Misunderstanding regarding locate_dofs_topological

Hello everyone,

I am currently having difficulties applying boundary conditions on a rectangle.

My question is quite simple: I don’t see why the nodes are not selected in the same way whether the mesh is made of triangular or quadrilateral elements. With the code shown below, I would expect that all nodes at x=0 would be considered, so there would be two nodes where u_{x}=0.1.

from ufl import VectorElement
from dolfinx.fem import Function, Constant, FunctionSpace, locate_dofs_topological, set_bc, Expression, dirichletbc
from dolfinx.mesh import locate_entities_boundary, create_unit_square, CellType
import numpy as np
from petsc4py.PETSc import ScalarType
from mpi4py import MPI

test = "tri"
if test == "quad":
    mesh = create_unit_square(MPI.COMM_WORLD, 1, 1, CellType.quadrilateral)
elif test == "tri":
    mesh = create_unit_square(MPI.COMM_WORLD, 1, 1)
U_e = VectorElement('CG', mesh.ufl_cell(), degree = 1)
V = FunctionSpace(mesh, U_e)
u = Function(V, name = "displacement")
fdim = mesh.topology.dim - 1

#Dirichlet BC
def l_boundary(x):
    return np.isclose(x[0], 0)

l_bc_facets = locate_entities_boundary(mesh, fdim, l_boundary)
final_loading = 1
l_magnitude = Constant(mesh, ScalarType(0.1))
loading = Constant(mesh, ScalarType(final_loading))
current_l_bc = Expression(l_magnitude*loading, V.element.interpolation_points())

l_bc_func = Function(V.sub(0).collapse()[0])
l_bc = dirichletbc(l_bc_func, locate_dofs_topological(V.sub(0), fdim, l_bc_facets))
l_bc_func.interpolate(current_l_bc)
load_steps = np.linspace(0, final_loading, 3)

for t in load_steps:
    loading.value = t
    l_bc_func.interpolate(current_l_bc)
    set_bc(u.vector, [l_bc])
    print("Displacement vector", u.vector.array)
    print("Loading parameter", loading.value)

For triangular element i get:

Displacement vector [0. 0. 0. 0. 0. 0. 0. 0.]
Loading parameter 0.0
Displacement vector [0.05 0.   0.   0.   0.   0.   0.   0.  ]
Loading parameter 0.5
Displacement vector [0.1 0.  0.  0.  0.  0.  0.  0. ]
Loading parameter 1.0

Whereas for CellType.quadrilateral, i get the correct answer:

Displacement vector [0. 0. 0. 0. 0. 0. 0. 0.]
Loading parameter 0.0
Displacement vector [0.05 0.   0.05 0.   0.   0.   0.   0.  ]
Loading parameter 0.5
Displacement vector [0.1 0.  0.1 0.  0.  0.  0.  0. ]
Loading parameter 1.0

Do you have a clue?

Best regards, Paul

You are not setting your boundary condition on the subspace correctly, i.e.

Should be:

V0, _ = V.sub(0).collapse()
l_bc_func = Function(V0)
l_bc = dirichletbc(l_bc_func, locate_dofs_topological(
    (V.sub(0), V0), fdim, l_bc_facets), V.sub(0))

Works perfectly, thank you!