Loading...
Searching...
No Matches
interpolate.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
#include <memory>
#include <vector>
#include <dune/common/exceptions.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/common/referencehelper.hh>
#include <dune/common/typetree/traversal.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/common/functionconcepts.hh>
#include <dune/functions/backends/concepts.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/functionspacebases/flatvectorview.hh>
#include <dune/functions/functionspacebases/hierarchicnodetorangemap.hh>
namespace Dune {
namespace Functions {
namespace Imp {
struct AllTrueBitSetVector
{
struct AllTrueBitSet
{
bool test(int) const { return true; }
} allTrue_;
operator bool() const
{
return true;
}
template<class I>
const AllTrueBitSetVector& operator[](const I&) const
{
return *this;
}
template<class SP>
void resize(const SP&) const
{}
};
// This helper function implements the restriction of some given function of type F.
// The restriction is a simple callback that is applied to the values of the
// function and the values of its derivative.
template<class F, class Restriction>
class ComponentFunction
{
public:
ComponentFunction(F f, Restriction restriction) :
f_(std::move(f)),
restriction_(std::move(restriction))
{}
template<class Domain>
auto operator()(const Domain& x) const
{
return restriction_(f_(x));
}
friend auto derivative(const ComponentFunction& cf)
{
// This provides support for capturing the derivative of the function by reference
// using forwardCapture for perfect forwarding capture. If the function caches its
// derivative, this saves a potentially costly copy.
auto&& df = derivative(Dune::resolveRef(cf.f_));
return [&, df=forwardCapture(std::forward<decltype(df)>(df))](auto&& x) {
return cf.restriction_(df.forward()(x));
};
}
private:
F f_;
Restriction restriction_;
};
// This helper function implements caching of the derivative for local functions.
// When using an algorithm that gets a LocalFunction and calls its derivative
// on each element, this leads to a costly call of derivative(f). E.g. for a
// DiscreteGlobalBasisFunction, this will allocate several buffer.
// To avoid this, this helper function caches the derivative and hands
// out the cached derivative by reference. To ensure that the handed out
// derivative is appropriately bound, binding the function will automatically
// bind the cached derivative.
//
// Notice that we cannot simply create the derivative in the constructor,
// because this may throw for functions that do not implement the derivative.
template<class F>
class CachedDerivativeLocalFunction
{
using Derivative = std::decay_t<decltype(derivative(Dune::resolveRef(std::declval<const F&>())))>;
public:
CachedDerivativeLocalFunction(F f) :
f_(f)
{}
template<class Element>
void bind(const Element& element)
{
Dune::resolveRef(f_).bind(element);
if (derivative_)
derivative_.value().bind(element);
}
template<class X>
auto operator()(const X& x) const
{
return f_(x);
}
friend const Derivative& derivative(const CachedDerivativeLocalFunction& cdlf)
{
if (not cdlf.derivative_)
{
auto&& lf = Dune::resolveRef(cdlf.f_);
cdlf.derivative_ = derivative(lf);
if (lf.bound())
cdlf.derivative_.value().bind(lf.localContext());
}
return cdlf.derivative_.value();
}
private:
F f_;
mutable std::optional<Derivative> derivative_;
};
template<class VectorBackend, class BitVectorBackend, class LocalFunction, class LocalView, class NodeToRangeEntry>
void interpolateLocal(VectorBackend& vector, const BitVectorBackend& bitVector, const LocalFunction& localF, const LocalView& localView, const NodeToRangeEntry& nodeToRangeEntry)
{
Dune::TypeTree::forEachLeafNode(localView.tree(), [&](auto&& node, auto&& treePath) {
if (node.empty())
return;
using Node = std::decay_t<decltype(node)>;
using FiniteElement = typename Node::FiniteElement;
using FiniteElementRangeField = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
auto&& fe = node.finiteElement();
auto localF_RE = ComponentFunction(std::cref(localF), [&](auto&& y) { return nodeToRangeEntry(node, treePath, y); });
fe.localInterpolation().interpolate(localF_RE, interpolationCoefficients);
for (size_t i=0; i<fe.localBasis().size(); ++i)
{
auto multiIndex = localView.index(node.localIndex(i));
if ( bitVector[multiIndex] )
vector[multiIndex] = interpolationCoefficients[i];
}
});
}
struct HasDerivative
{
template<class F>
auto require(F&& f) -> decltype(derivative(f));
};
} // namespace Imp
template <class B, class C, class F, class BV, class NTRE>
void interpolate(const B& basis, C&& coeff, const F& f, const BV& bv, const NTRE& nodeToRangeEntry)
{
using GridView = typename B::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalDomain = typename Element::Geometry::GlobalCoordinate;
static_assert(Dune::Functions::Concept::isCallable<F, GlobalDomain>(), "Function passed to interpolate does not model the Callable<GlobalCoordinate> concept");
auto&& gridView = basis.gridView();
// Small helper functions to wrap vectors using istlVectorBackend
// if they do not already satisfy the VectorBackend interface.
auto toVectorBackend = [&](auto& v) -> decltype(auto) {
if constexpr (models<Concept::VectorBackend<B>, decltype(v)>()) {
return v;
} else {
return istlVectorBackend(v);
}
};
auto toConstVectorBackend = [&](auto& v) -> decltype(auto) {
if constexpr (models<Concept::ConstVectorBackend<B>, decltype(v)>()) {
return v;
} else {
return istlVectorBackend(v);
}
};
auto&& bitVector = toConstVectorBackend(bv);
auto&& vector = toVectorBackend(coeff);
vector.resize(basis);
// Make a grid function supporting local evaluation out of f
auto gf = makeGridViewFunction(f, gridView);
// Obtain a local view of f
// To avoid costly reconstruction of the derivative on each element,
// we use the CachedDerivativeLocalFunction wrapper if the function
// is differentiable. This wrapper will handout
// a reference to a single cached derivative object.
auto localF = [&](){
if constexpr (models<Imp::HasDerivative, decltype(localFunction(gf))>())
return Imp::CachedDerivativeLocalFunction(localFunction(gf));
else
return localFunction(gf);
}();
auto localView = basis.localView();
for (const auto& e : elements(gridView))
{
localView.bind(e);
localF.bind(e);
Imp::interpolateLocal(vector, bitVector, localF, localView, nodeToRangeEntry);
}
}
template <class B, class C, class F, class BV>
void interpolate(const B& basis, C&& coeff, const F& f, const BV& bitVector)
{
interpolate(basis, coeff, f, bitVector, HierarchicNodeToRangeMap());
}
template <class B, class C, class F>
void interpolate(const B& basis, C&& coeff, const F& f)
{
interpolate (basis, coeff, f, Imp::AllTrueBitSetVector(), HierarchicNodeToRangeMap());
}
} // namespace Functions
} // namespace Dune
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
flatvectorview.hh
hierarchicnodetorangemap.hh
gridviewfunction.hh
istlvectorbackend.hh
functionconcepts.hh
Dune::Functions::derivative
TrigonometricFunction< K, -cosFactor, sinFactor > derivative(const TrigonometricFunction< K, sinFactor, cosFactor > &f)
Obtain derivative of TrigonometricFunction function.
Definition trigonometricfunction.hh:43
Dune::Functions::istlVectorBackend
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition istlvectorbackend.hh:350
Dune::Functions::Concept::isCallable
static constexpr auto isCallable()
Check if f is callable with given argument list.
Definition functionconcepts.hh:51
Dune
Definition monomialset.hh:19
Dune::Functions
Definition monomialset.hh:19
Dune::Functions::interpolate
void interpolate(const B &basis, C &&coeff, const F &f, const BV &bv, const NTRE &nodeToRangeEntry)
Interpolate given function in discrete function space.
Definition interpolate.hh:205
Dune::Functions::makeGridViewFunction
std::decay_t< F > makeGridViewFunction(F &&f, const GridView &gridView)
Construct a function modeling GridViewFunction from function and grid view.
Definition gridviewfunction.hh:72
Dune::Functions::forwardCapture
auto forwardCapture(T &&t)
Create a capture object for perfect forwarding.
Definition utility.hh:380
Dune::Functions::HierarchicNodeToRangeMap
A simple node to range map using the nested tree indices.
Definition hierarchicnodetorangemap.hh:34
concepts.hh