WIP: serialize subdomains & boundaries

[gdmcbain]

Fixes #261.

This is to explore some alternative ideas inspired by consideration of #680 that don't involve dealing with meshio.Mesh.cell_sets.

[gdmcbain]

This is with 59a474d.

[gdmcbain]

The idea for boundaries needs rethinking since a cell can have more than one facet on a boundary.

[gdmcbain]

Pressing ahead anyway for the moment, assuming that that doesn't happen (and checking for that)

https://github.com/gdmcbain/scikit-fem/blob/5bac8a3e469a89124b939ce9eb97ac91eca42fed/docs/examples/ex_261.py#L22-L23

[gdmcbain]

What about binary-encoding the facets of each cell that are on the boundary? 9e8cbe8

             (2 ** np.arange(mesh.t.shape[0])) @ np.isin(mesh.t2f, boundary)

[gdmcbain]

It seems to work on a MeshTet too. 6ec2cac

[kinnala]

I don't want us to invent a new representation because it won't be compatible with existing tools. E.g., ElmerFEM which I'm using in an upcoming project uses cell_data + explicit boundary mesh. I think code_aster also does it that way. Rather I'd fix the meshio gmsh exporter if .msh must be used, but I'm not very convinced of that format either, so I'll keep looking.

[gdmcbain]

A nice thing about this approach is that one can change

    meshio.write(Path(__file__).with_name(f"{name}.vtk"), mio)

from "vtk" to "xdmf" &c. and get XDMF which works the same but is much smaller (24 + 4K for beams.h5 & .xdmf, v. 60K for beams.vtk and 36K for the original meshes/beams.msh).

Doesn't work for "e" though; looks like meshio.write doesn't write cell_data for Exodus. :(

[gdmcbain]

Oops, didn't see your comment. Sorry. Yeah, O. K., I'll keep this approach for in-house.

[kinnala]

So basically the issue is that cell_data gets written for many formats but almost none of the formats (in meshio) support an explicit boundary mesh so you cannot use cell_data to define boundary conditions? Is this actually an issue of the formats or is it a lacking support in meshio?

How about point_data? Can tools figure out boundary facets based on those?

[gdmcbain]

So basically the issue is that cell_data gets written for many formats but almost none of the formats (in meshio) support an explicit boundary mesh so you cannot use cell_data to define boundary conditions?

I think the fairest way to answer this would be to say that meshio tries not impose any interpretation on .cell_data. Each entry has to have the same shape as .cells, i.e. be a list of one-dimensional np.ndarray with size equal to the number of cells in each block. I think the dtype has to be float too.

But no, meshio has no concept of subdomains or sub-boundaries. If any mesh file format does, meshio more or less ignores it, maybe treating it as .cell_data or maybe .cell_sets, or maybe point-data or sets.

There is "gmsh:physical" and "medit:ref" and maybe a couple of others. meshio-convert does turn "gmsh:physical" into "medit:ref" and a couple of others do something similar but it's kind of a special case, probably because "gmsh:physical" has been in meshio for a long time and Gmsh is a very popular mesh-generator. And it is true that a lot of external tools do understand Gmsh's MSH 2.2 and do interpret its physical entities and lower-dimensional cells as subdomains and boundaries. I think PETSc's DMPlex does, and a lot of other tools use that (libMesh, Firedrake, ...). It was because of that that I pushed to get something like "gmsh:physical" into meshio in the first place, but I have since been persuaded by

(I've been pondering a comment in nschloe/meshio#571: ‘The devil invented mixed topologies.’) (#261)

and the data-structures skfem.Mesh.subdomains and skfem.Mesh.boundaries that that's not so natural.

[kinnala]

What's going on in MSH 4.1? Is there some kind of implicit boundary representation?

[kinnala]

Let's move the discussion to #261.

[gdmcbain]

Reopening to include new commits demonstrating that this approach works automatically for "interfaces" #683 without needing to distinguish them from boundaries.

[gdmcbain]

Superseded by #680.