Hello.

It may sound silly, because the code runs, and it seems like a good

solution, but I want to know if this is the smart way of going from a 2D

discontinuous space to a 3D continuous representation with PyVista. The

following uses interpolation and needs to create an extra functional

space, which seems a bit much at first sight.

Thank you!

## Context

I am solving a problem in a 2D mesh with a vectorial variable

(`u`

in the code below; problem and solution not shown). This

variable with 2 components is then used to calculate a 2×2 tensor

(`strain_fun`

in the code below; mathematical formulation not

shown). To plot it and warp it (with PyVista’s viewer, or

ParaView–after saving it with PyVista’s VTK writer), the data needs to

be in 3D, even if it’s planar. Thus, the conversion of the vector from

1×2 to 1×3 and the tensor from 2×2 to 3×3.

The initial space of `strain_fun`

is discontinuous and the

one of `u`

is continuous by definition. Both of them need to

be plotted on the same `UnstructuredGrid`

. This means that

there are less number of points on the continuous space (elements share

nodes at the boundaries). I thought that the only way to go around it

was with an interpolation, but may be there is a smarter way of doing

it.

## MWE (minimal working example)

```
# Copyright (C) 2023 eDgar (fenicsproject.discourse.group)
#
# This program is free software: you can redistribute it
# and/or modify it under the terms of the GNU General Public
# License as published by the Free Software Foundation,
# version 3 of the License.
#
# This program is distributed in the hope that it will be
# useful, but WITHOUT ANY WARRANTY; without even the implied
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public
# License along with this program. If not, see
# <http://www.gnu.org/licenses/>.
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
from dolfinx import mesh
from dolfinx import plot
import basix
import pyvista as pv
import numpy as np
domain = mesh.create_unit_square(
MPI.COMM_WORLD, 2, 1, mesh.CellType.quadrilateral)
msh_cell_type = domain.basix_cell()
gdim = domain.geometry.dim
# Solution space (solved in gdim dimensions)
deg_u = 2
Ve = basix.ufl.element(
'Lagrange', msh_cell_type, degree=deg_u, shape=(gdim,))
V = fem.functionspace(domain, Ve)
u = fem.Function(V)
# Post-processing space (needs V, still in gdim dimensions)
deg_strain = deg_u
tensor2_shape = (domain.geometry.dim, domain.geometry.dim)
Ve_strain = basix.ufl.element(
'Discontinuous Lagrange', msh_cell_type, degree=deg_strain,
shape=tensor2_shape)
V_strain = fem.functionspace(domain, Ve_strain)
strain_fun = fem.Function(V_strain, name="strain_lagrange")
# Testing with interpolation function (from V_strain)
Ve_test = basix.ufl.element(
'Lagrange', msh_cell_type, degree=deg_strain,
shape=tensor2_shape)
V_test = fem.functionspace(domain, Ve_test)
test_fun = fem.Function(V_test, name="strain_lagrange")
# After solving
test_fun.interpolate(strain_fun)
cells, cell_types, geo_V = plot.vtk_mesh(V_test)
grid3d = pv.UnstructuredGrid(cells, cell_types, geo_V)
u_rows = u.function_space.dofmap.index_map.size_global
u_comp = u.function_space.dofmap.index_map_bs
assert u_comp == gdim
assert u_rows == geo_V.shape[0]
# 2D vector -> 3D vector by adding a column
grid3d["u"] = np.hstack((u.x.array.reshape(u_rows, u_comp),
np.zeros((u_rows, 1))))
# Add row and column
new_col_idx = gdim
new_row_idx = gdim
grid3d["test"] = np.insert(np.insert(
test_fun.x.array.reshape((geo_V.shape[0], gdim, gdim)),
new_row_idx, 0, axis=1), new_col_idx, 0, axis=2)
```

### Debugging code (not needed)

```
num_pts_per_cell = 9
print(f"{grid3d = }")
print(f"{grid3d.n_cells * num_pts_per_cell = }")
print(f"{geo_V.shape = }")
print(f"{V.tabulate_dof_coordinates().shape = }")
print(f"{strain_fun.x.array.shape = }")
print(f"{strain_fun.x.array.shape[0] / V.tabulate_dof_coordinates().shape[0] = }")
print(f"{strain_fun.x.array.shape[0] / V.tabulate_dof_coordinates().shape[0] / gdim = }")
print()
print(f"{test_fun.x.array.shape = }")
print(f"{test_fun.x.array.shape[0] / V.tabulate_dof_coordinates().shape[0] = }")
print(f"{test_fun.x.array.shape[0] / V.tabulate_dof_coordinates().shape[0] / gdim = }")
```

[EDIT: missing output]

```
grid3d = UnstructuredGrid (0x7f7900656d40)
N Cells: 2
N Points: 15
X Bounds: -4.163e-17, 1.000e+00
Y Bounds: -1.110e-16, 1.000e+00
Z Bounds: 0.000e+00, 0.000e+00
N Arrays: 2
grid3d.n_cells * num_pts_per_cell = 18
geo_V.shape = (15, 3)
V.tabulate_dof_coordinates().shape = (15, 3)
strain_fun.x.array.shape = (72,)
strain_fun.x.array.shape[0] / V.tabulate_dof_coordinates().shape[0] = 4.8
strain_fun.x.array.shape[0] / V.tabulate_dof_coordinates().shape[0] / gdim = 2.4
test_fun.x.array.shape = (60,)
test_fun.x.array.shape[0] / V.tabulate_dof_coordinates().shape[0] = 4.0
test_fun.x.array.shape[0] / V.tabulate_dof_coordinates().shape[0] / gdim = 2.0
```

## My system

```
- FEniCSx software
dolfinx: 0.8.0.dev0_r27624.095ef96-1
basix: 0.8.0.dev0_r975.15b915b-1
ufl: 2023.3.0.dev0_r3583.e8d3077-1
ffcx: 0.8.0.dev0_r7098.833967c-1
- Dependencies
python: 3.11.5-2
python-numpy: 1.26.0-1
petsc: 3.19.5.1172.g92f94a14170-1
hdf5-openmpi: 1.14.2-1
boost: 1.83.0-2
adios2: 2.8.3-5
scotch: 7.0.4-1
pybind11: 2.11.1-1
python-build: 1.0.1-1
python-cffi: 1.15.1-4
python-cppimport: 22.08.02.r6.g0849d17-1
python-imageio: 2.33.0-2
python-imageio-ffmpeg: 0.4.9-1
python-installer: 0.7.0-3
python-meshio: 5.3.4-2
python-mpi4py: 3.1.4-3
python-pytest: 7.4.2-1
python-scikit-build: 0.17.6-1
python-setuptools: 1:68.0.0-1
python-wheel: 0.40.0-3
xtensor: 0.24.0-2
xtensor-blas: 0.20.0-1
python-cattrs: 23.1.2-1
python-scikit-build: 0.17.6-1
nanobind: 1.8.0-1
- Operating system
Arch Linux: 6.5.4-arch2-1
```