Draft: Use a container with flexible key type as cache for quadrature rules

Summary

For a thread-safe cache of quadrature rules, currently a vector of vectors of vectors with several std::call_once flags is used. This makes reading the code very complicated and does not necessarily lead to any better performance than the proposed solution in this MR. Replace the nested vector with a simple std::map where the key allows to identify all quadrature rules and additionally makes it possible to extend quadrature rules properties in the future.

Discussion

1. Instead of a std::map container, any other associative container could be used. I have chose a simple map since comparison of three integers by < is very simple to implement. The alternative std::unordered_map would require a hash function. While this is possible, it requires more code. I do not expect a measurable performance difference.

2. The cache is implemented as thread_local. This allows to omit more complicated locking mechanisms and should be thread-safe by definition. It is, however, unclear to me, whether this has negative performance consequences it threads are closed and reopened multiple times, i.e., if this makes it necessary to re-"compute"/re-fill the cache or whether this preserves the previously cached values. An alternative implementation uses the static keyword:

// Container to store the cached values
static std::map<QuadratureKey, QuadratureRule> quadCache;

// mutex used to access the data in the container, necessary since
// access emplace is read-write.
using mutex_type = std::shared_timed_mutex;
static mutex_type access_mutex;

// define the quadrature key to be used to access the quadrature rule
const QuadratureKey key{t.id(),p,qt};

// first try to lock for read-only, if a quadratur rule for key is found, return it,
// if not, obtain a unique_lock to insert a new rule.
std::shared_lock<mutex_type> read_lock(access_mutex);
auto it = quadCache.find(key);
if (it != quadCache.end())
  return it->second;
else {
  read_lock.unlock();
  QuadratureRule rule = QuadratureRuleFactory<ctype,dim>::rule(t,p,qt);
  std::unique_lock<mutex_type> write_lock(access_mutex);
  auto new_it = quadCache.emplace(key, std::move(rule));
  return new_it.first->second;
}

dune/geometry/quadraturerules.hh

@@ -7,7 +7,7 @@
 #include <algorithm>
 #include <iostream>
 #include <limits>
-#include <mutex>
+#include <map>
 #include <utility>
 #include <vector>
 
@@ -195,98 +195,44 @@ namespace Dune {
       \ingroup Quadrature
    */
   template<typename ctype, int dim>
-  class QuadratureRules {
-
-    /** \brief Internal short-hand notation for the type of quadrature rules this container contains */
-    typedef Dune::QuadratureRule<ctype, dim> QuadratureRule;
-    //! \brief a quadrature rule (for each quadrature order, geometry type,
-    //!        and quadrature type)
-    static void initQuadratureRule(QuadratureRule *qr, QuadratureType::Enum qt,
-                                   const GeometryType &t, int p)
-    {
-      *qr = QuadratureRuleFactory<ctype,dim>::rule(t,p,qt);
-    }
-
-    typedef std::vector<std::pair<std::once_flag, QuadratureRule> >
-      QuadratureOrderVector; // indexed by quadrature order
-    //! \brief initialize the vector indexed by the quadrature order (for each
-    //!        geometry type and quadrature type)
-    static void initQuadratureOrderVector(QuadratureOrderVector *qov,
-                                          QuadratureType::Enum qt,
-                                          const GeometryType &t)
-    {
-      if(dim == 0)
-        // we only need one quadrature rule for points, not maxint
-        *qov = QuadratureOrderVector(1);
-      else
-        *qov = QuadratureOrderVector(QuadratureRuleFactory<ctype,dim>::maxOrder(t,qt)+1);
-    }
-
-    typedef std::vector<std::pair<std::once_flag, QuadratureOrderVector> >
-      GeometryTypeVector; // indexed by geometry type
-    //! \brief initialize the vector indexed by the geometry type (for each
-    //!        quadrature type)
-    static void initGeometryTypeVector(GeometryTypeVector *gtv)
-    {
-      *gtv = GeometryTypeVector(LocalGeometryTypeIndex::size(dim));
-    }
+  class QuadratureRules
+  {
+  public:
+    //! type of quadrature rule provided by this cache
+    using QuadratureRule = Dune::QuadratureRule<ctype, dim>;
 
-    //! real rule creator
-    DUNE_EXPORT const QuadratureRule& _rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
+    //! identifier for quadrature rules, used in the cache container
+    struct QuadratureKey
     {
-      assert(t.dim()==dim);
-
-      DUNE_ASSERT_CALL_ONCE();
+      unsigned int id;          // topologyId
+      int p;                    // order
+      QuadratureType::Enum qt;  // quadrature type
 
-      static std::vector<std::pair< // indexed by quadrature type
-        std::once_flag,
-        GeometryTypeVector
-        > > quadratureCache(QuadratureType::size);
-
-      auto & quadratureTypeLevel = quadratureCache[qt];
-      std::call_once(quadratureTypeLevel.first, initGeometryTypeVector,
-                     &quadratureTypeLevel.second);
-
-      auto & geometryTypeLevel =
-        quadratureTypeLevel.second[LocalGeometryTypeIndex::index(t)];
-      std::call_once(geometryTypeLevel.first, initQuadratureOrderVector,
-                     &geometryTypeLevel.second, qt, t);
-
-      // we only have one quadrature rule for points
-      auto & quadratureOrderLevel = geometryTypeLevel.second[dim == 0 ? 0 : p];
-      std::call_once(quadratureOrderLevel.first, initQuadratureRule,
-                     &quadratureOrderLevel.second, qt, t, p);
+      friend bool operator< (const QuadratureKey& lhs, const QuadratureKey& rhs)
+      {
+        return std::tie(lhs.id,lhs.p,lhs.qt) < std::tie(rhs.id,rhs.p,rhs.qt);
+      }
+    };
 
-      return quadratureOrderLevel.second;
-    }
-    //! singleton provider
-    DUNE_EXPORT static QuadratureRules& instance()
-    {
-      static QuadratureRules instance;
-      return instance;
-    }
-    //! private constructor
-    QuadratureRules () {}
   public:
     //! maximum quadrature order for given geometry type and quadrature type
-    static unsigned
-    maxOrder(const GeometryType& t,
-             QuadratureType::Enum qt=QuadratureType::GaussLegendre)
+    static unsigned maxOrder (const GeometryType& t,
+                              QuadratureType::Enum qt = QuadratureType::GaussLegendre)
     {
       return QuadratureRuleFactory<ctype,dim>::maxOrder(t,qt);
     }
 
     //! select the appropriate QuadratureRule for GeometryType t and order p
-    static const QuadratureRule& rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
+    static const QuadratureRule& rule (const GeometryType& t, int p,
+                                       QuadratureType::Enum qt = QuadratureType::GaussLegendre)
     {
-      return instance()._rule(t,p,qt);
-    }
+      // Container to store the cached quadrature rules
+      thread_local std::map<QuadratureKey, QuadratureRule> quadCache;
 
-    //! @copydoc rule
-    static const QuadratureRule& rule(const GeometryType::BasicType t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
-    {
-      GeometryType gt(t,dim);
-      return instance()._rule(gt,p,qt);
+      auto [it,inserted] = quadCache.try_emplace(QuadratureKey{t.id(),p,qt});
+      if (inserted)
+        it->second = QuadratureRuleFactory<ctype,dim>::rule(t,p,qt);
+      return it->second;
     }
   };