refine.hpp

Go to the documentation of this file.

// core — refinement criteria for a BlockOctree.
//
// Topology-mutation helpers that drive BlockOctree::refineIf from a geometric
// criterion. The headline one is SDF-driven: refine leaves the solid surface
// passes through (plus a band) down to a target level, leaving the interior /
// far field coarse — the usual AMR pattern for the suite's cut-cell IBM. Reuses
// the shared peclet::core::geom SDF (any callable returning a signed distance at a world
// point) and the AmrGeometry world mapping.
//
#ifndef PECLET_AMR_REFINE_HPP
#define PECLET_AMR_REFINE_HPP
 
#include <algorithm>
#include <cmath>
#include <utility>
#include <vector>
 
#include "peclet/amr/block_octree.hpp"
#include "peclet/amr/common.hpp"
#include "peclet/amr/leaf_field.hpp"
#include "peclet/core/common/types.hpp"
 
namespace peclet::amr {
 
template <int Dim, unsigned Bits, class SdfFn>
Index refineToSdf(BlockOctree<Dim, Bits>& t, const AmrGeometry<Dim>& geo, SdfFn&& sdf,
                  unsigned targetLevel, Real band = 1.0, bool balance = true) {
  using Code = typename BlockOctree<Dim, Bits>::Code;
  Index total = 0;
  for (;;) {
    std::vector<Code> toRefine;
    for (Index i = 0; i < t.numLeaves(); ++i) {
      const unsigned L = t.level(i);
      if (L <= targetLevel)
        continue;
      auto b = t.bounds(i);
      Vec<Dim> c = geo.center(b);
      // Phase 3 (`docs/amr_anisotropic.md` §6): a BOX leaf's half-diagonal is
      // 0.5*sqrt(sum_d w_d^2) and the band margin is a length, so it takes the COARSEST axis —
      // the conservative reading, which refines no less than the cubic rule did. Both reduce to
      // `0.5*sqrt(Dim)*w` and `band*h0` exactly when the three spacings are equal.
      const Vec<Dim> w = geo.leafSize(L);
      Real halfDiag2 = 0.0;
      for (int d = 0; d < Dim; ++d)
        halfDiag2 += w[d] * w[d];
      const Real halfDiag = 0.5 * std::sqrt(halfDiag2);
      if (std::fabs(static_cast<Real>(sdf(c))) <= halfDiag + band * geo.hMax())
        toRefine.push_back(t.code(i));
    }
    if (toRefine.empty())
      break;
    std::sort(toRefine.begin(), toRefine.end());
    toRefine.erase(std::unique(toRefine.begin(), toRefine.end()), toRefine.end());
    total += t.refineIf(
        [&](Code c, unsigned) { return std::binary_search(toRefine.begin(), toRefine.end(), c); });
  }
  if (balance)
    total += t.balance2to1();
  return total;
}
 
template <int Dim, unsigned Bits, class SdfFn, class TargetFn>
Index refineToSdfGraded(BlockOctree<Dim, Bits>& t, const AmrGeometry<Dim>& geo, SdfFn&& sdf,
                        TargetFn&& targetFn, Real band = 2.0, bool balance = true) {
  using Code = typename BlockOctree<Dim, Bits>::Code;
  const Real halfDiagFactor = 0.5 * std::sqrt(static_cast<Real>(Dim));
  Index total = 0;
  for (;;) {
    std::vector<Code> toRefine;
    for (Index i = 0; i < t.numLeaves(); ++i) {
      const unsigned L = t.level(i);
      if (L == 0)
        continue;
      auto b = t.bounds(i);
      Vec<Dim> c = geo.center(b);
      // Phase 3: the half-diagonal of a BOX leaf, and a band measured in cells of the level being
      // created — the coarsest axis's width, the conservative reading (§6 of the design note).
      const Vec<Dim> w = geo.leafSize(L);
      Real halfDiag2 = 0.0;
      Real wMax = w[0];
      for (int d = 0; d < Dim; ++d) {
        halfDiag2 += w[d] * w[d];
        wMax = w[d] > wMax ? w[d] : wMax;
      }
      const Real halfDiag = 0.5 * std::sqrt(halfDiag2);
      // Band predicate first: targetFn is the expensive one (a gap/medial-axis field, often a
      // Python callable) and only band cells can ever be refined.
      if (std::fabs(static_cast<Real>(sdf(c))) > halfDiag + band * 0.5 * wMax)
        continue;
      if (L > targetFn(c))
        toRefine.push_back(t.code(i));
    }
    if (toRefine.empty())
      break;
    std::sort(toRefine.begin(), toRefine.end());
    toRefine.erase(std::unique(toRefine.begin(), toRefine.end()), toRefine.end());
    total += t.refineIf(
        [&](Code c, unsigned) { return std::binary_search(toRefine.begin(), toRefine.end(), c); });
  }
  if (balance)
    total += t.balance2to1();
  return total;
}
 
template <int Dim, class GapFn>
auto gapFloorTarget(GapFn&& gapFn, Real h0, unsigned coarsestLevel, Real n = 4.0) {
  return [gapFn = std::forward<GapFn>(gapFn), h0, coarsestLevel, n](const Vec<Dim>& p) -> unsigned {
    const Real g = static_cast<Real>(gapFn(p));
    unsigned L = 0;
    while (L < coarsestLevel && h0 * static_cast<Real>(Index(1) << (L + 1)) * n <= g)
      ++L;
    return L;
  };
}
 
}  // namespace peclet::amr
 
#endif  // PECLET_AMR_REFINE_HPP