Merge branch 'issue/discretefunction_dofvectorview' into 'master'

DiscreteFunction vs. DOFVectorView

See merge request !125

src/amdis/DOFVector.hpp

@@ -10,7 +10,7 @@
 #include <amdis/Observer.hpp>
 #include <amdis/common/Concepts.hpp>
 #include <amdis/common/TypeTraits.hpp>
-#include <amdis/gridfunctions/GridFunction.hpp>
+#include <amdis/gridfunctions/DiscreteFunction.hpp>
 #include <amdis/typetree/TreePath.hpp>
 
 namespace Dune
@@ -88,7 +88,7 @@ namespace AMDiS
     auto child(TreePath const& path = {})
     {
       auto&& tp = makeTreePath(path);
-      return makeDOFVectorView(*this, tp);
+      return makeDiscreteFunction(*this, tp);
     }
 
     template <class TreePath = RootTreePath>
@@ -101,7 +101,7 @@ namespace AMDiS
 
     /// Interpolation of GridFunction to DOFVector, assuming that there is no
     /// reference to this DOFVector in the expression.
-    /// See \ref DOFVectorView::interpolate_noalias
+    /// See \ref DiscreteFunction::interpolate_noalias
     template <class Expr, class Tag = tag::average>
     void interpolate_noalias(Expr&& expr, Tag strategy)
     {
@@ -109,7 +109,7 @@ namespace AMDiS
     }
 
     /// Interpolation of GridFunction to DOFVector.
-    /// See \ref DOFVectorView::interpolate
+    /// See \ref DiscreteFunction::interpolate
     template <class Expr, class Tag = tag::average>
     void interpolate(Expr&& expr, Tag strategy)
     {
@@ -117,7 +117,7 @@ namespace AMDiS
     }
 
     /// Interpolation of GridFunction to DOFVector.
