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