Merge branch 'feature/operators_in_dofmatrix' into 'develop'

Feature/operators in dofmatrix

See merge request spraetor/dune-amdis!38

examples/init/stokes.dat.2d

@@ -6,7 +6,7 @@ stokesMesh->num cells: 4 4
 
 stokes->mesh: stokesMesh
 
-stokes->solver->name: bicgstab_ell
+stokes->solver->name: direct
 stokes->solver->ell: 3
 stokes->solver->max iteration: 1000
 stokes->solver->relative tolerance: 1e-8

src/amdis/Assembler.hpp

 #pragma once
 
-#include <memory>
-#include <tuple>
-
-#include <amdis/DirichletBC.hpp>
-#include <amdis/LocalAssemblerList.hpp>
-#include <amdis/common/Mpl.hpp>
-
 namespace AMDiS
 {
-  template <class Traits>
-  class Assembler
+  template <class GridView, class Element, class Operators, class ElementAssembler>
+  void assembleOperators(
+      GridView const& gridView,
+      Element const& element,
+      Operators& operators,
+      ElementAssembler const& localAssembler)
   {
-    using GlobalBasis = typename Traits::GlobalBasis;
-
-    /// The grid view the global FE basis lives on
-    using GridView = typename GlobalBasis::GridView;
-
-  public:
-
-    /// Constructor, stores a shared-pointer to the feSpaces
-    Assembler(GlobalBasis& globalBasis,
-              MatrixOperators<GlobalBasis>& matrixOperators,
-              VectorOperators<GlobalBasis>& rhsOperators,
-              Constraints<GlobalBasis>& constraints)
-      : globalBasis_(globalBasis)
-      , matrixOperators_(matrixOperators)
-      , rhsOperators_(rhsOperators)
-      , constraints_(constraints)
-    {}
-
-    /// Assemble the linear system
-    template <class SystemMatrixType, class SystemVectorType>
-    void assemble(
-        SystemMatrixType& matrix,
-        SystemVectorType& solution,
-        SystemVectorType& rhs,
-        bool asmMatrix, bool asmVector);
-
-
-  private:
-
-    /// Sets the system to zero and initializes all operators and boundary conditions
-    template <class SystemMatrixType, class SystemVectorType>
-    void initMatrixVector(
-        SystemMatrixType& matrix,
-        SystemVectorType& solution,
-        SystemVectorType& rhs,
-        bool asmMatrix, bool asmVector) const;
-
-    /// Assemble operators on an element, by passing the element/intersection to
-    /// `elementAssembler` functor.
-    template <class Element, class Operators, class ElementAssembler>
-    void assembleElementOperators(
-        Element const& element,
-        Operators& operators,
-        ElementAssembler const& elementAssembler) const;
-
-    /// Finish insertion into the matrix and assembles boundary conditions
-    /// Return the number of nonzeros assembled into the matrix
-    template <class SystemMatrixType, class SystemVectorType>
-    std::size_t finishMatrixVector(
-        SystemMatrixType& matrix,
-        SystemVectorType& solution,
-        SystemVectorType& rhs,
-        bool asmMatrix, bool asmVector) const;
-
-
-  private:
-
-    GlobalBasis& globalBasis_;
-    MatrixOperators<GlobalBasis>& matrixOperators_;
-    VectorOperators<GlobalBasis>& rhsOperators_;
-    Constraints<GlobalBasis>& constraints_;
-
-    //TODO: add caching of localBases
-  };
+    // assemble element operators
+    localAssembler(element, operators.element);
+
+    // assemble intersection operators
+    if (!operators.intersection.empty()
+        || (!operators.boundary.empty() && element.hasBoundaryIntersections()))
+    {
+      for (auto const& intersection : intersections(gridView, element)) {
+        if (intersection.boundary())
+          localAssembler(intersection, operators.boundary);
+        else
+          localAssembler(intersection, operators.intersection);
+      }
+    }
+  }
 
 } // end namespace AMDiS
-
-#include "Assembler.inc.hpp"

src/amdis/Assembler.inc.hpp

