Add functional assembler for Dune::Functions bases

dune/fufem/assemblers/CMakeLists.txt

@@ -5,6 +5,7 @@ install(FILES
     boundaryfunctionalassembler.hh
     boundaryoperatorassembler.hh
     dunefunctionsoperatorassembler.hh
+    dunefunctionsfunctionalassembler.hh
     functionalassembler.hh
     integraloperatorassembler.hh
     intersectionoperatorassembler.hh

dune/fufem/assemblers/dunefunctionsfunctionalassembler.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUFEM_ASSEMBERS_DUNEFUNCTIONSFUNCTIONALASSEMBLER_HH
+#define DUNE_FUFEM_ASSEMBERS_DUNEFUNCTIONSFUNCTIONALASSEMBLER_HH
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/matrixindexset.hh>
+
+#include <dune/matrix-vector/axpy.hh>
+
+#include "dune/fufem/functionspacebases/functionspacebasis.hh"
+
+
+namespace Dune {
+namespace Fufem {
+
+
+//! Generic global assembler for functionals on a gridview
+template <class TrialBasis>
+class DuneFunctionsFunctionalAssembler
+{
+  using GridView = typename TrialBasis::GridView;
+
+public:
+  //! create assembler for grid
+  DuneFunctionsFunctionalAssembler(const TrialBasis& tBasis) :
+      trialBasis_(tBasis)
+  {}
+
+
+
+  template <class VectorBackend, class LocalAssembler>
+  void assembleBulkEntries(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const
+  {
+    auto trialLocalView     = trialBasis_.localView();
+    auto trialLocalIndexSet = trialBasis_.localIndexSet();
+
+    using Field = std::decay_t<decltype(vectorBackend(trialLocalIndexSet.index(0)))>;
+    using LocalVector = Dune::BlockVector<Dune::FieldVector<Field,1>>;
+
+    auto localVector = LocalVector(trialLocalView.maxSize());
+
+    for (const auto& element : elements(trialBasis_.gridView()))
+    {
+      trialLocalView.bind(element);
+      trialLocalIndexSet.bind(trialLocalView);
+
+      localAssembler(element, localVector, trialLocalView);
+
+      // Add element stiffness matrix onto the global stiffness matrix
+      for (size_t localRow=0; localRow<trialLocalIndexSet.size(); ++localRow)
+      {
+        auto row = trialLocalIndexSet.index(localRow);
+        vectorBackend(row) += localVector[localRow];
+      }
+    }
+  }
+
+  template <class VectorBackend, class LocalAssembler>
+  void assembleBulk(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const
+  {
+    vectorBackend.resize(trialBasis_);
+    vectorBackend.vector() = 0.0;
+    assembleBulkEntries(vectorBackend, std::forward<LocalAssembler>(localAssembler));
+  }
+
+
+protected:
+  const TrialBasis& trialBasis_;
+};
+
+
+/**
+ * \brief Create DuneFunctionsFunctionalAssembler
+ */
+template <class TrialBasis>
+auto duneFunctionsFunctionalAssembler(const TrialBasis& trialBasis)
+{
+  return DuneFunctionsFunctionalAssembler<TrialBasis>(trialBasis);
+}
+
+
+
+} // namespace Fufem
+} // namespace Dune
+
+
+#endif // DUNE_FUFEM_ASSEMBERS_DUNEFUNCTIONSFUNCTIONALASSEMBLER_HH
+