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>
12 #include <dune/typetree/childextraction.hh>
23 #include <dune/typetree/traversal.hh>
24 #include <dune/typetree/visitor.hh>
31 struct AllTrueBitSetVector
35 bool test(
int i)
const {
return true; }
44 const AllTrueBitSetVector& operator[](
const I&)
const
50 void resize(
const SP&)
const
57 template <
class B,
class T,
class NTRE,
class HV,
class LF,
class HBV>
58 class 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& basis, 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& node, TreePath treePath)
99 template<
typename Node,
typename TreePath>
100 void post(Node& node, TreePath treePath)
103 template<
typename Node,
typename TreePath>
104 void leaf(Node& node, TreePath 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));
140 if (bitVectorBlock[j])
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_;
201 template <
class B,
class C,
class F,
class BV,
class NTRE>
202 void 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) {
225 auto toConstVectorBackend = [&](
auto& v) -> decltype(
auto) {
233 auto&& bitVector = toConstVectorBackend(bv);
234 auto&& vector = toVectorBackend(coeff);
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);
253 TypeTree::applyToTree(localView.tree(),localInterpolateVisitor);
273 template <
class B,
class C,
class F,
class BV>
274 void interpolate(
const B& basis, C&& coeff,
const F& f,
const BV& bitVector)
293 template <
class B,
class C,
class F>
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
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:202
std::decay< F >::type makeGridViewFunction(F &&f, const GridView &gridView)
Construct a function modeling GridViewFunction from function and grid view.
Definition: gridviewfunction.hh:68
auto flatVectorView(T &t)
Create flat vector view of passed mutable container.
Definition: flatvectorview.hh:159
SizeInfo< Basis > sizeInfo(const Basis &basis)
Definition: sizeinfo.hh:69
Definition: backends/concepts.hh:21
Definition: backends/concepts.hh:31
A simple node to range map using the nested tree indices.
Definition: hierarchicnodetorangemap.hh:30