TopoOffsetTetMesh.h

#pragma once
 
#include <array>
#include <cstdint>
#include <cstdlib>
#include <map>
#include <mutex>
#include <string>
 
#include <wmtk/TetMesh.h>
#include <wmtk/TetOptimizerMesh.h>
#include <wmtk/envelope/Envelope.hpp>
#include <wmtk/optimization/solver.hpp>
#include <wmtk/simplex/Simplex.hpp>
#include <wmtk/threading/enumerable_thread_specific.hpp>
#include "OffsetPotential.hpp"
#include "Parameters.h"
#include "SimplicialComplexBVH.hpp"
 
// clang-format off
#include <wmtk/utils/DisableWarnings.hpp>
#include <wmtk/utils/EnableWarnings.hpp>
// clang-format on
 
using CellTag = std::set<int64_t>;
 
 
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
 
class VertexExtra
{
public:
    int label = 0;
    size_t component_id = 0;
    bool m_is_on_region = false;
    bool m_is_on_offset = false; // is this vertex on the offset surface
 
    bool m_is_on_input_complex = false;
 
    uint64_t m_boundary_mask = 0;
 
    uint32_t m_born_epoch = 0;
};
 
 
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
};
 
 
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
};
 
 
class TopoOffsetTetMesh : public wmtk::TetOptimizerMesh
{
public:
    static constexpr int OFFSET_SURFACE_CLASS = 1;
 
    enum class OptPhase { A, B };
 
    OptPhase m_phase = OptPhase::A;
 
    enum class PhaseBSub { Offset, Background };
    PhaseBSub m_phase_b_sub = PhaseBSub::Offset;
 
    std::shared_ptr<SampleEnvelope> m_offset_envelope;
 
    void rebuild_offset_envelope();
 
    void check_no_vertex_on_both_surfaces(const char* when) const;
 
    void optimize_offset_alternating();
 
    size_t phase_b_smooth();
 
    bool smooth_interior_vertex_phase_b(const Tuple& t);
 
    double phase_b_band_gradient_linf();
 
    size_t update_band_sizing_from_tolerance();
 
public: // mode for splitting in marching tets
    enum class EdgeSplitMode {
        Midpoint = 0, // construction: simplicial embedding AND marching_tets
        Optimization = 5 // this is used in the optimization phase of the algorithm
    };
 
 
public:
    int m_vtu_counter = 0;
    std::array<size_t, 4> m_init_counts = {{0, 0, 0, 0}};
    size_t m_tags_count;
    SimplicialComplexBVH m_input_complex_bvh;
 
    std::shared_ptr<OffsetPotential3D> m_offset_potential;
 
    std::shared_ptr<SampleEnvelope> m_input_complex_envelope;
 
    std::map<int64_t, std::shared_ptr<SampleEnvelope>> m_tag_envelopes;
 
    std::map<int64_t, int> m_tag_bit;
 
    mutable std::map<uint64_t, std::shared_ptr<SampleEnvelope>> m_isect_cache;
    mutable std::mutex m_isect_mutex;
 
    EdgeSplitMode m_edge_split_mode = EdgeSplitMode::Midpoint;
 
    // tag map stuff
    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;
 
    // dont actually use, just for retaining in output
    bool m_has_envelope;
    std::vector<Vector3d> m_V_envelope;
    std::vector<Vector3i> m_F_envelope;
    // m_envelope itself lives on the base, which is what checks tracked-surface triangles
    // against it.
    double m_envelope_eps = -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 two are registered
    // alongside them in its attribute groups.
    VertexExtraCol m_vertex_extra;
    FaceExtraCol m_face_extra;
    EdgeAttCol m_edge_attribute;
    TetAttCol m_tet_attribute;
 
    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 set_vertex_position(const size_t vid, const Vector3d& p)
    {
        m_vertex_attribute[vid].m_posf = p;
        m_vertex_attribute[vid].m_pos = to_rational(p);
        m_vertex_attribute[vid].m_is_rounded = true;
    }
 
    bool cell_in_region(const size_t tid) const
    {
        const int l = m_tet_attribute[tid].label;
        return l == 1 || l == 2; // input complex, or the offset band
    }
 
    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;
    }
 
    // debug use
    std::atomic<int> cnt_split = 0, cnt_collapse = 0;
 
    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 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);
        p_edge_attrs = &m_edge_attribute;
        p_tet_attrs = &m_tet_attribute;
 
        m_collapse_check_link_condition = false;
        m_collapse_check_manifold = false;
 
        // TetWild parity, forced on for 3D regardless of the config key: Phase B places the
        // offset surface with its own root find, which never consults this flag, and Phase A is
        // pure TetWild quality work, where the veto is on. The spec default is false for the 2D
        // pipeline, whose offset-boundary placement the veto would freeze.
        _m_offset_params.smooth_quality_veto = true;
 
        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
    {
        // Spike tracker. Every operation that writes a quality comes through here, so this is
        // the one place that sees a tet cross from ordinary into absurd, and the spike's
        // position in the log names the pass that caused it. Only the crossing is reported, and
        // only the first few in full, or the flood would be the whole log.
        if (m_spike_track) {
            const double before = m_tet_attribute[tid].m_quality;
            const double a = spike_amips(q), b = spike_amips(before);
            // Two events, not a fixed bar -- a fixed bar spends the whole dump budget on the
            // mildest tets and never reaches the extremes.
            //   genesis: an ordinary tet ruined by one operation, i.e. the causal event.
            //   record:  a new global maximum, so the dumps end at the true worst and, since
            //            records only increase, the output stays bounded.
            if (b < m_spike_ordinary && a > m_spike_threshold) {
                log_quality_spike(tid, before, q, "genesis");
            }
            if (a > m_spike_record) {
                m_spike_record = a;
                log_quality_spike(tid, before, q, "record");
            }
            if (a > m_spike_threshold) ++m_spike_count;
        }
        m_tet_attribute[tid].m_quality = q;
    }
 
    static double spike_amips(const double q)
    {
        return std::isfinite(q) ? std::cbrt(q) : std::numeric_limits<double>::infinity();
    }
 
    void log_quality_spike(size_t tid, double before, double after, const char* kind) const;
 
    bool m_spike_track = true;
    double m_spike_threshold = 1e5;
    double m_spike_ordinary = 100.;
    mutable double m_spike_record = 0.;
    mutable std::atomic<int> m_spike_dumps{0};
    mutable std::atomic<int> m_spike_genesis_dumps{0};
    mutable std::atomic<int> m_spike_count{0};
    mutable std::atomic<int> m_spike_genesis{0};
    int m_spike_dump_budget = 24;
 
    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();
    }
 
    uint64_t tag_bits(const CellTag& tags) const
    {
        uint64_t bits = 0;
        for (const int64_t t : tags) {
            const auto it = m_tag_bit.find(t);
            if (it != m_tag_bit.end()) bits |= (uint64_t(1) << it->second);
        }
        return bits;
    }
 
    uint64_t vertex_boundary_mask(const size_t vid) const
    {
        return vertex_is_on_region(vid) ? m_vertex_extra[vid].m_boundary_mask : uint64_t(0);
    }
 
    uint64_t face_mask(const std::array<size_t, 3>& vids) const
    {
        return vertex_boundary_mask(vids[0]) & vertex_boundary_mask(vids[1]) &
               vertex_boundary_mask(vids[2]);
    }
 
    std::shared_ptr<SampleEnvelope> envelope_for_mask(uint64_t mask) const;
 
    std::shared_ptr<SampleEnvelope> surface_envelope_for_face(
        const std::array<size_t, 3>& vids) const override
    {
        // Boundary geometry first, in both phases: a face on any tag-region boundary may not
        // drift out of that boundary's tube, and a face on several is held in their
        // intersection. The mask carries the input complex too -- every complex simplex lies on
        // tag boundaries -- and E_t is built from the input mesh before construction touches it.
        if (const uint64_t mask = face_mask(vids)) {
            return envelope_for_mask(mask);
        }
        // Phase A holds the offset where Phase B left it; Phase B is what moves it, so it is
        // unconstrained there. Null when there is no offset envelope yet -- the construction
        // phase runs before the first one is built.
        bool all_offset = true;
        for (const size_t v : vids) {
            all_offset = all_offset && m_vertex_extra[v].m_is_on_offset;
        }
        if (all_offset && m_phase == OptPhase::A) return m_offset_envelope;
        return nullptr;
    }
 
    bool allow_surface_swap() const override { return true; }
    bool check_surface_topology() const override { return m_offset_params.perform_sanity_checks; }
 
    std::shared_ptr<SampleEnvelope> smoothing_energy_envelope(const size_t vid) const override
    {
        if (m_vertex_extra[vid].m_is_on_offset && !vertex_is_on_region(vid)) {
            return nullptr;
        }
        const uint64_t mask = vertex_boundary_mask(vid);
        if (mask == 0) {
            // Reachable only for a construction artefact (a wall-chord midpoint flagged
            // on-surface with disjoint endpoint masks); its containment is vacuous too, so no
            // pull is behavior-neutral. Not an error.
            return nullptr;
        }
        std::shared_ptr<SampleEnvelope> best;
        double worst_d2 = -1.;
        for (const auto& [tag, env] : m_tag_envelopes) {
            const auto it = m_tag_bit.find(tag);
            if (it == m_tag_bit.end() || !(mask & (uint64_t(1) << it->second))) continue;
            if (!best) {
                best = env;
                if ((mask & (mask - 1)) == 0) break; // single bit: no violation contest to run
                worst_d2 = env->squared_distance(m_vertex_attribute[vid].m_posf);
                continue;
            }
            const double d2 = env->squared_distance(m_vertex_attribute[vid].m_posf);
            if (d2 > worst_d2) {
                worst_d2 = d2;
                best = env;
            }
        }
        return best;
    }
 
    std::shared_ptr<SampleEnvelope> smoothing_containment_envelope(const size_t vid) const override
    {
        if (m_vertex_extra[vid].m_is_on_offset && !vertex_is_on_region(vid)) {
            return m_phase == OptPhase::A ? m_offset_envelope : nullptr;
        }
        // Boundary geometry is contained in the intersection of its tags' tubes, in BOTH
        // phases -- the caller only asks is_outside(triangle), which composites answer, so
        // unlike the pull this side may hand out an IntersectionEnvelope. Vertices with no
        // boundary bits (offset handled above; construction artefacts) have nothing to be
        // contained in.
        return envelope_for_mask(vertex_boundary_mask(vid));
    }
 
    // Do not re-add a smoothing_extra_energy() override; measured worse -- see git history of
    // this file. Phase B places the offset surface with its own root find and Phase A holds it
    // in m_offset_envelope, so the base's nullptr is correct in both phases.
 
    void write_optimization_debug_output(const std::string& path) override
    {
        // The per-pass debug_{N} series is opt-in: the driver's phase_{round}{A|B} frames are
        // the timeline the viewer shows, and the per-pass frames cost more time and disk than
        // the rest of the run put together.
        if (!m_offset_params.debug_output_per_pass && path.rfind("debug_", 0) == 0) {
            return;
        }
        write_vtu(m_offset_params.output_path + "_" + path);
    }
 
    size_t refine_sizing_around_worst(double) override;
 
    // optimization_bare_coarsen_passes() is deliberately not overridden: TetWild's opening and
    // closing unlimited-length collapse passes apply. Phase A rebuilds m_offset_envelope before
    // every mesh_improvement(), so the offset surface is held during them exactly as TetWild's
    // input surface is, and decimating to what the envelope tolerates is the contract.
 
    std::tuple<double, double> optimization_quality_stats() override;
 
    double optimization_stop_metric() const override
    {
        return m_phase == OptPhase::A ? wmtk::TetOptimizerMesh::optimization_stop_metric() : 1.;
    }
 
    double offset_residual_tolerance() const
    {
        return std::max(
            0.5 * m_offset_params.front_conv_rel * m_offset_params.target_distance,
            1e-16);
    }
 
    double offset_gradient_tolerance() const
    {
        // Normalized by the field's slope s (OffsetPotential::level_set_slope()), so
        // front_conv_rel means the same thing whichever field is in use: |Phi - c| is a field
        // difference, not a length, and s is the conversion. s == 1 for `euclidean`, where this
        // is a no-op; for `smooth` Phi is a barrier with s ~ 1/delta, growing as the offset gets
        // finer, so without the factor rel silently asks for a tiny fraction of what it names.
        //
        // At the field's reference slope the bound reduces to |Phi - c| / s <= (rel / 2) * delta,
        // so rel 0.2 means "within 10% of target_distance"; elsewhere it is that bound relaxed by
        // s / |grad Phi|. The relaxation is deliberate, not slack: it is what lets the criterion
        // be met at a stationary point of Phi, where |grad Phi| -> 0 and no residual bound is
        // achievable at all. Two bodies close enough that Phi's minimum contour between them
        // exceeds c have no level set in the gap; the sheets converge onto that contour from
        // either side and rel controls how close they get. Expect `converged: true` beside a
        // plateaued max_phi_residual there -- that geometry being reported, not a placement
        // failure.
        const double s = m_offset_potential ? m_offset_potential->level_set_slope() : 1.;
        return std::max(
            m_offset_params.front_conv_rel * m_offset_params.target_distance * s * s,
            1e-16);
    }
 
    void check_offset_within_support(const char* when) const;
 
    double band_vertex_residual(const size_t vid) const;
 
    bool band_vertex_is_reachable(const size_t vid) const
    {
        return m_vertex_attribute[vid].on_bbox_faces.empty();
    }
 
    std::vector<bool> band_vertex_mask() const;
 
    double max_band_vertex_distance() const;
 
    void write_phi_grid(const std::string& path, int n) const;
 
    struct FaceSamples
    {
        double max = 0.;
        double sum = 0.;
        size_t n = 0;
    };
 
    template <typename Visit>
    void for_each_offset_face_sample(const Tuple& f, Visit&& visit) const
    {
        const int k = m_offset_params.offset_residual_samples;
        if (k <= 0) return;
 
        const auto vs = get_face_vids(f);
        const Vector3d p0 = m_vertex_attribute[vs[0]].m_posf;
        const Vector3d p1 = m_vertex_attribute[vs[1]].m_posf;
        const Vector3d p2 = m_vertex_attribute[vs[2]].m_posf;
 
        const int n = k + 2;
        for (int i = 1; i < n; ++i) {
            for (int j = 1; j < n - i; ++j) {
                const int l = n - i - j;
                if (l < 1) continue;
                visit(Vector3d((double(i) * p0 + double(j) * p1 + double(l) * p2) / double(n)));
            }
        }
    }
 
    FaceSamples offset_face_samples(const Tuple& f) const;
 
    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;
 
    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 residual_split() 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;
        size_t worst_vid = static_cast<size_t>(-1);
    };
    GradientSplit gradient_split(bool include_face_samples = true) const;
    void report_outside_support(const char* when, const DistanceSplit& s) const;
 
    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 collapse_edge_before(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;
 
 
    struct OffsetCollapseRefusals
    {
        std::atomic<long> invariant{0}; 
        std::atomic<long> class_region{0}; 
        std::atomic<long> class_offset{0}; 
        std::atomic<long> order2{0}; 
        std::atomic<long> sublink{0}; 
        void reset() { invariant = class_region = class_offset = order2 = sublink = 0; }
    };
    mutable OffsetCollapseRefusals m_offset_collapse_refusals;
    // The shared collapse merges the sizing scalar with min(removed, survivor) -- see
    // TetOptimizerMeshCollapse.cpp -- making every collapse a one-way refinement ratchet. That
    // is deliberate: a collapse must never discard a finer request the field has already made.
 
    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); }
 
