Write cg1 and dg0 functions to VTKHDF

[mleoni-pf]

This PR introduces a convenience function to write CG1 and DG0 functions to VTKHDF.

[schnellerhase]

The decision to split up into different functions depending on function space looks problematic, it will not scale. A generic write_function would be favourable.

[mleoni-pf]

I was inspired by the Python interface having separate "write point data" and "write cell data" but I realize that that's not quite the same... I agree that a generic write_function would be better but I am not knowledgeable enough on the VTKHDF format to know how to write mixed or exotic elements...
An API-saving option could be to error out for unsupported function spaces; another one, used somewhere else in the repo though unsatisfactory, could be to interpolate whatever function is passed to CG1 or DG0 before saving. What are your thoughts?

[schnellerhase]

I think supporting subsets of elements is fine. Good to add tests, which check for a proper assertion being raised on such invocations.

[chrisrichardson]

I think in this case, raise an exception for unsupported function spaces.

[chrisrichardson]

This PR makes implicit assumptions about the dofmap for P1 and DG0 functions (i.e. that they are simply the same as the geometry).
I agree with @schnellerhase that there should be a common interface, and also there should be some tests that run with MPI.

[mleoni-pf]

This PR makes implicit assumptions about the dofmap for P1 and DG0 functions (i.e. that they are simply the same as the geometry). I agree with @schnellerhase that there should be a common interface, and also there should be some tests that run with MPI.

I'm interested in learning how to relax that assumption, if you can point me to whatever functions we have to map the geometry dofmap into the P1 and DG0 dofmaps.

As far as tests are concerned, I didn't know how to write a test for a writing without being able to read back but in #4027 I added a test that includes a full write-read-compare cycle.

[schnellerhase]

Have a look at add_function in xdmf_function.cpp, in particular how data_values is constructed.

[jorgensd]

This PR makes implicit assumptions about the dofmap for P1 and DG0 functions (i.e. that they are simply the same as the geometry). I agree with @schnellerhase that there should be a common interface, and also there should be some tests that run with MPI.

DG-0 functions are always the same as the mesh topology (i.e. dof index i corresponds to cell i) on any given process. This can of course be adapted by extracting the underlying dofmap of the DG-0 space, i.e. V.dofmap.list.

Of course, for other functions than P1, i.e. generalized geometry that is P-th order lagrangian, I would go down the route we use in VTKFile and VTXWriter: https://github.com/FEniCS/dolfinx/blob/main/cpp/dolfinx/io/ADIOS2Writers.h#L224-L285
which create a mesh geometry based on the dofmap of the underlying Lagrange space, and then writes the data to file on these.

I would not go down the route of what we do in XDMFFile, which is restricting ourselves to writing to the mesh nodes (which involves a non-trivial interpolation step if function space!=geometry space).

[mleoni-pf]

@jorgensd: thanks for the input, you are right that for higher-order P elements the other approach would be better but at the moment I think that the most critical use case is P1 functions because this functionality is mainly for interoperability between different simulations and different programs and, if needed, a higher-order function can be interpolated to a P1 function on a much finer mesh [and vice versa when reading it].

@chrisrichardson: in a pair programming session late last year, @schnellerhase and I adapted the implementation from xdmf_function and addressed the issues that you both raised. I requested a re-review!

[schnellerhase]

I would not go down the route of what we do in XDMFFile, which is restricting ourselves to writing to the mesh nodes (which involves a non-trivial interpolation step if function space!=geometry space).

VTKHDF supports arbitrary cell types - https://gitlab.kitware.com/keu-public/vtkhdf/vtkhdf-scripts/-/tree/main/tutorial?ref_type=heads#creating-a-basic-unstructured-grid. Therefore, if I understand correctly, we can aim for an actual exact high order output in the long run, i.e. supporting cell types such as VTK_HIGHER_ORDER_TETRAHEDRON. If one aims for this, the dofmap based construction is not exactly what we need for the higher order case, as we do not want to rely on locally linear interpolations, right?

[jorgensd]

I would not go down the route of what we do in XDMFFile, which is restricting ourselves to writing to the mesh nodes (which involves a non-trivial interpolation step if function space!=geometry space).

VTKHDF supports arbitrary cell types - https://gitlab.kitware.com/keu-public/vtkhdf/vtkhdf-scripts/-/tree/main/tutorial?ref_type=heads#creating-a-basic-unstructured-grid. Therefore, if I understand correctly, we can aim for an actual exact high order output in the long run, i.e. supporting cell types such as VTK_HIGHER_ORDER_TETRAHEDRON. If one aims for this, the dofmap based construction is not exactly what we need for the higher order case, as we do not want to rely on locally linear interpolations, right?

Vtkhdf supports the arbitrary Lagrange element we use within VTKFile and VTXWriter unfortunately it doesn’t currently support the DGCell construction (only exodus supports this), which would give us a split between geometry and function space.

for now I would aim for reading and writing point data, ie data from the geometry function space, meaning that if the geometry is based on a second order Lagrange element, we would write/read second order functions.