cleanup of source-tree, new directory common for commonly used functionality,...

cleanup of source-tree, new directory common for commonly used functionality, advanced addMatrixOperator syntax, OperatorTerm creation from Operators.

dune/amdis/CMakeLists.txt

@@ -62,37 +62,73 @@ install(FILES
     AdaptInstationary.hpp
     AdaptStationary.hpp
     AMDiS.hpp
-    Basic.hpp
+    CreatorInterface.hpp
+    CreatorMap.hpp
     DirichletBC.hpp
     DirichletBC.inc.hpp
+    FileWriter.hpp
     Flag.hpp
-    IndexSeq.hpp
     Initfile.hpp
+    LinearAlgebra.hpp
     Log.hpp
-    Loop.hpp
     Math.hpp
+    Mesh.hpp
+    OperatorEvaluation.hpp
     Operator.hpp
     Operator.inc.hpp
+    OperatorTermBase.hpp
     OperatorTerm.hpp
     ProblemInstatBase.hpp
     ProblemInstat.hpp
     ProblemInstat.inc.hpp
     ProblemInterationInterface.hpp
+    ProblemStatBase.hpp
     ProblemStat.hpp
     ProblemStat.inc.hpp
-    ProblemStatBase.hpp
     ProblemTimeInterface.hpp
     StandardProblemIteration.hpp
-    Timer.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)

dune/amdis/CreatorInterface.hpp

 #pragma once
 
+#include <memory>
 #include <string>
 
-#include <dune/amdis/Basic.hpp>
+#include <dune/amdis/Log.hpp>
 
 namespace AMDiS
 {
@@ -57,7 +58,14 @@ namespace AMDiS
     virtual std::shared_ptr<BaseClass> create(std::string) = 0;
   };
   
-  template <class T>
-  void addDefaultCreators(Type<T>);
+  
+  /// cast a ptr of CreatorInterface to CreatorInterfaceName
+  template <class BaseClass>
+  inline CreatorInterfaceName<BaseClass>* named(CreatorInterface<BaseClass>* ptr)
+  {
+    auto result = dynamic_cast<CreatorInterfaceName<BaseClass>*>(ptr);
+    AMDIS_TEST_EXIT_DBG(result, "Can not cast CreatorInterface to CreatorInterfaceName!");
+    return result;
+  }
 
 } // end namespace AMDiS

dune/amdis/CreatorMap.hpp

