Error when trying to interpolate a tensor function coming from a mixed element

Hello,
I am working on a formulation for linear viscoelasticity with a dissipation potential, where I solve for the viscous strain and the displacement (see Linear viscoelasticity — Numerical tours of continuum mechanics using FEniCS master documentation)
Therefore I have a mixed element with a vector and tensor element. Following this post I want to interpolate the functions after solving into a lower function space. It works fine for the vector element but for the tensor element I get the following error message:

  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/function.py", line 151, in __init__
    self._ufcx_expression, module, self._code = jit.ffcx_jit(
                                                ^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py", line 62, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py", line 216, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 305, in compile_expressions
    raise e
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 285, in compile_expressions
    impl = _compile_objects(
           ^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 330, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/compiler.py", line 108, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py", line 101, in analyze_ufl_objects
    processed_expression = _analyze_expression(original_expression, scalar_type)
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py", line 132, in _analyze_expression
    expression = ufl.algorithms.apply_function_pullbacks.apply_function_pullbacks(expression)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_function_pullbacks.py", line 61, in apply_function_pullbacks
    return map_integrand_dags(FunctionPullbackApplier(), expr)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 86, in map_integrand_dags
    return map_integrands(
           ^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 79, in map_integrands
    return function(integrand)
           ^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 87, in <lambda>
    lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
    (result,) = map_expr_dags(
                ^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 101, in map_expr_dags
    r = handlers[v._ufl_typecode_](v)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/multifunction.py", line 29, in _memoized_handler
    r = handler(self, o)
        ^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_function_pullbacks.py", line 43, in form_argument
    raise ValueError(
ValueError: Expecting pulled back expression with shape '(4,)', got '(2, 2)'

I am using DOLFINx v09 in a Docker container (image dolfinx/dolfinx:stable). Here is my MWE

import ufl
import dolfinx as dlf
import numpy as np
from mpi4py import MPI 
from basix.ufl import element, mixed_element


mesh = dlf.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dlf.mesh.CellType.triangle)
dim = mesh.topology.dim

### DEFINE FUNCTION SPACES
Ve_L3 = element('Lagrange', mesh.basix_cell(), 3, shape=(dim,))
Te_L2 = element('Lagrange', mesh.basix_cell(), 2, shape=(dim,dim))

W = dlf.fem.functionspace(mesh, mixed_element([Ve_L3, Te_L2]))

V, Vmap = W.sub(0).collapse()
T, Tmap = W.sub(1).collapse()

w = dlf.fem.Function(W)
u, eps = w.split()

# For saving the solution
Ve_L1 = element('Lagrange', mesh.basix_cell(), 1, shape=(dim,))
Te_L1 = element('Lagrange', mesh.basix_cell(), 1, shape=(dim,dim))

# interpolating into lower Function space 
XDMF_V = dlf.fem.functionspace(mesh, Ve_L1)
XDMF_T = dlf.fem.functionspace(mesh, Te_L1)

XDMF_u = dlf.fem.Function(XDMF_V, name='XDMF_u')
XDMF_eps = dlf.fem.Function(XDMF_T, name='XDMF_eps')

XDMF_u.interpolate(dlf.fem.Expression(u, XDMF_V.element.interpolation_points())) 
XDMF_eps.interpolate(dlf.fem.Expression(eps, XDMF_T.element.interpolation_points()))    

Even if I try to use the same Tensor element in the definition of the Function space

XDMF_T_L2 = dlf.fem.functionspace(mesh, Te_L2)
XDMF_eps_L2 = dlf.fem.Function(XDMF_T_L2, name='XDMF_eps_L2')
XDMF_eps_L2.interpolate(dlf.fem.Expression(eps, XDMF_T_L2.element.interpolation_points())) 

or collapsing the function space and using the subspace for the definition of a function

XDMF_eps_L2 = dlf.fem.Function(T, name='XDMF_eps_L2')
XDMF_eps_L2.interpolate(dlf.fem.Expression(eps, T.element.interpolation_points()))  

I get the same error. It seems, that there is a shape mismatch. How can I solve this?

I think you want to use:

u, eps = ufl.split(w)

rather than:

u, eps = w.split()

see this discussion for instance on the difference between the two: Function.split() vs Split(Function)

Also, in your particular problem, wouldn’t you want your viscous strain to be defined on a DG space ?

1 Like