Two questions on VTK output and Pyvista

I/O
1

In the following test code, I am trying to implement VTXWriter output (for Paraview, etc.) as well as Pyvista output. My first problem is that VTXWriter complains that ADIOS2 is not installed. How can I make this work? ADIOS appears to be installed on my system.

Second, I am trying to get Pyvista to work on multi-process MPI jobs. What am I missing?

from mpi4py import MPI
import numpy as np
from dolfinx import default_scalar_type, plot
from dolfinx.fem import Function, functionspace
from dolfinx.mesh import CellType, create_box, locate_entities, compute_midpoints
from dolfinx.io import XDMFFile

msh = create_box(MPI.COMM_WORLD, ((0, 0, 0), (1, 1, 1)), (25, 25, 25), CellType.tetrahedron)

V = functionspace(msh, ("Nedelec 1st kind H(curl)", 1))
u = Function(V, dtype=default_scalar_type)

tdim = msh.topology.dim
cells0 = locate_entities(msh, tdim, lambda x: x[0] <= 0.5)
cells1 = locate_entities(msh, tdim, lambda x: x[0] > 0.5)

u.interpolate(lambda x: np.vstack((x[0], x[1], x[2])), cells0)
u.interpolate(lambda x: np.vstack((x[0] + 1, x[1], x[2])), cells1)

gdim = msh.geometry.dim
V0 = functionspace(msh, ("Discontinuous Lagrange", 1, (gdim,)))
u0 = Function(V0, dtype=default_scalar_type)
u0.interpolate(u)

try:
    from dolfinx.io import VTXWriter
    with VTXWriter(msh.comm, "output_nedelec.bp", u0, "bp4") as f:
        f.write(0.0)
except ImportError:
    print("ADIOS2 required for VTX output")

topology, cell_types, geometry = plot.vtk_mesh(V0)
num_cells_local = msh.topology.index_map(msh.topology.dim).size_local
num_dofs_local = V0.dofmap.index_map.size_local*V0.dofmap.index_map_bs
num_dofs_per_cell = topology[0]
topology_dofs = (np.arange(len(topology)) % (num_dofs_per_cell+1)) != 0
global_dofs = V0.dofmap.index_map.local_to_global(topology[topology_dofs].copy())
topology[topology_dofs] = global_dofs

root = 0
global_topology = msh.comm.gather(topology[:(num_dofs_per_cell+1)*num_cells_local], root=root)
global_geometry = msh.comm.gather(geometry[:V0.dofmap.index_map.size_local,:], root=root)
global_ct = msh.comm.gather(cell_types[:num_cells_local])
global_vals = msh.comm.gather(u0.x.array[:num_dofs_local])

if(MPI.COMM_WORLD.rank == root):
    try:
        import pyvista
        root_geom = np.vstack(global_geometry)
        root_top = np.concatenate(global_topology)
        root_ct = np.concatenate(global_ct)
        root_vals = np.concatenate(global_vals)

        cells, types, x = plot.vtk_mesh(V0)
        grid = pyvista.UnstructuredGrid(root_top, root_ct, root_geom)
        grid.point_data["u"] = root_vals

        pl = pyvista.Plotter(shape=(2, 2))

        pl.subplot(0, 0)
        pl.add_text("magnitude", font_size=12, color="black", position="upper_edge")
        pl.add_mesh(grid.copy(), show_edges=True)

        pl.subplot(0, 1)
        glyphs = grid.glyph(orient="u", factor=0.08)
        pl.add_text("vector glyphs", font_size=12, color="black", position="upper_edge")
        pl.add_mesh(glyphs, show_scalar_bar=False)
        pl.add_mesh(grid.copy(), style="wireframe", line_width=2, color="black")

        pl.subplot(1, 0)
        pl.add_text("x-component", font_size=12, color="black", position="upper_edge")
        pl.add_mesh(grid.copy(), component=0, show_edges=True)
    
        pl.subplot(1, 1)
        pl.add_text("y-component", font_size=12, color="black", position="upper_edge")
        pl.add_mesh(grid.copy(), component=1, show_edges=True)

        pl.view_xy()
        pl.link_views()

    # If pyvista environment variable is set to off-screen (static)
    # plotting save png
        if pyvista.OFF_SCREEN:
            pyvista.start_xvfb(wait=0.1)
            pl.screenshot("uh.png")
        else:
            pl.show()
    except ModuleNotFoundError:
        print("pyvista is required to visualise the solution")

The returned error:

 File "/home/bill/Cluster/Fenics2020/Interplot.py", line 56, in <module>
    grid.point_data["u"] = root_vals
  File "/home/bill/.local/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 223, in __setitem__
    self.set_array(value, name=key)
  File "/home/bill/.local/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 623, in set_array
    vtk_arr = self._prepare_array(data, name, deep_copy)
  File "/home/bill/.local/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 810, in _prepare_array
    raise ValueError(f'data length of ({data.shape[0]}) != required length ({array_len})')
ValueError: data length of (1125000) != required length (375000)
--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------

