gridmanager.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2023 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#ifndef GRIDMANAGER_HH_
#define GRIDMANAGER_HH_
 
#include <ctime>
#include <iostream>
#include <mutex>
#include <numeric>
#include <utility>
 
#include <boost/signals2.hpp>
#include <boost/range/iterator_range.hpp>
 
#include <dune/common/config.h>
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/grid/io/file/dgfparser/gridptr.hh>
 
#include "fem/gridBasics.hh"
#include "utilities/threading.hh"
#include "utilities/timing.hh"
 
 
 
namespace Kaskade
{
 
  // forward declaration
  template<class G, class T> class CellData;
 
  // ----------------------------------------------------------------------------------------------

  struct GridSignals
  {
    enum { freeResources=0, allocResources=1 };

    enum Status { BeforeRefinement, BeforeRefinementComplete, AfterRefinement, TransferCompleted };
 

    boost::signals2::signal<void (Status)> informAboutRefinement;
  };
  
  // ----------------------------------------------------------------------------------------------

  template <class GridView>
  class CellRanges
  {
  public:
    typedef typename GridView::template Codim<0>::Iterator CellIterator;
    
    CellRanges() = default;
    
    CellRanges(CellRanges&&) = default;
    
    CellRanges(int nmax, GridView const& gridView): cellRanges(nmax+1)
    {
      // We know how many iterators we will insert - allocate space beforehand. 
      for (int n=1; n<=nmax; ++n)
        cellRanges[n].reserve(n+1);
      
      // number of cells
      size_t const size = gridView.size(0);
      
      // Step through all cells, storing the iterator in case it is at a boundary of cell ranges.
      // The iteration through the cells is monotone, hence the ranges we construct cannot overlap.
      // If few cells exist, some ranges may be empty, however.
      // The following postconditions hold: 
      // (i) cellRanges[n].size() <= n. 
      // Proof: Assume cellRanges[n].size()==n. Then there is no cell with j >= size
      //        (it is at most size-1), where a further push_back could happen.
      // (ii) cellRanges[n][0] == begin. 
      // Proof: In the first iterations j==0 and cellRanges[n].size()==0 holds,
      //        such that 0>=0 evaluates true and triggers the push_back(begin).
      // (iii) any range contains less than size/n + 1 cells.
      // Proof: Assume a range [r,s[. Then there is k with r>=k*size/n and t<(k+1)*size/n
      //        for all t<s. Thus s-1-r < ((k+1)-k)*size/n = size/n,
      //        such that the length is s-r < size/n + 1.
      size_t j = 0; // cell counter
      for (CellIterator i=gridView.template begin<0>(); i!=gridView.template end<0>(); ++i, ++j)
        for (int n=1; n<=nmax; ++n)
          if (j >= cellRanges[n].size()*size/n)
            cellRanges[n].push_back(i);
 
      // Now we have to make sure that every range vector contains the number of ranges
      // it's supposed to hold. Add empty ranges as needed.
      // As up to here the size of cellRanges[n] is at most n (not n+1), we automatically take care of
      // terminating the last range with the end iterator.
      for (int n=1; n<=nmax; ++n)
        while (cellRanges[n].size()<n+1)
          cellRanges[n].push_back(gridView.template end<0>());
    }
    
    CellRanges<GridView>& operator=(CellRanges<GridView>&&) = default;
    
    boost::iterator_range<CellIterator> range(int n, int k) const
    {
      assert(1<=n && n<cellRanges.size());
      assert(0<=k && k<n);
      return boost::iterator_range<CellIterator>(cellRanges[n][k],cellRanges[n][k+1]);
    }
    
    int maxRanges() const
    {
      return static_cast<int>(cellRanges.size())-1;
    }
    
  private:
    std::vector<std::vector<CellIterator>> cellRanges;
  };
  
  //----------------------------------------------------------------------------------------
 
  // Start of implementation for better numbering of entities. Very much in progress, still.
  template <class GridView>
  class EntityNumbering
  {
    using IndexSet = typename GridView::IndexSet;
    constexpr static int dim = GridView::dimension;
    
  public:
    using Index = typename IndexSet::IndexType;
    
    enum SortOrder { 
      SORTBYCODIM, 
      SEQUENTIAL, 
      CUTHILLMCKEE };
 
    EntityNumbering(GridView const& gv, SortOrder order)
    {
      IndexSet const& is = gv.indexSet();
      for (int codim=0; codim<=dim; ++codim)
        indices[codim].resize(is.size(codim));
      
      switch (order)
      {
        case SORTBYCODIM:
        {
          Index n = 0;
          for (int codim=0; codim<=dim; ++codim)
          {
            std::iota(begin(indices[codim]),end(indices[codim]),n);
            n += indices[codim].size();
          }
          break;
        }
        case SEQUENTIAL:
          abort();
          break;
        case CUTHILLMCKEE:
          abort();
          break;
      }
    }
    
