4
$\begingroup$

With the move away from launchpad I hope this is the right place for this question.

Is there a way to initialize the DOFs for a function? generated by

u = Function(V)

When I run a code which I will put in at the end I get the fenics error:

    *** 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** 

What is interesting about this error is that it does not occur when the mesh is generated by UnitCubeMesh() but when the same mesh is generated manually as bellow then it occurs (This is a example so that I didnt have to post my large main code, where I need to use mesh editor). It also runs fine on a single core but when using more than one the error occurs. Eg by

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)

Any ideas as to how to get this to work?

Thanks For any help

UPDATE: Thanks for your help Garth

With that sorted today I have been trying to get the next bit of my code to run in parallel. What this block does is takes coordinates from the mesh, deforms them by the u vector then recombines them into a new mesh (it adds some more elements as well). It works fine in serial but I cant work out how to link the thread specific dofs and interpolated coordinates back into a single vector to manually create the new mesh. Is this even doable, or is there a much better way of updating mesh coordinates in parallel?

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)

Again thanks for your help, hopefully this is trivial as well. Sorry about the questions, this is the first time I have ever done anything in parallel.

$\endgroup$
6
$\begingroup$

Your sample code will not work in parallel because you are creating the same vertices and cells on all processes.

If you add mesh entities only on process 0 and then call

MeshPartitioning.build_distributed_mesh(mesh)

on all processes, your sample code runs fine.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.