composible grid functions

dune/amdis/common/Mpl.hpp

@@ -110,15 +110,24 @@ namespace AMDiS
   template <bool... Bs>
   using all_of_t = std::is_same<Impl::all_helper<true, Bs...>, Impl::all_helper<Bs..., true>>;
 
+  template <bool... Bs>
+  constexpr bool all_of_v = all_of_t<Bs...>::value;
+
   template <bool... Bs>
   using and_t = all_of_t<Bs...>;
 
   template <bool... Bs>
   using none_of_t = std::is_same<Impl::all_helper<false, Bs...>, Impl::all_helper<Bs..., false>>;
 
+  template <bool... Bs>
+  constexpr bool none_of_v = none_of_t<Bs...>::value;
+
   template <bool... Bs>
   using any_of_t = bool_t<not none_of_t<Bs...>::value>;
 
+  template <bool... Bs>
+  constexpr bool any_of_v = any_of_t<Bs...>::value;
+
   template <bool... Bs>
   using or_t = any_of_t<Bs...>;
 

dune/amdis/common/Operations.hpp

@@ -32,6 +32,11 @@ namespace AMDiS
       {
         return lhs + rhs;
       }
+
+      friend int order(plus, int lhs, int rhs)
+      {
+        return std::max(lhs, rhs);
+      }
     };
 
     struct plus_assign
@@ -51,6 +56,11 @@ namespace AMDiS
       {
         return lhs - rhs;
       }
+
+      friend int order(minus, int lhs, int rhs)
+      {
+        return std::max(lhs, rhs);
+      }
     };
 
     struct minus_assign
@@ -70,6 +80,11 @@ namespace AMDiS
       {
         return lhs * rhs;
       }
+
+      friend int order(times, int lhs, int rhs)
+      {
+        return lhs + rhs;
+      }
     };
 
     struct times_assign
@@ -89,6 +104,12 @@ namespace AMDiS
       {
         return lhs / rhs;
       }
+
+      // NOTE: the order is given only approximately
+      friend int order(divides, int lhs, int rhs)
+      {
+        return lhs + rhs + 1;
+      }
     };
 
     struct divides_assign
@@ -108,6 +129,11 @@ namespace AMDiS
       {
         return -x;
       }
+
+      friend int order(divides, int x)
+      {
+        return x;
+      }
     };
 
     /// Operation that represents max(A,B)
@@ -118,6 +144,12 @@ namespace AMDiS
       {
         return Math::max(lhs, rhs);
       }
+
+      // NOTE: the order is given only approximately
+      friend int order(maximum, int lhs, int rhs)
+      {
+        return std::max(lhs, rhs);
+      }
     };
 
     /// Operation that represents min(A,B)
@@ -128,6 +160,12 @@ namespace AMDiS
       {
         return Math::min(lhs, rhs);
       }
+
+      // NOTE: the order is given only approximately
+      friend int order(minimum, int lhs, int rhs)
+      {
+        return std::max(lhs, rhs);
+      }
     };
 
     /// Operation that represents max(|A|,|B|)
@@ -138,6 +176,12 @@ namespace AMDiS
       {
         return Math::max(std::abs(lhs), std::abs(rhs));
       }
+
+      // NOTE: the order is given only approximately
+      friend int order(abs_max, int lhs, int rhs)
+      {
+        return std::max(lhs, rhs);
+      }
     };
 
     /// Operation that represents min(|A|,|B|)
@@ -148,6 +192,12 @@ namespace AMDiS
       {
         return Math::min(std::abs(lhs), std::abs(rhs));
       }
+
+      // NOTE: the order is given only approximately
+      friend int order(abs_min, int lhs, int rhs)
+      {
+        return std::max(lhs, rhs);
+      }
     };
 
     /// Operation that represents A*factor
@@ -166,6 +216,11 @@ namespace AMDiS
         return lhs * scalar;
       }
 
+      friend int order(times_scalar, int lhs)
+      {
+        return lhs;
+      }
+
     private:
       S scalar;
     };
@@ -186,6 +241,11 @@ namespace AMDiS
         return lhs / scalar;
       }
 
+      friend int order(divides_scalar, int lhs)
+      {
+        return lhs;
+      }
+
     private:
       S scalar;
     };
@@ -199,6 +259,11 @@ namespace AMDiS
       {
         return std::conj(x);
       }
