generalvertexorder.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_GEOMETRY_GENERALVERTEXORDER_HH
#define DUNE_GEOMETRY_GENERALVERTEXORDER_HH
 
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <iterator>
#include <vector>
 
#include <dune/common/iteratorfacades.hh>
 
#include "type.hh"
#include <dune/geometry/referenceelements.hh>
 
namespace Dune {
 
  template<class InIterator, class OutIterator>
  void reduceOrder(const InIterator& inBegin, const InIterator& inEnd,
                   OutIterator outIt)
  {
    for(InIterator inIt = inBegin; inIt != inEnd; ++inIt, ++outIt)
      *outIt = std::count_if(inBegin, inEnd, [&](const auto& v)
                             {
                               return v < *inIt;
                             });
  }
 
 
  template<std::size_t dim, class Index_ = std::size_t>
  class GeneralVertexOrder {
    typedef ReferenceElements<double, dim> RefElems;
    typedef typename RefElems::ReferenceElement RefElem;
 
    RefElem refelem;
    GeometryType gt;
    std::vector<Index_> vertexOrder;
 
  public:
    typedef Index_ Index;
 
    class iterator;
 
    static const std::size_t dimension = dim;
    const GeometryType &type() const { return gt; }
 
 
    template<class InIterator>
    GeneralVertexOrder(const GeometryType& gt_, const InIterator &inBegin,
                       const InIterator &inEnd) :
      refelem(RefElems::general(gt_)), gt(gt_),
      vertexOrder(refelem.size(dim))
    { reduceOrder(inBegin, inEnd, vertexOrder.begin()); }
 
 
    iterator begin(std::size_t codim, std::size_t subEntity) const
    { return iterator(*this, codim, subEntity); }
 
    iterator end(std::size_t codim, std::size_t subEntity) const {
      return iterator(*this, codim, subEntity,
                      refelem.size(subEntity, codim, dim));
    }
 
 
    void getReduced(std::size_t codim, std::size_t subEntity,
                    std::vector<Index>& order) const
    {
      order.resize(refelem.size(subEntity, codim, dim));
      reduceOrder(begin(codim, subEntity), end(codim, subEntity),
                  order.begin());
    }
  };
 
 
  template<std::size_t dim, class Index_>
  class GeneralVertexOrder<dim, Index_>::iterator :
    public Dune::RandomAccessIteratorFacade<iterator, const Index_>
  {
    const GeneralVertexOrder *order;
    std::size_t codim;
    std::size_t subEntity;
    std::size_t vertex;
 
    iterator(const GeneralVertexOrder &order_, std::size_t codim_,
             std::size_t subEntity_, std::size_t vertex_ = 0) :
      order(&order_), codim(codim_), subEntity(subEntity_), vertex(vertex_)
    { }
 
  public:
    const Index &dereference() const {
      return order->vertexOrder[order->refelem.subEntity(subEntity, codim,
                                                         vertex, dim)];
    }
    const Index &elementAt(std::ptrdiff_t n) const {
      return order->vertexOrder[order->refelem.subEntity(subEntity, codim,
                                                         vertex+n, dim)];
    }
    bool equals(const iterator &other) const {
      return order == other.order && codim == other.codim &&
             subEntity == other.subEntity && vertex == other.vertex;
    }
    void increment() { ++vertex; }
    void decrement() { --vertex; }
    void advance(std::ptrdiff_t n) { vertex += n; }
    std::ptrdiff_t distanceTo(const iterator &other) const {
      // make sure we reference the same container
      assert(order == other.order && codim == other.codim &&
             subEntity == other.subEntity);
      if(vertex < other.vertex) return other.vertex - vertex;
      else return -static_cast<std::ptrdiff_t>(vertex - other.vertex);
    }
 
    friend class GeneralVertexOrder<dim, Index>;
 
 
    iterator() { }
  };
} // namespace Dune
 
#endif // DUNE_GEOMETRY_GENERALVERTEXORDER_HH