some examples modified to make it running

dune/amdis/Assembler.hpp

@@ -6,6 +6,7 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
+#include <dune/amdis/DirichletBC.hpp>
 #include <dune/amdis/LinearAlgebra.hpp>
 #include <dune/amdis/Operator.hpp>
 #include <dune/amdis/common/Mpl.hpp>
@@ -23,10 +24,12 @@ namespace AMDiS
     /// Constructor, stores a shared-pointer to the feSpaces
     Assembler(GlobalBasis& globalBasis,
               MatrixOperators<GlobalBasis>& matrixOperators,
-              VectorOperators<GlobalBasis>& rhsOperators)
+              VectorOperators<GlobalBasis>& rhsOperators,
+              Constraints<GlobalBasis>& constraints)
       : globalBasis_(globalBasis)
       , matrixOperators_(matrixOperators)
       , rhsOperators_(rhsOperators)
+      , constraints_(constraints)
     {}
 
     void update(GridView const& gv)
@@ -88,6 +91,7 @@ namespace AMDiS
     GlobalBasis& globalBasis_;
     MatrixOperators<GlobalBasis>& matrixOperators_;
     VectorOperators<GlobalBasis>& rhsOperators_;
+    Constraints<GlobalBasis>& constraints_;
   };
 
 } // end namespace AMDiS

dune/amdis/Assembler.inc.hpp

 #pragma once
 
 #include <dune/functions/functionspacebases/subspacebasis.hh>
-#include <dune/typetree/treepath.hh>
+#include <dune/amdis/utility/TreePath.hpp>
 #include <dune/amdis/utility/Visitor.hpp>
 
 namespace AMDiS {
@@ -64,15 +64,19 @@ void Assembler<GlobalBasis>::assemble(
     for (std::size_t i = 0; i < localView.size(); ++i) {
       auto const row = localIndexSet.index(i);
       for (std::size_t j = 0; j < localView.size(); ++j) {
-        auto const col = localIndexSet.index(j);
-        matrix(row,col) += elementMatrix(i,j);
+        if (elementMatrix(i,j) != 0.0) {
+          auto const col = localIndexSet.index(j);
+          matrix(row,col) += elementMatrix(i,j);
+        }
       }
     }
 
     // add element-vector to system-vector
     for (std::size_t i = 0; i < localView.size(); ++i) {
-      auto const idx = localIndexSet.index(i);
-      rhs[idx] += elementVector[i];
+      if (elementVector[i] != 0.0) {
+        auto const idx = localIndexSet.index(i);
+        rhs[idx] += elementVector[i];
+      }
     }
 
     localIndexSet.unbind();
@@ -86,43 +90,6 @@ void Assembler<GlobalBasis>::assemble(
 }
 
 
-template <class GlobalBasis>
-  template <class SystemMatrixType, class SystemVectorType>
-void Assembler<GlobalBasis>::initMatrixVector(
-    SystemMatrixType& matrix,
-    SystemVectorType& solution,
-    SystemVectorType& rhs,
-    bool asmMatrix, bool asmVector) const
-{
-  matrix.init(asmMatrix);
-  solution.compress();
-  rhs.compress();
-
-  auto localView = globalBasis_.localView();
-  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
-  {
-    if (rowNode.isLeaf)
-      msg(0, " DOFs for Basis[", 1, "]"); // TODO: add right values
-
-    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
-    if (rhsOperators_[rowNode].assemble(asmVector))
-      rhsOperators_[rowNode].init(rowBasis);
-
-    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
-    {
-      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
-      if (matrixOperators_[rowNode][colNode].assemble(asmMatrix))
-        matrixOperators_[rowNode][colNode].init(rowBasis, colBasis);
-
-      // init boundary condition
-      // for (int c = 0; c < nComponents; ++c)
-      //   for (auto bc : matrixOperators[R][c].dirichlet)
-      //     bc->init(c == C, matrix(_r, _c), solution[_c], rhs[_r]);
-    });
-  });
-}
-
-
 template <class GlobalBasis>
   template <class ElementContainer, class Container, class Operators, class... LocalViews>
 void Assembler<GlobalBasis>::assembleElementOperators(
@@ -160,6 +127,44 @@ void Assembler<GlobalBasis>::assembleElementOperators(
 }
 
 
+template <class GlobalBasis>
+  template <class SystemMatrixType, class SystemVectorType>
+void Assembler<GlobalBasis>::initMatrixVector(
+    SystemMatrixType& matrix,
+    SystemVectorType& solution,
+    SystemVectorType& rhs,
+    bool asmMatrix, bool asmVector) const
+{
+  matrix.init(asmMatrix);
+  solution.compress();
+  rhs.compress();
+
+  auto localView = globalBasis_.localView();
+  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
+  {
+    if (rowNode.isLeaf)
+      msg(0, " DOFs for Basis[", to_string(rowTreePath), "]"); // TODO: add right values
+
+    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
+    if (rhsOperators_[rowNode].assemble(asmVector))
+      rhsOperators_[rowNode].init(rowBasis);
+
+    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
+    {
+      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
+      if (matrixOperators_[rowNode][colNode].assemble(asmMatrix))
+        matrixOperators_[rowNode][colNode].init(rowBasis, colBasis);
+
+      for (auto bc : constraints_[rowNode][colNode].scalar)
+        bc->init(matrix, solution, rhs, rowBasis, colBasis);
+      for (auto bc : constraints_[rowNode][colNode].vector)
+        bc->init(matrix, solution, rhs, rowBasis, colBasis);
+    });
+  });
+  msg(globalBasis_.dimension(), " total DOFs");
+}
+
+
 template <class GlobalBasis>
   template <class SystemMatrixType, class SystemVectorType>
 std::size_t Assembler<GlobalBasis>::finishMatrixVector(
@@ -171,27 +176,25 @@ std::size_t Assembler<GlobalBasis>::finishMatrixVector(
   matrix.finish();
 
   auto localView = globalBasis_.localView();
-  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto /*rowTreePath*/)
+  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
   {
+    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
     auto& rhsOp = rhsOperators_[rowNode];
     if (rhsOp.assemble(asmVector))
       rhsOp.assembled = true;
 
-    forEachNode(localView.tree(), [&,this](auto const& colNode, auto /*colTreePath*/)
+    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
     {
+      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
       auto& matOp = matrixOperators_[rowNode][colNode];
       if (matOp.assemble(asmMatrix))
         matOp.assembled = true;
 
       // finish boundary condition
-      // for (int c = 0; c < nComponents; ++c) {
-      //   for (int r = 0; r < nComponents; ++r) {
-      //     if (r != R && c != C)
-      //       continue;
-      //     for (auto bc : matrixOperators[r][c].dirichlet)
-      //       bc->finish(r == R, c == C, matrix(_r, _c), solution[_c], rhs[_r]);
-      //   }
-      // }
+      for (auto bc : constraints_[rowNode][colNode].scalar)
+        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
+      for (auto bc : constraints_[rowNode][colNode].vector)
+        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
     });
   });
 

