Dolfinx - Projection of a trial function into its tensor space in a mixed formulation

variational formulation
1

Hello everyone,

I am working with a stabilized mixed formulation with displacement and strain fields in FEniCSx, and I would like to project a trial function into its corresponding tensor space. (At the end of this post I’m attaching a minimal test code.)

The issue is that projection only seems to work if I introduce an auxiliary tensor field, but this changes the nature of my trial function (it becomes just a Function instead of remaining a TrialFunction).

On the other hand, if I try to directly obtain the tensor subspace from the mixed problem, it does not behave as expected — the tensor space is automatically converted into a vector space, and I lose the original tensorial structure.
Concretely, after collapsing W.sub(1) (the symmetric tensor component), the resulting space has shape (3,) in 2D (Voigt-like representation), so TestFunction there is a vector, while my strain trial has shape (2,2) — leading to a shape mismatch in ufl.inner.

It is worth mentioning that in FEniCS Legacy, if the subspace of the mixed space was collapsed, this procedure did work as I intended.

Has anyone faced this issue in FEniCSx? Is there a way to project a trial function while keeping its tensorial nature, without having to introduce an auxiliary field? Alternatively, is there a recommended tensor↔Voigt mapping workflow when working with collapsed subspaces?

Any suggestions or experiences would be very helpful!

Thanks in advance.

from mpi4py import MPI
import ufl
import dolfinx
import basix.ufl
import numpy as np
import dolfinx.fem.petsc

def project(v, V):
    """Proyección L2 de una expresión UFL `v` al espacio `V`."""
    bcs = []
    u = ufl.TrialFunction(V)
    w = ufl.TestFunction(V)
    a = ufl.inner(u, w)*ufl.dx
    L = ufl.inner(v, w)*ufl.dx
    problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=bcs,
                            petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    return problem.solve()

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

elem_e = basix.ufl.element("Lagrange", mesh.basix_cell(), 1, shape=(2,2), symmetry=True)
elem_u = basix.ufl.element("Lagrange", mesh.basix_cell(), 1, shape=(2,))
el_mixed = basix.ufl.mixed_element([elem_u, elem_e])

W = dolfinx.fem.functionspace(mesh, el_mixed)
u, eps = ufl.TrialFunctions(W)
v, eta = ufl.TestFunctions(W)

print(ufl.shape(u))
print(ufl.shape(eps))
print(ufl.shape(v))
print(ufl.shape(eta))

wh = dolfinx.fem.Function(W)
uh, epsh = ufl.split(wh)

print(ufl.shape(uh))
print(ufl.shape(epsh))

V0, _ = W.sub(1).collapse()

trial = dolfinx.fem.Function(V0)
print(ufl.shape(trial))

V1 = dolfinx.fem.functionspace(mesh, elem_e)

trial2 = dolfinx.fem.Function(V1)
print(ufl.shape(trial2))

This is what the prints gave:

(2,)
(2, 2)
(2,)
(2, 2)
(2,)
(2, 2)
(3,)
(2, 2)

Now I’m running tests with projections.


V_trial, _ = W.sub(0).collapse()

project(uh, V_trial)

#E_trial, _ = W.sub(1).collapse()

E_trial = dolfinx.fem.functionspace(mesh, elem_e)

project(epsh, E_trial)

E_trial, _ = W.sub(1).collapse()

project(epsh, E_trial)

The last one produces this error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[8], line 3
      1 E_trial, _ = W.sub(1).collapse()
----> 3 project(epsh, E_trial)

Cell In[2], line 29, in project(v, V)
     27 w = ufl.TestFunction(V)
     28 a = ufl.inner(u, w)*ufl.dx
---> 29 L = ufl.inner(v, w)*ufl.dx
     30 problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=bcs,
     31                         petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
     32 return problem.solve()

File ~/anaconda3/envs/fenicsx0.9/lib/python3.10/site-packages/ufl/operators.py:213, in inner(a, b)
    211 if a.ufl_shape == () and b.ufl_shape == ():
    212     return a * Conj(b)
--> 213 return Inner(a, b)

File ~/anaconda3/envs/fenicsx0.9/lib/python3.10/site-packages/ufl/tensoralgebra.py:166, in Inner.__new__(cls, a, b)
    164 ash, bsh = a.ufl_shape, b.ufl_shape
    165 if ash != bsh:
--> 166     raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
    168 # Simplification
    169 if isinstance(a, Zero) or isinstance(b, Zero):

ValueError: Shapes do not match: <ListTensor id=124046807197056> and <Argument id=124046785378048>
2

That was indeed unintended behavior of collapse! Fixed in Fix symmetric value shape by michalhabera · Pull Request #3535 · FEniCS/dolfinx · GitHub

Will be available in 0.10.0. For now, you can run your code with the nightly docker image to get your expected behavior, dolfinx/dolfinx:nightly.

3

Thank you very much for your reply. :smiley: