Read xdmf file to do pyvista visualizations

I/O
,
5

Thank you for your answer.

I managed to create a VTK-based workflow, which works well for static functions. For time-dependent functions, xdmf is much easier to use and does store the mesh only once, which is nice.

The option to read xdmf in pyvista was requested in this Github discussion, and then later implemented (see docs).

However i have not managed to make it work. Here is an MWE (based on the heat-equation example from the tutorial):

# HEAT EQUATION

import pyvista
import ufl
import numpy as np
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, io, plot
from dolfinx.fem import locate_dofs_topological,dirichletbc
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.mesh import *```

# sim parameters
t = 0  
T = 0.1 
num_steps = 20
dt = T / num_steps

# mesh and functions
domain = create_unit_square(MPI.COMM_WORLD, 10,10)
V = fem.FunctionSpace(domain, ("Lagrange", 1))
u_n, uh = fem.Function(V), fem.Function(V)
def initial_condition(x, a=5): return np.exp(-a * (x[0]**2 + x[1]**2 + x[2]**2))
u_n.interpolate(initial_condition)
uh.interpolate(initial_condition)

# BC
facets = locate_entities_boundary(domain, dim=1,
                                    marker=lambda x: np.logical_or.reduce((
                                        np.isclose(x[0], 0.0),
                                        np.isclose(x[0], 1.0),
                                        np.isclose(x[1], 0.0),
                                        np.isclose(x[1], 1.0))))
dofsV = locate_dofs_topological(V, entity_dim=1, entities=facets)
bc = dirichletbc(0.0,dofsV, V)

# problem
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
a = u * v * ufl.dx + dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = (u_n + dt * f) * v * ufl.dx
bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form)

# solver
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# xdmf file
from dolfinx import io
xdmf = io.XDMFFile(domain.comm, "test.xdmf", "w")
xdmf.write_mesh(domain)

for i in range(num_steps):
    t += dt

    # update right hand side, apply BC
    with b.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(b, linear_form)
    apply_lifting(b, [bilinear_form], [[bc]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, [bc])

    # Solve and update sol
    solver.solve(b, uh.vector)
    uh.x.scatter_forward()
    u_n.x.array[:] = uh.x.array

    # write solution to xdmf
    xdmf.write_function(uh, t) 

xdmf.close()

After running this and producing test.xdmf, I would like to read it in pyvista. However, neither

import pyvista as pv
xdmf_reader = pv.XdmfReader("test.xdmf")
mesh = xdmf_reader.read()
mesh.plot()

nor

import pyvista as pv
reader = pv.get_reader("test.xdmf")
mesh = reader.read()
mesh.plot()

work. They give a segmenation fault.
(These two code snippets are inspired by the github discussion and the pyvista documentation).

Im am very grateful for any advice.