From a4aeae1da5363195f6d07651daad4de63fa5c52c Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Mon, 22 Oct 2018 16:48:41 -0400
Subject: [PATCH] renamed DOFVectorConstView into DiscreteFunction and added
some more interpolation methods
---
src/amdis/FileWriter.hpp | 20 +-
src/amdis/ProblemStat.hpp | 3 +-
src/amdis/gridfunctions/CMakeLists.txt | 3 +-
src/amdis/gridfunctions/DOFVectorView.hpp | 332 +++---------------
.../gridfunctions/DerivativeGridFunction.hpp | 18 +-
src/amdis/gridfunctions/DiscreteFunction.hpp | 245 +++++++++++++
...rView.inc.hpp => DiscreteFunction.inc.hpp} | 24 +-
test/CMakeLists.txt | 4 +
test/DiscreteFunctionTest.cpp | 102 ++++++
9 files changed, 445 insertions(+), 306 deletions(-)
create mode 100644 src/amdis/gridfunctions/DiscreteFunction.hpp
rename src/amdis/gridfunctions/{DOFVectorView.inc.hpp => DiscreteFunction.inc.hpp} (83%)
create mode 100644 test/DiscreteFunctionTest.cpp
diff --git a/src/amdis/FileWriter.hpp b/src/amdis/FileWriter.hpp
index ff734622..a8d95f25 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>
@@ -49,8 +49,8 @@ namespace AMDiS
{
private: // typedefs and static constants
using GridView = typename GlobalBasis::GridView;
- using Vector = DOFVectorConstView<GlobalBasis,RangeType,TreePath>;
- using Range = typename Vector::Range;
+ using DiscreteFunction = AMDiS::DiscreteFunction<GlobalBasis,RangeType,TreePath>;
+ using Range = typename DiscreteFunction::Range;
/// Dimension of the mesh
static constexpr int dim = GridView::dimension;
@@ -61,9 +61,9 @@ namespace AMDiS
public:
/// Constructor.
FileWriter(std::string const& baseName,
- Vector const& dofvector)
+ DiscreteFunction const& discreteFct)
: FileWriterInterface(baseName)
- , dofvector_(dofvector)
+ , discreteFct_(discreteFct)
, animation_(false)
{
Parameters::get(baseName + "->ParaView animation", animation_);
@@ -86,7 +86,7 @@ namespace AMDiS
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 +105,11 @@ namespace AMDiS
protected:
GridView const& gridView() const
{
- return dofvector_.basis().gridView();
+ return discreteFct_.basis().gridView();
}
private:
- Vector dofvector_;
+ DiscreteFunction discreteFct_;
std::shared_ptr<Dune::VTKWriter<GridView>> vtkWriter_;
std::shared_ptr<Dune::VTKSequenceWriter<GridView>> vtkSeqWriter_;
@@ -125,9 +125,9 @@ 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)
+ DiscreteFunction<GlobalBasis,Range,TreePath> const& discreteFct)
{
- return std::make_shared<FileWriter<GlobalBasis,Range,TreePath>>(baseName, dofvector);
+ return std::make_shared<FileWriter<GlobalBasis,Range,TreePath>>(baseName, discreteFct);
}
} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 74eb5c10..a965675f 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 a2b00db0..328773db 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 86b52e7e..c39cf8df 100644
--- a/src/amdis/gridfunctions/DOFVectorView.hpp
+++ b/src/amdis/gridfunctions/DOFVectorView.hpp
@@ -1,299 +1,75 @@
#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
- * @{
- **/
-
+ /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
template <class GlobalBasisType, class RangeType, class TreePathType>
- class DOFVectorConstView
+ class DOFVectorView
+ : public DiscreteFunction<GlobalBasisType, RangeType, TreePathType>
{
- public:
+ using Self = DOFVectorView;
+ using Super = DiscreteFunction<GlobalBasisType, RangeType, TreePathType>;
+
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;
- };
-
-
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. Stores a pointer to the mutable `dofvector`.
+ DOFVectorView(DOFVector<GlobalBasis,RangeType>& dofVector, TreePath const& treePath)
+ : Super(dofVector, treePath)
+ , mutableDofVector_(&dofVector)
{}
- /// Evaluate the view on this DOFVector in global coordinates
- Range operator()(Domain const& x) const
+ public:
+ /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
+ /// reference to this DOFVector in the expression.
+ template <class Expr>
+ void interpolate_noalias(Expr&& expr)
{
- error_exit("Not implemented.");
- return Range(0);
- }
+ auto const& basis = this->basis();
+ auto const& treePath = this->treePath();
- /// \brief Create a local function for this view on the DOFVector. \relates LocalFunction
- friend LocalFunction localFunction(DOFVectorConstView const& self)
- {
- return LocalFunction{self};
- }
+ auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), basis.gridView());
- EntitySet const& entitySet() const
- {
- return entitySet_;
+ Dune::Functions::interpolate(basis, treePath, coefficients(), std::forward<decltype(gridFct)>(gridFct));
}
- public:
- /// Return global basis
- GlobalBasis const& basis() const
+ /// Interpolation of GridFunction to DOFVector
+ template <class Expr>
+ void interpolate(Expr&& expr)
{
- return dofVector_->basis();
- }
+ // create temporary copy of data
+ DOFVector<GlobalBasis,RangeType> tmp(coefficients());
+ Self tmpView{tmp, this->treePath()};
+ tmpView.interpolate_noalias(std::forward<Expr>(expr));
- /// Return treePath associated with this view
- TreePath const& treePath() const
- {
- return treePath_;
+ // move data from temporary vector into stored DOFVector
+ coefficients().vector() = std::move(tmp.vector());
}
- /// Return const coefficient vector
- DOFVector<GlobalBasis,RangeType> const& coefficients() const
+ /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
+ template <class Expr>
+ DOFVectorView& operator<<(Expr&& expr)
{
- return *dofVector_;
+ interpolate(expr);
+ return *this;
}
- 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)
- : Super(dofVector, treePath)
- , mutableDofVector_(&dofVector)
- {}
-
- public:
- /// Interpolation of GridFunction to DOFVector
template <class Expr>
- DOFVectorMutableView& interpolate(Expr&& expr)
+ DOFVectorView& operator+=(Expr&& expr)
{
- auto const& basis = Super::basis();
- auto const& treePath = Super::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));
-
- // move data from temporary vector into stored DOFVector
- mutableDofVector_->vector() = std::move(tmp.vector());
+ interpolate((*this) + expr);
return *this;
}
template <class Expr>
- DOFVectorMutableView& operator<<(Expr&& expr)
+ DOFVectorView& operator-=(Expr&& expr)
{
- return interpolate(expr);
+ interpolate((*this) - expr);
+ return *this;
}
@@ -307,42 +83,20 @@ namespace AMDiS
DOFVector<GlobalBasis,RangeType>* 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.
+ /// 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};
- }
-
-
- // A Generator for a const \ref DOFVectorView.
- template <class GlobalBasis, class RangeType>
- auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType> const& dofVector)
- {
- auto treePath = Dune::TypeTree::hybridTreePath();
- return DOFVectorConstView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
+ return DOFVectorView<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
}
- // A Generator for a mutable \ref DOFVectorView.
+ /// A Generator for a mutable \ref DOFVectorView
template <class GlobalBasis, class RangeType>
auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
- return DOFVectorMutableView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
+ return DOFVectorView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
}
-#endif
} // end namespace AMDiS
-
-#include "DOFVectorView.inc.hpp"
diff --git a/src/amdis/gridfunctions/DerivativeGridFunction.hpp b/src/amdis/gridfunctions/DerivativeGridFunction.hpp
index e4c0adc7..c9b3811d 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 00000000..396fe04c
--- /dev/null
+++ b/src/amdis/gridfunctions/DiscreteFunction.hpp
@@ -0,0 +1,245 @@
+#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
+ **/
+ template <class GlobalBasisType, class RangeType, class TreePathType>
+ class DiscreteFunction
+ {
+ 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;
+
+ 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;
+ };
+
+
+ public: // a local view on the values
+
+ /// A LocalFunction, i.e., an element local view on the DOFVector
+ class LocalFunction
+ {
+ public:
+ using Domain = typename DiscreteFunction::LocalDomain;
+ using Range = typename DiscreteFunction::Range;
+
+ enum { hasDerivative = true };
+
+ private:
+ using LocalView = typename GlobalBasis::LocalView;
+
+ public:
+ LocalFunction(DiscreteFunction const& globalFunction)
+ : globalFunction_(&globalFunction)
+ , localView_(globalFunction_->basis().localView())
+ , subTree_(&child(localView_.tree(), globalFunction_->treePath()))
+ {}
+
+ 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:
+ DiscreteFunction const* globalFunction_;
+ LocalView localView_;
+ SubTree const* subTree_;
+
+ bool bound_ = false;
+ };
+
+
+ public:
+ /// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
+ DiscreteFunction(DOFVector<GlobalBasis,RangeType> const& dofVector, TreePath const& treePath)
+ : dofVector_(&dofVector)
+ , treePath_(treePath)
+ , entitySet_(dofVector.basis().gridView())
+ , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(dofVector.basis(), treePath))
+ {}
+
+ /// 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};
+ }
+
+ 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 Generator for a \ref DiscreteFunction
+ template <class GlobalBasis, class RangeType, class TreePath>
+ auto makeDiscreteFunction(DOFVector<GlobalBasis, RangeType> const& dofVector, TreePath const& treePath)
+ {
+ return DiscreteFunction<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
+ }
+
+ /// A Generator for a \ref DiscreteFunction
+ template <class GlobalBasis, class RangeType>
+ auto makeDiscreteFunction(DOFVector<GlobalBasis, RangeType> const& dofVector)
+ {
+ auto treePath = Dune::TypeTree::hybridTreePath();
+ return DiscreteFunction<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
+ }
+
+} // end namespace AMDiS
+
+#include "DiscreteFunction.inc.hpp"
diff --git a/src/amdis/gridfunctions/DOFVectorView.inc.hpp b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
similarity index 83%
rename from src/amdis/gridfunctions/DOFVectorView.inc.hpp
rename to src/amdis/gridfunctions/DiscreteFunction.inc.hpp
index 939b89a0..fc6faa73 100644
--- a/src/amdis/gridfunctions/DOFVectorView.inc.hpp
+++ b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
@@ -3,10 +3,12 @@
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/LocalBasisCache.hpp>
+#include <dune/grid/utility/hierarchicsearch.hh>
+
namespace AMDiS {
template <class GB, class RT, class TP>
-typename DOFVectorConstView<GB,RT,TP>::Range DOFVectorConstView<GB,RT,TP>::
+typename DiscreteFunction<GB,RT,TP>::Range DiscreteFunction<GB,RT,TP>::
LocalFunction::operator()(LocalDomain const& x) const
{
assert( bound_ );
@@ -50,7 +52,7 @@ LocalFunction::operator()(LocalDomain const& x) const
template <class GB, class RT, class TP>
-typename DOFVectorConstView<GB,RT,TP>::DerivativeRange DOFVectorConstView<GB,RT,TP>::
+typename DiscreteFunction<GB,RT,TP>::DerivativeRange DiscreteFunction<GB,RT,TP>::
GradientLocalFunction::operator()(LocalDomain const& x) const
{
assert( bound_ );
@@ -106,4 +108,22 @@ GradientLocalFunction::operator()(LocalDomain const& x) const
return dy;
}
+
+template <class GB, class RT, class TP>
+typename DiscreteFunction<GB,RT,TP>::Range DiscreteFunction<GB,RT,TP>::
+operator()(Domain const& x) const
+{
+ using Grid = typename GB::GridView::Grid;
+ using IS = typename GB::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/test/CMakeLists.txt b/test/CMakeLists.txt
index a18847f2..363a35a6 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
CMD_ARGS "${CMAKE_SOURCE_DIR}/examples/init/ellipt.dat.2d")
diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp
new file mode 100644
index 00000000..bd3edbf4
--- /dev/null
+++ b/test/DiscreteFunctionTest.cpp
@@ -0,0 +1,102 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <iostream>
+
+#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();
+ }
+
+ AMDiS::finalize();
+ return 0;
+}
--
GitLab