GCC Code Coverage Report

Directory: ./
File: lib/geogram/mesh/mesh_reorder.cpp
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 236 380 62.1%
Functions: 153 291 52.6%
Branches: 123 311 39.5%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/mesh/mesh_reorder.h>
41 #include <geogram/mesh/mesh.h>
42 #include <geogram/mesh/mesh_geometry.h>
43 #include <geogram/mesh/mesh_repair.h>
44 #include <geogram/mesh/index.h>
45 #include <geogram/delaunay/periodic.h>
46 #include <geogram/basic/permutation.h>
47 #include <geogram/basic/process.h>
48 #include <geogram/basic/logger.h>
49 #include <geogram/basic/algorithm.h>
50 #include <geogram/bibliography/bibliography.h>
51
52 namespace {
53
54 using namespace GEO;
55
56 /**
57 * \brief Splits a sequence into two ordered halves.
58 * \details The algorithm shuffles the sequence and
59 * partitions its into two halves with the same number of elements
60 * and such that the elements of the first half are smaller
61 * than the elements of the second half.
62 * \param[in] begin an iterator to the first element
63 * \param[in] end an iterator one position past the last element
64 * \param[in] cmp the comparator object
65 * \return an iterator to the middle of the sequence that separates
66 * the two halves
67 */
68 template <class IT, class CMP>
69 1710350 inline IT reorder_split(
70 IT begin, IT end, CMP cmp
71 ) {
72
2/2
✓ Branch 0 taken 77824 times.
✓ Branch 1 taken 777351 times.
 1710350 if(begin >= end) {
73 155648 return begin;
74 }
75 1554702 IT middle = begin + (end - begin) / 2;
76
1/2
✓ Branch 1 taken 777351 times.
✗ Branch 2 not taken.
 1554702 std::nth_element(begin, middle, end, cmp);
77 1554702 return middle;
78 }
79
80 /************************************************************************/
81
82 /**
83 * \brief Used by VertexMesh.
84 * \details Exposes an interface compatible with the requirement
85 * of Hilbert sort templates for a raw array of vertices.
86 */
87 class VertexArray {
88 public:
89
90 /**
91 * \brief Constructs a new VertexArray.
92 * \param[in] base address of the points
93 * \param[in] stride number of doubles between
94 * two consecutive points
95 */
96 81 VertexArray(
97 index_t nb_vertices,
98 const double* base, index_t stride
99 81 ) :
100 81 base_(base),
101 81 stride_(stride) {
102 81 nb_vertices_ = nb_vertices;
103 81 }
104
105 /**
106 * \brief Gets a vertex by its index.
107 * \param[in] i the index of the point
108 * \return a const pointer to the coordinates of the vertex
109 */
110 18656396 const double* point_ptr(index_t i) const {
111
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 18656396 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 18656396 geo_debug_assert(i < nb_vertices_);
112 18656396 return base_ + i * stride_;
113 }
114
115 private:
116 const double* base_;
117 index_t stride_;
118 index_t nb_vertices_;
119 };
120
121
122 /**
123 * \brief Exposes an interface compatible with the requirement
124 * of Hilbert sort templates for a raw array of vertices.
125 */
126 class VertexMesh {
127 public:
128 /**
129 * \brief Constructs a new VertexMesh.
130 * \param[in] base address of the points
131 * \param[in] stride number of doubles between
132 * two consecutive points
133 */
134 81 VertexMesh(
135 index_t nb_vertices,
136 const double* base, index_t stride
137 81 ) : vertices(nb_vertices, base, stride) {
138 81 }
139 VertexArray vertices;
140 };
141
142 /************************************************************************/
143
144 /**
145 * \brief The generic comparator class for Hilbert vertex
146 * ordering.
147 * \tparam COORD the coordinate to compare
148 * \tparam UP if true, use direct order, else use reverse order
149 * \tparam MESH the class that represents meshes
150 */
151 template <int COORD, bool UP, class MESH>
152 struct Hilbert_vcmp {
153 };
154
155 /**
156 * \brief Specialization (UP=true) of the generic comparator class
157 * for Hilbert vertex ordering.
158 * \see Hilbert_vcmp
159 * \tparam COORD the coordinate to compare
160 * \tparam MESH the class that represents meshes
161 */
162 template <int COORD, class MESH>
163 struct Hilbert_vcmp<COORD, true, MESH> {
164
165 /**
166 * \brief Constructs a new Hilbert_vcmp.
167 * \param[in] mesh the mesh in which the compared
168 * points reside.
169 */
170 357218 Hilbert_vcmp(const MESH& mesh) :
171 357218 mesh_(mesh) {
172 357218 }
173
174 /**
175 * \brief Compares two points.
176 * \param[in] i1 index of the first point to compare
177 * \param[in] i2 index of the second point to compare
178 * \return true if point \p i1 is before point \p i2,
179 * false otherwise.
180 */
181 8322936 bool operator() (index_t i1, index_t i2) {
182 return
183 8322936 mesh_.vertices.point_ptr(i1)[COORD] <
184 8322936 mesh_.vertices.point_ptr(i2)[COORD];
185 }
186
187 const MESH& mesh_;
188 };
189
190 /**
191 * \brief Specialization (UP=false) of the generic comparator class
192 * for Hilbert vertex ordering.
193 * \see Hilbert_vcmp
194 * \tparam COORD the coordinate to compare
195 * \tparam MESH the class that represents meshes
196 */
197 template <int COORD, class MESH>
198 struct Hilbert_vcmp<COORD, false, MESH> {
199
200 /**
201 * \brief Constructs a new Hilbert_vcmp.
202 * \param[in] mesh the mesh in which the compared
203 * points reside.
204 */
205 359370 Hilbert_vcmp(const MESH& mesh) :
206 359370 mesh_(mesh) {
207 359370 }
208
209 /**
210 * \brief Compares two points.
211 * \param[in] i1 index of the first point to compare
212 * \param[in] i2 index of the second point to compare
213 * \return true if point \p i1 is before point \p i2,
214 * false otherwise.
215 */
216 10577524 bool operator() (index_t i1, index_t i2) {
217 return
218 10577524 mesh_.vertices.point_ptr(i1)[COORD] >
219 10577524 mesh_.vertices.point_ptr(i2)[COORD];
220 }
221
222 const MESH& mesh_;
223 };
224
225 /************************************************************************/
226
227 /**
228 * \brief Comparator class for Morton vertex
229 * ordering.
230 * \tparam COORD the coordinate to compare
231 * \tparam UP ignored in Morton order
232 * \tparam MESH the class that represents meshes
233 */
234 template <int COORD, bool UP, class MESH>
235 struct Morton_vcmp {
236
237 /**
238 * \brief Constructs a new Morton_vcmp.
239 * \param[in] mesh the mesh in which the compared
240 * points reside.
241 */
242 ✗ Morton_vcmp(const MESH& mesh) :
243 ✗ mesh_(mesh) {
244 ✗ }
245
246 /**
247 * \brief Compares two points.
248 * \param[in] i1 index of the first point to compare
249 * \param[in] i2 index of the second point to compare
250 * \return true if point \p i1 is before point \p i2,
251 * false otherwise.
252 */
253 ✗ bool operator() (index_t i1, index_t i2) {
254 return
255 ✗ mesh_.vertices.point_ptr(i1)[COORD] <
256 ✗ mesh_.vertices.point_ptr(i2)[COORD];
257 }
258
259 const MESH& mesh_;
260 };
261
262 /************************************************************************/
263
264 #ifndef GEOGRAM_PSM
265
266 /**
267 * \brief Base class for facets ordering.
268 * \tparam COORD the coordinate to compare
269 * \tparam MESH the class that represents meshes
270 */
271 template <int COORD, class MESH>
272 class Base_fcmp {
273 public:
274 /**
275 * \brief Constructs a new Base_vcmp.
276 * \param[in] mesh the mesh in which the compared
277 * facets reside.
278 */
279 989674 Base_fcmp(const MESH& mesh) :
280 989674 mesh_(mesh) {
281 989674 }
282
283 /**
284 * \brief Computes the compared coordinate from a facet index.
285 * \param[in] f the index of the facet
286 * \return the coordinate at the center of facet \p f
287 */
288 28450724 double center(index_t f) const {
289 28450724 double result = 0.0;
290 28450724 double s = 1.0 / double(mesh_.facets.nb_vertices(f));
291
3/4
✓ Branch 1 taken 14225362 times.
✗ Branch 2 not taken.
✓ Branch 7 taken 42683490 times.
✓ Branch 8 taken 14225362 times.
 113817704 for(index_t c: mesh_.facets.corners(f)) {
292
1/2
✓ Branch 1 taken 42683490 times.
✗ Branch 2 not taken.
 85366980 result += s*mesh_.vertices.point_ptr(
293
1/2
✓ Branch 1 taken 42683490 times.
✗ Branch 2 not taken.
 85366980 mesh_.facet_corners.vertex(c)
294 85366980 )[COORD];
295 }
296 28450724 return result;
297 }
298
299 private:
300 const MESH& mesh_;
301 };
302
303 /**
304 * \brief The generic comparator class for Hilbert facet
305 * ordering.
306 * \tparam COORD the coordinate to compare
307 * \tparam UP if true, use direct order, else use reverse order
308 * \tparam MESH the class that represents meshes
309 */
310 template <int COORD, bool UP, class MESH>
311 struct Hilbert_fcmp {
312 };
313
314 /**
315 * \brief Specialization (UP=true) of the generic comparator class
316 * for Hilbert vertex ordering.
317 * \see Hilbert_vcmp
318 * \tparam COORD the coordinate to compare
319 * \tparam MESH the class that represents meshes
320 */
321 template <int COORD, class MESH>
322 class Hilbert_fcmp<COORD, true, MESH> : public Base_fcmp<COORD, MESH> {
323 public:
324 /**
325 * \brief Constructs a new Hilbert_fcmp.
326 * \param[in] mesh the mesh in which the compared
327 * facets reside.
328 */
329 11278 Hilbert_fcmp(const MESH& mesh) :
330 11278 Base_fcmp<COORD, MESH>(mesh) {
331 11278 }
332
333 /**
334 * \brief Compares two facets.
335 * \param[in] f1 index of the first facet to compare
336 * \param[in] f2 index of the second facet to compare
337 * \return true if facet \p f1 is before facet \p f2,
338 * false otherwise.
339 */
340 224040 bool operator() (index_t f1, index_t f2) {
341 224040 return this->center(f1) < this->center(f2);
342 }
343 };
344
345 /**
346 * \brief Specialization (UP=false) of the generic comparator class
347 * for Hilbert vertex ordering.
348 * \see Hilbert_vcmp
349 * \tparam COORD the coordinate to compare
350 * \tparam MESH the class that represents meshes
351 */
352 template <int COORD, class MESH>
353 class Hilbert_fcmp<COORD, false, MESH> : public Base_fcmp<COORD, MESH> {
354 public:
355 /**
356 * \brief Constructs a new Hilbert_fcmp.
357 * \param[in] mesh the mesh in which the compared
358 * facets reside.
359 */
360 11430 Hilbert_fcmp(const MESH& mesh) :
361 11430 Base_fcmp<COORD, MESH>(mesh) {
362 11430 }
363
364 /**
365 * \brief Compares two facets.
366 * \param[in] f1 index of the first facet to compare
367 * \param[in] f2 index of the second facet to compare
368 * \return true if facet \p f1 is before facet \p f2,
369 * false otherwise.
370 */
371 285892 bool operator() (index_t f1, index_t f2) {
372 285892 return this->center(f1) > this->center(f2);
373 }
374 };
375
376 /**
377 * \brief Comparator class for Morton facet
378 * ordering.
379 * \tparam COORD the coordinate to compare
380 * \tparam UP ignored in Morton order
381 * \tparam MESH the class that represents meshes
382 */
383 template <int COORD, bool UP, class MESH>
384 class Morton_fcmp : public Base_fcmp<COORD, MESH> {
385 public:
386 /**
387 * \brief Constructs a new Morton_fcmp.
388 * \param[in] mesh the mesh in which the compared
389 * facets reside.
390 */
391 966966 Morton_fcmp(const MESH& mesh) :
392 966966 Base_fcmp<COORD, MESH>(mesh) {
393 966966 }
394
395 /**
396 * \brief Compares two facets.
397 * \param[in] f1 index of the first facet to compare
398 * \param[in] f2 index of the second facet to compare
399 * \return true if facet \p f1 is before facet \p f2,
400 * false otherwise.
401 */
402 13715430 bool operator() (index_t f1, index_t f2) {
403 13715430 return this->center(f1) < this->center(f2);
404 }
405 };
406
407 /************************************************************************/
408
409 /**
410 * \brief Base class for tetrahedra ordering.
411 * \tparam COORD the coordinate to compare
412 * \tparam MESH the class that represents meshes
413 */
414 template <int COORD, class MESH>
415 class Base_tcmp {
416 public:
417 /**
418 * \brief Constructs a new Base_tcmp.
419 * \param[in] mesh the mesh in which the compared
420 * tetrahedra reside.
421 */
422 4088 Base_tcmp(const MESH& mesh) :
423 4088 mesh_(mesh) {
424 4088 }
425
426 /**
427 * \brief Computes the compared coordinate from a tetra index.
428 * \param[in] t the index of the tetra
429 * \return the coordinate at the center of tetra \p f
430 */
431 135904 double center(index_t t) const {
432 135904 double result = 0.0;
433
2/2
✓ Branch 0 taken 271808 times.
✓ Branch 1 taken 67952 times.
 679520 for(index_t lv = 0; lv < 4; ++lv) {
434 543616 result += mesh_.vertices.point_ptr(
435 543616 mesh_.cells.vertex(t, lv)
436 543616 )[COORD];
437 }
438 135904 return result;
439 }
440
441 private:
442 const MESH& mesh_;
443 };
444
445 /**
446 * \brief The generic comparator class for Hilbert tetra
447 * ordering.
448 * \tparam COORD the coordinate to compare
449 * \tparam UP if true, use direct order, else use reverse order
450 * \tparam MESH the class that represents meshes
451 */
452 template <int COORD, bool UP, class MESH>
453 struct Hilbert_tcmp {
454 };
455
456 /**
457 * \brief Specialization (UP=true) of the generic comparator class
458 * for Hilbert tetra ordering.
459 * \see Hilbert_tcmp
460 * \tparam COORD the coordinate to compare
461 * \tparam MESH the class that represents meshes
462 */
463 template <int COORD, class MESH>
464 class Hilbert_tcmp<COORD, true, MESH> : public Base_tcmp<COORD, MESH> {
465 public:
466 /**
467 * \brief Constructs a new Hilbert_tcmp.
468 * \param[in] mesh the mesh in which the compared
469 * tetrahedra reside.
470 */
471 2016 Hilbert_tcmp(const MESH& mesh) :
472 2016 Base_tcmp<COORD, MESH>(mesh) {
473 2016 }
474
475 /**
476 * \brief Compares two tetrahedra.
477 * \param[in] t1 index of the first tetra to compare
478 * \param[in] t2 index of the second tetra to compare
479 * \return true if tetra \p t1 is before tetra \p t2,
480 * false otherwise.
481 */
482 28424 bool operator() (index_t t1, index_t t2) {
483 28424 return this->center(t1) < this->center(t2);
484 }
485 };
486
487 /**
488 * \brief Specialization (UP=false) of the generic comparator class
489 * for Hilbert tetra ordering.
490 * \see Hilbert_tcmp
491 * \tparam COORD the coordinate to compare
492 * \tparam MESH the class that represents meshes
493 */
494 template <int COORD, class MESH>
495 class Hilbert_tcmp<COORD, false, MESH> : public Base_tcmp<COORD, MESH> {
496 public:
497 /**
498 * \brief Constructs a new Hilbert_tcmp.
499 * \param[in] mesh the mesh in which the compared
500 * tetrahedra reside.
501 */
502 2072 Hilbert_tcmp(const MESH& mesh) :
503 2072 Base_tcmp<COORD, MESH>(mesh) {
504 2072 }
505
506 /**
507 * \brief Compares two tetrahedra.
508 * \param[in] t1 index of the first tetra to compare
509 * \param[in] t2 index of the second tetra to compare
510 * \return true if tetra \p t1 is before tetra \p t2,
511 * false otherwise.
512 */
513 39528 bool operator() (index_t t1, index_t t2) {
514 39528 return this->center(t1) > this->center(t2);
515 }
516 };
517
518 /**
519 * \brief Comparator class for Morton tetra
520 * ordering.
521 * \tparam COORD the coordinate to compare
522 * \tparam UP ignored in Morton order
523 * \tparam MESH the class that represents meshes
524 */
525 template <int COORD, bool UP, class MESH>
526 class Morton_tcmp : public Base_tcmp<COORD, MESH> {
527 public:
528 /**
529 * \brief Constructs a new Morton_tcmp.
530 * \param[in] mesh the mesh in which the compared
531 * tetrahedra reside.
532 */
533 ✗ Morton_tcmp(const MESH& mesh) :
534 ✗ Base_tcmp<COORD, MESH>(mesh) {
535 ✗ }
536
537 /**
538 * \brief Compares two tetrahedra.
539 * \param[in] t1 index of the first tetra to compare
540 * \param[in] t2 index of the second tetra to compare
541 * \return true if tetra \p t1 is before tetra \p t2,
542 * false otherwise.
543 */
544 ✗ bool operator() (index_t t1, index_t t2) {
545 ✗ return this->center(t1) < this->center(t2);
546 }
547 };
548
549 /************************************************************************/
550
551 /**
552 * \brief Base class for cells ordering.
553 * \tparam COORD the coordinate to compare
554 * \tparam MESH the class that represents meshes
555 */
556 template <int COORD, class MESH>
557 class Base_ccmp {
558 public:
559 /**
560 * \brief Constructs a new Base_ccmp.
561 * \param[in] mesh the mesh in which the compared
562 * cells reside.
563 */
564 ✗ Base_ccmp(const MESH& mesh) :
565 ✗ mesh_(mesh) {
566 ✗ }
567
568 /**
569 * \brief Computes the compared coordinate from a cell index.
570 * \param[in] c the index of the cell
571 * \return the coordinate at the center of cell \p c
572 */
573 ✗ double center(index_t c) const {
574 ✗ double result = 0.0;
575 ✗ for(index_t lv = 0; lv < mesh_.cells.nb_vertices(c); ++lv) {
576 ✗ result += mesh_.vertices.point_ptr(
577 ✗ mesh_.cells.vertex(c, lv)
578 ✗ )[COORD];
579 }
580 ✗ return result / double(mesh_.cells.nb_vertices(c));
581 }
582
583 private:
584 const MESH& mesh_;
585 };
586
587 /**
588 * \brief The generic comparator class for Hilbert cell
589 * ordering.
590 * \tparam COORD the coordinate to compare
591 * \tparam UP if true, use direct order, else use reverse order
592 * \tparam MESH the class that represents meshes
593 */
594 template <int COORD, bool UP, class MESH>
595 struct Hilbert_ccmp {
596 };
597
598 /**
599 * \brief Specialization (UP=true) of the generic comparator class
600 * for Hilbert cell ordering.
601 * \see Hilbert_ccmp
602 * \tparam COORD the coordinate to compare
603 * \tparam MESH the class that represents meshes
604 */
605 template <int COORD, class MESH>
606 class Hilbert_ccmp<COORD, true, MESH> : public Base_ccmp<COORD, MESH> {
607 public:
608 /**
609 * \brief Constructs a new Hilbert_ccmp.
610 * \param[in] mesh the mesh in which the compared
611 * cells reside.
612 */
613 ✗ Hilbert_ccmp(const MESH& mesh) :
614 ✗ Base_ccmp<COORD, MESH>(mesh) {
615 ✗ }
616
617 /**
618 * \brief Compares two cells
619 * \param[in] c1 index of the first cell to compare
620 * \param[in] c2 index of the second cell to compare
621 * \return true if cell \p c1 is before cell \p c2,
622 * false otherwise.
623 */
624 ✗ bool operator() (index_t c1, index_t c2) {
625 ✗ return this->center(c1) < this->center(c2);
626 }
627 };
628
629 /**
630 * \brief Specialization (UP=false) of the generic comparator class
631 * for Hilbert cell ordering.
632 * \see Hilbert_ccmp
633 * \tparam COORD the coordinate to compare
634 * \tparam MESH the class that represents meshes
635 */
636 template <int COORD, class MESH>
637 class Hilbert_ccmp<COORD, false, MESH> : public Base_ccmp<COORD, MESH> {
638 public:
639 /**
640 * \brief Constructs a new Hilbert_ccmp.
641 * \param[in] mesh the mesh in which the compared
642 * tetrahedra reside.
643 */
644 ✗ Hilbert_ccmp(const MESH& mesh) :
645 ✗ Base_ccmp<COORD, MESH>(mesh) {
646 ✗ }
647
648 /**
649 * \brief Compares two cells.
650 * \param[in] c1 index of the first cell to compare
651 * \param[in] c2 index of the second cell to compare
652 * \return true if cell \p c1 is before cell \p c2,
653 * false otherwise.
654 */
655 ✗ bool operator() (index_t c1, index_t c2) {
656 ✗ return this->center(c1) > this->center(c2);
657 }
658 };
659
660 /**
661 * \brief Comparator class for Morton cell
662 * ordering.
663 * \tparam COORD the coordinate to compare
664 * \tparam UP ignored in Morton order
665 * \tparam MESH the class that represents meshes
666 */
667 template <int COORD, bool UP, class MESH>
668 class Morton_ccmp : public Base_ccmp<COORD, MESH> {
669 public:
670 /**
671 * \brief Constructs a new Morton_ccmp.
672 * \param[in] mesh the mesh in which the compared
673 * cells reside.
674 */
675 ✗ Morton_ccmp(const MESH& mesh) :
676 ✗ Base_ccmp<COORD, MESH>(mesh) {
677 ✗ }
678
679 /**
680 * \brief Compares two tetrahedra.
681 * \param[in] c1 index of the first cell to compare
682 * \param[in] c2 index of the second cell to compare
683 * \return true if cell \p c1 is before cell \p c2,
684 * false otherwise.
685 */
686 ✗ bool operator() (index_t c1, index_t c2) {
687 ✗ return this->center(c1) < this->center(c2);
688 }
689 };
690
691 #endif
692
693 /************************************************************************/
694
695 /**
696 * \brief Generic class for sorting arbitrary elements in
697 * Hilbert and Morton orders in 3d.
698 * \details The implementation is inspired by:
699 * - Christophe Delage and Olivier Devillers. Spatial Sorting.
700 * In CGAL User and Reference Manual. CGAL Editorial Board,
701 * 3.9 edition, 2011
702 * \tparam CMP the comparator class for ordering the elements. CMP
703 * is itself a template parameterized by~:
704 * - COORD the coordinate along which elements should be
705 * sorted
706 * - UP a boolean that indicates whether direct or reverse
707 * order should be used
708 * - MESH the class that represents meshes
709 * \tparam MESH the class that represents meshes
710 */
711 template <template <int COORD, bool UP, class MESH> class CMP, class MESH>
712 struct HilbertSort3d {
713
714 /**
715 * \brief Low-level recursive spatial sorting function
716 * \details This function is recursive
717 * \param[in] M the mesh in which the elements reside
718 * \param[in] begin an iterator that points to the
719 * first element of the sequence
720 * \param[in] end an iterator that points one position past the
721 * last element of the sequence
722 * \param[in] limit subsequences smaller than limit are left unsorted
723 * \tparam COORDX the first coordinate, can be 0,1 or 2. The second
724 * and third coordinates are COORDX+1 modulo 3 and COORDX+2 modulo 3
725 * respectively
726 * \tparam UPX whether ordering along the first coordinate
727 * is direct or inverse
728 * \tparam UPY whether ordering along the second coordinate
729 * is direct or inverse
730 * \tparam UPZ whether ordering along the third coordinate
731 * is direct or inverse
732 */
733 template <int COORDX, bool UPX, bool UPY, bool UPZ, class IT>
734 1176192 static void sort(
735 const MESH& M, IT begin, IT end, index_t limit = 1
736 ) {
737 1176192 const int COORDY = (COORDX + 1) % 3, COORDZ = (COORDY + 1) % 3;
738
2/2
✓ Branch 0 taken 514626 times.
✓ Branch 1 taken 73470 times.
 1176192 if(end - begin <= signed_index_t(limit)) {
739 1029252 return;
740 }
741 146940 IT m0 = begin, m8 = end;
742
3/5
✓ Branch 1 taken 59587 times.
✓ Branch 2 taken 13883 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 59587 times.
✗ Branch 5 not taken.
 146940 IT m4 = reorder_split(m0, m8, CMP<COORDX, UPX, MESH>(M));
743
3/5
✓ Branch 1 taken 59138 times.
✓ Branch 2 taken 14332 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 59138 times.
✗ Branch 5 not taken.
 146940 IT m2 = reorder_split(m0, m4, CMP<COORDY, UPY, MESH>(M));
744
3/5
✓ Branch 1 taken 59106 times.
✓ Branch 2 taken 14364 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 59106 times.
✗ Branch 5 not taken.
 146940 IT m1 = reorder_split(m0, m2, CMP<COORDZ, UPZ, MESH>(M));
745
3/5
✓ Branch 1 taken 59115 times.
✓ Branch 2 taken 14355 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 59115 times.
✗ Branch 5 not taken.
 146940 IT m3 = reorder_split(m2, m4, CMP<COORDZ, !UPZ, MESH>(M));
746
3/5
✓ Branch 1 taken 59187 times.
✓ Branch 2 taken 14283 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 59187 times.
✗ Branch 5 not taken.
 146940 IT m6 = reorder_split(m4, m8, CMP<COORDY, !UPY, MESH>(M));
747
3/5
✓ Branch 1 taken 59106 times.
✓ Branch 2 taken 14364 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 59106 times.
✗ Branch 5 not taken.
 146940 IT m5 = reorder_split(m4, m6, CMP<COORDZ, UPZ, MESH>(M));
748
3/5
✓ Branch 1 taken 59115 times.
✓ Branch 2 taken 14355 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 59115 times.
✗ Branch 5 not taken.
 146940 IT m7 = reorder_split(m6, m8, CMP<COORDZ, !UPZ, MESH>(M));
749
1/2
✓ Branch 1 taken 73470 times.
✗ Branch 2 not taken.
 146940 sort<COORDZ, UPZ, UPX, UPY>(M, m0, m1);
750
1/2
✓ Branch 1 taken 73470 times.
✗ Branch 2 not taken.
 146940 sort<COORDY, UPY, UPZ, UPX>(M, m1, m2);
751
1/2
✓ Branch 1 taken 73470 times.
✗ Branch 2 not taken.
 146940 sort<COORDY, UPY, UPZ, UPX>(M, m2, m3);
752
1/2
✓ Branch 1 taken 73470 times.
✗ Branch 2 not taken.
 146940 sort<COORDX, UPX, !UPY, !UPZ>(M, m3, m4);
753
1/2
✓ Branch 1 taken 73470 times.
✗ Branch 2 not taken.
 146940 sort<COORDX, UPX, !UPY, !UPZ>(M, m4, m5);
754
1/2
✓ Branch 1 taken 73470 times.
✗ Branch 2 not taken.
 146940 sort<COORDY, !UPY, UPZ, !UPX>(M, m5, m6);
755
1/2
✓ Branch 1 taken 73470 times.
✗ Branch 2 not taken.
 146940 sort<COORDY, !UPY, UPZ, !UPX>(M, m6, m7);
756
1/2
✓ Branch 1 taken 73470 times.
✗ Branch 2 not taken.
 146940 sort<COORDZ, !UPZ, !UPX, UPY>(M, m7, m8);
757 }
758
759 /**
760 * \brief Sorts a sequence of elements spatially.
761 * \details This function does an indirect sort,
762 * in the sense that a sequence
763 * of indices that refer to the elements is sorted.
764 * This function uses a multithreaded implementation.
765 * \param[in] M the mesh in which the elements to sort reside
766 * \param[in] b an iterator to the first index to be sorted
767 * \param[in] e an iterator one position past the last index
768 * to be sorted
769 * \param[in] limit subsequences smaller than limit are left unsorted
770 */
771 238 HilbertSort3d(
772 const MESH& M,
773 vector<index_t>::iterator b,
774 vector<index_t>::iterator e,
775 index_t limit = 1
776 ) :
777 238 M_(M)
778 {
779
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 119 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 238 geo_debug_assert(e >= b);
780 238 geo_cite_with_info(
781 "WEB:SpatialSorting",
782 "The implementation of spatial sort is inspired by "
783 "the use of \\verb|std::nth_element()| and the recursive"
784 " template in the spatial sort package of CGAL"
785 );
786
787 // If the sequence is smaller than the limit, skip it
788
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 119 times.
 238 if(index_t(e - b) <= limit) {
789 ✗ return;
790 }
791
792 // If the sequence is smaller than 1024, use sequential sorting
793
2/2
✓ Branch 0 taken 88 times.
✓ Branch 1 taken 31 times.
 238 if(index_t(e - b) < 1024) {
794 176 sort<0, false, false, false>(M_, b, e);
795 176 return;
796 }
797
798 // Parallel sorting (2 then 4 then 8 sorts in parallel)
799
800 // Unfortunately we cannot access consts/constexprs for template
801 // arguments in lambdas in all compilers (gcc/clang OK but not
802 // MSVC) so I'm using macros here (it is ugly, but it is not a
803 // big drama), and I prefer that instead of hardwired constants
804 // that would make the code more difficult to read.
805
806 # define COORDX 0
807 # define COORDY 1
808 # define COORDZ 2
809 # define UPX false
810 # define UPY false
811 # define UPZ false
812
813 62 m0_ = b;
814 62 m8_ = e;
815
3/5
✓ Branch 1 taken 28 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
 62 m4_ = reorder_split(m0_, m8_, CMP<COORDX, UPX, MESH>(M));
816
817
818
1/2
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
 62 parallel(
819 155 [this]() { m2_ = reorder_split(m0_, m4_, CMP<COORDY, UPY, MESH>(M_)); },
820 155 [this]() { m6_ = reorder_split(m4_, m8_, CMP<COORDY, !UPY, MESH>(M_)); }
821 );
822
823
1/2
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
 62 parallel(
824 155 [this]() { m1_ = reorder_split(m0_, m2_, CMP<COORDZ, UPZ, MESH>(M_)); },
825 155 [this]() { m3_ = reorder_split(m2_, m4_, CMP<COORDZ, !UPZ, MESH>(M_)); },
826 155 [this]() { m5_ = reorder_split(m4_, m6_, CMP<COORDZ, UPZ, MESH>(M_)); },
827 155 [this]() { m7_ = reorder_split(m6_, m8_, CMP<COORDZ, !UPZ, MESH>(M_)); }
828 );
829
830
1/2
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
 62 parallel(
831 155 [this]() { sort<COORDZ, UPZ, UPX, UPY>(M_, m0_, m1_); },
832 155 [this]() { sort<COORDY, UPY, UPZ, UPX>(M_, m1_, m2_); },
833 155 [this]() { sort<COORDY, UPY, UPZ, UPX>(M_, m2_, m3_); },
834 155 [this]() { sort<COORDX, UPX, !UPY, !UPZ>(M_, m3_, m4_); },
835 155 [this]() { sort<COORDX, UPX, !UPY, !UPZ>(M_, m4_, m5_); },
836 155 [this]() { sort<COORDY, !UPY, UPZ, !UPX>(M_, m5_, m6_); },
837 155 [this]() { sort<COORDY, !UPY, UPZ, !UPX>(M_, m6_, m7_); },
838 155 [this]() { sort<COORDZ, !UPZ, !UPX, UPY>(M_, m7_, m8_); }
839 );
840
841 # undef COORDX
842 # undef COORDY
843 # undef COORDZ
844 # undef UPX
845 # undef UPY
846 # undef UPZ
847
848 }
849
850 private:
851 const MESH& M_;
852 vector<index_t>::iterator
853 m0_, m1_, m2_, m3_, m4_, m5_, m6_, m7_, m8_;
854 };
855
856 /************************************************************************/
857
858 /**
859 * \brief Generic class for sorting arbitrary elements in
860 * Hilbert and Morton orders in 3d.
861 * \details The implementation is inspired by:
862 * - Christophe Delage and Olivier Devillers. Spatial Sorting.
863 * In CGAL User and Reference Manual. CGAL Editorial Board,
864 * 3.9 edition, 2011
865 * \tparam CMP the comparator class for ordering the elements. CMP
866 * is itself a template parameterized by~:
867 * - COORD the coordinate along which elements should be
868 * sorted
869 * - UP a boolean that indicates whether direct or reverse
870 * order should be used
871 * - MESH the class that represents meshes
872 * \tparam MESH the class that represents meshes
873 */
874 template <template <int COORD, bool UP, class MESH> class CMP, class MESH>
875 struct HilbertSort2d {
876
877 /**
878 * \brief Low-level recursive spatial sorting function
879 * \details This function is recursive
880 * \param[in] M the mesh in which the elements reside
881 * \param[in] begin an iterator that points to the
882 * first element of the sequence
883 * \param[in] end an iterator that points one position past the
884 * last element of the sequence
885 * \param[in] limit subsequences smaller than limit are left unsorted
886 * \tparam COORDX the first coordinate, can be 0,1 or 2. The second
887 * coordinate is COORDX+1 modulo 2.
888 * \tparam UPX whether ordering along the first coordinate
889 * is direct or inverse
890 * \tparam UPY whether ordering along the second coordinate
891 * is direct or inverse
892 */
893 template <int COORDX, bool UPX, bool UPY, class IT>
894 908586 static void sort(
895 const MESH& M, IT begin, IT end, index_t limit = 1
896 ) {
897 908586 const int COORDY = (COORDX + 1) % 2;
898
2/2
✓ Branch 0 taken 340737 times.
✓ Branch 1 taken 113556 times.
 908586 if(end - begin <= signed_index_t(limit)) {
899 681474 return;
900 }
901 227112 IT m0 = begin, m4 = end;
902
903
1/2
✓ Branch 2 taken 113556 times.
✗ Branch 3 not taken.
 227112 IT m2 = reorder_split (m0, m4, CMP<COORDX, UPX, MESH>(M));
904
1/2
✓ Branch 2 taken 113556 times.
✗ Branch 3 not taken.
 227112 IT m1 = reorder_split (m0, m2, CMP<COORDY, UPY, MESH>(M));
905
1/2
✓ Branch 2 taken 113556 times.
✗ Branch 3 not taken.
 227112 IT m3 = reorder_split (m2, m4, CMP<COORDY, !UPY, MESH>(M));
906
907
1/2
✓ Branch 1 taken 113556 times.
✗ Branch 2 not taken.
 227112 sort<COORDY, UPY, UPX> (M, m0, m1);
908
1/2
✓ Branch 1 taken 113556 times.
✗ Branch 2 not taken.
 227112 sort<COORDX, UPX, UPY> (M, m1, m2);
909
1/2
✓ Branch 1 taken 113556 times.
✗ Branch 2 not taken.
 227112 sort<COORDX, UPX, UPY> (M, m2, m3);
910
1/2
✓ Branch 1 taken 113556 times.
✗ Branch 2 not taken.
 227112 sort<COORDY,!UPY,!UPX> (M, m3, m4);
911 }
912
913 /**
914 * \brief Sorts a sequence of elements spatially.
915 * \details This function does an indirect sort,
916 * in the sense that a sequence
917 * of indices that refer to the elements is sorted.
918 * This function uses a multithreaded implementation.
919 * \param[in] M the mesh in which the elements to sort reside
920 * \param[in] b an iterator to the first index to be sorted
921 * \param[in] e an iterator one position past the last index
922 * to be sorted
923 * \param[in] limit subsequences smaller than limit are left unsorted
924 */
925 69 HilbertSort2d(
926 const MESH& M,
927 vector<index_t>::iterator b,
928 vector<index_t>::iterator e,
929 index_t limit = 1
930 ) :
931 69 M_(M)
932 {
933
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 69 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 69 geo_debug_assert(e > b);
934 69 geo_cite_with_info(
935 "WEB:SpatialSorting",
936 "The implementation of spatial sort is inspired by "
937 "the use of \\verb|std::nth_element()| and the recursive"
938 " template in the spatial sort package of CGAL"
939 );
940
941 // If the sequence is smaller than the limit, skip it
942
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 69 times.
 69 if(index_t(e - b) <= limit) {
943 ✗ return;
944 }
945 69 sort<0, false, false>(M_, b, e);
946 }
947 private:
948 const MESH& M_;
949 };
950
951 /************************************************************************/
952
953 #ifndef GEOGRAM_PSM
954
955 /**
956 * \brief Sorts the vertices of a mesh according to the Hilbert ordering.
957 * \details The function does not change the mesh, it computes instead
958 * the permutation. The permutation can then be reused to order other
959 * arrays that may depend on the order of the vertices in the mesh (i.e.
960 * attributes).
961 * \param[in] M the mesh where the vertices to be sorted reside
962 * \param[out] sorted_indices the permutation to be applied
963 to the vertices
964 */
965 10 void hilbert_vsort_3d(
966 const Mesh& M, vector<index_t>& sorted_indices
967 ) {
968 10 sorted_indices.resize(M.vertices.nb());
969
2/2
✓ Branch 4 taken 4489 times.
✓ Branch 5 taken 10 times.
 4499 for(index_t i: M.vertices) {
970
1/2
✓ Branch 1 taken 4489 times.
✗ Branch 2 not taken.
 4489 sorted_indices[i] = i;
971 }
972 10 HilbertSort3d<Hilbert_vcmp, Mesh>(
973 M, sorted_indices.begin(), sorted_indices.end()
974
1/2
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
 10 );
975 10 }
976
977 /**
978 * \brief Sorts the facets of a mesh according to the Hilbert ordering.
979 * \details The function does not change the mesh, it computes instead
980 * the permutation. The permutation can then be reused to order other
981 * arrays that may depend on the order of the facets in the mesh (i.e.
982 * attributes).
983 * \param[in] M the mesh where the facets to be sorted reside
984 * \param[out] sorted_indices the permutation to be
985 * applied to the facets
986 */
987 10 void hilbert_fsort_3d(
988 const Mesh& M, vector<index_t>& sorted_indices
989 ) {
990 10 sorted_indices.resize(M.facets.nb());
991
2/2
✓ Branch 4 taken 8566 times.
✓ Branch 5 taken 10 times.
 8576 for(index_t i: M.facets) {
992
1/2
✓ Branch 1 taken 8566 times.
✗ Branch 2 not taken.
 8566 sorted_indices[i] = i;
993 }
994 10 HilbertSort3d<Hilbert_fcmp, Mesh>(
995 M, sorted_indices.begin(), sorted_indices.end()
996
1/2
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
 10 );
997 10 }
998
999 /**
1000 * \brief Sorts the cells of a mesh according to the Hilbert ordering.
1001 * \details The function does not change the mesh, it computes instead
1002 * the permutation. The permutation can then be reused to order other
1003 * arrays that may depend on the order of the tets in the mesh (i.e.
1004 * attributes).
1005 * \param[in] M the mesh where the cells to be sorted reside
1006 * \param[out] sorted_indices the permutation to be applied to the tets
1007 */
1008 4 void hilbert_csort_3d(
1009 const Mesh& M, vector<index_t>& sorted_indices
1010 ) {
1011 4 sorted_indices.resize(M.cells.nb());
1012
2/2
✓ Branch 4 taken 1536 times.
✓ Branch 5 taken 4 times.
 1540 for(index_t i: M.cells) {
1013
1/2
✓ Branch 1 taken 1536 times.
✗ Branch 2 not taken.
 1536 sorted_indices[i] = i;
1014 }
1015
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 if(M.cells.are_simplices()) {
1016 4 HilbertSort3d<Hilbert_tcmp, Mesh>(
1017 M, sorted_indices.begin(), sorted_indices.end()
1018
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 );
1019 } else {
1020 ✗ HilbertSort3d<Hilbert_ccmp, Mesh>(
1021 M, sorted_indices.begin(), sorted_indices.end()
1022 ✗ );
1023 }
1024 4 }
1025
1026 /**
1027 * \brief Sorts the vertices of a mesh according to the Morton ordering.
1028 * \details The function does not change the mesh, it computes instead
1029 * the permutation. The permutation can then be reused to order other
1030 * arrays that may depend on the order of the vertices in the mesh (i.e.
1031 * attributes).
1032 * \param[in] M the mesh where the vertices to be sorted reside
1033 * \param[out] sorted_indices the permutation to be applied to the vertices
1034 */
1035 ✗ void morton_vsort_3d(
1036 const Mesh& M, vector<index_t>& sorted_indices
1037 ) {
1038 ✗ sorted_indices.resize(M.vertices.nb());
1039 ✗ for(index_t i: M.vertices) {
1040 ✗ sorted_indices[i] = i;
1041 }
1042 ✗ HilbertSort3d<Morton_vcmp, Mesh>(
1043 M, sorted_indices.begin(), sorted_indices.end()
1044 ✗ );
1045 ✗ }
1046
1047 /**
1048 * \brief Sorts the facets of a mesh according to the Morton ordering.
1049 * \details The function does not change the mesh, it computes instead
1050 * the permutation. The permutation can then be reused to order other
1051 * arrays that may depend on the order of the facets in the mesh (i.e.
1052 * attributes).
1053 * \param[in] M the mesh where the facets to be sorted reside
1054 * \param[out] sorted_indices the permutation to be applied to the facets
1055 */
1056 83 void morton_fsort_3d(
1057 const Mesh& M, vector<index_t>& sorted_indices
1058 ) {
1059 83 sorted_indices.resize(M.facets.nb());
1060
2/2
✓ Branch 4 taken 192734 times.
✓ Branch 5 taken 83 times.
 192817 for(index_t i: M.facets) {
1061
1/2
✓ Branch 1 taken 192734 times.
✗ Branch 2 not taken.
 192734 sorted_indices[i] = i;
1062 }
1063 83 HilbertSort3d<Morton_fcmp, Mesh>(
1064 M, sorted_indices.begin(), sorted_indices.end()
1065
1/2
✓ Branch 3 taken 83 times.
✗ Branch 4 not taken.
 83 );
1066 83 }
1067
1068 /**
1069 * \brief Sorts the cells of a mesh according to the Morton ordering.
1070 * \details The function does not change the mesh, it computes instead
1071 * the permutation. The permutation can then be reused to order other
1072 * arrays that may depend on the order of the tets in the mesh (i.e.
1073 * attributes).
1074 * \param[in] M the mesh where the tets to be sorted reside
1075 * \param[out] sorted_indices the permutation to be applied to the tets
1076 */
1077 ✗ void morton_csort_3d(
1078 const Mesh& M, vector<index_t>& sorted_indices
1079 ) {
1080 ✗ sorted_indices.resize(M.cells.nb());
1081 ✗ for(index_t i: M.cells) {
1082 ✗ sorted_indices[i] = i;
1083 }
1084 ✗ if(M.cells.are_simplices()) {
1085 ✗ HilbertSort3d<Morton_tcmp, Mesh>(
1086 M, sorted_indices.begin(), sorted_indices.end()
1087 ✗ );
1088 } else {
1089 ✗ HilbertSort3d<Morton_ccmp, Mesh>(
1090 M, sorted_indices.begin(), sorted_indices.end()
1091 ✗ );
1092 }
1093 ✗ }
1094
1095 #endif
1096
1097 /**
1098 * \brief Computes the BRIO order for a set of 3D points.
1099 * \details Implementation of compute_BRIO_order().
1100 * It is used to accelerate incremental insertion in Delaunay triangulation
1101 * \param[in] nb_vertices number of vertices to sort
1102 * \param[in] vertices pointer to the coordinates of the vertices
1103 * \param[in] stride number of doubles between two consecutive vertices
1104 * \param[in,out] sorted_indices indices to sort
1105 * \param[in] b iterator to the first index to sort
1106 * \param[in] e iterator one position past the last index to sort
1107 * \param[in] threshold minimum size of interval to be sorted
1108 * \param[in] ratio splitting ratio between current interval and
1109 * the rest to be sorted
1110 * \param[in,out] depth iteration depth
1111 * \param[out] levels if non-null, bounds of each level
1112 */
1113 81 void compute_BRIO_order_recursive(
1114 index_t nb_vertices, const double* vertices,
1115 index_t dimension, index_t stride,
1116 vector<index_t>& sorted_indices,
1117 vector<index_t>::iterator b,
1118 vector<index_t>::iterator e,
1119 index_t threshold,
1120 double ratio,
1121 index_t& depth,
1122 vector<index_t>* levels
1123 ) {
1124
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 81 geo_debug_assert(e > b);
1125
1126 81 vector<index_t>::iterator m = b;
1127
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 41 times.
 81 if(index_t(e - b) > threshold) {
1128 40 ++depth;
1129 40 m = b + signed_index_t(double(e - b) * ratio);
1130
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 compute_BRIO_order_recursive(
1131 nb_vertices, vertices,
1132 dimension, stride,
1133 sorted_indices, b, m,
1134 threshold, ratio, depth,
1135 levels
1136 );
1137 }
1138
1139 81 VertexMesh M(nb_vertices, vertices, stride);
1140
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 69 times.
 81 if(dimension == 3) {
1141 12 HilbertSort3d<Hilbert_vcmp, VertexMesh>(
1142 M, m, e
1143
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 );
1144
1/2
✓ Branch 0 taken 69 times.
✗ Branch 1 not taken.
 69 } else if(dimension ==2) {
1145 69 HilbertSort2d<Hilbert_vcmp, VertexMesh>(
1146 M, m, e
1147
1/2
✓ Branch 1 taken 69 times.
✗ Branch 2 not taken.
 69 );
1148 } else {
1149 ✗ geo_assert_not_reached;
1150 }
1151
1152
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 72 times.
 81 if(levels != nullptr) {
1153
1/2
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
 18 levels->push_back(index_t(e - sorted_indices.begin()));
1154 }
1155 81 }
1156 }
1157
1158 /****************************************************************************/
1159
1160 namespace GEO {
1161
1162 #ifndef GEOGRAM_PSM
1163
1164 10 void mesh_reorder(Mesh& M, MeshOrder order, MeshElementsFlags elements) {
1165
1166
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10 geo_assert(M.vertices.dimension() >= 3);
1167
1168 // Step 1: reorder vertices
1169
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if((elements & MESH_VERTICES) != 0) {
1170 10 vector<index_t> sorted_indices;
1171
1/3
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 10 switch(order) {
1172 10 case MESH_ORDER_HILBERT:
1173
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 hilbert_vsort_3d(M, sorted_indices);
1174 10 break;
1175 ✗ case MESH_ORDER_MORTON:
1176 ✗ morton_vsort_3d(M, sorted_indices);
1177 ✗ break;
1178 }
1179
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 M.vertices.permute_elements(sorted_indices);
1180 10 }
1181
1182 // Step 2: reorder facets
1183
3/6
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
 10 if(((elements & MESH_FACETS) != 0) && (M.facets.nb() != 0)) {
1184 10 vector<index_t> sorted_indices;
1185
1/3
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 10 switch(order) {
1186 10 case MESH_ORDER_HILBERT:
1187
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 hilbert_fsort_3d(M, sorted_indices);
1188 10 break;
1189 ✗ case MESH_ORDER_MORTON:
1190 ✗ morton_fsort_3d(M, sorted_indices);
1191 ✗ break;
1192 }
1193
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 M.facets.permute_elements(sorted_indices);
1194 10 }
1195
1196 // Step 3: reorder cells
1197
5/6
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 6 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 6 times.
 10 if(((elements & MESH_CELLS) != 0) && (M.cells.nb() != 0)) {
1198 4 vector<index_t> sorted_indices;
1199
1/3
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 4 switch(order) {
1200 4 case MESH_ORDER_HILBERT:
1201
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 hilbert_csort_3d(M, sorted_indices);
1202 4 break;
1203 ✗ case MESH_ORDER_MORTON:
1204 ✗ morton_csort_3d(M, sorted_indices);
1205 ✗ break;
1206 }
1207
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 M.cells.permute_elements(sorted_indices);
1208 4 }
1209 10 }
1210
1211 83 void compute_mesh_elements_spatial_order(
1212 const Mesh& M, MeshElementsFlags elements,
1213 vector<index_t>& sorted_indices, MeshOrder order
1214 ) {
1215
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 83 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 83 geo_assert(M.vertices.dimension() >= 3);
1216
1/5
✗ Branch 0 not taken.
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 83 switch(elements) {
1217 ✗ case MESH_VERTICES: {
1218 ✗ switch(order) {
1219 ✗ case MESH_ORDER_HILBERT:
1220 ✗ hilbert_vsort_3d(M, sorted_indices);
1221 ✗ break;
1222 ✗ case MESH_ORDER_MORTON:
1223 ✗ morton_vsort_3d(M, sorted_indices);
1224 ✗ break;
1225 }
1226 ✗ } break;
1227 83 case MESH_FACETS: {
1228
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
 83 switch(order) {
1229 ✗ case MESH_ORDER_HILBERT:
1230 ✗ hilbert_fsort_3d(M, sorted_indices);
1231 ✗ break;
1232 83 case MESH_ORDER_MORTON:
1233 83 morton_fsort_3d(M, sorted_indices);
1234 83 break;
1235 }
1236 83 } break;
1237 ✗ case MESH_CELLS: {
1238 ✗ switch(order) {
1239 ✗ case MESH_ORDER_HILBERT:
1240 ✗ hilbert_csort_3d(M, sorted_indices);
1241 ✗ break;
1242 ✗ case MESH_ORDER_MORTON:
1243 ✗ morton_csort_3d(M, sorted_indices);
1244 ✗ break;
1245 }
1246 ✗ } break;
1247 ✗ case MESH_NONE:
1248 case MESH_EDGES:
1249 case MESH_ALL_ELEMENTS:
1250 case MESH_FACET_CORNERS:
1251 case MESH_CELL_CORNERS:
1252 case MESH_CELL_FACETS:
1253 case MESH_ALL_SUBELEMENTS:
1254 ✗ geo_assert_not_reached;
1255 }
1256 83 }
1257 #endif
1258
1259
1260 ✗ void compute_Hilbert_order(
1261 index_t total_nb_vertices, const double* vertices,
1262 vector<index_t>& sorted_indices,
1263 index_t first,
1264 index_t last,
1265 index_t dimension, index_t stride
1266 ) {
1267 ✗ geo_debug_assert(last > first);
1268 ✗ if(last - first <= 1) {
1269 ✗ return;
1270 }
1271 ✗ VertexMesh M(total_nb_vertices, vertices, stride);
1272 ✗ if(dimension == 3) {
1273 ✗ HilbertSort3d<Hilbert_vcmp, VertexMesh>(
1274 ✗ M, sorted_indices.begin() + int(first),
1275 ✗ sorted_indices.begin() + int(last)
1276 ✗ );
1277 ✗ } else if(dimension == 2) {
1278 ✗ HilbertSort2d<Hilbert_vcmp, VertexMesh>(
1279 ✗ M, sorted_indices.begin() + int(first),
1280 ✗ sorted_indices.begin() + int(last)
1281 ✗ );
1282 } else {
1283 ✗ geo_assert_not_reached;
1284 }
1285 }
1286
1287 41 void compute_BRIO_order(
1288 index_t nb_vertices, const double* vertices,
1289 vector<index_t>& sorted_indices,
1290 index_t dimension,
1291 index_t stride,
1292 index_t threshold,
1293 double ratio,
1294 vector<index_t>* levels
1295 ) {
1296
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 36 times.
 41 if(levels != nullptr) {
1297 5 levels->clear();
1298
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 levels->push_back(0);
1299 }
1300 41 index_t depth = 0;
1301
1/2
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
 41 sorted_indices.resize(nb_vertices);
1302
2/2
✓ Branch 0 taken 227234 times.
✓ Branch 1 taken 41 times.
 227275 for(index_t i = 0; i < nb_vertices; ++i) {
1303
1/2
✓ Branch 1 taken 227234 times.
✗ Branch 2 not taken.
 227234 sorted_indices[i] = i;
1304 }
1305
1306
1/2
✓ Branch 3 taken 41 times.
✗ Branch 4 not taken.
 41 GEO::random_shuffle(sorted_indices.begin(), sorted_indices.end());
1307
1308
1/2
✓ Branch 3 taken 41 times.
✗ Branch 4 not taken.
 41 compute_BRIO_order_recursive(
1309 nb_vertices, vertices,
1310 dimension, stride,
1311 sorted_indices,
1312 sorted_indices.begin(), sorted_indices.end(),
1313 threshold, ratio, depth, levels
1314 );
1315 41 }
1316 }
1317
1318 /**********************************************************************/
1319
1320 namespace {
1321 using namespace GEO;
1322
1323 // Same as in delaunay/periodic.cpp,
1324 // copied here for now because linker does not find
1325 // it under Android.
1326 int Periodic_translation[27][3] = {
1327 { 0, 0, 0}, //13 -> 0 + <-- zero displacement is first.
1328 { -1, -1, -1}, //0 -> 1 -
1329 { -1, -1, 0}, //1 -> 2 -
1330 { -1, -1, 1}, //2 -> 3 -
1331 { -1, 0, -1}, //3 -> 4 -
1332 { -1, 0, 0}, //4 -> 5 -
1333 { -1, 0, 1}, //5 -> 6 -
1334 { -1, 1, -1}, //6 -> 7 -
1335 { -1, 1, 0}, //7 -> 8 -
1336 { -1, 1, 1}, //8 -> 9 -
1337 { 0, -1, -1}, //9 -> 10 -
1338 { 0, -1, 0}, //10 -> 11 -
1339 { 0, -1, 1}, //11 -> 12 -
1340 { 0, 0, -1}, //12 -> 13 -
1341 // (zero displacement was there)
1342 { 0, 0, 1}, //14 -> 14 +
1343 { 0, 1, -1}, //15 -> 15 -
1344 { 0, 1, 0}, //16 -> 16 +
1345 { 0, 1, 1}, //17 -> 17 +
1346 { 1, -1, -1}, //18 -> 18 -
1347 { 1, -1, 0}, //19 -> 19 -
1348 { 1, -1, 1}, //20 -> 20 -
1349 { 1, 0, -1}, //21 -> 21 -
1350 { 1, 0, 0}, //22 -> 22 +
1351 { 1, 0, 1}, //23 -> 23 +
1352 { 1, 1, -1}, //24 -> 24 -
1353 { 1, 1, 0}, //25 -> 25 +
1354 { 1, 1, 1} //26 -> 26 +
1355 };
1356
1357
1358 /**
1359 * \details Exposes an interface compatible with the requirement
1360 * of Hilbert sort templates for a raw array of vertices.
1361 */
1362 class PeriodicVertexArray3d {
1363 public:
1364 /**
1365 * \brief Constructs a new PeriodicVertexArray.
1366 * \param[in] nb_vertices total number of vertices, including
1367 * the 27 copies.
1368 * \param[in] base address of the points.
1369 * \param[in] stride number of doubles between
1370 * two consecutive points.
1371 * \param[in] period the edge length of the periodic domain.
1372 */
1373 ✗ PeriodicVertexArray3d(
1374 index_t nb_vertices,
1375 const double* base, index_t stride,
1376 const vec3& period
1377 ✗ ) :
1378 ✗ base_(base),
1379 ✗ stride_(stride) {
1380 ✗ nb_vertices_ = nb_vertices;
1381 ✗ nb_real_vertices_ = nb_vertices_ / 27;
1382 ✗ geo_debug_assert(nb_vertices % 27 == 0);
1383
1384
1385 ✗ for(index_t i=0; i<27; ++i) {
1386 ✗ for(index_t j=0; j<3; ++j) {
1387 ✗ xlat_[i][j] = period[j] * double(Periodic_translation[i][j]);
1388 }
1389 }
1390 ✗ }
1391
1392 /**
1393 * \brief Gets a point coordinate by its index and coordinate.
1394 * \param[in] i the index of the point.
1395 * \param[in] coord the coordinate.
1396 * \return the value of the coordinate.
1397 */
1398 ✗ double point_coord(index_t i, index_t coord) const {
1399 ✗ index_t instance = i / nb_real_vertices_;
1400 ✗ i = i % nb_real_vertices_;
1401 ✗ return (base_ + i * stride_)[coord] + xlat_[instance][coord];
1402 }
1403
1404 private:
1405 const double* base_;
1406 index_t stride_;
1407 index_t nb_vertices_;
1408 index_t nb_real_vertices_;
1409 double xlat_[27][3];
1410 };
1411
1412 /**
1413 * \brief Exposes an interface compatible with the requirement
1414 * of Hilbert sort templates for a raw array of vertices.
1415 */
1416 class PeriodicVertexMesh3d {
1417 public:
1418 /**
1419 * \brief Constructs a new VertexMesh.
1420 * \param[in] nb_vertices total number of vertices, including
1421 * the 27 copies.
1422 * \param[in] base address of the points
1423 * \param[in] stride number of doubles between
1424 * two consecutive points
1425 */
1426 ✗ PeriodicVertexMesh3d(
1427 index_t nb_vertices,
1428 const double* base, index_t stride, const vec3& period
1429 ✗ ) : vertices(nb_vertices, base, stride, period) {
1430 ✗ }
1431
1432 PeriodicVertexArray3d vertices;
1433 };
1434
1435 /**
1436 * \brief Drop-in replacement of Hilbert_vcmp for
1437 * periodic vertices.
1438 * \details Needed because point_ptr() is not defined,
1439 * we need to use point_coord() instead.
1440 */
1441 template <int COORD, bool UP, class MESH>
1442 struct Hilbert_vcmp_periodic {
1443 };
1444
1445 /**
1446 * \brief Drop-in replacement of Hilbert_vcmp for
1447 * periodic vertices, specialization for UP.
1448 * \details Needed because point_ptr() is not defined,
1449 * we need to use point_coord() instead.
1450 */
1451 template <int COORD>
1452 struct Hilbert_vcmp_periodic<COORD, true, PeriodicVertexMesh3d> {
1453 ✗ Hilbert_vcmp_periodic(const PeriodicVertexMesh3d& mesh) :
1454 ✗ mesh_(mesh) {
1455 ✗ }
1456 ✗ bool operator() (index_t i1, index_t i2) {
1457 return
1458 ✗ mesh_.vertices.point_coord(i1,COORD) <
1459 ✗ mesh_.vertices.point_coord(i2,COORD);
1460 }
1461 const PeriodicVertexMesh3d& mesh_;
1462 };
1463
1464 /**
1465 * \brief Drop-in replacement of Hilbert_vcmp for
1466 * periodic vertices, specialization for !UP.
1467 * \details Needed because point_ptr() is not defined,
1468 * we need to use point_coord() instead.
1469 */
1470 template <int COORD>
1471 struct Hilbert_vcmp_periodic<COORD, false, PeriodicVertexMesh3d> {
1472 ✗ Hilbert_vcmp_periodic(const PeriodicVertexMesh3d& mesh) :
1473 ✗ mesh_(mesh) {
1474 ✗ }
1475 ✗ bool operator() (index_t i1, index_t i2) {
1476 return
1477 ✗ mesh_.vertices.point_coord(i1,COORD) >
1478 ✗ mesh_.vertices.point_coord(i2,COORD);
1479 }
1480 const PeriodicVertexMesh3d& mesh_;
1481 };
1482
1483 }
1484
1485 namespace GEO {
1486
1487 ✗ void Hilbert_sort_periodic(
1488 index_t nb_vertices, const double* vertices,
1489 vector<index_t>& sorted_indices,
1490 index_t dimension,
1491 index_t stride,
1492 vector<index_t>::iterator b,
1493 vector<index_t>::iterator e,
1494 const vec3& period
1495 ) {
1496 ✗ geo_assert(dimension == 3); // Only implemented for 3D.
1497 ✗ geo_argused(sorted_indices); // Accessed through b and e.
1498
1499 ✗ GEO::random_shuffle(b,e);
1500
1501 ✗ PeriodicVertexMesh3d M(nb_vertices, vertices, stride, period);
1502 ✗ HilbertSort3d<Hilbert_vcmp_periodic, PeriodicVertexMesh3d>(
1503 M, b, e
1504 ✗ );
1505 ✗ }
1506
1507 }
1508