+
+      friend int order(hermitian, int x)
+      {
+        return x;
+      }
     };
 
     template <class T>
@@ -208,6 +273,11 @@ namespace AMDiS
       {
         return x;
       }
+
+      friend int order(hermitian, int x)
+      {
+        return x;
+      }
     };
 
     /** @} **/

dune/amdis/gridfunctions/AnalyticGridFunction.hpp

+#pragma once
+
+#include <type_traits>
+
+#include <dune/common/std/optional.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+
+#include <dune/amdis/gridfunctions/GridFunctionConcepts.hpp>
+
+#define AMDIS_CALLABLE_DEFAULT_ORDER 1
+
+namespace AMDiS
+{
+  /**
+    * \addtogroup GridFunctions
+    * @{
+    **/
+
+
+  /// \brief A Gridfunction that applies a functor to the evaluated Gridfunctions
+  template <class Function, class GridView>
+  class AnalyticGridFunction
+  {
+  public:
+    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
+    using Element = typename EntitySet::Element;
+    using Geometry = typename Element::Geometry;
+
+    using Domain = typename Geometry::GlobalCoordinate;
+    using Range = typename std::result_of<Function(Domain)>::type;
+
+  public:
+
+    template <class Fct>
+    class LocalFunction
+    {
+    public:
+      using Range = typename std::result_of<Fct(Domain)>::type;
+      using Domain = typename Geometry::LocalCoordinate;
+
+    public:
+      LocalFunction(Fct const& fct)
+        : fct_(fct)
+      {}
+
+      void bind(Element const& element)
+      {
+        geometry_.emplace(element.geometry());
+      }
+
+      void unbind()
+      {
+        geometry_.reset();
+      }
+
+      Range operator()(Domain const& local) const
+      {
+        return fct_(geometry_.value().global(local));
+      }
+
+      template <class F,
+        REQUIRES(Concepts::HasDerivative<F>)>
+      friend auto derivative(LocalFunction<F> const& glf)
+      {
+        auto df = derivative(glf.fct_);
+        return LocalFunction<decltype(df)>{df};
+      }
+
+    private:
+      Fct fct_;
+      Dune::Std::optional<Geometry> geometry_;
+    };
+
+
+  public:
+
+    AnalyticGridFunction(Function const& fct, GridView const& gridView)
+      : fct_(fct)
+      , entitySet_(gridView)
+    {}
+
+    Range operator()(Domain const& x) const
+    {
+      return fct_(x);
+    }
+
+    template <class F, class GV>
+    friend auto localFunction(AnalyticGridFunction<F,GV> const& gf)
+      -> LocalFunction<F>
+    {
+      return LocalFunction<F>{gf.fct_};
+    }
+
+    EntitySet const& entitySet() const
+    {
+      return entitySet_;
+    }
+
+    Function const& fct() const { return fct_; }
+
+  private:
+    Function fct_;
+
+    EntitySet entitySet_;
+  };
+
+
+  template <class F, class GV>
+  auto derivative(AnalyticGridFunction<F,GV> const& gf,
+    REQUIRES(Concepts::HasDerivative<F>))
+  {
+    auto df = derivative(gf.fct());
+    return AnalyticGridFunction<decltype(df), GV>{df, gf.entitySet().gridView()};
+  }
+
+  template <class F, class GV>
+  int order(AnalyticGridFunction<F,GV> const& gf)
+  {
+    return Dune::Hybrid::ifElse( bool_<Concepts::HasOrder<F>>,
+      [&gf](auto id) { return id( order( gf.fct()) ); },
+      []   (auto id) { return AMDIS_CALLABLE_DEFAULT_ORDER; }
+    );
+  }
+
+
+  template <class Function>
+  struct AnalyticPreGridFunction
+  {
+    Function fct_;
+  };
+
+  namespace Concepts { namespace Definition
+  {
+    template <class Functor>
+    struct IsPreGridFunction<AnalyticPreGridFunction<Functor>>
+      : std::true_type {};
+  }}
+
+  template <class Function,
+    REQUIRES(Concepts::CallableDomain<Function>)>
+  auto eval_(Function const& f)
+  {
+    return AnalyticPreGridFunction<Function>{f};
+  }
+
+
+  /// Generator function for \ref FunctorGridFunction expressions
+  template <class Function, class GridView,
+    REQUIRES(Concepts::CallableDomain<Function>)>
+  auto makeGridFunction(Function const& f, GridView const& gridView, Dune::PriorityTag<3>)
+  {
+    return AnalyticGridFunction<Function, GridView>{f, gridView};
+  }
+
+  template <class Function, class GridView>
+  auto makeGridFunction(AnalyticPreGridFunction<Function> const& pre, GridView const& gridView, Dune::PriorityTag<3>)
+  {
+    return AnalyticGridFunction<Function, GridView>{pre.fct_, gridView};
+  }
+
+  /** @} **/
+
+} // end namespace AMDiS

