Rather than just on the type for points in the domain.
The point type is not enough: For example, if the density involves
a coefficient function given as a GridView function, then the
density has to be able to bind the local coefficient function
to the correct element, for which it has to know the element type.

Integration domains can be elements (i.e., Codim<0>::Entity)
or Intersections.

The actual bind method will follow in a separate commit.

Currently, densities only get the type used for coordinates
of the integration domain.  However, this is not enough:
The complete domain must be known, and this domain can be either
a grid element or a grid intersection.

As a first small step, this patch introduces ElementOrIntersection
as a template parameter of the CosseratShellDensity class.

This is one step towards turning SurfaceCosseratEnergy into a
generic SurfaceIntegralEnergy class.

In particular, there is a directory dune/elasticy/densities now.

They were expected to be the factors of a ProductManifold which is
also a template parameter.

Expecting the product and the factors is redundant.  We simply
extract the factors from the product now.

Previously, it expected two factor spaces, as separate template
parameters.  Taking a single ProductManifold type is more consistent
with GeodesicFEAssembler, and it allows to eventually generalize
the code to more factors, if desired.

... but disallow a plain pointer.  This increases consistency with the
GeodesicFEAssembler class.

dune-elasticity also has densities, but they implement a slightly
different interface.  Instead of having special code for such
densities in the LocalIntegralAssembler, add a wrapper that
makes densities from dune-elasticity look like dune-gfe densities.

Previously this wouldn't compile, because GlobalP2Mapper only supported
1d and 2d grids when used without dune-parmg.  This has now been fixed,
so we can enable filmonsubstratetest and filmonsubstratetest-mixed
(which involve a 3d grid) also for the situation without dune-parmg.

Strangely, the test without the -mixed produces different results
depending on whether the dune-parmg version of GlobalP2Mapper is used
or not. Minor differences are to be expected because the precise
behavior of multigrid methods (which are used by the TR solver) depends
on the dof ordering.  But the difference seen here is a bit large
for my taste.

I don't quite understand this, so this commit introduces some
subcasing to make CI pass, and document the current behavior.

Plus some minor cleanup

Rather than using preprocessor switches.  That's too fragile
and not extensible.

If the target manifold is a product space, then ProductManifold must be
used, and not a variadic list of the factor space types.

Together with the new test, this commit also contains a number
of fixes for the SurfaceCosseratEnergy class.  This class was
not unit-tested previously, and unintentionally I messed it up
badly during my recent modernization work.

Finally, SurfaceCosseratEnergy had apparently never been tested
without dune-curvedgeometry, which is however only an optional
dependency.  The commit therefore also adds the necessary
conditionals.

Then problem then is that the derivative of the outer normal
of an intersection is not available from the standard grid interface.
If dune-curvedgeometry is not available, this derivative is
simply set to a zero matrix.  Depending on the grid manager
this may not actually be correct, though.