AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
InterpolationDataTransfer.inc.hpp
1#pragma once
2
3#include <cmath>
4#include <functional>
5#include <limits>
6#include <map>
7#include <memory>
8#include <numeric>
9#include <type_traits>
10#include <unordered_map>
11#include <utility>
12#include <vector>
13
14#include <dune/common/fvector.hh>
15#include <dune/common/hash.hh>
16#include <dune/common/version.hh>
17
18#include <dune/grid/common/geometry.hh>
19#include <dune/grid/common/mcmgmapper.hh>
20#include <dune/grid/common/rangegenerators.hh>
21
22#include <amdis/Output.hpp>
23#include <amdis/common/ConcurrentCache.hpp>
24#include <amdis/functions/EntitySet.hpp>
25#include <amdis/operations/Assigner.hpp>
26#include <amdis/typetree/Traversal.hpp>
27
28namespace AMDiS {
29namespace Impl {
30
31// Hash function for cache container
32struct CoordHasher
33{
34 template <class LocalCoord>
35 std::size_t operator()(LocalCoord const& coord) const
36 {
37 std::size_t seed = 0;
38 for (std::size_t i = 0; i < coord.size(); ++i)
39 Dune::hash_combine(seed, coord[i]);
40 return seed;
41 }
42};
43
44} // end namespace Impl
45
46
47template <class B, class C>
49preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
50{
51 GridView gv = basis.gridView();
52 LocalView lv = basis.localView();
53 auto const& grid = gv.grid();
54 auto const& idSet = grid.localIdSet();
55
56 nodeDataTransfer_.resize(lv.tree());
57 Traversal::forEachLeafNode(lv.tree(), [&](auto const& node, auto const& tp) {
58 nodeDataTransfer_[tp].preAdaptInit(lv, coeff, node);
59 });
60
61 // Make persistent DoF container
62 persistentContainer_.clear(); // Redundant if postAdapt was correctly called last cycle
63 for (const auto& e : elements(gv))
64 {
65 auto it = persistentContainer_.emplace(idSet.id(e),
66 Dune::TypeTree::makeTreeContainer(lv.tree(), [](auto&& node) { return NodeElementData<TYPEOF(node)>{}; }));
67
68 lv.bind(e);
69 auto& treeContainer = it.first->second;
70 Traversal::forEachLeafNode(lv.tree(), [&](auto const& /*node*/, auto const& tp) {
71 nodeDataTransfer_[tp].cacheLocal(treeContainer[tp]);
72 });
73 }
74
75 if (!mightCoarsen)
76 return;
77
78 // Interpolate from possibly vanishing elements
79 auto maxLevel = grid.maxLevel();
80 using std::sqrt;
81 typename Grid::ctype const checkInsideTolerance = sqrt(std::numeric_limits<typename Grid::ctype>::epsilon());
82 using Seed = typename Grid::template Codim<0>::EntitySeed;
83 std::vector<Seed> seeds(maxLevel+1);
84
85 for (auto const& e : entitySet(basis))
86 {
87 auto father = e;
88 while (father.mightVanish() && father.hasFather())
89 {
90 father = father.father();
91 auto it = persistentContainer_.emplace(idSet.id(father),
92 Dune::TypeTree::makeTreeContainer(lv.tree(), [](auto&& node) { return NodeElementData<TYPEOF(node)>{}; }));
93 if (!it.second)
94 continue;
95
96 auto& treeContainer = it.first->second;
97 bool init = true; // init flag for first call on new father element
98 bool restrictLocalCompleted = false;
99 auto hItEnd = father.hend(maxLevel);
100 for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt)
101 {
102 seeds[hIt->level()] = hIt->seed(); // Save element in hierarchy to access geometryInFather at each step later
103
104 if (!hIt->isLeaf())
105 continue;
106
107 auto const& child = *hIt;
108 auto search = persistentContainer_.find(idSet.id(child));
109 assert(search != persistentContainer_.end());
110 auto const& childContainer = search->second;
111 lv.bind(child);
112
113 using BoolCoordPair = std::pair<bool, LocalCoordinate>;
114 using CacheImp = std::unordered_map<LocalCoordinate, BoolCoordPair, Impl::CoordHasher>;
115 using ChildCache = ConcurrentCache<LocalCoordinate, BoolCoordPair, ConsecutivePolicy, CacheImp>;
116
117 // Transfers input father-local point x into child-local point y
118 // Returns false if x is not inside the child
119 auto xInChild = [&](LocalCoordinate const& x) -> BoolCoordPair {
120 LocalCoordinate y = x;
121 for (int i = father.level()+1; i <= child.level(); ++i)
122 {
123 auto currentElement = grid.entity(seeds[i]);
124 auto const& geoInFather = currentElement.geometryInFather();
125 y = geoInFather.local(y);
126 // TODO(FM): Using an implementation detail as workaround for insufficient
127 // tolerance, see https://gitlab.dune-project.org/core/dune-grid/issues/84
128 bool isInside = Dune::Geo::Impl::checkInside(Dune::referenceElement(geoInFather).type().id(), Geometry::mydimension, y, checkInsideTolerance);
129 if (!isInside)
130 return BoolCoordPair(false, std::move(y));
131 }
132 return BoolCoordPair(true, std::move(y));
133 };
134 // Cache result of xInChild for subsequent calls from other basis nodes
135 // TODO(FM): Disable for single-node basis
136 ChildCache childCache;
137 auto xInChildCached = [&](LocalCoordinate const& x) -> BoolCoordPair {
138 return childCache.get(x, [&](LocalCoordinate const& x) { return xInChild(x); });
139 };
140
141 restrictLocalCompleted = true;
142 Traversal::forEachLeafNode(lv.tree(), [&](auto const& /*node*/, auto const& tp) {
143 restrictLocalCompleted &=
144 nodeDataTransfer_[tp].restrictLocal(father, treeContainer[tp], xInChildCached,
145 childContainer[tp], init);
146 });
147 init = false;
148 }
149 // test if restrictLocal was completed on all nodes
150 assert(restrictLocalCompleted);
151
152 } // end while (father.mightVanish)
153 } // end for (elements)
154}
155
156
157template <class B, class C>
158void InterpolationDataTransfer<B,C>::adapt(B const& basis, C& coeff)
159{
160 coeff.resize(basis);
161
162 // No data was saved before adapting the grid, make
163 // sure to call DataTransfer::preAdapt before calling adapt() on the grid
164 if (persistentContainer_.empty())
165 return;
166
167 GridView gv = basis.gridView();
168 LocalView lv = basis.localView();
169 auto const& idSet = gv.grid().localIdSet();
170 Traversal::forEachLeafNode(lv.tree(), [&](auto const& node, auto const& tp) {
171 nodeDataTransfer_[tp].adaptInit(lv, coeff, node);
172 });
173
174 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
175 Mapper mapper{gv, Dune::mcmgElementLayout()};
176
177 std::vector<bool> finished(mapper.size(), false);
178 for (const auto& e : entitySet(basis))
179 {
180 auto index = mapper.index(e);
181 if (finished[index])
182 continue;
183
184 auto it = persistentContainer_.find(idSet.id(e));
185
186 // Data already exists and no interpolation is required
187 if (it != persistentContainer_.end()) {
188 lv.bind(e);
189 auto const& treeContainer = it->second;
190 Traversal::forEachLeafNode(lv.tree(), [&](auto const& /*node*/, auto const& tp) {
191 nodeDataTransfer_[tp].copyLocal(treeContainer[tp]);
192 });
193 finished[index] = true;
194 continue;
195 }
196
197 // Data needs to be interpolated
198 auto father = e;
199 while (father.hasFather() && father.isNew())
200 father = father.father();
201
202 auto maxLevel = gv.grid().maxLevel();
203 bool init = true; // init flag for first call on new father element
204
205 auto father_it = persistentContainer_.find(idSet.id(father));
206 assert(father_it != persistentContainer_.end());
207 auto const& treeContainer = father_it->second;
208
209 auto hItEnd = father.hend(maxLevel);
210 for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt) {
211 if (!hIt->isLeaf())
212 continue;
213
214 auto const& child = *hIt;
215 lv.bind(child);
216
217 // coordinate transform from child to father element
218 auto xInFather = [&](LocalCoordinate const& x) -> LocalCoordinate
219 {
220 auto y = x;
221 auto currentElement = child;
222 while (currentElement.level() != father.level())
223 {
224 y = currentElement.geometryInFather().global(y);
225 currentElement = currentElement.father();
226 }
227 return y;
228 };
229
230 Traversal::forEachLeafNode(lv.tree(), [&](auto const& /*node*/, auto const& tp) {
231 nodeDataTransfer_[tp].prolongLocal(father, treeContainer[tp], xInFather, init);
232 });
233
234 finished[mapper.index(child)] = true;
235 init = false;
236 }
237 } // end for (elements)
238
239 coeff.finish();
240}
241
242
243template <class B, class C>
245{
246 persistentContainer_.clear();
247}
248
249
253template <class Node, class Container, class Basis>
255{
256 using T = typename Container::value_type;
257 using LocalView = typename Basis::LocalView;
258 using Element = typename Node::Element;
259
260 using LocalBasis = typename Node::FiniteElement::Traits::LocalBasisType;
261 using LBRangeType = typename LocalBasis::Traits::RangeType;
262 using LocalInterpolation = typename Node::FiniteElement::Traits::LocalBasisType;
263 using LIDomainType = typename LocalInterpolation::Traits::DomainType;
264 using LIRangeType = typename LocalInterpolation::Traits::RangeType;
265
266public:
267 using NodeElementData = std::vector<T>;
268
269public:
270 NodeDataTransfer() = default;
271
273 void preAdaptInit(LocalView const& lv, Container const& coeff, Node const& node)
274 {
275 lv_ = &lv;
276 node_ = &node;
277 fatherNode_ = std::make_shared<Node>(node);
278 constCoeff_ = &coeff;
279 }
280
282
286 // [[expects: preAdaptInit to be called before]]
287 void cacheLocal(NodeElementData& dofs) const
288 {
289 constCoeff_->gather(*lv_, *node_, dofs);
290 }
291
306 // [[expects: preAdaptInit to be called before]]
307 template <class Trafo>
308 bool restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& trafo,
309 NodeElementData const& childDOFs, bool init);
310
311
313 void adaptInit(LocalView const& lv, Container& coeff, Node const& node)
314 {
315 lv_ = &lv;
316 node_ = &node;
317 fatherNode_ = std::make_shared<Node>(node);
318 coeff_ = &coeff;
319 }
320
322 // [[expects: adaptInit to be called before]]
323 void copyLocal(NodeElementData const& dofs) const
324 {
325 coeff_->scatter(*lv_, *node_, dofs, Assigner::assign{});
326 }
327
339 // [[expects: adaptInit to be called before]]
340 template <class Trafo>
341 void prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
342 Trafo const& trafo, bool init);
343
344private:
345 LocalView const* lv_ = nullptr;
346 Node const* node_ = nullptr;
347 std::shared_ptr<Node> fatherNode_;
348 Container const* constCoeff_ = nullptr;
349 Container* coeff_ = nullptr;
350 std::vector<bool> finishedDOFs_;
351 NodeElementData fatherDOFsTemp_;
352};
353
354
355template <class N, class C, class B>
356 template <class Trafo>
358restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& trafo,
359 NodeElementData const& childDOFs, bool init)
360{
361 auto& fatherNode = *fatherNode_;
362 std::size_t currentDOF = 0;
363 if (init)
364 {
365 // TODO(FM): This is UB, see https://gitlab.com/amdis/amdis/-/issues/16 (case 2)
366 bindTree(fatherNode, father);
367 }
368 auto const& childNode = *node_;
369 auto const& childFE = childNode.finiteElement();
370 auto const& fatherFE = fatherNode.finiteElement();
371
372 if (init) {
373 finishedDOFs_.assign(fatherFE.size(), false);
374 fatherDOFsTemp_.assign(fatherFE.size(), 0);
375 }
376
377 auto evalLeaf = [&](LIDomainType const& x) -> LIRangeType {
378 if (!finishedDOFs_[currentDOF])
379 {
380 auto const& insideLocal = trafo(x);
381 bool isInside = insideLocal.first;
382 if (isInside)
383 {
384 auto const& local = insideLocal.second;
385 thread_local std::vector<LBRangeType> shapeValues;
386 childFE.localBasis().evaluateFunction(local, shapeValues);
387
388 assert(childDOFs.size() == shapeValues.size());
389
390 LIRangeType y(0);
391 for (std::size_t i = 0; i < shapeValues.size(); ++i)
392 y += shapeValues[i] * childDOFs[i];
393
394 fatherDOFsTemp_[currentDOF] = T(y);
395 finishedDOFs_[currentDOF++] = true;
396 return y;
397 }
398 }
399 return fatherDOFsTemp_[currentDOF++];
400 };
401
402 fatherFE.localInterpolation().interpolate(evalLeaf, fatherDOFs);
403
404 // Return true if all father DOFs have been evaluated
405 return std::accumulate(finishedDOFs_.begin(), finishedDOFs_.end(), true,
406 std::logical_and<bool>());
407}
408
409
410template <class N, class C, class B>
411 template <class Trafo>
413prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
414 Trafo const& trafo, bool init)
415{
416 auto& fatherNode = *fatherNode_;
417 if (init)
418 {
419 // TODO(FM): This is UB, see https://gitlab.com/amdis/amdis/-/issues/16 (case 1)
420 bindTree(fatherNode, father);
421 }
422 auto const& childNode = *node_;
423
424 // evaluate father in child coordinate x
425 auto evalFather = [&](LIDomainType const& x) -> LIRangeType
426 {
427 thread_local std::vector<LBRangeType> shapeValues;
428 fatherNode.finiteElement().localBasis().evaluateFunction(trafo(x), shapeValues);
429 assert(shapeValues.size() == fatherDOFs.size());
430
431 LIRangeType y(0);
432 for (std::size_t i = 0; i < shapeValues.size(); ++i)
433 y += shapeValues[i] * fatherDOFs[i];
434
435 return y;
436 };
437
438 auto const& childFE = childNode.finiteElement();
439 thread_local std::vector<T> childDOFs;
440 childFE.localInterpolation().interpolate(evalFather, childDOFs);
441
442 coeff_->scatter(*lv_, childNode, childDOFs, Assigner::assign{});
443}
444
445} // end namespace AMDiS
Definition: InterpolationDataTransfer.hpp:37
void preAdapt(Basis const &basis, Container const &coeff, bool mightCoarsen)
Definition: InterpolationDataTransfer.inc.hpp:49
The restriction of a finite element basis to a single element.
Definition: LocalView.hpp:16
Definition: InterpolationDataTransfer.inc.hpp:255
void copyLocal(NodeElementData const &dofs) const
Copy already existing data to element bound to node_.
Definition: InterpolationDataTransfer.inc.hpp:323
void cacheLocal(NodeElementData &dofs) const
Cache data on the element bound to node_.
Definition: InterpolationDataTransfer.inc.hpp:287
void adaptInit(LocalView const &lv, Container &coeff, Node const &node)
To be called once before copyLocal/prolongLocal are called within the adapt step.
Definition: InterpolationDataTransfer.inc.hpp:313
void preAdaptInit(LocalView const &lv, Container const &coeff, Node const &node)
To be called once before cacheLocal/restrictLocal are called within the preAdapt step.
Definition: InterpolationDataTransfer.inc.hpp:273
constexpr bool LocalView
A Dune::Functions::LocalView type.
Definition: Concepts.hpp:181
Definition: Assigner.hpp:8