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?