ygrid.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_YASPGRID_YGRID_HH
#define DUNE_GRID_YASPGRID_YGRID_HH
 
#include <array>
#include <vector>
#include <bitset>
#include <deque>
 
#include <dune/common/fvector.hh>
#include <dune/common/power.hh>
#include <dune/common/streamoperators.hh>
 
namespace Dune {
 
 namespace Yasp {
  template<int d, typename ct>
  std::array<int,d> sizeArray(const std::array<std::vector<ct>,d>& v)
  {
    std::array<int,d> tmp;
    for (int i=0; i<d; ++i)
      tmp[i] = v[i].size() - 1;
    return tmp;
  }
 } //namespace Yasp
 
  template<class Coordinates>
  class YGridComponent
  {
  public:
    //extract coordinate type and dimension from the coordinate container
    typedef typename Coordinates::ctype ct;
    static const int d = Coordinates::dimension;
 
    typedef std::array<int, d> iTupel;
    typedef FieldVector<ct,d> fTupel;
 
    YGridComponent () : _shift(0ULL)
    {
      std::fill(_origin.begin(), _origin.end(), 0);
      std::fill(_offset.begin(), _offset.end(), 0);
      std::fill(_size.begin(), _size.end(), 0);
    }
 
    YGridComponent(iTupel origin, iTupel size)
      : _origin(origin), _size(size)
    {}
 
    YGridComponent (iTupel origin, iTupel size, const YGridComponent<Coordinates>& enclosing)
      :  _origin(origin), _shift(enclosing.shift()), _coords(enclosing.getCoords()), _size(size), _supersize(enclosing.supersize())
    {
      for (int i=0; i<d; i++)
        _offset[i] = origin[i] - enclosing.origin(i) + enclosing.offset(i);
 
      // compute superincrements
      int inc = 1;
      for (int i=0; i<d; ++i)
        {
          _superincrement[i] = inc;
          inc *= _supersize[i];
        }
    }
 
    YGridComponent (iTupel origin, std::bitset<d> shift, Coordinates* coords, iTupel size, iTupel offset, iTupel supersize)
      : _origin(origin), _shift(shift), _coords(coords), _size(size), _offset(offset), _supersize(supersize)
    {
      // compute superincrements
      int inc = 1;
      for (int i=0; i<d; ++i)
        {
          _superincrement[i] = inc;
          inc *= _supersize[i];
        }
    }
 
    int origin (int i) const
    {
      return _origin[i];
    }
 
    const iTupel& origin () const
    {
      return _origin;
    }
 
    bool shift (int i) const
    {
      return _shift[i];
    }
 
    const std::bitset<d>& shift () const
    {
      return _shift;
    }
 
    Coordinates* getCoords() const
    {
      return _coords;
    }
 
    int offset (int i) const
    {
      return _offset[i];
    }
 
    const iTupel & offset () const
    {
      return _offset;
    }
 
    int supersize (int i) const
    {
      return _supersize[i];
    }
 
    const iTupel & supersize () const
    {
      return _supersize;
    }
 
    int size (int i) const
    {
      return _size[i];
    }
 
    iTupel size () const
    {
      return _size;
    }
 
    int totalsize () const
    {
      int s=1;
      for (int i=0; i<d; ++i)
        s *= size(i);
      return s;
    }
 
    int min (int i) const
    {
      return _origin[i];
    }
 
    int max (int i) const
    {
      return _origin[i] + size(i) - 1;
    }
 
    bool empty () const
    {
      for (int i=0; i<d; ++i)
      {
        if (size(i) == 0)
          return true;
      }
      return false;
    }
 
    bool inside (const iTupel& coord) const
    {
      for (int i=0; i<d; i++)
      {
        if ((coord[i]<_origin[i]) || (coord[i]>=_origin[i]+_size[i]))
          return false;
      }
      return true;
    }
 
    int index (const iTupel& coord) const
    {
      int index = (coord[d-1]-_origin[d-1]);
 
      for (int i=d-2; i>=0; i--)
        index = index*_size[i] + (coord[i]-_origin[i]);
 
      return index;
    }
 
