indicators.hpp

Go to the documentation of this file.

// core — solution-adaptive refinement indicators for a BlockOctree.
//
// Geometry-based refinement lives in refine.hpp (refineToSdf). This header adds
// *solution*-based error indicators: per-leaf scalars derived from a leaf field
// that flag where the mesh should be refined (large indicator) or may be coarsened
// (small indicator). They drive the dynamic-AMR cycle in adapt.hpp.
//
// The headline criterion is the Löhner (1987) indicator — a dimensionless,
// normalized second difference that is the standard robust feature detector for
// AMR (it tracks curvature / steep fronts and self-normalizes so a single
// threshold works across scales). A raw second-derivative magnitude is also
// provided for when an un-normalized curvature measure is wanted.
//
// Neighbour values come from BlockOctree::faceNeighbor (the covering same/coarser/
// corner-finer neighbour under 2:1 balance); a face with no neighbour (block /
// domain edge) contributes nothing on that axis.
//
// Header-only, guarded by PECLET_CORE_HAVE_MORTON.
#ifndef PECLET_CORE_AMR_INDICATORS_HPP
#define PECLET_CORE_AMR_INDICATORS_HPP
 
#ifdef PECLET_CORE_HAVE_MORTON
 
#include <cmath>
#include <vector>
 
#include "peclet/core/amr/block_octree.hpp"
#include "peclet/core/common/types.hpp"
 
namespace peclet::core::amr {
 
template <int Dim, unsigned Bits>
std::vector<double> lohnerIndicator(const BlockOctree<Dim, Bits>& t, const std::vector<double>& u,
                                    double eps = 0.01) {
  const Index n = t.numLeaves();
  std::vector<double> e(static_cast<std::size_t>(n), 0.0);
  for (Index i = 0; i < n; ++i) {
    const double ui = u[static_cast<std::size_t>(i)];
    double num2 = 0.0, den2 = 0.0;
    for (int axis = 0; axis < Dim; ++axis) {
      const Index jp = t.faceNeighbor(i, axis, +1);
      const Index jm = t.faceNeighbor(i, axis, -1);
      if (jp < 0 || jm < 0)
        continue;
      const double up = u[static_cast<std::size_t>(jp)];
      const double um = u[static_cast<std::size_t>(jm)];
      const double d2 = up - 2.0 * ui + um;
      const double nrm = std::fabs(up - ui) + std::fabs(ui - um) +
                         eps * (std::fabs(up) + 2.0 * std::fabs(ui) + std::fabs(um));
      num2 += d2 * d2;
      den2 += nrm * nrm;
    }
    e[static_cast<std::size_t>(i)] = (den2 > 0.0) ? std::sqrt(num2 / den2) : 0.0;
  }
  return e;
}
 
template <int Dim, unsigned Bits>
std::vector<double> secondDiffIndicator(const BlockOctree<Dim, Bits>& t,
                                        const std::vector<double>& u) {
  const Index n = t.numLeaves();
  std::vector<double> e(static_cast<std::size_t>(n), 0.0);
  for (Index i = 0; i < n; ++i) {
    const double ui = u[static_cast<std::size_t>(i)];
    double s = 0.0;
    for (int axis = 0; axis < Dim; ++axis) {
      const Index jp = t.faceNeighbor(i, axis, +1);
      const Index jm = t.faceNeighbor(i, axis, -1);
      if (jp < 0 || jm < 0)
        continue;
      const double d2 =
          u[static_cast<std::size_t>(jp)] - 2.0 * ui + u[static_cast<std::size_t>(jm)];
      s += d2 * d2;
    }
    e[static_cast<std::size_t>(i)] = std::sqrt(s);
  }
  return e;
}
 
}  // namespace peclet::core::amr
 
#endif  // PECLET_CORE_HAVE_MORTON
#endif  // PECLET_CORE_AMR_INDICATORS_HPP