DBCs on submesh - unexpected dof indices?

mesh
1

Dear Community,

For coupled problems (e.g. FSI), it can be beneficial to have a routine that removes Dirichlet BCs that are both set on a fluid and a solid mesh node (e.g. when two surfaces with DBCs meet at an edge), since the Lagrange multiplier field at these nodes does not have a proper solution.
Hence, my logic in doing so would be to go from the DBC dof indices - in my understanding local with respect to the submesh dofmap - to the global indices, intersect the arrays, remove common entries in one of them, and then go back to the local indices (w.r.t. the submesh).

However, I do not understand the logic of dof index oderings on the submeshes. The following MWE on the one hand selects the common dofs between two submeshes, and then tries to find these common dofs via the DBC dof indices of DBCs applied to the center facets of the main mesh - but at the submesh boundaries. (MWE with dolfinx nightly container from 6 Aug 2025)

#!/usr/bin/env python3
from dolfinx import mesh, fem
from mpi4py import MPI
import numpy as np

n = 2
# Demo assumes faces lie along x[0] == 0.5, so check n is even
assert n % 2 == 0
msh = mesh.create_unit_square(MPI.COMM_WORLD, n, n)
tdim = msh.topology.dim
fdim = tdim - 1

left_cells = mesh.locate_entities(msh, tdim, lambda x: x[0] <= 0.5)
left_sm, left_em, left_vm, left_gm = mesh.create_submesh(msh, tdim, left_cells)

right_cells = mesh.locate_entities(msh, tdim, lambda x: x[0] >= 0.5)
right_sm, right_em, right_vm, right_gm = mesh.create_submesh(msh, tdim, right_cells)

# function spaces and functions on the submeshes
V_l = fem.functionspace(left_sm, ("Lagrange", 1, (left_sm.geometry.dim,)))
f_l = fem.Function(V_l)
V_r = fem.functionspace(right_sm, ("Lagrange", 1, (right_sm.geometry.dim,)))
f_r = fem.Function(V_r)

centre_dofs_left = fem.locate_dofs_geometrical(V_l, lambda x: np.isclose(x[0], 0.5))
centre_dofs_right = fem.locate_dofs_geometrical(V_r, lambda x: np.isclose(x[0], 0.5))

dbcs_l, dbcs_r = [], []
dbcs_l.append( fem.dirichletbc(f_l, centre_dofs_left) )
dbcs_r.append( fem.dirichletbc(f_r, centre_dofs_right) )

# get dofs numbering with respect to main mesh - left side
num_vertices_l = V_l.dofmap.index_map.size_local + V_l.dofmap.index_map.num_ghosts
bs = V_l.dofmap.bs
dofs_left_ = np.empty(bs*num_vertices_l, dtype=np.int32)
for i, vert in enumerate(left_gm):
    for j in range(bs):
        dofs_left_[bs*i+j] = bs*vert+j

# get dofs numbering with respect to main mesh - right side
num_vertices_r = V_r.dofmap.index_map.size_local + V_r.dofmap.index_map.num_ghosts
bs = V_r.dofmap.bs
dofs_right_ = np.empty(bs*num_vertices_r, dtype=np.int32)
for i, vert in enumerate(right_gm):
    for j in range(bs):
        dofs_right_[bs*i+j] = bs*vert+j

common_dofs = np.intersect1d(dofs_left_,dofs_right_)

print("common_dofs btw. submeshes (=interface): ",common_dofs)

# get dofs of DBCs w.r.t. main mesh
dbc_dofs_l = dofs_left_[dbcs_l[0].dof_indices()[0]]
dbc_dofs_r = dofs_right_[dbcs_r[0].dof_indices()[0]]
print(dbc_dofs_l)
print(dbc_dofs_r)

common_dofs_dbcs = np.intersect1d(dbc_dofs_l,dbc_dofs_r)
# should yield the same set of dofs like common_dofs
print("common_dofs of DBCs: ",common_dofs_dbcs)

Here, I would expect that common_dofs_dbcs yields the same set of indices as common_dofs, right? But it does not. Would be great if someone has a hint on how dof orderings happen on submeshes.

Best,
Marc

2

I’m a bit confused as to how dofs that exist on two different sub-meshes, is expected to have the same indices.

