UFL tensors and matrices

UFL
1

Hi all!

I’m trying to modify this shell demo to the formulas (23), (24) and (25) of this paper that involve the contraction of a matrix and a 4th order tensor, and I get the following error

ValueError: Index to be removed (11) not part of indices ((13, 14)).

with the MWE code below. It’s probably due to a misunderstanding of the UFL types on my end, but I don’t really see where. Any help would be appreciated!
Thanks!

import math
from mpi4py import MPI
import ufl
from ufl import dx, grad, inner, split, sym, tr, Dx
from dolfinx.mesh import CellType, create_unit_square

# Parameters
E1, nu21 = 2.0e6, 0.45
nu12 = 0.45
E2 = nu21 * E1 / nu12
G = math.sqrt(E1 * E2) / (2 * (1 + math.sqrt(nu12 * nu21)))

# Geometry
mesh = create_unit_square(MPI.COMM_WORLD, 32, 32, CellType.triangle)

x = ufl.SpatialCoordinate(mesh)
phi0_ufl = ufl.as_vector([x[0], x[1], 0.])
a0_ufl = grad(phi0_ufl).T * grad(phi0_ufl)
a0_contra_ufl = ufl.inv(a0_ufl)
i, j, l, m = ufl.indices(4)  


# Elasticity tensor

den = (1 - nu12 * nu21)
K_tensor = ufl.as_matrix([ [ E1 / den, G / (1 - nu12) ], [ G / (1 - nu12), E2 / den ] ])
A_contra_ufl = ufl.as_tensor(
    (
        K_tensor[i, j] * (
            nu12 * a0_contra_ufl[i, j] * a0_contra_ufl[l, m]
            + 0.5 * (1 - nu12) * (
                a0_contra_ufl[i, l] * a0_contra_ufl[j, m]
                + a0_contra_ufl[i, m] * a0_contra_ufl[j, l])
        )
    ),
    [i, j, l, m],
)