generic_RVD_cell.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2000-2022 Inria
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *
 *  * Redistributions of source code must retain the above copyright notice,
 *  this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright notice,
 *  this list of conditions and the following disclaimer in the documentation
 *  and/or other materials provided with the distribution.
 *  * Neither the name of the ALICE Project-Team nor the names of its
 *  contributors may be used to endorse or promote products derived from this
 *  software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 *  Contact: Bruno Levy
 *
 *     https://www.inria.fr/fr/bruno-levy
 *
 *     Inria,
 *     Domaine de Voluceau,
 *     78150 Le Chesnay - Rocquencourt
 *     FRANCE
 *
 */
 
#ifndef GEOGRAM_VORONOI_GENERIC_RVD_CELL
#define GEOGRAM_VORONOI_GENERIC_RVD_CELL
 
#include <geogram/basic/common.h>
#include <geogram/voronoi/generic_RVD_vertex.h>
#include <geogram/basic/argused.h>
#include <geogram/basic/attributes.h>
#include <iosfwd>
#include <stack>
 
namespace GEOGen {
 
    using GEO::Mesh;
 
    class GEOGRAM_API ConvexCell {
 
        typedef ConvexCell thisclass;
 
    public:
 
    static const index_t NO_TRIANGLE = index_t(-1);
    static const index_t NO_VERTEX = index_t(-1);
    static const index_t END_OF_LIST = index_t(-1);
 
    enum TriangleStatus {
        TRI_IS_FREE = 0,
        TRI_IS_CONFLICT = 1,
        TRI_IS_USED = 2
    };
 
    struct Triangle {
 
        Triangle(
            index_t v0, index_t v1, index_t v2,
            index_t f0, index_t f1, index_t f2
        ) :
            next_(END_OF_LIST),
            status_(TRI_IS_FREE),
            id_(-1) {
            v[0] = v0;
            v[1] = v1;
            v[2] = v2;
            t[0] = f0;
            t[1] = f1;
            t[2] = f2;
        }
 
        Triangle() :
            next_(END_OF_LIST),
            status_(TRI_IS_FREE),
            id_(-1) {
            v[0] = NO_VERTEX;
            v[1] = NO_VERTEX;
            v[2] = NO_VERTEX;
            t[0] = NO_TRIANGLE;
            t[1] = NO_TRIANGLE;
            t[2] = NO_TRIANGLE;
        }
 
        GEOGen::Vertex dual_;
        index_t v[3];   // The 3 vertices of this triangle
        index_t t[3];   // The 3 triangles adjacent to this triangle
        index_t next_;  // Linked list management.
        TriangleStatus status_;
        // TRI_IS_FREE,TRI_IS_USED or TRI_IS_CONFLICT.
        signed_index_t id_;
    };
 
    class Vertex {
    public:
        Vertex() :
            t(-1),
            id_(-1) {
        }
 
        signed_index_t t;  // One triangle incident to this vertex
        signed_index_t id_;
    };
 
    ConvexCell(coord_index_t dim) :
    first_free_(END_OF_LIST),
    v_to_t_dirty_(false),
    intersections_(dim),
    symbolic_is_surface_(false),
    cell_id_(-1) {
    }
 
    void copy(const ConvexCell& rhs);
 
    coord_index_t dimension() const {
        return intersections_.dimension();
    }
 
    void clear() {
        first_free_ = END_OF_LIST;
        triangles_.resize(0);
        vertices_.resize(0);
        v_to_t_dirty_ = false;
        intersections_.clear();
    }
 
    void set_symbolic_is_surface(bool x) {
        symbolic_is_surface_ = x;
    }
 
    void initialize_from_mesh_tetrahedron(
        const Mesh* mesh, index_t t, bool symbolic,
        const GEO::Attribute<double>& vertex_weight
    );
 
 
    void initialize_from_surface_mesh(
        Mesh* mesh, bool symbolic
    );
 
 
    void convert_to_mesh(Mesh* mesh, bool copy_symbolic_info = false);
 
