foamgrid.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
#ifndef DUNE_FOAMGRID_HH
#define DUNE_FOAMGRID_HH
 
#include <array>
#include <list>
#include <set>
#include <map>
#include <tuple>
#include <utility>
#include <type_traits>
#include <functional>
 
#include <dune/common/parallel/communication.hh>
 
#include <dune/common/stdstreams.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/grid.hh>
 
// Implementation classes
#include "foamgrid/foamgridvertex.hh"
#include "foamgrid/foamgridedge.hh"
#include "foamgrid/foamgridelements.hh"
 
// The components of the FoamGrid interface
#include "foamgrid/foamgridgeometry.hh"
#include "foamgrid/foamgridentity.hh"
#include "foamgrid/foamgridentityseed.hh"
#include "foamgrid/foamgridintersectioniterators.hh"
#include "foamgrid/foamgridleveliterator.hh"
#include "foamgrid/foamgridleafiterator.hh"
#include "foamgrid/foamgridhierarchiciterator.hh"
#include "foamgrid/foamgridindexsets.hh"
 
namespace Dune {
 
// Forward declaration
template <int dimgrid, int dimworld, class ctype = double>
class FoamGrid;
 
 
template<int dimgrid, int dimworld, class ctype>
struct FoamGridFamily
{
    typedef GridTraits<
        dimgrid,   // dimgrid
        dimworld,   // dimworld
        Dune::FoamGrid<dimgrid, dimworld, ctype>,
        FoamGridGeometry,
        FoamGridEntity,
        FoamGridLevelIterator,
        FoamGridLeafIntersection,
        FoamGridLevelIntersection,
        FoamGridLeafIntersectionIterator,
        FoamGridLevelIntersectionIterator,
        FoamGridHierarchicIterator,
        FoamGridLeafIterator,
        FoamGridLevelIndexSet< const FoamGrid<dimgrid, dimworld, ctype> >,
        FoamGridLeafIndexSet< const FoamGrid<dimgrid, dimworld, ctype> >,
        FoamGridIdSet< const FoamGrid<dimgrid, dimworld, ctype> >,  // global IdSet
        unsigned int,   // global id type
        FoamGridIdSet< const FoamGrid<dimgrid, dimworld, ctype> >,  // local IdSet
        unsigned int,   // local id type
        Communication<Dune::No_Comm> ,
        DefaultLevelGridViewTraits,
        DefaultLeafGridViewTraits,
        FoamGridEntitySeed
            > Traits;
};
 
 
 
template <int dimgrid, int dimworld, class ct>
class FoamGrid :
        public GridDefaultImplementation  <dimgrid, dimworld, ct, FoamGridFamily<dimgrid, dimworld, ct> >
{
 
    friend class FoamGridLevelIndexSet<const FoamGrid >;
    friend class FoamGridLeafIndexSet<const FoamGrid >;
    friend class FoamGridIdSet<const FoamGrid >;
    friend class FoamGridHierarchicIterator<const FoamGrid >;
    friend class FoamGridLevelIntersectionIterator<const FoamGrid >;
    friend class FoamGridLeafIntersectionIterator<const FoamGrid >;
    friend class FoamGridLevelIntersection<const FoamGrid >;
 
    template<int codim, PartitionIteratorType pitype, class GridImp_>
    friend class FoamGridLevelIterator;
 
    template<int codim, PartitionIteratorType pitype, class GridImp_>
    friend class FoamGridLeafIterator;
 
    template <class GridType_>
    friend class GridFactory;
 
    template <int dimgrid_, int dimworld_, class ct_>
    friend class GridFactoryBase;
 
    template<int codim_, int dim_, class GridImp_>
    friend class FoamGridEntity;
 
public:
 
    static_assert(dimgrid==1 || dimgrid==2, "Use FoamGrid only for 1d and 2d grids!");
 
    //**********************************************************
    // The Interface Methods
    //**********************************************************
 
    typedef FoamGridFamily<dimgrid, dimworld, ct>  GridFamily;
 
    typedef typename FoamGridFamily<dimgrid, dimworld, ct>::Traits Traits;
 
    typedef ct ctype;
 
    FoamGrid()
    : entityImps_(makeEntityImps_())
    , levelIndexSets_(1) // we always have level 0 (even if it's empty)
    , leafIndexSet_(*this)
    , freeIdCounter_(0)
    , globalRefined(0)
    , numBoundarySegments_(0)
    , growing_(false)
    {}
 
        ~FoamGrid()
        {
            // Delete level index sets
            for (size_t i=0; i<levelIndexSets_.size(); i++)
                if (levelIndexSets_[i])
                    delete (levelIndexSets_[i]);
        }
 
 
        int maxLevel() const {
            return entityImps_.size()-1;
        }
 
 
        template<int codim>
        typename Traits::template Codim<codim>::LevelIterator lbegin (int level) const {
            if (level<0 || level>maxLevel())
                DUNE_THROW(Dune::GridError, "LevelIterator in nonexisting level " << level << " requested!");
 
            return Dune::FoamGridLevelIterator<codim,All_Partition, const Dune::FoamGrid<dimgrid, dimworld, ctype> >(std::get<dimgrid-codim>(entityImps_[level]).begin());
        }
 
 
        template<int codim>
        typename Traits::template Codim<codim>::LevelIterator lend (int level) const {
            if (level<0 || level>maxLevel())
                DUNE_THROW(GridError, "LevelIterator in nonexisting level " << level << " requested!");
 
            return Dune::FoamGridLevelIterator<codim,All_Partition, const Dune::FoamGrid<dimgrid, dimworld, ctype> >(std::get<dimgrid-codim>(entityImps_[level]).end());
        }
 
 
        template<int codim, PartitionIteratorType PiType>
        typename Traits::template Codim<codim>::template Partition<PiType>::LevelIterator lbegin (int level) const {
            if (level<0 || level>maxLevel())
                DUNE_THROW(Dune::GridError, "LevelIterator in nonexisting level " << level << " requested!");
 
            return Dune::FoamGridLevelIterator<codim,PiType, const Dune::FoamGrid<dimgrid, dimworld, ctype> >(std::get<dimgrid-codim>(entityImps_[level]).begin());
        }
 
 
        template<int codim, PartitionIteratorType PiType>
        typename Traits::template Codim<codim>::template Partition<PiType>::LevelIterator lend (int level) const {
            if (level<0 || level>maxLevel())
                DUNE_THROW(GridError, "LevelIterator in nonexisting level " << level << " requested!");
 
            return Dune::FoamGridLevelIterator<codim,PiType, const Dune::FoamGrid<dimgrid, dimworld, ctype> >(std::get<dimgrid-codim>(entityImps_[level]).end());
        }
 
 
        template<int codim>
        typename Traits::template Codim<codim>::LeafIterator leafbegin() const {
            return FoamGridLeafIterator<codim,All_Partition, const FoamGrid >(*this);
        }
 
 
        template<int codim>
        typename Traits::template Codim<codim>::LeafIterator leafend() const {
            return FoamGridLeafIterator<codim,All_Partition, const FoamGrid >();
        }
 
 
        template<int codim, PartitionIteratorType PiType>
        typename Traits::template Codim<codim>::template Partition<PiType>::LeafIterator leafbegin() const {
            return FoamGridLeafIterator<codim,PiType, const FoamGrid >(*this);
        }
 
 
        template<int codim, PartitionIteratorType PiType>
        typename Traits::template Codim<codim>::template Partition<PiType>::LeafIterator leafend() const {
            return FoamGridLeafIterator<codim,PiType, const FoamGrid >();
        }
 
 
        int size (int level, int codim) const {
 
            // Turn dynamic index into static index
            if ((codim==2 && dimgrid==2) || (codim==1 && dimgrid==1))
                return std::get<0>(entityImps_[level]).size();
            if ((codim==1 && dimgrid==2))
                return std::get<1>(entityImps_[level]).size();
            if (codim==0)
                return std::get<dimgrid>(entityImps_[level]).size();
 
            return 0;
        }
 
 
        int size (int codim) const{
            return leafIndexSet().size(codim);
        }
 
 
        int size (int level, GeometryType type) const {
            return this->levelIndexSet(level).size(type);
        }
 
 
        int size (GeometryType type) const
        {
            return this->leafIndexSet().size(type);
        }
 
        size_t numBoundarySegments() const
        {
            return numBoundarySegments_;
        }
 
        const typename Traits::GlobalIdSet& globalIdSet() const{
            return idSet_;
        }
 
 
        const typename Traits::LocalIdSet& localIdSet() const{
            return idSet_;
        }
 
 
        const typename Traits::LevelIndexSet& levelIndexSet(int level) const
        {
          if (level<0 || level>maxLevel())
            DUNE_THROW(GridError, "LevelIndexSet for nonexisting level " << level << " requested!");
 
          if (! levelIndexSets_[level])
          {
            levelIndexSets_[level] = new FoamGridLevelIndexSet<const FoamGrid>(*this, level);
            levelIndexSets_[level]->update();
          }
          return *levelIndexSets_[level];
        }
 
 
        const typename Traits::LeafIndexSet& leafIndexSet() const
        {
            return leafIndexSet_;
        }
 
        template <class EntitySeed>
        typename Traits::template Codim<EntitySeed::codimension>::Entity
        entity(const EntitySeed& seed) const
        {
          const int codim = EntitySeed::codimension;
          return FoamGridEntity<codim, dimgrid, const FoamGrid>(seed.impl().target());
        }
 
 
 
 
        void globalRefine (int refCount = 1);
 
        bool mark(int refCount, const typename Traits::template Codim<0>::Entity & e)
        {
            if (!e.isLeaf())
                return false;
 
            if (refCount>=1)
                const_cast<FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>*>(e.impl().target_)->markState_ = FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>::REFINE;
            else if (refCount<0)
                const_cast<FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>*>(e.impl().target_)->markState_ = FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>::COARSEN;
            else
                const_cast<FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>*>(e.impl().target_)->markState_ = FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>::DO_NOTHING;
 
            return true;
        }
 
        int getMark(const typename Traits::template Codim<0>::Entity & e) const
        {
            switch(e.impl().target_->markState_)
            {
              case FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>::DO_NOTHING:
              case FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>::IS_COARSENED:
                return 0;
              case FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>::REFINE:
                return 1;
              case FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>::COARSEN:
                return -1;
            }
            return 0;
        }
 
        bool preAdapt();
 
        bool adapt();
 
        void postAdapt();
 
        void setPosition(const typename Traits::template Codim<dimgrid>::Entity & e,
                         const FieldVector<ctype, dimworld>& pos);
 
 
        unsigned int insertVertex(const FieldVector<ctype,dimworld>& pos)
        {
          if(!growing_) initializeGrowth_();
 
          // the final level of the vertex will be the minimum common vertex level of the new element
          verticesToInsert_.push_back(FoamGridEntityImp<0, dimgrid, dimworld, ctype>(0, pos, -verticesToInsert_.size()-1)); // initialize with some invalid id
          FoamGridEntityImp<0, dimgrid, dimworld, ctype>& newVertex = verticesToInsert_.back();
          newVertex.isNew_ = true;
          // new vertices are numbered consecutively starting from
          // the highest available index in the leaf index set +1
          newVertex.growthInsertionIndex_ = this->leafGridView().size(dimgrid) - 1 + verticesToInsert_.size();
          return newVertex.growthInsertionIndex_;
        }
 
        unsigned int insertElement(const GeometryType& type,
                                   const std::vector<unsigned int>& vertices)
        {
          // foamgrid only supports simplices until now
          assert(type.isTriangle() || type.isLine());
          assert(vertices.size() == dimgrid + 1);
 
          // the final level of the element will be the minimum common vertex level
          elementsToInsert_.push_back(FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>(0, -elementsToInsert_.size()-1)); // initialize with some invalid id
          FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>& newElement = elementsToInsert_.back();
          assert(vertices.size() == newElement.vertex_.size());
 
          for(std::size_t i = 0; i < vertices.size(); i++)
          {
            if(int(vertices[i]) >= this->leafGridView().size(dimgrid))
            {
              // initialize with pointer to vertex in verticesToInsert_ vector, later overwrite with actual pointer
              auto vIt = verticesToInsert_.begin();
              std::advance(vIt, vertices[i] - this->leafGridView().size(dimgrid));
              newElement.vertex_[i] = &*vIt;
            }
            else
            {
              // make sure the index to vertex map has been initialized
              if(!growing_) initializeGrowth_();
              // the vertex already exists in the grid, initialize with leaf vertex, later overwrite with lowest level father
              assert(indexToVertexMap_[vertices[i]]->isLeaf());
              newElement.vertex_[i] = indexToVertexMap_[vertices[i]];
            }
          }
          newElement.isNew_ = true;
          newElement.growthInsertionIndex_ = elementsToInsert_.size()-1;
          return newElement.growthInsertionIndex_;
        }
 
        unsigned int insertElement(const GeometryType& type,
                                   const std::vector<unsigned int>& vertices,
                                   std::function<FieldVector<ctype,dimworld>(FieldVector<ctype,dimgrid>)> elementParametrization)
        {
          auto growthInsertionIndex = insertElement(type, vertices);
          // save the pointer to the element parametrization
          elementsToInsert_.back().elementParametrization_ = elementParametrization;
          return growthInsertionIndex;
        }
 
        void removeElement(const typename Traits::template Codim<0>::Entity & e)
        {
          // save entity for later, actual removal happens in grow()
          elementsToRemove_.push_back(const_cast<FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>*> (e.impl().target_));
        }
 
        bool preGrow();
 
        bool grow();
 
        void postGrow();
 
        unsigned int growthInsertionIndex(const typename Traits::template Codim<0>::Entity & e) const
        {
            int idx = e.impl().target_->growthInsertionIndex_;
            assert(idx >= 0);
            return static_cast<unsigned int>(idx);
        }
 
        unsigned int growthInsertionIndex(const typename Traits::template Codim<dimgrid>::Entity & e) const
        {
            int idx = e.impl().target_->growthInsertionIndex_;
            assert(idx >= 0);
            return static_cast<unsigned int>(idx);
        }
 
 
        unsigned int overlapSize(int codim) const {
            return 0;
        }
 
 
        unsigned int ghostSize(int codim) const {
            return 0;
        }
 
 
        unsigned int overlapSize(int level, int codim) const {
            return 0;
        }
 
 
        unsigned int ghostSize(int level, int codim) const {
            return 0;
        }
 
        template<class DataHandle>
        bool loadBalance(DataHandle& data)
        {
            return loadBalance();
        }
 
        bool loadBalance()
        {
            if (comm().size() > 1)
                DUNE_THROW(Dune::NotImplemented, "Load balancing not implemented. Foamgrid does not run in parallel yet!");
            return false;
        }
 
        template<class T, template<class> class P, int codim>
        void communicate (T& t, InterfaceType iftype, CommunicationDirection dir, int level) const
        {}
 
        template<class DataHandle>
        void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
        {}
 
        template<class DataHandle>
        void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const
        {}
 
        const typename Traits::Communication& comm () const
        {
            return ccobj_;
        }
        // **********************************************************
        // End of Interface Methods
        // **********************************************************
 
    private:
 
    void initializeGrowth_()
    {
      // update the index to vertex map
      indexToVertexMap_.resize(this->leafGridView().size(dimgrid));
      for (const auto& vertex : vertices(this->leafGridView()))
      {
        std::size_t index = leafIndexSet().index(vertex);
        indexToVertexMap_[index] = const_cast<FoamGridEntityImp<0, dimgrid ,dimworld, ctype>*>(vertex.impl().target_);
      }
 
      // tell the grid it's ready for growth
      growing_ = true;
    }
 
    void erasePointersToEntities(std::list<FoamGridEntityImp<dimgrid, dimgrid ,dimworld, ctype> >& elements);
 
    template<int i>
    void eraseVanishedEntities(std::list<FoamGridEntityImp<i, dimgrid, dimworld, ctype> >& levelEntities);
 
    void coarsenSimplexElement(FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>& element);
 
    void refineSimplexElement(FoamGridEntityImp<2, 2, dimworld, ctype>& element, int refCount);
    void refineSimplexElement(FoamGridEntityImp<1, 1, dimworld, ctype>& element, int refCount);
 
    bool removeSimplexElement(FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>& element);
 
    void overwriteFineLevelNeighbours(FoamGridEntityImp<dimgrid-1, dimgrid, dimworld, ctype>& edge,
                                      const FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>* son,
                                      const FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>* father);
 
    void addElementForFacet(const FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>* element,
                            FoamGridEntityImp<dimgrid-1, dimgrid, dimworld, ctype>* facet);
 
    void addNewFacet(FoamGridEntityImp<dimgrid-1, dimgrid, dimworld, ctype>* &facet,
                     std::array<FoamGridEntityImp<0, dimgrid, dimworld, ctype>*,dimgrid> vertexArray,
                     int level);
 
    void setIndices();
 
    unsigned int getNextFreeId()
    {
      return freeIdCounter_++;
    }
 
    void computeTwoZeroConnectivity()
    {
      // the elements vector for vertices is already available as they are facets in 1d grids
      for(int level = 0; level <= maxLevel(); ++level)
      {
        // do it everytime freshly so that we don't have to take care of 2-0 connectivity in adaptivity
        for(auto&& vertex : std::get<0>(entityImps_[level]))
          vertex.elements_.clear();
 
        for(auto eIt = std::get<dimgrid>(entityImps_[level]).begin(); eIt != std::get<dimgrid>(entityImps_[level]).end(); ++eIt)
          for(auto&& vertex : eIt->vertex_)
            vertex->elements_.push_back(&*eIt);
      }
    }
 
    typename Traits::Communication ccobj_;
 
    template<std::size_t... dimEntity>
    static auto makeEntityImpsImpl_(std::index_sequence<dimEntity...>, std::size_t numLevels)
    { return std::vector<std::tuple<std::list<FoamGridEntityImp<dimEntity, dimgrid, dimworld, ctype>>...>>(numLevels); }
 
    // Create the lists of vertices, edges, elements for each level
    static auto makeEntityImps_(std::size_t numLevels = 1)
    { return makeEntityImpsImpl_(std::make_index_sequence<dimgrid+1>{}, numLevels); }
 
    // The vector type for tuple of lists of vertices, edges, elements for each level
    using EntityImps = std::decay_t<decltype(makeEntityImps_())>;
 
    // The tuple type of lists of vertices, edges, elements
    using EntityTuple = typename EntityImps::value_type;
 
    EntityImps entityImps_;
 
    mutable std::vector<FoamGridLevelIndexSet<const FoamGrid>*> levelIndexSets_;
 
    FoamGridLeafIndexSet<const FoamGrid > leafIndexSet_;
 
    FoamGridIdSet<const FoamGrid > idSet_;
 
    unsigned int freeIdCounter_;
 
    int globalRefined;
 
    std::size_t numBoundarySegments_;
 
    // True if the last call to preadapt returned true
    bool willCoarsen;
 
    std::vector<FoamGridEntityImp<0, dimgrid, dimworld, ctype>* > indexToVertexMap_;
 
    std::vector<FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>* > elementsToRemove_;
 
    std::list<FoamGridEntityImp<0, dimgrid, dimworld, ctype> > verticesToInsert_;
 
    std::list<FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype> > elementsToInsert_;
 
    bool growing_;
 
}; // end Class FoamGrid
 
#include <dune/foamgrid/foamgrid/foamgrid.cc>
 
namespace Capabilities
{
    template<int dimgrid, int dimworld, class ctype, int codim>
    struct hasEntity< FoamGrid<dimgrid, dimworld, ctype>, codim>
    {
        static const bool v = true;
    };
 
    template <int dimgrid, int dimworld, class ctype>
    struct isLevelwiseConforming< FoamGrid<dimgrid, dimworld, ctype> >
    {
        static const bool v = false;
    };
 
    template <int dimgrid, int dimworld, class ctype>
    struct isLeafwiseConforming< FoamGrid<dimgrid, dimworld, ctype> >
    {
        static const bool v = false;
    };
 
    template<int dimgrid, int dimworld, class ctype>
    struct viewThreadSafe< FoamGrid<dimgrid, dimworld, ctype> > {
      static const bool v = true;
    };
}
 
} // namespace Dune
 
 
// The factory should be present whenever the user includes foamgrid.hh.
// However since the factory needs to know the grid the include directive
// comes here at the end.
#include <dune/foamgrid/foamgrid/foamgridfactory.hh>
 
#endif