mcmgmapper.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_GRID_COMMON_MCMGMAPPER_HH
#define DUNE_GRID_COMMON_MCMGMAPPER_HH
 
#include <functional>
#include <iostream>
 
#include <dune/common/deprecated.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/geometry/dimension.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/typeindex.hh>
 
#include "mapper.hh"
 
namespace Dune
{
  //
  //  Common Layout templates
  //
 
 
  template<int dimgrid> struct DUNE_DEPRECATED_MSG("The MCMG layout classes have been deprecated. Pass `mcmgElementLayout()` to the constructor instead")
  MCMGElementLayout {
    bool contains (Dune::GeometryType gt) const { return gt.dim()==dimgrid; }
  };
 
 
  template<int dim> struct DUNE_DEPRECATED_MSG("The MCMG layout classes have been deprecated. Pass `mcmgVertexLayout()` to the constructor instead")
  MCMGVertexLayout {
    bool contains (Dune::GeometryType gt) const { return gt.dim()==0; }
  };
 
  namespace Impl {
 
    /*
     * Dummy layout to be used as the default for
     * `MultipleCodimMultipleGeomTypeMapper`.  It should never be used, but
     * we need a default.
     *
     * This class can be removed once the `LayoutClass` template parameter
     * of `MultipleCodimMultipleGeomTypeMapper` is removed.
     */
    template<int dimgrid>
    struct MCMGFailLayout {
      MCMGFailLayout()
        { DUNE_THROW(Exception, "The default layout class cannot be used"); }
      bool contains(GeometryType gt) const
        { DUNE_THROW(Exception, "The default layout class cannot be used"); }
    };
 
  } /* namespace Impl */
 
  using MCMGLayout = std::function<size_t(GeometryType, int)>;
 
  template<int codim>
  MCMGLayout mcmgLayout(Codim<codim>)
  {
    return [](GeometryType gt, int dimgrid) {
      return dimgrid - gt.dim() == codim;
    };
  }
 
  template<int dim>
  MCMGLayout mcmgLayout(Dim<dim>)
  {
    return [](GeometryType gt, int) {
      return gt.dim() == dim;
    };
  }
 
  inline MCMGLayout mcmgElementLayout()
  {
    return mcmgLayout(Codim<0>());
  }
 
  inline MCMGLayout mcmgVertexLayout()
  {
    return mcmgLayout(Dim<0>());
  }
 
  //
  //  MultipleCodimMultipleGeomTypeMapper
  //
 
