[forms][bugfix] Fix GeometryType in QuadratureRuleKey exported by Forms::Coefficient

See merge request !261

[doc] Immprove documentation of AffineConstraints class

See merge request !262

So far a `Forms::Coefficient` exported a `QuadratureRuleKey` with
`GeometryType` being `none` with appropriate dimension. This works
when assembling linear and bilinear forms, because we always multiply
with at least one unary operator which then sets the correct `GeometryType`.

However, if we want to use the nullary operator `Forms::Coefficient`
standalone, e.g. to simply integrate it (it's a grid function after
all), we silently use the wrong `GeometryType`. Luckily we can always
derive the correct `GeometryType` from the element this function is
bound to allowing for an easy fix.

[examples] Add examples with Hermite-type bases

See merge request !238

This adds three examples:

* 2d Poisson equation with cubic Hermite-basis and weak enforcement of
  Dirichlet BC (Nitsches method).
* 2d Poisson equation with cubic Hermite-basis and strong enforcement
  of Dirichlet BC utilizing the `FunctionalDescriptor` interface.
* 2d biharmonic equation with nonconforming Morley basis and strong
  enforcement of clamped BC utilizing the `FunctionalDescriptor` interface.

[cmake] Do not add python flags unconditionally

See merge request !260

This add a patch from !239
and also removes some 2.7 backward compatibility code, since we only
support dune versions >= 2.9 which all require Python3.

Add utility function writeBasis()

See merge request !259

This plots writes all basis functions as vtk sequence
for viewing them with paraview. This is _very_ helpful
when investigating/debugging a dune-functions basis
and affine constraints e.g. for hanging nodes, Dirichlet
BC, ... .

Ideally this should be part of dune-functions, but the
latter does not depend on dune-vtk and furthermore does
not know about constraints so far.

Add constraints support for dune-functions basis

See merge request !258

The examples requires a fix that is not contained in 2.10
and makes composed grid functions work with bound unary
operators.

In order to demonstrate the (hanging-node) constraints
work with higher order bases, too, use 2nd order FEM
as ansatz space and modify the problem such that the
boundary data are not representable exactly in the FE
space.

As extension space we use 4th order finite elements.
By fixing certain DOFs we still get incremental
bubbles wrt. to the P2-basis.

This makes the rendering of Latex more consistent
with the rest of the documentation and improves readability.

Document what exactly `AffineConstraints` is supposed
to do and how it is intended to be used and rename
`constrainMatrixPattern()` to extendMatrixPattern()`
since this in fact enlarges the pattern and never
reduces it.

This avoids having to construct a boundary patch if we want
to process the whole boundary anyway.

If we have `markBoundaryPatchDofs()` there should also be
`markBoundaryDofs()` to avoid having to construct a `BoundaryPatch`
in the trivial cases.

* Adjust to updated interface of `AffineConstraints`.
* Separate construction of constraints object from
  computation of its content.

Before the class only supported linear constraints.
By adding support for constant offsets, this can now
be used e.g. for hanging node constraints (linear),
Dirichlet constraints (constant), or both at once.

The algorithm basically interpolates fine DOFs from
coarse basis functions of neighboring elements. This
requires some mechanism to interpolate only hanging fine
DOFs and not all.
Instead of checking if the fine interpolation point is in
within the coarse element, we now check if the interpolated
DOF is associated to a subentity of the intersection.
In most cases this is equivalent. But for hierarchical bases,
the former does probably not lead to the desired result (*).

Remark: (*)
In principle one would expect that the conforming hierarchical
basis consists of the conforming P1-basis plus conforming bubbles.
This is the case for both algorithms (except for some numerical
issues for the old one). However, it turnes out that there are
different possibilities to define the conforming bubbles and
it's not clear which one is desired. Ideally we would like to have:

* The bubbles are lagrangian, i.e. each one is zero at the interpolation
  point of the other ones.
* The bubbles are non-negative.

For hanging nodes you can't have both at once. The old algorithm
guarantees the former, the new one the latter property. I'm undecides
which one is 'correct'.

This demonstrates the usage of hanging node.
So far the refinement indicator is a very crude
quick and dirty hack.

This adds a class `AffineConstraints` for storing contraints
for general dune-functions bases and a method
`makeContinuityConstraints()` for computing the
constraints for the H^1-conforming subspace
for classical FE spaces on grid with hanging nodes.

This is experimental and the interface may change.
Also some features are still missing:

* The code is generic and should work with nested
  bases and subspace bases. But this has not been
  tested so far.
* One may want to generalize `makeContinuityConstraints()`
  such that one can extend an existing constraints
  object. This way one could add constraints for
  different branches in the tree one after the other.
* Currently constraints are stored in a simple
  `std::unordered_map<std::map<Index,Coeff>>`
  data structure. Maybe this could be more efficient
  in a vector-based data structure or a sparse matrix.
  I had a compress method in mind that transforms into
  a more efficient storage once the contraints are fully
  build. But in my tests copying into a vector-bases
  data structure did not improve significantly, probably
  because the number of entries is relatively small.
* The constraints class is called `AffineConstraints`
  because it is intended to also hold Dirichlet constraints.
  This feature is currently missing.

Add utility functions for recursive traversal of matrices and vectors

See merge request !257

This has a few advantages:
* It does not require the basis just to do algebraic transformations.
* It is also correct if the pattern is not exactly the one given by
  the local basis coupling. This is e.g. important for hanging nodes.
* It is faster, because we loop over the matrix entries just once
  and in consecutive order.

This allows

```cpp
recursiveForEachMatrixEntry(matrix, [](auto&& matrix_ij, auto&& i, auto&& j) {
  ...
});

recursiveForEachVectorEntry(vector, [](auto&& vector_i, auto&& i) {
  ...
});
```

Here the container is traversed hierachically and the row and column multi-indices
`i` and `j` are instances of `Dune::Functions::StaticMultiIndex` with appropriate
length depending on the nesting depth.

[examples] Fix matrix/vector types of examples

See merge request !256

The nested matrix and vector types should correspond
to the index tree of the basis. Here, this is not the
case since `Field*<T,*>` was used instead of `T`.
While this is in many cases no problem, it may fail
if the matrix is recursively traversed with `IsNumber`
as termination criterion.

[cleanup] Remove experimental header arithmetic.hh

See merge request !255

This header was added a long time ago to replace `staticmatrixtools.hh`
and was later itself superseded by dune-matrix-vector. Furthermore there
was also a copy of this header in dune-solvers which was deprecated
since 8 years. The header is not used anywhere and is thus removed.
Since it was still documented as experimental, removal does not require
a deprecation here.

Improve performance of assembleGlobalBasisTransferMatrix

See merge request !254