Each sub-mesh is in itself a proper mesh. it’s only relation to the parent mesh is through the entity_map, which tells you what cell in the sub-mesh corresponds to what cell in the parent grid.
You can in addition make a map from each facet of the sub-mesh to each facet of the parent mesh (if the meshes are of co-dimension 0).

The way assembly of systems with sub-meshes work, is that you either assemble into a “block” or “nest” matrix, which concatenates the dof indices of the individual submesh spaces into a single, monolithic matrix.

If you want to remove degrees of freedom that are “co-located”, I would probably find intersections through the geometrical coordinates of the degrees of freedom, as the orientation of the dofs from the two sides of the interfaces doesn’t necessarily align, ref Fig. 1 in: https://dl.acm.org/doi/pdf/10.1145/3524456

3

Thanks for the swift reply. So I forgot to emphasise, I am using the fourth return from create_submesh, which precisely should be the ordered list of node indices of the submesh with respect to the global mesh. This seems to be also working for the first intersection, returning exactly the mutual dofs between the two submeshes.

The only thing that wonders me is that it does not work when going via the dof indices from the DBC routine.

Maybe I’m missing something here, I will also consider another approach.

4

So what I don’t understand is how

are supposed to be on the “main” mesh.

The bcs are made on the sub-meshes, thus the dof indices have no relation to the parent mesh.

If one wants everything to “live” on the parent mesh, where one reduce the size of the dofmap, as opposed to having new meshes (which is usually fine, as long as the spaces on the restrictions are H^1 and continuous, not DG, H-div or H-curl) then multiphenicsx would give you these kind of “common” relations to the parent mesh.

5

Yes, but in this part, I get the dof indices with respect to the main mesh, by looping over the (ordered) list from left_gm (right_gm). And by means of these, intersecting also gives me the correct common dofs between the meshes.

6

This has some very strong assumption on the doflayout and ordering of vertices, which I am not sure is generally true (for anything that is not P1 or in parallel).

Since your example is P1, where one indeed can find a relation between the dofs and vertices, I would have a look at: Moving submesh diverges from its designated trajectory in a two-phase problem - #2 by dokken
as I try my best to explain how to map data from parent to sub-mesh, through geometry, topology and dof indices.

7

Thanks for pointing to this. I looked at the example and I believe that it does something very similar to what I do - just vice versa.

Nonetheless, I boiled my MWE down to a part where I believe the main problem arises. Consider

#!/usr/bin/env python3
from dolfinx import mesh, fem
from mpi4py import MPI

n = 2
# Demo assumes faces lie along x[0] == 0.5, so check n is even
assert n % 2 == 0
msh = mesh.create_unit_square(MPI.COMM_WORLD, n, n)
tdim = msh.topology.dim
fdim = tdim - 1

left_cells = mesh.locate_entities(msh, tdim, lambda x: x[0] <= 0.5)
left_sm, left_em, left_vm, left_gm = mesh.create_submesh(msh, tdim, left_cells)

right_cells = mesh.locate_entities(msh, tdim, lambda x: x[0] >= 0.5)
right_sm, right_em, right_vm, right_gm = mesh.create_submesh(msh, tdim, right_cells)

# function spaces and functions on the main mesh and the submeshes
V_G = fem.functionspace(msh, ("Lagrange", 1))
V_l = fem.functionspace(left_sm, ("Lagrange", 1))
V_r = fem.functionspace(right_sm, ("Lagrange", 1))
f_G = fem.Function(V_G)
f_l = fem.Function(V_l)
f_r = fem.Function(V_r)

coG = f_G.function_space.tabulate_dof_coordinates()
col = f_l.function_space.tabulate_dof_coordinates()
cor = f_r.function_space.tabulate_dof_coordinates()
print("coords main:")
print(coG)
print("coords left:")
print(col)
print("coords right:")
print(cor)

print("Left node ids w.r.t. main mesh:")
print(left_gm)
print("Right node ids w.r.t. main mesh:")
print(right_gm)

If I’m not mistaken, the node maps (left_gm, right_gm) store the node ids with respect to the main mesh, right? Comparing the printouts with the dof coordinates, however, it appears that only for the right mesh, the indices are correct (relating cor to coG). For the left mesh, it seems that this mapping is not correct anymore (e.g. submesh index 0 is not the same as main mesh index 0, cf. relation of col to coG).

Could this point to a potential bug in the node maps, or can these maps simply not be interpreted as mappings from sub to main mesh node IDs?

