diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index cee60873f8f84859292fdaa345ece6dfb496cecb..d7470373d98664d0443f98f159fa996c28e3be1e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -26,20 +26,35 @@ dune-2.6 ubuntu-18.04 clang-6-17: - dunecontrol --current make -j4 examples -dune-git debian-10 gcc-8-17: - image: registry.dune-project.org/docker/ci/dune:git-debian-10-gcc-8-17 +.dune-git: + before_script: + - . /duneci/bin/duneci-init-job + - duneci-install-module https://gitlab.dune-project.org/core/dune-common.git + - duneci-install-module https://gitlab.dune-project.org/core/dune-geometry.git + - duneci-install-module https://gitlab.dune-project.org/core/dune-localfunctions.git + - duneci-install-module https://gitlab.dune-project.org/staging/dune-uggrid.git + - duneci-install-module https://gitlab.dune-project.org/core/dune-grid.git + - duneci-install-module https://gitlab.dune-project.org/core/dune-istl.git + - duneci-install-module https://gitlab.dune-project.org/staging/dune-typetree.git + - duneci-install-module https://gitlab.dune-project.org/staging/dune-functions.git script: - duneci-standard-test - dunecontrol --current make -j4 examples +dune-git debian-10 gcc-8-17: + extends: .dune-git + image: registry.dune-project.org/docker/ci/debian:10 + variables: + DUNECI_TOOLCHAIN: gcc-8-17 + dune-git debian-9 gcc-6-14: - image: registry.dune-project.org/docker/ci/dune:git-debian-9-gcc-6-14 - script: - - duneci-standard-test - - dunecontrol --current make -j4 examples + extends: .dune-git + image: registry.dune-project.org/docker/ci/debian:9 + variables: + DUNECI_TOOLCHAIN: gcc-6-14 dune-git ubuntu-18.04 clang-6-17: - image: registry.dune-project.org/docker/ci/dune:git-ubuntu-18.04-clang-6-17 - script: - - duneci-standard-test - - dunecontrol --current make -j4 examples + extends: .dune-git + image: registry.dune-project.org/docker/ci/ubuntu:18.04 + variables: + DUNECI_TOOLCHAIN: clang-6-17 diff --git a/src/amdis/AdaptiveGrid.hpp b/src/amdis/AdaptiveGrid.hpp index e9303e8c10c75b446273b23620e3862b28fa392c..a40fcf91dc2b5552109c420a4fc03da1a61728a0 100644 --- a/src/amdis/AdaptiveGrid.hpp +++ b/src/amdis/AdaptiveGrid.hpp @@ -8,103 +8,132 @@ #include <type_traits> #include <utility> -#include <dune/common/parallel/mpihelper.hh> +#include <dune/common/hybridutilities.hh> #include <dune/common/timer.hh> +#include <dune/common/version.hh> +#include <dune/common/parallel/mpihelper.hh> #include <dune/geometry/type.hh> +#include <dune/grid/common/backuprestore.hh> +#include <dune/grid/common/capabilities.hh> +#include <dune/grid/common/grid.hh> +#include <dune/grid/utility/structuredgridfactory.hh> -#include <amdis/common/ConceptsBase.hpp> +#include <amdis/common/Concepts.hpp> +#include <amdis/common/DefaultGridView.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. + + // forward declaration + template <class HostGrid> + class AdaptiveGridFamily; + + + /// \brief Wrapper class for Dune-grids that allows automatic signalling of events + /// during grid adaptation. + /** * 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> + * + * \tparam HG Host grid to be wrapped. Must implement the dune grid interface. + **/ + template <class HG> class AdaptiveGrid - : public Signals<event::preAdapt, event::adapt, event::postAdapt> + : public Dune::GridDefaultImplementation< + HG::dimension, HG::dimensionworld, typename HG::ctype, AdaptiveGridFamily<HG> > + , public Signals<event::preAdapt, event::adapt, event::postAdapt> { - using Self = AdaptiveGrid<Grid>; - using Element = typename Grid::template Codim<0>::Entity; - struct HiddenStruct {}; + using Self = AdaptiveGrid<HG>; 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; + using HostGrid = HG; + using GridFamily = AdaptiveGridFamily<HG>; + using Traits = typename GridFamily::Traits; + using Element = typename Traits::template Codim<0>::Entity; + - /// Constructor that may only be called indirectly via the instance method - explicit AdaptiveGrid(std::shared_ptr<HostGrid> grid, HiddenStruct) - : hostGrid_(std::move(grid)) + public: + + template <class HostGrid_, + Dune::disableCopyMove<Self, HostGrid_> = 0> + explicit AdaptiveGrid(HostGrid_&& hostGrid) + : hostGrid_(wrap_or_share(FWD(hostGrid))) {} - /// 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."); - } + /// Return the underlying grid + std::shared_ptr<HostGrid> const& hostGrid() const { return hostGrid_; } + + public: + + /// Grid iterators + /// @{ + + /// Iterator to first entity of given codim on level + template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition> + auto lbegin(int level) const { return hostGrid_->levelGridView(level).template begin<codim,pt>(); } + + /// one past the end on this level + template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition> + auto lend(int level) const { return hostGrid_->levelGridView(level).template end<codim,pt>(); } + + + /// Iterator to first leaf entity of given codim + template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition> + auto leafbegin() const { return hostGrid_->leafGridView().template begin<codim,pt>(); } + + /// One past the end of the sequence of leaf entities + template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition> + auto leafend() const { return hostGrid_->leafGridView().template end<codim,pt>(); } + + + /// Obtain begin intersection iterator with respect to the level GridView + auto ilevelbegin(Element const& e) const { return hostGrid_->levelGridView(e.level()).ibegin(e); } + + /// Obtain end intersection iterator with respect to the level GridView + auto ilevelend(Element const& e) const { return hostGrid_->levelGridView(e.level()).iend(e); } + + /// Obtain begin intersection iterator with respect to the leaf GridView + auto ileafbegin(Element const& e) const { return hostGrid_->leafGridView().ibegin(e); } + + /// Obtain end intersection iterator with respect to the leaf GridView + auto ileafend(Element const& e) const { return hostGrid_->leafGridView().iend(e); } + + /// @} + + + /// Size methods + /// @{ + + /// Return maximum level defined in this grid. int maxLevel() const { return hostGrid_->maxLevel(); } + /// Number of grid entities per level and codim int size(int level, int codim) const { return hostGrid_->size(level, codim); } /// Return number of leaf entities of a given codim in this process @@ -116,24 +145,32 @@ namespace AMDiS /// 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(); } + /// Returns the number of boundary segments within the macro grid + std::size_t 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(); } + /// Access to index and id sets + /// @{ - /// return const reference to the host grids local id set - LocalIdSet const& localIdSet() const { return hostGrid_->localIdSet(); } + /// Return const reference to the grids global id set + auto const& globalIdSet() const { return hostGrid_->globalIdSet(); } - /// 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 local id set + auto const& localIdSet() const { return hostGrid_->localIdSet(); } - /// return const reference to the host grids leaf index set - LeafIndexSet const& leafIndexSet() const { return hostGrid_->leafIndexSet(); } + /// Return const reference to the host grids level index set for level level + auto const& levelIndexSet(int level) const { return hostGrid_->levelIndexSet(level); } + + /// Return const reference to the host grids leaf index set + auto const& leafIndexSet() const { return hostGrid_->leafIndexSet(); } + + /// @} + + + /// Adaptivity and grid refinement + /// @{ /// Refines all grid elements refCount times. void globalRefine(int refCount) @@ -148,7 +185,7 @@ namespace AMDiS } } - /// Mark entity for refinement + /// Marks an entity to be refined/coarsened in a subsequent adapt bool mark(int refCount, Element const& e) { return hostGrid_->mark(refCount, e); } /// Return refinement mark for entity @@ -161,7 +198,7 @@ namespace AMDiS mightCoarsen_ = hostGrid_->preAdapt(); Dune::MPIHelper::getCollectiveCommunication().max(&mightCoarsen_, 1); this->notify(event::preAdapt{mightCoarsen_}); - info(2,"preAdapt needed {} seconds", t.elapsed()); + info(2,"AdaptiveGrid::preAdapt needed {} seconds", t.elapsed()); return mightCoarsen_; } @@ -172,6 +209,7 @@ namespace AMDiS refined_ = hostGrid_->adapt(); Dune::MPIHelper::getCollectiveCommunication().max(&refined_, 1); this->notify(event::adapt{mightCoarsen_ || refined_}); + info(2,"AdaptiveGrid::adapt needed {} seconds", t.elapsed()); return refined_; } @@ -182,126 +220,340 @@ namespace AMDiS hostGrid_->postAdapt(); this->notify(event::postAdapt{}); changeIndex_++; - info(2,"postAdapt needed {} seconds", t.elapsed()); + info(2,"AdaptiveGrid::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() + /// Returns the grid change index, see \ref changeIndex. + unsigned long changeIndex() const { - return Impl::loadBalanceImpl(*hostGrid_); + return changeIndex_; } - /* 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. - */ + // @} + + + /// Parallel data distribution and communication + /// @{ + + /// Return const reference to a collective communication object. + auto const& comm() const { return hostGrid_->comm(); } + + /// Communicate data of level gridView template <class DataHandle> - bool loadBalance (DataHandle& handle) + void communicate(DataHandle& handle, Dune::InterfaceType iftype, + Dune::CommunicationDirection dir, int level) const { - return Impl::loadBalanceImpl(*hostGrid_, handle); + hostGrid_->levelGridView(level).communicate(handle,iftype,dir); } - /// Returns the grid change index, see \ref changeIndex. - unsigned long changeIndex() const + /// Communicate data of leaf gridView + template <class DataHandle> + void communicate(DataHandle& handle, Dune::InterfaceType iftype, + Dune::CommunicationDirection dir) const { - return changeIndex_; + hostGrid_->leafGridView().communicate(handle,iftype,dir); } - private: - template <class Grid_> - static std::shared_ptr<Self> - instanceImpl(std::shared_ptr<Grid_> grid) + /// \brief 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() { - 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(); - } + return Dune::Hybrid::ifElse( + /* if */ std::is_convertible<decltype(std::declval<HG>().loadBalance()), bool>{}, + [&](auto id) { return id(hostGrid_)->loadBalance(); }, + [&](auto id) { id(hostGrid_)->loadBalance(); return true; }); } - // 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) + /// \brief 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 instanceImpl(wrap_or_share(FWD(grid))); + return Dune::Hybrid::ifElse( + /* if */ std::is_convertible<decltype(std::declval<HG>().loadBalance(handle)), bool>{}, + [&](auto id) { return id(hostGrid_)->loadBalance(handle); }, + [&](auto id) { id(hostGrid_)->loadBalance(handle); return true; }); } - // 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_; } + /// Return size of the overlap region for a given codim on the level grid view. + int overlapSize(int level, int codim) const { return hostGrid_->levelGridView(level).overlapSize(codim); } + + /// Return size of the overlap region for a given codim on the leaf grid view. + int overlapSize(int codim) const { return hostGrid_->leafGridView().overlapSize(codim); } + + /// Return size of the ghost region for a given codim on the level grid view. + int ghostSize(int level, int codim) const { return hostGrid_->levelGridView(level).ghostSize(codim); } + + /// Return size of the ghost region for a given codim on the leaf grid view. + int ghostSize(int codim) const { return hostGrid_->leafGridView().ghostSize(codim); } + + /// @} + + + /// Obtain Entity from EntitySeed of the HostGrid. + template <class EntitySeed> + auto entity(EntitySeed const& seed) const { return hostGrid_->entity(seed); } + 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() + /// 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. + /// 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_; + + template <class HostGrid> + class AdaptiveGridFamily + { + public: + struct Traits : HostGrid::Traits + { + using Grid = AdaptiveGrid<HostGrid>; + using LeafGridView = Dune::GridView< AMDiS::DefaultLeafGridViewTraits<const Grid> >; + using LevelGridView = Dune::GridView< AMDiS::DefaultLevelGridViewTraits<const Grid> >; + }; + }; + + + /// Return a change index of a grid that is not an \ref AdaptiveGrid + template <class HostGrid> + unsigned long changeIndex(HostGrid const& /*hostGrid*/) + { + return 0; + } + + /// Return a change index of an \ref AdaptiveGrid + template <class HostGrid> + unsigned long changeIndex(AdaptiveGrid<HostGrid> const& grid) + { + return grid.changeIndex(); + } + + + namespace Impl + { + template <class HostGrid> + struct AdaptiveGridImpl + { + using type = AdaptiveGrid<HostGrid>; + }; + + template <class HostGrid> + struct AdaptiveGridImpl<AdaptiveGrid<HostGrid>> + { + using type = AdaptiveGrid<HostGrid>; + }; + } + + /// Always returning an \ref AdaptiveGrid. Returns the grid itself if it is + /// already an AdaptiveGrid. + template <class HostGrid> + using AdaptiveGrid_t = typename Impl::AdaptiveGridImpl<HostGrid>::type; + } // end namespace AMDiS + + +namespace Dune +{ + /// Specialization of a \ref GridFactory to \ref AdaptiveGrid. + /// Provide a generic factory class for unstructured grids. + template <class HostGrid> + class GridFactory<AMDiS::AdaptiveGrid<HostGrid> > + : public GridFactoryInterface<AMDiS::AdaptiveGrid<HostGrid> > + { + using Self = GridFactory; + using GridType = AMDiS::AdaptiveGrid<HostGrid>; + using HostGridFactory = GridFactory<HostGrid>; + + public: + + using ctype = typename GridType::ctype; + + enum { dim = GridType::dimension }; + enum { dimworld = GridType::dimensionworld }; + + template <class... Args, + Dune::disableCopyMove<Self, Args...> = 0> + GridFactory(Args&&... args) + : hostFactory_(FWD(args)...) + {} + + /// Insert a vertex into the coarse grid + void insertVertex(FieldVector<ctype,dimworld> const& pos) override + { + hostFactory_.insertVertex(pos); + } + + template <class F, class... Args> + using HasInsertElement = decltype(std::declval<F>().insertElement(std::declval<Args>()...)); + + /// Insert an element into the coarse grid + void insertElement(GeometryType const& type, + std::vector<unsigned int> const& vertices) override + { + hostFactory_.insertElement(type, vertices); + } + + using ElementParametrizationType = std::shared_ptr<VirtualFunction<FieldVector<ctype,dim>,FieldVector<ctype,dimworld> > >; + + /// Insert a parametrized element into the coarse grid + void insertElement(GeometryType const& type, + std::vector<unsigned int> const& vertices, + ElementParametrizationType const& elementParametrization) override + { + using A0 = GeometryType; + using A1 = std::vector<unsigned int>; + using A2 = ElementParametrizationType; + Hybrid::ifElse(Std::is_detected<HasInsertElement, HostGridFactory, A0,A1,A2>{}, + [&](auto id) { id(hostFactory_).insertElement(type, vertices, elementParametrization); }, + [&](auto id) { AMDiS::error_exit("insertElement() not implemented for HostGrid type."); }); + } + + template <class F, class... Args> + using HasInsertBoundarySegment = decltype(std::declval<F>().insertBoundarySegment(std::declval<Args>()...)); + + //// Insert a boundary segment + void insertBoundarySegment(std::vector<unsigned int> const& vertices) override + { + hostFactory_.insertBoundarySegment(vertices); + } + + using BoundarySegmentType = std::shared_ptr<BoundarySegment<dim,dimworld> >; + + /// Insert an arbitrarily shaped boundary segment + void insertBoundarySegment(std::vector<unsigned int> const& vertices, + BoundarySegmentType const& boundarySegment) override + { + using A0 = std::vector<unsigned int>; + using A1 = BoundarySegmentType; + Hybrid::ifElse(Std::is_detected<HasInsertBoundarySegment, HostGridFactory, A0,A1>{}, + [&](auto id) { id(hostFactory_).insertBoundarySegment(vertices, boundarySegment); }, + [&](auto id) { AMDiS::error_exit("insertBoundarySegment() not implemented for HostGrid type."); }); + } + + /// Finalize grid creation and hand over the grid +#if DUNE_VERSION_LT(DUNE_GRID,2,7) + GridType* createGrid() override + { + std::unique_ptr<HostGrid> hostGrid(hostFactory_.createGrid()); + return new GridType(std::move(hostGrid)); + } +#else + ToUniquePtr<GridType> createGrid() override + { + std::unique_ptr<HostGrid> hostGrid(hostFactory_.createGrid()); + return makeToUnique<GridType>(std::move(hostGrid)); + } +#endif + + private: + HostGridFactory hostFactory_; + }; + + + namespace Impl + { + template <class HostGrid, bool = Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v> + class BackupRestoreFacilityImpl {}; + + template <class HostGrid> + class BackupRestoreFacilityImpl<HostGrid,true> + { + using Grid = AMDiS::AdaptiveGrid<HostGrid>; + using HostBackupRestoreFacility = BackupRestoreFacility<HostGrid>; + + public: + + /// Backup the grid to file or stream + template <class Output> + static void backup(Grid const& grid, Output const& filename_or_stream) + { + HostBackupRestoreFacility::backup(*grid.hostGrid(), filename_or_stream); + } + + /// Restore the grid from file or stream + template <class Input> + static Grid* restore(Input const& filename_or_stream) + { + std::unique_ptr<HostGrid> hostGrid(HostBackupRestoreFacility::restore(filename_or_stream)); + return new Grid(std::move(hostGrid)); + } + }; + + } // end namespace Impl + + + /// Specialization of \ref BackupRestoreFacility to \ref AdaptiveGrid. + /// Facility for writing and reading grids. + template <class HostGrid> + class BackupRestoreFacility<AMDiS::AdaptiveGrid<HostGrid>> + : public Impl::BackupRestoreFacilityImpl<HostGrid> + { + public: + using Impl::BackupRestoreFacilityImpl<HostGrid>::BackupRestoreFacilityImpl; + }; + + + namespace Capabilities + { + template <class HostGrid, int codim> + struct hasEntity<AMDiS::AdaptiveGrid<HostGrid>, codim> + : hasEntity<HostGrid,codim>{}; + + template <class HostGrid, int codim> + struct hasEntityIterator<AMDiS::AdaptiveGrid<HostGrid>, codim> + : hasEntityIterator<HostGrid, codim> {}; + + template <class HostGrid> + struct isLevelwiseConforming<AMDiS::AdaptiveGrid<HostGrid> > + : isLevelwiseConforming<HostGrid> {}; + + template <class HostGrid> + struct isLeafwiseConforming<AMDiS::AdaptiveGrid<HostGrid> > + : isLeafwiseConforming<HostGrid> {}; + + template <class HostGrid> + struct hasSingleGeometryType<AMDiS::AdaptiveGrid<HostGrid> > + : hasSingleGeometryType<HostGrid> {}; + + template <class HostGrid, int codim > + struct canCommunicate<AMDiS::AdaptiveGrid<HostGrid>, codim> + : canCommunicate<HostGrid, codim> {}; + + template <class HostGrid> + struct hasBackupRestoreFacilities<AMDiS::AdaptiveGrid<HostGrid> > + : hasBackupRestoreFacilities<HostGrid> {}; + + template <class HostGrid> + struct threadSafe<AMDiS::AdaptiveGrid<HostGrid> > + : threadSafe<HostGrid> {}; + + template <class HostGrid> + struct viewThreadSafe<AMDiS::AdaptiveGrid<HostGrid> > + : viewThreadSafe<HostGrid> {}; + + } // end namespace Capabilities +} // end namespace Dune diff --git a/src/amdis/DOFVector.hpp b/src/amdis/DOFVector.hpp index 7b1e3235fb02749daf580c6bddfaa6788dabbbdc..02c80990d05678f0216193f585b29efbc276baf0 100644 --- a/src/amdis/DOFVector.hpp +++ b/src/amdis/DOFVector.hpp @@ -65,7 +65,7 @@ namespace AMDiS public: /// Constructor. Stores the shared_ptr of the basis and creates a new DataTransfer. - DOFVector(std::shared_ptr<GB> basis, + DOFVector(std::shared_ptr<GB> const& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE) : Super(basis) , Obs(basis) diff --git a/src/amdis/MeshCreator.hpp b/src/amdis/MeshCreator.hpp index 743ec9e01896399c720c618826ac5fe3ba6f40aa..432742d67497d37de570f3f1c3d24a0d38d7bdd9 100644 --- a/src/amdis/MeshCreator.hpp +++ b/src/amdis/MeshCreator.hpp @@ -43,10 +43,11 @@ namespace AMDiS enum { dimworld = Grid::dimensionworld }; using ctype = typename Grid::ctype; + using HostGrid = typename AdaptiveGrid_t<Grid>::HostGrid; /// Construct a new MeshCreator /** - * \param name The name of the mesh used in the initifile + * \param name The name of the mesh used in the initfile **/ MeshCreator(std::string const& name) : name_(name) @@ -64,11 +65,11 @@ namespace AMDiS * * Otherwise tries to create a grid depending on the grid type. **/ - std::shared_ptr<AdaptiveGrid<Grid>> create() const + std::shared_ptr<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; + std::unique_ptr<HostGrid> gridPtr; if (filename) { // read a macro file gridPtr = create_unstructured_grid(filename.value()); @@ -83,26 +84,26 @@ namespace AMDiS } } else { // decide by inspecting the grid type how to create the grid - gridPtr = create_by_gridtype<Grid>(Dune::PriorityTag<42>{}); + gridPtr = create_by_gridtype<HostGrid>(Dune::PriorityTag<42>{}); } - return AdaptiveGrid<Grid>::instance(gridPtr); + return construct(std::move(gridPtr)); } /// Create a structured cube grid - std::unique_ptr<Grid> create_cube_grid() const + std::unique_ptr<HostGrid> create_cube_grid() const { return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells) { - return Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, numCells); + return Dune::StructuredGridFactory<HostGrid>::createCubeGrid(lower, upper, numCells); }); } /// Create a structured simplex grid - std::unique_ptr<Grid> create_simplex_grid() const + std::unique_ptr<HostGrid> create_simplex_grid() const { return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells) { - return Dune::StructuredGridFactory<Grid>::createSimplexGrid(lower, upper, numCells); + return Dune::StructuredGridFactory<HostGrid>::createSimplexGrid(lower, upper, numCells); }); } @@ -119,9 +120,20 @@ namespace AMDiS } private: + static std::shared_ptr<Grid> construct(std::unique_ptr<Grid> hostGrid) + { + return std::move(hostGrid); + } + + template <class HG> + static std::shared_ptr<Grid> construct(std::unique_ptr<HG> hostGrid) + { + return std::make_shared<Grid>(std::move(hostGrid)); + } + // use the structured grid factory to create the grid template <class Size = unsigned int, class Factory> - std::unique_ptr<Grid> create_structured_grid(Factory factory) const + std::unique_ptr<HostGrid> create_structured_grid(Factory factory) const { // Lower left corner of the domain Dune::FieldVector<ctype,int(dimworld)> lower(0); @@ -137,16 +149,16 @@ namespace AMDiS } // read a filename from `[meshName]->macro file name` and determine from the extension the fileformat - std::unique_ptr<Grid> create_unstructured_grid(std::string const& filename) const + std::unique_ptr<HostGrid> create_unstructured_grid(std::string const& filename) const { filesystem::path fn(filename); auto ext = fn.extension(); if (ext == ".msh") { - return read_gmsh_file<Grid>(filename, Dune::PriorityTag<42>{}); + return read_gmsh_file<HostGrid>(filename, Dune::PriorityTag<42>{}); } else if (ext == ".1d" || ext == ".2d" || ext == ".3d" || ext == ".amc") { - return read_alberta_file<Grid>(filename, Dune::PriorityTag<42>{}); + return read_alberta_file<HostGrid>(filename, Dune::PriorityTag<42>{}); } else { error_exit("Cannot read grid file. Unknown file extension."); @@ -160,7 +172,7 @@ namespace AMDiS std::declval<std::shared_ptr<Dune::BoundarySegment<GridType::dimension, GridType::dimensionworld> >>()) ); // use GmshReader if GridFactory supports insertBoundarySegments - template <class GridType = Grid, + template <class GridType = HostGrid, REQUIRES(Dune::Std::is_detected<SupportsGmshReader, GridType>::value)> std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<1>) const { @@ -169,7 +181,7 @@ namespace AMDiS } // fallback if GmshReader cannot be used - template <class GridType = Grid> + template <class GridType = HostGrid> std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<0>) const { error_exit("Gmsh reader not supported for this grid type!"); @@ -183,7 +195,7 @@ namespace AMDiS using IsAlbertaGrid = decltype(std::declval<GridType>().alberta2dune(0,0)); // construct the albertagrid directly from a filename - template <class GridType = Grid, + template <class GridType = HostGrid, REQUIRES(GridType::dimensionworld == DIM_OF_WORLD), REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)> std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<3>) const @@ -192,7 +204,7 @@ namespace AMDiS } // use a gridfactory and the generic AlbertaReader - template <class GridType = Grid, + template <class GridType = HostGrid, REQUIRES(GridType::dimensionworld == DIM_OF_WORLD)> std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<2>) const { @@ -205,7 +217,7 @@ namespace AMDiS } // error if WORLDDIM not the same as Grid::dimensionworld - template <class GridType = Grid, + template <class GridType = HostGrid, REQUIRES(GridType::dimensionworld != DIM_OF_WORLD)> std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<1>) const { @@ -225,19 +237,19 @@ namespace AMDiS #if HAVE_ALBERTA // albertagrid -> simplex - template <class GridType = Grid, + template <class GridType = HostGrid, REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)> std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<3>) const { return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells) { - return MacroGridFactory<Grid>::createSimplexGrid(lower, upper, numCells); + return MacroGridFactory<GridType>::createSimplexGrid(lower, upper, numCells); }); } #endif // yasp grid -> cube - template <class GridType = Grid, + template <class GridType = HostGrid, class = typename GridType::YGrid> std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<2>) const { @@ -250,27 +262,27 @@ namespace AMDiS } template <int dim, class ct> - std::unique_ptr<Grid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantCoordinates<ct,dim>>>, + std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantCoordinates<ct,dim>>>, int overlap, std::bitset<dimension> const& per) const { return create_structured_grid<int>([&](auto&& /*lower*/, auto&& upper, std::array<int,dimension> const& numCells) { - return std::make_unique<Grid>(upper, numCells, per, overlap); + return std::make_unique<HostGrid>(upper, numCells, per, overlap); }); } template <int dim, class ct> - std::unique_ptr<Grid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<ct,dim>>>, + std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<ct,dim>>>, int overlap, std::bitset<dimension> const& per) const { return create_structured_grid<int>([&](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells) { - return std::make_unique<Grid>(lower, upper, numCells, per, overlap); + return std::make_unique<HostGrid>(lower, upper, numCells, per, overlap); }); } template <int dim, class ct> - std::unique_ptr<Grid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::TensorProductCoordinates<ct,dim>>>, + std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::TensorProductCoordinates<ct,dim>>>, int overlap, std::bitset<dimension> const& per) const { error_exit("MeshCreator cannot create YaspGrid with TensorProductCoordinates."); @@ -280,7 +292,7 @@ namespace AMDiS #if HAVE_DUNE_SPGRID // spgrid -> cube - template <class GridType = Grid, + template <class GridType = HostGrid, class = typename GridType::ReferenceCube, class = typename GridType::MultiIndex> std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<1>) const @@ -294,7 +306,7 @@ namespace AMDiS #endif // final fallback - template <class GridType = Grid> + template <class GridType = HostGrid> std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<0>) const { error_exit("Don't know how to create the grid."); diff --git a/src/amdis/Observer.hpp b/src/amdis/Observer.hpp index 62383a9f5fe7a5737abd733bffffd117388db050..3bf146e71fe00c5d751c99bec14e2e8ba5cafcfc 100644 --- a/src/amdis/Observer.hpp +++ b/src/amdis/Observer.hpp @@ -5,6 +5,8 @@ #include <memory> #include <type_traits> +#include <dune/common/typeutilities.hh> + #include <amdis/common/ConceptsBase.hpp> #include <amdis/Output.hpp> @@ -12,135 +14,128 @@ namespace AMDiS { namespace Impl { - // Forward declaration + // forward declaration template <class Event> - class SignalInterface; + class ObserverInterface; template <class Event> - class ObserverInterface + class SignalBase { 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) + /// Attaches an observer to this class. This method will be called by all observers with + /// with themselves as argument. + void attachObserver(ObserverInterface<Event>* o) const { - subject->attachObserver(this); + observers_.push_back(o); } - /// 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) + /// Detaches an observer to this class. This method will be called by all observers with + /// with themselves as argument. + void detachObserver(ObserverInterface<Event>* o) const { - subject->detachObserver(this); + auto it = std::find(observers_.begin(), observers_.end(), o); + if (it != observers_.end()) + observers_.erase(it); } - /// This method will be called by the subject when triggering the event. - virtual void update(Event const&) + /// Notify all observers that have called attachObserver but have not called detachObserver + void notify(Event const& e) const { - error_exit("Method must be overridden by derived class"); + for (ObserverInterface<Event>* o : observers_) + o->update(e); } - }; - - 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; + private: + /// List of observers that need to be notified in case of an event + // NOTE: this list is mutable, since the notification list itself is not part + // of the internal state of the object signaling the event. + mutable std::list<ObserverInterface<Event>*> observers_; }; + template <class Event> - class SignalsImpl - : virtual public SignalInterface<Event> + class ObserverInterface { public: - virtual ~SignalsImpl() = default; + virtual ~ObserverInterface() = default; - void attachObserver(ObserverInterface<Event>* o) override + /// 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<SignalBase<Event> const> const& subject) { - observers_.push_back(o); + if (bool(subject)) + subject->attachObserver(this); } - void detachObserver(ObserverInterface<Event>* o) override + /// 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<SignalBase<Event> const> const& subject) { - auto it = std::find(observers_.begin(), observers_.end(), o); - if (it != observers_.end()) - observers_.erase(it); + if (bool(subject)) + subject->detachObserver(this); } - void notify(Event const& e) const override + /// This method will be called by the subject when triggering the event. + virtual void update(Event const&) { - for (ObserverInterface<Event>* o : observers_) - o->update(e); + error_exit("Method must be overridden by derived class"); } - - 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 + class ObserverBase : virtual public ObserverInterface<Event> { - using Self = ObserverImpl; + using Self = ObserverBase; - ObserverImpl() - : subject_(std::make_shared<SignalsDummy<Event>>()) - {} + private: + ObserverBase() = default; - public: - ObserverImpl(std::shared_ptr<SignalInterface<Event>> subject) + // use subject if Event is handled directly + template <class S, + REQUIRES(std::is_base_of<SignalBase<Event>, S>::value)> + ObserverBase(std::shared_ptr<S const> subject, Dune::PriorityTag<2>) : 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()) + // get next subject in hierarchy if Event is not handled + template <class Observer, + class = void_t<decltype(std::declval<Observer>().subject_)> > + ObserverBase(std::shared_ptr<Observer const> const& o, Dune::PriorityTag<1>) + : ObserverBase(o->subject_) + {} + + // non-observable type + template <class T> + ObserverBase(std::shared_ptr<T const> const& other, Dune::PriorityTag<0>) + { + /* fallback implementation */ + } + + public: + template <class T> + ObserverBase(std::shared_ptr<T const> const& arg) + : ObserverBase(arg, Dune::PriorityTag<42>{}) {} /// Copy constructor - ObserverImpl(Self const& that) + ObserverBase(Self const& that) : subject_(that.subject_) { this->attach(subject_); } /// Move constructor - ObserverImpl(Self&& that) - : ObserverImpl() + ObserverBase(Self&& that) { swap(*this, that); } /// Destructor. Detaches observer from the subject - virtual ~ObserverImpl() + ~ObserverBase() override { this->detach(subject_); } @@ -165,11 +160,14 @@ namespace AMDiS private: /// The observed subject - std::shared_ptr<SignalInterface<Event>> subject_; + std::shared_ptr<SignalBase<Event> const> subject_ = nullptr; }; - } - /** Mixin for signaling of certain events. + } // end namespace Impl + + + /// \brief 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, @@ -179,54 +177,56 @@ namespace AMDiS */ template <class Event, class... Events> class Signals - : public Impl::SignalsImpl<Event> + : public Impl::SignalBase<Event> , public Signals<Events...> { public: - using Impl::SignalsImpl<Event>::notify; + using Impl::SignalBase<Event>::notify; using Signals<Events...>::notify; }; template <class Event> class Signals<Event> - : public Impl::SignalsImpl<Event> + : public Impl::SignalBase<Event> { public: - using Impl::SignalsImpl<Event>::notify; + using Impl::SignalBase<Event>::notify; }; - /** Mixin for reacting to certain events - * Derived classes T can react to events by implementing the functions update(Event& e) for + /// \brief Mixin for reacting to certain events. + /** + * Derived classes can react to events by implementing the functions `update(Event)` 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. + * `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. + * - the instance of class `Subject` passed to the constructor + * - the instance passed to the constructor of `Subject::Observer<S>` if Subject inherits from + * \ref Observer<S> + * - all other instances indirectly accessible in the above way + * + * For each event E the first object in the above hierarchy that implements \ref Signals<Es...> with E + * included in Es... will be observed by this class. */ template <class Subject, class... Events> class Observer - : public Impl::ObserverImpl<Events>... + : public Impl::ObserverBase<Events>... { template <class E> - friend class Impl::ObserverImpl; + friend class Impl::ObserverBase; public: - Observer(std::shared_ptr<Subject> s) - : Impl::ObserverImpl<Events>(s)... + Observer(std::shared_ptr<Subject const> const& s) + : Impl::ObserverBase<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_; } + Observer(std::shared_ptr<Subject> const& s) + : Observer(std::const_pointer_cast<Subject const>(s)) + {} - std::shared_ptr<Subject> subject_; + private: + std::shared_ptr<Subject const> subject_ = nullptr; }; } // end namespace AMDiS diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp index e44a627aa8cdd29ded5be325c862e205dd7ddf79..b50009b1bb24f000c20a981e662b8d69d710bb26 100644 --- a/src/amdis/ProblemStat.hpp +++ b/src/amdis/ProblemStat.hpp @@ -66,9 +66,7 @@ namespace AMDiS using GlobalBasis = typename Traits::GlobalBasis; using GridView = typename GlobalBasis::GridView; - // 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 Grid = AdaptiveGrid_t<typename GridView::Grid>; 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>; @@ -102,13 +100,13 @@ namespace AMDiS ProblemStat(std::string const& name, Grid_&& grid) : ProblemStat(name) { - adoptGrid(Grid::instance(FWD(grid))); + adoptGrid(wrap_or_share(FWD(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) + template <class Grid_, class Basis_> + ProblemStat(std::string const& name, Grid_&& grid, Basis_&& globalBasis) : ProblemStat(name, FWD(grid)) { adoptGlobalBasis(wrap_or_share(FWD(globalBasis))); @@ -161,7 +159,7 @@ namespace AMDiS * Adds an operator to the list of boundary operators to be assembled in * quadrature points on the boundary intersections. * - * \param b Boundary indentifier where to assemble this operator. Can be + * \param b Boundary identifier where to assemble this operator. Can be * constructed from an integer. \see BoundaryType * \param op A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator * \param row TreePath identifying the sub-basis in the global basis tree @@ -212,7 +210,7 @@ namespace AMDiS * Adds an operator to the list of boundary operators to be assembled in * quadrature points on the boundary intersections. * - * \param b Boundary indentifier where to assemble this operator. Can be + * \param b Boundary identifier where to assemble this operator. Can be * constructed from an integer. \see BoundaryType * \param op A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator * \param path TreePath identifying the sub-basis in the global basis tree @@ -373,60 +371,40 @@ namespace AMDiS public: // set-methods /// Set a new linear solver for the problem - void setSolver(std::shared_ptr<LinearSolverType> solver) + template <class Solver_> + void setSolver(Solver_&& solver) { - linearSolver_ = solver; - } - - /// Wrap the solver reference into a non-destroying shared_ptr, or move it - /// into a new shared_ptr if it is a temporary. - template <class S, - REQUIRES(std::is_base_of<LinearSolverType, remove_cvref_t<S>>::value)> - void setSolver(S&& solver) - { - setSolver(Dune::wrap_or_move(FWD(solver))); + linearSolver_ = wrap_or_share(FWD(solver)); } /// Set the grid. Stores pointer and initializes feSpaces /// matrices and vectors, as well as markers and file-writers. - void setGrid(std::shared_ptr<Grid> grid) + /// If grid is given as reference, wrap it into a non-destroying shared_ptr + template <class Grid_> + void setGrid(Grid_&& grid) { - adoptGrid(std::move(grid)); + adoptGrid(wrap_or_share(FWD(grid))); createGlobalBasis(); createMatricesAndVectors(); createMarker(); createFileWriter(); } - /// Wrap the grid into a non-destroying shared_ptr - template <class Grid_> - void setGrid(Grid_&& grid) - { - setGrid(Grid::instance(FWD(grid))); - } - /// Store the shared_ptr and the name of the marker in the problem /** * Note: multiple markers can be added but must have different names **/ - void addMarker(std::shared_ptr<Marker<Grid>> marker) + template <class Marker_> + void addMarker(Marker_&& m) { - auto it = marker_.emplace(marker->name(), std::move(marker)); + auto marker = wrap_or_share(FWD(m)); + auto it = marker_.emplace(marker->name(), marker); if (marker_.size() > 1) it.first->second->setMaximumMarking(true); } - /// Wrap the reference into a non-destroying shared_ptr or move it into a - /// new shared_ptr if it is a temporary. - template <class M, - REQUIRES(std::is_base_of<Marker<Grid>, remove_cvref_t<M>>::value)> - void addMarker(M&& marker) - { - addMarker(Dune::wrap_or_move(FWD(marker))); - } - /// Remove a marker with the given name from the problem void removeMarker(std::string name) { @@ -456,7 +434,7 @@ namespace AMDiS initGlobalBasis(); } - void adoptGrid(std::shared_ptr<Grid> grid, + void adoptGrid(std::shared_ptr<Grid> const& grid, std::shared_ptr<BoundaryManager<Grid>> const& boundaryManager) { grid_ = grid; @@ -464,8 +442,14 @@ namespace AMDiS Parameters::get(name_ + "->mesh", gridName_); } - void adoptGrid(std::shared_ptr<Grid> grid) + void adoptGrid(std::shared_ptr<Grid> const& grid) + { + adoptGrid(grid, std::make_shared<BoundaryManager<Grid>>(grid)); + } + + void adoptGrid(std::shared_ptr<typename Grid::HostGrid> const& hostGrid) { + auto grid = std::make_shared<Grid>(hostGrid); adoptGrid(grid, std::make_shared<BoundaryManager<Grid>>(grid)); } diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp index 163d9413c5fae1dee13723d7044b6d43f538e592..ab4699a61c3ff29589069d455a847933381b2bd0 100644 --- a/src/amdis/ProblemStat.inc.hpp +++ b/src/amdis/ProblemStat.inc.hpp @@ -136,11 +136,14 @@ restore(Flag initFlag) test_exit(filesystem::exists(grid_filename), "Restore file '{}' not found.", grid_filename); test_exit(filesystem::exists(solution_filename), "Restore file '{}' not found.", solution_filename); + // TODO(SP): implement BAckupRestore independent of wrapped grid + using HostGrid = typename Grid::HostGrid; + // restore grid from file if (Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v) - adoptGrid(Grid::instance(std::unique_ptr<HostGrid>(Dune::BackupRestoreFacility<HostGrid>::restore(grid_filename)))); + adoptGrid(std::shared_ptr<HostGrid>(Dune::BackupRestoreFacility<HostGrid>::restore(grid_filename))); else - adoptGrid(Grid::instance(std::unique_ptr<HostGrid>(BackupRestoreByGridFactory<HostGrid>::restore(grid_filename)))); + adoptGrid(std::shared_ptr<HostGrid>(BackupRestoreByGridFactory<HostGrid>::restore(grid_filename))); // create fespace if (initFlag.isSet(INIT_FE_SPACE) || initFlag.isSet(INIT_SYSTEM)) @@ -173,7 +176,8 @@ template <class Traits> void ProblemStat<Traits>::createGrid() { Parameters::get(name_ + "->mesh", gridName_); - MeshCreator<HostGrid> creator(gridName_); + + MeshCreator<Grid> creator(gridName_); grid_ = creator.create(); Dune::Timer t; @@ -211,8 +215,7 @@ template <class Traits> void ProblemStat<Traits>::createGlobalBasisImpl(std::true_type) { assert( bool(grid_) ); - // TODO(FM): Replace HostGrid -> Grid once AdaptiveGrid is full metagrid - static_assert(std::is_same<GridView, typename HostGrid::LeafGridView>::value, ""); + static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, ""); auto basis = Traits::create(name_, grid_->leafGridView()); globalBasis_ = std::make_shared<GlobalBasis>(std::move(basis)); } diff --git a/src/amdis/ProblemStatTraits.hpp b/src/amdis/ProblemStatTraits.hpp index 87393b725798540abd9fd9e02645d196588e86dc..b5a297f9d3741f10a7b1fad091bfe21edad9e911 100644 --- a/src/amdis/ProblemStatTraits.hpp +++ b/src/amdis/ProblemStatTraits.hpp @@ -67,10 +67,12 @@ namespace AMDiS }; - template <class Grid, class PreBasisCreator, class T = double> + template <class HostGrid, class PreBasisCreator, class T = double> struct DefaultBasisCreator { + using Grid = AdaptiveGrid<HostGrid>; using GridView = typename Grid::LeafGridView; + static auto create(std::string const& name, GridView const& gridView) { return makeGlobalBasis(name, gridView, PreBasisCreator::create()); diff --git a/src/amdis/common/CMakeLists.txt b/src/amdis/common/CMakeLists.txt index 32b15cb34d62007698ba72b1844924f8b0c2767d..8f345538998147c05d1ff8860127b303f2d565b4 100644 --- a/src/amdis/common/CMakeLists.txt +++ b/src/amdis/common/CMakeLists.txt @@ -10,6 +10,7 @@ install(FILES Concepts.hpp ConceptsBase.hpp ConcurrentCache.hpp + DefaultGridView.hpp DerivativeTraits.hpp FakeContainer.hpp FieldMatVec.hpp diff --git a/src/amdis/common/DefaultGridView.hpp b/src/amdis/common/DefaultGridView.hpp new file mode 100644 index 0000000000000000000000000000000000000000..a620e3ea78858de4213e416be578fdc5c08ad0bc --- /dev/null +++ b/src/amdis/common/DefaultGridView.hpp @@ -0,0 +1,318 @@ +#pragma once + +#include <dune/common/typetraits.hh> +#include <dune/common/exceptions.hh> + +#include <dune/grid/common/capabilities.hh> +#include <dune/grid/common/gridview.hh> + +namespace AMDiS +{ + /* + * NOTE: this is a modification of dune/grid/common/defaultgridview.hh since we + * want all implementation specific methods to be put into the grid class and not + * into the entity class. See ibegin(), iend(). + */ + + template <class GridImp> + class DefaultLevelGridView; + + template <class GridImp> + class DefaultLeafGridView; + + + template <class GridImp> + struct DefaultLevelGridViewTraits + { + using GridViewImp = DefaultLevelGridView<GridImp>; + + /// type of the grid + using Grid = std::remove_const_t<GridImp>; + + /// type of the index set + using IndexSet = typename Grid::Traits::LevelIndexSet; + + /// type of the intersection + using Intersection = typename Grid::Traits::LevelIntersection; + + /// type of the intersection iterator + using IntersectionIterator = typename Grid::Traits::LevelIntersectionIterator; + + /// type of the collective communication + using CollectiveCommunication = typename Grid::Traits::CollectiveCommunication; + + template <int cd> + struct Codim + { + using Entity = typename Grid::Traits::template Codim<cd>::Entity; + using Geometry = typename Grid::template Codim<cd>::Geometry; + using LocalGeometry = typename Grid::template Codim<cd>::LocalGeometry; + + /// Define types needed to iterate over entities of a given partition type + template <Dune::PartitionIteratorType pit> + struct Partition + { + /// iterator over a given codim and partition type + using Iterator = typename Grid::template Codim<cd>::template Partition<pit>::LevelIterator; + }; + + using Iterator = typename Partition<Dune::All_Partition>::Iterator; + }; + + enum { conforming = Dune::Capabilities::isLevelwiseConforming<Grid>::v }; + }; + + + template <class GridImp> + class DefaultLevelGridView + { + public: + using Traits = DefaultLevelGridViewTraits<GridImp>; + using Grid = typename Traits::Grid; + using IndexSet = typename Traits::IndexSet; + using Intersection = typename Traits::Intersection; + using IntersectionIterator = typename Traits::IntersectionIterator; + using CollectiveCommunication = typename Traits::CollectiveCommunication; + + template <int cd> + using Codim = typename Traits::template Codim<cd>; + + enum { conforming = Traits::conforming }; + + public: + DefaultLevelGridView(Grid const& grid, int level) + : grid_(&grid) + , level_(level) + {} + + /// Obtain a const reference to the underlying hierarchic grid + Grid const& grid() const + { + assert(grid_); + return *grid_; + } + + /// Obtain the index set + IndexSet const& indexSet() const + { + return grid().levelIndexSet(level_); + } + + /// Obtain number of entities in a given codimension + int size(int codim) const + { + return grid().size(level_, codim); + } + + /// Obtain number of entities with a given geometry type + int size(Dune::GeometryType const& type) const + { + return grid().size(level_, type); + } + + /// Obtain begin iterator for this view + template <int cd, Dune::PartitionIteratorType pit = Dune::All_Partition> + typename Codim<cd>::template Partition<pit>::Iterator begin() const + { + return grid().template lbegin<cd, pit>(level_); + } + + /// Obtain end iterator for this view + template <int cd, Dune::PartitionIteratorType pit = Dune::All_Partition> + typename Codim<cd>::template Partition<pit>::Iterator end() const + { + return grid().template lend<cd, pit>(level_); + } + + /// Obtain begin intersection iterator with respect to this view + IntersectionIterator ibegin(typename Codim<0>::Entity const& entity) const + { + return grid().ilevelbegin(entity); + } + + /// Obtain end intersection iterator with respect to this view + IntersectionIterator iend(typename Codim<0>::Entity const& entity) const + { + return grid().ilevelend(entity); + } + + /// Obtain collective communication object + CollectiveCommunication const& comm() const + { + return grid().comm(); + } + + /// Return size of the overlap region for a given codim on the grid view. + int overlapSize(int codim) const + { + return grid().overlapSize(level_, codim); + } + + /// Return size of the ghost region for a given codim on the grid view. + int ghostSize(int codim) const + { + return grid().ghostSize(level_, codim); + } + + /// Communicate data on this view + template <class DataHandleImp, class DataType> + void communicate(Dune::CommDataHandleIF<DataHandleImp, DataType>& data, + Dune::InterfaceType iftype, + Dune::CommunicationDirection dir) const + { + return grid().communicate(data, iftype, dir, level_); + } + + private: + Grid const* grid_; + int level_; + }; + + + template <class GridImp> + struct DefaultLeafGridViewTraits + { + using GridViewImp = DefaultLeafGridView<GridImp>; + + /// type of the grid + using Grid = std::remove_const_t<GridImp>; + + /// type of the index set + using IndexSet = typename Grid::Traits::LeafIndexSet; + + /// type of the intersection + using Intersection = typename Grid::Traits::LeafIntersection; + + /// type of the intersection iterator + using IntersectionIterator = typename Grid::Traits::LeafIntersectionIterator; + + /// type of the collective communication + using CollectiveCommunication = typename Grid::Traits::CollectiveCommunication; + + template <int cd> + struct Codim + { + using Entity = typename Grid::Traits::template Codim<cd>::Entity; + using Geometry = typename Grid::template Codim<cd>::Geometry; + using LocalGeometry = typename Grid::template Codim<cd>::LocalGeometry; + + /// Define types needed to iterate over entities of a given partition type + template <Dune::PartitionIteratorType pit> + struct Partition + { + /// iterator over a given codim and partition type + using Iterator = typename Grid::template Codim<cd>::template Partition<pit>::LeafIterator; + }; + + using Iterator = typename Partition<Dune::All_Partition>::Iterator; + }; + + enum { conforming = Dune::Capabilities::isLeafwiseConforming<Grid>::v }; + }; + + + template <class GridImp> + class DefaultLeafGridView + { + public: + using Traits = DefaultLeafGridViewTraits<GridImp>; + using Grid = typename Traits::Grid; + using IndexSet = typename Traits::IndexSet; + using Intersection = typename Traits::Intersection; + using IntersectionIterator = typename Traits::IntersectionIterator; + using CollectiveCommunication = typename Traits::CollectiveCommunication; + + template <int cd> + using Codim = typename Traits::template Codim<cd>; + + enum { conforming = Traits::conforming }; + + public: + DefaultLeafGridView(Grid const& grid) + : grid_(&grid) + {} + + /// obtain a const reference to the underlying hierarchic grid + Grid const& grid() const + { + assert(grid_); + return *grid_; + } + + /// Obtain the index set + IndexSet const& indexSet() const + { + return grid().leafIndexSet(); + } + + /// Obtain number of entities in a given codimension + int size(int codim) const + { + return grid().size(codim); + } + + /// Obtain number of entities with a given geometry type + int size(Dune::GeometryType const& type) const + { + return grid().size(type); + } + + + /// Obtain begin iterator for this view + template <int cd, Dune::PartitionIteratorType pit = Dune::All_Partition> + typename Codim<cd>::template Partition<pit>::Iterator begin() const + { + return grid().template leafbegin<cd, pit>(); + } + + /// Obtain end iterator for this view + template <int cd, Dune::PartitionIteratorType pit = Dune::All_Partition> + typename Codim<cd>::template Partition<pit>::Iterator end() const + { + return grid().template leafend<cd, pit>(); + } + + /// Obtain begin intersection iterator with respect to this view + IntersectionIterator ibegin(typename Codim<0>::Entity const& entity) const + { + return grid().ileafbegin(entity); + } + + /// Obtain end intersection iterator with respect to this view + IntersectionIterator iend(typename Codim<0>::Entity const& entity) const + { + return grid().ileafend(entity); + } + + /// Obtain collective communication object + CollectiveCommunication const& comm() const + { + return grid().comm(); + } + + /// Return size of the overlap region for a given codim on the grid view. + int overlapSize(int codim) const + { + return grid().overlapSize(codim); + } + + /// Return size of the ghost region for a given codim on the grid view. + int ghostSize(int codim) const + { + return grid().ghostSize(codim); + } + + /// Communicate data on this view + template <class DataHandleImp, class DataType> + void communicate(Dune::CommDataHandleIF<DataHandleImp, DataType>& data, + Dune::InterfaceType iftype, + Dune::CommunicationDirection dir) const + { + return grid().communicate(data, iftype, dir); + } + + private: + Grid const* grid_; + }; + +} // end namespace AMDiS diff --git a/src/amdis/functions/ParallelGlobalBasis.hpp b/src/amdis/functions/ParallelGlobalBasis.hpp index 81ca2ab80ec301e6f68da710681a33f69089bc26..47f278732d0223d55d56e775e3ae3478ebab9a96 100644 --- a/src/amdis/functions/ParallelGlobalBasis.hpp +++ b/src/amdis/functions/ParallelGlobalBasis.hpp @@ -50,11 +50,8 @@ namespace AMDiS template <class B> struct basisUpdate { - basisUpdate(B const& b_) - : basis(b_) - {} - - B const& basis; + using Basis = B; + Basis const& basis; }; } @@ -68,11 +65,11 @@ namespace AMDiS */ template <class PB> class ParallelGlobalBasis - : public Observer<AdaptiveGrid<typename PB::GridView::Grid>, event::adapt> - , public Signals<event::basisUpdate<ParallelGlobalBasis<PB>>> + : public Observer<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 Subject = typename PB::GridView::Grid; + using Super = Observer<Subject, event::adapt>; using Self = ParallelGlobalBasis<PB>; struct DummyImpl {}; @@ -87,8 +84,7 @@ namespace AMDiS /// 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; + using Grid = typename GridView::Grid; /// Type used for global numbering of the basis vectors using MultiIndex = typename PreBasis::MultiIndex; @@ -106,7 +102,7 @@ namespace AMDiS /// Node index set provided by PreBasis using NodeIndexSet = typename PreBasis::template IndexSet<PrefixPath>; - /// Type of local indixes set exported by localIndexSet() + /// Type of local index set exported by localIndexSet() using LocalIndexSet = Dune::Functions::DefaultLocalIndexSet<LocalView, NodeIndexSet>; #endif @@ -115,51 +111,48 @@ namespace AMDiS using ADH = Dune::AdaptDataHandle<Grid, DummyImpl>; + + public: + + /// Construct this global basis with given name and grid, and constructing a preBasis. /** - * \brief Constructor - * - * \tparam T Argument list for PreBasis - * \param s The AdaptiveGrid providing the GridView for this basis - * \param t Argument list for PreBasis + * \param grid The Grid providing the GridView for this basis + * \param args... 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_() + template <class... Args, + Dune::Functions::enableIfConstructible<PreBasis, Args...> = 0> + ParallelGlobalBasis(std::string const& name, Grid const& grid, Args&&... args) + : Super(Dune::stackobject_to_shared_ptr(grid)) + , preBasis_(FWD(args)...) { 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)...) + /// Construct this global basis with empty name + template <class... Args, + Dune::Functions::enableIfConstructible<PreBasis, Args...> = 0> + ParallelGlobalBasis(Grid const& grid, Args&&... args) + : ParallelGlobalBasis(std::string(""), grid, FWD(args)...) {} - /** Converting constructor from Dune::DefaultGlobalBasis<PB>. + /// 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); - } + : ParallelGlobalBasis(std::string(""), from.gridView().grid(), from.preBasis()) + {} + + + public: /// Obtain the grid view that the basis is defined on GridView const& gridView() const @@ -182,17 +175,15 @@ namespace AMDiS /// Updates the underlying basis when event::adapt is triggered by the observed grid void update(event::adapt const& e) override { - if (e.adapted) - { + if (e.adapted) { update(gridView()); event::basisUpdate<Self> eNew{*this}; this->notify(eNew); } } + /// \brief Update the stored grid view /** - * \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. */ @@ -216,7 +207,7 @@ namespace AMDiS } /// Return number of possible values for next position in multi index - size_type size(const SizePrefix& prefix) const + size_type size(SizePrefix const& prefix) const { return preBasis_.size(prefix); } @@ -239,12 +230,15 @@ namespace AMDiS return prefixPath_; } - /// Return the communicator. - /// This provides the means to communicate data associated to the basis with other processes. + /// \brief 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 + ADH globalRefineCallback() const { // TODO(FM): Implement error_exit("Not implemented: ParallelGlobalBasis::globalRefineCallback()"); @@ -253,41 +247,34 @@ namespace AMDiS protected: PreBasis preBasis_; - PrefixPath prefixPath_; + 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) + + template <class MultiIndex, class GV, class PBF> + auto makeGlobalBasis(std::string const& name, GV const& gridView, PBF&& 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)); + return ParallelGlobalBasis<TYPEOF(preBasis)>(name, gridView.grid(), std::move(preBasis)); } - template<class GridView, class PreBasisFactory> - auto makeGlobalBasis(std::string const& name, GridView const& gridView, - PreBasisFactory&& preBasisFactory) + template <class GV, class PBF> + auto makeGlobalBasis(std::string const& name, GV const& gridView, PBF&& preBasisFactory) { - using RawPreBasisFactory = std::decay_t<PreBasisFactory>; + using RawPreBasisFactory = remove_cvref_t<PBF>; 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)); + return makeGlobalBasis<MultiIndex, GV, PBF>(name, gridView, FWD(preBasisFactory)); } - template<class GridView, class PreBasisFactory> - auto makeGlobalBasis(GridView const& gridView, PreBasisFactory&& preBasisFactory) + template <class GV, class PBF> + auto makeGlobalBasis(GV const& gridView, PBF&& preBasisFactory) { - return makeGlobalBasis("", gridView, FWD(preBasisFactory)); + return makeGlobalBasis(std::string(""), gridView, FWD(preBasisFactory)); } } // end namespace AMDiS diff --git a/src/amdis/io/FileWriterCreator.hpp b/src/amdis/io/FileWriterCreator.hpp index 52f62913a21c5014bd55608d69de465409c1db73..1193f19b924e2c56e835c715592cd26212666fe1 100644 --- a/src/amdis/io/FileWriterCreator.hpp +++ b/src/amdis/io/FileWriterCreator.hpp @@ -24,8 +24,7 @@ namespace AMDiS { using GlobalBasis = typename SystemVector::GlobalBasis; using GridView = typename GlobalBasis::GridView; - using Grid = typename GlobalBasis::Grid; - //using Grid = typename GridView::Grid; // TODO(FM): Use once AdaptiveGrid is complete Metagrid + using Grid = typename GridView::Grid; public: /// Constructor. Stores the pointer to the systemVector and to the (optional) boundaryManager. diff --git a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp index 29e493b6c32f2dbf8048d00bbc883797735637cd..df56ecafd0ad12d7e6fd30d3bf37f4dc81b66be0 100644 --- a/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp +++ b/src/amdis/linearalgebra/petsc/MatrixNnzStructure.inc.hpp @@ -17,8 +17,7 @@ namespace AMDiS { template <class Basis> void MatrixNnzStructure::update(Basis const& basis) { - using Grid = AdaptiveGrid<typename Basis::GridView::Grid>; - unsigned long newChangeIndex = Grid::instance(basis.gridView().grid())->changeIndex(); + unsigned long newChangeIndex = changeIndex(basis.gridView().grid()); bool gridChanged = newChangeIndex > changeIndex_; changeIndex_ = newChangeIndex; diff --git a/test/AdaptiveGridTest.cpp b/test/AdaptiveGridTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4436e81e8820601cec2d8ea96304fc83aa95949c --- /dev/null +++ b/test/AdaptiveGridTest.cpp @@ -0,0 +1,38 @@ +#include "config.h" +#include <dune/grid/yaspgrid.hh> +#include <dune/grid/test/gridcheck.hh> + +#include <amdis/AdaptiveGrid.hpp> + +#include "Tests.hpp" + +using namespace AMDiS; + +// test IdentityGrid for given dimension +template <int dim> +void testDim() +{ + using HostGrid = Dune::YaspGrid<dim>; + std::array<int,dim> n; + std::fill(n.begin(), n.end(), 1 << (5 - dim)); + Dune::FieldVector<double,dim> extension(1.0); + + HostGrid grid(extension,n); + grid.globalRefine(1); + + AdaptiveGrid<HostGrid> adaptiveGrid(grid); + adaptiveGrid.globalRefine(1); + + gridcheck(adaptiveGrid); +} + + +int main(int argc, char** argv) +{ + Environment env(argc, argv); + + testDim<2>(); + testDim<3>(); + + return report_errors(); +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 046593ecb00fe01dff83394043b62ff0c251510b..ddbc3c4c80cae1ff2099720d9745595327f83053 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -6,6 +6,9 @@ dune_add_test(SOURCES AccessTest.cpp dune_add_test(SOURCES AdaptInfoTest.cpp LINK_LIBRARIES amdis) +dune_add_test(SOURCES AdaptiveGridTest.cpp + LINK_LIBRARIES amdis) + foreach(_GRID RANGE 7) dune_add_test(NAME "BackupRestoreTest_${_GRID}" SOURCES BackupRestoreTest.cpp diff --git a/test/DOFVectorTest.cpp b/test/DOFVectorTest.cpp index 7a448f56947f2da4182c3c53c4d205d0df3cbb40..e673e38b9423aa7956b64af26f6d76aa658ffd8e 100644 --- a/test/DOFVectorTest.cpp +++ b/test/DOFVectorTest.cpp @@ -60,8 +60,8 @@ int main(int argc, char** argv) Dune::FieldVector<double, 2> L; L = 1.0; auto s = Dune::filledArray<2>(1); HostGrid hostGrid(L, s); - auto grid = AdaptiveGrid<HostGrid>::instance(hostGrid); - auto const& gridView = grid->leafGridView(); + AdaptiveGrid<HostGrid> grid(hostGrid); + auto const& gridView = grid.leafGridView(); // create basis auto preBasis = composite(power<2>(lagrange<2>(), flatInterleaved()), @@ -99,10 +99,10 @@ int main(int argc, char** argv) test_dofvector(vec1); for (auto const& e : elements(gridView)) - grid->mark(1, e); - grid->preAdapt(); - grid->adapt(); - grid->postAdapt(); + 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); diff --git a/test/GradientTest.cpp b/test/GradientTest.cpp index 7d6d1ef3f2deff5a8c3a6428e1971e733ec1dcba..85023b41db963f41436db89d44048ac188df3ab2 100644 --- a/test/GradientTest.cpp +++ b/test/GradientTest.cpp @@ -21,8 +21,8 @@ void test() 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(); + AdaptiveGrid<HostGrid> grid(hostGrid); + auto const& gridView = grid.leafGridView(); auto basis = Basis1::create(gridView); diff --git a/test/ISTLCommTest.cpp b/test/ISTLCommTest.cpp index c2cdcb858b29d2a5ecc56c31a05730c546070d70..a8b5791202f91a8673e1c2ac0615e87a6ac17672 100644 --- a/test/ISTLCommTest.cpp +++ b/test/ISTLCommTest.cpp @@ -76,22 +76,22 @@ int main(int argc, char** argv) { YaspGrid2d hostGrid(lower2d, upper2d, nElements2d, std::bitset<2>(), ovlp); - auto grid = AdaptiveGrid<YaspGrid2d>::instance(hostGrid); + AdaptiveGrid<YaspGrid2d> grid(hostGrid); - auto l1 = LagrangeBasis<YaspGrid2d, 1>::create(grid->leafGridView()); + auto l1 = LagrangeBasis<YaspGrid2d, 1>::create(grid.leafGridView()); AMDIS_TEST(test(l1, "Yasp2d_L1_Ovlp")); - auto th = TaylorHoodBasis<YaspGrid2d, 1>::create(grid->leafGridView()); + auto th = TaylorHoodBasis<YaspGrid2d, 1>::create(grid.leafGridView()); AMDIS_TEST(test(th, "Yasp2d_TH_Ovlp")); } { YaspGrid3d hostGrid(lower3d, upper3d, nElements3d, std::bitset<3>(), ovlp); - auto grid = AdaptiveGrid<YaspGrid3d>::instance(hostGrid); + AdaptiveGrid<YaspGrid3d> grid(hostGrid); - auto l1 = LagrangeBasis<YaspGrid3d, 1>::create(grid->leafGridView()); + auto l1 = LagrangeBasis<YaspGrid3d, 1>::create(grid.leafGridView()); AMDIS_TEST(test(l1, "Yasp3d_L1_Ovlp")); - auto th = TaylorHoodBasis<YaspGrid3d, 1>::create(grid->leafGridView()); + 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 bc931654c48512dd8e9260c805342a4b15def61c..e2fa701ccea151ddc3c0fadf221bf35ff7e5b0a0 100644 --- a/test/MarkerTest.cpp +++ b/test/MarkerTest.cpp @@ -52,6 +52,7 @@ int main(int argc, char** argv) break; } + AMDIS_TEST_EQ(prob.solutionVector()->localSize(), prob.globalBasis()->dimension()); prob.solution().interpolate(markerFunc); prob.writeFiles(adaptInfo); diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp index 9d44993f3fc5d15f6383e224c6ea8d6a3cc1934f..4d99b534071c29c8699a7ff3fe13c412dd578564 100644 --- a/test/ProblemStatTest.cpp +++ b/test/ProblemStatTest.cpp @@ -17,14 +17,14 @@ void test() // use T as coordinate type using HostGrid = Dune::YaspGrid<2, Dune::EquidistantCoordinates<T,2>>; HostGrid hostGrid({T(1), T(1)}, {8,8}); - auto grid = AdaptiveGrid<HostGrid>::instance(hostGrid); + AdaptiveGrid<HostGrid> grid(hostGrid); // use T as range type for basis using namespace Dune::Functions::BasisFactory; #if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) - auto basis = makeGlobalBasis(grid->leafGridView(), power<1>(lagrange<1>(), flatLexicographic())); + auto basis = makeGlobalBasis(grid.leafGridView(), power<1>(lagrange<1>(), flatLexicographic())); #else - auto basis = makeGlobalBasis(grid->leafGridView(), power<1>(lagrange<1,T>(), flatLexicographic())); + auto basis = makeGlobalBasis(grid.leafGridView(), power<1>(lagrange<1,T>(), flatLexicographic())); #endif using Basis = decltype(basis);