DirichletBCs assignment for coupled vector field problem

dolfinx
1

Hi
I am trying to expand on the 3D hyperelastic demo problem but considering 2 coupled unknown vector fields, displacement vector u and orientation vector wf. I am currently stuck assigning the Dirichlet boundary conditions for both:

# ----------------------------------------------------------------------------------------------------------------------
# Create a box mesh of a beam
L = 20.0
H = 1.0
msh = mesh.create_box(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]), np.array([L, H, H])], [10, 1, 1],
                      mesh.CellType.hexahedron)
P1 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
P2 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
VY = fem.FunctionSpace(msh, ufl.MixedElement([P1, P2]))

# ----------------------------------------------------------------------------------------------------------------------
# Boundary conditions
def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], L)

# Dirichlet boundary application to affected dofs
u_left = np.array([0.0, 0.0, 0.0], dtype=PETSc.ScalarType)
u_right = np.array([0.0, 0.0, 1.0e-2], dtype=PETSc.ScalarType)
wf_left = np.array([0.0, 0.0, 0.0], dtype=PETSc.ScalarType)
wf_right = np.array([0.0, 0.0, 0.0], dtype=PETSc.ScalarType)

V0, submap = VY.sub(0).collapse()
left_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), left)
right_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), right)
bc1 = fem.dirichletbc(u_left, left_u_dofs[0], VY.sub(0))
bc2 = fem.dirichletbc(u_right, right_u_dofs[0], VY.sub(0))

V1, submap = VY.sub(1).collapse()
left_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), left)
right_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), right)
bc3 = fem.dirichletbc(wf_left, left_wf_dofs[0], VY.sub(1))
bc4 = fem.dirichletbc(wf_right, right_wf_dofs[0], VY.sub(1))

bcs = [bc1, bc2, bc3, bc4]

I am getting the following error:

Blockquote
RuntimeError: Creating a DirichletBC using a Constant is not supported when the Constant size is not equal to the block size of the constrained (sub-)space. Use a Function to create the DirichletBC.

Please advise.

2

The error message is quite clear.
Instead of using constants, you need to use a dolfinx.fem.Function to assign boundary conditions to your problem.
You should interpolate the constant value into the vector space.

3

Hi Dokken

Many thanks for your prompt response.

However, your suggestion seems to be a bit counter-intuitive, since the assigned boundary condition is a constant value, i.e. the function would be equally a constant. Besides, the constant BC assignment does work in the hyperelastic demo https://jorgensd.github.io/dolfinx-tutorial/chapter2/hyperelasticity.html:

u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

What also strikes me as odd, is the return value of the function fem.locate_dofs_topological in the hyperelastic demo which is a list of only 4 entries (the nodes?):

(Pdb) l
 63  	# Dirichlet boundary application to affected dofs
 64  	u_left = np.array([0.0, 0.0, 0.0], dtype=PETSc.ScalarType)
 65  	left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, np.where(facet_tag.values == 1))
 66  	bc1 = fem.dirichletbc(u_left, left_dofs, V)
 67  	pdb.set_trace()
 68  ->	u_right = np.array([0.0, 0.0, 1.0e-2], dtype=PETSc.ScalarType)
 69  	right_dofs = fem.locate_dofs_topological(V, facet_tag.dim, np.where(facet_tag.values == 2))
 70  	bc2 = fem.dirichletbc(u_right, right_dofs, V)
 71  	
 72  	bcs = [bc1, bc2]
 73  	
(Pdb) print (left_dofs)
[0 1 2 3]

whereas my equivalent two-field problem script gives a list of 12 entries:

(Pdb) l
 52  	wf_right = np.array([0.0, 0.0, 0.0], dtype=PETSc.ScalarType)
 53  	
 54  	V0, submap = VY.sub(0).collapse()
 55  	left_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), left)
 56  	pdb.set_trace()
 57  ->	right_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), right)
 58  	bc1 = fem.dirichletbc(u_left, left_u_dofs[0], VY.sub(0))
 59  	bc2 = fem.dirichletbc(u_right, right_u_dofs[0], VY.sub(0))
 60  	
 61  	V1, submap = VY.sub(1).collapse()
 62  	left_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), left)
(Pdb) print (left_u_dofs)
[array([ 0,  1,  2,  6,  7,  8, 12, 13, 14, 18, 19, 20], dtype=int32), array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11], dtype=int32)]

In both cases the Dirichlet boundary condition for the displacement field, u, is applied to a single quadrilateral element with four nodes x 3 displacement DOF = 12 displacement DOF for the element. I would have thought that the latter should be the return value (i.e. al list of 12 DOF) and the same for the single field and the two-field problem, respectively.

Could you please explain why the constant BC assignment works for the demo and what is the reason for DOF return value discrepancy between the demo and my own script. Thanks.

4

I implemented the constant bc in DOLFINx, and I know it’s limitations. The problem is not that the boundary condition is not spatially variying, but that you are supplying a fem.Constant and not a fem.Function.

It works in the hyperelasticity demo because you are not using a mixed space.

5

This is because dolfinx uses the fact that a vector space (when not part of a mixed space) has a block size, I.e. there are multiple dofs at the same spatial location.

6

Ok, thanks. Got it working:

def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], L)

def fixed_displacement_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def incremented_displacement_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.full(x.shape[1], 1.0)))

# Dirichlet boundary application to affected displacement dofs
V0, submap = VY.sub(0).collapse()
fixed_displacement = fem.Function(V0)
fixed_displacement.interpolate(fixed_displacement_expression)
incremented_displacement = fem.Function(V0)
incremented_displacement.interpolate(incremented_displacement_expression)
left_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), left)
right_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), right)
bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, V0)
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs,V0)