dune/amdis/gridfunctions/ConstantGridFunction.hpp

+#pragma once
+
+#include <functional>
+#include <type_traits>
+
+#include <dune/common/diagonalmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/identitymatrix.hh>
+#include <dune/functions/common/defaultderivativetraits.hh>
+
+#include <dune/amdis/common/Utility.hpp>
+#include <dune/amdis/gridfunctions/AnalyticGridFunction.hpp>
+#include <dune/amdis/gridfunctions/GridFunctionConcepts.hpp>
+
+namespace AMDiS
+{
+  /**
+    * \addtogroup GridFunctions
+    * @{
+    **/
+
+  template <class T>
+  struct ConstantFunction
+  {
+    using Range = Underlying_t<T>;
+
+    using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits<Range>::Range;
+    using Derivative = ConstantFunction<DerivativeRange>;
+
+    explicit ConstantFunction(T const& value)
+      : value_(value)
+    {}
+
+    template <class GlobalCoordinate>
+    Range const& operator()(GlobalCoordinate const& /*x*/) const
+    {
+      return value_;
+    }
+
+    friend Derivative derivative(ConstantFunction const& /*f*/)
+    {
+      return Derivative{DerivativeRange(0)};
+    }
+
+    friend int order(ConstantFunction const& /*f*/)
+    {
+      return 0;
+    }
+
+  private:
+    T value_;
+  };
+
+
+  /** @} **/
+
+  namespace Concepts
+  {
+    /** \addtogroup Concepts
+     *  @{
+     **/
+
+    namespace Definition
+    {
+      template <class T>
+      struct ConstantToGridFunction
+        : std::is_arithmetic<T> {};
+
+      template <class T>
+      struct ConstantToGridFunction<std::reference_wrapper<T>>
+        : ConstantToGridFunction<T> {};
+
+      template <class T, int N>
+      struct ConstantToGridFunction<Dune::FieldVector<T, N>>
+        : ConstantToGridFunction<T> {};
+
+      template <class T, int N, int M>
+      struct ConstantToGridFunction<Dune::FieldMatrix<T, N, M>>
+        : ConstantToGridFunction<T> {};
+
+      template <class T, int N>
+      struct ConstantToGridFunction<Dune::IdentityMatrix<T, N>>
+        : ConstantToGridFunction<T> {};
+
+      template <class T, int N>
+      struct ConstantToGridFunction<Dune::DiagonalMatrix<T, N>>
+        : ConstantToGridFunction<T> {};
+
+    } // end namespace Definition
+
+
+    template <class T>
+    constexpr bool ConstantToGridFunction =
+      Definition::ConstantToGridFunction<T>::value;
+
+    /** @} **/
+
+  } // end namespace Concepts
+
+
+  /**
+    * \addtogroup GridFunctions
+    * @{
+    **/
+
+  template <class T, class GridView,
+    REQUIRES(Concepts::ConstantToGridFunction<T>)>
+  auto makeGridFunction(T const& value, GridView const& gridView, Dune::PriorityTag<2>)
+  {
+    return makeGridFunction(ConstantFunction<T>{value}, gridView, Dune::PriorityTag<3>{});
+  }
+
+  /** @} **/
+
+} // end namespace AMDiS

dune/amdis/gridfunctions/CoordsGridFunction.hpp

