[BUG]: gmshio causes kernel crash if entities are tagged several times

[pierricmora]

Summarize the issue

It seems to me that gmshio requires that all facets are assigned into a single physical group. When facets get assigned several tags, I get a kernel crash in a jupyter lab.

I tracked the line causing kernel crash until create_mesh in model_to_mesh.

How to reproduce the bug

The following 2D example represents an ellipse inclusion inside a rectangular matrix. Another smaller rectangle represents a region where I would like to assign a source term, without being linked to different physical properties. The source region overlaps the inclusion; therefore the natural way to define tags is:

• (run example with mode=1) 1 for the matrix (in & out the source), 2 for the inclusion (in & out the source), and 3 for the source (matrix & inclusion). -> kernel crash

As for now, I can bypass the bug by working with more cumbersome tags:

• (run example with mode=2) 10 for the matrix inside the source, 11 for the matrix outside the source, 20 and 21 for the inclusion on & out the source.

Minimal Example (Python)

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

# inspired from https://github.com/FEniCS/dolfinx/blob/5d15d263f20a9d8ada051dd2ece0f35894d38d4a/python/demo/demo_gmsh.py
def gmsh_rect_ellipse(mode: int=1, show: bool=False) -> gmsh.model:
    """Create a Gmsh model of an ellipse inclusion inside a rectangle,
    with another smaller rectangle representing a source term

    Args:
        show: if True: pop the Gmsh viewer

    Returns:
        Gmsh model
    """
    model = gmsh.model()

    name = "rectangle-ellipse"
    model.add(name)
    model.setCurrent(name)

    # Phase 1; "the matrix"
    x, y, dx, dy = 0, 0, 2, 1
    rec1 = gmsh.model.occ.addRectangle(x, y, 0, dx, dy)

    # Phase 2; "the inclusion"
    x, y, big_radius, small_radius = 1, 0.5, 0.5, 0.25
    ell1 = gmsh.model.occ.addEllipse(x, y, 0, big_radius, small_radius)
    cl_ell1 = gmsh.model.occ.addCurveLoop([ell1])
    srf_ell1 = gmsh.model.occ.addPlaneSurface([cl_ell1])

    # Virtual rectangle that represents a source term; not a physical property
    x, y, dx, dy = 0.25, 0.25, 0.5, 0.5
    rec2 = gmsh.model.occ.addRectangle(x, y, 0, dx, dy)

    # Fragment: the ellipse inclusion gets split
    ov, ovv = gmsh.model.occ.fragment([(2, rec1)], [(2, rec2)] + [(2, srf_ell1)])
    model.occ.synchronize()

    gmsh.option.setNumber("Mesh.RecombineAll", 0)
    gmsh.option.setNumber("Mesh.Algorithm", 6)

    # #######
    # Find the labels before assigning physical groups
    srce = [e[1] for e in ovv[1]]  # source
    phs2 = [e[1] for e in ovv[2]]  # phase 2
    phs1 = [e[1] for e in ov if not(e[1] in phs2)]  # phase 1

    if mode == 1:  # kernel crash; desired way of defining tags
        model.add_physical_group(2, phs1, 1)
        model.add_physical_group(2, phs2, 2)
        model.add_physical_group(2, srce, 3)
        model.setPhysicalName(2, 1, "Phase 1")
        model.setPhysicalName(2, 2, "Phase 2")
        model.setPhysicalName(2, 3, "Source")

    else:  # works, but not very convenient
        phs1_in_srce = [e for e in phs1 if e in srce]
        phs1_out_srce = [e for e in phs1 if not(e in srce)]
        phs2_in_srce = [e for e in phs2 if e in srce]
        phs2_out_srce = [e for e in phs2 if not(e in srce)]

        model.add_physical_group(2, phs1_in_srce, 10)
        model.add_physical_group(2, phs1_out_srce, 11)
        model.add_physical_group(2, phs2_in_srce, 20)
        model.add_physical_group(2, phs2_out_srce, 21)
        model.setPhysicalName(2, 10, "Phase 1 - inside source")
        model.setPhysicalName(2, 11, "Phase 1 - outside source")
        model.setPhysicalName(2, 20, "Phase 2 - inside source")
        model.setPhysicalName(2, 21, "Phase 2 - outside source")

    model.mesh.generate(2)

    if show:
        gmsh.fltk.run()

    return model


gmsh.initialize()

# ##################
# First, this works:
# Create model
model = gmsh_rect_ellipse(mode=2, show=False)

# Import
msh, ct, ft = gmshio.model_to_mesh(model, MPI.COMM_WORLD, rank=0)


# ##################
# Then, this does not:
# Create model
model = gmsh_rect_ellipse(mode=1, show=False)

# Import
msh, ct, ft = gmshio.model_to_mesh(model, MPI.COMM_WORLD, rank=0)

Output (Python)

No response

Version

0.8.0

DOLFINx git commit

5a20e2b

Installation

Docker image in a ubuntu computer

Additional information

No response

[jorgensd]

