Only include the multiplication of a ScaledIdentityMatrix with another FieldMatrix if dune common is older or equal 2.8

This triggers bugs in the Rotation and RigidBodyMotion classes.

In particular, it triggers what Jonathan Youett fixed in the
never-merged merge request

  !2

These bugs are fixed in this commit, too:
* The RigidBodyMotion class gets a `log` method
* The Rotation<3>::log method now returns EmbeddedTangentVector
  instead of TangentVector.

It was TangentVector whereas it should have been
EmbeddedTangentVector.

Calling log(a,b) returned a-b instead of b-a.

Also, enable unit tests for the RealTuple class.  The current test
will not trigger the bug in the log method: that will come shortly.

Finally, this commit is the first that mentions its user-facing
changes in the newly introduced CHANGELOG.md file.

The compiler warns about &, and I think the author's intention
was && anyway.

This involves fixing the secondDerivativeOfDistanceSquaredWRTSecondArgument
method: Previously it returned a FieldMatrix, but the test requires
a SymmetricMatrix.

AlexanderMueller authored 3 years ago and Sander, Oliver committed 3 years ago

This commit adds the energy implementat for the nonlinear Reissner-Mindlin
shell model known as the Simo and Fox shell model.

Details may be found in

  "A consistent finite element formulation of the geometrically
   non-linear Reissner-Mindlin shell model [Müller, Bischoff]"

Move the polar decomposition to a separate file and add the algorithm by Higham and Noferini (from Robin Fraenzel)

In the test for the Higham and Noferini algorithm, we see:
Unfortunately, for matrices that are close to an orthogonal matrix, this algorithm is
about 2x slower. One can benefit from the Higham and Noferini algorithm only if the
matrix is quite far away from an orthogonal one.

Before, the calculation was terminated after exactly three iterations, now the
stopping criterion is fulfilled if the change in the norm of the polar factor is
small enough.

Correct the sign of the Neumann and volume forces to make it consistent throughout dune-gfe and dune-elasticity

In the energy file, the energies are *added*.
The Neumann and volume forces get the correct sign in the actual program file by mulitplying with (-homotopyParameter).

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 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.