Time-dependent mesh-output with higher order elements

I/O
, , ,
1

Hi everyone,

I want to use some ALE-method with remeshing, but I have problems writing my velocity field, which is composed of second order Lagrange-elements, in a single output file.

I have tried all 3 types of writers (dolfinx.io.VTKFile, dolfinx.io.VTXWriter, dolfinx.io.XDMFFile), but each of them has its drawbacks in this case:

  1. dolfinx.io.VTKFile is here practically the best one. The mesh geometry/topology can be updated every time step and you can write second order elements directly to the output files. It’s just annoying how many output files you produce, if you have many time steps. Would it be possible to move everything except the pvd-file (vtuand pvtu files) into a subdirectory of the output directory to avoid the clutter of files there? Further the simple vtu Format seems to be worse than the other two formats.
  2. dolfinx.io.VTXWriter seems to work, if you initialize the writer with the argument mesh_policy=VTXMeshPolicy.update. Also quadratic elements and changes in the mesh’s geometry are allowed, but you cannot change the topology of the mesh, e. g. applying remeshing or mesh refinment, (or paraview (Version 6.1.1) fails in the visualization of the topology change. I’m not sure since I don’t know the details of the file format). After each change in topology you have to create a new output file, which is not what I want for the visualization in paraview.
  3. dolfinx.io.XDMFFile also seems to work using the suggestion of the discussion in https://fenicsproject.discourse.group/t/update-mesh-over-time/13827, but it seems, that the xdmf-format is not ideal for changes in the geometry or topology of the mesh and the implementation is quite tricky (the solution in the post above suggest renaming the mesh every time step and setting mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']" as argument, when calling dolfinx.io.XDMFFile.write_function). Further the xdmf-format does not support quadratic Lagrange elements, so higher order elements must be interpolated or projected into a lower order space. This works, but then you loose some information.

To summarize this into a single question: Is there an output file format, that supports quadratic elements, changes in the mesh’s geometry and topology in time and do not blow up your output directory with files?

In my opinion the simplest solution would be using dolfinx.io.VTKFile as writer, but allow to create an additional subdirectory for all the vtuand pvtu files. You just need to add the (relative) path to this subdirectory to the filenames in your pvd-file.

Thanks in advance.

2

Long term, we hope that the optimal (and final format) would be VTKHDF) which combines the best of VTKFile and XDMFFile into one format.
I have for instance made writers and readers for this data at:

which can store point data and cell data.
Functionality is somewhat illustrated in:

If you want to have higher order fields, you can use dolfinx.fem.interpolate_geometry dolfinx.fem — DOLFINx Python 0.11.0.post0
to make the mesh higher order, and store data on it.
the idea would be:

from mpi4py import MPI
import dolfinx
import io4dolfinx
import numpy as np
import ufl
def ref(x):
    return x[0] + 2*x[1]


# create linear grid
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

# Create quadratic function
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
u = dolfinx.fem.Function(V)
u.interpolate(ref)

# Reconstruct second order mesh and interpolate function onto it
second_order_mesh = io4dolfinx.reconstruct_mesh(mesh, 2)
W = dolfinx.fem.functionspace(second_order_mesh, ("Lagrange", 2))
num_cells_local = second_order_mesh.topology.index_map(second_order_mesh.topology.dim).size_local + second_order_mesh.topology.index_map(second_order_mesh.topology.dim).num_ghosts
z = dolfinx.fem.Function(W)
z.interpolate(u, cells0=np.arange(num_cells_local, dtype=np.int32), cells1=np.arange(num_cells_local, dtype=np.int32))

and then store second_order_mesh and z to file.
Here you can replace io4dolfin.reconstruct mesh with interpolate_geometry as it was the prototype implementation of interpolate_geometry.

Time dependent topology: VTKHDF currently has a bug for time dependent topology, ref: Verifying connection...