+#pragma once
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/identitymatrix.hh>
+#include <dune/common/typeutilities.hh>
+
+namespace AMDiS
+{
+  /**
+    * \addtogroup GridFunctions
+    * @{
+    **/
+
+  struct CoordsFunction
+  {
+    struct Derivative
+    {
+      template <class T, int N>
+      Dune::IdentityMatrix<T, N> const& operator()(Dune::FieldVector<T, N> const& x) const
+      {
+        return Dune::IdentityMatrix<T,N>{};
+      }
+    };
+
+    template <class T, int N>
+    Dune::FieldVector<T, N> const& operator()(Dune::FieldVector<T, N> const& x) const
+    {
+      return x;
+    }
+
+    friend Derivative derivative(CoordsFunction const& /*f*/)
+    {
+      return Derivative{};
+    }
+
+    friend int order(CoordsFunction const& /*f*/)
+    {
+      return 1;
+    }
+  };
+
+
+  struct CoordsCompFunction
+  {
+    struct Derivative
+    {
+      explicit Derivative(int comp)
+        : comp_(comp)
+      {}
+
+      template <class T, int N>
+      Dune::FieldVector<T, N> const& operator()(Dune::FieldVector<T, N> const& x) const
+      {
+        Dune::FieldVector<T, N> result(0);
+        result[comp_] = T(1);
+
+        return result;
+      }
+
+    private:
+      int comp_;
+    };
+
+    explicit CoordsCompFunction(int comp)
+      : comp_(comp)
+    {}
+
+    template <class T, int N>
+    T const& operator()(Dune::FieldVector<T, N> const& x) const
+    {
+      return x[comp_];
+    }
+
+    friend Derivative derivative(CoordsCompFunction const& f)
+    {
+      return Derivative{f.comp_};
+    }
+
+    friend int order(CoordsCompFunction const& /*f*/)
+    {
+      return 1;
+    }
+
+  private:
+    int comp_;
+  };
+
+  namespace Concepts { namespace Definition
+  {
+    template <>
+    struct IsPreGridFunction<CoordsFunction>
+      : std::true_type {};
+
+    template <>
+    struct IsPreGridFunction<CoordsCompFunction>
+      : std::true_type {};
+  }}
+
+  /// Generator for \ref CoordsFunction
+  inline auto X()
+  {
+    return CoordsFunction{};
+  }
+
+  /// Generator for \ref CoordsCompFunction
+  inline auto X(int comp)
+  {
+    return CoordsCompFunction{comp};
+  }
+
+
+  template <class GridView>
+  auto makeGridFunction(CoordsFunction const& f, GridView const& gridView, Dune::PriorityTag<1>)
+  {
+    return makeAnalyticGridFunction(f, gridView);
+  }
+
+  template <class GridView>
+  auto makeGridFunction(CoordsCompFunction const& f, GridView const& gridView, Dune::PriorityTag<1>)
+  {
+    return makeAnalyticGridFunction(f, gridView);
+  }
+
+  /** @} **/
+
+} // end namespace AMDiS

dune/amdis/DOFVectorView.hpp → dune/amdis/gridfunctions/DOFVectorView.hpp

@@ -9,9 +9,6 @@
 #include <dune/functions/gridfunctions/gridviewentityset.hh>
 #include <dune/typetree/childextraction.hh>
 
-#include <dune/amdis/terms/GridViewExpression.hpp>
-#include <dune/amdis/terms/TermGenerator.hpp>
-
 namespace AMDiS
 {
   template <class GlobalBasisType, class TreePathType, bool isConst = true>
@@ -45,6 +42,7 @@ namespace AMDiS
     template <class Block>
     using Flat = Dune::Functions::FlatVectorBackend<Block>;
 
+
   public: // a local view on the gradients
 
     class GradientLocalFunction
@@ -209,10 +207,10 @@ namespace AMDiS
       return LocalFunction{self};
     }
 
-    /// Create a local function implementing the gradient of the DOFVector
-    friend GradientLocalFunction derivative(DOFVectorView const& self)
+    friend int order(DOFVectorView const& self)
     {
-      return GradientLocalFunction{self};
+      auto node = child(self.basis().localView().tree(), self.treePath());
+      return getPolynomialDegree(self.basis().gridView(), node);
     }
 
     EntitySet const& entitySet() const
@@ -273,31 +271,14 @@ namespace AMDiS
 
   public:
 
-    /// Interpolation of expr to DOFVector
-    template <class PreExpr>
-    DOFVectorView& operator<<(PreExpr&& preExpr)
+    /// Interpolation of GridFunction to DOFVector
+    template <class Expr>
+    DOFVectorView& interpolate(Expr&& expr)
     {
-      // create expression from passed parameter
-      typename ToTerm<PreExpr>::type expr = toTerm(std::forward<PreExpr>(preExpr));
-
       auto const& basis = DOFVectorConstView::basis();
       auto const& treePath = DOFVectorConstView::treePath();
 
-      auto gridViewExpr = makeGridViewExpression(expr, basis.gridView());
-      DOFVector<GlobalBasis> tmp(basis, "tmp");
-      Dune::Functions::interpolate(basis, treePath, tmp, gridViewExpr);
-
-      // move data from temporary vector into stored DOFVector
-      mutableDofVector_->getVector() = std::move(tmp.getVector());
-      return *this;
-    }
-
-    /// Interpolation of expr to DOFVector
-    template <class GridFct>
-    DOFVectorView& interpol(GridFct const& gridFct)
-    {
-      auto const& basis = DOFVectorConstView::basis();
-      auto const& treePath = DOFVectorConstView::treePath();
+      auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), basis.gridView(), Dune::PriorityTag<42>{});
 
       DOFVector<GlobalBasis> tmp(basis, "tmp");
       Dune::Functions::interpolate(basis, treePath, tmp, gridFct);

dune/amdis/DOFVectorView.inc.hpp → dune/amdis/gridfunctions/DOFVectorView.inc.hpp

@@ -55,6 +55,7 @@ LocalFunction::operator()(LocalDomain const& x) const
   return y;
 }
 
+
 template <class GlobalBasis, class TreePath>
 typename DOFVectorView<GlobalBasis, TreePath, true>::DerivativeRange DOFVectorView<GlobalBasis, TreePath, true>::
 GradientLocalFunction::operator()(LocalDomain const& x) const

dune/amdis/terms/DerivativeGridFunction.hpp → dune/amdis/gridfunctions/DerivativeGridFunction.hpp

@@ -4,10 +4,12 @@
 
 #include <dune/functions/common/defaultderivativetraits.hh>
 
+#include <dune/amdis/gridfunctions/GridFunctionConcepts.hpp>
+
 namespace AMDiS
 {
   /**
-    * \addtogroup OperatorTerms
+    * \addtogroup GridFunctions
     * @{
     **/
 
@@ -21,7 +23,7 @@ namespace AMDiS
   public:
     using Range = typename Dune::Functions::DefaultDerivativeTraits<GridFctRange(GridFctDomain)>::Range;
     using Domain = GridFctDomain;
-    using LocalFunction = std::decay_t<decltype(derivative(std::declval<GridFct>()))>;
+    using LocalFunction = std::decay_t<decltype(derivative(localFunction(std::declval<GridFct>())))>;
 
     using EntitySet = typename GridFct::EntitySet;
 
@@ -38,8 +40,12 @@ namespace AMDiS
 
     friend LocalFunction localFunction(DerivativeGridFunction const& gf)
     {
-      msg("DerivativeGridFunction::localFunction()\n");
-      return derivative(gf.gridFct_);
+      return derivative(localFunction(gf.gridFct_));
+    }
+
+    friend int order(DerivativeGridFunction const& gf)
+    {
+      return std::max(0, order(gf.gridFct_)-1);
     }
 
     EntitySet const& entitySet() const
@@ -53,10 +59,12 @@ namespace AMDiS
 
 
   /// Generator function for \ref DerivativeGridFunction expressions
-  template <class GridFct>
-  auto makeDerivativeGridFunction(GridFct&& gridFct)
+  template <class GridFct,
+    REQUIRES(not Concepts::HasDerivative<GridFct> &&
+                 Concepts::HasLocalFunctionDerivative<GridFct>)>
+  auto derivative(GridFct const& gridFct)
   {
-    return DerivativeGridFunction<std::decay_t<GridFct>>{std::forward<GridFct>(gridFct)};
+    return DerivativeGridFunction<GridFct>{gridFct};
   }
 
   /** @} **/

dune/amdis/terms/FunctorGridFunction.hpp → dune/amdis/gridfunctions/FunctorGridFunction.hpp

@@ -5,11 +5,15 @@
 
 #include <dune/amdis/common/IndexSeq.hpp>
 #include <dune/amdis/common/Loops.hpp>
+#include <dune/amdis/common/Mpl.hpp>
+#include <dune/amdis/gridfunctions/GridFunctionConcepts.hpp>
+