read multiple XDMF Grids

[gdmcbain]

This is for #568.

This is atop #621 and is the second part of the splitting up of #610.

[gdmcbain]

With this, one can successfully meshio.read('ref.xdmf'), with the example ref.xdmf from #620. That XDMF does not contain multiple grids but meshio thinks it does, thinking any XML child of a Domain must be a Grid.

One can also meshio.read multiple grids like the XDMF example from #568.

[nschloe]

So what does this do when multiple grids are present in a file? Merge them into one?

[gdmcbain]

Into one meshio.Mesh, yes. That the cells come from different grids is remembered in the .cell_sets attribute. This could be used either for different cell-types, to avoid TopologyType="Mixed", or for the same cell types to distinguish subdomains, which implements tagged subsets of elements #290 instead of relying on a MEDIT/MSH2-style cell_data tag.

[gdmcbain]

... That's for the cells. Multiple grids are implemented here for their Topology only, their cells. They are assumed to share a common Geometry, all cells being supported by a single set of points. This isn't a restriction of XDMF itself but unlike multiple sets of cells, I don't see much application for multiple sets of points.

[nschloe]

That the cells come from different grids is remembered in the .cell_sets attribute

Not sure if this is way we should do it. It seems we're mixing the notions of sets/tags and meshes here. If there are multiple meshes in a file, why not return multiple mesh objects, too?

[nschloe]

to avoid TopologyType="Mixed"

I wouldn't want to work around XDMF this way. If there are multiple cell types in a mesh, then it seems that -- unfortunately! -- the XDMF standards wants to see "Mixed" cells. That's slow to read in Python, but yeah, the standard isn't perfect. Better use something else then (e.g., HMF).

[nschloe]

I don't see much application for multiple sets of points.

Indeed, most codes divide the mesh into submeshes by cells. Some are point-centered though, point_sets make sense here as well.

[gdmcbain]

If there are multiple meshes in a file, why not return multiple mesh objects, too?

So that, in XDMF terms, multiple Topologies can share a Geometry; in meshio terms, Mesh.points is a single array but Mesh.cells is multiple.

That's for what I imagine is the more common use-case, when the different Grids have Topologies with different cell types. But then what should the reader do if the cell-types match? (Again with shared points.) The interpretation of the intent of such a XDMF Domain that I went with is subdomaina. I can't think of another interpretation, which would have the same cell-types and points but different connectivity.

[gdmcbain]

I wouldn't want to work around XDMF this way. If there are multiple cell types in a mesh, then it seems that -- unfortunately! -- the XDMF standards wants to see "Mixed" cells. That's slow to read in Python, but yeah, the standard isn't perfect.

No, it is the standard itself that suggests the present approach.

It is valid to have DataItem elements to be direct children of the Xdmf or Domain elements. This could be useful if several Grids share the same Geometry but have separate Topology:

http://xdmf.org/index.php/XDMF_Model_and_Format#Geometry

[gdmcbain]

Indeed, most codes divide the mesh into submeshes by cells. Some are point-centered though, point_sets make sense here as well.

This sounds interesting but different. Could you send an example? Aren't there points along the interface?

I had assumed that point_sets were for something like low-level implementation of essential boundary condition.

Anyway this doesn't seem to be what the XDMF authors are getting at with ‘several Grids share the same Geometry but have separate Topology‘.

[gdmcbain]

Not sure if this is way we should do it. It seems we're mixing the notions of sets/tags and meshes here. If there are multiple meshes in a file, why not return multiple mesh objects, too?

To get another idea of how to interpret an XDMF file containing several Grids sharing the same Geometry but with separate Topology, I opened the example from #568 in ParaView 5.7.0 and used 'Save Data' to export it as EXODUS. 568.zip

The ncdump of 568.e contains:

	num_nodes = 12 ;
	num_elem = 3 ;
	num_el_blk = 2 ;
	num_el_in_blk1 = 1 ;
	num_nod_per_el1 = 4 ;
	num_el_in_blk2 = 2 ;
	num_nod_per_el2 = 3 ;

That is: ParaView has doubled the points! Further:

 coordx = 0, 1, 2, 2, 1, 0, 0, 1, 2, 2, 1, 0 ;

 coordy = 1, 0, 0, 2, 2, 3, 1, 0, 0, 2, 2, 3 ;

 coordz = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;

and

 connect1 =
  2, 3, 4, 5 ;

 connect2 =
  7, 8, 11,
  7, 11, 12 ;

ParaView's behaviour in doubling the points is a bit disappointing, but it is exporting the Grids as EXODUS element blocks in a single EXODUS mesh rather than as multiple EXODUS meshes.

