TetMesh.h

#pragma once
 
#include <wmtk/utils/VectorUtils.h>
#include <wmtk/AttributeCollection.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 <array>
#include <atomic>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <limits>
#include <map>
#include <optional>
#include <vector>
 
namespace wmtk {
class TetMesh
{
private:
    static constexpr std::array<std::array<int, 2>, 6> m_local_edges = {
        {{{0, 1}}, {{1, 2}}, {{0, 2}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
 
    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_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;
 
        size_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(); }
    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 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::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);
 
    // 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;
        double c = std::ceil(m_preallocation_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; }
 
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)
    {
        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;
        }
    }
    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;
 
    bool lock_vertex_ball(const size_t* seeds, size_t n_seeds, int threadid, int n, size_t mark);
 
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;
 
    bool substructure_link_condition(const Tuple& e_tuple) const;
};
 
 
} // namespace wmtk