    Index operator()(int codim, Index i) const 
    {
      assert(0<=codim && codim<=dim);
      assert(0<=i && i<indices[codim].size());
      return indices[codim][i];
    }
    
  private:
    std::array<std::vector<Index>,dim+1> indices;
  };
  
  //----------------------------------------------------------------------------------------
  
  template<class Grd>
  class GridManagerBase
  {
  public:    
    using Grid = Grd;
    
    typedef GridManagerBase<Grid> Self;
    
    GridManagerBase(Dune::GridPtr<Grid> grid_, bool verbose_, bool enforceConcurrentReads = false)
    : GridManagerBase(grid_.release(),verbose_,enforceConcurrentReads)
    {}
 
    GridManagerBase(std::unique_ptr<Grid>&& grid_, bool verbose_, bool enforceConcurrentReads = false)
    : GridManagerBase(grid_.release(),verbose_,enforceConcurrentReads)
    {}
 
    GridManagerBase(Grid*&& grid_, bool verbose_, bool enforceConcurrentReads = false):
      markCount(0), gridptr(grid_),
      gridIsThreadSafe_(Dune::Capabilities::viewThreadSafe<Grid>::v), enforceConcurrentReads_(enforceConcurrentReads),
      levelRanges(gridptr->maxLevel()+1), verbose(verbose_)
    {}
 
    virtual ~GridManagerBase() {};
    

    Grid const& grid() const
    {
      return *gridptr;
    }
    
    std::shared_ptr<Grid const> gridShared() const
    {
      return gridptr;
    }

    
    
    bool mark(int refCount, typename Grid::template Codim<0>::Entity const& cell) 
    {
      if(refCount != 0)
      {
        ++markCount;
        return gridptr->mark(refCount,cell);
      }
      return false;
    }

    template<class S>
    void mark(CellData<Grid,S> const& cellData)
    {
      this->flushMarks();
      typename CellData<Grid,S>::CellDataVector::const_iterator imark;
      typename CellData<Grid,S>::CellDataVector::const_iterator imark_end = cellData.end();
      for(imark=cellData.begin(); imark != imark_end; ++imark)
      {
        auto ret = this->mark(imark->first,imark->second);
        ++markCount;
        if (!ret && imark->first!=0)
          std::cout << "GRIDMANAGER: Warning: Cannot mark element!" << '\n';
      }
    }

    size_t countMarked() const
    {
      return markCount;
    }

    bool preAdapt()
    {
      if (markCount > 0)
        return gridptr->preAdapt(); 
      else 
        return false; 
    };

    void postAdapt()
    {
      gridptr->postAdapt();
 
      if(verbose)
        std::cout << "after refinement: " << gridptr->size(0) << " cells, "
                  << gridptr->size(Grid::dimension) << " vertices.\n";
    };
 

    bool adaptAtOnce()
    {
      this->preAdapt();
      bool result = adapt();
      #if !defined(NDEBUG)
      if(!result)
        std::cout << "GRIDMANAGER: Nothing has been refined!" << '\n';
      #endif
      this->postAdapt();
      return result;
    }
    
    void flushMarks()
    {
      for (auto const& element: elements(gridptr->leafGridView()))
        mark(0,element);
 
      markCount = 0;
    }

    bool adapt()
    {
      bool result;
      if (this->markCount > 0)
      {
        Timings& timer = Timings::instance();
        ScopedTimingSection refineSection("grid refinement");
        
        timer.start("before refinement");
        this->beforeRefinement();
        timer.stop("before refinement");
        
        if (verbose)
          std::cout << "\nGRIDMANAGER: mesh is refined from " << this->gridptr->size(0) << " cells " << std::flush;
        
        timer.start("adapt");
        result = adaptGrid();
        update();
        timer.stop("adapt");
        
        if(verbose) 
          std::cout << "to " << this->gridptr->size(0) << " cells" << '\n';
        
        timer.start("after refinement");
        this->afterRefinement();
        timer.stop("after refinement");
        
        this->markCount = 0;
      } 
      else
        result = false;
      return result;
    }
    
    void globalRefine(int refCount) 
    {
      if (refCount<=0)
        return;
      
      Timings& timer = Timings::instance();
      ScopedTimingSection refineSection("grid refinement");
      timer.start("before refinement");
      this->beforeRefinement();
      timer.stop("before refinement");
      
      timer.start("adapt");
      refineGrid(refCount);
      update();
      timer.stop("adapt");
      
      timer.start("after refinement");
      this->afterRefinement();
      timer.stop("after refinement");
      
      this->markCount = 0;
    };
    
    virtual void update() = 0;
    
    
    Self& enforceConcurrentReads(bool enforceConcurrentReads) 
    { 
      enforceConcurrentReads_ = enforceConcurrentReads; 
      return *this; 
    }
 