Were ParaView able to realize that the second Geometry was the same (a second reference to the same Geometry), it should only have listed the first six points and subtracted six from the ints in connect2. Subtracting another 1 (for zero-based meshio & XDMF) would then give the same four (1 2 3 4) and six (0 1 4, 0 4 5) connectivity ints as in the two original Topology DataItems.

The interpretation of the present PR is: .points:

[[0. 1. 0.]
 [1. 0. 0.]
 [2. 0. 0.]
 [2. 2. 0.]
 [1. 2. 0.]
 [0. 3. 0.]]

and .cells:

{'quad': array([[1, 2, 3, 4]], dtype=int32), 
 'triangle': array([[0, 1, 4],
       [0, 4, 5]], dtype=int32)}

[nschloe]

I thought about it some more and I'm still having a hard time convincing myself this is right approach. When seeing multiple <Grid>s in an XDMF file, I can't help but reading this as multiple meshes. I'm sure they are supposed to be used this way. meshio as the shortcoming that, right now, it only supports one mesh for every file. I'd be perfectly happy with returning a list of meshes each containing the points and cells specified in the file, but not with merging multiple <Grid>s into one meshio Mesh.

[gdmcbain]

Multiple meshio meshes aren't connected; they don't share points ('geometry ' in XDMF). The purpose of the present approach is to support the case discussed in XDMF Model and format:

This could be useful if several Grids share the same Geometry but have separate Topology:

Practical use case: a boundary value problem. Say 2D. One defines the geometry in (py)gmsh and marks the domain as Physical Surface and the boundary as Physical Curve. On meshing this in Gmsh, one gets triangles and lines for the domain and boundary, respectively. These cells of different type ('topology' in XDMF) share the same points; their connectivity data indexes into a common array of points. This is naturally represented in meshio with a single Mesh with a single array of points but multiple arrays of cells. This is how meshio has always represented such oases. XDMF can do this too, as stated in the specification, but meshio can't read such XDMF without this PR.

[nschloe]

In the example you're linking, one has two <Grid>s which happen to share a <Geometry>. What if you provide a mesh where you just have two <Grid>s with different <Topology> and <Geometry>, say the output of a mesh optimizer? One would have to check if the <Geometry>s coincide. If they do, the two <Grid>s are merged down to one meshio.Mesh, otherwise they are two separate meshio.Meshes. That's what bothers me.

It'd be much clearer to say that one <Grid> is one meshio.Mesh.

It seems that the intended mixed-element mesh in XDMF is the <Grid> with Mixed element type. Unfortunately, Python one can't read that very fast.

[gdmcbain]

One would have to check if the s coincide. If they do, the two s are merged down to one

Yes, it should be checked that the Geometries coincide. This WIP doesn't do that yet. I think this should be easy though as the Geometry refers to a DataItem which has an XML name and these have to be unique so we only have to store and compare these nsmes rather then the numerical arrays of points. When I saw the recent use of the point_sets attribute in another meshio module, I thought that that looked like a natural home for these names.

otherwise they are two separate meshio.Meshes

Or one Mesh with point_sets. If the XDMF truly intended two Meshes, why wouldn't it use two Xdmf?

I don't get the mesh optimizer example. Could you elaborate? As in optimesh? I haven't used that yet and don't know what data structures are required. Without having done so, I thought it would be a function Mesh -> Mesh.

[gdmcbain]

It seems that the intended mixed-element mesh in XDMF is the with Mixed element type.

No, the spec explicitly refers to cases in which 'several Grids share the same Geometry but have separate Topology:'. The XDMF format is flexible, compared to say MEDIT, it enables more than one way of representing the same thing. The Mixed option is valid but not a good fit for meshio. Meshio should read it but should also feel free to avoid writing it whenever there is an alternative that is value XDMF. That's what I'm trying to get to. This might become clearer once I get to writing Mesh.cells as XDMF. This becomes easy using the idiom xdmf#15, except that without the present PR, meshio can't read it back.

[nschloe]

No, the spec explicitly refers to cases in which 'several Grids share the same Geometry but have separate Topology:'

I read this as "several meshes share the same Geometry...".

Or one Mesh with point_sets. If the XDMF truly intended two Meshes, why wouldn't it use two Xdmf?

optimesh indeed, or moving mesh methods. There are a bunch of applications where you have groups of meshes with different geometry and topology; XDMF allows you to store them in one file going like

<Grid>
  <Topo></Topo>
  <Geo></Geo>
</Grid>
<Grid>...</Grid>
<Grid>...</Grid>
...

So in this case, we'd probably want to return a list of meshes from meshio. Would you agree?