Physical labels incorrect in mesh generated with pygmsh and converted using meshio (gmsh or xdmf)

[cweickhmann]

mwe-files.zip Hi, here's a MWE for a problem generating and converting geometry/mesh using pygmsh and meshio. It is two problems really:

• Output of pygmsh.write being empty (gmsh format version 4.1) or wrong (gmsh format version 2.2)
• Meshio (potentially) jumbling up cell label order. I can post this as a separate issue in the meshio repo if you can confirm this.

What I am trying to do

I try to generate a geometry using pygmsh programatically, generate a mesh in a gmsh format (either 2.2 or 4.1) and try to convert it to XDMF using meshio.

Generate the geometry

#!/usr/bin/python3

import gmsh
import pygmsh as msh

w1, h1 = 10., 10.
w2, h2 = 5., 5.
meshsize = 2.

with msh.geo.Geometry() as geom:
    points = [
        [0., 0.],
        [-w1, 0.],
        [-w1, h1],
        [0., h1],
        [w2, 0.],
        [w2, h2],
        [0., h2],
    ]
    p = [geom.add_point(P[:2], mesh_size=meshsize) for P in points]

    n_pts = len(points)
    lines = [
        [0, 1],
        [1, 2],
        [2, 3],
        [3, 6],
        [6, 5],
        [5, 4],
        [4, 0],
        [6, 0],
    ]
    l = [geom.add_line(p[p1], p[p2]) for p1, p2 in lines]

    # Rect 1
    print("== Rect 1 "+"="*50)
    ls = [l[0], l[1], l[2], l[3], l[7]]
    rect1_cl = geom.add_curve_loop(ls)
    rect1_s = geom.add_plane_surface(rect1_cl)
    geom.add_physical(rect1_s, "1")

    # Rect 2
    print("== Rect 2 "+"="*50)
    ls = [-l[7], l[4], l[5], l[6]]
    rect2_cl = geom.add_curve_loop(ls)
    rect2_s = geom.add_plane_surface(rect2_cl)
    geom.add_physical(rect2_s, "2")

    # Make the two boundaries
    left = rect1_cl.curves[:-1]
    left = geom.add_physical(left, "10")
    right = rect2_cl.curves[1:]
    right = geom.add_physical(right, "11")

    mesh = geom.generate_mesh() # (verbose=True)
    name_root = "gmsh_mwe"

    # Variant 1: gmsh22
    mesh.write(name_root + ".pygmsh.22.msh", file_format="gmsh22")

    # Variant 2: gmsh 4.1
    # Results in WriteError: Specify entity information to deal with more than one cell type
    mesh.write(name_root + "pygmsh.41.msh", file_format="gmsh")

    # Variant 3: unroll geo file
    gmsh.write(name_root + ".geo_unrolled")

Running Variant 1, I obtain a msh file with the following output:

== Rect 1 ==================================================
== Rect 2 ==================================================
WARNING:root:Binary Gmsh needs 32-bit integers (got uint64). Converting.
WARNING:root:Binary Gmsh needs 32-bit integers (got uint64). Converting.
WARNING:root:Binary Gmsh needs 32-bit integers (got uint64). Converting.
WARNING:root:Appending zeros to replace the missing physical tag data.
WARNING:root:Appending zeros to replace the missing geometrical tag data.

Running Variant 2, I get a WriteError:

