AdaptiveGrid.hpp

#pragma once

#include <algorithm>
#include <list>
#include <memory>
#include <mutex>
#include <shared_mutex>
#include <type_traits>
#include <utility>

#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/Concepts.hpp>
#include <amdis/common/DefaultGridView.hpp>
#include <amdis/common/SharedPtr.hpp>
#include <amdis/Observer.hpp>
#include <amdis/Output.hpp>

namespace AMDiS
{
  // forward declaration
  template <class HostGrid>
  class AdaptiveGridFamily;



  template <class HG>
  class AdaptiveGrid
      : public Dune::GridDefaultImplementation<
          HG::dimension, HG::dimensionworld, typename HG::ctype, AdaptiveGridFamily<HG> >
      , public Notifier<event::preAdapt, event::adapt, event::postAdapt>
  {
    using Self = AdaptiveGrid<HG>;

  public:

    using HostGrid   = HG;
    using GridFamily = AdaptiveGridFamily<HG>;
    using Traits     = typename GridFamily::Traits;
    using Element    = typename Traits::template Codim<0>::Entity;


  public:

    template <class HostGrid_,
      Dune::disableCopyMove<Self, HostGrid_> = 0>
    explicit AdaptiveGrid(HostGrid_&& hostGrid)
      : hostGrid_(wrap_or_share(FWD(hostGrid)))
    {}

    std::shared_ptr<HostGrid> const& hostGrid() const { return hostGrid_; }


  public:


    template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
    auto lbegin(int level) const { return hostGrid_->levelGridView(level).template begin<codim,pt>(); }

    template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
    auto lend(int level) const { return hostGrid_->levelGridView(level).template end<codim,pt>(); }


    template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
    auto leafbegin() const { return hostGrid_->leafGridView().template begin<codim,pt>(); }

    template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
    auto leafend() const { return hostGrid_->leafGridView().template end<codim,pt>(); }


    auto ilevelbegin(Element const& e) const { return hostGrid_->levelGridView(e.level()).ibegin(e); }

    auto ilevelend(Element const& e) const { return hostGrid_->levelGridView(e.level()).iend(e); }

    auto ileafbegin(Element const& e) const { return hostGrid_->leafGridView().ibegin(e); }

    auto ileafend(Element const& e) const { return hostGrid_->leafGridView().iend(e); }




    int maxLevel() const { return hostGrid_->maxLevel(); }

    int size(int level, int codim) const { return hostGrid_->size(level, codim); }

    int size(int codim) const { return hostGrid_->size(codim); }

    int size(int level, Dune::GeometryType type) const { return hostGrid_->size(level, type); }

    int size(Dune::GeometryType type) const { return hostGrid_->size(type); }

    std::size_t numBoundarySegments() const { return hostGrid_->numBoundarySegments(); }




    auto const& globalIdSet() const { return hostGrid_->globalIdSet(); }

    auto const& localIdSet() const { return hostGrid_->localIdSet(); }

    auto const& levelIndexSet(int level) const { return hostGrid_->levelIndexSet(level); }

    auto const& leafIndexSet() const { return hostGrid_->leafIndexSet(); }




    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();
      }
    }

    bool mark(int refCount, Element const& e) { return hostGrid_->mark(refCount, e); }

    int getMark(Element const& e) const { return hostGrid_->getMark(e); }

    bool preAdapt()
    {
      Dune::Timer t;
      mightCoarsen_ = hostGrid_->preAdapt();
      Dune::MPIHelper::getCollectiveCommunication().max(&mightCoarsen_, 1);
      this->notify(event::preAdapt{mightCoarsen_});
      info(2,"AdaptiveGrid::preAdapt needed {} seconds", t.elapsed());
      return mightCoarsen_;
    }

    bool adapt()
    {
      Dune::Timer t;
      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_;
    }

    void postAdapt()
    {
      Dune::Timer t;
      hostGrid_->postAdapt();
      this->notify(event::postAdapt{});
      changeIndex_++;
      info(2,"AdaptiveGrid::postAdapt needed {} seconds", t.elapsed());
    }

    void update()
    {
      this->notify(event::adapt{true});
    }

    unsigned long changeIndex() const
    {
      return changeIndex_;
    }

    // @}



    auto const& comm() const { return hostGrid_->comm(); }

    template <class DataHandle>
    void communicate(DataHandle& handle, Dune::InterfaceType iftype,
                     Dune::CommunicationDirection dir, int level) const
    {
      hostGrid_->levelGridView(level).communicate(handle,iftype,dir);
    }

    template <class DataHandle>
    void communicate(DataHandle& handle, Dune::InterfaceType iftype,
                     Dune::CommunicationDirection dir) const
    {
      hostGrid_->leafGridView().communicate(handle,iftype,dir);
    }


    bool loadBalance()
    {
      if constexpr (std::is_convertible_v<decltype(std::declval<HG>().loadBalance()), bool>)
        return hostGrid_->loadBalance();
      else {
        hostGrid_->loadBalance();
        return true;
      }
    }


    template <class DataHandle>
    bool loadBalance(DataHandle& handle)
    {
      if constexpr (std::is_convertible_v<decltype(std::declval<HG>().loadBalance(handle)), bool>)
        return hostGrid_->loadBalance(handle);
      else {
        hostGrid_->loadBalance(handle);
        return true;
      }
    }


    int overlapSize(int level, int codim) const { return hostGrid_->levelGridView(level).overlapSize(codim); }

    int overlapSize(int codim) const { return hostGrid_->leafGridView().overlapSize(codim); }

    int ghostSize(int level, int codim) const { return hostGrid_->levelGridView(level).ghostSize(codim); }

    int ghostSize(int codim) const { return hostGrid_->leafGridView().ghostSize(codim); }



    template <class EntitySeed>
    auto entity(EntitySeed const& seed) const { return hostGrid_->entity(seed); }


  private:

    std::shared_ptr<HostGrid> hostGrid_;

    bool mightCoarsen_ = false;

    bool refined_ = false;

    unsigned long changeIndex_ = 0;
  };

  // deduction guide
  template <class HostGrid>
  AdaptiveGrid(HostGrid const&)
    -> AdaptiveGrid<HostGrid>;


  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> >;
    };
  };


  template <class HostGrid>
  unsigned long changeIndex(HostGrid const& /*hostGrid*/)
  {
    return 0;
  }

  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>;
    };
  }

  template <class HostGrid>
  using AdaptiveGrid_t = typename Impl::AdaptiveGridImpl<HostGrid>::type;


} // end namespace AMDiS


