TetMesh.h

#pragma once
 
#include <wmtk/utils/VectorUtils.h>
#include <wmtk/AttributeCollection.hpp>
#include <wmtk/OperationDryRun.hpp>
#include <wmtk/SlotPool.hpp>
#include <wmtk/Types.hpp>
#include <wmtk/simplex/Simplex.hpp>
#include <wmtk/simplex/SimplexCollection.hpp>
#include <wmtk/threading/enumerable_thread_specific.hpp>
#include <wmtk/threading/vertex_mutex.hpp>
#include <wmtk/utils/Logger.hpp>
 
#include <Eigen/Core>
 
#include <array>
#include <atomic>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <functional>
#include <limits>
#include <map>
#include <memory>
#include <optional>
#include <vector>
 
#include <wmtk/history/Fwd.hpp>
 
namespace wmtk {
class TetMesh
{
protected:
    static constexpr std::array<std::array<int, 2>, 6> m_local_edges = {
        {{{0, 1}}, {{1, 2}}, {{0, 2}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
 
private:
    static constexpr std::array<int, 4> m_map_vertex2edge = {{0, 0, 1, 3}};
    static constexpr std::array<int, 4> m_map_vertex2oppo_face = {{3, 1, 2, 0}};
    static constexpr std::array<int, 6> m_map_edge2face = {{0, 0, 0, 1, 2, 1}};
    static constexpr std::array<std::array<int, 3>, 4> m_local_faces = {
        {{{0, 1, 2}}, {{0, 2, 3}}, {{0, 1, 3}}, {{1, 2, 3}}}}; // sorted local vids
    static constexpr std::array<std::array<int, 3>, 4> m_local_edges_in_a_face = {
        {{{0, 1, 2}}, {{2, 5, 3}}, {{3, 4, 0}}, {{5, 1, 4}}}};
 
public:
    // Cell Tuple Navigator
    class Tuple
    {
        size_t m_global_vid = std::numeric_limits<size_t>::max();
        size_t m_local_eid = std::numeric_limits<size_t>::max();
        size_t m_local_fid = std::numeric_limits<size_t>::max();
        size_t m_global_tid = std::numeric_limits<size_t>::max();
 
        size_t m_hash = 0;
 
    private:
        Tuple(const TetMesh& m, size_t vid, size_t local_eid, size_t local_fid, size_t tid);
 
    public:
        Tuple() {}
 
        friend TetMesh;
 
        bool is_valid(const TetMesh& m) const;
        bool is_valid_unlocked(const TetMesh& m) const { return is_valid(m); }
        bool is_boundary_edge(const TetMesh& m) const;
        bool is_boundary_face(const TetMesh& m) const;
 
        void print_info() const;
 
        void print_info(const TetMesh& m) const;
 
        size_t vid(const TetMesh& m) const;
 
        size_t eid(const TetMesh& m) const;
 
        size_t fid(const TetMesh& m) const;
 
        size_t tid(const TetMesh& m) const;
 
        Tuple switch_vertex(const TetMesh& m) const;
        Tuple switch_edge(const TetMesh& m) const;
        Tuple switch_face(const TetMesh& m) const;
 
        std::optional<Tuple> switch_tetrahedron(const TetMesh& m) const;
        std::optional<Tuple> switch_tetrahedron_slow(const TetMesh& m) const;
 
 
 
        void check_validity(const TetMesh& m) const;
        friend bool operator==(const Tuple& a, const Tuple& t)
        {
            return (
                std::tie(a.m_global_vid, a.m_local_eid, a.m_local_fid, a.m_global_tid, a.m_hash) ==
                std::tie(t.m_global_vid, t.m_local_eid, t.m_local_fid, t.m_global_tid, t.m_hash));
        }
        friend bool operator<(const Tuple& a, const Tuple& t)
        {
            return (
                std::tie(a.m_global_vid, a.m_local_eid, a.m_local_fid, a.m_global_tid, a.m_hash) <
                std::tie(t.m_global_vid, t.m_local_eid, t.m_local_fid, t.m_global_tid, t.m_hash));
        }
    };
 
    class SmartTuple
    {
        const TetMesh& m_mesh;
        Tuple m_tuple;
 
    public:
        SmartTuple(const TetMesh& mesh, const Tuple& t)
            : m_mesh(mesh)
            , m_tuple(t)
        {}
 
        const Tuple& tuple() { return m_tuple; }
        const TetMesh& mesh() { return m_mesh; }
 
        SmartTuple& operator=(const SmartTuple& t)
        {
            m_tuple = t.m_tuple;
            return *this;
        }
 
        bool is_valid() const { return m_tuple.is_valid(m_mesh); }
        bool is_boundary_edge() const { return m_tuple.is_boundary_edge(m_mesh); }
        bool is_boundary_face() const { return m_tuple.is_boundary_face(m_mesh); }
        size_t vid() const { return m_tuple.vid(m_mesh); }
        size_t eid() const { return m_tuple.eid(m_mesh); }
        size_t fid() const { return m_tuple.fid(m_mesh); }
        size_t tid() const { return m_tuple.tid(m_mesh); }
        SmartTuple switch_vertex() const { return {m_mesh, m_tuple.switch_vertex(m_mesh)}; }
        SmartTuple switch_edge() const { return {m_mesh, m_tuple.switch_edge(m_mesh)}; }
        SmartTuple switch_face() const { return {m_mesh, m_tuple.switch_face(m_mesh)}; }
        std::optional<SmartTuple> switch_tetrahedron() const
        {
            const std::optional<Tuple> t = m_tuple.switch_tetrahedron(m_mesh);
            if (t) {
                return std::optional<SmartTuple>({m_mesh, t.value()});
            }
            return {};
        }
        void check_validity() const { return m_tuple.check_validity(m_mesh); }
    };
 
    class VertexConnectivity
    {
    public:
        std::vector<size_t> m_conn_tets; // todo: always keep it sorted
        bool m_is_removed = false;
 
        size_t& operator[](const size_t index)
        {
            assert(index < m_conn_tets.size());
            return m_conn_tets[index];
        }
 
        size_t operator[](const size_t index) const
        {
            assert(index < m_conn_tets.size());
            return m_conn_tets[index];
        }
 
        friend bool operator==(const VertexConnectivity& l, const VertexConnectivity& r)
        {
            return std::tie(l.m_conn_tets, l.m_is_removed) ==
                   std::tie(r.m_conn_tets, r.m_is_removed); // keep the same order
        }
 
        void print_info() {}
    };
 
    class TetrahedronConnectivity
    {
    public:
        std::array<size_t, 4> m_indices;
        bool m_is_removed = false;
 
        uint32_t hash = 0;
 
        size_t& operator[](size_t index)
        {
            assert(index < 4);
            return m_indices[index];
        }
 
        size_t operator[](size_t index) const
        {
            assert(index < 4);
            return m_indices[index];
        }
 
        int find(size_t v_id) const
        {
            for (int j = 0; j < 4; j++) {
                if (v_id == m_indices[j]) return j;
            }
            return -1;
        }
 
        int find_local_edge(size_t v1_id, size_t v2_id) const
        {
            std::array<int, 2> e;
            for (int j = 0; j < 4; j++) {
                if (v1_id == m_indices[j])
                    e[0] = j;
                else if (v2_id == m_indices[j])
                    e[1] = j;
            }
            if (e[0] > e[1]) std::swap(e[0], e[1]);
            int i =
                std::find(m_local_edges.begin(), m_local_edges.end(), e) - m_local_edges.begin();
            if (i >= m_local_edges.size()) return -1;
            return i;
        }
 
        int find_local_face(size_t v1_id, size_t v2_id, size_t v3_id) const
        {
            std::array<int, 3> f;
            for (int j = 0; j < 4; j++) {
                if (v1_id == m_indices[j])
                    f[0] = j;
                else if (v2_id == m_indices[j])
                    f[1] = j;
                else if (v3_id == m_indices[j])
                    f[2] = j;
            }
            std::sort(f.begin(), f.end());
            int i =
                std::find(m_local_faces.begin(), m_local_faces.end(), f) - m_local_faces.begin();
            if (i >= m_local_edges.size()) return -1;
            return i;
        }
 
        friend bool operator==(const TetrahedronConnectivity& l, const TetrahedronConnectivity& r)
        {
            return std::tie(l.m_indices, l.m_is_removed, l.hash) ==
                   std::tie(r.m_indices, r.m_is_removed, r.hash); // keep the same order
        }
 
        void print_info() {}
    };
 
    TetMesh();
    virtual ~TetMesh() = default;
    size_t vert_capacity() const { return m_vertex_connectivity.live(); }
    size_t tet_capacity() const { return m_tet_connectivity.live(); }
 
    static constexpr int EDGES_PER_CELL = 6;
    static constexpr int FACES_PER_CELL = 4;
    size_t cell_capacity() const { return tet_capacity(); }
    size_t cell_storage_capacity() const { return m_tet_connectivity.capacity(); }
    Tuple tuple_from_cell(size_t cid) const { return tuple_from_tet(cid); }
    size_t vertex_size() const
    {
        int cnt = 0;
        for (auto i = 0; i < vert_capacity(); i++) {
            if (!m_vertex_connectivity[i].m_is_removed) cnt++;
        }
        return cnt;
    }
    size_t tet_size() const
    {
        int cnt = 0;
        for (auto i = 0; i < tet_capacity(); i++) {
            if (!m_tet_connectivity[i].m_is_removed) cnt++;
        }
        return cnt;
    }
    void init(size_t n_vertices, const std::vector<std::array<size_t, 4>>& tets);
    void init_with_isolated_vertices(
        size_t n_vertices,
        const std::vector<std::array<size_t, 4>>& tets);
 
    void init(const MatrixXi& T);
 
    bool split_edge(const Tuple& t, std::vector<Tuple>& new_tets);
    virtual bool collapse_edge(const Tuple& t, std::vector<Tuple>& new_tets);
 
    bool link_condition(const Tuple& t);
 
    bool collapse_edge_conn(
        const Tuple& loc0,
        size_t& v1_id,
        Tuple& new_loc,
        std::map<size_t, wmtk::TetMesh::VertexConnectivity>& rollback_vert_conn,
        std::vector<size_t>& n1_t_ids_copy,
        std::vector<size_t>& new_tet_id,
        std::vector<TetrahedronConnectivity>& old_tets);
 
    bool collapse_edge_check_topology(const std::vector<size_t>& new_tet_id);
 
    void collapse_edge_rollback(
        size_t& v1_id,
        std::map<size_t, wmtk::TetMesh::VertexConnectivity>& rollback_vert_conn,
        std::vector<size_t>& n1_t_ids,
        std::vector<size_t>& new_tet_id,
        std::vector<TetrahedronConnectivity>& old_tets);
 
    bool swap_edge_56(const Tuple& t, std::vector<Tuple>& new_tets);
    bool swap_edge_44(const Tuple& t, std::vector<Tuple>& new_tets);
    bool swap_edge(const Tuple& t, std::vector<Tuple>& new_tets);
    bool swap_face(const Tuple& t, std::vector<Tuple>& new_tets);
    bool smooth_vertex(const Tuple& t);
    bool m_check_smoothing_stays_in_star = false;
 
    bool split_tet(const Tuple& t, std::vector<Tuple>& new_tets);
 
    bool split_face(const Tuple& t, std::vector<Tuple>& new_tets);
 
    void triangle_insertion(
        const std::vector<Tuple>& intersected_tets,
        const std::vector<Tuple>& intersected_edges,
        std::vector<size_t>& new_edge_vids,
        std::vector<size_t>& new_center_vids,
        std::vector<std::array<size_t, 4>>& center_split_tets);
 
    bool insert_point(const Tuple& t, std::vector<Tuple>& new_tets);
    virtual bool insert_point_before(const Tuple& t) { return true; };
    virtual bool insert_point_after(std::vector<Tuple>& new_tets) { return true; };
    void consolidate_mesh();
 
    std::shared_ptr<history::TrackedSurfaceRecorder> start_surface_history(
        std::function<Eigen::Vector3d(size_t)> position);
    void stop_surface_history();
    size_t resync_surface_history();
    std::shared_ptr<history::VolumeRecorder> start_volume_history(
        std::function<Eigen::Vector3d(size_t)> position);
    void stop_volume_history();
    virtual std::vector<int64_t> tet_tags(size_t tid) const { return {}; }
    std::vector<std::array<size_t, 3>> tracked_faces() const;
    virtual bool is_tracked_face(size_t tid, int j) const { return false; }
    virtual std::vector<std::array<size_t, 2>> tracked_feature_edges() const { return {}; }
    virtual std::array<size_t, 3> orient_tracked_face(
        const std::array<size_t, 3>& f,
        const std::vector<size_t>& tets) const
    {
        return f;
    }
    std::array<size_t, 3> orient_tracked_face(const std::array<size_t, 3>& f) const
    {
        return orient_tracked_face(f, tets_of_face(f));
    }
    std::array<size_t, 3> face_out_of(size_t tid, const std::array<size_t, 3>& face) const;
 
protected:
    std::vector<std::array<size_t, 3>> tracked_faces_around(const std::vector<size_t>& vs, bool all)
        const;
    struct VertexFace
    {
        std::array<size_t, 3> face;
        size_t fid;
        std::vector<size_t> tets;
    };
    std::vector<VertexFace> faces_at_vertex(size_t v) const;
    size_t face_id_in(const std::array<size_t, 3>& f, size_t tid) const;
 
    virtual bool checks_tracked_surface() const { return false; }
    virtual void tracked_surface_before(const std::vector<size_t>& vs, bool all) {}
    virtual bool tracked_surface_after(const std::vector<size_t>& vs, bool all) { return true; }
 
public:
    std::vector<size_t> tets_of_face(const std::array<size_t, 3>& f) const;
 
    std::vector<Tuple> get_edges() const;
    std::vector<Tuple> get_faces() const;
    std::vector<Tuple> get_vertices() const;
    std::vector<Tuple> get_tets() const;
    // virtual void for_each_edge(const std::function<void(const TetMesh::Tuple&)>&);
 
    // TODO: make this concurrent
    virtual void for_each_face(const std::function<void(const TetMesh::Tuple&)>&);
    // /**
    //  * @brief looping through all the unique vertices and perform the given function
    //  *
    //  */
    // virtual void for_each_vertex(const std::function<void(const TetMesh::Tuple&)>&);
    // /**
    //  * @brief looping through all the unique tet and perform the given function
    //  *
    //  */
    // virtual void for_each_tetra(const std::function<void(const TetMesh::Tuple&)>&);
 
public:
    template <typename T>
    using vector = std::vector<T>;
 
public:
    AbstractAttributeContainer* p_vertex_attrs = nullptr;
    AbstractAttributeContainer* p_edge_attrs = nullptr;
    AbstractAttributeContainer* p_face_attrs = nullptr;
    AbstractAttributeContainer* p_tet_attrs = nullptr;
    // AbstractAttributeContainer vertex_attrs, edge_attrs, face_attrs, tet_attrs;
 
 
public:
    void set_preallocation_factor(double factor)
    {
        if (factor >= 1.0) m_preallocation_factor = factor;
    }
    double preallocation_factor() const { return m_preallocation_factor; }
 
    // Atomically reserve `n` contiguous fresh tet/vertex slots. Returns the first
    // index of the block, or -1 if that would exceed the preallocated capacity (the
    // caller must then abort the operation before mutating any connectivity).
 
    size_t tet_storage_capacity() const { return m_tet_connectivity.capacity(); }
    size_t vert_storage_capacity() const { return m_vertex_connectivity.capacity(); }
    bool slots_exhausted() const
    {
        return m_vertex_connectivity.refused() || m_tet_connectivity.refused();
    }
    void clear_slots_exhausted()
    {
        m_vertex_connectivity.clear_refused();
        m_tet_connectivity.clear_refused();
    }
 
    size_t request_tet_slots(size_t n);
    size_t request_vert_slots(size_t n);
 
    size_t cell_slot_bound(const Tuple& t) const;
    void reserve_free_slots(size_t cells, size_t verts);
    virtual bool storage_grows_on_demand() const { return false; }
 
    struct SlotUsage
    {
        size_t cells = 0;
        size_t verts = 0;
    };
    static SlotUsage slot_usage_of_this_thread()
    {
        return {
            SlotPool<TetrahedronConnectivity>::taken_by_this_thread(),
            SlotPool<VertexConnectivity>::taken_by_this_thread()};
    }
protected:
    virtual void on_slot_storage_grown() {}
 
public:
    // Construction/insertion helpers (single-threaded ONLY): guarantee at least
    // `extra` free tet/vertex slots beyond the current used count, growing the
    // storage geometrically if necessary. Unlike request_*_slots these never fail;
    // they are used by the (legacy, non-concurrent) triangle-insertion path where
    // the mesh is being built rather than edited by concurrent operations.
    void ensure_free_tet_capacity(size_t extra)
    {
        const size_t need = m_tet_connectivity.live() + extra;
        if (m_tet_connectivity.capacity() >= need) return;
        const size_t newcap = std::max(need, m_tet_connectivity.capacity() * 2 + 1);
        m_tet_connectivity.resize(newcap);
        if (p_tet_attrs) p_tet_attrs->resize(newcap);
        if (p_face_attrs) p_face_attrs->resize(4 * newcap);
        if (p_edge_attrs) p_edge_attrs->resize(6 * newcap);
    }
    void ensure_free_vert_capacity(size_t extra)
    {
        const size_t need = m_vertex_connectivity.live() + extra;
        if (m_vertex_connectivity.capacity() >= need) return;
        const size_t newcap = std::max(need, m_vertex_connectivity.capacity() * 2 + 1);
        m_vertex_connectivity.resize(newcap);
        resize_vertex_mutex(newcap);
        if (p_vertex_attrs) p_vertex_attrs->resize(newcap);
    }
 
private:
    // capacity to reserve for a live element count
    size_t reserved_capacity(size_t live_count) const
    {
        const size_t floor = 64;
        // A mesh that grows on demand only needs room for the passes not to grow at once.
        const double factor = storage_grows_on_demand() ? 1.25 : m_preallocation_factor;
        double c = std::ceil(factor * static_cast<double>(live_count));
        size_t capacity = static_cast<size_t>(c);
        if (capacity < live_count) capacity = live_count;
        return capacity < floor ? floor : capacity;
    }
 
    // Stores the connectivity of the mesh
    SlotPool<VertexConnectivity> m_vertex_connectivity;
    SlotPool<TetrahedronConnectivity> m_tet_connectivity;
 
    double m_preallocation_factor = 6.0;
 
    int m_t_empty_slot = 0;
    int m_v_empty_slot = 0;
    size_t get_next_empty_slot_t();
    size_t get_next_empty_slot_v();
 
    // TODO: subdivide_tets function should not be in the TetMesh API.
    void subdivide_tets(
        const std::vector<size_t> t_ids,
        const std::vector<bool>& mark_surface,
        const std::map<std::array<size_t, 2>, size_t>& map_edge2vid,
        std::map<std::array<size_t, 3>, std::vector<std::array<size_t, 5>>>& new_face_vids,
        const std::vector<size_t>& new_vids,
        std::vector<size_t>& new_tids,
        std::vector<size_t>& new_center_vids,
        std::vector<std::array<size_t, 4>>& center_split_tets);
    void subdivide_a_tet(
        size_t t_id,
        const std::array<int, 6>& new_v_ids,
        bool mark_surface,
        std::map<std::array<size_t, 3>, std::vector<std::array<size_t, 5>>>& new_face_vids,
        std::vector<size_t>& new_tids,
        std::vector<size_t>& new_center_vids,
        std::vector<std::array<size_t, 4>>& center_split_tets);
 
public:
    virtual bool invariants(const std::vector<Tuple>&) { return true; }
 
    // ---- OPERATION ACCOUNTING (counting only) ----
    // Every call into split_edge, collapse_edge, swap_edge, swap_edge_44, swap_edge_56 and
    // swap_face reports `attempt` on entry and exactly ONE other event on the way out: `committed`,
    // or the place that refused it. The application's hooks are two of those places (before_hook,
    // after_hook); why a hook said no is the hook's own to count. So for each kind, attempt equals
    // the sum of the other events, and a subclass that counts can check both that identity and a
    // hook's refusals against the reasons the hook counted. Added because the refusals made here,
    // below every hook, were visible only as a difference: on the deliverable cube at
    // target_distance_rel 1e-2 / front_conv_rel 1e-4, 521142 of the 542123 edge swaps that
    // passed every before-hook in turn 5 were refused inside this class with no reason recorded.
    //
    // The default does nothing and no site reads a count: the mesh behaves identically with or
    // without an override.
    enum class OpKind : int { split, collapse, swap_32, swap_44, swap_56, swap_face, COUNT };
    enum class OpEvent : int {
        attempt,
        before_hook, // the application's *_before said no
        link_condition, // collapse: the classical link condition (m_collapse_check_link_condition)
        no_return_tet, // collapse: both probes for the return tuple hit a boundary face
        duplicate_tet, // a new tet repeats an existing one (collapse manifold check, 4-4, 5-6)
        topology, // collapse: m_collapse_check_topology found a face missing
        valence, // edge swap: the edge does not have this swap's number of tets around it
        boundary, // edge swap: the edge's ring of tets is open; face swap: a boundary face
        exists, // 3-2: the new face already exists; face swap: the new edge already exists
        no_case_allowed, // 4-4 / 5-6: every candidate was vetoed by *_accept_case
        no_better_case, // 4-4 / 5-6: no allowed candidate scored below the current cells
        out_of_slots, // the preallocated vertex or tet slots ran out
        after_hook, // the application's *_after said no, so the operation was rolled back
        invariants, // invariants() said no after the after-hook passed; rolled back
        committed,
        screened, // a dry run passed every check that precedes the change (see dry_run())
        COUNT
    };
    virtual void op_event(OpKind, OpEvent) const {}
    bool op_refused(const OpKind k, const OpEvent e) const
    {
        op_event(k, e);
        return false;
    }
 
    static bool dry_run() { return operation_dry_run(); }
    bool op_screened(const OpKind k) const
    {
        op_event(k, OpEvent::screened);
        return true;
    }
 
protected:
    virtual bool triangle_insertion_before(const std::vector<Tuple>& faces) { return true; }
    virtual bool triangle_insertion_after(const std::vector<std::vector<Tuple>>&) { return true; }
 
    // Checks if the split should be performed or not (user controlled)
    virtual bool split_edge_before(const Tuple& t) { return true; } // check edge condition
    // This function computes the attributes for the added simplices
    // if it returns false then the operation is undone
    virtual bool split_edge_after(const Tuple& t) { return true; } // check tet condition
 
    // Checks if the collapse should be performed or not (user controlled)
    virtual bool collapse_edge_before(const Tuple& t) { return true; }
    // If it returns false then the operation is undone (the tuple indexes a vertex and tet that
    // survived)
    virtual bool collapse_edge_after(const Tuple& t) { return true; }
    virtual bool swap_edge_44_before(const Tuple& t) { return true; }
    virtual double swap_edge_44_energy(
        const std::vector<std::array<size_t, 4>>& tets,
        const int op_case)
    {
        return -op_case;
    }
    virtual bool swap_edge_44_after(const Tuple& t) { return true; }
    virtual bool swap_edge_44_accept_case(const std::array<size_t, 2>& new_edge) { return true; }
    virtual bool swap_edge_56_before(const Tuple& t) { return true; }
    virtual double swap_edge_56_energy(
        const std::vector<std::array<size_t, 4>>& tets,
        const int op_case)
    {
        return -op_case;
    }
    virtual bool swap_edge_56_after(const Tuple& t) { return true; }
    virtual bool swap_edge_56_accept_case(const std::array<size_t, 3>& new_face) { return true; }
    virtual bool swap_edge_before(const Tuple& t) { return true; }
    virtual bool swap_edge_after(const Tuple& t) { return true; }
    virtual bool swap_face_before(const Tuple& t) { return true; }
    virtual bool swap_face_after(const Tuple& t) { return true; }
    virtual bool smooth_before(const Tuple& t) { return true; }
    virtual bool smooth_after(const Tuple& t) { return true; }
 
    virtual bool split_face_before(const Tuple& t) { return true; }
 
    virtual bool split_face_after(const Tuple& t) { return true; }
 
    virtual bool split_tet_before(const Tuple& t) { return true; }
 
    virtual bool split_tet_after(const Tuple& t) { return true; }
 
    // virtual void resize_vertex_mutex(size_t v) {}
 
public:
    Tuple tuple_from_edge(size_t tid, int local_eid) const;
    Tuple tuple_from_edge(const std::array<size_t, 2>& vids) const;
 
    Tuple tuple_from_face(size_t tid, int local_fid) const;
 
    std::tuple<Tuple, size_t> tuple_from_face(const std::array<size_t, 3>& vids) const;
    std::optional<std::tuple<Tuple, size_t>> try_tuple_from_face(
        const std::array<size_t, 3>& vids) const;
 
    size_t lowest_common_tet(size_t v0_id, size_t v1_id, size_t v2_id) const;
 
    size_t vertex_valence(const size_t vid) const
    {
        return m_vertex_connectivity[vid].m_conn_tets.size();
    }
    std::tuple<Tuple, size_t> tuple_from_face(const simplex::Face& f) const;
 
    Tuple tuple_from_vertex(size_t vid) const;
 
    Tuple tuple_from_tet(size_t tid) const;
 
    Tuple tuple_from_vids(size_t vid0, size_t vid1, size_t vid2, size_t vid3) const;
 
    simplex::Tet simplex_from_tet(const Tuple& t) const;
    simplex::Tet simplex_from_tet(const size_t tid) const;
    simplex::Face simplex_from_face(const Tuple& t) const;
    simplex::Edge simplex_from_edge(const Tuple& t) const;
 
 
    Tuple switch_vertex(const Tuple& t) const
    {
        auto loc = t.switch_vertex(*this);
        check_tuple_validity(loc);
        return loc;
    }
    Tuple switch_edge(const Tuple& t) const
    {
        auto loc = t.switch_edge(*this);
        check_tuple_validity(loc);
        return loc;
    }
    Tuple switch_face(const Tuple& t) const
    {
        auto loc = t.switch_face(*this);
        check_tuple_validity(loc);
        return loc;
    }
    std::optional<Tuple> switch_tetrahedron(const Tuple& t) const
    {
        auto loc = t.switch_tetrahedron(*this);
        if (loc.has_value()) check_tuple_validity(loc.value());
        return loc;
    }
 
    std::vector<Tuple> get_one_ring_tets_for_vertex(const Tuple& t) const;
    const std::vector<size_t>& get_one_ring_tids_for_vertex(const Tuple& t) const;
    const std::vector<size_t>& get_one_ring_tids_for_vertex(const size_t vid) const;
 
    std::vector<Tuple> get_one_ring_vertices_for_vertex(const Tuple& t) const;
    std::vector<size_t> get_one_ring_vids_for_vertex(size_t vid, std::vector<size_t>& cache);
    std::vector<size_t> get_one_ring_vids_for_vertex(size_t vid) const;
    std::vector<size_t> get_one_ring_vids_for_vertex_adj(size_t vid) const;
    std::vector<size_t> get_one_ring_vids_for_vertex_adj(size_t vid, std::vector<size_t>& cache);
 
    std::vector<Tuple> get_incident_tets_for_edge(const Tuple& t) const;
    std::vector<Tuple> get_incident_tets_for_edge(const size_t vid0, const size_t vid1) const;
 
    std::vector<size_t> get_incident_tids_for_edge(const Tuple& t) const;
    std::vector<size_t> get_incident_tids_for_edge(const size_t vid0, const size_t vid1) const;
 
    std::vector<Tuple> get_one_ring_tets_for_edge(const Tuple& t) const;
 
    std::vector<std::array<size_t, 3>> vertex_adjacent_boundary_faces(const Tuple& t) const;
    std::array<Tuple, 4> oriented_tet_vertices(const Tuple& t) const;
    std::array<size_t, 4> oriented_tet_vids(const Tuple& t) const;
    std::array<size_t, 4> oriented_tet_vids(const size_t tid) const;
    std::array<Tuple, 3> get_face_vertices(const Tuple& t) const;
    std::array<size_t, 3> get_face_vids(const Tuple& t) const;
 
    std::array<Tuple, 6> tet_edges(const Tuple& t) const;
    void check_tuple_validity(const Tuple& t) const { t.check_validity(*this); }
    bool check_mesh_connectivity_validity() const;
    void remove_tets_by_ids(const std::vector<size_t>& tids)
    {
        // A face or edge attribute lives at ONE canonical slot, (lowest incident tet) * 4
        // or 6 + local index. Deleting a tet that is a simplex's canonical owner silently
        // re-canonicalizes the simplex onto a surviving tet whose cell was never
        // maintained -- the operations only ever write the canonical slot -- so the
        // surviving cell holds stale garbage. Measured: after filtering, interior faces
        // read back m_is_surface_fs = true 0.12 (eps 0.09) from the real surface, and the
        // extracted "tracked surface" grows phantom triangles. Migrate every attribute
        // whose owner dies to the simplex's new owner BEFORE touching the connectivity.
        if (p_face_attrs != nullptr || p_edge_attrs != nullptr) {
            std::vector<bool> removed(tet_capacity(), false);
            for (const size_t tid : tids) {
                removed[tid] = true;
            }
            for (const size_t tid : tids) {
                const auto& tconn = m_tet_connectivity[tid];
                if (p_face_attrs != nullptr) {
                    for (int lf = 0; lf < 4; ++lf) {
                        std::array<size_t, 3> vids = {
                            {tconn[m_local_faces[lf][0]],
                             tconn[m_local_faces[lf][1]],
                             tconn[m_local_faces[lf][2]]}};
                        // The face's other incident tet, if any.
                        auto common = set_intersection(
                            m_vertex_connectivity[vids[0]].m_conn_tets,
                            m_vertex_connectivity[vids[1]].m_conn_tets);
                        common =
                            set_intersection(common, m_vertex_connectivity[vids[2]].m_conn_tets);
                        size_t other = std::numeric_limits<size_t>::max();
                        for (const size_t t : common) {
                            if (t != tid) {
                                other = t;
                            }
                        }
                        // Migrate only when this tet was the canonical owner (tid < other)
                        // and the other side survives; other < tid means the attribute
                        // already lives on the survivor.
                        if (other == std::numeric_limits<size_t>::max() || removed[other] ||
                            other < tid) {
                            continue;
                        }
                        std::array<size_t, 3> sorted = vids;
                        std::sort(sorted.begin(), sorted.end());
                        for (int lo = 0; lo < 4; ++lo) {
                            std::array<size_t, 3> ovids = {
                                {m_tet_connectivity[other][m_local_faces[lo][0]],
                                 m_tet_connectivity[other][m_local_faces[lo][1]],
                                 m_tet_connectivity[other][m_local_faces[lo][2]]}};
                            std::sort(ovids.begin(), ovids.end());
                            if (ovids == sorted) {
                                p_face_attrs->move(tid * 4 + lf, other * 4 + lo);
                                break;
                            }
                        }
                    }
                }
                if (p_edge_attrs != nullptr) {
                    for (int le = 0; le < 6; ++le) {
                        const size_t v0 = tconn[m_local_edges[le][0]];
                        const size_t v1 = tconn[m_local_edges[le][1]];
                        const auto fan = set_intersection(
                            m_vertex_connectivity[v0].m_conn_tets,
                            m_vertex_connectivity[v1].m_conn_tets);
                        // Only the canonical owner migrates, once, to the lowest survivor.
                        if (fan.empty() || fan.front() != tid) {
                            continue;
                        }
                        size_t new_owner = std::numeric_limits<size_t>::max();
                        for (const size_t t : fan) {
                            if (!removed[t]) {
                                new_owner = t; // fan is sorted ascending; first kept wins
                                break;
                            }
                        }
                        if (new_owner == std::numeric_limits<size_t>::max()) {
                            continue; // the edge dies with its whole fan
                        }
                        for (int lo = 0; lo < 6; ++lo) {
                            const size_t w0 = m_tet_connectivity[new_owner][m_local_edges[lo][0]];
                            const size_t w1 = m_tet_connectivity[new_owner][m_local_edges[lo][1]];
                            if ((w0 == v0 && w1 == v1) || (w0 == v1 && w1 == v0)) {
                                p_edge_attrs->move(tid * 6 + le, new_owner * 6 + lo);
                                break;
                            }
                        }
                    }
                }
            }
        }
 
        for (size_t tid : tids) {
            m_tet_connectivity[tid].m_is_removed = true;
            for (int j = 0; j < 4; j++)
                vector_erase(m_vertex_connectivity[m_tet_connectivity[tid][j]].m_conn_tets, tid);
        }
        for (auto& v : m_vertex_connectivity) {
            if (v.m_is_removed) continue;
            if (v.m_conn_tets.empty()) v.m_is_removed = true;
        }
        if (m_surface_history) resync_surface_history();
    }
    bool m_collapse_check_link_condition = true; // classical link condition
    bool m_collapse_check_topology = false; // sanity check
    bool m_collapse_check_manifold = true; // manifoldness check after collapse
 
private:
    // `ok` is set to false if the operation ran out of preallocated tet slots; in
    // that case no connectivity is mutated and the caller must abort the operation.
    std::map<size_t, VertexConnectivity> operation_update_connectivity_impl(
        std::vector<size_t>& affected_tid,
        const std::vector<std::array<size_t, 4>>& new_tet_conn,
        bool& ok);
    void operation_failure_rollback_imp(
        std::map<size_t, VertexConnectivity>& rollback_vert_conn,
        const std::vector<size_t>& affected,
        const std::vector<size_t>& new_tet_id,
        const std::vector<TetrahedronConnectivity>& old_tets);
    std::map<size_t, VertexConnectivity> operation_update_connectivity_impl(
        const std::vector<size_t>& remove_id,
        const std::vector<std::array<size_t, 4>>& new_tet_conn,
        std::vector<size_t>& allocate_id,
        bool& ok);
    static std::vector<TetrahedronConnectivity> record_old_tet_connectivity(
        const SlotPool<TetrahedronConnectivity>& conn,
        const std::vector<size_t>& tets)
    {
        std::vector<TetrahedronConnectivity> tet_conn;
        for (size_t i : tets) {
            tet_conn.push_back(conn[i]);
        }
        return tet_conn;
    }
 
public:
    void start_protect_attributes()
    {
        if (p_vertex_attrs) {
            p_vertex_attrs->begin_protect();
        }
        if (p_edge_attrs) {
            p_edge_attrs->begin_protect();
        }
        if (p_face_attrs) {
            p_face_attrs->begin_protect();
        }
        if (p_tet_attrs) {
            p_tet_attrs->begin_protect();
        }
    }
 
    void release_protect_attributes()
    {
        if (p_vertex_attrs) {
            p_vertex_attrs->end_protect();
        }
        if (p_edge_attrs) {
            p_edge_attrs->end_protect();
        }
        if (p_face_attrs) {
            p_face_attrs->end_protect();
        }
        if (p_tet_attrs) {
            p_tet_attrs->end_protect();
        }
    }
 
    void rollback_protected_attributes()
    {
        if (p_vertex_attrs) {
            p_vertex_attrs->rollback();
        }
        if (p_edge_attrs) {
            p_edge_attrs->rollback();
        }
        if (p_face_attrs) {
            p_face_attrs->rollback();
        }
        if (p_tet_attrs) {
            p_tet_attrs->rollback();
        }
    }
 
public:
    using VertexMutex = wmtk::threading::VertexMutex;
 
private:
    std::vector<VertexMutex> m_vertex_mutex;
 
    bool try_set_vertex_mutex(const Tuple& v, int threadid)
    {
        bool got = m_vertex_mutex[v.vid(*this)].trylock();
        if (got) m_vertex_mutex[v.vid(*this)].set_owner(threadid);
        return got;
    }
    bool try_set_vertex_mutex(size_t vid, int threadid)
    {
        bool got = m_vertex_mutex[vid].trylock();
        if (got) m_vertex_mutex[vid].set_owner(threadid);
        return got;
    }
 
    void unlock_vertex_mutex(const Tuple& v) { m_vertex_mutex[v.vid(*this)].unlock(); }
    void unlock_vertex_mutex(size_t vid) { m_vertex_mutex[vid].unlock(); }
 
    struct RingLockScratch
    {
        std::vector<uint32_t> stamp;
        uint32_t epoch = 0;
        std::vector<size_t> frontier;
        std::vector<size_t> next;
        std::vector<size_t> one_ring;
    };
    wmtk::threading::enumerable_thread_specific<RingLockScratch> m_ring_lock_scratch;
 
    std::shared_ptr<history::TrackedSurfaceRecorder> m_surface_history;
    wmtk::threading::enumerable_thread_specific<std::vector<std::array<size_t, 3>>>
        m_surface_history_before;
    void surface_history_before(const std::vector<size_t>& vs, bool all);
    std::array<size_t, 2> surface_history_commit(
        history::EditKind kind,
        std::vector<size_t> v,
        const std::vector<size_t>& vs_after,
        bool all);
 
    struct VolumeRecording;
    std::shared_ptr<VolumeRecording> m_volume_history;
    history::VolumeCapture volume_capture(const std::vector<size_t>& vs) const;
    void volume_history_before(const std::vector<size_t>& vs);
    void volume_history_commit(
        history::VolumeOp op,
        std::array<size_t, 2> surface_edits,
        const std::vector<size_t>& op_vertices,
        const std::vector<size_t>& vs_after);
    std::vector<size_t> vertices_of_tets(const std::vector<size_t>& tids) const;
 
    bool lock_vertex_ball(const size_t* seeds, size_t n_seeds, int threadid, int n, size_t mark);
 
    void check_smoothing_stayed_in_star(size_t vid);
 
protected:
    void resize_vertex_mutex(size_t v)
    {
        if (m_vertex_mutex.size() < v) m_vertex_mutex.resize(v);
    }
 
public:
    wmtk::threading::enumerable_thread_specific<std::vector<size_t>> mutex_release_stack;
    wmtk::threading::enumerable_thread_specific<std::vector<size_t>> get_one_ring_cache;
 
    // void init(size_t n_vertices, const std::vector<std::array<size_t, 4>>& tets);
    int release_vertex_mutex_in_stack();
    int release_vertex_mutex_to(size_t mark);
 
    bool try_set_vertex_mutex_n_ring(const Tuple& v, int threadid, int n);
    bool try_set_vertex_mutex_n_ring(size_t vid, int threadid, int n);
    bool try_set_edge_mutex_n_ring(const Tuple& e, int threadid, int n);
    bool try_set_face_mutex_n_ring(size_t v1, size_t v2, size_t v3, int threadid, int n);
 
    bool try_set_vertex_mutex_two_ring(const Tuple& v, int threadid);
    bool try_set_vertex_mutex_two_ring_vid(const Tuple& v, int threadid);
    bool try_set_vertex_mutex_two_ring_vid(size_t v, int threadid);
    bool try_set_edge_mutex_two_ring(const Tuple& e, int threadid = 0);
    bool try_set_face_mutex_two_ring(const Tuple& f, int threadid = 0);
    bool try_set_face_mutex_two_ring(
        const Tuple& v1,
        const Tuple& v2,
        const Tuple& v3,
        int threadid = 0);
    bool try_set_face_mutex_two_ring(size_t v1, size_t v2, size_t v3, int threadid = 0);
    bool try_set_vertex_mutex_one_ring(const Tuple& v, int threadid = 0);
public:
    void for_each_edge(const std::function<void(const TetMesh::Tuple&)>&);
    void for_each_vertex(const std::function<void(const TetMesh::Tuple&)>&);
    void for_each_tetra(const std::function<void(const TetMesh::Tuple&)>&);
    int NUM_THREADS = 0;
 
public:
    // substructure functionality
 
    virtual bool vertex_is_on_surface(const size_t vid) const { return false; }
 
    virtual bool face_is_on_surface(const size_t fid) const { return false; }
 
    simplex::SimplexCollection get_surface_faces_for_vertex(const size_t vid) const;
 
    simplex::SimplexCollection get_surface_faces_for_edge(const std::array<size_t, 2>& vids) const;
    size_t get_num_surface_faces_for_edge(const std::array<size_t, 2>& vids) const;
 
 
    size_t compute_vertex_order(const size_t vid) const;
 
    virtual size_t get_order_of_vertex(const size_t vid) const { return compute_vertex_order(vid); }
 
    size_t get_order_of_edge(const std::array<size_t, 2>& vids) const;
 
    virtual size_t substructure_order_of_edge(const std::array<size_t, 2>& vids) const
    {
        return get_order_of_edge(vids);
    }
    virtual size_t substructure_order_of_vertex(const size_t vid) const
    {
        return get_order_of_vertex(vid);
    }
    virtual void substructure_feature_neighbors(size_t, std::vector<size_t>&) const {}
 
    bool substructure_link_condition(const Tuple& e_tuple) const;
};
 
 
} // namespace wmtk