TopoOffsetTetMesh.h

#pragma once
 
#include <wmtk/MeshData.hpp>
 
#include <array>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <functional>
#include <future>
#include <limits>
#include <map>
#include <mutex>
#include <optional>
#include <set>
#include <string>
#include <type_traits>
 
#include <wmtk/TetMesh.h>
#include <wmtk/TetOptimizerMesh.h>
#include <wmtk/envelope/Envelope.hpp>
#include <wmtk/optimization/EnergySum.hpp>
#include <wmtk/optimization/solver.hpp>
#include <wmtk/simplex/Simplex.hpp>
#include <wmtk/threading/dynamic_parallel_for.hpp>
#include <wmtk/threading/enumerable_thread_specific.hpp>
#include "ExactBandIntegral.hpp"
#include "OffsetPotential.hpp"
#include "Parameters.h"
 
// clang-format off
#include <wmtk/utils/DisableWarnings.hpp>
#include <wmtk/utils/EnableWarnings.hpp>
// clang-format on
 
using CellTag = std::set<int64_t>;
 
namespace wmtk::utils {
class WindingNumberHierarchy;
}
 
 
namespace wmtk::components::topological_offset {
 
const int64_t TEMP_OFFSET_TET_TAG = -1;
const CellTag TEMP_OFFSET_TET_TAG_SET{TEMP_OFFSET_TET_TAG};
 
 
// for all attributes:
// label: 0=default, 1=input, 2=offset
 
inline int closed_label(const int a, const int b)
{
    if (a == 1 || b == 1) return 1;
    if (a == 2 || b == 2) return 2;
    return 0;
}
 
class VertexExtra
{
public:
    int label = 0;
    bool m_is_on_input = false; // on the input complex
    bool m_is_on_offset = false; // on the offset surface itself
    bool m_is_on_region = false; // on some OTHER tag region's boundary
    int64_t m_corr_input_vid = -1;
 
    uint32_t m_born_epoch = 0;
 
    double m_front_change = 0.;
    int m_front_change_group = -1;
};
 
 
class EdgeAttributes
{
public:
    int label = 0; // label: 0=default, 1=input, 2=offset
};
 
 
class FaceExtra
{
public:
    int label = 0; // label: 0=default, 1=input, 2=offset
    bool on_sheet = false;
 
    void merge(const FaceExtra& o)
    {
        label = closed_label(label, o.label);
        on_sheet = on_sheet || o.on_sheet;
    }
};
 
 
class TetAttributes
{
public:
    int label = 0; // label: 0=default, 1=input, 2=offset
    CellTag tag;
    double m_quality = 0; // AMIPS energy, kept up to date by smoothing
    bool rest_valid = false;
    std::array<Vector3d, 4> rest_pos;
};
 
 
class TopoOffsetTetMesh : public wmtk::TetOptimizerMesh
{
public: // mode for splitting in marching tets
    enum class EdgeSplitMode {
        Midpoint = 0, // construction: simplicial embedding AND marching_tets
        SphereTrace = 1, // marching_tets (construction_mode): sphere tracing along the edge to
                         // d(x) = m_construction_distance, midpoint when the trace leaves the
                         // edge
        Optimization = 5 // the optimization phase; the shared engine places the vertex
    };
 
public:
    std::array<size_t, 4> m_init_counts = {{0, 0, 0, 0}};
    size_t m_tags_count;
    int64_t m_sheet_tag = -1;
    std::vector<std::string> m_face_tag_names;
    std::map<int64_t, int> m_face_tag_bit;
    bool is_sheet_tag(int64_t tag) const { return tag == m_sheet_tag || m_face_tag_bit.count(tag); }
    MatrixXd m_sheet_V;
    MatrixXi m_sheet_F;
    void classify_sheet_faces();
    double input_complex_distance(const Vector3d& p) const;
    std::shared_ptr<const wmtk::utils::WindingNumberHierarchy> m_input_solid_winding;
 
    std::shared_ptr<OffsetPotential3D> m_offset_potential;
 
    int64_t front_face_band_cell(size_t a, size_t b, size_t c) const;
    double front_vertex_relative_residual(size_t vid) const;
    double front_vertex_residual_length(size_t vid) const;
    Vector3d front_vertex_field_gradient(size_t vid) const;
 
    std::shared_ptr<Envelope> m_input_complex_envelope;
    double m_input_complex_pad = 0.;
 
    int m_n_regions = 0; 
    std::vector<std::shared_ptr<OffsetPotential3D>> m_region_potentials; 
    std::vector<std::shared_ptr<Envelope>> m_region_envelopes;
    std::vector<int64_t> m_phi_vert_region; 
    std::vector<int64_t> m_phi_seg_region; 
    std::vector<int64_t> m_phi_face_region; 
    std::vector<int64_t> m_phi_point_region; 
    std::vector<int> m_cell_region; 
    std::vector<int> m_vertex_region; 
    void init_region_potentials(double delta, double effective_factor);
    void assign_band_regions(bool log = true);
    void log_front_profile(size_t vid);
    int vertex_region(const size_t vid) const
    {
        return vid < m_vertex_region.size() ? m_vertex_region[vid] : -1;
    }
    int edge_region(const size_t va, const size_t vb) const
    {
        const int a = vertex_region(va), b = vertex_region(vb);
        return (a >= 0 && a == b) ? a : -1;
    }
    const OffsetPotential3D& potential_for_region(const int region) const
    {
        return (region >= 0 && size_t(region) < m_region_potentials.size())
                   ? *m_region_potentials[size_t(region)]
                   : *m_offset_potential;
    }
    const OffsetPotential3D& potential_for(const size_t vid) const
    {
        return potential_for_region(vertex_region(vid));
    }
    std::shared_ptr<const OffsetPotential3D> potential_ptr_for(const size_t vid) const
    {
        const int r = vertex_region(vid);
        return (r >= 0 && size_t(r) < m_region_potentials.size()) ? m_region_potentials[size_t(r)]
                                                                  : m_offset_potential;
    }
    const OffsetPotential3D& potential_for_edge(const size_t va, const size_t vb) const
    {
        return potential_for_region(edge_region(va, vb));
    }
    const OffsetPotential3D& potential_for_face(const Tuple& f) const;
 
    bool face_is_held(const Tuple& f) const;
    bool vertex_is_held(size_t vid) const;
    bool face_is_complex_boundary(const Tuple& f) const;
    void build_envelopes();
    void build_final_envelopes();
    void release_final_envelopes();
    double tetwild_pass(const char* what, int max_its, bool coarsen = false);
    void optimize_construction_pass();
    double m_construction_envelope_eps = 0.;
    struct MeshSnapshot
    {
        size_t n_vertices = 0;
        std::vector<std::array<size_t, 4>> tets;
        std::vector<VertexAttributes> vertex;
        std::vector<VertexExtra> vertex_extra;
        std::vector<FaceAttributes> face;
        std::vector<FaceExtra> face_extra;
        std::vector<EdgeAttributes> edge;
        std::vector<TetAttributes> tet;
    };
    MeshSnapshot take_snapshot();
    void restore_snapshot(const MeshSnapshot& s);
    std::vector<Eigen::Vector3d> m_held_verts;
    std::vector<Eigen::Vector3i> m_held_tris;
    std::shared_ptr<Envelope> m_loop_envelope;
 
    bool m_tetwild_pass = false;
 
    enum class OpRules { Loop, TetWild };
    OpRules m_op_rules = OpRules::Loop;
    bool tetwild_rules() const { return m_op_rules == OpRules::TetWild; }
 
    std::shared_ptr<Envelope> m_offset_envelope;
    void build_offset_envelope();
 
    void check_no_vertex_on_both_surfaces(const char* when) const;
 
    bool optimize_offset_loop(
        bool final_stage = true,
        bool refine = true,
        const std::string& label = std::string(),
        int max_turns = 0,
        bool fold_free = false);
    void loop_group(const std::array<int, 4>& ops, const std::string& where);
    void front_pocket_split(const std::string& where);
    void label_cleanup(const std::string& where);
    double label_flip_surface_change(size_t tid) const;
    static constexpr int kEnergyDecreasePasses = 100;
    static constexpr int kPocketMaxRounds = 10;
    bool m_pocket_splitting = false;
    std::vector<std::array<size_t, 2>> m_settle_split_edges;
 