    YGridComponent<Coordinates> move (iTupel v) const
    {
      for (int i=0; i<d; i++)
        v[i] += _origin[i];
      return YGridComponent<Coordinates>(v,_size,*this);
    }
 
    YGridComponent<Coordinates> intersection (const YGridComponent<Coordinates>& r) const
    {
      for (int i=0; i<d; i++)
      {
        //empty coordinate vectors result in empty intersections
        if (empty() || r.empty())
          return YGridComponent<Coordinates>();
      }
 
      iTupel neworigin;
      iTupel newsize;
      for (int i=0; i<d; ++i)
      {
        neworigin[i] = std::max(origin(i),r.origin(i));
        newsize[i] = std::min(max(i),r.max(i)) - neworigin[i] + 1;
      }
 
      return YGridComponent<Coordinates>(neworigin,newsize,*this);
    }
 
 
    class Iterator {
    public:
      // default constructor
      Iterator () = default;
 
      Iterator (const YGridComponent<Coordinates>& r) : _grid(&r)
      {
        iTupel coord(r.origin());
        reinit(r,coord);
      }
 
      Iterator (const YGridComponent<Coordinates>& r, const iTupel& coord)
      {
        reinit(r,coord);
      }
 
      void reinit (const YGridComponent<Coordinates>& r, const iTupel& coord)
      {
        // initialize to given position in index set
        for (int i=0; i<d; ++i)
          _coord[i] = coord[i];
 
        // move superindex to first cell in subgrid
        _superindex = 0;
        for (int i=0; i<d; ++i)
          _superindex += (r.offset(i)+coord[i]-r.origin(i))*r.superincrement(i);
 
        _grid = &r;
      }
 
      bool operator== (const Iterator& i) const
      {
        return _superindex == i._superindex;
      }
 
      bool operator!= (const Iterator& i) const
      {
        return _superindex != i._superindex;
      }
 
      int superindex () const
      {
        return _superindex;
      }
 
      int coord (int i) const
      {
        return _coord[i];
      }
 
      const iTupel& coord () const
      {
        return _coord;
      }
 
      void move (int i, int dist)
      {
        _coord[i] += dist;
        _superindex += dist*_grid->superincrement(i);
      }
 
      void move (const iTupel& dist)
      {
        for (int i = 0; i < d; ++i)
          {
            _coord[i] += dist[i];
            _superindex += dist[i]*_grid->superincrement(i);
          }
      }
 
      Iterator& operator++ ()
      {
        for (int i=0; i<d; i++)         // check for wrap around
        {
          _superindex += _grid->superincrement(i);   // move on cell in direction i
          if (++_coord[i] <= _grid->max(i))
            return *this;
          else
          {
            _coord[i] = _grid->origin(i);         // move back to origin in direction i
            _superindex -= _grid->size(i) * _grid->superincrement(i);
          }
        }
        // if we wrapped around, back to to begin(), we must put the iterator to end()
        if (_coord == _grid->origin())
        {
          for (int i=0; i<d; i++)
            _superindex += (_grid->size(i)-1) * _grid->superincrement(i);
          _superindex += _grid->superincrement(0);
        }
        return *this;
      }
 
      ct lowerleft(int i) const
      {
        return _grid->getCoords()->coordinate(i,_coord[i]);
      }
 
      fTupel lowerleft() const
      {
        fTupel ll;
        for (int i=0; i<d; i++)
          ll[i] = lowerleft(i);
        return ll;
      }
 
      ct upperright(int i) const
      {
        int coord = _coord[i];
        if (shift(i))
          coord++;
        return _grid->getCoords()->coordinate(i,coord);
      }
 
      fTupel upperright() const
      {
        fTupel ur;
        for (int i=0; i<d; i++)
          ur[i] = upperright(i);
        return ur;
      }
 
      ct meshsize (int i) const
      {
        return _grid->getCoords()->meshsize(i,_coord[i]);
      }
 
