getSolution return view, DirichletBC implemented recursively

examples/ellipt.cc

@@ -38,8 +38,6 @@ int main(int argc, char** argv)
   auto dbcValues = [](auto const& x){ return 0.0; }; // set value
   prob.addDirichletBC(predicate, _0, _0, dbcValues);
 
-  *prob.getSolution() = 0.0; // maybe not necessary
-
   prob.buildAfterCoarsen(adaptInfo, Flag(0));
 
   // write matrix to file

examples/init/stokes.dat.2d

 dimension of world: 2
 
 stokesMesh->macro file name: ./macro/macro.stand.2d
-stokesMesh->global refinements: 6
+stokesMesh->global refinements: 0
 
 stokes->mesh: stokesMesh
 
 stokes->solver->name: bicgstab_ell
-stokes->solver->ell: 5
+stokes->solver->ell: 3
 stokes->solver->max iteration: 1000
 stokes->solver->reduction: 1e-8
 stokes->solver->restart: 50
@@ -26,7 +26,7 @@ stokes->output[1]->ParaView format: 1
 stokes->output[1]->ParaView mode: 1
 
 stokes->viscosity: 0.1
-stokes->viscosity: 1.0
+stokes->density: 1.0
 stokes->boundary velocity: 10.0
 
 adapt->timestep: 0.1

examples/vecellipt.cc

@@ -36,8 +36,6 @@ int main(int argc, char** argv)
   auto dbcValues = [](auto const& x){ return 0.0; }; // set value
   prob.addDirichletBC(predicate, 1, 1, dbcValues);
 
-  *prob.getSolution() = 0.0; // maybe not necessary
-
   prob.buildAfterCoarsen(adaptInfo, Flag(0));
 
   // write matrix to file

src/amdis/DirichletBC.hpp

@@ -49,8 +49,11 @@ namespace AMDiS
               RowBasis const& rowBasis,
               ColBasis const& /*colBasis*/)
     {
-      using RowNode = typename RowBasis::LocalView::Tree;
-      initImpl(rowBasis, typename RowNode::NodeTag{});
+      if (!initialized_) {
+        using RowNode = typename RowBasis::LocalView::Tree;
+        initImpl(rowBasis, typename RowNode::NodeTag{});
+        initialized_ = true;
+      }
     }
 
 
@@ -69,7 +72,14 @@ namespace AMDiS
                 RowBasis const& rowBasis,
                 ColBasis const& colBasis)
     {
+      test_exit_dbg(initialized_, "Boundary condition not initialized!");
+      auto columns = matrix.applyDirichletBC(dirichletNodes_);
+
       finishImpl(matrix, rhs, solution, rowBasis, colBasis, ValueCategory_t<Range>{});
+
+      // subtract columns of dirichlet nodes from rhs
+      for (auto const& triplet : columns)
+        rhs[triplet.row] -= triplet.value * solution[triplet.col];
     }
 
   protected:
@@ -79,6 +89,9 @@ namespace AMDiS
     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) {}
 

src/amdis/DirichletBC.inc.hpp

@@ -13,11 +13,7 @@ namespace AMDiS
       RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag)
   {
     using Dune::Functions::interpolate;
-
-    if (!initialized_) {
-      interpolate(rowBasis, hierarchicVectorWrapper(dirichletNodes_), predicate_);
-      initialized_ = true;
-    }
+    interpolate(rowBasis, hierarchicVectorWrapper(dirichletNodes_), predicate_);
   }
 
   template <class WorldVector, class Range>
@@ -25,16 +21,29 @@ namespace AMDiS
   void DirichletBC<WorldVector, Range>::initImpl(
       RowBasis const& rowBasis, Dune::TypeTree::PowerNodeTag)
   {
-    using Dune::Functions::interpolate;
+    using Dune::Functions::subspaceBasis;
+    using Node = typename RowBasis::LocalView::Tree;
+    using ChildNode = typename Node::template Child<0>::type;
 
-    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;
-    }
+    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{});
+    });
   }
 
 
@@ -45,18 +54,10 @@ namespace AMDiS
       RowBasis const& rowBasis, ColBasis const& colBasis, ValueCat)
   {
     using Dune::Functions::interpolate;
-
-    test_exit_dbg(initialized_, "Boundary condition not initialized!");
-    auto columns = matrix.applyDirichletBC(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];
   }
 
 

src/amdis/ProblemInstat.inc.hpp

@@ -59,7 +59,7 @@ template <class Traits>
 void ProblemInstat<Traits>::initTimestep(AdaptInfo&)
 {
   if (oldSolution)
-    oldSolution->copy(*problemStat.getSolution());
+    oldSolution->copy(*problemStat.getSolutionVec());
 }
 
 } // end namespace AMDiS

src/amdis/ProblemStat.hpp

@@ -207,15 +207,24 @@ namespace AMDiS
     auto getSystemMatrix() const { return systemMatrix; }
 
 