    bool m_init_optimize = false;
    std::string m_frame_prefix;
    void linear_growth();
    struct SmartGrowthSchedule
    {
        double base; 
        double delta; 
        int max_fails; 
        int j = 0; 
        int level = 1; 
        double target() const { return j == 0 ? delta : base + (delta - base) / std::ldexp(1., j); }
        enum class Next { Done, Advance, Retry, GiveUp };
        Next report(const bool converged)
        {
            if (converged) {
                if (j == 0) return Next::Done;
                base = target();
                j = 0;
                ++level;
                return Next::Advance;
            }
            ++j;
            return j >= max_fails ? Next::GiveUp : Next::Retry;
        }
    };
    void smart_growth();
    static double growth_stage_target(double d0, double delta, int k, int n);
 
    double front_gradient_linf();
    double m_front_gradient_reference = 0.;
 
    EdgeSplitMode m_edge_split_mode = EdgeSplitMode::Midpoint;
 
    // tag name maps
    std::map<std::string, int64_t> m_tag_name_to_id;
    std::map<int64_t, std::string> m_tag_id_to_name;
    CellTag m_offset_output_tag_ids;
 
    // if in 'singlebody' mode
    bool m_singlebody = false;
    int64_t m_single_tag;
 
    // just for retaining in output. dont actually use
    bool m_has_envelope = false;
    MatrixXd m_V_envelope;
    MatrixXi m_F_envelope;
    // m_envelope itself lives on the base, which is what checks tracked-surface triangles
    // against it.
    double m_envelope_eps = -1;
 
    static constexpr int INPUT_SURFACE_CLASS = 0;
    static constexpr int OFFSET_SURFACE_CLASS = 1;
 
    Parameters& m_offset_params;
 
    using VertexExtraCol = wmtk::AttributeCollection<VertexExtra>;
    using EdgeAttCol = wmtk::AttributeCollection<EdgeAttributes>;
    using FaceExtraCol = wmtk::AttributeCollection<FaceExtra>;
    using TetAttCol = wmtk::AttributeCollection<TetAttributes>;
    // m_vertex_attribute and m_face_attribute are the base's; these are registered alongside
    // them in its attribute groups.
    VertexExtraCol m_vertex_extra;
    FaceExtraCol m_face_extra;
    EdgeAttCol m_edge_attribute;
    TetAttCol m_tet_attribute;
 
    using FaceExtraMap = std::map<std::array<size_t, 3>, FaceExtra>;
    void face_extra_capture(const std::vector<size_t>& tids, FaceExtraMap& out) const;
    void face_extra_assign(const std::vector<size_t>& tids, const FaceExtraMap& captured);
    wmtk::threading::enumerable_thread_specific<FaceExtraMap> m_face_extra_split, m_face_extra_swap;
    wmtk::threading::enumerable_thread_specific<
        std::vector<std::pair<FaceExtra, std::array<size_t, 2>>>>
        m_face_extra_collapse;
    using EdgeLabelMap = std::map<std::array<size_t, 2>, int>;
    void edge_label_capture(
        const std::vector<size_t>& tids,
        EdgeLabelMap& out,
        size_t from = size_t(-1),
        size_t to = size_t(-1)) const;
    void edge_label_assign(const std::vector<size_t>& tids, const EdgeLabelMap& captured);
    wmtk::threading::enumerable_thread_specific<EdgeLabelMap> m_edge_label_split, m_edge_label_swap,
        m_edge_label_collapse;
    void close_labels(const std::vector<size_t>& tids);
 
    TopoOffsetTetMesh(Parameters& _m_offset_params, int _num_threads = 0)
        : wmtk::TetOptimizerMesh(_m_offset_params, nullptr)
        , m_offset_params(_m_offset_params)
    {
        NUM_THREADS = _num_threads;
        // The base owns the vertex, edge and face slots; register the offset's own data with
        // its groups so it is resized, protected and rolled back with them.
        m_vertex_attr_group.add(&m_vertex_extra);
        m_face_attr_group.add(&m_face_extra);
        add_edge_attributes(&m_edge_attribute);
        p_tet_attrs = &m_tet_attribute;
 
        m_collapse_check_link_condition = false;
        m_collapse_check_manifold = false;
 
        // As in 2D. The per-vertex Newton solver logs a line per smoothing attempt at info level,
        // which is one line per vertex per pass and buries the run's own output.
        optimization::deactivate_opt_logger();
    }
 
    ~TopoOffsetTetMesh() override = default;
 
 
    double cell_quality(const size_t tid) const override { return m_tet_attribute[tid].m_quality; }
    void set_cell_quality(const size_t tid, const double q) override
    {
        m_tet_attribute[tid].m_quality = q;
    }
 
    double tet_energy(size_t tid, bool band_term = true) const;
    bool band_term_cancels(const std::vector<size_t>& tids) const
    {
        // "exact_swept": every operation's sums leave the band term out, and its rule adds the
        // band term's exact change (EXACT SWEPT CHANGE below)
        if (swept_rule()) return true;
        if (m_offset_params.band_volume_rule != "exact" &&
            m_offset_params.band_volume_rule != "front_flux")
            return false;
        std::vector<size_t> vids;
        for (const size_t tid : tids) {
            if (!cell_is_offset_band(tid)) return false;
            for (const size_t v : oriented_tet_vids(tid)) vids.push_back(v);
        }
        wmtk::vector_unique(vids);
        std::vector<Vector3d> pts;
        pts.reserve(vids.size());
        for (const size_t v : vids) pts.push_back(m_vertex_attribute[v].m_posf);
        const int n = int(pts.size());
        if (m_band_exact->front_flux_flag(pts.data(), n)) return false;
        return m_band_exact->distance_bound(pts.data(), n) <= m_band_exact->front_flux_reach();
    }
    bool swept_rule() const { return m_offset_params.band_volume_rule == "exact_swept"; }
    double band_swept_change(
        size_t vid,
        const Vector3d& x0,
        const Vector3d& x1,
        ExactBandIntegral::Stats* st = nullptr) const;
    void band_swept_faces(
        size_t vid,
        const Vector3d& x0,
        std::vector<std::array<Vector3d, 2>>& faces) const;
    double band_swept_sum(
        const Vector3d& x0,
        const Vector3d& x1,
        const std::vector<std::array<Vector3d, 2>>& faces,
        ExactBandIntegral::Stats* st = nullptr) const;
    double band_chain_change(
        const std::vector<std::array<size_t, 4>>& old_band,
        const std::vector<std::array<size_t, 4>>& new_band,
        ExactBandIntegral::Stats* st = nullptr) const;
    double swap_swept_change(const std::vector<std::array<size_t, 4>>& new_band, bool flip);
    static constexpr int kSweptHalvings = 4;
    struct SweptCounters
    {
        std::atomic<size_t> solves{0}, refused{0};
        std::array<std::atomic<size_t>, kSweptHalvings + 1> accepted{};
    };
    bool swept_accept(
        size_t vid,
        const std::vector<size_t>& ring,
        const Vector3d& x0,
        double amips_before,
        SweptCounters& c);
    mutable SweptCounters m_swept_front, m_swept_nonfront;
    mutable BandVolumeEnergy3D::SweptDerivCounts m_swept_derivs;
    mutable wmtk::threading::enumerable_thread_specific<std::pair<size_t, Vector3d>> m_swept_stored;
    double smoothing_ring_energy(const std::vector<Tuple>& ring, const bool stored) const override
    {
        // TetWild's rules (m_op_rules): the base's max over the ring.
        if (tetwild_rules()) return TetOptimizerMesh::smoothing_ring_energy(ring, stored);
        double s = 0.;
        if (swept_rule()) {
            // "exact_swept": the ring's AMIPS terms plus the band term's change from the position
            // at the `stored` call (EXACT SWEPT CHANGE); the vertex is the one
            // smooth_nonfront_vertex() is smoothing, the same in both calls
            for (const Tuple& loc : ring) s += tet_energy(loc.tid(*this), false);
            auto& st = m_swept_stored.local();
            const size_t vid = st.first;
            const auto vs = oriented_tet_vids(ring.front().tid(*this));
            if (std::find(vs.begin(), vs.end(), vid) == vs.end()) {
                log_and_throw_error(
                    "exact_swept: the projected step's ring is not vertex {}'s",
                    vid);
            }
            if (stored) {
                st.second = m_vertex_attribute[vid].m_posf;
                return s;
            }
            return s + band_swept_change(vid, st.second, m_vertex_attribute[vid].m_posf);
        }
        for (const Tuple& loc : ring) s += tet_energy(loc.tid(*this));
        // EXPERIMENTAL_band_volume_rule "surface": the offset faces through the vertex move with
        // it (smooth_nonfront_vertex() names it).
        if (surface_energy()) {
            const size_t vid = m_smoothing_vid.local();
            if (vid != size_t(-1) && m_vertex_extra[vid].m_is_on_offset) s += surface_term_at(vid);
        }
        return s;
    }
    double vol_amips3(const std::array<size_t, 4>& vids, const double amips3) const
    {
        double tr = 0.;
        for (int i = 0; i < 4; ++i) {
            for (int j = i + 1; j < 4; ++j) {
                tr += (m_vertex_attribute[vids[size_t(i)]].m_posf -
                       m_vertex_attribute[vids[size_t(j)]].m_posf)
                          .squaredNorm();
            }
        }
        tr *= 0.5;
        return std::sqrt(2.) / 12. * tr * std::sqrt(tr) * std::sqrt(amips3);
    }
 
    double band_cell_term(const std::array<size_t, 4>& vids) const
    {
        if (m_n_regions > 1 && !m_offset_params.allow_mult_comp) {
            log_and_throw_error(
                "the energy's band term reads one field, and this input has {} regions",
                m_n_regions);
        }
        if (m_offset_params.band_volume_rule == "exact" ||
            m_offset_params.band_volume_rule == "front_flux" || swept_rule()) {
            // EXPERIMENTAL_band_volume_rule "exact" (and "front_flux", whose per-cell quantity
            // it is, see above): J_t / delta, J_t = int_t (d - delta) dV in closed form
            // (ExactBandIntegral), in place of the corner bound
            const std::array<Vector3d, 4> p = {
                {m_vertex_attribute[vids[0]].m_posf,
                 m_vertex_attribute[vids[1]].m_posf,
                 m_vertex_attribute[vids[2]].m_posf,
                 m_vertex_attribute[vids[3]].m_posf}};
            const double delta = m_offset_params.target_distance;
            return m_band_exact->value(p, delta) / delta;
        }
        int64_t best = -1;
        return band_cell_corner_bound(vids, best);
    }
    double band_cell_corner_bound(const std::array<size_t, 4>& vids, int64_t& best) const
    {
        const double delta = m_offset_params.target_distance;
        std::array<Vector3d, 4> p;
        for (int j = 0; j < 4; ++j) p[size_t(j)] = m_vertex_attribute[vids[size_t(j)]].m_posf;
        const double vol = std::abs((p[1] - p[0]).dot((p[2] - p[0]).cross(p[3] - p[0]))) / 6.;
        const double md =
            min_mean_distance(*m_input_complex_envelope, m_input_complex_pad, p.data(), 4, best);
        return vol * ((md - delta) / delta);
    }
    std::shared_ptr<const ExactBandIntegral> m_band_exact;
    double band_amips_weight() const
    {
        const double s = optimization_stop_metric();
        return m_offset_params.offset_amips_weight / (s * s * s);
    }
    double cell_energy(const std::array<size_t, 4>& vids, const double amips3, const bool band)
        const
    {
        if (amips3 >= MAX_ENERGY) return MAX_ENERGY;
        double a = band_amips_weight() * vol_amips3(vids, amips3);
        // EXPERIMENTAL_band_volume_rule "surface": no cell carries an offset term; it lives on the
        // front's faces (surface_face_term()), which every rule adds beside the cells' sums.
        if (band && !surface_energy()) a += band_cell_term(vids);
        if (m_repulsion_potential) a += repulsion_cell_term(vids);
        return a;
    }
    double tet_energy_sum(const std::vector<size_t>& tids, bool band_term = true) const;
    double candidate_energy(const std::array<size_t, 4>& vids);
    int candidate_label(const std::array<size_t, 4>& vids);
    mutable std::atomic<int> m_volume_reject_collapse{0}, m_volume_reject_swap{0};
    mutable wmtk::threading::enumerable_thread_specific<double> m_split_energy_before;
    mutable std::atomic<int> m_volume_reject_split{0};
 
    double global_energy() const;
    double offset_energy() const;
 
    bool surface_energy() const { return m_offset_params.band_volume_rule == "surface"; }
    double surface_term_weight() const;
    double surface_face_term(size_t a, size_t b, size_t c) const;
    double surface_term_at(size_t vid) const;
    double surface_term_total() const;
    double stencil_max_sq_error() const;
    double face_stencil_max_sq(const std::array<size_t, 3>& f) const;
    double surface_integral_raw() const;
    double settle_metric() const
    {
        if (!surface_energy()) return offset_energy();
        return m_offset_params.settle_test == "bar" ? surface_integral_raw()
                                                    : stencil_max_sq_error();
    }
    void live_offset_faces_at(size_t vid, std::vector<std::array<size_t, 3>>& out) const;
    mutable wmtk::threading::enumerable_thread_specific<size_t> m_smoothing_vid{size_t(-1)};
    double offset_term_weight() const
    {
        const double f = m_offset_params.front_conv_frac();
        return 1. / (f * f);
    }
 
    void set_vertex_position(const size_t vid, const Vector3d& p)
    {
        m_vertex_attribute[vid].m_posf = p;
        m_vertex_attribute[vid].set_pos_to_posf();
        m_vertex_attribute[vid].m_is_rounded = true;
    }
 
    bool face_is_offset(const size_t fid) const
    {
        return m_face_attribute[fid].m_is_surface_fs &&
               m_face_attribute[fid].m_surface_class == OFFSET_SURFACE_CLASS;
    }
    bool face_is_region(const size_t fid) const
    {
        return m_face_attribute[fid].m_is_surface_fs &&
               m_face_attribute[fid].m_surface_class != OFFSET_SURFACE_CLASS;
    }
 
    struct FaceSnapshot
    {
        FaceAttributes tags;
        FaceExtra extra;
    };
    FaceSnapshot face_snapshot(const size_t fid) const
    {
        return FaceSnapshot{m_face_attribute[fid], m_face_extra[fid]};
    }
    void restore_face(const size_t fid, const FaceSnapshot& s)
    {
        m_face_attribute[fid] = s.tags;
        m_face_extra[fid] = s.extra;
    }
 
    void label_offset_boundary();
 
    bool cell_in_region(const size_t tid) const
    {
        const int l = m_tet_attribute[tid].label;
        return l == 1 || l == 2;
    }
    bool cell_is_input_complex(const size_t tid) const { return m_tet_attribute[tid].label == 1; }
    bool cell_is_offset_band(const size_t tid) const { return m_tet_attribute[tid].label == 2; }
 
    // ------- the plastic medium (use_rest_pose) -------
 
    bool m_plastic_active = false;
    bool cell_is_plastic(const size_t tid) const
    {
        return m_plastic_active && !cell_in_region(tid);
    }
    bool plastic_rest(const size_t tid, Eigen::Matrix3d& R) const;
    void stamp_plastic_rests();
    void stamp_rest_cell(size_t tid);
    bool all_plastic(const std::vector<size_t>& tids) const;
    bool plastic_tie(const std::vector<size_t>& before, const std::vector<int>& after_labels) const;
 
    void optimize_offset(const std::filesystem::path& output_file);
 
    std::pair<double, double> compute_distance_deviation() const;
 
    mutable size_t m_worst_dist_vid = static_cast<size_t>(-1);
    void log_worst_dist_vertex() const;
 
