Explicit name for components of tensors exported in VTK

I/O
1

Hi,

I’m solving poromechanics problems and therefore I export tensor fields (eg. strain or Cauchy’s stress) in VTK files for postprocessing in Paraview.

I succeeded into projecting the tensor field, exporting them, and opening them in Paraview. However, the issue I have is that the components of the tensor are labeled only by an integer (between 0 and 8) in Paraview, and I don’t know the corresponding component (eg is Strain 1 the yy component of the strain or the xy component?).

Is there a way to have explicit name (like EXX, EXY or EYY) for these components in the VTK files?

I’m using dolfinx v0.9.0 installed on linux with conda.

Here is a working example (but any code exporting tensor field may have the same behavior) :

import matplotlib.pyplot as plt
import numpy as np
import gmsh
from mpi4py import MPI
import basix
import ufl
from dolfinx import fem, mesh, io, nls
from dolfinx.io.gmshio import model_to_mesh
import os
import pathlib

import dolfinx.fem.petsc
import dolfinx.nls.petsc

# Change working directory (allow to run this file from any where)
_dir_path = os.path.dirname(os.path.realpath(__file__))  # dir_path of this file
os.chdir(_dir_path)

res_folder = pathlib.Path("./results")
res_folder.mkdir(parents=True, exist_ok=True)


# Create mesh
L = 1.0
R = 0.1
N = 25

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)  # to disable meshing info
gdim = 2
model_rank = 0
gmsh.model.add("Model")

gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, L, L, tag=1)
gmsh.model.occ.addDisk(0.0, 0.0, 0.0, R, R, tag=2)
gmsh.model.occ.cut([(2, 1)], [(2, 2)])

gmsh.model.occ.synchronize()

gmsh.option.setNumber("Mesh.CharacteristicLengthMin", L / N)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", L / N)

volumes = gmsh.model.getEntities(gdim)
assert len(volumes) == 1
gmsh.model.addPhysicalGroup(
    gdim, [volumes[i][1] for i in range(len(volumes))], 1, name="Volume"
)

gmsh.model.addPhysicalGroup(gdim - 1, [1], 1, name="inner_circle")
gmsh.model.addPhysicalGroup(gdim - 1, [5], 2, name="bottom")
gmsh.model.addPhysicalGroup(gdim - 1, [2], 3, name="left")
gmsh.model.addPhysicalGroup(gdim - 1, [3, 4], 4, name="right_top")
gmsh.model.mesh.generate(gdim)

domain, _, facets = model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim=gdim)

gmsh.finalize()

# Define elements spaces
Vue = basix.ufl.element(
    "P", domain.basix_cell(), 2, shape=(2,)
)  # displacement finite element
Vpe = basix.ufl.element("P", domain.basix_cell(), 1)  # Pressure finite element
V = fem.functionspace(domain, basix.ufl.mixed_element([Vue, Vpe]))

V_ux, _ = V.sub(0).sub(0).collapse()  # used for Dirichlet BC
V_uy, _ = V.sub(0).sub(1).collapse()  # used for Dirichlet BC
V_p, _ = V.sub(1).collapse()  # used for Dirichlet BC


# --------------------------------------------------
# Define problem constants
# --------------------------------------------------
Pref_value = fem.Constant(domain, 0.0)
DPhole_value = 100e3
dt = fem.Constant(domain, 0.0)  # time step

DPhole = fem.Function(V_p)
DPhole.x.array[:] = DPhole_value

Pref = fem.Function(V_p)
Pref.x.array[:] = Pref_value

# --------------------------------------------------
# Dirichlet BC (u, p)
# --------------------------------------------------
inner_P_dofs = fem.locate_dofs_topological((V.sub(1), V_p), facets.dim, facets.find(1))
bottom_uy_dofs = fem.locate_dofs_topological(
    (V.sub(0).sub(1), V_uy), facets.dim, facets.find(2)
)
left_ux_dofs = fem.locate_dofs_topological(
    (V.sub(0).sub(0), V_ux), facets.dim, facets.find(3)
)
# used for post-processing
bottom_P_dofs = fem.locate_dofs_topological(
    (V.sub(1), V_p), facets.dim, facets.find(2)
)[1]


