3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
9#include <dune/common/exceptions.hh>
10#include <dune/common/bitsetvector.hh>
11#include <dune/common/hybridutilities.hh>
13#include <dune/typetree/childextraction.hh>
15#include <dune/localfunctions/common/virtualinterface.hh>
17#include <dune/functions/gridfunctions/gridviewfunction.hh>
18#include <dune/functions/common/functionfromcallable.hh>
19#include <dune/functions/common/functionconcepts.hh>
21#include <dune/functions/backends/concepts.hh>
22#include <dune/functions/backends/istlvectorbackend.hh>
23#include <dune/functions/functionspacebases/sizeinfo.hh>
24#include <dune/functions/functionspacebases/flatvectorview.hh>
25#include <dune/functions/functionspacebases/hierarchicnodetorangemap.hh>
27#include <dune/typetree/traversal.hh>
28#include <dune/typetree/visitor.hh>
35struct AllTrueBitSetVector
39 bool test(
int i)
const {
return true; }
48 const AllTrueBitSetVector& operator[](
const I&)
const
54 void resize(
const SP&)
const
61template <
class B,
class T,
class NTRE,
class HV,
class LF,
class HBV>
62class LocalInterpolateVisitor
63 :
public TypeTree::TreeVisitor
64 ,
public TypeTree::DynamicTraversal
70 using LocalView =
typename B::LocalView;
71 using MultiIndex =
typename LocalView::MultiIndex;
73 using LocalFunction = LF;
77 using VectorBackend = HV;
78 using BitVectorBackend = HBV;
80 using NodeToRangeEntry = NTRE;
82 using GridView =
typename Basis::GridView;
83 using Element =
typename GridView::template Codim<0>::Entity;
85 using LocalDomain =
typename Element::Geometry::LocalCoordinate;
87 using GlobalDomain =
typename Element::Geometry::GlobalCoordinate;
89 LocalInterpolateVisitor(
const B& basis, HV& coeff,
const HBV& bitVector,
const LF& localF,
const LocalView& localView,
const NodeToRangeEntry& nodeToRangeEntry) :
92 bitVector_(bitVector),
93 localView_(localView),
94 nodeToRangeEntry_(nodeToRangeEntry)
96 static_assert(Dune::Functions::Concept::isCallable<LocalFunction, LocalDomain>(),
"Function passed to LocalInterpolateVisitor does not model the Callable<LocalCoordinate> concept");
99 template<
typename Node,
typename TreePath>
100 void pre(Node& node, TreePath treePath)
103 template<
typename Node,
typename TreePath>
104 void post(Node& node, TreePath treePath)
107 template<
typename Node,
typename TreePath>
108 void leaf(Node& node, TreePath treePath)
110 using FiniteElement =
typename Node::FiniteElement;
111 using FiniteElementRange =
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
112 using FiniteElementRangeField =
typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
113 using FunctionBaseClass =
typename Dune::LocalFiniteElementFunctionBase<FiniteElement>::type;
120 auto localFj = [&](
const LocalDomain& x){
121 const auto& y = localF_(x);
122 return flatVectorView(nodeToRangeEntry_(node, treePath, y))[j];
125 using FunctionFromCallable =
typename Dune::Functions::FunctionFromCallable<FiniteElementRange(LocalDomain),
decltype(localFj), FunctionBaseClass>;
127 auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
129 auto&& fe = node.finiteElement();
134 auto blockSize = flatVectorView(vector_[localView_.index(0)]).size();
136 for(j=0; j<blockSize; ++j)
140 fe.localInterpolation().interpolate(FunctionFromCallable(localFj), interpolationCoefficients);
141 for (
size_t i=0; i<fe.localBasis().size(); ++i)
143 auto multiIndex = localView_.index(node.localIndex(i));
144 auto bitVectorBlock = flatVectorView(bitVector_[multiIndex]);
145 if (bitVectorBlock[j])
147 auto vectorBlock = flatVectorView(vector_[multiIndex]);
148 vectorBlock[j] = interpolationCoefficients[i];
157 VectorBackend& vector_;
158 const LocalFunction& localF_;
159 const BitVectorBackend& bitVector_;
160 const LocalView& localView_;
161 const NodeToRangeEntry& nodeToRangeEntry_;
187template <
class B,
class C,
class F,
class BV,
class NTRE>
188void interpolate(
const B& basis, C&& coeff,
const F& f,
const BV& bv,
const NTRE& nodeToRangeEntry)
190 using GridView =
typename B::GridView;
191 using Element =
typename GridView::template Codim<0>::Entity;
193 using Tree =
typename B::LocalView::Tree;
195 using GlobalDomain =
typename Element::Geometry::GlobalCoordinate;
197 static_assert(Dune::Functions::Concept::isCallable<F, GlobalDomain>(),
"Function passed to interpolate does not model the Callable<GlobalCoordinate> concept");
199 auto&& gridView = basis.gridView();
203 auto toVectorBackend = [&](
auto& v) ->
decltype(
auto) {
204 return Dune::Hybrid::ifElse(models<Concept::VectorBackend<B>,
decltype(v)>(),
205 [&](
auto id) ->
decltype(
auto) {
207 }, [&](
auto id) ->
decltype(
auto) {
211 auto toConstVectorBackend = [&](
auto& v) ->
decltype(
auto) {
212 return Dune::Hybrid::ifElse(models<Concept::ConstVectorBackend<B>,
decltype(v)>(),
213 [&](
auto id) ->
decltype(
auto) {
215 }, [&](
auto id) ->
decltype(
auto) {
220 auto&& bitVector = toConstVectorBackend(bv);
221 auto&& vector = toVectorBackend(coeff);
222 vector.resize(sizeInfo(basis));
227 auto gf = makeGridViewFunction(f, gridView);
232 auto localView = basis.localView();
234 for (
const auto& e : elements(gridView))
239 Imp::LocalInterpolateVisitor<B, Tree, NTRE,
decltype(vector),
decltype(localF),
decltype(bitVector)> localInterpolateVisitor(basis, vector, bitVector, localF, localView, nodeToRangeEntry);
240 TypeTree::applyToTree(localView.tree(),localInterpolateVisitor);
260template <
class B,
class C,
class F,
class BV>
261void interpolate(
const B& basis, C&& coeff,
const F& f,
const BV& bitVector)
263 interpolate(basis, coeff, f, bitVector, HierarchicNodeToRangeMap());
280template <
class B,
class C,
class F>
281void interpolate(
const B& basis, C&& coeff,
const F& f)
283 interpolate (basis, coeff, f, Imp::AllTrueBitSetVector(), HierarchicNodeToRangeMap());
void localFunction(DiscreteGlobalBasisFunction< TT... > &&t)=delete
Construction of local functions from a temporary DiscreteGlobalBasisFunction (forbidden)
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition: istlvectorbackend.hh:345
Definition: polynomial.hh:10