Let’s say, we want to solve a problem in linear elasticity with an anisotropic material:

\nabla\cdot S = 0\qquad +\text{bcs.}

Here, S denotes the stress, with S^{ij} = C^{ijkl}\varepsilon_{kl}, \varepsilon is the linearized Green-Lagrange strain \varepsilon(u) = \frac 1 2 \left(\nabla u + \nabla u^T\right) =: \operatorname{sym}\left(\nabla u\right), u is the displacement and C is the tensor of elasticity with rank four.

This tensor exhibits several symmetries: C^{ijkl} = C^{jikl} = C^{ijlk} = C^{klij}.

Additionally, it is constant.

The (symmetric) bilinear form then becomes

a(u, v) = \int_\Omega \left(C\varepsilon(u)\right):\varepsilon(v)\,\mathrm d\Omega.

The question is now, how to properly implement the tensor of elasticity.

The nicest thing would be, to be able to write the following:

```
a = inner(dot(C, sym(nabla_grad(u))), sym(nabla_grad(v)))
```

Additionally, it would be nice, if FEniCS was able to determine, that this expression is symmetric.

I found some old answer here: https://fenicsproject.org/qa/10309/small-question-on-defining-constant-tensors/

The solution there is to implement a subclass of `ufl.Coefficient`

and `cpp.Constant`

and make it a symmetric `TensorElement`

.

However, I found this old bug, that `symmetry=True`

is totally broken for `TensorElement`

: https://bitbucket.org/fenics-project/ffc/issues/78/tensorelement-symmetry-true-is-totally

The bug is quite old, but the discussion itself took more than two years and I cannot find, whether it is now fixed or not.

An alternative would be to use Voigt-Notation, as outlined in this tutorial here:

https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/orthotropic_elasticity.py.html

In a two dimensional setting, exploiting all the symmetries, the equation C\varepsilon(u) can be written as a matrix vector product:

\begin{bmatrix}S^{1}\\S^{2}\\S^{3}\end{bmatrix} = \begin{bmatrix}C^{11}&C^{12}&2C^{13}\\C^{21}&C^{22}&2C^{23}\\C^{31}&C^{32}&2C^{33}\end{bmatrix}\cdot\begin{bmatrix}u_{1}\\u_{2}\\u_{3}\end{bmatrix}

The indices 1, 2 and 3 correspond to xx, yy and xy respectively.

For three dimensions this becomes a 6\times 6 system.

A corresponding FEniCS implementation would use `as_vector`

to decompose the trial function into its components, multiply that vector with the matrix representation of C and afterwards use `as_tensor`

to reinterpret the result as a tensor.

However, while this works, the information about the symmetry is lost.

I know, it is possible to simply tell the linear solver, that the system is symmetric.

But I have read at other places, that is seems to be beneficial, when FEniCS already knows about the symmetry when assembling the system.

Is this true?

So: How would you implement a symmetric fourth rank tensor?