      fTupel meshsize () const
      {
        fTupel h;
        for (int i=0; i<d; i++)
          h[i] = meshsize(i);
        return h;
      }
 
      bool shift (int i) const
      {
        return _grid->shift(i);
      }
 
      std::bitset<d> shift() const
      {
        return _grid->shift();
      }
 
      Coordinates* coordCont() const
      {
        return _grid->getCoords();
      }
 
    protected:
      iTupel _coord;              
      int _superindex = 0;        
      const YGridComponent<Coordinates>* _grid = nullptr;
    };
 
 
    int superindex(iTupel coord) const
    {
      // move superindex to first cell in subgrid
      int si = 0;
      for (int i=0; i<d; ++i)
        si += (offset(i)+coord[i]-origin(i))*_superincrement[i];
      return si;
    }
 
    int superincrement(int i) const
    {
      return _superincrement[i];
    }
 
    Iterator begin () const
    {
      return Iterator(*this);
    }
 
    Iterator begin (const iTupel& co) const
    {
      return Iterator(*this,co);
    }
 
    Iterator end () const
    {
      iTupel last;
      for (int i=0; i<d; i++)
        last[i] = max(i);
      last[0] += 1;
      return Iterator(*this,last);
    }
 
  private:
    iTupel _origin;
    std::bitset<d> _shift;
    Coordinates* _coords;
    iTupel _size;
    iTupel _offset;    
    iTupel _supersize; 
    iTupel _superincrement; 
 
  };
 
 
  template <class Coordinates>
  inline std::ostream& operator<< (std::ostream& s, YGridComponent<Coordinates> e)
  {
    s << "Printing YGridComponent structure:" << std::endl;
    s << "Origin: " << e.origin() << std::endl;
    s << "Shift: " << e.shift() << std::endl;
    s << "Size: " << e.size() << std::endl;
    s << "Offset: " << e.offset() << std::endl;
    s << "Supersize: " << e.supersize() << std::endl;
    return s;
  }
 
  template <class Coordinates>
  inline std::ostream& operator<< (std::ostream& s, typename YGridComponent<Coordinates>::Iterator& e)
  {
    s << "Printing YGridComponent Iterator:" << std::endl << "Iterator at " << e.coord() << " (index ";
    s << e.index() << "), position " << e.position();
    return s;
  }
 
  template<class Coordinates>
  class YGrid
  {
    public:
    static const int dim = Coordinates::dimension;
 
    // define data array iterator
    typedef YGridComponent<Coordinates>* DAI;
 
    typedef typename std::array<int, dim> iTupel;
 
    void setBegin(DAI begin)
    {
      _begin = begin;
    }
 
    int shiftmapping(const std::bitset<dim>& shift) const
    {
      return _shiftmapping[shift.to_ulong()];
    }
 
    DAI dataBegin() const
    {
      return _begin;
    }
 
    DAI dataEnd() const
    {
      return _end;
    }
 
    bool inside(const iTupel& coord, const std::bitset<dim>& shift = std::bitset<dim>()) const
    {
      return (_begin+_shiftmapping[shift.to_ulong()])->inside(coord);
    }
 
    class Iterator
    {
      public:
 
      Iterator () = default;
 
      Iterator (const YGrid<Coordinates>& yg, const std::array<int,dim>& coords, int which = 0)
        : _which(which), _yg(&yg)
      {
        _it = typename YGridComponent<Coordinates>::Iterator(*(_yg->dataBegin()+which),coords);
      }
 
      Iterator (const YGrid<Coordinates>& yg, bool end=false) : _yg(&yg)
      {
        if (end)
        {
          _it = _yg->_itends.back();
          _which = _yg->_itends.size() - 1;
        }
        else
        {
          _it = _yg->_itbegins[0];
          _which = 0;
        }
      }
 
      void reinit(const YGrid<Coordinates>& yg, const std::array<int,dim>& coords, int which = 0)
      {
        _yg = &yg;
        _which = which;
        _it = typename YGridComponent<Coordinates>::Iterator(*(_yg->dataBegin()+which),coords);
      }
 
      int coord (int i) const
      {
        return _it.coord(i);
      }
 
      const std::array<int, dim>& coord () const
      {
        return _it.coord();
      }
 
      typename Coordinates::ctype lowerleft(int i) const
      {
        return _it.lowerleft(i);
      }
 
      Dune::FieldVector<typename Coordinates::ctype,dim> lowerleft() const
      {
        return _it.lowerleft();
      }
 
      typename Coordinates::ctype upperright(int i) const
      {
        return _it.upperright(i);
      }
 
      Dune::FieldVector<typename Coordinates::ctype,dim> upperright() const
      {
        return _it.upperright();
      }
 
      typename Coordinates::ctype meshsize (int i) const
      {
        return _it.meshsize(i);
      }
 
      Dune::FieldVector<typename Coordinates::ctype,dim> meshsize() const
      {
        return _it.meshsize();
      }
 
      bool shift (int i) const
      {
        return _it.shift(i);
      }
 
      std::bitset<dim> shift () const
      {
        return _it.shift();
      }
 
      int superindex() const
      {
        // the offset of the current component has to be taken into account
          return _yg->_indexOffset[_which] + _it.superindex();
      }
 
      Iterator& operator++ ()
      {
        if ((++_it == _yg->_itends[_which]) && (_which < _yg->_itends.size()-1))
          _it = _yg->_itbegins[++_which];
        return *this;
      }
 
      bool operator==(const Iterator& i) const
      {
        if (_which != i._which)
          return false;
        return _it == i._it;
      }
 
      bool operator!=(const Iterator& i) const
      {
        if (_it != i._it)
          return true;
        return _which != i._which;
      }
 
      int which() const
      {
        return _which;
      }
 
      void move(int i, int dist)
      {
        _it.move(i,dist);
      }
 
      void move(const iTupel& dist)
      {
        _it.move(dist);
      }
 
      Coordinates* coordCont() const
      {
        return _it.coordCont();
      }
 
 
      private:
      unsigned int _which = 0;
      const YGrid<Coordinates>* _yg = nullptr;
      typename YGridComponent<Coordinates>::Iterator _it;
    };
 
    Iterator begin() const
    {
      return Iterator(*this);
    }
 
    Iterator begin(const std::array<int, dim>& coord, int which = 0) const
    {
      return Iterator(*this, coord, which);
    }
 
    Iterator end() const
    {
      return Iterator(*this,true);
    }
 
    int superindex(const iTupel& coord, int which) const
    {
      return _indexOffset[which] + (dataBegin()+which)->superindex(coord);
    }
 
 
    // finalize the ygrid construction by storing component iterators
    void finalize(const DAI& end, int artificialOffset = 0)
    {
      // set the end iterator in the ygrid component array
      _end = end;
 
      _indexOffset.push_back(artificialOffset);
      int k = 0;
      for (DAI i=_begin; i != _end; ++i, ++k)
      {
        //store begin and end iterators
        _itbegins.push_back(i->begin());
        _itends.push_back(i->end());
 
        // store index offset
        _indexOffset.push_back(_indexOffset.back() + i->totalsize());
 
        // store shift to component mapping
        _shiftmapping[i->shift().to_ulong()] = k;
      }
      _indexOffset.resize(_itends.size());
    }
 
    private:
 
    friend class YGrid<Coordinates>::Iterator;
    DAI _begin;
    DAI _end;
    std::array<int,StaticPower<2,dim>::power> _shiftmapping;
    std::vector<typename YGridComponent<Coordinates>::Iterator> _itbegins;
    std::vector<typename YGridComponent<Coordinates>::Iterator> _itends;
    std::vector<int> _indexOffset;
  };
 
  template <class Coordinates>
  inline std::ostream& operator<< (std::ostream& s, const YGrid<Coordinates>& e)
  {
    s << "Printing YGrid structure:" << std::endl;
    for (auto it = e.dataBegin(); it != e.dataEnd(); ++it)
      s << *it << std::endl;
    return s;
  }
 
