Commit 276cfed9 authored by Praetorius, Simon's avatar Praetorius, Simon

removed nested LocalFunction and GradientLocalFunction from inline implementation

parent a4aeae1d
......@@ -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 DiscreteFunction = AMDiS::DiscreteFunction<GlobalBasis,RangeType,TreePath>;
using Range = typename DiscreteFunction::Range;
using Range = typename DiscreteFunction<GB,VT,TP>::Range;
/// Dimension of the mesh
static constexpr int dim = GridView::dimension;
......@@ -60,26 +63,25 @@ namespace AMDiS
public:
/// Constructor.
FileWriter(std::string const& baseName,
DiscreteFunction const& discreteFct)
: FileWriterInterface(baseName)
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());
......@@ -109,7 +111,7 @@ namespace AMDiS
}
private:
DiscreteFunction discreteFct_;
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,
DiscreteFunction<GlobalBasis,Range,TreePath> const& discreteFct)
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, discreteFct);
return std::make_shared<FileWriter<GlobalBasis,ValueType,TreePath>>(name, discreteFct);
}
} // end namespace AMDiS
......@@ -6,19 +6,24 @@
namespace AMDiS
{
/// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
template <class GlobalBasisType, class RangeType, class TreePathType>
template <class GB, class VT, class TP>
class DOFVectorView
: public DiscreteFunction<GlobalBasisType, RangeType, TreePathType>
: public DiscreteFunction<GB, VT, TP>
{
using Self = DOFVectorView;
using Super = DiscreteFunction<GlobalBasisType, RangeType, TreePathType>;
using Super = DiscreteFunction<GB, VT, TP>;
using GlobalBasis = GlobalBasisType;
using TreePath = TreePathType;
using GlobalBasis = GB;
using TreePath = TP;
public:
/// Constructor forwards to the treePath constructor, with empty TreePath
DOFVectorView(DOFVector<GB,VT>& dofVector)
: DOFVectorView{dofVector, Dune::TypeTree::hybridTreePath()}
{}
/// Constructor. Stores a pointer to the mutable `dofvector`.
DOFVectorView(DOFVector<GlobalBasis,RangeType>& dofVector, TreePath const& treePath)
DOFVectorView(DOFVector<GB,VT>& dofVector, TP const& treePath)
: Super(dofVector, treePath)
, mutableDofVector_(&dofVector)
{}
......@@ -26,6 +31,13 @@ namespace AMDiS
public:
/// \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>
void interpolate_noalias(Expr&& expr)
{
......@@ -34,15 +46,25 @@ namespace AMDiS
auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), basis.gridView());
Dune::Functions::interpolate(basis, treePath, coefficients(), std::forward<decltype(gridFct)>(gridFct));
Dune::Functions::interpolate(basis, treePath, coefficients(),
std::forward<decltype(gridFct)>(gridFct));
}
/// Interpolation of GridFunction to DOFVector
/// \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<GlobalBasis,RangeType> tmp(coefficients());
DOFVector<GB,VT> tmp(coefficients());
Self tmpView{tmp, this->treePath()};
tmpView.interpolate_noalias(std::forward<Expr>(expr));
......@@ -58,6 +80,7 @@ namespace AMDiS
return *this;
}
/// \brief interpolate `(*this) + expr` to DOFVector
template <class Expr>
DOFVectorView& operator+=(Expr&& expr)
{
......@@ -65,6 +88,7 @@ namespace AMDiS
return *this;
}
/// \brief interpolate `(*this) - expr` to DOFVector
template <class Expr>
DOFVectorView& operator-=(Expr&& expr)
{
......@@ -72,31 +96,37 @@ namespace AMDiS
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_;
};
#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 mutable \ref DOFVectorView
template <class GlobalBasis, class RangeType, class TreePath>
auto makeDOFVectorView(DOFVector<GlobalBasis, RangeType>& dofVector, TreePath const& treePath)
template <class GlobalBasis, class ValueType, class TreePath>
auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector, TreePath const& treePath)
{
return DOFVectorView<GlobalBasis, RangeType, 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)
template <class GlobalBasis, class ValueType>
auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
return DOFVectorView<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
return DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath};
}
} // end namespace AMDiS
......@@ -17,13 +17,19 @@ 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 GlobalBasisType, class RangeType, class TreePathType>
template <class GB, class VT, class TP>
class DiscreteFunction
{
public:
using GlobalBasis = GlobalBasisType;
using TreePath = TreePathType;
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>;
......@@ -32,159 +38,39 @@ namespace AMDiS
using GridView = typename GlobalBasis::GridView;
using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
public:
/// Global coordinates of the EntitySet
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;
/// 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 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:
/// 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<GlobalBasis,RangeType> const& dofVector, TreePath 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))
{}
/// Evaluate DiscreteFunction in global coordinates. NOTE: expensive
/// \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
......@@ -193,51 +79,59 @@ namespace AMDiS
return LocalFunction{self};
}
/// \brief Return a \ref Dune::Functions::GridViewEntitySet
EntitySet const& entitySet() const
{
return entitySet_;
}
public:
/// Return global basis
/// \brief Return global basis bound to the DOFVector
GlobalBasis const& basis() const
{
return dofVector_->basis();
}
/// Return treePath associated with this view
/// \brief Return treePath associated with this discrete function
TreePath const& treePath() const
{
return treePath_;
}
/// Return const coefficient vector
DOFVector<GlobalBasis,RangeType> const& coefficients() const
/// \brief Return const coefficient vector
DOFVector<GB,VT> const& coefficients() const
{
return *dofVector_;
}
protected:
DOFVector<GlobalBasis,RangeType> const* dofVector_;
TreePath const treePath_;
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 RangeType, class TreePath>
auto makeDiscreteFunction(DOFVector<GlobalBasis, RangeType> const& dofVector, TreePath const& treePath)
template <class GlobalBasis, class ValueType, class TreePath>
auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector, TreePath const& treePath)
{
return DiscreteFunction<GlobalBasis, RangeType, TreePath>{dofVector, treePath};
return DiscreteFunction<GlobalBasis, ValueType, TreePath>{dofVector, treePath};
}
/// A Generator for a \ref DiscreteFunction
template <class GlobalBasis, class RangeType>
auto makeDiscreteFunction(DOFVector<GlobalBasis, RangeType> const& dofVector)
template <class GlobalBasis, class ValueType>
auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
return DiscreteFunction<GlobalBasis, RangeType, decltype(treePath)>{dofVector, treePath};
return DiscreteFunction<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath};
}
} // end namespace AMDiS
......
......@@ -7,9 +7,140 @@
namespace AMDiS {
template <class GB, class RT, class TP>
typename DiscreteFunction<GB,RT,TP>::Range DiscreteFunction<GB,RT,TP>::
LocalFunction::operator()(LocalDomain const& x) const
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);
......@@ -51,9 +182,9 @@ LocalFunction::operator()(LocalDomain const& x) const
}
template <class GB, class RT, class TP>
typename DiscreteFunction<GB,RT,TP>::DerivativeRange DiscreteFunction<GB,RT,TP>::
GradientLocalFunction::operator()(LocalDomain const& x) const
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;
......@@ -109,12 +240,12 @@ GradientLocalFunction::operator()(LocalDomain const& x) const
}
template <class GB, class RT, class TP>
typename DiscreteFunction<GB,RT,TP>::Range DiscreteFunction<GB,RT,TP>::
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 GB::GridView::Grid;
using IS = typename GB::GridView::IndexSet;
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()};
......
......@@ -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