Skip to content
Snippets Groups Projects
Commit 2b9a6337 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

DiscreteFunction vs. DOFVectorView

parent bd634196
No related branches found
No related tags found
No related merge requests found
......@@ -10,7 +10,7 @@
#include <amdis/Observer.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/typetree/TreePath.hpp>
namespace Dune
......@@ -88,7 +88,7 @@ namespace AMDiS
auto child(TreePath const& path = {})
{
auto&& tp = makeTreePath(path);
return makeDOFVectorView(*this, tp);
return makeDiscreteFunction(*this, tp);
}
template <class TreePath = RootTreePath>
......@@ -101,7 +101,7 @@ namespace AMDiS
/// Interpolation of GridFunction to DOFVector, assuming that there is no
/// reference to this DOFVector in the expression.
/// See \ref DOFVectorView::interpolate_noalias
/// See \ref DiscreteFunction::interpolate_noalias
template <class Expr, class Tag = tag::average>
void interpolate_noalias(Expr&& expr, Tag strategy)
{
......@@ -109,7 +109,7 @@ namespace AMDiS
}
/// Interpolation of GridFunction to DOFVector.
/// See \ref DOFVectorView::interpolate
/// See \ref DiscreteFunction::interpolate
template <class Expr, class Tag = tag::average>
void interpolate(Expr&& expr, Tag strategy)
{
......@@ -117,7 +117,7 @@ namespace AMDiS
}
/// Interpolation of GridFunction to DOFVector.
/// See \ref DOFVectorView::interpolate
/// See \ref DiscreteFunction::interpolate
template <class Expr>
DOFVector& operator<<(Expr&& expr)
{
......
......@@ -40,7 +40,6 @@
#include <amdis/GridFunctions.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/gridfunctions/DOFVectorView.hpp>
#include <amdis/io/FileWriterBase.hpp>
......
......@@ -8,7 +8,7 @@ install(FILES
DerivativeGridFunction.hpp
DiscreteFunction.hpp
DiscreteFunction.inc.hpp
DOFVectorView.hpp
DiscreteLocalFunction.inc.hpp
FunctorGridFunction.hpp
GridFunction.hpp
OperationsGridFunction.hpp
......
#pragma once
#include <amdis/functions/Interpolate.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
namespace AMDiS
{
/// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
template <class GB, class VT, class TP>
class DOFVectorView
: public DiscreteFunction<GB, VT, TP>
{
using Self = DOFVectorView;
using Super = DiscreteFunction<GB, VT, TP>;
using GlobalBasis = GB;
using TreePath = TP;
public:
/// Constructor. Stores a pointer to the mutable `dofvector`.
template <class PreTreePath>
DOFVectorView(DOFVector<GB,VT>& dofVector, PreTreePath const& preTreePath)
: Super(dofVector, makeTreePath(preTreePath))
, mutableDofVector_(&dofVector)
{}
/// Constructor forwards to the treePath constructor, with empty TreePath
DOFVectorView(DOFVector<GB,VT>& dofVector)
: DOFVectorView(dofVector, Dune::TypeTree::hybridTreePath())
{}
public:
/// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
/// reference to this DOFVector in the expression.
/**
* **Example:**
* ```
* auto v = makeDOFVectorView(prob.solutionVector(),0);
* v.interpolate_noalias([](auto const& x) { return x[0]; });
* ```
**/
template <class Expr, class Tag = tag::average>
void interpolate_noalias(Expr&& expr, Tag strategy = {})
{
auto const& basis = *this->basis();
auto const& treePath = this->treePath();
auto&& gf = makeGridFunction(FWD(expr), basis.gridView());
if (std::is_same<Tag, tag::average>::value) {
auto counter = coefficients();
AMDiS::interpolate(basis, coefficients(), gf, treePath, counter);
coefficients().forEach([&counter](std::size_t dof, auto& coeff)
{
coeff /= std::max(double(counter.at(dof)), 1.0);
});
} else {
AMDiS::interpolate(basis, coefficients(), gf, treePath);
}
}
/// \brief Interpolation of GridFunction to DOFVector
/**
* **Example:**
* ```
* auto v = makeDOFVectorView(prob.solutionVector(),0);
* v.interpolate(v + [](auto const& x) { return x[0]; });
* ```
* Allows to have a reference to the DOFVector in the expression, e.g. as
* \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction.
**/
template <class Expr, class Tag = tag::average>
void interpolate(Expr&& expr, Tag strategy = {})
{
// create temporary copy of data
DOFVector<GB,VT> tmp(coefficients());
Self tmpView{tmp, this->treePath()};
tmpView.interpolate_noalias(FWD(expr), strategy);
// move data from temporary vector into stored DOFVector
coefficients().backend() = std::move(tmp.backend());
}
/// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
template <class Expr>
DOFVectorView& operator<<(Expr&& expr)
{
interpolate(FWD(expr));
return *this;
}
/// \brief interpolate `(*this) + expr` to DOFVector
template <class Expr>
DOFVectorView& operator+=(Expr&& expr)
{
interpolate((*this) + expr);
return *this;
}
/// \brief interpolate `(*this) - expr` to DOFVector
template <class Expr>
DOFVectorView& operator-=(Expr&& expr)
{
interpolate((*this) - expr);
return *this;
}
/// Return the mutable DOFVector
DOFVector<GB,VT>& coefficients() { return *mutableDofVector_; }
/// Return the const DOFVector
using Super::coefficients;
protected:
DOFVector<GB,VT>* mutableDofVector_;
};
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
// Deduction guide for DOFVectorView class
template <class GlobalBasis, class ValueType, class PreTreePath>
DOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector, PreTreePath const& preTreePath)
-> DOFVectorView<GlobalBasis, ValueType, TYPEOF(makeTreePath(preTreePath))>;
// Deduction guide for DOFVectorView class
template <class GlobalBasis, class ValueType>
DOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector)
-> DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>;
#endif
/// A Generator for a mutable \ref DOFVectorView
template <class GlobalBasis, class ValueType, class PreTreePath>
auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector, PreTreePath const& preTreePath)
{
auto treePath = makeTreePath(preTreePath);
return DOFVectorView<GlobalBasis, ValueType, decltype(treePath)>{dofVector, treePath};
}
/// A Generator for a mutable \ref DOFVectorView
template <class GlobalBasis, class ValueType>
auto makeDOFVectorView(DOFVector<GlobalBasis, ValueType>& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
return DOFVectorView<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath};
}
} // end namespace AMDiS
......@@ -16,23 +16,63 @@
namespace AMDiS
{
template <class GB, class VT>
class DOFVector;
/// \class DiscreteFunction
/// \brief A view on a subspace of a \ref DOFVector
/**
* \ingroup GridFunctions
*
* \tparam GB Type of the global basis
* \tparam VT Coefficient type of the DOFVector
* \tparam TP A realization of \ref Dune::TypeTree::HybridTreePath
*
* **Requirements:**
* - GB models \ref Dune::Functions::Concept::GlobalBasis
**/
template <class GB, class VT, class TP>
class DiscreteFunction
* \ingroup GridFunctions
*
* \tparam GB Type of the global basis
* \tparam VT Coefficient type of the DOFVector
* \tparam TP A realization of \ref Dune::TypeTree::HybridTreePath
* \tparam is_const Specifies whether a const or mutable view is implemented.
*
* **Requirements:**
* - GB models \ref Dune::Functions::Concept::GlobalBasis
**/
template <class GB, class VT, class TP, bool is_const>
class DiscreteFunction;
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
// Deduction guide for DiscreteFunction class
template <class GB, class VT>
DiscreteFunction(DOFVector<GB, VT> const& dofVector)
-> DiscreteFunction<GB, VT, Dune::TypeTree::HybridTreePath<>, true>;
template <class GB, class VT>
DiscreteFunction(DOFVector<GB, VT>& dofVector)
-> DiscreteFunction<GB, VT, Dune::TypeTree::HybridTreePath<>, false>;
#endif
/// A Generator for a const \ref DiscreteFunction
template <class GlobalBasis, class ValueType,
class PreTreePath = Dune::TypeTree::HybridTreePath<>>
auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector,
PreTreePath const& preTreePath = {})
{
// static_assert(Dune::IsNumber<VT>::value, "");
auto treePath = makeTreePath(preTreePath);
return DiscreteFunction<GlobalBasis, ValueType, decltype(treePath), true>{dofVector, treePath};
}
/// A Generator for a mutable \ref DiscreteFunction
template <class GlobalBasis, class ValueType,
class PreTreePath = Dune::TypeTree::HybridTreePath<>>
auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType>& dofVector,
PreTreePath const& preTreePath = {})
{
auto treePath = makeTreePath(preTreePath);
return DiscreteFunction<GlobalBasis, ValueType, decltype(treePath), false>{dofVector, treePath};
}
/// A Const DiscreteFunction
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP,true>
{
private:
using GlobalBasis = GB;
using TreePath = TP;
......@@ -71,18 +111,13 @@ namespace AMDiS
public:
/// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
DiscreteFunction(DOFVector<GB,VT> const& dofVector, TP const& treePath)
DiscreteFunction(DOFVector<GB,VT> const& dofVector, TP const& treePath = {})
: dofVector_(&dofVector)
, treePath_(treePath)
, entitySet_(dofVector.basis()->gridView())
, nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(*dofVector.basis(), treePath))
{}
/// Constructor forwards to the treePath constructor, with empty TreePath
DiscreteFunction(DOFVector<GB,VT> const& dofVector)
: DiscreteFunction(dofVector, Dune::TypeTree::hybridTreePath())
{}
/// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive
Range operator()(Domain const& x) const;
......@@ -98,7 +133,6 @@ namespace AMDiS
return entitySet_;
}
public:
/// \brief Return global basis bound to the DOFVector
std::shared_ptr<GlobalBasis const> basis() const
{
......@@ -125,29 +159,86 @@ namespace AMDiS
};
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
// Deduction guide for DiscreteFunction class
template <class GlobalBasis, class ValueType>
DiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector)
-> DiscreteFunction<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>;
#endif
/// A Generator for a \ref DiscreteFunction
template <class GlobalBasis, class ValueType, class PreTreePath>
auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector, PreTreePath const& preTreePath)
/// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP,false>
: public DiscreteFunction<GB, VT, TP, true>
{
auto treePath = makeTreePath(preTreePath);
return DiscreteFunction<GlobalBasis, ValueType, decltype(treePath)>{dofVector, treePath};
}
using Self = DiscreteFunction<GB, VT, TP, false>;
using Super = DiscreteFunction<GB, VT, TP, true>;
/// A Generator for a \ref DiscreteFunction
template <class GlobalBasis, class ValueType>
auto makeDiscreteFunction(DOFVector<GlobalBasis, ValueType> const& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
return DiscreteFunction<GlobalBasis, ValueType, Dune::TypeTree::HybridTreePath<>>{dofVector, treePath};
}
using GlobalBasis = GB;
using TreePath = TP;
public:
/// Constructor. Stores a pointer to the mutable `dofvector`.
DiscreteFunction(DOFVector<GB,VT>& dofVector, TP const& treePath = {})
: Super(dofVector, treePath)
, mutableDofVector_(&dofVector)
{}
public:
/// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
/// reference to this DOFVector in the expression.
/**
* **Example:**
* ```
* auto v = makeDiscreteFunction(prob.solutionVector(),0);
* v.interpolate_noalias([](auto const& x) { return x[0]; });
* ```
**/
template <class Expr, class Tag = tag::average>
void interpolate_noalias(Expr&& expr, Tag strategy = {});
/// \brief Interpolation of GridFunction to DOFVector
/**
* **Example:**
* ```
* auto v = makeDiscreteFunction(prob.solutionVector(),0);
* v.interpolate(v + [](auto const& x) { return x[0]; });
* ```
* Allows to have a reference to the DOFVector in the expression, e.g. as
* \ref DiscreteFunction or \ref gradientAtQP() of a DiscreteFunction.
**/
template <class Expr, class Tag = tag::average>
void interpolate(Expr&& expr, Tag strategy = {});
/// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
template <class Expr>
Self& operator<<(Expr&& expr)
{
interpolate(FWD(expr));
return *this;
}
/// \brief interpolate `(*this) + expr` to DOFVector
template <class Expr>
Self& operator+=(Expr&& expr)
{
interpolate((*this) + expr);
return *this;
}
/// \brief interpolate `(*this) - expr` to DOFVector
template <class Expr>
Self& operator-=(Expr&& expr)
{
interpolate((*this) - expr);
return *this;
}
/// Return the mutable DOFVector
DOFVector<GB,VT>& coefficients() { return *mutableDofVector_; }
/// Return the const DOFVector
using Super::coefficients;
protected:
DOFVector<GB,VT>* mutableDofVector_;
};
} // end namespace AMDiS
#include "DiscreteLocalFunction.inc.hpp"
#include "DiscreteFunction.inc.hpp"
#pragma once
#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/utility/LocalBasisCache.hpp>
#include <amdis/utility/LocalToGlobalAdapter.hpp>
#include <type_traits>
#include <dune/common/ftraits.hh>
#include <dune/grid/utility/hierarchicsearch.hh>
namespace AMDiS {
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP>::LocalFunction
{
public:
using Domain = typename EntitySet::LocalCoordinate;
using Range = typename DiscreteFunction::Range;
enum { hasDerivative = true };
private:
using LocalView = typename GlobalBasis::LocalView;
using Element = typename EntitySet::Element;
using Geometry = typename Element::Geometry;
public:
/// Constructor. Stores a copy of the DiscreteFunction.
LocalFunction(DiscreteFunction const& globalFunction)
: globalFunction_(globalFunction)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
/// Copy constructor.
LocalFunction(LocalFunction const& other)
: globalFunction_(other.globalFunction_)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
/// \brief Bind the LocalView to the element
void bind(Element const& element)
{
localView_.bind(element);
globalFunction_.coefficients().gather(localView_, localCoefficients_);
bound_ = true;
}
/// \brief Unbind the LocalView from the element
void unbind()
{
localView_.unbind();
bound_ = false;
}
/// \brief Evaluate LocalFunction at bound element in local coordinates
Range operator()(Domain const& x) const;
/// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
GradientLocalFunction makeDerivative(tag::gradient type) const
{
return GradientLocalFunction{globalFunction_, type};
}
DivergenceLocalFunction makeDerivative(tag::divergence type) const
{
return DivergenceLocalFunction{globalFunction_, type};
}
PartialLocalFunction makeDerivative(tag::partial type) const
{
return PartialLocalFunction{globalFunction_, type};
}
/// \brief The \ref polynomialDegree() of the LocalFunctions
int order() const
{
assert( bound_ );
return polynomialDegree(*subTree_);
}
/// \brief Return the bound element
Element const& localContext() const
{
assert( bound_ );
return localView_.element();
}
private:
DiscreteFunction globalFunction_;
LocalView localView_;
SubTree const* subTree_;
std::vector<VT> localCoefficients_;
bool bound_ = false;
};
#include <amdis/functions/Interpolate.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
namespace AMDiS {
// Evaluate DiscreteFunction in global coordinates
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
LocalFunction::operator()(Domain const& x) const
typename DiscreteFunction<GB,VT,TP,true>::Range DiscreteFunction<GB,VT,TP,true>::
operator()(Domain const& x) const
{
assert( bound_ );
Range y(0);
auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
{
auto localBasisCache = makeNodeCache(node);
auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(), x);
std::size_t size = node.finiteElement().size();
// Get range entry associated to this node
auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
for (std::size_t i = 0; i < size; ++i) {
// Get coefficient associated to i-th shape function
auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
// Get value of i-th shape function
auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
using Grid = typename GlobalBasis::GridView::Grid;
using IS = typename GlobalBasis::GridView::IndexSet;
std::size_t dimC = c.size();
std::size_t dimV = v.size();
assert(dimC*dimV == std::size_t(re.size()));
for(std::size_t j = 0; j < dimC; ++j) {
auto&& c_j = c[j];
for(std::size_t k = 0; k < dimV; ++k)
re[j*dimV + k] += c_j*v[k];
}
}
});
auto const& gv = this->basis()->gridView();
Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
return y;
auto element = hsearch.findEntity(x);
auto geometry = element.geometry();
auto localFct = localFunction(*this);
localFct.bind(element);
return localFct(geometry.local(x));
}
// Interpolation of GridFunction to DOFVector
template <class GB, class VT, class TP>
template <class Type>
class DiscreteFunction<GB,VT,TP>::DerivativeLocalFunctionBase
{
using R = typename DiscreteFunction::Range;
using D = typename DiscreteFunction::Domain;
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
using Traits = DerivativeTraits<RawSignature, Type>;
public:
using Domain = typename EntitySet::LocalCoordinate;
using Range = typename Traits::Range;
enum { hasDerivative = false };
private:
using LocalView = typename GlobalBasis::LocalView;
using Element = typename EntitySet::Element;
using Geometry = typename Element::Geometry;
public:
/// Constructor. Stores a copy of the DiscreteFunction.
DerivativeLocalFunctionBase(DiscreteFunction const& globalFunction, Type const& type)
: globalFunction_(globalFunction)
, type_(type)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
/// Copy constructor.
DerivativeLocalFunctionBase(DerivativeLocalFunctionBase const& other)
: globalFunction_(other.globalFunction_)
, type_(other.type_)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
void bind(Element const& element)
{
localView_.bind(element);
geometry_.emplace(element.geometry());
globalFunction_.coefficients().gather(localView_, localCoefficients_);
bound_ = true;
}
void unbind()
{
localView_.unbind();
geometry_.reset();
bound_ = false;
}
int order() const
{
assert( bound_ );
return std::max(0, polynomialDegree(*subTree_)-1);
}
/// Return the bound element
Element const& localContext() const
{
assert( bound_ );
return localView_.element();
}
Geometry const& geometry() const
{
assert( bound_ );
return geometry_.value();
}
LocalView const& localView() const
{
return localView_;
}
protected:
DiscreteFunction globalFunction_;
Type type_;
LocalView localView_;
SubTree const* subTree_;
Dune::Std::optional<Geometry> geometry_;
std::vector<VT> localCoefficients_;
bool bound_ = false;
};
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP>::GradientLocalFunction
: public DerivativeLocalFunctionBase<tag::gradient>
{
using Super = DerivativeLocalFunctionBase<tag::gradient>;
public:
using Range = typename Super::Range;
using Domain = typename Super::Domain;
// adopt constructor from base class
using DerivativeLocalFunctionBase<tag::gradient>::DerivativeLocalFunctionBase;
/// Evaluate Gradient at bound element in local coordinates
Range operator()(Domain const& x) const
{
assert( this->bound_ );
Range dy(0);
auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_;
for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp)
{
auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry());
auto const& gradients = localBasis.gradientsAt(x);
// 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) {
// Get coefficient associated to i-th shape function
auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]);
// 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;
}
using Super::localCoefficients_;
};
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP>::DivergenceLocalFunction
: public DerivativeLocalFunctionBase<tag::divergence>
{
using Super = DerivativeLocalFunctionBase<tag::divergence>;
public:
using Range = typename Super::Range;
using Domain = typename Super::Domain;
// adopt constructor from base class
using DerivativeLocalFunctionBase<tag::divergence>::DerivativeLocalFunctionBase;
/// Evaluate divergence at bound element in local coordinates
Range operator()(Domain const& x) const
{
return evaluate_node(x, bool_t<SubTree::isPower && SubTree::degree() == GridView::dimensionworld>{});
}
private:
Range evaluate_node(Domain const& x, std::false_type) const
{
error_exit("Cannot calculate divergence(node) where node is not a power node.");
return Range(0);
}
Range evaluate_node(Domain const& x, std::true_type) const
{
static_assert(static_size_v<Range> == 1, "Implemented for scalar coefficients only.");
assert( this->bound_ );
Range dy(0);
auto&& node = *this->subTree_;
auto localBasis = makeLocalToGlobalBasisAdapter(node.child(0), this->geometry());
auto const& gradients = localBasis.gradientsAt(x);
auto re = Dune::Functions::flatVectorView(dy);
assert(int(re.size()) == 1);
for (std::size_t i = 0; i < localBasis.size(); ++i) {
auto grad = Dune::Functions::flatVectorView(gradients[i]);
assert(int(grad.size()) == GridView::dimensionworld);
for (std::size_t j = 0; j < GridView::dimensionworld; ++j)
re[0] += localCoefficients_[node.child(j).localIndex(i)] * grad[j];
}
return dy;
}
using Super::localCoefficients_;
};
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP>::PartialLocalFunction
: public DerivativeLocalFunctionBase<tag::partial>
template <class Expr, class Tag>
void DiscreteFunction<GB,VT,TP,false>::
interpolate_noalias(Expr&& expr, Tag strategy)
{
using Super = DerivativeLocalFunctionBase<tag::partial>;
public:
using Range = typename Super::Range;
using Domain = typename Super::Domain;
auto const& basis = *this->basis();
auto const& treePath = this->treePath();
using DerivativeLocalFunctionBase<tag::partial>::DerivativeLocalFunctionBase;
auto&& gf = makeGridFunction(FWD(expr), basis.gridView());
/// Evaluate partial derivative at bound element in local coordinates
Range operator()(Domain const& x) const
{
assert( this->bound_ );
Range dy(0);
if (std::is_same<Tag, tag::average>::value) {
auto counter = coefficients();
AMDiS::interpolate(basis, coefficients(), gf, treePath, counter);
std::size_t comp = this->type_.comp;
auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_;
for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp)
coefficients().forEach([&counter](std::size_t dof, auto& coeff)
{
auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry());
auto const& partial = localBasis.partialsAt(x, comp);
// 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) {
// Get coefficient associated to i-th shape function
auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]);
// Get value of i-th transformed reference partial_derivative
auto d_comp = Dune::Functions::flatVectorView(partial[i]);
std::size_t dimC = c.size();
std::size_t dimV = d_comp.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*d_comp[k];
}
}
coeff /= std::max(double(counter.at(dof)), 1.0);
});
return dy;
} else {
AMDiS::interpolate(basis, coefficients(), gf, treePath);
}
using Super::localCoefficients_;
};
}
// Interpolation of GridFunction to DOFVector
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
operator()(Domain const& x) const
template <class Expr, class Tag>
void DiscreteFunction<GB,VT,TP,false>::
interpolate(Expr&& expr, Tag strategy)
{
using Grid = typename GlobalBasis::GridView::Grid;
using IS = typename GlobalBasis::GridView::IndexSet;
// create temporary copy of data
DOFVector<GB,VT> tmp(coefficients());
auto const& gv = this->basis()->gridView();
Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
Self tmpView{tmp, this->treePath()};
tmpView.interpolate_noalias(FWD(expr), strategy);
auto element = hsearch.findEntity(x);
auto geometry = element.geometry();
auto localFct = localFunction(*this);
localFct.bind(element);
return localFct(geometry.local(x));
// move data from temporary vector into stored DOFVector
coefficients().backend() = std::move(tmp.backend());
}
} // end namespace AMDiS
#pragma once
#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/utility/LocalBasisCache.hpp>
#include <amdis/utility/LocalToGlobalAdapter.hpp>
#include <dune/common/ftraits.hh>
namespace AMDiS {
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP,true>::LocalFunction
{
public:
using Domain = typename EntitySet::LocalCoordinate;
using Range = typename DiscreteFunction::Range;
enum { hasDerivative = true };
private:
using LocalView = typename GlobalBasis::LocalView;
using Element = typename EntitySet::Element;
using Geometry = typename Element::Geometry;
public:
/// Constructor. Stores a copy of the DiscreteFunction.
LocalFunction(DiscreteFunction const& globalFunction)
: globalFunction_(globalFunction)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
/// Copy constructor.
LocalFunction(LocalFunction const& other)
: globalFunction_(other.globalFunction_)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
/// \brief Bind the LocalView to the element
void bind(Element const& element)
{
localView_.bind(element);
globalFunction_.coefficients().gather(localView_, localCoefficients_);
bound_ = true;
}
/// \brief Unbind the LocalView from the element
void unbind()
{
localView_.unbind();
bound_ = false;
}
/// \brief Evaluate LocalFunction at bound element in local coordinates
Range operator()(Domain const& x) const
{
assert( bound_ );
Range y(0);
auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
{
auto localBasisCache = makeNodeCache(node);
auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(), x);
std::size_t size = node.finiteElement().size();
// Get range entry associated to this node
auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
for (std::size_t i = 0; i < size; ++i) {
// Get coefficient associated to i-th shape function
auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
// Get value of i-th shape function
auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
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;
}
/// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
GradientLocalFunction makeDerivative(tag::gradient type) const
{
return GradientLocalFunction{globalFunction_, type};
}
DivergenceLocalFunction makeDerivative(tag::divergence type) const
{
return DivergenceLocalFunction{globalFunction_, type};
}
PartialLocalFunction makeDerivative(tag::partial type) const
{
return PartialLocalFunction{globalFunction_, type};
}
/// \brief The \ref polynomialDegree() of the LocalFunctions
int order() const
{
assert( bound_ );
return polynomialDegree(*subTree_);
}
/// \brief Return the bound element
Element const& localContext() const
{
assert( bound_ );
return localView_.element();
}
private:
DiscreteFunction globalFunction_;
LocalView localView_;
SubTree const* subTree_;
std::vector<VT> localCoefficients_;
bool bound_ = false;
};
template <class GB, class VT, class TP>
template <class Type>
class DiscreteFunction<GB,VT,TP,true>::DerivativeLocalFunctionBase
{
using R = typename DiscreteFunction::Range;
using D = typename DiscreteFunction::Domain;
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
using Traits = DerivativeTraits<RawSignature, Type>;
public:
using Domain = typename EntitySet::LocalCoordinate;
using Range = typename Traits::Range;
enum { hasDerivative = false };
private:
using LocalView = typename GlobalBasis::LocalView;
using Element = typename EntitySet::Element;
using Geometry = typename Element::Geometry;
public:
/// Constructor. Stores a copy of the DiscreteFunction.
DerivativeLocalFunctionBase(DiscreteFunction const& globalFunction, Type const& type)
: globalFunction_(globalFunction)
, type_(type)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
/// Copy constructor.
DerivativeLocalFunctionBase(DerivativeLocalFunctionBase const& other)
: globalFunction_(other.globalFunction_)
, type_(other.type_)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
void bind(Element const& element)
{
localView_.bind(element);
geometry_.emplace(element.geometry());
globalFunction_.coefficients().gather(localView_, localCoefficients_);
bound_ = true;
}
void unbind()
{
localView_.unbind();
geometry_.reset();
bound_ = false;
}
int order() const
{
assert( bound_ );
return std::max(0, polynomialDegree(*subTree_)-1);
}
/// Return the bound element
Element const& localContext() const
{
assert( bound_ );
return localView_.element();
}
Geometry const& geometry() const
{
assert( bound_ );
return geometry_.value();
}
LocalView const& localView() const
{
return localView_;
}
protected:
DiscreteFunction globalFunction_;
Type type_;
LocalView localView_;
SubTree const* subTree_;
Dune::Std::optional<Geometry> geometry_;
std::vector<VT> localCoefficients_;
bool bound_ = false;
};
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP,true>::GradientLocalFunction
: public DerivativeLocalFunctionBase<tag::gradient>
{
using Super = DerivativeLocalFunctionBase<tag::gradient>;
public:
using Range = typename Super::Range;
using Domain = typename Super::Domain;
// adopt constructor from base class
using DerivativeLocalFunctionBase<tag::gradient>::DerivativeLocalFunctionBase;
/// Evaluate Gradient at bound element in local coordinates
Range operator()(Domain const& x) const
{
assert( this->bound_ );
Range dy(0);
auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_;
for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp)
{
auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry());
auto const& gradients = localBasis.gradientsAt(x);
// 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) {
// Get coefficient associated to i-th shape function
auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]);
// 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;
}
using Super::localCoefficients_;
};
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP,true>::DivergenceLocalFunction
: public DerivativeLocalFunctionBase<tag::divergence>
{
using Super = DerivativeLocalFunctionBase<tag::divergence>;
public:
using Range = typename Super::Range;
using Domain = typename Super::Domain;
// adopt constructor from base class
using DerivativeLocalFunctionBase<tag::divergence>::DerivativeLocalFunctionBase;
/// Evaluate divergence at bound element in local coordinates
Range operator()(Domain const& x) const
{
return evaluate_node(x, bool_t<SubTree::isPower && SubTree::degree() == GridView::dimensionworld>{});
}
private:
Range evaluate_node(Domain const& x, std::false_type) const
{
error_exit("Cannot calculate divergence(node) where node is not a power node.");
return Range(0);
}
Range evaluate_node(Domain const& x, std::true_type) const
{
static_assert(static_size_v<Range> == 1, "Implemented for scalar coefficients only.");
assert( this->bound_ );
Range dy(0);
auto&& node = *this->subTree_;
auto localBasis = makeLocalToGlobalBasisAdapter(node.child(0), this->geometry());
auto const& gradients = localBasis.gradientsAt(x);
auto re = Dune::Functions::flatVectorView(dy);
assert(int(re.size()) == 1);
for (std::size_t i = 0; i < localBasis.size(); ++i) {
auto grad = Dune::Functions::flatVectorView(gradients[i]);
assert(int(grad.size()) == GridView::dimensionworld);
for (std::size_t j = 0; j < GridView::dimensionworld; ++j)
re[0] += localCoefficients_[node.child(j).localIndex(i)] * grad[j];
}
return dy;
}
using Super::localCoefficients_;
};
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP,true>::PartialLocalFunction
: public DerivativeLocalFunctionBase<tag::partial>
{
using Super = DerivativeLocalFunctionBase<tag::partial>;
public:
using Range = typename Super::Range;
using Domain = typename Super::Domain;
using DerivativeLocalFunctionBase<tag::partial>::DerivativeLocalFunctionBase;
/// Evaluate partial derivative at bound element in local coordinates
Range operator()(Domain const& x) const
{
assert( this->bound_ );
Range dy(0);
std::size_t comp = this->type_.comp;
auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_;
for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp)
{
auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry());
auto const& partial = localBasis.partialsAt(x, comp);
// 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) {
// Get coefficient associated to i-th shape function
auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]);
// Get value of i-th transformed reference partial_derivative
auto d_comp = Dune::Functions::flatVectorView(partial[i]);
std::size_t dimC = c.size();
std::size_t dimV = d_comp.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*d_comp[k];
}
}
});
return dy;
}
using Super::localCoefficients_;
};
} // end namespace AMDiS
......@@ -53,7 +53,7 @@ namespace AMDiS
private:
template <class GB, class VT, class TP, class ValCat>
std::unique_ptr<FileWriterInterface>
create_impl(std::string type, std::string prefix, DiscreteFunction<GB,VT,TP> const& data, ValCat) const
create_impl(std::string type, std::string prefix, DiscreteFunction<GB,VT,TP,true> const& data, ValCat) const
{
// ParaView VTK format, writer from dune-grid
if (type == "vtk")
......@@ -92,7 +92,7 @@ namespace AMDiS
// The value-category is unknown, like a composite/hierarchic vector or any unknown type.
template <class GB, class VT, class TP>
std::unique_ptr<FileWriterInterface>
create_impl(std::string type, std::string prefix, DiscreteFunction<GB,VT,TP> const& /*data*/, tag::unknown) const
create_impl(std::string type, std::string prefix, DiscreteFunction<GB,VT,TP,true> const& /*data*/, tag::unknown) const
{
// Backup writer, writing the grid and the solution vector
if (type == "backup")
......
......@@ -10,6 +10,7 @@
#include <amdis/common/Filesystem.hpp>
#include <amdis/common/StaticSize.hpp>
#include <amdis/common/ValueCategory.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/io/FileWriterBase.hpp>
#include <amdis/io/VTKSequenceWriter.hpp>
......@@ -42,7 +43,9 @@ namespace AMDiS
using GridView = typename GB::GridView;
using Writer = Dune::VTKWriter<GridView>;
using SeqWriter = VTKSequenceWriter<GridView>;
using Range = typename DiscreteFunction<GB,VT,TP>::Range;
using Function = DiscreteFunction<GB,VT,TP,true>;
using Range = typename Function::Range;
template <class R>
static constexpr Dune::VTK::FieldInfo::Type VTKFieldType() {
......@@ -56,7 +59,7 @@ namespace AMDiS
public:
/// Constructor.
VTKWriter(std::string const& name, DiscreteFunction<GB,VT,TP> const& discreteFct)
VTKWriter(std::string const& name, Function const& discreteFct)
: FileWriterBase(name)
, discreteFct_(discreteFct)
{
......@@ -103,7 +106,7 @@ namespace AMDiS
}
private:
DiscreteFunction<GB,VT,TP> discreteFct_;
Function discreteFct_;
std::shared_ptr<Writer> vtkWriter_;
std::shared_ptr<SeqWriter> vtkSeqWriter_;
......
......@@ -7,7 +7,6 @@
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/gridfunctions/DOFVectorView.hpp>
#include <amdis/typetree/TreePath.hpp>
#include "Tests.hpp"
......@@ -55,9 +54,20 @@ int main(int argc, char** argv)
auto U1 = *prob.solutionVector();
auto U2 = *prob.solutionVector();
auto u0 = makeDOFVectorView(U0);
auto u1 = makeDOFVectorView(U1);
auto u2 = makeDOFVectorView(U2);
auto u0 = makeDiscreteFunction(U0);
auto u1 = makeDiscreteFunction(U1);
auto u2 = makeDiscreteFunction(U2);
// Test const and mutable construction
auto const& U_c = *prob.solutionVector();
auto& U_m = *prob.solutionVector();
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
DiscreteFunction u0_c(U_c);
DiscreteFunction u0_m(U_m);
#endif
auto u1_c = makeDiscreteFunction(U_c);
auto u1_m = makeDiscreteFunction(U_m);
using Range = typename decltype(u0)::Range;
......@@ -123,10 +133,10 @@ int main(int argc, char** argv)
}
auto V0 = makeDOFVector<double>(prob.globalBasis());
auto v0 = makeDOFVectorView(V0);
auto v0 = makeDiscreteFunction(V0);
v0 << expr1;
// test makeDiscreteFunction
// test discreteFunction
int preTP1 = 0;
std::integral_constant<std::size_t, 0> preTP2;
auto tp = treepath(preTP1);
......@@ -137,16 +147,16 @@ int main(int argc, char** argv)
auto w1 = makeDiscreteFunction(W1, preTP2);
auto w2 = makeDiscreteFunction(W2, tp);
// test makeDOFVectorView with (pre)treepath argument
// test discreteFunction with (pre)treepath argument
auto expr = [](auto const& x) { return 1 + x[0] + x[1]; };
auto W3 = *prob.solutionVector();
auto W4 = *prob.solutionVector();
auto W5 = *prob.solutionVector();
auto W7 = *prob.solutionVector();
auto W8 = *prob.solutionVector();
auto w3 = makeDOFVectorView(W3, preTP1);
auto w4 = makeDOFVectorView(W4, preTP2);
auto w5 = makeDOFVectorView(W5, tp);
auto w3 = makeDiscreteFunction(W3, preTP1);
auto w4 = makeDiscreteFunction(W4, preTP2);
auto w5 = makeDiscreteFunction(W5, tp);
auto w6 = prob.solution(tp);
auto& W6 = *prob.solutionVector();
w3 << expr;
......@@ -158,9 +168,9 @@ int main(int argc, char** argv)
AMDIS_TEST( comp(W3, W6) );
// test interpolation on subbasis
auto w7 = makeDOFVectorView(W7);
auto w8_0 = makeDOFVectorView(W8, 0);
auto w8_1 = makeDOFVectorView(W8, 1);
auto w7 = makeDiscreteFunction(W7);
auto w8_0 = makeDiscreteFunction(W8, 0);
auto w8_1 = makeDiscreteFunction(W8, 1);
w7 << expr;
w8_0 << expr;
w8_1 << expr;
......
......@@ -9,6 +9,7 @@
#include <amdis/LinearAlgebra.hpp>
#include <amdis/GridFunctions.hpp>
#include <amdis/Integrate.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include "Tests.hpp"
......@@ -28,9 +29,9 @@ void test()
auto uVector = makeDOFVector(basis);
auto u = makeDOFVectorView(uVector);
auto u_0 = makeDOFVectorView(uVector, 0);
auto u_1 = makeDOFVectorView(uVector, 1);
auto u = makeDiscreteFunction(uVector);
auto u_0 = makeDiscreteFunction(uVector, 0);
auto u_1 = makeDiscreteFunction(uVector, 1);
// eval a functor at global coordinates (at quadrature points)
u_0 << evalAtQP([](auto const& x) { return x[0] + x[1]; });
......@@ -65,9 +66,9 @@ void test()
(integrate(partialAtQP(u_0,0) + partialAtQP(u_1,1), gridView)));
auto vVector(uVector);
auto v = makeDOFVectorView(vVector);
auto v_0 = makeDOFVectorView(vVector, 0);
auto v_1 = makeDOFVectorView(vVector, 1);
auto v = makeDiscreteFunction(vVector);
auto v_0 = makeDiscreteFunction(vVector, 0);
auto v_1 = makeDiscreteFunction(vVector, 1);
// test gradient
v << evalAtQP([](auto const& x) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment