GCC Code Coverage Report

Directory: ./
File: delaunay/periodic_delaunay_3d.cpp
Date: 2026-09-27 03:10:11
Exec Total Coverage
Lines: 927 1193 77.7%
Functions: 52 66 78.8%
Branches: 689 1312 52.5%
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 #include <geogram/delaunay/periodic_delaunay_3d.h>
41 #include <geogram/delaunay/cavity.h>
42
43 #include <geogram/mesh/mesh_reorder.h>
44 #include <geogram/numerics/predicates.h>
45 #include <geogram/basic/geometry.h>
46 #include <geogram/basic/stopwatch.h>
47 #include <geogram/basic/command_line.h>
48 #include <geogram/basic/permutation.h>
49 #include <geogram/basic/algorithm.h>
50 #include <geogram/bibliography/bibliography.h>
51
52 #include <stack>
53 #include <algorithm>
54
55 #include <mutex>
56 #include <condition_variable>
57
58 // ParallelDelaunayThread class, declared locally, has
59 // no out-of-line virtual functions. It is not a
60 // problem since they are only visible from this translation
61 // unit, but clang will complain.
62 #ifdef __clang__
63 #pragma GCC diagnostic ignored "-Wweak-vtables"
64 #endif
65
66 namespace {
67
68 using namespace GEO;
69
70 /**
71 * \brief Generates a random integer.
72 * \return a random integer between 0 and \p choices - 1
73 * \param [in] choices_in number of possible choices for the
74 * random variable (maximum value + 1)
75 * \details The function is thread-safe, and uses one seed
76 * per thread.
77 */
78 5 GEO::index_t thread_safe_random_(GEO::index_t choices_in) {
79 #ifdef GARGANTUA
80 typedef Numeric::int64 Int;
81 #else
82 typedef long int Int;
83 #endif
84 5 GEO::signed_index_t choices = signed_index_t(choices_in);
85 static thread_local Int randomseed = 1l ;
86
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if (choices >= 714025l) {
87 ✗ Int newrandom = (randomseed * 1366l + 150889l) % 714025l;
88 ✗ randomseed = (newrandom * 1366l + 150889l) % 714025l;
89 ✗ newrandom = newrandom * (choices / 714025l) + randomseed;
90 ✗ if (newrandom >= choices) {
91 ✗ return GEO::index_t(newrandom - choices);
92 } else {
93 ✗ return GEO::index_t(newrandom);
94 }
95 } else {
96 5 randomseed = (randomseed * 1366l + 150889l) % 714025l;
97 5 return GEO::index_t(randomseed % choices);
98 }
99 }
100
101 /**
102 * \brief Generates a random integer between 0 and 3.
103 * \return a random integer between 0 and 3
104 * \details The function is thread-safe, and uses one seed
105 * per thread.
106 */
107 GEO::index_t thread_safe_random_4_() {
108 static thread_local long int randomseed = 1l ;
109 32558 randomseed = (randomseed * 1366l + 150889l) % 714025l;
110 32558 return GEO::index_t(randomseed % 4);
111 }
112
113 /**
114 * \brief Computes the number of bits set.
115 * \param[in] x an integer.
116 * \return the number of ones in the binary representation
117 * of the integer.
118 */
119 inline VBW::index_t pop_count(Numeric::uint32 x) {
120 #if defined(GEO_COMPILER_GCC_FAMILY)
121 ✗ return VBW::index_t(Numeric::uint32(__builtin_popcount(x)));
122 #elif defined(GEO_COMPILER_MSVC)
123 #if defined(_M_ARM64)
124 return VBW::index_t(_CountOneBits(x));
125 #else
126 return VBW::index_t(__popcnt(x));
127 #endif
128 #else
129 int result = 0;
130 for(int b=0; b<32; ++b) {
131 result += ((x & 1) != 0);
132 x >>= 1;
133 }
134 return VBW::index_t(result);
135 #endif
136 }
137
138 /************************************************************************/
139
140 // TODO: move these two functions to mesh_reorder.h
141
142 5 void compute_BRIO_order_periodic_recursive(
143 GEO::index_t nb_vertices, const double* vertices,
144 GEO::index_t dimension, GEO::index_t stride,
145 vector<GEO::index_t>& sorted_indices,
146 vector<GEO::index_t>::iterator b,
147 vector<GEO::index_t>::iterator e,
148 const vec3& period,
149 GEO::index_t threshold,
150 double ratio,
151 GEO::index_t& depth,
152 vector<GEO::index_t>* levels
153 ) {
154 geo_debug_assert(e > b);
155
156
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
 5 vector<GEO::index_t>::iterator m = b;
157
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
 5 if(GEO::index_t(e - b) > threshold) {
158 3 ++depth;
159 3 m = b + int(double(e - b) * ratio);
160 3 compute_BRIO_order_periodic_recursive(
161 nb_vertices, vertices,
162 dimension, stride,
163 sorted_indices, b, m,
164 period,
165 threshold, ratio, depth,
166 levels
167 );
168 }
169
170 5 Hilbert_sort_periodic(
171 nb_vertices, vertices,
172 sorted_indices,
173 dimension,
174 stride,
175 m,
176 e,
177 period
178 );
179
180
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if(levels != nullptr) {
181 5 levels->push_back(GEO::index_t(e - sorted_indices.begin()));
182 }
183 5 }
184
185 2 void compute_BRIO_order_periodic(
186 GEO::index_t nb_vertices, const double* vertices,
187 GEO::index_t dimension, GEO::index_t stride,
188 vector<GEO::index_t>& sorted_indices,
189 vector<GEO::index_t>::iterator first,
190 vector<GEO::index_t>::iterator last,
191 const vec3& period,
192 GEO::index_t threshold = 64,
193 double ratio = 0.125,
194 vector<GEO::index_t>* levels = nullptr
195 ) {
196
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(levels != nullptr) {
197 levels->clear();
198 2 levels->push_back(GEO::index_t(first - sorted_indices.begin()));
199 }
200 2 GEO::index_t depth = 0;
201 2 GEO::random_shuffle(first, last);
202
203 2 compute_BRIO_order_periodic_recursive(
204 nb_vertices, vertices,
205 dimension, stride,
206 sorted_indices,
207 first, last,
208 period, threshold, ratio, depth, levels
209 );
210 2 }
211
212 /************************************************************************/
213
214 2 void delaunay_citations() {
215 2 geo_cite_with_info(
216 "DBLP:journals/cj/Bowyer81",
217 "One of the two initial references to the algorithm, "
218 "discovered independently and simultaneously by Bowyer and Watson."
219 );
220 2 geo_cite_with_info(
221 "journals/cj/Watson81",
222 "One of the two initial references to the algorithm, "
223 "discovered independently and simultaneously by Bowyer and Watson."
224 );
225 2 geo_cite_with_info(
226 "DBLP:conf/compgeom/AmentaCR03",
227 "Using spatial sorting has a dramatic impact on the performances."
228 );
229 2 geo_cite_with_info(
230 "DBLP:journals/comgeo/FunkeMN05",
231 "Initializing \\verb|locate()| with a non-exact version "
232 " (structural filtering) gains (a bit of) performance."
233 );
234 2 geo_cite_with_info(
235 "DBLP:journals/comgeo/BoissonnatDPTY02",
236 "The idea of traversing the cavity from inside "
237 " used in GEOGRAM is inspired by the implementation of "
238 " \\verb|Delaunay_triangulation_3| in CGAL."
239 );
240 2 geo_cite_with_info(
241 "DBLP:conf/imr/Si06",
242 "The triangulation data structure used in GEOGRAM is inspired "
243 "by Tetgen."
244 );
245 2 geo_cite_with_info(
246 "DBLP:journals/ijfcs/DevillersPT02",
247 "Analysis of the different versions of the line walk algorithm "
248 " used by \\verb|locate()|."
249 );
250 2 }
251 }
252
253 /************************************************************************/
254
255 // These two functions are missing when compiling in PSM mode.
256 #ifdef GEOGRAM_PSM
257 namespace GEO {
258 namespace PCK {
259 inline Sign det_3d(const vec3& p0, const vec3& p1, const vec3& p2) {
260 return det_3d(p0.data(), p1.data(), p2.data());
261 }
262 inline Sign det_4d(
263 const vec4& p0, const vec4& p1, const vec4& p2, const vec4& p3
264 ) {
265 return det_4d(p0.data(), p1.data(), p2.data(), p3.data());
266 }
267 }
268 }
269 #endif
270
271 /************************************************************************/
272
273 namespace GEO {
274
275 /**
276 * \brief One of the threads of the multi-threaded
277 * 3d Delaunay implementation.
278 */
279 class PeriodicDelaunay3dThread : public GEO::Thread, public Periodic {
280 public:
281 friend class PeriodicDelaunay3d;
282
283 /**
284 * \brief Symbolic value for cell_status_[t] that
285 * indicates that no thread owns t.
286 */
287 static constexpr index_t NO_THREAD = CellStatusArray::FREE_CELL;
288
289 /**
290 * \brief Creates a new PeriodicDelaunay3dThread.
291 * \details Each PeriodicDelaunay3dThread has an affected working
292 * zone, i.e. a range of tetrahedra indices in which the
293 * thread is allowed to create tetrahedra.
294 * \param[in] master a pointer to the PeriodicDelaunay3d
295 * this thread belongs to
296 * \param[in] pool_begin first tetrahedron index of
297 * the working zone of this PeriodicDelaunay3dThread
298 * \param[in] pool_end one position past the last tetrahedron
299 * index of the working zone of this PeriodicDelaunay3dThread
300 */
301 8 PeriodicDelaunay3dThread(
302 PeriodicDelaunay3d* master,
303 index_t pool_begin,
304 index_t pool_end
305 8 ) :
306 8 master_(master),
307 8 periodic_(master->periodic_),
308 8 period_(master->period_),
309 8 cell_to_v_store_(master_->cell_to_v_store_),
310 8 cell_to_cell_store_(master_->cell_to_cell_store_),
311 8 cell_next_(master_->cell_next_),
312 8 cell_status_(master_->cell_status_),
313 8 abort_on_empty_cell_(master->abort_on_empty_cell_),
314
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
 8 has_empty_cells_(false) {
315
316
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 max_t_ = master_->cell_next_.size();
317
318 8 nb_vertices_ = master_->nb_vertices();
319 8 nb_vertices_non_periodic_ = master_->nb_vertices_non_periodic_;
320 8 vertices_ = master_->vertex_ptr(0);
321 8 weights_ = master_->weights_;
322 8 dimension_ = master_->dimension();
323 8 reorder_ = master_->reorder_.data();
324
325 8 b_hint_ = NO_TETRAHEDRON;
326 8 e_hint_ = NO_TETRAHEDRON;
327
328 8 nb_rollbacks_ = 0;
329 8 nb_failed_locate_ = 0;
330
331 8 set_pool(pool_begin, pool_end);
332 8 }
333
334 /**
335 * \brief Resets thread statistics
336 */
337 void reset_stats() {
338 16 nb_rollbacks_ = 0;
339 16 nb_failed_locate_ = 0;
340 }
341
342 /**
343 * \brief Initializes the pool of tetrahedra for this thread.
344 * \param[in] pool_begin first tetrahedron index of
345 * the working zone of this PeriodicDelaunay3dThread
346 * \param[in] pool_end one position past the last tetrahedron
347 * index of the working zone of this PeriodicDelaunay3dThread
348 */
349 8 void set_pool(index_t pool_begin, index_t pool_end) {
350 8 pool_begin_ = pool_begin;
351 8 pool_end_ = pool_end;
352 // Initialize free list in memory pool
353 8 first_free_ = pool_begin;
354
2/2
✓ Branch 0 taken 3395 times.
✓ Branch 1 taken 8 times.
 3403 for(index_t t=pool_begin; t<pool_end-1; ++t) {
355 3395 cell_next_[t] = t+1;
356 }
357 8 cell_next_[pool_end-1] = END_OF_LIST;
358 8 nb_free_ = pool_end - pool_begin;
359 8 memory_overflow_ = false;
360 8 work_begin_ = NO_INDEX;
361 8 work_rbegin_ = NO_INDEX;
362 8 finished_ = false;
363 8 direction_ = true;
364 #ifdef GEO_DEBUG
365 nb_acquired_tets_ = 0;
366 #endif
367 8 interfering_thread_ = NO_THREAD;
368 8 nb_tets_to_create_ = 0;
369 8 t_boundary_ = NO_TETRAHEDRON;
370 8 f_boundary_ = NO_INDEX;
371 8 used_tets_end_ = pool_begin;
372 8 }
373
374 /**
375 * \brief Tests whether this thread created empty cells.
376 * \retval true if this thread created empty cells.
377 * \retval false otherwise.
378 */
379 bool has_empty_cells() {
380 4 return has_empty_cells_;
381 }
382
383 /**
384 * \brief Picks a random tetrahedron in this thread's pool.
385 * \retval If no valid tet exists in this thread's pool,
386 * returns NO_TETRAHEDRON. It can be a real tetrahedron
387 * \retval Otherwise, returns a finite tetrahedon, an infinite
388 * tetrahedron or a tetrahedon in the free list. Caller needs to check.
389 */
390 index_t pick_random_tet() const {
391 // Shit happens [Forrest Gump]
392
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if(used_tets_end_ == pool_begin_) {
393 return NO_TETRAHEDRON;
394 }
395 5 return pool_begin_ + thread_safe_random_(
396 used_tets_end_ - pool_begin_
397 5 );
398 }
399
400 /**
401 * \brief Gets the number of rollbacks.
402 * \return the number of rollbacks
403 * \details rollbacks occur whenever a point
404 * could not be inserted, due to interferences
405 * from other threads
406 */
407 index_t nb_rollbacks() const {
408
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 return nb_rollbacks_;
409 }
410
411 /**
412 * \brief Gets the number of failed locate() calls.
413 * \return the number of failed locate() calls
414 * \details locate() can fail when it cannot acquire
415 * the tetrahedra that are traversed, due to
416 * interferences from another thread.
417 */
418 index_t nb_failed_locate() const {
419
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 return nb_failed_locate_;
420 }
421
422 /**
423 * \brief Gets the number of tetrahedra traversed by
424 * the latest locate() invocation.
425 */
426 index_t nb_traversed_tets() const {
427 return nb_traversed_tets_;
428 }
429
430 /**
431 * \brief Sets the point index sequence that
432 * should be processed by this thread.
433 * \param[in] b index of the first point to insert
434 * \param[in] e one position past the index of the
435 * last point to insert
436 */
437 void set_work(index_t b, index_t e) {
438 4 work_begin_ = b;
439 // e is one position past the last point index
440 // to insert. Internally we store the last point index
441 // to insert (like rbegin in STL containers). This is
442 // because we manipulate the point sequence to insert
443 // from both ends.
444 4 work_rbegin_ = e-1;
445 // reorder_ may have changed if new vertices were
446 // inserted into it
447 ✗ reorder_ = master_->reorder_.data();
448 }
449
450 /**
451 * \brief Gets the number of remaining points to
452 * be inserted.
453 * \return the number of points to be inserted by
454 * this thread
455 */
456 index_t work_size() const {
457 ✗ if(work_begin_ == NO_INDEX && work_rbegin_ == NO_INDEX) {
458 return 0;
459 }
460 geo_debug_assert(work_begin_ != NO_INDEX);
461 geo_debug_assert(work_rbegin_ != NO_INDEX);
462 ✗ return std::max(work_rbegin_ - work_begin_ + 1, index_t(0));
463 }
464
465 /**
466 * \brief Gets the number of threads.
467 * \return the number of threads created by
468 * the master PeriodicDelaunay3d of this thread.
469 */
470 index_t nb_threads() const {
471 return index_t(master_->threads_.size());
472 }
473
474 /**
475 * \brief Gets a thread by index
476 * \pre t < nb_threads()
477 * \param[in] t index of the thread
478 * \return a poiner to the \p t th thread
479 */
480 PeriodicDelaunay3dThread* thread(index_t t) {
481 return static_cast<PeriodicDelaunay3dThread*>(
482 ✗ master_->threads_[t].get()
483 );
484 }
485
486 /**
487 * \brief Inserts the point sequence allocated to
488 * this thread.
489 * \details The point sequence was previously defined
490 * by set_work().
491 */
492 4 void run() override {
493 4 has_empty_cells_ = false;
494 4 finished_ = false;
495
496
2/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 if(work_begin_ == NO_INDEX || work_rbegin_ == NO_INDEX) {
497 return ;
498 }
499
500 4 memory_overflow_ = false;
501
502 // Current hint associated with b
503 4 b_hint_ = NO_TETRAHEDRON;
504
505 // Current hint associated with e
506 4 e_hint_ = NO_TETRAHEDRON;
507
508 // If true, insert in b->e order,
509 // else insert in e->b order
510 4 direction_ = true;
511
512 4 while(
513 8042 work_rbegin_ >= work_begin_ &&
514
3/4
✓ Branch 0 taken 4019 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4019 times.
✗ Branch 3 not taken.
 4023 !memory_overflow_ && (
515
1/2
✓ Branch 0 taken 4019 times.
✗ Branch 1 not taken.
 4019 !abort_on_empty_cell_ || (
516 ✗ !has_empty_cells_ &&
517 ✗ !master_->has_empty_cells_
518 )
519 )
520 ) {
521
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4019 times.
 4019 index_t v = direction_ ? work_begin_ : work_rbegin_ ;
522
1/2
✓ Branch 0 taken 4019 times.
✗ Branch 1 not taken.
 4019 index_t& hint = direction_ ? b_hint_ : e_hint_ ;
523
524 // Try to insert v and update hint
525 4019 bool success = insert(reorder_[v],hint);
526
527 // Notify all threads that are waiting for
528 // this thread to release some tetrahedra.
529 send_event();
530
531
1/2
✓ Branch 0 taken 4019 times.
✗ Branch 1 not taken.
 4019 if(success) {
532
1/2
✓ Branch 0 taken 4019 times.
✗ Branch 1 not taken.
 4019 if(direction_) {
533 4019 ++work_begin_;
534 } else {
535 ✗ --work_rbegin_;
536 }
537 } else {
538 ✗ ++nb_rollbacks_;
539 ✗ if(interfering_thread_ != NO_THREAD) {
540 ✗ if(id() < interfering_thread_) {
541 // If this thread has a higher priority than
542 // the one that interfered, wait for the
543 // interfering thread to release the tets that
544 // it holds (then the loop will retry to insert
545 // the same vertex).
546 ✗ wait_for_event(interfering_thread_);
547 } else {
548 // If this thread has a lower priority than
549 // the interfering thread, try inserting
550 // from the other end of the points sequence.
551 ✗ direction_ = !direction_;
552 }
553 }
554 }
555 }
556
557
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if(has_empty_cells_) {
558 ✗ master_->has_empty_cells_ = true;
559 }
560
561 // Fix by Hiep Vu: wake up threads that potentially missed
562 // the previous wake ups.
563 4 mutex_.lock();
564 4 finished_ = true;
565 send_event();
566 mutex_.unlock();
567 }
568
569 /**
570 * \brief Symbolic constant for uninitialized hint.
571 * \details Locate functions can be accelerated by
572 * specifying a hint. This constant indicates that
573 * no hint is given.
574 */
575 static constexpr index_t NO_TETRAHEDRON = NO_INDEX;
576
577 /**
578 * \brief Symbolic value for a vertex of a
579 * tetrahedron that indicates a virtual tetrahedron.
580 * \details The three other vertices then correspond to a
581 * facet on the convex hull of the points.
582 */
583 static constexpr index_t VERTEX_AT_INFINITY = NO_INDEX;
584
585
586 /**
587 * \brief Maximum valid index for a tetrahedron.
588 * \return the maximum valid index for a tetrahedron
589 * This includes not only real tetrahedra, but also
590 * the virtual ones on the border, the conflict
591 * list and the free list.
592 */
593 index_t max_t() const {
594
1/4
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 3 return max_t_;
595 }
596
597 /**
598 * \brief Sets the maximum valid index for a tetrahedron.
599 * \details Needs to be called when starting the threads, and
600 * whenever memory allocation occured.
601 * \param[in] max_t the maximum valid index for a tetrahedron.
602 * This includes not only real tetrahedra,
603 * but also the virtual ones on the border, the conflict
604 * list and the free list.
605 */
606 void set_max_t(index_t max_t) {
607 ✗ max_t_ = max_t;
608 ✗ }
609
610 /**
611 * \brief Tests whether a given tetrahedron
612 * is a finite one.
613 * \details Infinite tetrahedra are the ones
614 * that are incident to the infinite vertex
615 * (index -1)
616 * \param[in] t the index of the tetrahedron
617 * \retval true if \p t is finite
618 * \retval false otherwise
619 */
620 63120 bool tet_is_finite(index_t t) const {
621 return
622
1/2
✓ Branch 0 taken 63120 times.
✗ Branch 1 not taken.
 63120 cell_to_v_store_[4 * t] != NO_INDEX &&
623
2/2
✓ Branch 0 taken 889 times.
✓ Branch 1 taken 62231 times.
 63120 cell_to_v_store_[4 * t + 1] != NO_INDEX &&
624
3/4
✓ Branch 0 taken 63120 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 484 times.
✓ Branch 3 taken 61747 times.
 125351 cell_to_v_store_[4 * t + 2] != NO_INDEX &&
625
2/2
✓ Branch 0 taken 798 times.
✓ Branch 1 taken 60949 times.
 61747 cell_to_v_store_[4 * t + 3] != NO_INDEX ;
626 }
627
628 /**
629 * \brief Tests whether a tetrahedron is
630 * a real one.
631 * \details Real tetrahedra are incident to
632 * four user-specified vertices (there are also
633 * virtual tetrahedra that are incident to the
634 * vertex at infinity, with index -1)
635 * \param[in] t index of the tetrahedron
636 * \retval true if tetrahedron \p t is a real one
637 * \retval false otherwise
638 */
639 bool tet_is_real(index_t t) const {
640
4/6
✓ Branch 0 taken 14247 times.
✓ Branch 1 taken 1414 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1199 times.
✓ Branch 4 taken 32558 times.
✗ Branch 5 not taken.
 49418 return !tet_is_free(t) && tet_is_finite(t);
641 }
642
643 /**
644 * \brief Tests whether a tetrahedron is
645 * a real one and has at least a vertex in the core.
646 * \details Real tetrahedra are incident to
647 * four user-specified vertices (there are also
648 * virtual tetrahedra that are incident to the
649 * vertex at infinity, with index -1)
650 * \param[in] t index of the tetrahedron
651 * \retval true if tetrahedron \p t is a real one
652 * \retval false otherwise
653 */
654 16764 bool tet_is_real_non_periodic(index_t t) const {
655
4/4
✓ Branch 0 taken 15228 times.
✓ Branch 1 taken 1536 times.
✓ Branch 2 taken 14244 times.
✓ Branch 3 taken 984 times.
 16764 return !tet_is_free(t) && tet_is_finite(t) && (
656
2/2
✓ Branch 0 taken 13781 times.
✓ Branch 1 taken 463 times.
 14244 finite_tet_vertex(t,0) < nb_vertices_non_periodic_ ||
657 12376 finite_tet_vertex(t,1) < nb_vertices_non_periodic_ ||
658 11866 finite_tet_vertex(t,2) < nb_vertices_non_periodic_ ||
659 16764 finite_tet_vertex(t,3) < nb_vertices_non_periodic_ ) ;
660 }
661
662 /**
663 * \brief Tests whether a tetrahedron is
664 * in the free list.
665 * \details Deleted tetrahedra are recycled
666 * in a free list.
667 * \param[in] t index of the tetrahedron
668 * \retval true if tetrahedron \p t is in
669 * the free list
670 * \retval false otherwise
671 */
672 bool tet_is_free(index_t t) const {
673 return tet_is_in_list(t);
674 }
675
676 /**
677 * \brief Tests whether a tetrahedron is contained
678 * by a given linked list.
679 * \details Used for debugging purposes.
680 * \param[in] t the tetrahedron
681 * \param[in] first the first element of the list
682 * \retval true if \p t is contained in the list starting
683 * at \p first
684 * \retval false otherwise
685 */
686 bool tet_is_in_list(index_t t, index_t first) const {
687 for(
688 index_t cur = first; cur != END_OF_LIST;
689 cur = tet_next(cur)
690 ) {
691 if(cur == t) {
692 return true;
693 }
694 }
695 return false;
696 }
697
698
699 /**
700 * \brief Finds in the pointset a set of four non-coplanar
701 * points and creates a tetrahedron that connects them.
702 * \details This function is used to initiate the incremental
703 * Delaunay construction, it should be called only once.
704 * \retval the index of the created tetrahedron
705 * \retval NO_TETRAHEDRON if all points were coplanar
706 */
707 2 index_t create_first_tetrahedron() {
708 index_t iv0,iv1,iv2,iv3;
709 2 if(nb_vertices() < 4) {
710 return NO_TETRAHEDRON;
711 }
712
713 iv0 = 0;
714
715 iv1 = 1;
716 while(
717
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 4 iv1 < nb_vertices_non_periodic_ &&
718 2 PCK::points_are_identical_3d(
719 non_periodic_vertex_ptr(iv0),
720 non_periodic_vertex_ptr(iv1)
721 )
722 ) {
723 ✗ ++iv1;
724 }
725 2 if(iv1 == nb_vertices()) {
726 return NO_TETRAHEDRON;
727 }
728
729 2 iv2 = iv1 + 1;
730 2 while(
731
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 4 iv2 < nb_vertices_non_periodic_ &&
732 2 PCK::points_are_colinear_3d(
733 non_periodic_vertex_ptr(iv0),
734 non_periodic_vertex_ptr(iv1),
735 non_periodic_vertex_ptr(iv2)
736 )
737 ) {
738 ✗ ++iv2;
739 }
740 2 if(iv2 == nb_vertices()) {
741 return NO_TETRAHEDRON;
742 }
743
744 2 iv3 = iv2 + 1;
745 Sign s = ZERO;
746 2 while(
747
3/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
 8 iv3 < nb_vertices_non_periodic_ &&
748 4 (s = PCK::orient_3d(
749 non_periodic_vertex_ptr(iv0),
750 non_periodic_vertex_ptr(iv1),
751 non_periodic_vertex_ptr(iv2),
752 non_periodic_vertex_ptr(iv3)
753 )) == ZERO
754 ) {
755 2 ++iv3;
756 }
757
758 2 if(iv3 == nb_vertices()) {
759 return NO_TETRAHEDRON;
760 }
761
762 geo_debug_assert(s != ZERO);
763
764
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(s == NEGATIVE) {
765 std::swap(iv2, iv3);
766 }
767
768 // Create the first tetrahedron
769 2 index_t t0 = new_tetrahedron(iv0, iv1, iv2, iv3);
770
771 // Create the first four virtual tetrahedra surrounding it
772 index_t t[4];
773
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for(index_t f = 0; f < 4; ++f) {
774 // In reverse order since it is an adjacent tetrahedron
775 index_t v1 = tet_vertex(t0, tet_facet_vertex(f,2));
776 index_t v2 = tet_vertex(t0, tet_facet_vertex(f,1));
777 index_t v3 = tet_vertex(t0, tet_facet_vertex(f,0));
778 8 t[f] = new_tetrahedron(VERTEX_AT_INFINITY, v1, v2, v3);
779 }
780
781 // Connect the virtual tetrahedra to the real one
782
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for(index_t f=0; f<4; ++f) {
783 8 set_tet_adjacent(t[f], 0, t0);
784 set_tet_adjacent(t0, f, t[f]);
785 }
786
787 // Interconnect the four virtual tetrahedra along their common
788 // faces
789
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for(index_t f = 0; f < 4; ++f) {
790 // In reverse order since it is an adjacent tetrahedron
791 index_t lv1 = tet_facet_vertex(f,2);
792 index_t lv2 = tet_facet_vertex(f,1);
793 index_t lv3 = tet_facet_vertex(f,0);
794 8 set_tet_adjacent(t[f], 1, t[lv1]);
795 8 set_tet_adjacent(t[f], 2, t[lv2]);
796 8 set_tet_adjacent(t[f], 3, t[lv3]);
797 }
798
799 2 release_tets();
800
801 return t0;
802 }
803
804
805 /**
806 * \brief Creates a star of tetrahedra filling the conflict zone.
807 * \param[in] v the index of the point to be inserted
808 * \details This function is used when the Cavity computed
809 * when traversing the conflict zone is OK, that is to say
810 * when its array sizes were not exceeded.
811 * \return the index of one the newly created tetrahedron
812 */
813 2742 index_t stellate_cavity(index_t v) {
814 index_t new_tet = NO_INDEX;
815
816
2/2
✓ Branch 0 taken 67740 times.
✓ Branch 1 taken 2742 times.
 70482 for(index_t f=0; f<cavity_.nb_facets(); ++f) {
817 index_t old_tet = cavity_.facet_tet(f);
818 index_t lf = cavity_.facet_facet(f);
819 index_t t_neigh = tet_adjacent(old_tet, lf);
820 index_t v1 = cavity_.facet_vertex(f,0);
821 index_t v2 = cavity_.facet_vertex(f,1);
822 index_t v3 = cavity_.facet_vertex(f,2);
823 67740 new_tet = new_tetrahedron(v, v1, v2, v3);
824 set_tet_adjacent(new_tet, 0, t_neigh);
825 set_tet_adjacent(
826 t_neigh, find_tet_adjacent(t_neigh,old_tet), new_tet
827 );
828 cavity_.set_facet_tet(f, new_tet);
829 }
830
831
2/2
✓ Branch 0 taken 67740 times.
✓ Branch 1 taken 2742 times.
 70482 for(index_t f=0; f<cavity_.nb_facets(); ++f) {
832 new_tet = cavity_.facet_tet(f);
833 index_t neigh1, neigh2, neigh3;
834 67740 cavity_.get_facet_neighbor_tets(f, neigh1, neigh2, neigh3);
835 67740 set_tet_adjacent(new_tet, 1, neigh1);
836 67740 set_tet_adjacent(new_tet, 2, neigh2);
837 67740 set_tet_adjacent(new_tet, 3, neigh3);
838 }
839
840 2742 return new_tet;
841 }
842
843
844 /**
845 * \brief Inserts a point in the triangulation.
846 * \param[in] v the index of the point to be inserted
847 * \param[in,out] hint the index of a tetrahedron as near as
848 * possible to \p v, or NO_TETRAHEDRON if unspecified. On
849 * exit, the index of one of the tetrahedra incident to
850 * point \p v
851 * \retval true if insertion was successful, that is, if there
852 * was no interference. This includes the situation where the
853 * point already exists (even if this does not create a new
854 * vertex)
855 * \retval false otherwise
856 */
857 4019 bool insert(index_t v, index_t& hint) {
858
859 4019 vec3 p = vertex(v);
860
861 Sign orient[4];
862 4019 index_t t = locate(v,p,hint,orient);
863
864 // locate() may fail due to tets already owned by
865 // other threads.
866
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4019 times.
 4019 if(t == NO_TETRAHEDRON) {
867 ✗ ++nb_failed_locate_;
868 geo_debug_assert(nb_acquired_tets_ == 0);
869 ✗ return false;
870 }
871
872 // At this point, t is a valid tetrahedron,
873 // and this thread acquired a lock on it.
874
875 // Test whether the point already exists in
876 // the triangulation. The point already exists
877 // if it's located on three faces of the
878 // tetrahedron returned by locate().
879 4019 int nb_zero =
880 4019 (orient[0] == ZERO) +
881 4019 (orient[1] == ZERO) +
882 4019 (orient[2] == ZERO) +
883 4019 (orient[3] == ZERO) ;
884
885
2/2
✓ Branch 0 taken 1244 times.
✓ Branch 1 taken 2775 times.
 4019 if(nb_zero >= 3) {
886 release_tet(t);
887 1244 return true;
888 }
889
890 geo_debug_assert(nb_acquired_tets_ == 1);
891
892 2775 index_t t_bndry = NO_TETRAHEDRON;
893 2775 index_t f_bndry = NO_INDEX;
894
895 2775 vec4 p_lifted = lifted_vertex(v,p);
896
897 cavity_.clear();
898
899 2775 bool ok = find_conflict_zone(v,p_lifted,t,t_bndry,f_bndry);
900
901 // When in multithreading mode, we cannot allocate memory
902 // dynamically and we use a fixed pool. If the fixed pool
903 // is full, then we exit the thread (and the missing points
904 // are inserted after, in sequential mode).
905 if(
906
3/4
✓ Branch 0 taken 948 times.
✓ Branch 1 taken 1827 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 948 times.
 3723 nb_tets_to_create_ > nb_free_ &&
907 948 Process::is_running_threads()
908 ) {
909 ✗ memory_overflow_ = true;
910 ok = false;
911 }
912
913
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2775 times.
 2775 if(!ok) {
914 // At this point, this thread did not successfully
915 // acquire all the tets in the conflict zone, so
916 // we need to rollback.
917 ✗ release_tets();
918 geo_debug_assert(nb_acquired_tets_ == 0);
919 ✗ return false;
920 }
921
922 // The conflict list can be empty if
923 // the triangulation is weighted and v is not visible
924
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2775 times.
 2775 if(tets_to_delete_.size() == 0) {
925 ✗ release_tets();
926 geo_debug_assert(nb_acquired_tets_ == 0);
927 ✗ has_empty_cells_ = true;
928 ✗ return true;
929 }
930
931
932 geo_debug_assert(
933 nb_acquired_tets_ ==
934 tets_to_delete_.size() + tets_to_release_.size()
935 );
936
937 #ifdef GEO_DEBUG
938 // Sanity check: make sure this threads owns all the tets
939 // in conflict and their neighbors.
940 for(index_t i=0; i<tets_to_delete_.size(); ++i) {
941 index_t tdel = tets_to_delete_[i];
942 geo_debug_assert(owns_tet(tdel));
943 for(index_t lf=0; lf<4; ++lf) {
944 geo_debug_assert(tet_adjacent(tdel,lf) != NO_INDEX);
945 geo_debug_assert(owns_tet(tet_adjacent(tdel,lf)));
946 }
947 }
948 #endif
949 geo_debug_assert(owns_tet(t_bndry));
950 geo_debug_assert(owns_tet(tet_adjacent(t_bndry,f_bndry)));
951 geo_debug_assert(
952 !tet_is_marked_as_conflict(tet_adjacent(t_bndry,f_bndry))
953 );
954
955 // At this point, this thread owns all the tets in conflict and
956 // their neighbors, therefore no other thread can interfere, and
957 // we can update the triangulation.
958
959 index_t new_tet = NO_INDEX;
960
2/2
✓ Branch 0 taken 2742 times.
✓ Branch 1 taken 33 times.
 2775 if(cavity_.OK()) {
961 2742 new_tet = stellate_cavity(v);
962 } else {
963 33 new_tet = stellate_conflict_zone_iterative(v,t_bndry,f_bndry);
964 }
965
966 // Recycle the tetrahedra of the conflict zone.
967
2/2
✓ Branch 0 taken 54010 times.
✓ Branch 1 taken 2775 times.
 59560 for(index_t i=0; i<tets_to_delete_.size()-1; ++i) {
968 54010 cell_next_[tets_to_delete_[i]] = tets_to_delete_[i+1];
969 }
970 2775 cell_next_[tets_to_delete_[tets_to_delete_.size()-1]] =
971 2775 first_free_;
972 2775 first_free_ = tets_to_delete_[0];
973 2775 nb_free_ += nb_tets_in_conflict();
974
975 // Reset deleted tet.
976 // This is needed because we update_v_to_cell() in
977 // a transient state.
978 // Note: update_v_to_cell() is overloaded here,
979 // with a check on nb_vertices_non_periodic_,
980 // this is why the VERTEX_OF_DELETED_TET (= -2)
981 // does not make everything crash.
982
2/2
✓ Branch 0 taken 56785 times.
✓ Branch 1 taken 2775 times.
 59560 for(index_t i=0; i<tets_to_delete_.size(); ++i) {
983 56785 index_t tdel = tets_to_delete_[i];
984 set_tet_vertex(tdel, 0, VERTEX_OF_DELETED_TET);
985 set_tet_vertex(tdel, 1, VERTEX_OF_DELETED_TET);
986 set_tet_vertex(tdel, 2, VERTEX_OF_DELETED_TET);
987 set_tet_vertex(tdel, 3, VERTEX_OF_DELETED_TET);
988 }
989
990 // Return one of the newly created tets
991 2775 hint=new_tet;
992 2775 release_tets();
993 geo_debug_assert(nb_acquired_tets_ == 0);
994 return true;
995 }
996
997 /**
998 * \brief Determines the list of tetrahedra in conflict
999 * with a given point.
1000 * \param[in] v the index of the point to be inserted
1001 * \param[in] t the index of a tetrahedron that contains
1002 * \p p, as returned by locate()
1003 * \param[out] t_bndry a tetrahedron adjacent to the
1004 * boundary of the conflict zone
1005 * \param[out] f_bndry the facet along which t_bndry is
1006 * adjacent to the boundary of the conflict zone
1007 * The other tetrahedra are linked, and can be traversed
1008 * from \p first by using tet_next() until \p last or END_OF_LIST
1009 * is reached.
1010 * The conflict zone can be empty under two circumstances:
1011 * - the vertex \p v already exists in the triangulation
1012 * - the triangulation is weighted and \p v is not visible
1013 * in either cases, both \p first and \p last contain END_OF_LIST
1014 * \retval true if all the tetrahedra of the conflict zone and their
1015 * neighbors could be acquired by this thread
1016 * \retval false otherwise
1017 */
1018 2775 bool find_conflict_zone(
1019 index_t v, const vec4& p, index_t t,
1020 index_t& t_bndry, index_t& f_bndry
1021 ) {
1022 2775 nb_tets_to_create_ = 0;
1023
1024 geo_debug_assert(t != NO_TETRAHEDRON);
1025 geo_debug_assert(owns_tet(t));
1026
1027 // Pointer to the coordinates of the point to be inserted
1028 //const double* p = vertex_ptr(v);
1029
1030 // Weighted triangulations can have dangling
1031 // vertices. Such vertices p are characterized by
1032 // the fact that p is not in conflict with the
1033 // tetrahedron returned by locate().
1034
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2775 times.
 2775 if(!tet_is_in_conflict(t,v,p)) {
1035 release_tet(t);
1036 ✗ return true;
1037 }
1038
1039 2775 mark_tet_as_conflict(t);
1040
1041 // Sanity check: the vertex to be inserted should
1042 // not correspond to one of the vertices of t.
1043 geo_debug_assert(v != tet_vertex(t,0));
1044 geo_debug_assert(v != tet_vertex(t,1));
1045 geo_debug_assert(v != tet_vertex(t,2));
1046 geo_debug_assert(v != tet_vertex(t,3));
1047
1048 // Note: points on edges and on facets are
1049 // handled by the way tet_is_in_conflict()
1050 // is implemented, that naturally inserts
1051 // the correct tetrahedra in the conflict list.
1052
1053 // Determine the conflict list by greedy propagation from t.
1054 2775 bool result = find_conflict_zone_iterative(v, p,t);
1055 2775 t_bndry = t_boundary_;
1056 2775 f_bndry = f_boundary_;
1057 2775 return result;
1058 }
1059
1060
1061 /**
1062 * \brief This function is used to implement find_conflict_zone.
1063 * \details This function detects the neighbors of \p t that are
1064 * in the conflict zone and calls itself recursively on them.
1065 * \param[in] v_in the index of the point to be inserted
1066 * \param[in] p_in the point to be inserted
1067 * \param[in] t_in index of a tetrahedron in the conflict zone
1068 * \pre The tetrahedron \p t was alredy marked as
1069 * conflict (tet_is_in_list(t))
1070 */
1071 2775 bool find_conflict_zone_iterative(
1072 index_t v_in, const vec4& p_in, index_t t_in
1073 ) {
1074 geo_debug_assert(owns_tet(t_in));
1075 5550 S_.push_back(SFrame(t_in, v_in, p_in));
1076
1077
2/2
✓ Branch 0 taken 56785 times.
✓ Branch 1 taken 2775 times.
 62335 while(S_.size() != 0) {
1078 56785 index_t t = S_.rbegin()->t;
1079 56785 index_t v = S_.rbegin()->v;
1080 56785 vec4 p = S_.rbegin()->p;
1081 S_.pop_back();
1082
1083 geo_debug_assert(owns_tet(t));
1084
1085
2/2
✓ Branch 0 taken 227140 times.
✓ Branch 1 taken 56785 times.
 283925 for(index_t lf = 0; lf < 4; ++lf) {
1086 index_t v2=v;
1087 227140 vec4 p2=p;
1088
1089 index_t t2 = tet_adjacent(t, lf);
1090
1091 // If t2 is already owned by current thread, then
1092 // its status was previously determined.
1093
2/2
✓ Branch 0 taken 104735 times.
✓ Branch 1 taken 122405 times.
 227140 if(owns_tet(t2)) {
1094
1095 geo_debug_assert(
1096 tet_is_marked_as_conflict(t2) ==
1097 tet_is_in_conflict(t2,v2,p2)
1098 );
1099
1100 // If t2 is not in conflict list, then t has a facet
1101 // on the border of the conflict zone, and there is
1102 // a tet to create.
1103
2/2
✓ Branch 0 taken 4693 times.
✓ Branch 1 taken 100042 times.
 104735 if(!tet_is_marked_as_conflict(t2)) {
1104 4693 ++nb_tets_to_create_;
1105 4693 cavity_.new_facet(
1106 t, lf,
1107 tet_vertex(t, tet_facet_vertex(lf,0)),
1108 tet_vertex(t, tet_facet_vertex(lf,1)),
1109 tet_vertex(t, tet_facet_vertex(lf,2))
1110 );
1111 }
1112 158745 continue;
1113 }
1114
1115 if(!acquire_tet(t2)) {
1116 ✗ S_.resize(0);
1117 ✗ return false;
1118 }
1119
1120 geo_debug_assert(owns_tet(t2));
1121
1122
3/4
✓ Branch 1 taken 122405 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 68395 times.
✓ Branch 4 taken 54010 times.
 122405 if(!tet_is_in_conflict(t2,v2,p2)) {
1123 68395 mark_tet_as_neighbor(t2);
1124 // If t2 is not in conflict list, then t has a facet
1125 // on the border of the conflict zone, and there is
1126 // a tet to create.
1127 68395 ++nb_tets_to_create_;
1128 } else {
1129
1/2
✓ Branch 1 taken 54010 times.
✗ Branch 2 not taken.
 54010 mark_tet_as_conflict(t2);
1130 geo_debug_assert(owns_tet(t2));
1131 54010 S_.push_back(SFrame(t2,v2,p2));
1132 54010 continue;
1133 }
1134
1135 // At this point, t is in conflict
1136 // and t2 is not in conflict.
1137 // We keep a reference to a tet on the boundary
1138 68395 t_boundary_ = t;
1139 68395 f_boundary_ = lf;
1140 68395 cavity_.new_facet(
1141 t, lf,
1142 tet_vertex(t, tet_facet_vertex(lf,0)),
1143 tet_vertex(t, tet_facet_vertex(lf,1)),
1144 tet_vertex(t, tet_facet_vertex(lf,2))
1145 );
1146 geo_debug_assert(tet_adjacent(t,lf) == t2);
1147 geo_debug_assert(owns_tet(t));
1148 geo_debug_assert(owns_tet(t2));
1149 }
1150 }
1151 return true;
1152 }
1153
1154 /**
1155 * \brief Gets a pointer to a vertex by its global index.
1156 * \param[in] i global index of the vertex
1157 * \return a pointer to vertex \p i
1158 */
1159 const double* non_periodic_vertex_ptr(index_t i) const {
1160 geo_debug_assert(i < nb_vertices_non_periodic_);
1161 8 return vertices_ + i*3;
1162 }
1163
1164 /**
1165 * \brief Gets the weight of a vertex by its global index.
1166 * \param[in] i global index of the vertex
1167 * \return the weight associated with vertex i
1168 */
1169 double non_periodic_weight(index_t i) const {
1170 geo_debug_assert(i < nb_vertices_non_periodic_);
1171
10/20
✗ Branch 0 not taken.
✓ Branch 1 taken 5281 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5281 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 5281 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5281 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 5281 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 691 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 691 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 691 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 691 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2775 times.
 8747 return (weights_ == nullptr) ? 0.0 : weights_[i];
1172 }
1173
1174
1175 /**
1176 * \brief gets a vertex.
1177 * \param[in] v the index of the vertex, in 0..nb_vertices_-1
1178 * \param[out] result a pointer to the 3d coordinates of the vertex.
1179 * \details In periodic mode, does vertex translation (v is a
1180 * virtual vertex index).
1181 */
1182 181603 void get_vertex(index_t v, double* result) const {
1183
2/2
✓ Branch 0 taken 15905 times.
✓ Branch 1 taken 165698 times.
 181603 if(!periodic_) {
1184 15905 result[0] = vertices_[3*v];
1185 15905 result[1] = vertices_[3*v+1];
1186 15905 result[2] = vertices_[3*v+2];
1187 15905 return;
1188 }
1189 index_t instance = periodic_vertex_instance(v);
1190 v = periodic_vertex_real(v);
1191 165698 result[0] = vertices_[3*v];
1192 165698 result[1] = vertices_[3*v+1];
1193 165698 result[2] = vertices_[3*v+2];
1194 165698 result[0] += double(translation[instance][0]) * period_.x;
1195 165698 result[1] += double(translation[instance][1]) * period_.y;
1196 165698 result[2] += double(translation[instance][2]) * period_.z;
1197 }
1198
1199 /**
1200 * \brief gets a vertex.
1201 * \param[in] v the index of the vertex, in 0..nb_vertices_-1
1202 * \details In periodic mode, does vertex translation (v is a
1203 * virtual vertex index).
1204 * \return the 3d vertex.
1205 */
1206 vec3 vertex(index_t v) const {
1207 vec3 result;
1208 36577 get_vertex(v, result.data());
1209 return result;
1210 }
1211
1212 /**
1213 * \brief gets a lifted vertex.
1214 * \param[in] v the index of the vertex, in 0..nb_vertices_-1
1215 * \param[out] result the 4d coordinates of the vertex,
1216 * shifted on the paraboloid and shifted by the weight.
1217 * \details In periodic mode, does vertex translation (v is a
1218 * virtual vertex index).
1219 */
1220 13040 void get_lifted_vertex(index_t v, double* result) const {
1221 index_t instance = 0;
1222
1/2
✓ Branch 0 taken 13040 times.
✗ Branch 1 not taken.
 13040 if(periodic_) {
1223 instance = periodic_vertex_instance(v);
1224 v = periodic_vertex_real(v);
1225 }
1226 13040 result[0] = vertices_[3*v];
1227 13040 result[1] = vertices_[3*v+1];
1228
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13040 times.
 13040 result[2] = vertices_[3*v+2];
1229 13040 result[3] = -non_periodic_weight(v);
1230
1/2
✓ Branch 0 taken 13040 times.
✗ Branch 1 not taken.
 13040 if(periodic_) {
1231 13040 result[0] += double(translation[instance][0]) * period_.x;
1232 13040 result[1] += double(translation[instance][1]) * period_.y;
1233 13040 result[2] += double(translation[instance][2]) * period_.z;
1234 }
1235 13040 result[3] +=
1236 13040 geo_sqr(result[0]) + geo_sqr(result[1]) + geo_sqr(result[2]);
1237 13040 }
1238
1239 /**
1240 * \brief gets a lifted vertex.
1241 * \param[in] v the index of the vertex, in 0..nb_vertices_-1
1242 * \details In periodic mode, does vertex translation (v is a
1243 * virtual vertex index).
1244 * \return the 4d vertex, lifted on the paraboloid, and shifted by
1245 * the weight.
1246 */
1247 vec4 lifted_vertex(index_t v) const {
1248 vec4 result;
1249 ✗ get_lifted_vertex(v,result.data());
1250 return result;
1251 }
1252
1253 /**
1254 * \brief gets a lifted vertex.
1255 * \param[in] v the index of the vertex, in 0..nb_vertices_-1
1256 * \param[in] p the geometric location of vertex v
1257 * \details In periodic mode, does vertex translation (v is a
1258 * virtual vertex index).
1259 * \return the 4d vertex, lifted on the paraboloid, and shifted by
1260 * the weight.
1261 */
1262 vec4 lifted_vertex(index_t v, const vec3& p) {
1263 return vec4(
1264 p.x, p.y, p.z,
1265
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2775 times.
 2775 geo_sqr(p.x) + geo_sqr(p.y) + geo_sqr(p.z)
1266 2775 - non_periodic_weight(periodic_vertex_real(v))
1267 );
1268 }
1269
1270 /**
1271 * \brief Orientation predicate with indices.
1272 * \param i , j , k , l the four indices of the four vertices.
1273 * \see PCK::orient_3d()
1274 */
1275 11838 Sign orient_3d(index_t i, index_t j, index_t k, index_t l) const {
1276 // No need to reorder since there is no SOS.
1277 double V[4][3];
1278 11838 get_vertex(i,V[0]);
1279 11838 get_vertex(j,V[1]);
1280 11838 get_vertex(k,V[2]);
1281 11838 get_vertex(l,V[3]);
1282 11838 return PCK::orient_3d(V[0], V[1], V[2], V[3]);
1283 }
1284
1285 /**
1286 * \brief Orientation predicate with indices.
1287 * \param i , j , k , l the four indices of the four vertices.
1288 * \see PCK::in_circle_3dlifted_SOS()
1289 */
1290 3951 Sign in_circle_3dlifted_SOS(
1291 index_t i, index_t j, index_t k, index_t l
1292 ) const {
1293
1294 // In non-periodic mode, directly access vertices.
1295
2/2
✓ Branch 0 taken 691 times.
✓ Branch 1 taken 3260 times.
 3951 if(!periodic_) {
1296 const double* pi = non_periodic_vertex_ptr(i);
1297 const double* pj = non_periodic_vertex_ptr(j);
1298 const double* pk = non_periodic_vertex_ptr(k);
1299 const double* pl = non_periodic_vertex_ptr(l);
1300
1301
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 691 times.
 691 double hi = geo_sqr(pi[0]) + geo_sqr(pi[1]) + geo_sqr(pi[2])
1302 691 - non_periodic_weight(i);
1303
1304
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 691 times.
 691 double hj = geo_sqr(pj[0]) + geo_sqr(pj[1]) + geo_sqr(pj[2])
1305 691 - non_periodic_weight(j);
1306
1307
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 691 times.
 691 double hk = geo_sqr(pk[0]) + geo_sqr(pk[1]) + geo_sqr(pk[2])
1308 691 - non_periodic_weight(k);
1309
1310
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 691 times.
 691 double hl = geo_sqr(pl[0]) + geo_sqr(pl[1]) + geo_sqr(pl[2])
1311 691 - non_periodic_weight(l);
1312
1313 691 return PCK::in_circle_3dlifted_SOS(
1314 pi, pj, pk, pl,
1315 hi, hj, hk, hl
1316 );
1317 }
1318
1319 // Note: in periodic mode, SOS mode is lexicographic.
1320 double V[4][4];
1321 3260 get_lifted_vertex(i,V[0]);
1322 3260 get_lifted_vertex(j,V[1]);
1323 3260 get_lifted_vertex(k,V[2]);
1324 3260 get_lifted_vertex(l,V[3]);
1325 3260 return PCK::in_circle_3dlifted_SOS(
1326 V[0], V[1], V[2], V[3],
1327 V[0][3], V[1][3], V[2][3], V[3][3]
1328 );
1329 }
1330
1331 /**
1332 * \brief 4D Orientation predicate with indices.
1333 * \param i , j , k , l , m the five indices of the four vertices.
1334 * \see PCK::orient_3dlifted_SOS()
1335 */
1336 113342 Sign orient_3dlifted_SOS(
1337 index_t i, index_t j, index_t k, index_t l, index_t m
1338 ) const {
1339
1340 // In non-periodic mode, directly access vertices.
1341
2/2
✓ Branch 0 taken 5281 times.
✓ Branch 1 taken 108061 times.
 113342 if(!periodic_) {
1342 const double* pi = non_periodic_vertex_ptr(i);
1343 const double* pj = non_periodic_vertex_ptr(j);
1344 const double* pk = non_periodic_vertex_ptr(k);
1345 const double* pl = non_periodic_vertex_ptr(l);
1346 const double* pm = non_periodic_vertex_ptr(m);
1347
1348
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5281 times.
 5281 double hi = geo_sqr(pi[0]) + geo_sqr(pi[1]) + geo_sqr(pi[2])
1349 5281 - non_periodic_weight(i);
1350
1351
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5281 times.
 5281 double hj = geo_sqr(pj[0]) + geo_sqr(pj[1]) + geo_sqr(pj[2])
1352 5281 - non_periodic_weight(j);
1353
1354
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5281 times.
 5281 double hk = geo_sqr(pk[0]) + geo_sqr(pk[1]) + geo_sqr(pk[2])
1355 5281 - non_periodic_weight(k);
1356
1357
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5281 times.
 5281 double hl = geo_sqr(pl[0]) + geo_sqr(pl[1]) + geo_sqr(pl[2])
1358 5281 - non_periodic_weight(l);
1359
1360
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5281 times.
 5281 double hm = geo_sqr(pm[0]) + geo_sqr(pm[1]) + geo_sqr(pm[2])
1361 5281 - non_periodic_weight(m);
1362
1363 5281 return PCK::orient_3dlifted_SOS(
1364 pi, pj, pk, pl, pm,
1365 hi, hj, hk, hl, hm
1366 );
1367 }
1368
1369 // Periodic mode.
1370 // Note: in periodic mode, SOS mode is lexicographic.
1371 double V[5][4];
1372
1373 /*
1374 The code below is an inlined version of:
1375 (gains a little bit)
1376 get_lifted_vertex(i,V[0]);
1377 get_lifted_vertex(j,V[1]);
1378 get_lifted_vertex(k,V[2]);
1379 get_lifted_vertex(l,V[3]);
1380 get_lifted_vertex(m,V[4]);
1381 */
1382
1383 index_t ii = periodic_vertex_instance(i);
1384 index_t ij = periodic_vertex_instance(j);
1385 index_t ik = periodic_vertex_instance(k);
1386 index_t il = periodic_vertex_instance(l);
1387 index_t im = periodic_vertex_instance(m);
1388
1389 i = periodic_vertex_real(i);
1390 j = periodic_vertex_real(j);
1391 k = periodic_vertex_real(k);
1392 l = periodic_vertex_real(l);
1393 m = periodic_vertex_real(m);
1394
1395 108061 V[0][0] = vertices_[3*i ] + double(translation[ii][0]) * period_.x;
1396 108061 V[0][1] = vertices_[3*i+1] + double(translation[ii][1]) * period_.y;
1397 108061 V[0][2] = vertices_[3*i+2] + double(translation[ii][2]) * period_.z;
1398 108061 V[1][0] = vertices_[3*j ] + double(translation[ij][0]) * period_.x;
1399 108061 V[1][1] = vertices_[3*j+1] + double(translation[ij][1]) * period_.y;
1400 108061 V[1][2] = vertices_[3*j+2] + double(translation[ij][2]) * period_.z;
1401 108061 V[2][0] = vertices_[3*k ] + double(translation[ik][0]) * period_.x;
1402 108061 V[2][1] = vertices_[3*k+1] + double(translation[ik][1]) * period_.y;
1403 108061 V[2][2] = vertices_[3*k+2] + double(translation[ik][2]) * period_.z;
1404 108061 V[3][0] = vertices_[3*l ] + double(translation[il][0]) * period_.x;
1405 108061 V[3][1] = vertices_[3*l+1] + double(translation[il][1]) * period_.y;
1406 108061 V[3][2] = vertices_[3*l+2] + double(translation[il][2]) * period_.z;
1407 108061 V[4][0] = vertices_[3*m ] + double(translation[im][0]) * period_.x;
1408 108061 V[4][1] = vertices_[3*m+1] + double(translation[im][1]) * period_.y;
1409
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 108061 times.
 108061 V[4][2] = vertices_[3*m+2] + double(translation[im][2]) * period_.z;
1410
1411 // Beware the parentheses, they are necessary to ensure that computations
1412 // | give always the same result.
1413 // |________________________________________________________________________
1414 // | |
1415 // v v
1416
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 108061 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 108061 times.
 108061 V[0][3] = -non_periodic_weight(i) + (geo_sqr(V[0][0]) + geo_sqr(V[0][1]) + geo_sqr(V[0][2]));
1417
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 108061 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 108061 times.
 108061 V[1][3] = -non_periodic_weight(j) + (geo_sqr(V[1][0]) + geo_sqr(V[1][1]) + geo_sqr(V[1][2]));
1418
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 108061 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 108061 times.
 108061 V[2][3] = -non_periodic_weight(k) + (geo_sqr(V[2][0]) + geo_sqr(V[2][1]) + geo_sqr(V[2][2]));
1419
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 108061 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 108061 times.
 108061 V[3][3] = -non_periodic_weight(l) + (geo_sqr(V[3][0]) + geo_sqr(V[3][1]) + geo_sqr(V[3][2]));
1420
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 108061 times.
 108061 V[4][3] = -non_periodic_weight(m) + (geo_sqr(V[4][0]) + geo_sqr(V[4][1]) + geo_sqr(V[4][2]));
1421
1422 108061 return PCK::orient_3dlifted_SOS(
1423 V[0], V[1], V[2], V[3], V[4],
1424 V[0][3], V[1][3], V[2][3], V[3][3], V[4][3]
1425 );
1426 }
1427
1428 /**
1429 * \brief Tests whether a given tetrahedron is in conflict with
1430 * a given 3d point.
1431 * \details A real tetrahedron is in conflict with a point whenever
1432 * the point is contained by its circumscribed sphere, and a
1433 * virtual tetrahedron is in conflict with a point whenever the
1434 * tetrahedron formed by its real face and with the point has
1435 * positive orientation.
1436 * \param[in] t the index of the tetrahedron
1437 * \param[in] v the index of the point
1438 * \param[in] p a pointer to the coordinates of the point
1439 * \retval true if point \p p is in conflict with tetrahedron \p t
1440 * \retval false otherwise
1441 */
1442 125180 bool tet_is_in_conflict(index_t t, index_t v, const vec4& p) const {
1443
1444 geo_argused(p);
1445
1446 // Lookup tetrahedron vertices
1447
1448 index_t iv[4];
1449
2/2
✓ Branch 0 taken 500720 times.
✓ Branch 1 taken 125180 times.
 625900 for(index_t i=0; i<4; ++i) {
1450 500720 iv[i] = tet_vertex(t,i);
1451 }
1452
1453
1454 // Check for virtual tetrahedra (then in_sphere()
1455 // is replaced with orient3d())
1456
2/2
✓ Branch 0 taken 488421 times.
✓ Branch 1 taken 113342 times.
 601763 for(index_t lf = 0; lf < 4; ++lf) {
1457
1458
2/2
✓ Branch 0 taken 11838 times.
✓ Branch 1 taken 476583 times.
 488421 if(iv[lf] == NO_INDEX) {
1459
1460 // Facet of a virtual tetrahedron opposite to
1461 // infinite vertex corresponds to
1462 // the triangle on the convex hull of the points.
1463 // Orientation is obtained by replacing vertex lf
1464 // with p.
1465 11838 iv[lf] = v;
1466
1467 // no SOS, we can directly use the PCK predicate
1468 11838 Sign sign = orient_3d(iv[0], iv[1], iv[2], iv[3]);
1469
1470
2/2
✓ Branch 0 taken 10070 times.
✓ Branch 1 taken 1768 times.
 11838 if(sign > 0) {
1471 return true;
1472 }
1473
1474
2/2
✓ Branch 0 taken 6825 times.
✓ Branch 1 taken 3245 times.
 10070 if(sign < 0) {
1475 return false;
1476 }
1477
1478 // If sign is zero, we check the real tetrahedron
1479 // adjacent to the facet on the convex hull.
1480 geo_debug_assert(tet_adjacent(t, lf) != NO_INDEX);
1481 index_t t2 = tet_adjacent(t, lf);
1482 geo_debug_assert(!tet_is_virtual(t2));
1483
1484 // If t2 was already visited by this thread, then
1485 // it is in conflict if it is already marked.
1486
2/2
✓ Branch 0 taken 2874 times.
✓ Branch 1 taken 3951 times.
 6825 if(owns_tet(t2)) {
1487 2874 return tet_is_marked_as_conflict(t2);
1488 }
1489
1490 // If t2 was not already visited, then we need to
1491 // switch to the in_circum_circle_3d() predicate.
1492
1493 3951 index_t iq0 = iv[(lf+1)%4];
1494 3951 index_t iq1 = iv[(lf+2)%4];
1495 3951 index_t iq2 = iv[(lf+3)%4];
1496
1497 3951 return (in_circle_3dlifted_SOS(iq0, iq1, iq2, v) > 0);
1498 }
1499 }
1500
1501 // If the tetrahedron is a finite one, it is in conflict
1502 // if its circumscribed sphere contains the point (this is
1503 // the standard case).
1504
1505 113342 return (orient_3dlifted_SOS(iv[0], iv[1], iv[2], iv[3], v) > 0);
1506 }
1507
1508
1509 /**
1510 * \brief Finds the tetrahedron that contains a point.
1511 * \details The tetrahedron is acquired by this thread. If the
1512 * tetrahedron could not be acquired, then NO_TETRAHEDRON is returned.
1513 * If the point is on a face, edge or vertex,
1514 * the function returns one of the tetrahedra incident
1515 * to that face, edge or vertex.
1516 * \param[in] v the index of the vertex to locate
1517 * \param[in] p the coordinates of the point that correspond to v
1518 * \param[out] orient a pointer to an array of four Sign%s
1519 * or nullptr. If non-nullptr, returns the orientation with respect
1520 * to the four facets of the tetrahedron that contains \p p.
1521 * \retval the index of a tetrahedron that contains \p p.
1522 * If the point is outside the convex hull of
1523 * the inserted so-far points, then the returned tetrahedron
1524 * is a virtual one (first vertex is the "vertex at infinity"
1525 * of index -1)
1526 * \retval NO_TETRAHEDRON if the tetrahedron could not be
1527 * acquired by this thread, or if the virtual tetrahedra
1528 * were previously removed
1529 */
1530 4019 index_t locate(
1531 index_t& v, vec3& p, index_t hint = NO_TETRAHEDRON,
1532 Sign* orient = nullptr
1533 ) {
1534 geo_argused(v);
1535 4019 nb_traversed_tets_ = 0;
1536
1537 // If no hint specified, find a tetrahedron randomly
1538
1539
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4015 times.
 4019 if(hint != NO_TETRAHEDRON) {
1540
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4015 times.
 4015 if(tet_is_free(hint)) {
1541 hint = NO_TETRAHEDRON;
1542 } else {
1543
1/2
✓ Branch 0 taken 4015 times.
✗ Branch 1 not taken.
 4015 if( !owns_tet(hint) && !acquire_tet(hint) ) {
1544 hint = NO_TETRAHEDRON;
1545 }
1546
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4015 times.
 4015 if((hint != NO_TETRAHEDRON) && tet_is_free(hint)) {
1547 release_tet(hint);
1548 hint = NO_TETRAHEDRON;
1549 }
1550 }
1551 }
1552
1553 do {
1554
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 4020 times.
 4025 while(hint == NO_TETRAHEDRON) {
1555
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 hint = master_->thread(0)->pick_random_tet();
1556 // we could also pick from a random thread,
1557 // but at initialization only thread0 has tets,
1558 // so let us keep thread0 for now
1559 }
1560 if(
1561
4/4
✓ Branch 0 taken 4019 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 4015 times.
 4020 tet_is_free(hint) || (!owns_tet(hint) && !acquire_tet(hint))
1562 ) {
1563
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if(owns_tet(hint)) {
1564 release_tet(hint);
1565 }
1566 hint = NO_TETRAHEDRON;
1567 } else {
1568
2/2
✓ Branch 0 taken 15588 times.
✓ Branch 1 taken 3543 times.
 19131 for(index_t f=0; f<4; ++f) {
1569
2/2
✓ Branch 0 taken 476 times.
✓ Branch 1 taken 15112 times.
 15588 if(tet_vertex(hint,f) == VERTEX_AT_INFINITY) {
1570 index_t new_hint = tet_adjacent(hint,f);
1571 if(
1572 476 tet_is_free(new_hint) ||
1573 !acquire_tet(new_hint)
1574 ) {
1575 new_hint = NO_TETRAHEDRON;
1576 }
1577 release_tet(hint);
1578 hint = new_hint;
1579 476 break;
1580 }
1581 }
1582 }
1583
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4019 times.
 4020 } while(hint == NO_TETRAHEDRON) ;
1584
1585 index_t t = hint;
1586 index_t t_pred = NO_TETRAHEDRON;
1587 Sign orient_local[4];
1588
1/2
✓ Branch 0 taken 4019 times.
✗ Branch 1 not taken.
 4019 if(orient == nullptr) {
1589 orient = orient_local;
1590 }
1591
1592
1593 4019 still_walking:
1594 {
1595
2/2
✓ Branch 0 taken 28539 times.
✓ Branch 1 taken 4019 times.
 32558 if(t_pred != NO_TETRAHEDRON) {
1596 release_tet(t_pred);
1597 }
1598
1599
1/2
✓ Branch 0 taken 32558 times.
✗ Branch 1 not taken.
 32558 if(tet_is_free(t)) {
1600 398 return NO_TETRAHEDRON;
1601 }
1602
1603
2/2
✓ Branch 0 taken 28539 times.
✓ Branch 1 taken 4019 times.
 32558 if(!owns_tet(t) && !acquire_tet(t)) {
1604 return NO_TETRAHEDRON;
1605 }
1606
1607
1608
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32558 times.
 32558 if(!tet_is_real(t)) {
1609 release_tet(t);
1610 ✗ return NO_TETRAHEDRON;
1611 }
1612
1613
2/2
✓ Branch 0 taken 130232 times.
✓ Branch 1 taken 32558 times.
 162790 vec3 pv[4];
1614 32558 pv[0] = vertex(finite_tet_vertex(t,0));
1615 32558 pv[1] = vertex(finite_tet_vertex(t,1));
1616 32558 pv[2] = vertex(finite_tet_vertex(t,2));
1617 32558 pv[3] = vertex(finite_tet_vertex(t,3));
1618
1619 // Start from a random facet
1620 index_t f0 = thread_safe_random_4_();
1621
2/2
✓ Branch 0 taken 78761 times.
✓ Branch 1 taken 3621 times.
 82382 for(index_t df = 0; df < 4; ++df) {
1622 78761 index_t f = (f0 + df) % 4;
1623
1624 index_t t_next = tet_adjacent(t,f);
1625
1626 // If the opposite tet is -1, then it means that
1627 // we are trying to locate() (e.g. called from
1628 // nearest_vertex) within a tetrahedralization
1629 // from which the infinite tets were removed.
1630
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 78761 times.
 78761 if(t_next == NO_INDEX) {
1631 release_tet(t);
1632 ✗ return NO_TETRAHEDRON;
1633 }
1634
1635 // If the candidate next tetrahedron is the
1636 // one we came from, then we know already that
1637 // the orientation is positive, thus we examine
1638 // the next candidate (or exit the loop if they
1639 // are exhausted).
1640
2/2
✓ Branch 0 taken 14322 times.
✓ Branch 1 taken 64439 times.
 78761 if(t_next == t_pred) {
1641 14322 orient[f] = POSITIVE ;
1642 14322 continue ;
1643 }
1644
1645 // To test the orientation of p w.r.t. the facet f of
1646 // t, we replace vertex number f with p in t (same
1647 // convention as in CGAL).
1648 // This is equivalent to tet_facet_point_orient3d(t,f,p)
1649 // (but less costly, saves a couple of lookups)
1650 64439 vec3 pv_bkp = pv[f];
1651 64439 pv[f] = p;
1652 64439 orient[f] = PCK::orient_3d(
1653 pv[0].data(), pv[1].data(), pv[2].data(), pv[3].data()
1654 );
1655
1656 // If the orientation is not negative, then we cannot
1657 // walk towards t_next, and examine the next candidate
1658 // (or exit the loop if they are exhausted).
1659
2/2
✓ Branch 0 taken 35502 times.
✓ Branch 1 taken 28937 times.
 64439 if(orient[f] != NEGATIVE) {
1660 35502 pv[f] = pv_bkp;
1661 35502 continue;
1662 }
1663
1664 // If the opposite tet is a virtual tet, then
1665 // the point has a positive orientation relative
1666 // to the facet on the border of the convex hull,
1667 // thus t_next is a tet in conflict and we are
1668 // done.
1669
2/2
✓ Branch 0 taken 398 times.
✓ Branch 1 taken 28539 times.
 28937 if(tet_is_virtual(t_next)) {
1670 release_tet(t);
1671 if(!acquire_tet(t_next)) {
1672 return NO_TETRAHEDRON;
1673 }
1674
2/2
✓ Branch 0 taken 1592 times.
✓ Branch 1 taken 398 times.
 1990 for(index_t lf = 0; lf < 4; ++lf) {
1675 1592 orient[lf] = POSITIVE;
1676 }
1677 return t_next;
1678 }
1679
1680 28539 ++nb_traversed_tets_;
1681
1682 // If we reach this point, then t_next is a valid
1683 // successor, thus we are still walking.
1684 t_pred = t;
1685 t = t_next;
1686 28539 goto still_walking;
1687 }
1688 }
1689
1690 // If we reach this point, we did not find a valid successor
1691 // for walking (a face for which p has negative orientation),
1692 // thus we reached the tet for which p has all positive
1693 // face orientations (i.e. the tet that contains p).
1694
1695 #ifdef GEO_DEBUG
1696 geo_debug_assert(tet_is_real(t));
1697
1698 vec3 pv[4];
1699 Sign signs[4];
1700 pv[0] = vertex(finite_tet_vertex(t,0));
1701 pv[1] = vertex(finite_tet_vertex(t,1));
1702 pv[2] = vertex(finite_tet_vertex(t,2));
1703 pv[3] = vertex(finite_tet_vertex(t,3));
1704 for(index_t f=0; f<4; ++f) {
1705 vec3 pv_bkp = pv[f];
1706 pv[f] = vec3(p.x, p.y, p.z);
1707 signs[f] = PCK::orient_3d(
1708 pv[0].data(), pv[1].data(), pv[2].data(), pv[3].data()
1709 );
1710 geo_debug_assert(signs[f] >= 0);
1711 pv[f] = pv_bkp;
1712 }
1713 #endif
1714
1715 3621 return t;
1716 }
1717
1718
1719 protected:
1720
1721 /**
1722 * \brief Tests whether a tetrahedron was marked as conflict.
1723 * \pre owns_tet(t)
1724 * \param[in] t the index of the tetrahedron to be tested
1725 * \retval true if \p t was marked as conflict
1726 * \retval false otherwise
1727 */
1728 bool tet_is_marked_as_conflict(index_t t) const {
1729 geo_debug_assert(owns_tet(t));
1730
4/4
✓ Branch 0 taken 10384 times.
✓ Branch 1 taken 8022 times.
✓ Branch 2 taken 4693 times.
✓ Branch 3 taken 100042 times.
 126015 return cell_status_.cell_is_marked_as_conflict(t);
1731 }
1732
1733
1734 /**
1735 * \brief Gets the number of tetrahedra in conflict.
1736 * \return the number of tetrahedra in conflict,
1737 * specified by mark_tet_as_conflict()
1738 */
1739 index_t nb_tets_in_conflict() const {
1740 return tets_to_delete_.size();
1741 }
1742
1743 /**
1744 * \brief Marks a tetrahedron as conflict.
1745 * \details The index of the tetrahedron is also
1746 * stored it in the list of conflict tetrahedra.
1747 * \param[in] t index of the tetrahedron to mark
1748 * \pre owns_tet(t)
1749 */
1750 56785 void mark_tet_as_conflict(index_t t) {
1751 geo_debug_assert(owns_tet(t));
1752
2/2
✓ Branch 0 taken 56767 times.
✓ Branch 1 taken 18 times.
 56785 tets_to_delete_.push_back(t);
1753 56785 cell_status_.mark_cell_as_conflict(t);
1754 geo_debug_assert(owns_tet(t));
1755 geo_debug_assert(tet_is_marked_as_conflict(t));
1756 56785 }
1757
1758 /**
1759 * \brief Marks a tetrahedron as neighbor of the conflict zone.
1760 * \details The index of the tetrahedron is also
1761 * stored it in the list of tetrahedra to release.
1762 * \param[in] t index of the tetrahedron to mark
1763 * \pre owns_tet(t)
1764 */
1765 void mark_tet_as_neighbor(index_t t) {
1766 // Note: nothing to change in cell_status_[t]
1767 // since MSB=0 means neigbhor tet.
1768
1/2
✓ Branch 0 taken 68395 times.
✗ Branch 1 not taken.
 68395 tets_to_release_.push_back(t);
1769 }
1770
1771 /**
1772 * \brief Acquires a lock on a tetrahedron and keep
1773 * it in the list of acquired tetrahedra.
1774 * \param[in] t index of the tetrahedron to acquire
1775 */
1776 73098 void acquire_and_mark_tet_as_created(index_t t) {
1777 // The tet was created in this thread's tet pool,
1778 // therefore there is no need to use sync
1779 // primitives to acquire a lock on it.
1780 geo_debug_assert(cell_status_.cell_thread(t) == NO_THREAD);
1781
2/2
✓ Branch 0 taken 73078 times.
✓ Branch 1 taken 20 times.
 73098 cell_status_.set_cell_status(
1782 t, CellStatusArray::thread_index_t(id())
1783 );
1784
1785 #ifdef GEO_DEBUG
1786 ++nb_acquired_tets_;
1787 #endif
1788
2/2
✓ Branch 0 taken 73078 times.
✓ Branch 1 taken 20 times.
 73098 tets_to_release_.push_back(t);
1789 73098 }
1790
1791
1792 /**
1793 * \brief Releases all the tetrahedron locks that were
1794 * acquired using mark_tet_as_neighbor(),
1795 * acquire_and_mark_tet_as_created() and mark_as_conflict().
1796 */
1797 2777 void release_tets() {
1798
2/2
✓ Branch 0 taken 141493 times.
✓ Branch 1 taken 2777 times.
 147047 for(index_t i=0; i<tets_to_release_.size(); ++i) {
1799 141493 release_tet(tets_to_release_[i]);
1800 }
1801 2777 tets_to_release_.resize(0);
1802
2/2
✓ Branch 0 taken 56785 times.
✓ Branch 1 taken 2777 times.
 62339 for(index_t i=0; i<tets_to_delete_.size(); ++i) {
1803 56785 release_tet(tets_to_delete_[i]);
1804 }
1805 2777 tets_to_delete_.resize(0);
1806 2777 }
1807
1808 /**
1809 * \brief Atomically acquires a lock on a tetrahedron.
1810 * \details When the lock could not be acquired, interfering_thread_
1811 * contains the id of the thread that owns the lock.
1812 * \param[in] t the index of the tetrahedron to acquire
1813 * \retval true if the lock was successfully acquired
1814 * \retval false otherwise
1815 */
1816 bool acquire_tet(index_t t) {
1817 geo_debug_assert(t < max_t());
1818 geo_debug_assert(!owns_tet(t));
1819
1820 155837 interfering_thread_ = cell_status_.acquire_cell(
1821 t, CellStatusArray::thread_index_t(id())
1822 );
1823
1824
6/12
✗ Branch 0 not taken.
✓ Branch 1 taken 4015 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 476 times.
✓ Branch 6 taken 28539 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 398 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 122405 times.
 155837 if(interfering_thread_ == NO_THREAD) {
1825 geo_debug_assert(t == first_free_ || !tet_is_in_list(t));
1826 #ifdef GEO_DEBUG
1827 ++nb_acquired_tets_;
1828 #endif
1829 return true;
1830 }
1831 return false;
1832 }
1833
1834 /**
1835 * \brief Releases a lock on a tetrahedron, making it
1836 * available to the other threads.
1837 */
1838 void release_tet(index_t t) {
1839 geo_debug_assert(t < max_t());
1840 geo_debug_assert(owns_tet(t));
1841 #ifdef GEO_DEBUG
1842 --nb_acquired_tets_;
1843 #endif
1844 228935 cell_status_.release_cell(t);
1845 28539 }
1846
1847
1848 /**
1849 * \brief Tests whether this thread owns a tetrahedron.
1850 * \param[in] t index of the tetrahedron
1851 * \retval true if this thread owns t
1852 * \retval false otherwise
1853 */
1854 bool owns_tet(index_t t) const {
1855 geo_debug_assert(t < max_t());
1856 return (
1857
10/12
✓ Branch 0 taken 4015 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 4015 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 28539 times.
✓ Branch 7 taken 4019 times.
✓ Branch 8 taken 2874 times.
✓ Branch 9 taken 3951 times.
✓ Branch 10 taken 104735 times.
✓ Branch 11 taken 122405 times.
 274558 cell_status_.cell_thread(t) ==
1858 CellStatusArray::thread_index_t(id())
1859 );
1860 }
1861
1862 /**
1863 * \brief Tests whether a tetrahedron is
1864 * a virtual one.
1865 * \details Virtual tetrahedra are tetrahedra
1866 * incident to the vertex at infinity.
1867 * \param[in] t index of the tetrahedron
1868 * \retval true if tetrahedron \p t is virtual
1869 * \retval false otherwise
1870 */
1871 28937 bool tet_is_virtual(index_t t) const {
1872 return
1873
3/4
✓ Branch 0 taken 28937 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 28935 times.
✓ Branch 3 taken 2 times.
 28937 !tet_is_free(t) && (
1874
2/2
✓ Branch 0 taken 28935 times.
✓ Branch 1 taken 2 times.
 28937 cell_to_v_store_[4 * t] == VERTEX_AT_INFINITY ||
1875
2/2
✓ Branch 0 taken 166 times.
✓ Branch 1 taken 28769 times.
 28935 cell_to_v_store_[4 * t + 1] == VERTEX_AT_INFINITY ||
1876
2/2
✓ Branch 0 taken 92 times.
✓ Branch 1 taken 28677 times.
 28769 cell_to_v_store_[4 * t + 2] == VERTEX_AT_INFINITY ||
1877
2/2
✓ Branch 0 taken 138 times.
✓ Branch 1 taken 28539 times.
 28677 cell_to_v_store_[4 * t + 3] == VERTEX_AT_INFINITY
1878 28937 ) ;
1879 }
1880
1881
1882 /**
1883 * \brief Returns the local index of a vertex by
1884 * facet and by local vertex index in the facet.
1885 * \details
1886 * tet facet vertex is such that the tetrahedron
1887 * formed with:
1888 * - vertex lv
1889 * - tet_facet_vertex(lv,0)
1890 * - tet_facet_vertex(lv,1)
1891 * - tet_facet_vertex(lv,2)
1892 * has the same orientation as the original tetrahedron for
1893 * any vertex lv.
1894 * \param[in] f local facet index, in (0,1,2,3)
1895 * \param[in] v local vertex index, in (0,1,2)
1896 * \return the local tetrahedron vertex index of
1897 * vertex \p v in facet \p f
1898 */
1899 static index_t tet_facet_vertex(index_t f, index_t v) {
1900 geo_debug_assert(f < 4);
1901 geo_debug_assert(v < 3);
1902 73104 return index_t(tet_facet_vertex_[f][v]);
1903 }
1904
1905 /**
1906 * \brief Gets the index of a vertex of a tetrahedron
1907 * \param[in] t index of the tetrahedron
1908 * \param[in] lv local vertex (0,1,2 or 3) index in \p t
1909 * \return the global index of the \p lv%th vertex of tetrahedron \p t
1910 * or -1 if the vertex is at infinity
1911 */
1912 index_t tet_vertex(index_t t, index_t lv) const {
1913 geo_debug_assert(t < max_t());
1914 geo_debug_assert(lv < 4);
1915
2/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 9 taken 476 times.
✓ Branch 10 taken 15112 times.
 589404 return cell_to_v_store_[4 * t + lv];
1916 }
1917
1918 /**
1919 * \brief Finds the index of the vertex in a tetrahedron.
1920 * \param[in] t the tetrahedron
1921 * \param[in] v the vertex
1922 * \return iv such that tet_vertex(t,v)==iv
1923 * \pre \p t is incident to \p v
1924 */
1925 index_t find_tet_vertex(index_t t, index_t v) const {
1926 geo_debug_assert(t < max_t());
1927 // Find local index of v in tetrahedron t vertices.
1928 const index_t* T = &(cell_to_v_store_[4 * t]);
1929 return find_4(T,v);
1930 }
1931
1932
1933 /**
1934 * \brief Gets the index of a vertex of a tetrahedron
1935 * \param[in] t index of the tetrahedron
1936 * \param[in] lv local vertex (0,1,2 or 3) index in \p t
1937 * \return the global index of the \p lv%th vertex of tetrahedron \p t
1938 * \pre Vertex \p lv of tetrahedron \p t is not at infinity
1939 */
1940 index_t finite_tet_vertex(index_t t, index_t lv) const {
1941 geo_debug_assert(t < max_t());
1942 geo_debug_assert(lv < 4);
1943 geo_debug_assert(cell_to_v_store_[4 * t + lv] != NO_INDEX);
1944
8/8
✓ Branch 4 taken 13781 times.
✓ Branch 5 taken 463 times.
✓ Branch 6 taken 1405 times.
✓ Branch 7 taken 12376 times.
✓ Branch 8 taken 510 times.
✓ Branch 9 taken 11866 times.
✓ Branch 10 taken 363 times.
✓ Branch 11 taken 11503 times.
 46802 return cell_to_v_store_[4 * t + lv];
1945 }
1946
1947 /**
1948 * \brief Sets a tetrahedron-to-vertex adjacency.
1949 * \param[in] t index of the tetrahedron
1950 * \param[in] lv local vertex index (0,1,2 or 3) in \p t
1951 * \param[in] v global index of the vertex
1952 */
1953 void set_tet_vertex(index_t t, index_t lv, index_t v) {
1954 geo_debug_assert(t < max_t());
1955 geo_debug_assert(lv < 4);
1956 geo_debug_assert(owns_tet(t));
1957 56785 cell_to_v_store_[4 * t + lv] = v;
1958 }
1959
1960 /**
1961 * \brief Gets the index of a tetrahedron adjacent to another one.
1962 * \param[in] t index of the tetrahedron
1963 * \param[in] lf local facet (0,1,2 or 3) index in \p t
1964 * \return the tetrahedron adjacent to \p t accorss facet \p lf
1965 */
1966 index_t tet_adjacent(index_t t, index_t lf) const {
1967 geo_debug_assert(t < max_t());
1968 geo_debug_assert(lf < 4);
1969
10/16
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4043 times.
✓ Branch 5 taken 1305 times.
✓ Branch 6 taken 8055 times.
✓ Branch 7 taken 8022 times.
✓ Branch 8 taken 476 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 78761 times.
✓ Branch 12 taken 2874 times.
✓ Branch 13 taken 3951 times.
✓ Branch 14 taken 104735 times.
✓ Branch 15 taken 122405 times.
 402367 index_t result = cell_to_cell_store_[4 * t + lf];
1970 return result;
1971 }
1972
1973 /**
1974 * \brief Sets a tetrahedron-to-tetrahedron adjacency.
1975 * \param[in] t1 index of the first tetrahedron
1976 * \param[in] lf1 local facet index (0,1,2 or 3) in t1
1977 * \param[in] t2 index of the tetrahedron
1978 * adjacent to \p t1 across \p lf1
1979 */
1980 void set_tet_adjacent(index_t t1, index_t lf1, index_t t2) {
1981 geo_debug_assert(t1 < max_t());
1982 geo_debug_assert(t2 < max_t());
1983 geo_debug_assert(lf1 < 4);
1984 geo_debug_assert(owns_tet(t1));
1985 geo_debug_assert(owns_tet(t2));
1986
2/2
✓ Branch 0 taken 49226 times.
✓ Branch 1 taken 18514 times.
 135496 cell_to_cell_store_[4 * t1 + lf1] = t2;
1987 8022 }
1988
1989 /**
1990 * \brief Finds the index of the facet across which t1 is
1991 * adjacent to t2.
1992 * \param[in] t1 first tetrahedron
1993 * \param[in] t2 second tetrahedron
1994 * \return f such that tet_adjacent(t1,f)==t2
1995 * \pre \p t1 and \p t2 are adjacent
1996 */
1997 index_t find_tet_adjacent(index_t t1, index_t t2) const {
1998 geo_debug_assert(t1 < max_t());
1999 geo_debug_assert(t2 < max_t());
2000 geo_debug_assert(t1 != t2);
2001
2002 // Find local index of t2 in tetrahedron t1 adajcent tets.
2003 const index_t* T = &(cell_to_cell_store_[4 * t1]);
2004 index_t result = find_4(T,t2);
2005
2006 // Sanity check: make sure that t1 is adjacent to t2
2007 // only once!
2008 geo_debug_assert(tet_adjacent(t1,(result+1)%4) != t2);
2009 geo_debug_assert(tet_adjacent(t1,(result+2)%4) != t2);
2010 geo_debug_assert(tet_adjacent(t1,(result+3)%4) != t2);
2011 return result;
2012 }
2013
2014 /**
2015 * Gets the local facet index incident to an
2016 * oriented halfedge.
2017 * \param[in] t index of the tetrahedron
2018 * \param[in] v1 global index of the first extremity
2019 * \param[in] v2 global index of the second extremity
2020 * \return the local index of the facet incident to
2021 * the oriented edge \p v1, \p v2.
2022 */
2023 10384 index_t get_facet_by_halfedge(index_t t, index_t v1, index_t v2) const {
2024 geo_debug_assert(t < max_t());
2025 geo_debug_assert(v1 != v2);
2026 // Find local index of v1 and v2 in tetrahedron t
2027
2/2
✓ Branch 0 taken 7941 times.
✓ Branch 1 taken 2443 times.
 10384 const index_t* T = &(cell_to_v_store_[4 * t]);
2028
2029 index_t lv1, lv2;
2030 lv1 = find_4(T,v1);
2031 lv2 = find_4(T,v2);
2032 geo_debug_assert(lv1 != lv2);
2033 10384 return index_t(halfedge_facet_[lv1][lv2]);
2034 }
2035
2036
2037 /**
2038 * Gets the local facet indices incident to an
2039 * oriented halfedge.
2040 * \param[in] t index of the tetrahedron
2041 * \param[in] v1 global index of the first extremity
2042 * \param[in] v2 global index of the second extremity
2043 * \param[out] f12 the local index of the facet
2044 * indicent to the halfedge [v1,v2]
2045 * \param[out] f21 the local index of the facet
2046 * indicent to the halfedge [v2,v1]
2047 */
2048 8022 void get_facets_by_halfedge(
2049 index_t t, index_t v1, index_t v2,
2050 index_t& f12, index_t& f21
2051 ) const {
2052 geo_debug_assert(t < max_t());
2053 geo_debug_assert(v1 != v2);
2054
2055 // Find local index of v1 and v2 in tetrahedron t
2056 // The following expression is 10% faster than using
2057 // if() statements (multiply by boolean result of test).
2058 // Thank to Laurent Alonso for this idea.
2059
2/2
✓ Branch 0 taken 5910 times.
✓ Branch 1 taken 2112 times.
 8022 const index_t* T = &(cell_to_v_store_[4 * t]);
2060
2061
4/4
✓ Branch 0 taken 5910 times.
✓ Branch 1 taken 2112 times.
✓ Branch 2 taken 5910 times.
✓ Branch 3 taken 2112 times.
 13932 index_t lv1 = index_t((T[1] == v1) | ((T[2] == v1) * 2) | ((T[3] == v1) * 3));
2062
4/4
✓ Branch 0 taken 6079 times.
✓ Branch 1 taken 1943 times.
✓ Branch 2 taken 6247 times.
✓ Branch 3 taken 1775 times.
 14101 index_t lv2 = index_t((T[1] == v2) | ((T[2] == v2) * 2) | ((T[3] == v2) * 3));
2063 geo_debug_assert(lv1 != 0 || T[0] == v1);
2064 geo_debug_assert(lv2 != 0 || T[0] == v2);
2065 geo_debug_assert(lv1 != lv2);
2066
2067 8022 f12 = index_t(halfedge_facet_[lv1][lv2]);
2068 8022 f21 = index_t(halfedge_facet_[lv2][lv1]);
2069 8022 }
2070
2071 /**
2072 * \brief Symbolic value of the cell_next_ field
2073 * that indicates the end of list in a linked
2074 * list of tetrahedra.
2075 */
2076 static constexpr index_t END_OF_LIST = NO_INDEX;
2077
2078 /**
2079 * \brief Symbolic value of the cell_next_ field
2080 * for a tetrahedron that is not in a list.
2081 */
2082 static constexpr index_t NOT_IN_LIST = index_t(-2);
2083
2084 /**
2085 * \brief Symbolic value for t2v_[] indicating a deleted tetrahedron.
2086 */
2087 static constexpr index_t VERTEX_OF_DELETED_TET = index_t(-2);
2088
2089 /**
2090 * \brief Gets the number of vertices.
2091 * \return the number of vertices in this Delaunay
2092 */
2093 index_t nb_vertices() const {
2094
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 8 return nb_vertices_;
2095 }
2096
2097 /**
2098 * \brief Tests whether a tetrahedron belongs to a linked
2099 * list.
2100 * \details Tetrahedra can be linked, it is used to manage
2101 * both the free list that recycles deleted tetrahedra,
2102 * the conflict region and the list of newly created
2103 * tetrahedra.
2104 * \param[in] t the index of the tetrahedron
2105 * \retval true if tetrahedron \p t belongs to a linked list
2106 * \retval false otherwise
2107 */
2108 bool tet_is_in_list(index_t t) const {
2109 geo_debug_assert(t < max_t());
2110
10/24
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 28937 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4015 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 4015 times.
✓ Branch 14 taken 4019 times.
✓ Branch 15 taken 1 times.
✓ Branch 16 taken 476 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 32558 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 32558 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 15228 times.
✓ Branch 23 taken 1536 times.
 123343 return (cell_next_[t] != NOT_IN_LIST);
2111 }
2112
2113 /**
2114 * \brief Gets the index of a successor of a tetrahedron.
2115 * \details Tetrahedra can be linked, it is used to manage
2116 * both the free list that recycles deleted tetrahedra.
2117 * \param[in] t the index of the tetrahedron
2118 * \retval END_OF_LIST if the end of the list is reached
2119 * \retval the index of the successor of
2120 * tetrahedron \t otherwise
2121 * \pre tet_is_in_list(t)
2122 */
2123 index_t tet_next(index_t t) const {
2124 geo_debug_assert(t < max_t());
2125 geo_debug_assert(tet_is_in_list(t));
2126 73098 return cell_next_[t];
2127 }
2128
2129
2130 index_t tet_thread(index_t t) const {
2131 geo_debug_assert(t < max_t());
2132 return cell_status_.cell_thread(t);
2133 }
2134
2135 /**
2136 * \brief Adds a tetrahedron to a linked list.
2137 * \details Tetrahedra can be linked, it is used to manage
2138 * the free list that recycles deleted tetrahedra.
2139 * \param[in] t the index of the tetrahedron
2140 * \param[in,out] first first item of the list or END_OF_LIST if
2141 * the list is empty
2142 * \param[in,out] last last item of the list or END_OF_LIST if
2143 * the list is empty
2144 */
2145 void add_tet_to_list(index_t t, index_t& first, index_t& last) {
2146 geo_debug_assert(t < max_t());
2147 geo_debug_assert(!tet_is_in_list(t));
2148 geo_debug_assert(owns_tet(t));
2149 if(last == END_OF_LIST) {
2150 geo_debug_assert(first == END_OF_LIST);
2151 first = last = t;
2152 cell_next_[t] = END_OF_LIST;
2153 } else {
2154 cell_next_[t] = first;
2155 first = t;
2156 }
2157 }
2158
2159 /**
2160 * \brief Removes a tetrahedron from the linked list it
2161 * belongs to.
2162 * \details Tetrahedra can be linked, it is used to manage
2163 * the free list that recycles deleted tetrahedra.
2164 * \param[in] t the index of the tetrahedron
2165 */
2166 void remove_tet_from_list(index_t t) {
2167 geo_debug_assert(t < max_t());
2168 geo_debug_assert(tet_is_in_list(t));
2169 geo_debug_assert(owns_tet(t));
2170 73098 cell_next_[t] = NOT_IN_LIST;
2171 }
2172
2173
2174 /**
2175 * \brief Creates a new tetrahedron.
2176 * \details Uses either a tetrahedron recycled
2177 * from the free list, or creates a new one by
2178 * expanding the two indices arrays.
2179 * \return the index of the newly created tetrahedron
2180 */
2181 73098 index_t new_tetrahedron() {
2182 // If the memory pool is full, then we expand it.
2183 // This cannot be done when running multiple threads.
2184
2/2
✓ Branch 0 taken 15645 times.
✓ Branch 1 taken 57453 times.
 73098 if(first_free_ == END_OF_LIST) {
2185 geo_debug_assert(!Process::is_running_threads());
2186
2187
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 15640 times.
 15645 if(
2188
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15645 times.
 15645 master_->cell_to_v_store_.size() ==
2189 master_->cell_to_v_store_.capacity()
2190 ) {
2191 5 master_->nb_reallocations_++;
2192 }
2193
2194 15645 master_->cell_to_v_store_.resize(
2195 master_->cell_to_v_store_.size() + 4, NO_INDEX
2196 );
2197
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15645 times.
 15645 master_->cell_to_cell_store_.resize(
2198 master_->cell_to_cell_store_.size() + 4, NO_INDEX
2199 );
2200
2/2
✓ Branch 0 taken 15640 times.
✓ Branch 1 taken 5 times.
 15645 master_->cell_next_.push_back(END_OF_LIST);
2201 15645 master_->cell_status_.grow();
2202 15645 ++nb_free_;
2203 15645 ++max_t_;
2204 15645 first_free_ = master_->cell_status_.size() - 1;
2205 }
2206
2207 73098 acquire_and_mark_tet_as_created(first_free_);
2208 73098 index_t result = first_free_;
2209
2210 73098 first_free_ = tet_next(first_free_);
2211 remove_tet_from_list(result);
2212
2213 73098 cell_to_cell_store_[4 * result] = NO_INDEX;
2214 73098 cell_to_cell_store_[4 * result + 1] = NO_INDEX;
2215 73098 cell_to_cell_store_[4 * result + 2] = NO_INDEX;
2216 73098 cell_to_cell_store_[4 * result + 3] = NO_INDEX;
2217
2218 73098 used_tets_end_ = std::max(used_tets_end_, result+1);
2219
2220 73098 --nb_free_;
2221 73098 return result;
2222 }
2223
2224 /**
2225 * \brief Creates a new tetrahedron.
2226 * \details Sets the vertices. Adjacent tetrahedra index are
2227 * left uninitialized. Uses either a tetrahedron recycled
2228 * from the free list, or creates a new one by
2229 * expanding the two indices arrays.
2230 * \param[in] v1 index of the first vertex
2231 * \param[in] v2 index of the second vertex
2232 * \param[in] v3 index of the third vertex
2233 * \param[in] v4 index of the fourth vertex
2234 * \return the index of the newly created tetrahedron
2235 */
2236 73098 index_t new_tetrahedron(
2237 index_t v1, index_t v2,
2238 index_t v3, index_t v4
2239 ) {
2240 73098 index_t result = new_tetrahedron();
2241 73098 cell_to_v_store_[4 * result] = v1;
2242 73098 cell_to_v_store_[4 * result + 1] = v2;
2243 73098 cell_to_v_store_[4 * result + 2] = v3;
2244 73098 cell_to_v_store_[4 * result + 3] = v4;
2245 73098 return result;
2246 }
2247
2248 /**
2249 * \brief Finds the index of an integer in an array of four integers.
2250 * \param[in] T a const pointer to an array of four integers
2251 * \param[in] v the integer to retrieve in \p T
2252 * \return the index (0,1,2 or 3) of \p v in \p T
2253 * \pre The four entries of \p T are different and one of them is
2254 * equal to \p v.
2255 */
2256 static index_t find_4(const index_t* T, index_t v) {
2257 // The following expression is 10% faster than using
2258 // if() statements. This uses the C++ norm, that
2259 // ensures that the 'true' boolean value converted to
2260 // an int is always 1. With most compilers, this avoids
2261 // generating branching instructions.
2262 // Thank to Laurent Alonso for this idea.
2263 // Note: Laurent also has this version:
2264 // (T[0] != v)+(T[2]==v)+2*(T[3]==v)
2265 // that avoids a *3 multiply, but it is not faster in
2266 // practice.
2267 111640 index_t result = index_t(
2268
22/22
✓ Branch 0 taken 7308 times.
✓ Branch 1 taken 2454 times.
✓ Branch 2 taken 6769 times.
✓ Branch 3 taken 2993 times.
✓ Branch 4 taken 5415 times.
✓ Branch 5 taken 2607 times.
✓ Branch 6 taken 4043 times.
✓ Branch 7 taken 1305 times.
✓ Branch 8 taken 4015 times.
✓ Branch 9 taken 1333 times.
✓ Branch 10 taken 7941 times.
✓ Branch 11 taken 2443 times.
✓ Branch 12 taken 7807 times.
✓ Branch 13 taken 2577 times.
✓ Branch 14 taken 7551 times.
✓ Branch 15 taken 2833 times.
✓ Branch 16 taken 8224 times.
✓ Branch 17 taken 2160 times.
✓ Branch 18 taken 49226 times.
✓ Branch 19 taken 18514 times.
✓ Branch 20 taken 50654 times.
✓ Branch 21 taken 17086 times.
 239801 (T[1] == v) | ((T[2] == v) * 2) | ((T[3] == v) * 3)
2269 );
2270 // Sanity check, important if it was T[0], not explicitly
2271 // tested (detects input that does not meet the precondition).
2272 geo_debug_assert(T[result] == v);
2273 return result;
2274 }
2275
2276
2277 /**
2278 * \brief Finds the index of an integer in an array of four integers.
2279 * \param[in] T a const pointer to an array of four integers
2280 * \param[in] v the real vertex to retrieve in \p T
2281 * \return the index (0,1,2 or 3) of \p v in \p T
2282 * \pre The four entries of \p T are different and one of them is
2283 * equal to \p v.
2284 */
2285 index_t find_4_periodic(const index_t* T, index_t v) const {
2286
2287 // v needs to be a real vertex.
2288 geo_debug_assert(periodic_vertex_instance(v) == 0);
2289
2290 geo_debug_assert(
2291 T[0] != NO_INDEX && T[1] != NO_INDEX &&
2292 T[2] != NO_INDEX && T[3] != NO_INDEX
2293 );
2294
2295 // The following expression is 10% faster than using
2296 // if() statements. This uses the C++ norm, that
2297 // ensures that the 'true' boolean value converted to
2298 // an int is always 1. With most compilers, this avoids
2299 // generating branching instructions.
2300 // Thank to Laurent Alonso for this idea.
2301 // Note: Laurent also has this version:
2302 // (T[0] != v)+(T[2]==v)+2*(T[3]==v)
2303 // that avoids a *3 multiply, but it is not faster in
2304 // practice.
2305 index_t result = index_t(
2306 ( periodic_vertex_real(index_t(T[1])) == v) |
2307 ((periodic_vertex_real(index_t(T[2])) == v) * 2) |
2308 ((periodic_vertex_real(index_t(T[3])) == v) * 3)
2309 );
2310
2311 // Sanity check, important if it was T[0], not explicitly
2312 // tested (detects input that does not meet the precondition).
2313 geo_debug_assert(periodic_vertex_real(index_t(T[result])) == v);
2314 return result;
2315 }
2316
2317
2318 /**
2319 * \brief Wakes up all the threads that are waiting for
2320 * this thread.
2321 */
2322 void send_event() {
2323 4023 cond_.notify_all();
2324 }
2325
2326 /**
2327 * \brief Waits for a thread.
2328 * \details Sleeps until thread \p t calls send_event().
2329 * \param[in] t index of the thread
2330 * \pre t < nb_threads()
2331 */
2332 ✗ void wait_for_event(index_t t) {
2333 // Fixed by Hiep Vu: enlarged critical section (contains
2334 // now the test (!thrd->finished)
2335 PeriodicDelaunay3dThread* thrd = thread(t);
2336 // RAII: ctor locks, dtor unlocks
2337 ✗ std::unique_lock<std::mutex> L(thrd->mutex_);
2338 ✗ if(!thrd->finished_) {
2339 ✗ thrd->cond_.wait(L);
2340 }
2341 ✗ }
2342
2343 /****** iterative stellate_conflict_zone *****************/
2344
2345 /**
2346 * \brief Used to represent the stack in the
2347 * (de-recursified) stellate_conflict_zone_iterative()
2348 * function.
2349 */
2350 class StellateConflictStack {
2351 public:
2352
2353 /**
2354 * \brief Pushes a new frame onto the stack.
2355 * \details This also creates the local variables (they are
2356 * left uninitialized).
2357 * \param[in] t1 index of a tetrahedron on the border of
2358 * the conflict zone
2359 * \param[in] t1fbord index of the facet of \p t1 that is
2360 * on the border of the conflict zone
2361 * \param[in] t1fprev index of the facet of \p t1 that we
2362 * come from, or NO_INDEX if \p t1 is the first tetrahedron
2363 */
2364 void push(index_t t1, index_t t1fbord, index_t t1fprev) {
2365 5348 store_.resize(store_.size()+1);
2366 5348 top().t1 = t1;
2367 5348 top().t1fbord = Numeric::uint8(t1fbord);
2368 5348 top().t1fprev = Numeric::uint8(t1fprev);
2369 33 }
2370
2371 /**
2372 * \brief Saves local variables into the current stack frame.
2373 * \param[in] new_t the index of the newly created tetrahedron
2374 * \param[in] t1ft2 the facet of t1 that is adjacent to t2
2375 * \param[in] t2ft1 the facet of t2 that is adjacent to t1
2376 */
2377 void save_locals(index_t new_t, index_t t1ft2, index_t t2ft1) {
2378 geo_debug_assert(!empty());
2379 5315 top().new_t = new_t;
2380 5315 top().t1ft2 = Numeric::uint8(t1ft2);
2381 5315 top().t2ft1 = Numeric::uint8(t2ft1);
2382 }
2383
2384 /**
2385 * \brief Gets the parameters from the current stack frame.
2386 * \param[out] t1 index of a tetrahedron on the border of
2387 * the conflict zone
2388 * \param[out] t1fbord index of the facet of \p t1 that is
2389 * on the border of the conflict zone
2390 * \param[out] t1fprev index of the facet of \p t1 that we
2391 * come from, or NO_INDEX if \p t1 is the first tetrahedron
2392 */
2393 void get_parameters(
2394 index_t& t1, index_t& t1fbord, index_t& t1fprev
2395 ) const {
2396 geo_debug_assert(!empty());
2397 10663 t1 = top().t1;
2398 10663 t1fbord = index_t(top().t1fbord);
2399 5348 t1fprev = index_t(top().t1fprev);
2400 }
2401
2402
2403 /**
2404 * \brief Gets the local variables from the current stack frame.
2405 * \param[out] new_t the index of the newly created tetrahedron
2406 * \param[out] t1ft2 the facet of t1 that is adjacent to t2
2407 * \param[out] t2ft1 the facet of t2 that is adjacent to t1
2408 */
2409 void get_locals(
2410 index_t& new_t, index_t& t1ft2, index_t& t2ft1
2411 ) const {
2412 geo_debug_assert(!empty());
2413 5315 new_t = top().new_t;
2414 5315 t1ft2 = index_t(top().t1ft2);
2415 5315 t2ft1 = index_t(top().t2ft1);
2416 }
2417
2418 /**
2419 * \brief Pops a stack frame.
2420 */
2421 void pop() {
2422 geo_debug_assert(!empty());
2423 store_.pop_back();
2424 }
2425
2426 /**
2427 * \brief Tests whether the stack is empty.
2428 * \retval true if the stack is empty
2429 * \retval false otherwise
2430 */
2431 bool empty() const {
2432 return store_.empty();
2433 }
2434
2435 private:
2436
2437 /**
2438 * \brief The parameters and local
2439 * variables stored in a stack frame.
2440 */
2441 struct Frame {
2442 // Parameters
2443 index_t t1;
2444 index_t new_t;
2445 Numeric::uint8 t1fbord ;
2446
2447 // Local variables
2448 Numeric::uint8 t1fprev ;
2449 Numeric::uint8 t1ft2 ;
2450 Numeric::uint8 t2ft1 ;
2451 };
2452
2453 /**
2454 * \brief Gets the top of the stack.
2455 * \return a modifiable reference to the Frame on
2456 * the top of the stack
2457 * \pre !empty()
2458 */
2459 Frame& top() {
2460 geo_debug_assert(!empty());
2461 return *store_.rbegin();
2462 }
2463
2464 /**
2465 * \brief Gets the top of the stack.
2466 * \return a const reference to the Frame on
2467 * the top of the stack
2468 * \pre !empty()
2469 */
2470 const Frame& top() const {
2471 geo_debug_assert(!empty());
2472 return *store_.rbegin();
2473 }
2474
2475 std::vector<Frame> store_;
2476 };
2477
2478 /**
2479 * \brief Creates a star of tetrahedra filling the conflict
2480 * zone.
2481 * \details For each tetrahedron facet on the border of the
2482 * conflict zone, a new tetrahedron is created, resting on
2483 * the facet and incident to vertex \p v. The function is
2484 * called recursively until the entire conflict zone is filled.
2485 * \param[in] v the index of the point to be inserted
2486 * \param[in] t1 index of a tetrahedron on the border
2487 * of the conflict zone.
2488 * \param[in] t1fbord index of the facet along which \p t_bndry
2489 * is incident to the border of the conflict zone
2490 * \param[in] t1fprev the facet of \p t_bndry connected to the
2491 * tetrahedron that \p t_bndry was reached from, or NO_INDEX
2492 * if it is the first tetrahedron.
2493 * \return the index of one the newly created tetrahedron
2494 */
2495
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 33 times.
 33 index_t stellate_conflict_zone_iterative(
2496 index_t v, index_t t1, index_t t1fbord,
2497 index_t t1fprev = NO_INDEX
2498 ) {
2499 // This function is de-recursified because some degenerate
2500 // inputs can cause stack overflow (system stack is limited to
2501 // a few megs). For instance, it can happen when a large number
2502 // of points are on the same sphere exactly.
2503
2504 // To de-recursify, it uses class StellateConflictStack
2505 // that emulates system's stack for storing functions's
2506 // parameters and local variables in all the nested stack
2507 // frames.
2508
2509 S2_.push(t1, t1fbord, t1fprev);
2510
2511 index_t new_t; // the newly created tetrahedron.
2512
2513 index_t t1ft2; // traverses the 4 facets of t1.
2514
2515 index_t t2; // the tetrahedron on the border of
2516 // the conflict zone that shares an
2517 // edge with t1 along t1ft2.
2518
2519 index_t t2fbord; // the facet of t2 on the border of
2520 // the conflict zone.
2521
2522 index_t t2ft1; // the facet of t2 that is incident to t1.
2523
2524 5348 entry_point:
2525 S2_.get_parameters(t1, t1fbord, t1fprev);
2526
2527
2528 geo_debug_assert(owns_tet(t1));
2529 geo_debug_assert(tet_adjacent(t1,t1fbord) != NO_INDEX);
2530 geo_debug_assert(owns_tet(tet_adjacent(t1,t1fbord)));
2531 geo_debug_assert(tet_is_marked_as_conflict(t1));
2532 geo_debug_assert(
2533 !tet_is_marked_as_conflict(tet_adjacent(t1,t1fbord))
2534 );
2535
2536 // Create new tetrahedron with same vertices as t_bndry
2537
2538 5348 new_t = new_tetrahedron(
2539 tet_vertex(t1,0),
2540 tet_vertex(t1,1),
2541 tet_vertex(t1,2),
2542 tet_vertex(t1,3)
2543 );
2544
2545 index_t tbord = tet_adjacent(t1,t1fbord);
2546
2547 // We generate the tetrahedron with the three vertices
2548 // of the tet outside the conflict zone and the newly
2549 // created vertex in the local frame of the tet outside
2550 // the conflict zone.
2551
2552 // Replace in new_t the vertex opposite to t1fbord with v
2553 set_tet_vertex(new_t, t1fbord, v);
2554
2555 {
2556 // Connect new_t with t1's neighbor across t1fbord
2557 set_tet_adjacent(new_t, t1fbord, tbord);
2558 set_tet_adjacent(tbord, find_tet_adjacent(tbord,t1), new_t);
2559 }
2560
2561 // Lookup new_t's neighbors across its three other
2562 // facets and connect them
2563
2/2
✓ Branch 0 taken 21392 times.
✓ Branch 1 taken 5348 times.
 26740 for(t1ft2=0; t1ft2<4; ++t1ft2) {
2564
2565
4/4
✓ Branch 0 taken 16077 times.
✓ Branch 1 taken 5315 times.
✓ Branch 2 taken 8055 times.
✓ Branch 3 taken 8022 times.
 21392 if(t1ft2 == t1fprev || tet_adjacent(new_t,t1ft2) != NO_INDEX) {
2566 13370 continue;
2567 }
2568
2569 // Get t1's neighbor along the border of the conflict zone
2570
2/2
✓ Branch 1 taken 5315 times.
✓ Branch 2 taken 2707 times.
 8022 if(!get_neighbor_along_conflict_zone_border(
2571 t1,t1fbord,t1ft2, t2,t2fbord,t2ft1
2572 )) {
2573 // If t1's neighbor is not a new tetrahedron,
2574 // create a new tetrahedron through a recursive call.
2575
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5315 times.
 5315 S2_.save_locals(new_t, t1ft2, t2ft1);
2576
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5315 times.
 5315 S2_.push(t2, t2fbord, t2ft1);
2577 5315 goto entry_point;
2578
2579 return_point:
2580 // This is the return value of the called function.
2581 index_t result = new_t;
2582 S2_.pop();
2583
2584 // Special case: we were in the outermost frame,
2585 // then we (truly) return from the function.
2586 5348 if(S2_.empty()) {
2587 33 return result;
2588 }
2589
2590 S2_.get_parameters(t1, t1fbord, t1fprev);
2591 S2_.get_locals(new_t, t1ft2, t2ft1);
2592 5315 t2 = result;
2593 }
2594
2595 8022 set_tet_adjacent(t2, t2ft1, new_t);
2596 set_tet_adjacent(new_t, t1ft2, t2);
2597 }
2598
2599 // Except for the initial call (see "Special case" above),
2600 // the nested calls all come from the same location,
2601 // thus there is only one possible return point
2602 // (no need to push any return address).
2603 5348 goto return_point;
2604 }
2605
2606 /**
2607 * \brief Finds the neighbor of a tetrahedron on the border of the
2608 * conflict zone.
2609 * \details This function is used by stellate_conflict_zone_iterative()
2610 * \param[in] t1 a tetrahedron on the border of the conflict zone
2611 * \param[in] t1fborder the local facet index of \p t1 along which it
2612 * is on the border of the conflict zone
2613 * \param[in] t1ft2 the local facet index of \p t1 that will be
2614 * traversed
2615 * \param[out] t2 a tetrahedron on the border of the conflict zone,
2616 * with an edge common to facets \p t1fborder and \p t1ft2 of
2617 * tetrahedron \p t1
2618 * \param[out] t2fborder the local facet index of \p t2 along which it
2619 * is on the border of the conflict zone
2620 * \param[out] t2ft1 the local index of the facet of \p t2 that has a
2621 * common edge with facets \p t1fborder and \p t1ft2 of tetrahedron
2622 * \p t1
2623 * \retval true if \p t2 is a newly created tetrahedron
2624 * \retval false if \p t2 is an old tetrahedron in conflict
2625 */
2626 8022 bool get_neighbor_along_conflict_zone_border(
2627 index_t t1,
2628 index_t t1fborder,
2629 index_t t1ft2,
2630 index_t& t2,
2631 index_t& t2fborder,
2632 index_t& t2ft1
2633 ) const {
2634
2635 // Note: this function is a bit long for an inline function,
2636 // but I observed a (modest) performance gain doing so.
2637
2638 // Find two vertices that are both on facets new_f and f1
2639 // (the edge around which we are turning)
2640 // This uses duality as follows:
2641 // Primal form (not used here):
2642 // halfedge_facet_[v1][v2] returns a facet that is incident
2643 // to both v1 and v2.
2644 // Dual form (used here):
2645 // halfedge_facet_[f1][f2] returns a vertex that both
2646 // f1 and f2 are incident to.
2647 index_t ev1 =
2648 8022 tet_vertex(t1, index_t(halfedge_facet_[t1ft2][t1fborder]));
2649 index_t ev2 =
2650 8022 tet_vertex(t1, index_t(halfedge_facet_[t1fborder][t1ft2]));
2651
2652 // Turn around edge [ev1,ev2] inside the conflict zone
2653 // until we reach again the boundary of the conflict zone.
2654 // Traversing inside the conflict zone is faster (as compared
2655 // to outside) since it traverses a smaller number of tets.
2656 index_t cur_t = t1;
2657 index_t cur_f = t1ft2;
2658 index_t next_t = tet_adjacent(cur_t,cur_f);
2659
2/2
✓ Branch 0 taken 10384 times.
✓ Branch 1 taken 8022 times.
 18406 while(tet_is_marked_as_conflict(next_t)) {
2660 geo_debug_assert(next_t != t1);
2661 cur_t = next_t;
2662 10384 cur_f = get_facet_by_halfedge(cur_t,ev1,ev2);
2663 next_t = tet_adjacent(cur_t, cur_f);
2664 }
2665
2666 // At this point, cur_t is in conflict zone and
2667 // next_t is outside the conflict zone.
2668 index_t f12,f21;
2669 8022 get_facets_by_halfedge(next_t, ev1, ev2, f12, f21);
2670
2/2
✓ Branch 0 taken 5783 times.
✓ Branch 1 taken 2239 times.
 8022 t2 = tet_adjacent(next_t,f21);
2671
2/2
✓ Branch 0 taken 5783 times.
✓ Branch 1 taken 2239 times.
 8022 index_t v_neigh_opposite = tet_vertex(next_t,f12);
2672 8022 t2ft1 = find_tet_vertex(t2, v_neigh_opposite);
2673 8022 t2fborder = cur_f;
2674
2675 // Test whether the found neighboring tet was created
2676 // (then return true) or is an old tet in conflict
2677 // (then return false).
2678 8022 return(t2 != cur_t);
2679 }
2680
2681 /**
2682 * \brief Used by the (de-recursified)
2683 * stellate_conflict_zone_iterative() function.
2684 */
2685 StellateConflictStack S2_;
2686
2687 /*************************** debugging ************************/
2688
2689 /**
2690 * \brief For debugging purposes, displays a tetrahedron adjacency.
2691 * \param[in] t index of the tetrahedron to display.
2692 * \param[in] lf local index (0,1,2 or 3) of the tetrahedron
2693 * facet adjacenty to display.
2694 */
2695 ✗ void show_tet_adjacent(index_t t, index_t lf) const {
2696 index_t adj = tet_adjacent(t, lf);
2697 ✗ if(adj != NO_INDEX) {
2698 ✗ std::cerr << (tet_is_in_list(adj) ? '*' : ' ');
2699 }
2700 std::cerr << adj;
2701 ✗ std::cerr << ' ';
2702 ✗ }
2703
2704
2705 /**
2706 * \brief For debugging purposes, displays a tetrahedron.
2707 * \param[in] t index of the tetrahedron to display.
2708 */
2709 ✗ void show_tet(index_t t) const {
2710 std::cerr << "tet"
2711 ✗ << (tet_is_in_list(t) ? '*' : ' ')
2712 << t
2713 << ", v=["
2714 << tet_vertex(t, 0)
2715 ✗ << ' '
2716 << tet_vertex(t, 1)
2717 ✗ << ' '
2718 << tet_vertex(t, 2)
2719 ✗ << ' '
2720 << tet_vertex(t, 3)
2721 ✗ << "] adj=[";
2722 ✗ show_tet_adjacent(t, 0);
2723 ✗ show_tet_adjacent(t, 1);
2724 ✗ show_tet_adjacent(t, 2);
2725 ✗ show_tet_adjacent(t, 3);
2726 ✗ std::cerr << "] ";
2727
2728 ✗ for(index_t f = 0; f < 4; ++f) {
2729 ✗ std::cerr << 'f' << f << ':';
2730 ✗ for(index_t v = 0; v < 3; ++v) {
2731 std::cerr << tet_vertex(t, tet_facet_vertex(f,v))
2732 ✗ << ',';
2733 }
2734 ✗ std::cerr << ' ';
2735 }
2736 std::cerr << std::endl;
2737 ✗ }
2738
2739 public:
2740
2741 /**
2742 * \brief For debugging purposes, tests some combinatorial properties.
2743 */
2744 ✗ void check_combinatorics(bool verbose) const {
2745 ✗ if(verbose) {
2746 std::cerr << std::endl;
2747 }
2748 bool ok = true;
2749 ✗ std::vector<bool> v_has_tet(nb_vertices(), false);
2750 ✗ for(index_t t = 0; t < max_t(); ++t) {
2751 ✗ if(tet_is_free(t)) {
2752 ✗ if(verbose) {
2753 ✗ std::cerr << "-Deleted tet: ";
2754 ✗ show_tet(t);
2755 }
2756 } else {
2757 ✗ if(verbose) {
2758 ✗ std::cerr << "Checking tet: ";
2759 ✗ show_tet(t);
2760 }
2761 ✗ for(index_t lf = 0; lf < 4; ++lf) {
2762 ✗ if(tet_adjacent(t, lf) == NO_INDEX) {
2763 std::cerr << lf << ":Missing adjacent tet"
2764 << std::endl;
2765 ok = false;
2766 ✗ } else if(tet_adjacent(t, lf) == t) {
2767 std::cerr << lf << ":Tet is adjacent to itself"
2768 << std::endl;
2769 ok = false;
2770 } else {
2771 index_t t2 = tet_adjacent(t, lf);
2772 bool found = false;
2773 ✗ for(index_t lf2 = 0; lf2 < 4; ++lf2) {
2774 ✗ if(tet_adjacent(t2, lf2) == t) {
2775 found = true;
2776 }
2777 }
2778 ✗ if(!found) {
2779 std::cerr
2780 << lf
2781 << ":Adjacent link is not bidirectional"
2782 << std::endl;
2783 ok = false;
2784 }
2785 }
2786 }
2787 index_t nb_infinite = 0;
2788 ✗ for(index_t lv = 0; lv < 4; ++lv) {
2789 ✗ if(tet_vertex(t, lv) == NO_INDEX) {
2790 ✗ ++nb_infinite;
2791 }
2792 }
2793 ✗ if(nb_infinite > 1) {
2794 ok = false;
2795 std::cerr << "More than one infinite vertex"
2796 << std::endl;
2797 }
2798 }
2799 ✗ for(index_t lv = 0; lv < 4; ++lv) {
2800 index_t v = tet_vertex(t, lv);
2801 ✗ if(v != NO_INDEX && v != NOT_IN_LIST) {
2802 v_has_tet[periodic_vertex_real(v)] = true;
2803 }
2804 }
2805 }
2806
2807 index_t nb_v = nb_vertices();
2808 ✗ for(index_t v = 0; v < nb_v; ++v) {
2809 ✗ if(!v_has_tet[v]) {
2810 ✗ if(verbose) {
2811 std::cerr << "Vertex " << v
2812 << " is isolated (duplicated ?)" << std::endl;
2813 }
2814 }
2815 }
2816 ✗ geo_assert(ok);
2817 ✗ if(verbose) {
2818 std::cerr << std::endl;
2819 }
2820 std::cerr << std::endl << "Delaunay Combi OK" << std::endl;
2821 ✗ }
2822
2823
2824 /**
2825 * \brief For debugging purposes, test some geometrical properties.
2826 */
2827 ✗ void check_geometry(bool verbose) const {
2828 bool ok = true;
2829 ✗ for(index_t t = 0; t < max_t(); ++t) {
2830 ✗ if(!tet_is_free(t)) {
2831 index_t v0 = tet_vertex(t, 0);
2832 index_t v1 = tet_vertex(t, 1);
2833 index_t v2 = tet_vertex(t, 2);
2834 index_t v3 = tet_vertex(t, 3);
2835 ✗ for(index_t v = 0; v < nb_vertices(); ++v) {
2836 vec4 p = lifted_vertex(v);
2837 ✗ if(v == v0 || v == v1 || v == v2 || v == v3) {
2838 ✗ continue;
2839 }
2840 ✗ if(tet_is_in_conflict(t, v, p)) {
2841 ok = false;
2842 ✗ if(verbose) {
2843 std::cerr << "Tet " << t <<
2844 " is in conflict with vertex " << v
2845 << std::endl;
2846
2847 ✗ std::cerr << " offending tet: ";
2848 ✗ show_tet(t);
2849 }
2850 }
2851 }
2852 }
2853 }
2854 ✗ geo_assert(ok);
2855 std::cerr << std::endl << "Delaunay Geo OK" << std::endl;
2856 ✗ }
2857
2858 private:
2859 PeriodicDelaunay3d* master_;
2860 bool periodic_;
2861 vec3 period_;
2862 index_t nb_vertices_;
2863 const double* vertices_;
2864 const double* weights_;
2865 index_t* reorder_;
2866 index_t dimension_;
2867 index_t pool_begin_;
2868 index_t pool_end_;
2869 index_t max_t_;
2870 index_t used_tets_end_;
2871
2872 vector<index_t>& cell_to_v_store_;
2873 vector<index_t>& cell_to_cell_store_;
2874 vector<index_t>& cell_next_;
2875 CellStatusArray& cell_status_;
2876
2877 index_t first_free_;
2878 index_t nb_free_;
2879 bool memory_overflow_;
2880
2881 /** \brief used by find_conflict_zone_iterative() */
2882 struct SFrame {
2883
2884 ✗ SFrame() {
2885 }
2886
2887 SFrame(
2888 index_t t_in,
2889 index_t v_in,
2890 const vec4& p_in
2891 56785 ) :
2892 56785 t(t_in),
2893 56785 v(v_in),
2894
1/2
✓ Branch 2 taken 54010 times.
✗ Branch 3 not taken.
 56785 p(p_in) {
2895 }
2896
2897 SFrame(
2898 const SFrame& rhs
2899 56975 ):
2900 56975 t(rhs.t),
2901 56975 v(rhs.v),
2902 56975 p(rhs.p) {
2903 }
2904
2905 SFrame& operator=(const SFrame& rhs) {
2906 t = rhs.t;
2907 v = rhs.v;
2908 p = rhs.p;
2909 return *this;
2910 }
2911
2912 index_t t;
2913 index_t v;
2914 vec4 p;
2915 };
2916
2917 vector<SFrame> S_;
2918 index_t nb_tets_to_create_;
2919 index_t t_boundary_; // index of a tet,facet on the bndry
2920 index_t f_boundary_; // of the conflict zone.
2921
2922 bool direction_;
2923 index_t work_begin_;
2924 index_t work_rbegin_;
2925 index_t b_hint_;
2926 index_t e_hint_;
2927 bool finished_;
2928
2929 // Whenever acquire_tet() is unsuccessful, contains
2930 // the index of the thread that was interfering
2931 // (shifted to the left by 1 !!)
2932 CellStatusArray::thread_index_t interfering_thread_;
2933
2934 #ifdef GEO_DEBUG
2935 index_t nb_acquired_tets_;
2936 #endif
2937
2938 vector<index_t> tets_to_delete_;
2939 vector<index_t> tets_to_release_;
2940
2941 index_t nb_rollbacks_;
2942 index_t nb_failed_locate_;
2943
2944 std::condition_variable cond_;
2945 std::mutex mutex_;
2946
2947 bool abort_on_empty_cell_;
2948 bool has_empty_cells_;
2949
2950 /**
2951 * \brief Gives the indexing of tetrahedron facet
2952 * vertices.
2953 * \details tet_facet_vertex[lf][lv] gives the
2954 * local vertex index (in 0,1,2,3) from a
2955 * local facet index lf (in 0,1,2,3) and a
2956 * local vertex index within the facet (in 0,1,2).
2957 */
2958 static char tet_facet_vertex_[4][3];
2959
2960 /**
2961 * \brief Gives a local facet index by
2962 * halfedge extremities local indices.
2963 */
2964 static char halfedge_facet_[4][4];
2965
2966 /**
2967 * \brief Optimized representation of triangles on
2968 * the border of the cavity, for fast generation of
2969 * tetrahedra.
2970 */
2971 Cavity cavity_;
2972
2973 /**
2974 * \brief Statistics for locate()
2975 */
2976 index_t nb_traversed_tets_;
2977 };
2978
2979
2980 char PeriodicDelaunay3dThread::halfedge_facet_[4][4] = {
2981 {4, 2, 3, 1},
2982 {3, 4, 0, 2},
2983 {1, 3, 4, 0},
2984 {2, 0, 1, 4}
2985 };
2986
2987 // tet facet vertex is such that the tetrahedron
2988 // formed with:
2989 // vertex lv
2990 // tet_facet_vertex[lv][0]
2991 // tet_facet_vertex[lv][1]
2992 // tet_facet_vertex[lv][2]
2993 // has the same orientation as the original tetrahedron for
2994 // any vertex lv.
2995
2996 char PeriodicDelaunay3dThread::tet_facet_vertex_[4][3] = {
2997 {1, 2, 3},
2998 {0, 3, 2},
2999 {3, 0, 1},
3000 {1, 0, 2}
3001 };
3002
3003
3004 /*************************************************************************/
3005
3006 2 PeriodicDelaunay3d::PeriodicDelaunay3d(
3007 bool periodic, double period
3008 2 ) :
3009 Delaunay(3),
3010 2 periodic_(periodic),
3011 period_(period,period,period),
3012
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 weights_(nullptr),
3013 2 update_periodic_v_to_cell_(false),
3014 2 abort_on_empty_cell_(false),
3015 2 has_empty_cells_(false),
3016
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 nb_reallocations_(0),
3017
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 convex_cell_exact_predicates_(true)
3018 {
3019
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
 2 debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay");
3020
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 2 verbose_debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay_verbose");
3021
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 debug_mode_ = (debug_mode_ || verbose_debug_mode_);
3022
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 2 benchmark_mode_ = CmdLine::get_arg_bool("dbg:delaunay_benchmark");
3023
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 detailed_benchmark_mode_ =
3024
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2 CmdLine::get_arg_bool("dbg:detailed_delaunay_benchmark");
3025 2 nb_vertices_non_periodic_ = 0;
3026
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 delaunay_citations();
3027 2 }
3028
3029 ✗ PeriodicDelaunay3d::PeriodicDelaunay3d(
3030 const vec3& period, bool periodic
3031 ✗ ) :
3032 Delaunay(3),
3033 ✗ periodic_(periodic),
3034 ✗ period_(period),
3035 ✗ weights_(nullptr),
3036 ✗ update_periodic_v_to_cell_(false),
3037 ✗ has_empty_cells_(false),
3038 ✗ nb_reallocations_(0),
3039 ✗ convex_cell_exact_predicates_(true)
3040 {
3041 ✗ debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay");
3042 ✗ verbose_debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay_verbose");
3043 ✗ debug_mode_ = (debug_mode_ || verbose_debug_mode_);
3044 ✗ benchmark_mode_ = CmdLine::get_arg_bool("dbg:delaunay_benchmark");
3045 ✗ detailed_benchmark_mode_ =
3046 ✗ CmdLine::get_arg_bool("dbg:detailed_delaunay_benchmark");
3047 ✗ nb_vertices_non_periodic_ = 0;
3048 ✗ delaunay_citations();
3049 ✗ }
3050
3051 2 void PeriodicDelaunay3d::set_vertices(
3052 index_t nb_vertices, const double* vertices
3053 ) {
3054 2 has_empty_cells_ = false;
3055
3056 #ifndef GARGANTUA
3057 {
3058 Numeric::uint64 expected_max_index =
3059 Numeric::uint64(nb_vertices) * 7 * 4;
3060 if(periodic_) {
3061 expected_max_index *= 2;
3062 }
3063 if(expected_max_index > Numeric::uint64(INT32_MAX)) {
3064 Logger::err("OTM") << "indices will overflow" << std::endl;
3065 Logger::err("OTM")
3066 << "recompile with -DGARGANTUA "
3067 << "to activate 64bit indices" << std::endl;
3068 exit(0);
3069 }
3070 }
3071 #endif
3072
3073
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(periodic_) {
3074 1 PCK::set_SOS_mode(PCK::SOS_LEXICO);
3075 }
3076
3077
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 Stopwatch W("BRIO", benchmark_mode_);
3078 2 nb_vertices_non_periodic_ = nb_vertices;
3079
3080
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Delaunay::set_vertices(nb_vertices, vertices);
3081 // Reorder the points
3082
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(do_reorder_) {
3083
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 compute_BRIO_order(
3084
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 nb_vertices, vertex_ptr(0), reorder_,
3085 3, dimension(),
3086 64, 0.125,
3087 &levels_
3088 );
3089 } else {
3090 ✗ reorder_.resize(nb_vertices);
3091 ✗ for(index_t i = 0; i < nb_vertices; ++i) {
3092 ✗ reorder_[i] = i;
3093 }
3094 geo_debug_assert(levels_[0] == 0);
3095 geo_debug_assert(levels_[levels_.size()-1] == nb_vertices);
3096 }
3097 2 }
3098
3099 ✗ void PeriodicDelaunay3d::set_weights(const double* weights) {
3100 ✗ has_empty_cells_ = false;
3101 ✗ weights_ = weights;
3102 ✗ }
3103
3104 2 void PeriodicDelaunay3d::compute() {
3105
3106 2 stats_.reset();
3107
3108
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 Stopwatch W_tot("total",false);
3109
3110 2 has_empty_cells_ = false;
3111
3112
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(periodic_) {
3113
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 reorder_.resize(nb_vertices_non_periodic_);
3114 }
3115
3116 {
3117
4/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 1 times.
 4 Stopwatch W("DelInternal", detailed_benchmark_mode_);
3118
3119 2 index_t expected_tetra = nb_vertices() * 7;
3120
3121 // Everything is allocated here, including for handling
3122 // periodic boundary conditions, much later. We need to
3123 // allocate sufficient space to have good chances of
3124 // inserting most of the additional points in parallel
3125 // (in insert_with_BRIO())
3126
3127
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(periodic_) {
3128 1 expected_tetra = index_t(double(expected_tetra)* 1.2);
3129 }
3130
3131 // Allocate the tetrahedra
3132
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_v_store_.assign(expected_tetra * 4, NO_INDEX);
3133
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_cell_store_.assign(expected_tetra * 4, NO_INDEX);
3134
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_next_.assign(expected_tetra,NO_INDEX);
3135
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_status_.resize(expected_tetra);
3136
3137 // Create the threads
3138 // The maximum number of threads is limited by the number
3139 // of bits used by cell_status_ (see delaunay_sync.h)
3140 index_t nb_threads = std::min(
3141
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Process::maximum_concurrent_threads(),
3142 CellStatusArray::MAX_THREADS
3143 );
3144 2 index_t pool_size = expected_tetra / nb_threads;
3145
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (pool_size == 0) {
3146 // There are more threads than expected_tetra
3147 pool_size = 1;
3148 nb_threads = expected_tetra;
3149 }
3150 index_t pool_begin = 0;
3151 2 threads_.clear();
3152
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for(index_t t=0; t<nb_threads; ++t) {
3153 6 index_t pool_end =
3154
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
 8 (t == nb_threads - 1) ? expected_tetra
3155 : pool_begin + pool_size;
3156 threads_.push_back(
3157
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 new PeriodicDelaunay3dThread(this, pool_begin, pool_end)
3158 );
3159 pool_begin = pool_end;
3160 }
3161
3162
3163 // Create first tetrahedron and triangulate first set of points
3164 // in sequential mode.
3165
3166 PeriodicDelaunay3dThread* thread0 = thread(0);
3167
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 thread0->create_first_tetrahedron();
3168 {
3169
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 Stopwatch Wmain("DelMain", detailed_benchmark_mode_);
3170
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 insert_vertices_with_BRIO("DelMain", levels_);
3171
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 stats_.phase_0_t_ = Wmain.elapsed_time();
3172 2 }
3173
3174
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if(abort_on_empty_cell_ && has_empty_cells_) {
3175 return;
3176 }
3177
3178
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(periodic_) {
3179
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 Stopwatch W12("DelPhaseI-II", detailed_benchmark_mode_);
3180
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 handle_periodic_boundaries();
3181 1 }
3182
3183
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if(abort_on_empty_cell_ && has_empty_cells_) {
3184 return;
3185 }
3186 2 }
3187
3188
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(debug_mode_) {
3189 ✗ for(index_t i=0; i<threads_.size(); ++i) {
3190 std::cerr << i << " : " <<
3191 static_cast<PeriodicDelaunay3dThread*>(threads_[i].get())
3192 ->max_t() << std::endl;
3193 }
3194
3195 ✗ thread(0)->check_combinatorics(verbose_debug_mode_);
3196 ✗ thread(0)->check_geometry(verbose_debug_mode_);
3197 }
3198
3199 index_t nb_tets = 0;
3200 {
3201
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 nb_tets = compress();
3202
3203
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 set_arrays(
3204 nb_tets,
3205 cell_to_v_store_.data(),
3206 cell_to_cell_store_.data()
3207 );
3208
3209 // We need v_to_cell even if CICL is not stored.
3210
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(!stores_cicl()) {
3211
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 update_v_to_cell();
3212 }
3213 }
3214
3215 if(periodic_) {
3216 #ifdef GEO_DEBUG
3217 FOR(v, nb_vertices_non_periodic_) {
3218 index_t t = v_to_cell_[v];
3219 geo_assert(t == NO_INDEX || t < nb_tets);
3220 }
3221 #endif
3222 }
3223
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 stats_.total_t_ = W_tot.elapsed_time();
3224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(benchmark_mode_) {
3225 ✗ Logger::out("Delaunay") << stats_.to_string() << std::endl;
3226 }
3227 2 }
3228
3229 2 index_t PeriodicDelaunay3d::compress(bool shrink) {
3230
3231
2/4
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 4 Stopwatch W("Compress",detailed_benchmark_mode_);
3232
3233 // Compress cell_to_v_store_ and cell_to_cell_store_
3234 // (remove free and virtual tetrahedra).
3235 // Since cell_next_ is not used at this point,
3236 // we reuse it for storing the conversion array that
3237 // maps old tet indices to new tet indices
3238 // Note: tet_is_real() uses the previous value of
3239 // cell_next(), but we are processing indices
3240 // in increasing order and since old2new[t] is always
3241 // smaller or equal to t, we never overwrite a value
3242 // before needing it.
3243
3244
3245 2 PeriodicDelaunay3dThread* thread0 = thread(0);
3246
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 vector<index_t>& old2new = cell_next_;
3247 index_t nb_tets = 0;
3248 index_t nb_tets_to_delete = 0;
3249
3250 {
3251 // Classify tets in parallel (on very large data sets,
3252 // >= 100M points, it gains a little bit of time)
3253
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 parallel_for(0, thread0->max_t(), [&,this](index_t t) {
3254 if(
3255
2/2
✓ Branch 0 taken 1199 times.
✓ Branch 1 taken 1085 times.
 2284 (keep_infinite_ && !thread0->tet_is_free(t)) ||
3256
8/8
✓ Branch 0 taken 2284 times.
✓ Branch 1 taken 16764 times.
✓ Branch 2 taken 16764 times.
✓ Branch 3 taken 1199 times.
✓ Branch 4 taken 14023 times.
✓ Branch 5 taken 2741 times.
✓ Branch 6 taken 1199 times.
✓ Branch 7 taken 14023 times.
 20247 (periodic_ && thread0->tet_is_real_non_periodic(t)) ||
3257
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1199 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1199 (!periodic_ && thread0->tet_is_real(t))
3258 ) {
3259 3826 old2new[t] = 0; // keep tetrahedron
3260 } else {
3261 15222 old2new[t] = NO_INDEX; // discard tetrahedron
3262 }
3263 19048 });
3264
3265 // Compress the tet array
3266
2/2
✓ Branch 0 taken 19048 times.
✓ Branch 1 taken 2 times.
 19050 for(index_t t = 0; t < thread0->max_t(); ++t) {
3267
2/2
✓ Branch 0 taken 3826 times.
✓ Branch 1 taken 15222 times.
 19048 if(old2new[t] != NO_INDEX) {
3268
2/2
✓ Branch 0 taken 3789 times.
✓ Branch 1 taken 37 times.
 3826 if(t != nb_tets) {
3269 Memory::copy(
3270 &cell_to_v_store_[nb_tets * 4],
3271 &cell_to_v_store_[t * 4],
3272 4 * sizeof(index_t)
3273 );
3274 Memory::copy(
3275 &cell_to_cell_store_[nb_tets * 4],
3276 &cell_to_cell_store_[t * 4],
3277 4 * sizeof(index_t)
3278 );
3279 }
3280 3826 old2new[t] = nb_tets;
3281 3826 ++nb_tets;
3282 } else {
3283 15222 ++nb_tets_to_delete;
3284 }
3285 }
3286
3287
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(shrink) {
3288
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_v_store_.resize(4 * nb_tets);
3289
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 cell_to_cell_store_.resize(4 * nb_tets);
3290 }
3291
3292 // Apply permutation to cell_to_cell_ array
3293
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 parallel_for(0, 4*nb_tets, [this, &old2new](index_t i) {
3294 15304 index_t t = cell_to_cell_store_[i];
3295 geo_debug_assert(t != NO_INDEX);
3296 15304 t = old2new[t];
3297 // Note: t can be equal to -1 when a real tet is
3298 // adjacent to a virtual one (and this is how the
3299 // rest of Vorpaline expects to see tets on the
3300 // border).
3301 geo_debug_assert(!(keep_infinite_ && (t == NO_INDEX)));
3302 15304 cell_to_cell_store_[i] = t;
3303 });
3304 }
3305
3306 // In "keep_infinite" mode, we reorder the cells in such
3307 // a way that finite cells have indices [0..nb_finite_cells_-1]
3308 // and infinite cells have indices [nb_finite_cells_ .. nb_cells_-1]
3309
3310
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(keep_infinite_) {
3311 1 nb_finite_cells_ = 0;
3312 index_t finite_ptr = 0;
3313 1 index_t infinite_ptr = nb_tets - 1;
3314 for(;;) {
3315
2/2
✓ Branch 0 taken 479 times.
✓ Branch 1 taken 253 times.
 732 while(thread0->tet_is_finite(finite_ptr)) {
3316 479 old2new[finite_ptr] = finite_ptr;
3317 479 ++finite_ptr;
3318 479 ++nb_finite_cells_;
3319 }
3320
2/2
✓ Branch 0 taken 102 times.
✓ Branch 1 taken 253 times.
 355 while(!thread0->tet_is_finite(infinite_ptr)) {
3321 102 old2new[infinite_ptr] = infinite_ptr;
3322 102 --infinite_ptr;
3323 }
3324
2/2
✓ Branch 0 taken 252 times.
✓ Branch 1 taken 1 times.
 253 if(finite_ptr > infinite_ptr) {
3325 break;
3326 }
3327 252 old2new[finite_ptr] = infinite_ptr;
3328 252 old2new[infinite_ptr] = finite_ptr;
3329 252 ++nb_finite_cells_;
3330
2/2
✓ Branch 0 taken 1008 times.
✓ Branch 1 taken 252 times.
 1260 for(index_t lf=0; lf<4; ++lf) {
3331 1008 std::swap(
3332 1008 cell_to_cell_store_[4*finite_ptr + lf],
3333 1008 cell_to_cell_store_[4*infinite_ptr + lf]
3334 );
3335 }
3336
2/2
✓ Branch 0 taken 1008 times.
✓ Branch 1 taken 252 times.
 1260 for(index_t lv=0; lv<4; ++lv) {
3337 1008 std::swap(
3338 1008 cell_to_v_store_[4*finite_ptr + lv],
3339 1008 cell_to_v_store_[4*infinite_ptr + lv]
3340 );
3341 }
3342 252 ++finite_ptr;
3343 252 --infinite_ptr;
3344 252 }
3345
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 parallel_for(0, 4*nb_tets, [this, &old2new](index_t i) {
3346 4340 index_t t = cell_to_cell_store_[i];
3347 geo_debug_assert(t != NO_INDEX);
3348 4340 t = old2new[t];
3349 geo_debug_assert(t != NO_INDEX);
3350 4340 cell_to_cell_store_[i] = t;
3351 });
3352 }
3353
3354
3355
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(detailed_benchmark_mode_) {
3356
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 Logger::out("DelCompress")
3357
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << "max tets " << thread0->max_t()
3358 << std::endl;
3359
3360
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 Logger::out("DelCompress")
3361 << "Final number of tets " << nb_tets
3362 << std::endl;
3363
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(keep_infinite_) {
3364
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 Logger::out("DelCompress")
3365 << "Removed " << nb_tets_to_delete
3366 << " tets (free list)"
3367 << " : "
3368
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 << double(nb_tets_to_delete)*100.0/double(nb_tets) << "%"
3369 << std::endl;
3370 } else {
3371
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 Logger::out("DelCompress")
3372 << "Removed " << nb_tets_to_delete
3373 << " tets (free list and infinite)"
3374 << " : "
3375
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 << double(nb_tets_to_delete)*100.0/double(nb_tets) << "%"
3376 << std::endl;
3377 }
3378 }
3379
3380
2/2
✓ Branch 0 taken 3826 times.
✓ Branch 1 taken 2 times.
 3828 for(index_t t=0; t<nb_tets; ++t) {
3381 3826 cell_next_[t] = PeriodicDelaunay3dThread::NOT_IN_LIST;
3382 }
3383
3384 // Disconnect tets that were connected to infinite tets
3385
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(periodic_) {
3386
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 parallel_for(0, 4*nb_tets, [this,nb_tets](index_t i) {
3387
2/2
✓ Branch 0 taken 904 times.
✓ Branch 1 taken 10060 times.
 10964 if(cell_to_cell_store_[i] >= nb_tets) {
3388 904 cell_to_cell_store_[i] = NO_INDEX;
3389 }
3390 });
3391 #ifdef GEO_DEBUG
3392 for(index_t i=0; i<4*nb_tets; ++i) {
3393 geo_debug_assert(cell_to_v_store_[i] != NO_INDEX);
3394 }
3395 #endif
3396 }
3397 2 return nb_tets;
3398 2 }
3399
3400 ✗ index_t PeriodicDelaunay3d::nearest_vertex(const double* p) const {
3401 // TODO
3402 ✗ return Delaunay::nearest_vertex(p);
3403 }
3404
3405 ✗ void PeriodicDelaunay3d::set_BRIO_levels(const vector<index_t>& levels) {
3406 levels_ = levels;
3407 ✗ }
3408
3409 3 void PeriodicDelaunay3d::update_v_to_cell() {
3410
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3 geo_assert(!is_locked_); // Not thread-safe
3411 3 is_locked_ = true;
3412
3413 // Optimized version for large scale optimal transport,
3414 // can be removed (all cases treated)
3415
3/4
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 1 times.
 3 if(!update_periodic_v_to_cell_ && !keeps_infinite()) {
3416 2 v_to_cell_.assign(nb_vertices(), NO_INDEX);
3417
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 parallel_for(0, nb_cells(), [this](index_t c) {
3418
2/2
✓ Branch 0 taken 21008 times.
✓ Branch 1 taken 5252 times.
 26260 for(index_t lv = 0; lv < 4; lv++) {
3419
2/2
✓ Branch 0 taken 9762 times.
✓ Branch 1 taken 11246 times.
 21008 index_t v = cell_vertex(c, lv);
3420 // discriminates both vertex at infinity (NO_INDEX)
3421 // and VERTEX_OF_DELETED_TET (index_t(-2)).
3422
2/2
✓ Branch 0 taken 9762 times.
✓ Branch 1 taken 11246 times.
 21008 if(v < nb_vertices_non_periodic_) {
3423 9762 v_to_cell_[v] = c;
3424 }
3425 }
3426 5252 });
3427 2 is_locked_ = false; // Do not forget to unlock !
3428 2 return;
3429 }
3430
3431
3432 // Note: if keeps_infinite is set, then infinite vertex
3433 // tet chaining is at t2v_[nb_vertices].
3434
3435 // Create periodic_v_to_cell_ structure in compressed row
3436 // storage format.
3437
3438 // It was used in previous version for handling periodic boundary
3439 // conditions based on ConvexCell, it is no longer the case, new
3440 // code solely uses tetrahedra. It is kept here for reference for
3441 // implementing the distributed version (using ConvexCell can save
3442 // points tranfers).
3443
3444
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if(update_periodic_v_to_cell_) {
3445 ✗ periodic_v_to_cell_rowptr_.resize(nb_vertices_non_periodic_ + 1);
3446 ✗ periodic_v_to_cell_rowptr_[0] = 0;
3447 index_t cur = 0;
3448 ✗ for(index_t v=0; v<nb_vertices_non_periodic_; ++v) {
3449 ✗ cur += pop_count(vertex_instances_[v])-1;
3450 ✗ periodic_v_to_cell_rowptr_[v+1] = cur;
3451 }
3452 ✗ periodic_v_to_cell_data_.assign(cur, NO_INDEX);
3453 }
3454
3455
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if(keeps_infinite()) {
3456
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1 geo_assert(!periodic_);
3457 1 v_to_cell_.assign(nb_vertices()+1, NO_INDEX);
3458
2/2
✓ Branch 0 taken 1085 times.
✓ Branch 1 taken 1 times.
 1086 for(index_t c = 0; c < nb_cells(); c++) {
3459
2/2
✓ Branch 0 taken 4340 times.
✓ Branch 1 taken 1085 times.
 5425 for(index_t lv = 0; lv < 4; lv++) {
3460 index_t v = cell_vertex(c, lv);
3461
2/2
✓ Branch 0 taken 354 times.
✓ Branch 1 taken 3986 times.
 4340 if(v == NO_INDEX) {
3462 v = nb_vertices();
3463 }
3464 4340 v_to_cell_[v] = c;
3465 }
3466 }
3467 } else {
3468 ✗ v_to_cell_.assign(nb_vertices(), NO_INDEX);
3469 ✗ for(index_t c = 0; c < nb_cells(); c++) {
3470 ✗ for(index_t lv = 0; lv < 4; lv++) {
3471 index_t v = cell_vertex(c, lv);
3472 ✗ if(v < nb_vertices_non_periodic_) {
3473 ✗ v_to_cell_[v] = c;
3474 ✗ } else if(
3475 ✗ update_periodic_v_to_cell_ &&
3476 ✗ v != NO_INDEX &&
3477 v != PeriodicDelaunay3dThread::VERTEX_OF_DELETED_TET
3478 ) {
3479 index_t v_real = periodic_vertex_real(v);
3480 index_t v_instance = periodic_vertex_instance(v);
3481
3482 geo_debug_assert(
3483 (vertex_instances_[v_real] & (1u << v_instance)) != 0
3484 );
3485
3486 ✗ index_t slot = pop_count(
3487 ✗ vertex_instances_[v_real] & ((1u << v_instance)-1)
3488 ✗ ) - 1;
3489
3490 periodic_v_to_cell_data_[
3491 periodic_v_to_cell_rowptr_[v_real] + slot
3492 ✗ ] = c;
3493 }
3494 }
3495 }
3496 }
3497
3498 1 is_locked_ = false;
3499 }
3500
3501 ✗ void PeriodicDelaunay3d::update_cicl() {
3502 // Note: updates CICL information only for the
3503 // corners that correspond to real vertices
3504 // (anyway, in the API we will be always
3505 // starting from a real vertex !)
3506
3507 ✗ geo_assert(!is_locked_); // Not thread-safe
3508 ✗ is_locked_ = true;
3509 ✗ cicl_.resize(4 * nb_cells());
3510
3511 ✗ for(index_t v = 0; v < nb_vertices_non_periodic_; ++v) {
3512 ✗ index_t t = v_to_cell_[v];
3513 ✗ if(t != NO_INDEX) {
3514 ✗ index_t lv = index(t, v);
3515 set_next_around_vertex(t, lv, t);
3516 }
3517 }
3518
3519 ✗ if(keeps_infinite()) {
3520
3521 {
3522 // Process the infinite vertex at index nb_vertices().
3523 ✗ index_t t = v_to_cell_[nb_vertices()];
3524 ✗ if(t != NO_INDEX) {
3525 ✗ index_t lv = index(t, NO_INDEX);
3526 set_next_around_vertex(t, lv, t);
3527 }
3528 }
3529
3530 ✗ for(index_t t = 0; t < nb_cells(); ++t) {
3531 ✗ for(index_t lv = 0; lv < 4; ++lv) {
3532 index_t v = cell_vertex(t, lv);
3533 ✗ index_t vv = (v == NO_INDEX) ? nb_vertices() : v;
3534 ✗ if(v_to_cell_[vv] != t) {
3535 index_t t1 = v_to_cell_[vv];
3536 ✗ index_t lv1 = index(t1, v);
3537 index_t t2 = next_around_vertex(t1, lv1);
3538 set_next_around_vertex(t1, lv1, t);
3539 set_next_around_vertex(t, lv, t2);
3540 }
3541 }
3542 }
3543
3544
3545 } else {
3546 ✗ for(index_t t = 0; t < nb_cells(); ++t) {
3547 ✗ for(index_t lv = 0; lv < 4; ++lv) {
3548 index_t v = cell_vertex(t, lv);
3549 ✗ if(v < nb_vertices_non_periodic_ && v_to_cell_[v] != t) {
3550 index_t t1 = v_to_cell_[v];
3551 ✗ index_t lv1 = index(t1, v);
3552 index_t t2 = next_around_vertex(t1, lv1);
3553 set_next_around_vertex(t1, lv1, t);
3554 set_next_around_vertex(t, lv, t2);
3555 }
3556 }
3557 }
3558 }
3559
3560 ✗ is_locked_ = false;
3561 ✗ }
3562
3563 442 void PeriodicDelaunay3d::get_incident_tets(
3564 index_t v, IncidentTetrahedra& W
3565 ) const {
3566
3567 geo_debug_assert(
3568 periodic_ || v < nb_vertices_non_periodic_
3569 );
3570
3571 W.clear_incident_tets();
3572
3573 index_t t = NO_INDEX;
3574
1/2
✓ Branch 0 taken 442 times.
✗ Branch 1 not taken.
 442 if(v < nb_vertices_non_periodic_) {
3575 442 t = v_to_cell_[v];
3576 } else {
3577 index_t v_real = periodic_vertex_real(v);
3578 index_t v_instance = periodic_vertex_instance(v);
3579
3580 geo_debug_assert(
3581 (vertex_instances_[v_real] & (1u << v_instance))!=0
3582 );
3583
3584 ✗ index_t slot = pop_count(
3585 ✗ vertex_instances_[v_real] & ((1u << v_instance)-1)
3586 ✗ ) - 1;
3587
3588 ✗ t = periodic_v_to_cell_data_[
3589 periodic_v_to_cell_rowptr_[v_real] + slot
3590 ✗ ];
3591 }
3592
3593 // Can happen: empty power cell.
3594
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
 442 if(t == NO_INDEX) {
3595 ✗ return;
3596 }
3597
3598 // TODO: different version if CICL is stored ?
3599 {
3600
1/2
✓ Branch 0 taken 442 times.
✗ Branch 1 not taken.
 442 W.add_incident_tet(t);
3601 W.S.push(t);
3602
2/2
✓ Branch 0 taken 9762 times.
✓ Branch 1 taken 442 times.
 10204 while(!W.S.empty()) {
3603
1/2
✓ Branch 0 taken 9762 times.
✗ Branch 1 not taken.
 9762 t = W.S.top();
3604 W.S.pop();
3605 const index_t* T = &(cell_to_v_store_[4 * t]);
3606 index_t lv = PeriodicDelaunay3dThread::find_4(T,v);
3607
1/2
✓ Branch 0 taken 9762 times.
✗ Branch 1 not taken.
 9762 index_t neigh = cell_to_cell_store_[4*t + (lv + 1)%4];
3608
3/4
✓ Branch 0 taken 9762 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2584 times.
✓ Branch 3 taken 7178 times.
 19524 if(neigh != NO_INDEX && !W.has_incident_tet(neigh)) {
3609
1/2
✓ Branch 0 taken 2584 times.
✗ Branch 1 not taken.
 2584 W.add_incident_tet(neigh);
3610 W.S.push(neigh);
3611 }
3612
1/2
✓ Branch 0 taken 9762 times.
✗ Branch 1 not taken.
 9762 neigh = cell_to_cell_store_[4*t + (lv + 2)%4];
3613
3/4
✓ Branch 0 taken 9762 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3047 times.
✓ Branch 3 taken 6715 times.
 19524 if(neigh != NO_INDEX && !W.has_incident_tet(neigh)) {
3614
1/2
✓ Branch 0 taken 3047 times.
✗ Branch 1 not taken.
 3047 W.add_incident_tet(neigh);
3615 W.S.push(neigh);
3616 }
3617
1/2
✓ Branch 0 taken 9762 times.
✗ Branch 1 not taken.
 9762 neigh = cell_to_cell_store_[4*t + (lv + 3)%4];
3618
3/4
✓ Branch 0 taken 9762 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3689 times.
✓ Branch 3 taken 6073 times.
 19524 if(neigh != NO_INDEX && !W.has_incident_tet(neigh)) {
3619
1/2
✓ Branch 0 taken 3689 times.
✗ Branch 1 not taken.
 3689 W.add_incident_tet(neigh);
3620 W.S.push(neigh);
3621 }
3622 }
3623 }
3624 }
3625
3626 inline double dist2(const double* p, const double* q) {
3627 return
3628 geo_sqr(p[0]-q[0]) +
3629 geo_sqr(p[1]-q[1]) +
3630 geo_sqr(p[2]-q[2]) ;
3631 }
3632
3633 /**
3634 * \brief Copies a Laguerre cell from the triangulation.
3635 * \param[in] i the index of the vertex of which the Laguerre cell
3636 * should be computed.
3637 * \param[out] C the Laguerre cell.
3638 * \param[out] W the vector of neighbor vertices indices.
3639 */
3640
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 440 times.
 442 void PeriodicDelaunay3d::copy_Laguerre_cell_from_Delaunay(
3641 GEO::index_t i,
3642 ConvexCell& C,
3643 IncidentTetrahedra& W
3644 ) const {
3645 // Create global vertex indices if not present.
3646 C.create_vglobal();
3647 442 C.clear();
3648
3649 // Create the vertex at infinity.
3650 442 C.create_vertex(vec4(0.0, 0.0, 0.0, 0.0), NO_INDEX);
3651
3652
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
 442 GEO::vec3 Pi = vertex(i);
3653 double wi = weight(i);
3654
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
 442 double Pi_len2 = Pi[0]*Pi[0] + Pi[1]*Pi[1] + Pi[2]*Pi[2];
3655
3656
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
 442 if(stores_cicl()) {
3657 // Get neighbors and initialize cell from the tetrahedra in
3658 // 1-ring neighborhood of vertex.
3659 GEO::index_t t = GEO::index_t(vertex_cell(i));
3660 // Special case: Laguerre cell is empty (vertex has
3661 // no incident tet).
3662 ✗ if(t == NO_INDEX) {
3663 ✗ return;
3664 }
3665 do {
3666 ✗ GEO::index_t f = copy_Laguerre_cell_facet_from_Delaunay(
3667 i, Pi, wi, Pi_len2, t, C, W
3668 );
3669 t = GEO::index_t(next_around_vertex(t,f));
3670 ✗ } while(t != GEO::index_t(vertex_cell(i)));
3671 } else {
3672 442 get_incident_tets(i,W);
3673
2/2
✓ Branch 0 taken 9762 times.
✓ Branch 1 taken 442 times.
 10204 for(index_t t: W) {
3674 9762 copy_Laguerre_cell_facet_from_Delaunay(
3675 i, Pi, wi, Pi_len2, t, C, W
3676 );
3677 }
3678 }
3679 442 C.connect_triangles();
3680 }
3681
3682 9762 GEO::index_t PeriodicDelaunay3d::copy_Laguerre_cell_facet_from_Delaunay(
3683 GEO::index_t i,
3684 const GEO::vec3& Pi,
3685 double wi,
3686 double Pi_len2,
3687 GEO::index_t t,
3688 ConvexCell& C,
3689 IncidentTetrahedra& W
3690 ) const {
3691 geo_argused(W);
3692
3693 // Local tet vertex indices from facet
3694 // and vertex in facet indices.
3695 static GEO::index_t fv[4][3] = {
3696 {2,3,1},
3697 {3,2,0},
3698 {0,1,3},
3699 {2,1,0}
3700 };
3701
3702 9762 GEO::index_t f = index(t,GEO::index_t(i));
3703 GEO::index_t jkl[3]; // Global index (in Delaunay) of triangle vertices
3704 VBW::index_t l_jkl[3];// Local index (in C) of triangle vertices
3705
3706
3707 // Find or create the three vertices of the facet.
3708
2/2
✓ Branch 0 taken 29286 times.
✓ Branch 1 taken 9762 times.
 39048 for(int lfv=0; lfv<3; ++lfv) {
3709
3710 29286 jkl[lfv] = GEO::index_t(cell_vertex(t, fv[f][lfv]));
3711 29286 l_jkl[lfv] = VBW::index_t(-1);
3712
3713 // Vertex already created in C (note:
3714 // also works for vertex at infinity)
3715
2/2
✓ Branch 0 taken 350285 times.
✓ Branch 1 taken 5586 times.
 355871 for(VBW::index_t u=0; u<C.nb_v(); ++u) {
3716
2/2
✓ Branch 0 taken 23700 times.
✓ Branch 1 taken 326585 times.
 350285 if(C.v_global_index(u) == jkl[lfv]) {
3717 23700 l_jkl[lfv] = VBW::index_t(u);
3718 23700 break;
3719 }
3720 }
3721
3722 // vertex not found, create vertex in C
3723
2/2
✓ Branch 0 taken 5586 times.
✓ Branch 1 taken 23700 times.
 29286 if(l_jkl[lfv] == VBW::index_t(-1)) {
3724 5586 l_jkl[lfv] = C.nb_v();
3725
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5586 times.
 5586 vec3 Pj = vertex(jkl[lfv]);
3726 double Pj_len2 = length2(Pj);
3727 double wj = weight(jkl[lfv]);
3728 5586 double a = 2.0 * (Pi[0] - Pj[0]);
3729 5586 double b = 2.0 * (Pi[1] - Pj[1]);
3730 5586 double c = 2.0 * (Pi[2] - Pj[2]);
3731 5586 double d = ((wi - Pi_len2) - (wj - Pj_len2));
3732 5586 C.create_vertex(vec4(a,b,c,d), jkl[lfv]);
3733 }
3734 }
3735
3736 9762 C.create_triangle(l_jkl[0], l_jkl[1], l_jkl[2]);
3737
3738 9762 return f;
3739 }
3740
3741 /*************************************************************************/
3742
3743 2 void PeriodicDelaunay3d::insert_vertices(
3744 const char* phase, index_t b, index_t e
3745 ) {
3746
3747
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 Stopwatch W(phase,detailed_benchmark_mode_);
3748
3749
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(detailed_benchmark_mode_) {
3750
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 6 Logger::out(phase) << "Inserting " << (e-b)
3751 << " additional vertices" << std::endl;
3752 }
3753
3754
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 has_empty_cells_ = false;
3755
3756
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 nb_vertices_ = reorder_.size();
3757 vector<index_t> levels;
3758
3759
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 compute_BRIO_order_periodic(
3760
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 nb_vertices_non_periodic_ * 27, // nb of possible periodic vertices
3761 vertex_ptr(0),
3762 3, dimension(),
3763
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 reorder_,
3764 reorder_.begin() + long(b),
3765 reorder_.begin() + long(e),
3766
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 period_,
3767 64, 0.125,
3768 &levels
3769 );
3770
3771
3772 #ifdef GEO_DEBUG
3773 // Check that the same vertex was not inserted twice
3774 for(index_t i=b; i+1<e; ++i) {
3775 geo_debug_assert(reorder_[i] != reorder_[i+1]);
3776 }
3777 #endif
3778
3779
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 insert_vertices_with_BRIO(phase, levels);
3780
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if(abort_on_empty_cell_ && has_empty_cells_) {
3781 return;
3782 }
3783 PeriodicDelaunay3dThread* thread0 = thread(0);
3784 2 nb_vertices_ = reorder_.size();
3785
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 set_arrays(
3786 thread0->max_t(),
3787 cell_to_v_store_.data(),
3788 cell_to_cell_store_.data()
3789 );
3790
3791
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(!strcmp(phase, "insert-I")) {
3792
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 stats_.phase_I_insert_t_ = W.elapsed_time();
3793 1 stats_.phase_I_insert_nb_ = e-b;
3794
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 } else if(!strcmp(phase, "insert-II")) {
3795
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 stats_.phase_II_insert_t_ = W.elapsed_time();
3796 1 stats_.phase_II_insert_nb_ = e-b;
3797 }
3798 2 }
3799
3800 4 void PeriodicDelaunay3d::insert_vertices_with_BRIO(
3801 const char* phase, const vector<index_t>& levels
3802 ) {
3803
3804
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 4 times.
 24 for(index_t t=0; t<threads_.size(); ++t) {
3805 thread(t)->reset_stats();
3806 }
3807
3808 PeriodicDelaunay3dThread* thread0 = thread(0);
3809
3810 index_t lvl = 1;
3811
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 13 while(lvl < (levels.size() - 1) && (levels[lvl] - levels[0]) < 1000) {
3812 5 ++lvl;
3813 }
3814
3815
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if(detailed_benchmark_mode_) {
3816
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 Logger::out(phase)
3817
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 << "Using " << levels.size()-1 << " levels" << std::endl;
3818
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 8 Logger::out(phase)
3819
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 << "Levels 0 - " << lvl-1
3820 << ": bootstraping with first levels in sequential mode"
3821 << std::endl;
3822 }
3823
3824 4 thread0->set_work(levels[0], levels[lvl]);
3825 4 thread0->run();
3826
3827
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if(thread0->has_empty_cells()) {
3828 ✗ has_empty_cells_ = true;
3829 ✗ if(abort_on_empty_cell_) {
3830 return;
3831 }
3832 }
3833
3834 index_t nb_sequential_points = 0;
3835 index_t first_lvl = lvl;
3836
3837 // Insert points in all BRIO levels
3838
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 for(; lvl<levels.size()-1; ++lvl) {
3839
3840 ✗ index_t lvl_b = levels[lvl];
3841 ✗ index_t lvl_e = levels[lvl+1];
3842
3843 ✗ if(detailed_benchmark_mode_) {
3844 ✗ Logger::out(phase) << "Level "
3845 << lvl << " : start "
3846 << " nbv = "
3847 ✗ << (lvl_e - lvl_b)
3848 << std::endl;
3849 }
3850
3851 ✗ index_t work_size = (lvl_e - lvl_b)/index_t(threads_.size());
3852
3853 // Initialize threads
3854 index_t b = lvl_b;
3855 ✗ for(index_t t=0; t<threads_.size(); ++t) {
3856 ✗ index_t e = (t == threads_.size()-1) ? lvl_e : b+work_size;
3857
3858 // Copy the indices of the first created tetrahedron
3859 // and the maximum valid tetrahedron index max_t_
3860 ✗ if(lvl == first_lvl && t!=0) {
3861 thread(t)->set_max_t(thread0->max_t());
3862 }
3863 thread(t)->set_work(b,e);
3864 b = e;
3865 }
3866
3867 ✗ check_max_t();
3868 ✗ Process::run_threads(threads_);
3869
3870 ✗ for(index_t t=0; t<this->nb_threads(); ++t) {
3871 ✗ if(thread(t)->has_empty_cells()) {
3872 ✗ has_empty_cells_ = true;
3873 ✗ if(abort_on_empty_cell_) {
3874 return;
3875 }
3876 }
3877 }
3878
3879 // Run threads sequentialy, to insert missing points if
3880 // memory overflow was encountered (in sequential mode,
3881 // dynamic memory growing works)
3882
3883 index_t this_level_nb_sequential_points = 0;
3884
3885 ✗ for(index_t t=0; t<threads_.size(); ++t) {
3886 PeriodicDelaunay3dThread* t1 = thread(t);
3887 ✗ this_level_nb_sequential_points += t1->work_size();
3888 ✗ if(t != 0) {
3889 // We need to copy max_t_ from previous thread,
3890 // since the memory pool may have grown.
3891 ✗ PeriodicDelaunay3dThread* t2 = thread(t-1);
3892 t1->set_max_t(t2->max_t());
3893 }
3894 ✗ t1->run();
3895 ✗ if(t1->has_empty_cells()) {
3896 ✗ has_empty_cells_ = true;
3897 ✗ if(abort_on_empty_cell_) {
3898 return;
3899 }
3900 }
3901 }
3902
3903 // If some tetrahedra were created in sequential mode, then
3904 // the maximum valid tetrahedron index was increased by all
3905 // the threads in increasing number, so we copy it from the
3906 // last thread into thread0 since we use thread0 afterwards
3907 // to get max_t()
3908
3909 ✗ if(this_level_nb_sequential_points != 0) {
3910 PeriodicDelaunay3dThread* t0 = thread(0);
3911 PeriodicDelaunay3dThread* tn = thread(
3912 this->nb_threads()-1
3913 );
3914 t0->set_max_t(tn->max_t());
3915 }
3916
3917 ✗ if(abort_on_empty_cell_ && has_empty_cells_) {
3918 return;
3919 }
3920
3921 ✗ nb_sequential_points += this_level_nb_sequential_points;
3922 }
3923
3924
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if(detailed_benchmark_mode_) {
3925 index_t tot_rollbacks = 0 ;
3926 index_t tot_failed_locate = 0 ;
3927
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 4 times.
 20 for(index_t t=0; t<threads_.size(); ++t) {
3928
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 32 Logger::out(phase)
3929 << "thread " << std::setw(3) << t << " : "
3930 << std::setw(3) << thread(t)->nb_rollbacks()
3931 << " rollbacks "
3932 << std::setw(3) << thread(t)->nb_failed_locate()
3933 << " restarted locate"
3934 << std::endl;
3935 16 tot_rollbacks += thread(t)->nb_rollbacks();
3936 16 tot_failed_locate += thread(t)->nb_failed_locate();
3937 }
3938
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 Logger::out(phase) << "------------------" << std::endl;
3939
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 Logger::out(phase) << "total: "
3940 << tot_rollbacks << " rollbacks "
3941 << tot_failed_locate << " restarted locate"
3942 << std::endl;
3943
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if(nb_sequential_points == 0) {
3944
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 8 Logger::out(phase) << "All points where inserted in parallel"
3945 << std::endl;
3946 } else {
3947 ✗ Logger::out(phase) << nb_sequential_points
3948 << " points inserted in sequential mode"
3949 << std::endl;
3950 }
3951 }
3952
3953 4 nb_vertices_ = reorder_.size();
3954 index_t nb_tets = thread0->max_t();
3955
3956
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 4 times.
 40 for(index_t i=0; i<nb_threads(); ++i) {
3957
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 16 geo_assert(thread(i)->max_t() <= nb_tets);
3958 }
3959 }
3960
3961 6510 bool PeriodicDelaunay3d::Laguerre_vertex_is_in_conflict_with_plane(
3962 index_t t, vec4 P
3963 ) const {
3964
3965 // Note: facet orientations and signs follow:
3966 // - ConvexCell
3967 // - copy_Laguerre_facet_from()
3968
3969 // Local tet vertex indices from facet and vertex in facet indices.
3970 // Carefully chosen in such a way that f[(lv+1)%4][2] == lv
3971 // This is used in the code block below, that handles tets with
3972 // a vertex at infinity
3973 static GEO::index_t fv[4][3] = {
3974 {1,2,3},
3975 {3,2,0},
3976 {3,0,1},
3977 {1,0,2}
3978 };
3979
3980 // Particular case, vertex at infinity
3981
2/2
✓ Branch 0 taken 23862 times.
✓ Branch 1 taken 4386 times.
 28248 for(index_t lv=0; lv<4; ++lv) {
3982
2/2
✓ Branch 0 taken 2124 times.
✓ Branch 1 taken 21738 times.
 23862 if(cell_vertex(t,lv) == NO_INDEX) {
3983 2124 index_t li = (lv + 1) % 4; // this vertex is not at infty
3984 // li is also the index of a facet with the vertex at infty
3985 // as the last vertex (fv[][] was constructed so), so we
3986 // are sure that vi, vj, vk are the not at infty
3987 index_t vi = cell_vertex(t, li);
3988 2124 index_t vj = cell_vertex(t, fv[li][0]);
3989 2124 index_t vk = cell_vertex(t, fv[li][1]);
3990 geo_debug_assert(fv[li][2] == lv); // we know that l == -1
3991 2124 vec3 pi = vertex(vi);
3992 2124 vec3 pj = vertex(vj);
3993 2124 vec3 pk = vertex(vk);
3994 Sign s = PCK::det_3d(
3995 2124 pi-pj,
3996 2124 pi-pk,
3997 2124 vec3(P.x,P.y,P.z)
3998 );
3999 2124 return (s <= 0);
4000 }
4001 }
4002
4003 index_t v1 = cell_vertex(t,0);
4004 index_t v2 = cell_vertex(t,1);
4005 index_t v3 = cell_vertex(t,2);
4006 index_t v4 = cell_vertex(t,3);
4007
4008 4386 vec3 p1 = vertex(v1);
4009 4386 vec3 p2 = vertex(v2);
4010 4386 vec3 p3 = vertex(v3);
4011
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4386 times.
 4386 vec3 p4 = vertex(v4);
4012
4013
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4386 times.
 4386 double l1 = weight(v1) - length2(p1);
4014
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4386 times.
 4386 double l2 = weight(v2) - length2(p2);
4015
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4386 times.
 4386 double l3 = weight(v3) - length2(p3);
4016 4386 double l4 = weight(v4) - length2(p4);
4017
4018 Sign s = PCK::det_4d(
4019 8772 vec4(2.0*(p1.x-p2.x), 2.0*(p1.y-p2.y), 2.0*(p1.z-p2.z), l1-l2),
4020 4386 vec4(2.0*(p1.x-p3.x), 2.0*(p1.y-p3.y), 2.0*(p1.z-p3.z), l1-l3),
4021 4386 vec4(2.0*(p1.x-p4.x), 2.0*(p1.y-p4.y), 2.0*(p1.z-p4.z), l1-l4),
4022 P
4023 );
4024
4025 4386 return (s >= 0);
4026 }
4027
4028
4029 1 void PeriodicDelaunay3d::handle_periodic_boundaries_phase_I() {
4030
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 Stopwatch W_classify_I("classify-I", detailed_benchmark_mode_);
4031 1 PeriodicDelaunay3dThread* thread0 = thread(0);
4032
4033 // Lag_cell_status:
4034 //
4035 // each bit k in 0..5 indicate that cell has at least one vertex
4036 // in conflict with plane Pk (outside central cube)
4037 // each bit k in 6..11 indicate that cell has all its vertices
4038 // in conflict with plane Pk (outside central cube)
4039 // status = 0 -> cell is contained by central cube
4040 // (status & conflict_mask) != 0 -> cell straddles bndr of central cube
4041 // (status & all_conflict_mask) != 0 -> cell is outside central cube
4042
4043 // static Numeric::uint16 conflict_mask = Numeric::uint16(63u);
4044 static Numeric::uint16 all_conflict_mask = Numeric::uint16(63u << 6);
4045
4046 std::atomic<Numeric::uint16>* Lag_cell_status
4047
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 = new std::atomic<Numeric::uint16>[nb_vertices_non_periodic_];
4048
4049
2/2
✓ Branch 0 taken 221 times.
✓ Branch 1 taken 1 times.
 222 for(index_t i=0; i<nb_vertices_non_periodic_; ++i) {
4050 221 Lag_cell_status[i] = all_conflict_mask;
4051 }
4052
4053 {
4054 vec4 cube_face[6] = {
4055 vec4( 1.0, 0.0, 0.0, 0.0),
4056 vec4(-1.0, 0.0, 0.0, period_.x),
4057 vec4( 0.0, 1.0, 0.0, 0.0),
4058 vec4( 0.0,-1.0, 0.0, period_.y),
4059 vec4( 0.0, 0.0, 1.0, 0.0),
4060 vec4( 0.0, 0.0,-1.0, period_.z),
4061 };
4062
4063
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 parallel_for(0, thread0->max_t(), [&](index_t t) {
4064
2/2
✓ Branch 0 taken 1085 times.
✓ Branch 1 taken 1426 times.
 2511 if(thread0->tet_is_free(t)) {
4065 return;
4066 }
4067
2/2
✓ Branch 0 taken 6510 times.
✓ Branch 1 taken 1085 times.
 7595 for(index_t k=0; k<6; ++k) {
4068 6510 bool conflict = Laguerre_vertex_is_in_conflict_with_plane(
4069 6510 t, cube_face[k]
4070 );
4071
2/2
✓ Branch 0 taken 26040 times.
✓ Branch 1 taken 6510 times.
 32550 for(index_t lv=0; lv<4; ++lv) {
4072
2/2
✓ Branch 0 taken 2124 times.
✓ Branch 1 taken 23916 times.
 26040 index_t v = cell_vertex(t,lv);
4073
2/2
✓ Branch 0 taken 2124 times.
✓ Branch 1 taken 23916 times.
 26040 if(v == NO_INDEX) {
4074 2124 continue;
4075 }
4076
2/2
✓ Branch 0 taken 5660 times.
✓ Branch 1 taken 18256 times.
 23916 if(conflict) {
4077 // set 'conflict' bit k
4078 5660 Lag_cell_status[v].fetch_or(
4079 5660 Numeric::uint16(1u << k),
4080 std::memory_order_relaxed
4081 );
4082 } else {
4083 // reset 'all conflict' bit k
4084 18256 Lag_cell_status[v].fetch_and(
4085 18256 Numeric::uint16(~(1u << (k+6))),
4086 std::memory_order_relaxed
4087 );
4088 }
4089 }
4090 }
4091 }
4092 );
4093 }
4094
4095 // Count cells inside, crossing, outside
4096 {
4097 1 stats_.phase_I_nb_inside_ = 0;
4098 1 stats_.phase_I_nb_cross_ = 0;
4099 1 stats_.phase_I_nb_outside_ = 0;
4100
2/2
✓ Branch 0 taken 221 times.
✓ Branch 1 taken 1 times.
 222 for(index_t v=0; v<nb_vertices_non_periodic_; ++v) {
4101
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 203 times.
 221 Numeric::uint16 status = Lag_cell_status[v];
4102
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 203 times.
 221 if(status == 0) {
4103 18 ++stats_.phase_I_nb_inside_;
4104
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 203 times.
 203 } else if((status & all_conflict_mask) != 0) {
4105 ✗ ++stats_.phase_I_nb_outside_;
4106 } else {
4107 203 ++stats_.phase_I_nb_cross_;
4108 }
4109 }
4110
4111
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if(detailed_benchmark_mode_) {
4112
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 Logger::out("classify-I") << "Nb cells inside cube: "
4113
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 << stats_.phase_I_nb_inside_
4114 << std::endl;
4115
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 Logger::out("classify-I") << "Nb cells on boundary: "
4116
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 << stats_.phase_I_nb_cross_
4117 << std::endl;
4118
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 Logger::out("classify-I") << "Nb cells outside cube: "
4119
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 << stats_.phase_I_nb_outside_
4120 << std::endl;
4121 }
4122 }
4123
4124 // Indicates for each real vertex the instances it has.
4125 // Each bit of vertex_instances_[v] indicates which instance
4126 // is used.
4127
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vertex_instances_.assign(nb_vertices_non_periodic_,1);
4128
4129
2/2
✓ Branch 0 taken 221 times.
✓ Branch 1 taken 1 times.
 222 for(index_t v=0; v<nb_vertices_non_periodic_; ++v) {
4130
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 203 times.
 221 Numeric::uint16 status = Lag_cell_status[v];
4131
4132 bool status_inside = (status == 0);
4133
4134 // In the distributed version, we might need to distinguish
4135 // also the following two cases:
4136 // bool status_outside = ((status & all_conflict_mask) != 0);
4137 // bool status_crossing = !status_inside && !status_outside;
4138
4139 // if cell is inside cube, no instance to create
4140
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 203 times.
 221 if(status_inside) {
4141 18 continue;
4142 }
4143
4144 // Integer translations associated with the six plane equations
4145 static int T[6][3]= {
4146 {-1, 0, 0},
4147 { 1, 0, 0},
4148 { 0,-1, 0},
4149 { 0, 1, 0},
4150 { 0, 0,-1},
4151 { 0, 0, 1}
4152 };
4153
4154 // Detect the bounds of the sub-(rubic's) cube overlapped
4155 // by the cell.
4156
4157 203 int TXmin = 2, TXmax = -2,
4158 203 TYmin = 2, TYmax = -2,
4159 203 TZmin = 2, TZmax = -2;
4160
4161
2/2
✓ Branch 0 taken 1218 times.
✓ Branch 1 taken 203 times.
 1421 FOR(i,6) {
4162
2/2
✓ Branch 0 taken 964 times.
✓ Branch 1 taken 254 times.
 1218 if((status & Numeric::uint8(1u << i)) != 0) {
4163 964 TXmin = std::min(TXmin, T[i][0]);
4164 964 TXmax = std::max(TXmax, T[i][0]);
4165 964 TYmin = std::min(TYmin, T[i][1]);
4166 964 TYmax = std::max(TYmax, T[i][1]);
4167 964 TZmin = std::min(TZmin, T[i][2]);
4168 964 TZmax = std::max(TZmax, T[i][2]);
4169 }
4170 }
4171
4172
2/2
✓ Branch 0 taken 525 times.
✓ Branch 1 taken 203 times.
 728 for(int TX = TXmin; TX <= TXmax; ++TX) {
4173
2/2
✓ Branch 0 taken 1229 times.
✓ Branch 1 taken 525 times.
 1754 for(int TY = TYmin; TY <= TYmax; ++TY) {
4174
2/2
✓ Branch 0 taken 3621 times.
✓ Branch 1 taken 1229 times.
 4850 for(int TZ = TZmin; TZ <= TZmax; ++TZ) {
4175 index_t instance = T_to_instance(-TX,-TY,-TZ);
4176 // Skip instance 0 (it was already inserted !)
4177
2/2
✓ Branch 0 taken 3437 times.
✓ Branch 1 taken 184 times.
 3621 if(instance != 0) {
4178 3437 vertex_instances_[v] |= (1u << instance);
4179
1/2
✓ Branch 1 taken 3437 times.
✗ Branch 2 not taken.
 3437 reorder_.push_back(make_periodic_vertex(v,instance));
4180 }
4181 }
4182 }
4183 }
4184 }
4185
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 delete[] Lag_cell_status;
4186
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 stats_.phase_I_classify_t_ = W_classify_I.elapsed_time();
4187 1 }
4188
4189 1 void PeriodicDelaunay3d::handle_periodic_boundaries_phase_II() {
4190
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 Stopwatch W_classify_II("classify-II", detailed_benchmark_mode_);
4191 1 PeriodicDelaunay3dThread* thread0 = thread(0);
4192
4193 // Computes translation_table[][]
4194 // translation_table[instance2][instance1] transforms
4195 // instance2 into the frame of instance1
4196 Numeric::int8 translation_table[27][27];
4197
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 1 times.
 28 for(index_t instance1=0; instance1<27; ++instance1) {
4198 27 int Tx1 = translation[instance1][0];
4199 27 int Ty1 = translation[instance1][1];
4200 27 int Tz1 = translation[instance1][2];
4201
2/2
✓ Branch 0 taken 729 times.
✓ Branch 1 taken 27 times.
 756 for(index_t instance2=0; instance2<27; ++instance2) {
4202 729 int Tx2 = translation[instance2][0];
4203 729 int Ty2 = translation[instance2][1];
4204 729 int Tz2 = translation[instance2][2];
4205
4206 729 translation_table[instance2][instance1] =
4207 Numeric::int8(instance2);
4208
4209
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 702 times.
 729 if(instance1 == instance2) {
4210 27 continue;
4211 }
4212
4213 1088 if(
4214
2/2
✓ Branch 0 taken 540 times.
✓ Branch 1 taken 162 times.
 702 std::abs(Tx2 - Tx1) >= 2 ||
4215
2/2
✓ Branch 0 taken 414 times.
✓ Branch 1 taken 126 times.
 540 std::abs(Ty2 - Ty1) >= 2 ||
4216
2/2
✓ Branch 0 taken 98 times.
✓ Branch 1 taken 316 times.
 414 std::abs(Tz2 - Tz1) >= 2
4217 ) {
4218 // Here we could use -1 to encode large displacements,
4219 // and issue an error message (for now they are ignored)
4220 386 continue;
4221 }
4222
4223 316 translation_table[instance2][instance1] =
4224 Numeric::int8(T_to_instance(Tx2-Tx1,Ty2-Ty1,Tz2-Tz1));
4225 }
4226 }
4227
4228 std::atomic<Numeric::uint32>* new_vertex_instances =
4229
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 new std::atomic<Numeric::uint32>[vertex_instances_.size()];
4230
4231
2/2
✓ Branch 0 taken 221 times.
✓ Branch 1 taken 1 times.
 223 for(index_t i=0; i<vertex_instances_.size(); ++i) {
4232 221 new_vertex_instances[i] = vertex_instances_[i];
4233 }
4234
4235
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Process::spinlock perio_lock = GEOGRAM_SPINLOCK_INIT;
4236
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 parallel_for(0, thread0->max_t(), [&,this](index_t t) {
4237
4/4
✓ Branch 0 taken 14247 times.
✓ Branch 1 taken 1414 times.
✓ Branch 2 taken 13415 times.
✓ Branch 3 taken 832 times.
 15661 if(!thread0->tet_is_real(t)) {
4238 return;
4239 }
4240 // Find the edges v1,v2 such that:
4241 // v1 is a vertex that was inserted in phase-I (instance != 0)
4242 // v2 is a vertex in an instance different from v1
4243
2/2
✓ Branch 0 taken 53660 times.
✓ Branch 1 taken 13415 times.
 67075 for(index_t lv=0; lv<4; ++lv) {
4244
2/2
✓ Branch 0 taken 5614 times.
✓ Branch 1 taken 48046 times.
 53660 index_t v1 = thread0->finite_tet_vertex(t, lv);
4245
2/2
✓ Branch 0 taken 5614 times.
✓ Branch 1 taken 48046 times.
 53660 index_t v1_instance = periodic_vertex_instance(v1);
4246
2/2
✓ Branch 0 taken 5614 times.
✓ Branch 1 taken 48046 times.
 53660 if(v1_instance == 0) {
4247 5614 continue;
4248 }
4249
2/2
✓ Branch 0 taken 144138 times.
✓ Branch 1 taken 48046 times.
 192184 for(index_t dlv=1; dlv<4; ++dlv) {
4250
2/2
✓ Branch 0 taken 58500 times.
✓ Branch 1 taken 85638 times.
 144138 index_t v2 = thread0->finite_tet_vertex(t, (lv + dlv)%4);
4251
2/2
✓ Branch 0 taken 58500 times.
✓ Branch 1 taken 85638 times.
 144138 index_t v2_real = periodic_vertex_real(v2);
4252 index_t v2_instance = periodic_vertex_instance(v2);
4253
4254 // transform v2_instance into the local frame of v1
4255 144138 v2_instance = index_t(
4256 144138 translation_table[v2_instance][v1_instance]
4257 );
4258
2/2
✓ Branch 0 taken 58500 times.
✓ Branch 1 taken 85638 times.
 144138 if(v2_instance == v1_instance) {
4259 58500 continue;
4260 }
4261
4262 // create the transformed v2_instance if it does not
4263 // already exist, and memorize it in the new list of
4264 // vertices to create
4265
4266 85638 Numeric::uint32 mask = (1u << v2_instance);
4267 Numeric::uint32 prev_instances =
4268
2/2
✓ Branch 0 taken 140 times.
✓ Branch 1 taken 85498 times.
 85638 new_vertex_instances[v2_real].fetch_or(
4269 mask, std::memory_order_relaxed // only need atomic
4270 );
4271
4272
2/2
✓ Branch 0 taken 140 times.
✓ Branch 1 taken 85498 times.
 85638 if((prev_instances & mask) == 0) {
4273 140 Process::acquire_spinlock(perio_lock);
4274 140 reorder_.push_back(
4275 140 make_periodic_vertex(v2_real, v2_instance)
4276 );
4277 140 Process::release_spinlock(perio_lock);
4278 }
4279 }
4280 }
4281 });
4282
2/2
✓ Branch 0 taken 221 times.
✓ Branch 1 taken 1 times.
 223 for(index_t i=0; i<vertex_instances_.size(); ++i) {
4283 221 vertex_instances_[i] = new_vertex_instances[i];
4284 }
4285
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 delete[] new_vertex_instances;
4286
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 stats_.phase_II_classify_t_ = W_classify_II.elapsed_time();
4287 1 }
4288
4289 1 void PeriodicDelaunay3d::handle_periodic_boundaries() {
4290
4291 // Update pointers so that queries function will work (even in our
4292 // "transient state").
4293 PeriodicDelaunay3dThread* thread0 = thread(0);
4294
4295 1 set_arrays(
4296 thread0->max_t(),
4297 cell_to_v_store_.data(),
4298 cell_to_cell_store_.data()
4299 );
4300
4301 // Test for empty cells
4302 // TODO: I tested, it really occurs that there is still an empty
4303 // cell here, why is it not detected before ? To be understood.
4304 1 update_v_to_cell();
4305
2/2
✓ Branch 0 taken 221 times.
✓ Branch 1 taken 1 times.
 222 for(index_t v=0; v<nb_vertices_non_periodic_; ++v) {
4306
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 221 times.
 221 if(v_to_cell_[v] == NO_INDEX) {
4307 ✗ has_empty_cells_ = true;
4308 ✗ if(abort_on_empty_cell_) {
4309 return;
4310 }
4311 }
4312 }
4313
4314 // Phase I: find the cells that intersect the boundaries, and
4315 // create periodic vertices for each possible translation.
4316 {
4317
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Stopwatch W_phase_I("phase-I",false);
4318
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 handle_periodic_boundaries_phase_I();
4319
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 insert_vertices(
4320 "insert-I", nb_vertices_non_periodic_, reorder_.size()
4321 );
4322
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 stats_.phase_I_t_ = W_phase_I.elapsed_time();
4323 1 }
4324
4325 // Phase II: Insert the real neighbors of the virtual vertices,
4326 // back-translated to the original position.
4327 {
4328 index_t nb_vertices_phase_I = reorder_.size();
4329
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Stopwatch W_phase_II("phase-II",false);
4330
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 handle_periodic_boundaries_phase_II();
4331
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 insert_vertices("insert-II", nb_vertices_phase_I, reorder_.size());
4332
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 stats_.phase_II_t_ = W_phase_II.elapsed_time();
4333 1 }
4334 }
4335
4336 ✗ void PeriodicDelaunay3d::check_volume() {
4337 ✗ ConvexCell C;
4338 ✗ C.use_exact_predicates(convex_cell_exact_predicates_);
4339
4340 ✗ Logger::out("Periodic") << "Checking total volume..." << std::endl;
4341 double sumV = 0;
4342 ✗ IncidentTetrahedra W;
4343
4344 ✗ FOR(v, nb_vertices_non_periodic_) {
4345 ✗ copy_Laguerre_cell_from_Delaunay(v, C, W);
4346 #ifdef GEO_DEBUG
4347 for(
4348 VBW::ushort t = C.first_triangle();
4349 t!=VBW::END_OF_LIST; t=C.next_triangle(t)
4350 ) {
4351 for(index_t lv=0; lv<3; ++lv) {
4352 // Make sure there is no vertex at infinity
4353 geo_debug_assert(
4354 C.triangle_v_local_index(t,VBW::index_t(lv)) != 0
4355 );
4356 }
4357 }
4358 #endif
4359 ✗ C.compute_geometry();
4360 ✗ sumV += C.volume();
4361 }
4362
4363 ✗ double expectedV = period_.x*period_.y*period_.z;
4364
4365 ✗ Logger::out("Periodic") << "Sum volumes = " << sumV << std::endl;
4366 ✗ Logger::out("Periodic") << " (expected " << expectedV << ")"
4367 << std::endl;
4368
4369 ✗ if(::fabs(sumV - expectedV) / expectedV >= 1.0 / 10000.0) {
4370 ✗ Logger::err("Periodic") << "FATAL, volume error is too large"
4371 << std::endl;
4372 ✗ exit(-1);
4373 }
4374
4375 ✗ }
4376
4377 2 void PeriodicDelaunay3d::save_cells(
4378 const std::string& basename, bool clipped
4379 ) {
4380 static int index = 1;
4381
4382 2 std::string index_string = String::to_string(index);
4383
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 while(index_string.length() < 3) {
4384
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 index_string = "0" + index_string;
4385 }
4386
4387
1/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2 std::ofstream out((basename+index_string+".obj").c_str());
4388 index_t v_off = 1;
4389
4390
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 ConvexCell C;
4391
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 C.use_exact_predicates(convex_cell_exact_predicates_);
4392
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 IncidentTetrahedra W;
4393
2/2
✓ Branch 0 taken 442 times.
✓ Branch 1 taken 2 times.
 444 for(index_t vv=0; vv<nb_vertices_non_periodic_; ++vv) {
4394
1/2
✓ Branch 1 taken 442 times.
✗ Branch 2 not taken.
 442 copy_Laguerre_cell_from_Delaunay(vv, C, W);
4395
2/2
✓ Branch 0 taken 221 times.
✓ Branch 1 taken 221 times.
 442 if(clipped) {
4396
2/4
✓ Branch 1 taken 221 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 221 times.
✗ Branch 5 not taken.
 221 C.clip_by_plane(vec4( 1.0, 0.0, 0.0, 0.0));
4397
2/4
✓ Branch 1 taken 221 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 221 times.
✗ Branch 5 not taken.
 221 C.clip_by_plane(vec4(-1.0, 0.0, 0.0, period_.x));
4398
2/4
✓ Branch 1 taken 221 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 221 times.
✗ Branch 5 not taken.
 221 C.clip_by_plane(vec4( 0.0, 1.0, 0.0, 0.0));
4399
2/4
✓ Branch 1 taken 221 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 221 times.
✗ Branch 5 not taken.
 221 C.clip_by_plane(vec4( 0.0,-1.0, 0.0, period_.y));
4400
2/4
✓ Branch 1 taken 221 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 221 times.
✗ Branch 5 not taken.
 221 C.clip_by_plane(vec4( 0.0, 0.0, 1.0, 0.0));
4401
1/2
✓ Branch 1 taken 221 times.
✗ Branch 2 not taken.
 221 C.clip_by_plane(vec4( 0.0, 0.0,-1.0, period_.z));
4402 }
4403
1/2
✓ Branch 1 taken 442 times.
✗ Branch 2 not taken.
 442 v_off += C.save(out, v_off, 0.1);
4404 }
4405 2 ++index;
4406 4 }
4407
4408 ✗ void PeriodicDelaunay3d::check_max_t() {
4409 ✗ index_t max_t = 0;
4410 ✗ for(index_t i=0; i<nb_threads(); ++i) {
4411 ✗ max_t = std::max(max_t, thread(i)->max_t());
4412 }
4413 ✗ for(index_t i=0; i<nb_threads(); ++i) {
4414 ✗ geo_assert(thread(i)->max_t() == max_t);
4415 }
4416 ✗ }
4417
4418 /***********************************************************/
4419
4420 2 PeriodicDelaunay3d::Stats::Stats() {
4421 2 reset();
4422 2 }
4423
4424 4 void PeriodicDelaunay3d::Stats::reset() {
4425 Memory::clear(this, sizeof(Stats));
4426
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 raw_ = CmdLine::get_arg_bool("dbg:raw_logs");
4427 4 }
4428
4429 ✗ std::string PeriodicDelaunay3d::Stats::to_string_raw() const {
4430 return String::format(
4431 "%.1f %.1f %.1f %.1f %d %d %d %.1f %d %.1f %.1f %.1f %d",
4432 ✗ total_t_,
4433
4434 ✗ phase_0_t_,
4435
4436 ✗ phase_I_t_, phase_I_classify_t_,
4437
4438 ✗ int(phase_I_nb_inside_), int(phase_I_nb_cross_),
4439 ✗ int(phase_I_nb_outside_),
4440
4441 ✗ phase_I_insert_t_, int(phase_I_insert_nb_),
4442
4443 ✗ phase_II_t_, phase_II_classify_t_, phase_II_insert_t_,
4444 ✗ int(phase_II_insert_nb_)
4445 ✗ );
4446 }
4447
4448
4449 ✗ std::string PeriodicDelaunay3d::Stats::to_string_pretty() const {
4450 return String::format(
4451 "total | t:%.1f\n"
4452 "phase0 | t:%.1f\n"
4453 "phaseI | t:%.1f t_cls:%.1f t_ins:%.1f nb_ins:%d\n"
4454 " | in:%d bndry:%d out:%d\n"
4455 "phaseII | t:%.1f t_cls:%.1f t_ins:%.1f nb_ins:%d",
4456 ✗ total_t_,
4457
4458 ✗ phase_0_t_,
4459
4460 ✗ phase_I_t_, phase_I_classify_t_,
4461 ✗ phase_I_insert_t_, int(phase_I_insert_nb_),
4462
4463 ✗ int(phase_I_nb_inside_), int(phase_I_nb_cross_),
4464 ✗ int(phase_I_nb_outside_),
4465
4466 ✗ phase_II_t_, phase_II_classify_t_,
4467 ✗ phase_II_insert_t_, int(phase_II_insert_nb_)
4468 ✗ );
4469 }
4470 }
4471