-#pragma once
-
-#include <dune/common/dynmatrix.hh>
-#include <dune/common/dynvector.hh>
-#include <dune/functions/functionspacebases/subspacebasis.hh>
-#include <amdis/utility/TreePath.hpp>
-#include <amdis/utility/Visitor.hpp>
-
-#include <amdis/common/Math.hpp>
-
-namespace AMDiS {
-
-template <class Traits>
-  template <class SystemMatrixType, class SystemVectorType>
-void Assembler<Traits>::assemble(
-    SystemMatrixType& matrix,
-    SystemVectorType& solution,
-    SystemVectorType& rhs,
-    bool asmMatrix, bool asmVector)
-{
-  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
-  initMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
-
-  auto localView = globalBasis_.localView();
-
-  // 2. create a local matrix and vector
-  std::size_t localSize = localView.maxSize();
-
-  Dune::DynamicMatrix<double> elementMatrix(localSize, localSize);
-  Dune::DynamicVector<double> elementVector(localSize);
-
-  // 3. traverse grid and assemble operators on the elements
-  for (auto const& element : elements(globalBasis_.gridView()))
-  {
-    elementMatrix = 0;
-    elementVector = 0;
-
-    localView.bind(element);
-    auto geometry = element.geometry();
-
-    // traverse type-tree of global-basis
-    forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
-    {
-      auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
-      auto rowLocalView = rowBasis.localView();
-      rowLocalView.bind(element);
-
-      auto& rhsOp = rhsOperators_[rowNode];
-      if (rhsOp.doAssemble(asmVector) && !rhsOp.empty()) {
-        rhsOp.bind(element, geometry);
-
-        auto vecAssembler = [&](auto const& context, auto& operator_list) {
-          for (auto scaled : operator_list)
-            scaled.op->assemble(context, rowLocalView.tree(), elementVector);
-        };
-
-        this->assembleElementOperators(element, rhsOp, vecAssembler);
-      }
-
-      forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
-      {
-        auto& matOp = matrixOperators_[rowNode][colNode];
-        if (matOp.doAssemble(asmMatrix) && !matOp.empty()) {
-          auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
-          auto colLocalView = colBasis.localView();
-          colLocalView.bind(element);
-
-          matOp.bind(element, geometry);
-
-          auto matAssembler = [&](auto const& context, auto& operator_list) {
-            for (auto scaled : operator_list)
-              scaled.op->assemble(context, rowLocalView.tree(), colLocalView.tree(), elementMatrix);
-          };
-
-          this->assembleElementOperators(element, matOp, matAssembler);
-        }
-      });
-    });
-
-    // add element-matrix to system-matrix and element-vector to rhs
-    matrix.insert(localView, localView, elementMatrix);
-    rhs.insert(localView, elementVector);
-
-    // unbind all operators
-    forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
-      rhsOperators_[rowNode].unbind();
-      forEachNode_(localView.tree(), [&,this](auto const& colNode, auto&&) {
-        matrixOperators_[rowNode][colNode].unbind();
-      });
-    });
-
-    localView.unbind();
-  }
-
-  // 4. finish matrix insertion and apply dirichlet boundary conditions
-  std::size_t nnz = finishMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
-
-  msg("fill-in of assembled matrix: {}", nnz);
-}
-
-
-template <class Traits>
-  template <class Element, class Operators, class ElementAssembler>
-void Assembler<Traits>::assembleElementOperators(
-    Element const& element,
-    Operators& operators,
-    ElementAssembler const& localAssembler) const
-{
-  // assemble element operators
-  localAssembler(element, operators.element);
-
-  // assemble intersection operators
-  if (!operators.intersection.empty()
-      || (!operators.boundary.empty() && element.hasBoundaryIntersections()))
-  {
-    for (auto const& intersection : intersections(globalBasis_.gridView(), element)) {
-      if (intersection.boundary())
-        localAssembler(intersection, operators.boundary);
-      else
-        localAssembler(intersection, operators.intersection);
-    }
-  }
-}
-
-
-template <class Traits>
-  template <class SystemMatrixType, class SystemVectorType>
-void Assembler<Traits>::initMatrixVector(
-    SystemMatrixType& matrix,
-    SystemVectorType& solution,
-    SystemVectorType& rhs,
-    bool asmMatrix, bool asmVector) const
-{
-  matrix.init(asmMatrix);
-  solution.compress();
-  rhs.compress();
-  if (asmVector)
-    rhs = 0;
-
-  auto localView = globalBasis_.localView();
-  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
-  {
-#ifdef HAVE_EXTENDED_DUNE_FUNCTIONS
-    if (rowNode.isLeaf)
-      msg("{} DOFs for Basis[{}]",globalBasis_.dimension(rowTreePath), to_string(rowTreePath));
-#endif
-
-    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
-
-    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
-    {
-      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
-
-      for (auto bc : constraints_[rowNode][colNode])
-        bc->init(matrix, solution, rhs, rowBasis, colBasis);
-    });
-  });
-  msg("{} total DOFs", globalBasis_.dimension());
-}
-
-
-template <class Traits>
-  template <class SystemMatrixType, class SystemVectorType>
-std::size_t Assembler<Traits>::finishMatrixVector(
-    SystemMatrixType& matrix,
-    SystemVectorType& solution,
-    SystemVectorType& rhs,
-    bool asmMatrix, bool asmVector) const
-{
-  matrix.finish();
-
-  auto localView = globalBasis_.localView();
-  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
-  {
-    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
-    auto& rhsOp = rhsOperators_[rowNode];
-    if (rhsOp.doAssemble(asmVector))
-      rhsOp.assembled = true;
-
-    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
-    {
-      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
-      auto& matOp = matrixOperators_[rowNode][colNode];
-      if (matOp.doAssemble(asmMatrix))
-        matOp.assembled = true;
-
-      // finish boundary condition
-      for (auto bc : constraints_[rowNode][colNode])
-        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
-    });
-  });
-
-  return matrix.nnz();
-}
-
-} // end namespace AMDiS