-    /// See \ref DOFVectorView::interpolate
+    /// See \ref DiscreteFunction::interpolate
     template <class Expr>
     DOFVector& operator<<(Expr&& expr)
     {

src/amdis/ProblemStat.hpp

@@ -40,7 +40,6 @@
 
 #include <amdis/GridFunctions.hpp>
 #include <amdis/gridfunctions/DiscreteFunction.hpp>
-#include <amdis/gridfunctions/DOFVectorView.hpp>
 
 #include <amdis/io/FileWriterBase.hpp>
 

src/amdis/gridfunctions/CMakeLists.txt

@@ -8,7 +8,7 @@ install(FILES
     DerivativeGridFunction.hpp
     DiscreteFunction.hpp
     DiscreteFunction.inc.hpp
-    DOFVectorView.hpp
+    DiscreteLocalFunction.inc.hpp
     FunctorGridFunction.hpp
     GridFunction.hpp
     OperationsGridFunction.hpp

src/amdis/gridfunctions/DOFVectorView.hpp

-#pragma once
-
-#include <amdis/functions/Interpolate.hpp>
-#include <amdis/gridfunctions/DiscreteFunction.hpp>
-#include <amdis/gridfunctions/GridFunction.hpp>
-
-namespace AMDiS
-{
-  /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
-  template <class GB, class VT, class TP>
-  class DOFVectorView
-      : public DiscreteFunction<GB, VT, TP>
-  {
-    using Self = DOFVectorView;
-    using Super = DiscreteFunction<GB, VT, TP>;
-
-    using GlobalBasis = GB;
-    using TreePath = TP;
-
-  public:
-    /// Constructor. Stores a pointer to the mutable `dofvector`.
-    template <class PreTreePath>
-    DOFVectorView(DOFVector<GB,VT>& dofVector, PreTreePath const& preTreePath)
-      : Super(dofVector, makeTreePath(preTreePath))
-      , mutableDofVector_(&dofVector)
-    {}
-
-    /// Constructor forwards to the treePath constructor, with empty TreePath
-    DOFVectorView(DOFVector<GB,VT>& dofVector)
-      : DOFVectorView(dofVector, Dune::TypeTree::hybridTreePath())
-    {}
-
-  public:
-    /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
-    /// reference to this DOFVector in the expression.
-    /**
-     * **Example:**
-     * ```
-     * auto v = makeDOFVectorView(prob.solutionVector(),0);
-     * v.interpolate_noalias([](auto const& x) { return x[0]; });
-     * ```
-     **/
-    template <class Expr, class Tag = tag::average>
-    void interpolate_noalias(Expr&& expr, Tag strategy = {})
-    {
-      auto const& basis = *this->basis();
-      auto const& treePath = this->treePath();
-
-      auto&& gf = makeGridFunction(FWD(expr), basis.gridView());
-
-      if (std::is_same<Tag, tag::average>::value) {
-        auto counter = coefficients();
-        AMDiS::interpolate(basis, coefficients(), gf, treePath, counter);
-
-        coefficients().forEach([&counter](std::size_t dof, auto& coeff)
-        {
-          coeff /= std::max(double(counter.at(dof)), 1.0);
-        });
-      } else {
-        AMDiS::interpolate(basis, coefficients(), gf, treePath);
-      }
-    }
-
-    /// \brief Interpolation of GridFunction to DOFVector
-    /**
-     * **Example:**
-     * ```
-     * auto v = makeDOFVectorView(prob.solutionVector(),0);
-     * v.interpolate(v + [](auto const& x) { return x[0]; });
-     * ```
-     * Allows to have a reference to the DOFVector in the expression, e.g. as
-     * \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction.
-     **/
-    template <class Expr, class Tag = tag::average>
-    void interpolate(Expr&& expr, Tag strategy = {})
-    {
-      // create temporary copy of data
-      DOFVector<GB,VT> tmp(coefficients());
-
-      Self tmpView{tmp, this->treePath()};
-      tmpView.interpolate_noalias(FWD(expr), strategy);
-
-      // move data from temporary vector into stored DOFVector
-      coefficients().backend() = std::move(tmp.backend());
-    }
-
-    /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
-    template <class Expr>
-    DOFVectorView& operator<<(Expr&& expr)
-    {
-      interpolate(FWD(expr));
-      return *this;
-    }
-
-    /// \brief interpolate `(*this) + expr` to DOFVector
-    template <class Expr>
-    DOFVectorView& operator+=(Expr&& expr)
-    {
-      interpolate((*this) + expr);
-      return *this;
-    }
-
-    /// \brief interpolate `(*this) - expr` to DOFVector
-    template <class Expr>
-    DOFVectorView& operator-=(Expr&& expr)
-    {
-      interpolate((*this) - expr);
-      return *this;
-    }
-
-    /// Return the mutable DOFVector
-    DOFVector<GB,VT>& coefficients() { return *mutableDofVector_; }
-
-    /// Return the const DOFVector
-    using Super::coefficients;
-
-  protected:
-    DOFVector<GB,VT>* mutableDofVector_;
-  };
-
-
-#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
-  // Deduction guide for DOFVectorView class
-  template <class GlobalBasis, class ValueType, class PreTreePath>
-  DOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector, PreTreePath const& preTreePath)
-    -> DOFVectorView<GlobalBasis, ValueType, TYPEOF(makeTreePath(preTreePath))>;
-
-  // Deduction guide for DOFVectorView class
-  template <class GlobalBasis, class ValueType>
-  DOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector)
-    -> DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>;
-#endif
-
-  /// A Generator for a mutable \ref DOFVectorView
-  template <class GlobalBasis, class ValueType, class PreTreePath>
-  auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector, PreTreePath const& preTreePath)
-  {
-    auto treePath = makeTreePath(preTreePath);
-    return DOFVectorView<GlobalBasis, ValueType, decltype(treePath)>{dofVector, treePath};
-  }
-
-  /// A Generator for a mutable \ref DOFVectorView
-  template <class GlobalBasis, class ValueType>
-  auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector)
-  {
-    auto treePath = Dune::TypeTree::hybridTreePath();
-    return DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath};
-  }
-
-} // end namespace AMDiS

src/amdis/gridfunctions/DiscreteFunction.hpp

