Add Geometry::jacobianInverse()

This adds a method Geometry::jacobianInverse(localCoordinate) and a corresponding typedef Geometry::JacobianInverse to the Geometry interface class. The method is default implemented using the TransposedMatrixWrapper from dune-common which internally captures (by value) the matrix returned by Geometry::jacobianInverseTransposed().

This is marked WIP, because a few things need to be discussed:

• Currently the method is default-implemented in the interface class Dune::Geometry. Hence it will work for all grid implementations. But one may argue (I tend to do so) that it should be implemented directly in the geometry implementations. Otherwise a Geometry provided by a grid provides a different interface than e.g. a MultilinearGeometry. To avoid that all grid implementations need to be modified, one could conditionally activate the default-implementation depending on presence of an implemenation in the grid implementation.
• Which interface do we want to provide? Currently you can do exactly one thing: Multiply the matrix from the right to a FieldMatrix using operator*. This is exactly what you need for the chain rule. This is in some sense the 'linear operator interface' proposed by @peter, but for the case where you want to chain the operator instead of directly applying it. Is this sufficient?
• Maybe one also wants a cast to a suitable FieldMatrix?
• Should Geometry::jacobian() also be implemented in a similar way? I currently don't have a use case for this, but having both, jacobian() and jacobianTransposed() seems to be reasonable. This would e.g. allow to write g.jacobianTransposed(x)*g.jacobian(x) to compute the metric tensor which is conceptually nice (but probably inefficient).
• Tests are missing (in dune-grid) due to the open questions. But the TransposedMatrixWrapper is already tested in dune-common and I also tested with the dune-functions examples, where this allows to replacing the ridiculous transpose(geometry.jacobianInverseTransposed(x) by a plain geometry.jacobianInverse(x).

This is related to several other MRs that should be merged in appropriate order:

• dune-common!1100 (merged) provides a utility needed for the subsequent extension of transpose().
• dune-common!1101 (closed) extends the transpose() function to support capturing matrices by value.
• dune-common!1138 (merged) extends support for transpose(matrix) and matrix.transposed()
• !577 (merged) (present MR) extends the (grid-)Geometry interface class by the new interface. This is conditionally default-implemented: If the actual grid geometry provides the new interface it is used, otherwise a default implementation based on transpose() is used.
• dune-geometry!193 (merged) extends the actual geometry check and the implementations provided in dune-geomtry by the new interface.

dune/grid/common/geometry.hh

@@ -11,6 +11,8 @@
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/typetraits.hh>
+#include <dune/common/transpose.hh>
+#include <dune/common/std/type_traits.hh>
 
 #include <dune/geometry/referenceelements.hh>
 
@@ -126,6 +128,42 @@ namespace Dune
      */
     typedef typename Implementation::JacobianTransposed JacobianTransposed;
 
+  private:
+
+    template<class Implementation_T>
+    using JacobianInverseOfImplementation = decltype(typename Implementation_T::JacobianInverse{std::declval<Implementation_T>().jacobianInverse(std::declval<LocalCoordinate>())});
+
+    using JacobianInverseDefault = decltype(transpose(std::declval<JacobianInverseTransposed>()));
+
+    template<class Implementation_T>
+    using JacobianOfImplementation = decltype(typename Implementation_T::Jacobian{std::declval<Implementation_T>().jacobian(std::declval<LocalCoordinate>())});
+
+    using JacobianDefault = decltype(transpose(std::declval<JacobianTransposed>()));
+
+  public:
+
+    /**
+     * \brief type of jacobian inverse
+     *
+     * The exact type is implementation-dependent.
+     * However, it is guaranteed to have the following properties:
+     * - You can multilply it from the right to a suitable FieldMatrix
+     * - It is copy constructible and copy assignable.
+     * .
+     */
+    using JacobianInverse = Std::detected_or_t<JacobianInverseDefault, JacobianInverseOfImplementation, Implementation>;
+
+    /**
+     * \brief type of jacobian
+     *
+     * The exact type is implementation-dependent.
+     * However, it is guaranteed to have the following properties:
+     * - You can multilply it from the right to a suitable FieldMatrix
+     * - It is copy constructible and copy assignable.
+     * .
+     */
+    using Jacobian = Std::detected_or_t<JacobianDefault, JacobianOfImplementation, Implementation>;
+
     /** \brief Return the type of the reference element. The type can
        be used to access the Dune::ReferenceElement.
      */
@@ -267,6 +305,55 @@ namespace Dune
     {
       return impl().jacobianInverseTransposed(local);
     }
+
+    /** \brief Return the Jacobian
+     *
+     *  The Jacobian is defined in the documentation of
+     *  \ref Dune::Geometry::integrationElement "integrationElement".
+     *
+     *  \param[in]  local  position \f$x\in D\f$
+     *
+     *  \return \f$J_g(x)\f$
+     *
+     *  \note The exact return type is implementation defined.
+     */
+    Jacobian jacobian ( const LocalCoordinate& local ) const
+    {
+      if constexpr(Std::is_detected_v<JacobianOfImplementation, Implementation>)
+        return impl().jacobian(local);
+      else
+        return transpose(jacobianTransposed(local));
+    }
+
+    /** \brief Return inverse of Jacobian
+     *
+     *  The Jacobian is defined in the documentation of
+     *  \ref Dune::Geometry::integrationElement "integrationElement".
+     *
+     *  \param[in]  local  position \f$x\in D\f$
+     *  \return \f$J_g^{-1}(x)\f$
+     *
+     *  The use of this function is to compute the jacobians of some function
+     *  \f$f : W \to \textbf{R}\f$ at some position \f$y=g(x)\f$, where
+     *  \f$x\in D\f$ and \f$g\f$ the transformation of the Geometry.
+     *  When we set \f$\hat{f}(x) = f(g(x))\f$ and apply the chain rule we obtain
+     *  \f[\nabla f(g(x)) = J_{\hat{f}}(x) J_g^{-1}(x).\f]
+     *
+     *  \note In the non-quadratic case \f$\textrm{cdim} \neq \textrm{mydim}\f$, the
+     *        pseudoinverse of \f$J_g(x)\f$ is returned.
+     *        This means that its transposed is inverse for all tangential vectors in
+     *        \f$g(x)\f$ while mapping all normal vectors to zero.
+     *
+     *  \note The exact return type is implementation defined.
+     */
+    JacobianInverse jacobianInverse ( const LocalCoordinate &local ) const
+    {
+      if constexpr(Std::is_detected_v<JacobianInverseOfImplementation, Implementation>)
+        return impl().jacobianInverse(local);
+      else
+        return transpose(jacobianInverseTransposed(local));
+    }
+
     //===========================================================
     /** @name Interface for grid implementers
      */