preferred format for static meshes

[gdmcbain]

Should docs/examples prefer a particular format for static meshes? Or rather demonstrate the ability to read from multiple formats?

The first split is between the internal JSON format and whatever can be read from meshio. The second split, in the latter, is between Gmsh's MSH, which can contain subdomains and boundaries (as (n − 1)-dimensional cells), other formats in which subdomains and boundary-facets have been encoded using #680, and other formats which are just plain meshes without subdomains or boundaries (e.g. which might have those defined subsequently using .with_boundaries, &c.). Thus:

• skfem.io.json.
• skfem.io.meshio
  • Gmsh's MSH
    • MSH 4.1 (see also ignore cell_sets "gmsh:*" inserted by meshio #695)
    • MSH 2.2
  • meshio-read formats with subdomains or boundary-facets encoded using Improvements to skfem.io.meshio #680
  • meshio-read formats of plain meshes without subdomains or boundaries

Does that cover it?

[gdmcbain]

Should the old internal JSON format be deprecated? It does have the advantage that it depends only on the Python Standard Library, not using meshio (which does evolve, sometimes by leaps & bounds). It's portable in that it is JSON but I don't think anything other than scikit-fem will understand the encoding within the JSON. There's also the slight disadvantage that the serialization of subdomains and boundaries in the JSON differs from that adopted in #680: it lists indices of subsets of cells and facets rather than enumerating the corresponding indicator-function for the subdomains or the indicator function multiplied by binary-encoding of subsets of facets for boundaries. It's probably not worth revising the JSON to match the new serialization, but what's the role of the JSON henceforth? I don't have an opinion on this, I probably won't use it much given the new options available in #680 (unless maybe a makeshift is required during a future upheaval of meshio).

I don't think there's anything that needs doing here, unless it is to replace the JSON in docs/examples/meshes should the decision be to deprecate it.

[gdmcbain]

It might be worthwhile indicating to the user what a good fast compact format for serializing meshes with subdomains and boundaries (and possibly other point- or cell-data).

I did have some findings on this in #680 (2021-08-10) and #693.

In the former:

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 her

In the latter, some steps were taken to replace the JSON but withdrawn as that made the pull-request as too composite. It was replacing

scikit-fem/docs/examples/ex27.py

Line 109 in e4a4b2b

self.mesh = skfem.io.json.from_file(Path(__file__).parent / 'meshes' / 'backward-facing_step.json')

with first MSH 4.1 dac551a and then VTU f799d23 for exx 24, 27. In that case, the MSH 4.1 was slightly smaller (178 K) than the JSON (193 K) and the VTU was much smaller (82.4 KB).

VTU was chosen after inspecting the charts in the meshio Performance comparision, passing over those formats don't support cell-data (Exodus) or have additional dependencies (Exodus, MED, XDMF), and formats that I just wasn't familiar with (MOAB).

I've been using XDMF heavily the last couple of years and it suits me well, but there might be situations where one mightn't want to depend on HDF5.

MSH 2.2 is a poor choice for large files because it's very slow to read. It's one of several formats that defines the type of each cell so that the cells have to be read one at a time. This is intrinsically slow in Python. The fast formats have blocks of cells of the same type, which can be usually be read quickly with NumPy. Thus when Gmsh updated to MSH 4.1, it became literally a thousand times faster to read in my real-world application.

Note that XDMF (which has quite a broad specification, only incompletely covered by meshio) writes inhomogeneous meshes
(e.g. triangles plus lines to represent boundary-facets) either as homogeneous blocks or as a ‘mixed topology’ but currently only the latter is supported by meshio (possibly pending nschloe/meshio#622) and I suspect that that would radically lower the position of XDMF in the Performance Comparison (which only compares homogeneous meshes, which are stored in XDMF by meshio as a homogeneous block). Since #680 doesn't use (n − 1)-dimenisonal cells to represent boundary-facets, it doesn't encounter the mixed-topology problem, originally raised in xdmf#15 multiple topology blocks; this is a distinct advantage of #680.

[kinnala]

What I think would be important is to have some serialization format which is directly compatible with other tools in the sense that if we export "boundaries" then some other FEM package can use them via their own importers. I don't yet know what is the universal format but I think it will include separate facet elements based on my brief overview.

[gdmcbain]

Ah, yes, of course. Interoperability is another very important criterion, at least equal to size and speed, perhaps moreso. It is unfortunate that at the moment these goals conflict.

At the moment, maybe MSH 2.2 is the most widely accepted input format? (Is that supported by your brief overview?) By format I mean format that includes subdomains and boundaries, and accepted and understood by finite element systems. If that's correct, perhaps we should have a method or option to save an skfem.Mesh that records our subdomains and boundaries in the way that they are encoded by Gmsh so that they can be reread by other FEM packages? This is not exactly the way that meshio does it. Currently we have

scikit-fem/skfem/io/meshio.py

Lines 182 to 189 in a1496b6

	cell_data.update(mesh._encode_cell_data())
	mio = meshio.Mesh(
	mesh.p.T,
	cells,
	point_data=point_data,
	cell_data=cell_data,
	)

When MSH 2.2 was the main Gmsh output format, which it was for a very long time, Gmsh encoded subdomains and boundaries as what it called ‘physical tags’ of cells, every cell having exactly one physical tag, using (n − 1)-dimensional cells to represent the boundary-facets. When meshio reads these, it encodes them as cell_data. Gmsh MSH 2.2 (& 4.1) also has a separate concept of ElementData which meshio also encodes as cell_data, so meshio had just one concept corresponding to two different ones in Gmsh MSH. This two-to-one on meshio-read causes problems on inverting: the cell_data has to be classified as either physical tag or element-data, so the special "gmsh:physical" name was invented and used in cell_data. Currently we do decode this in

scikit-fem/skfem/io/meshio.py

Lines 105 to 121 in a1496b6

	try:
	elements_tag = m.cell_data_dict['gmsh:physical'][meshio_type]
	subdomains = {}
	tags = np.unique(elements_tag)
	def find_tagname(tag):
	for key in m.field_data:
	if m.field_data[key][0] == tag:
	return key
	return None
	for tag in tags:
	t_set = np.nonzero(tag == elements_tag)[0]
	subdomains[find_tagname(tag)] = t_set
	# find tagged boundaries
	if bnd_type in m.cell_data_dict['gmsh:physical']:

but not re-encode it in skfem.io.meshio.to_meshio. Should this go in _encode_cell_data? That is, create a "gmsh:physical" list of arrays? (Of course this encounters a difficulty when subdomains or boundaries overlap, which we have dealt with differently in #680 but which Gmsh dealt with by deprecating MSH 2.2 in favour of MSH 4.1, and we don't know how to write MSH 4.1 with more than one type of cell because we don't know how to write the new MSH 4.1 ‘entities’.)

Amongst higher-end finite-element packages, perhaps Exodus is more popular? What was your finding in your brief overview? It does require netCDF4; not impossible, but an additional requirement. I'm not sure about the level of support of Exodus in meshio. It might be possible to lift that with some careful PRs, or perhaps the necessary work can be done prior to meshio in skfem.io.meshio.

A basic difficulty is that meshio and VTK have no real concept of subdomain or boundary. meshio does how have cell_sets which do naturally correspond to subdomains and can also be used for boundaries via (n − 1)-dimensional cells, but nschloe/meshio#619 ‘unification of point / cell sets among formats’ seems stalled with nschloe/meshio#691 ‘Unifying point and cell tags’ having been replaced by (the less ambitious? less prescriptive?) nschloe/meshio#704 ‘more streamlined handling of point and cell data’.

But really my question here wasn't so much about the ultimate ‘universal format’ (as desirable and worthwhile a goal as that is), but rather just the immediate consequences of the splitting up of #693: whether

scikit-fem/docs/examples/ex27.py

Line 109 in e4a4b2b

self.mesh = skfem.io.json.from_file(Path(__file__).parent / 'meshes' / 'backward-facing_step.json')

is to be retained for now or the JSON to be deprecated and here replaced by something else (MSH 2.2 or MSH 4.1 straight from Gmsh, something exploiting #680, ...).

I don't think it's urgent to deprecate the JSON, but I recall the comment in #556 (2021-03-01)

I plan to deprecate our in-house JSON format and try to use something else (HDF5?) instead.

So for exx 24, 27, I have prepared the pure MSH 4.1 and 2.2 files which could be used, also the #680-style VTU, or it can be left as JSON indefinitely till something decisively better is found. Should the file-formats used in docs/examples/meshes be overhauled? Do we want to demonstrate multiple approaches? Rather than pushing a particular line here, I'm more seeking an opinion, a house style guide for scikit-fem examples.

[kinnala]

I don't have an opinion on which format to use. Moving away from JSON in the examples makes sense because of the new developments. I don't think there is an urgent need to remove it though.

[kinnala]

The approach in #680 was chosen so that it's compatible with many different formats. However, if we want to have a separate approach for, e.g., MSH2.1 via gmsh:physical, I think it would be fine to have a variant of io.meshio like io.msh21 which writes explicit facet elements. Purpose of the variant would be to have better interoperability with other tools. However, I'd like to verify this first by checking what kind of formats other open source tools support.

[gdmcbain]

For the specific example of exx 24, 27, I was thinking of regenerating the mesh as a nonuniform MeshQuad.init_tensor (compressed along the x- and y-lines through the reentrant corner, perhaps with an arc tanh mapping), deleting the cells that would be inside the step and then identifying the boundary-facets geometrically. But haven't gotten around to it… Then I was thinking of in those examples demonstrating that one can .save a mesh with identified boundary-facets and preceding all that with a try Mesh.load, except generate as said. And no, none of this is urgent.

[gdmcbain]

Yes, a separate .io module for a Gmsh-based format focused on interoperability sounds good. It'd be MSH 2.2 rather than 2.1. I presume that there was a 2.1 but whereas MSH 2.2 continues as ‘legacy’, 2.1 has just disappeared.

Besides Gmsh's MSH 2.2, FreeFEM reads MEDIT, including int-labelled nonoverlapping subdomains and boundaries, the latter again represented as (n − 1)-dimensional simplices—FreeFEM only supports simplicial meshes, no quads or hexas. I'm not sure that I ever attempted to write a mesh with subdomains or boundaries from FreeFEM.

FEniCS also reads MSH 2.2, via meshio and https://github.com/floiseau/msh2xdmf.

Generally there's a lot more emphasis on reading than writing; one can't do much FEA without a mesh but serialization is more specialized. Often the output is only for visualization, e.g. to VTK or XDMF, and doesn't record the subdomains or boundaries.