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

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.

Use correct vector type for the combination of displacement and rotation in geodesicfeassemblerwrappertest

Add a simple test with a grid containing one element testing if the energy is invariant of the number of grid refinements.

This test assembles and solves a prototype Cosserat problem.
In then checks whether the required number of iterations and the
final energy match given (hard-coded) values.

This reflect the current purpose of the test much better.

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.

GeodesicFEAssembler can replace it without problems.

Only compile geodesicfeassemblerwrappertest if the version of dune-elasticity is greater or equal 2.7

Virtual Function is deprecated after dune 2.7.

Add a GeodesicFEAssemblerWrapper - it wraps a MixedGFEAssembler so it can be used like a GeodesicFEAssembler

The GeodesicFEAssemblerWrapper assembles the Gradient and the Hessian using the MixedGFEAssembler and then
resturctures them so they can be used with the normal RiemannianTRSolver.
This only works, if the FE spaces have the same order.

Previously, the method constructBoundaryDofs would accept scalar bases
together with blocked bit vectors, and would tacitly do
The Right Thing.

This has changed: Nowadays, the basis really has match the blocking
structure of the bit vectors.  This means that we have to introduce
a second type of power basis:  The one that already exists has the
dimension of the embedding space.  That is correct for sampling
initial configurations.  However, the Dirichlet values apply to
corrections, which live in the tangent space.  Therefore, their
'power-order' needs to be the dimension of that, i.e., lower.

Scalar bases used to work for that, but that feature got
removed a long time ago.

Do this by making RodAssembler more similar to the generic
GeodesicFEAssembler: the constructor argument is now a
general LocalGeodesicFEStiffness, and not something particular
to rod problems.

This is needed to test against the 2.6 core modules.

This new test computes a harmonic map from a square to a sphere.
In other words, rather than testing some specific aspect of some
class, it does a whole simulation.  It should therefore cover
much more of the dune-gfe code.

I wrote it without a clear idea of what is needed of such a
test function, and I never finished it.  Now it is unclear
what the code does or is supposed to do, so let's just
remove it.

Using LocalIndexSet objects was mandatory with dune-functions-2.5
and earlier.  However, we only support 2.6 and later anyway,
so let's get rid of the backward-compatibility code.

In the test, the grid dimension was defined twice:  Once globally
at the top of the file, and once through a chain of template parameters.
This patch removes the global variable.

Because it really is a functionality of rotation matrices in quaternion
representation.  Only the Cosserat energies use it at the moment, but
this is coincidental.

In the same process, the method is renamed to quaternionTangentToMatrixTangent.
That seems to be a more telling name.

In particular, we move the test for the computeDR method from
cosseratenergytest.cc to here.  It is a method that does not
have anything to do with Cosserat mechanics directly.

The test still seems overly complicated.  Do I really need
a LocalGeodesicFEFunction?  Maybe I can simplify it later.

* Do not test with configurations that result in zero-area elements
  (because the energy is inf in these cases).
* Test for the relative energy difference, rather than the absolute one,
  because the energy values can get very large.