interpolation of expressions added

dune/amdis/ContextGeometry.hpp

@@ -14,6 +14,14 @@ namespace AMDiS
     Geometry const& geometry;
     LocalGeometry const& localGeometry;
 
+    ContextGeometry(LocalContext const& localContext,
+                    Geometry const& geometry,
+                    LocalGeometry const& localGeometry)
+      : localContext(localContext)
+      , geometry(geometry)
+      , localGeometry(localGeometry)
+    {}
+
     /// Coordinate `p` given in `localGeometry`, transformed to coordinate in `geometry`.
     template <class Coordinate>
     decltype(auto) position(Coordinate const& p) const

dune/amdis/FileWriter.hpp

@@ -80,7 +80,7 @@ namespace AMDiS
     {
       using Dune::Functions::BasisBuilder::makeBasis;
       using Dune::Functions::BasisBuilder::lagrange;
-      auto fct = makeDiscreteFunction(basis_,treePath_,vector_);
+      auto fct = makeDiscreteFunctionPtr(basis_,treePath_,vector_);
 
       auto p1basis = makeBasis(basis_->gridView(), lagrange<1>());
       auto data = makeDOFVector(p1basis, name_);
@@ -102,7 +102,7 @@ namespace AMDiS
       using Dune::Functions::BasisBuilder::flatLexicographic;
 
       assert( C == dow );
-      auto fct = makeDiscreteFunction(basis_,treePath_,vector_);
+      auto fct = makeDiscreteFunctionPtr(basis_,treePath_,vector_);
 
       auto p1basis = makeBasis(basis_->gridView(), power<C>(lagrange<1>(), flatLexicographic()));
       auto data = makeDOFVector(p1basis, name_);

dune/amdis/ProblemStat.hpp

@@ -228,12 +228,12 @@ namespace AMDiS
     std::shared_ptr<SystemVector>       getSolution()       { return solution; }
     std::shared_ptr<SystemVector const> getSolution() const { return solution; }
 
-    /// Return a view to a solution component
+    /// Return a view to a solution component (as shared_ptr)
     template <class TreePath>
     auto getSolution(TreePath path) const
     {
       auto tp = makeTreePath(path);
-      return makeDiscreteFunction(globalBasis, tp, solution);
+      return makeDiscreteFunctionPtr(globalBasis, tp, solution);
     }
 
 

dune/amdis/common/IndexSeq.hpp

@@ -5,21 +5,6 @@
 // from http://stackoverflow.com/questions/17424477/implementation-c14-make-integer-sequence/17426611#17426611
 namespace AMDiS
 {
-
-
-  // template <std::size_t i, std::size_t j>
-  // constexpr index_<i+j> operator+(index_<i>, index_<j>) { return {}; }
-
-  // template <std::size_t i, std::size_t j>
-  // constexpr index_<i-j> operator-(index_<i>, index_<j>) { return {}; }
-
-  // template <std::size_t i, std::size_t j>
-  // constexpr index_<i*j> operator*(index_<i>, index_<j>) { return {}; }
-
-  // template <std::size_t i, std::size_t j>
-  // constexpr index_<i/j> operator/(index_<i>, index_<j>) { return {}; }
-
-
   namespace Impl
   {
     template <std::size_t s, class S>

dune/amdis/common/ScalarTypes.hpp

@@ -8,10 +8,12 @@ namespace AMDiS
   namespace traits
   {
     template <class T>
-    using is_integral = std::is_integral<std::decay_t<T>>;
+    struct is_integral
+        : public std::is_integral<std::decay_t<T>> {};
 
     template <class T>
-    using is_arithmetic = std::is_arithmetic<std::decay_t<T>>;
+    struct is_arithmetic
+        : public std::is_arithmetic<std::decay_t<T>> {};
 
   } // end namespace traits
 

dune/amdis/common/TypeDefs.hpp

@@ -2,6 +2,7 @@
 
 #include <type_traits>
 
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
 #include <boost/numeric/mtl/vector/dense_vector.hpp>
 #include <boost/numeric/mtl/matrix/dense2D.hpp>
 
@@ -17,26 +18,6 @@ namespace AMDiS
 
   } // end namespace Impl
 
-  enum FirstOrderType
-  {
-    GRD_PSI,
-    GRD_PHI
-  };
-
-  template <FirstOrderType type>
-  using FirstOrderType_t = std::integral_constant<FirstOrderType, type>;
-
-  template <FirstOrderType type>
-  constexpr FirstOrderType_t<type> FirstOrderType_ = {};
-
   struct BoundaryType { int b; };
 
-  namespace tag
-  {
-    struct zot {};
-    struct fot_grd_phi {};
-    struct fot_grd_psi {};
-    struct sot {};
-  }
-
 } // end namespace AMDiS

dune/amdis/common/Utility.hpp

@@ -21,16 +21,22 @@ namespace AMDiS
 
   namespace Impl
   {
-      // workaround for MSVC (problems with alias templates in pack expansion)
-      template <class, class T>
-      struct InvokeType { using type = T; };
-
-      template <class, class, class T>
-      struct InvokeType2 { using type = T; };
+    // workaround for MSVC (problems with alias templates in pack expansion)
+    template <class, class T>
+    struct InvokeType { using type = T; };
+
+    template <class, class, class T>
+    struct InvokeType2 { using type = T; };
+
+    template <class T>
+    struct ValueType
+    {
+      using type = typename InvokeType<T, typename T::value_type>::type;
+    };
   }
 
   template <class T>
-  using Value_t = typename Impl::InvokeType<T, typename T::value_type>::type;
+  using Value_t = typename Impl::ValueType<T>::type;
 
   template <class T>
   using Size_t = typename Impl::InvokeType<T, typename T::size_type>::type;
@@ -44,6 +50,24 @@ namespace AMDiS
   template <class T1, class T2>
   using Common_t = typename Impl::InvokeType2<T1, T2, typename std::common_type<T1,T2>::type>::type;
 
+  namespace Impl
+  {
+    template <class T>
+    struct UnderlyingType
+    {
+      using type = T;
+    };
+
+    template <class T>
+    struct UnderlyingType<std::reference_wrapper<T>>
+    {
+      using type = T;
+    };
+  }
+
+  template <class T>
+  using Underlying_t = typename Impl::UnderlyingType<T>::type;
+
 
   // ---------------------------------------------------------------------------
 

dune/amdis/common/ValueCategory.hpp

@@ -35,12 +35,6 @@ namespace AMDiS
     using type = tag::scalar;
   };
 