    bool face_is_offset_surface_live(const Tuple& f) const;
    static constexpr double kCos150 = -0.86602540378443864676;
    bool offset_face_normal_out_of_band(const Tuple& f, Vector3d& n) const;
    size_t count_real_folds(
        std::vector<char>* corners = nullptr,
        std::vector<std::array<size_t, 3>>* faces = nullptr) const;
    bool edge_is_offset_surface_live(size_t a, size_t b) const;
    std::vector<std::array<size_t, 3>> offset_surface_faces() const;
    std::vector<Tuple> offset_surface_faces_live_at(size_t vid) const;
    bool vertex_has_live_offset_face(size_t vid) const;
    template <typename F>
    bool for_each_live_offset_face_at(size_t vid, F&& fn) const;
    void refresh_offset_membership(size_t vid);
    std::pair<size_t, size_t> offset_membership_mismatches() const;
    void check_offset_membership(const char* when) const;
 
    mutable std::atomic<long long> m_offset_face_lookup_misses{0};
    mutable std::atomic<long long> m_offset_face_invalid_tuple{0};
    void report_offset_face_lookup_misses(const char* when) const;
 
    double front_face_field_angle(size_t a, size_t b, size_t c) const;
    bool front_face_out_normal(size_t a, size_t b, size_t c, Vector3d& n) const;
    bool front_face_out_normal_in(size_t a, size_t b, size_t c, size_t band_tid, Vector3d& n) const;
    double front_face_field_angle_in(size_t a, size_t b, size_t c, size_t band_tid) const;
    double front_face_angle_to_field(size_t a, size_t b, size_t c, const Vector3d& n) const;
    static double band_side_dihedral(const Vector3d& n1, const Vector3d& n2, const Vector3d& t)
    {
        return 180. - std::atan2(n1.cross(n2).dot(t), n1.dot(n2)) * 180. / M_PI;
    }
    static bool dihedral_is_foldover(const double dihedral, const double angle)
    {
        return dihedral < angle || dihedral > 360. - angle;
    }
    bool front_face_is_foldover(size_t a, size_t b, size_t c) const
    {
        return front_face_field_angle(a, b, c) > m_offset_params.foldover_angle;
    }
    struct FoldCensus
    {
        size_t n_fold = 0, n_faces = 0, n_unmeasured = 0;
        double worst = 0.;
        Vector3d worst_at = Vector3d::Zero();
    };
    FoldCensus fold_census() const;
    bool fold_reject_active() const
    {
        return m_offset_params.growth_mode == "individual_op_reject" && !m_tetwild_pass &&
               !m_repulsion_potential;
    }
    void folded_front_faces_at(size_t vid, std::set<std::array<size_t, 3>>& out) const;
    bool face_folded_in(size_t a, size_t b, size_t c, const std::vector<size_t>& cells) const;
    mutable wmtk::threading::enumerable_thread_specific<int> m_collapse_folds_before{-1},
        m_swap_folds_before{-1};
    mutable std::atomic<int> m_fold_reject_collapse{0}, m_fold_reject_swap{0},
        m_fold_reject_smooth{0};
 
    std::vector<std::array<double, 8>> optimization_metrics;
    int m_round = 0;
    mutable size_t m_debug_seq = 0;
    struct DebugTimeline
    {
        std::vector<std::string> labels;
        std::map<std::string, std::vector<size_t>> pvd_series;
    };
    DebugTimeline m_debug_timeline;
    std::future<DebugTimeline> m_debug_frame_write;
    static void write_debug_pvd(const std::string& out, const DebugTimeline& timeline);
    void wait_debug_frame_write();
    mutable int m_debug_pass = 0;
    mutable int m_debug_last_round = -1;
    mutable char m_debug_last_tag = '?';
    mutable std::string m_debug_last_prefix = "?";
    bool m_converged = false;
    bool m_quality_converged = true;
    double m_quality_max_amips = 0.;
 
    std::atomic<int> iter_cnt_split = 0, iter_cnt_collapse = 0, iter_cnt_swap = 0;
    std::atomic<int> iter_cnt_collapse_offset_removed{0};
    std::atomic<int> iter_cnt_collapse_offset_reject{0};
    std::atomic<int> iter_cnt_swap_offset_reject{0};
    std::atomic<long long> m_swap_scoring_checked{0};
    std::atomic<long long> m_swap_scoring_mismatch{0};
 
    mutable std::atomic<long long> funnel_offered{0};
    mutable std::array<std::atomic<long long>, 3> funnel_kind{};
    mutable std::atomic<long long> funnel_cases{0};
    mutable std::atomic<long long> funnel_case_inverted{0};
    mutable std::atomic<long long> funnel_case_not_better{0};
    mutable std::atomic<long long> funnel_case_better{0};
    mutable std::atomic<long long> funnel_after_cells{0};
    mutable std::atomic<long long> funnel_energy{0};
    mutable std::atomic<long long> funnel_committed{0};
    std::string flip_funnel_report() const;
    void flip_funnel_reset();
    std::atomic<int> iter_cnt_split_offset_before{0};
    std::atomic<int> iter_cnt_split_offset{0};
    std::atomic<long> m_split_order_waits{0};
    mutable std::atomic<long> m_split_off_longest{0};
    static bool& split_off_longest()
    {
        static thread_local bool off = false;
        return off;
    }
    void op_event(OpKind k, OpEvent e) const override;
    bool split_edge_is_due(const Tuple& e) const;
 
    struct OptSplitCache
    {
        bool is_edge_on_region = false;
        bool is_edge_on_offset = false;
        std::map<simplex::Edge, TetAttributes> tets;
        double parent_q_max = -1.;
        double parent_flatness = 1.;
    };
    wmtk::threading::enumerable_thread_specific<OptSplitCache> m_opt_split_cache;
 
    bool marching_split_edge_before(const Tuple& t);
    bool marching_split_edge_after(const Tuple& t);
    bool edge_split_sphere_trace(
        const Vector3d& p_in,
        const Vector3d& p_out,
        Vector3d& p_new,
        size_t& steps) const;
    size_t m_marching_root_splits = 0, m_marching_midpoint_splits = 0;
    size_t m_marching_trace_steps = 0, m_marching_trace_steps_max = 0;
    double m_construction_distance = 0.;
    double m_growth_start = 0.;
 
    bool collapse_edge_before(const Tuple& t) override;
    bool collapse_edge_after(const Tuple& t) override;
    bool collapse_before_vertex(size_t v1, size_t v2, double edge_length) override;
    bool collapse_after_connectivity(
        size_t v1,
        size_t v2,
        const std::vector<std::array<size_t, 2>>& boundary_edges) override;
    bool collapse_is_order_2_edge(const std::array<size_t, 2>& e) override
    {
        return is_order_2_edge(e);
    }
    void collapse_after_vertex(size_t v1, size_t v2) override;
 
    bool swap_before_interior(const std::vector<size_t>& tids) override;
    bool swap_before_surface(
        const std::vector<size_t>& tids,
        size_t a,
        size_t b,
        size_t c,
        size_t d) override;
    bool swap_after_cells(const std::vector<size_t>& tids, bool is_surface_flip) override;
 
    bool split_before_cells(const Tuple& edge, const std::vector<Tuple>& parents) override;
    bool split_after_cells(size_t v1, size_t v2, size_t v_new, const std::vector<Tuple>& children)
        override;
    bool split_adjust_position(size_t v_new, const std::vector<Tuple>& children) override;
    void split_after_vertex(size_t v_new, bool is_edge_open_boundary) override;
 
    bool is_open_boundary_edge(const Tuple& e) override { return is_order_2_edge(e); }
 
    bool smooth_before(const Tuple& t) override;
    bool smooth_after(const Tuple& t) override;
 
    bool vertex_is_on_domain_boundary(const size_t vid) const
    {
        return !m_vertex_attribute[vid].on_bbox_faces.empty();
    }
 
    void init_surfaces_and_boundaries();
 
    bool is_edge_on_region(const Tuple& loc);
    bool is_edge_on_offset(const Tuple& loc);
 
    void mark_input_complex_vertices();
 
    void warn_if_offset_reaches_domain_boundary() const;
 