dune/amdis/DirichletBC.hpp

@@ -4,9 +4,9 @@
 #include <type_traits>
 #include <vector>
 
-#include <dune/amdis/common/Concepts.hpp>
-
 #include <dune/amdis/Output.hpp>
+#include <dune/amdis/common/Concepts.hpp>
+#include <dune/amdis/utility/TreeData.hpp>
 
 namespace AMDiS
 {
@@ -25,13 +25,13 @@ namespace AMDiS
    * support this symmetric modification. As a result, this method returns a list
    * of columns values, that should be subtracted from the rhs.
    **/
-  template <class WorldVector>
+  template <class Domain, class Range = double>
   class DirichletBC
   {
   public:
     template <class Predicate, class Values,
-      REQUIRES(Concepts::Functor<Predicate, bool(WorldVector)> &&
-				       Concepts::Functor<Values,  double(WorldVector)>) >
+      REQUIRES(Concepts::Functor<Predicate, bool(Domain)> &&
+				       Concepts::Functor<Values,  Range(Domain)>) >
     DirichletBC(Predicate&& predicate, Values&& values)
       : predicate(std::forward<Predicate>(predicate))
       , values(std::forward<Values>(values))
@@ -41,11 +41,15 @@ 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>
-    void init(bool apply,
-              Matrix& matrix,
-              VectorX& rhs,
-              VectorB& solution);
+    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*/)
+    {
+      using RowNode = typename RowBasis::LocalView::Tree;
+      initImpl(rowBasis, typename RowNode::NodeTag{});
+    }
+
 
 
     /// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row
@@ -55,21 +59,51 @@ namespace AMDiS
       * 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>
-    void finish(bool apply_row,
-                bool apply_col,
-                Matrix& matrix,
+    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
+    void finish(Matrix& matrix,
                 VectorX& rhs,
-                VectorB& solution);
+                VectorB& solution,
+                RowBasis const& rowBasis,
+                ColBasis const& colBasis);
+
+  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 RB, class RowNodeTag>
+    void initImpl(RB const&, RowNodeTag) {}
 
   private:
-    std::function<bool(WorldVector)> predicate;
-    std::function<double(WorldVector)> values;
+    std::function<bool(Domain)> predicate;
+    std::function<Range(Domain)> values;
 
     bool initialized = false;
     std::vector<char> dirichletNodes;
   };
 
+
+  template <class GlobalBasis>
+  struct DirichletData
+  {
+    using WorldVector = typename GlobalBasis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
+
+    template <class Node>
+    struct type
+    {
+      std::list<std::shared_ptr<DirichletBC<WorldVector, double>>> scalar;
+      std::list<std::shared_ptr<DirichletBC<WorldVector, WorldVector>>> vector;
+
+      void push_back(std::shared_ptr<DirichletBC<WorldVector, double>> const& bc) { scalar.push_back(bc); }
+      void push_back(std::shared_ptr<DirichletBC<WorldVector, WorldVector>> const& bc) { vector.push_back(bc); }
+    };
+  };
+
+  template <class GlobalBasis>
+  using Constraints = MatrixData<GlobalBasis, DirichletData<GlobalBasis>::template type>;
+
 } // end namespace AMDiS
 
 #include "DirichletBC.inc.hpp"

dune/amdis/DirichletBC.inc.hpp

 #pragma once
 
 #include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
 #include <dune/amdis/LinearAlgebra.hpp>
+#include <dune/amdis/linear_algebra/HierarchicWrapper.hpp>
 
 namespace AMDiS
 {
-  template <class WorldVector>
-    template <class Matrix, class VectorX, class VectorB>
-  void DirichletBC<WorldVector>::init(bool /*apply*/,
-				                              Matrix& matrix,
-                                      VectorX& /*solution*/,
-				                              VectorB& /*rhs*/)
+  template <class WorldVector, class Range>
+    template <class RowBasis>
+  void DirichletBC<WorldVector, Range>::initImpl(
+      RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag)
   {
     using Dune::Functions::interpolate;
 
     if (!initialized) {
-      interpolate(matrix.getRowFeSpace(), dirichletNodes, predicate);
+      interpolate(rowBasis, hierarchicVectorWrapper(dirichletNodes), predicate);
       initialized = true;
     }
   }
 
+  template <class WorldVector, class Range>
+    template <class RowBasis>
+  void DirichletBC<WorldVector, Range>::initImpl(
+      RowBasis const& rowBasis, Dune::TypeTree::PowerNodeTag)
+  {
+    using Dune::Functions::interpolate;
 
-  template <class WorldVector>
-    template <class Matrix, class VectorX, class VectorB>
-  void DirichletBC<WorldVector>::finish(bool apply_row,
-                                        bool apply_col,
-					                              Matrix& matrix,
-                                        VectorX& solution,
-					                              VectorB& rhs)
+    if (!initialized) {
+      auto tp = rowBasis.prefixPath();
+      auto const& basis = rowBasis.rootBasis();
+      for (std::size_t i = 0; i < degree(rowBasis.localView().tree()); ++i)
+        interpolate(Dune::Functions::subspaceBasis(basis, push_back(tp,i)),
+          hierarchicVectorWrapper(dirichletNodes), predicate);
+      initialized = true;
+    }
+  }
+
+
+  template <class WorldVector, class Range>
+    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
+  void DirichletBC<WorldVector, Range>::finish(
+      Matrix& matrix, VectorX& solution, VectorB& rhs,
+      RowBasis const& rowBasis, ColBasis const& colBasis)
+      // Dune::TypeTree::LeafNodeTag, Dune::TypeTree::LeafNodeTag)
   {
     using Dune::Functions::interpolate;
 
     test_exit_dbg(initialized, "Boundary condition not initialized!");
-    auto columns = matrix.applyDirichletBC(dirichletNodes, apply_row, apply_col);
+    auto columns = matrix.applyDirichletBC(dirichletNodes);
 
-    if (apply_col) {
-      interpolate(matrix.getRowFeSpace(), wrapper(rhs.getVector()), values, dirichletNodes);
-      if (apply_row)
-        interpolate(matrix.getColFeSpace(), wrapper(solution.getVector()), values, dirichletNodes);
+    interpolate(rowBasis, hierarchicVectorWrapper(rhs.getVector()), values, hierarchicVectorWrapper(dirichletNodes));
+    interpolate(colBasis, hierarchicVectorWrapper(solution.getVector()), values, hierarchicVectorWrapper(dirichletNodes));
 
-      // subtract columns of dirichlet nodes from rhs
-      for (auto const& triplet : columns)
-        rhs[triplet.row] -= triplet.value * solution[triplet.col];
-    }
+    // subtract columns of dirichlet nodes from rhs
+    for (auto const& triplet : columns)
+      rhs[triplet.row] -= triplet.value * solution[triplet.col];
   }
 
 

dune/amdis/FirstOrderAssembler.hpp

@@ -60,19 +60,17 @@ namespace AMDiS
 
     template <class Operator, class RowView, class ColView>
     void calculateElementMatrix(Operator& op,
-                                RowView const& rowView,
-                                ColView const& colView,
+                                RowView const& rowView, ColView const& colView,
                                 ElementMatrix& elementMatrix,
                                 double fac,
-                                Dune::TypeTree::LeafNodeTag,
-                                Dune::TypeTree::LeafNodeTag,
+                                Dune::TypeTree::LeafNodeTag, Dune::TypeTree::LeafNodeTag,
                                 FirstOrderType_t<GRD_PHI>)
     {
       auto geometry = rowView.element().geometry();
       auto const& quad_geometry = Super::getGeometry();
 
-      auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
-      auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
+      auto const& rowLocalFE = rowView.tree().finiteElement();
+      auto const& colLocalFE = colView.tree().finiteElement();
       auto const& quad = Super::getQuadrature().getRule();
 
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
@@ -86,11 +84,11 @@ namespace AMDiS
         const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight() * fac;
 
         std::vector<Dune::FieldVector<double,1> > rowShapeValues;
-        rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
+        rowLocalFE.localBasis().evaluateFunction(quadPos, rowShapeValues);
 
         // The gradients of the shape functions on the reference element
         std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
-        colLocalBasis.evaluateJacobian(quadPos, colReferenceGradients);
+        colLocalFE.localBasis().evaluateJacobian(quadPos, colReferenceGradients);
 
         // Compute the shape function gradients on the real element
         std::vector<Dune::FieldVector<double,dow> > colGradients(colReferenceGradients.size());
@@ -98,8 +96,8 @@ namespace AMDiS
         for (std::size_t i = 0; i < colGradients.size(); ++i)
           jacobian.mv(colReferenceGradients[i][0], colGradients[i]);
 
-        for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
-          for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
+        for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
+          for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
             const int local_i = rowView.tree().localIndex(i);
             const int local_j = colView.tree().localIndex(j);
             op.evalFot1(elementMatrix[local_i][local_j], iq, factor, rowShapeValues[i], colGradients[j]);
@@ -110,19 +108,17 @@ namespace AMDiS
 
     template <class Operator, class RowView, class ColView>
     void calculateElementMatrix(Operator& op,
-                                RowView const& rowView,
-                                ColView const& colView,
+                                RowView const& rowView, ColView const& colView,
                                 ElementMatrix& elementMatrix,
                                 double fac,
-                                Dune::TypeTree::LeafNodeTag,
-                                Dune::TypeTree::LeafNodeTag,
+                                Dune::TypeTree::LeafNodeTag, Dune::TypeTree::LeafNodeTag,
                                 FirstOrderType_t<GRD_PSI>)
     {
       auto geometry = rowView.element().geometry();
       auto const& quad_geometry = Super::getGeometry();
 
-      auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
-      auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
+      auto const& rowLocalFE = rowView.tree().finiteElement();
+      auto const& colLocalFE = colView.tree().finiteElement();
       auto const& quad = Super::getQuadrature().getRule();
 
       for (std::size_t iq = 0; iq < quad.size(); ++iq) {
@@ -137,7 +133,7 @@ namespace AMDiS
 
         // The gradients of the shape functions on the reference element
         std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
-        rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);
+        rowLocalFE.localBasis().evaluateJacobian(quadPos, rowReferenceGradients);
 
         // Compute the shape function gradients on the real element
         std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());
@@ -146,10 +142,10 @@ namespace AMDiS
           jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
 
         std::vector<Dune::FieldVector<double,1> > colShapeValues;
-        colLocalBasis.evaluateFunction(quadPos, colShapeValues);
+        colLocalFE.localBasis().evaluateFunction(quadPos, colShapeValues);