    template <index_t DIM>
    signed_index_t clip_by_plane(
        const Mesh* mesh, const Delaunay* delaunay,
        index_t i, index_t j,
        bool exact, bool symbolic
    ) {
        index_t new_v = create_vertex();
        set_vertex_id(index_t(new_v), signed_index_t(j)+1);
        index_t conflict_begin, conflict_end;
 
        // Phase I: Determine the conflict zone and chain the triangles
        // Note: they are not immediately deleted, since we need the
        // geometric information in the triangles to compute the new
        // intersections.
        get_conflict_list<DIM>(
            mesh, delaunay, i, j, exact, conflict_begin, conflict_end
        );
 
        // Special case: the clipping plane did not clip anything.
        if(conflict_begin == END_OF_LIST) {
            return signed_index_t(new_v);
        }
 
        // Phase II: Find a triangle on the border of the conflict zone
        // (by traversing the conflict list).
        index_t first_conflict_t;
        index_t first_conflict_e;
        bool found_h = find_triangle_on_border(
            conflict_begin, conflict_end,
            first_conflict_t, first_conflict_e
        );
 
        // The clipping plane removed everything !
        // (note: cannot be empty conflict list, since this
        //  case was detected by previous test at the end of
        //  phase I)
        if(!found_h) {
            clear();
            return -1;
        }
 
        // Phase III: Triangulate hole.
        triangulate_hole<DIM>(
            delaunay, i, j, symbolic,
            first_conflict_t, first_conflict_e,
            new_v
        );
 
        // Phase IV: Merge the conflict zone into the free list.
        merge_into_free_list(conflict_begin, conflict_end);
        return signed_index_t(new_v);
    }
 
    index_t max_t() const {
        return triangles_.size();
    }
 
    index_t nb_t() const {
        index_t result = 0;
        for(index_t t = 0; t < max_t(); t++) {
            if(triangle_is_used(t)) {
                result++;
            }
        }
        return result;
    }
 
    index_t max_v() const {
        return vertices_.size();
    }
 
    bool triangle_is_free(index_t t) const {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        return triangles_[t].status_ == TRI_IS_FREE;
    }
 
    bool triangle_is_valid(index_t t) const {
        return !triangle_is_free(t);
    }
 
    bool triangle_is_used(index_t t) const {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        return triangles_[t].status_ == TRI_IS_USED;
    }
 
    bool triangle_is_conflict(index_t t) const {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        return triangles_[t].status_ == TRI_IS_CONFLICT;
    }
 
    index_t triangle_vertex(index_t t, index_t iv) const {
        geo_debug_assert(iv < 3);
        geo_debug_assert(triangle_is_valid(t));
        return triangles_[t].v[iv];
    }
 
    index_t triangle_adjacent(index_t t, index_t e) const {
        geo_debug_assert(e < 3);
        geo_debug_assert(triangle_is_valid(t));
        return triangles_[t].t[e];
    }
 
    void set_triangle_vertex(index_t t, index_t iv, index_t v) {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        geo_debug_assert(iv < 3);
        geo_debug_assert(v < max_v());
        triangles_[t].v[iv] = v;
    }
 
    void set_triangle_adjacent(
        index_t t, index_t e, index_t t2
    ) {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        geo_debug_assert(e < 3);
        geo_debug_assert(t2 < max_t());
        triangles_[t].t[e] = t2;
    }
 
    index_t find_triangle_vertex(index_t t, index_t v) const {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        geo_debug_assert(v < max_v());
 
        // The following expression is 10% faster than using
        // if() statements (multiply by boolean result of test).
        // Thank to Laurent Alonso for this idea.
        index_t result = index_t(
            (triangles_[t].v[1] == v) | ((triangles_[t].v[2] == v) * 2)
        );
        geo_debug_assert(triangles_[t].v[result] == v);
        return result;
    }
 
    index_t triangle_adjacent_index(
        index_t t1, index_t t2
    ) const {
        geo_debug_assert(t1 != NO_TRIANGLE);
        geo_debug_assert(t1 < max_t());
        geo_debug_assert(t2 != NO_TRIANGLE);
        geo_debug_assert(t2 < max_t());
 
        // The following expression is 10% faster than using
        // if() statements (multiply by boolean result of test).
        // Thank to Laurent Alonso for this idea.
        index_t result = index_t(
            (triangles_[t1].t[1] == t2) | ((triangles_[t1].t[2] == t2) * 2)
        );
 
        geo_debug_assert(triangles_[t1].t[result] == t2);
 
        return result;
    }
 
    signed_index_t vertex_triangle(index_t v) const {
        if(v_to_t_dirty_) {
            const_cast<ConvexCell*>(this)->init_v_to_t();
        }
        geo_debug_assert(v != NO_VERTEX);
        geo_debug_assert(v < max_v());
        return vertices_[v].t;
    }
 
    void set_vertex_triangle(index_t v, index_t t) {
        geo_debug_assert(v != NO_VERTEX);
        geo_debug_assert(index_t(v) < max_v());
        vertices_[v].t = signed_index_t(t);
    }
 
    const GEOGen::Vertex& triangle_dual(index_t t) const {
        geo_debug_assert(triangle_is_valid(t));
        return triangles_[t].dual_;
    }
 
    GEOGen::Vertex& triangle_dual(index_t t) {
        geo_debug_assert(triangle_is_valid(t));
        return triangles_[t].dual_;
    }
 
    signed_index_t cell_id() const {
        return cell_id_;
    }
 
    void set_cell_id(signed_index_t i) {
        cell_id_ = i;
    }
 
    signed_index_t triangle_id(index_t t) const {
        geo_debug_assert(triangle_is_valid(t));
        return triangles_[t].id_;
    }
 
    void set_triangle_id(index_t t, signed_index_t id) {
        geo_debug_assert(triangle_is_valid(t));
        triangles_[t].id_ = id;
    }
 
    signed_index_t vertex_id(index_t v) const {
        geo_debug_assert(v != NO_VERTEX);
        geo_debug_assert(v < max_v());
        return vertices_[v].id_;
    }
 
    void set_vertex_id(index_t v, signed_index_t id) {
        geo_debug_assert(v != NO_VERTEX);
        geo_debug_assert(v < max_v());
        vertices_[v].id_ = id;
    }
 
    class Corner {
    public:
        Corner() :
            t(NO_TRIANGLE),
            v(3) {
        }
 
        Corner(index_t t_in, index_t v_in) :
            t(t_in),
            v(v_in) {
        }
 
        bool operator== (const Corner& rhs) const {
            return t == rhs.t && v == rhs.v;
        }
 
        bool operator!= (const Corner& rhs) const {
            return t != rhs.t || v != rhs.v;
        }
 
        index_t t;
        index_t v;
    };
 
    void move_to_next_around_vertex(Corner& c) const {
        index_t t2 = triangle_adjacent(c.t, plus1mod3(c.v));
        index_t v = triangle_vertex(c.t, c.v);
        c.v = find_triangle_vertex(t2, v);
        c.t = t2;
    }
 
    void init_v_to_t() {
        v_to_t_dirty_ = false;
        for(index_t v = 0; v < max_v(); v++) {
            set_vertex_triangle(v, NO_TRIANGLE);
        }
        for(index_t t = 0; t < max_t(); t++) {
            if(triangle_is_used(t)) {
                for(index_t iv = 0; iv < 3; iv++) {
                    set_vertex_triangle(triangle_vertex(t, iv), t);
                }
            }
        }
    }
 
    index_t create_triangle() {
        if(first_free_ == END_OF_LIST) {
            grow();
        }
        index_t result = first_free_;
        first_free_ = next_triangle(first_free_);
        mark_as_used(result);
        set_triangle_id(result, -1);
        return result;
    }
 
    index_t create_triangle(index_t v0, index_t v1, index_t v2) {
        index_t t = create_triangle();
        triangles_[t].v[0] = v0;
        triangles_[t].v[1] = v1;
        triangles_[t].v[2] = v2;
        return t;
    }
 
    index_t create_triangle(
        index_t v0, index_t v1, index_t v2,
        index_t t0, index_t t1, index_t t2
    ) {
        index_t t = create_triangle();
        triangles_[t].v[0] = v0;
        triangles_[t].v[1] = v1;
        triangles_[t].v[2] = v2;
        triangles_[t].t[0] = t0;
        triangles_[t].t[1] = t1;
        triangles_[t].t[2] = t2;
        return t;
    }
 
    void set_triangle(
        index_t t,
        index_t v0, index_t v1, index_t v2,
        index_t t0, index_t t1, index_t t2
    ) {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        triangles_[t].v[0] = v0;
        triangles_[t].v[1] = v1;
        triangles_[t].v[2] = v2;
        triangles_[t].t[0] = t0;
        triangles_[t].t[1] = t1;
        triangles_[t].t[2] = t2;
    }
 
    index_t create_triangle(
        const double* p,
        double w,
        index_t v0, index_t v1, index_t v2,
        index_t t0, index_t t1, index_t t2
    ) {
        index_t t = create_triangle(v0, v1, v2, t0, t1, t2);
        triangle_dual(t).set_point(p);
        triangle_dual(t).set_weight(w);
        return t;
    }
 
    index_t create_triangle_copy(
        const double* p,
        index_t v0, index_t v1, index_t v2,
        index_t t0, index_t t1, index_t t2
    ) {
        index_t t = create_triangle(v0, v1, v2, t0, t1, t2);
        double* np = intersections_.new_item();
        for(coord_index_t c = 0; c < dimension(); ++c) {
            np[c] = p[c];
        }
        triangle_dual(t).set_point(np);
        return t;
    }
 
    index_t create_vertex() {
        v_to_t_dirty_ = true;
        vertices_.push_back(Vertex());
        return vertices_.size() - 1;
    }
 
    template <index_t DIM>
    index_t triangulate_hole(
        const Delaunay* delaunay,
        index_t i, index_t j, bool symbolic,
        index_t t1, index_t t1ebord,
        index_t v_in
    ) {
        index_t t = t1;
        index_t e = t1ebord;
        index_t t_adj = triangle_adjacent(t,e);
        geo_debug_assert(t_adj != index_t(-1));
 
        geo_debug_assert(triangle_is_conflict(t));
        geo_debug_assert(!triangle_is_conflict(t_adj));
 
        index_t new_t_first = index_t(-1);
        index_t new_t_prev  = index_t(-1);
 
        do {
 
            index_t v1 = triangle_vertex(t, plus1mod3(e));
            index_t v2 = triangle_vertex(t, minus1mod3(e));
 
            // Create new triangle
            index_t new_t = create_triangle(v_in, v1, v2);
 
            triangle_dual(new_t).intersect_geom<DIM>(
                intersections_,
                triangle_dual(t),
                triangle_dual(triangle_adjacent(t, e)),
                delaunay->vertex_ptr(i), delaunay->vertex_ptr(j)
            );
 
            if(symbolic) {
                triangle_dual(new_t).sym().intersect_symbolic(
                    triangle_dual(t).sym(),
                    triangle_dual(triangle_adjacent(t, e)).sym(),
                    j
                );
            }
 
            //   Connect new triangle to triangle on the other
            // side of the conflict zone.
            set_triangle_adjacent(new_t, 0, t_adj);
            index_t adj_e = triangle_adjacent_index(t_adj, t);
            set_triangle_adjacent(t_adj, adj_e, new_t);
 
 
            // Move to next triangle
            e = plus1mod3(e);
            t_adj = index_t(triangle_adjacent(t,e));
            while(triangle_is_conflict(t_adj)) {
                t = t_adj;
                e = minus1mod3(find_triangle_vertex(t,v2));
                t_adj = index_t(triangle_adjacent(t,e));
                geo_debug_assert(t_adj != index_t(-1));
            }
 
            if(new_t_prev == index_t(-1)) {
                new_t_first = new_t;
            } else {
                set_triangle_adjacent(new_t_prev, 1, new_t);
                set_triangle_adjacent(new_t, 2, new_t_prev);
            }
 
            new_t_prev = new_t;
 
        } while((t != t1) || (e != t1ebord));
 
        // Connect last triangle to first triangle
        set_triangle_adjacent(new_t_prev, 1, new_t_first);
        set_triangle_adjacent(new_t_first, 2, new_t_prev);
 
        return new_t_prev;
    }
 
    template <index_t DIM>
    void get_conflict_list(
        const Mesh* mesh, const Delaunay* delaunay,
        index_t i, index_t j, bool exact,
        index_t& conflict_begin, index_t& conflict_end
    ) {
        conflict_begin = END_OF_LIST;
        conflict_end = END_OF_LIST;
        if(exact) {
            // In exact mode, we classify each vertex
            // using the exact predicate. Note that
            // "climbing/walking" from a random vertex
            // would be more efficient, but it would require a
            // "comparison" exact predicate (with 8 different
            // versions according to the configuration of the
            // two vertices to be compared)
            for(index_t t = 0; t < max_t(); t++) {
                if(triangle_is_used(t)) {
                    Sign s = side<DIM>(
                        mesh, delaunay,
                        triangle_dual(t),
                        i, j,
                        exact
                    );
                    if(s == GEO::NEGATIVE) {
                        append_triangle_to_conflict_list(
                            t, conflict_begin, conflict_end
                        );
                    }
                }
            }
        } else {
            // In non-exact mode, we first detect the
            // vertex that is furthest away along the
            // normal vector of the clipping bisector,
            // and then we propagate using a flood-fill
            // algorithm. This strategy ensures that
            // the conflict zone remains connected, even
            // in the presence of numerical errors. Clearly
            // it does not ensure validity, but in practice
            // it improves resistance to degeneracies.
            index_t furthest_t =
                find_furthest_point_linear_scan<DIM>(
                    delaunay, i, j
                );
            propagate_conflict_list<DIM>(
                mesh, delaunay, furthest_t,
                i, j, exact,
                conflict_begin, conflict_end
            );
        }
    }
 
    template <index_t DIM>
    index_t find_furthest_point_linear_scan(
        const Delaunay* delaunay, index_t i, index_t j
    ) const {
        index_t result = NO_TRIANGLE;
        double furthest_dist = 0.0;
        for(index_t t = 0; t < max_t(); ++t) {
            if(triangle_is_used(t)) {
                double d = signed_bisector_distance<DIM>(
                    delaunay, i, j, triangle_dual(t).point()
                );
                if(d < furthest_dist) {
                    result = t;
                    furthest_dist = d;
                }
            }
        }
        return (furthest_dist < 0) ? result : NO_TRIANGLE;
    }
 
    template <index_t DIM>
    static double signed_bisector_distance(
        const Delaunay* delaunay, index_t i, index_t j, const double* q
    ) {
        const double* pi = delaunay->vertex_ptr(i);
        const double* pj = delaunay->vertex_ptr(j);
        double result = 0;
        for(coord_index_t c = 0; c < DIM; ++c) {
            result += GEO::geo_sqr(q[c] - pj[c]);
            result -= GEO::geo_sqr(q[c] - pi[c]);
        }
        return result;
    }
 
    template <index_t DIM>
    void propagate_conflict_list(
        const Mesh* mesh, const Delaunay* delaunay,
        index_t first_t,
        index_t i, index_t j, bool exact,
        index_t& conflict_begin, index_t& conflict_end
    ) {
        conflict_begin = END_OF_LIST;
        conflict_end = END_OF_LIST;
 
        // Special case, clipping plane does not clip anything
        if(first_t == NO_TRIANGLE) {
            return;
        }
 
        std::stack<index_t> S;
        S.push(first_t);
        append_triangle_to_conflict_list(
            first_t, conflict_begin, conflict_end
        );
        while(!S.empty()) {
            index_t t = S.top();
            S.pop();
            for(unsigned int e = 0; e < 3; e++) {
                index_t neigh = index_t(triangle_adjacent(t, e));
                if(!triangle_is_conflict(neigh)) {
                    if(
                        side<DIM>(
                            mesh, delaunay, triangle_dual(neigh),
                            i, j, exact
                        ) == GEO::NEGATIVE
                    ) {
                        S.push(neigh);
                        append_triangle_to_conflict_list(
                            neigh, conflict_begin, conflict_end
                        );
                    }
                }
            }
        }
    }
 
