Introduce preprocessor variable PROJECTED_INTERPOLATION to switch between geodesic FE and projection-based FE

This is a temporary solution; we need something smarter eventually.

Rationale
a) Code is much simpler without it
b) Code compiles faster without it
c) P3 functions are now downsampled to P2 functions, rather than to P1 functions.
   Hence there is less loss of information now.

Simply ask the first local finite element of the dune-functions basis for the order.

In particular, assemblers now hand LocalView objects around, rather than pairs of
elements and LocalFiniteElements.

Introduce a dedicated class LocalQuickAndDirtyFEFunction to produce initial iterates for the GFE minimization problems

Previously, a LocalProjectedFEFunction was used for this.  This turned out to be a bad idea for two reasons:
- Evaluating such a function can actually pretty expensive.  For example,
  to project onto SO(3), the polar decomposition needs to be evaluated by
  an iterative method.
- Screwing up the LocalProjectedFEFunction implementation can actually screw up
  the GFE implementation as well in very subtle ways.  Specifically, if the
  LocalProjectedFEFunction is not ADOL-C-differentiable, then the GFE function
  will not be differentiable, either.  This fact cost me a few hours searching...

Because these two classes are for the most part identical.

Currently, this merged version involves some trickery, including some template
specializations and multiple inheritances.  I hope to get rid of that, eventually.

This is how the new dune-functions bases should be used.  Also, it is yet
another step towards merging CosseratEnergyLocalStiffness and MixedCosseratEnergy.

This model was given to me by Patrizio Neff in private communication.

The code is a prototype.  It appears to work, but it has not been tested very
thoroughly, and quite likely there are some bugs left.  Also, the code contains
the shape of the reference configuration hard-wired (a sphere).

This patch adds several more implementations of operator+ etc for FieldMatrices.
In particular, you can now combine matrices of adoubles (the 'active' type from
the ADOL-C automatic differentiation library) with regular double types.

These new operators are needed for the NonplanarCosseratShellEnergy class,
which is coming in the next patch.

The old version kept was too big, and kept a useless row of zeros around.
Removing this row should lead to a tiny increase in efficiency.  More importantly,
it gets us closer to the CosseratEnergyLocalStiffness class again, which
is scheduled to replace MixedCosseratEnergy eventually.

This leads to some simplification

This is an important step towards making the 'mixed' Cosserat code more like
the regular one.  The mid-term goal is to merge the two implementations, because
there is way too much code duplication here.

The DR variable (the derivative of the microrotation R wrt space) had dimensions 3x3x3.
This is two big in a shell context, because there are actually only two space dimensions.
Hence the correct size is 3x3xgridDim.

The results without this patch were not wrong, though, only a bit wasteful.

It is needed for shell problems with non-planar reference configuration.

This makes 2nd-order orientation functions work again.

The implemented formula was actually for the projection of Stiefel manifolds, of which
the orthogonal matrices are only a special case.  For this special case, one term
drops out and doesn't have to be implemented.

... with values expected to be given in the parameter file.

Scaling the trust region is the key to success, because in general the midsurface
deformation and the rotations will have values in different ranges.

This marks the last step of porting the 'cosserat-continuum' application away
from the dune-fufem bases.  One must not(!!!) mix dune-fufem bases with
dune-functions bases, as they use different dof orderings at least in the P2 case.