DUNE PDELab (git)

interpolate.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4#ifndef DUNE_PDELAB_GRIDFUNCTIONSPACE_INTERPOLATE_HH
5#define DUNE_PDELAB_GRIDFUNCTIONSPACE_INTERPOLATE_HH
6
7#include <vector>
8#include <utility>
9
11
12#include <dune/localfunctions/common/interfaceswitch.hh>
13#include <dune/localfunctions/common/virtualinterface.hh>
14
15#include <dune/typetree/typetree.hh>
16#include <dune/typetree/pairtraversal.hh>
17
18#include <dune/pdelab/function/localfunction.hh>
19#include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
20#include <dune/pdelab/gridfunctionspace/lfsindexcache.hh>
21
22namespace Dune {
23 namespace PDELab {
24
28 namespace {
29
30 template<typename LF, typename XG>
31 struct InterpolateLeafFromScalarVisitor
32 : public TypeTree::TreeVisitor
33 , public TypeTree::DynamicTraversal
34 {
35
36 template<typename LFS, typename TreePath>
37 void leaf(const LFS& lfs, TreePath treePath) const
38 {
39 std::vector<typename XG::ElementType> xl(lfs.size());
40
41 // call interpolate for the basis
43 interpolate(lf, xl);
44
45 // write coefficients into local vector
46 xg.write_sub_container(lfs,xl);
47 }
48
49 InterpolateLeafFromScalarVisitor(const LF& lf_, XG& xg_)
50 : lf(lf_)
51 , xg(xg_)
52 {}
53
54 const LF& lf;
55 XG& xg;
56 };
57
58
59 template<typename LF, typename XG>
60 struct InterpolateLeafFromVectorVisitor
61 : public TypeTree::TreeVisitor
62 , public TypeTree::DynamicTraversal
63 {
64 using Domain = typename Functions::SignatureTraits<LF>::Domain;
65 using Range = typename Functions::SignatureTraits<LF>::Range;
66 using RangeField = typename FieldTraits<Range>::field_type;
67
68 template<typename LFS, typename TreePath>
69 void leaf(const LFS& lfs, TreePath treePath) const
70 {
71 std::vector<typename XG::ElementType> xl(lfs.size());
72
73 using LFSRange = typename LFS::Traits::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
74 static_assert(std::is_convertible<LFSRange, typename FieldTraits< LFSRange >::field_type>::value,
75 "only interpolation into scalar leaf function spaces is implemented");
76
77 // call interpolate for the basis
78 auto f = [&](const Domain& x) -> RangeField { return lf(x)[index]; };
79
81 interpolate(f, xl);
82
83 // write coefficients into local vector
84 xg.write_sub_container(lfs,xl);
85
86 // increment index
87 index++;
88 }
89
90 InterpolateLeafFromVectorVisitor(const LF& lf_, XG& xg_)
91 : lf(lf_)
92 , index(0)
93 , xg(xg_)
94 {}
95
96 const LF& lf;
97 mutable std::size_t index;
98 XG& xg;
99 };
100
101
102 template<typename E, typename XG>
103 struct InterpolateVisitor
104 : public TypeTree::TreePairVisitor
105 , public TypeTree::DynamicTraversal
106 {
107
108 template<typename F, typename LFS, typename TreePath>
109 typename std::enable_if<F::isLeaf && LFS::isLeaf>::type
110 leaf(const F& f, const LFS& lfs, TreePath treePath) const
111 {
112 if (lfs.size() == 0)
113 return;
114
115 std::vector<typename XG::ElementType> xl(lfs.size());
117 interpolation(lfs.finiteElement())
118 .interpolate(f, xl);
119
120 // write coefficients into local vector
121 xg.write_sub_container(lfs, xl);
122 }
123
124 // interpolate PowerLFS from scalar-valued function
125 template<typename F, typename LFS, typename TreePath,
126 typename Range = typename Functions::SignatureTraits<F>::Range>
127 typename std::enable_if<F::isLeaf &&
128 std::is_convertible<Range, typename FieldTraits< Range >::field_type>::value &&
129 (!LFS::isLeaf)>::type
130 leaf(const F& f, const LFS& lfs, TreePath treePath) const
131 {
132 // call interpolate for the basis
133 TypeTree::applyToTree(lfs,InterpolateLeafFromScalarVisitor<F,XG>(f, xg));
134 }
135
136 // interpolate PowerLFS from vector-valued function
137 template<typename F, typename LFS, typename TreePath,
138 typename Range = typename Functions::SignatureTraits<F>::Range>
139 typename std::enable_if<F::isLeaf &&
140 (!std::is_convertible<Range, typename FieldTraits< Range >::field_type>::value) &&
141 (!LFS::isLeaf)>::type
142 leaf(const F& f, const LFS& lfs, TreePath treePath) const
143 {
144 static_assert(TypeTree::TreeInfo<LFS>::leafCount == Range::dimension,
145 "Number of leaves and dimension of range type " \
146 "must match for automatic interpolation of " \
147 "vector-valued function");
148
149 TypeTree::applyToTree(lfs,InterpolateLeafFromVectorVisitor<F,XG>(f,xg));
150 }
151
152 InterpolateVisitor(const E& e_, XG& xg_)
153 : e(e_)
154 , xg(xg_)
155 {}
156
157 private:
158 const E& e;
159 XG& xg;
160 };
161
162 } // anonymous namespace
163
165
176 template<typename F, typename GFS, typename XG>
177 void interpolate (const F& f, const GFS& gfs, XG& xg)
178 {
179 // this is the leaf version now
180
181 // get some types
182 using EntitySet = typename GFS::Traits::EntitySet;
183 using Element = typename EntitySet::Element;
184
185 auto entity_set = gfs.entitySet();
186
187 // make local function
188 auto lf = makeLocalFunctionTree(f, gfs.gridView());
189
190 // make local function space
191 typedef LocalFunctionSpace<GFS> LFS;
192 LFS lfs(gfs);
193 typedef LFSIndexCache<LFS> LFSCache;
194 LFSCache lfs_cache(lfs);
195 typedef typename XG::template LocalView<LFSCache> XView;
196
197 XView x_view(xg);
198
199 // loop once over the grid
200 for (const auto& element : elements(entity_set))
201 {
202 // bind local function space to element
203 lf.bind(element);
204 lfs.bind(element);
205 lfs_cache.update();
206 x_view.bind(lfs_cache);
207
208 // call interpolate
209 TypeTree::applyToTreePair(lf,lfs,InterpolateVisitor<Element,XView>(element,x_view));
210
211 x_view.unbind();
212 lf.unbind();
213 }
214
215 x_view.detach();
216 }
217
219 } // namespace PDELab
220} // namespace Dune
221
222#endif // DUNE_PDELAB_GRIDFUNCTIONSPACE_INTERPOLATE_HH
A few common exception classes.
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
constexpr auto treePath(const T &... t)
Constructs a new HybridTreePath from the given indices.
Definition: treepath.hh:326
void applyToTreePair(Tree1 &&tree1, Tree2 &&tree2, Visitor &&visitor)
Apply visitor to a pair of TypeTrees.
Definition: pairtraversal.hh:128
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:239
Dune namespace.
Definition: alignedallocator.hh:13
static const Interpolation & interpolation(const FiniteElement &fe)
access interpolation
Definition: interfaceswitch.hh:42
static const std::size_t leafCount
The number of leaf nodes in the TypeTree.
Definition: utility.hh:136
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)