diff --git a/examples/periodic.cc b/examples/periodic.cc
index 02218df1a8db3d8c677f8ac94d56eb9287bab05b..8236d99c8bf05b7df63d1141fb3aaee34359360c 100644
--- a/examples/periodic.cc
+++ b/examples/periodic.cc
@@ -65,7 +65,7 @@ void run(Grid& grid)
using BC = PeriodicBC<FieldVector<double,2>, typename Traits::GlobalBasis::MultiIndex>;
BC periodicBC(prob.boundaryManagerPtr(),-1,{{{1.0,0.0}, {0.0,1.0}}, {1.0, 0.0}});
- periodicBC.init(prob.globalBasis(), prob.globalBasis());
+ periodicBC.init(*prob.globalBasis(), *prob.globalBasis());
std::cout << "periodicNodes:\n";
for (std::size_t i = 0; i < periodicBC.periodicNodes().size(); ++i)
std::cout << i << ": " << periodicBC.periodicNodes()[i] << "\n";
diff --git a/src/amdis/AdaptiveGrid.hpp b/src/amdis/AdaptiveGrid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e9303e8c10c75b446273b23620e3862b28fa392c
--- /dev/null
+++ b/src/amdis/AdaptiveGrid.hpp
@@ -0,0 +1,307 @@
+#pragma once
+
+#include <algorithm>
+#include <list>
+#include <memory>
+#include <mutex>
+#include <shared_mutex>
+#include <type_traits>
+#include <utility>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/geometry/type.hh>
+
+#include <amdis/common/ConceptsBase.hpp>
+#include <amdis/common/SharedPtr.hpp>
+#include <amdis/Observer.hpp>
+#include <amdis/Output.hpp>
+
+namespace AMDiS
+{
+ namespace Impl
+ {
+ template <class Grid, class... Args,
+ REQUIRES(std::is_convertible<decltype(std::declval<Grid>().loadBalance()), bool>::value)>
+ bool loadBalanceImpl(Grid& grid, Args&&... args)
+ {
+ return grid.loadBalance(FWD(args)...);
+ }
+
+ // Workaround for grids not adhering to the correct interface
+ template <class Grid, class... Args,
+ REQUIRES(not std::is_convertible<decltype(std::declval<Grid>().loadBalance()), bool>::value)>
+ bool loadBalanceImpl(Grid& grid, Args&&... args)
+ {
+ grid.loadBalance(FWD(args)...);
+ return true;
+ }
+ }
+
+ namespace event
+ {
+ /** Event generated from an AdaptiveGrid when calling preAdapt(). It contains the return value
+ * of preAdapt() as its 'mightCoarsen' member and is passed to registered observers after
+ * calling preAdapt on the underlying grid.
+ */
+ struct preAdapt { bool mightCoarsen; };
+
+ /** Event generated from an AdaptiveGrid when calling adapt(). Its 'adapted' member contains the
+ * value true if either preAdapt() or adapt() returned true. This event is passed to registered
+ * observers after calling adapt on the underlying grid.
+ */
+ struct adapt { bool adapted; };
+
+ /** Event generated from an AdaptiveGrid when calling postAdapt().This event is passed to
+ * registered observers after calling postAdapt on the underlying grid.
+ */
+ struct postAdapt {};
+ }
+
+ /** Wrapper class for Dune-grids that allows automatic signalling of events during grid adaptation
+ * This class needs to be created after construction of the associated grid by calling the
+ * instance method.
+ * Calls to grid.preAdapt(), grid.adapt() and grid.postAdapt() need to be replaced by calls to
+ * the corresponding methods of this class to use the automatic update functionality.
+ */
+ template <class Grid>
+ class AdaptiveGrid
+ : public Signals<event::preAdapt, event::adapt, event::postAdapt>
+ {
+ using Self = AdaptiveGrid<Grid>;
+ using Element = typename Grid::template Codim<0>::Entity;
+ struct HiddenStruct {};
+
+ public:
+ using HostGrid = Grid;
+ enum { dimension = HostGrid::dimension };
+ enum { dimensionworld = HostGrid::dimensionworld };
+ using LeafGridView = typename HostGrid::LeafGridView;
+ using LevelGridView = typename HostGrid::LevelGridView;
+
+ template <int cd>
+ struct Codim
+ {
+ using Geometry = typename HostGrid::template Codim<cd>::Geometry;
+ using Entity = typename HostGrid::template Codim<cd>::Entity;
+ };
+
+ using GlobalIdSet = typename HostGrid::GlobalIdSet;
+ using LocalIdSet = typename HostGrid::LocalIdSet;
+ using LevelIndexSet = typename HostGrid::LevelIndexSet;
+ using LeafIndexSet = typename HostGrid::LeafIndexSet;
+
+ /// Constructor that may only be called indirectly via the instance method
+ explicit AdaptiveGrid(std::shared_ptr<HostGrid> grid, HiddenStruct)
+ : hostGrid_(std::move(grid))
+ {}
+
+ /// Unreachable constructor required to compile the unreachable branch in instance method
+ explicit AdaptiveGrid(std::shared_ptr<HostGrid const> grid, HiddenStruct)
+ {
+ error_exit("Cannot create AdaptiveGrid from const Grid.");
+ }
+
+ int maxLevel() const { return hostGrid_->maxLevel(); }
+
+ int size(int level, int codim) const { return hostGrid_->size(level, codim); }
+
+ /// Return number of leaf entities of a given codim in this process
+ int size(int codim) const { return hostGrid_->size(codim); }
+
+ /// Return number of entities per level and geometry type in this process
+ int size(int level, Dune::GeometryType type) const { return hostGrid_->size(level, type); }
+
+ /// Return number of leaf entities per geometry type in this process
+ int size(Dune::GeometryType type) const { return hostGrid_->size(type); }
+
+ auto numBoundarySegments () const { return hostGrid_->numBoundarySegments(); }
+
+ /// View for a grid level for All_Partition
+ LevelGridView levelGridView(int l) const { return hostGrid_->levelGridView(l); }
+
+ /// View for the leaf grid for All_Partition
+ LeafGridView leafGridView() const { return hostGrid_->leafGridView(); }
+
+ GlobalIdSet const& globalIdSet() const { return hostGrid_->globalIdSet(); }
+
+ /// return const reference to the host grids local id set
+ LocalIdSet const& localIdSet() const { return hostGrid_->localIdSet(); }
+
+ /// return const reference to the host grids level index set for level level
+ LevelIndexSet const& levelIndexSet(int level) const { return hostGrid_->levelIndexSet(level); }
+
+ /// return const reference to the host grids leaf index set
+ LeafIndexSet const& leafIndexSet() const { return hostGrid_->leafIndexSet(); }
+
+ /// Refines all grid elements refCount times.
+ void globalRefine(int refCount)
+ {
+ for (int i = 0; i < refCount; ++i) {
+ // mark all entities for grid refinement
+ for (const auto& element : elements(hostGrid_->leafGridView()))
+ hostGrid_->mark(1, element);
+ preAdapt();
+ adapt();
+ postAdapt();
+ }
+ }
+
+ /// Mark entity for refinement
+ bool mark(int refCount, Element const& e) { return hostGrid_->mark(refCount, e); }
+
+ /// Return refinement mark for entity
+ int getMark(Element const& e) const { return hostGrid_->getMark(e); }
+
+ /// Prepare the grid for adaptation and notify observers of the preAdapt event
+ bool preAdapt()
+ {
+ Dune::Timer t;
+ mightCoarsen_ = hostGrid_->preAdapt();
+ Dune::MPIHelper::getCollectiveCommunication().max(&mightCoarsen_, 1);
+ this->notify(event::preAdapt{mightCoarsen_});
+ info(2,"preAdapt needed {} seconds", t.elapsed());
+ return mightCoarsen_;
+ }
+
+ /// Adapt the grid and notify observers of the adapt event
+ bool adapt()
+ {
+ Dune::Timer t;
+ refined_ = hostGrid_->adapt();
+ Dune::MPIHelper::getCollectiveCommunication().max(&refined_, 1);
+ this->notify(event::adapt{mightCoarsen_ || refined_});
+ return refined_;
+ }
+
+ /// Perform cleanup after grid adaptation and notify observers of the postAdapt event
+ void postAdapt()
+ {
+ Dune::Timer t;
+ hostGrid_->postAdapt();
+ this->notify(event::postAdapt{});
+ changeIndex_++;
+ info(2,"postAdapt needed {} seconds", t.elapsed());
+ }
+
+ /** Calls loadBalance on the underlying grid.
+ * Re-balances the load each process has to handle for a parallel grid.
+ * \return True if the grid has changed.
+ */
+ bool loadBalance()
+ {
+ return Impl::loadBalanceImpl(*hostGrid_);
+ }
+
+ /* Calls loadBalance(handle) on the underlying grid.
+ * Re-balances the load each process has to handle for a parallel grid and moves the data.
+ * \param data A data handle telling the method which data should be communicated
+ * and how. Has to adhere to the interface described by Dune::CommDataHandleIf.
+ * \return True if the grid has changed.
+ */
+ template <class DataHandle>
+ bool loadBalance (DataHandle& handle)
+ {
+ return Impl::loadBalanceImpl(*hostGrid_, handle);
+ }
+
+ /// Returns the grid change index, see \ref changeIndex.
+ unsigned long changeIndex() const
+ {
+ return changeIndex_;
+ }
+
+ private:
+ template <class Grid_>
+ static std::shared_ptr<Self>
+ instanceImpl(std::shared_ptr<Grid_> grid)
+ {
+ using mutex_type = std::shared_timed_mutex;
+ static mutex_type access_mutex;
+
+ // 1. Cleanup. Since the iterators may be invalidated only one accessor is allowed. Erases all
+ // expired weak pointers (all pointers that no longer have valid instances attached).
+ std::unique_lock<mutex_type> write_lock(access_mutex);
+ for (auto it = adaptiveGrids_.begin(); it != adaptiveGrids_.end(); ++it)
+ if (it->expired())
+ it = adaptiveGrids_.erase(it);
+ write_lock.unlock();
+
+ std::shared_lock<mutex_type> read_lock(access_mutex);
+ // 2. Find matching observed grid object. We obtain a lock here to avoid complications with
+ // insertions made by other threads in step 3.
+ auto it = find_if(adaptiveGrids_.begin(), adaptiveGrids_.end(),
+ [&](std::weak_ptr<Self>& wPtr) {
+ auto ag = wPtr.lock();
+ return ag->hostGrid_ == grid;
+ });
+
+ // 3. Register new object or return existing. In case of inserting a new instance we obtain a
+ // write lock.
+ if (it == adaptiveGrids_.end())
+ {
+ test_exit(!std::is_const<Grid_>::value,
+ "No existing AdaptiveGrid found and no mutable grid passed to create a new instance");
+ auto ptr = std::make_shared<Self>(std::move(grid), HiddenStruct{});
+ read_lock.unlock();
+ write_lock.lock();
+ adaptiveGrids_.emplace_back(ptr);
+ return std::move(ptr);
+ }
+ else
+ {
+ return it->lock();
+ }
+ }
+
+ // Do-nothing specialization if argument is already an AdaptiveGrid
+ static std::shared_ptr<Self> instanceImpl(std::shared_ptr<Self> grid) { return grid; }
+
+ public:
+ /** Returns the AdaptiveGrid associated to the grid passed.
+ * If no AdaptiveGrid exists, this method creates a new one if the passed grid is mutable,
+ * otherwise the call fails with an error.
+ */
+ template <class Grid_>
+ static std::shared_ptr<Self> instance(Grid_&& grid)
+ {
+ return instanceImpl(wrap_or_share(FWD(grid)));
+ }
+
+ // Test for equality of the grid pointers
+ bool operator==(AdaptiveGrid<HostGrid> const& that) const
+ {
+ return (hostGrid_ == that.hostGrid_);
+ }
+
+ /// Return the underlying grid
+ std::shared_ptr<HostGrid> hostGrid() const { return hostGrid_; }
+
+ private:
+ /// List of previously created instances, indexed by address of the HostGrid.
+ /**
+ * For each grid type Grid we maintain a list of instances handed out by the instance method.
+ * We use weak pointers here that are valid as long as there is at least one other place where
+ * the shared pointer to the instance is used. When the instance is no longer used the weak
+ * pointers here do not hinder the object's destruction.
+ */
+ static std::list<std::weak_ptr<Self>> adaptiveGrids_;
+
+ /// The underlying grid implementation
+ std::shared_ptr<HostGrid> hostGrid_;
+
+ /// Flag set during \ref preAdapt(), indicating whether any element might be coarsened in \ref adapt()
+ bool mightCoarsen_ = false;
+
+ /// Flag set during \ref adapt() indicating that at least one entity was refined
+ bool refined_ = false;
+
+ /// This index is incremented every time the grid is changed, e.g. by refinement or coarsening.
+ unsigned long changeIndex_ = 0;
+ };
+
+ template <class Grid>
+ std::list<std::weak_ptr<AdaptiveGrid<Grid>>> AdaptiveGrid<Grid>::adaptiveGrids_;
+
+} // end namespace AMDiS
diff --git a/src/amdis/BiLinearForm.hpp b/src/amdis/BiLinearForm.hpp
index ff1fda26299f625bd62f962bc52c530a2b0058ce..d3864d6d22c0c75472109f49c3b5ca43aaa22320 100644
--- a/src/amdis/BiLinearForm.hpp
+++ b/src/amdis/BiLinearForm.hpp
@@ -45,12 +45,11 @@ namespace AMDiS
using ElementMatrix = FlatMatrix<CoefficientType>;
public:
- /// Constructor. Wraps the reference into a non-destroying shared_ptr or moves
- /// the basis into a new shared_ptr. The additional arguments can be the communication
- /// object. If not provided, a new communication is created.
- template <class RB_, class CB_, class... Args>
- BiLinearForm(RB_&& rowBasis, CB_&& colBasis, Args&&... args)
- : Super(FWD(rowBasis), FWD(colBasis), FWD(args)...)
+ /// Constructor. Wraps the reference into a non-destroying shared_ptr or moves the basis into
+ /// a new shared_ptr.
+ template <class RB_, class CB_>
+ BiLinearForm(RB_&& rowBasis, CB_&& colBasis)
+ : Super(FWD(rowBasis), FWD(colBasis))
{
operators_.init(*Super::rowBasis(), *Super::colBasis());
}
@@ -123,16 +122,16 @@ namespace AMDiS
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
- template <class RB, class CB, class... Args>
- BiLinearForm(RB&&, CB&&, Args&&...)
+ template <class RB, class CB>
+ BiLinearForm(RB&&, CB&&)
-> BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;
#endif
template <class RB, class CB, class... Args>
- auto makeBiLinearForm(RB&& rowBasis, CB&& colBasis, Args&&... args)
+ auto makeBiLinearForm(RB&& rowBasis, CB&& colBasis)
{
using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;
- return BLF{FWD(rowBasis), FWD(colBasis), FWD(args)...};
+ return BLF{FWD(rowBasis), FWD(colBasis)};
}
} // end namespace AMDiS
diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt
index 25ca9c6c9d9959bd57412260359a4449e02c05fc..3527447c658906ad0ace05eee338eaa1f0a06e6d 100644
--- a/src/amdis/CMakeLists.txt
+++ b/src/amdis/CMakeLists.txt
@@ -18,6 +18,7 @@ install(FILES
AdaptInfo.hpp
AdaptInstationary.hpp
AdaptionInterface.hpp
+ AdaptiveGrid.hpp
AdaptStationary.hpp
AMDiS.hpp
Assembler.hpp
@@ -37,13 +38,10 @@ install(FILES
DirichletBC.inc.hpp
DOFVector.hpp
DOFVector.inc.hpp
- DOFVectorInterface.hpp
Environment.hpp
Flag.hpp
GridFunctionOperator.hpp
GridFunctions.hpp
- GridTransfer.hpp
- GridTransferManager.hpp
Initfile.hpp
InitfileParser.hpp
Integrate.hpp
@@ -55,6 +53,7 @@ install(FILES
Marker.hpp
Marker.inc.hpp
MeshCreator.hpp
+ Observer.hpp
Operations.hpp
OperatorList.hpp
Output.hpp
diff --git a/src/amdis/DOFVector.hpp b/src/amdis/DOFVector.hpp
index 9a606d9be2508c95844595e2a1f61866e99e6eb8..7b1e3235fb02749daf580c6bddfaa6788dabbbdc 100644
--- a/src/amdis/DOFVector.hpp
+++ b/src/amdis/DOFVector.hpp
@@ -6,16 +6,35 @@
#include <dune/common/shared_ptr.hh>
#include <amdis/DataTransfer.hpp>
-#include <amdis/DOFVectorInterface.hpp>
-#include <amdis/GridTransferManager.hpp>
#include <amdis/LinearAlgebra.hpp>
+#include <amdis/Observer.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
#include <amdis/typetree/TreePath.hpp>
+namespace Dune
+{
+ namespace Functions
+ {
+ template <class PB>
+ class DefaultGlobalBasis;
+ }
+}
+
namespace AMDiS
{
+ // Forward declarations
+ template <class PB>
+ class ParallelGlobalBasis;
+
+ namespace event
+ {
+ struct preAdapt;
+ struct postAdapt;
+ template <class B> struct basisUpdate;
+ }
+
/// \brief The basic container that stores a base vector and a corresponding basis
/**
* \tparam GB Basis of the vector
@@ -24,14 +43,14 @@ namespace AMDiS
template <class GB, class T = double>
class DOFVector
: public VectorBase<GB, VectorBackend<BackendTraits<GB,T>>>
- , public DOFVectorInterface
+ , public Observer<GB, event::preAdapt, event::basisUpdate<GB>, event::postAdapt>
{
using Self = DOFVector;
using Super = VectorBase<GB, VectorBackend<BackendTraits<GB,T>>>;
+ using Obs = Observer<GB, event::preAdapt, event::basisUpdate<GB>, event::postAdapt>;
public:
using Backend = VectorBackend<BackendTraits<GB,T>>;
- using Comm = typename Backend::Traits::Comm;
using GlobalBasis = GB;
@@ -41,89 +60,32 @@ namespace AMDiS
/// The type of the elements of the DOFVector
using value_type = typename Backend::value_type;
- /// Defines an interface to transfer the data during grid adaption
- using DataTransfer = DataTransferInterface<Self>;
-
- /// A creator for a concrete data transfer object, depending on \ref DataTransferOperation
- using DataTransferFactory = AMDiS::DataTransferFactory<Self>;
+ /// Wrapper for the implementation of the transfer of data during grid adaption
+ using DataTransfer = DataTransferWrapper<Self>;
public:
/// Constructor. Stores the shared_ptr of the basis and creates a new DataTransfer.
- DOFVector(std::shared_ptr<GlobalBasis> basis, std::shared_ptr<Comm> comm, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
- : Super(basis, std::move(comm))
- , dataTransfer_(DataTransferFactory::create(op, basis))
- {
- attachToGridTransfer();
- }
-
- /// Forwarding constructor that wraps the arguments into shared_ptr
- template <class GB_, class Comm_,
- REQUIRES(Concepts::Similar<GB_,GB>),
- REQUIRES(Concepts::Similar<Comm_,Comm>)>
- DOFVector(GB_&& basis, Comm_&& comm, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
- : DOFVector(Dune::wrap_or_move(FWD(basis)), Dune::wrap_or_move(FWD(comm)), op)
+ DOFVector(std::shared_ptr<GB> basis,
+ DataTransferOperation op = DataTransferOperation::INTERPOLATE)
+ : Super(basis)
+ , Obs(basis)
+ , dataTransfer_(op, basis)
+ , basis_(basis)
{}
+ /// Constructor creating a new basis from a Dune::Functions::DefaultGlobalBasis.
+ // TODO(FM): Replace explicit type with concept check
+ DOFVector(Dune::Functions::DefaultGlobalBasis<typename GB::PreBasis>&& basis,
+ DataTransferOperation op = DataTransferOperation::INTERPOLATE)
+ : DOFVector(std::make_shared<GB>(std::move(basis)), op)
+ {}
- /// Construct the DOFVector from a basis by first creating a comm object. This constructor
- /// Registers the basis and comm object into the \ref GridTransferManager so that it gets updated
- /// on grid changes.
- DOFVector(std::shared_ptr<GlobalBasis> basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
- : DOFVector(basis, CommunicationCreator<Comm>::create(*basis), op)
- {
- GridTransferManager::attach(this->basis()->gridView().grid(), *this->comm(), [c=this->comm(), b=this->basis()]() -> void {
- b->update(b->gridView());
- c->update(*b);
- });
- }
-
- /// Forwarding constructor that wraps the arguments into shared_ptr
template <class GB_,
REQUIRES(Concepts::Similar<GB_,GB>)>
DOFVector(GB_&& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
: DOFVector(Dune::wrap_or_move(FWD(basis)), op)
{}
- /// Copy constructor
- DOFVector(Self const& that)
- : Super(static_cast<Super const&>(that))
- , dataTransfer_(that.dataTransfer_)
- {
- attachToGridTransfer();
- }
-
- /// Move constructor
- DOFVector(Self&& that)
- : Super(static_cast<Super&&>(that))
- , dataTransfer_(std::move(that.dataTransfer_))
- {
- attachToGridTransfer();
- }
-
- /// Destructor
- ~DOFVector() override
- {
- detachFromGridTransfer();
- }
-
- /// Copy assignment operator
- Self& operator=(Self const& that)
- {
- detachFromGridTransfer();
- Super::operator=(that);
- dataTransfer_ = that.dataTransfer_;
- attachToGridTransfer();
- return *this;
- }
-
- /// Move assignment
- Self& operator=(Self&& that) = default;
-
- void resize() override
- {
- Super::resize();
- }
-
template <class TreePath = RootTreePath>
auto child(TreePath const& path = {})
@@ -175,71 +137,68 @@ namespace AMDiS
/// Return the associated DataTransfer object
- std::shared_ptr<DataTransfer const> dataTransfer() const
+ DataTransfer const& dataTransfer() const
{
return dataTransfer_;
}
/// Return the associated DataTransfer object
- std::shared_ptr<DataTransfer> dataTransfer()
+ DataTransfer& dataTransfer()
{
return dataTransfer_;
}
- /// Create a new DataTransfer object based on the operation type
- void setDataTransfer(DataTransferOperation op)
- {
- dataTransfer_ = DataTransferFactory::create(op, this->basis());
- }
-
/// Assign the DataTransfer object
- void setDataTransfer(std::shared_ptr<DataTransfer> const& dataTransfer)
+ void setDataTransfer(DataTransfer const& dataTransfer)
{
dataTransfer_ = dataTransfer;
}
+ void setDataTransfer(DataTransfer&& dataTransfer)
+ {
+ dataTransfer_ = std::move(dataTransfer);
+ }
- /// Implementation of \ref DOFVectorInterface::preAdapt
- /// Redirects to a \ref DataTransfer object.
- void preAdapt(bool mightCoarsen) override
+ /// Create a new DataTransfer object based on the operation type
+ void setDataTransfer(DataTransferOperation op)
{
- dataTransfer_->preAdapt(*this, mightCoarsen);
+ setDataTransfer(newDataTransfer(op, basis_));
}
- /// Implementation of \ref DOFVectorInterface::postAdapt
- /// Redirects to a \ref DataTransfer object.
- void postAdapt(bool refined) override
+ /// Override of Observer update(event::preAdapt) method. Redirects to preAdapt method of the
+ /// \ref DataTransfer object.
+ void update(event::preAdapt const& e) override
{
- dataTransfer_->postAdapt(*this, refined);
+ dataTransfer_.preAdapt(*this, e.mightCoarsen);
}
- private:
- // register this DOFVector and its basis to the DataTransfer
- void attachToGridTransfer()
+ /// Override of Observer update(event::adapt) method. Redirects to adapt method of the
+ /// \ref DataTransfer object.
+ void update(event::basisUpdate<GB> const&) override
{
- GridTransferManager::attach(this->basis()->gridView().grid(), *this);
+ this->resize();
+ dataTransfer_.adapt(*this);
}
- // deregister this DOFVector and its basis from the DataTransfer
- void detachFromGridTransfer()
+ /// Override of Observer update(event::postAdapt) method. Redirects to postAdapt method of the
+ /// \ref DataTransfer object.
+ void update(event::postAdapt const&) override
{
- GridTransferManager::detach(this->basis()->gridView().grid(), *this);
+ dataTransfer_.postAdapt(*this);
}
private:
/// Data interpolation when the grid changes, set by default
/// to \ref DataTransferOperation::INTERPOLATE.
- std::shared_ptr<DataTransfer> dataTransfer_;
+ DataTransfer dataTransfer_;
+
+ std::shared_ptr<GB> basis_;
};
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
- template <class GB, class C>
- DOFVector(GB&& basis, C&& comm, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
- -> DOFVector<Underlying_t<GB>>;
-
template <class GB>
DOFVector(GB&& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
- -> DOFVector<Underlying_t<GB>>;
+ -> DOFVector<ParallelGlobalBasis<typename Underlying_t<GB>::PreBasis>>;
#endif
/// \brief Create a DOFVector from a basis.
@@ -252,21 +211,13 @@ namespace AMDiS
* DOFVector. The default is interpolation of the data to the new grid. See
* \ref DataTransferOperation for more options.
**/
- template <class ValueType = double, class GlobalBasis, class Communication>
- DOFVector<Underlying_t<GlobalBasis>, ValueType>
- makeDOFVector(GlobalBasis&& basis, Communication&& comm, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
- {
- return {FWD(basis), FWD(comm), op};
- }
-
- template <class ValueType = double, class GlobalBasis>
- DOFVector<Underlying_t<GlobalBasis>, ValueType>
- makeDOFVector(GlobalBasis&& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
+ template <class ValueType = double, class GB>
+ DOFVector<ParallelGlobalBasis<typename Underlying_t<GB>::PreBasis>, ValueType>
+ makeDOFVector(GB&& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
{
return {FWD(basis), op};
}
-
} // end namespace AMDiS
#include <amdis/DOFVector.inc.hpp>
diff --git a/src/amdis/DOFVectorInterface.hpp b/src/amdis/DOFVectorInterface.hpp
deleted file mode 100644
index 316186b35067862b30eb41c2ef4f2a160b8afa56..0000000000000000000000000000000000000000
--- a/src/amdis/DOFVectorInterface.hpp
+++ /dev/null
@@ -1,28 +0,0 @@
-#pragma once
-
-namespace AMDiS
-{
- class DOFVectorInterface
- {
- public:
- /// Virtual destructor
- virtual ~DOFVectorInterface() = default;
-
- /// Change dimension of DOFVector to dimension of basis
- virtual void resize() = 0;
-
- /// \brief Prepare the data interpolation between grids
- /**
- * \param mightCoarsen Indicator that any element in the grid might be
- * coarsened in adapt()
- **/
- virtual void preAdapt(bool mightCoarsen) = 0;
-
- /// \brief Apply the actual interpolation
- /**
- * \param refined Indicator that at least one entity was refined
- **/
- virtual void postAdapt(bool refined) = 0;
- };
-
-} // end namespace AMDiS
diff --git a/src/amdis/DataTransfer.hpp b/src/amdis/DataTransfer.hpp
index 01c740550b5927f75bef8c424ab330e18ead4091..3ca2a27cb39ccba754ce4f904ad82970fd6dfb76 100644
--- a/src/amdis/DataTransfer.hpp
+++ b/src/amdis/DataTransfer.hpp
@@ -29,11 +29,17 @@ namespace AMDiS
/// Virtual destructor
virtual ~DataTransferInterface() = default;
+ /// Clone method
+ virtual std::unique_ptr<DataTransferInterface> clone() const = 0;
+
/// Collect data that is needed before grid adaption
virtual void preAdapt(Container const& container, bool mightCoarsen) = 0;
/// Interpolate data to new grid after grid adaption
- virtual void postAdapt(Container& container, bool refined) = 0;
+ virtual void adapt(Container& container) = 0;
+
+ /// Perform cleanup after grid adaption
+ virtual void postAdapt(Container& container) = 0;
};
@@ -45,13 +51,22 @@ namespace AMDiS
class NoDataTransfer
: public DataTransferInterface<Container>
{
+ using Interface = DataTransferInterface<Container>;
+
public:
- void preAdapt(Container const& container, bool) override {}
+ std::unique_ptr<Interface> clone() const override
+ {
+ return std::make_unique<NoDataTransfer>();
+ }
+
+ void preAdapt(Container const&, bool) override {}
- void postAdapt(Container& container, bool) override
+ void adapt(Container& container) override
{
container.resize();
}
+
+ void postAdapt(Container&) override {}
};
@@ -81,8 +96,15 @@ namespace AMDiS
using NodeElementData = typename NodeDataTransfer<Node, Container, Basis>::NodeElementData;
using ElementData = TYPEOF(makeTreeContainer<NodeElementData>(std::declval<const Tree&>()));
+ using Interface = DataTransferInterface<Container>;
+
public:
- DataTransfer(std::shared_ptr<Basis> basis);
+ DataTransfer(std::shared_ptr<Basis const> basis);
+
+ std::unique_ptr<Interface> clone() const override
+ {
+ return std::make_unique<DataTransfer<Container, Basis>>(basis_);
+ }
/** Saves data contained in coeff in the PersistentContainer
* To be called after grid.preAdapt() and before grid.adapt()
@@ -92,9 +114,14 @@ namespace AMDiS
/** Unpacks data from the PersistentContainer
* To be called after grid.adapt() and before grid.postAdapt()
*/
- // [[expects: basis is update in gridView]]
+ // [[expects: basis is updated in gridView]]
// [[expects: comm is updated in basis]]
- void postAdapt(Container& coeff, bool refined) override;
+ void adapt(Container& coeff) override;
+
+ /** Performs cleanup
+ * To be called after grid.postAdapt()
+ */
+ void postAdapt(Container& coeff) override;
private:
/// The global basis
@@ -139,6 +166,50 @@ namespace AMDiS
}
};
+
+ /** Wrapper class for (No)DataTransfer objects.
+ * This wrapper class may be used instead of a pointer type to DataTransferInterface. It wraps
+ * and redirects to the implementation, while providing deep copy and move functionality.
+ * A virtual clone() method in the interface and implementation classes is required to perform
+ * the deep copy of the implementation.
+ */
+ template <class Container>
+ class DataTransferWrapper
+ {
+ using Interface = DataTransferInterface<Container>;
+ using Factory = DataTransferFactory<Container>;
+ using Self = DataTransferWrapper<Container>;
+
+ public:
+ template <class Basis>
+ DataTransferWrapper(DataTransferOperation op, std::shared_ptr<Basis> basis)
+ : impl_(std::move(Factory::create(op, basis)))
+ {}
+
+ DataTransferWrapper(Self const& that)
+ : impl_(std::move(that.impl_->clone()))
+ {}
+
+ DataTransferWrapper(Self&& that) = default;
+
+ Self& operator=(Self const& that)
+ {
+ impl_ = std::move(that.impl_->clone());
+ return *this;
+ }
+
+ Self& operator=(Self&& that) = default;
+
+ void preAdapt(Container const& c, bool b) { impl_->preAdapt(c, b); }
+
+ void adapt(Container& c) { impl_->adapt(c); }
+
+ void postAdapt(Container& c) { impl_->postAdapt(c); }
+
+ private:
+ std::unique_ptr<Interface> impl_;
+ };
+
/// @}
} // end namespace AMDiS
diff --git a/src/amdis/DataTransfer.inc.hpp b/src/amdis/DataTransfer.inc.hpp
index 64319984308b8a81e43fab67ddcba599a4be8bcb..bf84ccfc60d08724c2860b28224903c94e603279 100644
--- a/src/amdis/DataTransfer.inc.hpp
+++ b/src/amdis/DataTransfer.inc.hpp
@@ -43,8 +43,7 @@ struct CoordHasher
template <class C, class B>
-DataTransfer<C,B>::
-DataTransfer(std::shared_ptr<B> basis)
+DataTransfer<C,B>::DataTransfer(std::shared_ptr<B const> basis)
: basis_(std::move(basis))
, mapper_(basis_->gridView().grid(), Dune::mcmgElementLayout())
, nodeDataTransfer_(makeTreeContainer<NDT>(basis_->localView().tree()))
@@ -64,7 +63,7 @@ preAdapt(C const& coeff, bool mightCoarsen)
});
// Make persistent DoF container
- persistentContainer_.clear();
+ persistentContainer_.clear(); // Redundant if postAdapt was correctly called last cycle
for (const auto& e : elements(gv))
{
auto it = persistentContainer_.emplace(idSet.id(e), makeTreeContainer<NodeElementData>(lv.tree()));
@@ -147,16 +146,18 @@ preAdapt(C const& coeff, bool mightCoarsen)
template <class C, class B>
-void DataTransfer<C,B>::
-postAdapt(C& coeff, bool refined)
+void DataTransfer<C,B>::adapt(C& coeff)
{
coeff.init(false);
+ test_exit(!persistentContainer_.empty(), "No data was saved before adapting the grid, make "
+ "sure to call DataTransfer::preAdapt before calling adapt() on the grid");
+
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);
+ nodeDataTransfer_[tp].adaptInit(lv, coeff, node);
});
mapper_.update();
@@ -221,6 +222,13 @@ postAdapt(C& coeff, bool refined)
}
+template <class C, class B>
+void DataTransfer<C,B>::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
*/
@@ -283,8 +291,8 @@ public:
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)
+ /// 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;
@@ -293,7 +301,7 @@ public:
}
/// \brief Copy already existing data to element bound to node_
- // [[expects: postAdaptInit to be called before]]
+ // [[expects: adaptInit to be called before]]
void copyLocal(NodeElementData const& dofs) const
{
coeff_->scatter(*lv_, *node_, dofs, Assigner::assign{});
@@ -310,7 +318,7 @@ public:
* coordinates in father element
* \param init Father element is visited for the first time
**/
- // [[expects: postAdaptInit to be called before]]
+ // [[expects: adaptInit to be called before]]
template <class Trafo>
void prolongLocal(Element const& father, NodeElementData const& fatherDOFs,
Trafo const& trafo, bool init);
diff --git a/src/amdis/GridTransfer.hpp b/src/amdis/GridTransfer.hpp
deleted file mode 100644
index 406245fbccf899b8095d562c42c229c1ce7dab62..0000000000000000000000000000000000000000
--- a/src/amdis/GridTransfer.hpp
+++ /dev/null
@@ -1,203 +0,0 @@
-#pragma once
-
-#include <functional>
-#include <list>
-
-#include <amdis/AdaptionInterface.hpp>
-#include <amdis/Output.hpp>
-
-namespace AMDiS
-{
- /**
- * \addtogroup Adaption
- * @{
- **/
-
- /// \brief Implementation of \ref AdaptionInterface for concrete Grid type.
- template <class Grid>
- class GridTransfer
- : public AdaptionInterface
- {
- using Self = GridTransfer;
- using Key = std::uintptr_t;
-
- struct InterpolateCallback
- {
- std::function<void(bool)> preAdapt;
- std::function<void(bool)> postAdapt;
- int count = 0;
- };
-
- struct UpdateCallback
- {
- std::function<void(void)> update;
- int count = 0;
- };
-
- public:
- /// Bind a (mutable) grid to this GridTransfer. Must be called before any xxxAdapt() method.
- // [[ensures: bound()]]
- void bind(Grid& grid)
- {
- grid_ = &grid;
- }
-
- /// Release the grid from this GridTransfer
- // [[ensures: not bound()]]
- void unbind()
- {
- grid_ = nullptr;
- }
-
- /// Returns whether this GridTransfer is bound to a grid.
- bool bound() const
- {
- return grid_ != nullptr;
- }
-
- /// Register data to GridTransfer.
- /**
- * This inserts a callback that is called in \ref preAdapt() and another that is
- * called in \ref postAdapt() step.
- *
- * *Requirement:*
- * - Data must be a model of \ref Concepts::InterpolateData
- **/
- template <class Data,
- REQUIRES(Concepts::InterpolateData<Data>)>
- void attach(Data& data)
- {
- auto it = interpolateCallbacks_.emplace(Key(&data), InterpolateCallback{
- [&data](bool mightCoarsen) -> void { data.preAdapt(mightCoarsen); },
- [&data](bool refined) -> void { data.postAdapt(refined); }});
-
- it.first->second.count++;
- }
-
- /// Register data to GridTransfer.
- /**
- * This inserts a callback that is called after the \ref adapt() step.
- * The callback calls `data.update(data.gridView())`.
- *
- * *Requirement:*
- * - Data must be a model of \ref Concepts::UpdateData
- **/
- template <class Data,
- REQUIRES(Concepts::UpdateData<Data>)>
- void attach(Data& data)
- {
- attach(data, [&data]() -> void { data.update(data.gridView()); });
- }
-
- /// Register callback to GridTransfer
- template <class Data, class Callback>
- void attach(Data const& data, Callback&& callback)
- {
- auto it = updateCallbacks_.emplace(Key(&data), UpdateCallback{FWD(callback)});
- it.first->second.count++;
- }
-
- /// Unregister data from GridTransfer. Data is identified by its address.
- template <class Data,
- REQUIRES(Concepts::InterpolateData<Data>)>
- void detach(Data const& data)
- {
- eraseCallback(interpolateCallbacks_, Key(&data));
- }
-
- /// Unregister data from GridTransfer. Data is identified by its address.
- template <class Data,
- REQUIRES(Concepts::UpdateData<Data>)>
- void detach(Data const& data)
- {
- eraseCallback(updateCallbacks_, Key(&data));
- }
-
- /// Implements \ref GridTransferInterface::preAdapt
- // [[expects: bound()]]
- bool preAdapt() override
- {
- assert(bound());
- Dune::Timer t;
- mightCoarsen_ = grid_->preAdapt();
- Dune::MPIHelper::getCollectiveCommunication().max(&mightCoarsen_, 1);
-
- for (auto&& data : interpolateCallbacks_)
- data.second.preAdapt(mightCoarsen_);
-
- info(2,"preAdapt needed {} seconds", t.elapsed());
- return mightCoarsen_;
- }
-
- /// Implements \ref GridTransferInterface::adapt
- // [[expects: bound()]]
- bool adapt() override
- {
- assert(bound());
- Dune::Timer t;
- refined_ = grid_->adapt();
- Dune::MPIHelper::getCollectiveCommunication().max(&refined_, 1);
-
- info(2,"adapt needed {} seconds", t.elapsed());
- return refined_;
- }
-
- /// Implements \ref GridTransferInterface::postAdapt
- // [[expects: bound()]]
- void postAdapt() override
- {
- assert(bound());
- Dune::Timer t;
- if (mightCoarsen_ || refined_) {
- for (auto&& data : updateCallbacks_)
- data.second.update();
- for (auto&& data : interpolateCallbacks_)
- data.second.postAdapt(refined_);
- }
-
- grid_->postAdapt();
- changeIndex_++;
- info(2,"postAdapt needed {} seconds", t.elapsed());
- }
-
- /// Returns the grid change index, see \ref changeIndex.
- unsigned long changeIndex() const
- {
- return changeIndex_;
- }
-
- private:
- // remove entry from map if count <= 1, otherwise reduce count by 1
- template <class Map>
- void eraseCallback(Map& map, Key key)
- {
- auto it = map.find(key);
- if (it == map.end())
- warning("Data to detach not found");
- else if (it->second.count > 1)
- it->second.count--;
- else
- map.erase(it);
- }
-
- private:
- /// A grid this GridTransfer is bound to in \ref bind()
- Grid* grid_ = nullptr;
-
- /// A list of data containers handled during grid adaption
- std::map<Key,InterpolateCallback> interpolateCallbacks_;
- std::map<Key,UpdateCallback> updateCallbacks_;
-
- /// Flag set during \ref preAdapt(), indicating whether any element might be coarsened in \ref adapt()
- bool mightCoarsen_ = false;
-
- /// Flag set during \ref adapt() indicating that at least one entity was refined
- bool refined_ = false;
-
- /// This index is incremented every time the grid is changed, e.g. by refinement or coarsening.
- unsigned long changeIndex_ = 0;
- };
-
- /// @}
-
-} // end namespace AMDiS
diff --git a/src/amdis/GridTransferManager.hpp b/src/amdis/GridTransferManager.hpp
deleted file mode 100644
index 45375be2805e81888131b0180d05855aecf6e273..0000000000000000000000000000000000000000
--- a/src/amdis/GridTransferManager.hpp
+++ /dev/null
@@ -1,106 +0,0 @@
-#pragma once
-
-#include <cstdint>
-#include <map>
-#include <memory>
-
-#include <amdis/AdaptionInterface.hpp>
-#include <amdis/GridTransfer.hpp>
-#include <amdis/common/ConcurrentCache.hpp>
-
-namespace AMDiS
-{
- /**
- * \addtogroup Adaption
- * @{
- **/
-
- /// Static administration class for automatic handling of DOFVectors during grid adaption
- class GridTransferManager
- {
- using Self = GridTransferManager;
- using Key = std::uintptr_t;
- using Data = std::shared_ptr<AdaptionInterface>;
-
- private:
- // Constructors. Since this is a static class you cannot call any constructor.
- GridTransferManager() = delete;
-
- public:
- /**
- * \brief Adapts the grid according to previously set marks and performs data transfer on
- * all DOFVectors associated to the given grid.
- *
- * This function retrieves the grid's stored GridTransfer and calls its adaptation methods.
- * During this process the grid will change according to previously set grid element markings.
- * All DOFVectors associated to the given grid will be interpolated to the new grid according to
- * their DataTransfer member object.
- *
- * Takes a grid as argument. Marking of the grid needs to be performed prior to calling this.
- * Returns true if the grid changed during adaptation.
- **/
- template <class Grid>
- static bool adapt(Grid& grid)
- {
- auto transfer = gridTransfer(grid);
- bool adapted = transfer->preAdapt();
- adapted |= transfer->adapt();
- transfer->postAdapt();
- return adapted;
- }
-
- /// Returns the GridTransfer bound to a Grid, to be used during the adaptation cycle
- template <class Grid>
- static std::shared_ptr<GridTransfer<Grid>> gridTransfer(Grid& grid)
- {
- auto transfer = get(grid);
- transfer->bind(grid);
- return transfer;
- }
-
- public:
- template <class Grid, class Data>
- static void attach(Grid const& grid, Data& data)
- {
- get(grid)->attach(data);
- }
-
- template <class Grid, class Data, class Callback>
- static void attach(Grid const& grid, Data const& data, Callback&& callback)
- {
- get(grid)->attach(data, FWD(callback));
- }
-
- template <class Grid, class Data>
- static void detach(Grid const& grid, Data const& data)
- {
- get(grid)->detach(data);
- }
-
-
- /// \brief Returns the GridTransfer class, used for attaching and detaching Data
- /**
- * NOTE: The returned GridTransfer is not yet bound to a grid, since a
- * mutable reference `Grid&` is needed for the binding.
- **/
- template <class Grid>
- static std::shared_ptr<GridTransfer<Grid>> get(Grid const& grid)
- {
- Key key = Key(&grid);
- GridTransferCache cache;
- auto& gridTransferInterface = cache.get(key, [&](Key const&)
- {
- return std::make_unique<GridTransfer<Grid>>();
- });
-
- return std::dynamic_pointer_cast<GridTransfer<Grid>>(gridTransferInterface);
- }
-
- private:
- // Associate to each Grid (represented as address) a GridTransfer
- using GridTransferCache = ConcurrentCache< Key, Data, StaticLockedPolicy, std::map<Key, Data> >;
- };
-
- /// @}
-
-} // end namespace AMDiS
diff --git a/src/amdis/LinearForm.hpp b/src/amdis/LinearForm.hpp
index a57ec0f1b0790cbee3bbd57f2382b80dea20bbbe..8601d9047ef01e60dac06b12b895e75f83d3906a 100644
--- a/src/amdis/LinearForm.hpp
+++ b/src/amdis/LinearForm.hpp
@@ -1,5 +1,7 @@
#pragma once
+#include <dune/common/typeutilities.hh>
+
#include <amdis/LinearAlgebra.hpp>
#include <amdis/OperatorList.hpp>
#include <amdis/common/FlatVector.hpp>
@@ -25,7 +27,6 @@ namespace AMDiS
public:
using Backend = VectorBackend<BackendTraits<GB,T>>;
- using Comm = typename Backend::Traits::Comm;
/// The type of the functionspace basis associated to this linearform
using GlobalBasis = GB;
@@ -42,12 +43,10 @@ namespace AMDiS
public:
/// Constructor. Stores the shared_ptr of the basis and creates a new DataTransfer.
- /// The additional arguments can be the communication object. If not provided,
- /// a new communication is created.
- template <class GB_, class... Args,
- Dune::disableCopyMove<Self, GB_, Args...> = 0>
- explicit LinearForm(GB_&& basis, Args&&... args)
- : Super(FWD(basis), FWD(args)...)
+ template <class GB_,
+ Dune::disableCopyMove<Self, GB_> = 0>
+ explicit LinearForm(GB_&& basis)
+ : Super(FWD(basis))
{
operators_.init(*Super::basis());
}
@@ -117,16 +116,16 @@ namespace AMDiS
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
- template <class GB, class... Args>
- LinearForm(GB&& basis, Args&&... args)
- -> LinearForm<Underlying_t<GB>, double>;
+ template <class GB>
+ LinearForm(GB&& basis)
+ -> LinearForm<Underlying_t<GB>>;
#endif
- template <class T = double, class GB, class... Args>
- auto makeLinearForm(GB&& basis, Args&&... args)
+ template <class T = double, class GB>
+ auto makeLinearForm(GB&& basis)
{
using LF = LinearForm<Underlying_t<GB>, T>;
- return LF{FWD(basis), FWD(args)...};
+ return LF{FWD(basis)};
}
} // end namespace AMDiS
diff --git a/src/amdis/MeshCreator.hpp b/src/amdis/MeshCreator.hpp
index 41d1d0b9c6a5a7d350865a011ea2713a0f48e330..743ec9e01896399c720c618826ac5fe3ba6f40aa 100644
--- a/src/amdis/MeshCreator.hpp
+++ b/src/amdis/MeshCreator.hpp
@@ -15,6 +15,7 @@
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
+#include <amdis/AdaptiveGrid.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/Concepts.hpp>
@@ -63,27 +64,28 @@ namespace AMDiS
*
* Otherwise tries to create a grid depending on the grid type.
**/
- std::unique_ptr<Grid> create() const
+ std::shared_ptr<AdaptiveGrid<Grid>> create() const
{
auto filename = Parameters::get<std::string>(name_ + "->macro file name");
auto structured = Parameters::get<std::string>(name_ + "->structured");
+ std::shared_ptr<Grid> gridPtr;
if (filename) {
// read a macro file
- return create_unstructured_grid(filename.value());
+ gridPtr = create_unstructured_grid(filename.value());
} else if (structured) {
// create structured grids
if (structured.value() == "cube") {
- return create_cube_grid();
+ gridPtr = create_cube_grid();
} else if (structured && structured.value() == "simplex") {
- return create_simplex_grid();
+ gridPtr = create_simplex_grid();
} else {
error_exit("Unknown structured grid type. Must be either `cube` or `simplex` in parameter [meshName]->structured.");
- return {};
}
} else {
// decide by inspecting the grid type how to create the grid
- return create_by_gridtype<Grid>(Dune::PriorityTag<42>{});
+ gridPtr = create_by_gridtype<Grid>(Dune::PriorityTag<42>{});
}
+ return AdaptiveGrid<Grid>::instance(gridPtr);
}
/// Create a structured cube grid
diff --git a/src/amdis/Observer.hpp b/src/amdis/Observer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..62383a9f5fe7a5737abd733bffffd117388db050
--- /dev/null
+++ b/src/amdis/Observer.hpp
@@ -0,0 +1,232 @@
+#pragma once
+
+#include <algorithm>
+#include <list>
+#include <memory>
+#include <type_traits>
+
+#include <amdis/common/ConceptsBase.hpp>
+#include <amdis/Output.hpp>
+
+namespace AMDiS
+{
+ namespace Impl
+ {
+ // Forward declaration
+ template <class Event>
+ class SignalInterface;
+
+ template <class Event>
+ class ObserverInterface
+ {
+ public:
+ virtual ~ObserverInterface() = default;
+
+ /// Attach the observer to a subject. It will then receive notifications from the subject when
+ /// an event of type Event it triggered.
+ void attach(std::shared_ptr<SignalInterface<Event>> const& subject)
+ {
+ subject->attachObserver(this);
+ }
+
+ /// Detach the observer from the subject. It will then no longer receive notifications. This
+ /// must be called before the observer is deleted.
+ void detach(std::shared_ptr<SignalInterface<Event>> const& subject)
+ {
+ subject->detachObserver(this);
+ }
+
+ /// This method will be called by the subject when triggering the event.
+ virtual void update(Event const&)
+ {
+ error_exit("Method must be overridden by derived class");
+ }
+ };
+
+ template <class Event>
+ class SignalInterface
+ {
+ public:
+ virtual ~SignalInterface() = default;
+
+ /// Attaches an observer to this class. This method will be called by all observers with
+ /// with themselves as argument.
+ virtual void attachObserver(ObserverInterface<Event>*) = 0;
+
+ /// Detaches an observer to this class. This method will be called by all observers with
+ /// with themselves as argument.
+ virtual void detachObserver(ObserverInterface<Event>*) = 0;
+
+ /// Notify all observers that have called attachObserver but have not called detachObserver
+ virtual void notify(Event const&) const = 0;
+ };
+
+ template <class Event>
+ class SignalsImpl
+ : virtual public SignalInterface<Event>
+ {
+ public:
+ virtual ~SignalsImpl() = default;
+
+ void attachObserver(ObserverInterface<Event>* o) override
+ {
+ observers_.push_back(o);
+ }
+
+ void detachObserver(ObserverInterface<Event>* o) override
+ {
+ auto it = std::find(observers_.begin(), observers_.end(), o);
+ if (it != observers_.end())
+ observers_.erase(it);
+ }
+
+ void notify(Event const& e) const override
+ {
+ for (ObserverInterface<Event>* o : observers_)
+ o->update(e);
+ }
+
+ private:
+ std::list<ObserverInterface<Event>*> observers_;
+ };
+
+ template <class Event>
+ class SignalsDummy
+ : virtual public SignalInterface<Event>
+ {
+ public:
+ virtual ~SignalsDummy() = default;
+ void attachObserver(ObserverInterface<Event>* o) override {}
+ void detachObserver(ObserverInterface<Event>* o) override {}
+ void notify(Event const& e) const override {}
+ };
+
+ template <class Event>
+ class ObserverImpl
+ : virtual public ObserverInterface<Event>
+ {
+ using Self = ObserverImpl;
+
+ ObserverImpl()
+ : subject_(std::make_shared<SignalsDummy<Event>>())
+ {}
+
+ public:
+ ObserverImpl(std::shared_ptr<SignalInterface<Event>> subject)
+ : subject_(std::move(subject))
+ {
+ this->attach(subject_);
+ }
+
+ // Get next subject in hierarchy if Event is not handled
+ template <class Subject,
+ REQUIRES(!std::is_convertible<std::shared_ptr<Subject>, std::shared_ptr<SignalInterface<Event>>>::value)>
+ ObserverImpl(std::shared_ptr<Subject> subject)
+ : Self(subject->getSubject())
+ {}
+
+ /// Copy constructor
+ ObserverImpl(Self const& that)
+ : subject_(that.subject_)
+ {
+ this->attach(subject_);
+ }
+
+ /// Move constructor
+ ObserverImpl(Self&& that)
+ : ObserverImpl()
+ {
+ swap(*this, that);
+ }
+
+ /// Destructor. Detaches observer from the subject
+ virtual ~ObserverImpl()
+ {
+ this->detach(subject_);
+ }
+
+ /// Copy/Move assignment
+ Self& operator=(Self that)
+ {
+ swap(*this, that);
+ return *this;
+ }
+
+ /// Swaps the contents of lhs and rhs
+ friend void swap(Self& lhs, Self& rhs)
+ {
+ using std::swap;
+ lhs.detach(lhs.subject_);
+ rhs.detach(rhs.subject_);
+ swap(lhs.subject_, rhs.subject_);
+ lhs.attach(lhs.subject_);
+ rhs.attach(rhs.subject_);
+ }
+
+ private:
+ /// The observed subject
+ std::shared_ptr<SignalInterface<Event>> subject_;
+ };
+ }
+
+ /** Mixin for signaling of certain events.
+ * Derived classes T can emit a signal e by calling notify(e). This will send the signal to all
+ * classes U... deriving from Observer<S, Events...> if
+ * - the type of the event is included in Events,
+ * - T = S or S is an Observer of T (directly or indirectly via other Observers),
+ * - U called the Observer constructor with an instance of S that has direct or indirect access
+ * to the instance of T (see Observer)
+ */
+ template <class Event, class... Events>
+ class Signals
+ : public Impl::SignalsImpl<Event>
+ , public Signals<Events...>
+ {
+ public:
+ using Impl::SignalsImpl<Event>::notify;
+ using Signals<Events...>::notify;
+ };
+
+ template <class Event>
+ class Signals<Event>
+ : public Impl::SignalsImpl<Event>
+ {
+ public:
+ using Impl::SignalsImpl<Event>::notify;
+ };
+
+
+ /** Mixin for reacting to certain events
+ * Derived classes T can react to events by implementing the functions update(Event& e) for
+ * all Events specified in the template parameter list, whose origin is an instance of class
+ * Subject specified in the constructor to Observer and any indirectly observed class instances.
+ * Observed instances may include:
+ * - the instance of class Subject passed to the constructor
+ * - the instance passed to the constructor of Subject::Observer<S> if Subject inherits from
+ * Observer<S>
+ * - all other instances indirectly accesible in the above way
+ * For each event E the first object in the above hierarchy that implements Signals<Es> with E
+ * included in Es will be observed by this class.
+ */
+ template <class Subject, class... Events>
+ class Observer
+ : public Impl::ObserverImpl<Events>...
+ {
+ template <class E>
+ friend class Impl::ObserverImpl;
+
+ public:
+ Observer(std::shared_ptr<Subject> s)
+ : Impl::ObserverImpl<Events>(s)...
+ , subject_(std::move(s))
+ {}
+
+ private:
+ /// Get the direct subject of this instance. Used by the implementation to move upwards in the
+ /// observer hierarchy.
+ std::shared_ptr<Subject> const& getSubject() { return subject_; }
+
+ std::shared_ptr<Subject> subject_;
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 02796a2d423f235f280c0f9a13ea93bea1efb978..e44a627aa8cdd29ded5be325c862e205dd7ddf79 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -10,8 +10,12 @@
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
+#include <dune/common/shared_ptr.hh>
+
#include <dune/grid/common/grid.hh>
+
#include <amdis/AdaptInfo.hpp>
+#include <amdis/AdaptiveGrid.hpp>
#include <amdis/BiLinearForm.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/CreatorMap.hpp>
@@ -22,14 +26,15 @@
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/LinearForm.hpp>
-#include <amdis/OperatorList.hpp>
#include <amdis/Marker.hpp>
#include <amdis/MeshCreator.hpp>
+#include <amdis/OperatorList.hpp>
#include <amdis/PeriodicBC.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>
+#include <amdis/common/SharedPtr.hpp>
#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/TypeTraits.hpp>
@@ -61,7 +66,9 @@ namespace AMDiS
using GlobalBasis = typename Traits::GlobalBasis;
using GridView = typename GlobalBasis::GridView;
- using Grid = typename GridView::Grid;
+ // TODO(FM): Make this GridView::Grid::HostGrid as soon as AdaptiveGrid is a complete metagrid
+ using HostGrid = typename GridView::Grid;
+ using Grid = AdaptiveGrid<HostGrid>;
using Element = typename GridView::template Codim<0>::Entity;
using WorldVector = typename Element::Geometry::GlobalCoordinate;
using WorldMatrix = FieldMatrix<typename WorldVector::field_type, WorldVector::dimension, WorldVector::dimension>;
@@ -89,20 +96,22 @@ namespace AMDiS
: name_(name)
{}
- /// Constructor taking additionally a reference to a grid that is used
+ /// Constructor taking additionally a grid that is used
/// instead of the default created grid, \ref ProblemStat
- ProblemStat(std::string const& name, Grid& grid)
+ template <class Grid_>
+ ProblemStat(std::string const& name, Grid_&& grid)
: ProblemStat(name)
{
- adoptGrid(Dune::stackobject_to_shared_ptr(grid));
+ adoptGrid(Grid::instance(FWD(grid)));
}
- /// \brief Constructor taking a grid reference and a basis reference.
- /// Stores pointers to both.
- ProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
- : ProblemStat(name, grid)
+ /// \brief Constructor taking a grid and basis
+ /// Wraps both in shared pointers.
+ template <class Grid_, class GB>
+ ProblemStat(std::string const& name, Grid_&& grid, GB&& globalBasis)
+ : ProblemStat(name, FWD(grid))
{
- adoptGlobalBasis(Dune::stackobject_to_shared_ptr(globalBasis));
+ adoptGlobalBasis(wrap_or_share(FWD(globalBasis)));
}
@@ -312,12 +321,11 @@ namespace AMDiS
/// Implementation of \ref ProblemStatBase::name
std::string const& name() const override { return name_; }
-
- /// Return a reference to the grid, \ref grid
+ /// Return the \ref grid_
std::shared_ptr<Grid> grid() { return grid_; }
std::shared_ptr<Grid const> grid() const { return grid_; }
- /// Return the gridView of the leaf-level
+ /// Return the gridView of the basis
GridView gridView() const { return globalBasis_->gridView(); }
/// Return the boundary manager to identify boundary segments
@@ -328,9 +336,6 @@ namespace AMDiS
std::shared_ptr<GlobalBasis> globalBasis() { return globalBasis_; }
std::shared_ptr<GlobalBasis const> globalBasis() const { return globalBasis_; }
- std::shared_ptr<CommunicationType> communication() { return communication_; }
- std::shared_ptr<CommunicationType const> communication() const { return communication_; }
-
/// Return a reference to the linear solver, \ref linearSolver
std::shared_ptr<LinearSolverType> solver() { return linearSolver_; }
std::shared_ptr<LinearSolverType const> solver() const { return linearSolver_; }
@@ -385,9 +390,9 @@ namespace AMDiS
/// Set the grid. Stores pointer and initializes feSpaces
/// matrices and vectors, as well as markers and file-writers.
- void setGrid(std::shared_ptr<Grid> const& grid)
+ void setGrid(std::shared_ptr<Grid> grid)
{
- adoptGrid(grid);
+ adoptGrid(std::move(grid));
createGlobalBasis();
createMatricesAndVectors();
createMarker();
@@ -395,9 +400,10 @@ namespace AMDiS
}
/// Wrap the grid into a non-destroying shared_ptr
- void setGrid(Grid& grid)
+ template <class Grid_>
+ void setGrid(Grid_&& grid)
{
- setGrid(Dune::stackobject_to_shared_ptr(grid));
+ setGrid(Grid::instance(FWD(grid)));
}
@@ -444,13 +450,13 @@ namespace AMDiS
void createMarker();
void createFileWriter();
- void adoptGlobalBasis(std::shared_ptr<GlobalBasis> const& globalBasis)
+ void adoptGlobalBasis(std::shared_ptr<GlobalBasis> globalBasis)
{
- globalBasis_ = globalBasis;
+ globalBasis_ = std::move(globalBasis);
initGlobalBasis();
}
- void adoptGrid(std::shared_ptr<Grid> const& grid,
+ void adoptGrid(std::shared_ptr<Grid> grid,
std::shared_ptr<BoundaryManager<Grid>> const& boundaryManager)
{
grid_ = grid;
@@ -458,7 +464,7 @@ namespace AMDiS
Parameters::get(name_ + "->mesh", gridName_);
}
- void adoptGrid(std::shared_ptr<Grid> const& grid)
+ void adoptGrid(std::shared_ptr<Grid> grid)
{
adoptGrid(grid, std::make_shared<BoundaryManager<Grid>>(grid));
}
@@ -483,7 +489,7 @@ namespace AMDiS
/// Management of boundary conditions
std::shared_ptr<BoundaryManager<Grid>> boundaryManager_;
- /// FE spaces of this problem.
+ /// FE space of this problem.
std::shared_ptr<GlobalBasis> globalBasis_;
/// A FileWriter object
@@ -508,8 +514,6 @@ namespace AMDiS
/// of the equation, filled during assembling
std::shared_ptr<SystemVector> rhs_;
- std::shared_ptr<CommunicationType> communication_;
-
/// A vector with the local element error estimates
/// for each node in the basis tree, indexed by [to_string(treePath)][element index]
std::map<std::string, std::vector<double>> estimates_;
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index c6e7fe9f2aa126b29beeb27050d21f3cdc5db6e3..b12361d2613a95e197e7213bb4bdbd6f12f0ddc8 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -14,7 +14,6 @@
#include <amdis/BackupRestore.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/GridFunctionOperator.hpp>
-#include <amdis/GridTransferManager.hpp>
#include <amdis/io/FileWriterCreator.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
@@ -57,12 +56,10 @@ void ProblemStat<Traits>::initialize(
bool loadBalance = false;
Parameters::get(gridName_ + "->load balance", loadBalance);
if (loadBalance)
- grid_->loadBalance();
+ loadBalance = grid_->loadBalance();
- if (globalBasis_ && (globalRefinements > 0 || loadBalance)) {
+ if (globalBasis_ && (globalRefinements > 0 || loadBalance))
globalBasis_->update(globalBasis_->gridView());
- communication_->update(*globalBasis_);
- }
}
// create fespace
@@ -84,7 +81,6 @@ void ProblemStat<Traits>::initialize(
if (!globalBasis_)
warning("no globalBasis created\n");
-
// create system
if (initFlag.isSet(INIT_SYSTEM))
createMatricesAndVectors();
@@ -141,10 +137,10 @@ restore(Flag initFlag)
test_exit(filesystem::exists(solution_filename), "Restore file '{}' not found.", solution_filename);
// restore grid from file
- if (Dune::Capabilities::hasBackupRestoreFacilities<Grid>::v)
- adoptGrid(std::unique_ptr<Grid>(Dune::BackupRestoreFacility<Grid>::restore(grid_filename)));
+ if (Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v)
+ adoptGrid(Grid::instance(Dune::BackupRestoreFacility<HostGrid>::restore(grid_filename)));
else
- adoptGrid(std::unique_ptr<Grid>(BackupRestoreByGridFactory<Grid>::restore(grid_filename)));
+ adoptGrid(Grid::instance(BackupRestoreByGridFactory<HostGrid>::restore(grid_filename)));
// create fespace
if (initFlag.isSet(INIT_FE_SPACE) || initFlag.isSet(INIT_SYSTEM))
@@ -177,12 +173,13 @@ template <class Traits>
void ProblemStat<Traits>::createGrid()
{
Parameters::get(name_ + "->mesh", gridName_);
- MeshCreator<Grid> creator(gridName_);
+ MeshCreator<HostGrid> creator(gridName_);
grid_ = creator.create();
Dune::Timer t;
- grid_->loadBalance();
- info(2,"load balance needed {} seconds", t.elapsed());
+ bool loadBalance = grid_->loadBalance();
+ if (loadBalance)
+ info(2,"load balance needed {} seconds", t.elapsed());
boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
if (!creator.boundaryIds().empty())
@@ -214,8 +211,10 @@ template <class Traits>
void ProblemStat<Traits>::createGlobalBasisImpl(std::true_type)
{
assert( bool(grid_) );
- static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
- globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
+ // TODO(FM): Replace HostGrid -> Grid once AdaptiveGrid is full metagrid
+ static_assert(std::is_same<GridView, typename HostGrid::LeafGridView>::value, "");
+ auto basis = Traits::create(name_, grid_->leafGridView());
+ globalBasis_ = std::make_shared<GlobalBasis>(std::move(basis));
}
@@ -231,24 +230,15 @@ void ProblemStat<Traits>::initGlobalBasis()
{
dirichletBCs_.init(*globalBasis_, *globalBasis_);
periodicBCs_.init(*globalBasis_, *globalBasis_);
-
- using CommCreator = CommunicationCreator<CommunicationType>;
- communication_ = CommCreator::create(*globalBasis_, name_ + "->solver");
-
- GridTransferManager::attach(*grid_, *communication_, [c=communication_, b=globalBasis_]() -> void
- {
- b->update(b->gridView());
- c->update(*b);
- });
}
template <class Traits>
void ProblemStat<Traits>::createMatricesAndVectors()
{
- systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_, communication_);
- solution_ = std::make_shared<CoefficientVector>(globalBasis_, communication_);
- rhs_ = std::make_shared<SystemVector>(globalBasis_, communication_);
+ systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
+ solution_ = std::make_shared<CoefficientVector>(globalBasis_);
+ rhs_ = std::make_shared<SystemVector>(globalBasis_);
auto localView = globalBasis_->localView();
for_each_node(localView.tree(), [&,this](auto const& node, auto treePath) -> void
@@ -338,7 +328,7 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val
auto i = child(localView.tree(), makeTreePath(row));
auto j = child(localView.tree(), makeTreePath(col));
- auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
+ auto valueGridFct = makeGridFunction(values, this->gridView());
using Range = RangeType_t<decltype(i)>;
using BcType = DirichletBC<WorldVector,Range>;
@@ -357,7 +347,7 @@ addDirichletBC(BoundaryType id, RowTreePath row, ColTreePath col, Values const&
auto i = child(localView.tree(), makeTreePath(row));
auto j = child(localView.tree(), makeTreePath(col));
- auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
+ auto valueGridFct = makeGridFunction(values, this->gridView());
using Range = RangeType_t<decltype(i)>;
using BcType = DirichletBC<WorldVector,Range>;
@@ -389,8 +379,8 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
solverInfo.setStoreMatrixData(storeMatrixData);
solution_->resize();
- linearSolver_->solve(systemMatrix_->backend().matrix(), solution_->backend().vector(), rhs_->backend().vector(),
- *communication_, solverInfo);
+ linearSolver_->solve(systemMatrix_->backend().matrix(), solution_->backend().vector(),
+ rhs_->backend().vector(), globalBasis_->communicator(), solverInfo);
if (solverInfo.info() > 0) {
msg("solution of discrete system needed {} seconds", t.elapsed());
@@ -439,12 +429,19 @@ globalCoarsen(int n)
{
Dune::Timer t;
bool adapted = false;
+ // TODO(FM): Find a less expensive alternative to the loop adaption
for (int i = 0; i < n; ++i) {
- // mark all entities for grid refinement
+ // mark all entities for coarsening
for (const auto& element : elements(grid_->leafGridView()))
grid_->mark(-1, element);
- adapted |= GridTransferManager::adapt(*grid_);
+ bool adaptedInLoop = grid_->preAdapt();
+ adaptedInLoop |= grid_->adapt();
+ grid_->postAdapt();
+ if (!adaptedInLoop)
+ break;
+ else
+ adapted = true;
}
msg("globalCoarsen needed {} seconds", t.elapsed());
@@ -454,31 +451,24 @@ globalCoarsen(int n)
// grid has globalRefine(int, AdaptDataHandleInterface&)
template <class G>
-using HasGlobalRefineADHI = decltype(std::declval<G>().globalRefine(1,std::declval<GridTransfer<G>&>()));
+using HasGlobalRefineADHI = decltype(
+ std::declval<G>().globalRefine(1,std::declval<typename G::ADHI&>()));
template <class Traits>
Flag ProblemStat<Traits>::
-globalRefine(int refCount)
+globalRefine(int n)
{
Dune::Timer t;
- bool adapted = false;
Dune::Hybrid::ifElse(Dune::Std::is_detected<HasGlobalRefineADHI, Grid>{},
/*then*/ [&](auto id) -> void {
- // TODO(FM): Add a way to pass a GridTransfer as ADH with *globalBasis_ argument
- id(grid_)->globalRefine(refCount, GridTransferManager::gridTransfer(*grid_));
+ id(grid_)->globalRefine(n, globalBasis_->globalRefineCallback());
},
/*else*/ [&](auto id) -> void {
- for (int i = 0; i < refCount; ++i) {
- // mark all entities for grid refinement
- for (const auto& element : elements(grid_->leafGridView()))
- grid_->mark(1, element);
-
- adapted |= GridTransferManager::adapt(*id(grid_));
- }
+ id(grid_)->globalRefine(n);
});
msg("globalRefine needed {} seconds", t.elapsed());
- return adapted ? MESH_ADAPTED : Flag(0);
+ return n > 0 ? MESH_ADAPTED : Flag(0);
}
@@ -488,7 +478,9 @@ adaptGrid(AdaptInfo& adaptInfo)
{
Dune::Timer t;
- bool adapted = GridTransferManager::adapt(*grid_);
+ bool adapted = grid_->preAdapt();
+ adapted |= grid_->adapt();
+ grid_->postAdapt();
msg("adaptGrid needed {} seconds", t.elapsed());
return adapted ? MESH_ADAPTED : Flag(0);
diff --git a/src/amdis/ProblemStatTraits.hpp b/src/amdis/ProblemStatTraits.hpp
index befa48b5c9bc13f18fb810d9dd31017cfe313253..87393b725798540abd9fd9e02645d196588e86dc 100644
--- a/src/amdis/ProblemStatTraits.hpp
+++ b/src/amdis/ProblemStatTraits.hpp
@@ -7,6 +7,8 @@
#include <dune/grid/yaspgrid.hh>
#include <amdis/common/Logical.hpp>
+#include <amdis/common/TypeTraits.hpp>
+#include <amdis/functions/ParallelGlobalBasis.hpp>
namespace AMDiS
{
@@ -69,13 +71,17 @@ namespace AMDiS
struct DefaultBasisCreator
{
using GridView = typename Grid::LeafGridView;
+ static auto create(std::string const& name, GridView const& gridView)
+ {
+ return makeGlobalBasis(name, gridView, PreBasisCreator::create());
+ }
+
static auto create(GridView const& gridView)
{
- using namespace Dune::Functions::BasisFactory;
- return makeBasis(gridView, PreBasisCreator::create());
+ return makeGlobalBasis(gridView, PreBasisCreator::create());
}
- using GlobalBasis = decltype(create(std::declval<GridView>()));
+ using GlobalBasis = Underlying_t<decltype(create(std::declval<GridView>()))>;
using CoefficientType = T;
};
diff --git a/src/amdis/common/CMakeLists.txt b/src/amdis/common/CMakeLists.txt
index f14c1c2468855aab773327cc6784d7b3a23de38c..32b15cb34d62007698ba72b1844924f8b0c2767d 100644
--- a/src/amdis/common/CMakeLists.txt
+++ b/src/amdis/common/CMakeLists.txt
@@ -28,6 +28,7 @@ install(FILES
Range.hpp
Resize.hpp
QuadMath.hpp
+ SharedPtr.hpp
StaticSize.hpp
String.hpp
SwitchCases.hpp
@@ -38,4 +39,4 @@ install(FILES
ValueCategory.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/common)
-add_subdirectory(parallel)
\ No newline at end of file
+add_subdirectory(parallel)
diff --git a/src/amdis/common/SharedPtr.hpp b/src/amdis/common/SharedPtr.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3cb9bffd796263f45007da801594ba893297657b
--- /dev/null
+++ b/src/amdis/common/SharedPtr.hpp
@@ -0,0 +1,62 @@
+#pragma once
+
+#include <memory>
+#include <utility>
+
+#include <dune/common/shared_ptr.hh>
+
+#include <amdis/common/TypeTraits.hpp>
+
+namespace AMDiS
+{
+ /**
+ * Helper function converting the argument into a shared_ptr. The following cases are handled:
+ * - objects passed as rvalue references are used to construct a new shared_ptr
+ * - objects passed as references are used to construct a shared_ptr with a null deleter that
+ * does nothing when the pointer is destroyed, using the stored base pointer in case of a
+ * unique_ptr argument
+ * - shared_ptr are simply forwarded
+ */
+ template <class T>
+ std::shared_ptr<T> wrap_or_share(T& t)
+ {
+ return Dune::stackobject_to_shared_ptr(t);
+ }
+
+ template<class T>
+ std::shared_ptr<T> wrap_or_share(T&& t)
+ {
+ return std::make_shared<T>(FWD(t));
+ }
+
+ template <class T>
+ std::shared_ptr<T> wrap_or_share(T*& t)
+ {
+ return std::shared_ptr<T>(t, Dune::null_deleter<T>());
+ }
+
+ template <class T>
+ std::shared_ptr<T> wrap_or_share(T*&& t)
+ {
+ return std::shared_ptr<T>(t);
+ }
+
+ template <class T>
+ std::shared_ptr<T> wrap_or_share(std::shared_ptr<T> t)
+ {
+ return std::move(t);
+ }
+
+ template <class T>
+ std::shared_ptr<T> wrap_or_share(std::unique_ptr<T>& t)
+ {
+ return std::shared_ptr<T>(t.get(), Dune::null_deleter<T>());
+ }
+
+ template <class T>
+ std::shared_ptr<T> wrap_or_share(std::unique_ptr<T>&& t)
+ {
+ return std::shared_ptr<T>(std::move(t));
+ }
+
+} // end namespace AMDiS
diff --git a/src/amdis/functions/CMakeLists.txt b/src/amdis/functions/CMakeLists.txt
index 4dd010be0334767552c3423aa93ffa4ff432075c..73697b6352cfe82f81c4258f7bf7ecc4ea1ef8af 100644
--- a/src/amdis/functions/CMakeLists.txt
+++ b/src/amdis/functions/CMakeLists.txt
@@ -5,4 +5,5 @@ install(FILES
Interpolate.hpp
NodeIndices.hpp
Nodes.hpp
+ ParallelGlobalBasis.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions)
diff --git a/src/amdis/functions/ParallelGlobalBasis.hpp b/src/amdis/functions/ParallelGlobalBasis.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..81ca2ab80ec301e6f68da710681a33f69089bc26
--- /dev/null
+++ b/src/amdis/functions/ParallelGlobalBasis.hpp
@@ -0,0 +1,293 @@
+#pragma once
+
+
+#include <algorithm>
+#include <list>
+#include <memory>
+#include <type_traits>
+#include <utility>
+
+#include <dune/common/concept.hh>
+#include <dune/common/reservedvector.hh>
+#include <dune/common/shared_ptr.hh>
+#include <dune/common/typeutilities.hh>
+#include <dune/common/version.hh>
+
+#include <dune/functions/common/type_traits.hh>
+#include <dune/functions/functionspacebases/concepts.hh>
+#include <dune/functions/functionspacebases/defaultlocalview.hh>
+#include <dune/functions/functionspacebases/flatmultiindex.hh>
+
+#include <dune/grid/common/adaptcallback.hh>
+
+#include <dune/typetree/treepath.hh>
+
+#include <amdis/common/Concepts.hpp>
+#include <amdis/common/TypeTraits.hpp>
+#include <amdis/linearalgebra/Traits.hpp>
+#include <amdis/AdaptiveGrid.hpp>
+#include <amdis/Observer.hpp>
+#include <amdis/Output.hpp>
+
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+#include <dune/functions/functionspacebases/defaultlocalindexset.hh>
+#endif
+
+namespace Dune
+{
+ namespace Functions
+ {
+ // Forward declaration
+ template <class PB>
+ class DefaultGlobalBasis;
+ }
+}
+
+namespace AMDiS
+{
+ namespace event
+ {
+ template <class B>
+ struct basisUpdate
+ {
+ basisUpdate(B const& b_)
+ : basis(b_)
+ {}
+
+ B const& basis;
+ };
+ }
+
+ /**
+ * \brief Parallel global basis defined on a (sequential) pre-basis
+ *
+ * This class is an expansion to Dune::Functions::DefaultGlobalBasis<PB>. It adds a communicator
+ * to use the basis in parallel as well as automatic update functionality.
+ *
+ * \tparam PB Pre-basis providing the implementation details
+ */
+ template <class PB>
+ class ParallelGlobalBasis
+ : public Observer<AdaptiveGrid<typename PB::GridView::Grid>, event::adapt>
+ , public Signals<event::basisUpdate<ParallelGlobalBasis<PB>>>
+ {
+ using Subject = AdaptiveGrid<typename PB::GridView::Grid>;
+ using Obs = Observer<Subject, event::adapt>;
+ using Self = ParallelGlobalBasis<PB>;
+
+ struct DummyImpl {};
+
+ public:
+
+ /// Pre-basis providing the implementation details
+ using PreBasis = PB;
+
+ /// The empty prefix path that identifies the root in the local ansatz tree
+ using PrefixPath = Dune::TypeTree::HybridTreePath<>;
+
+ /// The grid view that the FE space is defined on
+ using GridView = typename PreBasis::GridView;
+ // using Grid = typename GridView::Grid; // TODO(FM): Use once AdaptiveGrid is complete Metagrid
+ using Grid = Subject;
+
+ /// Type used for global numbering of the basis vectors
+ using MultiIndex = typename PreBasis::MultiIndex;
+
+ /// Type used for indices and size information
+ using size_type = std::size_t;
+
+ /// Type of the local view on the restriction of the basis to a single element
+ using LocalView = Dune::Functions::DefaultLocalView<Self>;
+
+ /// Type used for prefixes handed to the size() method
+ using SizePrefix = typename PreBasis::SizePrefix;
+
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+ /// Node index set provided by PreBasis
+ using NodeIndexSet = typename PreBasis::template IndexSet<PrefixPath>;
+
+ /// Type of local indixes set exported by localIndexSet()
+ using LocalIndexSet = Dune::Functions::DefaultLocalIndexSet<LocalView, NodeIndexSet>;
+#endif
+
+ /// Type of the communicator
+ using Comm = typename BackendTraits<Self>::Comm;
+
+ using ADH = Dune::AdaptDataHandle<Grid, DummyImpl>;
+
+ /**
+ * \brief Constructor
+ *
+ * \tparam T Argument list for PreBasis
+ * \param s The AdaptiveGrid providing the GridView for this basis
+ * \param t Argument list for PreBasis
+ *
+ * This will forward all arguments to the constructor of PreBasis
+ */
+ template<class... T,
+ Dune::Functions::enableIfConstructible<PreBasis, T...> = 0>
+ ParallelGlobalBasis(std::string const& name, std::shared_ptr<Subject> s, T&&... t)
+ : Obs(std::move(s))
+ , preBasis_(FWD(t)...)
+ , prefixPath_()
+ {
+ static_assert(Dune::models<Dune::Functions::Concept::PreBasis<GridView>, PreBasis>(),
+ "Type passed to ParallelGlobalBasis does not model the PreBasis concept.");
+ preBasis_.initializeIndices();
+ comm_ = CommunicationCreator<Comm>::create(*this, name + "->solver");
+ }
+
+ template <class... T>
+ ParallelGlobalBasis(std::string const& name, Subject& s, T&&... t)
+ : ParallelGlobalBasis(name, Dune::stackobject_to_shared_ptr(s), FWD(t)...)
+ {}
+
+ template <class S, class... T,
+ REQUIRES(Concepts::Same<Underlying_t<S>, Subject>)>
+ ParallelGlobalBasis(S&& s, T&&... t)
+ : ParallelGlobalBasis("", FWD(s), FWD(t)...)
+ {}
+
+ /** Converting constructor from Dune::DefaultGlobalBasis<PB>.
+ * This will create a new ParallellGlobalBasis. The pre-basis is copied from the constructor
+ * argument and a new communication object is built.
+ */
+ // TODO(FM): Replace explicit type with concept
+ ParallelGlobalBasis(Dune::Functions::DefaultGlobalBasis<PB>&& from)
+ : Obs(Subject::instance(from.gridView().grid()))
+ , preBasis_(from.preBasis())
+ , prefixPath_(from.prefixPath())
+ {
+ comm_ = CommunicationCreator<Comm>::create(*this);
+ }
+
+ /// Obtain the grid view that the basis is defined on
+ GridView const& gridView() const
+ {
+ return preBasis_.gridView();
+ }
+
+ /// Obtain the pre-basis providing the implementation details
+ PreBasis const& preBasis() const
+ {
+ return preBasis_;
+ }
+
+ /// Obtain the pre-basis providing the implementation details
+ PreBasis& preBasis()
+ {
+ return preBasis_;
+ }
+
+ /// Updates the underlying basis when event::adapt is triggered by the observed grid
+ void update(event::adapt const& e) override
+ {
+ if (e.adapted)
+ {
+ update(gridView());
+ event::basisUpdate<Self> eNew{*this};
+ this->notify(eNew);
+ }
+ }
+
+ /**
+ * \brief Update the stored grid view
+ *
+ * This will update the indexing information of the global basis as well as the communicator.
+ * It is called automatically if the grid has changed.
+ */
+ void update(GridView const& gv)
+ {
+ preBasis_.update(gv);
+ preBasis_.initializeIndices();
+ comm_->update(*this);
+ }
+
+ /// Get the total dimension of the space spanned by this basis
+ size_type dimension() const
+ {
+ return preBasis_.dimension();
+ }
+
+ /// Return number of possible values for next position in empty multi index
+ size_type size() const
+ {
+ return preBasis_.size();
+ }
+
+ /// Return number of possible values for next position in multi index
+ size_type size(const SizePrefix& prefix) const
+ {
+ return preBasis_.size(prefix);
+ }
+
+ /// Return local view for basis
+ LocalView localView() const
+ {
+ return LocalView(*this);
+ }
+
+ /// Return *this because we are not embedded in a larger basis
+ ParallelGlobalBasis const& rootBasis() const
+ {
+ return *this;
+ }
+
+ /// Return empty path, because this is the root in the local ansatz tree
+ PrefixPath const& prefixPath() const
+ {
+ return prefixPath_;
+ }
+
+ /// Return the communicator.
+ /// This provides the means to communicate data associated to the basis with other processes.
+ Comm const& communicator() const { return *comm_; }
+ Comm& communicator() { return *comm_; }
+
+ ADH& globalRefineCallback() const
+ {
+ // TODO(FM): Implement
+ error_exit("Not implemented: ParallelGlobalBasis::globalRefineCallback()");
+ return ADH{};
+ }
+
+ protected:
+ PreBasis preBasis_;
+ PrefixPath prefixPath_;
+ std::unique_ptr<Comm> comm_;
+ };
+
+ template<class MultiIndex, class GridView, class PreBasisFactory>
+ auto makeGlobalBasis(std::string const& name, GridView const& gridView,
+ PreBasisFactory&& preBasisFactory)
+ {
+ auto preBasis = preBasisFactory.template makePreBasis<MultiIndex>(gridView);
+ using PreBasis = std::decay_t<decltype(preBasis)>;
+
+ using Grid = typename GridView::Grid;
+ auto ag = AdaptiveGrid<Grid>::instance(gridView.grid());
+
+ return ParallelGlobalBasis<PreBasis>(name, std::move(ag), std::move(preBasis));
+ }
+
+ template<class GridView, class PreBasisFactory>
+ auto makeGlobalBasis(std::string const& name, GridView const& gridView,
+ PreBasisFactory&& preBasisFactory)
+ {
+ using RawPreBasisFactory = std::decay_t<PreBasisFactory>;
+ using MultiIndex = std::conditional_t<
+ (RawPreBasisFactory::requiredMultiIndexSize == 1),
+ Dune::Functions::FlatMultiIndex<std::size_t>,
+ Dune::ReservedVector<std::size_t, RawPreBasisFactory::requiredMultiIndexSize>>;
+
+ return makeGlobalBasis<MultiIndex, GridView, PreBasisFactory>(name, gridView,
+ FWD(preBasisFactory));
+ }
+
+ template<class GridView, class PreBasisFactory>
+ auto makeGlobalBasis(GridView const& gridView, PreBasisFactory&& preBasisFactory)
+ {
+ return makeGlobalBasis("", gridView, FWD(preBasisFactory));
+ }
+
+} // end namespace AMDiS
diff --git a/src/amdis/io/FileWriterCreator.hpp b/src/amdis/io/FileWriterCreator.hpp
index 1193f19b924e2c56e835c715592cd26212666fe1..52f62913a21c5014bd55608d69de465409c1db73 100644
--- a/src/amdis/io/FileWriterCreator.hpp
+++ b/src/amdis/io/FileWriterCreator.hpp
@@ -24,7 +24,8 @@ namespace AMDiS
{
using GlobalBasis = typename SystemVector::GlobalBasis;
using GridView = typename GlobalBasis::GridView;
- using Grid = typename GridView::Grid;
+ using Grid = typename GlobalBasis::Grid;
+ //using Grid = typename GridView::Grid; // TODO(FM): Use once AdaptiveGrid is complete Metagrid
public:
/// Constructor. Stores the pointer to the systemVector and to the (optional) boundaryManager.
diff --git a/src/amdis/linearalgebra/Communication.hpp b/src/amdis/linearalgebra/Communication.hpp
index ba5612e2ba5a788bc88ad957e7aad9da7780a65b..e08092386de8b16db09c0842f59116d0887313bb 100644
--- a/src/amdis/linearalgebra/Communication.hpp
+++ b/src/amdis/linearalgebra/Communication.hpp
@@ -16,8 +16,8 @@ namespace AMDiS
public:
using Impl = SequentialCommunication;
- template <class Comm>
- SequentialCommunication(Comm&&) {}
+ template <class T>
+ SequentialCommunication(T&&) {}
Impl const& get() const
{
diff --git a/src/amdis/linearalgebra/MatrixBase.hpp b/src/amdis/linearalgebra/MatrixBase.hpp
index 601330327a3f88fd5b8854ad1013bbbcb7308999..38481bda4705124bcb43500dbf4ed08e3422a28d 100644
--- a/src/amdis/linearalgebra/MatrixBase.hpp
+++ b/src/amdis/linearalgebra/MatrixBase.hpp
@@ -38,34 +38,16 @@ namespace AMDiS
/// The Linear-Algebra backend used to store the assembled coefficients
using Backend = B;
- using Comm = typename B::Traits::Comm;
public:
/// (1) Constructor. Stores the shared_ptr to the bases and the communication object.
- MatrixBase(std::shared_ptr<RB> rowBasis, std::shared_ptr<CB> colBasis,
- std::shared_ptr<Comm> comm)
+ MatrixBase(std::shared_ptr<RB> rowBasis, std::shared_ptr<CB> colBasis)
: rowBasis_(std::move(rowBasis))
, colBasis_(std::move(colBasis))
- , backend_(std::move(comm))
+ , backend_(*rowBasis_, *colBasis_)
{}
- /// (2) Constructor. Forwards to (1) by wraping reference into non-destroying
- /// shared_ptr, see \ref Dune::wrap_or_move.
- template <class RB_, class CB_, class Comm_,
- REQUIRES(Concepts::Similar<Types<RB_,CB_,Comm_>, Types<RB,CB,Comm>>)>
- MatrixBase(RB_&& rowBasis, CB_&& colBasis, Comm_&& comm)
- : MatrixBase(Dune::wrap_or_move(FWD(rowBasis)),
- Dune::wrap_or_move(FWD(colBasis)),
- Dune::wrap_or_move(FWD(comm)))
- {}
-
- /// (3) Constructor. Forwards to (1) by creating a new communicator
- MatrixBase(std::shared_ptr<RB> const& rowBasis, std::shared_ptr<CB> const& colBasis)
- : MatrixBase(rowBasis, colBasis,
- std::shared_ptr<Comm>(CommunicationCreator<Comm>::create(*rowBasis)))
- {}
-
- /// (4) Constructor. Forwards to (3) by wrapping references into non-destroying
+ /// (2) Constructor. Forwards to (1) by wrapping reference into non-destroying
/// shared_ptr, see \ref Dune::wrap_or_move.
template <class RB_, class CB_,
REQUIRES(Concepts::Similar<Types<RB_,CB_>, Types<RB,CB>>)>
diff --git a/src/amdis/linearalgebra/Traits.hpp b/src/amdis/linearalgebra/Traits.hpp
index 9f858ff210289820a86db56bfc9aedb3c631f771..30f08052c2bf61bd583a9ab464e82dc5d68d86d1 100644
--- a/src/amdis/linearalgebra/Traits.hpp
+++ b/src/amdis/linearalgebra/Traits.hpp
@@ -1,5 +1,18 @@
#pragma once
+#if HAVE_MTL
+#include <amdis/linearalgebra/mtl/Traits.hpp>
+
+#elif HAVE_EIGEN
+#include <amdis/linearalgebra/eigen/Traits.hpp>
+
+#elif HAVE_PETSC
+#include <amdis/linearalgebra/petsc/Traits.hpp>
+
+#else // ISTL
+#include <amdis/linearalgebra/istl/Traits.hpp>
+#endif
+
namespace AMDiS
{
#ifdef DOXYGEN
@@ -23,5 +36,4 @@ namespace AMDiS
using PartitionSet = Dune::Partitions::All; //< The dune partition set where to assemble operators
};
#endif
-
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/VectorBase.hpp b/src/amdis/linearalgebra/VectorBase.hpp
index 36a55ad16480943166b6f967d0f69b7512a5938b..25fc10edfe6b7a1b4d8fe92e6768f1b996749e99 100644
--- a/src/amdis/linearalgebra/VectorBase.hpp
+++ b/src/amdis/linearalgebra/VectorBase.hpp
@@ -44,7 +44,6 @@ namespace AMDiS
/// The Linear-Algebra backend used to store the assembled coefficients
using Backend = B;
- using Comm = typename B::Traits::Comm;
private:
// proxy class redirecting the mutable element access to the insertValue method
@@ -65,33 +64,18 @@ namespace AMDiS
std::integral_constant<VectorState, VectorState::insert_values>>;
public:
- /// (1) Constructor. Stores the shared_ptr of the basis and communication object.
- VectorBase(std::shared_ptr<GB> basis, std::shared_ptr<Comm> comm)
+ /// (1) Constructor. Stores the shared_ptr of the basis.
+ VectorBase(std::shared_ptr<GB> basis)
: basis_(std::move(basis))
- , comm_(std::move(comm))
- , backend_(comm_)
+ , backend_(*basis_)
{
resizeZero();
}
- /// (2) Constructor. Forwards to (1) by wraping reference into non-destroying
- /// shared_ptr.
- template <class GB_, class Comm_,
- REQUIRES(Concepts::Similar<Types<GB_,Comm_>, Types<GB,Comm>>)>
- VectorBase(GB_&& basis, Comm_&& comm)
- : VectorBase(Dune::wrap_or_move(FWD(basis)), Dune::wrap_or_move(FWD(comm)))
- {}
-
- /// (3) Constructor. Forwards to (1) by creating a new communicator
- explicit VectorBase(std::shared_ptr<GB> const& basis)
- : VectorBase(basis, std::shared_ptr<Comm>(CommunicationCreator<Comm>::create(*basis)))
- {}
-
- /// (4) Constructor. Forwards to (3) by wrapping references into non-destroying
- /// shared_ptr, see \ref Dune::wrap_or_move.
+ /// (2) Constructor. Forwards to (1) by wrapping reference into non-destroying shared_ptr.
template <class GB_,
REQUIRES(Concepts::Similar<GB_,GB>)>
- explicit VectorBase(GB_&& basis)
+ VectorBase(GB_&& basis)
: VectorBase(Dune::wrap_or_move(FWD(basis)))
{}
@@ -99,10 +83,6 @@ namespace AMDiS
std::shared_ptr<GlobalBasis> const& basis() const { return basis_; }
std::shared_ptr<GlobalBasis> const& basis() { return basis_; }
- /// Return the communication object \ref comm_
- std::shared_ptr<Comm> const& comm() const { return comm_; }
- std::shared_ptr<Comm> const& comm() { return comm_; }
-
/// Return the underlying linear algebra backend
Backend const& backend() const { return backend_; }
Backend& backend() { return backend_; }
@@ -423,9 +403,6 @@ namespace AMDiS
/// The finite element space / basis associated with the data vector
std::shared_ptr<GlobalBasis> basis_;
- /// A Communication object providing an indexSet, remoteIndices, and a dofMap.
- std::shared_ptr<Comm> comm_;
-
/// Data backend
Backend backend_;
diff --git a/src/amdis/linearalgebra/eigen/MatrixBackend.hpp b/src/amdis/linearalgebra/eigen/MatrixBackend.hpp
index 459d5203c2f408b2a0d87789f7e580b6f415d415..ed184958029c7b253225cf7ae28ea2dc18e5b494 100644
--- a/src/amdis/linearalgebra/eigen/MatrixBackend.hpp
+++ b/src/amdis/linearalgebra/eigen/MatrixBackend.hpp
@@ -31,7 +31,8 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseMatrix.
- MatrixBackend(std::shared_ptr<typename Traits::Comm> const&) {}
+ template <class Basis>
+ explicit MatrixBackend(Basis const&, Basis const&) {}
/// Return a reference to the data-matrix \ref matrix
BaseMatrix& matrix()
diff --git a/src/amdis/linearalgebra/eigen/VectorBackend.hpp b/src/amdis/linearalgebra/eigen/VectorBackend.hpp
index af9e0f6879079ec124b4ac00b58b968b5cd9ae16..5efb9506bab2b8a132d5eefd301270267032ecb5 100644
--- a/src/amdis/linearalgebra/eigen/VectorBackend.hpp
+++ b/src/amdis/linearalgebra/eigen/VectorBackend.hpp
@@ -33,7 +33,8 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseVector.
- VectorBackend(std::shared_ptr<typename Traits::Comm> const&) {}
+ template <class Basis>
+ explicit VectorBackend(Basis const&) {}
/// Return the data-vector \ref vector_
BaseVector const& vector() const
diff --git a/src/amdis/linearalgebra/istl/Communication.hpp b/src/amdis/linearalgebra/istl/Communication.hpp
index a4139075eb6e06870ad1d24c81945b22425d1174..08a1f2a79089ee582fe4c5d13d874d719ebbc2fd 100644
--- a/src/amdis/linearalgebra/istl/Communication.hpp
+++ b/src/amdis/linearalgebra/istl/Communication.hpp
@@ -4,6 +4,7 @@
#include <string>
#include <dune/common/std/optional.hh>
+
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/solvercategory.hh>
#include <dune/istl/paamg/pinfo.hh>
diff --git a/src/amdis/linearalgebra/istl/Communication.inc.hpp b/src/amdis/linearalgebra/istl/Communication.inc.hpp
index b58be7fd44e21a1763806df51b3d599e8a8c19eb..3886a9d35b8256e4845545e55fd2d82b888fff70 100644
--- a/src/amdis/linearalgebra/istl/Communication.inc.hpp
+++ b/src/amdis/linearalgebra/istl/Communication.inc.hpp
@@ -9,12 +9,18 @@ namespace AMDiS {
template <class Basis>
std::unique_ptr<ISTLCommunication<Basis>>
-CommunicationCreator<ISTLCommunication<Basis>>::
-create(Basis const& basis, std::string const& prefix)
+CommunicationCreator<ISTLCommunication<Basis>>
+ ::create(Basis const& basis, std::string const& prefix)
{
using SolverType = Dune::SolverCategory::Category;
- std::string category = "default";
+ std::string category = "";
Parameters::get(prefix + "->category", category);
+ if (category == "")
+ {
+ // check global parameter if no local one was found
+ category = "default";
+ Parameters::get("solver category", category);
+ }
SolverType cat;
auto const& gv = basis.gridView();
int mpiSize = gv.comm().size();
diff --git a/src/amdis/linearalgebra/istl/MatrixBackend.hpp b/src/amdis/linearalgebra/istl/MatrixBackend.hpp
index 4d0b4a65545871e7197bbed4ffb203658322e969..a62e1135208ab8752918ff58e4492188c92a6489 100644
--- a/src/amdis/linearalgebra/istl/MatrixBackend.hpp
+++ b/src/amdis/linearalgebra/istl/MatrixBackend.hpp
@@ -28,7 +28,8 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseVector.
- MatrixBackend(std::shared_ptr<typename Traits::Comm> const&) {}
+ template <class Basis>
+ explicit MatrixBackend(Basis const&, Basis const&) {}
/// Return the data-vector \ref vector
BaseMatrix const& matrix() const
diff --git a/src/amdis/linearalgebra/istl/VectorBackend.hpp b/src/amdis/linearalgebra/istl/VectorBackend.hpp
index 03bc10283c42a50fba5e0a3a54bceb100c650923..378248a630c4cada6d7c3ce077430f786fdf071c 100644
--- a/src/amdis/linearalgebra/istl/VectorBackend.hpp
+++ b/src/amdis/linearalgebra/istl/VectorBackend.hpp
@@ -28,7 +28,8 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseVector.
- VectorBackend(std::shared_ptr<typename Traits::Comm> const&) {}
+ template <class Basis>
+ explicit VectorBackend(Basis const&) {}
/// Return the data-vector \ref vector
BaseVector const& vector() const
diff --git a/src/amdis/linearalgebra/mtl/MatrixBackend.hpp b/src/amdis/linearalgebra/mtl/MatrixBackend.hpp
index c002b8e80f2a9a678e56c6e5e3ff6b11c12b9923..364bf7c97c21abd66a52ec010ba70305e2a67a52 100644
--- a/src/amdis/linearalgebra/mtl/MatrixBackend.hpp
+++ b/src/amdis/linearalgebra/mtl/MatrixBackend.hpp
@@ -35,7 +35,8 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseMatrix.
- MatrixBackend(std::shared_ptr<typename Traits::Comm> const&) {}
+ template <class Basis>
+ explicit MatrixBackend(Basis const&, Basis const&) {}
/// Return a reference to the data-matrix \ref matrix
BaseMatrix& matrix()
diff --git a/src/amdis/linearalgebra/mtl/VectorBackend.hpp b/src/amdis/linearalgebra/mtl/VectorBackend.hpp
index f9b3602aa38fea4978a9352ab3ff46b14d165ede..c01085222015f090517b32997fafae629de6d2f4 100644
--- a/src/amdis/linearalgebra/mtl/VectorBackend.hpp
+++ b/src/amdis/linearalgebra/mtl/VectorBackend.hpp
@@ -31,7 +31,8 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseVector.
- VectorBackend(std::shared_ptr<typename Traits::Comm> const&) {}
+ template <class Basis>
+ explicit VectorBackend(Basis const&) {}
/// Return the data-vector \ref vector_
BaseVector const& vector() const
diff --git a/src/amdis/linearalgebra/petsc/MatrixBackend.hpp b/src/amdis/linearalgebra/petsc/MatrixBackend.hpp
index 656fdbdc04bc54efa90eaca50c9487f1263999d5..2c70d3040f63ca17edfb15b7e408b8220b5bc7e1 100644
--- a/src/amdis/linearalgebra/petsc/MatrixBackend.hpp
+++ b/src/amdis/linearalgebra/petsc/MatrixBackend.hpp
@@ -34,8 +34,9 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseMatrix.
- MatrixBackend(std::shared_ptr<Communication> comm)
- : comm_(std::move(comm))
+ template <class Basis>
+ explicit MatrixBackend(Basis const& rowBasis, Basis const& /*colBasis*/)
+ : comm_(rowBasis.communicator())
{}
// disable copy and move operations
@@ -109,7 +110,7 @@ namespace AMDiS
Dune::Timer t;
// update the non-zeros structure of the matrix.
- pattern_.update(rowBasis, *comm_);
+ pattern_.update(rowBasis);
info(3, " update pattern needed {} seconds", t.elapsed());
t.reset();
@@ -173,7 +174,7 @@ namespace AMDiS
// The local-to-global map
auto const& dofMap() const
{
- return comm_->dofMap();
+ return comm_.dofMap();
}
// An MPI Communicator or PETSC_COMM_SELF
@@ -194,7 +195,7 @@ namespace AMDiS
}
private:
- std::shared_ptr<Communication> comm_;
+ Communication const& comm_;
/// The data-matrix
Mat matrix_;
diff --git a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.hpp b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.hpp
index ea504b285e19f1598d56f136185f19f4746bb829..29423da2538722e705112295b194f35058b09016 100644
--- a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.hpp
+++ b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.hpp
@@ -25,8 +25,8 @@ namespace AMDiS
/// Calculate the matrix pattern from the local pattern of a basis.
// NOTE: this assumes a square matrix with row-basis = col-basis.
- template <class Basis, class Communication>
- void update(Basis const& basis, Communication const& comm);
+ template <class Basis>
+ void update(Basis const& basis);
private:
std::vector<PetscInt> dnnz_; //< number of nonzeros in the diagonal part (owner part)
diff --git a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp
index 63793587b43834c2f4805b8441971f682001247f..29e493b6c32f2dbf8048d00bbc883797735637cd 100644
--- a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp
+++ b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp
@@ -6,7 +6,7 @@
#include <dune/common/timer.hh>
#include <dune/grid/common/partitionset.hh>
-#include <amdis/GridTransferManager.hpp>
+#include <amdis/AdaptiveGrid.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/parallel/Communicator.hpp>
#include <amdis/common/parallel/RequestOperations.hpp>
@@ -14,10 +14,11 @@
namespace AMDiS {
-template <class Basis, class Communication>
-void MatrixNnzStructure::update(Basis const& basis, Communication const& c)
+template <class Basis>
+void MatrixNnzStructure::update(Basis const& basis)
{
- unsigned long newChangeIndex = GridTransferManager::get(basis.gridView().grid())->changeIndex();
+ using Grid = AdaptiveGrid<typename Basis::GridView::Grid>;
+ unsigned long newChangeIndex = Grid::instance(basis.gridView().grid())->changeIndex();
bool gridChanged = newChangeIndex > changeIndex_;
changeIndex_ = newChangeIndex;
@@ -27,7 +28,7 @@ void MatrixNnzStructure::update(Basis const& basis, Communication const& c)
Dune::Timer t;
- auto const& dofMap = c.dofMap();
+ auto const& dofMap = basis.communicator().dofMap();
std::size_t localSize = dofMap.localSize();
test_exit(localSize > 0, "Local domain has no owned DOFs!");
@@ -84,7 +85,7 @@ void MatrixNnzStructure::update(Basis const& basis, Communication const& c)
std::vector<RemoteNz> sendList(world.size());
std::vector<std::size_t> receiveList(world.size(), 0);
- for (auto const& rim : c.remoteIndices()) {
+ for (auto const& rim : basis.communicator().remoteIndices()) {
int p = rim.first;
auto* sourceRemoteIndexList = rim.second.first;
@@ -187,4 +188,4 @@ void MatrixNnzStructure::update(Basis const& basis, Communication const& c)
msg("update MatrixNnzStructure need {} sec", t.elapsed());
}
-} // end namespace AMDiS
\ No newline at end of file
+} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/petsc/VectorBackend.hpp b/src/amdis/linearalgebra/petsc/VectorBackend.hpp
index 2148624169e6437c43f5a850ce967d163d8905aa..3634c607f3701c137d73e9adf8380575025cb7d3 100644
--- a/src/amdis/linearalgebra/petsc/VectorBackend.hpp
+++ b/src/amdis/linearalgebra/petsc/VectorBackend.hpp
@@ -3,6 +3,7 @@
#include <petscvec.h>
#include <dune/common/ftraits.hh>
+#include <dune/common/shared_ptr.hh>
#include <amdis/Output.hpp>
#include <amdis/common/FakeContainer.hpp>
@@ -36,8 +37,9 @@ namespace AMDiS
public:
/// Constructor. Constructs new BaseVector.
- VectorBackend(std::shared_ptr<Communication> comm)
- : comm_(std::move(comm))
+ template <class Basis>
+ explicit VectorBackend(Basis const& basis)
+ : comm_(Dune::stackobject_to_shared_ptr(basis.communicator()))
{}
/// Move constructor
@@ -305,7 +307,7 @@ namespace AMDiS
}
private:
- std::shared_ptr<Communication> comm_;
+ std::shared_ptr<Communication const> comm_;
/// The data-vector
mutable Vec vector_ = nullptr;
diff --git a/test/BackupRestoreTest.cpp b/test/BackupRestoreTest.cpp
index 4f6cc20f3201d2599f2742594324083f194cbb01..4749405fb7996085e556c9adf86ccf2fdfaffcdb 100644
--- a/test/BackupRestoreTest.cpp
+++ b/test/BackupRestoreTest.cpp
@@ -58,7 +58,7 @@ void test(Factory factory)
AMDIS_TEST_EQ(num_elements, prob.grid()->size(0));
AMDIS_TEST_EQ(num_vertices, prob.grid()->size(Grid::dimension));
- DOFVector<typename Param::GlobalBasis> vec(prob.globalBasis(), prob.communication());
+ DOFVector<typename Param::GlobalBasis> vec(prob.globalBasis());
vec.child(Dune::Indices::_0) << [](auto const& x) { return x; };
vec.child(Dune::Indices::_1) << [](auto const& x) { return two_norm(x); };
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index a1095005bfa09cb53384b413bc7492d46e045760..046593ecb00fe01dff83394043b62ff0c251510b 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -67,11 +67,13 @@ dune_add_test(SOURCES GradientTest.cpp
dune_add_test(SOURCES HybridSizeTest.cpp
LINK_LIBRARIES amdis)
-dune_add_test(SOURCES ISTLCommTest.cpp
- LINK_LIBRARIES amdis
- MPI_RANKS 2 4
- TIMEOUT 300
- CMAKE_GUARD "MPI_FOUND AND dune-istl_FOUND")
+if(BACKEND STREQUAL "ISTL")
+ dune_add_test(SOURCES ISTLCommTest.cpp
+ LINK_LIBRARIES amdis
+ MPI_RANKS 2 4
+ TIMEOUT 300
+ CMAKE_GUARD MPI_FOUND)
+endif()
dune_add_test(SOURCES IntegrateTest.cpp
LINK_LIBRARIES amdis)
@@ -94,6 +96,9 @@ dune_add_test(SOURCES MultiTypeMatrixTest.cpp
dune_add_test(SOURCES NodeIndicesTest.cpp
LINK_LIBRARIES amdis)
+dune_add_test(SOURCES ObserverTest.cpp
+ LINK_LIBRARIES amdis)
+
dune_add_test(SOURCES OperationsTest.cpp
LINK_LIBRARIES amdis)
diff --git a/test/DOFVectorTest.cpp b/test/DOFVectorTest.cpp
index b4258119dd093c19afb7945ccfd10ec0e9fbecec..7a448f56947f2da4182c3c53c4d205d0df3cbb40 100644
--- a/test/DOFVectorTest.cpp
+++ b/test/DOFVectorTest.cpp
@@ -1,15 +1,18 @@
#include "config.h"
+#include <memory>
+
#include <dune/common/filledarray.hh>
+#include <dune/common/shared_ptr.hh>
#include <dune/grid/yaspgrid.hh>
-#include <dune/functions/backends/concepts.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <amdis/AdaptiveGrid.hpp>
#include <amdis/AMDiS.hpp>
-#include <amdis/GridTransferManager.hpp>
#include <amdis/LinearAlgebra.hpp>
+#include <amdis/functions/ParallelGlobalBasis.hpp>
#include "Tests.hpp"
@@ -17,23 +20,35 @@ using namespace AMDiS;
using namespace Dune::Functions::BasisFactory;
template <class B, class T>
-void test_dofvector(B const& basis, DOFVector<B,T>& vec)
+void test_dofvector(DOFVector<B,T>& vec)
{
- // TODO:
- // AMDIS_TEST( vec.size() == basis.dimension() );
+ using size_type = typename DOFVector<B,T>::GlobalBasis::LocalView::size_type;
+ using value_type = typename DOFVector<B,T>::value_type;
+
+ auto const& basis = *vec.basis();
+ AMDIS_TEST( vec.localSize() == basis.dimension() );
vec.resizeZero();
-// vec = T(0);
-// vec[0] = T(1);
-// #if HAVE_MTL || HAVE_ISTL
-// auto m0 = std::min_element(std::begin(vec.vector()), std::end(vec.vector()));
-// auto m1 = std::max_element(std::begin(vec.vector()), std::end(vec.vector()));
-// AMDIS_TEST( *m0 == T(0) );
-// AMDIS_TEST( *m1 == T(1) );
-// #endif
-
-// // DOFVector models the concept of a VectorBackend in Dune::Functions
-// static_assert(Dune::models<Dune::Functions::Concept::VectorBackend<B>, decltype(vec)>(), "");
+ vec.init(false);
+ vec.insert(0, value_type(1));
+ vec.finish();
+
+ value_type min(10);
+ value_type max(-10);
+ auto lv = basis.localView();
+ auto const& gv = basis.gridView();
+ for (auto const& e : elements(gv))
+ {
+ lv.bind(e);
+ for (size_type i = 0; i < lv.size(); ++i)
+ {
+ auto value = vec.at(lv.index(i));
+ min = std::min(min, value);
+ max = std::max(max, value);
+ }
+ }
+ AMDIS_TEST( min == T(0) || vec.localSize() == 1 );
+ AMDIS_TEST( max == T(1) );
}
int main(int argc, char** argv)
@@ -41,52 +56,56 @@ int main(int argc, char** argv)
Environment env(argc, argv);
// create grid
+ using HostGrid = Dune::YaspGrid<2>;
Dune::FieldVector<double, 2> L; L = 1.0;
auto s = Dune::filledArray<2>(1);
- Dune::YaspGrid<2> grid(L, s);
- auto const& gridView = grid.leafGridView();
+ HostGrid hostGrid(L, s);
+ auto grid = AdaptiveGrid<HostGrid>::instance(hostGrid);
+ auto const& gridView = grid->leafGridView();
// create basis
- auto basis = makeBasis(gridView,
- composite(power<2>(lagrange<2>(), flatInterleaved()), lagrange<1>(), flatLexicographic()));
- using Basis = decltype(basis);
+ auto preBasis = composite(power<2>(lagrange<2>(), flatInterleaved()),
+ lagrange<1>(),
+ flatLexicographic());
+ auto basis = makeGlobalBasis(gridView, preBasis);
- using Traits = BackendTraits<Basis>;
- auto comm = CommunicationCreator<typename Traits::Comm>::create(basis);
-
- GridTransferManager::attach(grid, *comm, [&]()
- {
- basis.update(gridView);
- comm->update(basis);
- });
+ using Basis = decltype(basis);
{
- DOFVector<Basis> vec1(basis, *comm);
- test_dofvector(basis, vec1);
+ DOFVector<Basis> vec1(basis);
+ test_dofvector(vec1);
- DOFVector<Basis, double> vec2(basis, *comm);
- test_dofvector(basis, vec2);
+ DOFVector<Basis, double> vec2(basis);
+ test_dofvector(vec2);
- auto vec3 = makeDOFVector(basis, *comm);
- test_dofvector(basis, vec3);
+ auto vec3 = makeDOFVector(basis);
+ test_dofvector(vec3);
- auto vec4 = makeDOFVector<double>(basis, *comm);
- test_dofvector(basis, vec4);
+ auto vec4 = makeDOFVector<double>(basis);
+ test_dofvector(vec4);
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
- DOFVector vec5(basis, *comm);
- test_dofvector(basis, vec5);
+ DOFVector vec5(basis);
+ test_dofvector(vec5);
#endif
}
- // test GridTransferManager registration
+ // test automatic updates
{
- DOFVector<Basis> vec1(basis, *comm);
- test_dofvector(basis, vec1);
+ DOFVector<Basis> vec1(basis);
+
+ // Conversion from Dune::Functions::DefaultGlobalBasis
+ auto vec2 = makeDOFVector(makeBasis(gridView, preBasis));
+
+ test_dofvector(vec1);
for (auto const& e : elements(gridView))
- grid.mark(1, e);
- GridTransferManager::adapt(grid);
- // AMDIS_TEST_EQ(vec1.size(), basis.dimension());
+ grid->mark(1, e);
+ grid->preAdapt();
+ grid->adapt();
+ grid->postAdapt();
+ std::size_t newSize = basis.dimension();
+ AMDIS_TEST_EQ(vec1.localSize(), newSize);
+ AMDIS_TEST_EQ(vec2.localSize(), newSize);
}
report_errors();
diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp
index bf164eafa7cb6809bdd0faf7ec15cd4db016f503..f3ad5286e5980c38ddeb760ff6d42401d81346ce 100644
--- a/test/DiscreteFunctionTest.cpp
+++ b/test/DiscreteFunctionTest.cpp
@@ -14,21 +14,30 @@
using namespace AMDiS;
-using ElliptParam = YaspGridBasis<2, 2>;
+using ElliptParam = YaspGridBasis<2, 2, 2>;
using ElliptProblem = ProblemStat<ElliptParam>;
template <class GB, class T>
bool comp(DOFVector<GB,T> const& U, DOFVector<GB,T> const& V)
{
- // TODO:
- // if (U.size() != V.size())
- // return false;
-
- // using Int = typename DOFVector<GB,T>::size_type;
- // for (Int i = 0; i < U.size(); ++i) {
- // if (U[i] != V[i])
- // return false;
- // }
+ if (U.localSize() != V.localSize() || U.globalSize() != V.globalSize())
+ return false;
+ using size_type = typename DOFVector<GB,T>::GlobalBasis::LocalView::size_type;
+ auto const& basis = U.basis();
+ if (basis.get() != V.basis().get())
+ return false;
+ auto lv = basis->localView();
+ auto const& gv = basis->gridView();
+ for (auto const& e : elements(gv))
+ {
+ lv.bind(e);
+ for (size_type i = 0; i < lv.size(); ++i)
+ {
+ auto multiIndex = lv.index(i);
+ if (std::abs(U.at(multiIndex) - V.at(multiIndex)) > AMDIS_TEST_TOL)
+ return false;
+ }
+ }
return true;
}
@@ -50,15 +59,24 @@ int main(int argc, char** argv)
auto u1 = makeDOFVectorView(U1);
auto u2 = makeDOFVectorView(U2);
- auto expr = [](auto const& x) { return 1 + x[0] + x[1]; };
- u0.interpolate_noalias(expr);
- u1.interpolate(expr);
- u2 << expr;
+ using Range = typename decltype(u0)::Range;
+
+ auto expr1 = [](auto const& x) {
+ return Range{ 1 + x[0] + x[1],
+ 2 + x[0] + x[1]};
+ };
+
+ u0.interpolate_noalias(expr1);
+ u1.interpolate(expr1);
+ u2 << expr1;
AMDIS_TEST( comp(U0, U1) );
AMDIS_TEST( comp(U0, U2) );
- auto expr2 = [](auto const& x) { return 1 + 2*x[0] - 3*x[1]; };
+ auto expr2 = [](auto const& x) {
+ return Range{ 1 + 2*x[0] - 3*x[1],
+ 2 + 2*x[0] - 3*x[1]};
+ };
u0.interpolate_noalias(u2 + expr2);
u1.interpolate(u1 + expr2);
@@ -67,7 +85,10 @@ int main(int argc, char** argv)
AMDIS_TEST( comp(U0, U1) );
AMDIS_TEST( comp(U0, U2) );
- auto expr3 = [](auto const& x) { return 1 - 0.5*x[0] - 2*x[1]; };
+ auto expr3 = [](auto const& x) {
+ return Range{ 1 - 0.5*x[0] - 2*x[1],
+ 2 - 0.5*x[0] - 2*x[1]};
+ };
u0.interpolate_noalias(u2 - expr3);
u1.interpolate(u1 - expr3);
@@ -87,20 +108,23 @@ int main(int argc, char** argv)
auto y0 = localFct(local);
auto y1 = u0(geo.global(local));
- AMDIS_TEST(std::abs(y0 - y1) < 1.e-10);
+ for (auto it1 = y0.begin(), it2 = y1.begin(); it1 != y0.end(); ++it1, ++it2)
+ AMDIS_TEST(std::abs(*it1 - *it2) < 1.e-10);
+
localFct.unbind();
dlocalFct.bind(e);
auto g0 = dlocalFct(local);
auto g1 = du0(geo.global(local));
- AMDIS_TEST(two_norm(g0 - g1) < 1.e-10);
+ for (auto it1 = g0.begin(), it2 = g1.begin(); it1 != g0.end(); ++it1, ++it2)
+ AMDIS_TEST(two_norm(*it1 - *it2) < 1.e-10);
dlocalFct.unbind();
}
- auto V0 = makeDOFVector<double>(*prob.globalBasis());
+ auto V0 = makeDOFVector<double>(prob.globalBasis());
auto v0 = makeDOFVectorView(V0);
- v0 << expr;
+ v0 << expr1;
// test makeDiscreteFunction
int preTP1 = 0;
@@ -114,9 +138,12 @@ int main(int argc, char** argv)
auto w2 = makeDiscreteFunction(W2, tp);
// test makeDOFVectorView with (pre)treepath argument
+ auto expr = [](auto const& x) { return 1 + x[0] + x[1]; };
auto W3 = *prob.solutionVector();
auto W4 = *prob.solutionVector();
auto W5 = *prob.solutionVector();
+ auto W7 = *prob.solutionVector();
+ auto W8 = *prob.solutionVector();
auto w3 = makeDOFVectorView(W3, preTP1);
auto w4 = makeDOFVectorView(W4, preTP2);
auto w5 = makeDOFVectorView(W5, tp);
@@ -130,5 +157,14 @@ int main(int argc, char** argv)
AMDIS_TEST( comp(W3, W5) );
AMDIS_TEST( comp(W3, W6) );
+ // test interpolation on subbasis
+ auto w7 = makeDOFVectorView(W7);
+ auto w8_0 = makeDOFVectorView(W8, 0);
+ auto w8_1 = makeDOFVectorView(W8, 1);
+ w7 << expr;
+ w8_0 << expr;
+ w8_1 << expr;
+ AMDIS_TEST( comp(W7, W8) );
+
return 0;
}
diff --git a/test/GlobalIdSetTest.cpp b/test/GlobalIdSetTest.cpp
index 6cd937d1661648097de28882e540d499ea44816b..b92485e7eeac845e94e3891392c3b92f423be17b 100644
--- a/test/GlobalIdSetTest.cpp
+++ b/test/GlobalIdSetTest.cpp
@@ -42,7 +42,7 @@ template <class Grid>
void checkGrid(const Grid& grid)
{
// create basis
- using namespace Dune::Functions::BasisBuilder;
+ using namespace Dune::Functions::BasisFactory;
auto basis1 = makeBasis(grid.leafGridView(),
composite(
power<2>(lagrange<2>()),
diff --git a/test/GradientTest.cpp b/test/GradientTest.cpp
index 17b6228237e4086f7988579cbcdd0c739ba3e9c0..7d6d1ef3f2deff5a8c3a6428e1971e733ec1dcba 100644
--- a/test/GradientTest.cpp
+++ b/test/GradientTest.cpp
@@ -4,6 +4,8 @@
#include <iostream>
#include <amdis/AMDiS.hpp>
+
+#include <amdis/AdaptiveGrid.hpp>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/GridFunctions.hpp>
#include <amdis/Integrate.hpp>
@@ -15,11 +17,13 @@ using namespace AMDiS;
template <std::size_t p>
void test()
{
- using Grid = Dune::YaspGrid<2>;
- using Basis1 = LagrangeBasis<Grid,p,p>;
+ using HostGrid = Dune::YaspGrid<2>;
+ using Basis1 = LagrangeBasis<HostGrid,p,p>;
+
+ HostGrid hostGrid({1.0, 1.0}, {16,16});
+ auto grid = AdaptiveGrid<HostGrid>::instance(hostGrid);
+ auto const& gridView = grid->leafGridView();
- Grid grid({1.0, 1.0}, {16,16});
- auto gridView = grid.leafGridView();
auto basis = Basis1::create(gridView);
auto uVector = makeDOFVector(basis);
@@ -57,7 +61,8 @@ void test()
u_0 << evalAtQP([](auto const& x) { return std::sin(x[0]) + std::cos(x[1]); });
u_1 << evalAtQP([](auto const& x) { return -std::sin(x[0]) * std::cos(x[1]); });
- AMDIS_TEST_APPROX((integrate(divergenceAtQP(u), gridView)), (integrate(partialAtQP(u_0,0) + partialAtQP(u_1,1), gridView)));
+ AMDIS_TEST_APPROX((integrate(divergenceAtQP(u), gridView)),
+ (integrate(partialAtQP(u_0,0) + partialAtQP(u_1,1), gridView)));
auto vVector(uVector);
auto v = makeDOFVectorView(vVector);
@@ -65,7 +70,9 @@ void test()
auto v_1 = makeDOFVectorView(vVector, 1);
// test gradient
- v << evalAtQP([](auto const& x) { return Dune::FieldVector<double,2>{std::cos(x[0]), -std::sin(x[1])}; });
+ v << evalAtQP([](auto const& x) {
+ return Dune::FieldVector<double,2>{std::cos(x[0]), -std::sin(x[1])};
+ });
AMDIS_TEST((std::sqrt(integrate(unary_dot(v - gradientAtQP(u_0)), gridView))) < 0.05);
// test divergence
diff --git a/test/ISTLCommTest.cpp b/test/ISTLCommTest.cpp
index 60a25034b8be87a4648b6b299a05bcf96e8bae2e..1ef3d0e1644daa6eb956ac47381c85a7bae1e5ec 100644
--- a/test/ISTLCommTest.cpp
+++ b/test/ISTLCommTest.cpp
@@ -13,6 +13,7 @@
#include <dune/grid/yaspgrid.hh>
#include <amdis/linearalgebra/istl/Communication.hpp>
+#include <amdis/AdaptiveGrid.hpp>
#include <amdis/Environment.hpp>
#include <amdis/ProblemStatTraits.hpp>
@@ -23,7 +24,6 @@ using namespace AMDiS;
template <class Basis>
bool test(Basis& basis, std::string const& prefix)
{
- using Comm = ISTLCommunication<Basis>;
using GV = typename Basis::GridView;
using Grid = typename GV::Grid;
using Coord = typename Grid::template Codim<0>::Geometry::GlobalCoordinate;
@@ -32,17 +32,15 @@ bool test(Basis& basis, std::string const& prefix)
bool passed = true;
+#if HAVE_MPI
// Interpolate test function f onto dofs and reference
auto f = [&](const Coord& x) -> T { return std::pow(x[0],3) + x[1] + (d == 3 ? x[2] : 0.0); };
std::vector<T> ref(basis.size());
Dune::Functions::interpolate(basis, ref, f);
std::vector<T> dofs = ref;
- // Make communicator and exchange dofs
-#if HAVE_MPI
- using CommCreator = CommunicationCreator<Comm>;
- auto comm = CommCreator::create(basis, prefix);
- comm->get().copyOwnerToAll(dofs, dofs);
+ // Exchange dofs
+ basis.communicator().get().copyOwnerToAll(dofs, dofs);
// Compare post-exchange data with reference
for (std::size_t i = 0; i < dofs.size(); ++i)
@@ -77,16 +75,24 @@ int main(int argc, char** argv)
int ovlp = 1;
{
- auto grid = std::make_unique<YaspGrid2d>(lower2d, upper2d, nElements2d, std::bitset<2>(), ovlp);
+ auto hostGrid = std::make_unique<YaspGrid2d>(lower2d, upper2d, nElements2d, std::bitset<2>(),
+ ovlp);
+ auto grid = AdaptiveGrid<YaspGrid2d>::instance(hostGrid);
+
auto l1 = LagrangeBasis<YaspGrid2d, 1>::create(grid->leafGridView());
AMDIS_TEST(test(l1, "Yasp2d_L1_Ovlp"));
+
auto th = TaylorHoodBasis<YaspGrid2d, 1>::create(grid->leafGridView());
AMDIS_TEST(test(th, "Yasp2d_TH_Ovlp"));
}
{
- auto grid = std::make_unique<YaspGrid3d>(lower3d, upper3d, nElements3d, std::bitset<3>(), ovlp);
+ auto hostGrid = std::make_unique<YaspGrid3d>(lower3d, upper3d, nElements3d, std::bitset<3>(),
+ ovlp);
+ auto grid = AdaptiveGrid<YaspGrid3d>::instance(hostGrid);
+
auto l1 = LagrangeBasis<YaspGrid3d, 1>::create(grid->leafGridView());
AMDIS_TEST(test(l1, "Yasp3d_L1_Ovlp"));
+
auto th = TaylorHoodBasis<YaspGrid3d, 1>::create(grid->leafGridView());
AMDIS_TEST(test(th, "Yasp3d_TH_Ovlp"));
}
diff --git a/test/MarkerTest.cpp b/test/MarkerTest.cpp
index 6ae27b0def0ca4b167c257ce9bcf5428ebca2811..bc931654c48512dd8e9260c805342a4b15def61c 100644
--- a/test/MarkerTest.cpp
+++ b/test/MarkerTest.cpp
@@ -39,9 +39,9 @@ int main(int argc, char** argv)
};
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
- GridFunctionMarker marker("mymarker", grid, markerFunc);
+ GridFunctionMarker marker("mymarker", prob.grid(), markerFunc);
#else
- auto marker = makeGridFunctionMarker("mymarker", grid, markerFunc);
+ auto marker = makeGridFunctionMarker("mymarker", prob.grid(), markerFunc);
#endif
prob.addMarker(marker);
diff --git a/test/ObserverTest.cpp b/test/ObserverTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..77039feeb3e94f4dce53ace0d9e41dc9dcf404b6
--- /dev/null
+++ b/test/ObserverTest.cpp
@@ -0,0 +1,180 @@
+#include <iostream>
+#include <memory>
+#include <string>
+
+#include <amdis/Observer.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+namespace event {
+ struct fooEvent {};
+ struct barEvent {};
+ struct bazEvent {};
+}
+
+std::string testValue = "";
+
+bool checkAndResetTestValue(std::string testName, std::string ref, std::string refOpt = "X")
+{
+ bool result = testValue == ref || testValue == refOpt;
+ testValue.clear();
+ return result;
+}
+
+class A
+ : public Signals<event::fooEvent, event::barEvent, event::bazEvent>
+{
+public:
+ void foo()
+ {
+ std::cout << "A::foo()\n";
+ testValue += "A.foo ";
+ event::fooEvent e;
+ this->notify(e);
+ }
+
+ void bar()
+ {
+ std::cout << "A::bar()\n";
+ testValue += "A.bar ";
+ event::barEvent e;
+ this->notify(e);
+ }
+
+ void baz()
+ {
+ std::cout << "A::baz()\n";
+ testValue += "A.baz ";
+ event::bazEvent e;
+ this->notify(e);
+ }
+};
+
+class B
+ : public Observer<A, event::fooEvent, event::bazEvent>
+ , public Signals<event::bazEvent>
+{
+public:
+ B(std::shared_ptr<A> a)
+ : Observer<A, event::fooEvent, event::bazEvent>(a)
+ {}
+
+ void update(event::fooEvent const& e) override
+ {
+ std::cout << "B::update(foo)\n";
+ testValue += "B.foo ";
+ }
+
+ void update(event::bazEvent const& e) override
+ {
+ std::cout << "B::update(baz)\n";
+ testValue += "B.baz ";
+ this->notify(e);
+ }
+};
+
+class C
+ : public Observer<B, event::barEvent, event::bazEvent>
+{
+public:
+ C(std::shared_ptr<B> b)
+ : Observer<B, event::barEvent, event::bazEvent>(b)
+ {}
+
+ void update(event::barEvent const& e) override
+ {
+ std::cout << "C::update(bar)\n";
+ testValue += "C.bar ";
+ }
+
+ void update(event::bazEvent const& e) override
+ {
+ std::cout << "C::update(baz)\n";
+ testValue += "C.baz ";
+ }
+};
+
+class D
+ : public Signals<event::barEvent>
+{
+public:
+ void bar()
+ {
+ std::cout << "D::bar()\n";
+ testValue += "D.bar ";
+ event::barEvent e;
+ this->notify(e);
+ }
+};
+
+class E
+ : public Observer<A, event::fooEvent>
+ , public Observer<D, event::barEvent>
+{
+public:
+ E(std::shared_ptr<A> a, std::shared_ptr<D> d)
+ : Observer<A, event::fooEvent>(a)
+ , Observer<D, event::barEvent>(d)
+ {}
+
+ void update(event::fooEvent const& e) override
+ {
+ std::cout << "E::update(foo)\n";
+ testValue += "E.foo ";
+ }
+
+ void update(event::barEvent const& e) override
+ {
+ std::cout << "E::update(bar)\n";
+ testValue += "E.bar ";
+ }
+};
+
+/// Observing a class that inherits from Observer<T, Events...> more than once is not supported
+// class F
+// : public Observer<E, event::fooEvent>
+// {
+// public:
+// F(std::shared_ptr<E> e)
+// : Observer<E, event::fooEvent>(e) // Error
+// {}
+// };
+
+/// Test the observer hierarchy.
+/**
+ * There are 4 hierarchies tested.
+ * A::foo(): A -> B, E (multiple observers [B,E] of one event)
+ * A::bar(): A -> (B) -> C (indirect observer [C of A])
+ * A::baz(): A -> B -> C (simultanious observer and signaller of an event [B])
+ * D::bar(): D -> E (observer of multiple classes [E])
+ */
+int main()
+{
+ auto a = std::make_shared<A>();
+ auto b = std::make_shared<B>(a);
+ auto c = std::make_shared<C>(b);
+ auto d = std::make_shared<D>();
+ auto e = std::make_shared<E>(a, d);
+ std::string ref, refOpt;
+
+ a->foo();
+ ref = "A.foo B.foo E.foo ";
+ refOpt = "A.foo E.foo B.foo "; // Order of b->foo() and e->foo() unspecified
+ AMDIS_TEST(checkAndResetTestValue("A:foo()", ref, refOpt));
+
+ a->bar();
+ ref = "A.bar C.bar ";
+ AMDIS_TEST(checkAndResetTestValue("A:bar()", ref));
+
+ a->baz();
+ ref = "A.baz B.baz C.baz ";
+ AMDIS_TEST(checkAndResetTestValue("A:baz()", ref));
+
+ d->bar();
+ ref = "D.bar E.bar ";
+ AMDIS_TEST(checkAndResetTestValue("D:bar()", ref));
+
+ return report_errors();
+}
diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp
index c74a9cd69a7759992490c17125f49bb52b735bd2..9d44993f3fc5d15f6383e224c6ea8d6a3cc1934f 100644
--- a/test/ProblemStatTest.cpp
+++ b/test/ProblemStatTest.cpp
@@ -1,8 +1,10 @@
#include <amdis/AMDiS.hpp>
+#include <amdis/AdaptiveGrid.hpp>
#include <amdis/AdaptStationary.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/common/QuadMath.hpp>
+#include <amdis/functions/ParallelGlobalBasis.hpp>
#include <dune/common/version.hh>
#include <dune/grid/yaspgrid.hh>
@@ -13,15 +15,16 @@ template <class T>
void test()
{
// use T as coordinate type
- using Grid = Dune::YaspGrid<2, Dune::EquidistantCoordinates<T,2>>;
- Grid grid({T(1), T(1)}, {8,8});
+ using HostGrid = Dune::YaspGrid<2, Dune::EquidistantCoordinates<T,2>>;
+ HostGrid hostGrid({T(1), T(1)}, {8,8});
+ auto grid = AdaptiveGrid<HostGrid>::instance(hostGrid);
// use T as range type for basis
using namespace Dune::Functions::BasisFactory;
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
- auto basis = makeBasis(grid.leafGridView(), power<1>(lagrange<1>(), flatLexicographic()));
+ auto basis = makeGlobalBasis(grid->leafGridView(), power<1>(lagrange<1>(), flatLexicographic()));
#else
- auto basis = makeBasis(grid.leafGridView(), power<1>(lagrange<1,T>(), flatLexicographic()));
+ auto basis = makeGlobalBasis(grid->leafGridView(), power<1>(lagrange<1,T>(), flatLexicographic()));
#endif
using Basis = decltype(basis);