TTK
Loading...
Searching...
No Matches
ImplicitTriangulation.cpp
Go to the documentation of this file.
2
3#include <numeric>
4
5using namespace std;
6using namespace ttk;
7
8#define CASE_EDGE_POSITION_L_3D \
9 case EdgePosition::L_xnn_3D: \
10 case EdgePosition::L_xn0_3D: \
11 case EdgePosition::L_xnN_3D: \
12 case EdgePosition::L_x0n_3D: \
13 case EdgePosition::L_x00_3D: \
14 case EdgePosition::L_x0N_3D: \
15 case EdgePosition::L_xNn_3D: \
16 case EdgePosition::L_xN0_3D: \
17 case EdgePosition::L_xNN_3D
18#define CASE_EDGE_POSITION_H_3D \
19 case EdgePosition::H_nyn_3D: \
20 case EdgePosition::H_ny0_3D: \
21 case EdgePosition::H_nyN_3D: \
22 case EdgePosition::H_0yn_3D: \
23 case EdgePosition::H_0y0_3D: \
24 case EdgePosition::H_0yN_3D: \
25 case EdgePosition::H_Nyn_3D: \
26 case EdgePosition::H_Ny0_3D: \
27 case EdgePosition::H_NyN_3D
28#define CASE_EDGE_POSITION_P_3D \
29 case EdgePosition::P_nnz_3D: \
30 case EdgePosition::P_n0z_3D: \
31 case EdgePosition::P_nNz_3D: \
32 case EdgePosition::P_0nz_3D: \
33 case EdgePosition::P_00z_3D: \
34 case EdgePosition::P_0Nz_3D: \
35 case EdgePosition::P_Nnz_3D: \
36 case EdgePosition::P_N0z_3D: \
37 case EdgePosition::P_NNz_3D
38#define CASE_EDGE_POSITION_D1_3D \
39 case EdgePosition::D1_xyn_3D: \
40 case EdgePosition::D1_xy0_3D: \
41 case EdgePosition::D1_xyN_3D
42#define CASE_EDGE_POSITION_D2_3D \
43 case EdgePosition::D2_nyz_3D: \
44 case EdgePosition::D2_0yz_3D: \
45 case EdgePosition::D2_Nyz_3D
46#define CASE_EDGE_POSITION_D3_3D \
47 case EdgePosition::D3_xnz_3D: \
48 case EdgePosition::D3_x0z_3D: \
49 case EdgePosition::D3_xNz_3D
50#define CASE_EDGE_POSITION_L_2D \
51 case EdgePosition::L_xn_2D: \
52 case EdgePosition::L_x0_2D: \
53 case EdgePosition::L_xN_2D
54#define CASE_EDGE_POSITION_H_2D \
55 case EdgePosition::H_ny_2D: \
56 case EdgePosition::H_0y_2D: \
57 case EdgePosition::H_Ny_2D
58
60 : cellNumber_{}, vertexNumber_{}, edgeNumber_{}, triangleNumber_{},
61 tetrahedronNumber_{}, isAccelerated_{} {
62 setDebugMsgPrefix("ImplicitTriangulation");
63#ifdef TTK_ENABLE_MPI
64 this->hasPreconditionedDistributedEdges_ = true;
65 this->hasPreconditionedDistributedTriangles_ = true;
66#endif // TTK_ENABLE_MPI
67}
68
70
71int ImplicitTriangulation::setInputGrid(const float &xOrigin,
72 const float &yOrigin,
73 const float &zOrigin,
74 const float &xSpacing,
75 const float &ySpacing,
76 const float &zSpacing,
77 const SimplexId &xDim,
78 const SimplexId &yDim,
79 const SimplexId &zDim) {
80
81 // Dimensionality //
82 if(xDim < 1 or yDim < 1 or zDim < 1)
83 dimensionality_ = -1;
84 else if(xDim > 1 and yDim > 1 and zDim > 1)
86 else if((xDim > 1 and yDim > 1) or (yDim > 1 and zDim > 1)
87 or (xDim > 1 and zDim > 1))
89 else if(xDim > 1 or yDim > 1 or zDim > 1)
91 else
93
94 // Essentials //
95 origin_[0] = xOrigin;
96 origin_[1] = yOrigin;
97 origin_[2] = zOrigin;
98 spacing_[0] = xSpacing;
99 spacing_[1] = ySpacing;
100 spacing_[2] = zSpacing;
101 dimensions_[0] = xDim;
102 dimensions_[1] = yDim;
103 dimensions_[2] = zDim;
104 nbvoxels_[0] = xDim - 1;
105 nbvoxels_[1] = yDim - 1;
106 nbvoxels_[2] = zDim - 1;
107
108 if(dimensionality_ == 3) {
109 // VertexShift
110 vshift_[0] = xDim;
111 vshift_[1] = xDim * yDim;
112 // EdgeSetDimensions
113 esetdims_[0] = (xDim - 1) * yDim * zDim;
114 esetdims_[1] = xDim * (yDim - 1) * zDim;
115 esetdims_[2] = xDim * yDim * (zDim - 1);
116 esetdims_[3] = (xDim - 1) * (yDim - 1) * zDim;
117 esetdims_[4] = xDim * (yDim - 1) * (zDim - 1);
118 esetdims_[5] = (xDim - 1) * yDim * (zDim - 1);
119 esetdims_[6] = (xDim - 1) * (yDim - 1) * (zDim - 1);
120 // EdgeSetShift
121 esetshift_[0] = esetdims_[0];
122 for(int k = 1; k < 7; ++k)
123 esetshift_[k] = esetshift_[k - 1] + esetdims_[k];
124 // EdgeShift
125 eshift_[0] = xDim - 1;
126 eshift_[1] = (xDim - 1) * yDim;
127 eshift_[2] = xDim;
128 eshift_[3] = xDim * (yDim - 1);
129 eshift_[4] = xDim;
130 eshift_[5] = xDim * yDim;
131 eshift_[6] = xDim - 1;
132 eshift_[7] = (xDim - 1) * (yDim - 1);
133 eshift_[8] = xDim;
134 eshift_[9] = xDim * (yDim - 1);
135 eshift_[10] = xDim - 1;
136 eshift_[11] = (xDim - 1) * yDim;
137 eshift_[12] = xDim - 1;
138 eshift_[13] = (xDim - 1) * (yDim - 1);
139 // TriangleSetDimensions
140 tsetdims_[0] = (xDim - 1) * (yDim - 1) * zDim * 2;
141 tsetdims_[1] = (xDim - 1) * yDim * (zDim - 1) * 2;
142 tsetdims_[2] = xDim * (yDim - 1) * (zDim - 1) * 2;
143 tsetdims_[3] = (xDim - 1) * (yDim - 1) * (zDim - 1) * 2;
144 tsetdims_[4] = (xDim - 1) * (yDim - 1) * (zDim - 1) * 2;
145 tsetdims_[5] = (xDim - 1) * (yDim - 1) * (zDim - 1) * 2;
146 // TriangleSetShift
147 tsetshift_[0] = tsetdims_[0];
148 for(int k = 1; k < 6; ++k)
149 tsetshift_[k] = tsetshift_[k - 1] + tsetdims_[k];
150 // TriangleShift
151 tshift_[0] = (xDim - 1) * 2;
152 tshift_[1] = (xDim - 1) * (yDim - 1) * 2;
153 tshift_[2] = (xDim - 1) * 2;
154 tshift_[3] = (xDim - 1) * yDim * 2;
155 tshift_[4] = xDim * 2;
156 tshift_[5] = xDim * (yDim - 1) * 2;
157 tshift_[6] = (xDim - 1) * 2;
158 tshift_[7] = (xDim - 1) * (yDim - 1) * 2;
159 tshift_[8] = (xDim - 1) * 2;
160 tshift_[9] = (xDim - 1) * (yDim - 1) * 2;
161 tshift_[10] = (xDim - 1) * 2;
162 tshift_[11] = (xDim - 1) * (yDim - 1) * 2;
163 // TetrahedronShift
164 tetshift_[0] = (xDim - 1) * 6;
165 tetshift_[1] = (xDim - 1) * (yDim - 1) * 6;
166
167 // Numbers
168 vertexNumber_ = xDim * yDim * zDim;
169 edgeNumber_ = 0;
170 for(int k = 0; k < 7; ++k)
172 triangleNumber_ = 0;
173 for(int k = 0; k < 6; ++k)
175 tetrahedronNumber_ = (xDim - 1) * (yDim - 1) * (zDim - 1) * 6;
177
179 } else if(dimensionality_ == 2) {
180 // dimensions selectors
181 if(xDim == 1) {
182 Di_ = 1;
183 Dj_ = 2;
184 } else if(yDim == 1) {
185 Di_ = 0;
186 Dj_ = 2;
187 } else {
188 Di_ = 0;
189 Dj_ = 1;
190 }
191 // VertexShift
192 vshift_[0] = dimensions_[Di_];
193 // EdgeSetDimensions
194 esetdims_[0] = (dimensions_[Di_] - 1) * dimensions_[Dj_];
195 esetdims_[1] = dimensions_[Di_] * (dimensions_[Dj_] - 1);
196 esetdims_[2] = (dimensions_[Di_] - 1) * (dimensions_[Dj_] - 1);
197 // EdgeSetShift
198 esetshift_[0] = esetdims_[0];
199 for(int k = 1; k < 3; ++k)
200 esetshift_[k] = esetshift_[k - 1] + esetdims_[k];
201 // EdgeShift
202 eshift_[0] = dimensions_[Di_] - 1;
203 eshift_[2] = dimensions_[Di_];
204 eshift_[4] = dimensions_[Di_] - 1;
205 // TriangleShift
206 tshift_[0] = (dimensions_[Di_] - 1) * 2;
207
208 // Numbers
210 edgeNumber_ = 0;
211 for(int k = 0; k < 3; ++k)
213 triangleNumber_ = (dimensions_[Di_] - 1) * (dimensions_[Dj_] - 1) * 2;
215
217 } else if(dimensionality_ == 1) {
218 // dimensions selectors
219 for(int k = 0; k < 3; ++k) {
220 if(dimensions_[k] > 1) {
221 Di_ = k;
222 break;
223 }
224 }
225
226 // Numbers
230 }
231
232 return 0;
233}
234
236 isAccelerated_ = false;
237
238 unsigned long long int msb[3];
239 if(dimensionality_ == 3) {
240 bool allDimensionsArePowerOfTwo = true;
241 for(int k = 0; k < 3; ++k)
242 if(!isPowerOfTwo(dimensions_[k], msb[k]))
243 allDimensionsArePowerOfTwo = false;
244
245 if(allDimensionsArePowerOfTwo) {
246 mod_[0] = dimensions_[0] - 1;
247 mod_[1] = dimensions_[0] * dimensions_[1] - 1;
248 div_[0] = msb[0];
249 div_[1] = msb[0] + msb[1];
250 isAccelerated_ = true;
251 }
252 } else if(dimensionality_ == 2) {
253 bool const isDi = isPowerOfTwo(dimensions_[Di_], msb[Di_]);
254 bool const isDj = isPowerOfTwo(dimensions_[Dj_], msb[Dj_]);
255 bool const allDimensionsArePowerOfTwo = (isDi and isDj);
256
257 if(allDimensionsArePowerOfTwo) {
258 mod_[0] = dimensions_[Di_] - 1;
259 div_[0] = msb[Di_];
260 isAccelerated_ = true;
261 }
262 }
263
264 if(isAccelerated_) {
265 printMsg("Accelerated getVertex*() requests.", debug::Priority::INFO);
266 }
267
268 return 0;
269}
270
271bool ImplicitTriangulation::isPowerOfTwo(unsigned long long int v,
272 unsigned long long int &r) {
273 if(v && !(v & (v - 1))) {
274 r = 0;
275 while(v >>= 1)
276 r++;
277 return true;
278 }
279 return false;
280}
281
282template <typename Derived>
284 isVertexOnBoundary)(const SimplexId &vertexId) const {
285
286#ifdef TTK_ENABLE_MPI
287 if(this->metaGrid_ != nullptr) {
288 return this->isVertexOnGlobalBoundaryInternal(vertexId);
289 }
290#endif // TTK_ENABLE_MPI
291
292#ifndef TTK_ENABLE_KAMIKAZE
293 if(vertexId < 0 or vertexId >= vertexNumber_)
294 return false;
295#endif // !TTK_ENABLE_KAMIKAZE
296
297 switch(this->underlying().getVertexPosition(vertexId)) {
298 case VertexPosition::CENTER_3D:
299 case VertexPosition::CENTER_2D:
300 case VertexPosition::CENTER_1D:
301 return false;
302 default:
303 return true;
304 }
305}
306
307template <typename Derived>
309 isEdgeOnBoundary)(const SimplexId &edgeId) const {
310
311#ifdef TTK_ENABLE_MPI
312 if(this->metaGrid_ != nullptr) {
313 return this->isEdgeOnGlobalBoundaryInternal(edgeId);
314 }
315#endif // TTK_ENABLE_MPI
316
317#ifndef TTK_ENABLE_KAMIKAZE
318 if(edgeId < 0 or edgeId >= edgeNumber_)
319 return false;
320#endif // !TTK_ENABLE_KAMIKAZE
321
322 switch(this->underlying().getEdgePosition(edgeId)) {
323 case EdgePosition::L_xnn_3D:
324 case EdgePosition::H_nyn_3D:
325 case EdgePosition::P_nnz_3D:
326 case EdgePosition::D1_xyn_3D:
327 case EdgePosition::D2_nyz_3D:
328 case EdgePosition::D3_xnz_3D:
329 case EdgePosition::D4_3D:
330 case EdgePosition::L_xn_2D:
331 case EdgePosition::H_ny_2D:
332 case EdgePosition::D1_2D:
333 return false;
334 default:
335 break;
336 }
337 return true;
338}
339
340bool ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(isTriangleOnBoundary)(
341 const SimplexId &triangleId) const {
342
343#ifdef TTK_ENABLE_MPI
344 if(this->metaGrid_ != nullptr) {
345 return this->isTriangleOnGlobalBoundaryInternal(triangleId);
346 }
347#endif // TTK_ENABLE_MPI
348
349#ifndef TTK_ENABLE_KAMIKAZE
350 if(triangleId < 0 or triangleId >= triangleNumber_)
351 return false;
352#endif // !TTK_ENABLE_KAMIKAZE
353
354 if(dimensionality_ == 3)
355 return (TTK_TRIANGULATION_INTERNAL(getTriangleStarNumber)(triangleId) == 1);
356
357 return false;
358}
359
360const vector<vector<SimplexId>> *
361 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getVertexNeighbors)() {
362 if(vertexNeighborList_.empty()) {
363 Timer t;
364 vertexNeighborList_.resize(vertexNumber_);
365 for(SimplexId i = 0; i < vertexNumber_; ++i) {
366 vertexNeighborList_[i].resize(getVertexNeighborNumber(i));
367 for(SimplexId j = 0; j < (SimplexId)vertexNeighborList_[i].size(); ++j)
368 getVertexNeighbor(i, j, vertexNeighborList_[i][j]);
369 }
370
371 printMsg("Built " + to_string(vertexNumber_) + " vertex neighbors.", 1,
372 t.getElapsedTime(), 1);
373 }
374
375 return &vertexNeighborList_;
376}
377
382
383template <typename Derived>
385 const SimplexId &vertexId, const int &localEdgeId, SimplexId &edgeId) const {
386#ifndef TTK_ENABLE_KAMIKAZE
387 if(localEdgeId < 0 or localEdgeId >= getVertexEdgeNumberInternal(vertexId))
388 return -1;
389#endif
390 // e--------f
391 // /| /|
392 // / | / |
393 // a--g-----b--h
394 // | / | /
395 // |/ |/
396 // c--------d
397 //
398 // Classement des "Edges" et dans cet ordre:
399 // L: largeur (type ab)
400 // H: hauteur (type ac)
401 // P: profondeur (type ae)
402 // D1: diagonale1 (type bc)
403 // D2: diagonale2 (type ag)
404 // D3: diagonale3 (type be)
405 // D4: diagonale4 (type bg)
406
407 const auto &p = this->underlying().getVertexCoords(vertexId);
408
409 switch(this->underlying().getVertexPosition(vertexId)) {
410 case VertexPosition::CENTER_3D:
411 edgeId = getVertexEdgeABCDEFGH(p.data(), localEdgeId);
412 break;
413 case VertexPosition::FRONT_FACE_3D:
414 edgeId = getVertexEdgeABDC(p.data(), localEdgeId);
415 break;
416 case VertexPosition::BACK_FACE_3D:
417 edgeId = getVertexEdgeEFHG(p.data(), localEdgeId);
418 break;
419 case VertexPosition::TOP_FACE_3D:
420 edgeId = getVertexEdgeAEFB(p.data(), localEdgeId);
421 break;
422 case VertexPosition::BOTTOM_FACE_3D:
423 edgeId = getVertexEdgeGHDC(p.data(), localEdgeId);
424 break;
425 case VertexPosition::LEFT_FACE_3D:
426 edgeId = getVertexEdgeAEGC(p.data(), localEdgeId);
427 break;
428 case VertexPosition::RIGHT_FACE_3D:
429 edgeId = getVertexEdgeBFHD(p.data(), localEdgeId);
430 break;
431 case VertexPosition::TOP_FRONT_EDGE_3D: // ab
432 edgeId = getVertexEdgeAB(p.data(), localEdgeId);
433 break;
434 case VertexPosition::BOTTOM_FRONT_EDGE_3D: // cd
435 edgeId = getVertexEdgeCD(p.data(), localEdgeId);
436 break;
437 case VertexPosition::LEFT_FRONT_EDGE_3D: // ac
438 edgeId = getVertexEdgeAC(p.data(), localEdgeId);
439 break;
440 case VertexPosition::RIGHT_FRONT_EDGE_3D: // bd
441 edgeId = getVertexEdgeBD(p.data(), localEdgeId);
442 break;
443 case VertexPosition::TOP_BACK_EDGE_3D: // ef
444 edgeId = getVertexEdgeEF(p.data(), localEdgeId);
445 break;
446 case VertexPosition::BOTTOM_BACK_EDGE_3D: // gh
447 edgeId = getVertexEdgeGH(p.data(), localEdgeId);
448 break;
449 case VertexPosition::LEFT_BACK_EDGE_3D: // eg
450 edgeId = getVertexEdgeEG(p.data(), localEdgeId);
451 break;
452 case VertexPosition::RIGHT_BACK_EDGE_3D: // fh
453 edgeId = getVertexEdgeFH(p.data(), localEdgeId);
454 break;
455 case VertexPosition::TOP_LEFT_EDGE_3D: // ae
456 edgeId = getVertexEdgeAE(p.data(), localEdgeId);
457 break;
458 case VertexPosition::TOP_RIGHT_EDGE_3D: // bf
459 edgeId = getVertexEdgeBF(p.data(), localEdgeId);
460 break;
461 case VertexPosition::BOTTOM_LEFT_EDGE_3D: // cg
462 edgeId = getVertexEdgeCG(p.data(), localEdgeId);
463 break;
464 case VertexPosition::BOTTOM_RIGHT_EDGE_3D: // dh
465 edgeId = getVertexEdgeDH(p.data(), localEdgeId);
466 break;
467 case VertexPosition::TOP_LEFT_FRONT_CORNER_3D: // a
468 edgeId = getVertexEdgeA(p.data(), localEdgeId);
469 break;
470 case VertexPosition::TOP_RIGHT_FRONT_CORNER_3D: // b
471 edgeId = getVertexEdgeB(p.data(), localEdgeId);
472 break;
473 case VertexPosition::BOTTOM_LEFT_FRONT_CORNER_3D: // c
474 edgeId = getVertexEdgeC(p.data(), localEdgeId);
475 break;
476 case VertexPosition::BOTTOM_RIGHT_FRONT_CORNER_3D: // d
477 edgeId = getVertexEdgeD(p.data(), localEdgeId);
478 break;
479 case VertexPosition::TOP_LEFT_BACK_CORNER_3D: // e
480 edgeId = getVertexEdgeE(p.data(), localEdgeId);
481 break;
482 case VertexPosition::TOP_RIGHT_BACK_CORNER_3D: // f
483 edgeId = getVertexEdgeF(p.data(), localEdgeId);
484 break;
485 case VertexPosition::BOTTOM_LEFT_BACK_CORNER_3D: // g
486 edgeId = getVertexEdgeG(p.data(), localEdgeId);
487 break;
488 case VertexPosition::BOTTOM_RIGHT_BACK_CORNER_3D: // h
489 edgeId = getVertexEdgeH(p.data(), localEdgeId);
490 break;
491 case VertexPosition::CENTER_2D:
492 edgeId = getVertexEdge2dABCD(p.data(), localEdgeId);
493 break;
494 case VertexPosition::TOP_EDGE_2D:
495 edgeId = getVertexEdge2dAB(p.data(), localEdgeId);
496 break;
497 case VertexPosition::BOTTOM_EDGE_2D:
498 edgeId = getVertexEdge2dCD(p.data(), localEdgeId);
499 break;
500 case VertexPosition::LEFT_EDGE_2D:
501 edgeId = getVertexEdge2dAC(p.data(), localEdgeId);
502 break;
503 case VertexPosition::RIGHT_EDGE_2D:
504 edgeId = getVertexEdge2dBD(p.data(), localEdgeId);
505 break;
506 case VertexPosition::TOP_LEFT_CORNER_2D: // a
507 edgeId = getVertexEdge2dA(p.data(), localEdgeId);
508 break;
509 case VertexPosition::TOP_RIGHT_CORNER_2D: // b
510 edgeId = getVertexEdge2dB(p.data(), localEdgeId);
511 break;
512 case VertexPosition::BOTTOM_LEFT_CORNER_2D: // c
513 edgeId = getVertexEdge2dC(p.data(), localEdgeId);
514 break;
515 case VertexPosition::BOTTOM_RIGHT_CORNER_2D: // d
516 edgeId = getVertexEdge2dD(p.data(), localEdgeId);
517 break;
518 case VertexPosition::CENTER_1D:
519 edgeId = (localEdgeId == 0 ? vertexId : vertexId - 1);
520 break;
521 case VertexPosition::LEFT_CORNER_1D:
522 edgeId = vertexId;
523 break;
524 case VertexPosition::RIGHT_CORNER_1D:
525 edgeId = vertexId - 1;
526 break;
527 default:
528 edgeId = -1;
529 }
530
531 return 0;
532}
533
534const vector<vector<SimplexId>> *
536 if(vertexEdgeList_.empty()) {
537 Timer t;
538
540 for(SimplexId i = 0; i < vertexNumber_; ++i) {
542 for(SimplexId j = 0; j < (SimplexId)vertexEdgeList_[i].size(); ++j)
544 }
545
546 printMsg("Built " + to_string(vertexNumber_) + " vertex edges.", 1,
547 t.getElapsedTime(), 1);
548 }
549
550 return &vertexEdgeList_;
551}
552
553template <typename Derived>
555 const SimplexId &vertexId) const {
556#ifndef TTK_ENABLE_KAMIKAZE
557 if(vertexId < 0 or vertexId >= vertexNumber_)
558 return -1;
559#endif
560
561 switch(this->underlying().getVertexPosition(vertexId)) {
562 case VertexPosition::CENTER_3D:
563 return 36;
564 case VertexPosition::FRONT_FACE_3D:
565 case VertexPosition::BACK_FACE_3D:
566 case VertexPosition::TOP_FACE_3D:
567 case VertexPosition::BOTTOM_FACE_3D:
568 case VertexPosition::LEFT_FACE_3D:
569 case VertexPosition::RIGHT_FACE_3D:
570 return 21;
571 case VertexPosition::TOP_FRONT_EDGE_3D: // ab
572 case VertexPosition::RIGHT_FRONT_EDGE_3D: // bd
573 case VertexPosition::BOTTOM_BACK_EDGE_3D: // gh
574 case VertexPosition::LEFT_BACK_EDGE_3D: // eg
575 case VertexPosition::BOTTOM_LEFT_EDGE_3D: // cg
576 case VertexPosition::TOP_RIGHT_EDGE_3D: // bf
577 return 15;
578 case VertexPosition::TOP_RIGHT_FRONT_CORNER_3D: // b
579 case VertexPosition::BOTTOM_LEFT_BACK_CORNER_3D: // g
580 return 12;
581 case VertexPosition::TOP_BACK_EDGE_3D: // ef
582 case VertexPosition::BOTTOM_FRONT_EDGE_3D: // cd
583 case VertexPosition::LEFT_FRONT_EDGE_3D: // ac
584 case VertexPosition::TOP_LEFT_EDGE_3D: // ae
585 case VertexPosition::RIGHT_BACK_EDGE_3D: // fh
586 case VertexPosition::BOTTOM_RIGHT_EDGE_3D: // dh
587 return 9;
588 case VertexPosition::TOP_LEFT_FRONT_CORNER_3D: // a
589 case VertexPosition::BOTTOM_LEFT_FRONT_CORNER_3D: // c
590 case VertexPosition::BOTTOM_RIGHT_FRONT_CORNER_3D: // d
591 case VertexPosition::TOP_LEFT_BACK_CORNER_3D: // e
592 case VertexPosition::TOP_RIGHT_BACK_CORNER_3D: // f
593 case VertexPosition::BOTTOM_RIGHT_BACK_CORNER_3D: // h
594 return 5;
595 default: // 1D + 2D
596 break;
597 }
598
599 return 0;
600}
601
602template <typename Derived>
604 const SimplexId &vertexId,
605 const int &localTriangleId,
606 SimplexId &triangleId) const {
607#ifndef TTK_ENABLE_KAMIKAZE
608 if(localTriangleId < 0
609 or localTriangleId >= getVertexTriangleNumberInternal(vertexId))
610 return -1;
611#endif
612
613 const auto &p = this->underlying().getVertexCoords(vertexId);
614
615 switch(this->underlying().getVertexPosition(vertexId)) {
616 case VertexPosition::CENTER_3D:
617 triangleId = getVertexTriangleABCDEFGH(p.data(), localTriangleId);
618 break;
619 case VertexPosition::FRONT_FACE_3D:
620 triangleId = getVertexTriangleABDC(p.data(), localTriangleId);
621 break;
622 case VertexPosition::BACK_FACE_3D:
623 triangleId = getVertexTriangleEFHG(p.data(), localTriangleId);
624 break;
625 case VertexPosition::TOP_FACE_3D:
626 triangleId = getVertexTriangleAEFB(p.data(), localTriangleId);
627 break;
628 case VertexPosition::BOTTOM_FACE_3D:
629 triangleId = getVertexTriangleGHDC(p.data(), localTriangleId);
630 break;
631 case VertexPosition::LEFT_FACE_3D:
632 triangleId = getVertexTriangleAEGC(p.data(), localTriangleId);
633 break;
634 case VertexPosition::RIGHT_FACE_3D:
635 triangleId = getVertexTriangleBFHD(p.data(), localTriangleId);
636 break;
637 case VertexPosition::TOP_FRONT_EDGE_3D: // ab
638 triangleId = getVertexTriangleAB(p.data(), localTriangleId);
639 break;
640 case VertexPosition::BOTTOM_FRONT_EDGE_3D: // cd
641 triangleId = getVertexTriangleCD(p.data(), localTriangleId);
642 break;
643 case VertexPosition::LEFT_FRONT_EDGE_3D: // ac
644 triangleId = getVertexTriangleAC(p.data(), localTriangleId);
645 break;
646 case VertexPosition::RIGHT_FRONT_EDGE_3D: // bd
647 triangleId = getVertexTriangleBD(p.data(), localTriangleId);
648 break;
649 case VertexPosition::TOP_BACK_EDGE_3D: // ef
650 triangleId = getVertexTriangleEF(p.data(), localTriangleId);
651 break;
652 case VertexPosition::BOTTOM_BACK_EDGE_3D: // gh
653 triangleId = getVertexTriangleGH(p.data(), localTriangleId);
654 break;
655 case VertexPosition::LEFT_BACK_EDGE_3D: // eg
656 triangleId = getVertexTriangleEG(p.data(), localTriangleId);
657 break;
658 case VertexPosition::RIGHT_BACK_EDGE_3D: // fh
659 triangleId = getVertexTriangleFH(p.data(), localTriangleId);
660 break;
661 case VertexPosition::TOP_LEFT_EDGE_3D: // ae
662 triangleId = getVertexTriangleAE(p.data(), localTriangleId);
663 break;
664 case VertexPosition::TOP_RIGHT_EDGE_3D: // bf
665 triangleId = getVertexTriangleBF(p.data(), localTriangleId);
666 break;
667 case VertexPosition::BOTTOM_LEFT_EDGE_3D: // cg
668 triangleId = getVertexTriangleCG(p.data(), localTriangleId);
669 break;
670 case VertexPosition::BOTTOM_RIGHT_EDGE_3D: // dh
671 triangleId = getVertexTriangleDH(p.data(), localTriangleId);
672 break;
673 case VertexPosition::TOP_LEFT_FRONT_CORNER_3D: // a
674 triangleId = getVertexTriangleA(p.data(), localTriangleId);
675 break;
676 case VertexPosition::TOP_RIGHT_FRONT_CORNER_3D: // b
677 triangleId = getVertexTriangleB(p.data(), localTriangleId);
678 break;
679 case VertexPosition::BOTTOM_LEFT_FRONT_CORNER_3D: // c
680 triangleId = getVertexTriangleC(p.data(), localTriangleId);
681 break;
682 case VertexPosition::BOTTOM_RIGHT_FRONT_CORNER_3D: // d
683 triangleId = getVertexTriangleD(p.data(), localTriangleId);
684 break;
685 case VertexPosition::TOP_LEFT_BACK_CORNER_3D: // e
686 triangleId = getVertexTriangleE(p.data(), localTriangleId);
687 break;
688 case VertexPosition::TOP_RIGHT_BACK_CORNER_3D: // f
689 triangleId = getVertexTriangleF(p.data(), localTriangleId);
690 break;
691 case VertexPosition::BOTTOM_LEFT_BACK_CORNER_3D: // g
692 triangleId = getVertexTriangleG(p.data(), localTriangleId);
693 break;
694 case VertexPosition::BOTTOM_RIGHT_BACK_CORNER_3D: // h
695 triangleId = getVertexTriangleH(p.data(), localTriangleId);
696 break;
697 default: // 1D + 2D
698 triangleId = -1;
699 break;
700 }
701
702 return 0;
703}
704
705const vector<vector<SimplexId>> *
707 if(vertexTriangleList_.empty()) {
708 Timer t;
709
711 for(SimplexId i = 0; i < vertexNumber_; ++i) {
713 for(SimplexId j = 0; j < (SimplexId)vertexTriangleList_[i].size(); ++j)
715 }
716
717 printMsg("Built " + to_string(vertexNumber_) + " vertex triangles.", 1,
718 t.getElapsedTime(), 1);
719 }
720
721 return &vertexTriangleList_;
722}
723
724SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(
725 getVertexLinkNumber)(const SimplexId &vertexId) const {
727}
728
729template <typename Derived>
731 getVertexLink)(const SimplexId &vertexId,
732 const int &localLinkId,
733 SimplexId &linkId) const {
734
735#ifndef TTK_ENABLE_KAMIKAZE
736 if(localLinkId < 0 or localLinkId >= getVertexLinkNumber(vertexId))
737 return -1;
738#endif // !TTK_ENABLE_KAMIKAZE
739
740 const auto &p = this->underlying().getVertexCoords(vertexId);
741
742 switch(this->underlying().getVertexPosition(vertexId)) {
743 case VertexPosition::CENTER_3D:
744 linkId = getVertexLinkABCDEFGH(p.data(), localLinkId);
745 break;
746 case VertexPosition::FRONT_FACE_3D:
747 linkId = getVertexLinkABDC(p.data(), localLinkId);
748 break;
749 case VertexPosition::BACK_FACE_3D:
750 linkId = getVertexLinkEFHG(p.data(), localLinkId);
751 break;
752 case VertexPosition::TOP_FACE_3D:
753 linkId = getVertexLinkAEFB(p.data(), localLinkId);
754 break;
755 case VertexPosition::BOTTOM_FACE_3D:
756 linkId = getVertexLinkGHDC(p.data(), localLinkId);
757 break;
758 case VertexPosition::LEFT_FACE_3D:
759 linkId = getVertexLinkAEGC(p.data(), localLinkId);
760 break;
761 case VertexPosition::RIGHT_FACE_3D:
762 linkId = getVertexLinkBFHD(p.data(), localLinkId);
763 break;
764 case VertexPosition::TOP_FRONT_EDGE_3D: // ab
765 linkId = getVertexLinkAB(p.data(), localLinkId);
766 break;
767 case VertexPosition::BOTTOM_FRONT_EDGE_3D: // cd
768 linkId = getVertexLinkCD(p.data(), localLinkId);
769 break;
770 case VertexPosition::LEFT_FRONT_EDGE_3D: // ac
771 linkId = getVertexLinkAC(p.data(), localLinkId);
772 break;
773 case VertexPosition::RIGHT_FRONT_EDGE_3D: // bd
774 linkId = getVertexLinkBD(p.data(), localLinkId);
775 break;
776 case VertexPosition::TOP_BACK_EDGE_3D: // ef
777 linkId = getVertexLinkEF(p.data(), localLinkId);
778 break;
779 case VertexPosition::BOTTOM_BACK_EDGE_3D: // gh
780 linkId = getVertexLinkGH(p.data(), localLinkId);
781 break;
782 case VertexPosition::LEFT_BACK_EDGE_3D: // eg
783 linkId = getVertexLinkEG(p.data(), localLinkId);
784 break;
785 case VertexPosition::RIGHT_BACK_EDGE_3D: // fh
786 linkId = getVertexLinkFH(p.data(), localLinkId);
787 break;
788 case VertexPosition::TOP_LEFT_EDGE_3D: // ae
789 linkId = getVertexLinkAE(p.data(), localLinkId);
790 break;
791 case VertexPosition::TOP_RIGHT_EDGE_3D: // bf
792 linkId = getVertexLinkBF(p.data(), localLinkId);
793 break;
794 case VertexPosition::BOTTOM_LEFT_EDGE_3D: // cg
795 linkId = getVertexLinkCG(p.data(), localLinkId);
796 break;
797 case VertexPosition::BOTTOM_RIGHT_EDGE_3D: // dh
798 linkId = getVertexLinkDH(p.data(), localLinkId);
799 break;
800 case VertexPosition::TOP_LEFT_FRONT_CORNER_3D: // a
801 linkId = getVertexLinkA(p.data(), localLinkId);
802 break;
803 case VertexPosition::TOP_RIGHT_FRONT_CORNER_3D: // b
804 linkId = getVertexLinkB(p.data(), localLinkId);
805 break;
806 case VertexPosition::BOTTOM_LEFT_FRONT_CORNER_3D: // c
807 linkId = getVertexLinkC(p.data(), localLinkId);
808 break;
809 case VertexPosition::BOTTOM_RIGHT_FRONT_CORNER_3D: // d
810 linkId = getVertexLinkD(p.data(), localLinkId);
811 break;
812 case VertexPosition::TOP_LEFT_BACK_CORNER_3D: // e
813 linkId = getVertexLinkE(p.data(), localLinkId);
814 break;
815 case VertexPosition::TOP_RIGHT_BACK_CORNER_3D: // f
816 linkId = getVertexLinkF(p.data(), localLinkId);
817 break;
818 case VertexPosition::BOTTOM_LEFT_BACK_CORNER_3D: // g
819 linkId = getVertexLinkG(p.data(), localLinkId);
820 break;
821 case VertexPosition::BOTTOM_RIGHT_BACK_CORNER_3D: // h
822 linkId = getVertexLinkH(p.data(), localLinkId);
823 break;
824 case VertexPosition::CENTER_2D:
825 linkId = getVertexLink2dABCD(p.data(), localLinkId);
826 break;
827 case VertexPosition::TOP_EDGE_2D:
828 linkId = getVertexLink2dAB(p.data(), localLinkId);
829 break;
830 case VertexPosition::BOTTOM_EDGE_2D:
831 linkId = getVertexLink2dCD(p.data(), localLinkId);
832 break;
833 case VertexPosition::LEFT_EDGE_2D:
834 linkId = getVertexLink2dAC(p.data(), localLinkId);
835 break;
836 case VertexPosition::RIGHT_EDGE_2D:
837 linkId = getVertexLink2dBD(p.data(), localLinkId);
838 break;
839 case VertexPosition::TOP_LEFT_CORNER_2D: // a
840 linkId = getVertexLink2dA(p.data(), localLinkId);
841 break;
842 case VertexPosition::TOP_RIGHT_CORNER_2D: // b
843 linkId = getVertexLink2dB(p.data(), localLinkId);
844 break;
845 case VertexPosition::BOTTOM_LEFT_CORNER_2D: // c
846 linkId = getVertexLink2dC(p.data(), localLinkId);
847 break;
848 case VertexPosition::BOTTOM_RIGHT_CORNER_2D: // d
849 linkId = getVertexLink2dD(p.data(), localLinkId);
850 break;
851 default: // 1D
852 linkId = -1;
853 break;
854 };
855
856 return 0;
857}
858
859const vector<vector<SimplexId>> *
860 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getVertexLinks)() {
861 if(vertexLinkList_.empty()) {
862 Timer t;
863
864 vertexLinkList_.resize(vertexNumber_);
865 for(SimplexId i = 0; i < vertexNumber_; ++i) {
866 vertexLinkList_[i].resize(getVertexLinkNumber(i));
867 for(SimplexId j = 0; j < (SimplexId)vertexLinkList_[i].size(); ++j)
868 getVertexLink(i, j, vertexLinkList_[i][j]);
869 }
870
871 printMsg("Built " + to_string(vertexNumber_) + " vertex links.", 1,
872 t.getElapsedTime(), 1);
873 }
874
875 return &vertexLinkList_;
876}
877
878template <typename Derived>
880 getVertexStarNumber)(const SimplexId &vertexId) const {
881
882#ifndef TTK_ENABLE_KAMIKAZE
883 if(vertexId < 0 or vertexId >= vertexNumber_)
884 return -1;
885#endif // !TTK_ENABLE_KAMIKAZE
886
887 switch(this->underlying().getVertexPosition(vertexId)) {
888 case VertexPosition::CENTER_3D:
889 return 24;
890 case VertexPosition::FRONT_FACE_3D:
891 case VertexPosition::BACK_FACE_3D:
892 case VertexPosition::TOP_FACE_3D:
893 case VertexPosition::BOTTOM_FACE_3D:
894 case VertexPosition::LEFT_FACE_3D:
895 case VertexPosition::RIGHT_FACE_3D:
896 return 12;
897 case VertexPosition::TOP_FRONT_EDGE_3D: // ab
898 case VertexPosition::RIGHT_FRONT_EDGE_3D: // bd
899 case VertexPosition::BOTTOM_BACK_EDGE_3D: // gh
900 case VertexPosition::LEFT_BACK_EDGE_3D: // eg
901 case VertexPosition::BOTTOM_LEFT_EDGE_3D: // cg
902 case VertexPosition::TOP_RIGHT_EDGE_3D: // bf
903 return 8;
904 case VertexPosition::TOP_RIGHT_FRONT_CORNER_3D: // b
905 case VertexPosition::BOTTOM_LEFT_BACK_CORNER_3D: // g
906 case VertexPosition::CENTER_2D:
907 return 6;
908 case VertexPosition::TOP_BACK_EDGE_3D: // ef
909 case VertexPosition::BOTTOM_FRONT_EDGE_3D: // cd
910 case VertexPosition::LEFT_FRONT_EDGE_3D: // ac
911 case VertexPosition::TOP_LEFT_EDGE_3D: // ae
912 case VertexPosition::RIGHT_BACK_EDGE_3D: // fh
913 case VertexPosition::BOTTOM_RIGHT_EDGE_3D: // dh
914 return 4;
915 case VertexPosition::TOP_EDGE_2D: // ab
916 case VertexPosition::BOTTOM_EDGE_2D: // cd
917 case VertexPosition::LEFT_EDGE_2D: // ac
918 case VertexPosition::RIGHT_EDGE_2D: // bd
919 return 3;
920 case VertexPosition::TOP_LEFT_FRONT_CORNER_3D: // a
921 case VertexPosition::BOTTOM_LEFT_FRONT_CORNER_3D: // c
922 case VertexPosition::BOTTOM_RIGHT_FRONT_CORNER_3D: // d
923 case VertexPosition::TOP_LEFT_BACK_CORNER_3D: // e
924 case VertexPosition::TOP_RIGHT_BACK_CORNER_3D: // f
925 case VertexPosition::BOTTOM_RIGHT_BACK_CORNER_3D: // h
926 case VertexPosition::TOP_RIGHT_CORNER_2D: // b
927 case VertexPosition::BOTTOM_LEFT_CORNER_2D: // c
928 return 2;
929 case VertexPosition::TOP_LEFT_CORNER_2D: // a
930 case VertexPosition::BOTTOM_RIGHT_CORNER_2D: // d
931 return 1;
932 default: // 1D
933 break;
934 }
935
936 return 0;
937}
938
939template <typename Derived>
941 getVertexStar)(const SimplexId &vertexId,
942 const int &localStarId,
943 SimplexId &starId) const {
944
945#ifndef TTK_ENABLE_KAMIKAZE
946 if(localStarId < 0 or localStarId >= getVertexStarNumber(vertexId))
947 return -1;
948#endif // !TTK_ENABLE_KAMIKAZE
949
950 const auto &p = this->underlying().getVertexCoords(vertexId);
951
952 switch(this->underlying().getVertexPosition(vertexId)) {
953 case VertexPosition::CENTER_3D:
954 starId = getVertexStarABCDEFGH(p.data(), localStarId);
955 break;
956 case VertexPosition::FRONT_FACE_3D:
957 starId = getVertexStarABDC(p.data(), localStarId);
958 break;
959 case VertexPosition::BACK_FACE_3D:
960 starId = getVertexStarEFHG(p.data(), localStarId);
961 break;
962 case VertexPosition::TOP_FACE_3D:
963 starId = getVertexStarAEFB(p.data(), localStarId);
964 break;
965 case VertexPosition::BOTTOM_FACE_3D:
966 starId = getVertexStarGHDC(p.data(), localStarId);
967 break;
968 case VertexPosition::LEFT_FACE_3D:
969 starId = getVertexStarAEGC(p.data(), localStarId);
970 break;
971 case VertexPosition::RIGHT_FACE_3D:
972 starId = getVertexStarBFHD(p.data(), localStarId);
973 break;
974 case VertexPosition::TOP_FRONT_EDGE_3D: // ab
975 starId = getVertexStarAB(p.data(), localStarId);
976 break;
977 case VertexPosition::BOTTOM_FRONT_EDGE_3D: // cd
978 starId = getVertexStarCD(p.data(), localStarId);
979 break;
980 case VertexPosition::LEFT_FRONT_EDGE_3D: // ac
981 starId = getVertexStarAC(p.data(), localStarId);
982 break;
983 case VertexPosition::RIGHT_FRONT_EDGE_3D: // bd
984 starId = getVertexStarBD(p.data(), localStarId);
985 break;
986 case VertexPosition::TOP_BACK_EDGE_3D: // ef
987 starId = getVertexStarEF(p.data(), localStarId);
988 break;
989 case VertexPosition::BOTTOM_BACK_EDGE_3D: // gh
990 starId = getVertexStarGH(p.data(), localStarId);
991 break;
992 case VertexPosition::LEFT_BACK_EDGE_3D: // eg
993 starId = getVertexStarEG(p.data(), localStarId);
994 break;
995 case VertexPosition::RIGHT_BACK_EDGE_3D: // fh
996 starId = getVertexStarFH(p.data(), localStarId);
997 break;
998 case VertexPosition::TOP_LEFT_EDGE_3D: // ae
999 starId = getVertexStarAE(p.data(), localStarId);
1000 break;
1001 case VertexPosition::TOP_RIGHT_EDGE_3D: // bf
1002 starId = getVertexStarBF(p.data(), localStarId);
1003 break;
1004 case VertexPosition::BOTTOM_LEFT_EDGE_3D: // cg
1005 starId = getVertexStarCG(p.data(), localStarId);
1006 break;
1007 case VertexPosition::BOTTOM_RIGHT_EDGE_3D: // dh
1008 starId = getVertexStarDH(p.data(), localStarId);
1009 break;
1010 case VertexPosition::TOP_LEFT_FRONT_CORNER_3D: // a
1011 starId = getVertexStarA(p.data(), localStarId);
1012 break;
1013 case VertexPosition::TOP_RIGHT_FRONT_CORNER_3D: // b
1014 starId = getVertexStarB(p.data(), localStarId);
1015 break;
1016 case VertexPosition::BOTTOM_LEFT_FRONT_CORNER_3D: // c
1017 starId = getVertexStarC(p.data(), localStarId);
1018 break;
1019 case VertexPosition::BOTTOM_RIGHT_FRONT_CORNER_3D: // d
1020 starId = getVertexStarD(p.data(), localStarId);
1021 break;
1022 case VertexPosition::TOP_LEFT_BACK_CORNER_3D: // e
1023 starId = getVertexStarE(p.data(), localStarId);
1024 break;
1025 case VertexPosition::TOP_RIGHT_BACK_CORNER_3D: // f
1026 starId = getVertexStarF(p.data(), localStarId);
1027 break;
1028 case VertexPosition::BOTTOM_LEFT_BACK_CORNER_3D: // g
1029 starId = getVertexStarG(p.data(), localStarId);
1030 break;
1031 case VertexPosition::BOTTOM_RIGHT_BACK_CORNER_3D: // h
1032 starId = getVertexStarH(p.data(), localStarId);
1033 break;
1034 case VertexPosition::CENTER_2D:
1035 starId = getVertexStar2dABCD(p.data(), localStarId);
1036 break;
1037 case VertexPosition::TOP_EDGE_2D:
1038 starId = getVertexStar2dAB(p.data(), localStarId);
1039 break;
1040 case VertexPosition::BOTTOM_EDGE_2D:
1041 starId = getVertexStar2dCD(p.data(), localStarId);
1042 break;
1043 case VertexPosition::LEFT_EDGE_2D:
1044 starId = getVertexStar2dAC(p.data(), localStarId);
1045 break;
1046 case VertexPosition::RIGHT_EDGE_2D:
1047 starId = getVertexStar2dBD(p.data(), localStarId);
1048 break;
1049 case VertexPosition::TOP_LEFT_CORNER_2D: // a
1050 starId = getVertexStar2dA(p.data(), localStarId);
1051 break;
1052 case VertexPosition::TOP_RIGHT_CORNER_2D: // b
1053 starId = getVertexStar2dB(p.data(), localStarId);
1054 break;
1055 case VertexPosition::BOTTOM_LEFT_CORNER_2D: // c
1056 starId = getVertexStar2dC(p.data(), localStarId);
1057 break;
1058 case VertexPosition::BOTTOM_RIGHT_CORNER_2D: // d
1059 starId = getVertexStar2dD(p.data(), localStarId);
1060 break;
1061 default: // 1D
1062 starId = -1;
1063 break;
1065
1066 return 0;
1067}
1068
1069const vector<vector<SimplexId>> *
1070 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getVertexStars)() {
1071
1072 if(vertexStarList_.empty()) {
1073 Timer t;
1074 vertexStarList_.resize(vertexNumber_);
1075 for(SimplexId i = 0; i < vertexNumber_; ++i) {
1076 vertexStarList_[i].resize(getVertexStarNumber(i));
1077 for(SimplexId j = 0; j < (SimplexId)vertexStarList_[i].size(); ++j)
1078 getVertexStar(i, j, vertexStarList_[i][j]);
1079 }
1080
1081 printMsg("Built " + to_string(vertexNumber_) + " vertex stars.", 1,
1082 t.getElapsedTime(), 1);
1083 }
1084
1085 return &vertexStarList_;
1086}
1087
1088template <typename Derived>
1090 getVertexPoint)(const SimplexId &vertexId,
1091 float &x,
1092 float &y,
1093 float &z) const {
1094
1095 if(dimensionality_ == 3) {
1096 const auto &p = this->underlying().getVertexCoords(vertexId);
1097
1098 x = origin_[0] + spacing_[0] * p[0];
1099 y = origin_[1] + spacing_[1] * p[1];
1100 z = origin_[2] + spacing_[2] * p[2];
1101 } else if(dimensionality_ == 2) {
1102 const auto &p = this->underlying().getVertexCoords(vertexId);
1103
1104 if(dimensions_[0] > 1 and dimensions_[1] > 1) {
1105 x = origin_[0] + spacing_[0] * p[0];
1106 y = origin_[1] + spacing_[1] * p[1];
1107 z = origin_[2];
1108 } else if(dimensions_[1] > 1 and dimensions_[2] > 1) {
1109 x = origin_[0];
1110 y = origin_[1] + spacing_[1] * p[0];
1111 z = origin_[2] + spacing_[2] * p[1];
1112 } else if(dimensions_[0] > 1 and dimensions_[2] > 1) {
1113 x = origin_[0] + spacing_[0] * p[0];
1114 y = origin_[1];
1115 z = origin_[2] + spacing_[2] * p[1];
1117 } else if(dimensionality_ == 1) {
1118 if(dimensions_[0] > 1) {
1119 x = origin_[0] + spacing_[0] * vertexId;
1120 y = origin_[1];
1121 z = origin_[2];
1122 } else if(dimensions_[1] > 1) {
1123 x = origin_[0];
1124 y = origin_[1] + spacing_[1] * vertexId;
1125 z = origin_[2];
1126 } else if(dimensions_[2] > 1) {
1127 x = origin_[0];
1128 y = origin_[1];
1129 z = origin_[2] + spacing_[2] * vertexId;
1130 }
1131 }
1132
1133 return 0;
1134}
1135
1136template <typename Derived>
1138 const SimplexId &edgeId,
1139 const int &localVertexId,
1140 SimplexId &vertexId) const {
1141#ifndef TTK_ENABLE_KAMIKAZE
1142 if(edgeId < 0 or edgeId >= edgeNumber_)
1143 return -1;
1144 if(localVertexId < 0 or localVertexId >= 2)
1145 return -2;
1146#endif
1147
1148 const auto &p = this->underlying().getEdgeCoords(edgeId);
1149
1150 const auto helper3d = [&](const SimplexId a, const SimplexId b) -> SimplexId {
1151 if(isAccelerated_) {
1152 const auto tmp = p[0] + (p[1] << div_[0]) + (p[2] << div_[1]);
1153 return (localVertexId == 0) ? tmp + a : tmp + b;
1154 } else {
1155 const auto tmp = p[0] + (p[1] * vshift_[0]) + (p[2] * vshift_[1]);
1156 return (localVertexId == 0) ? tmp + a : tmp + b;
1157 }
1158 };
1159
1160 const auto helper2d = [&](const SimplexId a, const SimplexId b) -> SimplexId {
1161 if(isAccelerated_) {
1162 const auto tmp = p[0] + (p[1] << div_[0]);
1163 return localVertexId == 0 ? tmp + a : tmp + b;
1164 } else {
1165 const auto tmp = p[0] + (p[1] * vshift_[0]);
1166 return localVertexId == 0 ? tmp + a : tmp + b;
1167 }
1168 };
1169
1170 switch(this->underlying().getEdgePosition(edgeId)) {
1172 vertexId = helper3d(0, 1);
1173 break;
1175 vertexId = helper3d(0, vshift_[0]);
1176 break;
1178 vertexId = helper3d(0, vshift_[1]);
1179 break;
1181 vertexId = helper3d(1, vshift_[0]);
1182 break;
1184 vertexId = helper3d(0, vshift_[0] + vshift_[1]);
1185 break;
1187 vertexId = helper3d(1, vshift_[1]);
1188 break;
1189 case EdgePosition::D4_3D:
1190 vertexId = helper3d(1, vshift_[0] + vshift_[1]);
1191 break;
1192
1194 vertexId = helper2d(0, 1);
1195 break;
1197 vertexId = helper2d(0, vshift_[0]);
1198 break;
1199 case EdgePosition::D1_2D:
1200 vertexId = helper2d(1, vshift_[0]);
1201 break;
1202
1203 case EdgePosition::FIRST_EDGE_1D:
1204 vertexId = localVertexId == 0 ? 0 : 1;
1205 break;
1206 case EdgePosition::LAST_EDGE_1D:
1207 vertexId = localVertexId == 0 ? edgeNumber_ - 1 : edgeNumber_;
1208 break;
1209 case EdgePosition::CENTER_1D:
1210 vertexId = localVertexId == 0 ? edgeId : edgeId + 1;
1211 break;
1212 }
1213
1214 return 0;
1215}
1216
1217const vector<std::array<SimplexId, 2>> *
1218 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getEdges)() {
1219
1220 if(edgeList_.empty()) {
1221 Timer t;
1222
1223 edgeList_.resize(edgeNumber_);
1224 for(SimplexId i = 0; i < edgeNumber_; ++i) {
1225 SimplexId id0, id1;
1226 getEdgeVertexInternal(i, 0, id0);
1227 getEdgeVertexInternal(i, 1, id1);
1228 edgeList_[i] = {id0, id1};
1229 }
1230
1231 printMsg(
1232 "Built " + to_string(edgeNumber_) + " edges.", 1, t.getElapsedTime(), 1);
1233 }
1234
1235 return &edgeList_;
1236}
1237
1238template <typename Derived>
1240 const SimplexId &edgeId) const {
1241#ifndef TTK_ENABLE_KAMIKAZE
1242 if(edgeId < 0 or edgeId >= edgeNumber_)
1243 return -1;
1244#endif
1245
1246 switch(this->underlying().getEdgePosition(edgeId)) {
1247 case EdgePosition::L_xnn_3D:
1248 case EdgePosition::H_nyn_3D:
1249 case EdgePosition::P_nnz_3D:
1250 case EdgePosition::D4_3D:
1251 return 6;
1252 case EdgePosition::L_x0n_3D:
1253 case EdgePosition::L_xNn_3D:
1254 case EdgePosition::L_xn0_3D:
1255 case EdgePosition::L_xnN_3D:
1256 case EdgePosition::H_ny0_3D:
1257 case EdgePosition::H_nyN_3D:
1258 case EdgePosition::H_0yn_3D:
1259 case EdgePosition::H_Nyn_3D:
1260 case EdgePosition::P_n0z_3D:
1261 case EdgePosition::P_nNz_3D:
1262 case EdgePosition::P_0nz_3D:
1263 case EdgePosition::P_Nnz_3D:
1264 case EdgePosition::D1_xyn_3D:
1265 case EdgePosition::D2_nyz_3D:
1266 case EdgePosition::D3_xnz_3D:
1267 return 4;
1268 case EdgePosition::L_x00_3D:
1269 case EdgePosition::L_xNN_3D:
1270 case EdgePosition::H_0yN_3D:
1271 case EdgePosition::H_Ny0_3D:
1272 case EdgePosition::P_0Nz_3D:
1273 case EdgePosition::P_N0z_3D:
1274 case EdgePosition::D1_xy0_3D:
1275 case EdgePosition::D1_xyN_3D:
1276 case EdgePosition::D2_0yz_3D:
1277 case EdgePosition::D2_Nyz_3D:
1278 case EdgePosition::D3_x0z_3D:
1279 case EdgePosition::D3_xNz_3D:
1280 return 3;
1281 case EdgePosition::L_xN0_3D:
1282 case EdgePosition::L_x0N_3D:
1283 case EdgePosition::H_0y0_3D:
1284 case EdgePosition::H_NyN_3D:
1285 case EdgePosition::P_00z_3D:
1286 case EdgePosition::P_NNz_3D:
1287 case EdgePosition::L_xn_2D:
1288 case EdgePosition::H_ny_2D:
1289 case EdgePosition::D1_2D:
1290 return 2;
1291 case EdgePosition::L_x0_2D:
1292 case EdgePosition::L_xN_2D:
1293 case EdgePosition::H_0y_2D:
1294 case EdgePosition::H_Ny_2D:
1295 return 1;
1296
1297 default: // 1D
1298 break;
1299 }
1300
1301 return 0;
1302}
1303
1304template <typename Derived>
1306 const SimplexId &edgeId,
1307 const int &localTriangleId,
1308 SimplexId &triangleId) const {
1309#ifndef TTK_ENABLE_KAMIKAZE
1310 if(localTriangleId < 0
1311 or localTriangleId >= getEdgeTriangleNumberInternal(edgeId))
1312 return -1;
1313#endif
1314
1315 const auto &p = this->underlying().getEdgeCoords(edgeId);
1316
1317 switch(this->underlying().getEdgePosition(edgeId)) {
1318 case EdgePosition::L_xnn_3D:
1319 triangleId = getEdgeTriangleL_xnn(p.data(), localTriangleId);
1320 break;
1321 case EdgePosition::L_xn0_3D:
1322 triangleId = getEdgeTriangleL_xn0(p.data(), localTriangleId);
1323 break;
1324 case EdgePosition::L_xnN_3D:
1325 triangleId = getEdgeTriangleL_xnN(p.data(), localTriangleId);
1326 break;
1327 case EdgePosition::L_x0n_3D:
1328 triangleId = getEdgeTriangleL_x0n(p.data(), localTriangleId);
1329 break;
1330 case EdgePosition::L_x00_3D:
1331 triangleId = getEdgeTriangleL_x00(p.data(), localTriangleId);
1332 break;
1333 case EdgePosition::L_x0N_3D:
1334 triangleId = getEdgeTriangleL_x0N(p.data(), localTriangleId);
1335 break;
1336 case EdgePosition::L_xNn_3D:
1337 triangleId = getEdgeTriangleL_xNn(p.data(), localTriangleId);
1338 break;
1339 case EdgePosition::L_xN0_3D:
1340 triangleId = getEdgeTriangleL_xN0(p.data(), localTriangleId);
1341 break;
1342 case EdgePosition::L_xNN_3D:
1343 triangleId = getEdgeTriangleL_xNN(p.data(), localTriangleId);
1344 break;
1345 case EdgePosition::H_nyn_3D:
1346 triangleId = getEdgeTriangleH_nyn(p.data(), localTriangleId);
1347 break;
1348 case EdgePosition::H_ny0_3D:
1349 triangleId = getEdgeTriangleH_ny0(p.data(), localTriangleId);
1350 break;
1351 case EdgePosition::H_nyN_3D:
1352 triangleId = getEdgeTriangleH_nyN(p.data(), localTriangleId);
1353 break;
1354 case EdgePosition::H_0yn_3D:
1355 triangleId = getEdgeTriangleH_0yn(p.data(), localTriangleId);
1356 break;
1357 case EdgePosition::H_0y0_3D:
1358 triangleId = getEdgeTriangleH_0y0(p.data(), localTriangleId);
1359 break;
1360 case EdgePosition::H_0yN_3D:
1361 triangleId = getEdgeTriangleH_0yN(p.data(), localTriangleId);
1362 break;
1363 case EdgePosition::H_Nyn_3D:
1364 triangleId = getEdgeTriangleH_Nyn(p.data(), localTriangleId);
1365 break;
1366 case EdgePosition::H_Ny0_3D:
1367 triangleId = getEdgeTriangleH_Ny0(p.data(), localTriangleId);
1368 break;
1369 case EdgePosition::H_NyN_3D:
1370 triangleId = getEdgeTriangleH_NyN(p.data(), localTriangleId);
1371 break;
1372 case EdgePosition::P_nnz_3D:
1373 triangleId = getEdgeTriangleP_nnz(p.data(), localTriangleId);
1374 break;
1375 case EdgePosition::P_n0z_3D:
1376 triangleId = getEdgeTriangleP_n0z(p.data(), localTriangleId);
1377 break;
1378 case EdgePosition::P_nNz_3D:
1379 triangleId = getEdgeTriangleP_nNz(p.data(), localTriangleId);
1380 break;
1381 case EdgePosition::P_0nz_3D:
1382 triangleId = getEdgeTriangleP_0nz(p.data(), localTriangleId);
1383 break;
1384 case EdgePosition::P_00z_3D:
1385 triangleId = getEdgeTriangleP_00z(p.data(), localTriangleId);
1386 break;
1387 case EdgePosition::P_0Nz_3D:
1388 triangleId = getEdgeTriangleP_0Nz(p.data(), localTriangleId);
1389 break;
1390 case EdgePosition::P_Nnz_3D:
1391 triangleId = getEdgeTriangleP_Nnz(p.data(), localTriangleId);
1392 break;
1393 case EdgePosition::P_N0z_3D:
1394 triangleId = getEdgeTriangleP_N0z(p.data(), localTriangleId);
1395 break;
1396 case EdgePosition::P_NNz_3D:
1397 triangleId = getEdgeTriangleP_NNz(p.data(), localTriangleId);
1398 break;
1399 case EdgePosition::D1_xyn_3D:
1400 triangleId = getEdgeTriangleD1_xyn(p.data(), localTriangleId);
1401 break;
1402 case EdgePosition::D1_xy0_3D:
1403 triangleId = getEdgeTriangleD1_xy0(p.data(), localTriangleId);
1404 break;
1405 case EdgePosition::D1_xyN_3D:
1406 triangleId = getEdgeTriangleD1_xyN(p.data(), localTriangleId);
1407 break;
1408 case EdgePosition::D2_nyz_3D:
1409 triangleId = getEdgeTriangleD2_nyz(p.data(), localTriangleId);
1410 break;
1411 case EdgePosition::D2_0yz_3D:
1412 triangleId = getEdgeTriangleD2_0yz(p.data(), localTriangleId);
1413 break;
1414 case EdgePosition::D2_Nyz_3D:
1415 triangleId = getEdgeTriangleD2_Nyz(p.data(), localTriangleId);
1416 break;
1417 case EdgePosition::D3_xnz_3D:
1418 triangleId = getEdgeTriangleD3_xnz(p.data(), localTriangleId);
1419 break;
1420 case EdgePosition::D3_x0z_3D:
1421 triangleId = getEdgeTriangleD3_x0z(p.data(), localTriangleId);
1422 break;
1423 case EdgePosition::D3_xNz_3D:
1424 triangleId = getEdgeTriangleD3_xNz(p.data(), localTriangleId);
1425 break;
1426 case EdgePosition::D4_3D:
1427 triangleId = getEdgeTriangleD4_xyz(p.data(), localTriangleId);
1428 break;
1429
1430 case EdgePosition::L_xn_2D:
1431 triangleId = getEdgeTriangleL_xn(p.data(), localTriangleId);
1432 break;
1433 case EdgePosition::L_x0_2D:
1434 triangleId = getEdgeTriangleL_x0(p.data(), localTriangleId);
1435 break;
1436 case EdgePosition::L_xN_2D:
1437 triangleId = getEdgeTriangleL_xN(p.data(), localTriangleId);
1438 break;
1439 case EdgePosition::H_ny_2D:
1440 triangleId = getEdgeTriangleH_ny(p.data(), localTriangleId);
1441 break;
1442 case EdgePosition::H_0y_2D:
1443 triangleId = getEdgeTriangleH_0y(p.data(), localTriangleId);
1444 break;
1445 case EdgePosition::H_Ny_2D:
1446 triangleId = getEdgeTriangleH_Ny(p.data(), localTriangleId);
1447 break;
1448 case EdgePosition::D1_2D:
1449 triangleId = getEdgeTriangleD1_xy(p.data(), localTriangleId);
1450 break;
1451
1452 default: // 1D
1453 triangleId = -1;
1454 break;
1455 }
1456
1457 return 0;
1458}
1459
1460const vector<vector<SimplexId>> *
1462 if(edgeTriangleList_.empty()) {
1463 Timer t;
1464
1466 for(SimplexId i = 0; i < edgeNumber_; ++i) {
1468 for(SimplexId j = 0; j < (SimplexId)edgeTriangleList_[i].size(); ++j)
1470 }
1471
1472 printMsg("Built " + to_string(edgeNumber_) + " edge triangles.", 1,
1473 t.getElapsedTime(), 1);
1474 }
1475
1476 return &edgeTriangleList_;
1477}
1478
1479SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getEdgeLinkNumber)(
1480 const SimplexId &edgeId) const {
1482}
1483
1484template <typename Derived>
1486 const SimplexId &edgeId, const int &localLinkId, SimplexId &linkId) const {
1487
1488#ifndef TTK_ENABLE_KAMIKAZE
1489 if(localLinkId < 0 or localLinkId >= getEdgeLinkNumber(edgeId))
1490 return -1;
1491#endif
1492
1493 const auto &p = this->underlying().getEdgeCoords(edgeId);
1494
1495 switch(this->underlying().getEdgePosition(edgeId)) {
1497 linkId = getEdgeLinkL(p.data(), localLinkId);
1498 break;
1500 linkId = getEdgeLinkH(p.data(), localLinkId);
1501 break;
1503 linkId = getEdgeLinkP(p.data(), localLinkId);
1504 break;
1506 linkId = getEdgeLinkD1(p.data(), localLinkId);
1507 break;
1509 linkId = getEdgeLinkD2(p.data(), localLinkId);
1510 break;
1512 linkId = getEdgeLinkD3(p.data(), localLinkId);
1513 break;
1514 case EdgePosition::D4_3D:
1515 linkId = getEdgeLinkD4(p.data(), localLinkId);
1516 break;
1517
1519 linkId = getEdgeLink2dL(p.data(), localLinkId);
1520 break;
1522 linkId = getEdgeLink2dH(p.data(), localLinkId);
1523 break;
1524 case EdgePosition::D1_2D:
1525 linkId = getEdgeLink2dD1(p.data(), localLinkId);
1526 break;
1527
1528 default: // 1D
1529 linkId = -1;
1530 break;
1531 }
1532
1533 return 0;
1534}
1535
1536const vector<vector<SimplexId>> *
1537 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getEdgeLinks)() {
1538
1539 if(edgeLinkList_.empty()) {
1540 Timer t;
1541
1542 edgeLinkList_.resize(edgeNumber_);
1543 for(SimplexId i = 0; i < edgeNumber_; ++i) {
1544 edgeLinkList_[i].resize(getEdgeLinkNumber(i));
1545 for(SimplexId j = 0; j < (SimplexId)edgeLinkList_[i].size(); ++j)
1546 getEdgeLink(i, j, edgeLinkList_[i][j]);
1547 }
1548
1549 printMsg("Built " + to_string(edgeNumber_) + " edge links.", 1,
1550 t.getElapsedTime(), 1);
1551 }
1552
1553 return &edgeLinkList_;
1554}
1555
1556template <typename Derived>
1558 getEdgeStarNumber)(const SimplexId &edgeId) const {
1559
1560#ifndef TTK_ENABLE_KAMIKAZE
1561 if(edgeId < 0 or edgeId >= edgeNumber_)
1562 return -1;
1563#endif
1564
1565 switch(this->underlying().getEdgePosition(edgeId)) {
1566 case EdgePosition::L_xnn_3D:
1567 case EdgePosition::H_nyn_3D:
1568 case EdgePosition::P_nnz_3D:
1569 case EdgePosition::D4_3D:
1570 return 6;
1571 case EdgePosition::D1_xyn_3D:
1572 case EdgePosition::D2_nyz_3D:
1573 case EdgePosition::D3_xnz_3D:
1574 return 4;
1575 case EdgePosition::L_x0n_3D:
1576 case EdgePosition::L_xNn_3D:
1577 case EdgePosition::L_xn0_3D:
1578 case EdgePosition::L_xnN_3D:
1579 case EdgePosition::H_ny0_3D:
1580 case EdgePosition::H_nyN_3D:
1581 case EdgePosition::H_0yn_3D:
1582 case EdgePosition::H_Nyn_3D:
1583 case EdgePosition::P_n0z_3D:
1584 case EdgePosition::P_nNz_3D:
1585 case EdgePosition::P_0nz_3D:
1586 case EdgePosition::P_Nnz_3D:
1587 return 3;
1588 case EdgePosition::L_x00_3D:
1589 case EdgePosition::L_xNN_3D:
1590 case EdgePosition::H_0yN_3D:
1591 case EdgePosition::H_Ny0_3D:
1592 case EdgePosition::P_0Nz_3D:
1593 case EdgePosition::P_N0z_3D:
1594 case EdgePosition::D1_xy0_3D:
1595 case EdgePosition::D1_xyN_3D:
1596 case EdgePosition::D2_0yz_3D:
1597 case EdgePosition::D2_Nyz_3D:
1598 case EdgePosition::D3_x0z_3D:
1599 case EdgePosition::D3_xNz_3D:
1600 case EdgePosition::L_xn_2D:
1601 case EdgePosition::H_ny_2D:
1602 case EdgePosition::D1_2D:
1603 return 2;
1604 case EdgePosition::L_xN0_3D:
1605 case EdgePosition::L_x0N_3D:
1606 case EdgePosition::H_0y0_3D:
1607 case EdgePosition::H_NyN_3D:
1608 case EdgePosition::P_00z_3D:
1609 case EdgePosition::P_NNz_3D:
1610 case EdgePosition::L_x0_2D:
1611 case EdgePosition::L_xN_2D:
1612 case EdgePosition::H_0y_2D:
1613 case EdgePosition::H_Ny_2D:
1614 return 1;
1615
1616 default: // 1D
1617 break;
1618 }
1619
1620 return 0;
1621}
1622
1623template <typename Derived>
1625 const SimplexId &edgeId, const int &localStarId, SimplexId &starId) const {
1626
1627#ifndef TTK_ENABLE_KAMIKAZE
1628 if(localStarId < 0 or localStarId >= getEdgeStarNumber(edgeId))
1629 return -1;
1630#endif
1631
1632 const auto &p = this->underlying().getEdgeCoords(edgeId);
1633
1634 switch(this->underlying().getEdgePosition(edgeId)) {
1636 starId = getEdgeStarL(p.data(), localStarId);
1637 break;
1639 starId = getEdgeStarH(p.data(), localStarId);
1640 break;
1642 starId = getEdgeStarP(p.data(), localStarId);
1643 break;
1645 starId = getEdgeStarD1(p.data(), localStarId);
1646 break;
1648 starId = getEdgeStarD2(p.data(), localStarId);
1649 break;
1651 starId = getEdgeStarD3(p.data(), localStarId);
1652 break;
1653 case EdgePosition::D4_3D:
1654 starId
1655 = p[2] * tetshift_[1] + p[1] * tetshift_[0] + p[0] * 6 + localStarId;
1656 break;
1657
1659 starId = getEdgeStar2dL(p.data(), localStarId);
1660 break;
1662 starId = getEdgeStar2dH(p.data(), localStarId);
1663 break;
1664 case EdgePosition::D1_2D:
1665 starId = p[0] * 2 + p[1] * tshift_[0] + localStarId;
1666 break;
1667
1668 default: // 1D
1669 starId = -1;
1670 break;
1671 }
1672
1673 return 0;
1674}
1675
1676const vector<vector<SimplexId>> *
1677 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getEdgeStars)() {
1678
1679 if(edgeStarList_.empty()) {
1680 Timer t;
1681
1682 edgeStarList_.resize(edgeNumber_);
1683 for(SimplexId i = 0; i < edgeNumber_; ++i) {
1684 edgeStarList_[i].resize(getEdgeStarNumber(i));
1685 for(SimplexId j = 0; j < (SimplexId)edgeStarList_[i].size(); ++j)
1686 getEdgeStar(i, j, edgeStarList_[i][j]);
1687 }
1688
1689 printMsg("Built " + to_string(edgeNumber_) + " edge stars.", 1,
1690 t.getElapsedTime(), 1);
1691 }
1692
1693 return &edgeStarList_;
1694}
1695
1696template <typename Derived>
1698 const SimplexId &triangleId,
1699 const int &localVertexId,
1700 SimplexId &vertexId) const {
1701#ifndef TTK_ENABLE_KAMIKAZE
1702 if(triangleId < 0 or triangleId >= triangleNumber_)
1703 return -1;
1704 if(localVertexId < 0 or localVertexId >= 3)
1705 return -2;
1706#endif
1707
1708 // e--------f
1709 // /| /|
1710 // / | / |
1711 // a--g-----b--h
1712 // | / | /
1713 // |/ |/
1714 // c--------d
1715 //
1716 // Classement des "Triangles" et dans cet ordre:
1717 // F: face (type abc/bcd)
1718 // C: cote (type abe/bef)
1719 // H: haut (type acg/aeg)
1720 // D1: diagonale1 (type bdg/beg)
1721 // D2: diagonale2 (type abg/bgh)
1722 // D3: diagonale3 (type bcg/bfg)
1723
1724 const auto &p = this->underlying().getTriangleCoords(triangleId);
1725 vertexId = -1;
1726
1727 switch(this->underlying().getTrianglePosition(triangleId)) {
1728 case TrianglePosition::F_3D:
1729 vertexId = getTriangleVertexF(p.data(), localVertexId);
1730 break;
1731 case TrianglePosition::H_3D:
1732 vertexId = getTriangleVertexH(p.data(), localVertexId);
1733 break;
1734 case TrianglePosition::C_3D:
1735 vertexId = getTriangleVertexC(p.data(), localVertexId);
1736 break;
1737 case TrianglePosition::D1_3D:
1738 vertexId = getTriangleVertexD1(p.data(), localVertexId);
1739 break;
1740 case TrianglePosition::D2_3D:
1741 vertexId = getTriangleVertexD2(p.data(), localVertexId);
1742 break;
1743 case TrianglePosition::D3_3D:
1744 vertexId = getTriangleVertexD3(p.data(), localVertexId);
1745 break;
1746 case TrianglePosition::TOP_2D:
1747 switch(localVertexId) {
1748 break;
1749 case 0:
1750 vertexId = p[0] / 2 + p[1] * vshift_[0];
1751 break;
1752 case 1:
1753 vertexId = p[0] / 2 + p[1] * vshift_[0] + 1;
1754 break;
1755 case 2:
1756 vertexId = p[0] / 2 + p[1] * vshift_[0] + vshift_[0];
1757 break;
1758 }
1759 break;
1760 case TrianglePosition::BOTTOM_2D:
1761 switch(localVertexId) {
1762 break;
1763 case 0:
1764 vertexId = p[0] / 2 + p[1] * vshift_[0] + 1;
1765 break;
1766 case 1:
1767 vertexId = p[0] / 2 + p[1] * vshift_[0] + vshift_[0] + 1;
1768 break;
1769 case 2:
1770 vertexId = p[0] / 2 + p[1] * vshift_[0] + vshift_[0];
1771 break;
1772 }
1773 }
1774
1775 return 0;
1776}
1777
1778template <typename Derived>
1780 const SimplexId &triangleId,
1781 const int &localEdgeId,
1782 SimplexId &edgeId) const {
1783#ifndef TTK_ENABLE_KAMIKAZE
1784 if(triangleId < 0 or triangleId >= triangleNumber_)
1785 return -1;
1786 if(localEdgeId < 0 or localEdgeId >= 3)
1787 return -2;
1788#endif
1789
1790 const auto &p = this->underlying().getTriangleCoords(triangleId);
1791 const auto par = triangleId % 2;
1792 edgeId = -1;
1793
1794 switch(this->underlying().getTrianglePosition(triangleId)) {
1795 case TrianglePosition::F_3D:
1796 edgeId = (par == 1) ? getTriangleEdgeF_1(p.data(), localEdgeId)
1797 : getTriangleEdgeF_0(p.data(), localEdgeId);
1798 break;
1799 case TrianglePosition::H_3D:
1800 edgeId = (par == 1) ? getTriangleEdgeH_1(p.data(), localEdgeId)
1801 : getTriangleEdgeH_0(p.data(), localEdgeId);
1802 break;
1803 case TrianglePosition::C_3D:
1804 edgeId = (par == 1) ? getTriangleEdgeC_1(p.data(), localEdgeId)
1805 : getTriangleEdgeC_0(p.data(), localEdgeId);
1806 break;
1807 case TrianglePosition::D1_3D:
1808 edgeId = (par == 1) ? getTriangleEdgeD1_1(p.data(), localEdgeId)
1809 : getTriangleEdgeD1_0(p.data(), localEdgeId);
1810 break;
1811 case TrianglePosition::D2_3D:
1812 edgeId = (par == 1) ? getTriangleEdgeD2_1(p.data(), localEdgeId)
1813 : getTriangleEdgeD2_0(p.data(), localEdgeId);
1814 break;
1815 case TrianglePosition::D3_3D:
1816 edgeId = (par == 1) ? getTriangleEdgeD3_1(p.data(), localEdgeId)
1817 : getTriangleEdgeD3_0(p.data(), localEdgeId);
1818 break;
1819 case TrianglePosition::TOP_2D:
1820 switch(localEdgeId) {
1821 break;
1822 case 0:
1823 edgeId = p[0] / 2 + p[1] * eshift_[0];
1824 break;
1825 case 1:
1826 edgeId = esetshift_[0] + p[0] / 2 + p[1] * eshift_[2];
1827 break;
1828 case 2:
1829 edgeId = esetshift_[1] + p[0] / 2 + p[1] * eshift_[4];
1830 break;
1831 }
1832 break;
1833 case TrianglePosition::BOTTOM_2D:
1834 switch(localEdgeId) {
1835 break;
1836 case 0:
1837 edgeId = p[0] / 2 + (p[1] + 1) * eshift_[0];
1838 break;
1839 case 1:
1840 edgeId = esetshift_[0] + (p[0] + 1) / 2 + p[1] * eshift_[2];
1841 break;
1842 case 2:
1843 edgeId = esetshift_[1] + p[0] / 2 + p[1] * eshift_[4];
1844 break;
1845 }
1846 }
1847
1848 return 0;
1849}
1850
1852 vector<vector<SimplexId>> &edges) const {
1853 edges.resize(triangleNumber_);
1854 for(SimplexId i = 0; i < triangleNumber_; ++i) {
1855 edges[i].resize(3);
1856 for(int j = 0; j < 3; ++j)
1857 getTriangleEdgeInternal(i, j, edges[i][j]);
1858 }
1859 return 0;
1860}
1861
1862const vector<vector<SimplexId>> *
1864 if(triangleEdgeVector_.empty()) {
1865 Timer t;
1866
1868
1869 printMsg("Built " + to_string(triangleNumber_) + " triangle edges.", 1,
1870 t.getElapsedTime(), 1);
1871 }
1872
1873 return &triangleEdgeVector_;
1874}
1875
1876const vector<std::array<SimplexId, 3>> *
1877 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getTriangles)() {
1878
1879 if(triangleList_.empty()) {
1880 Timer t;
1881
1882 triangleList_.resize(triangleNumber_);
1883 for(SimplexId i = 0; i < triangleNumber_; ++i) {
1884 for(int j = 0; j < 3; ++j)
1885 getTriangleVertexInternal(i, j, triangleList_[i][j]);
1886 }
1887
1888 printMsg("Built " + to_string(triangleNumber_) + " triangles.", 1,
1889 t.getElapsedTime(), 1);
1890 }
1891
1892 return &triangleList_;
1893}
1894
1895template <typename Derived>
1897 getTriangleLink)(const SimplexId &triangleId,
1898 const int &localLinkId,
1899 SimplexId &linkId) const {
1900
1901#ifndef TTK_ENABLE_KAMIKAZE
1902 if(localLinkId < 0 or localLinkId >= getTriangleLinkNumber(triangleId))
1903 return -1;
1904#endif
1905
1906 const auto &p = this->underlying().getTriangleCoords(triangleId);
1907
1908 switch(this->underlying().getTrianglePosition(triangleId)) {
1909 case TrianglePosition::F_3D:
1910 linkId = getTriangleLinkF(p.data(), localLinkId);
1911 break;
1912 case TrianglePosition::H_3D:
1913 linkId = getTriangleLinkH(p.data(), localLinkId);
1914 break;
1915 case TrianglePosition::C_3D:
1916 linkId = getTriangleLinkC(p.data(), localLinkId);
1917 break;
1918 case TrianglePosition::D1_3D:
1919 linkId = getTriangleLinkD1(p.data(), localLinkId);
1920 break;
1921 case TrianglePosition::D2_3D:
1922 linkId = getTriangleLinkD2(p.data(), localLinkId);
1923 break;
1924 case TrianglePosition::D3_3D:
1925 linkId = getTriangleLinkD3(p.data(), localLinkId);
1926 break;
1927 default: // 2D
1928 linkId = -1;
1929 break;
1930 }
1931
1932 return 0;
1933}
1934
1935SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(
1936 getTriangleLinkNumber)(const SimplexId &triangleId) const {
1938}
1939
1940const vector<vector<SimplexId>> *
1941 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getTriangleLinks)() {
1942 if(triangleLinkList_.empty()) {
1943 Timer t;
1944
1945 triangleLinkList_.resize(triangleNumber_);
1946 for(SimplexId i = 0; i < triangleNumber_; ++i) {
1947 triangleLinkList_[i].resize(getTriangleLinkNumber(i));
1948 for(SimplexId j = 0; j < (SimplexId)triangleLinkList_[i].size(); ++j)
1949 getTriangleLink(i, j, triangleLinkList_[i][j]);
1950 }
1951
1952 printMsg("Built " + to_string(triangleNumber_) + " triangle links.", 1,
1953 t.getElapsedTime(), 1);
1954 }
1955 return &triangleLinkList_;
1956}
1957
1958template <typename Derived>
1960 getTriangleStarNumber)(const SimplexId &triangleId) const {
1961
1962#ifndef TTK_ENABLE_KAMIKAZE
1963 if(triangleId < 0 or triangleId >= triangleNumber_)
1964 return -1;
1965#endif
1966
1967 const auto &p = this->underlying().getTriangleCoords(triangleId);
1968
1969 switch(this->underlying().getTrianglePosition(triangleId)) {
1970 case TrianglePosition::F_3D:
1971 return (p[2] > 0 and p[2] < nbvoxels_[2]) ? 2 : 1;
1972 case TrianglePosition::H_3D:
1973 return (p[1] > 0 and p[1] < nbvoxels_[1]) ? 2 : 1;
1974 case TrianglePosition::C_3D:
1975 return (p[0] < 2 or p[0] >= (dimensions_[0] * 2 - 2)) ? 1 : 2;
1976
1977 case TrianglePosition::D1_3D:
1978 case TrianglePosition::D2_3D:
1979 case TrianglePosition::D3_3D:
1980 return 2;
1981 default: // 2D
1982 break;
1983 }
1984 return 0;
1985}
1986
1987template <typename Derived>
1989 getTriangleStar)(const SimplexId &triangleId,
1990 const int &localStarId,
1991 SimplexId &starId) const {
1992
1993#ifndef TTK_ENABLE_KAMIKAZE
1994 if(localStarId < 0 or localStarId >= getTriangleStarNumber(triangleId))
1995 return -1;
1996#endif
1997
1998 const auto &p = this->underlying().getTriangleCoords(triangleId);
1999
2000 switch(this->underlying().getTrianglePosition(triangleId)) {
2001 case TrianglePosition::F_3D:
2002 starId = getTriangleStarF(p.data(), localStarId);
2003 break;
2004 case TrianglePosition::H_3D:
2005 starId = getTriangleStarH(p.data(), localStarId);
2006 break;
2007 case TrianglePosition::C_3D:
2008 starId = getTriangleStarC(p.data(), localStarId);
2009 break;
2010 case TrianglePosition::D1_3D:
2011 starId = getTriangleStarD1(p.data(), localStarId);
2012 break;
2013 case TrianglePosition::D2_3D:
2014 starId = getTriangleStarD2(p.data(), localStarId);
2015 break;
2016 case TrianglePosition::D3_3D:
2017 starId = getTriangleStarD3(p.data(), localStarId);
2018 break;
2019 default: // 2D
2020 starId = -1;
2021 break;
2022 }
2023
2024 return 0;
2025}
2026
2027const vector<vector<SimplexId>> *
2028 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getTriangleStars)() {
2029
2030 if(triangleStarList_.empty()) {
2031 Timer t;
2032
2033 triangleStarList_.resize(triangleNumber_);
2034 for(SimplexId i = 0; i < triangleNumber_; ++i) {
2035 triangleStarList_[i].resize(getTriangleStarNumber(i));
2036 for(SimplexId j = 0; j < (SimplexId)triangleStarList_[i].size(); ++j)
2037 getTriangleStar(i, j, triangleStarList_[i][j]);
2038 }
2039
2040 printMsg("Built " + to_string(triangleNumber_) + " triangle stars.", 1,
2041 t.getElapsedTime(), 1);
2042 }
2043 return &triangleStarList_;
2044}
2045
2046template <typename Derived>
2048 const SimplexId &triangleId) const {
2049#ifndef TTK_ENABLE_KAMIKAZE
2050 if(triangleId < 0 or triangleId >= triangleNumber_)
2051 return -1;
2052#endif
2053
2054 if(dimensionality_ == 2) {
2055 const auto &p = this->underlying().getTriangleCoords(triangleId);
2056 const SimplexId id = triangleId % 2;
2057
2058 if(id) {
2059 if(p[0] / 2 == nbvoxels_[Di_] - 1 and p[1] == nbvoxels_[Dj_] - 1)
2060 return 1;
2061 else if(p[0] / 2 == nbvoxels_[Di_] - 1 or p[1] == nbvoxels_[Dj_] - 1)
2062 return 2;
2063 else
2064 return 3;
2065 } else {
2066 if(p[0] == 0 and p[1] == 0)
2067 return 1;
2068 else if(p[0] == 0 or p[1] == 0)
2069 return 2;
2070 else
2071 return 3;
2072 }
2073 }
2074
2075 return 0;
2076}
2077
2078template <typename Derived>
2080 const SimplexId &triangleId,
2081 const int &localNeighborId,
2082 SimplexId &neighborId) const {
2083#ifndef TTK_ENABLE_KAMIKAZE
2084 if(localNeighborId < 0
2085 or localNeighborId >= getTriangleNeighborNumber(triangleId))
2086 return -1;
2087#endif
2088
2089 neighborId = -1;
2090
2091 if(dimensionality_ == 2) {
2092 const auto &p = this->underlying().getTriangleCoords(triangleId);
2093 const SimplexId id = triangleId % 2;
2094
2095 if(id) {
2096 if(p[0] / 2 == nbvoxels_[Di_] - 1 and p[1] == nbvoxels_[Dj_] - 1)
2097 neighborId = triangleId - 1;
2098 else if(p[0] / 2 == nbvoxels_[Di_] - 1) {
2099 switch(localNeighborId) {
2100 case 0:
2101 neighborId = triangleId - 1;
2102 break;
2103 case 1:
2104 neighborId = triangleId + tshift_[0] - 1;
2105 break;
2106 }
2107 } else if(p[1] == nbvoxels_[Dj_] - 1) {
2108 switch(localNeighborId) {
2109 case 0:
2110 neighborId = triangleId - 1;
2111 break;
2112 case 1:
2113 neighborId = triangleId + 1;
2114 break;
2115 }
2116 } else {
2117 switch(localNeighborId) {
2118 case 0:
2119 neighborId = triangleId - 1;
2120 break;
2121 case 1:
2122 neighborId = triangleId + 1;
2123 break;
2124 case 2:
2125 neighborId = triangleId + tshift_[0] - 1;
2126 break;
2127 }
2128 }
2129 } else {
2130 if(p[0] == 0 and p[1] == 0)
2131 neighborId = triangleId + 1;
2132 else if(p[0] == 0) {
2133 switch(localNeighborId) {
2134 case 0:
2135 neighborId = triangleId + 1;
2136 break;
2137 case 1:
2138 neighborId = triangleId - tshift_[0] + 1;
2139 break;
2140 }
2141 } else if(p[1] == 0) {
2142 switch(localNeighborId) {
2143 case 0:
2144 neighborId = triangleId + 1;
2145 break;
2146 case 1:
2147 neighborId = triangleId - 1;
2148 break;
2149 }
2150 } else {
2151 switch(localNeighborId) {
2152 case 0:
2153 neighborId = triangleId + 1;
2154 break;
2155 case 1:
2156 neighborId = triangleId - 1;
2157 break;
2158 case 2:
2159 neighborId = triangleId - tshift_[0] + 1;
2160 break;
2161 }
2162 }
2163 }
2164 }
2165
2166 return 0;
2167}
2168
2170 vector<vector<SimplexId>> &neighbors) {
2171 neighbors.resize(triangleNumber_);
2172 for(SimplexId i = 0; i < triangleNumber_; ++i) {
2173 neighbors[i].resize(getTriangleNeighborNumber(i));
2174 for(SimplexId j = 0; j < (SimplexId)neighbors[i].size(); ++j)
2175 getTriangleNeighbor(i, j, neighbors[i][j]);
2176 }
2177 return 0;
2178}
2179
2180template <typename Derived>
2182 const SimplexId &tetId, const int &localVertexId, SimplexId &vertexId) const {
2183#ifndef TTK_ENABLE_KAMIKAZE
2184 if(tetId < 0 or tetId >= tetrahedronNumber_)
2185 return -1;
2186 if(localVertexId < 0 or localVertexId >= 4)
2187 return -2;
2188#endif
2189
2190 vertexId = -1;
2191
2192 if(dimensionality_ == 3) {
2193 const SimplexId id = tetId % 6;
2194 const auto &c = this->underlying().getTetrahedronCoords(tetId);
2195 const auto p{c.data()};
2196
2197 switch(id) {
2198 case 0:
2199 vertexId = getTetrahedronVertexABCG(p, localVertexId);
2200 break;
2201 case 1:
2202 vertexId = getTetrahedronVertexBCDG(p, localVertexId);
2203 break;
2204 case 2:
2205 vertexId = getTetrahedronVertexABEG(p, localVertexId);
2206 break;
2207 case 3:
2208 vertexId = getTetrahedronVertexBEFG(p, localVertexId);
2209 break;
2210 case 4:
2211 vertexId = getTetrahedronVertexBFGH(p, localVertexId);
2212 break;
2213 case 5:
2214 vertexId = getTetrahedronVertexBDGH(p, localVertexId);
2215 break;
2216 }
2217 }
2218 return 0;
2219}
2220
2221template <typename Derived>
2223 const SimplexId &tetId, const int &localEdgeId, SimplexId &edgeId) const {
2224#ifndef TTK_ENABLE_KAMIKAZE
2225 if(tetId < 0 or tetId >= tetrahedronNumber_)
2226 return -1;
2227 if(localEdgeId < 0 or localEdgeId >= 6)
2228 return -2;
2229#endif
2230
2231 edgeId = -1;
2232
2233 if(dimensionality_ == 3) {
2234 const SimplexId id = tetId % 6;
2235 const auto &c = this->underlying().getTetrahedronCoords(tetId);
2236 const auto p{c.data()};
2237
2238 switch(id) {
2239 case 0:
2240 edgeId = getTetrahedronEdgeABCG(p, localEdgeId);
2241 break;
2242 case 1:
2243 edgeId = getTetrahedronEdgeBCDG(p, localEdgeId);
2244 break;
2245 case 2:
2246 edgeId = getTetrahedronEdgeABEG(p, localEdgeId);
2247 break;
2248 case 3:
2249 edgeId = getTetrahedronEdgeBEFG(p, localEdgeId);
2250 break;
2251 case 4:
2252 edgeId = getTetrahedronEdgeBFGH(p, localEdgeId);
2253 break;
2254 case 5:
2255 edgeId = getTetrahedronEdgeBDGH(p, localEdgeId);
2256 break;
2257 }
2258 }
2259
2260 return 0;
2261}
2262
2264 vector<vector<SimplexId>> &edges) const {
2265 edges.resize(tetrahedronNumber_);
2266 for(SimplexId i = 0; i < tetrahedronNumber_; ++i) {
2267 edges[i].resize(6);
2268 for(int j = 0; j < 6; ++j)
2269 getTetrahedronEdge(i, j, edges[i][j]);
2270 }
2271
2272 return 0;
2273}
2274
2275template <typename Derived>
2277 const SimplexId &tetId,
2278 const int &localTriangleId,
2279 SimplexId &triangleId) const {
2280#ifndef TTK_ENABLE_KAMIKAZE
2281 if(tetId < 0 or tetId >= tetrahedronNumber_)
2282 return -1;
2283 if(localTriangleId < 0 or localTriangleId >= 4)
2284 return -2;
2285#endif
2286
2287 triangleId = -1;
2288
2289 if(dimensionality_ == 3) {
2290 const SimplexId id = tetId % 6;
2291 const auto &c = this->underlying().getTetrahedronCoords(tetId);
2292 const auto p{c.data()};
2293
2294 switch(id) {
2295 case 0:
2296 triangleId = getTetrahedronTriangleABCG(p, localTriangleId);
2297 break;
2298 case 1:
2299 triangleId = getTetrahedronTriangleBCDG(p, localTriangleId);
2300 break;
2301 case 2:
2302 triangleId = getTetrahedronTriangleABEG(p, localTriangleId);
2303 break;
2304 case 3:
2305 triangleId = getTetrahedronTriangleBEFG(p, localTriangleId);
2306 break;
2307 case 4:
2308 triangleId = getTetrahedronTriangleBFGH(p, localTriangleId);
2309 break;
2310 case 5:
2311 triangleId = getTetrahedronTriangleBDGH(p, localTriangleId);
2312 break;
2313 }
2314 }
2315
2316 return 0;
2317}
2318
2320 vector<vector<SimplexId>> &triangles) const {
2321 triangles.resize(tetrahedronNumber_);
2322 for(SimplexId i = 0; i < tetrahedronNumber_; ++i) {
2323 triangles[i].resize(4);
2324 for(int j = 0; j < 4; ++j)
2325 getTetrahedronTriangle(i, j, triangles[i][j]);
2326 }
2327
2328 return 0;
2329}
2330
2331template <typename Derived>
2333 const SimplexId &tetId) const {
2334#ifndef TTK_ENABLE_KAMIKAZE
2335 if(tetId < 0 or tetId >= tetrahedronNumber_)
2336 return -1;
2337#endif
2338
2339 if(dimensionality_ == 3) {
2340 const SimplexId id = tetId % 6;
2341 const auto &c = this->underlying().getTetrahedronCoords(tetId);
2342 const auto p{c.data()};
2343
2344 switch(id) {
2345 case 0: // ABCG
2346 if(p[0] == 0 and p[2] == 0)
2347 return 2;
2348 else if(p[0] == 0 or p[2] == 0)
2349 return 3;
2350 else
2351 return 4;
2352 break;
2353 case 1: // BCDG
2354 if(p[1] == nbvoxels_[1] - 1 and p[2] == 0)
2355 return 2;
2356 else if(p[1] == nbvoxels_[1] - 1 or p[2] == 0)
2357 return 3;
2358 else
2359 return 4;
2360 break;
2361 case 2: // ABEG
2362 if(p[0] == 0 and p[1] == 0)
2363 return 2;
2364 else if(p[0] == 0 or p[1] == 0)
2365 return 3;
2366 else
2367 return 4;
2368 break;
2369 case 3: // BEFG
2370 if(p[1] == 0 and p[2] == nbvoxels_[2] - 1)
2371 return 2;
2372 else if(p[1] == 0 or p[2] == nbvoxels_[2] - 1)
2373 return 3;
2374 else
2375 return 4;
2376 break;
2377 case 4: // BFGH
2378 if(p[0] == nbvoxels_[0] - 1 and p[2] == nbvoxels_[2] - 1)
2379 return 2;
2380 else if(p[0] == nbvoxels_[0] - 1 or p[2] == nbvoxels_[2] - 1)
2381 return 3;
2382 else
2383 return 4;
2384 break;
2385 case 5: // BDGH
2386 if(p[0] == nbvoxels_[0] - 1 and p[1] == nbvoxels_[1] - 1)
2387 return 2;
2388 else if(p[0] == nbvoxels_[0] - 1 or p[1] == nbvoxels_[1] - 1)
2389 return 3;
2390 else
2391 return 4;
2392 break;
2393 }
2394 }
2395
2396 return 0;
2397}
2398
2399template <typename Derived>
2401 const SimplexId &tetId,
2402 const int &localNeighborId,
2403 SimplexId &neighborId) const {
2404#ifndef TTK_ENABLE_KAMIKAZE
2405 if(localNeighborId < 0
2406 or localNeighborId >= getTetrahedronNeighborNumber(tetId))
2407 return -1;
2408#endif
2409
2410 neighborId = -1;
2411
2412 if(dimensionality_ == 3) {
2413 const SimplexId id = tetId % 6;
2414 const auto &c = this->underlying().getTetrahedronCoords(tetId);
2415 const auto p{c.data()};
2416
2417 switch(id) {
2418 case 0:
2419 neighborId = getTetrahedronNeighborABCG(tetId, p, localNeighborId);
2420 break;
2421 case 1:
2422 neighborId = getTetrahedronNeighborBCDG(tetId, p, localNeighborId);
2423 break;
2424 case 2:
2425 neighborId = getTetrahedronNeighborABEG(tetId, p, localNeighborId);
2426 break;
2427 case 3:
2428 neighborId = getTetrahedronNeighborBEFG(tetId, p, localNeighborId);
2429 break;
2430 case 4:
2431 neighborId = getTetrahedronNeighborBFGH(tetId, p, localNeighborId);
2432 break;
2433 case 5:
2434 neighborId = getTetrahedronNeighborBDGH(tetId, p, localNeighborId);
2435 break;
2436 }
2437 }
2438
2439 return 0;
2440}
2441
2443 vector<vector<SimplexId>> &neighbors) {
2444 neighbors.resize(tetrahedronNumber_);
2445 for(SimplexId i = 0; i < tetrahedronNumber_; ++i) {
2446 neighbors[i].resize(getTetrahedronNeighborNumber(i));
2447 for(SimplexId j = 0; j < (SimplexId)neighbors[i].size(); ++j)
2448 getTetrahedronNeighbor(i, j, neighbors[i][j]);
2449 }
2450
2451 return 0;
2452}
2453
2454SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(
2455 getCellVertexNumber)(const SimplexId & /*cellId*/) const {
2456 return dimensionality_ + 1;
2457}
2458
2459int ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getCellVertex)(
2460 const SimplexId &cellId,
2461 const int &localVertexId,
2462 SimplexId &vertexId) const {
2463
2464 if(dimensionality_ == 3)
2465 getTetrahedronVertex(cellId, localVertexId, vertexId);
2466 else if(dimensionality_ == 2)
2467 getTriangleVertexInternal(cellId, localVertexId, vertexId);
2468 else if(dimensionality_ == 1)
2469 getEdgeVertexInternal(cellId, localVertexId, vertexId);
2470
2471 return 0;
2472}
2473
2475 const SimplexId & /*cellId*/) const {
2476 if(dimensionality_ == 3)
2477 return 6;
2478 else if(dimensionality_ == 2)
2479 return 3;
2480
2481 return 0;
2482}
2483
2485 const int &localEdgeId,
2486 SimplexId &edgeId) const {
2487 if(dimensionality_ == 3)
2488 getTetrahedronEdge(cellId, localEdgeId, edgeId);
2489 else if(dimensionality_ == 2)
2490 getTriangleEdgeInternal(cellId, localEdgeId, edgeId);
2491 else if(dimensionality_ == 1)
2492 getCellNeighbor(cellId, localEdgeId, edgeId);
2493
2494 return 0;
2495}
2496
2497const vector<vector<SimplexId>> *ImplicitTriangulation::getCellEdgesInternal() {
2498 if(cellEdgeVector_.empty()) {
2499 Timer t;
2500
2501 if(dimensionality_ == 3)
2503 else if(dimensionality_ == 2)
2505
2506 printMsg("Built " + to_string(cellNumber_) + " cell edges.", 1,
2507 t.getElapsedTime(), 1);
2508 }
2509
2510 return &cellEdgeVector_;
2511}
2512
2514 const SimplexId &cellId,
2515 const int &localTriangleId,
2516 SimplexId &triangleId) const {
2517 if(dimensionality_ == 3)
2518 getTetrahedronTriangle(cellId, localTriangleId, triangleId);
2519
2520 return 0;
2521}
2522
2523const vector<vector<SimplexId>> *
2525 if(cellTriangleVector_.empty()) {
2526 Timer t;
2527
2528 if(dimensionality_ == 3)
2530
2531 printMsg("Built " + to_string(cellNumber_) + " cell triangles.", 1,
2532 t.getElapsedTime(), 1);
2533 }
2534
2535 return &cellTriangleVector_;
2536}
2537
2538SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(
2539 getCellNeighborNumber)(const SimplexId &cellId) const {
2540 if(dimensionality_ == 3)
2541 return getTetrahedronNeighborNumber(cellId);
2542 else if(dimensionality_ == 2)
2543 return getTriangleNeighborNumber(cellId);
2544 else if(dimensionality_ == 1) {
2545 printErr("getCellNeighborNumber() not implemented in 1D! (TODO)");
2546 return -1;
2547 }
2548
2549 return 0;
2550}
2551
2552int ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getCellNeighbor)(
2553 const SimplexId &cellId,
2554 const int &localNeighborId,
2555 SimplexId &neighborId) const {
2556 if(dimensionality_ == 3)
2557 getTetrahedronNeighbor(cellId, localNeighborId, neighborId);
2558 else if(dimensionality_ == 2)
2559 getTriangleNeighbor(cellId, localNeighborId, neighborId);
2560 else if(dimensionality_ == 1) {
2561 printErr("getCellNeighbor() not implemented in 1D! (TODO)");
2562 return -1;
2563 }
2564
2565 return 0;
2566}
2567
2568const vector<vector<SimplexId>> *
2569 ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL(getCellNeighbors)() {
2570 if(cellNeighborList_.empty()) {
2571 Timer t;
2572
2573 if(dimensionality_ == 3)
2574 getTetrahedronNeighbors(cellNeighborList_);
2575 else if(dimensionality_ == 2)
2576 getTriangleNeighbors(cellNeighborList_);
2577 else if(dimensionality_ == 1) {
2578 printErr("getCellNeighbors() not implemented in 1D! (TODO)");
2579 return nullptr;
2580 }
2581
2582 printMsg("Built " + to_string(cellNumber_) + " cell neighbors.", 1,
2583 t.getElapsedTime(), 1);
2584 }
2585
2586 return &cellNeighborList_;
2587}
2588
2590 // V(abcdefgh)=V(g)+V(d)::{g,h}+V(h)::{g}+V(b)::{c,d,g,h}
2591 this->vertexNeighborABCDEFGH_ = {
2592 -vshift_[0] - vshift_[1], // a
2593 1 - vshift_[0] - vshift_[1], // b
2594 -vshift_[1], // c
2595 1 - vshift_[1], // d
2596 -vshift_[0], // e
2597 1 - vshift_[0], // f
2598 1, // h
2599 -1 + vshift_[1], // V(d)::{g}
2600 vshift_[1], // V(d)::{h}
2601 -1, // V(h)::{g}
2602 -1 + vshift_[0], // V(b)::{c}
2603 vshift_[0], // V(b)::{d}
2604 -1 + vshift_[0] + vshift_[1], // V(b)::{g}
2605 vshift_[0] + vshift_[1] // V(b)::{h}
2606 };
2607
2608 // V(abdc)=V(b)+V(d)::{b}+V(c)::{b}+V(a)::{b}
2609 this->vertexNeighborABCD_ = {
2610 -1, // a
2611 -1 + vshift_[0], // c
2612 vshift_[0], // d
2613 -1 + vshift_[1], // e
2614 vshift_[1], // f
2615 -1 + vshift_[0] + vshift_[1], // g
2616 vshift_[0] + vshift_[1], // h
2617 -vshift_[0], // V(d)::{b}
2618 1 - vshift_[0], // V(c)::{b}
2619 1, // V(a)::{b}
2620 };
2621 // V(efhg)=V(g)+V(h)::{g}+V(f)::{g,h}
2622 this->vertexNeighborEFGH_ = {
2623 -vshift_[0] - vshift_[1], // a
2624 1 - vshift_[0] - vshift_[1], // b
2625 -vshift_[1], // c
2626 1 - vshift_[1], // d
2627 -vshift_[0], // e
2628 1 - vshift_[0], // f
2629 1, // h
2630 -1, // V(h)::{g}
2631 -1 + vshift_[0], // V(f)::{g}
2632 vshift_[0], // V(f)::{h}
2633 };
2634 // V(aefb)=V(b)+V(a)::{b}+V(e)::{b}+V(f)::{b}
2635 this->vertexNeighborAEFB_ = {
2636 -1, // a
2637 -1 + vshift_[0], // c
2638 vshift_[0], // d
2639 -1 + vshift_[1], // e
2640 vshift_[1], // f
2641 -1 + vshift_[0] + vshift_[1], // g
2642 vshift_[0] + vshift_[1], // h
2643 1, // V(a)::{b}
2644 1 - vshift_[1], // V(e)::{b}
2645 -vshift_[1], // V(f)::{b}
2646 };
2647 // V(ghdc)=V(g)+V(h)::{g}+V(d)::{g,h}
2648 this->vertexNeighborGHDC_ = {
2649 -vshift_[0] - vshift_[1], // a
2650 1 - vshift_[0] - vshift_[1], // b
2651 -vshift_[1], // c
2652 1 - vshift_[1], // d
2653 -vshift_[0], // e
2654 1 - vshift_[0], // f
2655 1, // h
2656 -1, // V(h)::{g}
2657 -1 + vshift_[1], // V(d)::{g}
2658 vshift_[1] // V(d)::{h}
2659 };
2660 // V(aegc)=V(g)+V(a)::{c,g}+V(c)::{g}
2661 this->vertexNeighborAEGC_ = {
2662 -vshift_[0] - vshift_[1], // a
2663 1 - vshift_[0] - vshift_[1], // b
2664 -vshift_[1], // c
2665 1 - vshift_[1], // d
2666 -vshift_[0], // e
2667 1 - vshift_[0], // f
2668 1, // h
2669 vshift_[0], // V(a)::{c}
2670 vshift_[0] + vshift_[1], // V(a)::{g}
2671 vshift_[1], // V(c)::{g}
2672 };
2673 // V(bfhd)=V(b)+V(f)::{b}+V(h)::{b}+V(d)::{b}
2674 this->vertexNeighborBFHD_ = {
2675 -1, // a
2676 -1 + vshift_[0], // c
2677 vshift_[0], // d
2678 -1 + vshift_[1], // e
2679 vshift_[1], // f
2680 -1 + vshift_[0] + vshift_[1], // g
2681 vshift_[0] + vshift_[1], // h
2682 -vshift_[1], // V(f)::{b}
2683 -vshift_[0] - vshift_[1], // V(h)::{b}
2684 -vshift_[0], // V(d)::{b}
2685 };
2686
2687 // V(ab)=V(b)+V(a)::{b}
2688 this->vertexNeighborAB_ = {
2689 -1, // a
2690 -1 + vshift_[0], // c
2691 vshift_[0], // d
2692 -1 + vshift_[1], // e
2693 vshift_[1], // f
2694 -1 + vshift_[0] + vshift_[1], // g
2695 vshift_[0] + vshift_[1], // h
2696 1, // V(a)::{b}
2697 };
2698 // V(bd)=V(b)+V(d)::{b}
2699 this->vertexNeighborBD_ = {
2700 -1, // a
2701 -1 + vshift_[0], // c
2702 vshift_[0], // d
2703 -1 + vshift_[1], // e
2704 vshift_[1], // f
2705 -1 + vshift_[0] + vshift_[1], // g
2706 vshift_[0] + vshift_[1], // h
2707 -vshift_[0], // V(d)::{b}
2708 };
2709 // V(gh)=V(g)+V(h)::{g}
2710 this->vertexNeighborGH_ = {
2711 -vshift_[0] - vshift_[1], // a
2712 1 - vshift_[0] - vshift_[1], // b
2713 -vshift_[1], // c
2714 1 - vshift_[1], // d
2715 -vshift_[0], // e
2716 1 - vshift_[0], // f
2717 1, // h
2718 -1, // V(h)::{g}
2719 };
2720 // V(eg)=V(g)+V(e)::{g}
2721 this->vertexNeighborEG_ = {
2722 -vshift_[0] - vshift_[1], // a
2723 1 - vshift_[0] - vshift_[1], // b
2724 -vshift_[1], // c
2725 1 - vshift_[1], // d
2726 -vshift_[0], // e
2727 1 - vshift_[0], // f
2728 1, // h
2729 vshift_[0], // V(e)::{g}
2730 };
2731 // V(cg)=V(g)+V(c)::{g}
2732 this->vertexNeighborCG_ = {
2733 -vshift_[0] - vshift_[1], // a
2734 1 - vshift_[0] - vshift_[1], // b
2735 -vshift_[1], // c
2736 1 - vshift_[1], // d
2737 -vshift_[0], // e
2738 1 - vshift_[0], // f
2739 1, // h
2740 vshift_[1], // V(c)::{g}
2741 };
2742 // V(bf)=V(b)+V(f)::{b}
2743 this->vertexNeighborBF_ = {
2744 -1, // a
2745 -1 + vshift_[0], // c
2746 vshift_[0], // d
2747 -1 + vshift_[1], // e
2748 vshift_[1], // f
2749 -1 + vshift_[0] + vshift_[1], // g
2750 vshift_[0] + vshift_[1], // h
2751 -vshift_[1], // V(f)::{b}
2752 };
2753
2754 // V(b)={a,c,d,e,f,g,h}
2755 this->vertexNeighborB_ = {
2756 -1, // a
2757 -1 + vshift_[0], // c
2758 vshift_[0], // d
2759 -1 + vshift_[1], // e
2760 vshift_[1], // f
2761 -1 + vshift_[0] + vshift_[1], // g
2762 vshift_[0] + vshift_[1], // h
2763 };
2764 // V(g)={a,b,c,d,e,f,h}
2765 this->vertexNeighborG_ = {
2766 -vshift_[0] - vshift_[1], // a
2767 1 - vshift_[0] - vshift_[1], // b
2768 -vshift_[1], // c
2769 1 - vshift_[1], // d
2770 -vshift_[0], // e
2771 1 - vshift_[0], // f
2772 1, // h
2773 };
2774
2775 // V(ef)=V(f)+V(e)::{b,f}
2776 this->vertexNeighborEF_ = {
2777 -vshift_[1], // b
2778 -1, // e
2779 -1 + vshift_[0], // g
2780 vshift_[0], // h
2781 1 - vshift_[1], // V(e)::{b}
2782 1, // V(e)::{f}
2783 };
2784 // V(cd)=V(d)+V(c)::{b,d}
2785 this->vertexNeighborCD_ = {
2786 -vshift_[0], // b
2787 -1, // c
2788 -1 + vshift_[1], // g
2789 vshift_[1], // h
2790 1 - vshift_[0], // V(c)::{b}
2791 1, // V(c)::{d}
2792 };
2793 // V(ac)=V(c)+V(a)::{c,g}
2794 this->vertexNeighborAC_ = {
2795 -vshift_[0], // a
2796 1 - vshift_[0], // b
2797 1, // d
2798 vshift_[1], // g
2799 vshift_[0], // V(a)::{c}
2800 vshift_[0] + vshift_[1], // V(a)::{c}
2801 };
2802 // V(ae)=V(a)+V(e)::{a,b}
2803 this->vertexNeighborAE_ = {
2804 1, // b
2805 vshift_[0], // c
2806 vshift_[1], // e
2807 vshift_[0] + vshift_[1], // g
2808 -vshift_[1], // V(e)::{a}
2809 1 - vshift_[1], // V(e)::{b}
2810 };
2811 // V(fh)=V(f)+V(h)::{b,f}
2812 this->vertexNeighborFH_ = {
2813 -vshift_[1], // b
2814 -1, // e
2815 -1 + vshift_[0], // g
2816 vshift_[0], // h
2817 -vshift_[0] - vshift_[1], // V(h)::{b}
2818 -vshift_[0], // V(h)::{f}
2819 };
2820 // V(dh)=V(d)+V(h)::{b,d}
2821 this->vertexNeighborDH_ = {
2822 -vshift_[0], // b
2823 -1, // c
2824 -1 + vshift_[1], // g
2825 vshift_[1], // h
2826 -vshift_[0] - vshift_[1], // V(h)::{b}
2827 -vshift_[1], // V(h)::{d}
2828 };
2829
2830 // V(a)={b,c,e,g}
2831 this->vertexNeighborA_ = {
2832 1, // b
2833 vshift_[0], // c
2834 vshift_[1], // e
2835 vshift_[0] + vshift_[1], // g
2836 };
2837 // V(c)={a,b,d,g}
2838 this->vertexNeighborC_ = {
2839 -vshift_[0], // a
2840 1 - vshift_[0], // b
2841 1, // d
2842 +vshift_[1], // g
2843 };
2844 // V(d)={b,c,g,h}
2845 this->vertexNeighborD_ = {
2846 -vshift_[0], // b
2847 -1, // c
2848 -1 + vshift_[1], // g
2849 vshift_[1], // h
2850 };
2851 // V(e)={a,b,f,g}
2852 this->vertexNeighborE_ = {
2853 -vshift_[1], // a
2854 1 - vshift_[1], // b
2855 1, // f
2856 +vshift_[0], // g
2857 };
2858 // V(f)={b,e,g,h}
2859 this->vertexNeighborF_ = {
2860 -vshift_[1], // b
2861 -1, // e
2862 -1 + vshift_[0], // g
2863 vshift_[0], // h
2864 };
2865 // V(h)={b,d,f,g}
2866 this->vertexNeighborH_ = {
2867 -vshift_[0] - vshift_[1], // b
2868 -vshift_[1], // d
2869 -vshift_[0], // f
2870 -1, // g
2871 };
2872
2873 // V(abcd)=V(d)::{b,c}+V(c)::{b,d}+V(a)::{b}+V(b)::{c}
2874 this->vertexNeighbor2dABCD_ = {
2875 -1, -vshift_[0], -vshift_[0] + 1, 1, vshift_[0], vshift_[0] - 1,
2876 };
2877 // V(ab)=V(b)::{a,c,d}+V(a)::{b}
2878 this->vertexNeighbor2dAB_ = {
2879 -1, // V(b)::a
2880 vshift_[0] - 1, // V(b)::c
2881 vshift_[0], // V(b)::d
2882 +1, // V(a)::b
2883 };
2884 // V(cd)=V(c)::{a,b,d}+V(d)::{c}
2885 this->vertexNeighbor2dCD_ = {
2886 -1, // V(d)::c
2887 -vshift_[0], // V(c)::a
2888 -vshift_[0] + 1, // V(c)::b
2889 1, // V(c)::d
2890 };
2891 // V(ac)=V(c)::{a,b,d}+V(a)::{c}
2892 this->vertexNeighbor2dAC_ = {
2893 -vshift_[0], // V(c)::{a}
2894 -vshift_[0] + 1, // V(c)::{b}
2895 1, // V(c)::{d}
2896 vshift_[0], // V(a)::{c}
2897 };
2898 // V(bd)=V(b)::{c,d}+V(d)::{b,c}
2899 this->vertexNeighbor2dBD_ = {
2900 vshift_[0] - 1, // V(b)::{c}
2901 vshift_[0], // V(b)::{d}
2902 -vshift_[0], // V(d)::{b}
2903 -1, // V(d)::{c}
2904 };
2905 // V(b)={a,c,d}
2906 this->vertexNeighbor2dB_ = {
2907 -1, // a
2908 vshift_[0], // d
2909 vshift_[0] - 1, // c
2910 };
2911 // V(c)={a,b,d}
2912 this->vertexNeighbor2dC_ = {
2913 1, // d
2914 -vshift_[0], // a
2915 -vshift_[0] + 1, // b
2916 };
2917 this->vertexNeighbor2dA_ = {};
2918 // V(a)={b,c}
2919 this->vertexNeighbor2dA_ = {
2920 1, // b
2921 vshift_[0] // c
2922 };
2923 // V(d)={c,b}
2924 this->vertexNeighbor2dD_ = {
2925 -1, // c
2926 -vshift_[0], // b
2927
2928 };
2929
2930 return 0;
2931}
2932
2933#ifdef TTK_ENABLE_MPI
2934
2935int ttk::ImplicitTriangulation::preconditionDistributedCells() {
2936 if(this->hasPreconditionedDistributedCells_) {
2937 return 0;
2938 }
2939 if(!ttk::hasInitializedMPI()) {
2940 return -1;
2941 }
2942 if(this->cellGhost_ == nullptr) {
2943 if(ttk::isRunningWithMPI()) {
2944 this->printErr("Missing cell ghost array!");
2945 }
2946 return -3;
2947 }
2948
2949 Timer tm{};
2950
2951 // number of local cells (with ghost cells...)
2952 const auto nLocCells{this->getNumberOfCells()};
2953
2954 // there are 6 tetrahedra per cubic cell (and 2 triangles per square)
2955 const int nTetraPerCube{this->dimensionality_ == 3 ? 6 : 2};
2956 std::vector<unsigned char> fillCells(nLocCells / nTetraPerCube);
2957
2958 this->ghostCellsPerOwner_.resize(ttk::MPIsize_);
2959
2960 this->neighborCellBBoxes_.resize(ttk::MPIsize_);
2961 auto &localBBox{this->neighborCellBBoxes_[ttk::MPIrank_]};
2962
2963 ttk::SimplexId localBBox_x_min{this->localGridOffset_[0]
2964 + this->dimensions_[0]},
2965 localBBox_y_min{this->localGridOffset_[1] + this->dimensions_[1]},
2966 localBBox_z_min{this->localGridOffset_[2] + this->dimensions_[2]};
2967 ttk::SimplexId localBBox_x_max{this->localGridOffset_[0]},
2968 localBBox_y_max{this->localGridOffset_[1]},
2969 localBBox_z_max{this->localGridOffset_[2]};
2970
2971#ifdef TTK_ENABLE_OPENMP
2972#pragma omp parallel for reduction( \
2973 min \
2974 : localBBox_x_min, localBBox_y_min, localBBox_z_min) \
2975 reduction(max \
2976 : localBBox_x_max, localBBox_y_max, localBBox_z_max)
2977#endif
2978 for(SimplexId lcid = 0; lcid < nLocCells; ++lcid) {
2979 // only keep non-ghost cells
2980 if(this->cellGhost_[lcid / nTetraPerCube] == 1) {
2981 continue;
2982 }
2983 // local vertex coordinates
2984 std::array<SimplexId, 3> p{};
2985 if(this->dimensionality_ == 3) {
2986 this->tetrahedronToPosition(lcid, p.data());
2987 } else if(this->dimensionality_ == 2) {
2988 this->triangleToPosition2d(lcid, p.data());
2989 // compatibility with tetrahedronToPosition; fix a bounding box
2990 // error in the first axis
2991 p[0] /= 2;
2992 }
2993
2994 // global vertex coordinates
2995 p[0] += this->localGridOffset_[0];
2996 p[1] += this->localGridOffset_[1];
2997 p[2] += this->localGridOffset_[2];
2998
2999 if(p[0] < localBBox_x_min) {
3000 localBBox_x_min = p[0];
3001 }
3002 if(p[0] > localBBox_x_max) {
3003 localBBox_x_max = p[0];
3004 }
3005 if(p[1] < localBBox_y_min) {
3006 localBBox_y_min = p[1];
3007 }
3008 if(p[1] > localBBox_y_max) {
3009 localBBox_y_max = p[1];
3010 }
3011 if(p[2] < localBBox_z_min) {
3012 localBBox_z_min = p[2];
3013 }
3014 if(p[2] > localBBox_z_max) {
3015 localBBox_z_max = p[2];
3016 }
3017 }
3018 localBBox_x_max++;
3019 localBBox_y_max++;
3020 localBBox_z_max++;
3021
3022 localBBox = {
3023 localBBox_x_min, localBBox_x_max, localBBox_y_min,
3024 localBBox_y_max, localBBox_z_min, localBBox_z_max,
3025 };
3026
3027 for(size_t i = 0; i < this->neighborRanks_.size(); ++i) {
3028 const auto neigh{this->neighborRanks_[i]};
3029 MPI_Sendrecv(this->neighborCellBBoxes_[ttk::MPIrank_].data(), 6,
3030 ttk::getMPIType(SimplexId{}), neigh, ttk::MPIrank_,
3031 this->neighborCellBBoxes_[neigh].data(), 6,
3032 ttk::getMPIType(SimplexId{}), neigh, neigh, ttk::MPIcomm_,
3033 MPI_STATUS_IGNORE);
3034 }
3035
3036 this->hasPreconditionedDistributedCells_ = true;
3037
3038 return 0;
3039}
3040
3041void ttk::ImplicitTriangulation::createMetaGrid(const double *const bounds) {
3042
3043 // only works with 2 processes or more
3044 if(!ttk::isRunningWithMPI()) {
3045 return;
3046 }
3047
3048 // no need to create it anew?
3049 if(this->metaGrid_ != nullptr) {
3050 return;
3051 }
3052
3053 // Reorganize bounds to only execute Allreduce twice
3054 std::array<double, 6> tempBounds = {
3055 bounds[0], bounds[2], bounds[4], bounds[1], bounds[3], bounds[5],
3056 };
3057 std::array<double, 6> tempGlobalBounds{};
3058
3059 // Compute and send to all processes the lower bounds of the data set
3060 MPI_Allreduce(tempBounds.data(), tempGlobalBounds.data(), 3, MPI_DOUBLE,
3061 MPI_MIN, ttk::MPIcomm_);
3062 // Compute and send to all processes the higher bounds of the data set
3063 MPI_Allreduce(&tempBounds[3], &tempGlobalBounds[3], 3, MPI_DOUBLE, MPI_MAX,
3064 ttk::MPIcomm_);
3065
3066 // re-order tempGlobalBounds
3067 std::array<double, 6> globalBounds{
3068 tempGlobalBounds[0], tempGlobalBounds[3], tempGlobalBounds[1],
3069 tempGlobalBounds[4], tempGlobalBounds[2], tempGlobalBounds[5],
3070 };
3071
3072 const std::array<int, 3> dimensions = {
3073 static_cast<int>(
3074 std::round((globalBounds[1] - globalBounds[0]) / this->spacing_[0]))
3075 + 1,
3076 static_cast<int>(
3077 std::round((globalBounds[3] - globalBounds[2]) / this->spacing_[1]))
3078 + 1,
3079 static_cast<int>(
3080 std::round((globalBounds[5] - globalBounds[4]) / this->spacing_[2]))
3081 + 1,
3082 };
3083
3084 this->localGridOffset_ = {
3085 static_cast<SimplexId>(
3086 std::round((this->origin_[0] - globalBounds[0]) / this->spacing_[0])),
3087 static_cast<SimplexId>(
3088 std::round((this->origin_[1] - globalBounds[2]) / this->spacing_[1])),
3089 static_cast<SimplexId>(
3090 std::round((this->origin_[2] - globalBounds[4]) / this->spacing_[2])),
3091 };
3092
3093 this->metaGrid_ = std::make_shared<ImplicitNoPreconditions>();
3094 this->metaGrid_->setInputGrid(globalBounds[0], globalBounds[1],
3095 globalBounds[2], this->spacing_[0],
3096 this->spacing_[1], this->spacing_[2],
3097 dimensions[0], dimensions[1], dimensions[2]);
3098 this->metaGrid_->preconditionBoundaryVertices();
3099 this->metaGrid_->preconditionBoundaryEdges();
3100 this->metaGrid_->preconditionBoundaryTriangles();
3101}
3102
3103std::array<SimplexId, 3>
3104 ttk::ImplicitTriangulation::getVertGlobalCoords(const SimplexId lvid) const {
3105
3106 // local vertex coordinates
3107 std::array<SimplexId, 3> p{};
3108 if(this->dimensionality_ == 3) {
3109 this->vertexToPosition(lvid, p.data());
3110 } else if(this->dimensionality_ == 2) {
3111 this->vertexToPosition2d(lvid, p.data());
3112 }
3113
3114 // global vertex coordinates
3115 p[0] += this->localGridOffset_[0];
3116 p[1] += this->localGridOffset_[1];
3117 p[2] += this->localGridOffset_[2];
3118
3119 return p;
3120}
3121
3122std::array<SimplexId, 3>
3123 ttk::ImplicitTriangulation::getVertLocalCoords(const SimplexId gvid) const {
3124
3125 // global vertex coordinates
3126 std::array<SimplexId, 3> p{};
3127 if(this->dimensionality_ == 3) {
3128 this->metaGrid_->vertexToPosition(gvid, p.data());
3129 } else if(this->dimensionality_ == 2) {
3130 this->metaGrid_->vertexToPosition2d(gvid, p.data());
3131 }
3132
3133 // local vertex coordinates
3134 p[0] -= this->localGridOffset_[0];
3135 p[1] -= this->localGridOffset_[1];
3136 p[2] -= this->localGridOffset_[2];
3137
3138 return p;
3139}
3140
3141bool ImplicitTriangulation::isVertexOnGlobalBoundaryInternal(
3142 const SimplexId lvid) const {
3143
3144 if(!ttk::isRunningWithMPI()) {
3145 return this->isVertexOnBoundary(lvid);
3146 }
3147
3148#ifndef TTK_ENABLE_KAMIKAZE
3149 if(lvid > this->TTK_TRIANGULATION_INTERNAL(getNumberOfVertices)() - 1
3150 || lvid < 0) {
3151 return false;
3152 }
3153 if(this->metaGrid_ == nullptr) {
3154 return false;
3155 }
3156#endif // TTK_ENABLE_KAMIKAZE
3157
3158 const auto gvid{this->getVertexGlobalIdInternal(lvid)};
3159 if(gvid == -1) {
3160 return false;
3161 }
3162 return this->metaGrid_->isVertexOnBoundary(gvid);
3163}
3164
3165bool ImplicitTriangulation::isEdgeOnGlobalBoundaryInternal(
3166 const SimplexId leid) const {
3167
3168 if(!ttk::isRunningWithMPI()) {
3169 return this->isEdgeOnBoundary(leid);
3170 }
3171
3172#ifndef TTK_ENABLE_KAMIKAZE
3173 if(leid > this->getNumberOfEdgesInternal() - 1 || leid < 0) {
3174 return false;
3175 }
3176 if(this->metaGrid_ == nullptr) {
3177 return false;
3178 }
3179#endif // TTK_ENABLE_KAMIKAZE
3180
3181 const auto geid{this->getEdgeGlobalIdInternal(leid)};
3182 if(geid == -1) {
3183 return false;
3184 }
3185 return this->metaGrid_->isEdgeOnBoundary(geid);
3186}
3187
3188bool ImplicitTriangulation::isTriangleOnGlobalBoundaryInternal(
3189 const SimplexId ltid) const {
3190
3191 if(!ttk::isRunningWithMPI()) {
3192 return this->isTriangleOnBoundary(ltid);
3193 }
3194
3195#ifndef TTK_ENABLE_KAMIKAZE
3196 if(ltid > this->getNumberOfTrianglesInternal() - 1 || ltid < 0) {
3197 return false;
3198 }
3199 if(this->metaGrid_ == nullptr) {
3200 return false;
3201 }
3202#endif // TTK_ENABLE_KAMIKAZE
3203
3204 const auto gtid{this->getTriangleGlobalIdInternal(ltid)};
3205 if(gtid == -1) {
3206 return false;
3207 }
3208 return this->metaGrid_->isTriangleOnBoundary(gtid);
3209}
3210
3211int ttk::ImplicitTriangulation::getCellRankInternal(
3212 const SimplexId lcid) const {
3213
3214 const int nTetraPerCube{this->dimensionality_ == 3 ? 6 : 2};
3215 const auto locCubeId{lcid / nTetraPerCube};
3216
3217 if(this->cellGhost_[locCubeId] == 0) {
3218 return ttk::MPIrank_;
3219 }
3220
3221#ifndef TTK_ENABLE_KAMIKAZE
3222 if(this->neighborRanks_.empty() && ttk::MPIsize_ > 1) {
3223 this->printErr("Empty neighborsRanks_!");
3224 return -1;
3225 }
3226#endif // TTK_ENABLE_KAMIKAZE
3227
3228 const auto nVertsCell{this->getCellVertexNumber(lcid)};
3229 std::vector<bool> inRank(nVertsCell);
3230 for(const auto neigh : this->neighborRanks_) {
3231 std::fill(inRank.begin(), inRank.end(), false);
3232 const auto &bbox{this->neighborCellBBoxes_[neigh]};
3233 for(SimplexId i = 0; i < nVertsCell; ++i) {
3234 SimplexId v{};
3235 this->getCellVertex(lcid, i, v);
3236 if(this->vertexGhost_[v] == 0) {
3237 inRank[i] = true;
3238 } else {
3239 const auto p{this->getVertGlobalCoords(v)};
3240 if(p[0] >= bbox[0] && p[0] <= bbox[1] && p[1] >= bbox[2]
3241 && p[1] <= bbox[3] && p[2] >= bbox[4] && p[2] <= bbox[5]) {
3242 inRank[i] = true;
3243 }
3244 }
3245 }
3246 if(std::all_of(
3247 inRank.begin(), inRank.end(), [](const bool v) { return v; })) {
3248 return neigh;
3249 }
3250 }
3251
3252 return -1;
3253}
3254
3255#endif // TTK_ENABLE_MPI
3256
3257// explicit instantiations
#define TTK_TRIANGULATION_INTERNAL(NAME)
bool ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() isTriangleOnBoundary(const SimplexId &triangleId) const
#define CASE_EDGE_POSITION_P_3D
int ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getEdgeStar(const SimplexId &edgeId, const int &localStarId, SimplexId &starId) const
const vector< vector< SimplexId > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getCellNeighbors()
SimplexId ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getTriangleStarNumber(const SimplexId &triangleId) const
int ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getCellNeighbor(const SimplexId &cellId, const int &localNeighborId, SimplexId &neighborId) const
int ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getVertexStar(const SimplexId &vertexId, const int &localStarId, SimplexId &starId) const
SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getTriangleLinkNumber(const SimplexId &triangleId) const
int ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getEdgeLink(const SimplexId &edgeId, const int &localLinkId, SimplexId &linkId) const
#define CASE_EDGE_POSITION_D2_3D
#define CASE_EDGE_POSITION_D1_3D
int ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getTriangleLink(const SimplexId &triangleId, const int &localLinkId, SimplexId &linkId) const
SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getCellVertexNumber(const SimplexId &) const
SimplexId ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getEdgeStarNumber(const SimplexId &edgeId) const
#define CASE_EDGE_POSITION_L_2D
#define CASE_EDGE_POSITION_H_2D
#define CASE_EDGE_POSITION_H_3D
const vector< vector< SimplexId > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getTriangleLinks()
#define CASE_EDGE_POSITION_L_3D
const vector< vector< SimplexId > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getTriangleStars()
SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getVertexLinkNumber(const SimplexId &vertexId) const
SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getCellNeighborNumber(const SimplexId &cellId) const
SimplexId ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getVertexStarNumber(const SimplexId &vertexId) const
const vector< vector< SimplexId > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getVertexLinks()
const vector< vector< SimplexId > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getEdgeStars()
#define CASE_EDGE_POSITION_D3_3D
bool ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() isEdgeOnBoundary(const SimplexId &edgeId) const
int ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getCellVertex(const SimplexId &cellId, const int &localVertexId, SimplexId &vertexId) const
SimplexId ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getEdgeLinkNumber(const SimplexId &edgeId) const
const vector< vector< SimplexId > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getVertexNeighbors()
int ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getTriangleStar(const SimplexId &triangleId, const int &localStarId, SimplexId &starId) const
const vector< std::array< SimplexId, 3 > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getTriangles()
const vector< vector< SimplexId > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getVertexStars()
int ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getVertexPoint(const SimplexId &vertexId, float &x, float &y, float &z) const
const vector< std::array< SimplexId, 2 > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getEdges()
bool ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() isVertexOnBoundary(const SimplexId &vertexId) const
int ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getVertexLink(const SimplexId &vertexId, const int &localLinkId, SimplexId &linkId) const
const vector< vector< SimplexId > > *ImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getEdgeLinks()
SimplexId PeriodicImplicitTriangulation::TTK_TRIANGULATION_INTERNAL() getVertexNeighborNumber(const SimplexId &vertexId) const
int PeriodicImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getVertexNeighbor(const SimplexId &vertexId, const int &localNeighborId, SimplexId &neighborId) const
virtual SimplexId getVertexTriangleNumberInternal(const SimplexId &ttkNotUsed(vertexId)) const
virtual SimplexId getVertexNeighborNumber(const SimplexId &vertexId) const
virtual bool isVertexOnBoundary(const SimplexId &vertexId) const
virtual int getVertexTriangleInternal(const SimplexId &ttkNotUsed(vertexId), const int &ttkNotUsed(localTriangleId), SimplexId &ttkNotUsed(triangleId)) const
virtual bool isEdgeOnBoundary(const SimplexId &edgeId) const
std::vector< std::vector< SimplexId > > vertexEdgeList_
std::vector< std::vector< SimplexId > > triangleEdgeVector_
virtual int getTriangleEdgeInternal(const SimplexId &ttkNotUsed(triangleId), const int &ttkNotUsed(localEdgeId), SimplexId &ttkNotUsed(edgeId)) const
virtual SimplexId getEdgeTriangleNumberInternal(const SimplexId &ttkNotUsed(edgeId)) const
virtual int getEdgeTriangleInternal(const SimplexId &ttkNotUsed(edgeId), const int &ttkNotUsed(localTriangleId), SimplexId &ttkNotUsed(triangleId)) const
virtual int getVertexEdgeInternal(const SimplexId &ttkNotUsed(vertexId), const int &ttkNotUsed(localEdgeId), SimplexId &ttkNotUsed(edgeId)) const
std::vector< std::vector< SimplexId > > vertexTriangleList_
std::vector< std::vector< SimplexId > > cellEdgeVector_
std::vector< std::vector< SimplexId > > cellTriangleVector_
std::vector< std::vector< SimplexId > > edgeTriangleList_
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
int getTriangleVertexInternal(const SimplexId &triangleId, const int &localVertexId, SimplexId &vertexId) const override
int getTetrahedronVertex(const SimplexId &tetId, const int &localVertexId, SimplexId &vertexId) const override
int getTetrahedronTriangle(const SimplexId &tetId, const int &id, SimplexId &triangleId) const override
SimplexId getEdgeTriangleNumberInternal(const SimplexId &edgeId) const override
int getVertexEdgeInternal(const SimplexId &vertexId, const int &id, SimplexId &edgeId) const override
int getEdgeVertexInternal(const SimplexId &edgeId, const int &localVertexId, SimplexId &vertexId) const override
SimplexId getVertexTriangleNumberInternal(const SimplexId &vertexId) const override
SimplexId getTriangleNeighborNumber(const SimplexId &triangleId) const override
int getEdgeTriangleInternal(const SimplexId &edgeId, const int &id, SimplexId &triangleId) const override
int getTetrahedronNeighbor(const SimplexId &tetId, const int &localNeighborId, SimplexId &neighborId) const override
int getVertexTriangleInternal(const SimplexId &vertexId, const int &id, SimplexId &triangleId) const override
SimplexId getTetrahedronNeighborNumber(const SimplexId &tetId) const override
int getTriangleNeighbor(const SimplexId &triangleId, const int &localNeighborId, SimplexId &neighborId) const override
int getTriangleEdgeInternal(const SimplexId &triangleId, const int &id, SimplexId &edgeId) const override
int getTetrahedronEdge(const SimplexId &tetId, const int &id, SimplexId &edgeId) const override
std::array< SimplexId, 8 > vertexNeighborGH_
std::array< SimplexId, 4 > vertexNeighbor2dAC_
std::array< SimplexId, 8 > vertexNeighborBF_
std::array< SimplexId, 8 > vertexNeighborBD_
std::array< SimplexId, 3 > vertexNeighbor2dB_
std::array< SimplexId, 7 > vertexNeighborB_
std::array< SimplexId, 6 > vertexNeighborEF_
std::array< SimplexId, 14 > vertexNeighborABCDEFGH_
const std::vector< std::vector< SimplexId > > * getTriangleEdgesInternal() override
SimplexId getCellEdgeNumberInternal(const SimplexId &cellId) const override
std::array< SimplexId, 8 > vertexNeighborEG_
int getCellTriangleInternal(const SimplexId &cellId, const int &id, SimplexId &triangleId) const override
std::array< SimplexId, 4 > vertexNeighborH_
std::array< SimplexId, 4 > vertexNeighborE_
std::array< SimplexId, 2 > vertexNeighbor2dD_
bool TTK_TRIANGULATION_INTERNAL() isTriangleOnBoundary(const SimplexId &triangleId) const override
const std::vector< std::vector< SimplexId > > * getVertexEdgesInternal() override
std::array< SimplexId, 10 > vertexNeighborAEFB_
virtual int getTriangleNeighbor(const SimplexId &triangleId, const int &localNeighborId, SimplexId &neighborId) const =0
std::array< SimplexId, 6 > vertexNeighbor2dABCD_
std::array< SimplexId, 6 > vertexNeighborFH_
std::array< SimplexId, 6 > vertexNeighborDH_
std::array< SimplexId, 6 > vertexNeighborCD_
const std::vector< std::vector< SimplexId > > * getCellEdgesInternal() override
std::array< SimplexId, 7 > vertexNeighborG_
int setInputGrid(const float &xOrigin, const float &yOrigin, const float &zOrigin, const float &xSpacing, const float &ySpacing, const float &zSpacing, const SimplexId &xDim, const SimplexId &yDim, const SimplexId &zDim) override
int getTetrahedronNeighbors(std::vector< std::vector< SimplexId > > &neighbors)
int TTK_TRIANGULATION_INTERNAL() getCellNeighbor(const SimplexId &cellId, const int &localNeighborId, SimplexId &neighborId) const override
const std::vector< std::vector< SimplexId > > * getVertexTrianglesInternal() override
virtual SimplexId getTetrahedronNeighborNumber(const SimplexId &tetId) const =0
virtual int getTetrahedronNeighbor(const SimplexId &tetId, const int &localNeighborId, SimplexId &neighborId) const =0
std::array< SimplexId, 10 > vertexNeighborAEGC_
virtual int getTetrahedronTriangle(const SimplexId &tetId, const int &id, SimplexId &triangleId) const =0
std::array< SimplexId, 10 > vertexNeighborEFGH_
std::array< SimplexId, 2 > vertexNeighbor2dA_
bool isPowerOfTwo(unsigned long long int v, unsigned long long int &r)
SimplexId TTK_TRIANGULATION_INTERNAL() getNumberOfVertices() const override
SimplexId getVertexEdgeNumberInternal(const SimplexId &vertexId) const override
std::array< SimplexId, 4 > vertexNeighbor2dBD_
std::array< SimplexId, 10 > vertexNeighborBFHD_
int getTriangleNeighbors(std::vector< std::vector< SimplexId > > &neighbors)
std::array< SimplexId, 4 > vertexNeighborD_
SimplexId getNumberOfEdgesInternal() const override
std::array< SimplexId, 4 > vertexNeighborC_
virtual int getTetrahedronEdge(const SimplexId &tetId, const int &id, SimplexId &edgeId) const =0
virtual SimplexId getTriangleNeighborNumber(const SimplexId &triangleId) const =0
int getTetrahedronEdges(std::vector< std::vector< SimplexId > > &edges) const
std::array< SimplexId, 4 > vertexNeighbor2dCD_
std::array< SimplexId, 8 > vertexNeighborCG_
std::array< SimplexId, 3 > vertexNeighbor2dC_
std::array< SimplexId, 6 > vertexNeighborAC_
const std::vector< std::vector< SimplexId > > * getEdgeTrianglesInternal() override
std::array< SimplexId, 10 > vertexNeighborABCD_
std::array< SimplexId, 10 > vertexNeighborGHDC_
std::array< SimplexId, 8 > vertexNeighborAB_
std::array< SimplexId, 6 > vertexNeighborAE_
int getCellEdgeInternal(const SimplexId &cellId, const int &id, SimplexId &edgeId) const override
std::array< SimplexId, 4 > vertexNeighborA_
SimplexId getNumberOfTrianglesInternal() const override
std::array< SimplexId, 4 > vertexNeighborF_
const std::vector< std::vector< SimplexId > > * getCellTrianglesInternal() override
std::array< SimplexId, 4 > vertexNeighbor2dAB_
int getTetrahedronTriangles(std::vector< std::vector< SimplexId > > &triangles) const
std::array< SimplexId, 3 > dimensions_
double getElapsedTime()
Definition Timer.h:15
std::string to_string(__int128)
Definition ripserpy.cpp:99
The Topology ToolKit.
COMMON_EXPORTS int MPIsize_
Definition BaseClass.cpp:10
COMMON_EXPORTS int MPIrank_
Definition BaseClass.cpp:9
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)