Orthant.h

Go to the documentation of this file.

/***************************************************************************
 *  ThunderEgg, a library for solvers on adaptively refined block-structured
 *  Cartesian grids.
 *
 *  Copyright (c) 2020-2021 Scott Aiton
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 ***************************************************************************/
 
#ifndef THUNDEREGG_ORTHANT_H
#define THUNDEREGG_ORTHANT_H
 
#include <ThunderEgg/Face.h>
#include <array>
#include <iostream>
#include <numeric>
#include <vector>
namespace ThunderEgg {
template<int D>
class Orthant
{
private:
  unsigned char val = num_orthants;
 
public:
  static constexpr size_t num_orthants = 1 << D;
  explicit Orthant(const unsigned char val_in)
    : val(val_in)
  {}
  Orthant() {}
  static Orthant<D> null() { return Orthant<D>(num_orthants); }
  static Orthant<1> lower() { return Orthant<1>(0b0); }
  static Orthant<1> upper() { return Orthant<1>(0b1); }
  static Orthant<2> sw() { return Orthant<2>(0b00); }
  static Orthant<2> se() { return Orthant<2>(0b01); }
  static Orthant<2> nw() { return Orthant<2>(0b10); }
  static Orthant<2> ne() { return Orthant<2>(0b11); }
  static Orthant<3> bsw() { return Orthant<3>(0b000); }
  static Orthant<3> bse() { return Orthant<3>(0b001); }
  static Orthant<3> bnw() { return Orthant<3>(0b010); }
  static Orthant<3> bne() { return Orthant<3>(0b011); }
  static Orthant<3> tsw() { return Orthant<3>(0b100); }
  static Orthant<3> tse() { return Orthant<3>(0b101); }
  static Orthant<3> tnw() { return Orthant<3>(0b110); }
  static Orthant<3> tne() { return Orthant<3>(0b111); }
  class Range
  {
  public:
    class Iterator : public std::input_iterator_tag
    {
    private:
      Orthant<D> o;
 
    public:
      explicit Iterator(Orthant<D> o_in)
        : o(o_in)
      {}
      const Orthant<D>& operator++()
      {
        ++o.val;
        return o;
      }
      const Orthant<D>& operator*() const { return o; }
      const Orthant<D>* operator->() const { return &o; }
      bool operator==(const Iterator& b) const { return o.val == b.o.val; }
      bool operator!=(const Iterator& b) const { return o.val != b.o.val; }
    };
    Iterator begin() { return Iterator(Orthant<D>(0)); }
    Iterator end() { return Iterator(null()); }
  };
  static Range getValues() { return Range(); }
  size_t getIndex() const { return val; }
  Orthant<D> getNbrOnSide(Side<D> s) const
  {
    Orthant<D> retval = *this;
    // flip the bit for that side
    retval.val ^= (0x1 << s.getAxisIndex());
    return retval;
  }
  std::array<Side<D>, D> getInteriorSides() const
  {
    std::array<Side<D>, D> retval;
    for (size_t i = 0; i < D; i++) {
      size_t side = 2 * i;
      if (!((1 << i) & val)) {
        side |= 1;
      }
      retval[i] = Side<D>(side);
    }
    return retval;
  }
  std::array<Side<D>, D> getExteriorSides() const
  {
    std::array<Side<D>, D> retval;
    for (size_t i = 0; i < D; i++) {
      size_t side = 2 * i;
      if ((1 << i) & val) {
        side |= 1;
      }
      retval[i] = Side<D>(side);
    }
    return retval;
  }
  bool isOnSide(Side<D> s) const
  {
    int idx = s.getIndex() / 2;
    int remainder = s.getIndex() % 2;
    bool is_bit_set = val & (0x1 << idx);
    return is_bit_set == remainder;
  }
  bool isHigherOnAxis(size_t axis) const { return val & (0b1 << axis); }
  bool isLowerOnAxis(size_t axis) const { return !(val & (0b1 << axis)); }
  Orthant<D - 1> collapseOnAxis(size_t axis) const
  {
    size_t upper_mask = (~0x0U) << axis;
    return Orthant<D - 1>(((val >> 1) & upper_mask) | (val & ~upper_mask));
  }
  static std::array<Orthant, num_orthants / 2> getValuesOnSide(Side<D> s)
  {
    unsigned int bit_to_insert = s.getAxisIndex();
    unsigned int set_bit = s.isLowerOnAxis() ? 0 : 1;
    unsigned int lower_mask = ~((~0x0U) << bit_to_insert);
    unsigned int upper_mask = (~0x0U) << (bit_to_insert + 1);
 
    std::array<Orthant<D>, Orthant<D>::num_orthants / 2> retval;
    for (size_t i = 0; i < Orthant<D>::num_orthants / 2; i++) {
      size_t value = (i << 1) & upper_mask;
      value |= i & lower_mask;
      value |= set_bit << bit_to_insert;
      retval[i] = Orthant<D>(value);
    }
    return retval;
  }
  bool operator==(const Orthant<D>& other) const { return val == other.val; }
  bool operator!=(const Orthant<D>& other) const { return val != other.val; }
  bool operator<(const Orthant<D>& other) const { return val < other.val; }
};
 
inline std::ostream&
operator<<(std::ostream& os, const Orthant<0>& o)
{
  if (o == Orthant<0>::null()) {
    os << "Orthant<0>::null()";
  } else {
    os << "Orthant<0> invalid value: " << o.getIndex();
  }
  return os;
}
inline std::ostream&
operator<<(std::ostream& os, const Orthant<1>& o)
{
  if (o == Orthant<1>::lower()) {
    os << "Orthant<1>::lower()";
  } else if (o == Orthant<1>::upper()) {
    os << "Orthant<1>::upper()";
  } else if (o == Orthant<1>::null()) {
    os << "Orthant<1>::null()";
  } else {
    os << "Orthant<1> invalid value: " << o.getIndex();
  }
  return os;
}
inline std::ostream&
operator<<(std::ostream& os, const Orthant<2>& o)
{
  if (o == Orthant<2>::sw()) {
    os << "Orthant<2>::sw()";
  } else if (o == Orthant<2>::se()) {
    os << "Orthant<2>::se()";
  } else if (o == Orthant<2>::nw()) {
    os << "Orthant<2>::nw()";
  } else if (o == Orthant<2>::ne()) {
    os << "Orthant<2>::ne()";
  } else if (o == Orthant<2>::null()) {
    os << "Orthant<2>::null()";
  } else {
    os << "Orthant<2> invalid value: " << o.getIndex();
  }
  return os;
}
inline std::ostream&
operator<<(std::ostream& os, const Orthant<3>& o)
{
  if (o == Orthant<3>::bsw()) {
    os << "Orthant<3>::bsw()";
  } else if (o == Orthant<3>::bse()) {
    os << "Orthant<3>::bse()";
  } else if (o == Orthant<3>::bnw()) {
    os << "Orthant<3>::bnw()";
  } else if (o == Orthant<3>::bne()) {
    os << "Orthant<3>::bne()";
  } else if (o == Orthant<3>::tsw()) {
    os << "Orthant<3>::tsw()";
  } else if (o == Orthant<3>::tse()) {
    os << "Orthant<3>::tse()";
  } else if (o == Orthant<3>::tnw()) {
    os << "Orthant<3>::tnw()";
  } else if (o == Orthant<3>::tne()) {
    os << "Orthant<3>::tne()";
  } else if (o == Orthant<3>::null()) {
    os << "Orthant<3>::null()";
  } else {
    os << "Orthant<3> invalid value: " << o.getIndex();
  }
  return os;
}
void
to_json(tpl::nlohmann::json& j, const Orthant<0>& o);
void
to_json(tpl::nlohmann::json& j, const Orthant<1>& o);
void
to_json(tpl::nlohmann::json& j, const Orthant<2>& o);
void
to_json(tpl::nlohmann::json& j, const Orthant<3>& o);
void
from_json(const tpl::nlohmann::json& j, Orthant<0>& o);
void
from_json(const tpl::nlohmann::json& j, Orthant<1>& o);
void
from_json(const tpl::nlohmann::json& j, Orthant<2>& o);
void
from_json(const tpl::nlohmann::json& j, Orthant<3>& o);
} // namespace ThunderEgg
#endif