src/amdis/Boundary.hpp

+#pragma once
+
+namespace AMDiS
+{
+  struct BoundaryType { int b; };
+
+} // end namespace AMDiS

src/amdis/CMakeLists.txt

@@ -20,12 +20,11 @@ install(FILES
     AdaptStationary.hpp
     AMDiS.hpp
     Assembler.hpp
-    Assembler.inc.hpp
+    Boundary.hpp
     ContextGeometry.hpp
     CreatorInterface.hpp
     CreatorMap.hpp
     DirichletBC.hpp
-    DirichletBC.inc.hpp
     FileWriter.hpp
     Flag.hpp
     GridFunctionOperator.hpp

src/amdis/DirichletBC.hpp

@@ -6,18 +6,20 @@
 #include <vector>
 
 #include <dune/common/hybridutilities.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
 #include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
 
+#include <amdis/Boundary.hpp>
 #include <amdis/Output.hpp>
 #include <amdis/common/Concepts.hpp>
 #include <amdis/common/ValueCategory.hpp>
 #include <amdis/utility/RangeType.hpp>
 #include <amdis/utility/TreeData.hpp>
+#include <amdis/utility/Visitor.hpp>
 
 namespace AMDiS
 {
-  struct BoundaryType { int b; };
-
   /// Implements a boundary condition of Dirichlet-type.
   /**
    * By calling the methods \ref init() and \ref finish before and after
@@ -49,66 +51,61 @@ namespace AMDiS
     /// Prepare the matrix, rhs and solution block for assembling dirichlet
     /// boundary conditions, e.g. collect the corresponding indices of boundary
     /// DOFS for later elimination in \ref finish.
-    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
-    void init(Matrix& /*matrix*/, VectorX& /*rhs*/, VectorB& /*solution*/,
-              RowBasis const& rowBasis,
-              ColBasis const& /*colBasis*/)
+    template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
+    void init(Matrix& /*matrix*/, VectorX& /*solution*/, VectorB& rhs,
+              RowNode const& rowNode, ColNode const&)
     {
       if (!initialized_) {
-        using RowNode = typename RowBasis::LocalView::Tree;
-        initImpl(rowBasis, typename RowNode::NodeTag{});
+        forEachLeafNode_(rowNode, [&](auto const& /*node*/, auto const& tp) {
+          auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), cat(rowNode.treePath(),tp));
+          this->initImpl(subBasis);
+        });
         initialized_ = true;
       }
     }
 
 
-
     /// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row
     /// to the matrix and optionally delete the corresponding matrix-column.
