leaf_field.hpp

Go to the documentation of this file.

// core — leaf-indexed fields and world geometry for a BlockOctree.
//
// A BlockOctree addresses leaves in integer *fine units* relative to the block
// origin. AmrGeometry maps those to world coordinates (origin + h0 * fine), and
// LeafField<T> is a value-per-leaf array in the octree's leaf (Z-order) slot
// order — the simulation field the device kernels and the VTU writer operate on.
// The on-device counterpart is simply a peclet::core::View<T> of length numLeaves(); this
// host container is the I/O / setup form.
//
#ifndef PECLET_AMR_LEAF_FIELD_HPP
#define PECLET_AMR_LEAF_FIELD_HPP
 
#include <array>
#include <cstddef>
#include <vector>
 
#include "peclet/amr/block_octree.hpp"
#include "peclet/amr/common.hpp"
#include "peclet/core/common/types.hpp"
 
namespace peclet::amr {
 
namespace detail {
template <int Dim>
inline Vec<Dim> filledVec(Real v) {
  Vec<Dim> x{};
  for (int d = 0; d < Dim; ++d)
    x[d] = v;
  return x;
}
}  // namespace detail
 
template <int Dim>
struct AmrGeometry {
  Vec<Dim> origin{};
  Vec<Dim> h0 = detail::filledVec<Dim>(1.0);
 
  static AmrGeometry isotropicAt(const Vec<Dim>& o, Real h) {
    AmrGeometry g;
    g.origin = o;
    g.h0 = detail::filledVec<Dim>(h);
    return g;
  }
  void setIsotropic(Real h) { h0 = detail::filledVec<Dim>(h); }
 
  Real hMin() const {
    Real m = h0[0];
    for (int d = 1; d < Dim; ++d)
      m = h0[d] < m ? h0[d] : m;
    return m;
  }
  Real hMax() const {
    Real m = h0[0];
    for (int d = 1; d < Dim; ++d)
      m = h0[d] > m ? h0[d] : m;
    return m;
  }
  bool isotropic() const {
    for (int d = 1; d < Dim; ++d)
      if (h0[d] != h0[0])
        return false;
    return true;
  }
  Real cellVolume() const {
    Real v = h0[0];
    for (int d = 1; d < Dim; ++d)
      v *= h0[d];
    return v;
  }
 
  Real leafSize(unsigned level, int d) const {
    return h0[d] * static_cast<Real>(Index(1) << level);
  }
  Vec<Dim> leafSize(unsigned level) const {
    Vec<Dim> s{};
    const Real f = static_cast<Real>(Index(1) << level);
    for (int d = 0; d < Dim; ++d)
      s[d] = h0[d] * f;
    return s;
  }
 
  template <class Coord>
  Vec<Dim> lowerCorner(const std::array<Coord, Dim>& lo) const {
    Vec<Dim> p{};
    for (int d = 0; d < Dim; ++d)
      p[d] = origin[d] + static_cast<Real>(lo[d]) * h0[d];
    return p;
  }
 
  template <class Coord>
  Vec<Dim> center(const std::array<std::array<Coord, Dim>, 2>& b) const {
    Vec<Dim> p{};
    for (int d = 0; d < Dim; ++d)
      p[d] =
          origin[d] +
          (static_cast<Real>(b[0][d]) + 0.5 * (static_cast<Real>(b[1][d] - b[0][d]) + 1.0)) * h0[d];
    return p;
  }
};
 
template <class T>
struct LeafField {
  std::vector<T> values;
 
  LeafField() = default;
  explicit LeafField(Index n, T init = T{}) : values(static_cast<std::size_t>(n), init) {}
 
  template <int Dim, unsigned Bits>
  explicit LeafField(const BlockOctree<Dim, Bits>& t, T init = T{})
      : values(static_cast<std::size_t>(t.numLeaves()), init) {}
 
  Index size() const { return static_cast<Index>(values.size()); }
  T& operator[](Index i) { return values[static_cast<std::size_t>(i)]; }
  const T& operator[](Index i) const { return values[static_cast<std::size_t>(i)]; }
  T* data() { return values.data(); }
  const T* data() const { return values.data(); }
};
 
}  // namespace peclet::amr
 
#endif  // PECLET_AMR_LEAF_FIELD_HPP