diff --git a/src/amdis/FileWriter.hpp b/src/amdis/FileWriter.hpp
index ff734622af2e36a5947b59f972ca8f3b2bdeb12f..6f6488af74291444aa7693d4a0935db3961f05ea 100644
--- a/src/amdis/FileWriter.hpp
+++ b/src/amdis/FileWriter.hpp
@@ -12,7 +12,7 @@
#include <amdis/Initfile.hpp>
#include <amdis/common/Size.hpp>
#include <amdis/common/ValueCategory.hpp>
-#include <amdis/gridfunctions/DOFVectorView.hpp>
+#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/io/FileWriterInterface.hpp>
#include <amdis/utility/Filesystem.hpp>
@@ -43,14 +43,17 @@ namespace AMDiS
constexpr std::size_t VTKFieldSize = Size<Range>;
- template <class GlobalBasis, class RangeType, class TreePath>
+ template <class GB, class VT, class TP>
class FileWriter
: public FileWriterInterface
{
- private: // typedefs and static constants
+ using GlobalBasis = GB;
+ using ValueType = VT;
+ using TreePath = TP;
+
+ private:
using GridView = typename GlobalBasis::GridView;
- using Vector = DOFVectorConstView<GlobalBasis,RangeType,TreePath>;
- using Range = typename Vector::Range;
+ using Range = typename DiscreteFunction<GB,VT,TP>::Range;
/// Dimension of the mesh
static constexpr int dim = GridView::dimension;
@@ -60,33 +63,32 @@ namespace AMDiS
public:
/// Constructor.
- FileWriter(std::string const& baseName,
- Vector const& dofvector)
- : FileWriterInterface(baseName)
- , dofvector_(dofvector)
+ FileWriter(std::string const& name, DiscreteFunction<GB,VT,TP> const& discreteFct)
+ : FileWriterInterface(name)
+ , discreteFct_(discreteFct)
, animation_(false)
{
- Parameters::get(baseName + "->ParaView animation", animation_);
+ Parameters::get(name + "->ParaView animation", animation_);
- int m = Parameters::get<int>(baseName + "->ParaView mode").value_or(0);
+ int m = Parameters::get<int>(name + "->ParaView mode").value_or(0);
mode_ = (m == 0 ? Dune::VTK::ascii : Dune::VTK::appendedraw);
- init(baseName, ValueCategory_t<Range>{});
+ init(name, ValueCategory_t<Range>{});
}
template <class ValueCat>
void init(std::string const& baseName, ValueCat)
{
- int subSampling = Parameters::get<int>(baseName + "->subsampling").value_or(0);
- if (subSampling > 0)
- vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(gridView(), subSampling);
+ auto subSampling = Parameters::get<int>(baseName + "->subsampling");
+ if (subSampling)
+ vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(gridView(), subSampling.value());
else
vtkWriter_ = std::make_shared<Dune::VTKWriter<GridView>>(gridView());
if (animation_)
vtkSeqWriter_ = std::make_shared<Dune::VTKSequenceWriter<GridView>>(vtkWriter_, filename_, dir_, "");
- vtkWriter_->addVertexData(dofvector_, Dune::VTK::FieldInfo(name_, VTKFieldType<Range>, VTKFieldSize<Range>));
+ vtkWriter_->addVertexData(discreteFct_, Dune::VTK::FieldInfo(name_, VTKFieldType<Range>, VTKFieldSize<Range>));
}
void init(std::string const&, tag::unknown) {}
@@ -105,11 +107,11 @@ namespace AMDiS
protected:
GridView const& gridView() const
{
- return dofvector_.basis().gridView();
+ return discreteFct_.basis().gridView();
}
private:
- Vector dofvector_;
+ DiscreteFunction<GB,VT,TP> discreteFct_;
std::shared_ptr<Dune::VTKWriter<GridView>> vtkWriter_;
std::shared_ptr<Dune::VTKSequenceWriter<GridView>> vtkSeqWriter_;
@@ -122,12 +124,12 @@ namespace AMDiS
};
- template <class GlobalBasis, class Range, class TreePath>
- std::shared_ptr<FileWriter<GlobalBasis,Range,TreePath>>
- makeFileWriterPtr(std::string baseName,
- DOFVectorConstView<GlobalBasis,Range,TreePath> const& dofvector)
+ template <class GlobalBasis, class ValueType, class TreePath>
+ std::shared_ptr<FileWriter<GlobalBasis,ValueType,TreePath>>
+ makeFileWriterPtr(std::string const& name,
+ DiscreteFunction<GlobalBasis,ValueType,TreePath> const& discreteFct)
{
- return std::make_shared<FileWriter<GlobalBasis,Range,TreePath>>(baseName, dofvector);
+ return std::make_shared<FileWriter<GlobalBasis,ValueType,TreePath>>(name, discreteFct);
}
} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 20c212d681e35ff0fa074066167796d893fed1ae..5a737a44d2aa2bc6c66bc6cfb8398455417b0aa2 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -32,6 +32,7 @@
#include <amdis/common/Utility.hpp>
#include <amdis/GridFunctions.hpp>
+#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/gridfunctions/DOFVectorView.hpp>
#include <amdis/io/FileWriterInterface.hpp>
@@ -198,7 +199,7 @@ namespace AMDiS
auto getSolution(TreePath const& path = {}) const
{
auto&& tp = makeTreePath(path);
- return makeDOFVectorView(*solution_, tp);
+ return makeDiscreteFunction(*solution_, tp);
}
diff --git a/src/amdis/gridfunctions/CMakeLists.txt b/src/amdis/gridfunctions/CMakeLists.txt
index a2b00db043f10c3680ffcfbb29aeb8a69fc6e4c7..328773dbb331542870e5923ecb608e05839620d7 100644
--- a/src/amdis/gridfunctions/CMakeLists.txt
+++ b/src/amdis/gridfunctions/CMakeLists.txt
@@ -5,8 +5,9 @@ install(FILES
ConstantGridFunction.hpp
CoordsGridFunction.hpp
DerivativeGridFunction.hpp
+ DiscreteFunction.hpp
+ DiscreteFunction.inc.hpp
DOFVectorView.hpp
- DOFVectorView.inc.hpp
FunctorGridFunction.hpp
GridFunctionConcepts.hpp
Integrate.hpp
diff --git a/src/amdis/gridfunctions/DOFVectorView.hpp b/src/amdis/gridfunctions/DOFVectorView.hpp
index 86b52e7e7771030973eeae3b3a83932e9672cd47..734c5e1ebc35037d2a13f252b35f9dae9b972e04 100644
--- a/src/amdis/gridfunctions/DOFVectorView.hpp
+++ b/src/amdis/gridfunctions/DOFVectorView.hpp
@@ -1,348 +1,132 @@
#pragma once
-#include <vector>
-
-#include <dune/common/std/optional.hh>
-#include <dune/functions/common/defaultderivativetraits.hh>
-#include <dune/functions/common/treedata.hh>
-#include <dune/functions/functionspacebases/defaultnodetorangemap.hh>
-#include <dune/functions/functionspacebases/flatvectorview.hh>
-#include <dune/functions/gridfunctions/gridviewentityset.hh>
-#include <dune/typetree/childextraction.hh>
-
#include <amdis/GridFunctions.hpp>
-#include <amdis/utility/FiniteElementType.hpp>
+#include <amdis/gridfunctions/DiscreteFunction.hpp>
namespace AMDiS
{
- /**
- * \addtogroup GridFunctions
- * @{
- **/
-
- template <class GlobalBasisType, class RangeType, class TreePathType>
- class DOFVectorConstView
+ /// 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>
{
- public:
- using GlobalBasis = GlobalBasisType;
- using TreePath = TreePathType;
-
- using Tree = typename GlobalBasis::LocalView::Tree;
- using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
- using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
-
- using GridView = typename GlobalBasis::GridView;
- using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
-
- using Domain = typename EntitySet::GlobalCoordinate;
- using Range = RangeType_t<SubTree>;
- static_assert(std::is_arithmetic<RangeType>::value, "");
- // Don't know how to determine Range with non-trivial RangeType
-
- using RawSignature = typename Dune::Functions::SignatureTraits<Range(Domain)>::RawSignature;
- using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
- using DerivativeRange = typename DerivativeTraits::Range;
-
- using LocalDomain = typename EntitySet::LocalCoordinate;
- using Element = typename EntitySet::Element;
- using Geometry = typename Element::Geometry;
-
- enum { hasDerivative = false };
-
- public: // a local view on the gradients
-
- /// A LocalFunction representing the derivative of the DOFVector
- class GradientLocalFunction
- {
- public:
- using Domain = LocalDomain;
- using Range = DerivativeRange;
-
- enum { hasDerivative = false };
-
- private:
- using LocalView = typename GlobalBasis::LocalView;
-
- template <class LeafNode>
- using LocalBasisJacobian = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
-
- template <class Node>
- using NodeData = typename std::vector<LocalBasisJacobian<Node>>;
-
- using ReferenceGradientContainer = Dune::Functions::TreeData<SubTree, NodeData, true>;
-
- public:
- GradientLocalFunction(DOFVectorConstView const& globalFunction)
- : globalFunction_(&globalFunction)
- , localView_(globalFunction_->basis().localView())
- , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_->treePath()))
- {
- referenceGradientContainer_.init(*subTree_);
- }
-
- void bind(Element const& element)
- {
- localView_.bind(element);
- geometry_.emplace(element.geometry());
- bound_ = true;
- }
-
- void unbind()
- {
- localView_.unbind();
- geometry_.reset();
- bound_ = false;
- }
-
- /// Evaluate Gradient at bound element in local coordinates
- Range operator()(Domain const& x) const;
-
- friend int order(GradientLocalFunction const& self)
- {
- assert( self.bound_ );
- return std::max(0, polynomialDegree(*self.subTree_)-1);
- }
-
- /// Return the bound element
- Element const& localContext() const
- {
- assert( bound_ );
- return localView_.element();
- }
-
- private:
- DOFVectorConstView const* globalFunction_;
- LocalView localView_;
-
- SubTree const* subTree_;
- mutable ReferenceGradientContainer referenceGradientContainer_;
-
- Dune::Std::optional<Geometry> geometry_;
- bool bound_ = false;
- };
-
-
- public: // a local view on the values
-
- /// A LocalFunction, i.e., an element local view on the DOFVector
- class LocalFunction
- {
- public:
- using Domain = typename DOFVectorConstView::LocalDomain;
- using Range = typename DOFVectorConstView::Range;
-
- enum { hasDerivative = true };
-
- private:
- using LocalView = typename GlobalBasis::LocalView;
-
- template <class LeafNode>
- using LocalBasisRange = RangeType_t<LeafNode>;
- // = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
-
- template <class Node>
- using NodeData = typename std::vector<LocalBasisRange<Node>>;
-
- using ShapeFunctionValueContainer = Dune::Functions::TreeData<SubTree, NodeData, true>;
-
- public:
- LocalFunction(DOFVectorConstView const& globalFunction)
- : globalFunction_(&globalFunction)
- , localView_(globalFunction_->basis().localView())
- , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_->treePath()))
- {
- shapeFunctionValueContainer_.init(*subTree_);
- }
-
- void bind(Element const& element)
- {
- localView_.bind(element);
- bound_ = true;
- }
-
- void unbind()
- {
- localView_.unbind();
- bound_ = false;
- }
-
- /// Evaluate LocalFunction at bound element in local coordinates
- Range operator()(Domain const& x) const;
-
- /// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
- friend GradientLocalFunction derivative(LocalFunction const& localFunction)
- {
- static_assert(isValidRange<DerivativeTraits>(),"Derivative of DOFVector not defined.");
- return GradientLocalFunction{*localFunction.globalFunction_};
- }
-
- friend int order(LocalFunction const& self)
- {
- assert( self.bound_ );
- return polynomialDegree(*self.subTree_);
- }
-
- /// Return the bound element
- Element const& localContext() const
- {
- assert( bound_ );
- return localView_.element();
- }
-
- private:
- DOFVectorConstView const* globalFunction_;
- LocalView localView_;
- SubTree const* subTree_;
-
- mutable ShapeFunctionValueContainer shapeFunctionValueContainer_;
-
- bool bound_ = false;
- };
+ using Self = DOFVectorView;
+ using Super = DiscreteFunction<GB, VT, TP>;
+ using GlobalBasis = GB;
+ using TreePath = TP;
public:
- /// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
- DOFVectorConstView(DOFVector<GlobalBasis,RangeType> const& dofVector, TreePath 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
+ DOFVectorView(DOFVector<GB,VT>& dofVector)
+ : DOFVectorView{dofVector, Dune::TypeTree::hybridTreePath()}
{}
- /// Evaluate the view on this DOFVector in global coordinates
- Range operator()(Domain const& x) const
- {
- error_exit("Not implemented.");
- return Range(0);
- }
-
- /// \brief Create a local function for this view on the DOFVector. \relates LocalFunction
- friend LocalFunction localFunction(DOFVectorConstView const& self)
- {
- return LocalFunction{self};
- }
-
- EntitySet const& entitySet() const
- {
- return entitySet_;
- }
-
- public:
- /// Return global basis
- GlobalBasis const& basis() const
- {
- return dofVector_->basis();
- }
-
- /// Return treePath associated with this view
- TreePath const& treePath() const
- {
- return treePath_;
- }
-
- /// Return const coefficient vector
- DOFVector<GlobalBasis,RangeType> const& coefficients() const
- {
- return *dofVector_;
- }
-
- protected:
- DOFVector<GlobalBasis,RangeType> const* dofVector_;
- TreePath const treePath_;
-
- EntitySet entitySet_;
- NodeToRangeEntry nodeToRangeEntry_;
- };
-
-
- // A mutable version of DOFVectorView
- template <class GlobalBasisType, class RangeType, class TreePathType>
- class DOFVectorMutableView
- : public DOFVectorConstView<GlobalBasisType, RangeType, TreePathType>
- {
- using Super = DOFVectorConstView<GlobalBasisType, RangeType, TreePathType>;
-
- using GlobalBasis = GlobalBasisType;
- using TreePath = TreePathType;
-
- public:
/// Constructor. Stores a pointer to the mutable `dofvector`.
- DOFVectorMutableView(DOFVector<GlobalBasis,RangeType>& dofVector, TreePath const& treePath)
+ DOFVectorView(DOFVector<GB,VT>& dofVector, TP const& treePath)
: Super(dofVector, treePath)
, mutableDofVector_(&dofVector)
{}
public:
- /// Interpolation of GridFunction to DOFVector
+ /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
+ /// reference to this DOFVector in the expression.
+ /**
+ * **Example:**
+ * ```
+ * auto v = makeDOFVectorView(prob.getSolutionVector(),0);
+ * v.interpolate([](auto const& x) { return x[0]; });
+ * ```
+ **/
template <class Expr>
- DOFVectorMutableView& interpolate(Expr&& expr)
+ void interpolate_noalias(Expr&& expr)
{
- auto const& basis = Super::basis();
- auto const& treePath = Super::treePath();
+ auto const& basis = this->basis();
+ auto const& treePath = this->treePath();
auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), basis.gridView());
- DOFVector<GlobalBasis,RangeType> tmp(*mutableDofVector_);
- Dune::Functions::interpolate(basis, treePath, tmp, std::forward<decltype(gridFct)>(gridFct));
+ Dune::Functions::interpolate(basis, treePath, coefficients(),
+ std::forward<decltype(gridFct)>(gridFct));
+ }
+
+ /// \brief Interpolation of GridFunction to DOFVector
+ /**
+ * **Example:**
+ * ```
+ * auto v = makeDOFVectorView(prob.getSolutionVector(),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>
+ void interpolate(Expr&& expr)
+ {
+ // create temporary copy of data
+ DOFVector<GB,VT> tmp(coefficients());
+ Self tmpView{tmp, this->treePath()};
+ tmpView.interpolate_noalias(std::forward<Expr>(expr));
// move data from temporary vector into stored DOFVector
- mutableDofVector_->vector() = std::move(tmp.vector());
+ coefficients().vector() = std::move(tmp.vector());
+ }
+
+ /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
+ template <class Expr>
+ DOFVectorView& operator<<(Expr&& expr)
+ {
+ interpolate(expr);
return *this;
}
+ /// \brief interpolate `(*this) + expr` to DOFVector
template <class Expr>
- DOFVectorMutableView& operator<<(Expr&& expr)
+ DOFVectorView& operator+=(Expr&& expr)
{
- return interpolate(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<GlobalBasis,RangeType>& coefficients() { return *mutableDofVector_; }
+ DOFVector<GB,VT>& coefficients() { return *mutableDofVector_; }
/// Return the const DOFVector
using Super::coefficients;
protected:
- DOFVector<GlobalBasis,RangeType>* mutableDofVector_;
+ DOFVector<GB,VT>* mutableDofVector_;
};
- /** @} **/
-
-
-#ifndef DOXYGEN
- // A Generator for a const \ref DOFVectorView.
- template <class GlobalBasis, class RangeType, class TreePath>
- auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType> const& dofVector, TreePath const& treePath)
- {
- return DOFVectorConstView<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
- }
-
- // A Generator for a mutable \ref DOFVectorView.
- template <class GlobalBasis, class RangeType, class TreePath>
- auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector, TreePath const& treePath)
- {
- return DOFVectorMutableView<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
- }
+#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
+ // 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 const \ref DOFVectorView.
- template <class GlobalBasis, class RangeType>
- auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType> const& dofVector)
+ /// A Generator for a mutable \ref DOFVectorView
+ template <class GlobalBasis, class ValueType, class TreePath>
+ auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector, TreePath const& treePath)
{
- auto treePath = Dune::TypeTree::hybridTreePath();
- return DOFVectorConstView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
+ return DOFVectorView<GlobalBasis, ValueType, TreePath>{dofVector, treePath};
}
- // A Generator for a mutable \ref DOFVectorView.
- template <class GlobalBasis, class RangeType>
- auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector)
+ /// A Generator for a mutable \ref DOFVectorView
+ template <class GlobalBasis, class ValueType>
+ auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
- return DOFVectorMutableView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
+ return DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath};
}
-#endif
} // end namespace AMDiS
-
-#include "DOFVectorView.inc.hpp"
diff --git a/src/amdis/gridfunctions/DOFVectorView.inc.hpp b/src/amdis/gridfunctions/DOFVectorView.inc.hpp
deleted file mode 100644
index 939b89a015e10a60317f97599ff5a834dcd38d8e..0000000000000000000000000000000000000000
--- a/src/amdis/gridfunctions/DOFVectorView.inc.hpp
+++ /dev/null
@@ -1,109 +0,0 @@
-#pragma once
-
-#include <amdis/common/FieldMatVec.hpp>
-#include <amdis/LocalBasisCache.hpp>
-
-namespace AMDiS {
-
-template <class GB, class RT, class TP>
-typename DOFVectorConstView<GB,RT,TP>::Range DOFVectorConstView<GB,RT,TP>::
-LocalFunction::operator()(LocalDomain const& x) const
-{
- assert( bound_ );
- auto y = Range(0);
-
- auto&& coefficients = *globalFunction_->dofVector_;
- auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
- forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
- {
- auto&& fe = node.finiteElement();
- auto&& localBasis = fe.localBasis();
-
- LocalBasisCache<std::decay_t<decltype(localBasis)>> localBasisCache(localBasis);
- auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element(),x);
-
- // Get range entry associated to this node
- auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
-
- for (std::size_t i = 0; i < localBasis.size(); ++i) {
- auto&& multiIndex = localView_.index(node.localIndex(i));
-
- // Get coefficient associated to i-th shape function
- auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
-
- // Get value of i-th shape function
- auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
-
- 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];
- }
- }
- });
-
- return y;
-}
-
-
-template <class GB, class RT, class TP>
-typename DOFVectorConstView<GB,RT,TP>::DerivativeRange DOFVectorConstView<GB,RT,TP>::
-GradientLocalFunction::operator()(LocalDomain const& x) const
-{
- assert( bound_ );
- Range dy;
- for (std::size_t j = 0; j < dy.size(); ++j)
- dy[j] = 0;
-
- auto&& coefficients = *globalFunction_->dofVector_;
- auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
- forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
- {
- // TODO: may DOFVectorView::Range to FieldVector type if necessary
- using LocalDerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
- using GradientBlock = typename LocalDerivativeTraits::Range;
-
- auto&& fe = node.finiteElement();
- auto&& localBasis = fe.localBasis();
-
- LocalBasisCache<std::decay_t<decltype(localBasis)>> localBasisCache(localBasis);
- auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element(),x);
-
- // The transposed inverse Jacobian of the map from the reference element to the element
- auto&& jacobian = geometry_.value().jacobianInverseTransposed(x);
-
- // Compute the shape function gradients on the real element
- std::vector<GradientBlock> gradients(referenceGradients.size());
- for (std::size_t i = 0; i < gradients.size(); ++i)
- multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad
-
- // Get range entry associated to this node
- auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
-
- for (std::size_t i = 0; i < localBasis.size(); ++i) {
- auto&& multiIndex = localView_.index(node.localIndex(i));
-
- // Get coefficient associated to i-th shape function
- auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
-
- // Get value of i-th transformed reference gradient
- auto grad = Dune::Functions::flatVectorView(gradients[i]);
-
- std::size_t dimC = c.size();
- std::size_t dimV = grad.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*grad[k];
- }
- }
- });
-
- return dy;
-}
-
-} // end namespace AMDiS
diff --git a/src/amdis/gridfunctions/DerivativeGridFunction.hpp b/src/amdis/gridfunctions/DerivativeGridFunction.hpp
index e4c0adc7f64304342260e00b5f0aa3d5f4bee75c..c9b3811d32947e7c0beaeac09980024874c58864 100644
--- a/src/amdis/gridfunctions/DerivativeGridFunction.hpp
+++ b/src/amdis/gridfunctions/DerivativeGridFunction.hpp
@@ -3,6 +3,7 @@
#include <type_traits>
#include <dune/functions/common/defaultderivativetraits.hh>
+#include <dune/grid/utility/hierarchicsearch.hh>
#include <amdis/gridfunctions/GridFunctionConcepts.hpp>
@@ -52,11 +53,22 @@ namespace AMDiS
static_assert(isValidRange<DerivativeTraits>(), "Derivative of GridFunction not defined");
}
- /// NOTE: no global derivative available
+ /// Evaluate derivative in global coordinates. NOTE: expensive
Range operator()(Domain const& x) const
{
- error_exit("Not implemented");
- return Range(0);
+ auto gv = entitySet().gridView();
+
+ using GridView = decltype(gv);
+ using Grid = typename GridView::Grid;
+ using IS = typename GridView::IndexSet;
+
+ Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
+
+ auto element = hsearch.findEntity(x);
+ auto geometry = element.geometry();
+ auto localFct = derivative(localFunction(gridFct_));
+ localFct.bind(element);
+ return localFct(geometry.local(x));
}
/// Return the derivative-localFunction of the GridFunction.
diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..71b2045e3fc8796c4ca50667e27731a62734d40f
--- /dev/null
+++ b/src/amdis/gridfunctions/DiscreteFunction.hpp
@@ -0,0 +1,141 @@
+#pragma once
+
+#include <vector>
+
+#include <dune/common/std/optional.hh>
+#include <dune/functions/common/defaultderivativetraits.hh>
+#include <dune/functions/functionspacebases/defaultnodetorangemap.hh>
+#include <dune/functions/functionspacebases/flatvectorview.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/typetree/childextraction.hh>
+
+#include <amdis/GridFunctions.hpp>
+#include <amdis/utility/FiniteElementType.hpp>
+
+namespace AMDiS
+{
+ /// \brief A view on a subspace of a \ref DOFVector
+ /**
+ * \ingroup GridFunctions
+ *
+ * \tparam GB GlobalBasis type that models \ref Dune::Functions::Concept::GlobalBasis
+ * \tparam VT Coefficient type of the DOFVector
+ * \tparam TP TreePath type a realization of \ref Dune::TypeTree::HybridTreePath
+ **/
+ template <class GB, class VT, class TP>
+ class DiscreteFunction
+ {
+ static_assert(std::is_arithmetic<VT>::value, "");
+
+ private:
+ using GlobalBasis = GB;
+ using TreePath = TP;
+
+ using Tree = typename GlobalBasis::LocalView::Tree;
+ using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
+ using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
+
+ using GridView = typename GlobalBasis::GridView;
+
+ public:
+ /// Set of entities the DiscreteFunction is defined on
+ using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
+
+ /// Global coordinates of the EntitySet
+ using Domain = typename EntitySet::GlobalCoordinate;
+
+ /// Range type of this DiscreteFunction
+ using Range = RangeType_t<SubTree,VT>;
+
+ /// \brief This GridFunction has no derivative function, it can be created
+ /// by \ref DiscreteGridFunction.
+ enum { hasDerivative = false };
+
+ public:
+ /// A LocalFunction representing the derivative of the DOFVector on a bound element
+ class GradientLocalFunction;
+
+ /// A LocalFunction representign the value the DOFVector on a bound element
+ class LocalFunction;
+
+ public:
+ /// Constructor forwards to the treePath constructor, with empty TreePath
+ DiscreteFunction(DOFVector<GB,VT> const& dofVector)
+ : DiscreteFunction{dofVector, Dune::TypeTree::hybridTreePath()}
+ {}
+
+ /// Constructor. Stores a pointer to the dofVector and a copy of the 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))
+ {}
+
+ /// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive
+ Range operator()(Domain const& x) const;
+
+ /// \brief Create a local function for this view on the DOFVector. \relates LocalFunction
+ friend LocalFunction localFunction(DiscreteFunction const& self)
+ {
+ return LocalFunction{self};
+ }
+
+ /// \brief Return a \ref Dune::Functions::GridViewEntitySet
+ EntitySet const& entitySet() const
+ {
+ return entitySet_;
+ }
+
+ public:
+ /// \brief Return global basis bound to the DOFVector
+ GlobalBasis const& basis() const
+ {
+ return dofVector_->basis();
+ }
+
+ /// \brief Return treePath associated with this discrete function
+ TreePath const& treePath() const
+ {
+ return treePath_;
+ }
+
+ /// \brief Return const coefficient vector
+ DOFVector<GB,VT> const& coefficients() const
+ {
+ return *dofVector_;
+ }
+
+ protected:
+ DOFVector<GB,VT> const* dofVector_;
+ TreePath treePath_;
+ EntitySet entitySet_;
+ NodeToRangeEntry nodeToRangeEntry_;
+ };
+
+
+#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 TreePath>
+ auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector, TreePath const& treePath)
+ {
+ return DiscreteFunction<GlobalBasis, ValueType, TreePath>{dofVector, treePath};
+ }
+
+ /// 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};
+ }
+
+} // end namespace AMDiS
+
+#include "DiscreteFunction.inc.hpp"
diff --git a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..730e4b82d30dd55733eb438efa376077ff4a83c4
--- /dev/null
+++ b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
@@ -0,0 +1,260 @@
+#pragma once
+
+#include <amdis/common/FieldMatVec.hpp>
+#include <amdis/LocalBasisCache.hpp>
+
+#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:
+ LocalFunction(DiscreteFunction const& globalFunction)
+ : globalFunction_(&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);
+ 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
+ friend GradientLocalFunction derivative(LocalFunction const& localFunction)
+ {
+ return GradientLocalFunction{*localFunction.globalFunction_};
+ }
+
+ /// \brief The \ref polynomialDegree() of the LocalFunctions
+ friend int order(LocalFunction const& self)
+ {
+ assert( self.bound_ );
+ return polynomialDegree(*self.subTree_);
+ }
+
+ /// \brief Return the bound element
+ Element const& localContext() const
+ {
+ assert( bound_ );
+ return localView_.element();
+ }
+
+private:
+ DiscreteFunction const* globalFunction_;
+ LocalView localView_;
+ SubTree const* subTree_;
+ bool bound_ = false;
+};
+
+
+template <class GB, class VT, class TP>
+class DiscreteFunction<GB,VT,TP>::GradientLocalFunction
+{
+ using R = typename DiscreteFunction::Range;
+ using D = typename DiscreteFunction::Domain;
+ using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
+ using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
+
+public:
+ using Domain = typename EntitySet::LocalCoordinate;
+ using Range = typename DerivativeTraits::Range;
+
+ enum { hasDerivative = false };
+
+private:
+ using LocalView = typename GlobalBasis::LocalView;
+ using Element = typename EntitySet::Element;
+ using Geometry = typename Element::Geometry;
+
+public:
+ GradientLocalFunction(DiscreteFunction const& globalFunction)
+ : globalFunction_(&globalFunction)
+ , localView_(globalFunction_->basis().localView())
+ , subTree_(&child(localView_.tree(), globalFunction_->treePath()))
+ {}
+
+ void bind(Element const& element)
+ {
+ localView_.bind(element);
+ geometry_.emplace(element.geometry());
+ bound_ = true;
+ }
+
+ void unbind()
+ {
+ localView_.unbind();
+ geometry_.reset();
+ bound_ = false;
+ }
+
+ /// Evaluate Gradient at bound element in local coordinates
+ Range operator()(Domain const& x) const;
+
+ friend int order(GradientLocalFunction const& self)
+ {
+ assert( self.bound_ );
+ return std::max(0, polynomialDegree(*self.subTree_)-1);
+ }
+
+ /// Return the bound element
+ Element const& localContext() const
+ {
+ assert( bound_ );
+ return localView_.element();
+ }
+
+private:
+ DiscreteFunction const* globalFunction_;
+ LocalView localView_;
+ SubTree const* subTree_;
+ Dune::Std::optional<Geometry> geometry_;
+ bool bound_ = false;
+};
+
+
+template <class GB, class VT, class TP>
+typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
+LocalFunction::operator()(Domain const& x) const
+{
+ assert( bound_ );
+ auto y = Range(0);
+
+ auto&& coefficients = *globalFunction_->dofVector_;
+ auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
+ forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
+ {
+ auto&& fe = node.finiteElement();
+ auto&& localBasis = fe.localBasis();
+
+ LocalBasisCache<std::decay_t<decltype(localBasis)>> localBasisCache(localBasis);
+ auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element(),x);
+
+ // Get range entry associated to this node
+ auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
+
+ for (std::size_t i = 0; i < localBasis.size(); ++i) {
+ auto&& multiIndex = localView_.index(node.localIndex(i));
+
+ // Get coefficient associated to i-th shape function
+ auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
+
+ // Get value of i-th shape function
+ auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
+
+ 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];
+ }
+ }
+ });
+
+ return y;
+}
+
+
+template <class GB, class VT, class TP>
+typename DiscreteFunction<GB,VT,TP>::GradientLocalFunction::Range DiscreteFunction<GB,VT,TP>::
+GradientLocalFunction::operator()(Domain const& x) const
+{
+ assert( bound_ );
+ Range dy;
+ for (std::size_t j = 0; j < dy.size(); ++j)
+ dy[j] = 0;
+
+ auto&& coefficients = *globalFunction_->dofVector_;
+ auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
+ forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
+ {
+ // TODO: may DOFVectorView::Range to FieldVector type if necessary
+ using LocalDerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
+ using GradientBlock = typename LocalDerivativeTraits::Range;
+
+ auto&& fe = node.finiteElement();
+ auto&& localBasis = fe.localBasis();
+
+ LocalBasisCache<std::decay_t<decltype(localBasis)>> localBasisCache(localBasis);
+ auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element(),x);
+
+ // The transposed inverse Jacobian of the map from the reference element to the element
+ auto&& jacobian = geometry_.value().jacobianInverseTransposed(x);
+
+ // Compute the shape function gradients on the real element
+ std::vector<GradientBlock> gradients(referenceGradients.size());
+ for (std::size_t i = 0; i < gradients.size(); ++i)
+ multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad
+
+ // Get range entry associated to this node
+ auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
+
+ for (std::size_t i = 0; i < localBasis.size(); ++i) {
+ auto&& multiIndex = localView_.index(node.localIndex(i));
+
+ // Get coefficient associated to i-th shape function
+ auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
+
+ // Get value of i-th transformed reference gradient
+ auto grad = Dune::Functions::flatVectorView(gradients[i]);
+
+ std::size_t dimC = c.size();
+ std::size_t dimV = grad.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*grad[k];
+ }
+ }
+ });
+
+ return dy;
+}
+
+
+template <class GB, class VT, class TP>
+typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
+operator()(Domain const& x) const
+{
+ using Grid = typename GlobalBasis::GridView::Grid;
+ using IS = typename GlobalBasis::GridView::IndexSet;
+
+ auto const& gv = this->basis().gridView();
+ Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
+
+ auto element = hsearch.findEntity(x);
+ auto geometry = element.geometry();
+ auto localFct = localFunction(*this);
+ localFct.bind(element);
+ return localFct(geometry.local(x));
+}
+
+} // end namespace AMDiS
diff --git a/src/amdis/utility/RangeType.hpp b/src/amdis/utility/RangeType.hpp
index cee8d959e970b86769290983aa38c2797ec0d86d..27116104982776775ae39e9369c5d117623c9536 100644
--- a/src/amdis/utility/RangeType.hpp
+++ b/src/amdis/utility/RangeType.hpp
@@ -11,31 +11,43 @@ namespace AMDiS
{
namespace Impl
{
- template <class Node, bool isLeaf, bool isPower, bool isComposite>
+ template <class Node, class R, bool isLeaf, bool isPower, bool isComposite>
struct RangeTypeImpl
{
static_assert( isLeaf || isPower || isComposite, "Unknown node-type for range definition" );
};
}
- template <class Node>
+ /// \brief Range type of a node in the basis tree, composed of the leaf basis
+ /// range types.
+ /**
+ * Generate the range type by recursively combining leaf range types to a
+ * hybrid node range type. Range types for PowerNodes are thereby constructed
+ * as Dune::FieldVector of child range types. CompositeNodes produce a
+ * MultiTypeVector of the difference child range types.
+ *
+ * \tparam Node Type of a basis-tree node
+ * \tparam R Coefficient type [double]
+ **/
+ template <class Node, class R = double>
using RangeType_t =
- typename Impl::RangeTypeImpl<Node, Node::isLeaf, Node::isPower, Node::isComposite>::type;
+ typename Impl::RangeTypeImpl<Node, R, Node::isLeaf, Node::isPower, Node::isComposite>::type;
namespace Impl
{
// Leaf node
- template <class Node>
- struct RangeTypeImpl<Node, true, false, false>
+ template <class Node, class R>
+ struct RangeTypeImpl<Node, R, true, false, false>
{
using LocalBasis = typename Node::FiniteElement::Traits::LocalBasisType;
- using type = typename LocalBasis::Traits::RangeType;
+ using T = typename LocalBasis::Traits::RangeType;
+ using type = std::decay_t<decltype(std::declval<R>() * std::declval<T>())>;
};
// Power node
- template <class Node>
- struct RangeTypeImpl<Node, false, true, false>
+ template <class Node, class R>
+ struct RangeTypeImpl<Node, R, false, true, false>
{
using ChildNode = typename Node::template Child<0>::type;
@@ -49,12 +61,12 @@ namespace AMDiS
using type = Dune::FieldVector<T, int(Node::CHILDREN)>;
};
- using type = typename FlatType<ChildNode::isLeaf, RangeType_t<ChildNode>>::type;
+ using type = typename FlatType<ChildNode::isLeaf, RangeType_t<ChildNode,R>>::type;
};
// Composite node
- template <class Node>
- struct RangeTypeImpl<Node, false, false, true>
+ template <class Node, class R>
+ struct RangeTypeImpl<Node, R, false, false, true>
{
template <class Idx> struct RangeTypeGenerator;
template <std::size_t... I>
@@ -62,7 +74,7 @@ namespace AMDiS
{
template <std::size_t J>
using ChildNode = typename Node::template Child<J>::type;
- using type = MultiTypeVector<RangeType_t<ChildNode<I>>...>;
+ using type = MultiTypeVector<RangeType_t<ChildNode<I>,R>...>;
};
using type = typename RangeTypeGenerator<MakeSeq_t<Node::CHILDREN>>::type;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 88b9c3065228b2ad046d74d12d152311399ec091..57754d5dc8a495c3869ad293c8482b29928803de 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -11,6 +11,10 @@ dune_add_test(SOURCES ConceptsTest.cpp
dune_add_test(SOURCES DOFVectorTest.cpp
LINK_LIBRARIES amdis)
+dune_add_test(SOURCES DiscreteFunctionTest.cpp
+ LINK_LIBRARIES amdis
+ CMD_ARGS "${CMAKE_SOURCE_DIR}/examples/init/ellipt.dat.2d")
+
dune_add_test(SOURCES ExpressionsTest.cpp
LINK_LIBRARIES amdis)
diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f390af6ea44eb63d875770f46b22bce15768193e
--- /dev/null
+++ b/test/DiscreteFunctionTest.cpp
@@ -0,0 +1,107 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <iostream>
+#include <type_traits>
+
+#include <amdis/AMDiS.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/Operators.hpp>
+#include <amdis/common/Literals.hpp>
+#include <amdis/common/FieldMatVec.hpp>
+#include <amdis/gridfunctions/Integrate.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+using ElliptParam = YaspGridBasis<2, 2>;
+using ElliptProblem = ProblemStat<ElliptParam>;
+
+template <class GB, class T>
+bool comp(DOFVector<GB,T> const& U, DOFVector<GB,T> const& V)
+{
+ if (U.size() != V.size())
+ return false;
+
+ using Int = typename DOFVector<GB,T>::size_type;
+ for (Int i = 0; i < U.size(); ++i) {
+ if (U[i] != V[i])
+ return false;
+ }
+ return true;
+}
+
+int main(int argc, char** argv)
+{
+ AMDiS::init(argc, argv);
+
+ using namespace Dune::Indices;
+
+ ElliptProblem prob("ellipt");
+ prob.initialize(INIT_ALL);
+
+ // create 3 copies of the solution vector
+ auto U0 = prob.getSolutionVector();
+ auto U1 = prob.getSolutionVector();
+ auto U2 = prob.getSolutionVector();
+
+ auto u0 = makeDOFVectorView(U0);
+ auto u1 = makeDOFVectorView(U1);
+ auto u2 = makeDOFVectorView(U2);
+
+ auto expr = [](auto const& x) { return x[0] + x[1]; };
+ u0.interpolate_noalias(expr);
+ u1.interpolate(expr);
+ u2 << expr;
+
+ AMDIS_TEST( comp(U0, U1) );
+ AMDIS_TEST( comp(U0, U2) );
+
+ auto expr2 = [](auto const& x) { return 2*x[0] - 3*x[1]; };
+
+ u0.interpolate_noalias(u2 + expr2);
+ u1.interpolate(u1 + expr2);
+ u2 += expr2;
+
+ AMDIS_TEST( comp(U0, U1) );
+ AMDIS_TEST( comp(U0, U2) );
+
+ auto expr3 = [](auto const& x) { return -0.5*x[0] - 2*x[1]; };
+
+ u0.interpolate_noalias(u2 - expr3);
+ u1.interpolate(u1 - expr3);
+ u2 -= expr3;
+
+ AMDIS_TEST( comp(U0, U1) );
+ AMDIS_TEST( comp(U0, U2) );
+
+ auto du0 = derivative(u0);
+ auto localFct = localFunction(u0);
+ auto dlocalFct = derivative(localFct);
+ for (auto const& e : elements(prob.gridView()))
+ {
+ localFct.bind(e);
+ auto geo = e.geometry();
+ auto local = referenceElement(geo).position(0,0);
+ auto y0 = localFct(local);
+ auto y1 = u0(geo.global(local));
+
+ AMDIS_TEST(std::abs(y0 - y1) < 1.e-10);
+ localFct.unbind();
+
+ dlocalFct.bind(e);
+ auto g0 = dlocalFct(local);
+ auto g1 = du0(geo.global(local));
+
+ AMDIS_TEST(two_norm(g0 - g1) < 1.e-10);
+ dlocalFct.unbind();
+ }
+
+ auto V0 = makeDOFVector<float>(prob.globalBasis());
+ auto v0 = makeDOFVectorView(V0);
+ v0 << expr;
+
+ AMDiS::finalize();
+ return 0;
+}