assembler clened up, functors added

CMakeLists.txt

 cmake_minimum_required(VERSION 3.1)
 project(dune-amdis CXX)
 
-set(CXX_MAX_STANDARD 14 CACHE BOOL "" FORCE)
+#set(CXX_MAX_STANDARD 14 CACHE BOOL "" FORCE)
 
 if(NOT (dune-common_DIR OR dune-common_ROOT OR
       "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
@@ -20,11 +20,34 @@ include(DuneMacros)
 
 # start a dune project with information from dune.module
 dune_project()
-dune_enable_all_packages()
+dune_enable_all_packages(MODULE_LIBRARIES amdis)
+
+include(AmdisMacros)
 add_subdirectory("src")
 add_subdirectory("dune")
 add_subdirectory("doc")
 add_subdirectory("cmake/modules")
 
+
+# some additional packages and flags
+add_dune_alberta_flags("amdis" OBJECT USE_GENERIC)
+target_compile_definitions("amdis" PUBLIC AMDIS_BACKEND_MTL=1)
+target_compile_options("amdis" PUBLIC -ftemplate-backtrace-limit=0 -Wall -pedantic -Wno-unused-parameter)
+
+find_package(MTL REQUIRED
+             PATHS /usr/local/lib/mtl4 /opt/sources/mtl4 /opt/development/mtl4)
+if (MTL_FOUND)
+  target_include_directories("amdis" PUBLIC ${MTL_INCLUDE_DIRS})
+
+  set (CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
+  foreach (feature ${CXX_ELEVEN_FEATURE_LIST})
+      target_compile_definitions("amdis" PUBLIC MTL_WITH_${feature})
+  endforeach ()
+
+  if (HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_FOUND)
+      target_compile_definitions("amdis" PUBLIC MTL_HAS_UMFPACK)
+  endif ()
+endif (MTL_FOUND)
+
 # finalize the dune project, e.g. generating config.h etc.
 finalize_dune_project(GENERATE_CONFIG_H_CMAKE)

cmake/modules/AMDiSCXXFeatures.cmake

+#include(CheckIncludeFileCXX)
+include(CheckCXXSourceCompiles)
+#include(CheckCXXSymbolExists)
+
+# fold expressions (a + ...)
+check_cxx_source_compiles("
+  template <class... Args>
+  int f(Args... args)
+  {
+    return (args + ...);
+  }
+  int main()
+  {
+    f(0,1,2,3,4,5);
+  }
+"  AMDIS_HAS_CXX_FOLD_EXPRESSION
+)
+
+check_cxx_source_compiles("
+  template <int I>
+  auto f()
+  {
+    if constexpr(I == 0)
+      return double{1.0};
+    else
+      return int{0};
+  }
+  int main()
+  {
+    return f<1>();
+  }
+"  AMDIS_HAS_CXX_CONSTEXPR_IF
+)
\ No newline at end of file

cmake/modules/AmdisMacros.cmake

 # File for module specific CMake tests.
+
+include(AMDiSCXXFeatures)
\ No newline at end of file

config.h.cmake

@@ -40,6 +40,12 @@
 /* Define to the revision of amdis */
 #define AMDIS_VERSION_REVISION @AMDIS_VERSION_REVISION@
 
-/* end amdis
+
+/* some detected compiler features may be used in AMDiS */
+#cmakedefine AMDIS_HAS_CXX_FOLD_EXPRESSIONS 1
+
+#cmakedefine AMDIS_HAS_CXX_CONSTEXPR_IF 1
+
+/* end dune-amdis
    Everything below here will be overwritten
 */

dune/amdis/AdaptInfo.hpp

@@ -9,7 +9,7 @@
 
 // AMDiS includes
 #include <dune/amdis/Output.hpp>
-#include <dune/amdis/Math.hpp>
+#include <dune/amdis/common/Math.hpp>
 
 namespace AMDiS
 {

dune/amdis/Assembler.inc.hpp

@@ -4,7 +4,7 @@
 #include <dune/amdis/utility/TreePath.hpp>
 #include <dune/amdis/utility/Visitor.hpp>
 
-#include <dune/amdis/Math.hpp>
+#include <dune/amdis/common/Math.hpp>
 
 namespace AMDiS {
 

dune/amdis/CMakeLists.txt

 #install headers
 
-dune_add_library("duneamdis" NO_EXPORT
+dune_library_add_sources(amdis SOURCES
     AdaptBase.cpp
     AdaptInfo.cpp
     AdaptInstationary.cpp
@@ -11,33 +11,7 @@ dune_add_library("duneamdis" NO_EXPORT
     # ProblemInstat.cpp
     ProblemStat.cpp
     StandardProblemIteration.cpp
-    #linear_algebra/istl/SystemMatrix.cpp
-    #linear_algebra/istl/SystemVector.cpp
-    linear_algebra/mtl/SystemMatrix.cpp
-    linear_algebra/mtl/SystemVector.cpp
-    utility/Filesystem.cpp
-    )
-add_dune_alberta_flags("duneamdis" OBJECT USE_GENERIC)
-target_compile_definitions("duneamdis" PUBLIC AMDIS_BACKEND_MTL=1)
-target_compile_options("duneamdis" PUBLIC -ftemplate-backtrace-limit=0 -Wall -pedantic -Wno-unused-parameter)
-
-find_package(MTL REQUIRED
-             PATHS /usr/local/lib/mtl4)
-if (MTL_FOUND)
-  target_include_directories("duneamdis" PUBLIC ${MTL_INCLUDE_DIRS})
-#   target_link_libraries("duneamdis" PUBLIC ${MTL_LIBRARIES})
-#   target_compile_options("duneamdis" PUBLIC ${MTL_CXX_DEFINITIONS})
-
-  set (CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
-  foreach (feature ${CXX_ELEVEN_FEATURE_LIST})
-      target_compile_definitions("duneamdis" PUBLIC MTL_WITH_${feature})
-  endforeach ()
-
-  if (HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_FOUND)
-      target_compile_definitions("duneamdis" PUBLIC MTL_HAS_UMFPACK)
-  endif ()
-endif (MTL_FOUND)
-
+)
 
 install(FILES
     AdaptBase.hpp
@@ -45,73 +19,40 @@ install(FILES
     AdaptInstationary.hpp
     AdaptStationary.hpp
     AMDiS.hpp
+    Assembler.hpp
+    Assembler.inc.hpp
+    ContextGeometry.hpp
     CreatorInterface.hpp
     CreatorMap.hpp
     DirichletBC.hpp
     DirichletBC.inc.hpp
     FileWriter.hpp
     Flag.hpp
+    GridFunctionOperator.hpp
+    GridFunctions.hpp
     Initfile.hpp
     LinearAlgebra.hpp
-    Log.hpp
-    Math.hpp
+    LocalAssembler.hpp
+    LocalAssemblerBase.hpp
     Mesh.hpp
-    OperatorEvaluation.hpp
-    Operator.hpp
-    Operator.inc.hpp
-    OperatorTermBase.hpp
-    OperatorTerm.hpp
-    ProblemInstatBase.hpp
+    Operators.hpp
+    Output.hpp
     ProblemInstat.hpp
     ProblemInstat.inc.hpp
+    ProblemInstatBase.hpp
     ProblemInterationInterface.hpp
-    ProblemStatBase.hpp
     ProblemStat.hpp
     ProblemStat.inc.hpp
+    ProblemStatBase.hpp
+    ProblemStatTraits.hpp
     ProblemTimeInterface.hpp
     StandardProblemIteration.hpp
-    TermGenerator.hpp
-    common/ConceptsBase.hpp
-    common/ClonablePtr.hpp
-    common/IndexSeq.hpp
-    common/Literals.hpp
-    common/Loops.hpp
-    common/Mpl.hpp
-    common/ScalarTypes.hpp
-    common/Timer.hpp
-    common/Traits.hpp
-    common/TupleUtility.hpp
-    common/Utility.hpp
-    common/ValueCategory.hpp
-    linear_algebra/LinearAlgebraBse.hpp
-    linear_algebra/LinearSolverInterface.hpp
-    linear_algebra/PreconditionerInterface.hpp
-    linear_algebra/RunnerInterface.hpp
-    linear_algebra/SolverInfo.hpp
-    linear_algebra/istl/DOFMatrix.hpp
-    linear_algebra/istl/DOFVector.hpp
-    linear_algebra/istl/ISTL_Preconditioner.hpp
-    linear_algebra/istl/ISTLRunner.hpp
-    linear_algebra/istl/ISTL_Solver.hpp
-    linear_algebra/istl/LinearSolver.hpp
-    linear_algebra/istl/Preconditioner.hpp
-    linear_algebra/istl/SystemMatrix.hpp
-    linear_algebra/istl/SystemVector.hpp
-    linear_algebra/mtl/BITL_Preconditioner.hpp
-    linear_algebra/mtl/BITL_Solver.hpp
-    linear_algebra/mtl/BlockMTLMatrix.hpp
-    linear_algebra/mtl/BlockMTLVector.hpp
-    linear_algebra/mtl/Copy.hpp
-    linear_algebra/mtl/DOFMatrix.hpp
-    linear_algebra/mtl/DOFVector.hpp
-    linear_algebra/mtl/ITL_Preconditioner.hpp
-    linear_algebra/mtl/ITL_Solver.hpp
-    linear_algebra/mtl/KrylovRunner.hpp
-    linear_algebra/mtl/LinearSolver.hpp
-    linear_algebra/mtl/Mapper.hpp
-    linear_algebra/mtl/MTLDenseVector.hpp
-    linear_algebra/mtl/Preconditioner.hpp
-    linear_algebra/mtl/SystemMatrix.hpp
-    linear_algebra/mtl/SystemVector.hpp
-    linear_algebra/mtl/UmfpackRunner.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis)
+
+add_subdirectory("assembler")
+add_subdirectory("common")
+add_subdirectory("gridfunctions")
+add_subdirectory("io")
+add_subdirectory("linear_algebra")
+add_subdirectory("operations")
+add_subdirectory("utility")
\ No newline at end of file

dune/amdis/OperatorBase.hpp → dune/amdis/GridFunctionOperator.hpp

 #pragma once
 
 #include <cassert>
+#include <type_traits>
 
-#include <dune/amdis/terms/TermGenerator.hpp>
+#include <dune/amdis/gridfunctions/GridFunctionConcepts.hpp>
 #include <dune/amdis/utility/GetDegree.hpp>
 
 namespace AMDiS
 {
-  template <class Expr>
-  class ExpressionOperatorBase
+  template <class GridFunction>
+  class GridFunctionOperatorBase
   {
+    using LocalFunction = decltype(localFunction(std::declval<GridFunction>()));
+    using Element = typename GridFunction::EntitySet::Element;
+    using Geometry = typename Element::Geometry;
+
+    using LocalCoordinate = typename GridFunction::EntitySet::LocalCoordinate;
+
   public:
     /// \brief Constructor. Stores a copy of `expr`.
     /**
@@ -17,8 +24,9 @@ namespace AMDiS
      * \ref ExpressionBase, and stores a copy. Additionally, it gets the
      * differentiation order, to calculate the quadrature degree in \ref getDegree.
      **/
-    ExpressionOperatorBase(Expr const& expr, int order)
-      : expr_(expr)
+    GridFunctionOperatorBase(GridFunction const& gridFct, int order)
+      : gridFct_(gridFct)
+      , localFct_(localFunction(gridFct_))
       , order_(order)
     {}
 
@@ -28,13 +36,12 @@ namespace AMDiS
      * Since all operators need geometry information, the `Element::Geometry`
      * object `geometry` is created once, during grid traversal, and passed in.
      *
-     * By default, it binds the \ref expr_ to the `element`.
+     * By default, it binds the \ref localFct_ to the `element`.
      **/
-    template <class Element, class Geometry>
     void bind(Element const& element, Geometry const& geometry)
     {
       if (!bound_) {
-        expr_.bind(element);
+        localFct_.bind(element);
         isSimplex_ = geometry.type().isSimplex();
         isAffine_ = geometry.affine();
         bound_ = true;
@@ -45,32 +52,14 @@ namespace AMDiS
     void unbind()
     {
       bound_ = false;
-      initSize_ = 0;
-      expr_.unbind();
-    }
-
-    /// \brief Evaluate expression at quadrature points.
-    /**
-     * Initialization at all points `quad[iq].position()` where
-     * `iq` in `[0, quad.size())`.
-     *
-     * \param context Container for localContext, localGeometry and geometry, see \ref LocalAssembler
-     * \param quad Liste of evaluations points. Follows interface of \ref Dune::QuadratureRule
-     **/
-    template <class Context, class Quad>
-    void init(Context const& context, Quad const& quad)
-    {
-      if (initSize_ != quad.size()) {
-        expr_.init(context, quad);
-        initSize_ = quad.size();
-      }
+      localFct_.unbind();
     }
 
     /// Return expressions iq'th value. Must be initialized in \ref init before.
-    auto eval(std::size_t iq) const
+    auto coefficient(LocalCoordinate const& local) const
     {
-      assert( initSize_ > 0u );
-      return expr_[iq];
+      assert( bound_ );
+      return localFct_(local);
     }
 
 
@@ -86,9 +75,9 @@ namespace AMDiS
       assert( bound_ );
 
       int psiDegree = getPolynomialDegree(node);
-      int termDegree = expr_.getDegree();
+      int coeffDegree = order(localFct_);
 
-      int degree = psiDegree + termDegree;
+      int degree = psiDegree + coeffDegree;
       if (isSimplex_)
         degree -= order_;
       if (isAffine_)
@@ -111,9 +100,9 @@ namespace AMDiS
 
       int psiDegree = getPolynomialDegree(rowNode);
       int phiDegree = getPolynomialDegree(colNode);
-      int termDegree = expr_.getDegree();
+      int coeffDegree = order(localFct_);
 
-      int degree = psiDegree + phiDegree + termDegree;
+      int degree = psiDegree + phiDegree + coeffDegree;
       if (isSimplex_)
         degree -= order_;
       if (isAffine_)
@@ -122,45 +111,31 @@ namespace AMDiS
       return degree;
     }
 
-
-  protected:
-
-    Expr& expr()
-    {
-      return expr_;
-    }
-
-    Expr const& expr() const
-    {
-      return expr_;
-    }
-
-
   private:
+    GridFunction gridFct_;
+    LocalFunction localFct_;
 
-    Expr expr_;   //< the stores expression, see \ref ExpressionBase
     int order_;   //< the derivative order of this operator
 
     bool isSimplex_ = false;    //< the bound element is a simplex
     bool isAffine_ = false;     //< the bound geometry is affine
     bool bound_ = false;        //< an element is bound to the operator
-    std::size_t initSize_ = 0;  //< number of quadrature points, the expression was evaluated at.
   };
 
 
-  /// A default implementation of an ExpressionOperator if no specialization is available.
+  /// A default implementation of an GridFunctionOperator if no specialization is available.
   /**
    * An operator must implement either \ref calculateElementVector, or
    * \ref calculateElementMatrix, if it is a vector or matrix operator,
    * respectively.
    **/
-  template <class Tag, class Expr>
-  class ExpressionOperator
-      : public ExpressionOperatorBase<Expr>
+  template <class Tag, class GridFct>
+  class GridFunctionOperator
+      : public GridFunctionOperatorBase<GridFct>
   {
   public:
-    ExpressionOperator(Tag, Expr const& expr)
-      : ExpressionOperatorBase<Expr>(expr, 0)
+    GridFunctionOperator(Tag, GridFct const& gridFct)
+      : GridFunctionOperatorBase<GridFct>(gridFct, 0)
     {}
 
     /// Assemble a local element vector on the element that is bound.
@@ -191,21 +166,50 @@ namespace AMDiS
     }
   };
 
+  template <class Tag, class Expr>
+  struct ExpressionPreOperator
+  {
+    Tag tag;
+    Expr expr;
+  };
+
+  /// Store tag and expression in struct
+  template <class Tag, class Expr>
+  auto makeOperator(Tag t, Expr const& expr)
+  {
+    return ExpressionPreOperator<Tag, Expr>{t, expr};
+  }
+
+  /// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
+  template <class Tag, class Expr, class GridView>
+  auto makeGridOperator(ExpressionPreOperator<Tag, Expr> const& op, GridView const& gridView)
+  {
+    auto gridFct = makeGridFunction(op.expr, gridView, Dune::PriorityTag<42>{});
+    return GridFunctionOperator<Tag, decltype(gridFct)>{op.tag, gridFct};
+  }
+
+  template <class Tag, class Expr, class GridView>
+  auto makeGridOperator(std::reference_wrapper<ExpressionPreOperator<Tag, Expr>> op, GridView const& gridView)
+  {
+    ExpressionPreOperator<Tag, Expr> const& op_ref = op;
+    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView, Dune::PriorityTag<42>{});
+    return GridFunctionOperator<Tag, decltype(gridFct)>{op_ref.tag, gridFct};
+  }
 
-  /// Generate an \ref ExpressionOperator with a given tag `t` and pre-expression `epxr`.
-  template <class Tag, class PreExpr>
-  auto makeOperator(Tag t, PreExpr const& expr)
+  /// Generate a shared_ptr to \ref GridFunctionOperator from a PreOperator (tag, expr).
+  template <class Tag, class Expr, class GridView>
+  auto makeGridOperatorPtr(ExpressionPreOperator<Tag, Expr> const& op, GridView const& gridView)
   {
-    using Expr = ToTerm_t<PreExpr>;
-    return ExpressionOperator<Tag, Expr>{t, toTerm(expr)};
+    auto gridFct = makeGridFunction(op.expr, gridView, Dune::PriorityTag<42>{});
+    return std::make_shared<GridFunctionOperator<Tag, decltype(gridFct)>>(op.tag, gridFct);
   }
 
-  /// Generate a shared_ptr to \ref ExpressionOperator with a given tag `t` and pre-expression `epxr`.
-  template <class Tag, class PreExpr>
-  auto makeOperatorPtr(Tag t, PreExpr const& expr)
+  template <class Tag, class Expr, class GridView>
+  auto makeGridOperatorPtr(std::reference_wrapper<ExpressionPreOperator<Tag, Expr>> op, GridView const& gridView)
   {
-    using Expr = ToTerm_t<PreExpr>;
-    return std::make_shared<ExpressionOperator<Tag, Expr>>(t, toTerm(expr));
+    ExpressionPreOperator<Tag, Expr> const& op_ref = op;
+    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView, Dune::PriorityTag<42>{});
+    return std::make_shared<GridFunctionOperator<Tag, decltype(gridFct)>>(op_ref.tag, gridFct);
   }
 
 } // end namespace AMDiS

dune/amdis/GridFunctions.hpp

+#pragma once
+
+#include <dune/amdis/gridfunctions/AnalyticGridFunction.hpp>
+#include <dune/amdis/gridfunctions/ConstantGridFunction.hpp>
+#include <dune/amdis/gridfunctions/CoordsGridFunction.hpp>
+#include <dune/amdis/gridfunctions/DerivativeGridFunction.hpp>
+#include <dune/amdis/gridfunctions/FunctorGridFunction.hpp>
+#include <dune/amdis/gridfunctions/OperationsGridFunction.hpp>

dune/amdis/Initfile.hpp

@@ -16,7 +16,7 @@
 #include <dune/common/fvector.hh>
 
 #include <dune/amdis/Output.hpp>
-#include <dune/amdis/Math.hpp>
+#include <dune/amdis/common/Math.hpp>
 
 namespace AMDiS
 {

dune/amdis/LocalAssembler.hpp

@@ -189,7 +189,7 @@ namespace AMDiS
 
   private:
 
-    std::unique_ptr<Operator> storage_;     //< the stored operator, implementing \ref ExpressionOperatorBase
+    std::unique_ptr<Operator> storage_;     //< the stored operator, implementing \ref GridFunctionOperatorBase
     Operator& op_;
 
     Element const* element_ = nullptr;

dune/amdis/Operations.hpp

+#pragma once
+
+#include <dune/amdis/operations/Arithmetic.hpp>
+#include <dune/amdis/operations/Basic.hpp>
+#include <dune/amdis/operations/CMath.hpp>
+#include <dune/amdis/operations/Composer.hpp>
+#include <dune/amdis/operations/FieldMatVec.hpp>
+#include <dune/amdis/operations/MaxMin.hpp>

dune/amdis/Operators.hpp

 #pragma once
 
 #include <dune/amdis/LocalAssembler.hpp>
-#include <dune/amdis/OperatorBase.hpp>
+#include <dune/amdis/GridFunctionOperator.hpp>
 
 /**
  * In the following comments we use the notation
@@ -15,25 +15,25 @@
  */
 
 // zero-order operators
-#include <dune/amdis/operators/ZeroOrderTest.hpp>             // <psi * c>
-#include <dune/amdis/operators/ZeroOrderTestTrial.hpp>        // <psi, c * phi>
-#include <dune/amdis/operators/ZeroOrderTestTrialvec.hpp>     // <psi, b * Phi>
-#include <dune/amdis/operators/ZeroOrderTestvec.hpp>          // <Psi * b>
-#include <dune/amdis/operators/ZeroOrderTestvecTrial.hpp>     // <Psi, b * phi>
-#include <dune/amdis/operators/ZeroOrderTestvecTrialvec.hpp>  // <Psi, A * Phi>, <Psi, c * Phi>
+#include <dune/amdis/assembler/ZeroOrderTest.hpp>             // <psi * c>
+#include <dune/amdis/assembler/ZeroOrderTestTrial.hpp>        // <psi, c * phi>
+#include <dune/amdis/assembler/ZeroOrderTestTrialvec.hpp>     // <psi, b * Phi>
+#include <dune/amdis/assembler/ZeroOrderTestvec.hpp>          // <Psi * b>
+#include <dune/amdis/assembler/ZeroOrderTestvecTrial.hpp>     // <Psi, b * phi>
+#include <dune/amdis/assembler/ZeroOrderTestvecTrialvec.hpp>  // <Psi, A * Phi>, <Psi, c * Phi>
 
 // first-order operators
-#include <dune/amdis/operators/FirstOrderDivTestvecTrial.hpp>   // <div(Psi), c * phi>
-#include <dune/amdis/operators/FirstOrderGradTestTrial.hpp>     // <grad(psi), b * phi>
-#include <dune/amdis/operators/FirstOrderGradTestTrialvec.hpp>  // <grad(psi), c * Phi>
-#include <dune/amdis/operators/FirstOrderPartialTestTrial.hpp>  // <d_i(psi), c * psi>
-#include <dune/amdis/operators/FirstOrderTestDivTrialvec.hpp>   // <psi, c * div(Phi)>
-#include <dune/amdis/operators/FirstOrderTestGradTrial.hpp>     // <psi, b * grad(phi)>
-#include <dune/amdis/operators/FirstOrderTestPartialTrial.hpp>  // <psi, c * d_i(phi)>
-#include <dune/amdis/operators/FirstOrderTestvecGradTrial.hpp>  // <Psi, c * grad(phi)>
+#include <dune/amdis/assembler/FirstOrderDivTestvecTrial.hpp>   // <div(Psi), c * phi>
+#include <dune/amdis/assembler/FirstOrderGradTestTrial.hpp>     // <grad(psi), b * phi>
+#include <dune/amdis/assembler/FirstOrderGradTestTrialvec.hpp>  // <grad(psi), c * Phi>
+#include <dune/amdis/assembler/FirstOrderPartialTestTrial.hpp>  // <d_i(psi), c * psi>
+#include <dune/amdis/assembler/FirstOrderTestDivTrialvec.hpp>   // <psi, c * div(Phi)>
+#include <dune/amdis/assembler/FirstOrderTestGradTrial.hpp>     // <psi, b * grad(phi)>
+#include <dune/amdis/assembler/FirstOrderTestPartialTrial.hpp>  // <psi, c * d_i(phi)>
+#include <dune/amdis/assembler/FirstOrderTestvecGradTrial.hpp>  // <Psi, c * grad(phi)>
 
 // second-order operators
-#include <dune/amdis/operators/SecondOrderDivTestvecDivTrialvec.hpp>    // <div(Psi), c * div(Phi)>
-#include <dune/amdis/operators/SecondOrderGradTestGradTrial.hpp>        // <grad(psi), A * grad(phi)>
-#include <dune/amdis/operators/SecondOrderPartialTestPartialTrial.hpp>  // <d_i(psi), c * d_j(phi)>
+#include <dune/amdis/assembler/SecondOrderDivTestvecDivTrialvec.hpp>    // <div(Psi), c * div(Phi)>
+#include <dune/amdis/assembler/SecondOrderGradTestGradTrial.hpp>        // <grad(psi), A * grad(phi)>
+#include <dune/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp>  // <d_i(psi), c * d_j(phi)>
 

dune/amdis/ProblemStat.hpp

@@ -31,6 +31,8 @@
 #include <dune/amdis/common/TypeDefs.hpp>
 #include <dune/amdis/common/Utility.hpp>
 
+#include <dune/amdis/gridfunctions/DOFVectorView.hpp>
+
 #include <dune/amdis/io/FileWriterInterface.hpp>
 
 #include <dune/amdis/utility/DiscreteFunction.hpp>
@@ -142,21 +144,6 @@ namespace AMDiS
                         RowTreePath row, ColTreePath col,
                         Values const& values);
 
-
-  private: // implementation detail
-
-    template <class Predicate, class RowNode, class ColNode>
-    void addDirichletBCImpl(Predicate const& predicate,
-                            RowNode const& rowNode, ColNode const& colNode,
-                            RangeType<RowNode> const& value);
-
-    template <class Predicate, class RowNode, class ColNode, class Values>
-      std::enable_if_t<not std::is_convertible<Values, RangeType<RowNode>>::value>
-    addDirichletBCImpl(Predicate const& predicate,
-                       RowNode const& rowNode, ColNode const& colNode,
-                       Values const& values);
-
-
   public:
 
     /// Implementation of \ref ProblemStatBase::solve
@@ -233,7 +220,7 @@ namespace AMDiS
     auto getSolution(TreePath path) const
     {
       auto tp = makeTreePath(path);
-      return makeDiscreteFunctionPtr(globalBasis, tp, solution);
+      return makeDOFVectorView(*solution, tp);
     }
 
 

dune/amdis/ProblemStat.inc.hpp

@@ -10,6 +10,7 @@
 #include <dune/amdis/Assembler.hpp>
 #include <dune/amdis/FileWriter.hpp>
 #include <dune/amdis/LocalAssembler.hpp>
+#include <dune/amdis/GridFunctionOperator.hpp>
 #include <dune/amdis/common/Loops.hpp>
 #include <dune/amdis/common/Timer.hpp>
 
@@ -132,13 +133,14 @@ void ProblemStat<Traits>::createFileWriter()
 template <class Traits>
   template <class Operator, class RowTreePath, class ColTreePath>
 void ProblemStat<Traits>::addMatrixOperator(
-    Operator const& op,
+    Operator const& preOp,
     RowTreePath row, ColTreePath col,
     double* factor, double* estFactor)
 {
   auto i = child(globalBasis->localView().tree(), makeTreePath(row));
   auto j = child(globalBasis->localView().tree(), makeTreePath(col));
 
+  auto op = makeGridOperator(preOp, globalBasis->gridView());
   auto localAssembler = makeLocalAssemblerPtr<Element>(op, i, j);
 
   matrixOperators[i][j].element.push_back({localAssembler, factor, estFactor});
@@ -150,7 +152,7 @@ template <class Traits>
   template <class Operator, class RowTreePath, class ColTreePath>
 void ProblemStat<Traits>::addMatrixOperator(
     BoundaryType b,
-    Operator const& op,
+    Operator const& preOp,
     RowTreePath row, ColTreePath col,
     double* factor, double* estFactor)
 {
@@ -158,6 +160,7 @@ void ProblemStat<Traits>::addMatrixOperator(
   auto j = child(globalBasis->localView().tree(), makeTreePath(col));
   using Intersection = typename GridView::Intersection;
 
+  auto op = makeGridOperator(preOp, globalBasis->gridView());
   auto localAssembler = makeLocalAssemblerPtr<Intersection>(op, i, j);
 
   matrixOperators[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
@@ -168,12 +171,13 @@ void ProblemStat<Traits>::addMatrixOperator(
 template <class Traits>
   template <class Operator, class TreePath>
 void ProblemStat<Traits>::addVectorOperator(
-    Operator const& op,
+    Operator const& preOp,
     TreePath path,
     double* factor, double* estFactor)
 {
   auto i = child(globalBasis->localView().tree(), makeTreePath(path));
 
+  auto op = makeGridOperator(preOp, globalBasis->gridView());
   auto localAssembler = makeLocalAssemblerPtr<Element>(op, i);
 
   rhsOperators[i].element.push_back({localAssembler, factor, estFactor});
@@ -185,13 +189,14 @@ template <class Traits>
   template <class Operator, class TreePath>
 void ProblemStat<Traits>::addVectorOperator(
     BoundaryType b,
-    Operator const& op,
+    Operator const& preOp,
     TreePath path,
     double* factor, double* estFactor)
 {
   auto i = child(globalBasis->localView().tree(), makeTreePath(path));
   using Intersection = typename GridView::Intersection;
 
+  auto op = makeGridOperator(preOp, globalBasis->gridView());
   auto localAssembler = makeLocalAssemblerPtr<Intersection>(op, i);
 
   rhsOperators[i].boundary.push_back({localAssembler, factor, estFactor, b});
@@ -211,33 +216,13 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val
   auto i = child(globalBasis->localView().tree(), makeTreePath(row));
   auto j = child(globalBasis->localView().tree(), makeTreePath(col));
 
-  addDirichletBCImpl(predicate, i, j, values);
-}
-
-template <class Traits>
-  template <class Predicate, class RowNode, class ColNode>
-void ProblemStat<Traits>::
-addDirichletBCImpl(Predicate const& predicate, RowNode const& row, ColNode const& col, RangeType<RowNode> const& value)
-{
-  using Range = RangeType<RowNode>;
-  using BcType = DirichletBC<WorldVector,Range>;
-  auto bc = std::make_shared<BcType>(predicate, [value](auto const&) { return value; });
-  constraints[row][col].push_back(bc);
-}
-
-template <class Traits>
-  template <class Predicate, class RowNode, class ColNode, class Values>
-    std::enable_if_t<not std::is_convertible<Values, RangeType<RowNode>>::value>
-ProblemStat<Traits>::
-addDirichletBCImpl(Predicate const& predicate, RowNode const& row, ColNode const& col, Values const& values)
-{
-  using Range = RangeType<RowNode>;