tpl_ca_ghost_lattice.H

Go to the documentation of this file.

/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
#ifndef TPL_CA_GHOST_LATTICE_H
#define TPL_CA_GHOST_LATTICE_H
 
#include <cstddef>
#include <type_traits>
#include <utility>
 
#include <ah-errors.H>
 
#include <ca-traits.H>
#include <tpl_ca_concepts.H>
#include <tpl_ca_lattice.H>  // for ca_lattice_detail::is_constant_boundary_v
 
namespace Aleph {
namespace CA {
 
namespace ca_ghost_detail {
 
template <std::size_t N>
inline std::array<ca_size_t, N> enlarge_extents(const std::array<ca_size_t, N> &user,
                                                const std::size_t halo) noexcept
{
  std::array<ca_size_t, N> e{};
  for (std::size_t d = 0; d < N; ++d)
    e[d] = user[d] == 0 ? 0 : user[d] + 2 * halo;
  return e;
}
 
inline ca_size_t positive_mod(ca_index_t v, const ca_size_t n) noexcept
{
  const ca_index_t m = static_cast<ca_index_t>(n);
  ca_index_t r = v % m;
  if (r < 0)
    r += m;
  return static_cast<ca_size_t>(r);
}
 
inline ca_size_t reflect_index(const ca_index_t v, const ca_size_t n) noexcept
{
  if (n == 1)
    return 0;
  const ca_index_t period = static_cast<ca_index_t>(2 * n);
  ca_index_t r = v % period;
  if (r < 0)
    r += period;
  const ca_index_t ni = static_cast<ca_index_t>(n);
  if (r >= ni)
    r = period - 1 - r;
  return static_cast<ca_size_t>(r);
}
 
inline ca_size_t clamp_index(const ca_index_t v, const ca_size_t n) noexcept
{
  if (v < 0)
    return 0;
  const ca_size_t u = static_cast<ca_size_t>(v);
  return u >= n ? n - 1 : u;
}
 
}  // namespace ca_ghost_detail
 
template <typename Storage, typename Boundary = OpenBoundary, std::size_t Halo = 1>
class Ghost_Lattice
{
  static_assert(Halo >= 1, "Ghost_Lattice requires Halo >= 1");
 
public:
  using storage_type = Storage;
  using boundary_type = Boundary;
  using state_type = typename Storage::state_type;
  using coord_type = typename Storage::coord_type;
  using extents_type = typename Storage::extents_type;
 
  static constexpr std::size_t rank = Storage::rank;
  static constexpr std::size_t halo_v = Halo;
 
private:
  extents_type user_ext_{};   
  extents_type store_ext_{};  
  Storage store_{};           
 
  [[nodiscard]] coord_type to_store(const coord_type &c) const noexcept
  {
    coord_type s = c;
    for (std::size_t d = 0; d < rank; ++d)
      s[d] = c[d] + static_cast<ca_index_t>(Halo);
    return s;
  }
 
  [[nodiscard]] state_type out_of_range_value() const
  {
    if constexpr (ca_lattice_detail::is_constant_boundary_v<Boundary>)
      return static_cast<state_type>(Boundary::value);
    else
      return state_type{};
  }
 
  [[nodiscard]] state_type resolve_via_policy(coord_type c) const
  {
    using namespace ca_ghost_detail;
    for (std::size_t d = 0; d < rank; ++d)
      {
        const ca_size_t n = user_ext_[d];
        if (n == 0)
          return out_of_range_value();
 
        if constexpr (std::is_same_v<Boundary, ToroidalBoundary>)
          {
            c[d] = static_cast<ca_index_t>(positive_mod(c[d], n));
          }
        else if constexpr (std::is_same_v<Boundary, ReflectiveBoundary>)
          {
            c[d] = static_cast<ca_index_t>(reflect_index(c[d], n));
          }
        else if constexpr (std::is_same_v<Boundary, NeumannBoundary>)
          {
            c[d] = static_cast<ca_index_t>(clamp_index(c[d], n));
          }
        else
          {
            if (c[d] < 0 or static_cast<ca_size_t>(c[d]) >= n)
              return out_of_range_value();
          }
      }
    return store_.at(to_store(c));
  }
 
