Pyvista for dolfinx

This package reintroduces simple plot and show functionality to FEniCSx, capable of plotting most common data-structures during both serial and parallel execution.

FEniCSx’ design principles prioritize parallel efficiency, and (in an effort to avoid opaque performance pitfalls) favor fine-grained control over a high-level interface. Consequently, simple, general-purpose plotting routines are no-longer in the standard library. Instead users are deferred to pyvista for their plotting needs. Unfortunately, quite a bit of boiler-plate code is required to interface dolfinx and pyvista. At a prototyping stage one often desires quick-and-dirty visualization with a simple interface. To facilitate this, pyvista4dolfinx provides a single plot function that can be used to plot most of dolfinx visualizable data-structures; scalar- and vector-valued Function, Mesh, FacetMarker, and even integration Measure. The function returns a pyvista.Plotter instance, such that the user still has full access to pyvista’s full range of capabilities.

While this package supports plotting during parallel execution, it should be understood that this comes at a performance penalty; data is gathered and maniputated on a single core. For more data visualization tasks that go beyond basic prototyping, users are encouraged to make use of the vtk export capabilities natively supported in dolfinx.

Demo

The following demo-script illustrates the ease of use. The pyvista4dolfinx-specific lines are highlighted.

import os
import sys

sys.path.insert(0, os.path.abspath("../"))  # Access to pyvista4dolfinx
from pyvista4dolfinx import plot, show, Plotter

from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
import matplotlib.pyplot as plt
from mpi4py import MPI
import numpy as np
import ufl

# Mesh and function space
domain = mesh.create_unit_square(MPI.COMM_WORLD, 20, 20, cell_type=mesh.CellType.quadrilateral)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))

# Quick visualization of mesh
plot(domain, show_partitioning=True, show=True)

# Exogeneous data
u_D = np.array([0, 0], dtype=default_scalar_type)
T = fem.Constant(domain, default_scalar_type((-0.1, 0)))
f = fem.Constant(domain, default_scalar_type((0, -5)))

# Boundary tags
facets_horizontal = mesh.locate_entities_boundary(domain, 
                                              domain.topology.dim - 1, 
                                              lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 1))
facets_bottom = mesh.locate_entities_boundary(domain, 
                                              domain.topology.dim - 1, 
                                              lambda x: np.isclose(x[1], 0))
facets_top = mesh.locate_entities_boundary(domain, 
                                           domain.topology.dim - 1, 
                                           lambda x: np.isclose(x[1], 1))
marked_facets = np.hstack([facets_bottom, facets_top, facets_horizontal])
marked_values = np.hstack([np.full_like(facets_bottom, 1), np.full_like(facets_top, 2), np.full_like(facets_horizontal, 99)])
sorting_map = np.argsort(marked_facets)
facet_tags = mesh.meshtags(domain,domain.topology.dim-1,
                           marked_facets[sorting_map],
                           marked_values[sorting_map])

# Volume tags
def volume_inclusion(x):
    return (x[0] <= 0.75) & (x[0] >= 0.25) & (x[1] <= 0.75) & (x[1] >= 0.25)
elements_inclusion = np.sort(mesh.locate_entities(domain,domain.topology.dim,volume_inclusion))
element_tags = mesh.meshtags(domain,domain.topology.dim,elements_inclusion,3)

# Quick visualization of tags
plot(domain)
plot(element_tags, domain, name="Element tags", cmap="cool", \
                    scalar_bar_args={"vertical": True,"position_x": 0.85,"position_y": 0.25})
plot(facet_tags, domain, name="Facet tags", show=True, 
                    scalar_bar_args={"position_x": 0.21,"position_y": 0.05})

# Measures from tags
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
dx = ufl.Measure("dx", domain=domain, subdomain_data=element_tags)

# Quick visualization of integration measure
plot(dx(3), domain, name="Volume measure", cmap="cool", \
                    scalar_bar_args={"vertical": True,"position_x": 0.85,"position_y": 0.25})
plot(ds(2), domain, name="Surface measure", show=True, 
                    scalar_bar_args={"position_x": 0.21,"position_y": 0.05})

# Clamping condition on the bottom
Vdofs_bottom = fem.locate_dofs_topological(V, domain.topology.dim - 1, facets_bottom)
bc = fem.dirichletbc(u_D, Vdofs_bottom, V)

# Linear elasticity formulation
def epsilon(u):
    return ufl.sym(ufl.grad(u))
def sigma(u):
    return ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * epsilon(u)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * dx
L = ufl.dot(f, v) * dx(3) + ufl.dot(T, v) * ds(2)

# Solve problem
problem = LinearProblem(
    a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
uh = problem.solve()
uh.name = "Displacement"

# Post process; von Mises stress
s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))
V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises, name='Stress')
stresses.interpolate(stress_expr)

# Final visualization
plotter = Plotter(shape=(1, 2))
plot(uh, plotter=plotter, show_mesh=False, subplot=(0,0), factor=0.5)
plot(stresses, warp=uh, show_mesh=False, subplot=(0,1))
show()

When executed on three cores, the following plot routines output the below graphics.

plot(domain, show_partitioning=True, show=True)

plot(domain)
plot(element_tags, domain, name="Element tags", cmap="cool", \
     scalar_bar_args={"vertical": True,"position_x": 0.85,"position_y": 0.25})
plot(facet_tags, domain, name="Facet tags", show=True,
     scalar_bar_args={"position_x": 0.21,"position_y": 0.05})

plot(dx(3), domain, name="Volume measure", cmap="cool", \
     scalar_bar_args={"vertical": True,"position_x": 0.85,"position_y": 0.25})
plot(ds(2), domain, name="Surface measure", show=True,
     scalar_bar_args={"position_x": 0.21,"position_y": 0.05})

plotter = Plotter(shape=(1, 2))
plot(uh, plotter=plotter, show_mesh=False, subplot=(0,0), factor=0.5)
plot(stresses, warp=uh, show_mesh=False, subplot=(0,1))
show()

Documentation

pyvista4dolfinx