-    /** If DBC is set for matrix block {R,C}, then \p apply_row means that
-      * the block-row R is assembled, and \p apply_col means that the block-col C
-      * is assembled. The \p matrix, \p rhs, and \p solution correspond to the
-      * currently visited blocks of system-matrix and system-vector.
-      **/
-    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
-    void finish(Matrix& matrix,
-                VectorX& rhs,
-                VectorB& solution,
-                RowBasis const& rowBasis,
-                ColBasis const& colBasis)
+    template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
+    void finish(Matrix& matrix, VectorX& solution, VectorB& rhs,
+                RowNode const& rowNode, ColNode const& colNode)
     {
       test_exit_dbg(initialized_, "Boundary condition not initialized!");
       auto columns = matrix.applyDirichletBC(dirichletNodes_);
 
-      using Dune::Functions::interpolate;
-      Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename RowBasis::LocalView::Tree>, Range>{},
+      Dune::Hybrid::ifElse(std::is_same<RangeType_t<RowNode>, Range>{},
         [&](auto id) {
-          interpolate(id(rowBasis), rhs, values_, dirichletNodes_);
+          auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), rowNode.treePath());
+          Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_);
         });
 
-      Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename ColBasis::LocalView::Tree>, Range>{},
+      Dune::Hybrid::ifElse(std::is_same<RangeType_t<ColNode>, Range>{},
         [&](auto id) {
-          interpolate(id(colBasis), solution, values_, dirichletNodes_);
+          auto subBasis = Dune::Functions::subspaceBasis(solution.basis(), colNode.treePath());
+          Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_);
         });
 
       // subtract columns of dirichlet nodes from rhs
-      // for (auto const& triplet : columns)
-      //   rhs[triplet.row] -= triplet.value * solution[triplet.col];
+      for (auto const& triplet : columns)
+        rhs[triplet.row] -= triplet.value * solution[triplet.col];
+
       initialized_ = false;
     }
 
   protected:
-    template <class RowBasis>
-    void initImpl(RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag);
-
-    template <class RowBasis>
-    void initImpl(RowBasis const& rowBasis, Dune::TypeTree::PowerNodeTag);
-
-    template <class RowBasis>
-    void initImpl(RowBasis const& rowBasis, Dune::TypeTree::CompositeNodeTag);
-
-    template <class RB, class RowNodeTag>
-    void initImpl(RB const&, RowNodeTag) {}
+    // fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the boundary or not.
+    template <class SubBasis>
+    void initImpl(SubBasis const& subBasis)
+    {
+      using Dune::Functions::forEachBoundaryDOF;
+      using LV = typename SubBasis::LocalView;
+      using I = typename SubBasis::GridView::Intersection;
+      dirichletNodes_.resize(subBasis.dimension());
+      forEachBoundaryDOF(subBasis, [&](int localIndex, LV const& localView, I const& intersection) {
+        dirichletNodes_[localView.index(localIndex)] = predicate_(intersection.geometry().center());
+      });
+    }
 
   private:
     std::function<bool(Domain)> predicate_;
@@ -119,18 +116,16 @@ namespace AMDiS
   };
 
 
-  template <class GlobalBasis>
+  template <class Basis>
   struct DirichletData
   {
-    using WorldVector = typename GlobalBasis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
+    using WorldVector = typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
 
     template <class RowNode, class ColNode>
     using type = std::list< std::shared_ptr<DirichletBC<WorldVector, RangeType_t<RowNode>>> >;
   };
 
-  template <class GlobalBasis>
-  using Constraints = MatrixData<GlobalBasis, DirichletData<GlobalBasis>::template type>;
+  template <class RowBasis, class ColBasis>
+  using Constraints = MatrixData<RowBasis, ColBasis, DirichletData<RowBasis>::template type>;
 
 } // end namespace AMDiS
-
-#include "DirichletBC.inc.hpp"

src/amdis/DirichletBC.inc.hpp