@@ -16,23 +16,63 @@
 
 namespace AMDiS
 {
+  template <class GB, class VT>
+  class DOFVector;
+  
   /// \class DiscreteFunction
   /// \brief A view on a subspace of a \ref DOFVector
   /**
-   * \ingroup GridFunctions
-   *
-   * \tparam GB  Type of the global basis
-   * \tparam VT  Coefficient type of the DOFVector
-   * \tparam TP  A realization of \ref Dune::TypeTree::HybridTreePath
-   *
-   * **Requirements:**
-   * - GB models \ref Dune::Functions::Concept::GlobalBasis
-   **/
-  template <class GB, class VT, class TP>
-  class DiscreteFunction
+  * \ingroup GridFunctions
+  *
+  * \tparam GB  Type of the global basis
+  * \tparam VT  Coefficient type of the DOFVector
+  * \tparam TP  A realization of \ref Dune::TypeTree::HybridTreePath
+  * \tparam is_const  Specifies whether a const or mutable view is implemented.
+  *
+  * **Requirements:**
+  * - GB models \ref Dune::Functions::Concept::GlobalBasis
+  **/
+  template <class GB, class VT, class TP, bool is_const>
+  class DiscreteFunction;
+
+#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
+  // Deduction guide for DiscreteFunction class
+  template <class GB, class VT>
+  DiscreteFunction(DOFVector<GB, VT> const& dofVector)
+    -> DiscreteFunction<GB, VT, Dune::TypeTree::HybridTreePath<>, true>;
+
+  template <class GB, class VT>
+  DiscreteFunction(DOFVector<GB, VT>& dofVector)
+    -> DiscreteFunction<GB, VT, Dune::TypeTree::HybridTreePath<>, false>;
+#endif
+
+
+  /// A Generator for a const \ref DiscreteFunction
+  template <class GlobalBasis, class ValueType,
+            class PreTreePath = Dune::TypeTree::HybridTreePath<>>
+  auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector,
+                            PreTreePath const& preTreePath = {})
   {
-    // static_assert(Dune::IsNumber<VT>::value, "");
+    auto treePath = makeTreePath(preTreePath);
+    return DiscreteFunction<GlobalBasis, ValueType, decltype(treePath), true>{dofVector, treePath};
+  }
 
+  /// A Generator for a mutable \ref DiscreteFunction
+  template <class GlobalBasis, class ValueType,
+            class PreTreePath = Dune::TypeTree::HybridTreePath<>>
+  auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType>& dofVector,
+                            PreTreePath const& preTreePath = {})
+  {
+    auto treePath = makeTreePath(preTreePath);
+    return DiscreteFunction<GlobalBasis, ValueType, decltype(treePath), false>{dofVector, treePath};
+  }
+
+
+
+  /// A Const DiscreteFunction
+  template <class GB, class VT, class TP>
+  class DiscreteFunction<GB,VT,TP,true>
+  {
   private:
     using GlobalBasis = GB;
     using TreePath = TP;
@@ -71,18 +111,13 @@ namespace AMDiS
 
   public:
     /// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
-    DiscreteFunction(DOFVector<GB,VT> const& dofVector, TP const& treePath)
+    DiscreteFunction(DOFVector<GB,VT> const& dofVector, TP const& treePath = {})
       : dofVector_(&dofVector)
       , treePath_(treePath)
       , entitySet_(dofVector.basis()->gridView())
       , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(*dofVector.basis(), treePath))
     {}
 
-    /// Constructor forwards to the treePath constructor, with empty TreePath
-    DiscreteFunction(DOFVector<GB,VT> const& dofVector)
-      : DiscreteFunction(dofVector, Dune::TypeTree::hybridTreePath())
-    {}
-
     /// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive
     Range operator()(Domain const& x) const;
 
@@ -98,7 +133,6 @@ namespace AMDiS
       return entitySet_;
     }
 
-  public:
     /// \brief Return global basis bound to the DOFVector
     std::shared_ptr<GlobalBasis const> basis() const
     {
@@ -125,29 +159,86 @@ namespace AMDiS
   };
 
 
-#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
-  // Deduction guide for DiscreteFunction class
-  template <class GlobalBasis, class ValueType>
-  DiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector)
-    -> DiscreteFunction<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>;
-#endif
 
-  /// A Generator for a \ref DiscreteFunction
-  template <class GlobalBasis, class ValueType, class PreTreePath>
-  auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector, PreTreePath const& preTreePath)
+  /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
+  template <class GB, class VT, class TP>
+  class DiscreteFunction<GB,VT,TP,false>
+      : public DiscreteFunction<GB, VT, TP, true>
   {
-    auto treePath = makeTreePath(preTreePath);
-    return DiscreteFunction<GlobalBasis, ValueType, decltype(treePath)>{dofVector, treePath};
-  }
+    using Self = DiscreteFunction<GB, VT, TP, false>;
+    using Super = DiscreteFunction<GB, VT, TP, true>;
 
