diff --git a/src/amdis/DataTransfer.hpp b/src/amdis/DataTransfer.hpp
index 362bca0d4e11051c33e7e8dfe2e5bdd8aa80e45d..e7da6447add188f909540e06968739658ab99fcf 100644
--- a/src/amdis/DataTransfer.hpp
+++ b/src/amdis/DataTransfer.hpp
@@ -1,8 +1,12 @@
#pragma once
+#include <map>
#include <memory>
+#include <dune/grid/common/mcmgmapper.hh>
+
#include <amdis/Output.hpp>
+#include <amdis/typetree/TreeContainer.hpp>
namespace AMDiS
{
@@ -51,8 +55,64 @@ namespace AMDiS
};
+ template <class Node, class Container, class Basis>
+ class NodeDataTransfer;
+
+
+ /** Data Transfer implementation for a single grid using interpolation
+ * Handles computations related to the geometric information of the grid and passes that to the
+ * underlying NodeDataTransfer classes
+ */
template <class Container, class Basis>
- class DataTransfer;
+ class DataTransfer
+ : public DataTransferInterface<Container>
+ {
+ using LocalView = typename Basis::LocalView;
+ using Tree = typename LocalView::Tree;
+ using GridView = typename Basis::GridView;
+ using Grid = typename GridView::Grid;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Geometry = typename Element::Geometry;
+ using LocalCoordinate = typename Geometry::LocalCoordinate;
+ using IdType = typename Grid::LocalIdSet::IdType;
+
+ template <class Node>
+ using NodeElementData = typename NodeDataTransfer<Node, Container, Basis>::NodeElementData;
+ using ElementData = TYPEOF(makeTreeContainer<NodeElementData>(std::declval<const Tree&>()));
+
+ public:
+ DataTransfer(std::shared_ptr<Basis> basis);
+
+ /** Saves data contained in coeff in the PersistentContainer
+ * To be called after grid.preAdapt() and before grid.adapt()
+ */
+ void preAdapt(Container const& coeff, bool mightCoarsen) override;
+
+ /** Unpacks data from the PersistentContainer
+ * To be called after grid.adapt() and before grid.postAdapt()
+ */
+ void postAdapt(Container& coeff, bool refined) override;
+
+ private:
+ /// The global basis
+ std::shared_ptr<Basis const> basis_;
+
+ /// Container with data that persists during grid adaptation
+ using PersistentContainer = std::map<IdType, ElementData>;
+ PersistentContainer persistentContainer_;
+
+ /// Map leaf entities to unique index
+ using Mapper = Dune::LeafMultipleCodimMultipleGeomTypeMapper<Grid>;
+ Mapper mapper_;
+
+ /// Data transfer on a single basis node
+ template <class Node>
+ using NDT = NodeDataTransfer<Node, Container, Basis>;
+ using NodeDataTransferContainer = TYPEOF(makeTreeContainer<NDT>(std::declval<const Tree&>()));
+ NodeDataTransferContainer nodeDataTransfer_;
+ };
+
/// Factory to create DataTransfer objects based on the \ref DataTransferOperation
template <class Container>
@@ -62,14 +122,14 @@ namespace AMDiS
public:
template <class Basis>
- static std::unique_ptr<Interface> create(DataTransferOperation op, Basis const& basis)
+ static std::unique_ptr<Interface> create(DataTransferOperation op, std::shared_ptr<Basis> basis)
{
switch (op)
{
case DataTransferOperation::NO_OPERATION:
return std::make_unique<NoDataTransfer<Container>>();
case DataTransferOperation::INTERPOLATE:
- return std::make_unique<DataTransfer<Container, Basis>>(basis);
+ return std::make_unique<DataTransfer<Container, Basis>>(std::move(basis));
default:
error_exit("Invalid data transfer\n");
return nullptr; // avoid warnings
diff --git a/src/amdis/DataTransfer.inc.hpp b/src/amdis/DataTransfer.inc.hpp
index ca5d200fe20737810cb41e2150720bac3b7d4c19..c5a8a1775c46e14346555bb470cd609d5fff72d3 100644
--- a/src/amdis/DataTransfer.inc.hpp
+++ b/src/amdis/DataTransfer.inc.hpp
@@ -15,454 +15,410 @@
#include <dune/common/hash.hh>
#include <dune/grid/common/geometry.hh>
-#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/functions/common/functionfromcallable.hh>
-#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <amdis/Output.hpp>
#include <amdis/common/ConcurrentCache.hpp>
#include <amdis/typetree/Traversal.hpp>
-#include <amdis/typetree/TreeContainer.hpp>
-
-namespace AMDiS
-{
- template <class Node, class Container, class Basis>
- class NodeDataTransfer;
+namespace AMDiS {
- /** Data Transfer implementation for a single grid
- * Handles computations related to the geometric information of the grid and passes that to the
- * underlying NodeDataTransfer classes
- */
- template <class Container, class Basis>
- class DataTransfer
- : public DataTransferInterface<Container>
+namespace Impl
+{
+ // Hash function for cache container
+ struct CoordHasher
{
- using LocalView = typename Basis::LocalView;
- using Tree = typename LocalView::Tree;
- using GridView = typename Basis::GridView;
- using Grid = typename GridView::Grid;
- using Mapper = Dune::LeafMultipleCodimMultipleGeomTypeMapper<Grid>;
-
- using IdSet = typename Grid::LocalIdSet;
- using IdType = typename IdSet::IdType;
-
- using Element = typename GridView::template Codim<0>::Entity;
- using Geometry = typename Element::Geometry;
- using LocalCoordinate = typename Geometry::LocalCoordinate;
- using BoolCoordPair = std::pair<bool, LocalCoordinate>;
-
- // Hash function for cache container
- struct Hasher
+ template <class LocalCoord>
+ std::size_t operator()(LocalCoord const& coord) const
{
- std::size_t operator()(LocalCoordinate 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;
- }
- };
+ std::size_t seed = 0;
+ for (std::size_t i = 0; i < coord.size(); ++i)
+ Dune::hash_combine(seed, coord[i]);
+ return seed;
+ }
+ };
- using CacheImp = std::unordered_map<LocalCoordinate, BoolCoordPair, Hasher>;
- using ChildCache = ConcurrentCache<LocalCoordinate, BoolCoordPair, ConsecutivePolicy, CacheImp>;
+} // end namespace Impl
- constexpr static int d = Geometry::coorddimension;
- using ctype = typename Geometry::ctype;
- // Returns the Node's NodeDataTransfer
- struct NDT
- {
- template<class Node>
- auto operator()(const Node& node)
- {
- using NDT = NodeDataTransfer<Node, Container, Basis>;
- return NDT{};
- }
- };
+template <class C, class B>
+DataTransfer<C,B>::DataTransfer(std::shared_ptr<B> basis)
+ : basis_(std::move(basis))
+ , mapper_(basis_->gridView().grid(), Dune::mcmgElementLayout())
+ , nodeDataTransfer_(makeTreeContainer<NDT>(basis_->localView().tree()))
+{}
- using NodeDataTransferContainer = TYPEOF(makeTreeContainer<Tree, NDT>(std::declval<const Tree&>(), NDT()));
- // Returns the Node's NodeElementData
- struct NodeElementData
- {
- template<class Node>
- auto operator()(const Node& node)
- {
- using NED = typename NodeDataTransfer<Node, Container, Basis>
- ::NodeElementData;
- return NED{};
- }
- };
-
- using ElementData = TYPEOF(makeTreeContainer<Tree, NodeElementData>(std::declval<const Tree&>(), NodeElementData()));
-
- public:
- // Container with data that persists during grid adaptation
- using PersistentContainer = std::map<IdType, ElementData>;
-
- public:
- DataTransfer(Basis const& basis)
- : basis_(&basis)
- , mapper_(basis.gridView().grid(), Dune::mcmgElementLayout())
- , finished_(mapper_.size(), false)
- , nodeDataTransfer_(makeTreeContainer<Tree, NDT>(basis_->localView().tree(), NDT()))
- {}
-
- /** Saves data contained in coeff in the PersistentContainer
- * To be called after grid.preAdapt() and before grid.adapt()
- */
- void preAdapt(Container const& coeff, bool mightCoarsen) override;
-
- /** Unpacks data from the PersistentContainer
- * To be called after grid.adapt() and before grid.postAdapt()
- */
- void postAdapt(Container& coeff, bool refined) override;
-
- private:
- Basis const* basis_;
- PersistentContainer persistentContainer_;
- Mapper mapper_;
- std::vector<bool> finished_;
- NodeDataTransferContainer nodeDataTransfer_;
- };
+template <class C, class B>
+void DataTransfer<C,B>::preAdapt(C const& coeff, bool mightCoarsen)
+{
+ GridView gv = basis_->gridView();
+ LocalView lv = basis_->localView();
+ auto const& idSet = gv.grid().localIdSet();
- template <class C, class B>
- void DataTransfer<C,B>::preAdapt(C const& coeff, bool mightCoarsen)
+ for_each_leaf_node(lv.tree(), [&](auto const& node, auto const& tp) {
+ nodeDataTransfer_[tp].preAdaptInit(lv, coeff, node);
+ });
+
+ // Make persistent DoF container
+ persistentContainer_.clear();
+ for (const auto& e : elements(gv))
{
- auto gv = basis_->gridView();
- auto lv = basis_->localView();
- auto const& idSet = gv.grid().localIdSet();
+ auto it = persistentContainer_.emplace(idSet.id(e), makeTreeContainer<NodeElementData>(lv.tree()));
+ lv.bind(e);
+ auto& treeContainer = it.first->second;
for_each_leaf_node(lv.tree(), [&](auto const& node, auto const& tp) {
- nodeDataTransfer_[tp].preAdaptInit(lv, coeff, node);
+ nodeDataTransfer_[tp].cacheLocal(treeContainer[tp]);
});
+ }
+
+ if (!mightCoarsen)
+ return;
- // Make persistent DoF container
- persistentContainer_.clear();
- for (const auto& e : elements(gv))
+ // Interpolate from possibly vanishing elements
+ auto maxLevel = gv.grid().maxLevel();
+ using std::sqrt;
+ typename Grid::ctype const checkInsideTolerance = sqrt(std::numeric_limits<typename Grid::ctype>::epsilon());
+ for (auto const& e : elements(gv))
+ {
+ auto father = e;
+ while (father.mightVanish() && father.hasFather())
{
- auto it = persistentContainer_.emplace(idSet.id(e),
- makeTreeContainer<Tree, NodeElementData>(lv.tree(), NodeElementData()));
+ father = father.father();
+ auto it = persistentContainer_.emplace(idSet.id(father), makeTreeContainer<NodeElementData>(lv.tree()));
+ if (!it.second)
+ continue;
- lv.bind(e);
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<bool, LocalCoordinate>;
+ using CacheImp = std::unordered_map<LocalCoordinate, BoolCoordPair, Impl::CoordHasher>;
+ using ChildCache = ConcurrentCache<LocalCoordinate, BoolCoordPair, ConsecutivePolicy, CacheImp>;
+
+ // 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 <class C, class B>
+void DataTransfer<C,B>::postAdapt(C& coeff, bool refined)
+{
+ coeff.resize(*basis_);
+ 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].postAdaptInit(lv, coeff, node);
+ });
+
+ mapper_.update();
+ std::vector<bool> finished(mapper_.size(), false);
+ for (const auto& e : elements(gv))
+ {
+ 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].cacheLocal(treeContainer[tp]);
+ nodeDataTransfer_[tp].copyLocal(treeContainer[tp]);
});
+ finished[index] = true;
+ continue;
}
- // Interpolate from possibly vanishing elements
- if (mightCoarsen) {
- auto maxLevel = gv.grid().maxLevel();
- using std::sqrt;
- ctype const checkInsideTolerance = sqrt(std::numeric_limits<ctype>::epsilon());
- for (auto const& e : elements(gv))
- {
- auto father = e;
- while (father.mightVanish() && father.hasFather())
- {
- father = father.father();
- auto it = persistentContainer_.emplace(idSet.id(father),
- makeTreeContainer<Tree, NodeElementData>(lv.tree(), NodeElementData()));
- if (it.second) {
- 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()) {
- 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 const& ref = Dune::referenceElement(childGeo);
- auto const& refTypeId = ref.type().id();
-
- // 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, d, 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);
- }
- }
- }
- }
- }
+ // Data needs to be interpolated
+ auto father = e;
+ while (father.hasFather() && father.isNew())
+ father = father.father();
- template <class C, class B>
- void DataTransfer<C,B>::postAdapt(C& coeff, bool refined)
- {
- coeff.resize(*basis_);
- auto gv = basis_->gridView();
- auto lv = basis_->localView();
- auto const& idSet = gv.grid().localIdSet();
- for_each_leaf_node(lv.tree(), [&](auto const& node, auto const& tp) {
- nodeDataTransfer_[tp].postAdaptInit(lv, coeff, node);
- });
+ auto maxLevel = gv.grid().maxLevel();
+ auto fatherGeo = father.geometry();
+ bool init = true; // init flag for first call on new father element
- mapper_.update();
- finished_.clear();
- finished_.resize(mapper_.size(), false);
- for (const auto& e : elements(gv))
- {
- auto index = mapper_.index(e);
- if (finished_[index])
+ 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 it = persistentContainer_.find(idSet.id(e));
+ auto const& child = *hIt;
+ lv.bind(child);
- // 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;
- }
- // Data needs to be interpolated
- else {
- 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 it = persistentContainer_.find(idSet.id(father));
- assert(it != persistentContainer_.end());
- auto const& treeContainer = it->second;
-
- auto hItEnd = father.hend(maxLevel);
- for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt) {
- if (hIt->isLeaf()) {
- 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;
- }
- }
- }
+ // 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)
+}
+
+
+/** Element-local data transfer on a single leaf node of the basis tree
+ * Handles computations related to the finite element basis node
+ */
+template <class Node, class Container, class Basis>
+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<T>;
+
+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>(node);
+ constCoeff_ = &coeff;
}
- /** Element-local data transfer on a single leaf node of the basis tree
- * Handles computations related to the finite element basis node
- */
- template <class Node, class Container, class Basis>
- class NodeDataTransfer
+ /// \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
{
- using T = typename Container::value_type;
- using LocalView = typename Basis::LocalView;
- using Element = typename Node::Element;
+ auto const& fe = node_->finiteElement();
+ auto feSize = fe.size();
+ dofs.resize(feSize);
+ for (std::size_t i = 0; i < feSize; ++i)
+ dofs[i] = (*constCoeff_)[lv_->index(node_->localIndex(i))];
+ }
- 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;
+ /** \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 <class Trafo>
+ 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 postAdapt step
+ void postAdaptInit(LocalView const& lv, Container& coeff, Node const& node)
+ {
+ lv_ = &lv;
+ node_ = &node;
+ fatherNode_ = std::make_unique<Node>(node);
+ coeff_ = &coeff;
+ }
- public:
- using NodeElementData = std::vector<T>;
+ /// \brief Copy already existing data to element bound to node_
+ // [[expects: postAdaptInit to be called before]]
+ void copyLocal(NodeElementData const& dofs) const
+ {
+ for (std::size_t i = 0; i < dofs.size(); ++i)
+ (*coeff_)[lv_->index(node_->localIndex(i))] = dofs[i];
+ }
- public:
- NodeDataTransfer() = default;
+ /** \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: postAdaptInit to be called before]]
+ template <class Trafo>
+ 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<Node> fatherNode_;
+ Container const* constCoeff_ = nullptr;
+ Container* coeff_ = nullptr;
+ std::vector<bool> finishedDOFs_;
+ NodeElementData fatherDOFsTemp_;
+};
+
+
+template <class N, class C, class B>
+ template <class Trafo>
+bool NodeDataTransfer<N,C,B>::
+ 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();
- /// 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;
- auto nodeCopy = node;
- fatherNode_ = std::make_unique<Node>(std::move(nodeCopy));
- constCoeff_ = &coeff;
- }
+ if (init) {
+ finishedDOFs_.assign(fatherFE.size(), false);
+ fatherDOFsTemp_.assign(fatherFE.size(), 0);
+ }
- /// To be called once before copyLocal/prolongLocal are called within the postAdapt step
- void postAdaptInit(LocalView const& lv, Container& coeff, Node const& node)
+ auto evalLeaf = [&](LIDomainType const& x) -> LIRangeType {
+ if (!finishedDOFs_[currentDOF])
{
- lv_ = &lv;
- node_ = &node;
- auto nodeCopy = node;
- fatherNode_ = std::make_unique<Node>(std::move(nodeCopy));
- coeff_ = &coeff;
- }
+ auto const& insideLocal = trafo(x);
+ bool isInside = insideLocal.first;
+ if (isInside)
+ {
+ auto const& local = insideLocal.second;
+ thread_local std::vector<LBRangeType> shapeValues;
+ childFE.localBasis().evaluateFunction(local, shapeValues);
- /// Cache data on the element bound to node_
- void cacheLocal(NodeElementData& dofs) const
- {
- auto const& fe = node_->finiteElement();
- auto feSize = fe.size();
- dofs.resize(feSize);
- for (std::size_t i = 0; i < feSize; ++i)
- dofs[i] = (*constCoeff_)[lv_->index(node_->localIndex(i))];
- }
+ assert(childDOFs.size() == shapeValues.size());
- /** Evaluate data on the child element bound to node_ and interpolate onto father entity
- * using the coordinate transformation trafo from father to child
- */
- template <class Trafo>
- bool restrictLocal(Element const& father, NodeElementData& fatherDOFs, Trafo const& trafo,
- NodeElementData const& childDOFs, bool init);
+ LIRangeType y(0);
+ for (std::size_t i = 0; i < shapeValues.size(); ++i)
+ y += shapeValues[i] * childDOFs[i];
- /// Copy already existing data to element bound to node_
- void copyLocal(NodeElementData const& dofs) const
- {
- for (std::size_t i = 0; i < dofs.size(); ++i)
- (*coeff_)[lv_->index(node_->localIndex(i))] = dofs[i];
+ fatherDOFsTemp_[currentDOF] = y;
+ finishedDOFs_[currentDOF++] = true;
+ return y;
+ }
}
-
- /** Interpolate data from father onto the child element bound to node_ using the transformation
- * trafo from child to father
- */
- template <class Trafo>
- void prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
- Trafo const& trafo, bool init);
-
- private:
- LocalView const* lv_;
- Node const* node_;
- std::unique_ptr<Node> fatherNode_;
- Container const* constCoeff_;
- Container* coeff_;
- std::vector<bool> finishedDOFs_;
- NodeElementData fatherDOFsTemp_;
+ return fatherDOFsTemp_[currentDOF++];
};
- template <class Node, class C, class B>
- template <class Trafo>
- bool NodeDataTransfer<Node,C,B>::
- 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);
- }
+ using FFC
+ = Dune::Functions::FunctionFromCallable<LIRangeType(LIDomainType), decltype(evalLeaf)>;
+ fatherFE.localInterpolation().interpolate(FFC(evalLeaf), fatherDOFs);
- 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<LBRangeType> 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] = y;
- finishedDOFs_[currentDOF++] = true;
- return y;
- }
- }
- return fatherDOFsTemp_[currentDOF++];
- };
+ // Return true if all father DOFs have been evaluated
+ return std::accumulate(finishedDOFs_.begin(), finishedDOFs_.end(), true,
+ std::logical_and<bool>());
+}
- using FFC
- = Dune::Functions::FunctionFromCallable<LIRangeType(LIDomainType), decltype(evalLeaf)>;
- fatherFE.localInterpolation().interpolate(FFC(evalLeaf), fatherDOFs);
- // Return true if all father DOFs have been evaluated
- return std::accumulate(finishedDOFs_.begin(), finishedDOFs_.end(), true,
- std::logical_and<bool>());
+template <class N, class C, class B>
+ template <class Trafo>
+void NodeDataTransfer<N,C,B>::
+ 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_;
- template <class Node, class C, class B>
- template <class Trafo>
- void NodeDataTransfer<Node,C,B>::
- prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
- Trafo const& trafo, bool init)
+ // evaluate father in child coordinate x
+ auto evalFather = [&](LIDomainType const& x) -> LIRangeType
{
- 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<LBRangeType> shapeValues;
- fatherNode.finiteElement().localBasis().evaluateFunction(trafo(x), shapeValues);
- assert(shapeValues.size() == fatherDOFs.size());
+ thread_local std::vector<LBRangeType> 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];
+ LIRangeType y(0);
+ for (std::size_t i = 0; i < shapeValues.size(); ++i)
+ y += shapeValues[i] * fatherDOFs[i];
- return y;
- };
+ return y;
+ };
- auto const& childFE = childNode.finiteElement();
- thread_local std::vector<T> childDOFs;
- using FFC = Dune::Functions
- ::FunctionFromCallable<LIRangeType(LIDomainType), decltype(evalFather)>;
- childFE.localInterpolation().interpolate(FFC(evalFather), childDOFs);
+ auto const& childFE = childNode.finiteElement();
+ thread_local std::vector<T> childDOFs;
+ using FFC = Dune::Functions
+ ::FunctionFromCallable<LIRangeType(LIDomainType), decltype(evalFather)>;
+ childFE.localInterpolation().interpolate(FFC(evalFather), childDOFs);
- for (std::size_t i = 0; i < childDOFs.size(); ++i)
- (*coeff_)[lv_->index(childNode.localIndex(i))] = childDOFs[i];
- }
+ for (std::size_t i = 0; i < childDOFs.size(); ++i)
+ (*coeff_)[lv_->index(childNode.localIndex(i))] = childDOFs[i];
+}
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/DOFVectorBase.hpp b/src/amdis/linearalgebra/DOFVectorBase.hpp
index c505124d17842927b239ee8e351e7dcb64d2b0e1..50f6d6a71ebb9f4f55fbe9e378fffabea0c33528 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.hpp
@@ -84,8 +84,8 @@ namespace AMDiS
public:
/// Constructor. Stores the shared_ptr of the basis and creates a new DataTransfer.
DOFVectorBase(std::shared_ptr<GlobalBasis> basis, DataTransferOperation op)
- : basis_(basis)
- , dataTransfer_(DataTransferFactory::create(op, *basis_))
+ : basis_(std::move(basis))
+ , dataTransfer_(DataTransferFactory::create(op, basis_))
{
compress();
attachToGridTransfer();
diff --git a/src/amdis/typetree/TreeContainer.hpp b/src/amdis/typetree/TreeContainer.hpp
index 4ecc70865e73012b5e9d36b36679bc83863e32d9..62be88589848d7c595f39408ec7182c53a906f8b 100644
--- a/src/amdis/typetree/TreeContainer.hpp
+++ b/src/amdis/typetree/TreeContainer.hpp
@@ -202,6 +202,13 @@ namespace AMDiS
return makeTreeContainer(tree, [](const auto&) {return Value{};});
}
+ template<template<class> class NodeData, class Tree>
+ auto makeTreeContainer(const Tree& tree)
+ {
+ return makeTreeContainer(tree, [](const auto& node) {return NodeData<TYPEOF(node)>{};});
+ }
+
+
/**
* \brief Alias to container type generated by makeTreeContainer for given value and tree type
*/