-  template <class T>
-  struct ValueCategory<std::reference_wrapper<T>, std::enable_if_t< std::is_arithmetic<T>::value >>
-  {
-    using type = tag::scalar;
-  };
-
   template <class K, int SIZE>
   struct ValueCategory< Dune::FieldVector<K, SIZE> >
   {
@@ -65,6 +59,12 @@ namespace AMDiS
     using type = tag::scalar;
   };
 
+  template <class T>
+  struct ValueCategory<std::reference_wrapper<T>>
+  {
+    using type = typename ValueCategory<T>::type;
+  };
+
   namespace Category
   {
     template <class T>

dune/amdis/linear_algebra/mtl/DOFVector.hpp

@@ -123,9 +123,9 @@ namespace AMDiS
     /// \brief interpolate a function \p f to the basis \ref feSpace and store the
     /// coefficients in \ref vector.
     template <class F>
-    void interpol(F const& f)
+    void interpol(F&& f)
     {
-      Dune::Functions::interpolate(feSpace, *this, f);
+      Dune::Functions::interpolate(feSpace, *this, std::forward<F>(f));
     }
 
     /// Scale each DOFVector by the factor \p s.

dune/amdis/terms/ConstantTerm.hpp

@@ -3,6 +3,7 @@
 #include <type_traits>
 
 #include <dune/amdis/common/ScalarTypes.hpp>
+#include <dune/amdis/common/Utility.hpp>
 #include <dune/amdis/terms/ExpressionBase.hpp>
 
 namespace AMDiS
@@ -19,44 +20,32 @@ namespace AMDiS
   {
   public:
     // TODO: possibly convert to plain type, in case of expression templates.
-    using value_type = T;
+    using value_type = Underlying_t<T>;
 
     explicit ConstantTerm(value_type const& value)
       : ExpressionBase(0)
       , value_(value)
     {}
 
-    value_type const& operator[](std::size_t iq) const
+    template <class Context, class LocalCoordinate>
+    value_type const& local(Context const& /*context*/, LocalCoordinate const& /*lambda*/) const
     {
       return value_;
     }
 
-  private:
-    value_type value_;
-  };
-
-
-  /// An expression representing a reference to a constant (arithmetic) value, \see constant
-  template <class T>
-  class ConstantTerm<std::reference_wrapper<T>>
-      : public ExpressionBase
-  {
-  public:
-    // TODO: possibly convert to plain type, in case of expression templates.
-    using value_type = T;
-
-    explicit ConstantTerm(std::reference_wrapper<value_type> const& ref)
-      : ExpressionBase(0)
-      , ref_(ref)
-    {}
+    template <class GlobalCoordinate>
+    value_type const& global(GlobalCoordinate const& /*x*/) const
+    {
+      return value_;
+    }
 
-    value_type const& operator[](std::size_t iq) const
+    value_type const& operator[](std::size_t /*iq*/) const
     {
-      return ref_;
+      return value_;
     }
 
   private:
-    std::reference_wrapper<value_type> ref_;
+    value_type value_;
   };
 
 

dune/amdis/terms/CoordsTerm.hpp

@@ -64,7 +64,19 @@ namespace AMDiS
     {
       values_.resize(points.size());
       for (std::size_t iq = 0; iq < points.size(); ++iq)
-        values_[iq] = fct_(context.global(points[iq].position()));
+        values_[iq] = local(context, points[iq].position());
+    }
+
+    template <class Context, class LocalCoordinate>
+    value_type local(Context const& context, LocalCoordinate const& lambda) const
+    {
+      return fct_(context.global(lambda));
+    }
+
+    template <class GlobalCoordinate>
+    value_type global(GlobalCoordinate const& x) const
+    {
+      return fct_(x);
     }
 
     value_type const& operator[](std::size_t iq) const
@@ -108,6 +120,18 @@ namespace AMDiS
         values_[iq] = context.global( points[iq].position() );
     }
 