-  /// A Generator for a \ref DiscreteFunction
-  template <class GlobalBasis, class ValueType>
-  auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector)
-  {
-    auto treePath = Dune::TypeTree::hybridTreePath();
-    return DiscreteFunction<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath};
-  }
+    using GlobalBasis = GB;
+    using TreePath = TP;
+
+  public:
+    /// Constructor. Stores a pointer to the mutable `dofvector`.
+    DiscreteFunction(DOFVector<GB,VT>& dofVector, TP const& treePath = {})
+      : Super(dofVector, treePath)
+      , mutableDofVector_(&dofVector)
+    {}
+
+  public:
+    /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
+    /// reference to this DOFVector in the expression.
+    /**
+     * **Example:**
+     * ```
+     * auto v = makeDiscreteFunction(prob.solutionVector(),0);
+     * v.interpolate_noalias([](auto const& x) { return x[0]; });
+     * ```
+     **/
+    template <class Expr, class Tag = tag::average>
+    void interpolate_noalias(Expr&& expr, Tag strategy = {});
+
+    /// \brief Interpolation of GridFunction to DOFVector
+    /**
+     * **Example:**
+     * ```
+     * auto v = makeDiscreteFunction(prob.solutionVector(),0);
+     * v.interpolate(v + [](auto const& x) { return x[0]; });
+     * ```
+     * Allows to have a reference to the DOFVector in the expression, e.g. as
+     * \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction.
+     **/
+    template <class Expr, class Tag = tag::average>
+    void interpolate(Expr&& expr, Tag strategy = {});
+
+    /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
+    template <class Expr>
+    Self& operator<<(Expr&& expr)
+    {
+      interpolate(FWD(expr));
+      return *this;
+    }
+
+    /// \brief interpolate `(*this) + expr` to DOFVector
+    template <class Expr>
+    Self& operator+=(Expr&& expr)
+    {
+      interpolate((*this) + expr);
+      return *this;
+    }
+
+    /// \brief interpolate `(*this) - expr` to DOFVector
+    template <class Expr>
+    Self& operator-=(Expr&& expr)
+    {
+      interpolate((*this) - expr);
+      return *this;
+    }
+
+    /// Return the mutable DOFVector
+    DOFVector<GB,VT>& coefficients() { return *mutableDofVector_; }
+
+    /// Return the const DOFVector
+    using Super::coefficients;
+
+  protected:
+    DOFVector<GB,VT>* mutableDofVector_;
+  };
 
 } // end namespace AMDiS
 
+#include "DiscreteLocalFunction.inc.hpp"
 #include "DiscreteFunction.inc.hpp"

src/amdis/gridfunctions/DiscreteFunction.inc.hpp

 #pragma once
 
-#include <amdis/common/DerivativeTraits.hpp>
-#include <amdis/common/FieldMatVec.hpp>
-#include <amdis/utility/LocalBasisCache.hpp>
-#include <amdis/utility/LocalToGlobalAdapter.hpp>
+#include <type_traits>
 
-#include <dune/common/ftraits.hh>
 #include <dune/grid/utility/hierarchicsearch.hh>
 