My output of this is:

coords main:
[[ 5.00000000e-01 -1.07711404e-17  0.00000000e+00]
 [ 1.00000000e+00 -1.07711404e-17  0.00000000e+00]
 [ 1.00000000e+00  5.00000000e-01  0.00000000e+00]
 [ 5.00000000e-01  5.00000000e-01  0.00000000e+00]
 [ 1.07711404e-17  1.07711404e-17  0.00000000e+00]
 [ 1.00000000e+00  1.00000000e+00  0.00000000e+00]
 [ 1.07711404e-17  5.00000000e-01  0.00000000e+00]
 [ 5.00000000e-01  1.00000000e+00  0.00000000e+00]
 [-1.07711404e-17  1.00000000e+00  0.00000000e+00]]
coords left:
[[ 1.07711404e-17  1.07711404e-17  0.00000000e+00]
 [ 5.00000000e-01 -1.07711404e-17  0.00000000e+00]
 [ 5.00000000e-01  5.00000000e-01  0.00000000e+00]
 [ 1.07711404e-17  5.00000000e-01  0.00000000e+00]
 [ 5.00000000e-01  1.00000000e+00  0.00000000e+00]
 [-1.07711404e-17  1.00000000e+00  0.00000000e+00]]
coords right:
[[ 5.00000000e-01  1.07711404e-17  0.00000000e+00]
 [ 1.00000000e+00 -1.07711404e-17  0.00000000e+00]
 [ 1.00000000e+00  5.00000000e-01  0.00000000e+00]
 [ 5.00000000e-01  5.00000000e-01  0.00000000e+00]
 [ 1.00000000e+00  1.00000000e+00  0.00000000e+00]
 [ 5.00000000e-01  1.00000000e+00  0.00000000e+00]]
Left node ids w.r.t. main mesh:
[0 3 4 6 7 8]
Right node ids w.r.t. main mesh:
[0 1 2 3 5 7]

Entries 0 3 4 of left mesh are wrong.

8

Ok, I now have found the solution, which indeed was related to different orderings of dof and node coordinates on submeshes. The following is the correct version to the MWE from my opening post, using Jorgen’s vertex_to_dofmap function:

#!/usr/bin/env python3
from dolfinx import mesh, fem, cpp
from mpi4py import MPI
import numpy as np

def vertex_to_dof_map_vectorized(V):
    """Create a map from the vertices of the mesh to the corresponding degree of freedom."""
    msh = V.mesh
    num_vertices_per_cell = cpp.mesh.cell_num_entities(
        msh.topology.cell_type, 0
    )

    dof_layout2 = np.empty((num_vertices_per_cell,), dtype=np.int32)
    for i in range(num_vertices_per_cell):
        var = V.dofmap.dof_layout.entity_dofs(0, i)
        assert len(var) == 1
        dof_layout2[i] = var[0]

    num_vertices = (
        msh.topology.index_map(0).size_local + msh.topology.index_map(0).num_ghosts
    )

    c_to_v = msh.topology.connectivity(msh.topology.dim, 0)
    assert (c_to_v.offsets[1:] - c_to_v.offsets[:-1] == c_to_v.offsets[1]).all(), (
        "Single cell type supported"
    )

    vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
    vertex_to_dof_map[c_to_v.array] = V.dofmap.list[:, dof_layout2].reshape(-1)
    return vertex_to_dof_map

n = 2
# Demo assumes faces lie along x[0] == 0.5, so check n is even
assert n % 2 == 0
msh = mesh.create_unit_square(MPI.COMM_WORLD, n, n)
tdim = msh.topology.dim
fdim = tdim - 1

left_cells = mesh.locate_entities(msh, tdim, lambda x: x[0] <= 0.5)
left_sm, left_em, left_vm, left_gm = mesh.create_submesh(msh, tdim, left_cells)

right_cells = mesh.locate_entities(msh, tdim, lambda x: x[0] >= 0.5)
right_sm, right_em, right_vm, right_gm = mesh.create_submesh(msh, tdim, right_cells)

V_G = fem.functionspace(msh, ("Lagrange", 1))
f_G = fem.Function(V_G)

# function spaces and functions on the submeshes
V_l = fem.functionspace(left_sm, ("Lagrange", 1, (left_sm.geometry.dim,)))
V_r = fem.functionspace(right_sm, ("Lagrange", 1, (right_sm.geometry.dim,)))
f_l = fem.Function(V_l)
f_r = fem.Function(V_r)