== Rect 1 ==================================================
== Rect 2 ==================================================
Traceback (most recent call last):
  File "gmsh_mwe.py", line 58, in <module>
    mesh.write(name_root + ".pygmsh.41.msh", file_format="gmsh")
  File "..../.local/lib/python3.8/site-packages/meshio/_mesh.py", line 178, in write
    write(path_or_buf, self, file_format, **kwargs)
  File "..../.local/lib/python3.8/site-packages/meshio/_helpers.py", line 147, in write
    return writer(filename, mesh, **kwargs)
  File "..../.local/lib/python3.8/site-packages/meshio/gmsh/main.py", line 111, in <lambda>
    "gmsh": lambda f, m, **kwargs: write(f, m, "4.1", **kwargs),
  File "..../.local/lib/python3.8/site-packages/meshio/gmsh/main.py", line 102, in write
    writer.write(filename, mesh, binary=binary, float_fmt=float_fmt)
  File "..../.local/lib/python3.8/site-packages/meshio/gmsh/_gmsh41.py", line 358, in write
    _write_nodes(fh, mesh.points, mesh.cells, mesh.point_data, float_fmt, binary)
  File "..../.local/lib/python3.8/site-packages/meshio/gmsh/_gmsh41.py", line 610, in _write_nodes
    raise WriteError(
meshio._exceptions.WriteError: Specify entity information to deal with more than one cell type

Using Variant 3, I obtain a geo file which I can run through gmsh on the command line:

gmsh -2 -format msh22 gmsh_mwe.geo_unrolled

or

gmsh -2 -format msh41 gmsh_mwe.geo_unrolled

Finally I convert

• gmsh_mwe.gmsh22.msh
• gmsh_mwe.gmsh41.msh
• gmsh_mwe.pygmsh.22.msh
• gmsh_mwe.pygmsh.41.msh to xdmf using meshio using the following function (courtesy dokken):

import meshio
import numpy

def create_entity_mesh(mesh, cell_type, prune_z=False,
                       remove_unused_points=False):
    """
    Given a meshio mesh, extract mesh and physical markers for a given entity.
    We assume that all unused points are at the end of the mesh.points
    (this happens when we use physical markers with pygmsh)
    """
    cells = mesh.get_cells_type(cell_type)
    try:
        # If mesh created with gmsh API it is simple to extract entity data
        cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    except KeyError:
        # If mehs created with pygmsh, we need to parse through cell sets and sort the data
        cell_entities = []
        cell_data = []
        cell_sets = mesh.cell_sets_dict
        for marker, set in cell_sets.items():
            for type, entities in set.items():
                if type == cell_type:
                    cell_entities.append(entities)
                    cell_data.append(np.full(len(entities), int(marker)))
        cell_entities = np.hstack(cell_entities)
        sorted = np.argsort(cell_entities)
        cell_data = np.hstack(cell_data)[sorted]
    if remove_unused_points:
        num_vertices = len(np.unique(cells.reshape(-1)))
        # We assume that the mesh has been created with physical tags,
        # then unused points are at the end of the array
        points = mesh.points[:num_vertices]
    else:
        points = mesh.points

    # Create output mesh
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells},
                           cell_data={"name_to_read": [cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

fn_root = "gmsh_mwe"
infile_mesh = fn_root + ".msh"
outfile_mesh = fn_root + "_mesh.xdmf"
outfile_boundary = fn_root + "_boundaries.xdmf"

inmsh = meshio.read(infile_mesh)

outmsh = create_entity_mesh(inmsh, "triangle", prune_z=True)
meshio.write(outfile_mesh, outmsh)

outboundary = create_entity_mesh(inmsh, "line", prune_z=True)
meshio.write(filename=outfile_boundary, mesh=outboundary)

What happens

Variant 1 / gmsh_mwe.pygmsh.22.msh

Converting the binary gmsh file to XDMF produces a mesh where the cell markers show values of 1e-38. Looks like a type conversion problem to me.

Variant 2 / gmsh_mwe.pygmsh.41.msh

Fails upon creation of the msh file with a WriteError (see above). The msh file is thus almost empty.

Variant 3a (gmsh 2.2 format via the terminal)

Successfully produces XDMF files with desired cell tags (screenshot below).

Variant 3b (gmsh 4.1 format via the terminal)

Produces cell tags with the correct value range, but all jumbled up (see screenshot, note the boundary mesh is offset to make the jumbling more easily visible - coordinates are fine!).

Versions

Machine/OS

• Architecture x86_64

• Windows 10 20H2 with a Ubuntu 20.04 in a WSL 1 (but reproduces same results on a pure Ubuntu 20.04 x86_64)

• Python Version 3.8.5

This was tested using the following packages:

• pygmsh 7.1.8
• meshio 4.3.10
• gmsh 4.8.0
• Paraview 5.6.2

Attachments

mwe-files.zip

[cweickhmann]

Looks like part 1 of my problem is related to #427.

[nschloe]

There have been some label fixes lately, perhaps this is fixed now.

[jkneer]

I still have this problem with pygmsh v7.1.17

extract_to_meshio() messes up the indices of physicals.