This implements element normals as (discontinuous) grid function.
When evaluated at a point `x` of an element `e`, it computes the
unit outward normal of the 'closest face' of this element.

* The 'closes face' is determined in local coordinates within the
  reference elements.
* The implementation also works for curved faces, since the normal
  of the real element is computed using the (potentially non-constant)
  covariant Piola transformation of the (constante) normal in the
  reference element.
* Notice that the function is an L^2 function in general. It is
  defined piecewise with non-unique value on intersections.
  On elements, with affine geometries, it is piecewise constant.
* This also adds a test checking for consistency with
  `intersection.unitOuterNormal()` when restricting the function
  to `intersection`. A grid with distorted cubes serves as test case
  for non-constant normals.
* This satisfied the standard grid function interface. Extending this
  to an interface for skeleton functions to be bound on intersections
  may be interesting. But such an interface was not defined so far.