# How to go from a 2D mesh with [2nd-degree, 2×2 (shape) tensorial, discontinuous Lagrange] elements to 3×3 PyVista UnstructuredGrid with continuous elements with vectors and tensors?

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
# 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
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(
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))))
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
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
``````

Instead of moving the discontinuous function to a continuous space, i would interpolate the continuous function into the discontinuous space.
This is because the continuous space is a sub space of the discontinuous space.

Regarding going from 2D to 3D, see Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial

Thank you very much, @dokken. Do you mean like this? (the diff is below
in case that is faster).

``````# 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
# 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
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
from dolfinx import mesh
from dolfinx import plot
from dolfinx import default_scalar_type as num_type
import basix
import pyvista as pv
import numpy as np

domain = mesh.create_unit_square(
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)
Ve_test = basix.ufl.element(
'Discontinuous Lagrange', msh_cell_type, degree=deg_u,
shape=(gdim,))
V_test = fem.functionspace(domain, Ve_test)
test_fun = fem.Function(V_test, name="test")

# Fake value
u.x.array[:5] = 1

# After solving
test_fun.interpolate(u)
assert np.any(np.isclose(test_fun.x.array, 1))

cells, cell_types, geo_V = plot.vtk_mesh(V_test)
grid3d = pv.UnstructuredGrid(cells, cell_types, geo_V)
assert np.array_equal(geo_V, V_test.tabulate_dof_coordinates())

u_rows = test_fun.function_space.dofmap.index_map.size_global
u_comp = test_fun.function_space.dofmap.index_map_bs
assert u_comp == gdim
assert u_rows == geo_V.shape[0]
# 2D vector -> 3D vector by partially filling up 1 × 3 array
values = np.zeros((u_rows, 3), dtype=num_type)
values[:, :-1] = test_fun.x.array.reshape(u_rows, u_comp)
grid3d["u"] = values
# Fill only gdim × gdim
values = np.zeros((u_rows, 3, 3), dtype=num_type)
values[:, :-1, :-1] = strain_fun.x.array.reshape(
u_rows, *tensor2_shape)
grid3d["test"] = values
``````

[EDIT: (above) from V_test → from V; ++ assert np.array_equal(… ]

May be this is easier to check the difference:

``````Copyright (C) 2023 eDgar (fenicsproject.discourse.group)
Last changed: 2023-12-04
--- /tmp/ediffv7uYJl  2023-12-04
+++ /tmp/ediffAsQf5L  2023-12-04
@@ -19,6 +19,7 @@
from dolfinx import fem
from dolfinx import mesh
from dolfinx import plot
+from dolfinx import default_scalar_type as num_type
import basix
import pyvista as pv
import numpy as np
@@ -44,29 +45,33 @@
V_strain = fem.functionspace(domain, Ve_strain)
strain_fun = fem.Function(V_strain, name="strain_lagrange")

-# Testing with interpolation function (from V_strain)
+# Testing with interpolation function (from V_test)
Ve_test = basix.ufl.element(
-    'Lagrange', msh_cell_type, degree=deg_strain,
-    shape=tensor2_shape)
+    'Discontinuous Lagrange', msh_cell_type, degree=deg_u,
+    shape=(gdim,))
V_test = fem.functionspace(domain, Ve_test)
-test_fun = fem.Function(V_test, name="strain_lagrange")
+test_fun = fem.Function(V_test, name="test")
+
+# Fake value
+u.x.array[:5] = 1

# After solving
-test_fun.interpolate(strain_fun)
+test_fun.interpolate(u)
+assert np.any(np.isclose(test_fun.x.array, 1))

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
+u_rows = test_fun.function_space.dofmap.index_map.size_global
+u_comp = test_fun.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))))
-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)
+# 2D vector -> 3D vector by partially filling up 1 × 3 array
+values = np.zeros((u_rows, 3), dtype=num_type)
+values[:, :-1] = test_fun.x.array.reshape(u_rows, u_comp)
+grid3d["u"] = values
+# Fill only gdim × gdim
+values = np.zeros((u_rows, 3, 3), dtype=num_type)
+values[:, :-1, :-1] = strain_fun.x.array.reshape(
+    u_rows, *tensor2_shape)
+grid3d["test"] = values
``````