+    template <class Context, class LocalCoordinate>
+    value_type const& local(Context const& context, LocalCoordinate const& lambda) const
+    {
+      return context.global(lambda);
+    }
+
+    template <class GlobalCoordinate>
+    value_type const& global(GlobalCoordinate const& x) const
+    {
+      return x;
+    }
+
     value_type const& operator[](std::size_t iq) const
     {
       return values_[iq];

dune/amdis/terms/DOFVectorTerm.hpp

@@ -54,11 +54,24 @@ namespace AMDiS
     }
 
     template <class Context, class PointList>
-    void init(Context const& /*context*/, PointList const& points)
+    void init(Context const& context, PointList const& points)
     {
       values_.resize(points.size());
       for (std::size_t iq = 0; iq < points.size(); ++iq)
-        values_[iq] = factor_ * localFunction_(points[iq].position());
+        values_[iq] = factor_ * local(context, points[iq].position());
+    }
+
+    template <class Context, class LocalCoordinate>
+    value_type local(Context const& /*context*/, LocalCoordinate const& lambda) const
+    {
+      return localFunction_(lambda);
+    }
+
+    template <class GlobalCoordinate>
+    value_type global(GlobalCoordinate const& x) const
+    {
+      error_exit("Not implemented.");
+      return value_type(0);
     }
 
     value_type const& operator[](std::size_t iq) const
@@ -86,10 +99,12 @@ namespace AMDiS
    *  an optional factor \p factor the value of the DOFVector at the quadrature
    *  points are scaled with.
    **/
-  template <class DiscreteGlobalBasisFct>
-  auto valueOf(std::shared_ptr<DiscreteGlobalBasisFct> const& vector, double factor = 1.0)
+  template <class... Args>
+  auto valueOf(std::shared_ptr<Dune::Functions::DiscreteGlobalBasisFunction<Args...>> const& dofvector,
+               double factor = 1.0)
   {
-    return DOFVectorTerm<DiscreteGlobalBasisFct>{vector, factor};
+    using DiscreteGlobalBasisFct = Dune::Functions::DiscreteGlobalBasisFunction<Args...>;
+    return DOFVectorTerm<DiscreteGlobalBasisFct>{dofvector, factor};
   }
 
 