It looks like there is some problem with passing the three vector components to Pyvista.

2

How did you install ADIOS2, and how did you install DOLFINx? Was adios2 installed prior to DOLFINx? DOLFINx needs to build with ADIOS2 to be able to use the ADIOS2 interfaces

You do not have point data in your mesh, but vector valued quantities. See for instance: Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial on how to get vector valued quantities into a pyvista grid. (thus replace grid.point_data with grid["u"]= and make sure that root_vals has the expected shape

3

Hi Dokken,
Thanks! Pyvista seems to work making the modifications as you suggested - viz.

from mpi4py import MPI
import numpy as np
from dolfinx import default_scalar_type, plot
from dolfinx.fem import Function, functionspace
from dolfinx.mesh import CellType, create_box, locate_entities, compute_midpoints
from dolfinx.io import XDMFFile

msh = create_box(MPI.COMM_WORLD, ((0, 0, 0), (1, 1, 1)), (25, 25, 25), CellType.tetrahedron)

V = functionspace(msh, ("Nedelec 1st kind H(curl)", 1))
u = Function(V, dtype=default_scalar_type)

tdim = msh.topology.dim
cells0 = locate_entities(msh, tdim, lambda x: x[0] <= 0.5)
cells1 = locate_entities(msh, tdim, lambda x: x[0] > 0.5)

u.interpolate(lambda x: np.vstack((x[0], x[1], x[2])), cells0)
u.interpolate(lambda x: np.vstack((x[0] + 1, x[1], x[2])), cells1)

gdim = msh.geometry.dim
V0 = functionspace(msh, ("Discontinuous Lagrange", 1, (gdim,)))
u0 = Function(V0, dtype=default_scalar_type)
u0.interpolate(u)

try:
    from dolfinx.io import VTXWriter
    with VTXWriter(msh.comm, "output_nedelec.bp", u0, "bp4") as f:
        f.write(0.0)
except ImportError:
    print("ADIOS2 required for VTX output")

topology, cell_types, geometry = plot.vtk_mesh(V0)
num_cells_local = msh.topology.index_map(msh.topology.dim).size_local
num_dofs_local = V0.dofmap.index_map.size_local*V0.dofmap.index_map_bs
num_dofs_per_cell = topology[0]
topology_dofs = (np.arange(len(topology)) % (num_dofs_per_cell+1)) != 0
global_dofs = V0.dofmap.index_map.local_to_global(topology[topology_dofs].copy())
topology[topology_dofs] = global_dofs

root = 0
global_topology = msh.comm.gather(topology[:(num_dofs_per_cell+1)*num_cells_local], root=root)
global_geometry = msh.comm.gather(geometry[:V0.dofmap.index_map.size_local,:], root=root)
global_ct = msh.comm.gather(cell_types[:num_cells_local])
global_vals = msh.comm.gather(u0.x.array[:num_dofs_local])
if(MPI.COMM_WORLD.rank == root):
    try:
        import pyvista
        root_geom = np.vstack(global_geometry)
        root_top = np.concatenate(global_topology)
        root_ct = np.concatenate(global_ct)
        root_vals = np.concatenate(global_vals)
        vals = np.zeros((root_geom.shape[0], 3), dtype=np.float64)
        for p in range(root_geom.shape[0]):
            vals[p][0] = np.real(root_vals[3*p])
            vals[p][1] = np.real(root_vals[3*p+1])
            vals[p][2] = np.real(root_vals[3*p+2])

        print(root_vals, vals)
        grid = pyvista.UnstructuredGrid(root_top, root_ct, root_geom)
        grid["u"] = vals

        pl = pyvista.Plotter(shape=(2, 2))

        pl.subplot(0, 0)
        pl.add_text("magnitude", font_size=12, color="black", position="upper_edge")
        pl.add_mesh(grid.copy(), show_edges=True)

        pl.subplot(0, 1)
        glyphs = grid.glyph(orient="u", factor=0.08)
        pl.add_text("vector glyphs", font_size=12, color="black", position="upper_edge")
        pl.add_mesh(glyphs, show_scalar_bar=False)
        pl.add_mesh(grid.copy(), style="wireframe", line_width=2, color="black")

        pl.subplot(1, 0)
        pl.add_text("x-component", font_size=12, color="black", position="upper_edge")
        pl.add_mesh(grid.copy(), component=0, show_edges=True)
    
        pl.subplot(1, 1)
        pl.add_text("y-component", font_size=12, color="black", position="upper_edge")
        pl.add_mesh(grid.copy(), component=1, show_edges=True)

        pl.view_xy()
        pl.link_views()

    # If pyvista environment variable is set to off-screen (static)
    # plotting save png
        if pyvista.OFF_SCREEN:
            pyvista.start_xvfb(wait=0.1)
            pl.screenshot("uh.png")
        else:
            pl.show()
    except ModuleNotFoundError:
        print("pyvista is required to visualise the solution")

Unfortunately, when I did the DOLFINx install, Adios2 was not installed. I presume I would need to do a reinstall of DOLFINx to make this work. :slightly_frowning_face: