Merge branch 'feature/jacobianinverse' into 'master'

Extend jacobian interface of geometry classes

See merge request !193

CHANGELOG.md

 # Master (will become release 2.9)
 
+- The `Geometry` interface was extended by methods `jacobian(local)` and `jacobianInverse(local)`
+  and corresponding typedefs `Jacobian` and `JacobianInverse`. This is implemented by all geometry
+  implementations provided by dune-geometry. But external implementations need to be adjusted
+  to pass the interface check provided by `checkgeometry.hh`.
+
 - The `Geometry::integrationElement` now needs to return the type `Volume`
   instead of `ctype`. Note that this may be different from `ctype` if the
   geometry supports typed dimensions. In such case, `ctype` is a length, and not

dune/geometry/affinegeometry.hh

@@ -506,6 +506,12 @@ namespace Dune
     /** \brief Type for the transposed inverse Jacobian matrix */
     typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
 
+    /** \brief Type for the Jacobian matrix */
+    typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
+
+    /** \brief Type for the inverse Jacobian matrix */
+    typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse;
+
   private:
     //! type of reference element
     typedef Geo::ReferenceElement< Geo::ReferenceElementImplementation< ctype, mydimension > > ReferenceElement;
@@ -640,6 +646,28 @@ namespace Dune
       return jacobianInverseTransposed_;
     }
 
