From 1f6515fec8d94c22a00a783d7c9b35d113314858 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Wed, 4 Nov 2020 10:59:13 +0100
Subject: [PATCH] Allow storage of GridFunctions in DirichletBC

---
 amdis/DirichletBC.hpp     | 43 +++++++++++++++++----------------------
 amdis/DirichletBC.inc.hpp | 14 ++++++-------
 2 files changed, 25 insertions(+), 32 deletions(-)

diff --git a/amdis/DirichletBC.hpp b/amdis/DirichletBC.hpp
index 413fbe33..91a1eb24 100644
--- a/amdis/DirichletBC.hpp
+++ b/amdis/DirichletBC.hpp
@@ -40,17 +40,18 @@ namespace AMDiS
    * \tparam CB   Basis of the column FE space
    * \tparam CTP  Treepath to the column node this constraint applies to
    **/
-  template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP>
+  template <class Mat, class Sol, class Rhs, class RowSubBasis, class ColSubBasis, class ValueGridFct>
   class DirichletBC
       : public BoundaryCondition<Mat, Sol, Rhs>
   {
     using Super = BoundaryCondition<Mat, Sol, Rhs>;
-    using RowSubBasis = Dune::Functions::SubspaceBasis<RB, RTP>;
-    using ColSubBasis = Dune::Functions::SubspaceBasis<CB, CTP>;
+
     // We assume CB's GridView to be the same as RB's
-    using Domain = typename RB::GridView::template Codim<0>::Geometry::GlobalCoordinate;
-    using Intersection = typename RB::GridView::Intersection;
-    using Range = RangeType_t<typename ColSubBasis::LocalView::Tree>;
+    using GridView = typename RowSubBasis::GridView;
+    using Intersection = typename GridView::Intersection;
+
+    using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+    using Range = TYPEOF(std::declval<ValueGridFct>()(std::declval<Domain>()));
 
   public:
     /// Constructor accepting subspacebases.
@@ -64,17 +65,6 @@ namespace AMDiS
       , values_(FWD(values))
     {}
 
-    /// Constructor accepting global bases and treepaths.
-    template <class Values,
-      REQUIRES(Concepts::Functor<Values, Range(Domain)>)>
-    DirichletBC(RB const& rowBasis, RTP const& rowTreePath,
-                CB const& colBasis, CTP const& colTreePath,
-                BoundarySubset<Intersection> boundarySubset, Values&& values)
-      : DirichletBC(Dune::Functions::subspaceBasis(rowBasis, rowTreePath),
-                    Dune::Functions::subspaceBasis(colBasis, colTreePath),
-                    std::move(boundarySubset), FWD(values))
-    {}
-
     // fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the boundary or not.
     /**
      * \see BoundaryCondition::init
@@ -92,7 +82,8 @@ namespace AMDiS
     RowSubBasis rowBasis_;
     ColSubBasis colBasis_;
     BoundarySubset<Intersection> boundarySubset_;
-    std::function<Range(Domain)> values_;
+    ValueGridFct values_;
+
     std::vector<bool> dirichletNodes_;
   };
 
@@ -104,7 +95,9 @@ namespace AMDiS
                        BoundarySubset<typename RB::GridView::Intersection> boundarySubset,
                        Values&& values)
   {
-    using BcType = DirichletBC<Mat, Sol, Rhs, RB, RTP, CB, CTP>;
+    using RowSubBasis = Dune::Functions::SubspaceBasis<RB, RTP>;
+    using ColSubBasis = Dune::Functions::SubspaceBasis<CB, CTP>;
+    using BcType = DirichletBC<Mat, Sol, Rhs, RowSubBasis, ColSubBasis, TYPEOF(values)>;
     return BcType(std::move(rowBasis), std::move(colBasis), std::move(boundarySubset), FWD(values));
   }
 
@@ -117,9 +110,10 @@ namespace AMDiS
                        BoundarySubset<typename RB::GridView::Intersection> boundarySubset,
                        Values&& values)
   {
-    using BcType = DirichletBC<Mat, Sol, Rhs, RB, RTP, CB, CTP>;
-    return BcType(rowBasis, rowTreePath, colBasis, colTreePath,
-                  std::move(boundarySubset), FWD(values));
+    return makeDirichletBC<Mat,Sol,Rhs>(
+        Dune::Functions::subspaceBasis(rowBasis, rowTreePath),
+        Dune::Functions::subspaceBasis(colBasis, colTreePath),
+        std::move(boundarySubset), FWD(values));
   }
 
   /// Make a DirichletBC from a global row- and colbasis
@@ -130,8 +124,9 @@ namespace AMDiS
                        BoundarySubset<typename RB::GridView::Intersection> boundarySubset,
                        Values&& values)
   {
-    return makeDirichletBC<Mat, Sol, Rhs>(rowBasis, makeTreePath(), colBasis, makeTreePath(),
-                                          std::move(boundarySubset), FWD(values));
+    return makeDirichletBC<Mat, Sol, Rhs>(
+        rowBasis, makeTreePath(), colBasis, makeTreePath(),
+        std::move(boundarySubset), FWD(values));
   }
 
 } // end namespace AMDiS
diff --git a/amdis/DirichletBC.inc.hpp b/amdis/DirichletBC.inc.hpp
index be722778..fe99a52e 100644
--- a/amdis/DirichletBC.inc.hpp
+++ b/amdis/DirichletBC.inc.hpp
@@ -33,11 +33,11 @@ void forEachInteriorBoundaryDOF(Basis const& basis, F&& f)
 } // end namespace Impl
 
 
-template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP>
-void DirichletBC<Mat, Sol, Rhs, RB, RTP, CB, CTP>::
+template <class Mat, class Sol, class Rhs, class RB, class CB, class V>
+void DirichletBC<Mat, Sol, Rhs, RB, CB, V>::
 init()
 {
-  using LV = typename ColSubBasis::LocalView;
+  using LV = typename CB::LocalView;
   dirichletNodes_.assign(colBasis_.dimension(), false);
   Impl::forEachInteriorBoundaryDOF(colBasis_,
     [&](int localIndex, LV const& localView, Intersection const& intersection) {
@@ -47,13 +47,11 @@ init()
 }
 
 
-template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP>
-void DirichletBC<Mat, Sol, Rhs, RB, RTP, CB, CTP>::
+template <class Mat, class Sol, class Rhs, class RB, class CB, class V>
+void DirichletBC<Mat, Sol, Rhs, RB, CB, V>::
 apply(Mat& matrix, Sol& solution, Rhs& rhs)
 {
-  auto&& gf = makeGridFunction(values_, colBasis_.gridView());
-  AMDiS::interpolate(colBasis_, solution, gf, tag::defaulted{}, tag::defaulted{}, dirichletNodes_);
-
+  AMDiS::interpolate(colBasis_, solution, values_, tag::defaulted{}, tag::defaulted{}, dirichletNodes_);
   dirichletBC(matrix, solution, rhs, dirichletNodes_);
 }
 
-- 
GitLab