load auxiliary data with meshes

[gdmcbain]

This is kind of the opposite of #268, the removal of the old .external attribute of a loaded Mesh, which served as a catch-all for whatever else had been read from the file and wasn't immediately interpretable as part of the mesh. That attribute was no longer populated by Mesh.load and so was dropped #269.

However, there are cases in which it makes sense to load auxiliary data along with a mesh, e.g. fields that are defined on the points or cells of the mesh; e.g.

• reloaded a saved solution for continuation (transient evolution, numerical continuation, &c.)
• importing a solution from another code
• continuously variable coefficients (not so suited to modelling with subdomains as in exx 17, 26, 28)

or subsets of the points or cells; e.g.

• subdomains or boundaries Mesh.save subdomains and boundaries #261

meshio.Mesh does have attributes point_data, cell_data, and field_data and has more recently added:

• node_sets and element_sets for Exodus (only)
• point_tags and cell_tags for MED (only)

Ignoring these latter for the moment (I would hope that they'd be standardized in meshio eventually), should Mesh.load return both a mesh and some kind of container of auxiliary data? Maybe return Tuple[Mesh, MeshData] where MeshData is an object like a meshio.Mesh without the .points or .cells? I guess the .cell_data and .field_data would have the 'gmsh:physical' items popped.

I'm not sure that auxiliary data belongs in an object with the mesh (which should just be geometric and topological data), although it does make sense to have data that is associated with a particular mesh, and does make sense to be able to store mesh and auxiliary data together in a file (as many formats allow).

[kinnala]

When you have one value per node you also implicitly think of piecewise linear interpolation between the nodes which corresponds to ElementTriP1 or similar.

I think it makes sense to read this kind of information from the files. However, there is a similar kind of "vagueness" than in the solution vector returned from solve: it doesn't really mean anything unless you know which Element and which Mesh it corresponds to.

I'm trying to think what kind of container would be sufficient but simple enough.

[gdmcbain]

However, there is a similar kind of "vagueness" than in the solution vector returned from solve: it doesn't really mean anything unless you know which Element and which Mesh it corresponds to.

Yes. Or which Basis; that contains a mesh and an element.

I had suggested something like this in nschloe/meshio#630 (comment) but it wasn't well received

the project name should be changed! It is currently mesh-io and not fem-io.

I do like that idea myself, though I'm also still

trying to think what kind of container would be sufficient but simple enough.

Maybe besides a str for a name, an Enum for the basis? (Or rather Element which with the mesh almost implies a basis.) That is awkward because then the data can no longer be held as Dict[str, ndarray]…

Besides degrees of freedom, the other kind of vector which is passed around in finite element analysis is the DiscreteField; i.e. values at quadrature points. I think this is popular for stress in elasticity.

[kinnala]

This is now pretty topical w.r.t #556 . For now I decided to make the new mesh classes as backwards-compatible as possible while still being "immutable".

However, I'm now also thinking that subdomains and boundaries should be inside another object and not inside Mesh. I'm willing to remove those from Mesh at some point but I think it should probably be 4.0.0 and done in smaller steps.

First we could simply introduce the new objects, be it MeshData or similar, and start filling those when calling skfem.io directly and not through Mesh.load

[gdmcbain]

However, I'm now also thinking that subdomains and boundaries should be inside another object and not inside Mesh

Yes, as discussed 2021-01-15:

In the same way as class Mesh might be divided up into class Topology and class Geometry (or hopefully with better names), perhaps the boundaries attribute (and similarly subdomains) might be separated from the object so that another new class (‘Domain’?), again immutable, might have an immutable Mesh as an attribute and also boundaries and subdomains, which might be implemented as namedtuples rather than dicts so that they're less mutable too.

Yes, I think this a good idea and the small steps towards 4.0.0 too, yes.

[gdmcbain]

The name ‘MeshData’ suggests to me more data that is defined on a mesh (functions from points or cells or facets or edges to ℝ or ℝ^k or whatever), like the point_data and cell_data of meshio. Subdomains can be defined in this way (a function from cells to bool), and maybe boundaries (at a stretch, but I don't have an immediate suggestion, since a cell might have zero, one, or more than one facet on a boundary), but more natural I think is subsets of the indices of cells and facets, as is done now in skfem.Mesh.

On mesh data in the sense of point_data and cell_data, I'm still thinking (long term) about the exchange of 2020-01-09 in nschloe/meshio#630:

@gdmcbain:

But, I wonder, are we indeed only interested in their geometrical shapes? What about the data that's currently stored in node_data and cell_data? Isn't that more appropriately associated with the degrees of freedom of a finite element space?

@tianyikillua:

...Yes...but if we pursue further in this direction

1. we will be reinventing a new finite element data format...
2. the project name should be changed! It is currently mesh-io and not fem-io.

Those two points are well taken but I do still sometimes think that I do want a femio, I do want to be able to store the values of the degrees of freedom of a finite element function as such and not as point_data or cell_data. I guess this would currently have to go in field_data in meshio for the general case with the corresponding element-type stored along side it so that the basis can be reconstructed from the mesh and the element. Of course I would very much like to avoid inventing a new finite element data format and so I continue searching for an existing one like this.