namespace Dune
{
  template <class HostGrid>
  class GridFactory<AMDiS::AdaptiveGrid<HostGrid> >
      : public GridFactoryInterface<AMDiS::AdaptiveGrid<HostGrid> >
  {
    using Self = GridFactory;
    using Super = GridFactoryInterface<AMDiS::AdaptiveGrid<HostGrid> >;
    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)...)
    {}

    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>()...));

    void insertElement(GeometryType const& type,
                       std::vector<unsigned int> const& vertices) override
    {
      hostFactory_.insertElement(type, vertices);
    }

#if DUNE_VERSION_LT(DUNE_GRID,2,8)
    using ElementParametrizationType = std::shared_ptr<VirtualFunction<FieldVector<ctype,dim>,FieldVector<ctype,dimworld> > >;

    void insertElement(GeometryType const& type,
                       std::vector<unsigned int> const& vertices,
                       ElementParametrizationType const& elementParametrization) override
#else
    using ElementParametrizationType
      = std::function<FieldVector<ctype,dimworld>(FieldVector<ctype,dim>)>;

    void insertElement(GeometryType const& type,
                       std::vector<unsigned int> const& vertices,
                       ElementParametrizationType elementParametrization) override
#endif
    {
      using A0 = GeometryType;
      using A1 = std::vector<unsigned int>;
      using A2 = ElementParametrizationType;
      if constexpr (Std::is_detected<HasInsertElement, HostGridFactory, A0,A1,A2>::value)
        hostFactory_.insertElement(type, vertices, elementParametrization);
      else
        AMDiS::error_exit("insertElement() not implemented for HostGrid type.");
    }
    using Super::insertElement;

    template <class F, class... Args>
    using HasInsertBoundarySegment = decltype(std::declval<F>().insertBoundarySegment(std::declval<Args>()...));

    void insertBoundarySegment(std::vector<unsigned int> const& vertices) override
    {
      hostFactory_.insertBoundarySegment(vertices);
    }

    using BoundarySegmentType = std::shared_ptr<BoundarySegment<dim,dimworld> >;

    void insertBoundarySegment(std::vector<unsigned int> const& vertices,
                               BoundarySegmentType const& boundarySegment) override
    {
      using A0 = std::vector<unsigned int>;
      using A1 = BoundarySegmentType;
      if constexpr (Std::is_detected<HasInsertBoundarySegment, HostGridFactory, A0,A1>::value)
        hostFactory_.insertBoundarySegment(vertices, boundarySegment);
      else
        AMDiS::error_exit("insertBoundarySegment() not implemented for HostGrid type.");
    }

#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:

      template <class Output>
      static void backup(Grid const& grid, Output const& filename_or_stream)
      {
        HostBackupRestoreFacility::backup(*grid.hostGrid(), filename_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


  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