Improvements to skfem.io.meshio #680

kinnala · 2021-08-03T15:51:20Z

Write/read boundaries/subdomains to/from external file formats
Optionally, read mesh data from external formats

Fixes #158
Fixes #261
Fixes #271
Fixes #683

Example 1: Export boundaries to external formats.

m = MeshHex().with_boundaries({'left': lambda x: x[0] == 0, 'bottom': lambda x: x[2] == 0})
m.save('test.msh', file_format='gmsh22')

Example 2: Load the exported boundaries

In [1]: from skfem import *                          
In [2]: m = MeshHex().with_boundaries({'left': lambda x: x[0] == 0, 'bottom': lambda x: x[2] == 0})       
In [3]: m.save('test.vtk')                           
In [4]: M = MeshHex.load('test.vtk')                 
In [5]: M.boundaries                                 
Out[5]: {'left': array([0]), 'bottom': array([2])}
In [6]: m.boundaries                                 
Out[6]: {'left': array([0]), 'bottom': array([2])}

Example 3: Load raw data from meshio

In [1]: from skfem import *                          
In [2]: out = ['cell_data']                          
In [3]: m = Mesh.load('docs/examples/meshes/beams.msh', out=out)                                         
In [4]: out[0].keys()                               
Out[4]: dict_keys(['gmsh:physical', 'gmsh:geometrical'])

gdmcbain · 2021-08-03T23:18:53Z

Good morning. Serialization #158 is back on the agenda, excellent: I believe this is an important long-term goal. Would this fix
#261 ‘Mesh.save boundaries and subdomains’ ?

There's quite a history of efforts on that here and even more in meshio. I'll have to refresh myself on that to see what the status is. Prior to that, my impressions are that:

meshio.Mesh.cell_sets and equivalents (!) have not been implemented consistently with respect to different external file formats
it is not possible to write a nontrivial mesh, i.e. with subdomains or boundaries, in MSH 4.1 using meshio

I do however routinely save meshes with subdomains and boundaries. I can outline my current practice, in case that's helpful, or share code—basically I write MSH 2.2 using meshio using cell_data rather than cell_sets. Are you trying to solve a specific problem? Or rather to generally extend the capabilities of skfem.io? I can also share my current thoughts on strategy for #261, if the current PR doesn't cover it.

While there may not be documentation for meshio.Mesh.cell_sets, ‘cell_sets’ is mentioned in 29 issues ! A lot of discussion !

gdmcbain · 2021-08-04T03:36:16Z

it is not possible to write a nontrivial mesh, i.e. with subdomains or boundaries, in MSH 4.1 using meshio

What I mean by this is that if just before

scikit-fem/docs/examples/ex28.py

Line 145 in 80e39e2

return from_meshio(geom.generate_mesh(dim=2))

I put

        mesh = geom.generate_mesh(dim=2)
        mesh.write(Path(__file__).with_suffix(".msh"), binary=False)

and then in the script

make_mesh(length, halfheight, thickness)

I get

