diff --git a/src/amdis/DataTransfer.hpp b/src/amdis/DataTransfer.hpp
index 539628cc4de5eb52c0f1ab856148aed969c025b9..3f9f8ca6df4ece3ea9e365157d72169e9da85f27 100644
--- a/src/amdis/DataTransfer.hpp
+++ b/src/amdis/DataTransfer.hpp
@@ -1,14 +1,8 @@
#pragma once
#include <memory>
-#include <type_traits>
-#include <utility>
-
-#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <amdis/Output.hpp>
-#include <amdis/utility/TreeData.hpp>
-#include <amdis/utility/Visitor.hpp>
namespace AMDiS
{
@@ -17,17 +11,6 @@ namespace AMDiS
INTERPOLATE = 1
} DataTransferOperation;
- template <class Node, class RangeType>
- class NodeDataTransfer
- {
- public:
- template <class Basis, class Container>
- void preAdapt(Basis const& basis, Container const& coeff, bool mightCoarsen) {}
-
- template <class Basis, class Container>
- void postAdapt(Basis const& basis, Container& coeff, bool refined) const {}
- };
-
template <class Container>
class DataTransferInterface
@@ -40,7 +23,7 @@ namespace AMDiS
virtual void preAdapt(Container const& container, bool mightCoarsen) = 0;
/// Interpolate data to new grid after grid adaption
- virtual void postAdapt(Container& container, bool refined) const = 0;
+ virtual void postAdapt(Container& container, bool refined) = 0;
};
@@ -51,9 +34,9 @@ namespace AMDiS
: public DataTransferInterface<Container>
{
public:
- virtual void preAdapt(Container const& container, bool) override {}
+ void preAdapt(Container const& container, bool) override {}
- virtual void postAdapt(Container& container, bool) const override
+ void postAdapt(Container& container, bool) override
{
container.compress();
}
@@ -61,45 +44,7 @@ namespace AMDiS
template <class Container, class Basis>
- class DataTransfer
- : public DataTransferInterface<Container>
- {
- template <class Node>
- using NDT = NodeDataTransfer<std::decay_t<Node>, typename Container::value_type>;
-
- public:
- DataTransfer(Basis const& basis)
- : basis_(&basis)
- {}
-
- /// Calls \ref NodeDataTransfer::preAdapt() on each basis node
- virtual void preAdapt(Container const& container, bool mightCoarsen) override
- {
- nodeDataTransfer_.init(*basis_);
- AMDiS::forEachLeafNode_(basis_->localView().tree(), [&](auto const& node, auto const& treePath)
- {
- auto subBasis = Dune::Functions::subspaceBasis(*basis_, treePath);
- nodeDataTransfer_[node].preAdapt(subBasis, container.vector(), mightCoarsen);
- });
- }
-
- /// Calls \ref NodeDataTransfer::postAdapt() on each basis node after compressing the
- /// Container the dimension of the basis
- virtual void postAdapt(Container& container, bool refined) const override
- {
- container.compress();
- AMDiS::forEachLeafNode_(basis_->localView().tree(), [&](auto const& node, auto const& treePath)
- {
- auto subBasis = Dune::Functions::subspaceBasis(*basis_, treePath);
- nodeDataTransfer_[node].postAdapt(subBasis, container.vector(), refined);
- });
- }
-
- private:
- Basis const* basis_;
- TreeData<Basis, NDT, true> nodeDataTransfer_;
- };
-
+ class DataTransfer;
/// Factory to create DataTransfer objects based on the \ref DataTransferOperation
template <class Container>
@@ -116,7 +61,7 @@ namespace AMDiS
case NO_OPERATION:
return std::make_unique<NoDataTransfer<Container>>();
case INTERPOLATE:
- return std::make_unique<DataTransfer<Container, Basis>>(basis);
+ return std::make_unique< DataTransfer<Container, Basis> >(basis);
default:
error_exit("Invalid data transfer\n");
return nullptr; // avoid warnings
@@ -125,3 +70,5 @@ namespace AMDiS
};
} // end namespace AMDiS
+
+#include "DataTransfer.inc.hpp"
diff --git a/src/amdis/DataTransfer.inc.hpp b/src/amdis/DataTransfer.inc.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9040221aad2fa15da2ee23cf2dcbee74502497e9
--- /dev/null
+++ b/src/amdis/DataTransfer.inc.hpp
@@ -0,0 +1,470 @@
+#pragma once
+
+#include <cmath>
+#include <functional>
+#include <limits>
+#include <map>
+#include <memory>
+#include <numeric>
+#include <type_traits>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/hash.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/utility/ConcurrentCache.hpp>
+#include <amdis/utility/TreeContainer.hpp>
+#include <amdis/utility/Visitor.hpp>
+
+namespace AMDiS
+{
+ template <class Node, class Container, class Basis>
+ class NodeDataTransfer {};
+
+
+ /** 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>
+ {
+ 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
+ {
+ 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;
+ }
+ };
+
+ using CacheImp = std::unordered_map<LocalCoordinate, BoolCoordPair, Hasher>;
+ using ChildCache = ConcurrentCache<LocalCoordinate, BoolCoordPair, ConsecutivePolicy, CacheImp>;
+
+ 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{};
+ }
+ };
+
+ using NodeDataTransferContainer = std::decay_t<decltype(
+ 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 = std::decay_t<decltype(
+ 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 Container, class Basis>
+ void DataTransfer<Container, Basis>::preAdapt(Container const& coeff, bool mightCoarsen)
+ {
+ auto gv = basis_->gridView();
+ auto lv = basis_->localView();
+ auto const& idSet = gv.grid().localIdSet();
+
+ forEachLeafNode_(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 it = persistentContainer_.emplace(idSet.id(e),
+ makeTreeContainer<Tree, NodeElementData>(lv.tree(), NodeElementData()));
+
+ lv.bind(e);
+ auto& treeContainer = it.first->second;
+ forEachLeafNode_(lv.tree(), [&](auto const& node, auto const& tp) {
+ nodeDataTransfer_[tp].cacheLocal(treeContainer[tp]);
+ });
+ }
+
+ // Interpolate from possibly vanishing elements
+ if (mightCoarsen) {
+ auto maxLevel = gv.grid().maxLevel();
+ ctype const checkInsideTolerance = std::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::template checkInside<ctype,d>
+ (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;
+ forEachLeafNode_(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);
+ }
+ }
+ }
+ }
+ }
+
+ template <class Container, class Basis>
+ void DataTransfer<Container, Basis>::postAdapt(Container& coeff, bool refined)
+ {
+ coeff.resize(*basis_);
+ auto gv = basis_->gridView();
+ auto lv = basis_->localView();
+ auto const& idSet = gv.grid().localIdSet();
+ forEachLeafNode_(lv.tree(), [&](auto const& node, auto const& tp) {
+ nodeDataTransfer_[tp].postAdaptInit(lv, coeff, node);
+ });
+
+ mapper_.update();
+ finished_.clear();
+ finished_.resize(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;
+ forEachLeafNode_(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));
+ };
+
+ forEachLeafNode_(lv.tree(), [&](auto const& node, auto const& tp) {
+ nodeDataTransfer_[tp].prolongLocal(father, treeContainer[tp], xInFather, init);
+ });
+
+ finished_[mapper_.index(child)] = true;
+ init = false;
+ }
+ }
+ }
+ }
+ }
+
+ /** Element-local data transfer on a single leaf node of the basis tree
+ * Handles computations related to the finite element basis node
+ */
+ template <class GV, int k, class TP, class Container, class Basis>
+ class NodeDataTransfer<Dune::Functions::LagrangeNode<GV,k,TP>, Container, Basis>
+ {
+ using T = typename Container::value_type;
+ using LocalView = typename Basis::LocalView;
+ using Element = typename GV::template Codim<0>::Entity;
+
+ using Node = typename Dune::Functions::LagrangeNode<GV,k,TP>;
+ 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;
+ auto nodeCopy = node;
+ fatherNode_ = std::make_unique<Node>(std::move(nodeCopy));
+ constCoeff_ = &coeff;
+ }
+
+ /// 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;
+ auto nodeCopy = node;
+ fatherNode_ = std::make_unique<Node>(std::move(nodeCopy));
+ coeff_ = &coeff;
+ }
+
+ /// 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))];
+ }
+
+ /** 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);
+
+ /// 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];
+ }
+
+ /** 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_;
+ };
+
+ template <class GV, int k, class TP, class Container, class Basis>
+ template <class Trafo>
+ bool NodeDataTransfer<Dune::Functions::LagrangeNode<GV,k,TP>, Container, Basis>::
+ 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<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++];
+ };
+
+ 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 GV, int k, class TP, class Container, class Basis>
+ template <class Trafo>
+ void NodeDataTransfer<Dune::Functions::LagrangeNode<GV,k,TP>, Container, Basis>::
+ 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<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];
+
+ 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);
+
+ 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/common/Apply.hpp b/src/amdis/common/Apply.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..05bab973c292df65934fb2f3e4237032d93684b8
--- /dev/null
+++ b/src/amdis/common/Apply.hpp
@@ -0,0 +1,59 @@
+#pragma once
+
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+#include <amdis/common/Mpl.hpp>
+
+namespace AMDiS
+{
+ namespace Tools
+ {
+ namespace Impl_
+ {
+ template <class F, class Tuple, std::size_t... I>
+ constexpr decltype(auto) apply_impl(F&& f, Tuple&& t, std::index_sequence<I...>)
+ {
+ return f(std::get<I>(std::forward<Tuple>(t))...);
+ }
+
+ template <class F, std::size_t I0, std::size_t... I>
+ constexpr decltype(auto) apply_indices_impl(F&& f, index_t<I0>, std::index_sequence<I...>)
+ {
+ return f(index_t<I0+I>{}...);
+ }
+ } // namespace Impl_
+
+ template <class F, class Tuple>
+ constexpr decltype(auto) apply(F&& f, Tuple&& t)
+ {
+ return Impl_::apply_impl(
+ std::forward<F>(f), std::forward<Tuple>(t),
+ std::make_index_sequence<std::tuple_size<std::remove_reference_t<Tuple>>::value>{});
+ }
+
+ template <class F, class... Args>
+ constexpr decltype(auto) apply_variadic(F&& f, Args&&... args)
+ {
+ return apply(std::forward<F>(f), std::forward_as_tuple(args...));
+ }
+
+ template <class F, std::size_t N>
+ constexpr decltype(auto) apply_indices(F&& f, index_t<N>)
+ {
+ return Impl_::apply_indices_impl(
+ std::forward<F>(f),
+ index_t<0>{},
+ std::make_index_sequence<N>{});
+ }
+
+ template <class F, std::size_t I0, std::size_t I1>
+ constexpr decltype(auto) apply_indices(F&& f, index_t<I0>, index_t<I1>)
+ {
+ return Impl_::apply_indices_impl(
+ std::forward<F>(f), index_t<I0>{},
+ std::make_index_sequence<I1-I0>{});
+ }
+ }
+} // end namespace AMDiS
diff --git a/src/amdis/utility/TreeContainer.hpp b/src/amdis/utility/TreeContainer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cd013a67158ab5ef5033b0d4921fb5c4898151bf
--- /dev/null
+++ b/src/amdis/utility/TreeContainer.hpp
@@ -0,0 +1,199 @@
+#pragma once
+
+#include <array>
+#include <functional>
+#include <type_traits>
+#include <utility>
+
+#include <dune/common/indices.hh>
+#include <dune/common/tuplevector.hh>
+
+#include <dune/typetree/treepath.hh>
+
+#include <amdis/common/Apply.hpp>
+
+// NOTE: backport of dune/typetree/treecontainer.hh
+
+namespace AMDiS
+{
+ namespace Impl
+ {
+ /*
+ * \brief A factory class creating a hybrid container compatible with a type tree
+ *
+ * This class allows to create a nested hybrid container having the same structure
+ * as a given type tree. Power nodes are represented as std::array's while composite
+ * nodes are represented as Dune::TupleVector's. The stored values for the leaf nodes
+ * are creating using a given predicate. Once created, the factory provides an
+ * operator() creating the container for the tree given as argument.
+ *
+ * \tparam LeafToValue Type of a predicate that determines the stored values at the leafs
+ */
+ template<class LeafToValue>
+ class ContainerFactory
+ {
+ public:
+
+ /**
+ * \brief Create ContainerFactory
+ *
+ * The given predicate will be stored by value.
+ *
+ * \param A predicate used to generate the stored values for the leaves
+ */
+ ContainerFactory(LeafToValue leafToValue) :
+ leafToValue_(leafToValue)
+ {}
+
+ template<class Node,
+ std::enable_if_t<Node::isLeaf, int> = 0>
+ auto operator()(const Node& node)
+ {
+ return leafToValue_(node);
+ }
+
+ template<class Node,
+ std::enable_if_t<Node::isPower, int> = 0>
+ auto operator()(const Node& node)
+ {
+ using TransformedChild = decltype((*this)(node.child(0)));
+ return std::array<TransformedChild, Node::degree()>();
+ }
+
+ template<class Node,
+ std::enable_if_t<Node::isComposite, int> = 0>
+ auto operator()(const Node& node)
+ {
+ return Tools::apply_indices([&](auto... indices) { return Dune::makeTupleVector((*this)(node.child(indices))...); },
+ index_t<Node::degree()>{});
+ }
+
+ private:
+ LeafToValue leafToValue_;
+ };
+
+
+ /*
+ * \brief Wrap nested container to provide a VectorBackend
+ */
+ template<class Container>
+ class TreeContainerVectorBackend
+ {
+ template<class C>
+ static constexpr decltype(auto) accessByTreePath(C&& container, const Dune::TypeTree::HybridTreePath<>& path)
+ {
+ return container;
+ }
+
+ template<class C, class... T>
+ static constexpr decltype(auto) accessByTreePath(C&& container, const Dune::TypeTree::HybridTreePath<T...>& path)
+ {
+ auto head = Dune::TypeTree::treePathEntry(path,Dune::Indices::_0);
+ auto tailPath = Tools::apply_indices([&](auto... i)
+ {
+ using namespace Dune::TypeTree;
+ return hybridTreePath(treePathEntry(path,index_t<i.value+1>{})...);
+ }, index_t<sizeof...(T)-1>{});
+ return accessByTreePath(container[head], tailPath);
+ }
+
+ public:
+ TreeContainerVectorBackend(Container&& container) :
+ container_(std::move(container))
+ {}
+
+ TreeContainerVectorBackend(TreeContainerVectorBackend&& other) :
+ container_(std::move(other.container_))
+ {}
+
+ template<class... T>
+ decltype(auto) operator[](const Dune::TypeTree::HybridTreePath<T...>& path) const
+ {
+ return accessByTreePath(container_, path);
+ }
+
+ template<class... T>
+ decltype(auto) operator[](const Dune::TypeTree::HybridTreePath<T...>& path)
+ {
+ return accessByTreePath(container_, path);
+ }
+
+ const Container& data() const
+ {
+ return container_;
+ }
+
+ Container& data()
+ {
+ return container_;
+ }
+
+ private:
+ Container container_;
+ };
+
+ template<class Container>
+ auto makeTreeContainerVectorBackend(Container&& container)
+ {
+ return TreeContainerVectorBackend<std::decay_t<Container>>(std::forward<Container>(container));
+ }
+
+ } // end namespace Impl
+
+ /** \addtogroup TypeTree
+ * \{
+ */
+
+ /**
+ * \brief Create container havin the same structure as the given tree
+ *
+ * This class allows to create a nested hybrid container having the same structure
+ * as a given type tree. Power nodes are represented as std::array's while composite
+ * nodes are represented as Dune::TupleVector's. The stored values for the leaf nodes
+ * are creating using a given predicate. For convenience the created container is
+ * not returned directly. Instead, the returned object stores the container and
+ * provides operator[] access using a HybridTreePath.
+ *
+ * \param tree The tree which should be mapper to a container
+ * \param leafToValue A predicate used to generate the stored values for the leaves
+ *
+ * \returns A container matching the tree structure
+ */
+ template<class Tree, class LeafToValue>
+ auto makeTreeContainer(const Tree& tree, LeafToValue&& leafToValue)
+ {
+ auto f = std::ref(leafToValue);
+ auto factory = Impl::ContainerFactory<decltype(f)>(f);
+ return Impl::makeTreeContainerVectorBackend(factory(tree));
+ }
+
+ /**
+ * \brief Create container havin the same structure as the given tree
+ *
+ * This class allows to create a nested hybrid container having the same structure
+ * as a given type tree. Power nodes are represented as std::array's while composite
+ * nodes are represented as Dune::TupleVector's. The stored values for the leaf nodes
+ * are of the given type Value. For convenience the created container is
+ * not returned directly. Instead, the returned object stores the container and
+ * provides operator[] access using a HybridTreePath.
+ *
+ * \tparam Value Type of the values to be stored for the leafs. Should be default constructible.
+ * \param leafToValue A predicate used to generate the stored values for the leaves
+ *
+ * \returns A container matching the tree structure
+ */
+ template<class Value, class Tree>
+ auto makeTreeContainer(const Tree& tree)
+ {
+ return makeTreeContainer(tree, [](const auto&) {return Value{};});
+ }
+
+ /**
+ * \brief Alias to container type generated by makeTreeContainer for given value and tree type
+ */
+ template<class Value, class Tree>
+ using TreeContainer = std::decay_t<decltype(makeTreeContainer<Value>(std::declval<const Tree&>()))>;
+
+ //! \} group TypeTree
+
+} //namespace AMDiS
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f7819bc79529290f40abf4a17a9b798a99d0f011..bb4e86298bb60fb70e9d89fff6c5d804e2706a93 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -8,6 +8,12 @@ dune_add_test(SOURCES ClonablePtrTest.cpp
dune_add_test(SOURCES ConceptsTest.cpp
LINK_LIBRARIES amdis)
+dune_add_test(SOURCES DataTransferTest2d.cpp
+ LINK_LIBRARIES amdis)
+
+dune_add_test(SOURCES DataTransferTest3d.cpp
+ LINK_LIBRARIES amdis)
+
dune_add_test(SOURCES DOFVectorTest.cpp
LINK_LIBRARIES amdis)
diff --git a/test/DOFVectorTest.cpp b/test/DOFVectorTest.cpp
index 2f883436d854db6ecce3ba64af3706064542b831..61721ef05a6c2e1fd3bf95b396b4c35090a5a5db 100644
--- a/test/DOFVectorTest.cpp
+++ b/test/DOFVectorTest.cpp
@@ -53,13 +53,13 @@ int main(int argc, char** argv)
DOFVector<Basis> vec1(basis);
test_dofvector(basis, vec1);
- DOFVector<Basis, float> vec2(basis);
+ DOFVector<Basis, double> vec2(basis);
test_dofvector(basis, vec2);
auto vec3 = makeDOFVector(basis);
test_dofvector(basis, vec3);
- auto vec4 = makeDOFVector<float>(basis);
+ auto vec4 = makeDOFVector<double>(basis);
test_dofvector(basis, vec4);
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
diff --git a/test/DataTransferTest.hpp b/test/DataTransferTest.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7006d3cd9156d2f1c0c260b955f4c8da79eeacc7
--- /dev/null
+++ b/test/DataTransferTest.hpp
@@ -0,0 +1,166 @@
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <array>
+#include <cmath>
+#include <functional>
+#include <memory>
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#ifdef HAVE_DUNE_UGGRID
+#include <dune/grid/uggrid.hh>
+#else
+#include <dune/grid/yaspgrid.hh>
+#endif
+
+#include <amdis/ProblemStat.hpp>
+#include <amdis/ProblemStatTraits.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+template <class BasisCreator>
+auto makeGrid(bool simplex, int globalRefines = 0)
+{
+ using Grid = typename BasisCreator::GlobalBasis::GridView::Grid;
+ static constexpr int d = Grid::dimension; // problem dimension
+ using DomainType = typename Dune::FieldVector<double, d>;
+
+ // constants
+ DomainType lowerLeft; lowerLeft = 0.0; // lower left grid corner
+ DomainType upperRight; upperRight = 1.0; // upper right grid corner
+ std::array<unsigned int, d> s; s.fill(1); // number of elements on each axis
+
+ // make grid
+ std::unique_ptr<Grid> gridPtr;
+ if (simplex)
+ {
+ gridPtr = std::unique_ptr<Grid>{
+ Dune::StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, s)};
+ }
+ else
+ {
+ gridPtr = std::unique_ptr<Grid>{
+ Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, s)};
+ }
+ gridPtr->globalRefine(globalRefines);
+
+ return gridPtr;
+}
+
+template <class BasisCreator, class Fcts>
+auto makeProblem(typename BasisCreator::GlobalBasis::GridView::Grid& grid, Fcts const& funcs)
+{
+ using Problem = ProblemStat<BasisCreator>;
+
+ // make problem
+ Problem prob("test", grid);
+ prob.initialize(INIT_ALL);
+
+ auto& globalBasis = prob.globalBasis();
+ auto localView = globalBasis.localView();
+
+ // interpolate given function to initial grid
+ int k = 0;
+ AMDiS::forEachLeafNode_(localView.tree(), [&](auto const& node, auto tp)
+ {
+ interpolate(globalBasis, tp, prob.solution(tp).coefficients(), funcs[k]);
+ k++;
+ });
+
+ return prob;
+}
+
+template <class Problem, class Fcts>
+double calcError(Problem const& prob, Fcts const& funcs)
+{
+ auto& globalBasis = prob.globalBasis();
+ auto localView = globalBasis.localView();
+ auto sol = prob.solution().coefficients();
+ std::vector<double> ref;
+ ref.resize(globalBasis.size());
+ double error = 0;
+ double maxError = 0;
+ int k = 0;
+
+ // interpolate given function onto reference vector
+ AMDiS::forEachLeafNode_(localView.tree(), [&](auto const& node, auto tp)
+ {
+ interpolate(globalBasis, tp, ref, funcs[k]);
+ k++;
+ });
+
+ // compare the solution with the reference
+ for (std::size_t i = 0; i < ref.size(); ++i)
+ {
+ error = std::abs(ref[i]-sol[i]);
+ maxError = std::max(maxError, error);
+ }
+ return maxError;
+}
+
+template <class BasisCreator>
+ using Fcts = std::vector<std::function<double(
+ Dune::FieldVector<double, BasisCreator::GlobalBasis::GridView::Grid::dimension>)> >;
+
+/// Test data transfer for the case where no grid changes are made
+template <class BasisCreator>
+bool unchanged_test(Fcts<BasisCreator> const& funcs, bool simplex = true)
+{
+ using Grid = typename BasisCreator::GlobalBasis::GridView::Grid;
+ static constexpr int d = Grid::dimension; // problem dimension
+
+ auto gridPtr = makeGrid<BasisCreator>(simplex, (d > 2 ? 3 : 5));
+ auto prob = makeProblem<BasisCreator, Fcts<BasisCreator>>(*gridPtr, funcs);
+ // mark a single element for coarsening -> no adaptation is done
+ auto e = *gridPtr->leafGridView().template begin<0>();
+ gridPtr->mark(-1, e);
+ AdaptInfo adaptInfo("adapt");
+ prob.adaptGrid(adaptInfo);
+ auto error = calcError(prob, funcs);
+
+ return error < AMDIS_TEST_TOL;
+}
+
+/// Test data transfer for the case of global coarsening
+template <class BasisCreator>
+bool coarsen_test(Fcts<BasisCreator> const& funcs, bool simplex = true)
+{
+ using Grid = typename BasisCreator::GlobalBasis::GridView::Grid;
+ static constexpr int d = Grid::dimension; // problem dimension
+
+ auto gridPtr = makeGrid<BasisCreator>(simplex, (d > 2 ? 2 : 4));
+ auto prob = makeProblem<BasisCreator, Fcts<BasisCreator>>(*gridPtr, funcs);
+ prob.globalCoarsen(1);
+ auto error = calcError(prob, funcs);
+
+ return error < AMDIS_TEST_TOL;
+}
+
+/// Test data transfer for the case of global refinement
+template <class BasisCreator>
+bool refine_test(Fcts<BasisCreator> const& funcs, bool simplex = true)
+{
+ using Grid = typename BasisCreator::GlobalBasis::GridView::Grid;
+ static constexpr int d = Grid::dimension; // problem dimension
+
+ auto gridPtr = makeGrid<BasisCreator>(simplex, (d > 2 ? 1 : 3));
+ auto prob = makeProblem<BasisCreator, Fcts<BasisCreator>>(*gridPtr, funcs);
+ prob.globalRefine(1);
+ auto error = calcError(prob, funcs);
+
+ return error < AMDIS_TEST_TOL;
+}
+
+template <class Grid>
+ using Lagrange3 = LagrangeBasis<typename Grid::LeafGridView, 3>;
+template <class Grid>
+ using TaylorHood = TaylorHoodBasis<typename Grid::LeafGridView>;
+
+ constexpr double pi = 3.141592653589793238463;
diff --git a/test/DataTransferTest2d.cpp b/test/DataTransferTest2d.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2f4e92b2f8393fb63a4f46610ff8f68d6fb773e3
--- /dev/null
+++ b/test/DataTransferTest2d.cpp
@@ -0,0 +1,37 @@
+
+#include "DataTransferTest.hpp"
+
+int main(int argc, char** argv)
+{
+#ifdef HAVE_DUNE_UGGRID
+ using Grid = Dune::UGGrid<2>;
+#else
+ using Grid = Dune::YaspGrid<2>;
+#endif
+
+ Dune::MPIHelper::instance(argc, argv);
+
+ using Domain = typename Dune::FieldVector<double, 2>;
+
+ // polynomial of order 1
+ auto p1 = [](const Domain& x) -> double { return {0.5 + 0.25*(x[0]+x[1])}; };
+
+ // polynomial of order 2
+ auto p2 = [](const Domain& x) -> double { return {0.5 + 0.25*(2*std::pow(x[0],2) + x[0]*x[1] - std::pow(x[1],2))}; };
+
+ // polynomial of order 3
+ auto p3 = [](const Domain& x) -> double { return {0.5 + 0.25*(2*std::pow(x[0],3) + x[0]*x[1] - std::pow(x[1],3))}; };
+
+ // analytic function
+ auto f = [](const Domain& x) -> double { return {0.5 + 0.25*(std::sin(2*pi*x[0]) + 0.25*std::cos(6*pi*x[1]))}; };
+
+ AMDIS_TEST( (unchanged_test< Lagrange3<Grid> >({f})) );
+ AMDIS_TEST( (coarsen_test< Lagrange3<Grid> >({f})) );
+ AMDIS_TEST( (refine_test< Lagrange3<Grid> >({p3})) );
+
+ AMDIS_TEST( (unchanged_test< TaylorHood<Grid> >({f,f,f})) );
+ AMDIS_TEST( (coarsen_test< TaylorHood<Grid> >({f,f,f})) );
+ AMDIS_TEST( (refine_test< TaylorHood<Grid> >({p2,p2,p1})) );
+
+ return report_errors();
+}
\ No newline at end of file
diff --git a/test/DataTransferTest3d.cpp b/test/DataTransferTest3d.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e1b818c7e6d20564419f6a6e81b27c724ee121ba
--- /dev/null
+++ b/test/DataTransferTest3d.cpp
@@ -0,0 +1,37 @@
+
+#include "DataTransferTest.hpp"
+
+int main(int argc, char** argv)
+{
+#ifdef HAVE_DUNE_UGGRID
+ using Grid = Dune::UGGrid<3>;
+#else
+ using Grid = Dune::YaspGrid<3>;
+#endif
+
+ Dune::MPIHelper::instance(argc, argv);
+
+ using Domain = typename Dune::FieldVector<double, 3>;
+
+ // polynomial of order 1
+ auto p1 = [](const Domain& x) -> double { return {0.5 + 0.25*(x[0]+x[2])}; };
+
+ // polynomial of order 2
+ auto p2 = [](const Domain& x) -> double { return {0.5 + 0.25*(2*std::pow(x[0],2) + x[0]*x[1] - std::pow(x[2],2))}; };
+
+ // polynomial of order 3
+ auto p3 = [](const Domain& x) -> double { return {0.5 + 0.25*(2*std::pow(x[0],3) + x[0]*x[1] - std::pow(x[2],3))}; };
+
+ // analytic function
+ auto f = [](const Domain& x) -> double { return {0.5 + 0.25*(std::sin(2*pi*x[0]) + 0.25*std::cos(6*pi*x[2]))}; };
+
+ AMDIS_TEST( (unchanged_test< Lagrange3<Grid> >({f})) );
+ AMDIS_TEST( (coarsen_test< Lagrange3<Grid> >({f})) );
+ AMDIS_TEST( (refine_test< Lagrange3<Grid> >({p3})) );
+
+ AMDIS_TEST( (unchanged_test< TaylorHood<Grid> >({f,f,f,f})) );
+ AMDIS_TEST( (coarsen_test< TaylorHood<Grid> >({f,f,f,f})) );
+ AMDIS_TEST( (refine_test< TaylorHood<Grid> >({p2,p2,p2,p1})) );
+
+ return report_errors();
+}
\ No newline at end of file
diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp
index 0862d4ff73649eeea0814577fc15b4d5b35765db..14419ce8a50a1991aaa93cd52f6829a65c16b104 100644
--- a/test/DiscreteFunctionTest.cpp
+++ b/test/DiscreteFunctionTest.cpp
@@ -97,7 +97,7 @@ int main(int argc, char** argv)
dlocalFct.unbind();
}
- auto V0 = makeDOFVector<float>(prob.globalBasis());
+ auto V0 = makeDOFVector<double>(prob.globalBasis());
auto v0 = makeDOFVectorView(V0);
v0 << expr;