3D cross section of pyvista video/gif

I/O
,
1

I would like to use Pyvista for creating a gif/video of a time-dependent FEM solution on a 3D domain. To see the interior, i would like to create a cross section.

I managed to create

  • a video of the solution without cross-section
  • a static snapshot plot of the solution with cross-section

but not the video with cross section. Here is my MWE code (the code is based on this tutorial):


import matplotlib as mpl
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 = 5
dt = T / num_steps

# mesh and functions
domain = create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [1, 1, 1]], [4,4,4], CellType.hexahedron)
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),
                                        np.isclose(x[2], 0.0),
                                        np.isclose(x[2], 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)

# pyvista
pyvista.start_xvfb()
grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V))
plotter = pyvista.Plotter()
plotter.open_gif("gif.gif", fps=10)
grid.point_data["uh"] = uh.x.array
warped = grid.warp_by_scalar("uh", factor=1)
viridis = mpl.colormaps.get_cmap("viridis").resampled(25)
sargs = dict(title_font_size=25, label_font_size=20, fmt="%.2e", color="black",
             position_x=0.1, position_y=0.8, width=0.8, height=0.1)
renderer = plotter.add_mesh(warped, show_edges=True, lighting=False,
                            cmap=viridis, scalar_bar_args=sargs,
                            clim=[0, max(uh.x.array)])

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

    # 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
    solver.solve(b, uh.vector)
    uh.x.scatter_forward()

    # Update solution
    u_n.x.array[:] = uh.x.array

    # Update gif
    warped = grid.warp_by_scalar("uh", factor=1)
    plotter.update_coordinates(warped.points.copy(), render=False)
    plotter.update_scalars(uh.x.array, render=False)
    plotter.write_frame()
plotter.close()


# plot
cells, types, x = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = uh.x.array
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
clipped = grid.clip("y", invert=True)
plotter.add_mesh(grid, style='wireframe')
plotter.add_mesh(clipped)
maxval=1
plotter.save_graphic("test.pdf")

Related posts, that do not fully answer the question:

Could you help me how to create a video with cross-section in pyvista?

2

If you swap out warped (the grid you use in the .gif making) with clipped (the grid you create at the end of your script which I assume is a clipped cross section), does this produce the solution you seek?

Cheers,
Halvor

3

It sounds correct, but I have not managed to make it work. I dont see how I should edit the code inside the time loop. Do you have a working example?

4

Ok I managed to solve the problem.

As @hherlyng already mentioned, it is of course clear that one has to use clipping with a surface to show the inside of the domain (as e.g. described in the pyvista documentation).

However, it was not obvious to me that one has to add the clipped object only inside the time-loop, and that the methods update_coordinates() as well as update_coordinates() are deprecated. One way to make it work is to call add_mesh() in every time step.

Here is an MWE in case anybody needs it:

import matplotlib as mpl
import pyvista
import ufl
import numpy as np
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, 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 create_box, CellType, locate_entities_boundary

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

# mesh and functions
domain = create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [1, 1, 1]], [4, 4, 4], CellType.hexahedron)
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=2,
                                  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),
                                      np.isclose(x[2], 0.0),
                                      np.isclose(x[2], 1.0))))
dofsV = locate_dofs_topological(V, entity_dim=2, 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)

# pyvista
pyvista.start_xvfb()
grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V))
plotter = pyvista.Plotter()
plotter.open_gif("gif.gif", fps=10)
grid.point_data["uh"] = uh.x.array

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

    # 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
    solver.solve(b, uh.vector)
    uh.x.scatter_forward()

    # Update solution
    u_n.x.array[:] = uh.x.array

    # Update gif with cross-section
    grid.point_data["uh"] = uh.x.array
    grad_clipped = grid.clip("y", invert=True)
    plotter.clear()
    plotter.add_mesh(grid, style='wireframe',clim=[0, 0.1])
    plotter.add_mesh(grad_clipped, clim=[0, 0.1])
    plotter.write_frame()

plotter.close()

# plot
cells, types, x = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = uh.x.array
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
grad_clipped = grid.clip("y", invert=True)
plotter.add_mesh(grid, style='wireframe')
plotter.add_mesh(grad_clipped)
plotter.save_graphic("test.pdf")
1 Like