private:
    struct OptSplitCache
    {
        bool is_edge_on_region = false;
        bool is_edge_on_offset = false;
        std::map<simplex::Edge, TetAttributes> tets;
    };
    wmtk::threading::enumerable_thread_specific<OptSplitCache> m_opt_split_cache;
 
public:
private:
public:
private:
    wmtk::threading::enumerable_thread_specific<double> m_collapse_offset_rel_before;
    wmtk::threading::enumerable_thread_specific<double> m_swap_offset_rel_before;
 
    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;
 
public:
    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);
 
    void init_surfaces_and_boundaries();
 
    bool is_edge_on_region(const Tuple& loc);
    bool is_edge_on_offset(const Tuple& loc);
 
    bool ambient_assert();
 
    void label_input_complex();
 
    void mark_input_complex_vertices();
 
    bool empty_input_complex();
 
    void init_input_complex_bvh();
 
    void init_offset_potential();
 
    MatrixXd m_phi_V;
    MatrixXi m_phi_E;
    MatrixXi m_phi_F;
    std::vector<int> m_phi_P;
 
 
    size_t flood_fill();
 
    std::vector<std::array<size_t, 3>> get_faces_by_condition(
        std::function<bool(const FaceAttributes&)> cond) const;
 
    bool split_edge_before(const Tuple& t) override;
    bool split_edge_after(const Tuple& t) override;
    bool marching_split_edge_before(const Tuple& t);
    bool marching_split_edge_after(const Tuple& t);
    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;
    bool smooth_before(const Tuple& t) override;
    bool smooth_after(const Tuple& t) override;
 
    bool smooth_offset_vertex_backtracking(const Tuple& t);
 
    void execute_offset(const std::filesystem::path& output_file);
 
    void optimize_offset(const std::filesystem::path& output_file);
 
 
    bool is_offset_face(const Tuple& f) const;
    bool is_offset_face(const size_t fid) const;
 
    std::vector<Tuple> get_offset_surface_faces_for_vertex(const Tuple& t) const;
 
    struct SmoothTrace
    {
        std::atomic<int> attempted{0}; 
        std::atomic<int> before_bbox{0}; 
        std::atomic<int> before_unrounded{0}; 
        std::atomic<int> before_on_region{0}; 
        std::atomic<int> before_phase_b_not_offset{0}; 
        std::atomic<int> before_phase_b_enveloped_background{
            0}; 
        std::atomic<int> offset_attempted{0}; 
        std::atomic<int> interior_attempted{0}; 
 
        void reset()
        {
            for (std::atomic<int>* c :
                 {&attempted,
                  &before_bbox,
                  &before_unrounded,
                  &before_on_region,
                  &before_phase_b_not_offset,
                  &before_phase_b_enveloped_background,
                  &offset_attempted,
                  &interior_attempted}) {
                c->store(0);
            }
        }
    };
    mutable SmoothTrace m_smooth_trace;
 
    void log_smooth_trace() const;
 
    mutable std::atomic<int> m_phase_b_constrained{0};
 
    mutable wmtk::threading::enumerable_thread_specific<
        std::unique_ptr<polysolve::nonlinear::Solver>>
        m_phase_b_solver;
 
    enum class VClass { Bbox = 0, Offset = 1, Input = 2, Interior = 3, Count = 4 };
 
    VClass vertex_class(const size_t vid) const
    {
        if (!m_vertex_attribute[vid].on_bbox_faces.empty()) return VClass::Bbox;
        if (m_vertex_extra[vid].m_is_on_offset) return VClass::Offset;
        if (m_vertex_extra[vid].m_is_on_region) return VClass::Input;
        return VClass::Interior;
    }
 
    static const char* vclass_name(VClass c)
    {
        switch (c) {
        case VClass::Bbox: return "bbox";
        case VClass::Offset: return "offset";
        case VClass::Input: return "input";
        default: return "interior";
        }
    }
 
    struct MoveStats
    {
        std::atomic<int> attempted{0};
        std::atomic<int> refused_before{0};
        std::atomic<int> accepted{0};
        std::atomic<double> sum_disp{0.};
        std::atomic<double> max_disp{0.};
 
        void add(double d)
        {
            ++accepted;
            for (double cur = sum_disp.load(); !sum_disp.compare_exchange_weak(cur, cur + d);) {
            }
            for (double cur = max_disp.load();
                 d > cur && !max_disp.compare_exchange_weak(cur, d);) {
            }
        }
    };
    mutable std::array<MoveStats, size_t(VClass::Count)> m_move_stats;
 
    struct WallMoveStats
    {
        std::atomic<int> accepted{0}; 
        std::atomic<int> accepted_zero_tracked{0}; 
        std::atomic<int> added_deviation{0}; 
        std::atomic<int> added_dev_zero_tracked{0}; 
        std::atomic<double> max_single_step{0.}; 
 
        void note(bool zero_tracked, double dev_before, double dev_after)
        {
            ++accepted;
            if (zero_tracked) ++accepted_zero_tracked;
            const double d = dev_after - dev_before;
            if (d > 0.) {
                ++added_deviation;
                if (zero_tracked) ++added_dev_zero_tracked;
                for (double cur = max_single_step.load();
                     d > cur && !max_single_step.compare_exchange_weak(cur, d);) {
                }
            }
        }
    };
    mutable WallMoveStats m_wall_moves;
 
    double wall_offplane_deviation(const size_t vid) const;
 
    void log_vertex_movement(const char* when) const;
 
 
 
 
    void init_offset_sizing_field();
 
    std::pair<double, double> compute_distance_deviation() const;
 
    bool face_is_offset_surface_live(const Tuple& f) const;
    void diag_offset_bands(const char* tag) const;
 
    bool cell_is_offset_band(const size_t tid) const { return m_tet_attribute[tid].label == 2; }
 
    bool cell_is_input_complex(const size_t tid) const { return m_tet_attribute[tid].label == 1; }
 
    void warn_if_offset_reaches_domain_boundary() const;
 
    mutable size_t m_worst_dist_vid = static_cast<size_t>(-1);
    void log_worst_dist_vertex() const;
 
    void log_worst_tet(const char* when) const;
 
    std::vector<std::array<double, 8>> optimization_metrics;
    std::vector<std::array<int, 3>> churn_counts;
    std::vector<std::array<int, 3>> op_counts;
    bool m_converged = false;
 
    std::atomic<int> iter_cnt_split_born{0};
    std::atomic<int> iter_cnt_recollapsed{0};
    std::atomic<int> iter_cnt_recollapsed_same_pass{0};
    std::atomic<int> iter_cnt_split{0};
    std::atomic<int> iter_cnt_collapse{0};
    std::atomic<int> 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<int> iter_cnt_split_offset_before{0};
    std::atomic<int> iter_cnt_split_offset_endpoints{0};
    std::atomic<int> iter_cnt_split_offset_base_reject{0};
    std::atomic<int> iter_cnt_split_offset_tried{0};
    std::atomic<int> iter_cnt_split_offset{0};
    std::atomic<int> iter_cnt_split_offset_reject{0};
 
 
 
 
    // The sizing field's bounds are user-configurable: see Parameters::min_sizing_scalar /
    // max_sizing_scalar / sizing_gradation. Refinement stops at min_sizing_scalar; coarsening
    // never exceeds max_sizing_scalar, since a sizing scalar above 1 would be coarser than the
    // target edge length the user asked for.
 
    static double mean_ratio_metric(const Vector3d& p0, const Vector3d& p1, const Vector3d& p2);
 
 
    void marching_tets();
 
 
 
    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_input_complex(const std::string& path); // write components labeled to be offset
    void write_vtu(const std::string& path);
    // void write_msh(const std::string& file);
    void write_msh_groups(const std::string& file);
 
