Use cholmod if the dune-solvers version >= 2.8, use umfpack for older versions

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.

Add multiplication of a ScaledIdentityMatrix with another FieldMatrix to the collection in linearalgebra.hh

The tests failed since the construction of values of ProductManifold<> in ValueFactory uses random entries between [0.9..1.1]. These are then used for the tests and are projected onto the manifold. 

To pass this tests it is needed to adjust several Taylor expansions and thresholds. For example this commit increases the threshold from `1e-4` to `1e-2` and adds more terms to the Taylor expansion. The reason for this change is explained in the following. 

The problem is, even if the tolerance `1e-4` is sufficient to have a correct function value within machine precision, it is not always sufficient to get correct derivatives since here we lose orders of correctness.

For  example of for sinc(x) we need to put the threshold at `1e-4` to get the correct function value if we use  `1.0-x*x/6.0` as approximation formula.

If we then use this formula within automatic differentiation or finite differences, e.g. the derivative algorithms  "sees" only the following formulas of the first and second derivative:

- Function value: `1.0-x*x/6.0`   This function is implemented
- First deriv: `-x/3.0`      This sees the derivative algorithms as first derivative
- Second: `-1.0/3.0`       This sees the derivative algorithms as second derivative

Obviously, this is the case if the function value is inside the region where the Taylor expansion is used.

If we use these functions to test the exactness of the derivatives we need to set the threshold to `x<1e-6` to get exact second order derivatives where the error is within machine precision. Therefore, for larger values `x>1e-6` the exact formula has to be used to get correct results. Unfortunately, in this range the exact derivative are already unstable.
E.g. the first derivative formula behaves already strange near `x=1e-4`.

Therefore,
to get a correct derivative value we need to switch to the Taylor expansion earlier (`1e-4`) to prevent using the unstable exact formula. But in this range the Taylor expansion is unable to reproduce an approximation error within machine precision.

I think the only way to fix this problem is to add more terms to the Taylor expansions. Since even if they seem to be sufficient in terms of function value, usually they are not sufficient in terms of derivatives.

For `sin(x)/x` a test ist at https://godbolt.org/z/T995hGec3 and for `acos(x)^2` https://godbolt.org/z/TG9E15jjf.

The StressFreeStateGridFunction is the actual parametrization from the grid to the
nonplanar Cosserat shell in stress-free state.
The geometries of the StressFreeStateGridFunction will then be used for all calculations
instead of the linear geometries of the piecewise linear grid.

Otherwise the Cosserat rod code will not compile without dune-parmg.

So far, the Cosserat rod energy implementation hat a first-order
finite element space hardcoded.  This patch removes that restriction.
As for the other models in this Dune module, the finite element basis
is now a template parameter of the model energy, and can be set to
any reasonable basis.

It is not a 'stiffness' anymore, but really only implements the
energy.

Previously, the rod3d code used a hand-implemented first derivative
of the rod energy, and FD for the second derivatives.  This patch
replaces this by ADOL-C for both first and second derivatives.
Advantages are:
* The code is much shorter, and easier to understand.
* The previous code contained interpolation formulae for
  first-order 1d geodesic finite elements.  This restricted
  rod problems to first-order GFE.  A follow-up patch will
  allow higher-order approximations.
* The new code runs roughly twice as fast.

GeodesicFEAssembler can replace it without problems.

As a first step towards getting rid of the RodAssembler class,
this patch moves the getStrain, getStress, and getResultantForce
methods to the RodLocalStiffness class.  I am not 100% convinced
that that is the best place for them, but I can't think of a
better one right now.

ADOL-C implements most standard math functions for its
'adouble' type, and they are not in the namespace 'std'.
dune-fufem contains a file adolcnamespaceinjections.hh
which imports some of these methods into the 'std'
namespace, but that is not the proper way to do it.
Rather, they should be found by argument-dependent
lookup (ADL), that is

  using std::sin;
  auto v = sin(x);

When dune-curvedgeometry is not found, the normal Derivative is set to zero.

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

Add localintegralenergy.hh and neumannenergy.hh, they assemble the energy for a single element and work with a CompositeBasis

These classes work similarly to Dune::Elasticity::LocalIntegralEnergy and Dune::Elasticity::NeumannEnergy,
where the ones in Dune::Elasticity with a power basis (for the displacement function) and
Dune::GFE::LocalIntegralEnergy / Dune::GFE::NeumannEnergy work with a CompositeBasis (for the displacement AND the rotation),
so the indices to access the different parts of the basis are different.

Add sumenergy.hh, this class assembles the a sum of energies for a single element by summing up the energies of each GFE::LocalEnergy.

This class works similarly to the class Dune::Elasticity::SumEnergy, where Dune::Elasticity::SumEnergy extends
Dune::Elasticity::LocalEnergy and Dune::GFE::SumEnergy extends Dune::GFE::LocalEnergy.