@pierricmora I don’t think we can support the behavior that you are suggesting without adding some very expensive checks to the gmshio readers. One would have to create the intersection of Physical Surfaces tagged by Gmsh.

You can do this by modifiying: https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/io/gmshio.py#L278 to be the unique set of rows in cells.
However, i am not certain the behavior you are suggesting would work for MeshTags of cells, https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/io/gmshio.py#L297
as there is an implicit assumption that a MeshTags object should have a unique set of indices, which is not true in your overlapping case (

dolfinx/cpp/dolfinx/mesh/MeshTags.h

Line 35 in 8e4122b

/// @brief Create a MeshTag from entities of given dimension on a

)

[pierricmora]

I think that the problem comes from the fact that the cells are read based on the physical groups, rather than by directly asking gmsh to return them. I mean: in extract_topology_and_markers, the for loop first gets the dimensions and tags in model.getPhysicalGroups(), and then only look for cells that share the tag with (entity_types, entity_tags, entity_topologies) = model.mesh.getElements(dim, tag=entity). So if a cell is tagged twice, it is returned twice. But also, if a cell is not tagged, then it is never returned.

Could the for loop be reversed, by first reading all cells regardless of tags with model.mesh.getElements(dim, tag=-1), and then looking for markers?

[jorgensd]

I think that the problem comes from the fact that the cells are read based on the physical groups, rather than by directly asking gmsh to return them. I mean: in extract_topology_and_markers, the for loop first gets the dimensions and tags in model.getPhysicalGroups(), and then only look for cells that share the tag with (entity_types, entity_tags, entity_topologies) = model.mesh.getElements(dim, tag=entity). So if a cell is tagged twice, it is returned twice. But also, if a cell is not tagged, then it is never returned.

Could the for loop be reversed, by first reading all cells regardless of tags with model.mesh.getElements(dim, tag=-1), and then looking for markers?

You can write your own wrapper for this if you would like to. As it is using the GMSH Python API you can create your own custom reader.

Adding all cells regardless of tags is not necessarily a good idea, as GMSH usually only generates cells/facets/edges for Physical entities, to avoid duplicated cells (as in your case).

From section 1.2.3 of https://gmsh.info/doc/texinfo/gmsh.html

1.2.3 Elementary entities vs. physical groups
It is usually convenient to combine elementary geometrical entities into more meaningful groups, e.g. to define some mathematical (“domain”, “boundary with Neumann condition”), functional (“left wing”, “fuselage”) or material (“steel”, “carbon”) properties. Such grouping is done in Gmsh’s geometry module (see Geometry module) through the definition of “physical groups”.
By default in the native Gmsh MSH mesh file format (see Gmsh file formats), as well as in most other mesh formats, if physical groups are defined, the output mesh only contains those elements that belong to at least one physical group. (Different mesh file formats treat physical groups in slightly different ways, depending on their capability to define groups.) To save all mesh elements whether or not physical groups are defined, use the Mesh.SaveAll option (see Mesh options) or specify -save_all on the command line. In some formats (e.g. MSH2), setting Mesh.SaveAll will however discard all physical group definitions.

As I previously stated, our MeshTags objects will not be able to read in duplicate elements in a sensible way, so rather having a slightly more verbose GMSH tagging routine, as described in

 phs1_in_srce = [e for e in phs1 if e in srce]
        phs1_out_srce = [e for e in phs1 if not(e in srce)]
        phs2_in_srce = [e for e in phs2 if e in srce]
        phs2_out_srce = [e for e in phs2 if not(e in srce)]

        model.add_physical_group(2, phs1_in_srce, 10)
        model.add_physical_group(2, phs1_out_srce, 11)
        model.add_physical_group(2, phs2_in_srce, 20)
        model.add_physical_group(2, phs2_out_srce, 21)
        model.setPhysicalName(2, 10, "Phase 1 - inside source")
        model.setPhysicalName(2, 11, "Phase 1 - outside source")
        model.setPhysicalName(2, 20, "Phase 2 - inside source")
        model.setPhysicalName(2, 21, "Phase 2 - outside source")

would be preferable from my point of view.

[pierricmora]

Ok, yes it makes sense from Gmsh doc. In my case

e = model.mesh.getElements(dim=2, tag=-1)
print(e[1][0].shape)  # (154,): the correct number of cells

does produce the same output regardless the way I have defined - or not - physical groups, so I thought it could be the stable way to go. But indeed, given what's stated in the doc, maybe not, so I think I'll follow your advice. This is a pity, because I find it so much more convenient to work with overlapping or incomplete tags.

[pierricmora]

@jorgensd Then, what about adding these simple checks and error messages to avoid crashes and give the user a hint on what is expected?