    struct SmoothTrace
    {
        std::atomic<int> attempted{0}; 
        std::atomic<int> before_bbox{0}; 
        std::atomic<int> before_unrounded{0}; 
        std::atomic<int> offset_attempted{0}; 
        std::atomic<int> offset_accepted{0}; 
        std::atomic<int> interior_attempted{0}; 
        std::atomic<int> region_attempted{0}; 
 
        void reset()
        {
            for (std::atomic<int>* c :
                 {&attempted,
                  &before_bbox,
                  &before_unrounded,
                  &offset_attempted,
                  &offset_accepted,
                  &interior_attempted,
                  &region_attempted}) {
                c->store(0);
            }
        }
    };
    SmoothTrace m_smooth_trace;
 
    optimization::NewtonCounters m_newton_front;
    struct FrontSolveRecord
    {
        size_t vid;
        int iterations;
        int status; 
    };
    std::vector<FrontSolveRecord> m_front_solve_log;
    std::mutex m_front_solve_log_mutex;
    static constexpr double kSmoothRelGradNormTol = 1e-6;
    polysolve::nonlinear::Solver& smoothing_solver();
    std::atomic<size_t> m_front_veto_asked{0}, m_front_veto_fired{0};
    std::atomic<size_t> m_damping_moves{0}, m_damping_bounded{0};
    std::shared_ptr<const OffsetPotential3D> m_repulsion_potential;
    optimization::NewtonCounters m_newton_repulsion;
    std::atomic<size_t> m_repulsion_veto_asked{0}, m_repulsion_veto_fired{0};
    bool is_marched_edge(const size_t a, const size_t b) const
    {
        return (m_vertex_extra[a].label == 0) != (m_vertex_extra[b].label == 0);
    }
    bool is_repulsion_vertex(const size_t vid) const
    {
        if (!m_repulsion_potential || m_vertex_extra[vid].label != 0 || vertex_is_held(vid)) {
            return false;
        }
        for (const size_t u : get_one_ring_vids_for_vertex(vid)) {
            if (m_vertex_extra[u].label != 0) return true;
        }
        return false;
    }
    std::shared_ptr<OffsetEnergy3D> m_repulsion_term;
    double repulsion_term_at(const size_t vid) const
    {
        Eigen::VectorXd x = m_vertex_attribute[vid].m_posf;
        return m_repulsion_term->value(x);
    }
    double repulsion_cell_term(const std::array<size_t, 4>& vids) const
    {
        if (!m_repulsion_potential) return 0.;
        bool has_complex = false;
        for (const size_t v : vids) has_complex = has_complex || m_vertex_extra[v].label != 0;
        if (!has_complex) return 0.;
        double e = 0.;
        for (const size_t v : vids) {
            if (m_vertex_extra[v].label == 0 && !vertex_is_held(v)) {
                e += repulsion_term_at(v);
            }
        }
        return e;
    }
    size_t repulsion_cells_at(const size_t vid) const
    {
        size_t n = 0;
        for (const size_t tid : get_one_ring_tids_for_vertex(vid)) {
            for (const size_t u : oriented_tet_vids(tid)) {
                if (m_vertex_extra[u].label != 0) {
                    ++n;
                    break;
                }
            }
        }
        return n;
    }
    std::shared_ptr<OffsetEnergy3D> repulsion_energy(const size_t vid) const
    {
        return std::make_shared<OffsetEnergy3D>(
            m_repulsion_potential,
            double(repulsion_cells_at(vid)) * 4. * offset_term_weight(),
            true,
            true,
            /*one_sided=*/true);
    }
    bool repulsion_embedding_kept(const std::vector<size_t>& tids) const;
    bool simplex_in_input_complex(const std::vector<size_t>& c) const;
    std::atomic<size_t> m_repulsion_embed_refused{0};
    static constexpr int kGradBins = 22;
    std::array<std::atomic<size_t>, kGradBins> m_front_grad_abs{}, m_front_grad_rel{};
    void log_smoothing_pass_accounting() override;
    void smoothing_pass_done() override
    {
        if (m_energy_step_turn > 0) {
            log_energy_step(m_energy_step_first ? m_energy_step_group : "smooth");
            m_energy_step_first = false;
        }
        stamp_plastic_rests();
        if (m_energy_step_turn > 0 && m_plastic_active) log_energy_step("stamp");
    }
    void log_energy_step(const char* step) const;
    std::pair<double, double> energy_parts() const;
    template <typename F>
    void diagnostic_parallel_for(const size_t n, const size_t chunk, F&& fn) const
    {
        wmtk::threading::dynamic_parallel_for(n, std::max(1, NUM_THREADS), chunk, fn);
    }
    static constexpr size_t kDiagnosticChunk = 64;
    int m_energy_step_turn = 0;
    const char* m_energy_step_group = "";
    bool m_energy_step_first = false;
 
    std::vector<double> m_cross_ring;
    std::vector<Vector3d> m_cross_pos;
    bool m_cross_valid = false;
    std::atomic<size_t> m_cross_own{0}, m_cross_neighbour{0};
    std::vector<double> front_ring_measures() const;
    double ring_measure_at(size_t vid) const;
    void crossing_snapshot(const std::string& pass, bool compare, bool match_positions);
    void update_attributes() override;
 
    mutable wmtk::threading::enumerable_thread_specific<std::pair<double, Vector3d>> m_needle_pre;
    mutable std::atomic<size_t> m_needle_smooth_offered{0};
    mutable std::atomic<size_t> m_needle_smooth_reached{0};
    mutable std::atomic<size_t> m_needle_smooth_fixed{0};
    mutable std::atomic<size_t> m_needle_smooth_stationary{0};
    mutable std::atomic<size_t> m_needle_smooth_reports{0};
 
    double ring_max_quality(size_t vid) const;
    double tet_amips(const size_t tid) const { return std::cbrt(cell_quality(tid)); }
 
    double tet_flatness(size_t tid) const;
 
    void needle_forensics() const;
 
    void record_flatness(const char* op, double parent_flat, size_t child_tid) const;
    mutable std::atomic<size_t> m_flat_created_split{0};
    mutable std::atomic<size_t> m_flat_created_collapse{0};
    mutable std::atomic<size_t> m_flat_worsened_split{0};
    mutable std::atomic<size_t> m_flat_genesis_reports{0};
    static constexpr double kFlatThreshold = 1e-3;
    mutable wmtk::threading::enumerable_thread_specific<double> m_collapse_parent_flatness;
    mutable wmtk::threading::enumerable_thread_specific<double> m_collapse_survivor_sizing;
    struct CollapseVolume
    {
        double before = 0.;
        std::vector<size_t> v2_only;
        bool tie = false;
        bool band_cancels = false;
        double surface_before = 0.;
        std::vector<std::array<Vector3d, 2>> swept_faces;
        Vector3d swept_x1 = Vector3d::Zero(), swept_x2 = Vector3d::Zero();
    };
    mutable wmtk::threading::enumerable_thread_specific<CollapseVolume> m_collapse_volume;
    mutable wmtk::threading::enumerable_thread_specific<std::vector<size_t>> m_collapse_edge_link;
    void log_smooth_trace() const;
 
    void audit_surface_containment(const std::string& when) const;
 
 
    bool vertex_is_on_region(const size_t vid) const
    {
        return m_vertex_extra[vid].m_is_on_region || !m_vertex_attribute[vid].on_bbox_faces.empty();
    }
 
    std::shared_ptr<Envelope> surface_envelope_for_face(
        const std::array<size_t, 3>& vids,
        const FaceAttributes& /*attr*/) const override
    {
        const auto found = try_tuple_from_face(vids);
        if (!found) return m_envelope; // not reached: every caller asks about an existing face
        const Tuple& f = std::get<0>(*found);
        if (f.switch_tetrahedron(*this) && face_is_offset_surface_live(f)) {
            return m_tetwild_pass ? m_offset_envelope : nullptr;
        }
        return face_is_held(f) ? m_envelope : nullptr;
    }
 
    bool allow_surface_swap() const override { return true; }
 
