GCC Code Coverage Report

Directory: ./
File: lib/geogram/voronoi/RVD.cpp
Date: 2026-09-07 02:25:23
Exec Total Coverage
Lines: 341 855 39.9%
Functions: 78 376 20.7%
Branches: 219 1080 20.3%
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/voronoi/RVD.h>
41 #include <geogram/voronoi/generic_RVD.h>
42 #include <geogram/voronoi/RVD_mesh_builder.h>
43 #include <geogram/voronoi/integration_simplex.h>
44 #include <geogram/voronoi/RVD_callback.h>
45 #include <geogram/mesh/mesh_partition.h>
46 #include <geogram/mesh/mesh_sampling.h>
47 #include <geogram/mesh/mesh_repair.h>
48 #include <geogram/mesh/mesh_AABB.h>
49 #include <geogram/delaunay/delaunay.h>
50 #include <geogram/basic/geometry_nd.h>
51 #include <geogram/basic/process.h>
52 #include <geogram/basic/command_line.h>
53 #include <geogram/basic/argused.h>
54 #include <geogram/basic/algorithm.h>
55 #include <geogram/bibliography/bibliography.h>
56
57 /*
58 * There are three levels of implementation:
59 * Level 1: RestrictedVoronoiDiagram is the abstract API seen from client code
60 * Level 2: RVD_Nd_Impl<DIM> implements RestrictedVoronoiDiagram
61 * Level 3: RVD_Nd_Impl<DIM>::GenRestrictedVoronoiDiagram is
62 * an instantiation of GEOGen::RestrictedVoronoiDiagram (from generic_RVD.h)
63 *
64 * Warning: there are approx. 1000 lines of boring code ahead.
65 */
66
67 namespace {
68
69 using namespace GEO;
70
71 /**
72 * \brief Generic implementation of RestrictedVoronoiDiagram.
73 * \tparam DIM dimension
74 */
75 template <unsigned int DIM>
76 class RVD_Nd_Impl : public GEO::RestrictedVoronoiDiagram {
77
78 /** \brief This class type */
79 typedef RVD_Nd_Impl<DIM> thisclass;
80
81 /** \brief The base class of this class */
82 typedef RestrictedVoronoiDiagram baseclass;
83
84 public:
85 /** \brief Implementation based on the generic version. */
86 typedef GEOGen::RestrictedVoronoiDiagram<DIM>
87 GenRestrictedVoronoiDiagram;
88
89 /** \brief Representation of points. */
90 typedef vecng<DIM, double> Point;
91
92 /** \brief Representation of vectors. */
93 typedef vecng<DIM, double> Vector;
94
95 /** \brief Represents a point and its symbolic information. */
96 typedef typename GenRestrictedVoronoiDiagram::Vertex Vertex;
97
98 /**
99 * \brief Specifies the computation done by the threads.
100 */
101 enum ThreadMode {
102 MT_NONE, /**< uninitialized */
103 MT_LLOYD, /**< Lloyd iteration */
104 MT_NEWTON, /**< Newton optimization */
105 MT_INT_SMPLX, /**< Newton with integration simplex */
106 MT_POLYG, /**< Polygon callback */
107 MT_POLYH /**< Polyhedron callback */
108 };
109
110 /**
111 * \brief Gets a mesh vertex from its index.
112 * \param[in] v index of the vertex
113 * \return a const reference to a Point
114 */
115 const Point& mesh_vertex(index_t v) {
116 ✗ return *(const Point*) mesh_->vertices.point_ptr(v);
117 }
118
119
120 /**
121 * \brief Creates a RVD_Nd_Impl.
122 *
123 * \details The dimension is determined by \p mesh->dimension().
124 * \param[in] delaunay the Delaunay triangulation
125 * \param[in] mesh the input mesh
126 * \param[in] R3_embedding gives for each vertex
127 * its mapping in 3D space.
128 * \param[in] R3_embedding_stride gives the stride between
129 * two consecutive vertices in R3_embedding
130 */
131 92 RVD_Nd_Impl(
132 Delaunay* delaunay, Mesh* mesh,
133 const double* R3_embedding, index_t R3_embedding_stride
134 ) :
135 RestrictedVoronoiDiagram(
136 delaunay, mesh, R3_embedding, R3_embedding_stride
137 ),
138 92 RVD_(delaunay, mesh) {
139 92 use_exact_projection_ = false;
140 92 is_slave_ = false;
141 92 master_ = nullptr;
142
1/2
✓ Branch 1 taken 46 times.
✗ Branch 2 not taken.
 92 has_weights_ = false;
143
2/4
✓ Branch 1 taken 46 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 46 times.
 184 if(mesh->vertices.attributes().is_defined("weight")) {
144 ✗ vertex_weight_.bind(mesh->vertices.attributes(), "weight");
145 ✗ has_weights_ = true;
146 }
147 92 parts_ = nullptr;
148 92 nb_parts_ = 0;
149 92 funcval_ = 0.0;
150 92 simplex_func_ = nullptr;
151 92 polygon_callback_ = nullptr;
152 92 polyhedron_callback_ = nullptr;
153 92 arg_vectors_ = nullptr;
154 92 arg_scalars_ = nullptr;
155 92 thread_mode_ = MT_NONE;
156 92 nb_triangles_ = 0;
157 92 }
158
159 /**
160 * \brief Constructor for parts, used in multithreading mode.
161 */
162 80 RVD_Nd_Impl() :
163 RestrictedVoronoiDiagram(nullptr, nullptr, nullptr, 0),
164 80 RVD_(nullptr, nullptr) {
165 80 use_exact_projection_ = false;
166 80 is_slave_ = true;
167 80 master_ = nullptr;
168 80 mesh_ = nullptr;
169 80 parts_ = nullptr;
170 80 nb_parts_ = 0;
171 80 facets_begin_ = NO_INDEX;
172 80 facets_end_ = NO_INDEX;
173 80 funcval_ = 0.0;
174 80 simplex_func_ = nullptr;
175 80 polygon_callback_ = nullptr;
176 80 polyhedron_callback_ = nullptr;
177 80 arg_vectors_ = nullptr;
178 80 arg_scalars_ = nullptr;
179 80 thread_mode_ = MT_NONE;
180 80 nb_triangles_ = 0;
181 80 }
182
183 80 void set_delaunay(Delaunay* delaunay) override {
184 80 baseclass::set_delaunay(delaunay);
185 RVD_.set_delaunay(delaunay);
186
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
 80 for(index_t p = 0; p < nb_parts_; ++p) {
187 ✗ parts_[p].set_delaunay(delaunay);
188 }
189 80 }
190
191 152 void set_check_SR(bool x) override {
192 RVD_.set_check_SR(x);
193
6/14
✓ Branch 0 taken 56 times.
✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 304 for(index_t p = 0; p < nb_parts_; ++p) {
194 112 parts_[p].set_check_SR(x);
195 }
196 152 }
197
198
2/2
✓ Branch 0 taken 55 times.
✓ Branch 1 taken 40 times.
 190 void set_exact_predicates(bool x) override {
199 RVD_.set_exact_predicates(x);
200
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 95 times.
 222 for(index_t p = 0; p < nb_parts_; ++p) {
201 32 parts_[p].set_exact_predicates(x);
202 }
203 190 }
204
205 ✗ bool exact_predicates() const override {
206 ✗ return RVD_.exact_predicates();
207 }
208
209 /********************************************************************/
210
211 /**
212 * \brief Place holder, "no locking" policy.
213 * \details NoLocks is used by algorithms templated
214 * by locking policy, for the single-threaded instances
215 * that do not need synchronization. The multi-threaded
216 * instances are parameterized by SpinLockArray.
217 */
218 class NoLocks {
219 public:
220 /**
221 * \brief Acquires a spinlock.
222 * \details Does nothing in this version
223 * \param[in] i index of the spinlock to acquire
224 */
225 void acquire_spinlock(index_t i) {
226 geo_argused(i);
227 }
228
229 /**
230 * \brief Releases a spinlock.
231 * \details Does nothing in this version
232 * \param[in] i index of the spinlock to release
233 */
234 void release_spinlock(index_t i) {
235 geo_argused(i);
236 }
237 };
238
239 /******************************************************************/
240
241 /**
242 * \brief Implementation class for surfacic Lloyd relaxation.
243 * \details To be used as a template argument
244 * to RVD::for_each_triangle().
245 * This version ignores the weights.
246 *
247 * Computes for each RVD cell:
248 * - mg[v] (v's Voronoi cell's total area times centroid)
249 * - m[v] (v's total area)
250 * \tparam LOCKS locking policy
251 * (can be one of Process::SpinLockArray, NoLocks)
252 */
253 template <class LOCKS>
254 class ComputeCentroids {
255 public:
256 /**
257 * \brief Constructs a ComputeCentroids.
258 * \param[out] mg where to store the centroids
259 * \param[out] m where to store the masses
260 * \param[in] locks the array of locks
261 * (or NoLocks in single thread mode)
262 */
263 960 ComputeCentroids(
264 double* mg,
265 double* m,
266 LOCKS& locks
267 ) :
268 960 mg_(mg),
269 960 m_(m),
270 960 locks_(locks) {
271 }
272
273 /**
274 * \brief The callback called for each integration simplex.
275 * \param[in] v index of current center vertex
276 * \param[in] p1 first vertex of current integration simplex
277 * \param[in] p2 second vertex of current integration simplex
278 * \param[in] p3 third vertex of current integration simplex
279 */
280 16224392 void operator() (
281 index_t v,
282 const double* p1,
283 const double* p2,
284 const double* p3
285 ) const {
286 16224392 double cur_m = Geom::triangle_area(p1, p2, p3, DIM);
287 16224392 double s = cur_m / 3.0;
288 16224392 locks_.acquire_spinlock(v);
289 16224392 m_[v] += cur_m;
290 16224392 double* cur_mg_out = mg_ + v * DIM;
291
2/2
✓ Branch 0 taken 41805768 times.
✓ Branch 1 taken 8112196 times.
 99835928 for(coord_index_t coord = 0; coord < DIM; coord++) {
292 83611536 cur_mg_out[coord] +=
293 83611536 s * (p1[coord] + p2[coord] + p3[coord]);
294 }
295 16224392 locks_.release_spinlock(v);
296 16224392 }
297
298 private:
299 double* mg_;
300 double* m_;
301 LOCKS& locks_;
302 };
303
304 /**
305 * \brief Implementation class for surfacic Lloyd relaxation.
306 * \details To be used as a template
307 * argument to RVD::for_each_triangle().
308 * This version takes the weights into account.
309 *
310 * Computes for each RVD cell:
311 * - mg[v] (v's Voronoi cell's total area times centroid)
312 * - m[v] (v's total area)
313 * \tparam LOCKS locking policy
314 * (can be one of Process::SpinLockArray, NoLocks)
315 */
316 template <class LOCKS>
317 class ComputeCentroidsWeighted {
318 public:
319 /**
320 * \brief Constructs a ComputeCentroidsWeighted.
321 * \param[out] mg where to store the centroids
322 * \param[out] m where to store the masses
323 * \param[in] locks the array of locks
324 * (or NoLocks in single thread mode)
325 */
326 ✗ ComputeCentroidsWeighted(
327 double* mg,
328 double* m,
329 LOCKS& locks
330 ) :
331 ✗ mg_(mg),
332 ✗ m_(m),
333 ✗ locks_(locks) {
334 }
335
336 /**
337 * \brief The callback called for each integration simplex.
338 * \param[in] v index of current center vertex
339 * \param[in] v1 first vertex of current integration simplex
340 * \param[in] v2 second vertex of current integration simplex
341 * \param[in] v3 third vertex of current integration simplex
342 */
343 ✗ void operator() (
344 index_t v,
345 const Vertex& v1,
346 const Vertex& v2,
347 const Vertex& v3
348 ) const {
349 double cur_m;
350 double cur_Vg[DIM];
351 ✗ Geom::triangle_centroid(
352 v1.point(), v2.point(), v3.point(),
353 v1.weight(), v2.weight(), v3.weight(),
354 cur_Vg, cur_m, DIM
355 );
356 ✗ locks_.acquire_spinlock(v);
357 ✗ m_[v] += cur_m;
358 ✗ double* cur_mg_out = mg_ + v * DIM;
359 ✗ for(coord_index_t coord = 0; coord < DIM; coord++) {
360 ✗ cur_mg_out[coord] += cur_Vg[coord];
361 }
362 ✗ locks_.release_spinlock(v);
363 ✗ }
364
365 private:
366 double* mg_;
367 double* m_;
368 LOCKS& locks_;
369 };
370
371 2400 void compute_centroids_on_surface(double* mg, double* m) override {
372 2400 create_threads();
373 2400 if(nb_parts() == 0) {
374
1/2
✓ Branch 0 taken 960 times.
✗ Branch 1 not taken.
 1920 if(master_ != nullptr) {
375
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 960 times.
 1920 if(has_weights_) {
376 ✗ RVD_.for_each_triangle(
377 ✗ ComputeCentroidsWeighted<Process::SpinLockArray>(
378 ✗ mg, m, master_->spinlocks_
379 )
380 );
381 } else {
382 1920 RVD_.for_each_triangle(
383 1920 ComputeCentroids<Process::SpinLockArray>(
384 1920 mg, m, master_->spinlocks_
385 )
386 );
387 }
388 } else {
389 NoLocks nolocks;
390 ✗ if(has_weights_) {
391 ✗ RVD_.for_each_triangle(
392 ✗ ComputeCentroidsWeighted<NoLocks>(
393 mg, m, nolocks
394 )
395 );
396 } else {
397 ✗ RVD_.for_each_triangle(
398 ✗ ComputeCentroids<NoLocks>(mg, m, nolocks)
399 );
400 }
401 }
402 } else {
403 480 thread_mode_ = MT_LLOYD;
404 480 arg_vectors_ = mg;
405 480 arg_scalars_ = m;
406 480 spinlocks_.resize(delaunay_->nb_vertices());
407 480 parallel_for(
408 0, nb_parts(),
409 1440 [this](index_t i) { run_thread(i); }
410 );
411 }
412 2400 }
413
414 /********************************************************************/
415
416 /**
417 * \brief Implementation class for surfacic Lloyd relaxation.
418 * \details To be used as a template argument
419 * to RVD::for_each_volumetric_integration_simplex().
420 * This version ignores the weights.
421 *
422 * Computes for each RVD cell:
423 * - mg[v] (v's Voronoi cell's total area times centroid)
424 * - m[v] (v's total area)
425 * \tparam LOCKS locking policy
426 * (can be one of Process::SpinLockArray, NoLocks)
427 */
428 template <class LOCKS>
429 class ComputeCentroidsVolumetric {
430 public:
431 /**
432 * \brief Constructs a ComputeCentroidsVolumetric.
433 * \param[out] mg where to store the centroids
434 * \param[out] m where to store the masses
435 * \param[in] delaunay the Delaunay triangulation
436 * \param[in] locks the array of locks
437 * (or NoLocks in single thread mode)
438 */
439 640 ComputeCentroidsVolumetric(
440 double* mg,
441 double* m,
442 const Delaunay* delaunay,
443 LOCKS& locks
444 ) :
445 640 mg_(mg),
446 640 m_(m),
447 640 delaunay_(delaunay),
448 640 locks_(locks) {
449 }
450
451 /**
452 * \brief The callback called for each integration simplex.
453 * \param[in] v index of current center vertex
454 * \param[in] v_adj (unused here) is the index of the Voronoi cell
455 * adjacent to t accros facet (\p v1, \p v2, \p v3) or
456 * NO_INDEX if it does not exists
457 * \param[in] t (unused here) is the index of the current
458 * tetrahedron
459 * \param[in] t_adj (unused here) is the index of the
460 * tetrahedron adjacent to t accros facet (\p v1, \p v2, \p v3)
461 * or NO_INDEX if it does not exists
462 * \param[in] p0 first vertex of current integration simplex
463 * \param[in] p1 second vertex of current integration simplex
464 * \param[in] p2 third vertex of current integration simplex
465 * \param[in] p3 fourth vertex of current integration simplex
466 */
467 11666782 void operator() (
468 index_t v, index_t v_adj,
469 index_t t, index_t t_adj,
470 const double* p0,
471 const double* p1,
472 const double* p2,
473 const double* p3
474 ) const {
475 geo_argused(v_adj);
476 geo_argused(t);
477 geo_argused(t_adj);
478 9389046 double cur_m = Geom::tetra_volume<DIM>(
479 p0, p1, p2, p3
480 );
481 11666782 double s = cur_m / 4.0;
482 11666782 locks_.acquire_spinlock(v);
483 11666782 m_[v] += cur_m;
484 11666782 double* cur_mg_out = mg_ + v * DIM;
485
2/2
✓ Branch 0 taken 32227042 times.
✓ Branch 1 taken 5833391 times.
 76120866 for(coord_index_t coord = 0; coord < DIM; coord++) {
486 64454084 cur_mg_out[coord] += s * (
487 64454084 p0[coord] + p1[coord] + p2[coord] + p3[coord]
488 );
489 }
490 11666782 locks_.release_spinlock(v);
491 11666782 }
492
493 private:
494 double* mg_;
495 double* m_;
496 const Delaunay* delaunay_;
497 LOCKS& locks_;
498 };
499
500 1600 void compute_centroids_in_volume(double* mg, double* m) override {
501 1600 create_threads();
502 1600 if(nb_parts() == 0) {
503
1/2
✓ Branch 0 taken 640 times.
✗ Branch 1 not taken.
 1280 if(master_ != nullptr) {
504 1280 RVD_.for_each_tetrahedron(
505 1280 ComputeCentroidsVolumetric<Process::SpinLockArray>(
506 1280 mg, m, RVD_.delaunay(), master_->spinlocks_
507 )
508 );
509 } else {
510 NoLocks nolocks;
511 ✗ RVD_.for_each_tetrahedron(
512 ✗ ComputeCentroidsVolumetric<NoLocks>(
513 mg, m, RVD_.delaunay(), nolocks
514 )
515 );
516 }
517 } else {
518 320 thread_mode_ = MT_LLOYD;
519 320 arg_vectors_ = mg;
520 320 arg_scalars_ = m;
521 320 spinlocks_.resize(delaunay_->nb_vertices());
522 320 parallel_for(
523 0, nb_parts(),
524 960 [this](index_t i) { run_thread(i); }
525 );
526 }
527 1600 }
528
529 /********************************************************************/
530
531 /**
532 * \brief Implementation class for Newton-based restricted CVT.
533 * \details To be used as a template argument
534 * to RVD::for_each_triangle().
535 * This version ignores the weights.
536 *
537 * Computes for each RVD cell:
538 * - g (gradient)
539 * - f (CVT energy)
540 * \tparam LOCKS locking policy
541 * (can be one of Process::SpinLockArray, NoLocks)
542 */
543 template <class LOCKS>
544 class ComputeCVTFuncGrad {
545 public:
546 /**
547 * \brief Constructs a ComputeCVTFuncGrad.
548 * \param[in] RVD the restricted Voronoi diagram
549 * \param[out] f the computed function value
550 * \param[out] g the computed gradient of f,
551 * allocated by caller, and managed
552 * by caller
553 * \param[in] locks the array of locks
554 * (or NoLocks in single thread mode)
555 */
556 ✗ ComputeCVTFuncGrad(
557 const GenRestrictedVoronoiDiagram& RVD,
558 double& f,
559 double* g,
560 LOCKS& locks
561 ) :
562 ✗ f_(f),
563 ✗ g_(g),
564 ✗ locks_(locks),
565 ✗ RVD_(RVD) {
566 }
567
568 /**
569 * \brief The callback called for each integration simplex.
570 * \param[in] v index of current center vertex
571 * \param[in] p1 first vertex of current integration simplex
572 * \param[in] p2 second vertex of current integration simplex
573 * \param[in] p3 third vertex of current integration simplex
574 */
575 ✗ void operator() (
576 index_t v,
577 const double* p1,
578 const double* p2,
579 const double* p3
580 ) const {
581
582 ✗ const double* p0 = RVD_.delaunay()->vertex_ptr(v);
583
584 ✗ double t_area = Geom::triangle_area(p1, p2, p3, DIM);
585
586 double cur_f = 0.0;
587 ✗ for(index_t c = 0; c < DIM; c++) {
588 ✗ double u0 = p0[c] - p1[c];
589 ✗ double u1 = p0[c] - p2[c];
590 ✗ double u2 = p0[c] - p3[c];
591 ✗ cur_f += u0 * u0;
592 ✗ cur_f += u1 * (u0 + u1);
593 ✗ cur_f += u2 * (u0 + u1 + u2);
594 }
595
596 ✗ f_ += t_area * cur_f / 6.0;
597
598 ✗ locks_.acquire_spinlock(v);
599 ✗ for(index_t c = 0; c < DIM; c++) {
600 ✗ double Gc = (1.0 / 3.0) * (p1[c] + p2[c] + p3[c]);
601 ✗ g_[DIM * v + c] += (2.0 * t_area) * (p0[c] - Gc);
602 }
603 ✗ locks_.release_spinlock(v);
604 ✗ }
605
606 double& f_;
607 double* g_;
608 LOCKS& locks_;
609 const GenRestrictedVoronoiDiagram& RVD_;
610 };
611
612 /**
613 * \brief Implementation class for Newton-based restricted CVT.
614 * \details To be used as a template argument
615 * to RVD::for_each_triangle().
616 * This version takes the weights into account.
617 *
618 * Computes for each RVD cell:
619 * - g (gradient)
620 * - f (CVT energy)
621 * \tparam LOCKS locking policy
622 * (can be one of Process::SpinLockArray, NoLocks)
623 */
624 template <class LOCKS>
625 class ComputeCVTFuncGradWeighted {
626 public:
627 /**
628 * \brief Constructs a ComputeCVTFuncGradWeighted.
629 * \param[in] RVD the restricted Voronoi diagram
630 * \param[out] f the computed function value
631 * \param[out] g the computed gradient of f,
632 * allocated by caller, and managed by caller
633 * \param[in] locks the array of locks
634 * (or NoLocks in single thread mode)
635 */
636 ✗ ComputeCVTFuncGradWeighted(
637 const GenRestrictedVoronoiDiagram& RVD,
638 double& f,
639 double* g,
640 LOCKS& locks
641 ) :
642 ✗ f_(f),
643 ✗ g_(g),
644 ✗ locks_(locks),
645 ✗ RVD_(RVD) {
646 }
647
648 /**
649 * \brief The callback called for each integration simplex.
650 * \param[in] v index of current center vertex
651 * \param[in] v1 first vertex of current integration simplex
652 * \param[in] v2 second vertex of current integration simplex
653 * \param[in] v3 third vertex of current integration simplex
654 */
655 ✗ void operator() (
656 index_t v,
657 const Vertex& v1,
658 const Vertex& v2,
659 const Vertex& v3
660 ) const {
661
662 ✗ const double* p0 = RVD_.delaunay()->vertex_ptr(v);
663
664 const double* p1 = v1.point();
665 const double* p2 = v2.point();
666 const double* p3 = v3.point();
667
668 ✗ double t_area = Geom::triangle_area(p1, p2, p3, DIM);
669
670 ✗ double Sp = v1.weight() + v2.weight() + v3.weight();
671 double rho[3], alpha[3];
672 rho[0] = v1.weight();
673 rho[1] = v2.weight();
674 rho[2] = v3.weight();
675 ✗ alpha[0] = Sp + rho[0];
676 ✗ alpha[1] = Sp + rho[1];
677 ✗ alpha[2] = Sp + rho[2];
678
679 double dotprod_00 = 0.0;
680 double dotprod_10 = 0.0;
681 double dotprod_11 = 0.0;
682 double dotprod_20 = 0.0;
683 double dotprod_21 = 0.0;
684 double dotprod_22 = 0.0;
685 ✗ for(unsigned int c = 0; c < DIM; c++) {
686 ✗ double sp0 = p0[c] - p1[c];
687 ✗ double sp1 = p0[c] - p2[c];
688 ✗ double sp2 = p0[c] - p3[c];
689 ✗ dotprod_00 += sp0 * sp0;
690 ✗ dotprod_10 += sp1 * sp0;
691 ✗ dotprod_11 += sp1 * sp1;
692 ✗ dotprod_20 += sp2 * sp0;
693 ✗ dotprod_21 += sp2 * sp1;
694 ✗ dotprod_22 += sp2 * sp2;
695 }
696
697 double cur_f = 0.0;
698 ✗ cur_f += (alpha[0] + rho[0]) * dotprod_00; // 0 0
699 ✗ cur_f += (alpha[1] + rho[0]) * dotprod_10; // 1 0
700 ✗ cur_f += (alpha[1] + rho[1]) * dotprod_11; // 1 1
701 ✗ cur_f += (alpha[2] + rho[0]) * dotprod_20; // 2 0
702 ✗ cur_f += (alpha[2] + rho[1]) * dotprod_21; // 2 1
703 ✗ cur_f += (alpha[2] + rho[2]) * dotprod_22; // 2 2
704
705 ✗ f_ += t_area * cur_f / 30.0;
706 ✗ double* g_out = g_ + v * DIM;
707 ✗ locks_.acquire_spinlock(v);
708 ✗ for(index_t c = 0; c < DIM; c++) {
709 ✗ g_out[c] += (t_area / 6.0) * (
710 ✗ 4.0 * Sp * p0[c] - (
711 ✗ alpha[0] * p1[c] +
712 ✗ alpha[1] * p2[c] +
713 ✗ alpha[2] * p3[c]
714 )
715 );
716 }
717 ✗ locks_.release_spinlock(v);
718 ✗ }
719
720 double& f_;
721 double* g_;
722 LOCKS& locks_;
723 const GenRestrictedVoronoiDiagram& RVD_;
724 };
725
726 ✗ void compute_CVT_func_grad_on_surface(double& f, double* g) override {
727 ✗ create_threads();
728 ✗ if(nb_parts() == 0) {
729 ✗ if(master_ != nullptr) {
730 ✗ if(has_weights_) {
731 ✗ RVD_.for_each_triangle(
732 ✗ ComputeCVTFuncGradWeighted<Process::SpinLockArray>(
733 ✗ RVD_, f, g, master_->spinlocks_
734 )
735 );
736 } else {
737 ✗ RVD_.for_each_triangle(
738 ✗ ComputeCVTFuncGrad<Process::SpinLockArray>(
739 ✗ RVD_, f, g, master_->spinlocks_
740 )
741 );
742 }
743 } else {
744 NoLocks nolocks;
745 ✗ if(has_weights_) {
746 ✗ RVD_.for_each_triangle(
747 ✗ ComputeCVTFuncGradWeighted<NoLocks>(
748 RVD_, f, g, nolocks
749 )
750 );
751 } else {
752 ✗ RVD_.for_each_triangle(
753 ✗ ComputeCVTFuncGrad<NoLocks>(
754 RVD_, f, g, nolocks
755 )
756 );
757 }
758 }
759 } else {
760 ✗ thread_mode_ = MT_NEWTON;
761 ✗ arg_vectors_ = g;
762 ✗ spinlocks_.resize(delaunay_->nb_vertices());
763 ✗ for(index_t t = 0; t < nb_parts(); t++) {
764 ✗ part(t).funcval_ = 0.0;
765 }
766 ✗ parallel_for(
767 0, nb_parts(),
768 ✗ [this](index_t i) { run_thread(i); }
769 );
770 ✗ for(index_t t = 0; t < nb_parts(); t++) {
771 ✗ f += part(t).funcval_;
772 }
773 }
774 ✗ }
775
776 /********************************************************************/
777
778 /**
779 * \brief Implementation class for Newton-based restricted CVT
780 * in volume.
781 * \details To be used as a template argument
782 * to RVD::for_each_volumetric_integration_simplex().
783 * This version ignores the weights.
784 *
785 * Computes for each RVD cell:
786 * - g (gradient)
787 * - f (CVT energy)
788 * \tparam LOCKS locking policy
789 * (can be one of Process::SpinLockArray, NoLocks)
790 */
791 template <class LOCKS>
792 class ComputeCVTFuncGradVolumetric {
793 public:
794 /**
795 * \brief Constructs a ComputeCentroidsFuncGradVolumetric.
796 * \param[in] RVD the restricted Voronoi diagram
797 * \param[out] f the computed function value
798 * \param[out] g the computed gradient of f,
799 * allocated by caller, and managed
800 * by caller
801 * \param[in] locks the array of locks
802 * (or NoLocks in single thread mode)
803 */
804 ✗ ComputeCVTFuncGradVolumetric(
805 const GenRestrictedVoronoiDiagram& RVD,
806 double& f,
807 double* g,
808 LOCKS& locks
809 ) :
810 ✗ f_(f),
811 ✗ g_(g),
812 ✗ locks_(locks),
813 ✗ RVD_(RVD) {
814 }
815
816 /**
817 * \brief The callback called for each integration simplex.
818 * \param[in] v index of current center vertex
819 * \param[in] v_adj (unused here) is the index of the Voronoi cell
820 * adjacent to t accros facet (\p v1, \p v2, \p v3) or
821 * NO_INDEX if it does not exists
822 * \param[in] t (unused here) is the index of the current
823 * tetrahedron
824 * \param[in] t_adj (unused here) is the index of the
825 * tetrahedron adjacent to t accros facet (\p v1, \p v2, \p v3)
826 * or NO_INDEX if it does not exists
827 * \param[in] p1 first vertex of current integration simplex
828 * \param[in] p2 second vertex of current integration simplex
829 * \param[in] p3 third vertex of current integration simplex
830 */
831 ✗ void operator() (
832 index_t v,
833 index_t v_adj,
834 index_t t,
835 index_t t_adj,
836 const double* p1,
837 const double* p2,
838 const double* p3
839 ) const {
840 geo_argused(v_adj);
841 geo_argused(t);
842 geo_argused(t_adj);
843 ✗ const double* p0 = RVD_.delaunay()->vertex_ptr(v);
844
845 ✗ double mi = Geom::tetra_volume<DIM>(p0, p1, p2, p3);
846
847 // fi = (mi/10)*(U.U + V.V + W.W + U.V + V.W + W.U)
848 // where: U = p1-p0 ; V = p2-p0 and W=p3-p0
849 double fi = 0.0;
850 ✗ for(coord_index_t c = 0; c < DIM; ++c) {
851 ✗ double Uc = p1[c] - p0[c];
852 ✗ double Vc = p2[c] - p0[c];
853 ✗ double Wc = p3[c] - p0[c];
854 ✗ fi += geo_sqr(Uc) + geo_sqr(Vc) + geo_sqr(Wc);
855 ✗ fi += (Uc * Vc + Vc * Wc + Wc * Uc);
856 }
857 ✗ fi *= (mi / 10.0);
858 ✗ f_ += fi;
859
860 // gi = 2*mi(p0 - 1/4(p0 + p1 + p2 + p3))
861 ✗ double* g_out = g_ + v * DIM;
862 ✗ locks_.acquire_spinlock(v);
863 ✗ for(coord_index_t c = 0; c < DIM; ++c) {
864 ✗ g_out[c] += 2.0 * mi * (
865 ✗ 0.75 * p0[c]
866 ✗ - 0.25 * p1[c] - 0.25 * p2[c] - 0.25 * p3[c]
867 );
868 }
869 ✗ locks_.release_spinlock(v);
870 ✗ }
871
872 double& f_;
873 double* g_;
874 LOCKS& locks_;
875 const GenRestrictedVoronoiDiagram& RVD_;
876 };
877
878 ✗ void compute_CVT_func_grad_in_volume(double& f, double* g) override {
879 ✗ create_threads();
880 ✗ if(nb_parts() == 0) {
881 ✗ if(master_ != nullptr) {
882 ✗ RVD_.for_each_volumetric_integration_simplex(
883 ✗ ComputeCVTFuncGradVolumetric<Process::SpinLockArray>(
884 ✗ RVD_, f, g, master_->spinlocks_
885 )
886 );
887 } else {
888 NoLocks nolocks;
889 ✗ RVD_.for_each_volumetric_integration_simplex(
890 ✗ ComputeCVTFuncGradVolumetric<NoLocks>(
891 RVD_, f, g, nolocks
892 )
893 );
894 }
895 } else {
896 ✗ thread_mode_ = MT_NEWTON;
897 ✗ arg_vectors_ = g;
898 ✗ spinlocks_.resize(delaunay_->nb_vertices());
899 ✗ for(index_t t = 0; t < nb_parts(); t++) {
900 ✗ part(t).funcval_ = 0.0;
901 }
902 ✗ parallel_for(
903 0, nb_parts(),
904 ✗ [this](index_t i) { run_thread(i); }
905 );
906 ✗ for(index_t t = 0; t < nb_parts(); t++) {
907 ✗ f += part(t).funcval_;
908 }
909 }
910 ✗ }
911
912 /********************************************************************/
913
914 /**
915 * \brief Implementation class for computing function integrals
916 * and gradients over integration simplices.
917 * \details To be used as a template argument
918 * to RVD::for_each_triangle() and
919 * RVD::for_each_volumetric_integration_simplex()
920 */
921 class ComputeCVTFuncGradIntegrationSimplex {
922 public:
923 /**
924 * \brief Constructs a ComputeCVTFuncGradIntegrationSimplex.
925 * \param[in] RVD the Restricted Voronoi Diagram
926 * \param[in] F the IntegrationSimplex
927 */
928 ✗ ComputeCVTFuncGradIntegrationSimplex(
929 const GenRestrictedVoronoiDiagram& RVD,
930 IntegrationSimplex* F
931 ) :
932 ✗ f_(0.0),
933 ✗ RVD_(RVD),
934 ✗ simplex_func_(F) {
935 ✗ simplex_func_->reset_thread_local_storage();
936 }
937
938 /**
939 * \brief The callback called for each surfacic integration simplex.
940 * \param[in] i index of current center vertex
941 * \param[in] v1 first vertex of current integration simplex
942 * \param[in] v2 second vertex of current integration simplex
943 * \param[in] v3 third vertex of current integration simplex
944 */
945 void operator() (
946 index_t i,
947 const Vertex& v1,
948 const Vertex& v2,
949 const Vertex& v3
950 ) {
951 ✗ f_ += simplex_func_->eval(
952 ✗ i,v1,v2,v3,RVD_.current_facet()
953 );
954 }
955
956 /**
957 * \brief The callback called for each volumetric
958 * integration simplex.
959 * \param[in] v index of current center vertex
960 * \param[in] v_adj index of the Voronoi cell adjacent to t accros
961 * facet (\p v1, \p v2, \p v3) or NO_INDEX if it does not exists
962 * \param[in] t index of the current tetrahedron
963 * \param[in] t_adj index of the tetrahedron adjacent to t accros
964 * facet (\p v1, \p v2, \p v3) or NO_INDEX if it does not exists
965 * \param[in] v1 first vertex of current integration simplex
966 * \param[in] v2 second vertex of current integration simplex
967 * \param[in] v3 third vertex of current integration simplex
968 */
969 void operator() (
970 index_t v,
971 index_t v_adj,
972 index_t t,
973 index_t t_adj,
974 const Vertex& v1,
975 const Vertex& v2,
976 const Vertex& v3
977 ) {
978 geo_argused(v_adj);
979 geo_argused(t_adj);
980 ✗ f_ += simplex_func_->eval(v,v1,v2,v3,t,t_adj,v_adj);
981 }
982
983 /**
984 * \brief Gets the function value.
985 * \return The function value accumulated so far.
986 */
987 double f() const {
988 ✗ return f_;
989 }
990
991 private:
992 double f_;
993 const GenRestrictedVoronoiDiagram& RVD_;
994 IntegrationSimplex* simplex_func_;
995 };
996
997 ✗ void compute_integration_simplex_func_grad(
998 double& f, double* g, IntegrationSimplex* F
999 ) override {
1000 ✗ create_threads();
1001 ✗ if(nb_parts() == 0) {
1002 ✗ if(master_ == nullptr) {
1003 F->set_points_and_gradient(
1004 delaunay()->dimension(),
1005 delaunay()->nb_vertices(),
1006 delaunay()->vertex_ptr(0),
1007 g
1008 );
1009 }
1010 ✗ ComputeCVTFuncGradIntegrationSimplex C(RVD_,F);
1011 bool sym = RVD_.symbolic();
1012 RVD_.set_symbolic(true);
1013 ✗ if(F->volumetric()) {
1014 RVD_.for_each_volumetric_integration_simplex(
1015 C,
1016 F->background_mesh_has_varying_attribute(),
1017 false /* Coherent triangles */
1018 );
1019 } else {
1020 RVD_.for_each_triangle(C);
1021 }
1022 RVD_.set_symbolic(sym);
1023 ✗ funcval_ = C.f();
1024 ✗ f = C.f();
1025 } else {
1026 ✗ spinlocks_.resize(delaunay_->nb_vertices());
1027 F->set_points_and_gradient(
1028 delaunay()->dimension(),
1029 delaunay()->nb_vertices(),
1030 delaunay()->vertex_ptr(0),
1031 g,
1032 &spinlocks_
1033 );
1034 ✗ thread_mode_ = MT_INT_SMPLX;
1035 ✗ arg_vectors_ = g;
1036 ✗ simplex_func_ = F;
1037 ✗ funcval_ = 0.0;
1038 ✗ for(index_t t = 0; t < nb_parts(); t++) {
1039 ✗ part(t).arg_vectors_ = g;
1040 ✗ part(t).simplex_func_ = F;
1041 ✗ part(t).funcval_ = 0.0;
1042 }
1043
1044 ✗ parallel_for(
1045 0, nb_parts(),
1046 ✗ [this](index_t i) { run_thread(i); }
1047 );
1048
1049 ✗ f = 0.0;
1050 ✗ for(index_t t = 0; t < nb_parts(); t++) {
1051 ✗ f += part(t).funcval_;
1052 }
1053 }
1054 ✗ }
1055
1056 /********************************************************************/
1057
1058 /**
1059 * \brief Adapter class used internally to implement for_each_polygon()
1060 * \details Gets the current triangle from the RVD and passes it back
1061 * to the callback. It is needed because GenericRVD::for_each_polygon()
1062 * does not pass the current triangle.
1063 */
1064 // TODO: pass it through all the callbacks, because it is ridiculous:
1065 // we pass it through the first levels, then throw it, then retrieve it
1066 // (see GenRVD)
1067 class PolygonCallbackAction {
1068 public:
1069 /**
1070 * \brief PolygonCallbackAction constructor
1071 * \param[in] RVD a pointer to the restricted Voronoi diagram
1072 * \param[in] callback a pointer to the PolygonCallback
1073 */
1074 ✗ PolygonCallbackAction(
1075 GenRestrictedVoronoiDiagram& RVD,
1076 GEO::RVDPolygonCallback& callback
1077 ) :
1078 ✗ RVD_(RVD),
1079 ✗ callback_(callback) {
1080 }
1081
1082 /**
1083 * \brief Callback called for each polygon.
1084 * \details Routes the callback to the wrapped user action class.
1085 * \param[in] v index of current Delaunay seed
1086 * \param[in] P intersection between current mesh facet
1087 * and the Voronoi cell of \p v
1088 */
1089 void operator() (
1090 index_t v,
1091 const GEOGen::Polygon& P
1092 ) const {
1093 ✗ callback_(v, RVD_.current_facet(), P);
1094 }
1095
1096 protected:
1097 GenRestrictedVoronoiDiagram& RVD_;
1098 GEO::RVDPolygonCallback& callback_;
1099 };
1100
1101
1102 ✗ virtual void compute_with_polygon_callback(
1103 GEO::RVDPolygonCallback& polygon_callback
1104 ) {
1105 ✗ create_threads();
1106 ✗ if(nb_parts() == 0) {
1107 ✗ PolygonCallbackAction action(RVD_,polygon_callback);
1108 RVD_.for_each_polygon(action);
1109 } else {
1110 ✗ for(index_t t = 0; t < nb_parts(); t++) {
1111 part(t).RVD_.set_symbolic(RVD_.symbolic());
1112 part(t).RVD_.set_connected_components_priority(
1113 RVD_.connected_components_priority()
1114 );
1115 }
1116 ✗ spinlocks_.resize(delaunay_->nb_vertices());
1117 ✗ thread_mode_ = MT_POLYG;
1118 ✗ polygon_callback_ = &polygon_callback;
1119 polygon_callback_->set_spinlocks(&spinlocks_);
1120 // Note: callback begin()/end() is called in for_each_polygon()
1121 ✗ parallel_for(
1122 0, nb_parts(),
1123 ✗ [this](index_t i) { run_thread(i); }
1124 );
1125 ✗ polygon_callback_->set_spinlocks(nullptr);
1126 }
1127 ✗ }
1128
1129 ✗ virtual void compute_with_polyhedron_callback(
1130 GEO::RVDPolyhedronCallback& polyhedron_callback
1131 ) {
1132 ✗ create_threads();
1133 ✗ if(nb_parts() == 0) {
1134 ✗ RVD_.for_each_polyhedron(polyhedron_callback);
1135 } else {
1136 ✗ for(index_t t = 0; t < nb_parts(); t++) {
1137 part(t).RVD_.set_symbolic(RVD_.symbolic());
1138 part(t).RVD_.set_connected_components_priority(
1139 RVD_.connected_components_priority()
1140 );
1141 }
1142 ✗ spinlocks_.resize(delaunay_->nb_vertices());
1143 ✗ thread_mode_ = MT_POLYH;
1144 ✗ polyhedron_callback_ = &polyhedron_callback;
1145 polyhedron_callback_->set_spinlocks(&spinlocks_);
1146 // Note: callback begin()/end() is
1147 // called in for_each_polyhedron()
1148 ✗ parallel_for(
1149 0, nb_parts(),
1150 ✗ [this](index_t i) { run_thread(i); }
1151 );
1152 ✗ polyhedron_callback_->set_spinlocks(nullptr);
1153 }
1154 ✗ }
1155
1156 /********************************************************************/
1157
1158 /**
1159 * \brief Implementation class for explicitly constructing
1160 * a surfacic mesh that corresponds to the surfacic
1161 * restricted Voronoi diagram.
1162 * \details To be used as a template argument
1163 * to RVD::for_each_polygon(). The current Vornoi cell is
1164 * reported in facet region.
1165 * \tparam BUILDER a class that implements iterative mesh building,
1166 * e.g., MeshBuilder.
1167 */
1168 template <class BUILDER>
1169 class BuildRVD {
1170 public:
1171 /**
1172 * \brief Constructs a new BuildRVD.
1173 * \param[in] RVD_in the restricted Voronoi diagram
1174 * \param[in] builder the lesh builder
1175 */
1176 36 BuildRVD(
1177 const GenRestrictedVoronoiDiagram& RVD_in,
1178 BUILDER& builder
1179 ) :
1180 36 RVD(RVD_in),
1181 36 builder_(builder),
1182 36 current_facet_(NO_INDEX) {
1183 36 builder_.begin_surface();
1184 36 }
1185
1186 /**
1187 * \brief The destructor
1188 * \details Terminates the current facet
1189 * and the current surface.
1190 */
1191 ~BuildRVD() {
1192 if(current_facet_ != NO_INDEX) {
1193 builder_.end_reference_facet();
1194 }
1195 36 builder_.end_surface();
1196 ✗ }
1197
1198 /**
1199 * \brief The callback called for each restricted Voronoi cell.
1200 * \param[in] v index of current center vertex
1201 * \param[in] P current restricted Voronoi cell
1202 */
1203 420818 void operator() (
1204 index_t v,
1205 const typename GenRestrictedVoronoiDiagram::Polygon& P
1206 ) {
1207
2/2
✓ Branch 0 taken 60502 times.
✓ Branch 1 taken 149907 times.
 420818 index_t f = RVD.current_facet();
1208
2/2
✓ Branch 0 taken 60502 times.
✓ Branch 1 taken 149907 times.
 420818 if(f != current_facet_) {
1209 if(current_facet_ != NO_INDEX) {
1210 builder_.end_reference_facet();
1211 }
1212 121004 current_facet_ = f;
1213 builder_.begin_reference_facet(f);
1214 }
1215 420818 builder_.begin_facet(v);
1216
2/2
✓ Branch 0 taken 841011 times.
✓ Branch 1 taken 210409 times.
 4205680 for(index_t i = 0; i < P.nb_vertices(); i++) {
1217 const Vertex& ve = P.vertex(i);
1218 1682022 builder_.add_vertex_to_facet(ve.point(), ve.sym());
1219 }
1220 420818 builder_.end_facet();
1221 420818 }
1222
1223 private:
1224 const GenRestrictedVoronoiDiagram& RVD;
1225 BUILDER& builder_;
1226 index_t current_facet_;
1227 };
1228
1229 /**
1230 * \brief Implementation class for explicitly constructing
1231 * a volumetric mesh that corresponds to the volumetric
1232 * restricted Voronoi diagram.
1233 * \details To be used as a template argument
1234 * to RVD::for_each_volumetric_integration_simplex().
1235 * The current Voronoi cell is reported in tetrahedron region.
1236 * \note For the moment, vertices are duplicated (will be fixed
1237 * in a future version).
1238 */
1239 ✗ class BuildVolumetricRVD {
1240 public:
1241 /**
1242 * Constructs a new BuildVolumetricRVD.
1243 * \param[in] RVD the volumetric restricted Voronoi diagram
1244 * \param[in] dim dimension of the points (can be smaller
1245 * than actual dimension of the RVD).
1246 * \param[out] vertices coordinates of the generated vertices
1247 * \param[out] triangle_vertex_indices generated triangles, as
1248 * vertex indices triplets
1249 * \param[out] tet_vertex_indices generated tetrahedra, as
1250 * vertex indices 4-uples
1251 * \param[out] triangle_regions each generated triangle has
1252 * a region index, that corresponds to the index of the
1253 * Voronoi cell the triangle belongs to
1254 * \param[out] tet_regions each generated tetrahedron has
1255 * a region index, that corresponds to the index of the
1256 * Voronoi cell the tetrahedron belongs to
1257 * \param[in] cell_borders_only if true, only the surfacic
1258 * borders of the volumetric cells are saved in the mesh, else
1259 * volumetric cells are tetrahedralized.
1260 * \pre dim <= delaunay->dimension()
1261 */
1262 ✗ BuildVolumetricRVD(
1263 GenRestrictedVoronoiDiagram& RVD,
1264 coord_index_t dim,
1265 vector<double>& vertices,
1266 vector<index_t>& triangle_vertex_indices,
1267 vector<index_t>& tet_vertex_indices,
1268 vector<index_t>& triangle_regions,
1269 vector<index_t>& tet_regions,
1270 bool cell_borders_only
1271 ) :
1272 ✗ delaunay_(RVD.delaunay()),
1273 ✗ mesh_(RVD.mesh()),
1274 ✗ dim_(dim),
1275 ✗ vertices_(vertices),
1276 ✗ triangle_vertex_indices_(triangle_vertex_indices),
1277 ✗ tet_vertex_indices_(tet_vertex_indices),
1278 ✗ triangle_regions_(triangle_regions),
1279 ✗ tet_regions_(tet_regions),
1280 ✗ cell_borders_only_(cell_borders_only)
1281 {
1282 ✗ vertices_.clear();
1283 ✗ triangle_vertex_indices_.clear();
1284 ✗ tet_vertex_indices_.clear();
1285 ✗ triangle_regions_.clear();
1286 ✗ tet_regions_.clear();
1287
1288 // The first vertices are copied from Delaunay,
1289 // the other ones will be created during the traversal
1290 ✗ nb_vertices_ = delaunay_->nb_vertices();
1291 ✗ vertices_.resize(nb_vertices_ * dim);
1292 ✗ for(index_t v = 0; v < delaunay_->nb_vertices(); ++v) {
1293 ✗ for(coord_index_t c = 0; c < dim; ++c) {
1294 ✗ vertices_[v * dim + c] = delaunay_->vertex_ptr(v)[c];
1295 }
1296 }
1297 ✗ vertex_map_.set_first_vertex_index(nb_vertices_);
1298 ✗ }
1299
1300 /**
1301 * \brief The callback called for each integration simplex.
1302 * \param[in] v index of current center vertex
1303 * \param[in] v_adj (unused here) is the index of the Voronoi cell
1304 * adjacent to t accros facet (\p v1, \p v2, \p v3) or
1305 * NO_INDEX if it does not exists
1306 * \param[in] t (unused here) is the index of the current
1307 * tetrahedron
1308 * \param[in] t_adj (unused here) is the index of the
1309 * tetrahedron adjacent to t accros facet (\p v1, \p v2, \p v3)
1310 * or NO_INDEX if it does not exists
1311 * \param[in] v1 first vertex of current integration simplex
1312 * \param[in] v2 second vertex of current integration simplex
1313 * \param[in] v3 third vertex of current integration simplex
1314 */
1315 ✗ void operator() (
1316 index_t v, index_t v_adj,
1317 index_t t, index_t t_adj,
1318 const Vertex& v1, const Vertex& v2, const Vertex& v3
1319 ) {
1320 geo_argused(v_adj);
1321 geo_argused(t);
1322
1323 ✗ if(cell_borders_only_) {
1324 ✗ if(v > v_adj) {
1325 ✗ index_t iv2 = find_or_create_vertex(v, v1);
1326 ✗ index_t iv3 = find_or_create_vertex(v, v2);
1327 ✗ index_t iv4 = find_or_create_vertex(v, v3);
1328 ✗ triangle_vertex_indices_.push_back(iv4);
1329 ✗ triangle_vertex_indices_.push_back(iv3);
1330 ✗ triangle_vertex_indices_.push_back(iv2);
1331 ✗ triangle_regions_.push_back(v);
1332 }
1333 } else {
1334 ✗ index_t iv1 = v;
1335 ✗ index_t iv2 = find_or_create_vertex(v, v1);
1336 ✗ index_t iv3 = find_or_create_vertex(v, v2);
1337 ✗ index_t iv4 = find_or_create_vertex(v, v3);
1338
1339 // Triangle v1,v2,v3 is on border if there is
1340 // no adjacent seed and no adjacent tet.
1341 ✗ if(v_adj == NO_INDEX && t_adj == NO_INDEX) {
1342 ✗ triangle_vertex_indices_.push_back(iv4);
1343 ✗ triangle_vertex_indices_.push_back(iv3);
1344 ✗ triangle_vertex_indices_.push_back(iv2);
1345 ✗ triangle_regions_.push_back(v);
1346 }
1347
1348 ✗ tet_vertex_indices_.push_back(iv1);
1349 ✗ tet_vertex_indices_.push_back(iv2);
1350 ✗ tet_vertex_indices_.push_back(iv3);
1351 ✗ tet_vertex_indices_.push_back(iv4);
1352 ✗ tet_regions_.push_back(v);
1353 }
1354 ✗ }
1355
1356 /**
1357 * \brief The callback called for each tetrahedron
1358 * \param[in] v index of current center vertex
1359 * \param[in] v_adj (unused here) is the index of the Voronoi cell
1360 * adjacent to t accros facet (\p v1, \p v2, \p v3) or
1361 * NO_INDEX if it does not exists
1362 * \param[in] t (unused here) is the index of the current
1363 * tetrahedron
1364 * \param[in] t_adj (unused here) is the index of the
1365 * tetrahedron adjacent to t accros facet (\p v1, \p v2, \p v3)
1366 * or NO_INDEX if it does not exists
1367 * \param[in] v1 first vertex of current tetrahedron
1368 * \param[in] v2 second vertex of current tetrahedron
1369 * \param[in] v3 third vertex of current tetrahedron
1370 * \param[in] v4 fourth vertex of current tetrahedron
1371 */
1372 ✗ void operator() (
1373 index_t v, index_t v_adj,
1374 index_t t, index_t t_adj,
1375 const Vertex& v1, const Vertex& v2,
1376 const Vertex& v3, const Vertex& v4
1377 ) {
1378 geo_argused(v_adj);
1379 geo_argused(t);
1380 geo_argused(t_adj);
1381 ✗ index_t iv1 = vertices_.size() / dim_;
1382 ✗ for(index_t c = 0; c < dim_; ++c) {
1383 ✗ vertices_.push_back(v1.point()[c]);
1384 }
1385 ✗ index_t iv2 = vertices_.size() / dim_;
1386 ✗ for(index_t c = 0; c < dim_; ++c) {
1387 ✗ vertices_.push_back(v2.point()[c]);
1388 }
1389 ✗ index_t iv3 = vertices_.size() / dim_;
1390 ✗ for(index_t c = 0; c < dim_; ++c) {
1391 ✗ vertices_.push_back(v3.point()[c]);
1392 }
1393 ✗ index_t iv4 = vertices_.size() / dim_;
1394 ✗ for(index_t c = 0; c < dim_; ++c) {
1395 ✗ vertices_.push_back(v4.point()[c]);
1396 }
1397 ✗ tet_vertex_indices_.push_back(iv1);
1398 ✗ tet_vertex_indices_.push_back(iv2);
1399 ✗ tet_vertex_indices_.push_back(iv3);
1400 ✗ tet_vertex_indices_.push_back(iv4);
1401 ✗ tet_regions_.push_back(v);
1402 ✗ }
1403
1404 protected:
1405 /**
1406 * \brief Retrieves the index of a vertex given its symbolic
1407 * representation.
1408 * \param[in] center_vertex_id index of current Voronoi seed
1409 * \param[in] v symbolic and geometric representation of the vertex
1410 * \return the index of the vertex
1411 */
1412 ✗ index_t find_or_create_vertex(
1413 index_t center_vertex_id, const Vertex& v
1414 ) {
1415 ✗ index_t result = vertex_map_.find_or_create_vertex(
1416 center_vertex_id, v.sym()
1417 );
1418 ✗ if(result >= nb_vertices_) {
1419 ✗ geo_assert(result == nb_vertices_);
1420 ✗ nb_vertices_ = result + 1;
1421 ✗ for(coord_index_t c = 0; c < dim_; ++c) {
1422 ✗ vertices_.push_back(v.point()[c]);
1423 }
1424 }
1425 ✗ return result;
1426 }
1427
1428 private:
1429 const Delaunay* delaunay_;
1430 const Mesh* mesh_;
1431 coord_index_t dim_;
1432 vector<double>& vertices_;
1433 vector<index_t>& triangle_vertex_indices_;
1434 vector<index_t>& tet_vertex_indices_;
1435 vector<index_t>& triangle_regions_;
1436 vector<index_t>& tet_regions_;
1437 RVDVertexMap vertex_map_;
1438 index_t nb_vertices_;
1439 bool cell_borders_only_;
1440 };
1441
1442
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
 72 void compute_RVD(
1443 Mesh& M, coord_index_t dim, bool cell_borders_only,
1444 bool integration_simplices
1445 ) override {
1446 bool sym = RVD_.symbolic();
1447 RVD_.set_symbolic(true);
1448
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
 72 if(volumetric_) {
1449 ✗ if(dim == 0) {
1450 dim = dimension();
1451 }
1452 vector<double> vertices;
1453 vector<index_t> triangle_vertices;
1454 vector<index_t> tet_vertices;
1455 vector<index_t> triangle_regions;
1456 vector<index_t> tet_regions;
1457 ✗ if(cell_borders_only) {
1458 ✗ RVD_.for_each_volumetric_integration_simplex(
1459 ✗ BuildVolumetricRVD(
1460 RVD_, dim,
1461 vertices,
1462 triangle_vertices,
1463 tet_vertices,
1464 triangle_regions,
1465 tet_regions,
1466 cell_borders_only
1467 ),
1468 false, // Do not visit inner tetrahedra.
1469 true // Ensure that polygonal facets are triangulated
1470 // coherently.
1471 );
1472 } else {
1473 ✗ if(integration_simplices) {
1474 ✗ RVD_.for_each_volumetric_integration_simplex(
1475 ✗ BuildVolumetricRVD(
1476 RVD_, dim,
1477 vertices,
1478 triangle_vertices,
1479 tet_vertices,
1480 triangle_regions,
1481 tet_regions,
1482 cell_borders_only
1483 ),
1484 false, // Do not visit inner tetrahedra.
1485 true // Ensure that polygonal facets are
1486 // triangulated coherently.
1487 );
1488 } else {
1489 ✗ RVD_.for_each_tetrahedron(
1490 ✗ BuildVolumetricRVD(
1491 RVD_, dim,
1492 vertices,
1493 triangle_vertices,
1494 tet_vertices,
1495 triangle_regions,
1496 tet_regions,
1497 cell_borders_only
1498 )
1499 );
1500 }
1501 }
1502
1503 ✗ M.clear(true); // keep attributes
1504
1505 ✗ M.vertices.assign_points(vertices,dim,true);
1506 ✗ M.facets.assign_triangle_mesh(triangle_vertices, true);
1507 ✗ M.cells.assign_tet_mesh(tet_vertices, true);
1508
1509 // TODO: use Attribute::assign(vector, steal_args)
1510 // when it is there...
1511
1512 ✗ if(M.facets.nb() != 0) {
1513 Attribute<index_t> facet_region_attr(
1514 ✗ M.facets.attributes(), "region"
1515 );
1516 ✗ for(index_t f=0; f<M.facets.nb(); ++f) {
1517 ✗ facet_region_attr[f] = triangle_regions[f];
1518 }
1519 }
1520
1521 ✗ if(M.cells.nb() != 0) {
1522 Attribute<index_t> cell_region_attr(
1523 ✗ M.cells.attributes(), "region"
1524 );
1525 ✗ for(index_t c=0; c<M.cells.nb(); ++c) {
1526 ✗ cell_region_attr[c] = tet_regions[c];
1527 }
1528 }
1529
1530 ✗ if(cell_borders_only) {
1531 ✗ mesh_repair(M, MESH_REPAIR_TOPOLOGY);
1532 } else {
1533 ✗ M.facets.connect();
1534 }
1535 } else {
1536 72 RVDMeshBuilder builder(
1537 &M, mesh_, delaunay_
1538 );
1539 if(dim != 0) {
1540 builder.set_dimension(dim);
1541 }
1542
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 72 RVD_.for_each_polygon(
1543 72 BuildRVD<RVDMeshBuilder>(RVD_, builder)
1544 );
1545 72 }
1546 RVD_.set_symbolic(sym);
1547
1/2
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
 72 M.show_stats("RVD");
1548 72 }
1549
1550 /********************************************************************/
1551
1552 ✗ void compute_RVC(
1553 index_t i,
1554 Mesh& M,
1555 Mesh& result,
1556 bool copy_symbolic_info
1557 ) override {
1558 ✗ Mesh* tmp_mesh = mesh_;
1559 ✗ mesh_ = &M;
1560 RVD_.set_mesh(&M);
1561 typename GenRestrictedVoronoiDiagram::Polyhedron Cell(dimension());
1562 ✗ Cell.initialize_from_surface_mesh(&M, RVD_.symbolic());
1563 ✗ RVD_.intersect_cell_cell(i, Cell);
1564 ✗ Cell.convert_to_mesh(&result, copy_symbolic_info);
1565 ✗ mesh_ = tmp_mesh;
1566 RVD_.set_mesh(tmp_mesh);
1567 ✗ }
1568
1569 /********************************************************************/
1570
1571 ✗ void for_each_polygon(
1572 GEO::RVDPolygonCallback& callback,
1573 bool symbolic,
1574 bool connected_comp_priority,
1575 bool parallel
1576 ) override {
1577 bool sym_backup = RVD_.symbolic();
1578 RVD_.set_symbolic(symbolic);
1579 RVD_.set_connected_components_priority(connected_comp_priority);
1580 ✗ callback.begin();
1581 ✗ if(parallel) {
1582 ✗ compute_with_polygon_callback(callback);
1583 } else {
1584 ✗ PolygonCallbackAction action(RVD_,callback);
1585 RVD_.for_each_polygon(action);
1586 }
1587 ✗ callback.end();
1588 RVD_.set_symbolic(sym_backup);
1589 RVD_.set_connected_components_priority(false);
1590 ✗ }
1591
1592 /********************************************************************/
1593
1594
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 8 void for_each_polyhedron(
1595 GEO::RVDPolyhedronCallback& callback,
1596 bool symbolic,
1597 bool connected_comp_priority,
1598 bool parallel
1599 ) override {
1600 bool sym_backup = RVD_.symbolic();
1601 RVD_.set_symbolic(symbolic);
1602 RVD_.set_connected_components_priority(connected_comp_priority);
1603 callback.set_dimension(RVD_.mesh()->vertices.dimension());
1604 8 callback.begin();
1605
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 8 if(parallel) {
1606 ✗ compute_with_polyhedron_callback(callback);
1607 } else {
1608 8 RVD_.for_each_polyhedron(callback);
1609 }
1610 8 callback.end();
1611 RVD_.set_symbolic(sym_backup);
1612 RVD_.set_connected_components_priority(false);
1613 8 }
1614
1615 /********************************************************************/
1616
1617 /**
1618 * \brief Does the actual computation for a specific part
1619 * in multithread mode.
1620 * \param[in] t the index of the part.
1621 * \pre \p t < nb_parts()
1622 */
1623 3200 void run_thread(index_t t) {
1624 ✗ geo_assert(t < nb_parts());
1625 thisclass& T = part(t);
1626 3200 switch(thread_mode_) {
1627 3200 case MT_LLOYD:
1628 {
1629
2/2
✓ Branch 0 taken 640 times.
✓ Branch 1 taken 960 times.
 3200 T.compute_centroids(arg_vectors_, arg_scalars_);
1630 } break;
1631 ✗ case MT_NEWTON:
1632 {
1633 ✗ T.compute_CVT_func_grad(T.funcval_, arg_vectors_);
1634 } break;
1635 ✗ case MT_INT_SMPLX:
1636 {
1637 ✗ T.compute_integration_simplex_func_grad(
1638 ✗ T.funcval_, arg_vectors_, simplex_func_
1639 );
1640 ✗ } break;
1641 ✗ case MT_POLYG:
1642 {
1643 ✗ T.compute_with_polygon_callback(
1644 ✗ *polygon_callback_
1645 );
1646 ✗ } break;
1647 ✗ case MT_POLYH:
1648 {
1649 ✗ T.compute_with_polyhedron_callback(
1650 ✗ *polyhedron_callback_
1651 );
1652 ✗ } break;
1653 case MT_NONE:
1654 ✗ geo_assert_not_reached;
1655 }
1656 3200 }
1657
1658 12 bool compute_initial_sampling_on_surface(
1659 double* p, index_t nb_points, bool verbose
1660 ) override {
1661
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 12 geo_assert(mesh_->facets.are_simplices());
1662
1663 // We do that here, since this triggers partitioning,
1664 // that improves data locality. Then data locality is
1665 // inherited by the generated points.
1666 12 create_threads();
1667
1668
3/6
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 12 if(verbose && facets_begin_ == NO_INDEX && facets_end_ == NO_INDEX) {
1669
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 24 Logger::out("RVD")
1670 << "Computing initial sampling on surface, using dimension="
1671
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 << index_t(dimension_) << std::endl;
1672 }
1673
1674 12 return mesh_generate_random_samples_on_surface<DIM>(
1675 12 *mesh_, p, nb_points, vertex_weight_, facets_begin_, facets_end_
1676 12 );
1677 }
1678
1679 8 bool compute_initial_sampling_in_volume(
1680 double* p, index_t nb_points, bool verbose
1681 ) override {
1682
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 8 geo_assert(mesh_->cells.nb() != 0);
1683
1684 // We do that here, since this triggers partitioning,
1685 // that improves data locality. Then data locality is
1686 // inherited by the generated points.
1687 8 create_threads();
1688
1689
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 8 if(verbose && tets_begin_ == NO_INDEX && tets_end_ == NO_INDEX) {
1690 ✗ Logger::out("RVD")
1691 << "Computing initial sampling in volume, using dimension="
1692 ✗ << index_t(dimension_) << std::endl;
1693 }
1694
1695 8 return mesh_generate_random_samples_in_volume<DIM>(
1696 8 *mesh_, p, nb_points, vertex_weight_, tets_begin_, tets_end_
1697 8 );
1698 }
1699
1700 /**
1701 * \brief Creates the data structures for fast projection.
1702 * \details It decomposes the surface into triangles, and
1703 * stores the vertices in a KdTree.
1704 */
1705 ✗ void prepare_projection() {
1706 ✗ if(!mesh_vertices_.is_null()) {
1707 ✗ return;
1708 }
1709
1710 // Step 1: get triangles
1711 ✗ for(index_t f = 0; f < mesh_->facets.nb(); f++) {
1712 index_t i = mesh_->facets.corners_begin(f);
1713 ✗ for(
1714 ✗ index_t j = i + 1;
1715 ✗ j + 1 < mesh_->facets.corners_end(f); j++
1716 ) {
1717 ✗ triangles_.push_back(mesh_->facet_corners.vertex(i));
1718 ✗ triangles_.push_back(mesh_->facet_corners.vertex(j));
1719 ✗ triangles_.push_back(mesh_->facet_corners.vertex(j + 1));
1720 }
1721 }
1722 ✗ nb_triangles_ = index_t(triangles_.size() / 3);
1723
1724 // Step 2: get vertices stars
1725 // Step 2.1: get one-ring neighborhood
1726 vector<vector<index_t> > stars2(mesh_->vertices.nb());
1727 ✗ for(index_t t = 0; t < nb_triangles_; t++) {
1728 ✗ stars2[triangles_[3 * t]].push_back(t);
1729 ✗ stars2[triangles_[3 * t + 1]].push_back(t);
1730 ✗ stars2[triangles_[3 * t + 2]].push_back(t);
1731 }
1732
1733 // Step 2.2: get two-ring neighborhood
1734 ✗ stars_.resize(mesh_->vertices.nb());
1735 ✗ for(index_t i = 0; i < stars2.size(); i++) {
1736 vector<index_t> Ni;
1737 ✗ for(index_t j = 0; j < stars2[i].size(); j++) {
1738 ✗ index_t t = stars2[i][j];
1739 ✗ for(index_t iv = 0; iv < 3; iv++) {
1740 ✗ index_t v = triangles_[3 * t + iv];
1741 ✗ if(v != i) {
1742 Ni.push_back(v);
1743 }
1744 }
1745 }
1746 ✗ sort_unique(Ni);
1747 ✗ for(index_t j = 0; j < Ni.size(); j++) {
1748 ✗ index_t k = Ni[j];
1749 ✗ stars_[i].insert(
1750 stars_[i].end(), stars2[k].begin(), stars2[k].end()
1751 );
1752 }
1753 ✗ sort_unique(stars_[i]);
1754 }
1755
1756 // Step 3: create search structure
1757 ✗ mesh_vertices_ = Delaunay::create(dimension_, "NN");
1758 ✗ index_t nb_vertices = mesh_->vertices.nb();
1759
1760 // TODO: BUG !! mesh_vertices_ keeps a ref. to mesh_vertices
1761 // that is destroyed when leaving this function.
1762 ✗ vector<double> mesh_vertices(nb_vertices * dimension_);
1763 ✗ for(index_t i = 0; i < nb_vertices; i++) {
1764 ✗ for(index_t coord = 0; coord < dimension_; coord++) {
1765 ✗ mesh_vertices[i * dimension_ + coord] =
1766 ✗ mesh_->vertices.point_ptr(i)[coord];
1767 }
1768 }
1769 ✗ mesh_vertices_->set_vertices(nb_vertices, mesh_vertices.data());
1770 }
1771
1772 ✗ void project_points_on_surface(
1773 index_t nb_points, double* points, vec3* nearest, bool do_project
1774 ) override {
1775
1776 ✗ prepare_projection();
1777
1778 ✗ if(use_exact_projection_) {
1779 ✗ for(index_t p = 0; p < nb_points; p++) {
1780 ✗ Point P(points + p * dimension_);
1781 double d2 = Numeric::max_float64();
1782 ✗ for(index_t t = 0; t < nb_triangles_; t++) {
1783 ✗ const Point& p1 = mesh_vertex(triangles_[3 * t]);
1784 ✗ const Point& p2 = mesh_vertex(triangles_[3 * t + 1]);
1785 ✗ const Point& p3 = mesh_vertex(triangles_[3 * t + 2]);
1786
1787 double l1, l2, l3;
1788 Point nearestP;
1789
1790 ✗ double cur_d2 = Geom::point_triangle_squared_distance(
1791 P, p1, p2, p3, nearestP, l1, l2, l3
1792 );
1793
1794 ✗ if(cur_d2 < d2) {
1795 d2 = cur_d2;
1796 const vec3& p1_R3 =
1797 R3_embedding(triangles_[3 * t]);
1798 const vec3& p2_R3 =
1799 R3_embedding(triangles_[3 * t + 1]);
1800 const vec3& p3_R3 =
1801 R3_embedding(triangles_[3 * t + 2]);
1802 ✗ nearest[p] = l1 * p1_R3 + l2 * p2_R3 + l3 * p3_R3;
1803 ✗ if(do_project) {
1804 ✗ for(coord_index_t
1805 ✗ coord = 0; coord < dimension_; coord++) {
1806 ✗ (points + p * dimension_)[coord] =
1807 ✗ nearestP[coord];
1808 }
1809 }
1810 }
1811 }
1812 }
1813 return;
1814 }
1815
1816 // find nearest point on surface in star of nearest vertex
1817 ✗ for(index_t p = 0; p < nb_points; p++) {
1818 ✗ Point P(points + p * dimension_);
1819 ✗ index_t v = mesh_vertices_->nearest_vertex(
1820 ✗ points + p * dimension_
1821 );
1822 double d2 = Numeric::max_float64();
1823 ✗ nearest[p] = R3_embedding(v);
1824 ✗ for(index_t i = 0; i < stars_[v].size(); i++) {
1825 ✗ index_t t = stars_[v][i];
1826 ✗ const Point& p1 = mesh_vertex(triangles_[3 * t]);
1827 ✗ const Point& p2 = mesh_vertex(triangles_[3 * t + 1]);
1828 ✗ const Point& p3 = mesh_vertex(triangles_[3 * t + 2]);
1829
1830 double l1, l2, l3;
1831 Point nearestP;
1832 ✗ double cur_d2 = Geom::point_triangle_squared_distance(
1833 P, p1, p2, p3, nearestP, l1, l2, l3
1834 );
1835 ✗ if(cur_d2 < d2) {
1836 d2 = cur_d2;
1837 const vec3& p1_R3 = R3_embedding(triangles_[3 * t]);
1838 const vec3& p2_R3 = R3_embedding(triangles_[3 * t + 1]);
1839 const vec3& p3_R3 = R3_embedding(triangles_[3 * t + 2]);
1840 ✗ nearest[p] = l1 * p1_R3 + l2 * p2_R3 + l3 * p3_R3;
1841 ✗ if(do_project) {
1842 ✗ for(coord_index_t coord = 0;
1843 ✗ coord < dimension_; coord++
1844 ) {
1845 ✗ (points + p * dimension_)[coord] =
1846 ✗ nearestP[coord];
1847 }
1848 }
1849 }
1850 }
1851 }
1852 }
1853
1854 /********************************************************************/
1855
1856 /**
1857 * \brief Implementation class for computing the restricted Delaunay
1858 * triangulation.
1859 * \details To be used as a template argument
1860 * to RVD::for_each_primal_triangle().
1861 */
1862 class GetPrimalTriangles {
1863 public:
1864 /**
1865 * \brief Creates a new GetPrimalTriangles.
1866 * \param[out] triangles where to store the triangles
1867 */
1868 ✗ GetPrimalTriangles(
1869 vector<index_t>& triangles
1870 ) :
1871 ✗ triangles_(triangles) {
1872 }
1873
1874 /**
1875 * \brief The callback called for each primal triangle.
1876 * \param[in] v1 index of the first vertex
1877 * \param[in] v2 index of the second vertex
1878 * \param[in] v3 index of the third vertex
1879 */
1880 ✗ void operator() (index_t v1, index_t v2, index_t v3) {
1881 ✗ triangles_.push_back(v1);
1882 ✗ triangles_.push_back(v2);
1883 ✗ triangles_.push_back(v3);
1884 ✗ }
1885
1886 private:
1887 vector<index_t>& triangles_;
1888 };
1889
1890 /**
1891 * \brief Implementation class for computing the restricted Delaunay
1892 * triangulation of the connected components.
1893 *
1894 * \details The difference with GetPrimalTriangles is that when a
1895 * restricted Voronoi cell has multiple connected components,
1896 * more triangles are generated to account for the topology.
1897 * To be used as a template argument to RVD::for_each_polygon().
1898 * The RestrictedVoronoiDiagram needs to be in connected-components
1899 * priority mode.
1900 */
1901 class GetConnectedComponentsPrimalTriangles {
1902 public:
1903 /** Internal representation of the polygons. */
1904 typedef typename GenRestrictedVoronoiDiagram::Polygon Polygon;
1905
1906 /** Internal representation of the vertices. */
1907 typedef typename GenRestrictedVoronoiDiagram::Vertex Vertex;
1908
1909 static constexpr index_t UNINITIALIZED = index_t(-1);
1910 static constexpr index_t MULTI_COMP = index_t(-2);
1911 static constexpr index_t ON_BORDER = index_t(-3);
1912
1913 /**
1914 * \brief Constructs a new GetConnectedComponentsPrimalTriangles.
1915 * \param[in] RVD the restricted Voronoi diagram
1916 * \param[out] triangles where to store the triangles
1917 * \param[out] vertices where to store the vertices
1918 * \param[in] dimension dimension of the restricted Voronoi diagram
1919 * \param[in] mode a combination of constants defined in RDTMode
1920 * \param[in] seed_is_locked specifies for each seed whether it
1921 * can be moved (to RVC centroid or projected on surface). If
1922 * left uninitialized, all the seeds can be moved.
1923 * \param[in] AABB an axis-aligned bounding box tree defined on
1924 * the input surface. It is used if one of (RDT_SELECT_NEAREST,
1925 * RDT_PROJECT_ON_SURFACE) is set in \p mode. If needed and not
1926 * specified, then a new one is created locally.
1927 */
1928 12 GetConnectedComponentsPrimalTriangles(
1929 const GenRestrictedVoronoiDiagram& RVD,
1930 vector<index_t>& triangles,
1931 vector<double>& vertices,
1932 coord_index_t dimension,
1933 RDTMode mode,
1934 const std::vector<bool>& seed_is_locked,
1935 MeshFacetsAABB* AABB = nullptr
1936 ) :
1937 12 RVD_(RVD),
1938 12 dimension_(dimension),
1939 12 triangles_(triangles),
1940 12 vertices_(vertices),
1941 12 m_(0.0),
1942 12 cur_seed_(NO_INDEX),
1943 12 cur_vertex_(0),
1944 12 use_RVC_centroids_((mode & RDT_RVC_CENTROIDS) != 0),
1945 12 select_nearest_((mode & RDT_SELECT_NEAREST) != 0),
1946 12 project_on_surface_((mode & RDT_PROJECT_ON_SURFACE) != 0),
1947 12 seed_is_locked_(seed_is_locked),
1948
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 12 prefer_seeds_((mode & RDT_PREFER_SEEDS) != 0),
1949 12 AABB_(AABB)
1950 {
1951
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 12 if(prefer_seeds_) {
1952
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 seed_to_vertex_.assign(
1953 RVD.delaunay()->nb_vertices(), UNINITIALIZED
1954 );
1955 }
1956 12 }
1957
1958 /**
1959 * \brief The callback called for each restricted Voronoi cell.
1960 * \param[in] s1 index of current center vertex
1961 * \param[in] P current restricted Voronoi cell
1962 */
1963 150052 void operator() (index_t s1, const Polygon& P) {
1964
2/2
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 45026 times.
 150052 if(RVD_.connected_component_changed()) {
1965
2/2
✓ Branch 0 taken 29994 times.
✓ Branch 1 taken 6 times.
 60000 if(cur_seed_ != NO_INDEX) {
1966 59988 end_connected_component();
1967 }
1968 begin_connected_component(s1);
1969 }
1970
1971
4/4
✓ Branch 0 taken 56283 times.
✓ Branch 1 taken 18743 times.
✓ Branch 2 taken 55994 times.
✓ Branch 3 taken 289 times.
 150052 if(prefer_seeds_ && !component_on_border_) {
1972 // NOTE: there is one vertex shift between
1973 // adjacent facet and adjacent seed, there
1974 // must be something wrong in the way the
1975 // combinatorial information is initialized,
1976 // to be checked !!! (should be the i index
1977 // for both)
1978
2/2
✓ Branch 0 taken 256242 times.
✓ Branch 1 taken 55800 times.
 624084 for(index_t i=0; i<P.nb_vertices(); ++i) {
1979
2/2
✓ Branch 0 taken 177030 times.
✓ Branch 1 taken 79212 times.
 512484 index_t j = (i+1) % P.nb_vertices();
1980 if(
1981
4/4
✓ Branch 0 taken 177030 times.
✓ Branch 1 taken 79212 times.
✓ Branch 2 taken 194 times.
✓ Branch 3 taken 176836 times.
 512484 P.vertex(i).adjacent_facet() == -1 &&
1982 P.vertex(j).adjacent_seed() == -1
1983 ) {
1984 388 component_on_border_ = true;
1985 388 break;
1986 }
1987 }
1988 }
1989
1990 // Accumulate mass and barycenter
1991 150052 index_t vbase = cur_vertex_ * dimension_;
1992
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 277084 times.
✓ Branch 2 taken 202058 times.
✓ Branch 3 taken 75026 times.
 554168 for(index_t i = 1; i + 1 < P.nb_vertices(); ++i) {
1993 404116 double cur_m = Geom::triangle_area(
1994 P.vertex(0).point(),
1995 P.vertex(i).point(),
1996 P.vertex(i + 1).point(), dimension_
1997 );
1998
2/2
✓ Branch 0 taken 1041264 times.
✓ Branch 1 taken 202058 times.
 2486644 for(coord_index_t c = 0; c < dimension_; ++c) {
1999 2082528 vertices_[vbase + c] += cur_m / 3.0 * (
2000 2082528 P.vertex(0).point()[c] +
2001 2082528 P.vertex(i).point()[c] +
2002 2082528 P.vertex(i + 1).point()[c]
2003 );
2004 }
2005 404116 m_ += cur_m;
2006 }
2007
2008 // Detect Voronoi vertices and generate
2009 // triangles.
2010 // Note: they can be generated several times,
2011 // since we cannot know in advance whether
2012 // the other instances of the Voronoi vertex
2013 // were finalized or not (i.e. have their
2014 // three vertices ready).
2015 // Duplicate triangles are then filtered-out
2016 // by client code.
2017
2/2
✓ Branch 0 taken 352110 times.
✓ Branch 1 taken 75026 times.
 1558492 for(index_t i = 0; i < P.nb_vertices(); ++i) {
2018 const Vertex& V = P.vertex(i);
2019
2/2
✓ Branch 0 taken 179412 times.
✓ Branch 1 taken 172698 times.
 704220 if(V.sym().nb_bisectors() == 2) {
2020 index_t s2 = V.sym().bisector(0);
2021 index_t s3 = V.sym().bisector(1);
2022 index_t f = V.sym().boundary_facet(0);
2023
2024
2/2
✓ Branch 0 taken 89012 times.
✓ Branch 1 taken 90400 times.
 358824 index_t v1 = RVD_.current_connected_component();
2025 358824 index_t v2 = index_t(
2026 RVD_.get_facet_seed_connected_component(f,s2)
2027 );
2028 358824 index_t v3 = index_t(
2029 RVD_.get_facet_seed_connected_component(f,s3)
2030 );
2031
2032
4/4
✓ Branch 0 taken 89012 times.
✓ Branch 1 taken 90400 times.
✓ Branch 2 taken 59804 times.
✓ Branch 3 taken 29208 times.
 358824 if(v2 != NO_INDEX && v3 != NO_INDEX) {
2033
2/2
✓ Branch 0 taken 59798 times.
✓ Branch 1 taken 6 times.
 119608 triangles_.push_back(v1);
2034
2/2
✓ Branch 0 taken 59756 times.
✓ Branch 1 taken 48 times.
 119608 triangles_.push_back(v2);
2035
2/2
✓ Branch 0 taken 59762 times.
✓ Branch 1 taken 42 times.
 119608 triangles_.push_back(v3);
2036 }
2037 }
2038 }
2039 150052 }
2040
2041 /**
2042 * \brief The destructor
2043 */
2044 12 ~GetConnectedComponentsPrimalTriangles() {
2045
2046 bool owns_AABB = false;
2047
2/4
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
 12 if(select_nearest_ || project_on_surface_) {
2048 ✗ if(AABB_ == nullptr) {
2049 // Construct an axis-aligned bounding box tree,
2050 // do not reorder the mesh (needs to be pre-reordered)
2051 ✗ AABB_ = new MeshFacetsAABB(
2052 ✗ *const_cast<Mesh*>(RVD_.mesh()), AABB_NOREORDER
2053 );
2054 owns_AABB = true;
2055 }
2056 }
2057
2058
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 12 if(cur_seed_ != NO_INDEX) {
2059 12 end_connected_component();
2060 }
2061
2062
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 12 if(prefer_seeds_) {
2063
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 8 if(select_nearest_) {
2064 ✗ for(index_t s=0; s<seed_to_vertex_.size(); ++s) {
2065 if(
2066 ✗ seed_to_vertex_[s] != MULTI_COMP &&
2067 seed_to_vertex_[s] != UNINITIALIZED &&
2068 seed_to_vertex_[s] != ON_BORDER
2069 ) {
2070 ✗ index_t vbase = seed_to_vertex_[s] * dimension_;
2071
2072 const double* seed_ptr =
2073 ✗ RVD_.delaunay()->vertex_ptr(s);
2074
2075
2076 // At this step, vertex_ptr contains
2077 // the centroid of the connected component
2078 // of the RVC, we now determine whether it
2079 // should be replaced by the
2080 // seed (or by a projection onto the surface).
2081
2082 ✗ const double* vertex_ptr = &(vertices_[vbase]);
2083
2084 // If the seed is nearer to the surface
2085 // than the
2086 // centroid of the connected component of the
2087 // restricted Voronoi cell, then use the seed.
2088
2089 double seed_dist;
2090 vec3 seed_projection;
2091 double vertex_dist;
2092 vec3 vertex_projection;
2093
2094 ✗ AABB_->nearest_facet(
2095 ✗ vec3(seed_ptr), seed_projection, seed_dist
2096 );
2097 ✗ AABB_->nearest_facet(
2098 ✗ vec3(vertex_ptr),
2099 vertex_projection, vertex_dist
2100 );
2101
2102 ✗ if(seed_dist < vertex_dist) {
2103 ✗ if(project_on_surface_) {
2104 ✗ for(
2105 coord_index_t c = 0;
2106 ✗ c < dimension_; ++c
2107 ) {
2108 ✗ vertices_[vbase + c] =
2109 seed_projection[c];
2110 }
2111 } else {
2112 ✗ for(coord_index_t c = 0;
2113 ✗ c < dimension_; ++c
2114 ) {
2115 ✗ vertices_[vbase + c] = seed_ptr[c];
2116 }
2117 }
2118 } else {
2119 ✗ if(project_on_surface_) {
2120 ✗ for(
2121 coord_index_t c = 0;
2122 ✗ c < dimension_; ++c
2123 ) {
2124 ✗ vertices_[vbase + c] =
2125 vertex_projection[c];
2126 }
2127 }
2128 }
2129 }
2130 }
2131 } else {
2132
2133 // Current mode: prefer seeds and not select nearest
2134 // Replace all points with the seeds (provided that
2135 // they do not correspond to multiple
2136 // connected components).
2137
2138
2/2
✓ Branch 0 taken 20000 times.
✓ Branch 1 taken 4 times.
 80008 for(index_t s=0; s<seed_to_vertex_.size(); ++s) {
2139 if(
2140
2/2
✓ Branch 0 taken 19806 times.
✓ Branch 1 taken 194 times.
 40000 seed_to_vertex_[s] != MULTI_COMP &&
2141 seed_to_vertex_[s] != UNINITIALIZED &&
2142 seed_to_vertex_[s] != ON_BORDER
2143 ) {
2144 39612 index_t vbase = seed_to_vertex_[s] * dimension_;
2145
2146 const double* seed_ptr =
2147 39612 RVD_.delaunay()->vertex_ptr(s);
2148
2149
2/2
✓ Branch 0 taken 118836 times.
✓ Branch 1 taken 19806 times.
 277284 for(coord_index_t c = 0; c < dimension_; ++c) {
2150 237672 vertices_[vbase + c] = seed_ptr[c];
2151 }
2152 }
2153 }
2154 }
2155 }
2156
2157 12 if(
2158
4/6
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
 12 (!prefer_seeds_ || !select_nearest_) && project_on_surface_
2159 ) {
2160 ✗ for(index_t v=0; v<vertices_.size()/3; ++v) {
2161 vec3 p(
2162 vertices_[3*v], vertices_[3*v+1], vertices_[3*v+2]
2163 );
2164 vec3 q;
2165 double sq_dist;
2166 ✗ AABB_->nearest_facet(p,q,sq_dist);
2167 ✗ vertices_[3*v ] = q.x;
2168 ✗ vertices_[3*v+1] = q.y;
2169 ✗ vertices_[3*v+2] = q.z;
2170 }
2171 }
2172
2173
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 12 if(owns_AABB) {
2174 ✗ delete AABB_;
2175 AABB_ = nullptr;
2176 }
2177 12 }
2178
2179 protected:
2180 /**
2181 * \brief Tests whether a given seed is locked.
2182 */
2183 bool seed_is_locked(index_t s) {
2184 return
2185
1/28
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 20000 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
 20000 seed_is_locked_.size() > 0 &&
2186 seed_is_locked_[s]
2187 ;
2188 }
2189
2190 /**
2191 * \brief Starts a new connected component.
2192 * \param[in] s the seed the connected component
2193 * is associated with.
2194 */
2195 void begin_connected_component(index_t s) {
2196 30000 cur_seed_ = s;
2197
4/14
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 30000 times.
✓ Branch 3 taken 10000 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 120000 times.
✓ Branch 7 taken 20000 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 180000 for(coord_index_t c = 0; c < dimension_; ++c) {
2198 150000 vertices_.push_back(0.0);
2199 }
2200 30000 m_ = 0.0;
2201 30000 component_on_border_ = false;
2202 30000 }
2203
2204 /**
2205 * \brief Terminates the current connected component.
2206 */
2207 60000 void end_connected_component() {
2208
2209 if(
2210 60000 !use_RVC_centroids_ ||
2211
3/4
✓ Branch 0 taken 20000 times.
✓ Branch 1 taken 10000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 20000 times.
 60000 seed_is_locked(cur_seed_) ||
2212
2/2
✓ Branch 0 taken 19806 times.
✓ Branch 1 taken 194 times.
 40000 component_on_border_
2213 ) {
2214 // Copy seed
2215 20388 index_t vbase = cur_vertex_ * dimension_;
2216 const double* seed_ptr =
2217 20388 RVD_.delaunay()->vertex_ptr(cur_seed_);
2218
2/2
✓ Branch 0 taken 31164 times.
✓ Branch 1 taken 10194 times.
 82716 for(coord_index_t c = 0; c < dimension_; ++c) {
2219 62328 vertices_[vbase + c] = seed_ptr[c];
2220 }
2221 } else {
2222 // Use restricted Voronoi
2223 // cell component's centroid.
2224
1/2
✓ Branch 0 taken 19806 times.
✗ Branch 1 not taken.
 39612 double scal = (m_ < 1e-30 ? 0.0 : 1.0 / m_);
2225 39612 index_t vbase = cur_vertex_ * dimension_;
2226
2/2
✓ Branch 0 taken 118836 times.
✓ Branch 1 taken 19806 times.
 277284 for(coord_index_t c = 0; c < dimension_; ++c) {
2227 237672 vertices_[vbase + c] *= scal;
2228 }
2229 }
2230
2/2
✓ Branch 0 taken 20000 times.
✓ Branch 1 taken 10000 times.
 60000 if(prefer_seeds_) {
2231
2/2
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 19806 times.
 40000 if(component_on_border_) {
2232 388 seed_to_vertex_[cur_seed_] = ON_BORDER;
2233 }
2234
2/3
✓ Branch 0 taken 19806 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 194 times.
 40000 switch(seed_to_vertex_[cur_seed_]) {
2235 39612 case UNINITIALIZED:
2236 39612 seed_to_vertex_[cur_seed_] = cur_vertex_;
2237 39612 break;
2238 case ON_BORDER:
2239 break;
2240 ✗ default:
2241 ✗ seed_to_vertex_[cur_seed_] = MULTI_COMP;
2242 ✗ break;
2243 }
2244 }
2245 60000 ++cur_vertex_;
2246 60000 }
2247
2248 private:
2249 const GenRestrictedVoronoiDiagram& RVD_;
2250 coord_index_t dimension_;
2251 vector<index_t>& triangles_;
2252 vector<double>& vertices_;
2253 double m_;
2254 index_t cur_seed_;
2255 index_t cur_vertex_;
2256 bool use_RVC_centroids_;
2257 bool select_nearest_;
2258 bool project_on_surface_;
2259 const std::vector<bool>& seed_is_locked_;
2260 bool prefer_seeds_;
2261 vector<index_t> seed_to_vertex_;
2262 bool component_on_border_;
2263 MeshFacetsAABB* AABB_;
2264 };
2265
2266 /**
2267 * \brief Implementation class for computing the restricted Delaunay
2268 * triangulation in volume mode.
2269 * \details To be used as a template argument
2270 * to RVD::for_each_primal_tetrahedron().
2271 */
2272 class GetPrimalTetrahedra {
2273 public:
2274 /**
2275 * \brief Creates a new GetPrimalTetrahedra.
2276 * \param[out] tetrahedra where to store the tetrahedra
2277 */
2278 ✗ GetPrimalTetrahedra(
2279 vector<index_t>& tetrahedra
2280 ) :
2281 ✗ tetrahedra_(tetrahedra) {
2282 }
2283
2284 /**
2285 * \brief The callback called for each primal tetrahedron.
2286 * \param[in] v1 index of the first vertex
2287 * \param[in] v2 index of the second vertex
2288 * \param[in] v3 index of the third vertex
2289 * \param[in] v4 index of the fourth vertex
2290 */
2291 ✗ void operator() (index_t v1, index_t v2, index_t v3, index_t v4) {
2292 ✗ tetrahedra_.push_back(v1);
2293 ✗ tetrahedra_.push_back(v2);
2294 ✗ tetrahedra_.push_back(v3);
2295 ✗ tetrahedra_.push_back(v4);
2296 ✗ }
2297
2298 private:
2299 vector<index_t>& tetrahedra_;
2300 };
2301
2302 12 void compute_RDT(
2303 vector<index_t>& simplices,
2304 vector<double>& embedding,
2305 RDTMode mode,
2306 const vector<bool>& seed_is_locked,
2307 MeshFacetsAABB* AABB
2308 ) override {
2309
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 12 if(volumetric_) {
2310 // For the moment, only simple mode is supported
2311 simplices.clear();
2312 ✗ RVD_.for_each_primal_tetrahedron(
2313 ✗ GetPrimalTetrahedra(simplices)
2314 );
2315 // Reorient the tetrahedra
2316 ✗ index_t nb_tetrahedra = simplices.size() / 4;
2317 ✗ for(index_t t = 0; t < nb_tetrahedra; ++t) {
2318 const double* p1 =
2319 ✗ delaunay()->vertex_ptr(simplices[4 * t]);
2320 const double* p2 =
2321 ✗ delaunay()->vertex_ptr(simplices[4 * t + 1]);
2322 const double*
2323 ✗ p3 = delaunay()->vertex_ptr(simplices[4 * t + 2]);
2324 const double*
2325 ✗ p4 = delaunay()->vertex_ptr(simplices[4 * t + 3]);
2326 ✗ if(PCK::orient_3d(p1, p2, p3, p4) < 0) {
2327 std::swap(simplices[4 * t], simplices[4 * t + 1]);
2328 }
2329 }
2330 embedding.clear();
2331 ✗ embedding.reserve(dimension_ * delaunay_->nb_vertices());
2332 ✗ for(index_t i = 0; i < delaunay_->nb_vertices(); i++) {
2333 ✗ for(coord_index_t coord = 0; coord < dimension_; coord++) {
2334 ✗ embedding.push_back(delaunay_->vertex_ptr(i)[coord]);
2335 }
2336 }
2337 } else {
2338
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 12 if((mode & RDT_MULTINERVE) != 0) {
2339 simplices.clear();
2340 embedding.clear();
2341 bool sym = RVD_.symbolic();
2342 RVD_.set_symbolic(true);
2343 RVD_.set_connected_components_priority(true);
2344 12 RVD_.for_each_polygon(
2345
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
 24 GetConnectedComponentsPrimalTriangles(
2346 RVD_, simplices, embedding, RVD_.dimension(),
2347 mode, seed_is_locked, AABB
2348 )
2349 );
2350 RVD_.set_symbolic(sym);
2351 RVD_.set_connected_components_priority(false);
2352 } else {
2353 // Simple mode: compute RDT, without any post-processing
2354 simplices.clear();
2355 ✗ RVD_.for_each_primal_triangle(
2356 ✗ GetPrimalTriangles(simplices)
2357 );
2358 embedding.clear();
2359 ✗ embedding.reserve(dimension_ * delaunay_->nb_vertices());
2360 ✗ for(index_t i = 0; i < delaunay_->nb_vertices(); i++) {
2361 ✗ for(
2362 coord_index_t coord = 0;
2363 ✗ coord < dimension_; ++coord
2364 ){
2365 embedding.push_back(
2366 ✗ delaunay_->vertex_ptr(i)[coord]
2367 );
2368 }
2369 }
2370 }
2371 }
2372 12 }
2373
2374 4020 void create_threads() override {
2375 // TODO: check if number of facets is not smaller than
2376 // number of threads
2377 // TODO: create parts even if facets range is specified
2378 // (and subdivide facets range)
2379 4020 if(
2380
2/2
✓ Branch 0 taken 410 times.
✓ Branch 1 taken 1600 times.
 4020 is_slave_ ||
2381
2/4
✓ Branch 0 taken 410 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 410 times.
✗ Branch 3 not taken.
 820 facets_begin_ != NO_INDEX || facets_end_ != NO_INDEX
2382 ) {
2383 return;
2384 }
2385 820 index_t nb_parts_in = Process::maximum_concurrent_threads();
2386 820 if(nb_parts() != nb_parts_in) {
2387
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 20 if(nb_parts_in == 1) {
2388 ✗ delete_threads();
2389 } else {
2390 vector<index_t> facet_ptr;
2391 vector<index_t> tet_ptr;
2392 20 mesh_partition(
2393
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 20 *mesh_, MESH_PARTITION_HILBERT,
2394 facet_ptr, tet_ptr, nb_parts_in
2395 );
2396 20 delete_threads();
2397
5/10
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 40 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 40 times.
✓ Branch 9 taken 10 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 100 parts_ = new thisclass[nb_parts_in];
2398 20 nb_parts_ = nb_parts_in;
2399 20 for(index_t i = 0; i < nb_parts(); ++i) {
2400 80 part(i).mesh_ = mesh_;
2401 80 part(i).set_delaunay(delaunay_);
2402 80 part(i).R3_embedding_base_ = R3_embedding_base_;
2403 80 part(i).R3_embedding_stride_ = R3_embedding_stride_;
2404 80 part(i).has_weights_ = has_weights_;
2405 80 part(i).master_ = this;
2406 80 part(i).RVD_.set_mesh(mesh_);
2407 80 part(i).set_facets_range(
2408 facet_ptr[i], facet_ptr[i + 1]
2409 );
2410 80 part(i).set_exact_predicates(RVD_.exact_predicates());
2411 part(i).set_volumetric(volumetric());
2412 part(i).set_check_SR(RVD_.check_SR());
2413 }
2414
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 6 times.
 20 if(mesh_->cells.nb() != 0) {
2415
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 4 times.
 40 for(index_t i = 0; i < nb_parts(); ++i) {
2416 32 part(i).set_tetrahedra_range(
2417 tet_ptr[i], tet_ptr[i + 1]
2418 );
2419 }
2420 }
2421
2/12
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 20 geo_assert(!Process::is_running_threads());
2422 }
2423 }
2424 }
2425
2426 80 void set_volumetric(bool x) override {
2427 120 volumetric_ = x;
2428 120 for(index_t i = 0; i < nb_parts(); ++i) {
2429 ✗ part(i).set_volumetric(x);
2430 }
2431 80 }
2432
2433 ✗ void set_facets_range(
2434 index_t facets_begin, index_t facets_end
2435 ) override {
2436 RVD_.set_facets_range(facets_begin, facets_end);
2437 40 facets_begin_ = facets_begin;
2438 40 facets_end_ = facets_end;
2439 ✗ }
2440
2441 ✗ void set_tetrahedra_range(
2442 index_t tets_begin, index_t tets_end
2443 ) override {
2444 RVD_.set_tetrahedra_range(tets_begin, tets_end);
2445 16 tets_begin_ = tets_begin;
2446 16 tets_end_ = tets_end;
2447 ✗ }
2448
2449 204 void delete_threads() override {
2450
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 92 times.
✓ Branch 2 taken 40 times.
✓ Branch 3 taken 10 times.
 284 delete[] parts_;
2451 204 parts_ = nullptr;
2452 204 nb_parts_ = 0;
2453 204 }
2454
2455 /**
2456 * \brief Gets the number of parts (or number of threads).
2457 */
2458 index_t nb_parts() const {
2459
46/308
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 480 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 160 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 800 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 160 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✓ Branch 70 taken 3 times.
✓ Branch 71 taken 120 times.
✗ Branch 72 not taken.
✓ Branch 73 taken 12 times.
✓ Branch 74 taken 12 times.
✓ Branch 75 taken 3 times.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✓ Branch 94 taken 160 times.
✓ Branch 95 taken 40 times.
✓ Branch 97 taken 40 times.
✗ Branch 98 not taken.
✓ Branch 99 taken 320 times.
✓ Branch 100 taken 80 times.
✓ Branch 102 taken 80 times.
✗ Branch 103 not taken.
✗ Branch 104 not taken.
✓ Branch 105 taken 37 times.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 109 not taken.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✗ Branch 114 not taken.
✗ Branch 115 not taken.
✓ Branch 116 taken 1 times.
✓ Branch 117 taken 40 times.
✗ Branch 118 not taken.
✓ Branch 119 taken 4 times.
✓ Branch 120 taken 4 times.
✓ Branch 121 taken 1 times.
✗ Branch 122 not taken.
✗ Branch 123 not taken.
✗ Branch 124 not taken.
✗ Branch 125 not taken.
✗ Branch 126 not taken.
✗ Branch 127 not taken.
✗ Branch 128 not taken.
✗ Branch 129 not taken.
✗ Branch 130 not taken.
✗ Branch 131 not taken.
✗ Branch 132 not taken.
✗ Branch 133 not taken.
✗ Branch 134 not taken.
✗ Branch 135 not taken.
✗ Branch 136 not taken.
✗ Branch 137 not taken.
✗ Branch 138 not taken.
✗ Branch 139 not taken.
✓ Branch 140 taken 160 times.
✓ Branch 141 taken 40 times.
✓ Branch 143 taken 40 times.
✗ Branch 144 not taken.
✗ Branch 145 not taken.
✗ Branch 146 not taken.
✗ Branch 148 not taken.
✗ Branch 149 not taken.
✗ Branch 150 not taken.
✓ Branch 151 taken 1 times.
✗ Branch 152 not taken.
✗ Branch 153 not taken.
✗ Branch 155 not taken.
✗ Branch 156 not taken.
✗ Branch 157 not taken.
✗ Branch 158 not taken.
✗ Branch 160 not taken.
✗ Branch 161 not taken.
✓ Branch 162 taken 5 times.
✓ Branch 163 taken 200 times.
✗ Branch 164 not taken.
✓ Branch 165 taken 20 times.
✓ Branch 166 taken 20 times.
✓ Branch 167 taken 5 times.
✗ Branch 168 not taken.
✗ Branch 169 not taken.
✗ Branch 170 not taken.
✗ Branch 171 not taken.
✗ Branch 172 not taken.
✗ Branch 173 not taken.
✗ Branch 174 not taken.
✗ Branch 175 not taken.
✗ Branch 176 not taken.
✗ Branch 177 not taken.
✗ Branch 178 not taken.
✗ Branch 179 not taken.
✗ Branch 180 not taken.
✗ Branch 181 not taken.
✗ Branch 182 not taken.
✗ Branch 183 not taken.
✗ Branch 184 not taken.
✗ Branch 185 not taken.
✓ Branch 186 taken 160 times.
✓ Branch 187 taken 40 times.
✓ Branch 189 taken 40 times.
✗ Branch 190 not taken.
✓ Branch 191 taken 640 times.
✓ Branch 192 taken 160 times.
✓ Branch 194 taken 160 times.
✗ Branch 195 not taken.
✗ Branch 196 not taken.
✓ Branch 197 taken 1 times.
✗ Branch 198 not taken.
✗ Branch 199 not taken.
✗ Branch 201 not taken.
✗ Branch 202 not taken.
✗ Branch 203 not taken.
✗ Branch 204 not taken.
✗ Branch 206 not taken.
✗ Branch 207 not taken.
✓ Branch 208 taken 1 times.
✓ Branch 209 taken 40 times.
✗ Branch 210 not taken.
✓ Branch 211 taken 4 times.
✓ Branch 212 taken 4 times.
✓ Branch 213 taken 1 times.
✗ Branch 214 not taken.
✗ Branch 215 not taken.
✗ Branch 216 not taken.
✗ Branch 217 not taken.
✗ Branch 218 not taken.
✗ Branch 219 not taken.
✗ Branch 220 not taken.
✗ Branch 221 not taken.
✗ Branch 222 not taken.
✗ Branch 223 not taken.
✗ Branch 224 not taken.
✗ Branch 225 not taken.
✗ Branch 226 not taken.
✗ Branch 227 not taken.
✗ Branch 228 not taken.
✗ Branch 229 not taken.
✗ Branch 230 not taken.
✗ Branch 231 not taken.
✓ Branch 232 taken 160 times.
✓ Branch 233 taken 40 times.
✓ Branch 235 taken 40 times.
✗ Branch 236 not taken.
✗ Branch 237 not taken.
✗ Branch 238 not taken.
✗ Branch 240 not taken.
✗ Branch 241 not taken.
✗ Branch 242 not taken.
✓ Branch 243 taken 1 times.
✗ Branch 244 not taken.
✗ Branch 245 not taken.
✗ Branch 247 not taken.
✗ Branch 248 not taken.
✗ Branch 249 not taken.
✗ Branch 250 not taken.
✗ Branch 252 not taken.
✗ Branch 253 not taken.
✗ Branch 254 not taken.
✗ Branch 255 not taken.
✗ Branch 256 not taken.
✗ Branch 257 not taken.
✗ Branch 258 not taken.
✗ Branch 259 not taken.
✗ Branch 260 not taken.
✗ Branch 261 not taken.
✗ Branch 262 not taken.
✗ Branch 263 not taken.
✗ Branch 264 not taken.
✗ Branch 265 not taken.
✗ Branch 266 not taken.
✗ Branch 267 not taken.
✗ Branch 268 not taken.
✗ Branch 269 not taken.
✗ Branch 270 not taken.
✗ Branch 271 not taken.
✗ Branch 272 not taken.
✗ Branch 273 not taken.
✗ Branch 274 not taken.
✗ Branch 275 not taken.
✗ Branch 276 not taken.
✗ Branch 277 not taken.
✗ Branch 278 not taken.
✗ Branch 279 not taken.
✗ Branch 281 not taken.
✗ Branch 282 not taken.
✗ Branch 283 not taken.
✗ Branch 284 not taken.
✗ Branch 286 not taken.
✗ Branch 287 not taken.
✗ Branch 288 not taken.
✗ Branch 289 not taken.
✗ Branch 290 not taken.
✗ Branch 291 not taken.
✗ Branch 293 not taken.
✗ Branch 294 not taken.
✗ Branch 295 not taken.
✗ Branch 296 not taken.
✗ Branch 298 not taken.
✗ Branch 299 not taken.
✗ Branch 300 not taken.
✗ Branch 301 not taken.
✗ Branch 302 not taken.
✗ Branch 303 not taken.
✗ Branch 304 not taken.
✗ Branch 305 not taken.
✗ Branch 306 not taken.
✗ Branch 307 not taken.
✗ Branch 308 not taken.
✗ Branch 309 not taken.
✗ Branch 310 not taken.
✗ Branch 311 not taken.
✗ Branch 312 not taken.
✗ Branch 313 not taken.
✗ Branch 314 not taken.
✗ Branch 315 not taken.
✗ Branch 316 not taken.
✗ Branch 317 not taken.
✗ Branch 318 not taken.
✗ Branch 319 not taken.
✗ Branch 320 not taken.
✗ Branch 321 not taken.
✗ Branch 322 not taken.
✗ Branch 323 not taken.
✗ Branch 324 not taken.
✗ Branch 325 not taken.
✗ Branch 327 not taken.
✗ Branch 328 not taken.
✗ Branch 329 not taken.
✗ Branch 330 not taken.
✗ Branch 332 not taken.
✗ Branch 333 not taken.
✗ Branch 334 not taken.
✗ Branch 335 not taken.
 4540 return nb_parts_;
2460 }
2461
2462 /**
2463 * \brief Gets a given part from its index.
2464 * \param[in] i index of the part
2465 * \pre \p i < nb_parts()
2466 */
2467 thisclass& part(index_t i) {
2468 geo_debug_assert(i < nb_parts());
2469
8/91
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 480 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 160 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 800 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 160 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✓ Branch 64 taken 12 times.
✗ Branch 65 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✓ Branch 74 taken 4 times.
✗ Branch 75 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✓ Branch 84 taken 20 times.
✗ Branch 85 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✓ Branch 94 taken 4 times.
✗ Branch 95 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 102 not taken.
✗ Branch 104 not taken.
✗ Branch 105 not taken.
✗ Branch 109 not taken.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✗ Branch 114 not taken.
✗ Branch 115 not taken.
 1656 return parts_[i];
2470 }
2471
2472 /**
2473 * \copydoc RestrictedVoronoiDiagram::point_allocator()
2474 */
2475 ✗ GEOGen::PointAllocator* point_allocator() override {
2476 ✗ return RVD_.point_allocator();
2477 }
2478
2479
2480 protected:
2481
2482 GenRestrictedVoronoiDiagram RVD_;
2483
2484 // For projection
2485 bool use_exact_projection_;
2486 index_t nb_triangles_;
2487 vector<index_t> triangles_;
2488 vector<vector<index_t> > stars_;
2489 Delaunay_var mesh_vertices_;
2490
2491 // One of MT_NONE, MT_LLOYD, MT_NEWTON
2492 ThreadMode thread_mode_;
2493
2494 bool is_slave_;
2495
2496 // Variables for 'master' in multithreading mode
2497 thisclass* parts_;
2498 index_t nb_parts_;
2499 Process::SpinLockArray spinlocks_;
2500
2501 // Newton mode with int. simplex
2502 IntegrationSimplex* simplex_func_;
2503
2504 // PolygonCallback mode.
2505 RVDPolygonCallback* polygon_callback_;
2506
2507 // PolyhedronCallback mode.
2508 RVDPolyhedronCallback* polyhedron_callback_;
2509
2510 // master stores argument for compute_centroids() and
2511 // compute_CVT_func_grad() to pass it to the parts.
2512 double* arg_vectors_;
2513 double* arg_scalars_;
2514
2515 // Variables for 'slaves' in multithreading mode
2516 thisclass* master_;
2517 double funcval_; // Newton mode: function value
2518
2519 protected:
2520 /**
2521 * \brief Destructor
2522 */
2523 264 ~RVD_Nd_Impl() override {
2524 172 delete_threads();
2525 436 }
2526
2527 private:
2528 /** \brief Forbids construction by copy. */
2529 RVD_Nd_Impl(const thisclass&);
2530
2531 /** \brief Forbids assignment. */
2532 thisclass& operator= (const thisclass&);
2533 };
2534 }
2535
2536 /****************************************************************************/
2537
2538 namespace GEO {
2539
2540 46 RestrictedVoronoiDiagram* RestrictedVoronoiDiagram::create(
2541 Delaunay* delaunay, Mesh* mesh,
2542 const double* R3_embedding, index_t R3_embedding_stride
2543 ) {
2544
2545 46 geo_cite("DBLP:journals/tog/EdelsbrunnerM90");
2546 46 geo_cite("DBLP:conf/compgeom/Shewchuk96");
2547 46 geo_cite("meyer:inria-00344297");
2548 46 geo_cite("DBLP:conf/gmp/YanWLL10");
2549 46 geo_cite("DBLP:journals/cad/YanWLL13");
2550 46 geo_cite("DBLP:journals/cad/Levy16");
2551
2552 delaunay->set_stores_neighbors(true);
2553 RestrictedVoronoiDiagram* result = nullptr;
2554 geo_assert(delaunay != nullptr);
2555 coord_index_t dim = delaunay->dimension();
2556
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 46 switch(dim) {
2557 ✗ case 2:
2558 result = new RVD_Nd_Impl<2>(
2559 delaunay, mesh, R3_embedding, R3_embedding_stride
2560 ✗ );
2561 break;
2562 39 case 3:
2563 result = new RVD_Nd_Impl<3>(
2564 delaunay, mesh, R3_embedding, R3_embedding_stride
2565
1/2
✓ Branch 2 taken 39 times.
✗ Branch 3 not taken.
 39 );
2566 break;
2567 1 case 4:
2568 result = new RVD_Nd_Impl<4>(
2569 delaunay, mesh, R3_embedding, R3_embedding_stride
2570
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 );
2571 break;
2572 5 case 6:
2573 result = new RVD_Nd_Impl<6>(
2574 delaunay, mesh, R3_embedding, R3_embedding_stride
2575
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 );
2576 break;
2577 1 case 8:
2578 result = new RVD_Nd_Impl<8>(
2579 delaunay, mesh, R3_embedding, R3_embedding_stride
2580
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 );
2581 break;
2582 ✗ case 20:
2583 result = new RVD_Nd_Impl<20>(
2584 delaunay, mesh, R3_embedding, R3_embedding_stride
2585 ✗ );
2586 break;
2587 ✗ case 100:
2588 result = new RVD_Nd_Impl<100>(
2589 delaunay, mesh, R3_embedding, R3_embedding_stride
2590 ✗ );
2591 break;
2592 default:
2593 ✗ geo_assert_not_reached;
2594 }
2595
4/6
✓ Branch 2 taken 46 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 46 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 35 times.
✓ Branch 7 taken 11 times.
 92 if(CmdLine::get_arg("algo:predicates") == "exact") {
2596 35 result->set_exact_predicates(true);
2597 }
2598 46 return result;
2599 }
2600
2601 126 void RestrictedVoronoiDiagram::set_delaunay(Delaunay* delaunay) {
2602 126 delaunay_ = delaunay;
2603
2/2
✓ Branch 0 taken 86 times.
✓ Branch 1 taken 40 times.
 126 if(delaunay_ != nullptr) {
2604 86 dimension_ = delaunay->dimension();
2605 } else {
2606 40 dimension_ = 0;
2607 }
2608 126 }
2609
2610
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 86 times.
 172 RestrictedVoronoiDiagram::~RestrictedVoronoiDiagram() {
2611 172 }
2612
2613 86 RestrictedVoronoiDiagram::RestrictedVoronoiDiagram(
2614 Delaunay* delaunay, Mesh* mesh,
2615 const double* R3_embedding, index_t R3_embedding_stride
2616
1/2
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
 86 ) :
2617 86 dimension_(0),
2618 86 mesh_(mesh),
2619 86 R3_embedding_base_(R3_embedding),
2620
1/2
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
 86 R3_embedding_stride_(R3_embedding_stride) {
2621
1/2
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
 86 set_delaunay(delaunay);
2622 86 has_weights_ = false;
2623 86 facets_begin_ = NO_INDEX;
2624 86 facets_end_ = NO_INDEX;
2625 86 tets_begin_ = NO_INDEX;
2626 86 tets_end_ = NO_INDEX;
2627 86 volumetric_ = false;
2628 86 }
2629
2630
2631 ✗ void RestrictedVoronoiDiagram::compute_RDT(
2632 Mesh& RDT,
2633 RDTMode mode,
2634 const vector<bool>& seed_is_locked,
2635 MeshFacetsAABB* AABB
2636 ) {
2637 vector<index_t> simplices;
2638 vector<double> embedding;
2639 ✗ compute_RDT(
2640 simplices, embedding,
2641 mode, seed_is_locked,
2642 AABB
2643 );
2644 ✗ if(volumetric()) {
2645 ✗ RDT.cells.assign_tet_mesh(dimension(),embedding,simplices,true);
2646 } else {
2647 ✗ RDT.facets.assign_triangle_mesh(
2648 dimension(),embedding,simplices,true
2649 );
2650 ✗ if((mode & RDT_DONT_REPAIR) == 0) {
2651 ✗ mesh_repair(RDT); // Needed to reorient triangles
2652 }
2653 }
2654 ✗ }
2655
2656
2657 }
2658