  template<class Coordinates>
  class YGridList
  {
    public:
    static const int dim = Coordinates::dimension;
 
    struct Intersection
    {
      YGridComponent<Coordinates> grid;
      int rank;
      int distance;
      YGrid<Coordinates> yg;
    };
 
    // define data array iterator type
    typedef typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator DAI;
 
    // iterator that allows to iterate over a concatenation of deques. namely those
    // that belong to the same codimension.
    class Iterator
    {
      public:
 
        Iterator(const YGridList<Coordinates>& ygl, bool end=false) : _end(ygl.dataEnd()), _which(ygl.dataBegin())
      {
        _it = _which->begin();
 
        // advance the iterator to the first element that exists.
        // some deques might be empty and should be skipped
        while ((_which != _end) && (_it == _which->end()))
        {
          ++_which;
          if (_which != _end)
            _it = _which->begin();
        }
        // the iterator is at the end if and only if _which==_end
        if (end)
        {
          _which = _end;
        }
      }
 
      Iterator& operator++ ()
      {
        ++_it;
        // advance the iterator to the next element that exists.
        // some deques might be empty and should be skipped
        while ((_which != _end) && (_it == _which->end()))
        {
          ++_which;
          if (_which != _end)
            _it = _which->begin();
        }
        return *this;
      }
 
      typename std::deque<Intersection>::iterator  operator->() const
      {
        return _it;
      }
 
      typename std::deque<Intersection>::iterator  operator*() const
      {
        return _it;
      }
 
      bool operator== (const Iterator& i) const
      {
        if (_which != i._which)
          return false;
        if (_which == _end)
          return true;
        return _it == i._it;
      }
 
      bool operator!= (const Iterator& i) const
      {
        if (_which != i._which)
          return true;
        if (_which == _end)
          return false;
        return _it != i._it;
      }
 
      private:
      typename std::deque<Intersection>::iterator _it;
      DAI _end;
      DAI _which;
    };
 
    Iterator begin() const
    {
      return Iterator(*this);
    }
 
    Iterator end() const
    {
      return Iterator(*this,true);
    }
 
    void setBegin(typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator begin)
    {
      _begin = begin;
    }
 
    DAI dataBegin() const
    {
      return _begin;
    }
 
    DAI dataEnd() const
    {
      return _end;
    }
 
    int size() const
    {
      int count = 0;
      for (DAI it = _begin; it != _end; ++it)
        count += it->size();
      return count;
    }
 
    void finalize(DAI end, const YGrid<Coordinates>& ygrid)
    {
      // Instead of directly iterating over the intersection deques, this code
      // iterates over the components of an associated ygrid and works its way
      // through the list of intersection deques in parallel.
      // The reason for this convoluted iteration technique is that there are not
      // necessarily intersections for all possible shifts, but we have to make
      // sure that we stop at each shift to update the per-component index shift.
      // This is ensured by iterating over the ygrid, which is guaranteed to have
      // a component for each shift vector.
 
      // set end iterator in the data array
      _end = end;
 
      int offset = 0;
 
      DAI i = _begin;
 
      // make sure that we have a valid deque (i.e. a non-empty one)
      while (i != _end && i->begin() == i->end())
        ++i;
 
      for (auto yit = ygrid.dataBegin(); yit != ygrid.dataEnd(); ++yit)
      {
        if (i == _end)
          break;
        auto it = i->begin();
        if (it->grid.shift() == yit->shift())
        {
          // iterate over the intersections in the deque and set the offset
          for (; it != i->end(); ++it)
          {
            it->yg.setBegin(&(it->grid));
            it->yg.finalize(&(it->grid)+1, offset);
          }
 
          // advance to next non-empty deque
          ++i;
          while (i != _end && i->begin() == i->end())
            ++i;
        }
 
        // update the offset from the ygrid component
        int add = 1;
        for (int j=0; j<dim; j++)
          add *= yit->supersize(j);
        offset += add;
      }
      assert (i == end);
    }
 
    private:
    DAI _begin;
    DAI _end;
  };
 
} // namespace Dune
 
#endif