Plotting quantities on-the-fly with VTXWriter as is done with XDMF

dolfinx
1

I am trying to plot several quantities in a domain using Paraview. This works using the XDMFFile writer, but all quantities are stored with separate coordinates, even though they are on the same domain. This removes some functionality from Paraview, such as first warping a space by displacements, and then viewing the strain on that warped domain.
This problem was addressed in this forum question, where the solution is to switch to VTXWriter.

I want to do this, but don’t know how to make it work for my setup. At the moment of creating my file_writer, I don’t know yet which quantities will be plotted. Instead, I pass the writer around to a bunch of classes, which decide what to plot. I have provided a MWE example using XDMFFile below. How do I make this work with VTXWriter?

import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import fem, mesh, nls, io
import petsc4py.PETSc as PETSc
import dolfinx.fem.petsc
import dolfinx.nls.petsc
from petsc4py.PETSc import ScalarType
from dolfinx import nls

BC_type = "displacement"

L = 5.0
domain = mesh.create_rectangle(
    MPI.COMM_WORLD,
    [(0.0, -1 / 2), (L, 1 / 2)],
    (20, 5),
    diagonal=mesh.DiagonalType.crossed,
)

E = fem.Constant(domain, 1e5)
nu = fem.Constant(domain, 0.3)
mu = E / 2 / (1 + nu)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)


def eps(v):
    return ufl.sym(ufl.grad(v))


def sigma(v):
    return lmbda * ufl.tr(eps(v)) * ufl.Identity(2) + 2.0 * mu * eps(v)


V = fem.functionspace(domain, ("P", 1, (2,)))
u = fem.Function(V, name="Displacement")
u_ = ufl.TestFunction(V)

# Dirichlet BC
def left(x):
    return np.isclose(x[0], 0.0)

def right(x):
    return np.isclose(x[0], L)

fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

left_dofs = fem.locate_dofs_geometrical(V, left)
u_bc = fem.Function(V)
u_bc.x.array[:] = 0.0

# Right boundary
right_dofs = fem.locate_dofs_topological(V.sub(0), fdim, right_facets)
u_right = fem.Constant(domain, PETSc.ScalarType(1.))
bcs = [fem.dirichletbc(u_bc, left_dofs), fem.dirichletbc(u_right, right_dofs, V.sub(0))]
U_x = np.linspace(0, 0.2 , 11)[1:]

# Problem
residual = ufl.inner(sigma(u), eps(u_)) * ufl.dx
problem = fem.petsc.NonlinearProblem(residual, u, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)

class arbitrary_class:
    def __init__(self, domain):
        self.V_scalar = fem.functionspace(domain, ("P", 1))
        self.strain_xx = fem.Function(self.V_scalar, name="eps_xx")

    def plot(self, plotter, strain, i):
        self.strain_xx.interpolate(fem.Expression(strain[0, 0], self.V_scalar.element.interpolation_points()))
        plotter.write_function(self.strain_xx, i)

other_class = arbitrary_class(domain)

# Working example with XDMFFile
file_results = io.XDMFFile(
    domain.comm,
    f"plot_MWE_XDMF.xdmf",
    "w",
)

# # How to do this with VTXWriter or some other writer?
# file_results = io.VTXWriter(
#     domain.comm,
#     f"plot_MWE_VTK.bp",
#     ??????
# )

file_results.write_mesh(domain)

for i in range(len(U_x)):
    u_right.value = U_x[i]
    its, converged = solver.solve(u)
    assert converged
    u.x.scatter_forward()

    # Plot displacement
    file_results.write_function(u, i)

    # Plot some other quantity
    other_class.plot(file_results, sigma(u), i)
1 Like
2

That simply won’t work with VTXWriter.
This is due to how the VTX format is specifed at: Visualizing Data — ADIOS2 2.10.1 documentation

I guess you might be able to do it with the VTKFile.