Shape-mismatch between mixed-space mesh and the trial function

Hi Folks,

I’m new to FEniCSx and I’ve been trying learn to implement a coupled system of time-dependent equations using a mixed function space – in this case a 3D vector space which will represent a displacement of a mirror undergoing thermoelastic deformation, and a scalar finite element space which will represent the temperature of mirror.

I think I understand how to write a solver for my problem, and I have been able to run a solver end-to-end without issue. The problem I am facing arises when I try to make a movie, or plot a frame of one of the solutions following an example from Jørgen Dokken’s FEniCSx tutorial here using pyvista. Ultimately, I run into a shape mismatch between the size of the 3D mesh on one of the subspaces and the size of the solution array. Below is a MWE which just highlights the shape mismatch (I’ve removed the solver/solution because these are not necessary to reproduce the issue).

``````from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI
from dolfinx.io import gmshio, XDMFFile
import dolfinx
from dolfinx import mesh, fem, plot, io, nls
import numpy as np
import petsc4py.PETSc as PETSc
from petsc4py.PETSc import ScalarType
import ufl
import matplotlib.pyplot as plt
import pyvista as pv
import os

import gmsh
gmsh.initialize()

R = 0.34/2      # radius [m]
H = 0.2         # height [m]

mirror = gmsh.model.occ.addCylinder(0, 0, 0 ,0, 0, H, R)
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
surfaces = gmsh.model.occ.getEntities(dim=2)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",0.02)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",0.02)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(3)

model_rank = 0
mesh, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank)

# Define the function spaces

v_uel = ufl.VectorElement('Lagrange', mesh.ufl_cell(), 2)
v_tel = ufl.FiniteElement('Lagrange', mesh.ufl_cell(), 2)
mel = ufl.MixedElement([v_uel, v_tel])
Vue = fem.FunctionSpace(mesh, v_uel)
Vte = fem.FunctionSpace(mesh, v_tel)

#the mixed element space
Y = fem.FunctionSpace(mesh, v_uel* v_tel)

#trial functions
U = fem.Function(Y)  # current solution

# pyvista stuff
os.environ["PATH"] += os.pathsep + '/opt/X11/bin' # THIS is for xvfb (requires xQuartz on Mac)
pv.start_xvfb()

# Create plotter and pyvista grid
p = pv.Plotter()

#extract the subspace for the vector solution V0
V0, dofs = Y.sub(0).collapse()

topology, cell_types, geometry = plot.create_vtk_mesh(V0)

grid = pv.UnstructuredGrid(topology, cell_types, geometry)

# Attach vector values to grid and warp grid by vector
grid["u"] = U.sub(0).x.array.reshape((geometry.shape[0], 3))
``````

This fails with the message

``````---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[6], line 14
11 grid = pv.UnstructuredGrid(topology, cell_types, geometry)
13 # Attach vector values to grid and warp grid by vector
---> 14 grid["u"] = U.sub(0).x.array.reshape((geometry.shape[0], 3))

ValueError: cannot reshape array of size 68744 into shape (17186,3)
``````

The code versions are

dolfinx: ‘0.6.0’
ufl: ‘2023.1.1.post0’
pyvista: ‘0.38.2’
gmsh: ‘4.11.1’
OS: MacOS Monterey 12.2.1

Thus, you would have to call `U.sub(0).collapse().x.array.reshape(...)` to get the contiguous array you want.