Unexpected Behaviour with .sub in fenicsx

Hey there,
The following code generates a line mesh from (0,0) to (1,0) and solves a simple FEM-problem that can be interpreted as a simple bar of length, cross-section and Youngs Modulus = 1 where the point (0,0) is fixed and a force of magnitude 0.1 pulls on the right-hand side (at (1,0)) in positive x-direction.

import ufl
from dolfinx import io, fem, mesh
from mpi4py import MPI
from petsc4py import PETSc
import gmsh
import numpy as np

gmsh.initialize()

meshSize = 0.05
point1 = gmsh.model.geo.add_point(0.,0.,0.,meshSize)
point2 = gmsh.model.geo.add_point(1.,0.,0.,meshSize)

edge1 = gmsh.model.geo.add_line(point1,point2)

gmsh.model.geo.synchronize()

edges = [edge1]
model1 = gmsh.model.add_physical_group(1,edges)
gmsh.model.mesh.generate(1)

domain = io.ufl_mesh_from_gmsh(1,2)
geometry_data = io.extract_gmsh_geometry(gmsh.model)
topology_data = io.extract_gmsh_topology_and_markers(gmsh.model)
cells = topology_data[1]["topology"]
nodes = geometry_data[:,:2]

struct_mesh = mesh.create_mesh(MPI.COMM_WORLD, cells, nodes, domain)

Vu = fem.VectorFunctionSpace(struct_mesh,("CG",1))

du = ufl.TrialFunction(Vu)
u_ = ufl.TestFunction(Vu)

# bilinear form
aM = ufl.inner(ufl.grad(du),ufl.grad(u_))*ufl.dx

# "fixing" node at (0,0)
b_dofs_l = fem.locate_dofs_geometrical(Vu, lambda x : (np.isclose(x[0],0.)) & (np.isclose(x[1],0.)))
bc_left = fem.dirichletbc(fem.Constant(struct_mesh,np.array([0.,0.])),b_dofs_l,Vu)

# fixing y-dof at (1,0)
b_dofs_r = fem.locate_dofs_geometrical(Vu, lambda x : (np.isclose(x[0],1.) & (np.isclose(x[1],0.))))
bc_right = fem.dirichletbc(fem.Constant(struct_mesh,0.),b_dofs_r,Vu.sub(1))

bcs = [bc_left]#, bc_right]

# assemble lhs
A = fem.petsc.assemble_matrix(fem.form(aM),bcs=bcs)
A.assemble()

# create solver
solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setOperators(A)

# create empty rhs-vector
LM = ufl.inner(u_,fem.Constant(struct_mesh,np.array([0.,0.])))*ufl.dx
b = fem.petsc.create_vector(fem.form(LM))

# applying a force to the node at the very right, at (1,0)
xForce = 0.1
b.setValue(b_dofs_r[0]*2,xForce) # accessing x-component of node b_dofs_r
fem.petsc.set_bc(b,bcs=bcs)
b.assemble()

# create a function to write the solution into
u = fem.Function(Vu)
solver.solve(b,u.vector)

If I run the code like this, I get the expected result: The value of u increases linearly from 0. to .1 - exactly like the (trivial) analytical solution.

If I change the line bcs = [bc_left]#, bc_right] to bcs = [bc_left, bc_right], I expect that the boundary condition bc_right is applied too. This BC would constrain the y-dof at the end node ((1,0)). In this case this shouldn’t have an effect on the solution, but sets my engineering mind (which is concerned with static/kinematic determinacy considerations) to peace.
However, the result is not as expected at all. The u-function is constant (value 0.) between (0,0) and (0.5,0) and then increases linearly to 0.05 up to the point (1,0). It seems like instead of constraining the y-dof of (1,0), bc_right constrained the x-dof of (0.5,0).

How can this be explained?

PS: Another question: I noticed that b_dofs_r is 20 in this case. However, node 20 in the array nodes is at (0.95,0) - while node 1 (the second one) is at (1.,0.). Are the nodes reshuffled when the mesh is imported from gmsh to fenicsx?

You are not using the correct dofs in your boundary condition.
Consider the following modification:

from IPython import embed
import ufl
from dolfinx import io, fem, mesh
from dolfinx.io import gmshio
from mpi4py import MPI
from petsc4py import PETSc
import gmsh
import numpy as np

gmsh.initialize()

meshSize = 0.05
point1 = gmsh.model.geo.add_point(0., 0., 0., meshSize)
point2 = gmsh.model.geo.add_point(1., 0., 0., meshSize)

edge1 = gmsh.model.geo.add_line(point1, point2)

gmsh.model.geo.synchronize()

edges = [edge1]
model1 = gmsh.model.add_physical_group(1, edges)
gmsh.model.mesh.generate(1)


struct_mesh, _, _ = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)

Vu = fem.VectorFunctionSpace(struct_mesh, ("CG", 1))

du = ufl.TrialFunction(Vu)
u_ = ufl.TestFunction(Vu)

# bilinear form
aM = ufl.inner(ufl.grad(du), ufl.grad(u_))*ufl.dx

# "fixing" node at (0,0)
b_dofs_l = fem.locate_dofs_geometrical(Vu, lambda x: (
    np.isclose(x[0], 0.)) & (np.isclose(x[1], 0.)))
bc_left = fem.dirichletbc(fem.Constant(
    struct_mesh, np.array([0., 0.])), b_dofs_l, Vu)

# fixing y-dof at (1,0)
V_1, _ = Vu.sub(1).collapse()
b_dofs_r = fem.locate_dofs_geometrical((Vu.sub(1), V_1), lambda x: (
    np.isclose(x[0], 1.) & (np.isclose(x[1], 0.))))
bc_right = fem.dirichletbc(fem.Constant(
    struct_mesh, 0.), b_dofs_r[0], Vu.sub(1))

bcs = [bc_left, bc_right]

# assemble lhs
A = fem.petsc.assemble_matrix(fem.form(aM), bcs=bcs)
A.assemble()

# create solver
solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setOperators(A)

# create empty rhs-vector
LM = ufl.inner(u_, fem.Constant(struct_mesh, np.array([0., 0.])))*ufl.dx
b = fem.petsc.create_vector(fem.form(LM))

# applying a force to the node at the very right, at (1,0)
xForce = 0.1
b.setValue(b_dofs_r[1]*2, xForce)  # accessing x-component of node b_dofs_r
fem.petsc.set_bc(b, bcs=bcs)
b.assemble()

# create a function to write the solution into
u = fem.Function(Vu)
solver.solve(b, u.vector)
with io.XDMFFile(struct_mesh.comm, "u.xdmf", "w") as xdmf:
    xdmf.write_mesh(struct_mesh)
    xdmf.write_function(u)

Ah yes, like this it makes perfect sense! Thank you very much!