-#pragma once
-
-#include <dune/functions/functionspacebases/boundarydofs.hh>
-#include <dune/functions/functionspacebases/interpolate.hh>
-#include <dune/functions/functionspacebases/subspacebasis.hh>
-
-namespace AMDiS
-{
-  template <class WorldVector, class Range>
-    template <class RowBasis>
-  void DirichletBC<WorldVector, Range>::initImpl(
-      RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag)
-  {
-    using Dune::Functions::forEachBoundaryDOF;
-    using LV = typename RowBasis::LocalView;
-    using I = typename RowBasis::GridView::Intersection;
-    dirichletNodes_.resize(rowBasis.dimension());
-    forEachBoundaryDOF(rowBasis, [&](int localIndex, LV const& localView, I const& intersection) {
-      dirichletNodes_[localView.index(localIndex)] = predicate_(intersection.geometry().center());
-    });
-  }
-
-  template <class WorldVector, class Range>
-    template <class RowBasis>
-  void DirichletBC<WorldVector, Range>::initImpl(
-      RowBasis const& rowBasis, Dune::TypeTree::PowerNodeTag)
-  {
-    using Dune::Functions::subspaceBasis;
-    using Node = typename RowBasis::LocalView::Tree;
-    using ChildNode = typename Node::template Child<0>::type;
-
-    auto tp = rowBasis.prefixPath();
-    auto const& basis = rowBasis.rootBasis();
-    for (std::size_t i = 0; i < degree(rowBasis.localView().tree()); ++i)
-      initImpl(subspaceBasis(basis, push_back(tp,i)), typename ChildNode::NodeTag{});
-  }
-
-  template <class WorldVector, class Range>
-    template <class RowBasis>
-  void DirichletBC<WorldVector, Range>::initImpl(
-      RowBasis const& rowBasis, Dune::TypeTree::CompositeNodeTag)
-  {
-    using Dune::Functions::subspaceBasis;
-    using Node = typename RowBasis::LocalView::Tree;
-    auto tp = rowBasis.prefixPath();
-    auto const& basis = rowBasis.rootBasis();
-    forEach(range_<0,Node::CHILDREN>, [&,this](auto const _i) {
-      using ChildNode = typename Node::template Child<decltype(_i)::value>::type;
-      this->initImpl(subspaceBasis(basis, push_back(tp,_i)), typename ChildNode::NodeTag{});
-    });
-  }
-
-} // end namespace AMDiS

src/amdis/LocalAssemblerList.hpp

@@ -3,11 +3,19 @@
 #include <list>
 #include <memory>
 
+#include <amdis/Boundary.hpp>
 #include <amdis/LocalAssemblerBase.hpp>
 #include <amdis/utility/TreeData.hpp>
 
 namespace AMDiS
 {
+  namespace tag
+  {
+    template <class E> struct element_operator { using type = E; };
+    template <class I> struct intersection_operator { using type = I; };
+    template <class I> struct boundary_operator { using type = I; BoundaryType id; };
+  }
+
   template <class GridView>
   class OperatorLists
   {
@@ -18,8 +26,6 @@ namespace AMDiS
     struct DataElement
     {
       std::shared_ptr<OperatorType> op;
-      double* factor = nullptr;
-      double* estFactor = nullptr;
       BoundaryType b = {0};
     };
 
@@ -39,9 +45,9 @@ namespace AMDiS
       }
 
       /// Test whether to assemble on the node
-      bool doAssemble(bool flag) const
+      bool doAssemble() const
       {
-        return flag && (!assembled || changing);
+        return !assembled || changing;
       }
 
       /// Bind all operators to the grid element and geometry
@@ -61,6 +67,26 @@ namespace AMDiS
         for (auto& scaled : intersection) scaled.op->unbind();
       }
 
+
+      template <class... Args>
+      void push(tag::element_operator<Element>, Args&&... args)
+      {
+        element.push_back({std::forward<Args>(args)...});
+      }
+
+      template <class... Args>
+      void push(tag::boundary_operator<Intersection> b, Args&&... args)
+      {
+        boundary.push_back({std::forward<Args>(args)..., b.id});
+      }
+
+      template <class... Args>
+      void push(tag::intersection_operator<Intersection>, Args&&... args)
+      {
+        intersection.push_back({std::forward<Args>(args)...});
+      }
+
+
       /// List of operators to be assembled on grid elements
       std::list<DataElement<ElementOperator>> element;
       /// List of operators to be assembled on boundary intersections
@@ -72,7 +98,7 @@ namespace AMDiS
       bool assembled = false;
 
       /// if true, or \ref assembled false, do reassemble of all operators
-      bool changing = false;
+      bool changing = true;
     };
 
   public:
@@ -87,8 +113,8 @@ namespace AMDiS
   };