Reading functions from VTU files

[c-abird]

I'm trying to read meshes and functions (CG1) from VTU files using meshio (https://github.com/nschloe/meshio).

Reading 3d meshes works fine with the following:

import firedrake as fd
import meshio

mdata = meshio.read("some_file.vtu")
cells = np.asarray(mdata.cells_dict["tetra"], dtype=np.int32)
coords = np.asarray(mdata.points, dtype=np.double)
dm = PETSc.DMPlex().createFromCellList(3, cells, coords, comm=fd.COMM_WORLD)
mesh = fd.Mesh(dm, reorder=False)

However, I'm not able to reliably read function data from VTUs since Firedrake seems to either renumber the mesh vertices or the DOFs of the CG spaces. The following might or might not work:

V = fd.FunctionSpace(mesh, "CG", 1)
f = Function(V)
f_data = list(mdata.point_data.values())[0]
f.vector()[:] = fdata

Is there some way to get a map of the reordering in order to properly assign the function vector?

[wence-]

You would access the point data through the vtus cell to node map, and the firedrake function through its cell_node_map. Note that even with reordering off this won't work right in parallel. The generic way of doing this is via interpolation from external data (there's a section in the docs). I guess you have the vtu data from somewhere else...

[c-abird]

Thanks. I thought that interpolation would be the most robust way, but I'm trying to avoid installing the complete VTK suite which would AFAIK be required in order to do point evaluation on the VTU data.

Do I understand correctly that the DOFs in Firedrake are renumbered such that they appear in ascending order in the cell_node_map? I'm now deriving a mapping from the VTU cell to node map as follows

# get unique vertex ids in their order of appearance of the VTU cell to node map
_, idx = np.unique(cells.reshape(-1), return_index=True)
vertex_map = cells.reshape(-1)[np.sort(idx)]
f.vector()[:] = fdata[vertex_map]

which seems to work quite well.

[pefarrell]

In case it's of interest, here's how to use vtktools.py (from Fluidity) to read in data from a VTU file:

import numpy
from firedrake import *
import vtktools

def make_vtu():
    mesh = UnitSquareMesh(3, 3, quadrilateral=True)
    V = FunctionSpace(mesh, "CG", 1)
    (x, y) = SpatialCoordinate(mesh)
    u = interpolate(x + y, V)
    u.rename("Solution")
    File("output.pvd").write(u)

make_vtu()

mesh = UnitSquareMesh(5, 5, diagonal="crossed")
V = FunctionSpace(mesh, "CG", 2)

W = VectorFunctionSpace(mesh, "CG", 2)
X = interpolate(SpatialCoordinate(mesh), W)

vtu = vtktools.vtu("output_0.vtu")
reader = lambda X: vtu.ProbeData(numpy.c_[X, numpy.zeros(X.shape[0])], "Solution").reshape((-1,))

u = Function(V)
u.dat.data[:] = reader(X.dat.data_ro)

(x, y) = SpatialCoordinate(mesh)
w = interpolate(x + y, V)
print("||u - w||: ", norm(u - w))

vtktools.py.gz

[c-abird]

Thank you @pefarrell