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>
27 #include <dune/typetree/traversal.hh>
28 #include <dune/typetree/visitor.hh>
35 struct AllTrueBitSetVector
39 bool test(
int i)
const {
return true; }
48 const AllTrueBitSetVector& operator[](
const I&)
const
54 void resize(
const SP&)
const
61 template <
class B,
class T,
class NTRE,
class HV,
class LF,
class HBV>
62 class 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);
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));
145 if (bitVectorBlock[j])
148 vectorBlock[j] = interpolationCoefficients[i];
157 VectorBackend& vector_;
158 const LocalFunction& localF_;
159 const BitVectorBackend& bitVector_;
160 const LocalView& localView_;
161 const NodeToRangeEntry& nodeToRangeEntry_;
187 template <
class B,
class C,
class F,
class BV,
class NTRE>
188 void 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) {
205 [&](
auto id) -> decltype(
auto) {
207 }, [&](
auto id) -> decltype(
auto) {
211 auto toConstVectorBackend = [&](
auto& v) -> decltype(
auto) {
213 [&](
auto id) -> decltype(
auto) {
215 }, [&](
auto id) -> decltype(
auto) {
220 auto&& bitVector = toConstVectorBackend(bv);
221 auto&& vector = toVectorBackend(coeff);
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);
260 template <
class B,
class C,
class F,
class BV>
261 void interpolate(
const B& basis, C&& coeff,
const F& f,
const BV& bitVector)
280 template <
class B,
class C,
class F>
289 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH