Mixed formulation: Solution data length inconsistent with the mesh -> impossible to visualise

dolfinx
1

I’m solving a fourth order problem with a mixed formulation approach, and although the computation runs without crashing, I am not able to visualise the result due to inconsistent shape of the solution with the mesh (ValueError: data length of (195770) != required length (19729)).

The only thing I know is that when I solve a 2nd order problem without a mixed formulation but with similar code, it works.

Formulation

I need to solve the equation \left(\Delta^2 - \Delta + \kappa^2\right)u=0 on a disk domain \Omega with three circular holes which I call “channels”. The boundaries are noted \Gamma_i where i=0 on the external boundary and i=1,2,3 on the channels.

I use a two-fields formulation:
\begin{cases} \sigma - \Delta u = 0 \text{ on } \Omega\\ \Delta \sigma - \sigma + \kappa^2u = 0 \text{ on } \Omega\\ u = u_i \text{ on } \Gamma_i\\ \sigma = 0 \text{ on } \Gamma_i\\ \end{cases}

The variational formulation is
\int_\Omega \left( -\nabla\sigma\nabla v - \sigma v + \kappa^2uv + \sigma\tau + \nabla u\nabla\tau \right)\mathrm{d}\Omega = 0
with (\tau,v) being the test functions.

Implementation

import numpy as np

import gmsh
from basix.ufl import element, mixed_element
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
from dolfinx import fem, plot, default_scalar_type
from dolfinx.io import gmshio
from dolfinx.mesh import locate_entities, meshtags

import pyvista
from ufl import TrialFunctions, TestFunctions


domain_radius=30.0
channel_radius=0.500
channels_locations=[(-1.8, 0), (1.8, 0), (0, 1)]
boundary_conditions=[-0.5, -0.5, 1, ]
kappa=1.0

# Mesh #

gmsh.initialize()
gmsh.option.setNumber('General.Verbosity', 1)

## Domain
ext_boundary = gmsh.model.occ.addCircle(0, 0, 0, domain_radius)
ext_boundary_loop = gmsh.model.occ.addCurveLoop([ext_boundary])

## Channels
channels_loops = []
for (x, y) in channels_locations:
    channel = gmsh.model.occ.addCircle(x, y, 0, channel_radius)
    channels_loops.append(gmsh.model.occ.addCurveLoop([channel]))

membrane = gmsh.model.occ.addPlaneSurface([ext_boundary_loop] + channels_loops)
gmsh.model.occ.synchronize()

gdim = 2
gmsh.model.addPhysicalGroup(gdim, [membrane], 1)

## Mesh generation
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", channel_radius/32)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", domain_radius/8)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 128)

gmsh.model.mesh.generate(gdim)

gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
mesh_model = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
gmsh.finalize()

domain = mesh_model[0]

# Boundaries #

def external_boundary_locator(x) -> bool:
    return np.isclose(x[0]**2 + x[1]**2, domain_radius**2)

def channel_i_boundary_locator(x, i) -> bool:
    return np.isclose((x[0] - channels_locations[i][0])**2 + (x[1] - channels_locations[i][1])**2, channel_radius**2)

def boundary_locator(x) -> bool:
    return np.logical_or.reduce([external_boundary_locator(x)] + [channel_i_boundary_locator(x, i) for i in range(len(channels_locations))])

boundaries = list(enumerate([external_boundary_locator] + [lambda x, _i=i: channel_i_boundary_locator(x, _i) for i in range(len(channels_locations))]))

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1

for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))

facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# Solution #

k = 3
mixed_el = mixed_element([
element("Lagrange", domain.basix_cell(), k),
element("Lagrange", domain.basix_cell(), k - 2),
])
V = fem.functionspace(domain, mixed_el)
V_sigma, V_u = V.sub(0), V.sub(1)
V_sigma0, _ = V_sigma.collapse()
V_u0, _ = V_u.collapse()

## Dirichlet conditions

zero_constant = fem.Constant(domain, default_scalar_type(0))

### External boundary
_, locator_0 = boundaries[0]

dofs_dirichlet_sigma = fem.locate_dofs_geometrical((V_sigma, V_sigma0), boundary_locator)
dofs_dirichlet_u_0 = fem.locate_dofs_geometrical((V_u, V_u0), locator_0)

bcs_dirichlet = [
fem.dirichletbc(zero_constant, dofs_dirichlet_sigma[0], V_sigma),
fem.dirichletbc(zero_constant, dofs_dirichlet_u_0[0], V_u),
]

### Channels
for i, (marker, locator) in enumerate(boundaries[1:]):
    dofs_dirichlet_u = fem.locate_dofs_geometrical((V_u, V_u0), locator)
    value_dirichlet_u = fem.Constant(domain, default_scalar_type(boundary_conditions[i]))
    bc_u = fem.dirichletbc(value_dirichlet_u, dofs_dirichlet_u[0], V_u)
    bcs_dirichlet.append(bc_u)


## Variational formulation

(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)

a = (-ufl.dot(ufl.grad(sigma), ufl.grad(tau)) - sigma * tau + u * tau + sigma * v + ufl.dot(ufl.grad(u), ufl.grad(v))) * ufl.dx
L = zero_constant * v * ufl.dx