    bool swap_quality_allowed(double after, double before, bool is_surface_flip) const override
    {
        // TetWild's rules (m_op_rules): the base's, after < before.
        return !tetwild_rules() ||
               TetOptimizerMesh::swap_quality_allowed(after, before, is_surface_flip);
    }
    double swap_edge_44_energy(const std::vector<std::array<size_t, 4>>& tets, const int op_case)
        override;
    double swap_edge_56_energy(const std::vector<std::array<size_t, 4>>& tets, const int op_case)
        override;
    bool swap_face_before(const Tuple& t) override;
    bool check_surface_topology() const override { return m_offset_params.perform_sanity_checks; }
 
    std::shared_ptr<Envelope> smoothing_energy_envelope(const size_t vid) const override
    {
        return smoothing_envelope(vid);
    }
    std::shared_ptr<Envelope> smoothing_containment_envelope(const size_t vid) const
    {
        return smoothing_envelope(vid);
    }
    std::shared_ptr<Envelope> smoothing_envelope(const size_t vid) const
    {
        if (vertex_is_held(vid)) return m_envelope;
        if (m_tetwild_pass && m_vertex_extra[vid].m_is_on_offset) return m_offset_envelope;
        return nullptr;
    }
    std::shared_ptr<Envelope> smoothing_containment_envelope(
        const size_t vid,
        const std::array<size_t, 3>& /*face*/,
        const FaceAttributes& /*attr*/) const override
    {
        return smoothing_containment_envelope(vid);
    }
 
    bool smooth_front_vertex(const Tuple& t);
    bool smooth_repulsion_vertex(const Tuple& t);
    bool smooth_nonfront_vertex(const Tuple& t);
    double damping_step(size_t vid) const;
    bool damping_active() const
    {
        return m_offset_params.damping && !m_tetwild_pass && !m_repulsion_potential;
    }
    void count_damped_move(size_t vid, const Vector3d& x0, double max_step);
    std::shared_ptr<polysolve::nonlinear::Problem> amips3_energy(size_t vid) const;
    double front_vertex_normal_gradient(size_t vid) const;
    Vector3d front_vertex_move_direction(size_t vid) const;
    double front_move_alignment(size_t vid) const;
    double front_vertex_conv_ratio(size_t vid) const;
    int m_smooth_group = 0;
    bool front_placed_by_ratio(const double ratio) const
    {
        return std::isfinite(ratio) && ratio <= 1.;
    }
    double face_offset_term(size_t a, size_t b, size_t c) const;
    double face_offset_term(
        const OffsetPotential3D& pot,
        const Vector3d& pa,
        const Vector3d& pb,
        const Vector3d& pc) const;
    Vector3d front_vertex_normal(size_t vid) const;
    std::shared_ptr<polysolve::nonlinear::Problem> front_objective(size_t vid) const;
    bool vertex_carries_offset_term(const size_t vid) const
    {
        return !m_tetwild_pass && m_offset_potential && m_vertex_extra[vid].m_is_on_offset &&
               !vertex_is_held(vid);
    }
    std::shared_ptr<polysolve::nonlinear::Problem> smoothing_extra_energy(
        const size_t vid) const override
    {
        // TetWild's rules (m_op_rules): no extra energy, so the engine's smoother minimises the
        // raw sum of AMIPS over the one-ring (AMIPSEnergy3D), as TetWild's does.
        if (tetwild_rules()) return nullptr;
        auto sum = std::make_shared<optimization::EnergySum>();
        sum->add_energy(amips3_energy(vid));
        // Before the march (repulsion_smoothing()) there is no front.
        if (is_repulsion_vertex(vid)) {
            sum->add_energy(repulsion_energy(vid));
            return sum;
        }
        if (vertex_carries_offset_term(vid)) {
            sum->add_energy(front_energy(vid));
        } else if (!vertex_is_held(vid)) {
            // E's band term moves with this vertex too -- except a vertex held on the input.
            // There d = 0 at the vertex and d has a kink: the term's slope along the input's
            // normal jumps across the surface, and its curvature there is 0 (or ~1e10 near an
            // input edge), so the Newton step along the normal is set by AMIPS alone and lands
            // far off; the projected step then discards the normal motion anyway. MEASURED
            // 2026-10-09, cube, serial, 3 turns: with the term these solves took 24.6 Newton
            // iterations and ~20 line-search trials per iteration (72% of the smoother's value
            // calls), without it 3.9 and 1; 19.3 s -> 11.3 s. E still guards the move:
            // smooth_nonfront_vertex() vetoes on the ring's sum of tet_energy().
            // Under "exact", not where it cancels (BAND TERM THAT CANCELS, tet_energy()): a ring
            // of band cells clear of the input's bisectors tiles the star of vid, so the term is a
            // constant, and its cells' exact integrals would add only work and their rounding to
            // every Newton step. The veto's own test (cancels_before in smooth_nonfront_vertex()).
            // Constant wherever the move can end, not only here: vid is not held, so not on the
            // domain wall (face_is_held()), and its ring closes around it; the engine keeps no
            // position that inverts a ring cell, so vid stays inside the convex hull of its link;
            // and the test reads the points only through their box, their least side of each
            // bisector polytope's planes, and min_P max d_P (d_P convex), none of which a point
            // inside that hull changes. So the test here certifies every position the smoother
            // can keep. "front_flux" and "exact_swept" leave the constant out by their own forms
            // (band_volume_energy_at()). MEASURED (cube 1e-2 / 1e-4, 15 rounds, serial, two runs
            // each): 659 of 46,594 non-front solves had such a ring; at every one's end the
            // veto's test found it clear again, and every one was accepted; smoothing 7.7-7.9 s ->
            // 7.4-7.5 s.
            if (m_offset_params.band_volume_rule == "exact" &&
                band_term_cancels(get_one_ring_tids_for_vertex(vid)))
                return sum;
            if (const auto bv = band_volume_energy_at(vid)) sum->add_energy(bv);
        }
        return sum;
    }
 
    std::shared_ptr<polysolve::nonlinear::Problem> front_energy(size_t vid) const;
    std::shared_ptr<polysolve::nonlinear::Problem> band_volume_energy_at(size_t vid) const;
 
 
    int stencil_order() const { return m_offset_params.stencil_order; }
    const std::vector<std::array<double, 3>>& face_stencil() const
    {
        const std::string& type = m_offset_params.stencil_type;
        const int k = m_offset_params.stencil_order;
        return type == "fem"        ? fem_face_stencil(k)
               : type == "interior" ? interior_face_stencil(k)
                                    : subdivision_face_stencil(k);
    }
    const std::vector<double>* face_stencil_weights() const
    {
        return m_offset_params.stencil_type == "fem"
                   ? &fem_face_stencil_weights(m_offset_params.stencil_order)
                   : nullptr;
    }
    int stencil_points_per_face() const
    {
        const int k = m_offset_params.stencil_order;
        if (m_offset_params.stencil_type == "subdivision") {
            if (k < 0) return 0;
            if (k == 0) return 3;
            const int n = 1 << (k - 1);
            return (n + 1) * (n + 2) / 2 + n * n;
        }
        if (k < 1) return 0;
        if (m_offset_params.stencil_type == "fem") return (k + 1) * (k + 2) / 2;
        return k * (k + 1) / 2;
    }
 
    void check_offset_within_support(const char* when) const;
 
    struct DistanceSplit
    {
        double max_reachable = 0., avg_reachable = 0.;
        double max_pinned = 0.;
        size_t n_reachable = 0, n_pinned = 0;
        double max_at_vertex = 0., max_in_face = 0.;
        size_t n_outside_support = 0;
        size_t worst_outside_vid = static_cast<size_t>(-1);
        double worst_outside_dist = 0.;
    };
    DistanceSplit distance_deviation_split() const;
 
    DistanceSplit residual_split() const;
    std::vector<bool> band_vertex_mask() const;
 
    double max_band_vertex_distance() const;
    double band_vertex_distance_error(const size_t vid) const;
 
    double band_vertex_residual(const size_t vid) const;
 
    struct FaceSamples
    {
        double max = 0.;
        double sum = 0.;
        size_t n = 0;
    };
 
    template <typename Visit>
    void for_each_face_sample(
        const Vector3d& p0,
        const Vector3d& p1,
        const Vector3d& p2,
        Visit&& visit) const
    {
        for (const auto& w : face_stencil()) {
            const Vector3d q(w[0] * p0 + w[1] * p1 + w[2] * p2);
            visit(q, w[0], w[1], w[2]);
        }
    }
 
