Fenics,初始化函数的自由度向量

计算科学 Python 芬尼克斯
2021-12-04 12:19:43

随着离开启动板,我希望这是解决这个问题的正确地方。

有没有办法为函数初始化自由度?由产生

u = 函数(V)

当我运行我将在最后输入的代码时,我得到了 fenics 错误:

    *** Error:   Unable to initialize vector of degrees of freedom for function.
***Reason:  Cannot be created from subspace. Consider collapsing the function space. 
****Where:   This error was encountered inside Function.cpp. 
****Process: 0** 

这个错误的有趣之处在于它不会在 UnitCubeMesh() 生成网格时发生,但是当手动生成相同的网格时,它会发生(这是一个示例,因此我不必发布我的大型主代码,我需要使用网格编辑器)。它在单个内核上也可以正常运行,但是当使用多个内核时会发生错误。例如通过

mpirun -n 4 python demo_hyperelasticity.py

    from dolfin import *


mesh = Mesh();
editor = MeshEditor();
editor.open(mesh, 'tetrahedron', 3, 3)
editor.init_vertices(8)
editor.init_cells(6)

editor.add_vertex(0, 0, 0, 0)
editor.add_vertex(1, 1, 0, 0)
editor.add_vertex(2, 0, 1, 0)
editor.add_vertex(3, 1, 1, 0)
editor.add_vertex(4, 0, 0, 1)
editor.add_vertex(5, 1, 0, 1)
editor.add_vertex(6, 0, 1, 1)
editor.add_vertex(7, 1, 1, 1)


editor.add_cell(0, 0, 1, 2, 4)
editor.add_cell(1, 1, 2, 3, 5)
editor.add_cell(2, 1, 2, 4, 5)

editor.add_cell(3, 4, 5, 6, 2)
editor.add_cell(4, 5, 6, 7, 3)
editor.add_cell(5, 5, 6, 3, 2)


editor.close()

#mesh = UnitCubeMesh(1, 1, 1)

#plot(mesh1, interactive=True)
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# Create mesh and define function space
V = VectorFunctionSpace(mesh, "Lagrange", 1)


u  = Function(V) ###### ERROR HERE#####

def boundary_Z0(x):
    tol = 1E-15
    return x[2] < tol
u0 = Expression(('0','0','0'))
bc = DirichletBC(V, u0, boundary_Z0)
# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function

B  = Constant((0.0, -0.5, 0.0))  # Body force per unit volume
T  = Constant((0.0,  0.0, 0.0))  # Traction force on the boundary

# Kinematics
I = Identity(V.cell().d)    # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Solve variational problem
solve(F == 0, u, bc, J=J,
      form_compiler_parameters=ffc_options)

# Plot and hold solution
plot(u, mode = "displacement", interactive = True)

关于如何让它发挥作用的任何想法?

谢谢你的帮助

更新:感谢您的帮助加思

今天排序完毕后,我一直试图让我的下一部分代码并行运行。这个块所做的是从网格中获取坐标,通过 u 向量对其进行变形,然后将它们重新组合成一个新的网格(它还添加了更多元素)。它在串行中工作正常,但我无法弄清楚如何将线程特定的自由度和插值坐标链接回单个向量以手动创建新网格。这甚至是可行的,还是有更好的方法来并行更新网格坐标?

coor_int = interpolate(Expression(("x[0]", "x[1]", "x[2]"), value_shape=3), V).vector().array()
    #turing displacement into numpy array
u_vec = u.vector().array()
    #calc new coordinates by subtaracting the displacement at each interploated
    #node from the oringonal mesh
new_coor = coor_int+u_vec

    #initalising matries
dofs_to_vert = np.zeros(num_vert, dtype=np.uintp)
vectordofs_to_vert = np.zeros((num_vert*3), dtype=np.uintp)
vectordofs_to_subvert = np.zeros((num_vert*3),dtype=np.uintp)
cellinds = np.zeros((num_cells,4), dtype=np.uintp)
    #creating degrees of freedom maps
dm = VV.dofmap()
dms = [V.sub(i).dofmap() for i in range(3)]

for cell in cells(mesh):
    cell_ind = cell.index()
    vert_inds = cell.entities(0)
    cellinds[cell_ind,:] = vert_inds #retuns the local indacies for each cell
    dofs_to_vert[dm.cell_dofs(cell_ind)] = vert_inds #dofs to vertcies map
    for i, (dms_i, dmcs_i) in enumerate(zip(dms, dms)):
        vectordofs_to_vert[dms_i.cell_dofs(cell_ind)] = vert_inds #gives map for coords to cells, not sencitive to xyz position
        vectordofs_to_subvert[dms_i.cell_dofs(cell_ind)] = i #gives map for each x,y,z to to particular cells for map above
#            #initalise more arrays       
map_mat = np.zeros((3*num_vert,3), dtype=float)
coorxyz = np.zeros((num_vert,3), dtype=float)
#    #make matrix of maps above
map_mat[:,0] = vectordofs_to_vert
map_mat[:,1] = vectordofs_to_subvert
map_mat[:,2] = new_coor

    #sort entries using maps matrix above
for ij in range(0,len(map_mat)):
            if map_mat[ij,1] == 0:
                coorxyz[map_mat[ij,0],0] = map_mat[ij,2]
            if map_mat[ij,1] == 1:
                coorxyz[map_mat[ij,0],1] = map_mat[ij,2]
            if map_mat[ij,1] == 2:
                coorxyz[map_mat[ij,0],2] = map_mat[ij,2] 


mesh1 = Mesh()

me = MeshEditor()
me.open(mesh1,'tetrahedron', 3, 3)
me.init_vertices(10)
me.init_cells(9)

for ii in range(0, len(coorxyz)):
    me.add_vertex(ii, coorxyz[ii,:])
for jj in range(0, len(cellinds)):
    me.add_cell(jj, cellinds[jj,:])

me.add_vertex(8, 2, 1, 0)
me.add_vertex(9, 2, 0, 0)

me.add_cell(6, 1, 3, 8, 5)
me.add_cell(7, 3, 8, 7, 5)
me.add_cell(8, 1, 5, 9, 8)

me.close()

plot(mesh1, interactive=True)

再次感谢您的帮助,希望这也是微不足道的。抱歉这些问题,这是我第一次并行做任何事情。

1个回答

您的示例代码不会并行工作,因为您在所有进程上创建相同的顶点和单元格。

如果您仅在进程 0 上添加网格实体,然后调用

MeshPartitioning.build_distributed_mesh(mesh)

在所有进程中,您的示例代码运行良好。