Remove tracking of processed entries in assembleGlobalBasisTransferMatrix

To avoid duplicate work, `assembleGlobalBasisTransferMatrix()` used
a container to store all processed fine basis functions. The used
container is a `std::unordered_set<MultiIndex>`.

Unfortunately this container is really slow. With the improved
implemention it turned that checking, filling, clearning this
container is in fact more expensive and the assembly is faster
if the checks are removed. This is based on the test cases
of the corresponding test, that checks various combinations
of bases on 2d and 3d grids.

Using a suitable bit vector to store flags would indeed be faster.
However, it is very hard to derive a suitable container type
generically based on the current basis interface. (This will
hopefully iprove soon.)

dune/fufem/assemblers/basisinterpolationmatrixassembler.hh

@@ -8,6 +8,7 @@
 #include <unordered_set>
 
 #include <dune/common/bitsetvector.hh>
+#include <dune/common/timer.hh>
 
 #include <dune/istl/scaledidmatrix.hh>
 
@@ -128,7 +129,7 @@ namespace Impl{
  *       some assembly time.
  */
 
-template <class CLV, class FLV, class G, class B, class M, class I>
+template <class CLV, class FLV, class G, class M, class I>
 class LeafNodeInterpolationVisitor
   : public Dune::TypeTree::TreePairVisitor
   , public Dune::TypeTree::DynamicTraversal
@@ -138,7 +139,6 @@ public:
   using CoarseLocalView = CLV;
   using FineLocalView = FLV;
   using Geometry = G;
-  using BitSet = B;
   using Matrix = M;
   using IndexSet = I;
 
@@ -147,18 +147,16 @@ public:
    * \param [in] clv (coarseLocalView) coarse element information and localBasis
    * \param [in] flv (fineLocalView) fine element information and localBasis
    * \param [in] g (geometry) mapping fine -> coarse
-   * \param [in] p (processed) set of already processed dof's
    * \param [in] t (tolerance) threshold for considering an interpolated value as zero
    * \param [out] m (matrix) resulting tranferoperator matrix (coarse basis) -> (fine basis) wrapped in ISTLBackend
    * \param [out] i (indices) index set of the matrix
    * \param [in] e (exportIndices) boolean switch: true  = add indices to the index set
    *                                               false = increase the matrix entries
    */
-  LeafNodeInterpolationVisitor(const CLV& clv, const FLV& flv, const G& g, const B& p, double t, M& m, I& i, bool e)
+  LeafNodeInterpolationVisitor(const CLV& clv, const FLV& flv, const G& g, double t, M& m, I& i, bool e)
   : coarseLocalView_(clv)
   , fineLocalView_(flv)
   , geometry_(g)
-  , processed_(p)
   , tolerance_(t)
   , matrix_(m)
   , indices_(i)
@@ -233,9 +231,6 @@ public:
         // get the matrix row index
         const auto& globalFineIndex = fineLocalView_.index( fineNode.localIndex( i ) );
 
-        if (processed_.find(globalFineIndex) != processed_.end() )
-          continue;
-
         // Either insert the index or add the value.
         if (exportIndices_)
           indices_.insertEntry(globalFineIndex, globalCoarseIndex);
@@ -250,7 +245,6 @@ private:
   const CoarseLocalView& coarseLocalView_;
   const FineLocalView& fineLocalView_;
   const Geometry& geometry_;
-  const BitSet& processed_;
   double tolerance_;
   Matrix& matrix_;
   IndexSet& indices_;
@@ -277,6 +271,7 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
                                               const FineBasis& fineBasis,
                                               const double tolerance = 1e-12)
 {
+  auto timer = Dune::Timer();
   const auto& coarseGridView = coarseBasis.gridView();
   const auto& fineGridView   = fineBasis.gridView();
 
@@ -296,9 +291,6 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
   // Determine the indices present in the 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;
 
@@ -333,23 +325,16 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
 
     // visit all children of the coarse node and interpolate with the fine basis functions
     ::Impl::LeafNodeInterpolationVisitor visitor( coarseLocalView, fineLocalView, geometryInFathers,
-                                                processed, tolerance, matrixBackend, indices,
+                                                tolerance, matrixBackend, indices,
                                                 true /* = set indexSet */ );
 
     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++)
-      processed.insert(fineLocalView.index(i));
   }
 
   // apply the index set to the matrix
   indices.setupMatrix();
   matrix = 0;
 
-  // reset all dof's
-  processed.clear();
-
   //////////////////////
   // Now set the values
   //////////////////////
@@ -380,15 +365,12 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
 
     // visit all children of the coarse node and interpolate with the fine basis functions
     ::Impl::LeafNodeInterpolationVisitor visitor( coarseLocalView, fineLocalView, geometryInFathers,
-                                                processed, tolerance, matrixBackend, indices,
+                                                tolerance, matrixBackend, indices,
                                                 false /* = increase matrix entries */ );
 
     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++)
-      processed.insert(fineLocalView.index(i));
   }
+  std::cout << "Assembly took " << timer.elapsed() << std::endl;
 }