10#include <unordered_map>
14#include <dune/common/fvector.hh>
15#include <dune/common/hash.hh>
16#include <dune/common/version.hh>
18#include <dune/grid/common/geometry.hh>
19#include <dune/grid/common/mcmgmapper.hh>
20#include <dune/grid/common/rangegenerators.hh>
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>
34 template <
class LocalCoord>
35 std::size_t operator()(LocalCoord
const& coord)
const
38 for (std::size_t i = 0; i < coord.size(); ++i)
39 Dune::hash_combine(seed, coord[i]);
47template <
class B,
class C>
49preAdapt(B
const& basis, C
const& coeff,
bool mightCoarsen)
51 GridView gv = basis.gridView();
52 LocalView lv = basis.localView();
53 auto const& grid = gv.grid();
54 auto const& idSet = grid.localIdSet();
56 nodeDataTransfer_.resize(lv.tree());
57 Traversal::forEachLeafNode(lv.tree(), [&](
auto const& node,
auto const& tp) {
58 nodeDataTransfer_[tp].preAdaptInit(lv, coeff, node);
62 persistentContainer_.clear();
63 for (
const auto& e : elements(gv))
65 auto it = persistentContainer_.emplace(idSet.id(e),
66 Dune::TypeTree::makeTreeContainer(lv.tree(), [](
auto&& node) { return NodeElementData<TYPEOF(node)>{}; }));
69 auto& treeContainer = it.first->second;
70 Traversal::forEachLeafNode(lv.tree(), [&](
auto const& ,
auto const& tp) {
71 nodeDataTransfer_[tp].cacheLocal(treeContainer[tp]);
79 auto maxLevel = grid.maxLevel();
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);
85 for (
auto const& e : entitySet(basis))
88 while (father.mightVanish() && father.hasFather())
90 father = father.father();
91 auto it = persistentContainer_.emplace(idSet.id(father),
92 Dune::TypeTree::makeTreeContainer(lv.tree(), [](
auto&& node) { return NodeElementData<TYPEOF(node)>{}; }));
96 auto& treeContainer = it.first->second;
98 bool restrictLocalCompleted =
false;
99 auto hItEnd = father.hend(maxLevel);
100 for (
auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt)
102 seeds[hIt->level()] = hIt->seed();
107 auto const& child = *hIt;
108 auto search = persistentContainer_.find(idSet.id(child));
109 assert(search != persistentContainer_.end());
110 auto const& childContainer = search->second;
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>;
119 auto xInChild = [&](LocalCoordinate
const& x) -> BoolCoordPair {
120 LocalCoordinate y = x;
121 for (
int i = father.level()+1; i <= child.level(); ++i)
123 auto currentElement = grid.entity(seeds[i]);
124 auto const& geoInFather = currentElement.geometryInFather();
125 y = geoInFather.local(y);
128 bool isInside = Dune::Geo::Impl::checkInside(Dune::referenceElement(geoInFather).type().
id(), Geometry::mydimension, y, checkInsideTolerance);
130 return BoolCoordPair(
false, std::move(y));
132 return BoolCoordPair(
true, std::move(y));
136 ChildCache childCache;
137 auto xInChildCached = [&](LocalCoordinate
const& x) -> BoolCoordPair {
138 return childCache.get(x, [&](LocalCoordinate
const& x) {
return xInChild(x); });
141 restrictLocalCompleted =
true;
142 Traversal::forEachLeafNode(lv.tree(), [&](
auto const& ,
auto const& tp) {
143 restrictLocalCompleted &=
144 nodeDataTransfer_[tp].restrictLocal(father, treeContainer[tp], xInChildCached,
145 childContainer[tp], init);
150 assert(restrictLocalCompleted);
157template <
class B,
class C>
164 if (persistentContainer_.empty())
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);
174 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
175 Mapper mapper{gv, Dune::mcmgElementLayout()};
177 std::vector<bool> finished(mapper.size(),
false);
178 for (
const auto& e : entitySet(basis))
180 auto index = mapper.index(e);
184 auto it = persistentContainer_.find(idSet.id(e));
187 if (it != persistentContainer_.end()) {
189 auto const& treeContainer = it->second;
190 Traversal::forEachLeafNode(lv.tree(), [&](
auto const& ,
auto const& tp) {
191 nodeDataTransfer_[tp].copyLocal(treeContainer[tp]);
193 finished[index] =
true;
199 while (father.hasFather() && father.isNew())
200 father = father.father();
202 auto maxLevel = gv.grid().maxLevel();
205 auto father_it = persistentContainer_.find(idSet.id(father));
206 assert(father_it != persistentContainer_.end());
207 auto const& treeContainer = father_it->second;
209 auto hItEnd = father.hend(maxLevel);
210 for (
auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt) {
214 auto const& child = *hIt;
218 auto xInFather = [&](LocalCoordinate
const& x) -> LocalCoordinate
221 auto currentElement = child;
222 while (currentElement.level() != father.level())
224 y = currentElement.geometryInFather().global(y);
225 currentElement = currentElement.father();
230 Traversal::forEachLeafNode(lv.tree(), [&](
auto const& ,
auto const& tp) {
231 nodeDataTransfer_[tp].prolongLocal(father, treeContainer[tp], xInFather, init);
234 finished[mapper.index(child)] =
true;
243template <
class B,
class C>
246 persistentContainer_.clear();
253template <
class Node,
class Container,
class Basis>
256 using T =
typename Container::value_type;
258 using Element =
typename Node::Element;
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;
267 using NodeElementData = std::vector<T>;
273 void preAdaptInit(LocalView
const& lv, Container
const& coeff, Node
const& node)
277 fatherNode_ = std::make_shared<Node>(node);
278 constCoeff_ = &coeff;
289 constCoeff_->gather(*lv_, *node_, dofs);
307 template <
class Trafo>
308 bool restrictLocal(Element
const& father, NodeElementData& fatherDOFs, Trafo
const& trafo,
309 NodeElementData
const& childDOFs,
bool init);
313 void adaptInit(LocalView
const& lv, Container& coeff, Node
const& node)
317 fatherNode_ = std::make_shared<Node>(node);
340 template <
class Trafo>
341 void prolongLocal(Element
const& father, NodeElementData
const& fatherDOFs,
342 Trafo
const& trafo,
bool init);
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_;
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)
361 auto& fatherNode = *fatherNode_;
362 std::size_t currentDOF = 0;
366 bindTree(fatherNode, father);
368 auto const& childNode = *node_;
369 auto const& childFE = childNode.finiteElement();
370 auto const& fatherFE = fatherNode.finiteElement();
373 finishedDOFs_.assign(fatherFE.size(),
false);
374 fatherDOFsTemp_.assign(fatherFE.size(), 0);
377 auto evalLeaf = [&](LIDomainType
const& x) -> LIRangeType {
378 if (!finishedDOFs_[currentDOF])
380 auto const& insideLocal = trafo(x);
381 bool isInside = insideLocal.first;
384 auto const& local = insideLocal.second;
385 thread_local std::vector<LBRangeType> shapeValues;
386 childFE.localBasis().evaluateFunction(local, shapeValues);
388 assert(childDOFs.size() == shapeValues.size());
391 for (std::size_t i = 0; i < shapeValues.size(); ++i)
392 y += shapeValues[i] * childDOFs[i];
394 fatherDOFsTemp_[currentDOF] = T(y);
395 finishedDOFs_[currentDOF++] =
true;
399 return fatherDOFsTemp_[currentDOF++];
402 fatherFE.localInterpolation().interpolate(evalLeaf, fatherDOFs);
405 return std::accumulate(finishedDOFs_.begin(), finishedDOFs_.end(),
true,
406 std::logical_and<bool>());
410template <
class N,
class C,
class B>
411 template <
class Trafo>
413prolongLocal(Element
const& father, NodeElementData
const& fatherDOFs,
414 Trafo
const& trafo,
bool init)
416 auto& fatherNode = *fatherNode_;
420 bindTree(fatherNode, father);
422 auto const& childNode = *node_;
425 auto evalFather = [&](LIDomainType
const& x) -> LIRangeType
427 thread_local std::vector<LBRangeType> shapeValues;
428 fatherNode.finiteElement().localBasis().evaluateFunction(trafo(x), shapeValues);
429 assert(shapeValues.size() == fatherDOFs.size());
432 for (std::size_t i = 0; i < shapeValues.size(); ++i)
433 y += shapeValues[i] * fatherDOFs[i];
438 auto const& childFE = childNode.finiteElement();
439 thread_local std::vector<T> childDOFs;
440 childFE.localInterpolation().interpolate(evalFather, childDOFs);
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