problem = LinearProblem(a, L, bcs=bcs_dirichlet, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
wh = problem.solve()
_, uh = wh.split()


# Visualisation #

mesh = uh.function_space.mesh
cells, types, x = plot.vtk_mesh(mesh)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = uh.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
if pyvista.OFF_SCREEN:
    pyvista.start_xvfb(wait=0.1)
    plotter.screenshot("figures/uh_membrane.png")
else:
    plotter.show()

Traceback

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 136
    134 cells, types, x = plot.vtk_mesh(mesh)
    135 grid = pyvista.UnstructuredGrid(cells, types, x)
--> 136 grid.point_data["u"] = uh.x.array.real
    137 grid.set_active_scalars("u")
    138 plotter = pyvista.Plotter()

File ~/miniforge3/envs/fenicsx-env/lib/python3.12/site-packages/pyvista/core/datasetattributes.py:227, in DataSetAttributes.__setitem__(self, key, value)
    224     raise TypeError('Only strings are valid keys for DataSetAttributes.')
    226 has_arr = key in self
--> 227 self.set_array(value, name=key)
    229 # do not make array active if it already exists.  This covers
    230 # an inplace update like self.point_data[key] += 1
    231 if has_arr:

File ~/miniforge3/envs/fenicsx-env/lib/python3.12/site-packages/pyvista/core/datasetattributes.py:579, in DataSetAttributes.set_array(self, data, name, deep_copy)
    576 if not isinstance(name, str):
    577     raise TypeError('`name` must be a string')
--> 579 vtk_arr = self._prepare_array(data, name, deep_copy)
    580 self.VTKObject.AddArray(vtk_arr)
    581 self.VTKObject.Modified()

File ~/miniforge3/envs/fenicsx-env/lib/python3.12/site-packages/pyvista/core/datasetattributes.py:737, in DataSetAttributes._prepare_array(self, data, name, deep_copy)
    735     data = tmparray
    736 if data.shape[0] != array_len:
--> 737     raise ValueError(f'data length of ({data.shape[0]}) != required length ({array_len})')
    739 # attempt to reuse the existing pointer to underlying VTK data
    740 if isinstance(data, pyvista_ndarray):
    741     # pyvista_ndarray already contains the reference to the vtk object
    742     # pyvista needs to use the copy of this object rather than wrapping
    743     # the array (which leaves a C++ pointer uncollected.

ValueError: data length of (195770) != required length (19729)

Code that works with another equation:

import numpy as np

import gmsh
from basix.ufl import element, mixed_element
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
from dolfinx import fem, plot, default_scalar_type
from dolfinx.io import gmshio
from dolfinx.mesh import locate_entities, meshtags

import pyvista
from ufl import TrialFunctions, TestFunctions


domain_radius=30.0
channel_radius=0.500
channels_locations=[(-1.8, 0), (1.8, 0), (0, 1)]
boundary_conditions=[-0.5, -0.5, 1, ]

# Mesh #

gmsh.initialize()
gmsh.option.setNumber('General.Verbosity', 1)

## Domain
ext_boundary = gmsh.model.occ.addCircle(0, 0, 0, domain_radius)
ext_boundary_loop = gmsh.model.occ.addCurveLoop([ext_boundary])

## Channels
channels_loops = []
for (x, y) in channels_locations:
    channel = gmsh.model.occ.addCircle(x, y, 0, channel_radius)
    channels_loops.append(gmsh.model.occ.addCurveLoop([channel]))

membrane = gmsh.model.occ.addPlaneSurface([ext_boundary_loop] + channels_loops)
gmsh.model.occ.synchronize()

gdim = 2
gmsh.model.addPhysicalGroup(gdim, [membrane], 1)

## Mesh generation
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", channel_radius/32)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", domain_radius/8)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 128)

gmsh.model.mesh.generate(gdim)

gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
mesh_model = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
gmsh.finalize()

domain = mesh_model[0]

# Boundaries #

def external_boundary_locator(x) -> bool:
    return np.isclose(x[0]**2 + x[1]**2, domain_radius**2)

def channel_i_boundary_locator(x, i) -> bool:
    return np.isclose((x[0] - channels_locations[i][0])**2 + (x[1] - channels_locations[i][1])**2, channel_radius**2)

def boundary_locator(x) -> bool:
    return np.logical_or.reduce([external_boundary_locator(x)] + [channel_i_boundary_locator(x, i) for i in range(len(channels_locations))])

boundaries = list(enumerate([external_boundary_locator] + [lambda x, _i=i: channel_i_boundary_locator(x, _i) for i in range(len(channels_locations))]))

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1

for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))

facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# Solution #

## Function space
V = fem.functionspace(domain, ("Lagrange", 1, ))

## Dirichlet conditions

bcs_dirichlet = []
for i, (marker, locator) in enumerate(boundaries[1:]):
    dofs_dirichlet = fem.locate_dofs_geometrical(V, locator)
    value_dirichlet = fem.Constant(domain, default_scalar_type(boundary_conditions[i]))
    bc = fem.dirichletbc(value_dirichlet, dofs_dirichlet, V)
    bcs_dirichlet.append(bc)

## Variational formulation

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.dot(u, v) * ds(0) + ufl.dot(u, v) * ufl.dx
L = fem.Constant(v.ufl_function_space().mesh, default_scalar_type(0)) * v * ufl.dx

problem = LinearProblem(a, L, bcs=bcs_dirichlet, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()


# Visualisation #

mesh = uh.function_space.mesh
cells, types, x = plot.vtk_mesh(mesh)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = uh.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
if pyvista.OFF_SCREEN:
    pyvista.start_xvfb(wait=0.1)
    plotter.screenshot("figures/uh_membrane.png")
else:
    plotter.show()
2

You need to do a collapse here, as split is not sufficient (as uh resides in the mixed space, with the array being that of the mixed function space: dolfinx.fem — DOLFINx 0.9.0 documentation)
i.e. call uh = wh.split()[0].collapse() or uh = wh.sub(0).collapse()