# DirichletBCs assignment for coupled vector field problem

**URL:** <https://fenicsproject.discourse.group/t/dirichletbcs-assignment-for-coupled-vector-field-problem/9050>\
**Category:** dolfinx\
**Created:** [August 18, 2022, 9:24pm UTC](https://fenicsproject.discourse.group/t/dirichletbcs-assignment-for-coupled-vector-field-problem/9050 "2022-08-18T21:24:55Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![skatulla](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/skatulla/32/4089_2.png) [@skatulla](https://fenicsproject.discourse.group/u/skatulla)\
**Post date:** [August 18, 2022, 9:24pm UTC](https://fenicsproject.discourse.group/t/dirichletbcs-assignment-for-coupled-vector-field-problem/9050/1 "2022-08-18T21:24:56Z")

</div>

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:

```auto
# ----------------------------------------------------------------------------------------------------------------------
# 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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 18, 2022, 9:51pm UTC](https://fenicsproject.discourse.group/t/dirichletbcs-assignment-for-coupled-vector-field-problem/9050/2 "2022-08-18T21:51:06Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![skatulla](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/skatulla/32/4089_2.png) [@skatulla](https://fenicsproject.discourse.group/u/skatulla)\
**Post date:** [August 19, 2022, 7:50am UTC](https://fenicsproject.discourse.group/t/dirichletbcs-assignment-for-coupled-vector-field-problem/9050/3 "2022-08-19T07:50:49Z")

</div>

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](https://jorgensd.github.io/dolfinx-tutorial/chapter2/hyperelasticity.html):

```auto
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?):

```auto
(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:

```auto
(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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 19, 2022, 8:25am UTC](https://fenicsproject.discourse.group/t/dirichletbcs-assignment-for-coupled-vector-field-problem/9050/4 "2022-08-19T08:25:57Z")

</div>

> [@skatulla](#):
>
> since the assigned boundary condition is a constant value, i.e. the function would be equally a constant.

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`.

> [@skatulla](#):
>
> erelastic demo [https://jorgensd.github.io/dolfinx-tutorial/chapter2/hyperelasticity.html:](https://jorgensd.github.io/dolfinx-tutorial/chapter2/hyperelasticity.html:)
> 
> ```auto
> 
> ```

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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 19, 2022, 8:35am UTC](https://fenicsproject.discourse.group/t/dirichletbcs-assignment-for-coupled-vector-field-problem/9050/5 "2022-08-19T08:35:17Z")

</div>

> [@skatulla](#):
>
> 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?):
> 
> ```auto
> 
> ```

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.

---

<div class="post-metadata">

**Author:** ![skatulla](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/skatulla/32/4089_2.png) [@skatulla](https://fenicsproject.discourse.group/u/skatulla)\
**Post date:** [August 19, 2022, 10:02am UTC](https://fenicsproject.discourse.group/t/dirichletbcs-assignment-for-coupled-vector-field-problem/9050/6 "2022-08-19T10:02:40Z")

</div>

Ok, thanks. Got it working:

```auto
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)

```