7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
13#include <dune/common/exceptions.hh>
14#include <dune/common/bitsetvector.hh>
15#include <dune/common/referencehelper.hh>
17#include <dune/typetree/traversal.hh>
19#include <dune/functions/gridfunctions/gridviewfunction.hh>
20#include <dune/functions/common/functionconcepts.hh>
22#include <dune/functions/backends/concepts.hh>
23#include <dune/functions/backends/istlvectorbackend.hh>
24#include <dune/functions/functionspacebases/flatvectorview.hh>
25#include <dune/functions/functionspacebases/hierarchicnodetorangemap.hh>
32struct AllTrueBitSetVector
36 bool test(
int)
const {
return true; }
45 const AllTrueBitSetVector& operator[](
const I&)
const
51 void resize(
const SP&)
const
61template<
class F,
class Restriction>
62class ComponentFunction
66 ComponentFunction(F f, Restriction restriction) :
68 restriction_(std::move(restriction))
71 template<
class Domain>
72 auto operator()(
const Domain& x)
const
74 return restriction_(f_(x));
77 friend auto derivative(
const ComponentFunction& cf)
82 auto&& df =
derivative(Dune::resolveRef(cf.f_));
83 return [&, df=forwardCapture(std::forward<
decltype(df)>(df))](
auto&& x) {
84 return cf.restriction_(df.forward()(x));
90 Restriction restriction_;
108class CachedDerivativeLocalFunction
110 using Derivative = std::decay_t<decltype(derivative(Dune::resolveRef(std::declval<const F&>())))>;
114 CachedDerivativeLocalFunction(F f) :
118 template<
class Element>
119 void bind(
const Element& element)
121 Dune::resolveRef(f_).bind(element);
123 derivative_.value().bind(element);
127 auto operator()(
const X& x)
const
132 friend const Derivative&
derivative(
const CachedDerivativeLocalFunction& cdlf)
134 if (not cdlf.derivative_)
136 auto&& lf = Dune::resolveRef(cdlf.f_);
139 cdlf.derivative_.value().bind(lf.localContext());
141 return cdlf.derivative_.value();
146 mutable std::optional<Derivative> derivative_;
151template<
class VectorBackend,
class BitVectorBackend,
class LocalFunction,
class LocalView,
class NodeToRangeEntry>
152void interpolateLocal(VectorBackend& vector,
const BitVectorBackend& bitVector,
const LocalFunction& localF,
const LocalView& localView,
const NodeToRangeEntry& nodeToRangeEntry)
154 Dune::TypeTree::forEachLeafNode(localView.tree(), [&](
auto&& node,
auto&& treePath) {
155 using Node = std::decay_t<decltype(node)>;
156 using FiniteElement = typename Node::FiniteElement;
157 using FiniteElementRangeField = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
159 auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
160 auto&& fe = node.finiteElement();
161 auto localF_RE = ComponentFunction(std::cref(localF), [&](auto&& y) { return nodeToRangeEntry(node, treePath, y); });
163 fe.localInterpolation().interpolate(localF_RE, interpolationCoefficients);
164 for (
size_t i=0; i<fe.localBasis().size(); ++i)
166 auto multiIndex = localView.index(node.localIndex(i));
167 if ( bitVector[multiIndex] )
168 vector[multiIndex] = interpolationCoefficients[i];
177 auto require(F&& f) ->
decltype(
derivative(f));
202template <
class B,
class C,
class F,
class BV,
class NTRE>
203void interpolate(
const B& basis, C&& coeff,
const F& f,
const BV& bv,
const NTRE& nodeToRangeEntry)
205 using GridView =
typename B::GridView;
206 using Element =
typename GridView::template Codim<0>::Entity;
207 using GlobalDomain =
typename Element::Geometry::GlobalCoordinate;
209 static_assert(Dune::Functions::Concept::isCallable<F, GlobalDomain>(),
"Function passed to interpolate does not model the Callable<GlobalCoordinate> concept");
211 auto&& gridView = basis.gridView();
215 auto toVectorBackend = [&](
auto& v) ->
decltype(
auto) {
216 if constexpr (models<Concept::VectorBackend<B>,
decltype(v)>()) {
223 auto toConstVectorBackend = [&](
auto& v) ->
decltype(
auto) {
224 if constexpr (models<Concept::ConstVectorBackend<B>,
decltype(v)>()) {
231 auto&& bitVector = toConstVectorBackend(bv);
232 auto&& vector = toVectorBackend(coeff);
233 vector.resize(basis);
236 auto gf = makeGridViewFunction(f, gridView);
244 if constexpr (models<Imp::HasDerivative, decltype(localFunction(gf))>())
245 return Imp::CachedDerivativeLocalFunction(localFunction(gf));
247 return localFunction(gf);
250 auto localView = basis.localView();
252 for (
const auto& e : elements(gridView))
256 Imp::interpolateLocal(vector, bitVector, localF, localView, nodeToRangeEntry);
276template <
class B,
class C,
class F,
class BV>
277void interpolate(
const B& basis, C&& coeff,
const F& f,
const BV& bitVector)
279 interpolate(basis, coeff, f, bitVector, HierarchicNodeToRangeMap());
296template <
class B,
class C,
class F>
297void interpolate(
const B& basis, C&& coeff,
const F& f)
299 interpolate (basis, coeff, f, Imp::AllTrueBitSetVector(), HierarchicNodeToRangeMap());
TrigonometricFunction< K, -cosFactor, sinFactor > derivative(const TrigonometricFunction< K, sinFactor, cosFactor > &f)
Obtain derivative of TrigonometricFunction function.
Definition: trigonometricfunction.hh:43
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition: istlvectorbackend.hh:350
Definition: polynomial.hh:17