Improve performance of assembleGlobalBasisTransferMatrix

The old approach used a `std::map` and compared coordinates
in the inner most loop to avoid duplicate evaluations
of `LocalBasis`. Instead the new approach determines the
interpolation points first and caches the values in advaced
based on the evaluation order. This improved the performance
significantly.

dune/fufem/assemblers/basisinterpolationmatrixassembler.hh

@@ -183,16 +183,39 @@ public:
     using CoarseFEType = std::decay_t<decltype(coarseFiniteElement)>;
     using FunctionTraits = typename CoarseFEType::Traits::LocalBasisType::Traits;
     using LocalBasisWrapper = LocalBasisComponentWrapper<typename CoarseFEType::Traits::LocalBasisType, FunctionTraits>;
-
-    LocalBasisWrapper basisFctj(coarseLocalBasis,0);
+    using CoarseRange = typename std::decay_t<decltype(coarseLocalBasis)>::Traits::RangeType;
+
+    // We first compute the values at all interpolation points.
+    // To this end a dummy call to interpolate() is used to deduce
+    // the interpolation points and evaluate the coarse bases.
+    // Later the, we can simply return the pre-computed values
+    // based on the evaluation order.
+    // This requires that interpolate() does the evaluation at
+    // a deterministic point set in a deterministic order and
+    // thus rules out e.g. adaptive quadrature.
+    // If the interpolation is given by simple point evaluations,
+    // the cache computed here is essentially the local interpolation
+    // matrix. However, it is also fine, if interpolate() does
+    // additional computations based on the values, like, e.g.
+    // a Piola-transformation.
+    auto coarseValues = std::vector<std::vector<CoarseRange>>(fineLocalBasis.size());
+    auto k = 0;
+    auto trackingCoarseFunction = [&](const auto& x) {
+      coarseValues.resize(k+1);
+      coarseLocalBasis.evaluateFunction(geometry_.global(x), coarseValues[k]);
+      return coarseValues[k++][0];
+    };
+    fineFiniteElement.localInterpolation().interpolate(trackingCoarseFunction, coefficients);
 
     for (size_t j=0; j<coarseLocalBasis.size(); j++)
     {
-      // wrap each local basis function as a local function.
-      basisFctj.setIndex(j);
-
-      // transform the local fine function to a local coarse function
-      const auto baseFctj = OtherGeometryWrapper(basisFctj, geometry_);
+      // baseFctj behaves like the j-th coarse bases function if
+      // interpolate() evaluates it at the same points in the same
+      // order as above.
+      auto k = 0;
+      auto baseFctj = [&](const auto& x) {
+        return coarseValues[k++][j];
+      };
 
       // Interpolate j^th base function by the fine basis
       fineFiniteElement.localInterpolation().interpolate(baseFctj, coefficients);
@@ -213,16 +236,11 @@ public:
         if (processed_.find(globalFineIndex) != processed_.end() )
           continue;
 
-        if ( exportIndices_ )
-        {
-          // add only the index of the top level
+        // Either insert the index or add the value.
+        if (exportIndices_)
           indices_.insertEntry(globalFineIndex, globalCoarseIndex);
-        }
         else
-        {
-          // use ISTLBackend wrapper with ()-operator
-          matrix_(globalFineIndex,globalCoarseIndex) += coefficients[i];
-        }
+          matrix_(globalFineIndex,globalCoarseIndex) = coefficients[i];
       }
     }
   }
@@ -281,14 +299,19 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
   // Only handle every fine dof once
   std::unordered_set<typename FineLocalView::MultiIndex> processed;
 
+  using GeometryInFather = std::decay_t<decltype(elements(fineGridView).begin()->geometryInFather())>;
+  std::vector<GeometryInFather> geometryInFathersVector;
+
   // loop over all fine grid elements
   for ( const auto& fineElement : elements(fineGridView) )
   {
+    fineLocalView.bind(fineElement);
+
     // find a coarse element that is the parent of the fine element
     // start in the fine grid and track the geometry mapping
     auto fE = fineElement;
     // collect all geometryInFather's on the way
-    std::vector<decltype(fineElement.geometryInFather())> geometryInFathersVector;
+    geometryInFathersVector.clear();
     while ( not coarseGridView.contains(fE) and fE.hasFather() )
     {
       // add the geometry to the container
@@ -307,7 +330,6 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
     const auto geometryInFathers = GeometryInFathers(geometryInFathersVector);
 
     coarseLocalView.bind(coarseElement);
-    fineLocalView.bind(fineElement);
 
     // visit all children of the coarse node and interpolate with the fine basis functions
     ::Impl::LeafNodeInterpolationVisitor visitor( coarseLocalView, fineLocalView, geometryInFathers,
@@ -316,13 +338,9 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
 
     Dune::TypeTree::applyToTreePair(coarseLocalView.tree(),fineLocalView.tree(),visitor);
 
-
     // now the all dof's of this fine element to 'processed'
     for (size_t i=0; i<fineLocalView.size(); i++)
-    {
-      const auto& globalFineIndex = fineLocalView.index( i );
-      processed.insert(globalFineIndex);
-    }
+      processed.insert(fineLocalView.index(i));
   }
 
   // apply the index set to the matrix
@@ -338,11 +356,13 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
 
   for ( const auto& fineElement : elements(fineGridView) )
   {
+    fineLocalView.bind(fineElement);
+
     // find a coarse element that is the parent of the fine element
     // start in the fine grid and track the geometry mapping
     auto fE = fineElement;
     // collect all geometryInFather's on the way
-    std::vector<decltype(fineElement.geometryInFather())> geometryInFathersVector;
+    geometryInFathersVector.clear();
     while ( not coarseGridView.contains(fE) and fE.hasFather() )
     {
       // add the geometry to the container
@@ -357,7 +377,6 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
     const auto geometryInFathers = GeometryInFathers(geometryInFathersVector);
 
     coarseLocalView.bind(coarseElement);
-    fineLocalView.bind(fineElement);
 
     // visit all children of the coarse node and interpolate with the fine basis functions
     ::Impl::LeafNodeInterpolationVisitor visitor( coarseLocalView, fineLocalView, geometryInFathers,
@@ -368,10 +387,7 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
 
     // now the all dof's of this fine element to 'processed'
     for (size_t i=0; i<fineLocalView.size(); i++)
-    {
-      const auto& globalFineIndex = fineLocalView.index( i );
-      processed.insert(globalFineIndex);
-    }
+      processed.insert(fineLocalView.index(i));
   }
 }