GCC Code Coverage Report

Directory: ./
File: lib/geogram/voronoi/generic_RVD_cell.h
Date: 2026-09-07 02:37:58
Exec Total Coverage
Lines: 353 353 100.0%
Functions: 74 95 77.9%
Branches: 155 450 34.4%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #ifndef GEOGRAM_VORONOI_GENERIC_RVD_CELL
41 #define GEOGRAM_VORONOI_GENERIC_RVD_CELL
42
43 #include <geogram/basic/common.h>
44 #include <geogram/voronoi/generic_RVD_vertex.h>
45 #include <geogram/basic/argused.h>
46 #include <geogram/basic/attributes.h>
47 #include <iosfwd>
48 #include <stack>
49
50 /**
51 * \file geogram/voronoi/generic_RVD_cell.h
52 * \brief Internal representation of polyhedra for GEO::GenericVoronoiDiagram.
53 * \note This file contains functions and classes used by the
54 * internal implementation of GEO::GenericVoronoiDiagram.
55 * They are not meant to be used directly by client code.
56 * Users who whant similar functionalities may use GEO::ConvexCell instead.
57 */
58
59 namespace GEOGen {
60
61 using GEO::Mesh;
62
63 /**
64 * \brief Computes the intersection between a set of halfspaces.
65 * \note This is an internal implementation class used by
66 * GEO::RestrictedVoronoiDiagram. It is not meant to be
67 * used directly by client code.
68 */
69 class GEOGRAM_API ConvexCell {
70
71 /** \brief This class type */
72 typedef ConvexCell thisclass;
73
74 public:
75
76 static constexpr index_t NO_TRIANGLE = index_t(-1);
77 static constexpr index_t NO_VERTEX = index_t(-1);
78 static constexpr index_t END_OF_LIST = index_t(-1);
79
80 /**
81 * \brief Represents the current state of a triangle.
82 */
83 enum TriangleStatus {
84 TRI_IS_FREE = 0,
85 TRI_IS_CONFLICT = 1,
86 TRI_IS_USED = 2
87 };
88
89 /**
90 * \brief Represents a vertex of this ConvexCell in dual form.
91 * \details Each vertex of a ConvexCell is represented
92 * combinatorially as a Triangle, in dual form
93 * (in primal parlance, each vertex is
94 * of degree 3, and can be represented by a triangle).
95 * A Triangle knows its tree vertices, its tree adjacent triangles
96 * and the geometric point it corresponds to.
97 */
98
99 struct Triangle {
100
101 /**
102 * \brief Creates a new triangle.
103 * \param[in] v0 index of the first vertex
104 * \param[in] v1 index of the second vertex
105 * \param[in] v2 index of the third vertex
106 * \param[in] f0 index of the triangle opposite to \p v0
107 * \param[in] f1 index of the triangle opposite to \p v1
108 * \param[in] f2 index of the triangle opposite to \p v2
109 */
110 Triangle(
111 index_t v0, index_t v1, index_t v2,
112 index_t f0, index_t f1, index_t f2
113 ) :
114 next_(END_OF_LIST),
115 status_(TRI_IS_FREE),
116 id_(-1) {
117 v[0] = v0;
118 v[1] = v1;
119 v[2] = v2;
120 t[0] = f0;
121 t[1] = f1;
122 t[2] = f2;
123 }
124
125 /**
126 * \brief Creates a new uninitialized Triangle.
127 */
128 10002420 Triangle() :
129 10002420 next_(END_OF_LIST),
130 10002420 status_(TRI_IS_FREE),
131 10002420 id_(-1) {
132 10002420 v[0] = NO_VERTEX;
133 10002420 v[1] = NO_VERTEX;
134 10002420 v[2] = NO_VERTEX;
135 10002420 t[0] = NO_TRIANGLE;
136 10002420 t[1] = NO_TRIANGLE;
137 10002420 t[2] = NO_TRIANGLE;
138 10002420 }
139
140 GEOGen::Vertex dual_;
141 index_t v[3]; // The 3 vertices of this triangle
142 index_t t[3]; // The 3 triangles adjacent to this triangle
143 index_t next_; // Linked list management.
144 TriangleStatus status_;
145 // TRI_IS_FREE,TRI_IS_USED or TRI_IS_CONFLICT.
146 signed_index_t id_;
147 };
148
149 /**
150 * \brief Represents a facet of this ConvexCell in dual form.
151 * \details Each facet of a ConvexCell is represented
152 * combinatorially as a Vertex (dual form).
153 * Each vertex v knows an incident triangle (v.t).
154 */
155 class Vertex {
156 public:
157 /**
158 * \brief Creates a new uninitialized Vertex.
159 */
160 13373954 Vertex() :
161 13373954 t(-1),
162 13373954 id_(-1) {
163 13373954 }
164
165 signed_index_t t; // One triangle incident to this vertex
166 signed_index_t id_;
167 };
168
169 /**
170 * \brief ConvexCell constructor.
171 * \param[in] dim dimension of the ConvexCell, e.g. 3 for 3d
172 */
173 650 ConvexCell(coord_index_t dim) :
174 650 first_free_(END_OF_LIST),
175 650 v_to_t_dirty_(false),
176 650 intersections_(dim),
177 650 symbolic_is_surface_(false),
178 650 cell_id_(-1) {
179 650 }
180
181 /**
182 * \brief Copies a ConvexCell.
183 * \details The allocated vertices are shared with \p rhs, thus
184 * \p rhs should not be deleted before this ConvexCell.
185 * \param[in] rhs a const reference to the ConvexCell to be
186 * copied.
187 */
188 void copy(const ConvexCell& rhs);
189
190 /**
191 * \brief Gets the dimension of this ConvexCell.
192 * \return the dimension of this ConvexCell, e.g. 3 for 3d
193 */
194 coord_index_t dimension() const {
195 return intersections_.dimension();
196 }
197
198 /**
199 * \brief Clears this ConvexCell.
200 */
201 792613 void clear() {
202 792613 first_free_ = END_OF_LIST;
203 792613 triangles_.resize(0);
204 792613 vertices_.resize(0);
205 792613 v_to_t_dirty_ = false;
206 792613 intersections_.clear();
207 792613 }
208
209 /**
210 * \brief Specifies that symbolic information is
211 * relative to surfacic mesh (rather than volumetric mesh).
212 * \details Affects the behavior of side_exact().
213 * \param[in] x true if symbolic information is relative
214 * to surfacic mesh (triangles), false if symbolic information
215 * is relative to volumetric mesh (tetrahedra).
216 */
217 18 void set_symbolic_is_surface(bool x) {
218 18 symbolic_is_surface_ = x;
219 18 }
220
221 /**
222 * \brief Assigns a mesh tetrahedron to this ConvexCell
223 * \details The tetrahedron from the initial mesh is converted into
224 * the internal geometric/symbolic representation.
225 * \param[in] mesh the mesh from which the tetrahedron is copied
226 * \param[in] t the index of the tetrahedron in \p mesh
227 * \param[in] symbolic if true, symbolic information is copied
228 * \param[in] vertex_weight if bound, an attribute that gives
229 * the weight of each vertex in \p mesh.
230 */
231 void initialize_from_mesh_tetrahedron(
232 const Mesh* mesh, index_t t, bool symbolic,
233 const GEO::Attribute<double>& vertex_weight
234 );
235
236
237 /**
238 * \brief Copies a Mesh into a ConvexCell
239 * \details The surface mesh in \p mesh represents the boundary
240 * of the ConvexCell.
241 * \param[in] mesh a pointer to the input Mesh
242 * \param[in] symbolic if true, symbolic information is copied
243 */
244 void initialize_from_surface_mesh(
245 Mesh* mesh, bool symbolic
246 );
247
248
249 /**
250 * \brief Copies a ConvexCell into a Mesh
251 * \details On exit, the output mesh is a surfacic
252 * mesh with the boundary of the convex cell.
253 * \param[out] mesh a pointer to the target mesh
254 * \param[in] copy_symbolic_info if true, symbolic
255 * information is copied. An attribute "id" is attached
256 * to the facets. The value of id[f] is either 1 + the index of
257 * the Voronoi vertex that generated with \p i the bisector that
258 * created the facet, or -1-g if the facet was an original facet
259 * of mesh \p mesh, where g is the index of the original
260 * facet in \p mesh.
261 */
262 void convert_to_mesh(Mesh* mesh, bool copy_symbolic_info = false);
263
264 /**
265 * \brief Clips this ConvexCell with a plane.
266 * \details The plane is specified as a bisector
267 * in a Delaunay triangulation.
268 * \param[in] mesh input mesh, used by exact predicates
269 * \param[in] delaunay the Delaunay triangulation
270 * \param[in] i index of the first extremity of the bisector
271 * \param[in] j index of the second extremity of the bisector
272 * \param[in] exact if true, uses exact predicates (implies symbolic)
273 * \param[in] symbolic if true, computes symbolic information
274 * \return the index of the newly created vertex that corresponds to
275 * the dual of the new face or -1 if the input cell was completely
276 * on the negative side (removed everything)
277 * \tparam DIM dimension (specified as a template parameter for
278 * efficiency reasons).
279 * \pre DIM == dimension()
280 */
281 template <index_t DIM>
282 20391200 signed_index_t clip_by_plane(
283 const Mesh* mesh, const Delaunay* delaunay,
284 index_t i, index_t j,
285 bool exact, bool symbolic
286 ) {
287
1/2
✓ Branch 1 taken 10208926 times.
✗ Branch 2 not taken.
 20391200 index_t new_v = create_vertex();
288
1/2
✓ Branch 1 taken 10208926 times.
✗ Branch 2 not taken.
 20391200 set_vertex_id(index_t(new_v), signed_index_t(j)+1);
289 index_t conflict_begin, conflict_end;
290
291 // Phase I: Determine the conflict zone and chain the triangles
292 // Note: they are not immediately deleted, since we need the
293 // geometric information in the triangles to compute the new
294 // intersections.
295
1/2
✓ Branch 1 taken 10208926 times.
✗ Branch 2 not taken.
 20391200 get_conflict_list<DIM>(
296 mesh, delaunay, i, j, exact, conflict_begin, conflict_end
297 );
298
299 // Special case: the clipping plane did not clip anything.
300
2/2
✓ Branch 0 taken 6240951 times.
✓ Branch 1 taken 3967975 times.
 20391200 if(conflict_begin == END_OF_LIST) {
301 12461901 return signed_index_t(new_v);
302 }
303
304 // Phase II: Find a triangle on the border of the conflict zone
305 // (by traversing the conflict list).
306 index_t first_conflict_t;
307 index_t first_conflict_e;
308
1/2
✓ Branch 1 taken 3967975 times.
✗ Branch 2 not taken.
 7929299 bool found_h = find_triangle_on_border(
309 conflict_begin, conflict_end,
310 first_conflict_t, first_conflict_e
311 );
312
313 // The clipping plane removed everything !
314 // (note: cannot be empty conflict list, since this
315 // case was detected by previous test at the end of
316 // phase I)
317
2/2
✓ Branch 0 taken 1365 times.
✓ Branch 1 taken 3966610 times.
 7929299 if(!found_h) {
318
1/2
✓ Branch 1 taken 1365 times.
✗ Branch 2 not taken.
 2730 clear();
319 2730 return -1;
320 }
321
322 // Phase III: Triangulate hole.
323
1/2
✓ Branch 1 taken 3966610 times.
✗ Branch 2 not taken.
 7926569 triangulate_hole<DIM>(
324 delaunay, i, j, symbolic,
325 first_conflict_t, first_conflict_e,
326 new_v
327 );
328
329 // Phase IV: Merge the conflict zone into the free list.
330
1/2
✓ Branch 1 taken 3966610 times.
✗ Branch 2 not taken.
 7926569 merge_into_free_list(conflict_begin, conflict_end);
331 7926569 return signed_index_t(new_v);
332 }
333
334 /**
335 * \brief Gets the maximum valid triangle index plus one.
336 * \details May be greater than nb_t() if this ConvexCell has
337 * some free (unused) triangles.
338 */
339 1383861564 index_t max_t() const {
340 1383861564 return triangles_.size();
341 }
342
343 /**
344 * \brief Gets the number of used triangles.
345 */
346 index_t nb_t() const {
347 index_t result = 0;
348 for(index_t t = 0; t < max_t(); t++) {
349 if(triangle_is_used(t)) {
350 result++;
351 }
352 }
353 return result;
354 }
355
356 /**
357 * \brief Gets the maximum valid vertex index plus one.
358 */
359 155501278 index_t max_v() const {
360 155501278 return vertices_.size();
361 }
362
363 /**
364 * \brief Tests whether a given triangle is free.
365 */
366 469283675 bool triangle_is_free(index_t t) const {
367
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 469283675 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 469283675 geo_debug_assert(t != NO_TRIANGLE);
368
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 469283675 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 469283675 geo_debug_assert(t < max_t());
369 469283675 return triangles_[t].status_ == TRI_IS_FREE;
370 }
371
372 /**
373 * \brief Tests whether a given triangle is valid.
374 * \details A "valid" triangle is a triangle from which
375 * we can query information (vertices, adjacent triangles,
376 * embedding), therefore it is a "used" or "conflict" triangle.
377 * \param[in] t index of the triangle
378 * \pre t < max_t()
379 */
380 469283675 bool triangle_is_valid(index_t t) const {
381 469283675 return !triangle_is_free(t);
382 }
383
384 /**
385 * \brief Tests whether a given triangle is used.
386 * \param[in] t index of the triangle
387 * \pre t < max_t()
388 */
389 299247678 bool triangle_is_used(index_t t) const {
390
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 299247678 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 299247678 geo_debug_assert(t != NO_TRIANGLE);
391
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 299247678 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 299247678 geo_debug_assert(t < max_t());
392 299247678 return triangles_[t].status_ == TRI_IS_USED;
393 }
394
395 /**
396 * \brief Tests whether a given triangle belongs to
397 * the conflict zone.
398 * \param[in] t index of the triangle
399 * \pre t < max_t()
400 */
401 63884216 bool triangle_is_conflict(index_t t) const {
402
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 63884216 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 63884216 geo_debug_assert(t != NO_TRIANGLE);
403
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 63884216 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 63884216 geo_debug_assert(t < max_t());
404 63884216 return triangles_[t].status_ == TRI_IS_CONFLICT;
405 }
406
407 /**
408 * \brief Gets the index of a triangle vertex.
409 * \param[in] t the triangle index
410 * \param[in] iv local vertex index (0,1 or 2) in \p t
411 * \return the index of triangle's \p iv th vertex
412 * \pre t < max_t()
413 */
414 64785394 index_t triangle_vertex(index_t t, index_t iv) const {
415
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 64785394 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 64785394 geo_debug_assert(iv < 3);
416
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 64785394 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 64785394 geo_debug_assert(triangle_is_valid(t));
417 64785394 return triangles_[t].v[iv];
418 }
419
420 /**
421 * \brief Gets the index of a triangle adjacent to another one.
422 * \param[in] t the triangle index
423 * \param[in] e local edge index (0,1 or 2) in \p t
424 * \return the index of the triangle adjacent to \p t on edge \p e
425 */
426 94969794 index_t triangle_adjacent(index_t t, index_t e) const {
427
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 94969794 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 94969794 geo_debug_assert(e < 3);
428
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 94969794 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 94969794 geo_debug_assert(triangle_is_valid(t));
429 94969794 return triangles_[t].t[e];
430 }
431
432 /**
433 * \brief Sets a vertex of a triangle.
434 * \param[in] t the triangle index
435 * \param[in] iv local vertex index (0,1 or 2) in \p t
436 * \param[in] v global vertex index
437 */
438 void set_triangle_vertex(index_t t, index_t iv, index_t v) {
439 geo_debug_assert(t != NO_TRIANGLE);
440 geo_debug_assert(t < max_t());
441 geo_debug_assert(iv < 3);
442 geo_debug_assert(v < max_v());
443 triangles_[t].v[iv] = v;
444 }
445
446 /**
447 * \brief Sets a triangle adjacency.
448 * \param[in] t the triangle index
449 * \param[in] e local edge index (0,1 or 2)
450 * \param[in] t2 global triangle index
451 */
452 58142476 void set_triangle_adjacent(
453 index_t t, index_t e, index_t t2
454 ) {
455
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 58142476 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 58142476 geo_debug_assert(t != NO_TRIANGLE);
456
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 58142476 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 58142476 geo_debug_assert(t < max_t());
457
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 58142476 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 58142476 geo_debug_assert(e < 3);
458
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 58142476 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 58142476 geo_debug_assert(t2 < max_t());
459 58142476 triangles_[t].t[e] = t2;
460 58142476 }
461
462 /**
463 * \brief Finds the local index of a triangle vertex
464 * \param[in] t the triangle
465 * \param[in] v global vertex index
466 * \return local vertex index (0,1 or 2) of \p v in \p t
467 * \pre \p t is incident to \p v
468 */
469 28270374 index_t find_triangle_vertex(index_t t, index_t v) const {
470
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 28270374 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 28270374 geo_debug_assert(t != NO_TRIANGLE);
471
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 28270374 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 28270374 geo_debug_assert(t < max_t());
472
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 28270374 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 28270374 geo_debug_assert(v < max_v());
473
474 // The following expression is 10% faster than using
475 // if() statements (multiply by boolean result of test).
476 // Thank to Laurent Alonso for this idea.
477 index_t result = index_t(
478
2/2
✓ Branch 2 taken 7454430 times.
✓ Branch 3 taken 20815944 times.
 28270374 (triangles_[t].v[1] == v) | ((triangles_[t].v[2] == v) * 2)
479 28270374 );
480
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 28270374 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 28270374 geo_debug_assert(triangles_[t].v[result] == v);
481 28270374 return result;
482 }
483
484 /**
485 * \brief Finds the edge along which two triangles are adjacent
486 * \param[in] t1 first triangle
487 * \param[in] t2 second triangle
488 * \return the local edge index (0,1 or 2) in \p t1 along which
489 * \p t2 is adjacent to \p t1
490 * \pre \p t1 and \p t2 are adjacent
491 */
492 14535619 index_t triangle_adjacent_index(
493 index_t t1, index_t t2
494 ) const {
495
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 14535619 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 14535619 geo_debug_assert(t1 != NO_TRIANGLE);
496
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 14535619 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 14535619 geo_debug_assert(t1 < max_t());
497
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 14535619 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 14535619 geo_debug_assert(t2 != NO_TRIANGLE);
498
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 14535619 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 14535619 geo_debug_assert(t2 < max_t());
499
500 // The following expression is 10% faster than using
501 // if() statements (multiply by boolean result of test).
502 // Thank to Laurent Alonso for this idea.
503 index_t result = index_t(
504
2/2
✓ Branch 2 taken 4452447 times.
✓ Branch 3 taken 10083172 times.
 14535619 (triangles_[t1].t[1] == t2) | ((triangles_[t1].t[2] == t2) * 2)
505 14535619 );
506
507
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 14535619 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 14535619 geo_debug_assert(triangles_[t1].t[result] == t2);
508
509 14535619 return result;
510 }
511
512 /**
513 * \brief Gets one of the triangles incident to a vertex.
514 * \details The information is cached and reconstructed
515 * whenever the v_to_t_dirty_ flag is positioned.
516 * \param[in] v index of the vertex
517 * \return the index of a triangle incident to \p v
518 */
519 26683816 signed_index_t vertex_triangle(index_t v) const {
520
2/2
✓ Branch 0 taken 789883 times.
✓ Branch 1 taken 25893933 times.
 26683816 if(v_to_t_dirty_) {
521 789883 const_cast<ConvexCell*>(this)->init_v_to_t();
522 }
523
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 26683816 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 26683816 geo_debug_assert(v != NO_VERTEX);
524
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 26683816 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 26683816 geo_debug_assert(v < max_v());
525 26683816 return vertices_[v].t;
526 }
527
528 /**
529 * \brief Stores in a vertex the index of one
530 * the triangles incident to it.
531 * \param[in] v index of the vertex
532 * \param[in] t index of one of the triangles incident to \p v
533 */
534 34551452 void set_vertex_triangle(index_t v, index_t t) {
535
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 34551452 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 34551452 geo_debug_assert(v != NO_VERTEX);
536
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 34551452 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 34551452 geo_debug_assert(index_t(v) < max_v());
537 34551452 vertices_[v].t = signed_index_t(t);
538 34551452 }
539
540 /**
541 * \brief Gets the dual vertex that corresponds to a triangle.
542 * \details Each triangle corresponds to a vertex of the ConvexCell
543 * (combinatorics are stored in dual form).
544 * \param[in] t index of the triangle
545 * \return a const reference to the GEOGen::Vertex that corresponds
546 * to the triangle, with both geometrical and combinatorial
547 * representations
548 */
549 93312249 const GEOGen::Vertex& triangle_dual(index_t t) const {
550
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 93312249 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 93312249 geo_debug_assert(triangle_is_valid(t));
551 93312249 return triangles_[t].dual_;
552 }
553
554 /**
555 * \brief Gets the dual vertex that corresponds to a triangle.
556 * \details Each triangle corresponds to a vertex of the ConvexCell
557 * (combinatorics are stored in dual form).
558 * \param[in] t index of the triangle
559 * \return a reference to the GEOGen::Vertex that corresponds
560 * to the triangle, with both geometrical and
561 * combinatorial representations
562 */
563 198502178 GEOGen::Vertex& triangle_dual(index_t t) {
564
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 198502178 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 198502178 geo_debug_assert(triangle_is_valid(t));
565 198502178 return triangles_[t].dual_;
566 }
567
568 /**
569 * \brief Gets the id of this ConvexCell.
570 * \details The id can be used to maintain the
571 * correspondence with a tetrahedron in the Mesh.
572 * \return the id of this ConvexCell. Can be a
573 * negative number.
574 */
575 136346 signed_index_t cell_id() const {
576 136346 return cell_id_;
577 }
578
579 /**
580 * \brief Sets the id of this ConvexCell.
581 * \details The id can be used to maintain the
582 * correspondence with a tetrahedron in the Mesh.
583 * \param[in] i the id of this ConvexCell. Can be
584 * a negative number.
585 */
586 791230 void set_cell_id(signed_index_t i) {
587 791230 cell_id_ = i;
588 791230 }
589
590 /**
591 * \brief Gets the id of a triangle.
592 * \details Each triangle of a ConvexCell has an id, that can
593 * be used to maintain the correspondence with a vertex in the Mesh.
594 * \return The id of triangle \p t.
595 */
596 signed_index_t triangle_id(index_t t) const {
597 geo_debug_assert(triangle_is_valid(t));
598 return triangles_[t].id_;
599 }
600
601 /**
602 * \brief Sets the id of a triangle.
603 * \details Each triangle of a ConvexCell has an id, that can
604 * be used to maintain the correspondence with a vertex in the Mesh.
605 */
606 17700683 void set_triangle_id(index_t t, signed_index_t id) {
607
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 17700683 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 17700683 geo_debug_assert(triangle_is_valid(t));
608 17700683 triangles_[t].id_ = id;
609 17700683 }
610
611 /**
612 * \brief Gets the id of a vertex.
613 * \details Each vertex of a ConvexCell has an id, that can
614 * be used to maintain the correspondence with a facet in the Mesh.
615 */
616 10211582 signed_index_t vertex_id(index_t v) const {
617
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10211582 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 10211582 geo_debug_assert(v != NO_VERTEX);
618
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10211582 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10211582 geo_debug_assert(v < max_v());
619 10211582 return vertices_[v].id_;
620 }
621
622 /**
623 * \brief Sets the id of a vertex.
624 * \details Each vertex of a ConvexCell has an id, that can
625 * be used to maintain the correspondence with a facet in the Mesh.
626 */
627 13373954 void set_vertex_id(index_t v, signed_index_t id) {
628
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 13373954 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 13373954 geo_debug_assert(v != NO_VERTEX);
629
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 13373954 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 13373954 geo_debug_assert(v < max_v());
630 13373954 vertices_[v].id_ = id;
631 13373954 }
632
633 /**
634 * \brief A Corner corresponds to a vertex seen from a triangle.
635 * \details Corner has helper functions that facilitate traversing
636 * the vertices of a facet in dual form.
637 */
638 class Corner {
639 public:
640 /**
641 * \brief Creates an uninitialized Corner.
642 */
643 Corner() :
644 t(NO_TRIANGLE),
645 v(3) {
646 }
647
648 /**
649 * \brief Creates a Corner from a triangle index and
650 * local vertex index.
651 * \param[in] t_in index of the triangle
652 * \param[in] v_in local index (0,1 or 2) in triangle \p t_in.
653 */
654 3399890 Corner(index_t t_in, index_t v_in) :
655 3399890 t(t_in),
656 3399890 v(v_in) {
657 3399890 }
658
659 /**
660 * \brief Compares two corners.
661 * \param[in] rhs the right hand side
662 * \retval true if this Corner and \p rhs correspond to the same
663 * triangle and the same vertex
664 * \retval false otherwise
665 */
666 bool operator== (const Corner& rhs) const {
667 return t == rhs.t && v == rhs.v;
668 }
669
670 /**
671 * \brief Compares two corners.
672 * \param[in] rhs the right hand side
673 * \retval true if this Corner and \p rhs correspond to different
674 * triangles or different vertices
675 * \retval false otherwise
676 */
677 14517992 bool operator!= (const Corner& rhs) const {
678
3/4
✓ Branch 0 taken 3449259 times.
✓ Branch 1 taken 11068733 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3449259 times.
 14517992 return t != rhs.t || v != rhs.v;
679 }
680
681 index_t t;
682 index_t v;
683 };
684
685 /**
686 * \brief Replaces a corner by the next corner obtained by turing around
687 * the vertex.
688 * \param[in,out] c the corner
689 */
690 14517992 void move_to_next_around_vertex(Corner& c) const {
691 14517992 index_t t2 = triangle_adjacent(c.t, plus1mod3(c.v));
692 14517992 index_t v = triangle_vertex(c.t, c.v);
693 14517992 c.v = find_triangle_vertex(t2, v);
694 14517992 c.t = t2;
695 14517992 }
696
697 /**
698 * \brief Updates the cache that stores for each vertex a triangle
699 * incident to it.
700 */
701 789883 void init_v_to_t() {
702 789883 v_to_t_dirty_ = false;
703
2/2
✓ Branch 1 taken 13355288 times.
✓ Branch 2 taken 789883 times.
 14145171 for(index_t v = 0; v < max_v(); v++) {
704 13355288 set_vertex_triangle(v, NO_TRIANGLE);
705 }
706
2/2
✓ Branch 1 taken 9987685 times.
✓ Branch 2 taken 789883 times.
 10777568 for(index_t t = 0; t < max_t(); t++) {
707
2/2
✓ Branch 1 taken 7065388 times.
✓ Branch 2 taken 2922297 times.
 9987685 if(triangle_is_used(t)) {
708
2/2
✓ Branch 0 taken 21196164 times.
✓ Branch 1 taken 7065388 times.
 28261552 for(index_t iv = 0; iv < 3; iv++) {
709 21196164 set_vertex_triangle(triangle_vertex(t, iv), t);
710 }
711 }
712 }
713 789883 }
714
715 /**
716 * \brief Creates a new uninitialized triangle.
717 * \details The created triangle is marked as used.
718 * \return The index of the new triangle.
719 */
720 17700683 index_t create_triangle() {
721
2/2
✓ Branch 0 taken 10002420 times.
✓ Branch 1 taken 7698263 times.
 17700683 if(first_free_ == END_OF_LIST) {
722 10002420 grow();
723 }
724 17700683 index_t result = first_free_;
725 17700683 first_free_ = next_triangle(first_free_);
726 17700683 mark_as_used(result);
727 17700683 set_triangle_id(result, -1);
728 17700683 return result;
729 }
730
731 /**
732 * \brief Creates a new triangle with specified vertices.
733 * \details The created triangle is marked as used. Adjacent
734 * triangles are left uninitialized.
735 * \param[in] v0 index of first vertex
736 * \param[in] v1 index of second vertex
737 * \param[in] v2 index of third vertex
738 * \return the index of the new triangle
739 */
740 14535619 index_t create_triangle(index_t v0, index_t v1, index_t v2) {
741 14535619 index_t t = create_triangle();
742 14535619 triangles_[t].v[0] = v0;
743 14535619 triangles_[t].v[1] = v1;
744 14535619 triangles_[t].v[2] = v2;
745 14535619 return t;
746 }
747
748 /**
749 * \brief Creates a new triangles with specified vertices and
750 * adjacent triangles.
751 * \details The created triangle is marked as used.
752 * \param[in] v0 index of first vertex
753 * \param[in] v1 index of second vertex
754 * \param[in] v2 index of third vertex
755 * \param[in] t0 index of adjacent triangle opposite to \p v0
756 * \param[in] t1 index of adjacent triangle opposite to \p v1
757 * \param[in] t2 index of adjacent triangle opposite to \p v2
758 * \return the index of the new triangle
759 */
760 3165064 index_t create_triangle(
761 index_t v0, index_t v1, index_t v2,
762 index_t t0, index_t t1, index_t t2
763 ) {
764 3165064 index_t t = create_triangle();
765 3165064 triangles_[t].v[0] = v0;
766 3165064 triangles_[t].v[1] = v1;
767 3165064 triangles_[t].v[2] = v2;
768 3165064 triangles_[t].t[0] = t0;
769 3165064 triangles_[t].t[1] = t1;
770 3165064 triangles_[t].t[2] = t2;
771 3165064 return t;
772 }
773
774 /**
775 * \brief Sets the vertices and adjacent triangles of a triangle.
776 * \param[in] t index of the triangle
777 * \param[in] v0 index of first vertex
778 * \param[in] v1 index of second vertex
779 * \param[in] v2 index of third vertex
780 * \param[in] t0 index of adjacent triangle opposite to \p v0
781 * \param[in] t1 index of adjacent triangle opposite to \p v1
782 * \param[in] t2 index of adjacent triangle opposite to \p v2
783 */
784 void set_triangle(
785 index_t t,
786 index_t v0, index_t v1, index_t v2,
787 index_t t0, index_t t1, index_t t2
788 ) {
789 geo_debug_assert(t != NO_TRIANGLE);
790 geo_debug_assert(t < max_t());
791 triangles_[t].v[0] = v0;
792 triangles_[t].v[1] = v1;
793 triangles_[t].v[2] = v2;
794 triangles_[t].t[0] = t0;
795 triangles_[t].t[1] = t1;
796 triangles_[t].t[2] = t2;
797 }
798
799 /**
800 * \brief Creates a new triangles with specified vertices,
801 * adjacent triangles and geometric location at the dual
802 * vertex.
803 * \details The vertex is shared with caller.
804 * The created triangle is marked as used.
805 * \param[in] p geometric location at the dual vertex, shared with
806 * caller. Caller remains responsible for memory management.
807 * \param[in] w the weight associated with point \p p
808 * \param[in] v0 index of first vertex
809 * \param[in] v1 index of second vertex
810 * \param[in] v2 index of third vertex
811 * \param[in] t0 index of adjacent triangle opposite to \p v0
812 * \param[in] t1 index of adjacent triangle opposite to \p v1
813 * \param[in] t2 index of adjacent triangle opposite to \p v2
814 * \return the index of the new triangle
815 */
816 3165064 index_t create_triangle(
817 const double* p,
818 double w,
819 index_t v0, index_t v1, index_t v2,
820 index_t t0, index_t t1, index_t t2
821 ) {
822 3165064 index_t t = create_triangle(v0, v1, v2, t0, t1, t2);
823 3165064 triangle_dual(t).set_point(p);
824 3165064 triangle_dual(t).set_weight(w);
825 3165064 return t;
826 }
827
828 /**
829 * \brief Creates a new triangles with specified vertices,
830 * adjacent triangles and geometric location at the dual
831 * vertex.
832 * \details The vertex is copied into local storage.
833 * The created triangle is marked as used.
834 * \param[in] p geometric location at the dual vertex. Vertex
835 * is copied into local storage.
836 * \param[in] v0 index of first vertex
837 * \param[in] v1 index of second vertex
838 * \param[in] v2 index of third vertex
839 * \param[in] t0 index of adjacent triangle opposite to \p v0
840 * \param[in] t1 index of adjacent triangle opposite to \p v1
841 * \param[in] t2 index of adjacent triangle opposite to \p v2
842 * \return the index of the new triangle
843 */
844 index_t create_triangle_copy(
845 const double* p,
846 index_t v0, index_t v1, index_t v2,
847 index_t t0, index_t t1, index_t t2
848 ) {
849 index_t t = create_triangle(v0, v1, v2, t0, t1, t2);
850 double* np = intersections_.new_item();
851 for(coord_index_t c = 0; c < dimension(); ++c) {
852 np[c] = p[c];
853 }
854 triangle_dual(t).set_point(np);
855 return t;
856 }
857
858 /**
859 * \brief Creates a new vertex.
860 * \return the index of the new vertex
861 */
862 13373954 index_t create_vertex() {
863 13373954 v_to_t_dirty_ = true;
864
1/2
✓ Branch 2 taken 13373954 times.
✗ Branch 3 not taken.
 13373954 vertices_.push_back(Vertex());
865 13373954 return vertices_.size() - 1;
866 }
867
868 /**
869 * \brief Triangulates the conflict zone.
870 * \details Creates the triangles radiating from \p new_v and
871 * attached to the border of the conflict zone, indicated by \p t and \p e.
872 * \param[in] delaunay the Delaunay triangulation
873 * \param[in] i index of the first extremity of
874 * the bisector in \p delaunay.
875 * \param[in] j index of the first extremity of
876 * the bisector in \p delaunay.
877 * \param[in] symbolic if true, symbolic representation of the
878 * vertices is generated.
879 * \param[in] t1 a triangle adjacent to the border of the conflict zone from
880 * inside.
881 * \param[in] t1ebord the edge along which t is adjacent to the border of the conflict
882 * zone.
883 * \param[in] v_in index of the new vertex
884 * \return one of the created triangles
885 */
886 template <index_t DIM>
887 7926569 index_t triangulate_hole(
888 const Delaunay* delaunay,
889 index_t i, index_t j, bool symbolic,
890 index_t t1, index_t t1ebord,
891 index_t v_in
892 ) {
893 7926569 index_t t = t1;
894 7926569 index_t e = t1ebord;
895 7926569 index_t t_adj = triangle_adjacent(t,e);
896
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3966610 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 7926569 geo_debug_assert(t_adj != index_t(-1));
897
898
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3966610 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 7926569 geo_debug_assert(triangle_is_conflict(t));
899
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3966610 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 7926569 geo_debug_assert(!triangle_is_conflict(t_adj));
900
901 7926569 index_t new_t_first = index_t(-1);
902 7926569 index_t new_t_prev = index_t(-1);
903
904 do {
905
906 29034890 index_t v1 = triangle_vertex(t, plus1mod3(e));
907 29034890 index_t v2 = triangle_vertex(t, minus1mod3(e));
908
909 // Create new triangle
910 29034890 index_t new_t = create_triangle(v_in, v1, v2);
911
912 29034890 triangle_dual(new_t).intersect_geom<DIM>(
913 29034890 intersections_,
914 29034890 triangle_dual(t),
915 29034890 triangle_dual(triangle_adjacent(t, e)),
916 delaunay->vertex_ptr(i), delaunay->vertex_ptr(j)
917 );
918
919
2/2
✓ Branch 0 taken 385708 times.
✓ Branch 1 taken 14149911 times.
 29034890 if(symbolic) {
920 1494368 triangle_dual(new_t).sym().intersect_symbolic(
921 747184 triangle_dual(t).sym(),
922 747184 triangle_dual(triangle_adjacent(t, e)).sym(),
923 j
924 );
925 }
926
927 // Connect new triangle to triangle on the other
928 // side of the conflict zone.
929 29034890 set_triangle_adjacent(new_t, 0, t_adj);
930 29034890 index_t adj_e = triangle_adjacent_index(t_adj, t);
931 29034890 set_triangle_adjacent(t_adj, adj_e, new_t);
932
933
934 // Move to next triangle
935 29034890 e = plus1mod3(e);
936 29034890 t_adj = index_t(triangle_adjacent(t,e));
937
2/2
✓ Branch 1 taken 10352492 times.
✓ Branch 2 taken 14535619 times.
 49706904 while(triangle_is_conflict(t_adj)) {
938 20672014 t = t_adj;
939 20672014 e = minus1mod3(find_triangle_vertex(t,v2));
940 20672014 t_adj = index_t(triangle_adjacent(t,e));
941
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10352492 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 20672014 geo_debug_assert(t_adj != index_t(-1));
942 }
943
944
2/2
✓ Branch 0 taken 3966610 times.
✓ Branch 1 taken 10569009 times.
 29034890 if(new_t_prev == index_t(-1)) {
945 7926569 new_t_first = new_t;
946 } else {
947 21108321 set_triangle_adjacent(new_t_prev, 1, new_t);
948 21108321 set_triangle_adjacent(new_t, 2, new_t_prev);
949 }
950
951 29034890 new_t_prev = new_t;
952
953
4/4
✓ Branch 0 taken 6993547 times.
✓ Branch 1 taken 7542072 times.
✓ Branch 2 taken 3575462 times.
✓ Branch 3 taken 3966610 times.
 29034890 } while((t != t1) || (e != t1ebord));
954
955 // Connect last triangle to first triangle
956 7926569 set_triangle_adjacent(new_t_prev, 1, new_t_first);
957 7926569 set_triangle_adjacent(new_t_first, 2, new_t_prev);
958
959 7926569 return new_t_prev;
960 }
961
962 /**
963 * \brief Determines the conflict zone.
964 * \details The conflict zone corresponds to the set of triangles
965 * that have their dual vertices on the negative side of a
966 * bisector.
967 * \param[in] mesh the input mesh
968 * \param[in] delaunay the Delaunay triangulation
969 * \param[in] i index of the first extremity
970 * of the bisector in \p delaunay
971 * \param[in] j index of the second extremity
972 * of the bisector in \p delaunay
973 * \param[in] exact if true, exact predicates are used
974 * \param[out] conflict_begin
975 * index of the first triangle in conflict list
976 * \param[out] conflict_end one position past index of the
977 * last triangle in conflict list
978 */
979 template <index_t DIM>
980 20391200 void get_conflict_list(
981 const Mesh* mesh, const Delaunay* delaunay,
982 index_t i, index_t j, bool exact,
983 index_t& conflict_begin, index_t& conflict_end
984 ) {
985 20391200 conflict_begin = END_OF_LIST;
986 20391200 conflict_end = END_OF_LIST;
987
2/2
✓ Branch 0 taken 258485 times.
✓ Branch 1 taken 9950441 times.
 20391200 if(exact) {
988 // In exact mode, we classify each vertex
989 // using the exact predicate. Note that
990 // "climbing/walking" from a random vertex
991 // would be more efficient, but it would require a
992 // "comparison" exact predicate (with 8 different
993 // versions according to the configuration of the
994 // two vertices to be compared)
995
2/2
✓ Branch 1 taken 42869613 times.
✓ Branch 2 taken 258485 times.
 45794304 for(index_t t = 0; t < max_t(); t++) {
996
2/2
✓ Branch 1 taken 42079484 times.
✓ Branch 2 taken 790129 times.
 45301765 if(triangle_is_used(t)) {
997 131469978 Sign s = side<DIM>(
998 mesh, delaunay,
999 43823326 triangle_dual(t),
1000 i, j,
1001 exact
1002 );
1003
2/2
✓ Branch 0 taken 281000 times.
✓ Branch 1 taken 41798484 times.
 43823326 if(s == GEO::NEGATIVE) {
1004 546516 append_triangle_to_conflict_list(
1005 t, conflict_begin, conflict_end
1006 );
1007 }
1008 }
1009 }
1010 } else {
1011 // In non-exact mode, we first detect the
1012 // vertex that is furthest away along the
1013 // normal vector of the clipping bisector,
1014 // and then we propagate using a flood-fill
1015 // algorithm. This strategy ensures that
1016 // the conflict zone remains connected, even
1017 // in the presence of numerical errors. Clearly
1018 // it does not ensure validity, but in practice
1019 // it improves resistance to degeneracies.
1020 index_t furthest_t =
1021 19898661 find_furthest_point_linear_scan<DIM>(
1022 delaunay, i, j
1023 );
1024 19898661 propagate_conflict_list<DIM>(
1025 mesh, delaunay, furthest_t,
1026 i, j, exact,
1027 conflict_begin, conflict_end
1028 );
1029 }
1030 20391200 }
1031
1032 /**
1033 * \brief Finds the index of the vertex furthest away
1034 * on the negative side of a bisector.
1035 * \param[in] delaunay the Delaunay triangulation
1036 * \param[in] i index of the first extremity of the bisector
1037 * in \p delaunay
1038 * \param[in] j index of the second extremity of the bisector
1039 * in \p delaunay
1040 * \return the index of the vertex furthest away on \p j%'s side,
1041 * or -1 if all the vertices are on \p i%'s side.
1042 */
1043 template <index_t DIM>
1044 19898661 index_t find_furthest_point_linear_scan(
1045 const Delaunay* delaunay, index_t i, index_t j
1046 ) const {
1047 19898661 index_t result = NO_TRIANGLE;
1048 19898661 double furthest_dist = 0.0;
1049
3/4
✓ Branch 1 taken 117385485 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 107435044 times.
✓ Branch 4 taken 9950441 times.
 232743613 for(index_t t = 0; t < max_t(); ++t) {
1050
3/4
✓ Branch 1 taken 107435044 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 77847558 times.
✓ Branch 4 taken 29587486 times.
 212844952 if(triangle_is_used(t)) {
1051
1/2
✓ Branch 1 taken 77847558 times.
✗ Branch 2 not taken.
 153678918 double d = signed_bisector_distance<DIM>(
1052
1/2
✓ Branch 1 taken 77847558 times.
✗ Branch 2 not taken.
 153678918 delaunay, i, j, triangle_dual(t).point()
1053 );
1054
2/2
✓ Branch 0 taken 6129240 times.
✓ Branch 1 taken 71718318 times.
 153678918 if(d < furthest_dist) {
1055 12254397 result = t;
1056 12254397 furthest_dist = d;
1057 }
1058 }
1059 }
1060
2/2
✓ Branch 0 taken 3864729 times.
✓ Branch 1 taken 6085712 times.
 39797322 return (furthest_dist < 0) ? result : NO_TRIANGLE;
1061 }
1062
1063 /**
1064 * \brief Evaluates the equation of a bisector at a given point.
1065 * \details Positive side corresponds to vertex \p i and negative
1066 * side to vertex \p j.
1067 * \param[in] delaunay the Delaunay triangulation
1068 * \param[in] i index of the first extremity of the bisector
1069 * in \p delaunay
1070 * \param[in] j index of the second extremity of the bisector
1071 * in \p delaunay
1072 * \param[in] q the query point
1073 */
1074 template <index_t DIM>
1075 153678918 static double signed_bisector_distance(
1076 const Delaunay* delaunay, index_t i, index_t j, const double* q
1077 ) {
1078 153678918 const double* pi = delaunay->vertex_ptr(i);
1079 153678918 const double* pj = delaunay->vertex_ptr(j);
1080 153678918 double result = 0;
1081
2/2
✓ Branch 0 taken 423405664 times.
✓ Branch 1 taken 77847558 times.
 994441652 for(coord_index_t c = 0; c < DIM; ++c) {
1082 840762734 result += GEO::geo_sqr(q[c] - pj[c]);
1083 840762734 result -= GEO::geo_sqr(q[c] - pi[c]);
1084 }
1085 153678918 return result;
1086 }
1087
1088 /**
1089 * \brief Computes the conflict list by propagation from
1090 * a conflict triangle.
1091 * \param[in] mesh the input mesh
1092 * \param[in] delaunay the Delaunay triangulation
1093 * \param[in] first_t a triangle in the conflict zone
1094 * \param[in] i index of the first extremity of the bisector
1095 * in \p delaunay
1096 * \param[in] j index of the second extremity of the bisector
1097 * in \p delaunay
1098 * \param[in] exact if true, exact predicates are used
1099 * \param[out] conflict_begin
1100 * index of the first triangle in conflict list
1101 * \param[out] conflict_end one position past index of the
1102 * last triangle in conflict list
1103 */
1104 template <index_t DIM>
1105 19898661 void propagate_conflict_list(
1106 const Mesh* mesh, const Delaunay* delaunay,
1107 index_t first_t,
1108 index_t i, index_t j, bool exact,
1109 index_t& conflict_begin, index_t& conflict_end
1110 ) {
1111 19898661 conflict_begin = END_OF_LIST;
1112 19898661 conflict_end = END_OF_LIST;
1113
1114 // Special case, clipping plane does not clip anything
1115
2/2
✓ Branch 0 taken 6085712 times.
✓ Branch 1 taken 3864729 times.
 19898661 if(first_t == NO_TRIANGLE) {
1116 12171420 return;
1117 }
1118
1119
1/2
✓ Branch 1 taken 3864729 times.
✗ Branch 2 not taken.
 7727241 std::stack<index_t> S;
1120
1/2
✓ Branch 1 taken 3864729 times.
✗ Branch 2 not taken.
 7727241 S.push(first_t);
1121
1/2
✓ Branch 1 taken 3864729 times.
✗ Branch 2 not taken.
 7727241 append_triangle_to_conflict_list(
1122 first_t, conflict_begin, conflict_end
1123 );
1124
2/2
✓ Branch 1 taken 10354295 times.
✓ Branch 2 taken 3864729 times.
 28428089 while(!S.empty()) {
1125 20700848 index_t t = S.top();
1126 20700848 S.pop();
1127
2/2
✓ Branch 0 taken 31062885 times.
✓ Branch 1 taken 10354295 times.
 82803392 for(unsigned int e = 0; e < 3; e++) {
1128
1/2
✓ Branch 1 taken 31062885 times.
✗ Branch 2 not taken.
 62102544 index_t neigh = index_t(triangle_adjacent(t, e));
1129
3/4
✓ Branch 1 taken 31062885 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 20639477 times.
✓ Branch 4 taken 10423408 times.
 62102544 if(!triangle_is_conflict(neigh)) {
1130 41261313 if(
1131 123783939 side<DIM>(
1132
2/4
✓ Branch 1 taken 20639477 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20639477 times.
✗ Branch 5 not taken.
 41261313 mesh, delaunay, triangle_dual(neigh),
1133 i, j, exact
1134
2/2
✓ Branch 0 taken 6489566 times.
✓ Branch 1 taken 14149911 times.
 41261313 ) == GEO::NEGATIVE
1135 ) {
1136
1/2
✓ Branch 1 taken 6489566 times.
✗ Branch 2 not taken.
 12973607 S.push(neigh);
1137
1/2
✓ Branch 1 taken 6489566 times.
✗ Branch 2 not taken.
 12973607 append_triangle_to_conflict_list(
1138 neigh, conflict_begin, conflict_end
1139 );
1140 }
1141 }
1142 }
1143 }
1144 7727241 }
1145
1146 /**
1147 * \brief Tests on which side a vertex is relative to
1148 * a bisector.
1149 * \param[in] mesh the input mesh
1150 * \param[in] delaunay the Delaunay triangulation
1151 * \param[in] v the query vertex
1152 * \param[in] i index of the first extremity of the bisector
1153 * in \p delaunay
1154 * \param[in] j index of the second extremity of the bisector
1155 * in \p delaunay
1156 * \param[in] exact if true, exact predicates are used
1157 * \return POSITIVE if \p v is on vertex \p i%'s side,
1158 * NEGATIVE otherwise. ZERO is never returned since
1159 * globally coherent symbolic perturbations are used
1160 * in exact mode.
1161 * \tparam DIM dimension, specified as a template
1162 * parameter for efficiency considerations.
1163 */
1164 template <index_t DIM>
1165 85084639 Sign side(
1166 const Mesh* mesh, const Delaunay* delaunay,
1167 const GEOGen::Vertex& v,
1168 index_t i, index_t j, bool exact
1169 ) const {
1170 85084639 Sign result = GEO::ZERO;
1171
2/2
✓ Branch 0 taken 42079484 times.
✓ Branch 1 taken 20639477 times.
 85084639 if(exact) {
1172 43823326 result = side_exact(
1173 mesh, delaunay, v,
1174 delaunay->vertex_ptr(i),
1175 delaunay->vertex_ptr(j),
1176 DIM,
1177 43823326 symbolic_is_surface_
1178 );
1179 } else {
1180 41261313 result = v.side_fast<DIM>(
1181 delaunay->vertex_ptr(i),
1182 delaunay->vertex_ptr(j)
1183 );
1184 }
1185 85084639 return result;
1186 }
1187
1188 /**
1189 * \brief Tests on which side a vertex is relative to
1190 * a bisector using exact predicates.
1191 * \param[in] mesh the input mesh
1192 * \param[in] delaunay the Delaunay triangulation
1193 * \param[in] v the query vertex
1194 * \param[in] pi first extremity of the bisector
1195 * \param[in] pj second extremity of the bisector
1196 * \param[in] dim dimension of the points
1197 * \param[in] symbolic_is_surface if true, then symbolic
1198 * information is relative to a surface mesh (facets)
1199 * rather than volumetric mesh (tetrahedra).
1200 * \return POSITIVE if \p v is on vertex \p i%'s side,
1201 * NEGATIVE otherwise. ZERO is never returned since
1202 * globally coherent symbolic perturbations are used
1203 * in exact mode.
1204 * \note Only dimension=3 is implemented for now
1205 */
1206 Sign side_exact(
1207 const Mesh* mesh, const Delaunay* delaunay,
1208 const GEOGen::Vertex& v,
1209 const double* pi, const double* pj,
1210 coord_index_t dim,
1211 bool symbolic_is_surface = false
1212 ) const;
1213
1214 /**
1215 * \brief Gets a triangle and an edge on the internal border of the conflict zone.
1216 * \details The returned triangle touches the conflict zone from inside.
1217 * \param[in] conflict_begin first triangle of the conflict zone
1218 * \param[in] conflict_end one element past the last triangle of
1219 * the conflict zone
1220 * \param[out] t a triangle in the conflict zone adjacent to the border of the
1221 * conflict zone.
1222 * \param[out] e the edge along which \p t is adjacent to the border of the
1223 * conflict zone.
1224 * \return true if a triangle on the border was found, false otherwise.
1225 */
1226 3967975 bool find_triangle_on_border(
1227 index_t conflict_begin, index_t conflict_end,
1228 index_t& t, index_t& e
1229 ) const {
1230 3967975 GEO::geo_argused(conflict_end);
1231 3967975 t = conflict_begin;
1232 do {
1233
2/2
✓ Branch 0 taken 5612869 times.
✓ Branch 1 taken 39539 times.
 5652408 for(e = 0; e < 3; ++e) {
1234 5612869 index_t nt = triangle_adjacent(t, e);
1235
2/2
✓ Branch 1 taken 3966610 times.
✓ Branch 2 taken 1646259 times.
 5612869 if(triangle_is_used(nt)) {
1236 3966610 return true;
1237 }
1238 }
1239 39539 t = next_triangle(t);
1240
2/2
✓ Branch 0 taken 38174 times.
✓ Branch 1 taken 1365 times.
 39539 } while(t != END_OF_LIST);
1241 1365 return false;
1242 }
1243
1244 /**
1245 * \brief Gets the successor of a triangle.
1246 * \details Triangles are linked, for instance to represent
1247 * the conflict zone.
1248 */
1249 24400189 index_t next_triangle(index_t t) const {
1250
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 24400189 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 24400189 geo_debug_assert(t != NO_TRIANGLE);
1251
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 24400189 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 24400189 geo_debug_assert(t < max_t());
1252 24400189 return triangles_[t].next_;
1253 }
1254
1255 /**
1256 * \brief Sets the successor of a triangle.
1257 * \details Triangles are linked, for instance to represent
1258 * the conflict zone.
1259 * \param[in] t index of the triangle
1260 * \param[in] t2 index of the successor
1261 */
1262 14601905 void set_next_triangle(index_t t, index_t t2) {
1263
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 14601905 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 14601905 geo_debug_assert(t != NO_TRIANGLE);
1264
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 14601905 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 14601905 geo_debug_assert(t < max_t());
1265 14601905 triangles_[t].next_ = t2;
1266 14601905 }
1267
1268 /**
1269 * \brief Specify that a triangle is free.
1270 * \details A free triangle can be reused by subsequent
1271 * triangle creations.
1272 */
1273 10626577 void mark_as_free(index_t t) {
1274
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10626577 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 10626577 geo_debug_assert(t != NO_TRIANGLE);
1275
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10626577 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10626577 geo_debug_assert(t < max_t());
1276 10626577 triangles_[t].status_ = TRI_IS_FREE;
1277 10626577 }
1278
1279 /**
1280 * \brief Specify that a triangle belongs to the conflict zone.
1281 */
1282 10635295 void mark_as_conflict(index_t t) {
1283
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10635295 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 10635295 geo_debug_assert(t != NO_TRIANGLE);
1284
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10635295 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10635295 geo_debug_assert(t < max_t());
1285 10635295 triangles_[t].status_ = TRI_IS_CONFLICT;
1286 10635295 }
1287
1288 /**
1289 * \brief Specify that a triangle is used.
1290 */
1291 17700683 void mark_as_used(index_t t) {
1292
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 17700683 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 17700683 geo_debug_assert(t != NO_TRIANGLE);
1293
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 17700683 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 17700683 geo_debug_assert(t < max_t());
1294 17700683 triangles_[t].status_ = TRI_IS_USED;
1295 17700683 }
1296
1297 /**
1298 * \brief Appends a triangle to the conflict list.
1299 * \details The triangle is marked as conflict and
1300 * linked to the conflict list.
1301 * \param[in] t the triangle
1302 * \param[in,out] conflict_begin first triangle in the conflict list
1303 * \param[in,out] conflict_end one position past the last triangle
1304 * in the conflict list
1305 */
1306 10635295 void append_triangle_to_conflict_list(
1307 index_t t, index_t& conflict_begin, index_t& conflict_end
1308 ) {
1309
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10635295 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10635295 geo_debug_assert(triangle_is_used(t));
1310 10635295 set_next_triangle(t, conflict_begin);
1311 10635295 mark_as_conflict(t);
1312 10635295 conflict_begin = t;
1313
2/2
✓ Branch 0 taken 3967975 times.
✓ Branch 1 taken 6667320 times.
 10635295 if(conflict_end == END_OF_LIST) {
1314 3967975 conflict_end = t;
1315 }
1316 10635295 }
1317
1318 /**
1319 * \brief Merges a list of triangles into the free list.
1320 * \param[in] list_begin first triangle in the list to be freed
1321 * \param[in] list_end one position past the last triangle of
1322 * the list to be freed
1323 */
1324 3966610 void merge_into_free_list(index_t list_begin, index_t list_end) {
1325
1/2
✓ Branch 0 taken 3966610 times.
✗ Branch 1 not taken.
 3966610 if(list_begin != END_OF_LIST) {
1326
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3966610 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 3966610 geo_debug_assert(list_end != END_OF_LIST);
1327
1328 3966610 index_t cur = list_begin;
1329
2/2
✓ Branch 0 taken 6659967 times.
✓ Branch 1 taken 3966610 times.
 10626577 while(cur != list_end) {
1330 6659967 mark_as_free(cur);
1331 6659967 cur = next_triangle(cur);
1332 }
1333 3966610 mark_as_free(list_end);
1334 3966610 set_next_triangle(list_end, first_free_);
1335 3966610 first_free_ = list_begin;
1336 }
1337 3966610 }
1338
1339 /**
1340 * \brief Allocates a new triangle.
1341 * \details This function is called whenever a triangle needs
1342 * to be created and the free list is empty.
1343 */
1344 10002420 void grow() {
1345
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10002420 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 10002420 geo_debug_assert(first_free_ == END_OF_LIST);
1346 10002420 first_free_ = triangles_.size();
1347
2/4
✓ Branch 1 taken 10002420 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10002420 times.
✗ Branch 5 not taken.
 10002420 triangles_.push_back(Triangle());
1348 10002420 }
1349
1350 /**
1351 * \brief Displays the number of free,used,conflict triangles.
1352 * \details For debugging purposes.
1353 */
1354 std::ostream& show_stats(std::ostream& os) const;
1355
1356 /**
1357 * \brief Computes a unique global facet id from a mesh
1358 * tetrahedron and local facet index.
1359 * \details If a facet is shared by two tetrahedra t1 and t2,
1360 * its global index determined from t1 and from t2 is the same.
1361 * \param[in] mesh the mesh
1362 * \param[in] t the index of the tetrahedron
1363 * \param[in] lf the local facet index (0,1,2 or 3) in tetrahedron \p t
1364 * \return an index that uniquely identifies the facet
1365 * in the tetrahedron
1366 */
1367 78668 static index_t global_facet_id(
1368 const Mesh* mesh, index_t t, index_t lf
1369 ) {
1370 78668 index_t t2 = mesh->cells.tet_adjacent(t, lf);
1371
4/4
✓ Branch 0 taken 70965 times.
✓ Branch 1 taken 7703 times.
✓ Branch 2 taken 35598 times.
✓ Branch 3 taken 35367 times.
 78668 if(t2 != GEO::NO_CELL && t2 > t) {
1372 35598 index_t lf2 = mesh->cells.find_tet_adjacent(
1373 t2, t
1374 );
1375
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 35598 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 35598 geo_debug_assert(lf2 != GEO::NO_FACET);
1376 35598 return index_t(4 * t2 + lf2);
1377 }
1378 43070 return 4 * t + lf;
1379 }
1380
1381 private:
1382 GEO::vector<Triangle> triangles_;
1383 GEO::vector<Vertex> vertices_;
1384 index_t first_free_;
1385 bool v_to_t_dirty_;
1386 PointAllocator intersections_;
1387 bool symbolic_is_surface_;
1388 signed_index_t cell_id_;
1389
1390 static index_t plus1mod3_[3];
1391 static index_t minus1mod3_[3];
1392
1393 /**
1394 * \brief Gets the modulo-3 successor of an index.
1395 * \param[in] i the index
1396 * \return \p i plus 1 modulo 3
1397 */
1398 43589230 static index_t plus1mod3(index_t i) {
1399
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 43589230 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 43589230 geo_debug_assert(i < 3);
1400 43589230 return plus1mod3_[i];
1401 }
1402
1403 /**
1404 * \brief Gets the modulo-3 predecessor of an index.
1405 * \param[in] i the index
1406 * \return \p i minus 1 modulo 3
1407 */
1408 24888111 static index_t minus1mod3(index_t i) {
1409
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 24888111 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 24888111 geo_debug_assert(i < 3);
1410 24888111 return minus1mod3_[i];
1411 }
1412 };
1413 }
1414
1415 #endif
1416