@@ -8,6 +8,13 @@
 
 namespace AMDiS
 {  
+  // forward declaration. 
+  // All classes that need creators must specialize this class and implement
+  // a static void init() method.
+  template <class BaseClass>
+  struct DefaultCreators;
+  
+  
   /** \ingroup Common
    * \brief
    * A CreatorMap is used to construct objects, which types depends on key words
@@ -54,7 +61,7 @@ namespace AMDiS
     {
       if (!initialized) {
         initialized = true;
-        addDefaultCreators(Type<BaseClass>{});
+        DefaultCreators<BaseClass>::init();
       }
     }
 

dune/amdis/DirichletBC.hpp

@@ -4,7 +4,7 @@
 #include <type_traits>
 #include <vector>
 
-#include <dune/functions/common/functionconcepts.hh>
+#include <dune/amdis/common/ConceptsBase.hpp>
 
 #include "Log.hpp"
 
@@ -15,8 +15,8 @@ namespace AMDiS
   {
   public:
     template <class Predicate, class Values,
-      class = std::enable_if_t< Dune::Functions::Concept::isFunction<Predicate, bool(WorldVector)>() &&
-				Dune::Functions::Concept::isFunction<Values,  double(WorldVector)>() > >
+      class = std::enable_if_t< concepts::Functor<Predicate, bool(WorldVector)> &&
+				concepts::Functor<Values,  double(WorldVector)> > >
     DirichletBC(Predicate&& predicate, Values&& values)
       : predicate(std::forward<Predicate>(predicate))
       , values(std::forward<Values>(values))

dune/amdis/FileWriter.hpp

@@ -9,8 +9,8 @@
 #include <dune/grid/io/file/vtk/vtksequencewriter.hh>
 #include <dune/geometry/genericgeometry/referenceelements.hh>
 
-#include "Loops.hpp"
-#include "Initfile.hpp"
+#include <dune/amdis/Initfile.hpp>
+#include <dune/amdis/common/Loops.hpp>
 
 namespace AMDiS
 {
@@ -52,7 +52,7 @@ namespace AMDiS
     
     /// default write method for time-depended data
     template <class SystemVectorType>
-    void write(double time, SystemVectorType&& solutions)
+    void write(double time, SystemVectorType const& solutions)
     {
       vtkWriter->clear();
       // copy dofvector to vertex data
@@ -67,7 +67,7 @@ namespace AMDiS
 
     /// default write method for stationary data
     template <class SystemVectorType>
-    void write(SystemVectorType&& solutions)
+    void write(SystemVectorType const& solutions)
     {
       vtkWriter->clear();
       // copy dofvector to vertex data

dune/amdis/LinearAlgebra.hpp

 #pragma once
 
-#include "linear_algebra/LinearSolverInterface.hpp"
-#include "linear_algebra/SolverInfo.hpp"
+#include <dune/amdis/linear_algebra/LinearSolverInterface.hpp>
+#include <dune/amdis/linear_algebra/SolverInfo.hpp>
 
 #if defined(AMDIS_BACKEND_ISTL)
 
-#include "linear_algebra/istl/SystemVector.hpp"
-#include "linear_algebra/istl/SystemMatrix.hpp"
-#include "linear_algebra/istl/LinearSolver.hpp"
+#include <dune/amdis/linear_algebra/istl/SystemVector.hpp>
+#include <dune/amdis/linear_algebra/istl/SystemMatrix.hpp>
+#include <dune/amdis/linear_algebra/istl/LinearSolver.hpp>
 
 #elif defined(AMDIS_BACKEND_MTL)
 
-#include "linear_algebra/mtl/SystemVector.hpp"
-#include "linear_algebra/mtl/SystemMatrix.hpp"
-#include "linear_algebra/mtl/LinearSolver.hpp"
-#include "linear_algebra/mtl/ITL_Solver.hpp"
-#include "linear_algebra/mtl/BITL_Solver.hpp"
+#include <dune/amdis/linear_algebra/mtl/SystemVector.hpp>
+#include <dune/amdis/linear_algebra/mtl/SystemMatrix.hpp>
+#include <dune/amdis/linear_algebra/mtl/LinearSolver.hpp>
+#include <dune/amdis/linear_algebra/mtl/ITL_Solver.hpp>
+#include <dune/amdis/linear_algebra/mtl/BITL_Solver.hpp>
 
 #elif defined(AMDIS_BACKEND_PETSC)
 
-#include "linear_algebra/petsc/SystemVector.hpp"
-#include "linear_algebra/petsc/SystemMatrix.hpp"
-#include "linear_algebra/petsc/LinearSolver.hpp"
+#include <dune/amdis/linear_algebra/petsc/SystemVector.hpp>
+#include <dune/amdis/linear_algebra/petsc/SystemMatrix.hpp>
+#include <dune/amdis/linear_algebra/petsc/LinearSolver.hpp>
 
 #else
 

dune/amdis/Operator.hpp

@@ -29,6 +29,14 @@ namespace AMDiS
       return addZOTImpl(toTerm(c));
     }
     
+    template <class... Args>
+    static std::shared_ptr<Self> zot(Args&&... args)
+    {
+      auto op = std::make_shared<Self>();
+      op->addZOT(std::forward<Args>(args)...);
+      return op;
+    }
+    
     
     /// Add coefficients for first-order operator < psi, term * grad(phi) > or 
     /// < grad(psi), term * phi >, where the first corresponds to 
@@ -50,6 +58,14 @@ namespace AMDiS
       return addFOTImpl(toTerm(b), _i, firstOrderType);
     }
     
+    template <class... Args>
+    static std::shared_ptr<Self> fot(Args&&... args)
+    {
+      auto op = std::make_shared<Self>();
+      op->addFOT(std::forward<Args>(args)...);
+      return op;
+    }
+    
     
     /// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
     /// where \p A can be a matrix or a scalar expression. 
@@ -69,10 +85,21 @@ namespace AMDiS
       return addSOTImpl(toTerm(A), _i, _j);
     }
     
+    template <class... Args>
+    static std::shared_ptr<Self> sot(Args&&... args)
+    {
+      auto op = std::make_shared<Self>();
+      op->addSOT(std::forward<Args>(args)...);
+      return op;
+    }
     
-    static std::shared_ptr<Self> create()
+    
+    /// Calls \ref zot(), \ref for() or \ref sot(), depending on template 
+    /// parameter \p Order.
+    template <size_t Order, class... Args>
+    static std::shared_ptr<Self> create(Args&&... args)
     {
-      return std::make_shared<Self>();
+      return create(index_<Order>{}, std::forward<Args>(args)...);
     }
     
     
@@ -152,6 +179,25 @@ namespace AMDiS
     template <class Term, size_t I, size_t J>
     Self& addSOTImpl(Term const& term, const index_<I>, const index_<J>);
     
+    template <class... Args>
+    static std::shared_ptr<Self> create(index_<0>, Args&&... args)
+    {
+      return zot(std::forward<Args>(args)...);
+    }
+    
+    template <class... Args>
+    static std::shared_ptr<Self> create(index_<1>, Args&&... args)
+    {
+      return fot(std::forward<Args>(args)...);
+    }
+    
+    template <class... Args>
+    static std::shared_ptr<Self> create(index_<2>, Args&&... args)
+    {
+      return sot(std::forward<Args>(args)...);
+    }
+      
+      
   private:
     /// List of all second order terms
     std::list<OperatorTermType*> secondOrder;

dune/amdis/OperatorEvaluation.hpp

+#pragma once
+
+#include <cassert>
+
+#include <dune/amdis/common/ValueCategory.hpp>
+
+namespace AMDiS {
+
+
+template <class K, int SIZE>
+K sum(Dune::FieldVector<K, SIZE> const& vector)
+{
+  K result = 0;
+  for (int i = 0; i < SIZE; ++i)
+    result += vector[i];
+  return result;
+}
+    
+    
+
+class OperatorEvaluation
+{
+  using Scalar = Dune::FieldVector<double, 1>;
+    
+protected:
+  template <class... Args>
+  auto evalZotImpl(Args...) const { assert( false ); return 0; }
+      
+  template <class ValueType>
+  auto evalZotImpl(tag::scalar, tag::none, ValueType const& v, 
+                    Scalar const& test, Scalar const& trial) const
+  {
+    return v * test * trial;
+  }
+  
+  template <class... Args>
+  auto evalFotImpl(Args...) const { assert( false ); return 0; }
+      
+  template <size_t I, class ValueType, class Gradient>
+  auto evalFotImpl(tag::scalar, VectorComponent<I>, ValueType const& v, 
+                    Gradient const& grad_test, Scalar const& trial) const
+  {
+    return v * grad_test[I] * trial;
+  }
+      
+  template <class ValueType, class Gradient>
+  auto evalFotImpl(tag::vector, tag::none, ValueType const& v, 
+                    Gradient const& grad_test, Scalar const& trial) const
+  {
+    return (v * grad_test) * trial;
+  }
+      
+  template <class ValueType, class Gradient>
+  auto evalFotImpl(tag::scalar, tag::none, ValueType const& v, 
+                    Gradient const& grad_test, Scalar const& trial) const
+  {
+    return v * sum(grad_test) * trial;
+  }
+  
+  
+  template <class... Args>
+  auto evalSotImpl(Args...) const { assert( false ); return 0; }
+      
+  template <size_t R, size_t C, class ValueType, class Gradient>
+  auto evalSotImpl(tag::scalar, MatrixComponent<R, C>, ValueType const& v, 
+                    Gradient const& grad_test, Gradient const& grad_trial) const
+  {
+    return v * (grad_test[R] * grad_trial[C]);
+  }
+      
+  template <class ValueType, class Gradient>
+  auto evalSotImpl(tag::matrix, tag::none, ValueType const& M, 
+                    Gradient const& grad_test, Gradient const& grad_trial) const
+  {
+    return grad_test * (M * grad_trial);
+  }
+      
+  template <class ValueType, class Gradient>
+  auto evalSotImpl(tag::scalar, tag::none, ValueType const& v, 
+                    Gradient const& grad_test, Gradient const& grad_trial) const
+  {
+    return v * (grad_test * grad_trial);
+  }
+};
+  
+
+} // end namespace AMDiS
\ No newline at end of file

dune/amdis/OperatorTerm.hpp

@@ -4,125 +4,15 @@
 #include <vector>
 #include <type_traits>
 
-#include <dune/geometry/quadraturerules.hh>
 #include <dune/istl/bvector.hh>
 #include <dune/functions/common/functionconcepts.hh>
 #include <dune/functions/functionspacebases/pqknodalbasis.hh>
 
-#include "OperatorTermBase.hpp"
+#include <dune/amdis/OperatorTermBase.hpp>
+#include <dune/amdis/common/Traits.hpp>
 
 namespace AMDiS 
 {
-    
-  template <class MeshView>
-  class OperatorTerm
-  {
-  protected:
-    using Codim0  = typename MeshView::template Codim<0>;
-    using Element = typename Codim0::Entity;
-    
-    static constexpr int dim = Element::dimension;
-    
-    using QuadratureRule = Dune::QuadratureRule<double, dim>;
-    using PointList = std::vector<Dune::QuadraturePoint<double, dim>>;
-      
-  public:
-    virtual void init(Element const& element, 
-                      PointList const& points) = 0;
-                      
-    virtual double evalZot(size_t iq, 
-        Dune::FieldVector<double,1> const& test, 
-        Dune::FieldVector<double,1> const trial = 1.0) const = 0;
-        
-    virtual double evalFot1(size_t iq, 
-        Dune::FieldVector<double,1> const& test, 
-        Dune::FieldVector<double,dim> const& grad_trial) const = 0;
-        
-    virtual double evalFot2(size_t iq, 
-        Dune::FieldVector<double,dim> const& grad_test, 
-        Dune::FieldVector<double,1> const trial = 1.0) const = 0;
-        
-    virtual double evalSot(size_t iq, 
-        Dune::FieldVector<double,dim> const& grad_test, 
-        Dune::FieldVector<double,dim> const& grad_trial) const = 0;
-                          
-    virtual int getDegree() const = 0;
-  };
-  
-  
-  
-  /// Base class for all operators based on expressions
-  template <class MeshView, class Term, class Traits = _none>
-  class GenericOperatorTerm 
-      : public OperatorTerm<MeshView>
-      , public OperatorEvaluation
-  {
-    using Super   = OperatorTerm<MeshView>;
-    using Element = typename Super::Element;
-    using PointList = typename Super::PointList;
-    
-    static constexpr int dim = Element::dimension;
-
-  public:
-    GenericOperatorTerm(Term const& term)
-      : term(term)
-    {}
-
-    virtual void init(Element const& element, 
-                      PointList const& points) override
-    {
-      term.init(element, points);
-      
-      // cache term evaluation
-      values.resize(points.size());
-      for (size_t iq = 0; iq < points.size(); ++iq)
-        values[iq] = term[iq];
-    }
-                      
-    virtual double evalZot(size_t iq, 
-        Dune::FieldVector<double,1> const& test, 
-        Dune::FieldVector<double,1> const trial = 1.0) const override
-    {
-      return this->evalZotImpl(_cat{}, _traits{}, values[iq], test, trial);
-    }
-        
-    virtual double evalFot1(size_t iq, 
-        Dune::FieldVector<double,1> const& test, 
-        Dune::FieldVector<double,dim> const& grad_trial) const override
-    {
-      return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_trial, test);
-    }
-        
-    virtual double evalFot2(size_t iq, 
-        Dune::FieldVector<double,dim> const& grad_test, 
-        Dune::FieldVector<double,1> const trial = 1.0) const override
-    {
-      return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_test, trial);
-    }
-        
-    virtual double evalSot(size_t iq, 
-        Dune::FieldVector<double,dim> const& grad_test, 
-        Dune::FieldVector<double,dim> const& grad_trial) const override
-    {
-      return this->evalSotImpl(_cat{}, _traits{}, values[iq], grad_test, grad_trial);
-    }
-                            
-    virtual int getDegree() const override
-    {
-      return term.getDegree();
-    }
-
-  private:
-    Term term;
-    
-    using value_type = std::decay_t< decltype( std::declval<Term>()[std::declval<size_t>()] ) >;
-    using _cat    = ValueCategory_t<value_type>;
-    using _traits = Traits;
-    
-    std::vector<value_type> values;
-  };
-  
-
 // some example terms
 // -----------------------------------------------------------------------------
 
@@ -158,7 +48,7 @@ namespace AMDiS
     
   /// generator function for \ref ConstantTerm expressions
   template <class T>
-  std::enable_if_t< std::is_arithmetic<T>::value, ConstantTerm<T> >
+  std::enable_if_t< concepts::Arithmetic<T>, ConstantTerm<T> >
   constant(T value) { return {value}; }
   
   
@@ -172,7 +62,7 @@ namespace AMDiS
       
     template <class F, int dow>
     static constexpr bool isCallableDow()
-    { return Dune::Functions::Concept::isCallable<F, Dune::FieldVector<double, dow>>(); }
+    { return concepts::Callable<F(Dune::FieldVector<double, dow>)>; }
       
     template <class F>
     using CoordsFunctorResult_t = typename std::conditional_t<
@@ -192,7 +82,7 @@ namespace AMDiS
     using value_type = Impl::CoordsFunctorResult_t<Functor>;
     
     template <class F,
-      class = std::enable_if_t<std::is_same<Functor, std::decay_t<F>>::value> >
+      class = std::enable_if_t<concepts::Compatible<Functor, F>> >
     CoordsTerm(F&& f, int degree = 1)
       : fct(std::forward<F>(f))
       , degree(degree)
@@ -246,7 +136,7 @@ namespace AMDiS
   class DOFVectorTerm
   {      
   public:
-    using value_type = typename DOFVectorType::value_type;
+    using value_type = Value_t<DOFVectorType>;
     
     DOFVectorTerm(DOFVectorType const& dofvector, double factor = 1.0)
       : vector(dofvector.getVector())
@@ -326,7 +216,7 @@ namespace AMDiS
   class DOFVectorFuncTerm
   {      
   public:
-    using value_type = typename DOFVectorType::value_type;
+    using value_type = Value_t<DOFVectorType>;
     
     template <class F_>
     DOFVectorFuncTerm(DOFVectorType const& dofvector, F_&& f, int f_deg)
@@ -414,7 +304,7 @@ namespace AMDiS
     
   public:    
     using value_type = Jacobian;
-    using field_type = typename DOFVectorType::value_type;
+    using field_type = Value_t<DOFVectorType>;
     
     GradientTerm(DOFVectorType const& dofvector, double factor = 1.0)
       : vector(dofvector.getVector())
@@ -476,7 +366,7 @@ namespace AMDiS
   };
 
   
-  /// generator function for \ref DOFVector expressions
+  /// generator function for an expression of gradients of \ref DOFVectors
   template <class DOFVectorType>
   GradientTerm<std::decay_t<DOFVectorType>>
   gradientOf(DOFVectorType&& vector, double factor = 1.0)
@@ -496,7 +386,7 @@ namespace AMDiS
     
   public:
     using value_type = Derivative;
-    using field_type = typename DOFVectorType::value_type;
+    using field_type = Value_t<DOFVectorType>;
     
     DerivativeTerm(DOFVectorType const& dofvector, int component, double factor = 1.0)
       : vector(dofvector.getVector())
@@ -566,7 +456,7 @@ namespace AMDiS
   };
 
   
-  /// generator function for \ref DOFVector expressions
+  /// generator function for an expression of partial derivatives of \ref DOFVectors
   template <class DOFVectorType>
   DerivativeTerm<std::decay_t<DOFVectorType>>
   derivativeOf(DOFVectorType&& vector, int component, double factor = 1.0)

dune/amdis/OperatorTermBase.hpp

 #pragma once
 
-#include "Traits.hpp"
+#include <vector>
+#include <type_traits>
 
-namespace AMDiS {
+#include <dune/geometry/quadraturerules.hh>
 
+#include <dune/amdis/OperatorEvaluation.hpp>
+#include <dune/amdis/common/ValueCategory.hpp>
 
-template <class K, int SIZE>
-K sum(Dune::FieldVector<K, SIZE> const& vector)
+namespace AMDiS 
 {
-  K result = 0;
-  for (int i = 0; i < SIZE; ++i)
-    result += vector[i];
-  return result;
-}