vert_to_dof_l = vertex_to_dof_map_vectorized(V_l)
vert_to_dof_r = vertex_to_dof_map_vectorized(V_r)

centre_dofs_left = fem.locate_dofs_geometrical(V_l, lambda x: np.isclose(x[0], 0.5))
centre_dofs_right = fem.locate_dofs_geometrical(V_r, lambda x: np.isclose(x[0], 0.5))

dbcs_l, dbcs_r = [], []
dbcs_l.append( fem.dirichletbc(f_l, centre_dofs_left) )
dbcs_r.append( fem.dirichletbc(f_r, centre_dofs_right) )

# get dofs numbering with respect to main mesh - left side
num_vertices_l = V_l.dofmap.index_map.size_local + V_l.dofmap.index_map.num_ghosts
bs = V_l.dofmap.bs
dofs_left_ = np.empty(bs*num_vertices_l, dtype=np.int32)
map_l = np.empty(num_vertices_l, dtype=np.int32)
map_l[vert_to_dof_l] = left_gm
for i, vert in enumerate(map_l):
    for j in range(bs):
        dofs_left_[bs*i+j] = bs*vert+j

# get dofs numbering with respect to main mesh - right side
num_vertices_r = V_r.dofmap.index_map.size_local + V_r.dofmap.index_map.num_ghosts
bs = V_r.dofmap.bs
dofs_right_ = np.empty(bs*num_vertices_r, dtype=np.int32)
map_r = np.empty(num_vertices_r, dtype=np.int32)
map_r[vert_to_dof_r] = right_gm
for i, vert in enumerate(map_r):
    for j in range(bs):
        dofs_right_[bs*i+j] = bs*vert+j

common_dofs = np.intersect1d(dofs_left_,dofs_right_)

print("common_dofs btw. submeshes (=interface): ",common_dofs)

# get dofs of DBCs w.r.t. main mesh
dbc_dofs_l = dofs_left_[dbcs_l[0].dof_indices()[0]]
dbc_dofs_r = dofs_right_[dbcs_r[0].dof_indices()[0]]

print(dbcs_l[0].dof_indices()[0])
print("dofs_left_(g): ", dofs_left_)
print("dbc_dofs_l(g): ", dbc_dofs_l)
print("###")
print(dbcs_r[0].dof_indices()[0])
print("dofs_right_(g): ", dofs_right_)
print("dbc_dofs_r(g): ", dbc_dofs_r)

common_dofs_dbcs = np.intersect1d(dbc_dofs_l,dbc_dofs_r)
# should yield the same set of dofs like common_dofs
print("common_dofs of DBCs: ",common_dofs_dbcs)
9

@marchirschvogel you are getting closer, but you are still not quite right.
As you are using right_gm, you are talking about the physical nodes in mesh.geometry.x, not the vertices of mesh.topology.
You should introduce a third map, so that you can move between:
dofs<->vertices<->mesh nodes

Alternatively you can go from the sub_dofs<->sub_vertices<->parent_vertices<->parent_dofs
through either right_vm or left_vm.
Note the clear distinction in the output of the submesh.

Here is a code that does the way through the geometry

#!/usr/bin/env python3
from dolfinx import mesh, fem, cpp
import numpy as np
from mpi4py import MPI

n = 2
# Demo assumes faces lie along x[0] == 0.5, so check n is even
assert n % 2 == 0
msh = mesh.create_unit_square(MPI.COMM_WORLD, n, n)
tdim = msh.topology.dim
fdim = tdim - 1

left_cells = mesh.locate_entities(msh, tdim, lambda x: x[0] <= 0.5)
left_sm, left_em, left_vm, left_gm = mesh.create_submesh(msh, tdim, left_cells)

right_cells = mesh.locate_entities(msh, tdim, lambda x: x[0] >= 0.5)
right_sm, right_em, right_vm, right_gm = mesh.create_submesh(msh, tdim, right_cells)

# function spaces and functions on the main mesh and the submeshes
V_G = fem.functionspace(msh, ("Lagrange", 1))
V_l = fem.functionspace(left_sm, ("Lagrange", 1))
V_r = fem.functionspace(right_sm, ("Lagrange", 1))
f_G = fem.Function(V_G)
f_l = fem.Function(V_l)
f_r = fem.Function(V_r)

