[BUG]: gmshio causes kernel crash if entities are tagged several times about dolfinx HOT 8 CLOSED

@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 (

from dolfinx.

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?

from dolfinx.

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.

from dolfinx.

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.

from dolfinx.

@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

from dolfinx.

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.

from dolfinx.

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%

from dolfinx.

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.

from dolfinx.