The following seems to work on my side. In model_to_mesh, between lines 258 and 260, add:

        # Sort elements by ascending dimension
        perm_sort = np.argsort(cell_dimensions)

        # Check that all cells are tagged once
        _d = model.getDimension()
        assert _d in topologies.keys(), 'All cells are expected to be tagged once; none found'
        _elementTypes, _elementTags, _nodeTags = model.mesh.getElements(dim=_d, tag=-1)
        # assert len(_elementTypes) == 1  # assert only one type of elements; already checked in extract_topology_and_markers
        nbcells = len(_elementTags[0])
        nbcells_tagged = len(topologies[_d]["topology"])
        nbcells_tagged_once = len(np.unique(topologies[_d]['topology'], axis=0))
        assert nbcells == nbcells_tagged, f'All cells are expected to be tagged once; found: {nbcells_tagged}, expected: {nbcells}'
        assert nbcells_tagged == nbcells_tagged_once, 'All cells are expected to be tagged once; found duplicates'

        # Broadcast cell type data and geometric dimension

[jorgensd]

The only thing that worries me here is that the call to nbcells_tagged_once can be quite expensive, as illustrated in the following mwe:

import numpy as np
from time import perf_counter

N = 1_000_000
M = 8

a = np.arange(N*M).reshape(N, M).astype(np.int64)

b = np.vstack([a, a, a]).astype(np.int64)

start = perf_counter()
unique_b = np.unique(b, axis=0)
end = perf_counter()
print(end-start)
print(f"{unique_b.shape=}, {a.shape=}, {b.shape=}")

yielding

1.6669321240005956
unique_b.shape=(1000000, 8), a.shape=(1000000, 8), b.shape=(3000000, 8)

It would be interesting to see how this modification changes the runtime of "standard" problems. Feel free to make a pull request, and test some standard meshes.
I guess worst case scenario, one has to run Python with -O to turn of the assert.

[pierricmora]

By playing with N on your snippet I get a linear behaviour for the time spent in np.unique.

Then, I extended my example by adding a mesh_size_max input parameter for model generation, and a mute_assert bool in my modified model_to_mesh to turn on/off the assertions, including np.unique. I ran a few tests representing from 41k to 12M cells (mesh_size_max = 0.01 to 0.01 / 16) and I obtained a relatively constant +15% increase, unless for the smallest mesh which resulted in a +30% increase (the trend is a slow decrease: largest mesh got a +10%). Is it acceptable? If yes I'll make the PR.

Below is the code I run.

mesh_size_max = 0.01

# Create model
start = perf_counter()
model = gmsh_rect_ellipse(mode=3, mesh_size_max=mesh_size_max, show=False)
end = perf_counter()
gen_model_time = end - start

# model_to_mesh, original
start = perf_counter()
mesh, ct, ft = model_to_mesh(model, MPI.COMM_WORLD, rank=0, gdim=2, mute_assert=True)
end = perf_counter()
model_to_mesh_woutAssert_time = end - start

# model_to_mesh with checks
start = perf_counter()
mesh, ct, ft = model_to_mesh(model, MPI.COMM_WORLD, rank=0, gdim=2, mute_assert=False)
end = perf_counter()
model_to_mesh_withAssert_time = end - start

# say hello
print(f'mesh_size_max: {mesh_size_max},\n' + \
        f'number of cells: {len(model.mesh.getElements(dim=2, tag=-1)[1][0])},\n' + \
        f'model generation time (Gmsh): {gen_model_time:.3f}s,\n' + \
        f'model_to_mesh time, original version: {model_to_mesh_woutAssert_time:.3f}s,\n' + \
        f'model_to_mesh time, with checks: {model_to_mesh_withAssert_time:.3f}s,\n' + \
        f'overcost due to checks: {(model_to_mesh_withAssert_time - model_to_mesh_woutAssert_time) / model_to_mesh_woutAssert_time * 100:.1f}%')

Outputs:
mesh_size_max: 0.01,
number of cells: 47142,
model generation time (Gmsh): 0.843s,
model_to_mesh time, original version: 0.175s,
model_to_mesh time, with checks: 0.223s,
overcost due to checks: 26.9%

mesh_size_max: 0.005,
number of cells: 186674,
model generation time (Gmsh): 3.652s,
model_to_mesh time, original version: 0.913s,
model_to_mesh time, with checks: 1.035s,
overcost due to checks: 13.4%

mesh_size_max: 0.0025,
number of cells: 743524,
model generation time (Gmsh): 15.690s,
model_to_mesh time, original version: 4.133s,
model_to_mesh time, with checks: 4.747s,
overcost due to checks: 14.9%

mesh_size_max: 0.00125,
number of cells: 2963030,
model generation time (Gmsh): 69.363s,
model_to_mesh time, original version: 21.390s,
model_to_mesh time, with checks: 24.306s,
overcost due to checks: 13.6%

mesh_size_max: 0.000625,
number of cells: 11841660,
model generation time (Gmsh): 329.698s,
model_to_mesh time, original version: 122.245s,
model_to_mesh time, with checks: 134.240s,
overcost due to checks: 9.8%

[jorgensd]

I would be fine with this. 10% overhead isn’t to bad, as one should preferrably only use this wrapper once and not on hpc systems, as Gmsh only writes serial files and Gmsh.model doesn’t do mpi.