  template <typename GV, template<int> class LayoutClass = Impl::MCMGFailLayout>
  class MultipleCodimMultipleGeomTypeMapper :
    public Mapper<typename GV::Grid,MultipleCodimMultipleGeomTypeMapper<GV,LayoutClass>, typename GV::IndexSet::IndexType >
  {
  public:
 
    typedef GV GridView;
 
    typedef typename GV::IndexSet::IndexType Index;
 
    using size_type = decltype(std::declval<typename GV::IndexSet>().size(0));
 
    MultipleCodimMultipleGeomTypeMapper(const GV& gridView, const LayoutClass<GV::dimension> layout = {})
      DUNE_DEPRECATED_MSG("Use the constructor taking a `MCMGLayout` functional instead")
      : MultipleCodimMultipleGeomTypeMapper(gridView, wrapLayoutClass(layout))
    {}
 
    MultipleCodimMultipleGeomTypeMapper(const GV& gridView, const MCMGLayout& layout)
      : gridView_(gridView)
      , is(gridView_.indexSet())
      , layout_(layout)
    {
      update();
    }
 
    template<class EntityType>
    Index index (const EntityType& e) const
    {
      const GeometryType gt = e.type();
      assert(offset(gt) != invalidOffset);
      return is.index(e)*blockSize(gt) + offset(gt);
    }
 
    Index subIndex (const typename GV::template Codim<0>::Entity& e, int i, unsigned int codim) const
    {
      const GeometryType eType = e.type();
      GeometryType gt = eType.isNone() ?
        GeometryTypes::none( GV::dimension - codim ) :
        ReferenceElements<double,GV::dimension>::general(eType).type(i,codim) ;
      //GeometryType gt=ReferenceElements<double,GV::dimension>::general(e.type()).type(i,codim);
      assert(offset(gt) != invalidOffset);
      return is.subIndex(e, i, codim)*blockSize(gt) + offset(gt);
    }
 
    size_type size () const
    {
      return n;
    }
 
    size_type size(GeometryType gt) const
    {
      return blockSize(gt);
    }
 
    const std::vector< GeometryType >& types ( int codim ) const
    {
      return myTypes_[ codim ];
    }
 
    template<class EntityType>
    IntegralRange<Index> indices (const EntityType& e) const
    {
      if(!is.contains(e) || offset(e.type()) == invalidOffset)
        return {0,0};
      Index start = index(e);
      return {start, start+blockSize(e.type())};
    }
 
    IntegralRange<Index> indices (const typename GV::template Codim<0>::Entity& e, int i, int cc) const
    {
      const GeometryType eType = e.type();
      const GeometryType gt = eType.isNone() ?
        GeometryTypes::none(GV::dimension - cc) :
        ReferenceElements<double,GV::dimension>::general(eType).type(i,cc) ;
      if (offset(gt) == invalidOffset)
        return {0,0};
      else
      {
        Index start = subIndex(e,i,cc);
        return {start, start+blockSize(gt)};
      }
    }
 
    template<class EntityType>
    bool contains (const EntityType& e, Index& result) const
    {
      if(!is.contains(e) || offset(e.type()) == invalidOffset)
      {
        result = 0;
        return false;
      }
      result = index(e);
      return true;
    }
 
    bool contains (const typename GV::template Codim<0>::Entity& e, int i, int cc, Index& result) const
    {
      const GeometryType eType = e.type();
      const GeometryType gt = eType.isNone() ?
        GeometryTypes::none( GV::dimension - cc ) :
        ReferenceElements<double,GV::dimension>::general(eType).type(i,cc) ;
      //GeometryType gt=ReferenceElements<double,GV::dimension>::general(e.type()).type(i,cc);
      if (offset(gt) == invalidOffset)
        return false;
      result = is.subIndex(e, i, cc) + offset(gt);
      return true;
    }
 
    void update ()
    {
      n = 0;
 
      std::fill(offsets.begin(),offsets.end(),Index(0));
      std::fill(blocks.begin(),blocks.end(),Index(0));
 
      for (unsigned int codim = 0; codim <= GV::dimension; ++codim)
      {
        // walk over all geometry types in the codimension
        for (const GeometryType& gt : is.types(codim)) {
          Index offset;
          size_t block = layout()(gt, GV::Grid::dimension);
 
          // if the geometry type is contained in the layout, increment offset
          // and store geometry type
          if (block) {
            offset = n;
            n += is.size(gt) * block;
            myTypes_[codim].push_back(gt);
          }
          else {
            offset = invalidOffset;
          }
 
          offsets[GlobalGeometryTypeIndex::index(gt)] = offset;
          blocks[GlobalGeometryTypeIndex::index(gt)] = block;
        }
      }
    }
 
    const MCMGLayout &layout () const { return layout_; }
    const GridView &gridView () const { return gridView_; }
 
  private:
    Index offset(GeometryType gt) const
      { return offsets[GlobalGeometryTypeIndex::index(gt)]; }
    Index blockSize(GeometryType gt) const
      { return blocks[GlobalGeometryTypeIndex::index(gt)]; }
 
    static const Index invalidOffset = std::numeric_limits<Index>::max();
 
    // number of data elements required
    unsigned int n;
    // GridView is needed to keep the IndexSet valid
    const GV gridView_;
    const typename GV::IndexSet& is;
    // provide an array for the offsets
    std::array<Index, GlobalGeometryTypeIndex::size(GV::dimension)> offsets;
    std::array<Index, GlobalGeometryTypeIndex::size(GV::dimension)> blocks;
    const MCMGLayout layout_;     // get layout object
    std::vector<GeometryType> myTypes_[GV::dimension+1];
 
  protected:
    static MCMGLayout wrapLayoutClass(const LayoutClass<GV::dimension>& layout)
    {
      /* `mutable` as the `contains()` method is not required to be const */
      return [layout = layout](GeometryType gt, int) mutable {
        return layout.contains(gt);
      };
    }
  };
 
  //
  //  Leaf and level mapper
  //
 
  template <typename G, template<int> class LayoutClass = Impl::MCMGFailLayout>
  class LeafMultipleCodimMultipleGeomTypeMapper
    : public MultipleCodimMultipleGeomTypeMapper<typename G::LeafGridView,LayoutClass>
  {
    typedef MultipleCodimMultipleGeomTypeMapper<typename G::LeafGridView,
        LayoutClass> Base;
  public:
    LeafMultipleCodimMultipleGeomTypeMapper (const G& grid, const LayoutClass<G::dimension> layout = {})
      DUNE_DEPRECATED_MSG("Use the constructor taking a `MCMGLayout` functional instead")
      : LeafMultipleCodimMultipleGeomTypeMapper(grid, Base::wrapLayoutClass(layout))
    {}
 
    LeafMultipleCodimMultipleGeomTypeMapper (const G& grid, const MCMGLayout& layout)
      : Base(grid.leafGridView(), layout)
    {}
  };
 
  template <typename G, template<int> class LayoutClass = Impl::MCMGFailLayout>
  class LevelMultipleCodimMultipleGeomTypeMapper
    : public MultipleCodimMultipleGeomTypeMapper<typename G::LevelGridView,LayoutClass> {
    typedef MultipleCodimMultipleGeomTypeMapper<typename G::LevelGridView,
        LayoutClass> Base;
  public:
    LevelMultipleCodimMultipleGeomTypeMapper (const G& grid, int level, const LayoutClass<G::dimension> layout = {})
      DUNE_DEPRECATED_MSG("Use the constructor taking a `MCMGLayout` functional instead")
      : LevelMultipleCodimMultipleGeomTypeMapper(grid, level, Base::wrapLayoutClass(layout))
    {}
 
    LevelMultipleCodimMultipleGeomTypeMapper (const G& grid, int level, const MCMGLayout& layout)
      : Base(grid.levelGridView(level),layout)
    {}
  };
 
}
#endif