-namespace AMDiS {
-
-template <class GB, class VT, class TP>
-class DiscreteFunction<GB,VT,TP>::LocalFunction
-{
-public:
-  using Domain = typename EntitySet::LocalCoordinate;
-  using Range = typename DiscreteFunction::Range;
-
-  enum { hasDerivative = true };
-
-private:
-  using LocalView = typename GlobalBasis::LocalView;
-  using Element = typename EntitySet::Element;
-  using Geometry = typename Element::Geometry;
-
-public:
-  /// Constructor. Stores a copy of the DiscreteFunction.
-  LocalFunction(DiscreteFunction const& globalFunction)
-    : globalFunction_(globalFunction)
-    , localView_(globalFunction_.basis()->localView())
-    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
-  {}
-
-  /// Copy constructor.
-  LocalFunction(LocalFunction const& other)
-    : globalFunction_(other.globalFunction_)
-    , localView_(globalFunction_.basis()->localView())
-    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
-  {}
-
-  /// \brief Bind the LocalView to the element
-  void bind(Element const& element)
-  {
-    localView_.bind(element);
-
-    globalFunction_.coefficients().gather(localView_, localCoefficients_);
-    bound_ = true;
-  }
-
-  /// \brief Unbind the LocalView from the element
-  void unbind()
-  {
-    localView_.unbind();
-    bound_ = false;
-  }
-
-  /// \brief Evaluate LocalFunction at bound element in local coordinates
-  Range operator()(Domain const& x) const;
-
-  /// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
-  GradientLocalFunction makeDerivative(tag::gradient type) const
-  {
-    return GradientLocalFunction{globalFunction_, type};
-  }
-
-  DivergenceLocalFunction makeDerivative(tag::divergence type) const
-  {
-    return DivergenceLocalFunction{globalFunction_, type};
-  }
-
-  PartialLocalFunction makeDerivative(tag::partial type) const
-  {
-    return PartialLocalFunction{globalFunction_, type};
-  }
-
-  /// \brief The \ref polynomialDegree() of the LocalFunctions
-  int order() const
-  {
-    assert( bound_ );
-    return polynomialDegree(*subTree_);
-  }
-
-  /// \brief Return the bound element
-  Element const& localContext() const
-  {
-    assert( bound_ );
-    return localView_.element();
-  }
-
-private:
-  DiscreteFunction globalFunction_;
-  LocalView localView_;
-  SubTree const* subTree_;
-
-  std::vector<VT> localCoefficients_;
-  bool bound_ = false;
-};
+#include <amdis/functions/Interpolate.hpp>
+#include <amdis/gridfunctions/GridFunction.hpp>
 
+namespace AMDiS {
 
+// Evaluate DiscreteFunction in global coordinates
 template <class GB, class VT, class TP>
-typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
-LocalFunction::operator()(Domain const& x) const
+typename DiscreteFunction<GB,VT,TP,true>::Range DiscreteFunction<GB,VT,TP,true>::
+  operator()(Domain const& x) const
 {
-  assert( bound_ );
-  Range y(0);
-
-  auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
-  for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
-  {
-    auto localBasisCache = makeNodeCache(node);
-    auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(), x);
-    std::size_t size = node.finiteElement().size();
-
-    // Get range entry associated to this node
-    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
-
-    for (std::size_t i = 0; i < size; ++i) {
-      // Get coefficient associated to i-th shape function
-      auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
-
-      // Get value of i-th shape function
-      auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
+  using Grid = typename GlobalBasis::GridView::Grid;
+  using IS = typename GlobalBasis::GridView::IndexSet;
 
-      std::size_t dimC = c.size();
-      std::size_t dimV = v.size();
-      assert(dimC*dimV == std::size_t(re.size()));
-      for(std::size_t j = 0; j < dimC; ++j) {
-        auto&& c_j = c[j];
-        for(std::size_t k = 0; k < dimV; ++k)
-          re[j*dimV + k] += c_j*v[k];
-      }
-    }
-  });
+  auto const& gv = this->basis()->gridView();
+  Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
 
-  return y;
+  auto element = hsearch.findEntity(x);
+  auto geometry = element.geometry();
+  auto localFct = localFunction(*this);
+  localFct.bind(element);
+  return localFct(geometry.local(x));
 }
 
 
+// Interpolation of GridFunction to DOFVector
 template <class GB, class VT, class TP>
-  template <class Type>
-class DiscreteFunction<GB,VT,TP>::DerivativeLocalFunctionBase
-{
-  using R = typename DiscreteFunction::Range;
-  using D = typename DiscreteFunction::Domain;
-  using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
-  using Traits = DerivativeTraits<RawSignature, Type>;
-
-public:
-  using Domain = typename EntitySet::LocalCoordinate;
-  using Range = typename Traits::Range;
-
-  enum { hasDerivative = false };
-
-private:
-  using LocalView = typename GlobalBasis::LocalView;
-  using Element = typename EntitySet::Element;
-  using Geometry = typename Element::Geometry;
-
-public:
-  /// Constructor. Stores a copy of the DiscreteFunction.
-  DerivativeLocalFunctionBase(DiscreteFunction const& globalFunction, Type const& type)
-    : globalFunction_(globalFunction)
-    , type_(type)