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:

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

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

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!