-    /// Return a pointer to the solution, \ref solution
-    std::shared_ptr<SystemVector>       getSolution()       { return solution; }
-    std::shared_ptr<SystemVector const> getSolution() const { return solution; }
+    std::shared_ptr<SystemVector> getSolutionVec() const
+    {
+      return solution;
+    }
+
+    /// Return a mutable view to a solution component
+    template <class TreePath = RootTreePath>
+    auto getSolution(TreePath const& path = {})
+    {
+      auto&& tp = makeTreePath(path);
+      return makeDOFVectorView(*solution, tp);
+    }
 
-    /// Return a view to a solution component (as shared_ptr)
-    template <class TreePath>
-    auto getSolution(TreePath path) const
+    /// Return a const view to a solution component
+    template <class TreePath = RootTreePath>
+    auto getSolution(TreePath const& path = {}) const
     {
-      auto tp = makeTreePath(path);
+      auto&& tp = makeTreePath(path);
       return makeDOFVectorView(*solution, tp);
     }
 
@@ -343,7 +352,7 @@ namespace AMDiS
     std::vector<Grid*> componentGrids;
 
     /// FE spaces of this problem.
-    std::shared_ptr<GlobalBasis> globalBasis; // eventuell const
+    std::shared_ptr<GlobalBasis> globalBasis;
     std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
 
     /// A FileWriter object

src/amdis/ProblemStat.inc.hpp

@@ -65,7 +65,7 @@ void ProblemStat<Traits>::initialize(
 
     if (adoptProblem &&
         (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
-      globalBasis = adoptProblem->getGlobalBasis();
+      globalBasis = adoptProblem->globalBasis;
     }
   }
 
@@ -78,9 +78,9 @@ void ProblemStat<Traits>::initialize(
     createMatricesAndVectors();
 
   if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
-    solution = adoptProblem->getSolution();
-    rhs = adoptProblem->getRhs();
-    systemMatrix = adoptProblem->getSystemMatrix();
+    solution = adoptProblem->solution;
+    rhs = adoptProblem->rhs;
+    systemMatrix = adoptProblem->systemMatrix;
   }
 
 
@@ -94,7 +94,7 @@ void ProblemStat<Traits>::initialize(
 
     if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
       test_exit(!linearSolver, "solver already created\n");
-      linearSolver = adoptProblem->getSolver();
+      linearSolver = adoptProblem->linearSolver;
     }
   }
 

src/amdis/utility/TreePath.hpp

@@ -10,6 +10,8 @@
 
 namespace AMDiS
 {
+  struct RootTreePath {};
+
   namespace Concepts
   {
     /** \addtogroup Concepts
@@ -38,6 +40,11 @@ namespace AMDiS
         : std::true_type
       {};
 
+      template <>
+      struct IsPreTreePath<RootTreePath>
+        : std::true_type
+      {};
+
     } // end namespace Definition
 
     template <class TP>
@@ -72,26 +79,39 @@ namespace AMDiS
 
 #else // DOXYGEN
 
-  auto makeTreePath(int i) { return Dune::TypeTree::hybridTreePath(std::size_t(i)); }
+  auto makeTreePath(int i)         { return Dune::TypeTree::hybridTreePath(std::size_t(i)); }
   auto makeTreePath(std::size_t i) { return Dune::TypeTree::hybridTreePath(i); }
+  auto makeTreePath(RootTreePath)  { return Dune::TypeTree::hybridTreePath(); }
+
 
   template <int I>
-  auto makeTreePath(int_t<I>) { return Dune::TypeTree::hybridTreePath(index_<std::size_t(I)>); }
+  auto makeTreePath(int_t<I>)      { return Dune::TypeTree::hybridTreePath(index_<std::size_t(I)>); }
+
+  template <int... I>
+  auto makeTreePath(Ints<I...>)    { return Dune::TypeTree::hybridTreePath(index_<std::size_t(I)>...); }
 
   template <std::size_t I>
   auto makeTreePath(index_t<I> _i) { return Dune::TypeTree::hybridTreePath(_i); }
 
-  template <class... T>
-  auto const& makeTreePath(Dune::TypeTree::HybridTreePath<T...> const& tp) { return tp; }
-
   template <std::size_t... I>
   auto makeTreePath(Indices<I...>) { return Dune::TypeTree::hybridTreePath(index_<I>...); }
 
+
+  template <class... T>
+  auto const& makeTreePath(Dune::TypeTree::HybridTreePath<T...> const& tp)
+  {
+    return tp;
+  }
+
   template <std::size_t... I>
-  auto makeTreePath(Dune::TypeTree::TreePath<I...>) { return Dune::TypeTree::hybridTreePath(index_<I>...); }
+  auto makeTreePath(Dune::TypeTree::TreePath<I...>)
+  {
+    return Dune::TypeTree::hybridTreePath(index_<I>...);
+  }
 
   template <class TP>
-  auto makeTreePath(TP const&) {
+  auto makeTreePath(TP const&)
+  {
     static_assert( Concepts::PreTreePath<TP>,
       "Argument must be a valid treepath, or an integer/index-constant");
     return Dune::TypeTree::hybridTreePath();