Visualize mesh, vector field in paraview; export/import vector field with a parameter in DOLFINx

Hello I am using paraview to visualize my mesh and vector field solution with a parameter. Now I can export the mesh to a .xdmf file and import using dolfinx.io. But I can’t visualize it in paraview. This post seems to have the same issue as I am and solved by extract blocks. While I cannot click this option (it’s not clickable) when I open my xdmf file in paraview. I wonder what’s the solution in this case.

I have a parameter in my solution as a vector field, say it’s time. And I obtain a sequence of solution for different time. The way in my code to export this field could generate some file. But I don’t know how to visualize it in paraview or import it back to DOLFINx. What is the best way to do this?

import basix.ufl
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.nls.petsc
import dolfinx.log
import dolfinx.plot
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
import matplotlib.pyplot as plt
from petsc4py.PETSc import ScalarType
import pyvista
from pathlib import Path

# Mesh is only parametrized by its domain size 'a' and discretized size 'lcar'.
# First generate mesh

a = 5.0
lcar = 0.3

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, lcar)
p1 = gmsh.model.geo.addPoint(0.0, a, 0.0, lcar)
p2 = gmsh.model.geo.addPoint(0.0, -a, 0.0, lcar)
c0 = gmsh.model.geo.addCircleArc(p1, p0, p2)
c1 = gmsh.model.geo.addCircleArc(p2, p0, p1)
boundary = gmsh.model.geo.addCurveLoop([c0, c1])

domain = gmsh.model.geo.addPlaneSurface([boundary])

p_center = gmsh.model.geo.addPoint(0, 0, 0, lcar)

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [c0, c1], 1)
gmsh.model.addPhysicalGroup(0, [p_center], 2)
gmsh.model.addPhysicalGroup(2, [boundary], 0)

gmsh.model.mesh.embed(0, [p_center], 2, domain)

gmsh.model.mesh.generate(2)

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()

boundary_facet = boundaries.indices[boundaries.values == 1]
bb_tree = dolfinx.geometry.bb_tree(mesh, mesh.topology.dim)

V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 2))
uxuy = dolfinx.fem.Function(V)
uxuy.interpolate(lambda x: [x[0] + 1, (x[1])**2 - x[0]])

with dolfinx.io.XDMFFile(mesh.comm, "results/mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

with dolfinx.io.XDMFFile(mpi4py.MPI.COMM_WORLD, "results/mesh.xdmf", "r") as xdmf:
       mesh1 = xdmf.read_mesh(name="mesh")

vtx_u = dolfinx.io.VTXWriter(mesh.comm, "results/function/poiseuille_u.bp", uxuy, engine="BP4")
vtx_u.write(5.1)

Hello so far I haven’t found anyway to resolve the issue. However it could be bypassed, if I don’t use paraview for visualization and either print or append the array of the solution to another array, and eventually collect all solutions or print out the big array. So if you have any idea to do it in a better way, please leave your message.

I can open all files (mesh.xdmf and poiseuille_u.bp) in paraview 5.11.2 on Debian.

Can you try upgrading paraview to the newest version you have access to?

Hello Francesco thank you very much for your reply.

Could you tell me how exactly you open them? I am using 5.11.2 on Mac.

For mesh there are two separate files mesh.h5 and mesh.xdmf and options for readers in paraview for both. I tried multiple combinations and didn’t work out.

For the solution there is a folder called poiseuile_u.bp. If I drag the whole folder into paraview and it usually crashed…

I open paraview, click on “File → Open” and:

  • for the mesh, open the xdmf file (not the h5 one) and choose Xdmf3ReaderS (the one I get suggested by default, but the one with T would work too)
  • for the solution, select the name.bp folder and click “OK” (do not go inside the folder) and choose ADIOS2VTXReader (the one which is NOT suggested by default)

Hello thank you very much for your reply now its not crashing but only a disk appears…

This is for the mesh one:


I cannot click extract block either before or after I click apply.

Then I opened the folder while keep the mesh. Now I can click extract block. Nothing happens.

I delete the mesh.

This is so mysterious.

See the following for writing a solution with a varying parameter. And it can be played in paraview by clicking the play button. I still don’t know how to read the solution. There are old related posts 1 and 2.

vtx_u = dolfinx.io.VTXWriter(mesh.comm, "results/function/poiseuille_u.bp", uxuy, engine="BP4")
t = 1e-10
while t < 5.05:
    uxuy.interpolate(lambda x: [x[0] * t, x[1]])
    vtx_u.write(t)
    t += 0.1

Then

# Read function from file???

Here’s a screenshot of successful visualization.

I am lost: what’s the actual question? It seems to me that you are able to import the files, and (upon changing from “Solid color” to “f”) you are also able to visualize the solution.

Oh sorry. I mean to import it in dolfinx. Suppose I have obtained the solution and saved it in a file. Later I use dolfinx and want to read the file to import it in my program, maybe for some future use. (Not in paraview, but read the saved function in dolfinx)

For that operation (which we typically call “checkpointing”) see GitHub - jorgensd/adios4dolfinx: Interface of ADIOS2 for DOLFINx

Thank you very much!