Merge branch 'speedup-uggrid-geometry' into 'master'

Speed up `UGGridGeometry`

See merge request !750

CHANGELOG.md

@@ -30,6 +30,14 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
 - `MultipleCodimMultipleGeomTypeMapper` is assignable.
 
+- The method `UGGridGeometry::integrationElement` has been simplified, and is faster now.
+
+- The method `UGGridGeometry::affine` has been reimplemented, and is much faster now.
+
+- For simplices, `UGGridGeometry` now caches the values of `jacobianTransposed`,
+  `jacobianInverseTransposed`, etc. This may speed up the code when these methods
+  are called many times per element.
+
 ## Python
 
 - Improve pickling support (GridViews and some GridFunction objects can now be pickled).

dune/grid/uggrid/uggridgeometry.cc

@@ -62,13 +62,8 @@ FieldVector<typename GridImp::ctype, coorddim> UGGridGeometry<mydim,coorddim,Gri
 corner(int i) const
 {
   // This geometry is a vertex
-  if (mydim==0) {
-    FieldVector<typename GridImp::ctype, coorddim> result;
-    for (size_t j=0; j<coorddim; j++)
-      // The cast is here to make the code compile even when target_ is not a node
-      result[j] = ((typename UG_NS<coorddim>::Node*)target_)->myvertex->iv.x[j];
-    return result;
-  }
+  if constexpr (mydim==0)
+    return target_->myvertex->iv.x;
 
   // ////////////////////////////////
   //  This geometry is an element
@@ -77,11 +72,8 @@ corner(int i) const
 
   i = UGGridRenumberer<mydim>::verticesDUNEtoUG(i,type());
 
-  FieldVector<typename GridImp::ctype, coorddim> result;
-  for (size_t j=0; j<coorddim; j++)
-    // The cast is here to make the code compile even when target_ is not an element
-    result[j] = UG_NS<coorddim>::Corner(((typename UG_NS<coorddim>::Element*)target_),i)->myvertex->iv.x[j];
-  return result;
+  if constexpr (mydim==coorddim)
+    return UG_NS<coorddim>::Corner(target_,i)->myvertex->iv.x;
 }
 
 template< int mydim, int coorddim, class GridImp>
@@ -133,9 +125,16 @@ integrationElement (const FieldVector<typename GridImp::ctype, mydim>& local) co
 {
   if (mydim==0)
     return 1;
-  else
-    /** \todo No need to recompute the determinant every time on a simplex */
-    return std::abs(1/jacobianInverseTransposed(local).determinant());
+
+  // Geometry is not affine --> the integrationElement cannot be cached
+  if (!affine())
+    return std::abs(jacobianTransposed(local).determinant());
+
+  // Geometry is affine --> try to use the integrationElement from the cache
+  if (!cachedIntegrationElement_)
+    cachedIntegrationElement_.emplace(std::abs(jacobianTransposed(local).determinant()));
+
+  return *cachedIntegrationElement_;
 }
 
 
@@ -143,33 +142,40 @@ template< int mydim, int coorddim, class GridImp>
 FieldMatrix<typename GridImp::ctype, coorddim,mydim> UGGridGeometry<mydim,coorddim, GridImp>::
 jacobianInverseTransposed (const FieldVector<typename GridImp::ctype, mydim>& local) const
 {
-  FieldMatrix<UGCtype,coorddim,mydim> jIT;
+  if (!affine() || !cachedJacobianInverseTransposed_)
+  {
+    cachedJacobianInverseTransposed_.emplace();
 
-  // compile array of pointers to corner coordinates
-  // coorddim*coorddim is an upper bound for the number of vertices
-  UGCtype* cornerCoords[coorddim*coorddim];
-  int n = UG_NS<coorddim>::Corner_Coordinates(target_, cornerCoords);
+    // compile array of pointers to corner coordinates
+    // coorddim*coorddim is an upper bound for the number of vertices
+    UGCtype* cornerCoords[coorddim*coorddim];
+    const int n = UG_NS<coorddim>::Corner_Coordinates(target_, cornerCoords);
 
-  // compute the transformation onto the reference element (or vice versa?)
-  UG_NS<coorddim>::Transformation(n, cornerCoords, local, jIT);
+    // compute the transformation onto the reference element (or vice versa?)
+    UG_NS<coorddim>::JacobianInverseTransformation(n, cornerCoords, local, *cachedJacobianInverseTransposed_);
+  }
 
-  return jIT;
+  return *cachedJacobianInverseTransposed_;
 }
