TopoOffsetTetMesh.h

#pragma once
 
#include <array>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <functional>
#include <map>
#include <mutex>
#include <optional>
#include <set>
#include <string>
 
#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/enumerable_thread_specific.hpp>
#include "OffsetPotential.hpp"
#include "Parameters.h"
#include "SimplicialComplexBVH.hpp"
#include "TagEnvelopes.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_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
    Vector3d m_turn_start = Vector3d::Zero();
    bool m_turn_start_valid = 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
    bool on_sheet = false;
};
 
 
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 under sphere_trace_initialization: sphere tracing
                         // along the edge to d(x) = target_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;
    MatrixXd m_sheet_V;
    MatrixXi m_sheet_F;
    void classify_sheet_faces();
    std::shared_ptr<SimplicialComplexBVH> m_input_complex_bvh;
 
    std::shared_ptr<OffsetPotential3D> m_offset_potential;
 
    std::shared_ptr<SampleEnvelope> m_input_complex_envelope;
 
    int m_n_regions = 0; 
    std::vector<std::shared_ptr<OffsetPotential3D>> m_region_potentials; 
    std::vector<std::shared_ptr<SimplicialComplexBVH>> m_region_bvhs;
    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;
 
    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;
 
    mutable std::map<uint64_t, std::shared_ptr<SampleEnvelope>> m_offset_isect_cache;
 
    std::shared_ptr<SampleEnvelope> containment_for(uint64_t region_mask, bool on_offset) const;
 
    bool project_into_containment(size_t vid, Vector3d& x) const;
 
    enum class OptPhase { A, B, Single };
 
    bool phase_places_front() const { return m_phase != OptPhase::A; }
 
    OptPhase m_phase = OptPhase::A;
 
    bool m_freeze_front = false;
 
    enum class EnvelopeSetup { PerTag, WallComplex };
    EnvelopeSetup envelope_setup() const
    {
        return m_offset_params.deform_others ? EnvelopeSetup::WallComplex : EnvelopeSetup::PerTag;
    }
    static constexpr int64_t m_wall_tag = -2; 
    static constexpr int64_t m_complex_tag = -3; 
    std::string envelope_key_name(int64_t tag) const;
    bool face_is_complex_boundary(const Tuple& f) const;
    void build_boundary_envelopes(const char* when, EnvelopeSetup setup);
 
    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_single_phase();
 
    double phase_b_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;
 
    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;
 
        // 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;
    }
 
    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 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; }
 
    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;
    bool edge_is_offset_surface_live(size_t a, size_t b) const;
    std::vector<std::array<size_t, 2>> offset_surface_edges() 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;
    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;
 
    std::vector<char> offset_surface_foldover_labels() 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;
    int m_ab_round = 0;
    mutable size_t m_debug_seq = 0;
    mutable std::vector<std::string> m_debug_frame_labels;
    mutable std::map<std::string, std::vector<size_t>> m_debug_pvd_series;
    void write_debug_pvd() const;
    mutable int m_debug_pass = 0;
    mutable int m_debug_last_round = -1;
    mutable char m_debug_last_phase = '?';
    double m_gradient_reference = 0.;
    bool m_converged = false;
    bool m_quality_converged = true;
    double m_quality_max_amips = 0.;
 
    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, 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<int> iter_cnt_collapse_guard_reject{0};
    std::atomic<int> iter_cnt_swap_guard_reject{0};
 
    static constexpr long long kFlipTraceEvery = 20000;
    std::atomic<long long> flip_trace_n{0}; 
    std::atomic<long long> flip_trace_nonmono{
        0}; 
    std::atomic<long long> flip_trace_bar{0}; 
    std::array<std::atomic<long long>, 8> flip_trace_dec{};
    void flip_trace_record(double before, double after);
 
    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_quality{0};
    mutable std::atomic<long long> funnel_quality_ok{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;
 
    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();
    }
    bool face_is_on_domain_boundary(const size_t fid) const
    {
        return m_face_attribute[fid].m_is_bbox_fs >= 0;
    }
 
    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> before_phase_b_not_offset{
            0}; 
        std::atomic<int> before_phase_b_enveloped_background{0}; 
        std::atomic<int> before_phase_b_enveloped_offset{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}; 
        std::atomic<long long> res_before_nano{0};
        std::atomic<long long> res_after_nano{0};
        std::atomic<long long> res_max_before_nano{0};
        std::atomic<long long> res_max_after_nano{0};
 
        void reset()
        {
            for (std::atomic<int>* c :
                 {&attempted,
                  &before_bbox,
                  &before_unrounded,
                  &before_phase_b_not_offset,
                  &before_phase_b_enveloped_background,
                  &before_phase_b_enveloped_offset,
                  &offset_attempted,
                  &offset_accepted,
                  &interior_attempted,
                  &region_attempted}) {
                c->store(0);
            }
            for (std::atomic<long long>* c :
                 {&res_before_nano, &res_after_nano, &res_max_before_nano, &res_max_after_nano}) {
                c->store(0);
            }
        }
    };
    SmoothTrace m_smooth_trace;
 
    optimization::NewtonCounters m_newton_front;
    optimization::NewtonCounters m_newton_plastic;
    void log_smoothing_pass_accounting() 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;
    mutable wmtk::threading::enumerable_thread_specific<std::vector<size_t>> m_collapse_edge_link;
    double face_resolution_or_inf(size_t a, size_t b, size_t c) const;
    std::pair<double, double> ops_guard_measures(
        const std::vector<std::array<size_t, 3>>& before,
        const std::vector<std::array<size_t, 3>>& after,
        size_t removed) const;
    double vertex_ring_measure(const std::vector<std::array<size_t, 3>>& ring) const;
    double ring_face_area(size_t a, size_t b, size_t c) const;
    bool ops_guard_refuses_collapse(size_t v1, size_t v2) const;
    std::array<size_t, 3> face_vids(const Tuple& f) const;
    void log_smooth_trace() const;
 
    void log_region_face_mask_health(const std::string& when) const;
 
    void audit_surface_containment(const std::string& when) const;
 
    mutable std::atomic<int> m_placement_env_entry_outside{0};
    mutable std::atomic<int> m_placement_projected{0};
    mutable std::atomic<int> m_placement_tangential{0};
 
 
    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]);
    }
 
    uint64_t face_boundary_bits(const Tuple& f) const
    {
        const std::optional<Tuple> opp = f.switch_tetrahedron(*this);
        if (!opp) {
            return tag_bits(m_tet_attribute[f.tid(*this)].tag); // domain wall
        }
        const auto& t0 = m_tet_attribute[f.tid(*this)].tag;
        const auto& t1 = m_tet_attribute[opp->tid(*this)].tag;
        CellTag diff;
        std::set_symmetric_difference(
            t0.begin(),
            t0.end(),
            t1.begin(),
            t1.end(),
            std::inserter(diff, diff.begin()));
        return tag_bits(diff);
    }
 
    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
    {
        uint64_t mask = face_mask(vids);
        bool all_offset = true;
        for (const size_t v : vids) {
            all_offset = all_offset && m_vertex_extra[v].m_is_on_offset;
        }
        // The ambiguous case: all corners can be on region boundaries and on the offset surface
        // at once. The corner-mask AND is then necessary but not sufficient for the face lying
        // on a shared boundary; ask the face's own class, the only record that distinguishes a
        // chord from a boundary. Reading the slot is safe here: a split child never reaches
        // this branch, and the m_is_surface_fs guard leaves both constraints standing when a
        // slot is illegible.
        if (mask != 0 && all_offset) {
            if (const auto found = try_tuple_from_face(vids)) {
                const size_t fid = std::get<1>(*found);
                if (m_face_attribute[fid].m_is_surface_fs) {
                    if (face_is_offset(fid)) {
                        mask = 0; // an offset face lies on no region boundary
                    } else {
                        all_offset = false; // a region face is not the offset surface
                    }
                }
            }
        }
        const std::shared_ptr<SampleEnvelope> base = containment_for(mask, all_offset);
        if (base || m_deform_tags.empty()) return base;
        // deform_others' ops-only tube: a released boundary is held by no mask -- its vertices
        // were freed so smoothing can carry the object -- which would leave the operations free
        // to decimate and reposition it. A face the masks and the offset class do not claim, but
        // which lies on a released boundary by its incident tets' current tags, is held to the
        // tube around the boundary's current shape.
        if (const auto found = try_tuple_from_face(vids)) {
            if (face_borders_released_boundary(std::get<0>(*found))) return released_envelope();
        }
        return nullptr;
    }
 
    bool allow_surface_swap() const override { return true; }
 
    bool swap_quality_allowed(const double after, const double before, const bool is_surface_flip)
        const override
    {
        const bool ok = after < before;
        if (is_surface_flip && m_swap_sides.local().worthwhile) {
            ++funnel_quality;
            if (ok) ++funnel_quality_ok;
        }
        return ok;
    }
    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 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) {
            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
    {
        return containment_for(vertex_boundary_mask(vid), m_vertex_extra[vid].m_is_on_offset);
    }
 
    bool smooth_front_vertex_phase_b(const Tuple& t);
    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;
    bool front_vertex_alignment_traps_1d_solve(size_t vid) const;
    double front_vertex_conv_ratio(size_t vid) const;
    bool front_vertex_placed(size_t vid) const;
    bool front_placed_by_ratio(const double ratio) const
    {
        return std::isfinite(ratio) && ratio <= 1.;
    }
    double edge_conv_ratio(size_t a, size_t b) const;
    double face_conv_ratio(size_t a, size_t b, size_t c) const;
    mutable size_t m_front_gradient_worst_vid =
        static_cast<size_t>(-1); 
    Vector3d front_vertex_normal(size_t vid) const;
    std::shared_ptr<polysolve::nonlinear::Problem> phase_b_front_objective(
        size_t vid,
        const Vector3d& x) const;
    std::shared_ptr<polysolve::nonlinear::Problem> smoothing_extra_energy(
        const size_t vid) const override
    {
        std::shared_ptr<polysolve::nonlinear::Problem> front;
        if (phase_places_front() && m_offset_potential && m_vertex_extra[vid].m_is_on_offset &&
            vertex_boundary_mask(vid) == 0) {
            front = phase_b_front_energy(vid, potential_ptr_for(vid));
        }
        const std::shared_ptr<polysolve::nonlinear::Problem> rest = rest_energy_for_vertex(vid);
        if (!front) return rest;
        if (!rest) return front;
        auto sum = std::make_shared<optimization::EnergySum>();
        sum->add_energy(front);
        sum->add_energy(rest);
        return sum;
    }
 
    // ------- deform_others: other input regions deform instead of being envelope-held -------
 
    std::set<int64_t> m_deform_tags;
    std::set<int64_t> m_source_tags;
    mutable std::shared_ptr<SampleEnvelope> m_released_envelope;
    mutable std::atomic<bool> m_released_tube_dirty{false};
    mutable std::mutex m_released_mutex;
    std::shared_ptr<SampleEnvelope> released_envelope() const;
    bool face_borders_released_boundary(const Tuple& f) const;
    bool cell_is_deformable(size_t tid) const;
    bool m_plastic_active = false; 
    bool cell_is_plastic(size_t tid) const
    {
        // Everything outside the band: ambient, the other objects and the input complex's
        // interior alike -- one material. The complex's boundary is what its tube holds.
        return m_plastic_active && m_tet_attribute[tid].label != 2;
    }
    void stamp_plastic_rests();
    bool smooth_plastic_vertex(const Tuple& t);
    bool cell_is_released_band(size_t tid) const;
    void stamp_rest_cell(size_t tid);
    void release_deformable_regions();
    std::shared_ptr<polysolve::nonlinear::Problem> rest_energy_for_vertex(size_t vid) const;
    std::shared_ptr<polysolve::nonlinear::Problem> phase_b_front_energy(
        size_t vid,
        const std::shared_ptr<const OffsetPotential3D>& pot) const;
 
    std::tuple<double, double> optimization_quality_stats() override;
 
    double optimization_stop_metric() const override
    {
        return m_phase != OptPhase::B ? wmtk::TetOptimizerMesh::optimization_stop_metric() : 1.;
    }
 
    int stencil_order() const { return m_offset_params.stencil_order; }
    int stencil_points_per_face() const
    {
        const int k = m_offset_params.stencil_order;
        if (k < 0) return 0;
        if (k == 0) return 3;
        const int n = 1 << (k - 1);
        return (n + 1) * (n + 2) / 2 + n * n;
    }
 
    double offset_residual_tolerance() const
    {
        const double s = m_offset_potential ? m_offset_potential->level_set_slope() : 1.;
        return std::max(0.5 * offset_gradient_tolerance() / (s * s), 1e-16);
    }
 
    double offset_gradient_tolerance() const
    {
        return std::max(m_offset_params.front_conv_frac() * m_gradient_reference, 1e-16);
    }
 
    double gradient_reference() const { return m_gradient_reference; }
 
    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
    {
        const int k = m_offset_params.stencil_order;
        if (k < 0) return;
 
        const auto emit = [&](const double wa, const double wb, const double wc) {
            visit(Vector3d(wa * p0 + wb * p1 + wc * p2), wa, wb, wc);
        };
        // Order 0 is the three CORNERS alone. That is the whole point of including them: the
        // measure sampled here is a distance to the level set, which at a corner is exactly that
        // vertex's own placement error, so one stencil covers what used to be two criteria.
        if (k == 0) {
            emit(1., 0., 0.);
            emit(0., 1., 0.);
            emit(0., 0., 1.);
            return;
        }
        // Order k >= 1: the vertices of the triangle subdivided k-1 times by 4-way midpoint
        // refinement, plus the centroid of each of its 4^(k-1) sub-triangles. With n = 2^(k-1)
        // segments per side that is (n+1)(n+2)/2 + n^2 points: 4, 10, 31, 109, ...
        const int n = 1 << (k - 1);
        const double dn = double(n);
        for (int i = n; i >= 0; --i) {
            for (int j = n - i; j >= 0; --j) {
                emit(double(i) / dn, double(j) / dn, double(n - i - j) / dn);
            }
        }
        // The sub-triangles, in integer barycentric coordinates over 3n. "Up" triangles have
        // corners (i+1,j,l), (i,j+1,l), (i,j,l+1) for i+j+l = n-1, so centroid (3i+1, 3j+1,
        // 3l+1); "down" triangles (i+1,j+1,l), (i,j+1,l+1), (i+1,j,l+1) for i+j+l = n-2, so
        // centroid (3i+2, 3j+2, 3l+2). n(n+1)/2 + n(n-1)/2 = n^2 of them.
        const double d3n = 3. * dn;
        for (int i = n - 1; i >= 0; --i) {
            for (int j = n - 1 - i; j >= 0; --j) {
                const int l = n - 1 - i - j;
                emit((3. * i + 1.) / d3n, (3. * j + 1.) / d3n, (3. * l + 1.) / d3n);
            }
        }
        for (int i = n - 2; i >= 0; --i) {
            for (int j = n - 2 - i; j >= 0; --j) {
                const int l = n - 2 - i - j;
                emit((3. * i + 2.) / d3n, (3. * j + 2.) / d3n, (3. * l + 2.) / d3n);
            }
        }
    }
 
    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.;
        double max_in_face_pinned = 0.;
        double max_normal_aligned = 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;
 
    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);
        Vector3d worst_face_centroid = Vector3d::Zero();
        double worst_face_len = 0.; 
        size_t n_faces_over = 0, n_faces_over_placed = 0;
        double max_face_placed = 0.;
        Vector3d worst_placed_centroid = Vector3d::Zero();
        double tube = 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 "area-weighted 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 tube) 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);
 
    void grade_sizing(double grade, const std::vector<size_t>& seeds);
    size_t grade_sizing_by_distance(const std::vector<size_t>& seeds);
 
    struct SmoothingProgress
    {
        double front_max_ratio = 0.;
        size_t front_worst_vid = static_cast<size_t>(-1);
        size_t n_front = 0; 
        size_t n_front_unmeasurable = 0; 
        double front_max_step = 0.; 
        double background_max_step = 0.;
        size_t background_worst_vid = static_cast<size_t>(-1);
        size_t n_background = 0;
    };
    SmoothingProgress smoothing_progress(const std::vector<Vector3d>& before);
    void smooth_group_to_convergence(const char* group_name);
    std::optional<EnergyCriterion> m_energy_verdict;
    double edge_interpolation_residual(size_t a, size_t b) const;
 
    Vector3d offset_vertex_normal(const size_t vid) const;
 
    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_boundary_mask(vid) != 0) 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;
 
    mutable std::atomic<size_t> m_deg_split_created{0};
    mutable std::atomic<size_t> m_deg_collapse_offered{0};
    mutable std::atomic<size_t> m_deg_collapse_allowed{0};
    mutable std::atomic<size_t> m_deg_collapse_by_ringmax{0};
    mutable std::atomic<size_t> m_deg_collapse_by_stop{0};
    mutable std::atomic<size_t> m_deg_collapse_by_unrounded{0};
    std::array<size_t, 6> m_deg_prev_counts{{0, 0, 0, 0, 0, 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 false; }
 
    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;
 
    void write_optimization_debug_output(const std::string& path) override
    {
        const char ph = (m_phase == OptPhase::A) ? 'A' : (m_phase == OptPhase::B ? 'B' : 'S');
        if (m_ab_round != m_debug_last_round || ph != m_debug_last_phase) {
            m_debug_last_round = m_ab_round;
            m_debug_last_phase = ph;
            m_debug_pass = 0;
        }
        std::string label = path;
        if (path.rfind("debug_", 0) == 0) {
            label = fmt::format(
                "r{}{}{}{}",
                m_ab_round,
                ph,
                ++m_debug_pass,
                m_debug_pass_name.empty() ? std::string() : "_" + m_debug_pass_name);
        } else if (path.rfind("phase_", 0) == 0) {
            label = fmt::format("r{}{}_end", m_ab_round, ph);
        }
        write_debug_frame(label);
    }
    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 = "");
 
    bool ambient_assert();
 
    void label_input_complex();
 
    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 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 execute_offset(const std::filesystem::path& output_file);
 
    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_phi_grid(const std::string& path, int n) const;
 
    void write_input_complex(const std::string& path);
    void write_vtu(const std::string& path);
    void write_msh_groups(const std::string& file);
 
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 sag_measured = false;
        double sag_before = 0.0;
        double sag_after = 0.0;
        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);
    mutable wmtk::threading::enumerable_thread_specific<SwapSurfaceSides> m_swap_sides;
 
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<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