TTK
Loading...
Searching...
No Matches
DiscreteGradient_Template.h
Go to the documentation of this file.
1
14
15#pragma once
16
17#include <DiscreteGradient.h>
18
19using ttk::SimplexId;
20using ttk::dcg::Cell;
23
24template <typename dataType, typename triangulationType>
26 const Cell &up,
27 const Cell &down,
28 const dataType *const scalars,
29 const triangulationType &triangulation) const {
30
31 return scalars[getCellGreaterVertex(up, triangulation)]
32 - scalars[getCellLowerVertex(down, triangulation)];
33}
34
35template <typename triangulationType>
36int DiscreteGradient::buildGradient(const triangulationType &triangulation,
37 bool bypassCache) {
38
39 auto &cacheHandler = *triangulation.getGradientCacheHandler();
40 const auto findGradient
41 = [this, &cacheHandler]() -> AbstractTriangulation::gradientType * {
42 if(this->inputScalarField_.first == nullptr) {
43 return {};
44 }
45 return cacheHandler.get(this->inputScalarField_);
46 };
47
48#ifdef TTK_ENABLE_OPENMP
49 if(!bypassCache && omp_in_parallel()) {
50 this->printWrn(
51 "buildGradient() called inside a parallel region, disabling cache...");
52 bypassCache = true;
53 }
54#endif // TTK_ENABLE_OPENMP
55
56 // set member variables at each buildGradient() call
57 this->dimensionality_ = triangulation.getCellVertexNumber(0) - 1;
58 this->numberOfVertices_ = triangulation.getNumberOfVertices();
59
60 this->gradient_ = bypassCache ? &this->localGradient_ : findGradient();
61 if(this->gradient_ == nullptr || bypassCache) {
62
63 if(!bypassCache) {
64 // add new cache entry
65 cacheHandler.insert(this->inputScalarField_, {});
66 this->gradient_ = cacheHandler.get(this->inputScalarField_);
67 }
68
69 // allocate gradient memory
70 this->initMemory(triangulation);
71
72 Timer tm{};
73 // compute gradient pairs
74 this->processLowerStars(this->inputOffsets_, triangulation);
75
76 this->printMsg(
77 "Built discrete gradient", 1.0, tm.getElapsedTime(), this->threadNumber_);
78 } else {
79 this->printMsg("Fetched cached discrete gradient");
80 }
81
82 return 0;
83}
84
85template <typename triangulationType>
87 const std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
88 std::vector<std::array<float, 3>> &points,
89 std::vector<char> &cellDimensions,
90 std::vector<SimplexId> &cellIds,
91 std::vector<char> &isOnBoundary,
92 std::vector<SimplexId> &PLVertexIdentifiers,
93 const triangulationType &triangulation) const {
94
95 std::array<size_t, 5> partSums{};
96 for(size_t i = 0; i < criticalCellsByDim.size(); ++i) {
97 partSums[i + 1] = partSums[i] + criticalCellsByDim[i].size();
98 }
99
100 const auto nCritPoints = partSums.back();
101
102 points.resize(nCritPoints);
103 cellDimensions.resize(nCritPoints);
104 cellIds.resize(nCritPoints);
105 isOnBoundary.resize(nCritPoints);
106 PLVertexIdentifiers.resize(nCritPoints);
107
108 for(size_t i = 0; i < criticalCellsByDim.size(); ++i) {
109#ifdef TTK_ENABLE_OPENMP
110#pragma omp parallel for num_threads(threadNumber_)
111#endif // TTK_ENABLE_OPENMP
112 for(size_t j = 0; j < criticalCellsByDim[i].size(); ++j) {
113 const SimplexId cellId = criticalCellsByDim[i][j];
114 const int cellDim = i;
115 const auto o{partSums[i] + j};
116
117 triangulation.getCellIncenter(cellId, i, points[o].data());
118 cellDimensions[o] = cellDim;
119#ifdef TTK_ENABLE_MPI
120 ttk::SimplexId globalId{-1};
121 triangulation.getDistributedGlobalCellId(cellId, cellDim, globalId);
122 cellIds[o] = globalId;
123#else
124 cellIds[o] = cellId;
125#endif // TTK_ENABLE_MPI
126 const Cell cell{static_cast<int>(i), cellId};
127 isOnBoundary[o] = this->isBoundary(cell, triangulation);
128 PLVertexIdentifiers[o] = this->getCellGreaterVertex(cell, triangulation);
129 }
130 }
131
132 std::vector<std::vector<std::string>> rows(this->dimensionality_ + 1);
133 for(int i = 0; i < this->dimensionality_ + 1; ++i) {
134 rows[i]
135 = std::vector<std::string>{"#" + std::to_string(i) + "-cell(s)",
136 std::to_string(criticalCellsByDim[i].size())};
137 }
138 this->printMsg(rows);
139
140 return 0;
141}
142
143template <typename triangulationType>
145 std::vector<std::array<float, 3>> &points,
146 std::vector<char> &cellDimensions,
147 std::vector<SimplexId> &cellIds,
148 std::vector<char> &isOnBoundary,
149 std::vector<SimplexId> &PLVertexIdentifiers,
150 const triangulationType &triangulation) const {
151
152 std::array<std::vector<SimplexId>, 4> criticalCellsByDim;
153 getCriticalPoints(criticalCellsByDim, triangulation);
154 setCriticalPoints(criticalCellsByDim, points, cellDimensions, cellIds,
155 isOnBoundary, PLVertexIdentifiers, triangulation);
156
157 return 0;
158}
159
160template <typename triangulationType>
162 std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
163 const triangulationType &triangulation) const {
164
165 const auto dims{this->getNumberOfDimensions()};
166 for(int i = 0; i < dims; ++i) {
167
168 // map: store critical cell per dimension per thread
169 std::vector<std::vector<SimplexId>> critCellsPerThread(this->threadNumber_);
170 const auto numberOfCells{this->getNumberOfCells(i, triangulation)};
171
172 // use static scheduling to ensure that critical cells
173 // are sorted by id
174
175#ifdef TTK_ENABLE_OPENMP
176#pragma omp parallel for num_threads(this->threadNumber_) schedule(static)
177#endif // TTK_ENABLE_OPENMP
178 for(SimplexId j = 0; j < numberOfCells; ++j) {
179#ifdef TTK_ENABLE_OPENMP
180 const auto tid = omp_get_thread_num();
181#else
182 const auto tid = 0;
183#endif // TTK_ENABLE_OPENMP
184 if(this->isCellCritical(i, j)) {
185 critCellsPerThread[tid].emplace_back(j);
186 }
187 }
188
189 // reduce: aggregate critical cells per thread
190 criticalCellsByDim[i] = std::move(critCellsPerThread[0]);
191 for(size_t j = 1; j < critCellsPerThread.size(); ++j) {
192 const auto &vec{critCellsPerThread[j]};
193 criticalCellsByDim[i].insert(
194 criticalCellsByDim[i].end(), vec.begin(), vec.end());
195 }
196 }
197
198 return 0;
199}
200
201template <typename triangulationType>
203 const int dimension, const triangulationType &triangulation) const {
204
205 if(dimension > this->dimensionality_ || dimension < 0) {
206 return -1;
207 }
208
209 switch(dimension) {
210 case 0:
211 return triangulation.getNumberOfVertices();
212 break;
213
214 case 1:
215 return triangulation.getNumberOfEdges();
216 break;
217
218 case 2:
219 return triangulation.getNumberOfTriangles();
220 break;
221
222 case 3:
223 return triangulation.getNumberOfCells();
224 break;
225 }
226
227 return -1;
228}
229
230template <typename triangulationType>
231inline void
232 DiscreteGradient::lowerStar(lowerStarType &ls,
233 const SimplexId a,
234 const SimplexId *const offsets,
235 const triangulationType &triangulation) const {
236
237 // make sure that ls is cleared
238 for(auto &vec : ls) {
239 vec.clear();
240 }
241
242 // a belongs to its lower star
243 ls[0].emplace_back(CellExt{0, a});
244
245 // store lower edges
246 const auto nedges = triangulation.getVertexEdgeNumber(a);
247 ls[1].reserve(nedges);
248 for(SimplexId i = 0; i < nedges; i++) {
249 SimplexId edgeId;
250 triangulation.getVertexEdge(a, i, edgeId);
251 SimplexId vertexId;
252 triangulation.getEdgeVertex(edgeId, 0, vertexId);
253 if(vertexId == a) {
254 triangulation.getEdgeVertex(edgeId, 1, vertexId);
255 }
256 if(offsets[vertexId] < offsets[a]) {
257 ls[1].emplace_back(CellExt{1, edgeId, {offsets[vertexId], -1, -1}, {}});
258 }
259 }
260
261 if(ls[1].size() < 2) {
262 // at least two edges in the lower star for one triangle
263 return;
264 }
265
266 const auto processTriangle
267 = [&](const SimplexId triangleId, const SimplexId v0, const SimplexId v1,
268 const SimplexId v2) {
269 std::array<SimplexId, 3> lowVerts{-1, -1, -1};
270 if(v0 == a) {
271 lowVerts[0] = offsets[v1];
272 lowVerts[1] = offsets[v2];
273 } else if(v1 == a) {
274 lowVerts[0] = offsets[v0];
275 lowVerts[1] = offsets[v2];
276 } else if(v2 == a) {
277 lowVerts[0] = offsets[v0];
278 lowVerts[1] = offsets[v1];
279 }
280 // higher order vertex first
281 if(lowVerts[0] < lowVerts[1]) {
282 std::swap(lowVerts[0], lowVerts[1]);
283 }
284 if(offsets[a] > lowVerts[0]) { // triangle in lowerStar
285 uint8_t j{}, k{};
286 // store edges indices of current triangle
287 std::array<uint8_t, 3> faces{};
288 for(const auto &e : ls[1]) {
289 if(e.lowVerts_[0] == lowVerts[0] || e.lowVerts_[0] == lowVerts[1]) {
290 faces[k++] = j;
291 }
292 j++;
293 }
294 ls[2].emplace_back(CellExt{2, triangleId, lowVerts, faces});
295 }
296 };
297
298 if(dimensionality_ == 2) {
299 // store lower triangles
300
301 // use optimised triangulation methods:
302 // getVertexStar instead of getVertexTriangle
303 // getCellVertex instead of getTriangleVertex
304 const auto ncells = triangulation.getVertexStarNumber(a);
305 ls[2].reserve(ncells);
306 for(SimplexId i = 0; i < ncells; ++i) {
307 SimplexId cellId;
308 triangulation.getVertexStar(a, i, cellId);
309 SimplexId v0{}, v1{}, v2{};
310 triangulation.getCellVertex(cellId, 0, v0);
311 triangulation.getCellVertex(cellId, 1, v1);
312 triangulation.getCellVertex(cellId, 2, v2);
313 processTriangle(cellId, v0, v1, v2);
314 }
315 } else if(dimensionality_ == 3) {
316 // store lower triangles
317 const auto ntri = triangulation.getVertexTriangleNumber(a);
318 ls[2].reserve(ntri);
319 for(SimplexId i = 0; i < ntri; i++) {
320 SimplexId triangleId;
321 triangulation.getVertexTriangle(a, i, triangleId);
322 SimplexId v0{}, v1{}, v2{};
323 triangulation.getTriangleVertex(triangleId, 0, v0);
324 triangulation.getTriangleVertex(triangleId, 1, v1);
325 triangulation.getTriangleVertex(triangleId, 2, v2);
326 processTriangle(triangleId, v0, v1, v2);
327 }
328
329 // at least three triangles in the lower star for one tetra
330 if(ls[2].size() >= 3) {
331 // store lower tetra
332 const auto ncells = triangulation.getVertexStarNumber(a);
333 ls[3].reserve(ncells);
334 for(SimplexId i = 0; i < ncells; ++i) {
335 SimplexId cellId;
336 triangulation.getVertexStar(a, i, cellId);
337 std::array<SimplexId, 3> lowVerts{-1, -1, -1};
338 SimplexId v0{}, v1{}, v2{}, v3{};
339 triangulation.getCellVertex(cellId, 0, v0);
340 triangulation.getCellVertex(cellId, 1, v1);
341 triangulation.getCellVertex(cellId, 2, v2);
342 triangulation.getCellVertex(cellId, 3, v3);
343 if(v0 == a) {
344 lowVerts[0] = offsets[v1];
345 lowVerts[1] = offsets[v2];
346 lowVerts[2] = offsets[v3];
347 } else if(v1 == a) {
348 lowVerts[0] = offsets[v0];
349 lowVerts[1] = offsets[v2];
350 lowVerts[2] = offsets[v3];
351 } else if(v2 == a) {
352 lowVerts[0] = offsets[v0];
353 lowVerts[1] = offsets[v1];
354 lowVerts[2] = offsets[v3];
355 } else if(v3 == a) {
356 lowVerts[0] = offsets[v0];
357 lowVerts[1] = offsets[v1];
358 lowVerts[2] = offsets[v2];
359 }
360 if(offsets[a] > *std::max_element(
361 lowVerts.begin(), lowVerts.end())) { // tetra in lowerStar
362
363 // higher order vertex first
364 std::sort(lowVerts.rbegin(), lowVerts.rend());
365
366 uint8_t j{}, k{};
367 // store triangles indices of current tetra
368 std::array<uint8_t, 3> faces{};
369 for(const auto &t : ls[2]) {
370 // lowVerts & t.lowVerts are ordered, no need to check if
371 // t.lowVerts[0] == lowVerts[2] or t.lowVerts[1] == lowVerts[0]
372 if((t.lowVerts_[0] == lowVerts[0]
373 && (t.lowVerts_[1] == lowVerts[1]
374 || t.lowVerts_[1] == lowVerts[2]))
375 || (t.lowVerts_[0] == lowVerts[1]
376 && t.lowVerts_[1] == lowVerts[2])) {
377 faces[k++] = j;
378 }
379 j++;
380 }
381
382 ls[3].emplace_back(CellExt{3, cellId, lowVerts, faces});
383 }
384 }
385 }
386 }
387}
388
389template <typename triangulationType>
390inline void DiscreteGradient::pairCells(
391 CellExt &alpha, CellExt &beta, const triangulationType &triangulation) {
392#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
393 char localBId{0}, localAId{0};
394 SimplexId a{}, b{};
395
396 if(beta.dim_ == 1) {
397
398 for(SimplexId i = 0; i < 2; ++i) {
399 triangulation.getEdgeVertex(beta.id_, i, a);
400 if(a == alpha.id_) {
401 localAId = i;
402 break;
403 }
404 }
405 const auto nedges = triangulation.getVertexEdgeNumber(alpha.id_);
406 for(SimplexId i = 0; i < nedges; ++i) {
407 triangulation.getVertexEdge(alpha.id_, i, b);
408 if(b == beta.id_) {
409 localBId = i;
410 break;
411 }
412 }
413 } else if(beta.dim_ == 2) {
414 for(SimplexId i = 0; i < 3; ++i) {
415 triangulation.getTriangleEdge(beta.id_, i, a);
416 if(a == alpha.id_) {
417 localAId = i;
418 break;
419 }
420 }
421 const auto ntri = triangulation.getEdgeTriangleNumber(alpha.id_);
422 for(SimplexId i = 0; i < ntri; ++i) {
423 triangulation.getEdgeTriangle(alpha.id_, i, b);
424 if(b == beta.id_) {
425 localBId = i;
426 break;
427 }
428 }
429 } else {
430 for(SimplexId i = 0; i < 4; ++i) {
431 triangulation.getCellTriangle(beta.id_, i, a);
432 if(a == alpha.id_) {
433 localAId = i;
434 break;
435 }
436 }
437 const auto ntetra = triangulation.getTriangleStarNumber(alpha.id_);
438 for(SimplexId i = 0; i < ntetra; ++i) {
439 triangulation.getTriangleStar(alpha.id_, i, b);
440 if(b == beta.id_) {
441 localBId = i;
442 break;
443 }
444 }
445 }
446 (*gradient_)[2 * alpha.dim_][alpha.id_] = localBId;
447 (*gradient_)[2 * alpha.dim_ + 1][beta.id_] = localAId;
448#else
449 TTK_FORCE_USE(triangulation);
450 (*gradient_)[2 * alpha.dim_][alpha.id_] = beta.id_;
451 (*gradient_)[2 * alpha.dim_ + 1][beta.id_] = alpha.id_;
452#endif // TTK_ENABLE_DCG_OPTIMIZE_MEMORY
453 alpha.paired_ = true;
454 beta.paired_ = true;
455}
456
457template <typename triangulationType>
458int DiscreteGradient::processLowerStars(
459 const SimplexId *const offsets, const triangulationType &triangulation) {
460
461 /* Compute gradient */
462
463 auto nverts = triangulation.getNumberOfVertices();
464
465 // Comparison function for Cells inside priority queues
466 const auto orderCells = [&](const CellExt &a, const CellExt &b) -> bool {
467 return a.lowVerts_ > b.lowVerts_;
468 };
469
470 // Type alias for priority queues
471 using pqType
472 = std::priority_queue<std::reference_wrapper<CellExt>,
473 std::vector<std::reference_wrapper<CellExt>>,
474 decltype(orderCells)>;
475
476 // To reduce allocations, priority queues and lowerStar objects are
477 // cleaned & reused between iterations.
478
479 // Priority queues are pushed at the beginning and popped at the
480 // end. To pop the minimum, elements should be sorted in a
481 // decreasing order.
482 pqType pqZero{orderCells}, pqOne{orderCells};
483
484 // store lower star structure
485 lowerStarType Lx;
486
487#ifdef TTK_ENABLE_OPENMP
488#pragma omp parallel for num_threads(threadNumber_) \
489 firstprivate(Lx, pqZero, pqOne)
490#endif // TTK_ENABLE_OPENMP
491 for(SimplexId x = 0; x < nverts; x++) {
492
493 // clear priority queues (they should be empty at the end of the
494 // previous iteration)
495 while(!pqZero.empty()) {
496 pqZero.pop();
497 }
498 while(!pqOne.empty()) {
499 pqOne.pop();
500 }
501
502 // Insert into pqOne cofacets of cell c_alpha such as numUnpairedFaces == 1
503 const auto insertCofacets = [&](const CellExt &ca, lowerStarType &ls) {
504 if(ca.dim_ == 1) {
505 for(auto &beta : ls[2]) {
506 if(ls[1][beta.faces_[0]].id_ == ca.id_
507 || ls[1][beta.faces_[1]].id_ == ca.id_) {
508 // edge ca belongs to triangle beta
509 if(numUnpairedFacesTriangle(beta, ls).first == 1) {
510 pqOne.push(beta);
511 }
512 }
513 }
514
515 } else if(ca.dim_ == 2) {
516 for(auto &beta : ls[3]) {
517 if(ls[2][beta.faces_[0]].id_ == ca.id_
518 || ls[2][beta.faces_[1]].id_ == ca.id_
519 || ls[2][beta.faces_[2]].id_ == ca.id_) {
520 // triangle ca belongs to tetra beta
521 if(numUnpairedFacesTetra(beta, ls).first == 1) {
522 pqOne.push(beta);
523 }
524 }
525 }
526 }
527 };
528
529 lowerStar(Lx, x, offsets, triangulation);
530 // In case the vertex is a ghost, the gradient of the
531 // simplices of its star is set to GHOST_GRADIENT
532#ifdef TTK_ENABLE_MPI
533 if(ttk::isRunningWithMPI()
534 && triangulation.getVertexRank(x) != ttk::MPIrank_) {
535 int sizeDim = Lx.size();
536 for(int i = 0; i < sizeDim; i++) {
537 int nCells = Lx[i].size();
538 for(int j = 0; j < nCells; j++) {
539 setCellToGhost(Lx[i][j].dim_, Lx[i][j].id_);
540 }
541 }
542 } else
543#endif // TTK_ENABLE_MPI
544
545 {
546 // Lx[1] empty => x is a local minimum
547 if(!Lx[1].empty()) {
548 // get delta: 1-cell (edge) with minimal G value (steeper gradient)
549 size_t minId = 0;
550 for(size_t i = 1; i < Lx[1].size(); ++i) {
551 const auto &a = Lx[1][minId].lowVerts_[0];
552 const auto &b = Lx[1][i].lowVerts_[0];
553 if(a > b) {
554 // edge[i] < edge[0]
555 minId = i;
556 }
557 }
558
559 auto &c_delta = Lx[1][minId];
560
561 // store x (0-cell) -> delta (1-cell) V-path
562 pairCells(Lx[0][0], c_delta, triangulation);
563
564 // push every 1-cell in Lx that is not delta into pqZero
565 for(auto &alpha : Lx[1]) {
566 if(alpha.id_ != c_delta.id_) {
567 pqZero.push(alpha);
568 }
569 }
570
571 // push into pqOne every coface of delta in Lx (2-cells only,
572 // 3-cells have not any facet paired yet) such that
573 // numUnpairedFaces == 1
574 insertCofacets(c_delta, Lx);
575
576 while(!pqOne.empty() || !pqZero.empty()) {
577 while(!pqOne.empty()) {
578 auto &c_alpha = pqOne.top().get();
579 pqOne.pop();
580 auto unpairedFaces = numUnpairedFaces(c_alpha, Lx);
581 if(unpairedFaces.first == 0) {
582 pqZero.push(c_alpha);
583 } else {
584 auto &c_pair_alpha = Lx[c_alpha.dim_ - 1][unpairedFaces.second];
585
586 // store (pair_alpha) -> (alpha) V-path
587 pairCells(c_pair_alpha, c_alpha, triangulation);
588
589 // add cofaces of c_alpha and c_pair_alpha to pqOne
590 insertCofacets(c_alpha, Lx);
591 insertCofacets(c_pair_alpha, Lx);
592 }
593 }
594
595 // skip pair_alpha from pqZero:
596 // cells in pqZero are not critical if already paired
597 while(!pqZero.empty() && pqZero.top().get().paired_) {
598 pqZero.pop();
599 }
600
601 if(!pqZero.empty()) {
602 auto &c_gamma = pqZero.top().get();
603 pqZero.pop();
604
605 // gamma is a critical cell
606 // mark gamma as paired
607 c_gamma.paired_ = true;
608
609 // add cofacets of c_gamma to pqOne
610 insertCofacets(c_gamma, Lx);
611 }
612 }
613 }
614 }
615 }
616
617 return 0;
618}
619
620template <typename triangulationType>
622 const Cell &cell, const triangulationType &triangulation) const {
623
624 if(cell.dim_ > this->dimensionality_ || cell.dim_ < 0) {
625 return false;
626 }
627
628 const auto vert{this->getCellGreaterVertex(cell, triangulation)};
629 return triangulation.isVertexOnBoundary(vert);
630}
631
632template <typename triangulationType>
635 const triangulationType &triangulation,
636 bool isReverse) const {
637
638 // ensure that getPairedCell(Cell, boolean) calls are rejected
639 static_assert(
640 std::is_base_of<AbstractTriangulation, triangulationType>(),
641 "triangulationType should be an AbstractTriangulation derivative");
642
643#ifndef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
644 TTK_FORCE_USE(triangulation);
645#endif // !TTK_ENABLE_DCG_OPTIMIZE_MEMORY
646
647 if((cell.dim_ > this->dimensionality_ - 1 && !isReverse)
648 || (cell.dim_ > this->dimensionality_ && isReverse) || cell.dim_ < 0) {
649 return -1;
650 }
651
652 SimplexId id{-1};
653
654 if(cell.dim_ == 0) {
655 if(!isReverse) {
656#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
657 const auto locId{(*gradient_)[0][cell.id_]};
658 if(locId != -1) {
659 triangulation.getVertexEdge(cell.id_, locId, id);
660 }
661#else
662 id = (*gradient_)[0][cell.id_];
663#endif
664 }
665 }
666
667 else if(cell.dim_ == 1) {
668 if(isReverse) {
669#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
670 const auto locId{(*gradient_)[1][cell.id_]};
671 if(locId != -1) {
672 triangulation.getEdgeVertex(cell.id_, locId, id);
673 }
674#else
675 id = (*gradient_)[1][cell.id_];
676#endif
677 } else {
678#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
679 const auto locId{(*gradient_)[2][cell.id_]};
680 if(locId != -1) {
681 triangulation.getEdgeTriangle(cell.id_, locId, id);
682 }
683#else
684 id = (*gradient_)[2][cell.id_];
685#endif
686 }
687 }
688
689 else if(cell.dim_ == 2) {
690 if(isReverse) {
691#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
692 const auto locId{(*gradient_)[3][cell.id_]};
693 if(locId != -1) {
694 triangulation.getTriangleEdge(cell.id_, locId, id);
695 }
696#else
697 id = (*gradient_)[3][cell.id_];
698#endif
699 } else {
700#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
701 const auto locId{(*gradient_)[4][cell.id_]};
702 if(locId != -1) {
703 triangulation.getTriangleStar(cell.id_, locId, id);
704 }
705#else
706 id = (*gradient_)[4][cell.id_];
707#endif
708 }
709 }
710
711 else if(cell.dim_ == 3) {
712 if(isReverse) {
713#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
714 const auto locId{(*gradient_)[5][cell.id_]};
715 if(locId != -1) {
716 triangulation.getCellTriangle(cell.id_, locId, id);
717 }
718#else
719 id = (*gradient_)[5][cell.id_];
720#endif
721 }
722 }
723
724 return id;
725}
726
727template <typename triangulationType>
729 const Cell &cell,
730 std::vector<Cell> &vpath,
731 const triangulationType &triangulation) const {
732
733 if(cell.dim_ == 0) {
734 // assume that cellId is a vertex
735 SimplexId currentId = cell.id_;
736 SimplexId connectedEdgeId;
737 do {
738 // add a vertex
739 const Cell vertex(0, currentId);
740 vpath.push_back(vertex);
741
742 if(isCellCritical(vertex)) {
743 break;
744 }
745
746 connectedEdgeId = getPairedCell(vertex, triangulation);
747 if(connectedEdgeId == -1) {
748 break;
749 }
750
751 // add an edge
752 const Cell edge(1, connectedEdgeId);
753 vpath.push_back(edge);
754
755 if(isCellCritical(edge)) {
756 break;
757 }
758
759 for(int i = 0; i < 2; ++i) {
760 SimplexId vertexId;
761 triangulation.getEdgeVertex(connectedEdgeId, i, vertexId);
762
763 if(vertexId != currentId) {
764 currentId = vertexId;
765 break;
766 }
767 }
768
769 } while(connectedEdgeId != -1);
770 }
771
772 return 0;
773}
774
775template <typename triangulationType>
777 const Cell &saddle2,
778 const Cell &saddle1,
779 const std::vector<bool> &isVisited,
780 std::vector<Cell> *const vpath,
781 const triangulationType &triangulation,
782 const bool stopIfMultiConnected,
783 const bool enableCycleDetector) const {
784
785 // debug
786 const SimplexId numberOfEdges = triangulation.getNumberOfEdges();
787 std::vector<char> isCycle;
788 if(enableCycleDetector) {
789 isCycle.resize(numberOfEdges, 0);
790 }
791
792 if(dimensionality_ == 3) {
793 // add the 2-saddle to the path
794 if(vpath != nullptr) {
795 vpath->push_back(saddle2);
796 }
797
798 SimplexId currentId = -1;
799 {
800 int nconnections = 0;
801 for(int i = 0; i < 3; ++i) {
802 SimplexId edgeId;
803 triangulation.getTriangleEdge(saddle2.id_, i, edgeId);
804 if(isVisited[edgeId]) {
805 // saddle2 can be adjacent to saddle1 on the wall
806 if(isSaddle1(Cell(1, edgeId))) {
807 if(vpath != nullptr) {
808 vpath->push_back(Cell(1, edgeId));
809 }
810 return false;
811 }
812
813 currentId = edgeId;
814 ++nconnections;
815 }
816 }
817 if(stopIfMultiConnected && nconnections > 1) {
818 return true;
819 }
820 }
821
822 int oldId;
823 do {
824
825 // debug
826 if(enableCycleDetector) {
827 if(isCycle[currentId] == 0) {
828 isCycle[currentId] = 1;
829 } else {
830 this->printErr("Cycle detected on the wall of 1-saddle "
831 + std::to_string(saddle1.id_));
832 break;
833 }
834 }
835
836 oldId = currentId;
837
838 // add an edge
839 const Cell edge(1, currentId);
840 if(vpath != nullptr) {
841 vpath->push_back(edge);
842 }
843
844 if(isCellCritical(edge)) {
845 break;
846 }
847
848 const SimplexId connectedTriangleId = getPairedCell(edge, triangulation);
849
850 // add a triangle
851 const Cell triangle(2, connectedTriangleId);
852 if(vpath != nullptr) {
853 vpath->push_back(triangle);
854 }
855
856 if(isCellCritical(triangle)) {
857 break;
858 }
859
860 int nconnections = 0;
861 for(int i = 0; i < 3; ++i) {
862 SimplexId edgeId;
863 triangulation.getTriangleEdge(connectedTriangleId, i, edgeId);
864
865 if(isVisited[edgeId] and edgeId != oldId) {
866 currentId = edgeId;
867 ++nconnections;
868 }
869 }
870 if(stopIfMultiConnected && nconnections > 1) {
871 return true;
872 }
873
874 // stop at convergence caused by boundary effect
875 } while(currentId != oldId);
876 }
877
878 return false;
879}
880
881template <typename triangulationType>
883 std::vector<Cell> &vpath,
884 const triangulationType &triangulation,
885 const bool enableCycleDetector) const {
886
887 const SimplexId numberOfCells = triangulation.getNumberOfCells();
888 std::vector<char> isCycle;
889 if(enableCycleDetector) {
890 isCycle.resize(numberOfCells, 0);
891 }
892
893 if(dimensionality_ == 2) {
894 if(cell.dim_ == 2) {
895 // assume that cellId is a triangle
896 SimplexId currentId = cell.id_;
897 SimplexId oldId;
898 do {
899 oldId = currentId;
900
901 // add a triangle
902 const Cell triangle(2, currentId);
903 vpath.push_back(triangle);
904
905 if(isCellCritical(triangle)) {
906 break;
907 }
908
909 const SimplexId connectedEdgeId
910 = getPairedCell(triangle, triangulation, true);
911 if(connectedEdgeId == -1) {
912 break;
913 }
914
915 // add an edge
916 const Cell edge(1, connectedEdgeId);
917 vpath.push_back(edge);
918
919 if(isCellCritical(edge)) {
920 break;
921 }
922
923 const SimplexId starNumber
924 = triangulation.getEdgeStarNumber(connectedEdgeId);
925 for(SimplexId i = 0; i < starNumber; ++i) {
926 SimplexId starId;
927 triangulation.getEdgeStar(connectedEdgeId, i, starId);
928
929 if(starId != currentId) {
930 currentId = starId;
931 break;
932 }
933 }
934
935 // stop at convergence caused by boundary effect
936 } while(currentId != oldId);
937 }
938 } else if(dimensionality_ == 3) {
939 if(cell.dim_ == 3) {
940 // assume that cellId is a tetra
941 SimplexId currentId = cell.id_;
942 SimplexId oldId;
943 do {
944
945 // debug
946 if(enableCycleDetector) {
947 if(isCycle[currentId] == 0) {
948 isCycle[currentId] = 1;
949 } else {
950 this->printErr("cycle detected in the path from tetra "
951 + std::to_string(cell.id_));
952 break;
953 }
954 }
955
956 oldId = currentId;
957
958 // add a tetra
959 const Cell tetra(3, currentId);
960 vpath.push_back(tetra);
961
962 if(isCellCritical(tetra)) {
963 break;
964 }
965
966 const SimplexId connectedTriangleId
967 = getPairedCell(tetra, triangulation, true);
968 if(connectedTriangleId == -1) {
969 break;
970 }
971
972 // add a triangle
973 const Cell triangle(2, connectedTriangleId);
974 vpath.push_back(triangle);
975
976 if(isCellCritical(triangle)) {
977 break;
978 }
979
980 const SimplexId starNumber
981 = triangulation.getTriangleStarNumber(connectedTriangleId);
982 for(SimplexId i = 0; i < starNumber; ++i) {
983 SimplexId starId;
984 triangulation.getTriangleStar(connectedTriangleId, i, starId);
985
986 if(starId != currentId) {
987 currentId = starId;
988 break;
989 }
990 }
991
992 // stop at convergence caused by boundary effect
993 } while(currentId != oldId);
994 }
995 }
996
997 return 0;
998}
999
1000template <typename triangulationType>
1002 const Cell &saddle1,
1003 const Cell &saddle2,
1004 const std::vector<bool> &isVisited,
1005 std::vector<Cell> *const vpath,
1006 const triangulationType &triangulation,
1007 const bool stopIfMultiConnected,
1008 const bool enableCycleDetector) const {
1009
1010 // debug
1011 const SimplexId numberOfTriangles = triangulation.getNumberOfTriangles();
1012 std::vector<char> isCycle;
1013 if(enableCycleDetector) {
1014 isCycle.resize(numberOfTriangles, 0);
1015 }
1016
1017 if(dimensionality_ == 3) {
1018 // add the 1-saddle to the path
1019 if(vpath != nullptr) {
1020 vpath->push_back(saddle1);
1021 }
1022
1023 SimplexId currentId = -1;
1024 {
1025 int nconnections = 0;
1026 const SimplexId triangleNumber
1027 = triangulation.getEdgeTriangleNumber(saddle1.id_);
1028 for(SimplexId i = 0; i < triangleNumber; ++i) {
1029 SimplexId triangleId;
1030 triangulation.getEdgeTriangle(saddle1.id_, i, triangleId);
1031 if(isVisited[triangleId]) {
1032 // saddle1 can be adjacent to saddle2 on the wall
1033 if(isSaddle2(Cell(2, triangleId))) {
1034 if(vpath != nullptr) {
1035 vpath->push_back(Cell(2, triangleId));
1036 }
1037 return false;
1038 }
1039
1040 currentId = triangleId;
1041 ++nconnections;
1042 }
1043 }
1044 if(stopIfMultiConnected && nconnections > 1) {
1045 return true;
1046 }
1047 }
1048
1049 if(currentId == -1) {
1050 return true;
1051 }
1052
1053 SimplexId oldId;
1054 do {
1055
1056 // debug
1057 if(enableCycleDetector) {
1058 if(isCycle[currentId] == 0) {
1059 isCycle[currentId] = 1;
1060 } else {
1061 this->printErr("Cycle detected on the wall of 2-saddle "
1062 + std::to_string(saddle2.id_));
1063 break;
1064 }
1065 }
1066
1067 oldId = currentId;
1068
1069 // add a triangle
1070 const Cell triangle(2, currentId);
1071 if(vpath != nullptr) {
1072 vpath->push_back(triangle);
1073 }
1074
1075 if(isCellCritical(triangle)) {
1076 break;
1077 }
1078
1079 const SimplexId connectedEdgeId
1080 = getPairedCell(triangle, triangulation, true);
1081
1082 // add an edge
1083 const Cell edge(1, connectedEdgeId);
1084 if(vpath != nullptr) {
1085 vpath->push_back(edge);
1086 }
1087
1088 if(isCellCritical(edge)) {
1089 break;
1090 }
1091
1092 int nconnections = 0;
1093 const SimplexId triangleNumber
1094 = triangulation.getEdgeTriangleNumber(connectedEdgeId);
1095 for(SimplexId i = 0; i < triangleNumber; ++i) {
1096 SimplexId triangleId;
1097 triangulation.getEdgeTriangle(connectedEdgeId, i, triangleId);
1098
1099 if(isVisited[triangleId] and triangleId != oldId) {
1100 currentId = triangleId;
1101 ++nconnections;
1102 }
1103 }
1104 if(stopIfMultiConnected && nconnections > 1) {
1105 return true;
1106 }
1107
1108 // stop at convergence caused by boundary effect
1109 } while(currentId != oldId);
1110 }
1111
1112 return false;
1113}
1114
1115template <typename triangulationType>
1117 const Cell &cell,
1118 VisitedMask &mask,
1119 const triangulationType &triangulation,
1120 std::vector<Cell> *const wall,
1121 std::vector<SimplexId> *const saddles) const {
1122
1123 if(saddles != nullptr) {
1124 saddles->clear();
1125 }
1126
1127 if(dimensionality_ == 3) {
1128 if(cell.dim_ == 2) {
1129 // assume that cellId is a triangle
1130 const SimplexId originId = cell.id_;
1131
1132 std::queue<SimplexId> bfs;
1133 bfs.push(originId);
1134
1135 // BFS traversal
1136 while(!bfs.empty()) {
1137 const SimplexId triangleId = bfs.front();
1138 bfs.pop();
1139
1140 if(!mask.isVisited_[triangleId]) {
1141 mask.isVisited_[triangleId] = true;
1142 mask.visitedIds_.emplace_back(triangleId);
1143
1144 // add the triangle
1145 if(wall != nullptr) {
1146 wall->push_back(Cell(2, triangleId));
1147 }
1148
1149 for(int j = 0; j < 3; ++j) {
1150 SimplexId edgeId;
1151 triangulation.getTriangleEdge(triangleId, j, edgeId);
1152
1153 if((saddles != nullptr) and isSaddle1(Cell(1, edgeId))) {
1154 saddles->emplace_back(edgeId);
1155 }
1156
1157 const SimplexId pairedCellId
1158 = getPairedCell(Cell(1, edgeId), triangulation);
1159
1160 if(pairedCellId != -1 and pairedCellId != triangleId) {
1161 bfs.push(pairedCellId);
1162 }
1163 }
1164 }
1165 }
1166
1167 if(saddles != nullptr && saddles->size() > 1) {
1168 std::sort(saddles->begin(), saddles->end());
1169 const auto last = std::unique(saddles->begin(), saddles->end());
1170 saddles->erase(last, saddles->end());
1171 }
1172 }
1173 }
1174
1175 return 0;
1176}
1177
1178template <typename triangulationType>
1180 const Cell &cell,
1181 VisitedMask &mask,
1182 const triangulationType &triangulation,
1183 std::vector<Cell> *const wall,
1184 std::vector<SimplexId> *const saddles) const {
1185
1186 if(saddles != nullptr) {
1187 saddles->clear();
1188 }
1189
1190 if(dimensionality_ == 3) {
1191 if(cell.dim_ == 1) {
1192 // assume that cellId is an edge
1193 const SimplexId originId = cell.id_;
1194
1195 std::queue<SimplexId> bfs;
1196 bfs.push(originId);
1197
1198 // BFS traversal
1199 while(!bfs.empty()) {
1200 const SimplexId edgeId = bfs.front();
1201 bfs.pop();
1202
1203 if(!mask.isVisited_[edgeId]) {
1204 mask.isVisited_[edgeId] = true;
1205 mask.visitedIds_.emplace_back(edgeId);
1206
1207 // add the edge
1208 if(wall != nullptr) {
1209 wall->push_back(Cell(1, edgeId));
1210 }
1211
1212 const SimplexId triangleNumber
1213 = triangulation.getEdgeTriangleNumber(edgeId);
1214 for(SimplexId j = 0; j < triangleNumber; ++j) {
1215 SimplexId triangleId;
1216 triangulation.getEdgeTriangle(edgeId, j, triangleId);
1217
1218 if((saddles != nullptr) and isSaddle2(Cell(2, triangleId))) {
1219 saddles->emplace_back(triangleId);
1220 }
1221
1222 const SimplexId pairedCellId
1223 = getPairedCell(Cell(2, triangleId), triangulation, true);
1224
1225 if(pairedCellId != -1 and pairedCellId != edgeId) {
1226 bfs.push(pairedCellId);
1227 }
1228 }
1229 }
1230 }
1231
1232 if(saddles != nullptr && saddles->size() > 1) {
1233 std::sort(saddles->begin(), saddles->end());
1234 const auto last = std::unique(saddles->begin(), saddles->end());
1235 saddles->erase(last, saddles->end());
1236 }
1237 }
1238 }
1239
1240 return 0;
1241}
1242
1243template <typename triangulationType>
1245 const std::vector<Cell> &vpath,
1246 const triangulationType &triangulation) const {
1247
1248 if(dimensionality_ == 2) {
1249 // assume that the first cell is an edge
1250 const SimplexId numberOfCellsInPath = vpath.size();
1251 for(SimplexId i = 0; i < numberOfCellsInPath; i += 2) {
1252 const SimplexId edgeId = vpath[i].id_;
1253 const SimplexId triangleId = vpath[i + 1].id_;
1254
1255#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1256 for(int k = 0; k < 3; ++k) {
1257 SimplexId tmp;
1258 triangulation.getCellEdge(triangleId, k, tmp);
1259 if(tmp == edgeId) {
1260 (*gradient_)[3][triangleId] = k;
1261 break;
1262 }
1263 }
1264 for(int k = 0; k < triangulation.getEdgeStarNumber(edgeId); ++k) {
1265 SimplexId tmp;
1266 triangulation.getEdgeStar(edgeId, k, tmp);
1267 if(tmp == triangleId) {
1268 (*gradient_)[2][edgeId] = k;
1269 break;
1270 }
1271 }
1272#else
1273 TTK_FORCE_USE(triangulation);
1274 (*gradient_)[3][triangleId] = edgeId;
1275 (*gradient_)[2][edgeId] = triangleId;
1276#endif
1277 }
1278 } else if(dimensionality_ == 3) {
1279 // assume that the first cell is a triangle
1280 const SimplexId numberOfCellsInPath = vpath.size();
1281 for(SimplexId i = 0; i < numberOfCellsInPath; i += 2) {
1282 const SimplexId triangleId = vpath[i].id_;
1283 const SimplexId tetraId = vpath[i + 1].id_;
1284
1285#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1286 for(int k = 0; k < 4; ++k) {
1287 SimplexId tmp;
1288 triangulation.getCellTriangle(tetraId, k, tmp);
1289 if(tmp == triangleId) {
1290 (*gradient_)[5][tetraId] = k;
1291 break;
1292 }
1293 }
1294 for(int k = 0; k < triangulation.getTriangleStarNumber(triangleId); ++k) {
1295 SimplexId tmp;
1296 triangulation.getTriangleStar(triangleId, k, tmp);
1297 if(tmp == tetraId) {
1298 (*gradient_)[4][triangleId] = k;
1299 break;
1300 }
1301 }
1302#else
1303 (*gradient_)[5][tetraId] = triangleId;
1304 (*gradient_)[4][triangleId] = tetraId;
1305#endif
1306 }
1307 }
1308
1309 return 0;
1310}
1311
1312template <typename triangulationType>
1314 const std::vector<Cell> &vpath,
1315 const triangulationType &triangulation) const {
1316
1317 // assume that the first cell is an edge
1318 for(size_t i = 0; i < vpath.size(); i += 2) {
1319 const SimplexId edgeId = vpath[i].id_;
1320 const SimplexId vertId = vpath[i + 1].id_;
1321
1322#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1323 const auto nneighs = triangulation.getVertexEdgeNumber();
1324 for(int k = 0; k < nneighs; ++k) {
1325 SimplexId tmp;
1326 triangulation.getVertexEdge(vertId, k, tmp);
1327 if(tmp == edgeId) {
1328 (*gradient_)[0][vertId] = k;
1329 break;
1330 }
1331 }
1332 const auto nverts = triangulation.getEdgeStarNumber(edgeId);
1333 for(int k = 0; k < nverts; ++k) {
1334 SimplexId tmp;
1335 triangulation.getEdgeVertex(edgeId, k, tmp);
1336 if(tmp == vertId) {
1337 (*gradient_)[1][edgeId] = k;
1338 break;
1339 }
1340 }
1341#else
1342 TTK_FORCE_USE(triangulation);
1343 (*gradient_)[0][vertId] = edgeId;
1344 (*gradient_)[1][edgeId] = vertId;
1345#endif
1346 }
1347
1348 return 0;
1349}
1350
1351template <typename triangulationType>
1353 const std::vector<Cell> &vpath,
1354 const triangulationType &triangulation) const {
1355
1356 if(dimensionality_ == 3) {
1357 // assume that the first cell is an edge
1358 const SimplexId numberOfCellsInPath = vpath.size();
1359 for(SimplexId i = 0; i < numberOfCellsInPath; i += 2) {
1360 const SimplexId edgeId = vpath[i].id_;
1361 const SimplexId triangleId = vpath[i + 1].id_;
1362
1363#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1364 for(int k = 0; k < 3; ++k) {
1365 SimplexId tmp;
1366 triangulation.getTriangleEdge(triangleId, k, tmp);
1367 if(tmp == edgeId) {
1368 (*gradient_)[3][triangleId] = k;
1369 break;
1370 }
1371 }
1372 for(int k = 0; k < triangulation.getEdgeTriangleNumber(edgeId); ++k) {
1373 SimplexId tmp;
1374 triangulation.getEdgeTriangle(edgeId, k, tmp);
1375 if(tmp == triangleId) {
1376 (*gradient_)[2][edgeId] = k;
1377 break;
1378 }
1379 }
1380#else
1381 TTK_FORCE_USE(triangulation);
1382 (*gradient_)[3][triangleId] = edgeId;
1383 (*gradient_)[2][edgeId] = triangleId;
1384#endif
1385 }
1386 }
1387
1388 return 0;
1389}
1390
1391template <typename triangulationType>
1393 const std::vector<Cell> &vpath,
1394 const triangulationType &triangulation) const {
1395
1396 if(dimensionality_ == 3) {
1397 // assume that the first cell is a triangle
1398 const SimplexId numberOfCellsInPath = vpath.size();
1399 for(SimplexId i = 0; i < numberOfCellsInPath; i += 2) {
1400 const SimplexId triangleId = vpath[i].id_;
1401 const SimplexId edgeId = vpath[i + 1].id_;
1402
1403#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1404 for(int k = 0; k < 3; ++k) {
1405 SimplexId tmp;
1406 triangulation.getTriangleEdge(triangleId, k, tmp);
1407 if(tmp == edgeId) {
1408 (*gradient_)[3][triangleId] = k;
1409 break;
1410 }
1411 }
1412 for(int k = 0; k < triangulation.getEdgeTriangleNumber(edgeId); ++k) {
1413 SimplexId tmp;
1414 triangulation.getEdgeTriangle(edgeId, k, tmp);
1415 if(tmp == triangleId) {
1416 (*gradient_)[2][edgeId] = k;
1417 break;
1418 }
1419 }
1420#else
1421 TTK_FORCE_USE(triangulation);
1422 (*gradient_)[2][edgeId] = triangleId;
1423 (*gradient_)[3][triangleId] = edgeId;
1424#endif
1425 }
1426 }
1427
1428 return 0;
1429}
1430
1431template <typename triangulationType>
1433 const Cell c, const triangulationType &triangulation) const {
1434
1435 const auto offsets = this->inputOffsets_;
1436
1437 auto cellDim = c.dim_;
1438 auto cellId = c.id_;
1439
1440 SimplexId vertexId = -1;
1441 if(cellDim == 0) {
1442 vertexId = cellId;
1443 }
1444
1445 else if(cellDim == 1) {
1446 SimplexId v0;
1447 SimplexId v1;
1448 triangulation.getEdgeVertex(cellId, 0, v0);
1449 triangulation.getEdgeVertex(cellId, 1, v1);
1450
1451 if(offsets[v0] > offsets[v1]) {
1452 vertexId = v0;
1453 } else {
1454 vertexId = v1;
1455 }
1456 }
1457
1458 else if(cellDim == 2) {
1459 SimplexId v0{}, v1{}, v2{};
1460 triangulation.getTriangleVertex(cellId, 0, v0);
1461 triangulation.getTriangleVertex(cellId, 1, v1);
1462 triangulation.getTriangleVertex(cellId, 2, v2);
1463 if(offsets[v0] > offsets[v1] && offsets[v0] > offsets[v2]) {
1464 vertexId = v0;
1465 } else if(offsets[v1] > offsets[v0] && offsets[v1] > offsets[v2]) {
1466 vertexId = v1;
1467 } else {
1468 vertexId = v2;
1469 }
1470 }
1471
1472 else if(cellDim == 3) {
1473 SimplexId v0{}, v1{}, v2{}, v3{};
1474 triangulation.getCellVertex(cellId, 0, v0);
1475 triangulation.getCellVertex(cellId, 1, v1);
1476 triangulation.getCellVertex(cellId, 2, v2);
1477 triangulation.getCellVertex(cellId, 3, v3);
1478 if(offsets[v0] > offsets[v1] && offsets[v0] > offsets[v2]
1479 && offsets[v0] > offsets[v3]) {
1480 vertexId = v0;
1481 } else if(offsets[v1] > offsets[v0] && offsets[v1] > offsets[v2]
1482 && offsets[v1] > offsets[v3]) {
1483 vertexId = v1;
1484 } else if(offsets[v2] > offsets[v0] && offsets[v2] > offsets[v1]
1485 && offsets[v2] > offsets[v3]) {
1486 vertexId = v2;
1487 } else {
1488 vertexId = v3;
1489 }
1490 }
1491 return vertexId;
1492}
1493
1494template <typename triangulationType>
1496 const Cell c, const triangulationType &triangulation) const {
1497
1498 const auto offsets = this->inputOffsets_;
1499
1500 auto cellDim = c.dim_;
1501 auto cellId = c.id_;
1502
1503 SimplexId vertexId = -1;
1504 if(cellDim == 0) {
1505 vertexId = cellId;
1506 }
1507
1508 else if(cellDim == 1) {
1509 SimplexId v0;
1510 SimplexId v1;
1511 triangulation.getEdgeVertex(cellId, 0, v0);
1512 triangulation.getEdgeVertex(cellId, 1, v1);
1513
1514 if(offsets[v0] < offsets[v1]) {
1515 vertexId = v0;
1516 } else {
1517 vertexId = v1;
1518 }
1519 }
1520
1521 else if(cellDim == 2) {
1522 SimplexId v0{}, v1{}, v2{};
1523 triangulation.getTriangleVertex(cellId, 0, v0);
1524 triangulation.getTriangleVertex(cellId, 1, v1);
1525 triangulation.getTriangleVertex(cellId, 2, v2);
1526 if(offsets[v0] < offsets[v1] && offsets[v0] < offsets[v2]) {
1527 vertexId = v0;
1528 } else if(offsets[v1] < offsets[v0] && offsets[v1] < offsets[v2]) {
1529 vertexId = v1;
1530 } else {
1531 vertexId = v2;
1532 }
1533 }
1534
1535 else if(cellDim == 3) {
1536 SimplexId v0{}, v1{}, v2{}, v3{};
1537 triangulation.getCellVertex(cellId, 0, v0);
1538 triangulation.getCellVertex(cellId, 1, v1);
1539 triangulation.getCellVertex(cellId, 2, v2);
1540 triangulation.getCellVertex(cellId, 3, v3);
1541 if(offsets[v0] < offsets[v1] && offsets[v0] < offsets[v2]
1542 && offsets[v0] < offsets[v3]) {
1543 vertexId = v0;
1544 } else if(offsets[v1] < offsets[v0] && offsets[v1] < offsets[v2]
1545 && offsets[v1] < offsets[v3]) {
1546 vertexId = v1;
1547 } else if(offsets[v2] < offsets[v0] && offsets[v2] < offsets[v1]
1548 && offsets[v2] < offsets[v3]) {
1549 vertexId = v2;
1550 } else {
1551 vertexId = v3;
1552 }
1553 }
1554 return vertexId;
1555}
1556
1557template <typename triangulationType>
1559 std::vector<std::array<float, 3>> &points,
1560 std::vector<char> &points_pairOrigins,
1561 std::vector<char> &cells_pairTypes,
1562 std::vector<SimplexId> &cellIds,
1563 std::vector<char> &cellDimensions,
1564 const triangulationType &triangulation) const {
1565
1566 const auto nDims = this->getNumberOfDimensions();
1567
1568 // number of glyphs per dimension
1569 std::vector<size_t> nGlyphsPerDim(nDims);
1570
1571#ifdef TTK_ENABLE_OPENMP
1572#pragma omp parallel for num_threads(threadNumber_)
1573#endif // TTK_ENABLE_OPENMP
1574 for(int i = 0; i < nDims - 1; ++i) {
1575 const auto nCells = this->getNumberOfCells(i, triangulation);
1576 for(SimplexId j = 0; j < nCells; ++j) {
1577 if(this->getPairedCell(Cell{i, j}, triangulation) > -1) {
1578 nGlyphsPerDim[i]++;
1579 }
1580 }
1581 }
1582
1583 // partial sum of number of gradient glyphs
1584 std::vector<size_t> offsets(nDims + 1);
1585 for(SimplexId i = 0; i < nDims; ++i) {
1586 offsets[i + 1] = offsets[i] + nGlyphsPerDim[i];
1587 }
1588
1589 // total number of glyphs
1590 const auto nGlyphs = offsets.back();
1591
1592 // resize arrays accordingly
1593 points.resize(2 * nGlyphs);
1594 points_pairOrigins.resize(2 * nGlyphs);
1595 cells_pairTypes.resize(nGlyphs);
1596 cellIds.resize(2 * nGlyphs);
1597 cellDimensions.resize(2 * nGlyphs);
1598
1599#ifdef TTK_ENABLE_OPENMP
1600#pragma omp parallel for num_threads(threadNumber_)
1601#endif // TTK_ENABLE_OPENMP
1602 for(int i = 0; i < nDims - 1; ++i) {
1603 const SimplexId nCells = getNumberOfCells(i, triangulation);
1604 size_t nProcessedGlyphs{offsets[i]};
1605 for(SimplexId j = 0; j < nCells; ++j) {
1606 const Cell c{i, j};
1607 const auto pcid = this->getPairedCell(c, triangulation);
1608 if(pcid > -1) {
1609 const Cell pc{i + 1, pcid};
1610 triangulation.getCellIncenter(
1611 c.id_, c.dim_, points[2 * nProcessedGlyphs].data());
1612 triangulation.getCellIncenter(
1613 pc.id_, pc.dim_, points[2 * nProcessedGlyphs + 1].data());
1614 points_pairOrigins[2 * nProcessedGlyphs] = 0;
1615 points_pairOrigins[2 * nProcessedGlyphs + 1] = 1;
1616 cells_pairTypes[nProcessedGlyphs] = i;
1617#ifdef TTK_ENABLE_MPI
1618 ttk::SimplexId globalId{-1};
1619 triangulation.getDistributedGlobalCellId(j, i, globalId);
1620 cellIds[2 * nProcessedGlyphs + 0] = globalId;
1621 triangulation.getDistributedGlobalCellId(pcid, i + 1, globalId);
1622 cellIds[2 * nProcessedGlyphs + 1] = globalId;
1623#else
1624 cellIds[2 * nProcessedGlyphs + 0] = j;
1625 cellIds[2 * nProcessedGlyphs + 1] = pcid;
1626#endif // TTK_ENABLE_MPI
1627 cellDimensions[2 * nProcessedGlyphs + 0] = i;
1628 cellDimensions[2 * nProcessedGlyphs + 1] = i + 1;
1629 nProcessedGlyphs++;
1630 }
1631 }
1632 }
1633
1634 return 0;
1635}
#define TTK_FORCE_USE(x)
Force the compiler to use the function/method parameter.
Definition: BaseClass.h:57
std::array< std::vector< gradIdType >, 6 > gradientType
Discrete gradient internal struct.
int threadNumber_
Definition: BaseClass.h:95
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:159
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:149
TTK discreteGradient processing package.
SimplexId getCellLowerVertex(const Cell c, const triangulationType &triangulation) const
int reverseAscendingPath(const std::vector< Cell > &vpath, const triangulationType &triangulation) const
int reverseDescendingPath(const std::vector< Cell > &vpath, const triangulationType &triangulation) const
int setGradientGlyphs(std::vector< std::array< float, 3 > > &points, std::vector< char > &points_pairOrigins, std::vector< char > &cells_pairTypes, std::vector< SimplexId > &cellsIds, std::vector< char > &cellsDimensions, const triangulationType &triangulation) const
bool getDescendingPathThroughWall(const Cell &saddle2, const Cell &saddle1, const std::vector< bool > &isVisited, std::vector< Cell > *const vpath, const triangulationType &triangulation, const bool stopIfMultiConnected=false, const bool enableCycleDetector=false) const
bool isSaddle1(const Cell &cell) const
bool getAscendingPathThroughWall(const Cell &saddle1, const Cell &saddle2, const std::vector< bool > &isVisited, std::vector< Cell > *const vpath, const triangulationType &triangulation, const bool stopIfMultiConnected=false, const bool enableCycleDetector=false) const
bool isSaddle2(const Cell &cell) const
AbstractTriangulation::gradientKeyType inputScalarField_
SimplexId getNumberOfCells(const int dimension, const triangulationType &triangulation) const
int reverseAscendingPathOnWall(const std::vector< Cell > &vpath, const triangulationType &triangulation) const
bool isBoundary(const Cell &cell, const triangulationType &triangulation) const
bool isCellCritical(const int cellDim, const SimplexId cellId) const
SimplexId getPairedCell(const Cell &cell, const triangulationType &triangulation, bool isReverse=false) const
int getAscendingWall(const Cell &cell, VisitedMask &mask, const triangulationType &triangulation, std::vector< Cell > *const wall=nullptr, std::vector< SimplexId > *const saddles=nullptr) const
int reverseDescendingPathOnWall(const std::vector< Cell > &vpath, const triangulationType &triangulation) const
SimplexId getCellGreaterVertex(const Cell c, const triangulationType &triangulation) const
int getCriticalPoints(std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, const triangulationType &triangulation) const
dataType getPersistence(const Cell &up, const Cell &down, const dataType *const scalars, const triangulationType &triangulation) const
int buildGradient(const triangulationType &triangulation, bool bypassCache=false)
int setCriticalPoints(const std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, std::vector< std::array< float, 3 > > &points, std::vector< char > &cellDimensions, std::vector< SimplexId > &cellIds, std::vector< char > &isOnBoundary, std::vector< SimplexId > &PLVertexIdentifiers, const triangulationType &triangulation) const
int getDescendingPath(const Cell &cell, std::vector< Cell > &vpath, const triangulationType &triangulation) const
int getAscendingPath(const Cell &cell, std::vector< Cell > &vpath, const triangulationType &triangulation, const bool enableCycleDetector=false) const
AbstractTriangulation::gradientType localGradient_
int getDescendingWall(const Cell &cell, VisitedMask &mask, const triangulationType &triangulation, std::vector< Cell > *const wall=nullptr, std::vector< SimplexId > *const saddles=nullptr) const
AbstractTriangulation::gradientType * gradient_
const SimplexId * inputOffsets_
COMMON_EXPORTS int MPIrank_
Definition: BaseClass.cpp:9
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
Auto-cleaning re-usable graph propagations data structure.
Definition: VisitedMask.h:27
std::vector< bool > & isVisited_
Definition: VisitedMask.h:28
std::vector< SimplexId > & visitedIds_
Definition: VisitedMask.h:29
Extended Cell structure for processLowerStars.
const std::array< SimplexId, 3 > lowVerts_
const std::array< uint8_t, 3 > faces_