    template <index_t DIM>
    Sign side(
        const Mesh* mesh, const Delaunay* delaunay,
        const GEOGen::Vertex& v,
        index_t i, index_t j, bool exact
    ) const {
        Sign result = GEO::ZERO;
        if(exact) {
            result = side_exact(
                mesh, delaunay, v,
                delaunay->vertex_ptr(i),
                delaunay->vertex_ptr(j),
                DIM,
                symbolic_is_surface_
            );
        } else {
            result = v.side_fast<DIM>(
                delaunay->vertex_ptr(i),
                delaunay->vertex_ptr(j)
            );
        }
        return result;
    }
 
    Sign side_exact(
        const Mesh* mesh, const Delaunay* delaunay,
        const GEOGen::Vertex& v,
        const double* pi, const double* pj,
        coord_index_t dim,
        bool symbolic_is_surface = false
    ) const;
 
    bool find_triangle_on_border(
        index_t conflict_begin, index_t conflict_end,
        index_t& t, index_t& e
    ) const {
        GEO::geo_argused(conflict_end);
        t = conflict_begin;
        do {
            for(e = 0; e < 3; ++e) {
                index_t nt = triangle_adjacent(t, e);
                if(triangle_is_used(nt)) {
                    return true;
                }
            }
            t = next_triangle(t);
        } while(t != END_OF_LIST);
        return false;
    }
 
    index_t next_triangle(index_t t) const {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        return triangles_[t].next_;
    }
 
    void set_next_triangle(index_t t, index_t t2) {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        triangles_[t].next_ = t2;
    }
 
    void mark_as_free(index_t t) {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        triangles_[t].status_ = TRI_IS_FREE;
    }
 
    void mark_as_conflict(index_t t) {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        triangles_[t].status_ = TRI_IS_CONFLICT;
    }
 
    void mark_as_used(index_t t) {
        geo_debug_assert(t != NO_TRIANGLE);
        geo_debug_assert(t < max_t());
        triangles_[t].status_ = TRI_IS_USED;
    }
 
    void append_triangle_to_conflict_list(
        index_t t, index_t& conflict_begin, index_t& conflict_end
    ) {
        geo_debug_assert(triangle_is_used(t));
        set_next_triangle(t, conflict_begin);
        mark_as_conflict(t);
        conflict_begin = t;
        if(conflict_end == END_OF_LIST) {
            conflict_end = t;
        }
    }
 
    void merge_into_free_list(index_t list_begin, index_t list_end) {
        if(list_begin != END_OF_LIST) {
            geo_debug_assert(list_end != END_OF_LIST);
 
            index_t cur = list_begin;
            while(cur != list_end) {
                mark_as_free(cur);
                cur = next_triangle(cur);
            }
            mark_as_free(list_end);
            set_next_triangle(list_end, first_free_);
            first_free_ = list_begin;
        }
    }
 
    void grow() {
        geo_debug_assert(first_free_ == END_OF_LIST);
        first_free_ = triangles_.size();
        triangles_.push_back(Triangle());
    }
 
    std::ostream& show_stats(std::ostream& os) const;
 
    static index_t global_facet_id(
        const Mesh* mesh, index_t t, index_t lf
    ) {
        index_t t2 = mesh->cells.tet_adjacent(t, lf);
        if(t2 != GEO::NO_CELL && t2 > t) {
            index_t lf2 = mesh->cells.find_tet_adjacent(
                t2, t
            );
            geo_debug_assert(lf2 != GEO::NO_FACET);
            return index_t(4 * t2 + lf2);
        }
        return 4 * t + lf;
    }
 
    private:
    GEO::vector<Triangle> triangles_;
    GEO::vector<Vertex> vertices_;
    index_t first_free_;
    bool v_to_t_dirty_;
    PointAllocator intersections_;
    bool symbolic_is_surface_;
    signed_index_t cell_id_;
 
    static index_t plus1mod3_[3];
    static index_t minus1mod3_[3];
 
    static index_t plus1mod3(index_t i) {
        geo_debug_assert(i < 3);
        return plus1mod3_[i];
    }
 
    static index_t minus1mod3(index_t i) {
        geo_debug_assert(i < 3);
        return minus1mod3_[i];
    }
    };
}
 
#endif