I/O from XDMF/HDF5 files in dolfin-x

dolfinx
32

So the issue here is that you are making certain assumptions when you use a built in mesh.

A “built in” mesh is created as any other mesh, as set of points and connectivities as two separate arrays

These cells and points are sent to a mesh partitioner of choice (default is SCOTCH), which partitions the mesh. The mesh keeps track of its input index (under mesh.topology.original_cell_index).
These are in turn used to read in MeshTags of various dimensions.

However, when a mesh is written to file, we do not revert it to its original indices. It is then written in a continuous fashion, as this is the most efficient.

This means that whenever you read in that mesh from file, the cells are shuffled before they are sent to the mesh partitioner.

This means that there is no direct mapping between the mesh generated with create_unit_square and the one read from file.

You can see that xdmf and adios2 treats the mesh similarly, by calling:


from mpi4py import MPI
import dolfinx
from dolfinx.fem import (Function, FunctionSpace, Expression)
from dolfinx.io import XDMFFile
import ufl
from petsc4py import PETSc
import adios4dolfinx

comm = MPI.COMM_WORLD
mesh = dolfinx.mesh.create_unit_square(comm, 10, 10, cell_type=dolfinx.mesh.CellType.triangle)
filename = f"./{comm.size}proc.bp"

### write ###
v_cg2 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, v_cg2)
u = Function(V)
x = ufl.SpatialCoordinate(mesh)
dfnx_expr = Expression(x, V.element.interpolation_points())
u.interpolate(dfnx_expr)
adios4dolfinx.write_mesh(mesh, "./mesh/square.bp")
adios4dolfinx.write_function(u, filename)


### write mesh ###
with dolfinx.io.XDMFFile(comm, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    top = xdmf.read_topology_data()
### read ###

with dolfinx.io.XDMFFile(comm, "mesh.xdmf", "r") as xdmf:
    mesh_xdmf = xdmf.read_mesh()

mesh_new = adios4dolfinx.read_mesh(comm, "./mesh/square.bp", "BP4", dolfinx.mesh.GhostMode.none)


print(comm.rank,
      "mesh input index", mesh.topology.original_cell_index, "\n",
      mesh_new.topology.original_cell_index, "\n",
      "\n",
      mesh_xdmf.topology.original_cell_index)


v_new = ufl.VectorElement("Lagrange", mesh_new.ufl_cell(), 2)


W = FunctionSpace(mesh_new, v_new)
w = Function(W)
adios4dolfinx.read_function(w, filename)
#adios4dolfinx.read_function(v, filename)
with XDMFFile(comm, f"./0_{comm.size}proc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh_new)
    xdmf.write_function(w)
with XDMFFile(comm, f"./1_{comm.size}proc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(u)
with XDMFFile(comm, f"./2_{comm.size}proc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(w)

with XDMFFile(comm, f"./3_{comm.size}proc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh_xdmf)
    xdmf.write_function(w)

Here, the third example (the mesh read from XDMF) gives the correct output.

2 Likes