Cache a few quantities if the element is a simplex

In particular, we cache integrationElement, jacobianTransposed
and jacobianInverseTransposed.  This leads to speedup
for high-order methods on simplex grids.

CHANGELOG.md

@@ -32,6 +32,10 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
 - 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

@@ -125,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 */
+
+  // 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_;
 }
 
 
@@ -135,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>::JacobianInverseTransformation(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

@@ -121,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_;
   };