coG = f_G.function_space.tabulate_dof_coordinates()
col = f_l.function_space.tabulate_dof_coordinates()
cor = f_r.function_space.tabulate_dof_coordinates()


def vertex_to_dof_map_vectorized(V):
    mesh = V.mesh
    num_vertices_per_cell = cpp.mesh.cell_num_entities(mesh.topology.cell_type, 0)

    dof_layout2 = np.empty((num_vertices_per_cell,), dtype=np.int32)
    for i in range(num_vertices_per_cell):
        var = V.dofmap.dof_layout.entity_dofs(0, i)
        assert len(var) == 1
        dof_layout2[i] = var[0]

    num_vertices = (
        mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
    )

    c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
    assert (c_to_v.offsets[1:] - c_to_v.offsets[:-1] == c_to_v.offsets[1]).all(), (
        "Single cell type supported"
    )

    vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
    vertex_to_dof_map[c_to_v.array] = V.dofmap.list[:, dof_layout2].reshape(-1)
    return vertex_to_dof_map


# NOTE: Assume that there is a 1-1 relationship between dofs and vertices and nodes
G_v_to_d = vertex_to_dof_map_vectorized(V_G)
l_v_to_d = vertex_to_dof_map_vectorized(V_l)
r_v_to_d = vertex_to_dof_map_vectorized(V_r)
l_d_to_v = np.argsort(l_v_to_d)
r_d_to_v = np.argsort(r_v_to_d)
V_l.mesh.topology.create_connectivity(0, tdim)
l_v_to_n = cpp.mesh.entities_to_geometry(
    V_l.mesh._cpp_object, 0, np.arange(len(l_v_to_d), dtype=np.int32), False
)
V_r.mesh.topology.create_connectivity(0, tdim)
r_v_to_n = cpp.mesh.entities_to_geometry(
    V_r.mesh._cpp_object, 0, np.arange(len(r_v_to_d), dtype=np.int32), False
)

centre_dofs_left = fem.locate_dofs_geometrical(V_l, lambda x: np.isclose(x[0], 0.5))
centre_dofs_right = fem.locate_dofs_geometrical(V_r, lambda x: np.isclose(x[0], 0.5))

center_vertices_left = l_d_to_v[centre_dofs_left]
center_nodes_left = l_v_to_n[center_vertices_left]
parent_center_nodes_left = left_gm[center_nodes_left]

center_vertices_right = r_d_to_v[centre_dofs_right]
center_nodes_right = r_v_to_n[center_vertices_right]
parent_center_nodes_right = right_gm[center_nodes_right]
print(
    V_G.mesh.geometry.x[parent_center_nodes_left],
    V_G.mesh.geometry.x[parent_center_nodes_right],
)

# Map geometry nodes of parent mesh to dofs
V_G.mesh.topology.create_connectivity(0, tdim)
G_v_to_n = cpp.mesh.entities_to_geometry(
    V_G.mesh._cpp_object, 0, np.arange(len(G_v_to_d), dtype=np.int32), False
)
G_n_to_v = np.argsort(G_v_to_n.flatten())
xG_coords = V_G.tabulate_dof_coordinates()

print(xG_coords[G_v_to_d[G_n_to_v[parent_center_nodes_left]]])
10

Thanks for the swift reply. Ok, I think I understand the way to go through the different maps.
One question about the vertex_to_dof_map_vectorized function: if I run this with a natively quadratic mesh (where the higher order nodes are already defined in the input), I still get

print(msh.topology.cell_type)
print(num_vertices_per_cell)
CellType.quadrilateral
4

instead of

CellType.quadrilateral_9
9

Is that expected?

11

This is one of the core concepts of the dolfinx.mesh.Topology and dolfinx.mesh.Geometry split.

A topology will only encode vertices of cells, and the entities that can be made up by these (cells, facets, ridges/edges and peaks/vertices) and the connectivity between them.

The geometry will encode the curvature of a cell (if it is not affine, and contains it’s own index map for all the nodes that make up the mesh, and a coordinate element to move from an affine reference element to the curved element.
This split makes it possible to for instance share the topology between a linear and curved grid, as for instance done in: Adaptive mesh refinement with NetGen and DOLFINx — FEniCSx tutorial