dune/amdis/terms/FunctorTerm.hpp

@@ -60,16 +60,31 @@ namespace AMDiS
       });
     }
 
+    template <class Context, class LocalCoordinate>
+    value_type local(Context const& context, LocalCoordinate const& lambda) const
+    {
+      return eval([&context,&lambda](auto const& term) { return term.local(context, lambda); },
+        MakeSeq_t<sizeof...(Terms)>{});
+    }
+
+    template <class GlobalCoordinate>
+    value_type global(GlobalCoordinate const& x) const
+    {
+      return eval([&x](auto const& term) { return term.global(x); },
+        MakeSeq_t<sizeof...(Terms)>{});
+    }
+
     auto operator[](std::size_t iq) const
     {
-      return eval(iq, MakeSeq_t<sizeof...(Terms)>{});
+      return eval([iq](auto const& term) { return term[iq]; },
+        MakeSeq_t<sizeof...(Terms)>{});
     }
 
   private:
-    template <std::size_t... Is>
-    auto eval(std::size_t iq, Indices<Is...>) const
+    template <class Op, std::size_t... Is>
+    auto eval(Op op, Indices<Is...>) const
     {
-      return fct_(std::get<Is>(terms_)[iq]...);
+      return fct_(op(std::get<Is>(terms_))...);
     }
 
   private:

dune/amdis/terms/GridViewExpression.hpp

+#pragma once
+
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+
+#include <dune/amdis/Output.hpp>
+#include <dune/amdis/terms/LocalExpression.hpp>
+
+namespace AMDiS
+{
+  template <class Expr, class GridView>
+  class GridViewExpression
+  {
+    using Element = typename GridView::template Codim<0>::Entity;
+    using Domain = typename Element::Geometry::GlobalCoordinate;
+    using Range = typename Expr::value_type;
+
+    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
+
+    using LocalExpr = LocalExpression<Expr, GridView>;
+
+  public:
+    GridViewExpression(Expr const& expr, GridView const& gridView)
+      : expr_(expr)
+      , entitySet_(gridView)
+    {}
+
+    Range operator()(Domain const& x) const
+    {
+      return expr_.global(x);
+    }
+
+    friend LocalExpr localFunction(GridViewExpression const& self)
+    {
+      return LocalExpr{self.expr_};
+    }
+
+    EntitySet const& entitySet() const
+    {
+      return entitySet_;
+    }
+
+
+  private:
+
+    Expr expr_;
+    EntitySet entitySet_;
+  };
+
+  template <class Expr, class GridView,
+    std::enable_if_t<std::is_base_of<ExpressionBase, Expr>::value, int> = 0>
+  auto makeGridViewExpression(Expr const& expr, GridView const& gridView)
+  {
+    return GridViewExpression<Expr, GridView>{expr, gridView};
+  }
+
+} // end namespace AMDiS

dune/amdis/terms/Interpolation.hpp

+#pragma once
+
+#include <memory>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/amdis/terms/GridViewExpression.hpp>
+#include <dune/amdis/terms/TermGenerator.hpp>
+
+namespace AMDiS
+{
+  template <class... Args, class PreExpr>
+  auto operator<<(std::shared_ptr<Dune::Functions::DiscreteGlobalBasisFunction<Args...>> const& dofvector,
+                  PreExpr&& preExpr)
+  {
+    using DiscreteGlobalBasisFct = Dune::Functions::DiscreteGlobalBasisFunction<Args...>;
+    typename ToTerm<PreExpr>::type expr = toTerm(std::forward<PreExpr>(preExpr));
+
+    using Vector = typename DiscreteGlobalBasisFct::Vector;
+    Vector vec(dofvector->basis(), "tmp");
+
+    auto gridViewExpr = makeGridViewExpression(expr, dofvector->basis().gridView());
+    Dune::Functions::interpolate(dofvector->basis(), dofvector->treePath(), vec, gridViewExpr);
+
+    const_cast<Vector&>(dofvector->dofs()).getVector() = vec.getVector();
+    return dofvector;
+  }
+
+  template <class... Args, class PreExpr>
+  auto& operator<<(Dune::Functions::DiscreteGlobalBasisFunction<Args...>& dofvector,
+                   PreExpr& preExpr)
+  {
+    using DiscreteGlobalBasisFct = Dune::Functions::DiscreteGlobalBasisFunction<Args...>;
+    typename ToTerm<PreExpr>::type expr = toTerm(std::forward<PreExpr>(preExpr));
+
+    using Vector = typename DiscreteGlobalBasisFct::Vector;
+    Vector vec(dofvector.basis(), "tmp");
+
+    auto gridViewExpr = makeGridViewExpression(expr, dofvector.basis().gridView());
+    Dune::Functions::interpolate(dofvector.basis(), dofvector.treePath(), vec, gridViewExpr);
+
+    const_cast<Vector&>(dofvector.dofs()).getVector() = vec.getVector();
+    return dofvector;
+  }
+
+} // end namespace AMDiS

