Mesh topology and geometry (built-in vs. Gmsh)

dolfinx
1

Hi all,

I am trying to understand what the local ordering of cell entities (vertices, edges, faces, …) used by dolfinx is and how the mesh.topology can be related to the mesh.geometry.

Some cells are defined in cpp/dolfinx/io/cells.h:

"""
quadrilateral using local vertex order as in Basix

2--------3
|        |
|        |
0--------1
"""

However, this is the input-output module and I did not find a similar clear definition of the cells somewhere under cpp/dolfinx/mesh/.
After experimenting with a mesh using quadrilateral cells it seems that the local order of the vertices in dolfinx is

"""
dolfinx (built-in mesh)

1--------3
|        |
|        |
0--------2
"""

I can also confirm the above with the following MWE:

import tempfile
import gmsh
import dolfinx
from dolfinx.io import gmshio
from mpi4py import MPI


def create_rectangle_grid(
    xmin,
    xmax,
    ymin,
    ymax,
    z=0.0,
    lc=0.1,
    num_cells=None,
    recombine=False,
    out_file=None,
):
    """TODO docstring"""
    gmsh.initialize()
    gmsh.model.add("rectangle")

    p0 = gmsh.model.geo.addPoint(xmin, ymin, z, lc)
    p1 = gmsh.model.geo.addPoint(xmax, ymin, z, lc)
    p2 = gmsh.model.geo.addPoint(xmax, ymax, z, lc)
    p3 = gmsh.model.geo.addPoint(xmin, ymax, z, lc)

    l0 = gmsh.model.geo.addLine(p0, p1)
    l1 = gmsh.model.geo.addLine(p1, p2)
    l2 = gmsh.model.geo.addLine(p2, p3)
    l3 = gmsh.model.geo.addLine(p3, p0)

    curve_loop = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3])
    surface = gmsh.model.geo.addPlaneSurface([curve_loop])

    if num_cells is not None:
        try:
            nx, ny = num_cells
        except TypeError:
            nx = ny = num_cells
        gmsh.model.geo.mesh.setTransfiniteCurve(l0, nx + 1)
        gmsh.model.geo.mesh.setTransfiniteCurve(l2, nx + 1)
        gmsh.model.geo.mesh.setTransfiniteCurve(l1, ny + 1)
        gmsh.model.geo.mesh.setTransfiniteCurve(l3, ny + 1)

        gmsh.model.geo.mesh.setTransfiniteSurface(surface, "Left")

        if recombine:
            # setRecombine(dim, tag, angle=45.0)
            gmsh.model.geo.mesh.setRecombine(2, surface)

    gmsh.model.geo.synchronize()
    gmsh.model.addPhysicalGroup(2, [surface])

    filepath = out_file or "./rectangle.msh"
    gmsh.model.mesh.generate(2)
    gmsh.write(filepath)
    gmsh.finalize()


if __name__ == "__main__":
    with tempfile.NamedTemporaryFile(suffix=".msh") as tf:
        create_rectangle_grid(
            0.0, 1.0, 0.0, 1.0, num_cells=(2, 1), recombine=True, out_file=tf.name
        )
        gmsh_grid, _, _ = gmshio.read_from_msh(tf.name, MPI.COMM_WORLD, gdim=2)


    df_grid = dolfinx.mesh.create_unit_square(
            MPI.COMM_WORLD, 2, 1, dolfinx.mesh.CellType.quadrilateral
            )


    def do_test(grid):
        verts = grid.topology.connectivity(2, 0).links(0)
        print(f"{verts=}")
        x = dolfinx.mesh.compute_midpoints(grid, 0, verts)
        print(f"{x=}")
        print(f"{grid.geometry.x=}")


    do_test(gmsh_grid)
    do_test(df_grid)

Output:

verts=array([0, 1, 2, 3], dtype=int32)
x=array([[0. , 0. , 0. ],
       [0.5, 0. , 0. ],
       [0. , 1. , 0. ],
       [0.5, 1. , 0. ]])
grid.geometry.x=array([[0. , 0. , 0. ],
       [0.5, 0. , 0. ],
       [0. , 1. , 0. ],
       [0.5, 1. , 0. ],
       [1. , 0. , 0. ],
       [1. , 1. , 0. ]])
verts=array([0, 1, 2, 3], dtype=int32)
x=array([[0. , 0. , 0. ],
       [0. , 1. , 0. ],
       [0.5, 0. , 0. ],
       [0.5, 1. , 0. ]])
grid.geometry.x=array([[0. , 0. , 0. ],
       [0. , 1. , 0. ],
       [0.5, 0. , 0. ],
       [0.5, 1. , 0. ],
       [1. , 0. , 0. ],
       [1. , 1. , 0. ]])

Questions/Remarks:

  1. The second and third vertex (indices 1 and 2) are swapped as the ordering of the quadrilateral in /cpp/dolfinx/io/cells.h suggests. However, as a user I would have expected the same mesh whether this was constructed with dolfinx or with Gmsh (why does gmshio.cells_perm_array(dolfinx.mesh.CellType.quadrilateral, 4) return np.array([0, 1, 3, 2]) and not np.array([0, 3, 1, 2])?)
  2. I was of the impression that in the above code grid.geometry.x[verts] should give the correct coordinates corresponding to the global vertex indices, that is the vertex tags for the topology and geometry are the same. Is this correct or not true in general?
  3. There is also grid.geometry.dofmap, but df_grid.geometry.x[df_grid.geometry.dofmap.links(0)] and gmsh_grid.geometry.x[gmsh_grid.geometry.dofmap.links(0)] give the same result as above with vertices 2 and 3 swapped.

I feel I am overlooking something.

Thanks,
Philipp

2

It does not really matter if the ordering is:

or

It just means that it is a tensor product representation of the cell (as opposed to going counter clockwise).

The geometry and topology of a mesh has separate ordering,as the topology is always first order, only containing vertices, while the geometry can contain more nodes, as the mesh can be higher order.
The function entities_to_geometry creates the map from the topology to the geometry (for linear/first order meshes): https://github.com/FEniCS/dolfinx/blob/f2b660a7119345bae7c136c775323f10e7903d96/python/dolfinx/wrappers/mesh.cpp#L427-L440
See use-case: https://github.com/FEniCS/dolfinx/blob/04110c14ac58e5a9d556b229c4c96e50ef9f0da9/python/test/unit/geometry/test_bounding_box_tree.py#L57-L63

3

thank you @dokken for the clarification :slight_smile: