csys.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2025 VTT Technical Research Centre of Finland Ltd
// SPDX-License-Identifier: AGPL-3.0-or-later
 
#pragma once
 
#include "types.hpp"
 
namespace pfc {
namespace csys {
 
using pfc::types::Bool3;
using pfc::types::Int3;
using pfc::types::Real3;
 
template <typename Tag> struct CoordinateSystem;
 
template <typename CoordTag>
struct CoordinateSystemDefaults; // intentionally left undefined
 
// Coordinate system tags.
// struct LineTag {};
// struct PlaneTag {};
// struct PolarTag {};
// struct CylindricalTag {};
// struct SphericalTag {};
// struct Polar2DTag {};
// struct Toroidal2DTag {};
// struct LogPolar3DTag {};
 
struct CartesianTag {};
 
template <> struct CoordinateSystemDefaults<CartesianTag> {
  static constexpr Real3 offset = {0.0, 0.0, 0.0};
  static constexpr Real3 spacing = {1.0, 1.0, 1.0};
  static constexpr Bool3 periodic = {true, true, true};
  std::size_t dimensions = 3; 
};
 
template <> struct CoordinateSystem<CartesianTag> {
  const Real3 m_offset;   
  const Real3 m_spacing;  
  const Bool3 m_periodic; 
 
  CoordinateSystem(
      const Real3 &offset = CoordinateSystemDefaults<CartesianTag>::offset,
      const Real3 &spacing = CoordinateSystemDefaults<CartesianTag>::spacing,
      const Bool3 &periodic = CoordinateSystemDefaults<CartesianTag>::periodic)
      : m_offset(offset), m_spacing(spacing), m_periodic(periodic) {
    for (std::size_t i = 0; i < 3; ++i) {
      if (spacing[i] <= 0.0) {
        throw std::invalid_argument("Spacing must be positive.");
      }
    };
  };
 
  bool operator==(const CoordinateSystem<CartesianTag> &other) const {
    return m_offset == other.m_offset && m_spacing == other.m_spacing &&
           m_periodic == other.m_periodic;
  };
 
  bool operator!=(const CoordinateSystem<CartesianTag> &other) const {
    return !(*this == other);
  };
};
 
using CartesianCS = CoordinateSystem<CartesianTag>;
 
inline const Real3 &get_offset(const CartesianCS &cs) noexcept {
  return cs.m_offset;
};
 
inline double get_offset(const CartesianCS &cs, int i) { return cs.m_offset.at(i); }
 
inline const Real3 &get_spacing(const CartesianCS &cs) noexcept {
  return cs.m_spacing;
};
 
inline double get_spacing(const CartesianCS &cs, int i) {
  return cs.m_spacing.at(i);
}
 
inline const Bool3 &get_periodic(const CartesianCS &cs) noexcept {
  return cs.m_periodic;
};
 
inline bool is_periodic(const CartesianCS &cs, int i) { return cs.m_periodic.at(i); }
 
inline const Real3 to_coords(const CartesianCS &cs, const Int3 &idx) noexcept {
  Real3 xyz;
  const auto &offset = get_offset(cs);
  const auto &spacing = get_spacing(cs);
  for (int i = 0; i < 3; ++i) {
    xyz[i] = offset[i] + idx[i] * spacing[i];
  }
  return xyz;
}
 
inline const Int3 to_index(const CartesianCS &cs, const Real3 &xyz) noexcept {
  Int3 idx;
  const auto &offset = get_offset(cs);
  const auto &spacing = get_spacing(cs);
  for (int i = 0; i < 3; ++i) {
    idx[i] = static_cast<int>((xyz[i] - offset[i]) / spacing[i]);
  }
  return idx;
}
 
} // namespace csys
} // namespace pfc