Voigt-Notation distinguishes in the transformation from Matrix to Vector between stresses and strains.
The transformation for strains features an additional factor 2 for the non-diagonal entries. In order to avoid
the use of different data structures for both stresses & strains we use the same Matrix-to-Vector
mapping ('matrixToVoigt') and incorporate the factors in suitable places. namely:
  - The Stiffness matrix of the constitutive relation gets scaled by a factor of 2 in the last three columns
  - The 'voigtScalarProduct' scales the last three products by a factor of 2

Since all element geometries are identical, the element stiffness
matrices only depend on the phases at the quadrature nodes.
We use those to access a cache.  This leads to tremendous time savings.

This will be the cache key.

That's the first step towards caching the element matrices.

That's cheaper than first computing the symmetric part as a matrix,
and then converting that to Voigt notation.

This again saves a few conversions.

This avoids many transformations from symmetric matrices to
Voigt vectors and back.  In my (limited) testing, this reduces
the time to assemble the stiffness matrix by about 25%.

This patch also introduces a custom scalar product method for
Voigt vectors, which reproduces the Frobenius scalar product
in matrix space.  That way, the potentially confusing distinction
between stress-like Voigt vectors and strain-like Voigt vectors
can be avoided.

There is no need for this.  At least, I don't see any.

Yes, it is debatable whether that change makes the code more readable.

Use the corresponding member method of FieldMatrix.

That's Dune convention.

This saves another 15% run-time, at least in debug builds.

This saves about 20% run time.

dune-grid 2.9 introduced the Geometry::jacobianInverse method,
which  allows to transform shape function Jacobians from
reference element coordinates to world coordinates in a much
simpler way.

That's easy for the element stiffness matrix.  For the other matrix
and vector types it has to wait, because not everything that we use
from dune-istl support BCRSMatrix<double> yet.

They are not actually used.

Computing the material phase at a particular point may involve
a call to a Python function.  That is usually quite costly.
Therefore, make sure that this call is really made only once
per quadrature point.

Because such statements bleed into all the files that include
the headers, and this may lead to compilation errors later.