The contravariant base vectors were not calculated correctly.
The contravariant base vectors are the *columns* of the inverse of the covariant matrix, not the rows.
To fix this, take the rows of the transpose of inverse of the covariant matrix.

Use GFE::LocalEnergy, GFE::LocalIntegralEnergy, GFE::SumEnergy (and delete GFE::SumCosseratEnergy).

Several energy density implementations contained code for matrix methods
like sym, skew, transpose etc.  This patch moves them all to the file
linearalgebra.hh, to get rid of the redundancy.

Reduce quadrature order, before 16 quadrature points were used, an amount that's sufficient to integrate over a 3D-element, but as the integral is on the boundary, now 9 quadrature points, which are enough for 2D, are used.

Add option to read in a displacement file that will be used for the stress-free configuration of the shell elements
When startFromFile is set to true, then the stress-free configuration is read from the given file - this file needs to match the *finest* level of the grid
When startFromFile is set to false, the stress-free configuration can be given in as a function

In principle, a serious bug: but for dimworld=3d (i.e., the standard)
this loop was never executed anyway.  Hence there were no serious
consequences.

Previously, they all inherited from LocalGeodesicFEStiffness,
which was weird, because none of them implemented first or
second derivatives as you would expect from something called
'stiffness'.

Added program to solve a combined minimization problem:
Minmize the sum of an elastic energy (representing the substrate)
and a cosserat energy (representing a thin film on the substrate).