    template <typename Visit>
    void for_each_offset_face_sample(const Tuple& f, Visit&& visit) const
    {
        const auto vs = get_face_vids(f);
        for_each_face_sample(
            m_vertex_attribute[vs[0]].m_posf,
            m_vertex_attribute[vs[1]].m_posf,
            m_vertex_attribute[vs[2]].m_posf,
            std::forward<Visit>(visit));
    }
 
    FaceSamples offset_face_samples(const Tuple& f) const;
 
    struct GradientSplit
    {
        double max_reachable = 0., avg_reachable = 0.;
        double max_pinned = 0.;
        size_t n_reachable = 0, n_pinned = 0;
        double max_at_vertex = 0., max_in_face = 0.;
        size_t n_face_samples = 0;
        size_t n_skipped_inverted = 0, n_skipped_unrounded = 0;
    };
    GradientSplit gradient_split(bool include_face_samples = true) const;
 
    struct EnergyCriterion
    {
        double max_vertex = 0., max_face = 0.; 
        double sum_vertex = 0., sum_face = 0.;
        double bar = 1.;
        size_t n_vertices = 0, n_faces = 0, n_unmeasurable = 0;
        size_t worst_vid = static_cast<size_t>(-1);
        size_t n_faces_over = 0, n_faces_over_placed = 0;
        double max_face_placed = 0.;
        Vector3d worst_placed_centroid = Vector3d::Zero();
        double bar_length = 0.; 
        size_t n_at_floor = 0;
        double max_face_at_floor = 0.; 
        Vector3d worst_at_floor_centroid = Vector3d::Zero();
        double worst_at_floor_scalar = 0.; 
        size_t n_corners_at_floor = 0;
        double max_face_corners_at_floor = 0.; 
        Vector3d worst_corners_at_floor_centroid = Vector3d::Zero();
        double floor_scalar = 0.;
        bool floor_from_min_edge_length = false;
        size_t n_unplaced = 0; 
        struct Refinable
        {
            size_t a, b, c;
            double measure, len;
        };
        std::vector<Refinable> refinable;
        bool ring_exit = false; 
        static const char* ring_name() { return "ring measure"; }
        double max_ring = 0., sum_ring = 0.; 
        size_t n_rings = 0, n_rings_unmeasurable = 0;
        size_t worst_ring_vid = static_cast<size_t>(-1);
        size_t n_rings_over = 0; 
        std::vector<size_t> refinable_vertices;
        size_t n_rings_at_floor = 0;
        double max_ring_at_floor = 0.; 
        Vector3d worst_ring_at_floor_pos = Vector3d::Zero();
        bool rings_ok() const { return max_ring <= bar; }
        double avg_ring() const { return n_rings ? sum_ring / double(n_rings) : 0.; }
        bool vertices_ok() const { return n_unplaced == 0; }
        bool faces_ok() const { return max_face <= bar; }
        bool converged() const
        {
            return (ring_exit ? rings_ok() : faces_ok()) && n_unmeasurable == 0;
        }
        std::string sizing_floor_fact() const;
        double ratio() const { return bar > 0. ? std::max(max_vertex, max_face) / bar : 0.; }
        double avg_vertex() const { return n_vertices ? sum_vertex / double(n_vertices) : 0.; }
        double avg_face() const { return n_faces ? sum_face / double(n_faces) : 0.; }
    };
    EnergyCriterion energy_criterion();
    double front_chord_target(size_t va, size_t vb, double len, double sag, double bar) const;
 
    size_t refine_front_by_halving(const std::vector<EnergyCriterion::Refinable>& faces);
    size_t refine_front_by_halving(const std::vector<size_t>& vertices);
    bool settle_refinement_turn(int turn, double& e_prev, bool refine_now);
 
    void grade_sizing(double grade, const std::vector<size_t>& seeds);
    size_t grade_sizing_by_distance(const std::vector<size_t>& seeds);
 
    std::optional<EnergyCriterion> m_energy_verdict;
 
    void report_outside_support(const char* when, const DistanceSplit& s) const;
    bool band_vertex_is_reachable(const size_t vid) const
    {
        if (m_vertex_extra[vid].m_is_on_offset && vertex_is_held(vid)) return false;
        return !vertex_is_on_domain_boundary(vid);
    }
 
    size_t refine_sizing_around_worst(double max_metric) override;
 
    void log_stuck_refine_census(double max_metric, double filter_energy);
 
    void log_refine_block_census(const std::string& when, double filter_energy) const;
 
    bool collapse_quality_allowed(size_t v1, double q, double ring_max) const override
    {
        // TetWild's rules (m_op_rules): the base's, q <= ring_max for a rounded v1.
        return !tetwild_rules() || TetOptimizerMesh::collapse_quality_allowed(v1, q, ring_max);
    }
 
    mutable std::atomic<size_t> m_deg_split_created{0};
    size_t m_deg_prev_split_created = 0;
 
    void report_needle(const char* op, size_t tid, double parent_q) const;
    static constexpr size_t kNeedleReports = 12;
    static constexpr double kNeedleQuality = 1e6;
    mutable std::atomic<size_t> m_needle_reports{0};
 
    void needle_scan(const char* when) const;
 
    std::set<std::tuple<long, long, long>> m_stuck_prev_cells;
    size_t m_stuck_calls = 0;
 
    bool optimization_bare_coarsen_passes() const override { return tetwild_rules(); }
 
    bool coarsen_no_rise() const { return m_offset_params.coarsen_guard == "no_rise"; }
 
    double face_criterion_rel(const Tuple& f) const;
    double cell_quality_rel(const size_t tid) const;
    double amips_rel_at_face(const Tuple& f) const;
 
    double quality_rel(const size_t tid) const override
    {
        if (!coarsen_no_rise()) return TetOptimizerMesh::quality_rel(tid);
        return std::max(1., std::cbrt(cell_quality(tid)) / m_params.stop_energy);
    }
    std::atomic<size_t> m_coarsen_reject_amips{0}, m_coarsen_reject_offset{0};
    struct CoarsenCollapse
    {
        double amips3_before = 0.;
        std::vector<size_t> v2_only;
        std::vector<std::array<size_t, 3>> faces_before;
    };
    mutable wmtk::threading::enumerable_thread_specific<CoarsenCollapse> m_coarsen_collapse;
 
    void write_optimization_debug_output(const std::string& path) override
    {
        // DEBUG_remote_output: no per-operation or per-smoothing-pass frame (remote_frames()).
        if (remote_frames() && path.rfind("debug_", 0) == 0) return;
        const char ph = m_tetwild_pass ? 'F' : 'S'; // a TetWild pass, or the loop
        if (m_round != m_debug_last_round || ph != m_debug_last_tag ||
            m_frame_prefix != m_debug_last_prefix) {
            m_debug_last_round = m_round;
            m_debug_last_tag = ph;
            m_debug_last_prefix = m_frame_prefix;
            m_debug_pass = 0;
        }
        std::string label = path;
        // The prefix leads, since the loop after the init_optimize loop restarts at turn 1.
        const std::string& g = m_frame_prefix;
        if (path.rfind("debug_", 0) == 0) {
            label = fmt::format(
                "{}r{}{}{}{}",
                g,
                m_round,
                ph,
                ++m_debug_pass,
                m_debug_pass_name.empty() ? std::string() : "_" + m_debug_pass_name);
        } else if (path.rfind("end_", 0) == 0) {
            label = fmt::format("{}r{}{}_end", g, m_round, ph);
        }
        write_debug_frame(label);
    }
    bool remote_frames() const { return m_offset_params.debug_remote_output; }
    bool stage_frames() const { return m_offset_params.debug_output && !remote_frames(); }
    void optimization_iteration_done(const int it) override
    {
        if (remote_frames()) {
            write_debug_frame(fmt::format("{}r{}F_it{}", m_frame_prefix, m_round, it + 1));
        }
    }
    void append_frame_label(size_t idx, const std::string& label) const;
    void write_debug_frame(const std::string& label);
 