u0x = fem.Function(V_ux)
u0y = fem.Function(V_uy)
bcs = [
    fem.dirichletbc(DPhole, inner_P_dofs, V.sub(1)),
    fem.dirichletbc(u0y, bottom_uy_dofs, V.sub(0).sub(1)),
    fem.dirichletbc(u0x, left_ux_dofs, V.sub(0).sub(0)),
]

# --------------------------------------------------
# Variational problem definition
# --------------------------------------------------
v = fem.Function(V)
(u, p) = ufl.split(v)
v_ = ufl.TestFunction(V)
(u_, p_) = ufl.split(v_)
dv = ufl.TrialFunction(V)

V_aux = fem.functionspace(domain, ("DG", 1))
phi_old = fem.Function(V_aux, name="Previous_porosity")


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


K = fem.Constant(domain, 100e6)  # [Pa]
mu = fem.Constant(domain, 60e6)  # [Pa]
b = fem.Constant(domain, 1.0)  # [1]
N = fem.Constant(domain, 1e7)  # [Pa]
rho = fem.Constant(domain, 1e3)  # [kg/m^3]
# perm = fem.Constant(domain, 1e-9 / 1e-3)  # [m^3.s/kg] cf. Coussy p.47
perm = fem.Constant(domain, 1e-9)  # [m^3.s/kg] cf. Coussy p.47
grav = fem.Constant(domain, [0.0, 0.0])  # [m/s^2]

lam = K - 2.0 * mu / 3.0  # [Pa]

sig = (
    lam * ufl.tr(eps(u)) * ufl.Identity(gdim)
    + 2 * mu * eps(u)
    - p * b * ufl.Identity(gdim)
)

phi = b * ufl.tr(eps(u)) + p / N
phi_expr = fem.Expression(phi, V_aux.element.interpolation_points())

mech_res = (ufl.inner(sig, eps(u_)) + rho * ufl.dot(grav, u_)) * ufl.dx
poro_res = (
    (phi - phi_old) / dt * p_ + perm * ufl.dot(ufl.grad(p), ufl.grad(p_))
) * ufl.dx

Res = mech_res + poro_res
Jac = ufl.derivative(Res, v, dv)

problem = fem.petsc.NonlinearProblem(Res, v, bcs=bcs, J=Jac)

newton = nls.petsc.NewtonSolver(domain.comm, problem)
newton.rtol = 1e-8
newton.atol = 1e-8
newton.convergence_criterion = "incremental"
newton.report = True
newton.max_it = 10
ksp = newton.krylov_solver

# --------------------------------------------------
# Resolution
# --------------------------------------------------
Nincr = 50
t = np.logspace(1, 2, Nincr + 1)
x = V_p.tabulate_dof_coordinates()[bottom_P_dofs, 0]  # x position of dofs

vtk_d = io.VTKFile(domain.comm, res_folder / "displacement.pvd", "w")

for i, dti in enumerate(np.diff(t)):
    dt.value = dti

    num_its, converged = newton.solve(v)
    assert converged

    phi_old.interpolate(phi_expr)

    u_out = v.sub(0).collapse()
    u_out.name = "Displacement"

    V_strain_out = fem.functionspace(domain, ("Lagrange", 1, (2, 2)))
    strain_expr = fem.Expression(
        ufl.sym(ufl.grad(u_out)), V_strain_out.element.interpolation_points()
    )
    strain_out = fem.Function(V_strain_out)
    strain_out.name = "Strain"
    strain_out.interpolate(strain_expr)

    vtk_d.write_function(u_out, i)
    vtk_d.write_function(strain_out, i)


vtk_d.close()

Thanks in advance for help.

2

Paraview internally uses the following convention for tensors:
[[0, 1, 2],
[3, 4, 5],
[6, 7, 8]]
it does not care about symmetry whatsoever

2 Likes