Skip to content
Snippets Groups Projects
Commit ddd1060b authored by Happel, Lea's avatar Happel, Lea
Browse files

actually add RefinedGridFunction

parent 1afde968
No related tags found
No related merge requests found
#pragma once
#include <functional>
#include <type_traits>
#include <dune/common/referencehelper.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
// backport of dune/functions/gridfunctions/refinedgridfunction.hh
namespace AMDiS {
/**
* \brief A grid function defining data on a coarser entity level than bound to
*
* \ingroup FunctionImplementations
*
* This wrapper function allows to bind the wrapped grid function to elements on a finer
* grid than this grid function was defined on. It works by traversing from the refined element
* through the hierarchy until a coarser element is found that is in the entity-set of the
* bound grid function.
*
* \tparam GV The GridView the wrapper is bound to
* \tparam GF A GridFunction that is bound to a coarser GridView
*/
template <class GV, class GF>
class RefinedGridFunction
{
public:
/** \brief The set of entities this function can be evaluated on */
using EntitySet = Dune::Functions::GridViewEntitySet<GV,0>;
/** \brief The global coordinates */
using Domain = typename EntitySet::GlobalCoordinate;
/** \brief The result type of the function */
using Range = std::invoke_result_t<GF, Domain>;
/** \brief The local coordinates */
using LocalDomain = typename EntitySet::LocalCoordinate;
/** \brief Type of the grid element the LocalFunction can be bound to */
using Element = typename EntitySet::Element;
template <class ES, class LF>
class RefinedLocalFunction
{
public:
using Domain = LocalDomain;
using Range = std::invoke_result_t<LF, LocalDomain>;
using Element = RefinedGridFunction::Element;
public:
/** \brief Constructor. Stores the `localFct` by value. */
RefinedLocalFunction (const ES& entitySet, const LF& localFct)
: entitySet_{entitySet}
, localFct_{localFct}
{}
RefinedLocalFunction (const ES& entitySet, LF&& localFct)
: entitySet_{entitySet}
, localFct_{std::move(localFct)}
{}
RefinedLocalFunction (const RefinedLocalFunction& that) = default;
RefinedLocalFunction (RefinedLocalFunction&& that) = default;
RefinedLocalFunction& operator= (const RefinedLocalFunction& that) = default;
RefinedLocalFunction& operator= (RefinedLocalFunction&& that) = default;
/**
* \brief Bind the wrapped local-function to grid element.
*
* The given `element` might not be the element the wrapped local-function can be
* bound to. It has to be iterated in the hierarchy up (to coarser elements) until the
* entity-set of the local function contains the found father element.
*/
void bind (const Element& element)
{
element_ = element;
father_.emplace(element_);
local_ = [](const Domain& x) { return x; };
while (!entitySet_.contains(*father_) && father_->hasFather()) {
// store the coordinate transform
local_ = [local=local_, geo=father_->geometryInFather()](const Domain& x) {
return geo.global(local(x));
};
father_.emplace(father_->father());
}
assert(entitySet_.contains(*father_));
localFct_.bind(*father_);
bound_ = true;
}
/** \brief Unbind the wrapped local-function */
void unbind ()
{
localFct_.unbind();
bound_ = false;
}
/** \brief Check whether the local function is bound to an element */
bool bound () const
{
return bound_;
}
/** \brief Evaluate LocalFunction at bound element. */
Range operator() (const Domain& x) const
{
assert(bound());
return localFct_(local_(x));
}
template <class Entity>
void print (Entity const& entity) const
{
auto geo = entity.geometry();
std::cout << "element (l=" << entity.level() << ") [ ";
for (int i = 0; i < geo.corners(); ++i)
{
auto x = geo.corner(i);
std::cout << "(" << x[0] << ", " << x[1] << ") ";
}
std::cout << "]\n";
}
/** \brief Return the element this LocalFunction is bound to */
const Element& localContext () const
{
assert(bound());
return element_;
}
/** \brief Construct a derivative by wrapping the derivative of the wrapped local-function */
template <class LF_ = LF>
auto makeDerivative () const
-> RefinedLocalFunction<ES,
std::decay_t<decltype(derivative(std::declval<const LF_&>()))>>
{
auto dLocalFct = derivative(localFct_);
auto dRefinedLocalFct = RefinedLocalFunction<ES,decltype(dLocalFct)>{entitySet_, std::move(dLocalFct)};
if (bound())
dRefinedLocalFct.bind(localContext());
return dRefinedLocalFct;
}
private:
ES entitySet_;
LF localFct_;
// the following members are initialized during bind
Element element_;
std::optional<Element> father_ = std::nullopt;
std::function<Domain(const Domain&)> local_ = nullptr;
bool bound_ = false;
};
/**
* \brief Constructor.
*
* Stores the GridFunction by value. Use a reference wrapper if copying the GridFunction
* is an expensive operation!
**/
RefinedGridFunction (const GV& gridView, const GF& gridFct)
: gridView_{gridView}
, entitySet_{gridView}
, gridFct_{gridFct}
{}
RefinedGridFunction (const GV& gridView, GF&& gridFct)
: gridView_{gridView}
, entitySet_{gridView}
, gridFct_{std::move(gridFct)}
{}
/** \brief Evaluate the wrapped grid-function. */
Range operator() (const Domain& x) const
{
return gridFct_(x);
}
/** \brief Construct a derivative by wrapping the derivative of the wrapped grid-function */
template <class GF_ = GF>
auto makeDerivative () const
-> RefinedGridFunction<GV,
std::decay_t<decltype(derivative(Dune::resolveRef(std::declval<const GF_&>())))>>
{
return {gridView_, derivative(Dune::resolveRef(gridFct_))};
}
/**
* \brief Construct local function from a DiscreteGlobalBasisFunction
* \ingroup FunctionImplementations
*/
friend auto localFunction (const RefinedGridFunction& t)
{
using ES = std::decay_t<decltype(t.gridFct_.entitySet())>;
using LF = std::decay_t<decltype(localFunction(Dune::resolveRef(t.gridFct_)))>;
using LocalFunction = RefinedLocalFunction<ES,LF>;
return LocalFunction{t.gridFct_.entitySet(),
localFunction(Dune::resolveRef(t.gridFct_))};
}
/** \brief Get associated EntitySet */
const EntitySet& entitySet () const
{
return entitySet_;
}
/** \brief Get associated grid function */
const GF& gridFunction () const
{
return gridFct_;
}
GF& gridFunction ()
{
return gridFct_;
}
private:
GV gridView_;
EntitySet entitySet_;
GF gridFct_;
};
template <class F>
auto derivative (F const& f)
-> decltype(f.makeDerivative())
{
return f.makeDerivative();
}
} // end namespace AMDiS
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