dune/amdis/terms/LocalExpression.hpp

+#pragma once
+
+#include <cassert>
+
+#include <dune/common/std/optional.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+
+#include <dune/amdis/ContextGeometry.hpp>
+
+namespace AMDiS
+{
+  template <class Expr, class GridView>
+  class LocalExpression
+  {
+    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
+
+    using Element = typename EntitySet::Element;
+    using Geometry = typename Element::Geometry;
+
+    using Domain = typename EntitySet::LocalCoordinate;
+    using Range = typename Expr::value_type;
+
+    using LocalContext = Element;
+    using LocalGeometry = Geometry;
+    using GeometryWrapper = ContextGeometry<LocalContext, Geometry, LocalGeometry>; // TODO: allow intersections!
+
+  public:
+
+    /// Constructor. Stores a copy to the expression.
+    explicit LocalExpression(Expr const& expr)
+      : expr_(expr)
+    {}
+
+
+    /// Evaluate expression `expr_` at local coordinate `lambda`
+    /// [expects: LocalExpression is bound to element]
+    Range operator()(Domain const& lambda) const
+    {
+      assert( bound_ );
+      return expr_.local(context_.value(), lambda);
+    }
+
+    /// Bind expression to a local context. Must be called before you can evaluate it.
+    /// [expects: local function is not bound to other element]
+    void bind(Element const& element)
+    {
+      assert( not bound_ ); // NOTE: interpolateTreeSubset(...) does not unbind the element!
+      geometry_.emplace(element.geometry());
+      context_.emplace(element, *geometry_, *geometry_);
+
+      expr_.bind(element);
+      bound_ = true;
+    }
+
+    /// Unbind from local context
+    void unbind()
+    {
+      context_.reset();
+      geometry_.reset();
+      bound_ = false;
+    }
+
+    /// Obtain local contex this LocalFunction is bound to
+    LocalContext const& localContext() const
+    {
+      assert( bound_ );
+      return context_.value().localContext;
+    }
+
+  private:
+    Expr expr_;
+
+    Dune::Std::optional<GeometryWrapper> context_;
+    Dune::Std::optional<Geometry> geometry_;
+    bool bound_ = false;
+  };
+
+} // end namespace AMDiS

dune/amdis/utility/DiscreteFunction.hpp

@@ -26,9 +26,9 @@ namespace AMDiS
 
 
   template <class GlobalBasis, class TreePath, class Vector>
-  auto makeDiscreteFunction(std::shared_ptr<GlobalBasis> const& basis,
-                            TreePath const& tp,
-                            std::shared_ptr<Vector> const& vec)
+  auto makeDiscreteFunctionPtr(std::shared_ptr<GlobalBasis> const& basis,
+                               TreePath const& tp,
+                               std::shared_ptr<Vector> const& vec)
   {
     using NTRM = typename Impl::DiscreteFunction<GlobalBasis,TreePath,Vector>::NTRM;
     using DiscreteFunction = typename Impl::DiscreteFunction<GlobalBasis,TreePath,Vector>::type;

src/CMakeLists.txt

 
-set(projects2d "ellipt" "heat" "vecellipt" "stokes1") # "stokes2" "test") #"pfc"  "navier_stokes" "gridview")
+set(projects2d "ellipt" "heat" "vecellipt" "expressions" "stokes1") # "stokes2" "test") #"pfc"  "navier_stokes" "gridview")
 
 foreach(project ${projects2d})
     add_