Reimplementation of DiscreteGlobalBasisFunction

dune/amdis/DOFVectorView.hpp

+#pragma once
+
+namespace AMDiS
+{
+  template <class GlobalBasisType, class TreePathType>
+  class DOFVectorView
+  {
+  public:
+    using GlobalBasis = GlobalBasisType;
+    using TreePath = TreePathType;
+    using Vector = DOFVector<GlobalBasis>;
+
+    using Tree = typename GlobalBasis::LocalView::Tree;
+    using SubTree = typename TypeTree::ChildForTreePath<Tree, TreePath>;
+    using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>
+
+    using GridView = typename GlobalBasis::GridView;
+    using EntitySet = GridViewEntitySet<GridView, 0>;
+
+    using Domain = typename EntitySet::GlobalCoordinate;
+    using Range = typename Vector::value_type;
+
+    using LocalDomain = typename EntitySet::LocalCoordinate;
+    using Element = typename EntitySet::Element;
+
+    class LocalFunction
+    {
+    public:
+      using Domain = typename DOFVectorView::LocalDomain;
+      using Range = typename DOFVectorView::Range;
+      using Element = typename DOFVectorView::Element;
+
+      using LocalBasisView = typename GlobalBasis::LocalView;
+      using LocalIndexSet = typename GlobalBasis::LocalIndexSet;
+
+      template <class LeafNode>
+      using LocalBasisRange = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
+
+      template <class Node>
+      using NodeData = typename std::vector<LocalBasisRange<Node>>;
+
+      template <class Block>
+      using Flat = Dune::Functions::FlatVectorBackend<Block>;
+
+      using ShapeFunctionValueContainer = Dune::Functions::TreeData<SubTree, NodeData, true>;
+
+    public:
+      LocalFunction(DOFVectorView const& globalFunction)
+        : globalFunction_(&globalFunction)
+        , localBasisView_(globalFunction_->basis().localView())
+        , localIndexSet_(globalFunction_->basis().localIndexSet())
+        , subTree_(TypeTree::child(localBasisView_.tree(), globalFunction_->treePath())
+      {
+        shapeFunctionValueContainer_.init(subTree_);
+      }
+
+      void bind(Element const& element)
+      {
+        localBasisView_.bind(element);
+        localIndexSet_.bind(localBasisView_);
+        bound_ = true;
+      }
+
+      void unbind()
+      {
+        localIndexSet_.unbind();
+        localBasisView_.unbind();
+        bound_ = false;
+      }
+
+      /**
+      * \brief Evaluate LocalFunction at bound element.
+      *
+      * The result of this method is undefined if you did
+      * not call bind() beforehand or changed the coefficient
+      * vector after the last call to bind(). In the latter case
+      * you have to call bind() again in order to make operator()
+      * usable.
+      */
+      Range operator()(Domain const& x) const;
+
+      Element const& localContext() const
+      {
+        return localBasisView_.element();
+      }
+
+    private:
+      DOFVectorView const* globalFunction_;
+      LocalBasisView localBasisView_;
+      LocalIndexSet localIndexSet_;
+
+      mutable ShapeFunctionValueContainer shapeFunctionValueContainer_;
+      SubTree const& subTree_;
+
+      bool bound_ = false;
+    };
+
+  public:
+
+    DOFVectorView(DOFVector<GlobalBasis>& dofVector, TreePath const& treePath)
+      : basis_(basis)
+      , treePath_(treePath)
+      , coefficients_(coefficients)
+      , entitySet_(basis.gridView())
+      , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(basis, treePath))
+    {}
+
+    Range operator()(Domain const& x) const
+    {
+      error_exit("Not implemented.");
+      return Range(0);
+    }
+
+    friend LocalFunction localFunction(DOFVectorView const& self)
+    {
+      return LocalExpr{self};
+    }
+
+    EntitySet const& entitySet() const
+    {
+      return entitySet_;
+    }
+
+
+  public:
+
+    template <class PreExpr>
+    DOFVectorView& operator<<(PreExpr const& preExpr)
+    {
+      typename ToTerm<PreExpr>::type expr = toTerm(std::forward<PreExpr>(preExpr));
+
+      auto gridViewExpr = makeGridViewExpression(expr, basis_.gridView());
+      DOFVector<GlobalBasis> tmp(basis(), "tmp");
+      Dune::Functions::interpolate(basis_, treePath_, tmp, gridViewExpr);
+
+      coefficients_.getVector() = std::move(vec.getVector());
+      return *this;
+    }
+
+
+  public:
+
+    GlobalBasis const& basis() const { return dofVector_.getFeSpace(); }
+    TreePath const& treePath() const { return treePath_; }
+
+    DOFVector<GlobalBasis>& coefficients() { return dofVector_; }
+    DOFVector<GlobalBasis> const& coefficients() const { return dofVector_; }
+
+
+  private:
+
+    DOFVector<GlobalBasis>& dofVector_;
+    TreePath const treePath_;
+
+    EntitySet entitySet_;
+    NodeToRangeEntry nodeToRangeEntry_;
+  };
+
+} // end namespace AMDiS
+
+#include "DOFVectorView.inc.hpp"

dune/amdis/DOFVectorView.inc.hpp

+#pragma once
+
+namespace AMDiS
+{
+
+      {
+        auto y = Range(0);
+
+        auto const& coefficients = globalFunction_.coefficients();
+        auto const& nodeToRangeEntry = globalFunction_.nodeToRangeEntry();
+        forEachLeafNode(subTree_, [&,this](auto const& node, auto)
+        {
+          using Node = std::decay_t<decltype(node)>;
+          using LocalBasisRange = typename LocalFunction::template LocalBasisRange<Node>;
+          using MultiIndex = typename LocalIndexSet::MultiIndex;
+          using CoefficientBlock = Range;
+          using RangeBlock = std::decay_t<decltype(nodeToRangeEntry(node, y))>;
+
+          auto&& fe = node.finiteElement();
+          auto&& localBasis = fe.localBasis();
+
+          auto&& shapeFunctionValues = shapeFunctionValueContainer_[node];
+          localBasis.evaluateFunction(x, shapeFunctionValues);
+
+          // Get range entry associated to this node
+          auto&& re = nodeToRangeEntry(node, y);
+
+          for (std::size_t i = 0; i < localBasis.size(); ++i) {
+            auto&& multiIndex = localIndexSet_.index(node.localIndex(i));
+
+            // Get coefficient associated to i-th shape function
+            auto&& c = coefficients[multiIndex];
+
+            // Get value of i-th shape function
+            auto&& v = shapeFunctionValues[i];
+
+            // Notice that the range entry re, the coefficient c, and the shape functions
+            // value v may all be scalar, vector, matrix, or general container valued.
+            auto dimC = Flat<CoefficientBlock>::size(c);
+            auto dimV = Flat<LocalBasisRange>::size(v);
+            assert(dimC*dimV == Flat<RangeBlock>::size(re));
+            for(std::size_t j = 0; j < dimC; ++j) {
+              auto&& c_j = Flat<CoefficientBlock>::getEntry(c, j);
+              for(std::size_t k = 0; k < dimV; ++k) {
+                auto&& v_k = Flat<LocalBasisRange>::getEntry(v, k);
+                Flat<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k;
+              }
+            }
+          }
+        });
+
+        return y;
+      }
+
+
+} // end namespace AMDiS