+
 template< int mydim, int coorddim, class GridImp>
 FieldMatrix<typename GridImp::ctype, mydim,coorddim> UGGridGeometry<mydim,coorddim, GridImp>::
 jacobianTransposed (const FieldVector<typename GridImp::ctype, mydim>& local) const
 {
-  FieldMatrix<UGCtype,mydim,coorddim> jac;
+  if (!affine() || !cachedJacobianTransposed_)
+  {
+    cachedJacobianTransposed_.emplace();
 
-  // compile array of pointers to corner coordinates
-  // coorddim*coorddim is an upper bound for the number of vertices
-  UGCtype* cornerCoords[coorddim*coorddim];
-  int n = UG_NS<coorddim>::Corner_Coordinates(target_, cornerCoords);
+    // compile array of pointers to corner coordinates
+    // coorddim*coorddim is an upper bound for the number of vertices
+    UGCtype* cornerCoords[coorddim*coorddim];
+    const int n = UG_NS<coorddim>::Corner_Coordinates(target_, cornerCoords);
 
-  // compute the transformation onto the reference element (or vice versa?)
-  UG_NS<coorddim>::JacobianTransformation(n, cornerCoords, local, jac);
+    // compute the transformation onto the reference element (or vice versa?)
+    UG_NS<coorddim>::JacobianTransformation(n, cornerCoords, local, *cachedJacobianTransposed_);
+  }
 
-  return jac;
+  return *cachedJacobianTransposed_;
 }
 
 

dune/grid/uggrid/uggridgeometry.hh

@@ -49,7 +49,15 @@ namespace Dune {
     GeometryType type () const;
 
     //! returns true if type is simplex, false otherwise
-    bool affine() const { return type().isSimplex(); }
+    bool affine() const
+    {
+      if constexpr (mydim==0 || mydim==1)  // Vertices and edges are always affine
+        return true;
+      else if constexpr (mydim==2)
+        return UG_NS<coorddim>::Tag(target_)==UG::D2::TRIANGLE;
+      else
+        return UG_NS<coorddim>::Tag(target_)==UG::D3::TETRAHEDRON;
+    }
 
     //! return the number of corners of this element.
     int corners () const {
@@ -113,12 +121,21 @@ namespace Dune {
     void setToTarget(typename UG_NS<coorddim>::template Entity<coorddim-mydim>::T* target)
     {
       target_ = target;
+
+      cachedIntegrationElement_.reset();
+      cachedJacobianTransposed_.reset();
+      cachedJacobianInverseTransposed_.reset();
     }
 
     // in element mode this points to the element we map to
     // in coord_mode this is the element whose reference element is mapped into the father's one
     typename UG_NS<coorddim>::template Entity<coorddim-mydim>::T* target_;
 
+    // If the element is affine, then the geometry Jacobian is constant, and only needs to be
+    // computed once per element.  Therefore, keep them in a cache.
+    mutable std::optional<UGCtype> cachedIntegrationElement_;
+    mutable std::optional<FieldMatrix<UGCtype,mydim,coorddim> > cachedJacobianTransposed_;
+    mutable std::optional<FieldMatrix<UGCtype,coorddim,mydim> > cachedJacobianInverseTransposed_;
   };
 
 

dune/grid/uggrid/ugwrapper.hh

@@ -926,8 +926,8 @@ namespace Dune {
      * \return The return type is int because the macro INVERSE_TRANSFORMATION
      *  returns 1 on failure.
      */
-    static int Transformation(int n, double** x,
-                              const FieldVector<double, UG_DIM>& local, FieldMatrix<double,UG_DIM,UG_DIM>& mat) {
+    static int JacobianInverseTransformation(int n, double** x,
+                                             const FieldVector<double, UG_DIM>& local, FieldMatrix<double,UG_DIM,UG_DIM>& mat) {
       using UG_NAMESPACE ::DOUBLE_VECTOR;
       using UG::DOUBLE;
       double det;
@@ -939,8 +939,8 @@ namespace Dune {
     }
 
     /** \brief Dummy method for vertices */
-    static int Transformation(int n, double** x,
-                              const FieldVector<double, 0>& local, FieldMatrix<double,UG_DIM,0>& mat) {
+    static int JacobianInverseTransformation(int n, double** x,
+                                             const FieldVector<double, 0>& local, FieldMatrix<double,UG_DIM,0>& mat) {
       return 0;
     }
 

dune/python/test/test-ug.py

@@ -15,7 +15,7 @@ except ImportError as exception:
 base = "../../../doc/grids/gmsh/"
 reader = (dune.grid.reader.gmsh, base+"circle2ndorder.msh")
 grid = dune.grid.ugGrid(reader, dimgrid=2)
-# we either have UG, ot this test should not be run at all 
+# we either have UG, or this test should not be run at all
 assert(grid)
 
 exact = 6.27593206157460
@@ -44,7 +44,7 @@ for l in range(5):
 if math.isclose(eoc,2.0,abs_tol=0.02):
     print("Geometry approximation: 2d order convergence")
 else:
-    print("Geometry approximation does not show expected 2d order convergence")
+    print("Geometry approximation does not show expected 2nd order convergence")
     sys.exit(-1)
 
 sys.exit(0)