Non-matching coordinate ordering for CG1 functions defined on submeshes

[rsz7ncf24p]

I am using the docker container of v0.9.0.
DOLFINx version: 0.9.0 based on GIT commit: 4116cca8cf91f1a7e877d38039519349b3e08620

Hello,

I am used to rely on the 1-to-1 correspondance between the DOF coordinates of CG1 functions and mesh coordinates.

When extracting two submeshes from a mesh and defining function spaces on them, however, this relation seems to be uncertain. Consider the following MWE.

The function create_box creates a square mesh which is divided by a line in the middle from top to bottom.

from dolfinx import fem, mesh, io
import numpy as np
import basix
import gmsh
from mpi4py import MPI

def create_box(h):
    gmsh.initialize()
    mesh_comm = MPI.COMM_WORLD
    model_rank = 0
    if mesh_comm.rank == model_rank:
        gmsh.model.occ.addPoint(0, 0, 0, 1, 1)
        gmsh.model.occ.addPoint(0.5, 0, 0, 1, 2)
        gmsh.model.occ.addPoint(1, 0, 0, 1, 3)
        gmsh.model.occ.addPoint(1, 1, 0, 1, 4)
        gmsh.model.occ.addPoint(0.5, 1, 0, 1, 5)
        gmsh.model.occ.addPoint(0, 1, 0, 1, 6)
        gmsh.model.occ.addLine(1, 2, 1)
        gmsh.model.occ.addLine(2, 3, 2)
        gmsh.model.occ.addLine(3, 4, 3)
        gmsh.model.occ.addLine(4, 5, 4)
        gmsh.model.occ.addLine(5, 6, 5)
        gmsh.model.occ.addLine(6, 1, 6)
        gmsh.model.occ.addLine(2, 5, 7)
        gmsh.model.occ.addCurveLoop([1, 7, 5, 6], 1)
        gmsh.model.occ.addCurveLoop([2, 3, 4, 7], 2)
        gmsh.model.occ.addPlaneSurface([1], 1)
        gmsh.model.occ.addPlaneSurface([2], 2)
        gmsh.model.occ.synchronize()
        gmsh.model.addPhysicalGroup(dim=1, tags=[7], tag=1, name="interface")
        gmsh.model.addPhysicalGroup(dim=2, tags=[1], tag=1, name="left")
        gmsh.model.addPhysicalGroup(dim=2, tags=[2], tag=2, name="right")
        gmsh.model.mesh.setOrder(1)
        gmsh.option.setNumber("Mesh.MeshSizeMin", h)
        gmsh.option.setNumber("Mesh.MeshSizeMax", h)
        gmsh.model.mesh.generate(2)
    domain, cell_markers, facet_markers = io.gmshio.model_to_mesh(gmsh.model,
                                                               mesh_comm,
                                                               model_rank,
                                                               gdim=2)
    gmsh.finalize()
    return domain, cell_markers, facet_markers

Then, this domain is split into two submeshes and function spaces (CG1) are defined on them. The coordinates of the mesh (domain_X.geometry.x) and the DOF coordinates (V_X.tabulate_dof_coordinates()) are then compared. For the global and first submesh the 1-to-1 relation still holds, but it is violated for the second submesh. I added some printouts to visualize.

domain, cell_markers, facet_markers = create_box(0.5)
ddim = domain.geometry.dim
cell_name = domain.topology.cell_name()

P1 = basix.ufl.element("Lagrange", cell_name, 1)

# Define the subdomains
domain_1, entity_map_1, vertex_map_1, node_map_1 = mesh.create_submesh(domain, ddim, cell_markers.find(1))
domain_2, entity_map_2, vertex_map_2, node_map_2 = mesh.create_submesh(domain, ddim, cell_markers.find(2))
V = fem.functionspace(domain, P1)
V_1 = fem.functionspace(domain_1, P1)
V_2 = fem.functionspace(domain_2, P1)
np.set_printoptions(formatter={'float': lambda x: "{0:0.2f}".format(x)})
print("V DOF coords == domain mesh coords?")
print(np.isclose((V.tabulate_dof_coordinates() - domain.geometry.x), 0))
print()
print("V1 DOF coords == domain_1 mesh coords?")
print(np.isclose((V_1.tabulate_dof_coordinates() - domain_1.geometry.x), 0))
print()
print("V2 DOF coords == domain_2 mesh coords?")
print(np.isclose((V_2.tabulate_dof_coordinates() - domain_2.geometry.x), 0))
print("V2 DOF coords:")
print(V_2.tabulate_dof_coordinates())
print("domain_2 mesh coords:")
print(domain_2.geometry.x)

For other cases (different mesh sizes in create_box(h) for example), the first submesh might be fine and the second one not, but always one is fine and one is not fine).
Am I overlooking something or is this intended behavior? If so, why, and can I, in an inexpensive post-processing step, reorder the DOFs such that they are aligned with the mesh coordinates again?

Thank you.

[dokken]

@rsz7ncf24p

Am I overlooking something or is this intended behavior? If so, why, and can I, in an inexpensive post-processing step, reorder the DOFs such that they are aligned with the mesh coordinates again?

Thank you.

It is intended behavior as:

1. the function is designed to work with arbitrary order function spaces and meshes, you cannot rely on this 1-1 correspondence.
2. It is supposed to work in parallel (on data partitioned with mpi), which introduces a whole other set of considerations

Why do you rely on this ordering being the same?
There are several functions to map data from degrees of freedom to vertices, and even the mesh geometry nodes, as discussed in:

Application of point forces, mapping vertex indices to corresponding DOFs

Here is a vectorized version of the code above (credit to @dokken) def vertex_to_dof_map_original(mesh, V): num_vertices = ( mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts ) vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32) num_cells = ( mesh.topology.index_map(mesh.topology.dim).size_local + mesh.topology.index_map(mesh.topology.dim).num_ghosts ) c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0) …

How to apply a boundary condition to a range of DOFs?

I don’t have time to read that user manual. Assuming that what is in NSET describes the original nodes in the mesh you would like to constrain, you need to: Create a map from the input node index to a vertex on your process, this is simply a matter of inverting mesh.geometry.input_global_indices. Once you have a map from the global indices (in your case a map from 97, 96, 95, ... to what is local on the process), you remove all the results that are -1. You then map from geometry to topology …

(And to some extent in: Renumbered .msh imported data giving wrong subdomains - #4 by dokken)

I would start with mapping dofs in the submesh to the nodes of the geometry in the submesh, and then use the last output argument of ‘create_submesh` to map these geometry indices from the submesh to the parent mesh.
https://docs.fenicsproject.org/dolfinx/v0.9.0/python/generated/dolfinx.mesh.html#dolfinx.mesh.create_submesh

[rsz7ncf24p]

Thank you for the explanation. I was curious and will probably now adapt the existing code such that I do not rely on the mapping, which is a limitation anyways.

@dokken

Why do you rely on this ordering being the same?

Inherited burdens of existing code… For geometric computations and an extrapolation method. Of course, all these can be adapted.