  // -----------------------------------------------------------------
  // Halo fillers (one per boundary policy). Each walks the halo cells
  // and writes the right value into the backing storage. They all use
  // the same helper `for_each_halo_cell` because the structure of the
  // halo is identical across policies.
  // -----------------------------------------------------------------
 
  template <typename F>
  void for_each_halo_cell(F &&f) const
  {
    for (std::size_t d = 0; d < rank; ++d)
      if (user_ext_[d] == 0)
        return;
 
    coord_type u{};
    const ca_index_t h = static_cast<ca_index_t>(Halo);
 
    // Walk the enlarged box as a disjoint shell: until a halo band is
    // selected, recurse through interior plus the two halo bands; after
    // that, the remaining dimensions may span their full enlarged range.
    auto enumerate = [&](auto &&self_fn, std::size_t d, bool in_halo) -> void
    {
      if (d == rank)
        {
          if (in_halo)
            f(u);
          return;
        }
      const ca_index_t n = static_cast<ca_index_t>(user_ext_[d]);
      if (n == 0)
        return;
 
      if (in_halo)
        {
          for (ca_index_t i = -h; i < n + h; ++i)
            {
              u[d] = i;
              self_fn(self_fn, d + 1, true);
            }
          return;
        }
 
      if (d + 1 < rank)
        {
          for (ca_index_t i = 0; i < n; ++i)
            {
              u[d] = i;
              self_fn(self_fn, d + 1, false);
            }
        }
      
      if (h > 0)
        {
          for (ca_index_t i = -h; i < 0; ++i)
            {
              u[d] = i;
              self_fn(self_fn, d + 1, true);
            }
          for (ca_index_t i = n; i < n + h; ++i)
            {
              u[d] = i;
              self_fn(self_fn, d + 1, true);
            }
        }
    };
    enumerate(enumerate, 0, false);
  }
 
  void store_set(const coord_type &user_c, const state_type &value)
  {
    store_.set(to_store(user_c), value);
  }
 
  [[nodiscard]] state_type store_get(const coord_type &user_c) const
  {
    return store_.at(to_store(user_c));
  }
 
public:
  Ghost_Lattice() = default;
 
  explicit Ghost_Lattice(const extents_type &user_extents,
                         const state_type &init = state_type{})
    : user_ext_(user_extents)
    , store_ext_(ca_ghost_detail::enlarge_extents<rank>(user_extents, Halo))
    , store_(store_ext_, init)
  {}
 
  // ------------------------------------------------------------------
  // LatticeLike surface (see `tpl_ca_concepts.H`).
  // ------------------------------------------------------------------
 
  [[nodiscard]] static constexpr std::size_t dimension() noexcept { return rank; }
 
  [[nodiscard]] ca_size_t size() const noexcept
  {
    ca_size_t p = 1;
    for (std::size_t d = 0; d < rank; ++d)
      p *= user_ext_[d];
    return p;
  }
 
  [[nodiscard]] ca_size_t size(std::size_t d) const
  {
    ah_out_of_range_error_if(d >= rank)
      << "Ghost_Lattice::size: axis " << d << " out of range";
    return user_ext_[d];
  }
 
  [[nodiscard]] const extents_type &extents() const noexcept { return user_ext_; }
 
  [[nodiscard]] const extents_type &store_extents() const noexcept { return store_ext_; }
 
  [[nodiscard]] static constexpr std::size_t halo_radius() noexcept { return Halo; }
 
  [[nodiscard]] state_type at(const coord_type &c) const
  {
    for (std::size_t d = 0; d < rank; ++d)
      ah_out_of_range_error_if(c[d] < 0 or static_cast<ca_size_t>(c[d]) >= user_ext_[d])
        << "Ghost_Lattice::at: coord[" << d << "]=" << c[d] << " out of [0, "
        << user_ext_[d] << ")";
    return store_get(c);
  }
 
  void set(const coord_type &c, const state_type &v)
  {
    for (std::size_t d = 0; d < rank; ++d)
      ah_out_of_range_error_if(c[d] < 0 or static_cast<ca_size_t>(c[d]) >= user_ext_[d])
        << "Ghost_Lattice::set: coord[" << d << "]=" << c[d] << " out of [0, "
        << user_ext_[d] << ")";
    store_set(c, v);
  }
 
  [[nodiscard]] state_type at_safe(const coord_type &c) const
  {
    const ca_index_t h = static_cast<ca_index_t>(Halo);
    bool in_halo_range = true;
    for (std::size_t d = 0; d < rank; ++d)
      {
        if (user_ext_[d] == 0)
          return out_of_range_value();
        const ca_index_t n = static_cast<ca_index_t>(user_ext_[d]);
        if (c[d] < -h or c[d] >= n + h)
          {
            in_halo_range = false;
            break;
          }
      }
    if (in_halo_range)
      return store_get(c);
    return resolve_via_policy(c);
  }
 
  // ------------------------------------------------------------------
  // Halo refresh — the Phase 4 raison d'être.
  // ------------------------------------------------------------------
 
  void refresh_halo()
  {
    if constexpr (std::is_same_v<Boundary, OpenBoundary>)
      {
        for_each_halo_cell([&](const coord_type &u) { store_set(u, state_type{}); });
      }
    else if constexpr (ca_lattice_detail::is_constant_boundary_v<Boundary>)
      {
        const state_type v = static_cast<state_type>(Boundary::value);
        for_each_halo_cell([&](const coord_type &u) { store_set(u, v); });
      }
    else if constexpr (std::is_same_v<Boundary, NeumannBoundary>)
      {
        for_each_halo_cell([&](const coord_type &u)
        {
          coord_type src = u;
          for (std::size_t d = 0; d < rank; ++d)
            src[d] = static_cast<ca_index_t>(
              ca_ghost_detail::clamp_index(u[d], user_ext_[d]));
          store_set(u, store_get(src));
        });
      }
    else if constexpr (std::is_same_v<Boundary, ToroidalBoundary>)
      {
        for_each_halo_cell([&](const coord_type &u)
        {
          coord_type src = u;
          for (std::size_t d = 0; d < rank; ++d)
            src[d] = static_cast<ca_index_t>(
              ca_ghost_detail::positive_mod(u[d], user_ext_[d]));
          store_set(u, store_get(src));
        });
      }
    else if constexpr (std::is_same_v<Boundary, ReflectiveBoundary>)
      {
        for_each_halo_cell([&](const coord_type &u)
        {
          coord_type src = u;
          for (std::size_t d = 0; d < rank; ++d)
            src[d] = static_cast<ca_index_t>(
              ca_ghost_detail::reflect_index(u[d], user_ext_[d]));
          store_set(u, store_get(src));
        });
      }
    else
      {
        static_assert(std::is_same_v<Boundary, OpenBoundary>
                        or std::is_same_v<Boundary, ToroidalBoundary>
                        or std::is_same_v<Boundary, ReflectiveBoundary>
                        or std::is_same_v<Boundary, NeumannBoundary>
                        or ca_lattice_detail::is_constant_boundary_v<Boundary>,
                      "Ghost_Lattice: unsupported Boundary policy");
      }
  }
 
  // ------------------------------------------------------------------
  // Utilities shared with `Lattice`.
  // ------------------------------------------------------------------
 
  [[nodiscard]] const Storage &storage() const noexcept { return store_; }
 
  [[nodiscard]] Storage &storage() noexcept { return store_; }
 
  void fill(const state_type &value)
  {
    coord_type u{};
    auto go = [&](auto &&self_fn, std::size_t d) -> void
    {
      if (d == rank)
        {
          store_set(u, value);
          return;
        }
      const ca_index_t n = static_cast<ca_index_t>(user_ext_[d]);
      for (ca_index_t i = 0; i < n; ++i)
        {
          u[d] = i;
          self_fn(self_fn, d + 1);
        }
    };
    go(go, 0);
  }
 
  void swap(Ghost_Lattice &other) noexcept(noexcept(store_.swap(other.store_)))
  {
    using std::swap;
    swap(user_ext_, other.user_ext_);
    swap(store_ext_, other.store_ext_);
    store_.swap(other.store_);
  }
};
 
template <typename Storage, typename Boundary, std::size_t Halo>
inline void swap(Ghost_Lattice<Storage, Boundary, Halo> &a,
                 Ghost_Lattice<Storage, Boundary, Halo> &b) noexcept(noexcept(a.swap(b)))
{
  a.swap(b);
}
 
}  // namespace CA
}  // namespace Aleph
 
#endif  // TPL_CA_GHOST_LATTICE_H