Add FaceNormalGridFunction

This implements face 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.