3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
12#include <dune/typetree/childextraction.hh>
14#include <dune/functions/gridfunctions/gridviewfunction.hh>
15#include <dune/functions/common/functionconcepts.hh>
17#include <dune/functions/backends/concepts.hh>
18#include <dune/functions/backends/istlvectorbackend.hh>
19#include <dune/functions/functionspacebases/sizeinfo.hh>
20#include <dune/functions/functionspacebases/flatvectorview.hh>
21#include <dune/functions/functionspacebases/hierarchicnodetorangemap.hh>
23#include <dune/typetree/traversal.hh>
24#include <dune/typetree/visitor.hh>
31struct AllTrueBitSetVector
35 bool test(
int)
const {
return true; }
44 const AllTrueBitSetVector& operator[](
const I&)
const
50 void resize(
const SP&)
const
57template <
class B,
class T,
class NTRE,
class HV,
class LF,
class HBV>
58class LocalInterpolateVisitor
59 :
public TypeTree::TreeVisitor
60 ,
public TypeTree::DynamicTraversal
66 using LocalView =
typename B::LocalView;
67 using MultiIndex =
typename LocalView::MultiIndex;
69 using LocalFunction = LF;
73 using VectorBackend = HV;
74 using BitVectorBackend = HBV;
76 using NodeToRangeEntry = NTRE;
78 using GridView =
typename Basis::GridView;
79 using Element =
typename GridView::template Codim<0>::Entity;
81 using LocalDomain =
typename Element::Geometry::LocalCoordinate;
83 using GlobalDomain =
typename Element::Geometry::GlobalCoordinate;
85 LocalInterpolateVisitor(
const B& , HV& coeff,
const HBV& bitVector,
const LF& localF,
const LocalView& localView,
const NodeToRangeEntry& nodeToRangeEntry) :
88 bitVector_(bitVector),
89 localView_(localView),
90 nodeToRangeEntry_(nodeToRangeEntry)
92 static_assert(Dune::Functions::Concept::isCallable<LocalFunction, LocalDomain>(),
"Function passed to LocalInterpolateVisitor does not model the Callable<LocalCoordinate> concept");
95 template<
typename Node,
typename TreePath>
96 void pre(Node&, TreePath)
99 template<
typename Node,
typename TreePath>
100 void post(Node&, TreePath)
103 template<
typename Node,
typename TreePath>
106 using FiniteElement =
typename Node::FiniteElement;
107 using FiniteElementRange =
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
108 using FiniteElementRangeField =
typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
110 auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
111 auto&& fe = node.finiteElement();
116 if constexpr ( FiniteElement::Traits::LocalBasisType::Traits::dimRange == 1 )
123 auto localFj = [&](
const LocalDomain& x){
124 const auto& y = localF_(x);
125 return FiniteElementRange(flatVectorView(nodeToRangeEntry_(node,
treePath, y))[j]);
131 auto blockSize = flatVectorView(vector_[localView_.index(0)]).size();
133 for(j=0; j<blockSize; ++j)
135 fe.localInterpolation().interpolate(localFj, interpolationCoefficients);
136 for (
size_t i=0; i<fe.localBasis().size(); ++i)
138 auto multiIndex = localView_.index(node.localIndex(i));
139 auto bitVectorBlock = flatVectorView(bitVector_[multiIndex]);
140 if (bitVectorBlock[j])
142 auto vectorBlock = flatVectorView(vector_[multiIndex]);
143 vectorBlock[j] = interpolationCoefficients[i];
151 auto localF = [&](
const LocalDomain& x){
152 const auto& y = localF_(x);
153 return FiniteElementRange(nodeToRangeEntry_(node,
treePath, y));
156 fe.localInterpolation().interpolate(localF, interpolationCoefficients);
157 for (
size_t i=0; i<fe.localBasis().size(); ++i)
159 auto multiIndex = localView_.index(node.localIndex(i));
160 if ( bitVector_[multiIndex] )
162 vector_[multiIndex] = interpolationCoefficients[i];
171 VectorBackend& vector_;
172 const LocalFunction& localF_;
173 const BitVectorBackend& bitVector_;
174 const LocalView& localView_;
175 const NodeToRangeEntry& nodeToRangeEntry_;
201template <
class B,
class C,
class F,
class BV,
class NTRE>
202void interpolate(
const B& basis, C&& coeff,
const F& f,
const BV& bv,
const NTRE& nodeToRangeEntry)
204 using GridView =
typename B::GridView;
205 using Element =
typename GridView::template Codim<0>::Entity;
207 using Tree =
typename B::LocalView::Tree;
209 using GlobalDomain =
typename Element::Geometry::GlobalCoordinate;
211 static_assert(Dune::Functions::Concept::isCallable<F, GlobalDomain>(),
"Function passed to interpolate does not model the Callable<GlobalCoordinate> concept");
213 auto&& gridView = basis.gridView();
217 auto toVectorBackend = [&](
auto& v) ->
decltype(
auto) {
218 if constexpr (models<Concept::VectorBackend<B>,
decltype(v)>()) {
225 auto toConstVectorBackend = [&](
auto& v) ->
decltype(
auto) {
226 if constexpr (models<Concept::ConstVectorBackend<B>,
decltype(v)>()) {
233 auto&& bitVector = toConstVectorBackend(bv);
234 auto&& vector = toVectorBackend(coeff);
235 vector.resize(sizeInfo(basis));
240 auto gf = makeGridViewFunction(f, gridView);
243 auto localF = localFunction(gf);
245 auto localView = basis.localView();
247 for (
const auto& e : elements(gridView))
252 Imp::LocalInterpolateVisitor<B, Tree, NTRE,
decltype(vector),
decltype(localF),
decltype(bitVector)> localInterpolateVisitor(basis, vector, bitVector, localF, localView, nodeToRangeEntry);
273template <
class B,
class C,
class F,
class BV>
274void interpolate(
const B& basis, C&& coeff,
const F& f,
const BV& bitVector)
276 interpolate(basis, coeff, f, bitVector, HierarchicNodeToRangeMap());
293template <
class B,
class C,
class F>
294void interpolate(
const B& basis, C&& coeff,
const F& f)
296 interpolate (basis, coeff, f, Imp::AllTrueBitSetVector(), HierarchicNodeToRangeMap());
Efficient implementation of a dynamic array of static arrays of booleans.
A few common exception classes.
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition: istlvectorbackend.hh:346
constexpr HybridTreePath< T... > treePath(const T &... t)
Constructs a new HybridTreePath from the given indices.
Definition: treepath.hh:191
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:237
Dune namespace.
Definition: alignedallocator.hh:13
void post(T &&, TreePath) const
Method for postfix tree traversal.
Definition: visitor.hh:81
void leaf(T &&, TreePath) const
Method for leaf traversal.
Definition: visitor.hh:91
void pre(T &&, TreePath) const
Method for prefix tree traversal.
Definition: visitor.hh:58