mcmgmapper.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- 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/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
  //

  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>
  class MultipleCodimMultipleGeomTypeMapper :
    public Mapper<typename GV::Grid,MultipleCodimMultipleGeomTypeMapper<GV>, 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 MCMGLayout& layout)
      : gridView_(gridView)
      , indexSet_(&gridView_.indexSet())
      , layout_(layout)
    {
      update(gridView);
    }

    template<class EntityType>
    Index index (const EntityType& e) const
    {
      const GeometryType gt = e.type();
      assert(offset(gt) != invalidOffset);
      return indexSet_->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 indexSet_->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(!indexSet_->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(!indexSet_->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) ;
      if (offset(gt) == invalidOffset)
        return false;
      result = indexSet_->subIndex(e, i, cc)*blockSize(gt) + offset(gt);
      return true;
    }

    void update (const GV& gridView)
    {
      gridView_ = gridView;
      indexSet_ = &gridView_.indexSet();
      update_();
    }

    void update (GV&& gridView)
    {
      gridView_ = std::move(gridView);
      indexSet_ = &gridView_.indexSet();
      update_();
    }

    const MCMGLayout &layout () const { return layout_; }
    const GridView &gridView () const { return gridView_; }

  private:
    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 : indexSet_->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 += indexSet_->size(gt) * block;
            myTypes_[codim].push_back(gt);
          }
          else {
            offset = invalidOffset;
          }

          offsets[GlobalGeometryTypeIndex::index(gt)] = offset;
          blocks[GlobalGeometryTypeIndex::index(gt)] = block;
        }
      }
    }

    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
    GV gridView_;
    const typename GV::IndexSet* indexSet_;
    // provide an array for the offsets
    std::array<Index, GlobalGeometryTypeIndex::size(GV::dimension)> offsets;
    std::array<Index, GlobalGeometryTypeIndex::size(GV::dimension)> blocks;
    MCMGLayout layout_;     // get layout object
    std::vector<GeometryType> myTypes_[GV::dimension+1];
  };

}
#endif