#pragma once #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace AMDiS { namespace Impl { // Hash function for cache container struct CoordHasher { template std::size_t operator()(LocalCoord const& coord) const { std::size_t seed = 0; for (std::size_t i = 0; i < coord.size(); ++i) Dune::hash_combine(seed, coord[i]); return seed; } }; } // end namespace Impl template DataTransfer::DataTransfer(std::shared_ptr basis) : basis_(std::move(basis)) , mapper_(basis_->gridView().grid(), Dune::mcmgElementLayout()) , nodeDataTransfer_() {} template void DataTransfer:: preAdapt(C const& coeff, bool mightCoarsen) { GridView gv = basis_->gridView(); LocalView lv = basis_->localView(); auto const& idSet = gv.grid().localIdSet(); for_each_leaf_node(lv.tree(), [&](auto const& node, auto const& tp) { nodeDataTransfer_[tp].preAdaptInit(lv, coeff, node); }); // Make persistent DoF container persistentContainer_.clear(); // Redundant if postAdapt was correctly called last cycle for (const auto& e : elements(gv)) { auto it = persistentContainer_.emplace(idSet.id(e), TypeTree::treeContainer(lv.tree())); lv.bind(e); auto& treeContainer = it.first->second; for_each_leaf_node(lv.tree(), [&](auto const& node, auto const& tp) { nodeDataTransfer_[tp].cacheLocal(treeContainer[tp]); }); } if (!mightCoarsen) return; // Interpolate from possibly vanishing elements auto maxLevel = gv.grid().maxLevel(); using std::sqrt; typename Grid::ctype const checkInsideTolerance = sqrt(std::numeric_limits::epsilon()); for (auto const& e : elements(gv, typename C::Traits::PartitionSet{})) { auto father = e; while (father.mightVanish() && father.hasFather()) { father = father.father(); auto it = persistentContainer_.emplace(idSet.id(father), TypeTree::treeContainer(lv.tree())); if (!it.second) continue; auto& treeContainer = it.first->second; auto geo = father.geometry(); bool init = true; // init flag for first call on new father element bool restrictLocalCompleted = false; auto hItEnd = father.hend(maxLevel); for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt) { if (!hIt->isLeaf()) continue; auto const& child = *hIt; auto search = persistentContainer_.find(idSet.id(child)); assert(search != persistentContainer_.end()); auto const& childContainer = search->second; lv.bind(child); auto const& childGeo = child.geometry(); auto refTypeId = Dune::referenceElement(childGeo).type().id(); using BoolCoordPair = std::pair; using CacheImp = std::unordered_map; using ChildCache = ConcurrentCache; // Transfers input father-local point x into child-local point y // Returns false if x is not inside the child auto xInChild = [&](LocalCoordinate const& x) -> BoolCoordPair { LocalCoordinate local = childGeo.local(geo.global(x)); // TODO(FM): Using an implementation detail as workaround for insufficient // tolerance, see https://gitlab.dune-project.org/core/dune-grid/issues/84 bool isInside = Dune::Geo::Impl::checkInside(refTypeId, Geometry::coorddimension, local, checkInsideTolerance); return BoolCoordPair(isInside, std::move(local)); }; // TODO(FM): Disable for single-node basis ChildCache childCache; auto xInChildCached = [&](LocalCoordinate const& x) -> BoolCoordPair { return childCache.get(x, [&](LocalCoordinate const& x) { return xInChild(x); }); }; restrictLocalCompleted = true; for_each_leaf_node(lv.tree(), [&](auto const& node, auto const& tp) { restrictLocalCompleted &= nodeDataTransfer_[tp].restrictLocal(father, treeContainer[tp], xInChildCached, childContainer[tp], init); }); init = false; } // test if restrictLocal was completed on all nodes assert(restrictLocalCompleted); } // end while (father.mightVanish) } // end for (elements) } template void DataTransfer::adapt(C& coeff) { coeff.resize(); // No data was saved before adapting the grid, make // sure to call DataTransfer::preAdapt before calling adapt() on the grid if (persistentContainer_.empty()) return; GridView gv = basis_->gridView(); LocalView lv = basis_->localView(); auto const& idSet = gv.grid().localIdSet(); for_each_leaf_node(lv.tree(), [&](auto const& node, auto const& tp) { nodeDataTransfer_[tp].adaptInit(lv, coeff, node); }); mapper_.update(); std::vector finished(mapper_.size(), false); for (const auto& e : elements(gv, typename C::Traits::PartitionSet{})) { auto index = mapper_.index(e); if (finished[index]) continue; auto it = persistentContainer_.find(idSet.id(e)); // Data already exists and no interpolation is required if (it != persistentContainer_.end()) { lv.bind(e); auto const& treeContainer = it->second; for_each_leaf_node(lv.tree(), [&](auto const& node, auto const& tp) { nodeDataTransfer_[tp].copyLocal(treeContainer[tp]); }); finished[index] = true; continue; } // Data needs to be interpolated auto father = e; while (father.hasFather() && father.isNew()) father = father.father(); auto maxLevel = gv.grid().maxLevel(); auto fatherGeo = father.geometry(); bool init = true; // init flag for first call on new father element auto father_it = persistentContainer_.find(idSet.id(father)); assert(father_it != persistentContainer_.end()); auto const& treeContainer = father_it->second; auto hItEnd = father.hend(maxLevel); for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt) { if (!hIt->isLeaf()) continue; auto const& child = *hIt; lv.bind(child); // coordinate transform from child to father element auto xInFather = [&fatherGeo, childGeo = child.geometry()] (LocalCoordinate const& x) -> LocalCoordinate { return fatherGeo.local(childGeo.global(x)); }; for_each_leaf_node(lv.tree(), [&](auto const& node, auto const& tp) { nodeDataTransfer_[tp].prolongLocal(father, treeContainer[tp], xInFather, init); }); finished[mapper_.index(child)] = true; init = false; } } // end for (elements) coeff.finish(); } template void DataTransfer::postAdapt(C&) { persistentContainer_.clear(); } /** Element-local data transfer on a single leaf node of the basis tree * Handles computations related to the finite element basis node */ template class NodeDataTransfer { using T = typename Container::value_type; using LocalView = typename Basis::LocalView; using Element = typename Node::Element; using LocalBasis = typename Node::FiniteElement::Traits::LocalBasisType; using LBRangeType = typename LocalBasis::Traits::RangeType; using LocalInterpolation = typename Node::FiniteElement::Traits::LocalBasisType; using LIDomainType = typename LocalInterpolation::Traits::DomainType; using LIRangeType = typename LocalInterpolation::Traits::RangeType; public: using NodeElementData = std::vector; public: NodeDataTransfer() = default; /// To be called once before cacheLocal/restrictLocal are called within the preAdapt step void preAdaptInit(LocalView const& lv, Container const& coeff, Node const& node) { lv_ = &lv; node_ = &node; fatherNode_ = std::make_unique(node); constCoeff_ = &coeff; } /// \brief Cache data on the element bound to node_ /** * This functions is used whenever the element does not vanish and thus the * data can trivially be transferred to the new element **/ // [[expects: preAdaptInit to be called before]] void cacheLocal(NodeElementData& dofs) const { constCoeff_->gather(*lv_, *node_, dofs); } /** \brief Evaluate data on the child element bound to node_ and interpolate onto * father entity using the coordinate transformation trafo from father to child. * * Stores cached data in the NodeElementData argument. After grid adaption the * data is copied by \ref copyLocal or \ref prolongLocal to the target element * in the new grid. * * \param father The father element to interpolate to * \param fatherDOFs Container to store the interpolated DOFs * \param trafo Coordinate transform from local coordinates in father to local * coordinates in child element * \param childDOFs DOF values from the child element * \param init The father element is visited for the first time **/ // [[expects: preAdaptInit to be called before]] template bool restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& trafo, NodeElementData const& childDOFs, bool init); /// To be called once before copyLocal/prolongLocal are called within the adapt step void adaptInit(LocalView const& lv, Container& coeff, Node const& node) { lv_ = &lv; node_ = &node; fatherNode_ = std::make_unique(node); coeff_ = &coeff; } /// \brief Copy already existing data to element bound to node_ // [[expects: adaptInit to be called before]] void copyLocal(NodeElementData const& dofs) const { coeff_->scatter(*lv_, *node_, dofs, Assigner::assign{}); } /** \brief Interpolate data from father onto the child element bound to node_ using * the transformation trafo from child to father * * Stores the interpolated data from father to child in the container \ref coeff_. * * \param father The father element * \param fatherDOFs DOF values cached on the father element before adapt * \param trafo Coordinate transform from local coordinates in child to local * coordinates in father element * \param init Father element is visited for the first time **/ // [[expects: adaptInit to be called before]] template void prolongLocal(Element const& father, NodeElementData const& fatherDOFs, Trafo const& trafo, bool init); private: LocalView const* lv_ = nullptr; Node const* node_ = nullptr; std::unique_ptr fatherNode_; Container const* constCoeff_ = nullptr; Container* coeff_ = nullptr; std::vector finishedDOFs_; NodeElementData fatherDOFsTemp_; }; template template bool NodeDataTransfer:: restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& trafo, NodeElementData const& childDOFs, bool init) { auto& fatherNode = *fatherNode_; std::size_t currentDOF = 0; if (init) { // TODO(FM): This is UB, replace with FE cache for father bindTree(fatherNode, father); } auto const& childNode = *node_; auto const& childFE = childNode.finiteElement(); auto const& fatherFE = fatherNode.finiteElement(); if (init) { finishedDOFs_.assign(fatherFE.size(), false); fatherDOFsTemp_.assign(fatherFE.size(), 0); } auto evalLeaf = [&](LIDomainType const& x) -> LIRangeType { if (!finishedDOFs_[currentDOF]) { auto const& insideLocal = trafo(x); bool isInside = insideLocal.first; if (isInside) { auto const& local = insideLocal.second; thread_local std::vector shapeValues; childFE.localBasis().evaluateFunction(local, shapeValues); assert(childDOFs.size() == shapeValues.size()); LIRangeType y(0); for (std::size_t i = 0; i < shapeValues.size(); ++i) y += shapeValues[i] * childDOFs[i]; fatherDOFsTemp_[currentDOF] = T(y); finishedDOFs_[currentDOF++] = true; return y; } } return fatherDOFsTemp_[currentDOF++]; }; auto evalLeafFct = functionFromCallable(evalLeaf); fatherFE.localInterpolation().interpolate(evalLeafFct, fatherDOFs); // Return true if all father DOFs have been evaluated return std::accumulate(finishedDOFs_.begin(), finishedDOFs_.end(), true, std::logical_and()); } template template void NodeDataTransfer:: prolongLocal(Element const& father, NodeElementData const& fatherDOFs, Trafo const& trafo, bool init) { auto& fatherNode = *fatherNode_; if (init) { // TODO(FM): This is UB, replace with FE cache for father bindTree(fatherNode, father); } auto const& childNode = *node_; // evaluate father in child coordinate x auto evalFather = [&](LIDomainType const& x) -> LIRangeType { thread_local std::vector shapeValues; fatherNode.finiteElement().localBasis().evaluateFunction(trafo(x), shapeValues); assert(shapeValues.size() == fatherDOFs.size()); LIRangeType y(0); for (std::size_t i = 0; i < shapeValues.size(); ++i) y += shapeValues[i] * fatherDOFs[i]; return y; }; auto const& childFE = childNode.finiteElement(); thread_local std::vector childDOFs; auto evalFatherFct = functionFromCallable(evalFather); childFE.localInterpolation().interpolate(evalFatherFct, childDOFs); coeff_->scatter(*lv_, childNode, childDOFs, Assigner::assign{}); } } // end namespace AMDiS