    Self& setVerbosity(const bool verbosity) 
    { 
      verbose = verbosity; return *this; 
    }

    bool gridIsThreadSafe() const
    {
#ifdef NDEBUG
      return gridIsThreadSafe_ || enforceConcurrentReads_;
#else
      return false;
#endif
    }
 

    CellRanges<typename Grid::LeafGridView> const& cellRanges(typename Grid::LeafGridView const&) const
    {
      std::lock_guard<std::mutex> lock(cellRangeMutex);
      if (!leafRanges)
        leafRanges.reset(new CellRanges<typename Grid::LeafGridView>(2*NumaThreadPool::instance().cpus(),grid().leafGridView()));
      
      return *leafRanges;
    }
    
    CellRanges<typename Grid::LevelGridView> const& cellRanges(typename Grid::LevelGridView const& gridView) const
    {
      int const level = gridView.template begin<0>()->level(); // cells have the level of their level view (inferred from the Dune tutorial)
      assert(level<=grid().maxLevel());
      
      std::lock_guard<std::mutex> lock(cellRangeMutex);
      if (levelRanges[level].maxRanges()<0)
        levelRanges[level] = CellRanges<typename Grid::LevelGridView>(2*NumaThreadPool::instance().cpus(),grid().levelGridView(level));
        // we could have used gridView directly here. But can we make sure this is a level view of the grid we keep?
      return levelRanges[level];
    }
    
    EntityNumbering<typename Grid::LeafGridView> const& entityNumbering(typename Grid::LeafGridView const&) const
    {
      using Numbering = EntityNumbering<typename Grid::LeafGridView>;
      
      std::lock_guard<std::mutex> lock(cellRangeMutex); // maybe we should use our own mutex here?
      if (!leafEntityNumbers)
        leafEntityNumbers = std::make_unique<Numbering>(grid().leafGridView(),Numbering::SORTBYCODIM);
      
      return *leafEntityNumbers;
    }
    
    
  private:
    // TODO: docme!
    virtual void refineGrid(int refcount) = 0;
    virtual bool adaptGrid() = 0;
    
    size_t markCount; // number of cells marked for refinement/coarsening
 
  protected:
    void beforeRefinement()
    {
      // trigger creation of transfer data structures
      signals.informAboutRefinement(GridSignals::BeforeRefinement);
      signals.informAboutRefinement(GridSignals::BeforeRefinementComplete);
 
      // invalidate the cell range data structures
      leafRanges.reset();
      levelRanges.clear();
      leafEntityNumbers.reset();
    }
 
    void afterRefinement()
    {
      levelRanges.resize(grid().maxLevel()+1);
      signals.informAboutRefinement(GridSignals::AfterRefinement);
      signals.informAboutRefinement(GridSignals::TransferCompleted);
    }
 
    template<class LocalToGlobalMapper> friend class FEFunctionSpace;
    template<class GridMan, class ErrEst> friend class AdaptiveGrid;

    std::shared_ptr<Grid> gridptr;
    
  public:
    GridSignals signals;
 
  private:
    const bool gridIsThreadSafe_;
    
    bool enforceConcurrentReads_;
    
    // cell ranges and entity numberings are constructed on demand
    mutable std::mutex cellRangeMutex;
    mutable std::unique_ptr<CellRanges<typename Grid::LeafGridView>> leafRanges;
    mutable std::vector<CellRanges<typename Grid::LevelGridView>> levelRanges;
    mutable std::unique_ptr<EntityNumbering<typename Grid::LeafGridView>> leafEntityNumbers;
 
  protected:
    bool verbose;
  };
 
  // -------------------------------------------------------------------
  
  template <class Grd>
  class GridManager : public GridManagerBase<Grd>
  {
  public:
    typedef Grd Grid;

    explicit GridManager(Dune::GridPtr<Grid> grid_, bool verbose=false) 
    : GridManagerBase<Grid>(grid_,verbose)
    {}


    explicit GridManager(std::unique_ptr<Grid>&& grid_, bool verbose=false) 
    : GridManagerBase<Grid>( std::move(grid_),verbose)
    {}


    explicit GridManager(Grid*&& grid_, bool verbose=false) 
    : GridManagerBase<Grid>( std::move(grid_),verbose)
    {}
 
 
    virtual void update() {}
 
  private:
    virtual void refineGrid(int refCount)
    {
      this->gridptr->globalRefine(refCount);
    };
 
 
    virtual bool adaptGrid()
    {
      return this->gridptr->adapt();
    };
  };
 
 
  //---------------------------------------------------------------------

  template <class Grid>
  void refineAtBoundary(GridManagerBase<Grid>& gridManager)
  {
    for (auto const& cell: elements(gridManager.grid().leafGridView()))
      if (cell.hasBoundaryIntersections())
        gridManager.mark(1,cell);
    gridManager.adaptAtOnce();
  }
}
 
#endif /* GRIDMANAGER_HH_ */