Traceback (most recent call last):
  File "/home/gdmcbain/src/scikit-fem/docs/examples/ex28_make_mesh.py", line 149, in <module>
    make_mesh(length, halfheight, thickness)
  File "/home/gdmcbain/src/scikit-fem/docs/examples/ex28_make_mesh.py", line 146, in make_mesh
    mesh.write(Path(__file__).with_suffix(".msh"), binary=False)
  File "/home/gdmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/meshio/_mesh.py", line 219, in write
    write(path_or_buf, self, file_format, **kwargs)
  File "/home/gdmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/meshio/_helpers.py", line 141, in write
    return writer(filename, mesh, **kwargs)
  File "/home/gdmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/meshio/gmsh/main.py", line 111, in <lambda>
    "gmsh": lambda f, m, **kwargs: write(f, m, "4.1", **kwargs),
  File "/home/gdmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/meshio/gmsh/main.py", line 102, in write
    writer.write(filename, mesh, binary=binary, float_fmt=float_fmt)
  File "/home/gdmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/meshio/gmsh/_gmsh41.py", line 358, in write
    _write_nodes(fh, mesh.points, mesh.cells, mesh.point_data, float_fmt, binary)
  File "/home/gdmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/meshio/gmsh/_gmsh41.py", line 607, in _write_nodes
    raise WriteError(
meshio._exceptions.WriteError: Specify entity information (gmsh:dim_tags in point_data) to deal with more than one cell type.

I'm pretty sure that this is a meshio issue.

However, if I insert

file_format="gmsh22",

into the call to meshio.write, all goes well : ex28_make_mesh.msh.txt

kinnala · 2021-08-04T09:04:29Z

So I'm planning to write meshio.Mesh.cell_sets inside skfem.io.meshio.to_meshio in order to call meshio.Mesh.sets_to_int_data and then use the resulting meshio.Mesh.cell_data for further purposes. But I think I'm not properly understanding the correct format of cell_sets.

kinnala · 2021-08-04T09:12:36Z

I might have now got it correctly. Pushing some more commits soon.

gdmcbain · 2021-08-04T10:38:13Z

O. K.

The best way to see what meshio.Mesh.cell_sets should look like is to meshio.read a genuine MSH 4.1 file (genuine as in generated by Gmsh itself, e.g. from GEO code); meshio reads MSH 4.1 O. K., it's the writing that's tricky. I'm not sure it's even conceptually possible to write 4.1 for an arbitrary mesh that is acceptable to Gmsh.

kinnala · 2021-08-04T11:03:57Z

tests/test_mesh.py

+    from tempfile import NamedTemporaryFile
+    with NamedTemporaryFile() as f:
+        m.save(f.name + ".inp", write_boundaries=True)
+        m2 = Mesh.load(f.name + ".inp")


This is funny that Abaqus inp format works properly for write/read cycle but almost nothing else does without sorcery.

kinnala · 2021-08-04T11:07:04Z

This allows now doing stuff like this:

In [1]: from skfem import *                          

In [2]: m = MeshTri().with_boundaries({'left': lambda x: x[0] == 0})                                      

In [3]: m.save('testing.msh', write_boundaries=True, sets_to_int_data=True, file_format='gmsh22')         
WARNING:root:msh2 requires 3D points, but 2D points given. Appending 0 third component.
WARNING:root:Binary Gmsh needs 32-bit integers (got int64). Converting.
WARNING:root:Binary Gmsh needs 32-bit integers (got int64). Converting.
WARNING:root:Appending zeros to replace the missing physical tag data.
WARNING:root:Appending zeros to replace the missing geometrical tag data.

In [4]: M = MeshTri.load('testing.msh', int_data_to_sets=True)                                            

In [5]: M.boundaries                                 
Out[5]: 
{'set-1.0': array([0, 3, 4]),
 'set0.0': array([1]),
 'gmsh:physical': array([0, 1, 3, 4]),
 'gmsh:geometrical': array([0, 1, 3, 4])}

In [6]: m.boundaries                                 
Out[6]: {'left': array([1])}

As you can see, the tag names are lost due to the sets_to_int_data/int_data_to_sets cycle, but the correct set is there in the end.

gdmcbain · 2021-08-04T20:15:53Z

Fantastic!

As you can see, the tag names are lost due to the sets_to_int_data/int_data_to_sets cycle, but the correct set is there in the end.

I wonder whether the names can be stashed in field_data or similar. (Or in the names of the cell data?)

That's weird and unfortunate about Abaqus .inp. Is that a viable format? Can one do anything with it without access to the Abaqus program? It might just be that better work has been done on that submodule of meshio.

MSH 2.2 is not ideal, it's very slow. In my first industrial application of scikit-fem, reading MSH 2.2 was actually the bottleneck.

gdmcbain

Oh, no, I see that we've gone from having one named boundary

{'left': array([1])}

to four, all misnamed, and three spurious

{'set-1.0': array([0, 3, 4]), 'set0.0': array([1]), 'gmsh:physical': array([0, 1, 3, 4]), 'gmsh:geometrical': array([0, 1, 3, 4])}

Yeah, no, that's not so good. Hmm…Some further thoughts in #261.

skfem/io/meshio.py

tests/test_mesh.py

skfem/io/meshio.py

kinnala · 2021-08-07T11:44:20Z

MSH 4.1 will likely require some workaround, e.g., writing gmsh:dim_tags. I need to look into that further.

kinnala · 2021-08-07T17:20:21Z

If I remove WriteError line from meshio, then everything works just fine. However, I cannot figure out what is this gmsh:dim_tags.

gdmcbain · 2021-08-07T19:53:56Z

MSH 4.1 will likely require some workaround,

Personally i wouldn't see writing MSH 4.1 as a requirement. Is that needed for interoperation with Elmer?

It might (for my purposes #261) be worth working out a format that meshio writes and reads reasonably quickly. It does read MSH 4.1 much faster than MSH 2.2; this is an intrinsic benefit of the redesign from the Gmsh tesm, avoiding mixed topologies.

Unfortunately, the meshio speed tests don't test mixed topologies, and so while XDMF looks fast, that's only true if only one kind of cell is written. The current meshio XDMF writer, if given domain and boundary cells, will revert to a mixed topology which i fear will be slow to reread—each cell will have to be processed individually to see whst kind it is. This is not intrinsic to XDMF but a shortcoming of meshio; a fix for that is proposed in nschloe/meshio#622.

This was part of the motivation for writing boundary facets as such in #681 rather than as lower-dimensional cells (besides the main motivation of that being a more direct representation). But neither of those are overriding considerations ; good interoperation is.

I'm thinking that VTK will probably be a reasonable default for using this PR to fix #261 ; I.e. write and reread by skfem.

If Elmer doesn't use meshio to read MSH 2.2 but rather has a C++ reader, it might not be slow, so that could be a workable procedure there.

kinnala · 2021-08-09T09:13:05Z

Now user can pass the optional parameter out to Mesh.load. If it is a dictionary, it gets populated with cell_data, point_data and field_data.

kinnala · 2021-08-09T10:37:50Z

I noticed that in ex28.py, the lines depending on pygmsh are not being run at all. So I removed the unused lines and, thus, reverted the license to BSD-3-Clause. I suppose that's fine by you @gdmcbain?

kinnala · 2021-08-09T10:41:01Z

Note that I'm still planning to do #474 at some point, which will hopefully include stuff like how to use pygmsh.

gdmcbain · 2021-08-09T11:41:56Z

Yes, I noticed that too in ex28 and was puzzled. Thanks, it's better now.

gdmcbain · 2021-08-10T04:22:46Z

I have just fetched this branch and taken it for a spin. It's awesome. I wish we had done this years ago! So far everything just works, very easily, no surprises; very nice to use.

One slight disappointment is that the new docs/examples/ex28.msh is 156K whereas the old JSON was only 52K. A bit of a shame since the latter is not compressed at all and was hardly designed for efficiency. VTK is a bit better at 128K and parity is restored with XDMF at 48+4=52K, with the additional requirement of h5py (not difficult but it is an additional requirement). I suppose all our simple boolean arrays of indicator functions for subdomains and uints for boundaries are getting stored as big floats (float64?). I don't think this matters too much; in particular I don't think it is too much to pay for the ease of being able to immediately open the VTK or XDMF in ParaView and immediately see each subdomain and boundary, along with e.g. point_data={"temperature": temperature[basis["heat"].nodal_dofs.flatten()]}. Conclusion: the indicator functions are verbose but portable, which is a fair compromise and probably the right one to make here.

The pinked visualization of the boundaries is a bit zany, but it kind of makes sense, I think; I wonder what others will think of it.

Is this ready for a review? Anyway, I'm going to begin using it for everything immediately.

kinnala · 2021-08-10T07:25:01Z

Yes, I think it's ready.

kinnala · 2021-08-10T07:29:41Z

I updated the first message in the PR with some examples.

gdmcbain · 2021-08-11T04:13:41Z

tests/test_mesh.py

 class Loading(TestCase):
    """Check that Mesh.load works properly."""

    def runTest(self):
        # submeshes
-        path = Path(__file__).parents[1] / 'docs' / 'examples' / 'meshes'
-        m = MeshTet.load(str(path / 'box.msh'))
+        m = MeshTet.load(str(MESH_PATH / 'box.msh'))


I don't think we need to cast to str; see, e.g.,

scikit-fem/docs/examples/ex21.py

Line 70 in 96fa97a

m = MeshTet.load(Path(__file__).parent / 'meshes' / 'beams.msh')

scikit-fem/docs/examples/ex28.py

Line 147 in 96fa97a

mesh = from_file(Path(__file__).parent / 'meshes' / 'ex28.json')

just load Path

Add support for exporting cell_sets

1a57b1f

Add flags for writing facets and converting to/from integer data

c8d0610

kinnala commented Aug 4, 2021

View reviewed changes

kinnala added 2 commits August 4, 2021 14:42

Fix flake8

f68c697

Workaround for a bug in meshio

3c46d6e

gdmcbain mentioned this pull request Aug 5, 2021

Mesh.save subdomains and boundaries #261

Closed

gdmcbain reviewed Aug 5, 2021

View reviewed changes

skfem/io/meshio.py Outdated Show resolved Hide resolved

gdmcbain mentioned this pull request Aug 5, 2021

WIP: serialize subdomains & boundaries #681

Closed

Simplify and test using gmsh22

d09de69

kinnala commented Aug 6, 2021

View reviewed changes

tests/test_mesh.py Outdated Show resolved Hide resolved

Encode subdomains/boundaries as cell_data during export

ba563dd

gdmcbain reviewed Aug 6, 2021

View reviewed changes

skfem/io/meshio.py Outdated Show resolved Hide resolved

kinnala changed the title ~~Add support for exporting cell_sets~~ Export subdomains/boundaries Aug 6, 2021

kinnala added 7 commits August 6, 2021 14:16

Write tags using point_data for Paraview visualization

c62c81c

Fix flake8

507b310

Add a new test dependency

84e136b

Simplifications

679bfde

Do not encode subdomains as vertex data; it's already in cell data

44a16a4

Test point_data export

d0a8de6

Add a changelog entry

8f107f6

kinnala added 2 commits August 9, 2021 11:45

Fix the case with no boundaries/subdomains and add tests

9878700

Optionally export mesh data

d7d5d40

kinnala changed the title ~~Export subdomains/boundaries~~ Improvements to skfem.io.meshio Aug 9, 2021

Add a changelog entry

f7220bc

kinnala added 2 commits August 9, 2021 13:43

Remove unused lines from ex28 and modify the license accordingly

5540c68

Move subdomain/boundary encoding to Mesh class

07d9486

kinnala added 2 commits August 9, 2021 15:53

Fix test

503eeb5

Fix changelog entry

90fe988

kinnala added 2 commits August 10, 2021 09:56

Update docstrings

6c49c78

Change dict to list

2297b2f

Fix bug

509db9a

just load Path

79f5729

gdmcbain mentioned this pull request Aug 11, 2021

just load Path #685

Merged

gdmcbain approved these changes Aug 11, 2021

View reviewed changes

Merge pull request #685 from gdmcbain/load-path-no-need-to-cast-to-str

0fef727

just load Path

kinnala merged commit e7196f6 into master Aug 11, 2021

gdmcbain added a commit to gdmcbain/OpenFOAM-9 that referenced this pull request Aug 13, 2021

rename SkFem from ToMeshio following kinnala/scikit-fem#680

d4afc0e

This was referenced Aug 13, 2021

Examples/snippets demonstrating the integrations #474

Closed

from_meshio fails on old ex27 #691

Closed

This was referenced Aug 23, 2021

take care reading MSH 4.1 with meshio #693

Closed

preferred format for static meshes #696

Closed

gdmcbain mentioned this pull request Aug 30, 2021

periodic boundary conditions #688

Closed

kinnala deleted the improvements-to-meshio-support branch November 4, 2021 07:09

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Improvements to skfem.io.meshio #680

Improvements to skfem.io.meshio #680

kinnala commented Aug 3, 2021 •

edited

gdmcbain commented Aug 3, 2021

gdmcbain commented Aug 4, 2021

kinnala commented Aug 4, 2021

kinnala commented Aug 4, 2021

gdmcbain commented Aug 4, 2021

kinnala Aug 4, 2021

kinnala commented Aug 4, 2021

gdmcbain commented Aug 4, 2021

gdmcbain left a comment

kinnala commented Aug 7, 2021

kinnala commented Aug 7, 2021

gdmcbain commented Aug 7, 2021

kinnala commented Aug 9, 2021

kinnala commented Aug 9, 2021

kinnala commented Aug 9, 2021

gdmcbain commented Aug 9, 2021

gdmcbain commented Aug 10, 2021

kinnala commented Aug 10, 2021

kinnala commented Aug 10, 2021

gdmcbain Aug 11, 2021

Improvements to skfem.io.meshio #680

Improvements to skfem.io.meshio #680

Conversation

kinnala commented Aug 3, 2021 • edited

gdmcbain commented Aug 3, 2021

gdmcbain commented Aug 4, 2021

kinnala commented Aug 4, 2021

kinnala commented Aug 4, 2021

gdmcbain commented Aug 4, 2021

kinnala Aug 4, 2021

Choose a reason for hiding this comment

kinnala commented Aug 4, 2021

gdmcbain commented Aug 4, 2021

gdmcbain left a comment

Choose a reason for hiding this comment

kinnala commented Aug 7, 2021

kinnala commented Aug 7, 2021

gdmcbain commented Aug 7, 2021

kinnala commented Aug 9, 2021

kinnala commented Aug 9, 2021

kinnala commented Aug 9, 2021

gdmcbain commented Aug 9, 2021

gdmcbain commented Aug 10, 2021

kinnala commented Aug 10, 2021

kinnala commented Aug 10, 2021

gdmcbain Aug 11, 2021

Choose a reason for hiding this comment

kinnala commented Aug 3, 2021 •

edited