    void init_from_image(
        const MatrixXd& V,
        const MatrixXi& T,
        const MatrixSi& T_tags,
        const MatrixXd& V_env,
        const MatrixXi F_env,
        const std::vector<std::string>& tag_names,
        const std::string& sheet_name = "",
        const MatrixXi& F_tagged = MatrixXi(),
        const MatrixSi& F_tags = MatrixSi(),
        const std::vector<std::string>& face_tag_names = {});
 
    bool ambient_assert();
 
    void label_input_complex();
 
    bool empty_input_complex();
 
    void init_input_complex_bvh();
 
    void init_offset_potential();
 
    void init_band_exact();
 
    MatrixXd m_phi_V;
    MatrixXi m_phi_E;
    MatrixXi m_phi_F;
    std::vector<int8_t> m_phi_F_inner;
    std::vector<int> m_phi_P;
 
    size_t flood_fill();
 
    bool split_edge_before(const Tuple& t) override;
    bool split_edge_after(const Tuple& t) override;
    bool split_face_before(const Tuple& t) override;
    bool split_face_after(const Tuple& t) override;
    bool split_tet_before(const Tuple& t) override;
    bool split_tet_after(const Tuple& t) override;
    bool invariants(const std::vector<Tuple>& tets) override;
 
    void construct_offset(const std::filesystem::path& output_file);
 
    void marching_tets();
 
    void repulsion_smoothing();
    void relabel_input_complex();
 
    bool is_simplicially_embedded() const;
    bool tet_is_simp_emb(const Tuple& t) const;
    void simplicial_embedding();
 
    void set_offset_tet_tags();
 
    bool offset_is_manifold();
 
    void write_phi_grid(const std::string& path, int n) const;
 
    void write_input_complex(const std::string& path);
    void write_vtu(const std::string& path, MeshMap* output = nullptr);
    struct VtuFrame
    {
        std::string path;
        MatrixXd V;
        MatrixXi T, F_in, F_off, E;
        // per cell
        VectorXd amips, tet_e, amips_term, band_term, band_gap, plastic;
        std::vector<MatrixXd> tags;
        std::vector<std::string> tag_names;
        // per vertex
        VectorXd labels, order, vid, sizing;
        VectorXd in_dist, dist_ratio, front_err, front_change;
        // the _off file's: per offset face, and the ring measure per vertex
        VectorXd f_err, f_len, f_inside, f_outside, f_fold, f_field_fold, f_field_angle, f_dihedral,
            f_fold_dihedral, ring;
        // the _stencil file's: one point per stencil point of every offset face (S, a vertex
        // cell each in S_cells), and its fields
        MatrixXd S;
        MatrixXi S_cells;
        VectorXd s_err, s_over, s_signed, s_dist_ratio, s_sq, s_face, s_point;
    };
    VtuFrame gather_vtu_frame(const std::string& path);
    static void write_vtu_frame(const VtuFrame& frame, MeshMap* output = nullptr);
 
private:
    struct EdgeSplitCache
    {
        size_t v1_id;
        size_t v2_id;
        Vector3d new_v_pos;
        VertexExtra new_v_extra;
 
        bool is_edge_on_region = false;
        bool is_edge_on_offset = false;
        bool is_edge_open_boundary = false;
 
        std::vector<std::pair<FaceAttributes, std::array<size_t, 3>>> changed_faces;
 
        // cache edge attributes
        EdgeAttributes split_e;
        std::map<size_t, EdgeAttributes> internal_e;
        std::map<simplex::Edge, EdgeAttributes> external_e; // edge is boundary edge (not link)
        std::map<simplex::Edge, EdgeAttributes> link_e; // link edge around splitted edge
 
        // cache face attributes
        std::map<size_t, FaceSnapshot> split_f; // splitted faces
        std::map<simplex::Edge, FaceSnapshot> internal_f; // new faces created by split
        std::map<std::pair<simplex::Edge, size_t>, FaceSnapshot>
            external_f; // closed star boundary faces of splitted edge
 
        // cache tet attributes
        std::map<simplex::Edge, TetAttributes> tets;
    };
    wmtk::threading::enumerable_thread_specific<EdgeSplitCache> edge_split_cache;
 
    struct FaceSplitCache
    {
        size_t v1_id;
        size_t v2_id;
        size_t v3_id;
        std::map<simplex::Edge, EdgeAttributes> existing_e;
        std::map<simplex::Face, FaceSnapshot> existing_f;
        int splitf_label;
        std::map<size_t, TetAttributes> tets;
    };
    wmtk::threading::enumerable_thread_specific<FaceSplitCache> face_split_cache;
 
    struct TetSplitCache
    {
        std::array<size_t, 4> v_ids;
        std::map<simplex::Edge, EdgeAttributes> existing_e;
        std::map<simplex::Face, FaceSnapshot> existing_f;
        TetAttributes tet;
    };
    wmtk::threading::enumerable_thread_specific<TetSplitCache> tet_split_cache;
 
    bool swap_capture_tag(const std::vector<size_t>& tids);
    wmtk::threading::enumerable_thread_specific<CellTag> m_swap_tag;
    wmtk::threading::enumerable_thread_specific<int> m_swap_label;
 
    struct SwapSurfaceSides
    {
        std::map<size_t, std::pair<CellTag, int>> by_vertex;
        std::array<size_t, 4> abcd{};
        bool worthwhile = false;
        double case0_energy = 0.0;
        bool saw_case = false;
        int kind = 0;
    };
    bool swap_capture_surface_sides(
        const std::vector<size_t>& tids,
        size_t a,
        size_t b,
        size_t c,
        size_t d);
    wmtk::threading::enumerable_thread_specific<SwapSurfaceSides> m_swap_sides;
    mutable wmtk::threading::enumerable_thread_specific<double> m_swap_energy_before;
    mutable std::atomic<int> m_plastic_tie_swap{0}, m_plastic_tie_collapse{0};
    struct SwapRecord
    {
        bool scored = false;
        double scored_energy = std::numeric_limits<double>::max();
        bool flip = false;
        bool band_cancels = false;
        std::vector<std::array<size_t, 4>> old_band;
        void clear()
        {
            flip = false;
            band_cancels = false;
            old_band.clear();
            scored = false;
            scored_energy = std::numeric_limits<double>::max();
        }
    };
    mutable wmtk::threading::enumerable_thread_specific<SwapRecord> m_swap_record;
    void swap_record_fill(bool flip);
    void swap_scoring_check(const std::vector<size_t>& tids, double after);
 
public:
    // substructure functions
 
    bool is_order_2_edge(const Tuple& e) const;
    bool is_order_2_edge(const std::array<size_t, 2>& e) const;
 
    bool vertex_is_on_surface(const size_t vid) const override;
 
    bool face_is_on_surface(const size_t fid) const override;
 
    size_t get_order_of_vertex(const size_t vid) const override;
    void init_vertex_order();
 
private: // helpers
    void sort_edges_by_length(std::vector<simplex::Edge>& edges)
    {
        std::sort(
            edges.begin(),
            edges.end(),
            [this](const simplex::Edge& e1, const simplex::Edge& e2) {
                double len1 = (m_vertex_attribute[e1.vertices()[0]].m_posf -
                               m_vertex_attribute[e1.vertices()[1]].m_posf)
                                  .norm();
                double len2 = (m_vertex_attribute[e2.vertices()[0]].m_posf -
                               m_vertex_attribute[e2.vertices()[1]].m_posf)
                                  .norm();
                return len1 > len2;
            });
    }
 
public: // helpers
    std::vector<size_t> connected_components_helper(const size_t& v_id)
    {
        auto onering_v_ids = get_one_ring_vids_for_vertex(v_id);
        std::vector<size_t> ret_v_ids;
        for (const size_t& other_v_id : onering_v_ids) {
            size_t e_id = tuple_from_edge({{v_id, other_v_id}}).eid(*this);
            if (m_edge_attribute[e_id].label != 0) { // edge labelled 1 or 2
                ret_v_ids.push_back(other_v_id);
            }
        }
        return ret_v_ids;
    }
};
 
 
} // namespace wmtk::components::topological_offset