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 <iostream>
 
#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 MCMGElementLayout {
    bool contains (Dune::GeometryType gt) { return gt.dim()==dimgrid; }
  };
 
 
  template<int dim> struct MCMGVertexLayout {
    bool contains (Dune::GeometryType gt) { return gt.dim()==0; }
  };
 
  //
  //  MultipleCodimMultipleGeomTypeMapper
  //
 
  template <typename GV, template<int> class Layout>
  class MultipleCodimMultipleGeomTypeMapper :
    public Mapper<typename GV::Grid,MultipleCodimMultipleGeomTypeMapper<GV,Layout>, typename GV::IndexSet::IndexType >
  {
  public:
 
    typedef typename GV::IndexSet::IndexType Index;
 
    MultipleCodimMultipleGeomTypeMapper (const GV& gridView_, const Layout<GV::dimension> layout)
      : gridView(gridView_),
        is(gridView.indexSet()),
        offset(GlobalGeometryTypeIndex::size(GV::dimension)),
        layout(layout)
    {
      update();
    }
 
    MultipleCodimMultipleGeomTypeMapper (const GV& gridView_)
      : gridView(gridView_),
        is(gridView.indexSet()),
        offset(GlobalGeometryTypeIndex::size(GV::dimension))
    {
      update();
    }
 
    template<class EntityType>
    Index index (const EntityType& e) const
    {
      const GeometryType gt = e.type();
      assert(layout.contains(gt));
      return is.index(e) + offset[GlobalGeometryTypeIndex::index(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() ?
        GeometryType( GeometryType::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(layout.contains(gt));
      return is.subIndex(e, i, codim) + offset[GlobalGeometryTypeIndex::index(gt)];
    }
 
    int size () const
    {
      return n;
    }
 
    template<class EntityType>
    bool contains (const EntityType& e, Index& result) const
    {
      if(!is.contains(e) || !layout.contains(e.type()))
      {
        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() ?
        GeometryType( GeometryType::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 (not layout.contains(gt))
        return false;
      result = is.subIndex(e, i, cc) + offset[GlobalGeometryTypeIndex::index(gt)];
      return true;
    }
 
    void update ()
    {
      n = 0;
 
      for (unsigned int codim = 0; codim <= GV::dimension; ++codim)
      {
        // walk over all geometry types in the codimension
        typedef typename GV::IndexSet::Types GTV;
        GTV gtv = is.types(codim);
        for (typename GTV::const_iterator it = gtv.begin(); it != gtv.end(); ++it)
        {
          // if the geometry type is contained in the layout, increment offset
          if (layout.contains(*it))
          {
            offset[GlobalGeometryTypeIndex::index(*it)] = n;
            n += is.size(*it);
          }
        }
      }
    }
 
  private:
    // 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::vector<int> offset;
    mutable Layout<GV::dimension> layout;     // get layout object
  };
 
  //
  //  Leaf and level mapper
  //
 
  template <typename G, template<int> class Layout>
  class LeafMultipleCodimMultipleGeomTypeMapper
    : public MultipleCodimMultipleGeomTypeMapper<typename G::LeafGridView,Layout>
  {
    typedef MultipleCodimMultipleGeomTypeMapper<typename G::LeafGridView,
        Layout> Base;
  public:
    LeafMultipleCodimMultipleGeomTypeMapper (const G& grid)
      : Base(grid.leafGridView())
    {}
 
    LeafMultipleCodimMultipleGeomTypeMapper (const G& grid, const Layout<G::dimension> layout)
      : Base(grid.leafGridView(),layout)
    {}
 
  };
 
  template <typename G, template<int> class Layout>
  class LevelMultipleCodimMultipleGeomTypeMapper
    : public MultipleCodimMultipleGeomTypeMapper<typename G::LevelGridView,Layout> {
    typedef MultipleCodimMultipleGeomTypeMapper<typename G::LevelGridView,
        Layout> Base;
  public:
    LevelMultipleCodimMultipleGeomTypeMapper (const G& grid, int level)
      : Base(grid.levelGridView(level))
    {}
 
    LevelMultipleCodimMultipleGeomTypeMapper (const G& grid, int level, const Layout<G::dimension> layout)
      : Base(grid.levelGridView(level),layout)
    {}
 
  };
 
}
#endif