Previously, LocalEnergy would accept an arbitrary number of
template arguments, and (with the exception of the first one)
they would be interpreted as the factors of a product target space.
This patch replaces this template list by a single template
parameter, which has to be ProductManifold if a product space
is desired.

This has implications throughout the code.  In particular, there are
now two energy methods: One that still accepts coefficients sets
in the form

  std::vector<TargetSpace>

and a second one which accepts

  TupleVector<std::vector<Factors>...>

The second one really only makes sense for product manifolds.
However, as the 'energy' methods are pure virtual, they cannot be
disabled by SFINAE or C++20 concepts.  Therefore, the second 'energy'
method exists always, accepting

  TupleVector<std::vector<TargetSpace> >

if TargetSpace is not a product.

In the long run this is cleaner than setting up the type in
the LocalStiffness class itself.

This will later simplify merging the mixed and non-mixed versions
of the assembler code.

The rest of the HarmonicEnergy class is now independent of
the density.  It is functionally equivalent to the LocalIntegralEnergy
class, and shall be replaced by that eventually.  However,
LocalIntegralEnergy still has the deformation/rotation-split
hard-coded in a few places, and that needs more work.

Currently, Cosserat materials are hard-wired in the interface.
This patch changes the interface a little bit towards being
more general.

That's another step closer towards merging the
MixedLocalGeodesicFEStiffness and the LocalGeodesicFEStiffness.
Also, it is necessary to be later able to generalize the code
for product spaces with more than two factors.

It is deprecated, and not used anymore.

... rather than two separate spaces.  This is one step towards
getting rid of it altogether.

... regarding ambiguous logical expressions.

The ProductManifold class generalizes RigidBodyMotion, and can do
everything that the RigidBodyMotion class can.  Therefore there is
no point in keeping RigidBodyMotion any longer.  Having two
implementations for the same thing will just confuse people.

It contained code for mixed-dimensional domain decomposition algorithms
involving Cosserat materials.  This code has moved to a separate modul
'dune-gfe-coupling' a long time ago, and it is not used from anywhere
else in dune-gfe.

They should be installed, but since dune-gfe doesn't have any
downstream dependencies, nobody has missed the files yet.

It was essentially a reimplementation of Rotation::orthonormalFrame,
just with a different interface.  It was not directly used by any
code.

Divide the regularization parameter through the scaling parameter in the PN solver instead of multiplying

In the PN solver, we add identity*regularization/scaling to the stiffness matrix.
This has two effects:
- make the system matrix positive definite
- penalize large corrections to ensure an energy decrease in each step

The regularization parameter is adapted if either the system matrix is not positive definite
if there is an energy increase. This works similar to the trust-region radius in the TR solver.
The correction must be inside the trust-region, this ensures there is an energy decrease
in each step.
For each component of the correction, the trust-region radius as well as the regularization parameter
is multiplied with a scaling factor, because the correction for the displacement works depends
on the grid size where the correction for the rotation does not.

To align the effect of the scaling parameter, we need to divide through the scaling parameter in
the PN solver and multiply with it it in the TR solver.

... by LagrangeLocalFiniteElementCache.

It was apparently not used.

If you should ever need something like this in the future,
please use periodicbasis.hh from dune-functions.

This used to work A Long Time Ago, but the feature was removed
eventually.  Unfortunately, the code still used to build,
but failed to produce the correct result.  More recently,
some other commit in dune-functions apparently now leads to
a build failure when trying to call 'interpolate' for a
vector-valued function and a scalar basis.  This is nice!
It allows to catch errors earlier.

The old style has been deprecated for a long time, and is not
supported by dune-fufem anymore.

(In particular, the code discretizationerror.hh now requires
functions to satisfy the dune-functions function interface.)

In the wake of this, compute-disc-error.cc gets a number
of changes.  Let's hope I didn't screw up anything.

There are two important changes:  First of all, the new class implements
the `dune-functions` interface rather than the deprecated one based
on inheritance from `VirtualGridViewFunction`. Secondly, the new class
does not hard-wire geodesic interpolation anymore.  Rather, you can
give it interpolation classes which then govern how interpolation is done.

This change was long overdue.  It was now needed to make
gradient-flow.cc compile with the post 2.9 version of dune-fufem again.

In particular:
- Use the CosseratVTKWriter also for MIXED_SPACE = 1, only print out the deformation though
- Use a SubsamplingVTKWriter only for dim==gridDim==3