CellIndex.h

#pragma once
 
#include <iostream>
#include <iterator>
#include <string_view>
#include <type_traits>
 
#include "tarch/la/Vector.h"
 
namespace coupling {
namespace indexing {

enum class IndexTrait { vector, local, md2macro, noGhost };
namespace TraitOperations {

template <IndexTrait t1> constexpr bool is_same(const IndexTrait& t2) { return t1 == t2; }

template <IndexTrait t1, IndexTrait t2, IndexTrait... rest> constexpr bool is_ordered() {
  if constexpr (sizeof...(rest) == 0) {
    return t1 < t2;
  } else {
    return t1 < t2 and is_ordered<t2, rest...>();
  }
}
 
} // namespace TraitOperations

template <unsigned int dim, IndexTrait... traits> class CellIndex;

template <unsigned int dim> using BaseIndex = CellIndex<dim, IndexTrait::vector>;
 
// TODO: refactor as member functions
template <unsigned int dim, IndexTrait... traits> unsigned int convertToScalar(const CellIndex<dim, traits...>&);
 
template <unsigned int dim, IndexTrait... traits> tarch::la::Vector<dim, int> convertToVector(const CellIndex<dim, traits...>&);
} // namespace indexing
} // namespace coupling
 
template <unsigned int dim, coupling::indexing::IndexTrait... traits> class coupling::indexing::CellIndex {
public:
  static constexpr bool checkIndexTraitOrder() {
    if constexpr (sizeof...(traits) > 1) {
      return coupling::indexing::TraitOperations::is_ordered<traits...>();
    } else {
      return true;
    }
  }
  static_assert(checkIndexTraitOrder(), "Invalid order in IndexTrait parameter pack! Correct oder: "
                                        "IndexTrait::vector < IndexTrait::local < IndexTrait::md2macro "
                                        "< IndexTrait::noGhost");

  using value_T = std::conditional_t<(coupling::indexing::TraitOperations::is_same<coupling::indexing::IndexTrait::vector>(traits) or ...),
                                     tarch::la::Vector<dim, int>, unsigned int>;

  CellIndex() = default;
  CellIndex(const CellIndex& ci) = default;
  explicit CellIndex(const value_T& i) : _index(i) {}

  template <coupling::indexing::IndexTrait... converted_traits> operator CellIndex<dim, converted_traits...>() const;

  value_T get() const { return _index; } // TODO: why (value_T) cast?

  CellIndex& operator++() {
    if constexpr (std::is_same_v<value_T, tarch::la::Vector<dim, int>>) {
      ++_index[0];
      if (_index[0] == (int)numberCellsInDomain[0]) {
        _index[0] = 0;
        ++_index[1];
        if (_index[1] == (int)numberCellsInDomain[1]) {
          _index[1] = 0;
          if constexpr (dim == 3)
            ++_index[2];
        }
      }
    } else
      ++_index;
 
    return *this;
  }

  CellIndex& operator--() {
    if constexpr (std::is_same_v<value_T, tarch::la::Vector<dim, int>>) {
      --_index[0];
      if (_index[0] < 0) {
        _index[0] = numberCellsInDomain[0] - 1;
        --_index[1];
        if (_index[1] < 0) {
          _index[1] = numberCellsInDomain[1] - 1;
          if constexpr (dim == 3)
            --_index[2];
        }
      }
    } else
      --_index;
 
    return *this;
  }
 
  /*
   * Any two indices fulfill some relation iff the unsigned integers underlying
   * their CellIndex<dim, {}> equivalents fulfill that relation. Note that this
   * is different than what the static member 'contains(...)' does.
   *
   * @param CellIndex index to compare this index to
   */
  bool operator==(const CellIndex& i) const { return _index == i.get(); }
  bool operator!=(const CellIndex& i) const { return not(i == *this); }
  bool operator<(const CellIndex& i) const { return convertToScalar<dim, traits...>(*this) < convertToScalar<dim, traits...>(i); };
  bool operator<=(const CellIndex& i) const { return convertToScalar<dim, traits...>(*this) <= convertToScalar<dim, traits...>(i); };
  bool operator>(const CellIndex& i) const { return convertToScalar<dim, traits...>(*this) > convertToScalar<dim, traits...>(i); };
  bool operator>=(const CellIndex& i) const { return convertToScalar<dim, traits...>(*this) >= convertToScalar<dim, traits...>(i); };

  static void setDomainParameters() {
    auto cellnum = upperBoundary.get() - lowerBoundary.get() + tarch::la::Vector<dim, int>{1};
    for (unsigned int d = 0; d < dim; d++)
      cellnum[d] = std::max(0, cellnum[d]);
    numberCellsInDomain = tarch::la::Vector<dim, unsigned int>{cellnum};
 
    linearNumberCellsInDomain = 1;
    for (unsigned int d = 0; d < dim; d++)
      linearNumberCellsInDomain *= numberCellsInDomain[d];
 
    tarch::la::Vector<dim, unsigned int> divFactor{1};
    for (unsigned int d = 1; d < dim; d++)
      divFactor[d] = divFactor[d - 1] * (numberCellsInDomain[d - 1]);
    divisionFactor = divFactor;
  }

  static bool contains(const coupling::indexing::BaseIndex<dim>& index) {
    for (unsigned int d = 0; d < dim; d++) {
      if (index.get()[d] < CellIndex::lowerBoundary.get()[d])
        return false;
      if (index.get()[d] > CellIndex::upperBoundary.get()[d])
        return false;
    }
 
    return true;
  }
 
  tarch::la::Vector<dim, double> getCellMidPoint() const {
    BaseIndex<dim> globalIndex{*this};
    tarch::la::Vector<dim, double> cellMidPoint(IndexingService<dim>::getInstance().getGlobalMDDomainOffset() -
                                                0.5 * IndexingService<dim>::getInstance().getCouplingCellSize());
    for (unsigned int d = 0; d < dim; d++)
      cellMidPoint[d] += ((double)(globalIndex.get()[d])) * IndexingService<dim>::getInstance().getCouplingCellSize()[d];
    return cellMidPoint;
  }

  static BaseIndex<dim> lowerBoundary;
  static BaseIndex<dim> upperBoundary;

  static tarch::la::Vector<dim, unsigned int> numberCellsInDomain;

  static unsigned int linearNumberCellsInDomain;

  static tarch::la::Vector<dim, unsigned int> divisionFactor;
 
  class IndexIterator {
  public:
    IndexIterator(CellIndex x) : _idx(x) {}
    IndexIterator(const IndexIterator& a) : _idx(a._idx) {}
 
    const CellIndex& operator*() const { return _idx; }
    const CellIndex* operator->() const { return &_idx; }
 
    // Prefix increment
    IndexIterator& operator++() {
      ++_idx;
      return *this;
    }
 
    // Postfix increment
    IndexIterator operator++(int) {
      IndexIterator tmp = *this;
      ++(*this);
      return tmp;
    }
 
    friend bool operator==(const IndexIterator& a, const IndexIterator& b) { return a._idx == b._idx; };
    friend bool operator!=(const IndexIterator& a, const IndexIterator& b) { return a._idx != b._idx; };
 
  private:
    CellIndex _idx;
  };
  static IndexIterator begin() {
    // if constexpr is neccessary here (otherwise benchmark performance breaks, due to runtime if)
    if constexpr (std::is_same_v<unsigned int, value_T>) {
      // If this CellIndex is scalar and if the domain is empty, then lowerBoundary can not be converted to CellIndex
      // Thus we have to return something else as an iterator for the empty domain
      if (linearNumberCellsInDomain == 0)
        return IndexIterator(CellIndex{});
      else
        return IndexIterator(lowerBoundary);
    } else
      return IndexIterator(lowerBoundary);
  }
  static IndexIterator end() {
    // Todo: This if-else branching is quite slow? (~factor 2)
    if (linearNumberCellsInDomain == 0)
      return begin();
    else
      return ++IndexIterator(upperBoundary);
  }
 
  static const char TNAME[];
 
private:
  value_T _index;
 
  template <coupling::indexing::IndexTrait T1, coupling::indexing::IndexTrait... other_traits> CellIndex<dim, other_traits...> getScalarCellIndex(int value) {
    static_assert(T1 == coupling::indexing::IndexTrait::vector);
    return CellIndex<dim, other_traits...>(value);
  }
};
 
// Include implementation
#include "CellIndex.cpph"