+    /** \brief Obtain the Jacobian
+     *
+     *  \param[in]  local  local coordinate to evaluate Jacobian in
+     *
+     *  \returns a copy of the transposed of the Jacobian
+     */
+    Jacobian jacobian ([[maybe_unused]] const LocalCoordinate &local) const
+    {
+      return jacobianTransposed_.transposed();
+    }
+
+    /** \brief Obtain the Jacobian's inverse
+     *
+     *  The Jacobian's inverse is defined as a pseudo-inverse. If we denote
+     *  the Jacobian by \f$J(x)\f$, the following condition holds:
+     *  \f[J^{-1}(x) J(x) = I.\f]
+     */
+    JacobianInverse jacobianInverse ([[maybe_unused]] const LocalCoordinate &local) const
+    {
+      return jacobianInverseTransposed_.transposed();
+    }
+
     friend ReferenceElement referenceElement ( const AffineGeometry &geometry )
     {
       return geometry.refElement_;

dune/geometry/axisalignedcubegeometry.hh

@@ -88,6 +88,24 @@ namespace Dune {
         DiagonalMatrix<ctype,dim>,
         FieldMatrix<ctype,coorddim,dim> >::type JacobianInverseTransposed;
 
+    /**
+     * \brief Return type of jacobian
+     *
+     * This is a fast DiagonalMatrix if dim==coorddim, and a FieldMatrix otherwise.
+     * The FieldMatrix will never contain more than one entry per row,
+     * hence it could be replaced by something more efficient.
+     */
+    using Jacobian = std::conditional_t<dim==coorddim, DiagonalMatrix<ctype,dim>, FieldMatrix<ctype,coorddim,dim> >;
+
+    /**
+     * \brief Return type of jacobianInverse
+     *
+     * This is a fast DiagonalMatrix if dim==coorddim, and a FieldMatrix otherwise.
+     * The FieldMatrix will never contain more than one entry per row,
+     * hence it could be replaced by something more efficient.
+     */
+    using JacobianInverse = std::conditional_t<dim==coorddim, DiagonalMatrix<ctype,dim>, FieldMatrix<ctype,dim,coorddim> >;
+
     /** \brief Constructor from a lower left and an upper right corner
 
         \note Only for dim==coorddim
@@ -184,7 +202,7 @@ namespace Dune {
       return result;
     }
 
-    /** \brief Jacobian transposed of the transformation from local to global coordinates */
+    /** \brief Inverse Jacobian transposed of the transformation from local to global coordinates */
     JacobianInverseTransposed jacobianInverseTransposed([[maybe_unused]] const LocalCoordinate& local) const
     {
       JacobianInverseTransposed result;
@@ -196,6 +214,18 @@ namespace Dune {
       return result;
     }
 
+    /** \brief Jacobian of the transformation from local to global coordinates */
+    Jacobian jacobian([[maybe_unused]] const LocalCoordinate& local) const
+    {
+      return jacobianTransposed(local).transposed();
+    }
+
+    /** \brief Inverse Jacobian of the transformation from local to global coordinates */
+    JacobianInverse jacobianInverse([[maybe_unused]] const LocalCoordinate& local) const
+    {
+      return jacobianInverseTransposed(local).transposed();
+    }
+
     /** \brief Return the integration element, i.e., the determinant term in the integral
                transformation formula
      */

dune/geometry/multilineargeometry.hh

@@ -201,6 +201,12 @@ namespace Dune
     //! type of jacobian inverse transposed
     class JacobianInverseTransposed;
 
+    /** \brief Type for the Jacobian matrix */
+    typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
+
+    /** \brief Type for the inverse Jacobian matrix */
+    typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse;
+
   protected:
 
     typedef Dune::ReferenceElements< ctype, mydimension > ReferenceElements;
@@ -389,6 +395,29 @@ namespace Dune
       return geometry.refElement();
     }
 
+
+    /** \brief Obtain the Jacobian
+     *
+     *  \param[in]  local  local coordinate to evaluate Jacobian in
+     *
+     *  \returns a copy of the transposed of the Jacobian
+     */
+    Jacobian jacobian (const LocalCoordinate &local) const
+    {
+      return jacobianTransposed(local).transposed();
+    }
+
+    /** \brief Obtain the Jacobian's inverse
+     *
+     *  The Jacobian's inverse is defined as a pseudo-inverse. If we denote
+     *  the Jacobian by \f$J(x)\f$, the following condition holds:
+     *  \f[J^{-1}(x) J(x) = I.\f]
+     */
+    JacobianInverse jacobianInverse (const LocalCoordinate &local) const
+    {
+      return jacobianInverseTransposed(local).transposed();
+    }
+
   protected:
 
     ReferenceElement refElement () const
@@ -513,6 +542,8 @@ namespace Dune
 
     typedef typename Base::JacobianTransposed JacobianTransposed;
     typedef typename Base::JacobianInverseTransposed JacobianInverseTransposed;
+    typedef typename Base::Jacobian Jacobian;
+    typedef typename Base::JacobianInverse JacobianInverse;
 
     template< class CornerStorage >
     CachedMultiLinearGeometry ( const ReferenceElement &referenceElement, const CornerStorage &cornerStorage )
@@ -660,6 +691,28 @@ namespace Dune
         return Base::jacobianInverseTransposed( local );
     }
 
+    /** \brief Obtain the Jacobian
+     *
+     *  \param[in]  local  local coordinate to evaluate Jacobian in
+     *
+     *  \returns a copy of the transposed of the Jacobian
+     */
+    Jacobian jacobian (const LocalCoordinate &local) const
+    {
+      return jacobianTransposed(local).transposed();
+    }
+
+    /** \brief Obtain the Jacobian's inverse
+     *
+     *  The Jacobian's inverse is defined as a pseudo-inverse. If we denote
+     *  the Jacobian by \f$J(x)\f$, the following condition holds:
+     *  \f[J^{-1}(x) J(x) = I.\f]
+     */
+    JacobianInverse jacobianInverse (const LocalCoordinate &local) const
+    {
+      return jacobianInverseTransposed(local).transposed();
+    }
+
   protected:
     using Base::refElement;
 

dune/geometry/test/checkgeometry.hh

@@ -17,6 +17,39 @@
 namespace Dune
 {
 
+  namespace Impl
+  {
+
+    // Helper function to convert a matrix to a FieldMatrix
+    // using only the mv method.
+    template<class ctype, int rows, int cols, class M>
+    Dune::FieldMatrix<ctype, rows, cols> toFieldMatrix(const M& m){
+      Dune::FieldMatrix<ctype, rows, cols> result;
+      for (int j=0; j<cols; j++) {
+        FieldVector<ctype,cols> idColumn(0);
+        idColumn[j] = 1;
+        FieldVector<ctype,rows> column;
+        m.mv(idColumn,column);
+        for (int k=0; k<rows; k++)
+          result[k][j] = column[k];
+      }
+      return result;
+    }
+
+    // Helper function to check if FieldMatrix is an identity matrix
+    template<class ctype, int rows, int cols>
+    bool isIdentity(const Dune::FieldMatrix<ctype, rows, cols>& m, double tolerance){
+      if (rows != cols)
+        return false;
+      for(int i=0; i<rows; ++i)
+        for(int j=0; j<rows; ++j)
+          if (std::abs(m[i][j] - (i == j ? 1 : 0)) >= tolerance)
+            return false;
+      return true;
+    }
+
+  }
+
   /**
    * \brief Static and dynamic checks for all features of a Geometry
    *
@@ -60,6 +93,12 @@ namespace Dune
     // Matrix-like type for the return value of the jacobianInverseTransposed method
     typedef typename TestGeometry::JacobianInverseTransposed JacobianInverseTransposed;
 
+    // Matrix-like type for the return value of the jacobian method
+    typedef typename TestGeometry::Jacobian Jacobian;
+
+    // Matrix-like type for the return value of the jacobianInverse method
+    typedef typename TestGeometry::JacobianInverse JacobianInverse;
+
     const ctype tolerance = std::sqrt( std::numeric_limits< ctype >::epsilon() );
 
     ////////////////////////////////////////////////////////////////////////
@@ -161,60 +200,65 @@ namespace Dune
         pass = false;
       }
 
-      // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed'
-      // return matrices that are inverse to each other.
-      const JacobianTransposed &jt = geometry.jacobianTransposed( x );
-      const JacobianInverseTransposed &jit = geometry.jacobianInverseTransposed( x );
+      const JacobianTransposed &Jt = geometry.jacobianTransposed( x );
+      const JacobianInverseTransposed &Jit = geometry.jacobianInverseTransposed( x );
+      const Jacobian &J = geometry.jacobian( x );
+      const JacobianInverse &Ji = geometry.jacobianInverse( x );
 
       // Transform to FieldMatrix, so we can have coefficent access and other goodies
-      // We need some black magic for the transformation, because there is no
-      // official easy way yet.
-      // The following code does the transformation by multiplying jt and jit from
-      // the right by identity matrices.  That way, only the mv method is used.
-      FieldMatrix< ctype, mydim, coorddim > jtAsFieldMatrix;
-      for (int j=0; j<coorddim; j++) {
-
-        FieldVector<ctype,coorddim> idColumn(0);
-        idColumn[j] = 1;
+      auto JtAsFieldMatrix = Impl::toFieldMatrix< ctype, mydim, coorddim >(Jt);
+      auto JitAsFieldMatrix = Impl::toFieldMatrix< ctype, coorddim, mydim >(Jit);
+      auto JAsFieldMatrix = Impl::toFieldMatrix< ctype, coorddim, mydim >(J);
+      auto JiAsFieldMatrix = Impl::toFieldMatrix< ctype, mydim, coorddim >(Ji);
 
-        FieldVector<ctype,mydim> column;
-        jt.mv(idColumn,column);
-
-        for (int k=0; k<mydim; k++)
-          jtAsFieldMatrix[k][j] = column[k];
 
+      // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed'
+      // return matrices that are inverse to each other.
+      {
+        FieldMatrix< ctype, mydim, mydim > id = JtAsFieldMatrix * JitAsFieldMatrix;
+        if(not Impl::isIdentity(id, tolerance))
+        {
+          std::cerr << "Error: jacobianTransposed and jacobianInverseTransposed are not inverse to each other." << std::endl;
+          std::cout << "       id != [ ";
+          for( int j = 0; j < mydim; ++j )
+            std::cout << (j > 0 ? " | " : "") << id[ j ];
+          std::cout << " ]" << std::endl;
+          pass = false;
+        }
       }
 
-      FieldMatrix< ctype, coorddim, mydim > jitAsFieldMatrix;
-      for (int j=0; j<mydim; j++) {
-
-        FieldVector<ctype,mydim> idColumn(0);
-        idColumn[j] = 1;
-
-        FieldVector<ctype,coorddim> column;
-        jit.mv(idColumn,column);
-
-        for (int k=0; k<coorddim; k++)
-          jitAsFieldMatrix[k][j] = column[k];
-
+      // Test whether the methods 'jacobian' and 'jacobianInverse'
+      // return matrices that are inverse to each other.
+      {
+        FieldMatrix< ctype, mydim, mydim > id = JiAsFieldMatrix * JAsFieldMatrix;
+        if(not Impl::isIdentity(id, tolerance))
+        {
+          std::cerr << "Error: jacobian and jacobianInverse are not inverse to each other." << std::endl;
+          std::cout << "       id != [ ";
+          for( int j = 0; j < mydim; ++j )
+            std::cout << (j > 0 ? " | " : "") << id[ j ];
+          std::cout << " ]" << std::endl;
+          pass = false;
+        }
       }
 
-
-      FieldMatrix< ctype, mydim, mydim > id;
-      FMatrixHelp::multMatrix( jtAsFieldMatrix, jitAsFieldMatrix, id );
-      bool isId = true;
-      for( int j = 0; j < mydim; ++j )
-        for( int k = 0; k < mydim; ++k )
-          isId &= (std::abs( id[ j ][ k ] - (j == k ? 1 : 0) ) < tolerance);
-      if( !isId)
+      // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed'
+      // are the transposed of 'jacobian' and 'jacobianInverse', respectively.
       {
-        std::cerr << "Error: jacobianTransposed and jacobianInverseTransposed are not inverse to each other." << std::endl;
-        std::cout << "       id != [ ";
-        for( int j = 0; j < mydim; ++j )
-          std::cout << (j > 0 ? " | " : "") << id[ j ];
-        std::cout << " ]" << std::endl;
-        pass = false;
+        if( (JtAsFieldMatrix - JAsFieldMatrix.transposed()).infinity_norm() != 0 )
+        {
+          std::cerr << "Error: jacobian and jacobianTransposed are not transposed to each other." << std::endl;
+          pass = false;
+        }
       }
+      {
+        if( (JitAsFieldMatrix - JiAsFieldMatrix.transposed()).infinity_norm() != 0 )
+        {
+          std::cerr << "Error: jacobianInverse and jacobianInverseTransposed are not transposed to each other." << std::endl;
+          pass = false;
+        }
+      }
+
 
       // Test whether integrationElement returns something nonnegative
       if( geometry.integrationElement( x ) < 0 ) {
@@ -226,7 +270,7 @@ namespace Dune
       for( int i = 0; i < mydim; ++i )
         for( int j = 0; j < mydim; ++j )
           for( int k = 0; k < coorddim; ++k )
-            jtj[ i ][ j ] += jtAsFieldMatrix[ i ][ k ] * jtAsFieldMatrix[ j ][ k ];
+            jtj[ i ][ j ] += JtAsFieldMatrix[ i ][ k ] * JtAsFieldMatrix[ j ][ k ];
 
       if( abs( sqrt( jtj.determinant() ) - geometry.integrationElement( x ) ) > tolerance ) {
         std::cerr << "Error: integrationElement is not consistent with jacobianTransposed." << std::endl;