fix physical tags when generating the mesh

[capitalaslash]

Closes #397.

[capitalaslash]

I not super-happy with the way I fixed this issue, but it works for my applications.
I based it on the way meshio sets up the physical tags when reading a msh file.

[capitalaslash]

Do not merge. The geometrical tag must also be set to get this properly working, I am adding it.

[capitalaslash]

The branch is ready for review.
The geometrical tag is set reusing the same code that is used for the physical tag, so it is set only where there are physical tags set.
When no physical tags are set, those 2 tags will not be present (meshio will print 0s in msh format, as they are always present in the file format).
I believe the default behavior in gmsh when no physical tags are present is to set them to the geometrical one? Not sure it would be useful.

[nschloe]

I not super-happy with the way I fixed this issue,

Neither. I don't understand why this PR switches away from cell_sets to cell_data + field_data. I've now fixed the corresponding bug simply by deprecating empty labels. Is there anything else in the PR that's worth merging?

[capitalaslash]

I not super-happy with the way I fixed this issue,

Neither. I don't understand why this PR switches away from cell_sets to cell_data + field_data. I've now fixed the corresponding bug simply by deprecating empty labels. Is there anything else in the PR that's worth merging?

The problem goes beyond the referenced bug. A mesh constructed in pygmsh cannot be printed in gmsh41 format, spitting out errors about

  File "/home/aslash/.local/lib/python3.9/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

nor in gmsh22 format, since that one would not use cell_sets.
Other formats (such as vtu) ignore cell_sets.

My problem is simple: I want to print out a mesh that I can open in gmsh directly.

So I see 2 ways out: meshio._gmsh22 adds support for physical tags stored as cell_sets or pygmsh creates a mesh that can be fed to meshio._gmsh41.
Let me know what a suitable strategy would be to get over this problem.

[nschloe]

My problem is simple: I want to print out a mesh that I can open in gmsh directly.

Yeah, that should be made possible. Note that you can always fall back to gmsh like

geom.generate_mesh()
gmsh.write("mesh3D.msh")

I'm wondering how gmsh itself writes the file, and if the issue should not be fixed in meshio after all.

[capitalaslash]

geom.generate_mesh()
gmsh.write("mesh3D.msh")

That works. Why not add that to the pygmsh interface?

I'm wondering how gmsh itself writes the file, and if the issue should not be fixed in meshio after all.

Physical and geometrical (in pygmsh, elementary in gmsh) tags in gmsh are special: they are always present, that is why this branch moves that information to cell_data.
For any gmsh compatible mesh those flags must be there. If the user did not specify the physical tags, they correspond to the elementary ones.
This information cannot be generated at the meshio stage since there is no gmsh.model available there, so they must be stored somewhere. cell_sets can be adequate for physical tags, but they are awkward for elementary tags: there is always a geometrical object behind a mesh element so they cannot be defined on a partition of the mesh as we can do for the physical tags.

I'll point out that even meshio._gmsh41 when writing expects that that information is stored in cell_data

https://github.com/nschloe/meshio/blob/b5229a8181745bd63d819810a7f0dbcb4ace224a/meshio/gmsh/_gmsh41.py#L338

[nschloe]

Why not add that to the pygmsh interface?

Yeah, why not actually. Just added it as pygmsh.write() (#426).

[nschloe]

Physical and geometrical (in pygmsh, elementary in gmsh) tags in gmsh are special: they are always present, that is why this branch moves that information to cell_data.

A yes, right. Every cell is guaranteed to have both of those tags so we can treat it as a function (aka cell_data), not a cell_set.

I can't quite read the PR, still. Perhaps it could be made clearer, I don't know. For example, instead of

 cell_data["gmsh:physical"] = [
    numpy.array([0 for _ in range(len(cells[i][1]))]) for i in range(len(cells))
)

one could go

 cell_data["gmsh:physical"] = [
    numpy.array([0 for _ in range(len(cell[1]))]) for cell in cells
)

but I still don't understand what signficance cell[1] has, or why it's initialized with 0. Isn't it overridden afterwards? Perhaps not, who knows what values elem_tags - 1 can take.

Considering the gmsh error your describe, I think we need a "gmsh:dim_tags" point data entry.

[capitalaslash]

I can't quite read the PR, still. Perhaps it could be made clearer, I don't know. For example, instead of

 cell_data["gmsh:physical"] = [
    numpy.array([0 for _ in range(len(cells[i][1]))]) for i in range(len(cells))
)

one could go

 cell_data["gmsh:physical"] = [
    numpy.array([0 for _ in range(len(cell[1]))]) for cell in cells
)

I improved there with numpy.full.

but I still don't understand what signficance cell[1] has, or why it's initialized with 0. Isn't it overridden afterwards? Perhaps not, who knows what values elem_tags - 1 can take.

cell[1] stores the size of the group of elements. This guarantees that the array are properly sized, since cell_data must be complete as discussed above. Zero is an arbitrary initialization value (that gmsh does not use so it's easy when debugging).

After some thoughts, I believe the best way would be to always fill in a gmsh:geometrical cell_data based on the elementary tags. The gmsh:physical is then added only if the model has some physical tag, otherwise is the same as the first one, similarly to how gmsh fills them.

I will pursue this if you deem it mergeable.

Considering the gmsh error your describe, I think we need a "gmsh:dim_tags" point data entry.

I will add a separate bug report (and hopefully a merge request) about this.

[nschloe]

I improved there with numpy.full.

It should be zeros or empty I think.

I will pursue this if you deem it mergeable.

Let's first get the "gmsh:dim_tags" issue out of the way.

[nschloe]

This might have been fixed by #486.

[capitalaslash]

This might have been fixed by #486.

unfortunately not fixed yet.

[nschloe]

Okay, perhaps this becomes clearer with a test. Could you add one?

[capitalaslash]

Okay, perhaps this becomes clearer with a test. Could you add one?

see #488

[nschloe]

Thanks for the test! We need one that fails before this PR, though, and passes with the PR.

[capitalaslash]

Thanks for the test! We need one that fails before this PR, though, and passes with the PR.

test_physical.py should have failed with the other merge, but it does not in pytest, I'm trying to understand why...

[capitalaslash]

my bad, I asserted the read mesh in the main function instead of the test one.
see #489