Export field from triangular mesh to numpy

mesh
1

Hello,

I’m working with Mixed Poisson implementation and I want to export the pressure field solution to numpy and use it as a benchmark for other methods.

The function space is based on a triangular mesh with BDM*DG structure and for a 64x64 mesh, if I export the solution with

u_topology, u_cell_types, u_geometry = plot.create_vtk_mesh(BDM_out)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)

u_grid.cell_data["u"] = u_out.x.array

the size of the vector u_out.x.array is actually (8192,), which is 64*64*2.

The actual solution looks like this:

dpm_fem
dpm_fem744×657 197 KB

But if I try to plot it with numpy and matplotlib, apart from doubling one axes’ dimension, I get:

sol = u_out.x.array
np_u = sol.reshape(mesh_size, mesh_size*2)
fig, ax = plt.subplots()

pos = ax.imshow(u, interpolation = "bilinear")
fig.colorbar(pos, ax = ax)
fig.show()

Screenshot from 2023-03-08 17-18-21

Now, I’m aware of the counting method for Basix triangles
image

But still I’m not able to replicate the first plot.

How should I refactor the solution in order to visualize it correctly outside pyvista?

Thanks!

2

First of all, you should provide all definitions and code snippts to be able to reproduce your issue.
Note that matplotlib can only plot First order triangles (CG 1 functions) or the nodal values of your function.
This is the reason why we use pyvista in DOLFINx, as it can plot unstructured grids in 2D and 3D.

You need to relate the degrees of freedom (entries of u_out.x.array) with a physical coordinate, i.e.
BDM_out.tabulate_dof_coordinates() as shown below:

import matplotlib.pyplot as plt
from mpi4py import MPI
import dolfinx

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
BDM_out = dolfinx.fem.FunctionSpace(mesh, ("DG", 2))
u_out = dolfinx.fem.Function(BDM_out)
u_out.interpolate(lambda x: x[0] + 2 * x[1])

X = BDM_out.tabulate_dof_coordinates()


plt.scatter(X[:, 0], X[:, 1], c=u_out.x.array, cmap="viridis")
plt.savefig("figure.png")

I s there a specific reason for wanting to use matplotlib to plot your data?

Note that

are numpy arrays, dolfinx.plot — DOLFINx 0.6.0 documentation
which are input to pyvista:
Redirecting to https://docs.pyvista.org/version/stable/api/core/_autosummary/pyvista.UnstructuredGrid.html
See for instance their example:

cells = [4, 0, 1, 2, 3]
celltypes = [pyvista.CellType.TETRA]
points = [
    [1.0, 1.0, 1.0],
    [1.0, -1.0, -1.0],
    [-1.0, 1.0, -1.0],
    [-1.0, -1.0, 1.0],
]
grid = pyvista.UnstructuredGrid(cells, celltypes, points)
grid.plot(show_edges=True)
3

Hey there, this wasn’t my question but this answer really helped me. And I do have an answer to your question, though it may be different than nicolot’s.

While PyVista has some really cool features, generally, I am already quite proficient with matplotlib, and there is an incredible amount of documentation for it. I think Fenics’ use of PyVista is actually really important and nice, but often, I just want to pump out a simple 2D image of some simulation result, and that is extremely straightforward with matplotlib. Still not sure how to do that with PyVista!