GCC Code Coverage Report

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