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

**URL:** <https://fenicsproject.discourse.group/t/mesh-topology-and-geometry-built-in-vs-gmsh/9492>\
**Category:** dolfinx\
**Created:** [October 12, 2022, 1:55pm UTC](https://fenicsproject.discourse.group/t/mesh-topology-and-geometry-built-in-vs-gmsh/9492 "2022-10-12T13:55:17Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![pdiercks](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/pdiercks/32/225_2.png) [@pdiercks](https://fenicsproject.discourse.group/u/pdiercks)\
**Post date:** [October 12, 2022, 1:55pm UTC](https://fenicsproject.discourse.group/t/mesh-topology-and-geometry-built-in-vs-gmsh/9492/1 "2022-10-12T13:55:17Z")

</div>

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`:

```python
"""
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

```python
"""
dolfinx (built-in mesh)

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

```

I can also confirm the above with the following MWE:

```python
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:

```auto
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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 13, 2022, 10:29am UTC](https://fenicsproject.discourse.group/t/mesh-topology-and-geometry-built-in-vs-gmsh/9492/2 "2022-10-13T10:29:17Z")

</div>

> [@pdiercks](#):
>
> 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])`?)

It does not really matter if the ordering is:

> [@pdiercks](#):
>
> ```auto
> 2--------3
> | |
> | |
> 0--------1
> 
> ```

or

> [@pdiercks](#):
>
> ```auto
> 1--------3
> | |
> | |
> 0--------2
> 
> ```

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](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](https://github.com/FEniCS/dolfinx/blob/04110c14ac58e5a9d556b229c4c96e50ef9f0da9/python/test/unit/geometry/test_bounding_box_tree.py#L57-L63)

---

<div class="post-metadata">

**Author:** ![pdiercks](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/pdiercks/32/225_2.png) [@pdiercks](https://fenicsproject.discourse.group/u/pdiercks)\
**Post date:** [October 13, 2022, 11:10am UTC](https://fenicsproject.discourse.group/t/mesh-topology-and-geometry-built-in-vs-gmsh/9492/3 "2022-10-13T11:10:19Z")

</div>

thank you @dokken for the clarification 🙂