private:
    struct EdgeSplitCache
    {
        size_t v1_id;
        size_t v2_id;
        // The marching-tets placement and the offset's own data for the new vertex. The
        // shared attributes (position, rounding, bbox, order) are written in `after`.
        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;
 
        // cache edge attributes
        std::map<simplex::Edge, EdgeAttributes> existing_e;
 
        // cache face attributes
        std::map<simplex::Face, FaceSnapshot> existing_f;
        int splitf_label;
 
        // cache tet attributes
        std::map<size_t, TetAttributes> tets;
    };
    wmtk::threading::enumerable_thread_specific<FaceSplitCache> face_split_cache;
 
    struct TetSplitCache
    {
        std::array<size_t, 4> v_ids;
 
        // cache retained edge attributes
        std::map<simplex::Edge, EdgeAttributes> existing_e;
 
        // cache retained face attributes
        std::map<simplex::Face, FaceSnapshot> existing_f;
 
        // cache tet attribute
        TetAttributes tet;
    };
    wmtk::threading::enumerable_thread_specific<TetSplitCache> tet_split_cache;
 
 
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
    bool any_tag_present(const CellTag& tag1, const CellTag& tag2)
    {
        if (tag2.empty()) {
            return tag1.empty();
        }
        if (tag1.empty()) { // tag1 is ambient, tag2 is not
            return false;
        }
 
        for (const int64_t& i : tag1) {
            if (tag2.find(i) != tag2.end()) {
                return true;
            }
        }
        return false;
    }
 
    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
    void compute_vertex_partition();
 
    size_t get_partition_id(const Tuple& loc) const
    {
        return m_vertex_attribute[loc.vid(*this)].partition_id;
    }
 
    std::vector<Tuple> get_face_adjacent_tets(const Tuple& t) const
    {
        std::vector<Tuple> adj_tets;
        auto tet_1 = t.switch_tetrahedron(*this);
        if (tet_1) {
            adj_tets.push_back(tet_1.value());
        }
        auto tet_2 = t.switch_face(*this).switch_tetrahedron(*this);
        if (tet_2) {
            adj_tets.push_back(tet_2.value());
        }
        auto tet_3 = t.switch_edge(*this).switch_face(*this).switch_tetrahedron(*this);
        if (tet_3) {
            adj_tets.push_back(tet_3.value());
        }
        auto tet_4 =
            t.switch_vertex(*this).switch_edge(*this).switch_face(*this).switch_tetrahedron(*this);
        if (tet_4) {
            adj_tets.push_back(tet_4.value());
        }
        return adj_tets;
    }
 
    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;
    }
 
    void reset_connected_components()
    {
        auto verts = get_vertices();
        for (const Tuple& v : verts) {
            size_t v_id = v.vid(*this);
            m_vertex_extra[v_id].component_id = 0;
        }
    }
};
 
 
} // namespace wmtk::components::topological_offset