TTK
Loading...
Searching...
No Matches
IntegralLines.h
Go to the documentation of this file.
1
20
21#pragma once
22
23// base code includes
24#include <ArrayLinkedList.h>
25#include <Geometry.h>
26#include <ScalarFieldCriticalPoints.h>
27#include <Triangulation.h>
28// std includes
29#include <algorithm>
30#include <iterator>
31#include <limits>
32#include <numeric>
33#include <string>
34#include <unordered_set>
35
36#define INTEGRAL_LINE_TABULAR_SIZE 50
37#ifdef TTK_ENABLE_MPI
38#define INTEGRAL_LINE_IS_ELEMENT_TO_PROCESS 0
39#define INTEGRAL_LINE_IS_MESSAGE_SIZE 1
40#endif
41
42namespace ttk {
43 namespace intgl {
44
54 struct IntegralLine {
55 std::vector<ttk::SimplexId> trajectory;
56 std::vector<double> distanceFromSeed;
57 std::vector<ttk::SimplexId> localVertexIdentifier;
58 ttk::SimplexId seedIdentifier;
59 ttk::SimplexId forkIdentifier = -1;
60 };
61
62#ifdef TTK_ENABLE_MPI
70 struct GhostElementsToSort {
71 ttk::SimplexId ownedGlobalId;
72 ttk::SimplexId minLocalVertexId;
73 ttk::SimplexId seedIdentifier;
74 ttk::SimplexId forkIdentifier{-1};
75 ttk::SimplexId globalEdgeId{-1};
76 ttk::SimplexId ghostVertexLocalId;
77 ttk::SimplexId ghostEdgeLocalId;
78 };
79
80 inline bool operator<(const GhostElementsToSort &left,
81 const GhostElementsToSort &right) {
82 if(left.seedIdentifier != right.seedIdentifier) {
83 return left.seedIdentifier < right.seedIdentifier;
84 }
85 if(left.forkIdentifier != right.forkIdentifier) {
86 return left.forkIdentifier < right.forkIdentifier;
87 }
88 return left.minLocalVertexId < right.minLocalVertexId;
89 };
90
91 /* Counters for keeping tracks of number of integral lines
92 * left to compute. Each counter is aligned on a cache line
93 * to prevent false-sharing
94 */
95 static int finishedElement_ __attribute__((aligned(64)));
96 static int addedElement_ __attribute__((aligned(64)));
97
98 /*
99 * For each integral line continuing on another process, we send one layer
100 * of ghost cells and we keep one layers for this process, meaning that
101 * vertices Id1 belongs to this process and Id2 belongs to a neighboring
102 * process.
103 */
104 struct ElementToBeSent {
105 ttk::SimplexId Id1;
106 ttk::SimplexId Id2;
107 double DistanceFromSeed1;
108 double DistanceFromSeed2;
109 ttk::SimplexId LocalVertexIdentifier1;
110 ttk::SimplexId LocalVertexIdentifier2;
111 ttk::SimplexId SeedIdentifier;
112 ttk::SimplexId ForkIdentifier;
113 };
114#endif
115 } // namespace intgl
116
117 enum Direction { Forward = 0, Backward };
118
119 class IntegralLines : virtual public Debug {
120
121 public:
122 IntegralLines();
123 ~IntegralLines() override;
124
125 template <class triangulationType = ttk::AbstractTriangulation>
126 int execute(triangulationType *triangulation);
127
137 template <class triangulationType = ttk::AbstractTriangulation>
138 void computeIntegralLine(const triangulationType *triangulation,
139 ttk::intgl::IntegralLine *integralLine,
140 const ttk::SimplexId *offsets) const;
141
152 template <class triangulationType>
153 void createTask(const triangulationType *triangulation,
154 std::vector<ttk::intgl::IntegralLine *> &chunkIntegralLine,
155 const ttk::SimplexId *offsets,
156 int nbElement) const;
171 template <class triangulationType>
172 void
173 prepareForTask(const triangulationType *triangulation,
174 std::vector<ttk::intgl::IntegralLine *> &chunkIntegralLine,
175 int startingIndex,
176 int nbElement,
177 std::vector<SimplexId> *seeds) const;
178
179 inline void setVertexNumber(const SimplexId &vertexNumber) {
180 vertexNumber_ = vertexNumber;
181 }
182
183 inline void setSeedNumber(const SimplexId &seedNumber) {
184 seedNumber_ = seedNumber;
185 }
186
187 inline void setDirection(int direction) {
188 direction_ = direction;
189 }
190
200 void findNextVertex(ttk::SimplexId &vnext,
201 ttk::SimplexId &fnext,
202 std::vector<ttk::SimplexId> &component,
203 const SimplexId *offsets) const;
204
205#ifdef TTK_ENABLE_MPI
206
217 template <typename triangulationType>
218 int getGlobalIdentifiers(
219 std::vector<ttk::SimplexId> &globalVertexId,
220 std::vector<ttk::SimplexId> &globalCellId,
221 const std::vector<ttk::ArrayLinkedList<ttk::intgl::IntegralLine,
222 INTEGRAL_LINE_TABULAR_SIZE>>
223 &integralLines,
224 const triangulationType *triangulation);
225
236 int exchangeGhosts(
237 const std::vector<std::vector<ttk::intgl::GhostElementsToSort>>
238 &unmatchedGhosts,
239 std::vector<ttk::SimplexId> &globalVertexId,
240 std::vector<ttk::SimplexId> &globalCellId);
241
253 template <class triangulationType>
254 void storeToSendIfNecessary(const triangulationType *triangulation,
255 ttk::intgl::IntegralLine *integralLine,
256 bool &isMax) const;
257
273 template <class triangulationType>
274 void
275 receiveElement(const triangulationType *triangulation,
276 intgl::ElementToBeSent &element,
277 std::vector<ttk::intgl::IntegralLine *> &chunkIntegralLine,
278 int &index,
279 int taskSize,
280 const ttk::SimplexId *offsets);
281
282 inline void setToSend(
283 std::vector<std::vector<std::vector<intgl::ElementToBeSent>>> *send) {
284 toSend_ = send;
285 }
286
287 inline void setGlobalElementCounter(int counter) {
288 globalElementCounter_ = counter;
289 }
290
291 inline void setNeighbors(const std::vector<int> &neighbors) {
292 neighbors_ = neighbors;
293 neighborNumber_ = neighbors_.size();
294 int idx = 0;
295 for(int neighbor : neighbors) {
296 neighborsToId_[neighbor] = idx;
297 idx++;
298 }
299 }
305 void createMessageType() {
306 ttk::SimplexId id = 0;
307 MPI_Datatype types[]
308 = {getMPIType(id), getMPIType(id), MPI_DOUBLE, MPI_DOUBLE,
309 getMPIType(id), getMPIType(id), getMPIType(id), getMPIType(id)};
310 int lengths[] = {1, 1, 1, 1, 1, 1, 1, 1};
311 const long int mpi_offsets[]
312 = {offsetof(intgl::ElementToBeSent, Id1),
313 offsetof(intgl::ElementToBeSent, Id2),
314 offsetof(intgl::ElementToBeSent, DistanceFromSeed1),
315 offsetof(intgl::ElementToBeSent, DistanceFromSeed2),
316 offsetof(intgl::ElementToBeSent, LocalVertexIdentifier1),
317 offsetof(intgl::ElementToBeSent, LocalVertexIdentifier2),
318 offsetof(intgl::ElementToBeSent, SeedIdentifier),
319 offsetof(intgl::ElementToBeSent, ForkIdentifier)};
320 MPI_Type_create_struct(
321 8, lengths, mpi_offsets, types, &(this->MessageType));
322 MPI_Type_commit(&(this->MessageType));
323 }
324#endif
325
326 int preconditionTriangulation(ttk::AbstractTriangulation *triangulation) {
327 int status = triangulation->preconditionVertexNeighbors();
328 // For critical points computation
329 status += triangulation->preconditionVertexStars();
330 status += triangulation->preconditionBoundaryVertices();
331 this->scalarFieldCriticalPoints_.preconditionTriangulation(triangulation);
332 return status;
333 }
334
335 inline void setInputScalarField(void *data) {
336 inputScalarField_ = data;
337 }
338
347 inline void setInputOffsets(const SimplexId *const data) {
348 inputOffsets_ = data;
349 }
350
351 inline void
352 setVertexIdentifierScalarField(std::vector<SimplexId> *const data) {
353 vertexIdentifierScalarField_ = data;
354 }
355
356 inline void setOutputIntegralLines(
357 std::vector<ttk::ArrayLinkedList<ttk::intgl::IntegralLine,
358 INTEGRAL_LINE_TABULAR_SIZE>>
359 *integralLines) {
360 this->outputIntegralLines_ = integralLines;
361 }
362
363 inline void setChunkSize(int size) {
364 chunkSize_ = size;
365 }
366
370
371 protected:
372 ttk::SimplexId vertexNumber_;
373 ttk::SimplexId seedNumber_;
374 ttk::SimplexId chunkSize_;
375 ttk::SimplexId direction_;
376 void *inputScalarField_;
377 const ttk::SimplexId *inputOffsets_;
378 std::vector<ttk::SimplexId> *vertexIdentifierScalarField_;
379 std::vector<ttk::ArrayLinkedList<ttk::intgl::IntegralLine,
380 INTEGRAL_LINE_TABULAR_SIZE>>
381 *outputIntegralLines_;
382 ttk::ScalarFieldCriticalPoints scalarFieldCriticalPoints_;
383 bool EnableForking{false};
384
385#ifdef TTK_ENABLE_MPI
386 std::vector<std::vector<std::vector<intgl::ElementToBeSent>>> *toSend_{
387 nullptr};
388 int neighborNumber_;
389 std::unordered_map<int, int> neighborsToId_;
390 std::vector<int> neighbors_;
391 SimplexId keepWorking_;
392 SimplexId globalElementCounter_;
393 MPI_Datatype MessageType;
394#endif
395 };
396} // namespace ttk
397
398#ifdef TTK_ENABLE_MPI
399
400template <class triangulationType>
401void ttk::IntegralLines::receiveElement(
402 const triangulationType *triangulation,
403 intgl::ElementToBeSent &element,
404 std::vector<ttk::intgl::IntegralLine *> &chunkIntegralLine,
405 int &index,
406 int taskSize,
407 const ttk::SimplexId *offsets) {
408
409 // Create integral line object on this process
410 int threadNum{0};
411#ifdef TTK_ENABLE_OPENMP4
412 threadNum = omp_get_thread_num();
413#endif // TTK_ENABLE_OPENMP4
414
415 ttk::intgl::IntegralLine *integralLine
416 = outputIntegralLines_->at(threadNum).addArrayElement(
417 ttk::intgl::IntegralLine{
418 std::vector<ttk::SimplexId>(
419 {triangulation->getVertexLocalId(element.Id1),
420 triangulation->getVertexLocalId(element.Id2)}),
421 std::vector<double>(
422 {element.DistanceFromSeed1, element.DistanceFromSeed2}),
423 std::vector<ttk::SimplexId>(
424 {element.LocalVertexIdentifier1, element.LocalVertexIdentifier2}),
425 element.SeedIdentifier, element.ForkIdentifier});
426
427 // Add to chunks for task granularity
428 chunkIntegralLine[index] = integralLine;
429
430 // Start task if necessary
431 if(index == taskSize - 1) {
432 this->createTask<triangulationType>(
433 triangulation, chunkIntegralLine, offsets, taskSize);
434 index = 0;
435 } else {
436 index++;
437 }
438}
439
440template <class triangulationType>
441void ttk::IntegralLines::storeToSendIfNecessary(
442 const triangulationType *triangulation,
443 ttk::intgl::IntegralLine *integralLine,
444 bool &isMax) const {
445#ifdef TTK_ENABLE_MPI
446 if(ttk::isRunningWithMPI()) {
447 int size = integralLine->trajectory.size();
448 if(size > 1) {
449 int rankArray
450 = triangulation->getVertexRank(integralLine->trajectory.back());
451 if(rankArray != ttk::MPIrank_) {
452 intgl::ElementToBeSent element
453 = intgl::ElementToBeSent{-1,
454 -1,
455 0,
456 0,
457 -1,
458 -1,
459 integralLine->seedIdentifier,
460 integralLine->forkIdentifier};
461 element.Id2
462 = triangulation->getVertexGlobalId(integralLine->trajectory.back());
463 element.Id1 = triangulation->getVertexGlobalId(
464 integralLine->trajectory.at(size - 2));
465 element.DistanceFromSeed2 = integralLine->distanceFromSeed.back();
466 element.DistanceFromSeed1 = integralLine->distanceFromSeed.at(size - 2);
467 element.LocalVertexIdentifier2
468 = integralLine->localVertexIdentifier.back();
469 element.LocalVertexIdentifier1
470 = integralLine->localVertexIdentifier.at(size - 2);
471#ifdef TTK_ENABLE_OPENMP4
472 toSend_
473 ->at(neighborsToId_.find(rankArray)->second)[omp_get_thread_num()]
474 .push_back(element);
475#else
476 toSend_->at(neighborsToId_.find(rankArray)->second)[0].push_back(
477 element);
478#endif // TTK_ENABLE_OPENMP4
479 isMax = true;
480 }
481 }
482 }
483#endif
484}
485#endif
486
487inline void
488 ttk::IntegralLines::findNextVertex(ttk::SimplexId &vnext,
489 ttk::SimplexId &fnext,
490 std::vector<ttk::SimplexId> &component,
491 const SimplexId *offsets) const {
492 ttk::SimplexId const elementInComponentNumber = component.size();
493 for(ttk::SimplexId k = 0; k < elementInComponentNumber; ++k) {
494 if(direction_ == static_cast<int>(Direction::Forward)) {
495 if(fnext < offsets[component[k]]) {
496 vnext = component[k];
497 fnext = offsets[component[k]];
498 }
499 } else {
500 if(fnext > offsets[component[k]]) {
501 vnext = component[k];
502 fnext = offsets[component[k]];
503 }
504 }
505 }
506}
507
508template <class triangulationType>
509void ttk::IntegralLines::computeIntegralLine(
510 const triangulationType *triangulation,
511 ttk::intgl::IntegralLine *integralLine,
512 const SimplexId *offsets) const {
513 double distance = integralLine->distanceFromSeed.back();
514 ttk::SimplexId v = integralLine->trajectory.back();
515 float p0[3];
516 float p1[3];
517 triangulation->getVertexPoint(v, p0[0], p0[1], p0[2]);
518 bool isMax{};
519 std::vector<std::vector<ttk::SimplexId>> *components;
520 while(!isMax) {
521 std::vector<std::vector<ttk::SimplexId>> upperComponents;
522 std::vector<std::vector<ttk::SimplexId>> lowerComponents;
523 // Computation of the critical type
524 char const criticalType
525 = this->scalarFieldCriticalPoints_.getCriticalType<triangulationType>(
526 v, offsets, triangulation, &upperComponents, &lowerComponents);
527 // End of integral line if an appropriate maxima is reached
528 if((criticalType == (char)CriticalType::Local_maximum
529 && direction_ == static_cast<int>(Direction::Forward))
530 || (criticalType == (char)CriticalType::Local_minimum
531 && direction_ != static_cast<int>(Direction::Forward))) {
532 isMax = true;
533#ifdef TTK_ENABLE_MPI
534 if(ttk::isRunningWithMPI()
535 && triangulation->getVertexRank(v) == ttk::MPIrank_) {
536#ifdef TTK_ENABLE_OPENMP4
537#pragma omp atomic update seq_cst
538#endif // TTK_ENABLE_OPENMP4
539 ttk::intgl::finishedElement_++;
540 } else {
541 this->storeToSendIfNecessary<triangulationType>(
542 triangulation, integralLine, isMax);
543 }
544#endif
545 } else {
546 // If the integral line continues, the computation continues using lower
547 // and upper components
548 components = &(lowerComponents);
549 if(direction_ == static_cast<int>(Direction::Forward)) {
550 components = &(upperComponents);
551 }
552 if((criticalType == (char)CriticalType::Saddle1
553 || criticalType == (char)CriticalType::Saddle2
554 || criticalType == (char)CriticalType::Degenerate)
555 && EnableForking) {
556 // For each connected components, the max (or the min) is computed
557 // and a task is created to further the computation of the integral line
558 ttk::SimplexId const numberOfComponents = components->size();
559#ifdef TTK_ENABLE_MPI
560#ifdef TTK_ENABLE_OPENMP4
561#pragma omp atomic update seq_cst
562#endif // TTK_ENABLE_OPENMP4
563 ttk::intgl::finishedElement_++;
564#ifdef TTK_ENABLE_OPENMP4
565#pragma omp atomic update seq_cst
566#endif // TTK_ENABLE_OPENMP4
567 ttk::intgl::addedElement_ += numberOfComponents;
568#endif
569 isMax = true;
570 for(int i = 0; i < numberOfComponents; i++) {
571 SimplexId vnext = -1;
572 ttk::SimplexId fnext = offsets[v];
573 this->findNextVertex(vnext, fnext, components->at(i), offsets);
574 // The task is created
575#ifdef TTK_ENABLE_MPI
576 ttk::SimplexId forkIdentifier
577 = triangulation->getVertexGlobalId(vnext);
578#else
579 ttk::SimplexId const forkIdentifier = vnext;
580#endif
581 triangulation->getVertexPoint(vnext, p1[0], p1[1], p1[2]);
582 double const distanceFork = Geometry::distance(p0, p1, 3);
583 int threadNum{0};
584#ifdef TTK_ENABLE_OPENMP4
585 threadNum = omp_get_thread_num();
586#endif // TTK_ENABLE_OPENMP4
587 ttk::intgl::IntegralLine *integralLineFork
588 = outputIntegralLines_->at(threadNum).addArrayElement(
589 ttk::intgl::IntegralLine{
590 std::vector<ttk::SimplexId>({v, vnext}),
591 std::vector<double>({distance, distance + distanceFork}),
592 std::vector<ttk::SimplexId>(
593 {integralLine->localVertexIdentifier.back(),
594 integralLine->localVertexIdentifier.back() + 1}),
595 integralLine->seedIdentifier, forkIdentifier});
596
597#ifdef TTK_ENABLE_OPENMP4
598#pragma omp task firstprivate(integralLineFork)
599 {
600#endif // TTK_ENABLE_OPENMP4
601#ifdef TTK_ENABLE_MPI
602 bool hasBeenSent = false;
603 this->storeToSendIfNecessary<triangulationType>(
604 triangulation, integralLineFork, hasBeenSent);
605 if(!hasBeenSent) {
606#endif
607 this->computeIntegralLine<triangulationType>(
608 triangulation, integralLineFork, offsets);
609#ifdef TTK_ENABLE_MPI
610 }
611#endif
612#ifdef TTK_ENABLE_OPENMP4
613 }
614#endif // TTK_ENABLE_OPENMP4
615 }
616 } else {
617 // In case the vertex is not a saddle point, all neighbor vertices
618 // are used for the computation
619 components->clear();
620 ttk::SimplexId const neighborNumber
621 = triangulation->getVertexNeighborNumber(v);
622 components->push_back(std::vector<ttk::SimplexId>());
623 ttk::SimplexId id;
624 for(ttk::SimplexId i = 0; i < neighborNumber; i++) {
625 triangulation->getVertexNeighbor(v, i, id);
626 components->at(0).push_back(id);
627 }
628 SimplexId vnext = -1;
629 ttk::SimplexId fnext = offsets[v];
630 this->findNextVertex(vnext, fnext, components->at(0), offsets);
631 triangulation->getVertexPoint(vnext, p1[0], p1[1], p1[2]);
632 distance += Geometry::distance(p0, p1, 3);
633 integralLine->trajectory.push_back(vnext);
634 p0[0] = p1[0];
635 p0[1] = p1[1];
636 p0[2] = p1[2];
637 integralLine->distanceFromSeed.push_back(distance);
638 integralLine->localVertexIdentifier.push_back(
639 integralLine->localVertexIdentifier.back() + 1);
640 v = vnext;
641#ifdef TTK_ENABLE_MPI
642 this->storeToSendIfNecessary<triangulationType>(
643 triangulation, integralLine, isMax);
644#endif
645 }
646 }
647 }
648}
649
650template <class triangulationType>
651void ttk::IntegralLines::prepareForTask(
652#ifdef TTK_ENABLE_MPI
653 const triangulationType *triangulation,
654#else
655 const triangulationType *ttkNotUsed(triangulation),
656#endif
657 std::vector<ttk::intgl::IntegralLine *> &chunkIntegralLine,
658 int startingIndex,
659 int nbElement,
660 std::vector<SimplexId> *seeds) const {
661
662 for(SimplexId j = 0; j < nbElement; j++) {
663 SimplexId const v{seeds->at(j + startingIndex)};
664 int seedIdentifier;
665#ifdef TTK_ENABLE_MPI
666 seedIdentifier = triangulation->getVertexGlobalId(v);
667#else
668 seedIdentifier = v;
669#endif
670 int threadNum{0};
671#ifdef TTK_ENABLE_OPENMP4
672 threadNum = omp_get_thread_num();
673#endif // TTK_ENABLE_OPENMP4
674 chunkIntegralLine[j] = outputIntegralLines_->at(threadNum).addArrayElement(
675 ttk::intgl::IntegralLine{
676 std::vector<ttk::SimplexId>(1, v), std::vector<double>(1, 0),
677 std::vector<ttk::SimplexId>(1, 0), seedIdentifier, -1});
678 }
679}
680
681template <class triangulationType>
682void ttk::IntegralLines::createTask(
683 const triangulationType *triangulation,
684 std::vector<ttk::intgl::IntegralLine *> &chunkIntegralLine,
685 const ttk::SimplexId *offsets,
686 int nbElement) const {
687#ifdef TTK_ENABLE_OPENMP4
688#pragma omp task firstprivate(chunkIntegralLine)
689 {
690#endif // TTK_ENABLE_OPENMP4
691 for(int j = 0; j < nbElement; j++) {
692 this->computeIntegralLine<triangulationType>(
693 triangulation, chunkIntegralLine[j], offsets);
694 }
695#ifdef TTK_ENABLE_OPENMP4
696 }
697#endif // TTK_ENABLE_OPENMP4
698}
699
700template <class triangulationType>
701int ttk::IntegralLines::execute(triangulationType *triangulation) {
702
703#ifdef TTK_ENABLE_MPI
704 keepWorking_ = 1;
705 ttk::intgl::finishedElement_ = 0;
706 ttk::intgl::addedElement_ = 0;
707#endif
708 const SimplexId *offsets = inputOffsets_;
709 std::vector<SimplexId> *seeds = vertexIdentifierScalarField_;
710 Timer t;
711
712 std::vector<ttk::intgl::IntegralLine *> chunkIntegralLine(chunkSize_);
713 int const taskNumber = (int)seedNumber_ / chunkSize_;
714#ifdef TTK_ENABLE_OPENMP4
715#ifdef TTK_ENABLE_MPI
716#pragma omp parallel shared( \
717 ttk::intgl::finishedElement_, toSend_, ttk::intgl::addedElement_) \
718 num_threads(threadNumber_)
719 {
720#else
721#pragma omp parallel num_threads(threadNumber_)
722 {
723#endif // TTK_ENABLE_MPI
724#pragma omp master
725 {
726#endif // TTK_ENABLE_OPENMP4
727 for(SimplexId i = 0; i < taskNumber; ++i) {
728 this->prepareForTask<triangulationType>(
729 triangulation, chunkIntegralLine, i * chunkSize_, chunkSize_, seeds);
730 this->createTask<triangulationType>(
731 triangulation, chunkIntegralLine, offsets, chunkSize_);
732 }
733 int const rest = seedNumber_ % chunkSize_;
734 if(rest > 0) {
735 this->prepareForTask<triangulationType>(
736 triangulation, chunkIntegralLine, taskNumber * chunkSize_, rest,
737 seeds);
738 this->createTask<triangulationType>(
739 triangulation, chunkIntegralLine, offsets, rest);
740 }
741#ifdef TTK_ENABLE_OPENMP4
742 }
743 }
744#endif
745#ifdef TTK_ENABLE_MPI
746 if(ttk::isRunningWithMPI()) {
747 int i;
748 int finishedElementReceived = 0;
749 std::vector<int> sendMessageSize(neighborNumber_);
750 std::vector<int> recvMessageSize(neighborNumber_);
751 std::vector<std::vector<intgl::ElementToBeSent>> send_buf(neighborNumber_);
752 std::vector<std::vector<intgl::ElementToBeSent>> recv_buf(neighborNumber_);
753 for(i = 0; i < neighborNumber_; i++) {
754 send_buf.reserve((int)seedNumber_ * 0.005);
755 recv_buf.reserve((int)seedNumber_ * 0.005);
756 }
757 std::vector<MPI_Request> requests(2 * neighborNumber_, MPI_REQUEST_NULL);
758 std::vector<MPI_Status> statuses(4 * neighborNumber_);
759 int taskSize;
760 int index;
761 int totalMessageSize;
762 while(keepWorking_) {
763 ttk::intgl::finishedElement_ -= ttk::intgl::addedElement_;
764 ttk::intgl::addedElement_ = 0;
765 // Exchange of the number of integral lines finished on all processes
766 MPI_Allreduce(&ttk::intgl::finishedElement_, &finishedElementReceived, 1,
767 MPI_INTEGER, MPI_SUM, ttk::MPIcomm_);
768 ttk::intgl::finishedElement_ = 0;
769 // Update the number of integral lines left to compute
770 globalElementCounter_ -= finishedElementReceived;
771 // Stop working in case there are no more computation to be done
772 if(globalElementCounter_ == 0) {
773 keepWorking_ = 0;
774 }
775 if(keepWorking_) {
776 totalMessageSize = 0;
777 // Preparation of the send buffers and exchange of the size of messages
778 // to be sent
779 for(i = 0; i < neighborNumber_; i++) {
780 for(int j = 0; j < threadNumber_; j++) {
781 send_buf[i].insert(send_buf[i].end(), toSend_->at(i)[j].begin(),
782 toSend_->at(i)[j].end());
783 toSend_->at(i)[j].clear();
784 }
785 sendMessageSize[i] = (int)send_buf[i].size();
786 MPI_Isend(&sendMessageSize[i], 1, MPI_INTEGER, neighbors_.at(i),
787 INTEGRAL_LINE_IS_MESSAGE_SIZE, ttk::MPIcomm_,
788 &requests[2 * i]);
789 MPI_Irecv(&recvMessageSize[i], 1, MPI_INTEGER, neighbors_.at(i),
790 INTEGRAL_LINE_IS_MESSAGE_SIZE, ttk::MPIcomm_,
791 &requests[2 * i + 1]);
792 }
793 MPI_Waitall(2 * neighborNumber_, requests.data(), MPI_STATUSES_IGNORE);
794 // Exchange of the data
795 for(i = 0; i < neighborNumber_; i++) {
796 if(recv_buf[i].size() < (size_t)recvMessageSize[i]) {
797 recv_buf[i].resize(recvMessageSize[i]);
798 }
799 if(recvMessageSize[i] > 0) {
800 MPI_Irecv(recv_buf[i].data(), recvMessageSize[i], this->MessageType,
801 neighbors_.at(i), INTEGRAL_LINE_IS_ELEMENT_TO_PROCESS,
802 ttk::MPIcomm_, &requests[2 * i]);
803 totalMessageSize += recvMessageSize[i];
804 }
805
806 if(sendMessageSize[i] > 0) {
807 MPI_Isend(send_buf[i].data(), sendMessageSize[i], this->MessageType,
808 neighbors_.at(i), INTEGRAL_LINE_IS_ELEMENT_TO_PROCESS,
809 ttk::MPIcomm_, &requests[2 * i + 1]);
810 }
811 }
812 MPI_Waitall(2 * neighborNumber_, requests.data(), MPI_STATUSES_IGNORE);
813 for(i = 0; i < neighborNumber_; i++) {
814 send_buf[i].clear();
815 }
816 // Extraction of the received data and creation of the tasks
817#ifdef TTK_ENABLE_OPENMP4
818#pragma omp parallel shared(ttk::intgl::finishedElement_, toSend_) \
819 num_threads(threadNumber_)
820 {
821#pragma omp master
822 {
823#endif // TTK_ENABLE_OPENMP4
824 index = 0;
825 taskSize = std::min(
826 (ttk::SimplexId)std::max(totalMessageSize / (threadNumber_ * 100),
827 std::min(totalMessageSize, 50)),
828 chunkSize_);
829 chunkIntegralLine.resize(taskSize);
830 for(i = 0; i < neighborNumber_; i++) {
831 for(int j = 0; j < recvMessageSize[i]; j++) {
832 this->receiveElement<triangulationType>(
833 triangulation, recv_buf[i][j], chunkIntegralLine, index,
834 taskSize, offsets);
835 }
836 }
837 if(index > 0) {
838 this->createTask<triangulationType>(
839 triangulation, chunkIntegralLine, offsets, index);
840 }
841#ifdef TTK_ENABLE_OPENMP4
842 }
843 }
844#endif // TTK_ENABLE_OPENMP4
845 }
846 }
847 }
848#endif // TTK_ENABLE_MPI
849 {
850 std::stringstream msg;
851 msg << "Processed " << vertexNumber_ << " points";
852 this->printMsg(msg.str(), 1, t.getElapsedTime(), threadNumber_);
853 }
854 return 0;
855}
856
857#ifdef TTK_ENABLE_MPI
858template <typename triangulationType>
859int ttk::IntegralLines::getGlobalIdentifiers(
860 std::vector<ttk::SimplexId> &globalVertexId,
861 std::vector<ttk::SimplexId> &globalCellId,
862 const std::vector<
863 ttk::ArrayLinkedList<ttk::intgl::IntegralLine, INTEGRAL_LINE_TABULAR_SIZE>>
864 &integralLines,
865 const triangulationType *triangulation) {
866 ttk::SimplexId outputVertexNumber = 0;
867 ttk::SimplexId outputCellNumber = 0;
868 ttk::SimplexId realVertexNumber = 0;
869 ttk::SimplexId realCellNumber = 0;
870 ttk::SimplexId intervalSize;
871 // Counts vertices and edges number (with and without ghosts)
872#ifdef TTK_ENABLE_OPENMP4
873#pragma omp parallel for reduction(+ : outputVertexNumber, outputCellNumber, \
874 realCellNumber, realVertexNumber) \
875 schedule(static, 1) private(intervalSize)
876#endif
877 for(int thread = 0; thread < threadNumber_; thread++) {
878 std::list<std::array<ttk::intgl::IntegralLine,
879 INTEGRAL_LINE_TABULAR_SIZE>>::const_iterator
880 integralLine
881 = integralLines[thread].list_.begin();
882 while(integralLine != integralLines[thread].list_.end()) {
883 for(int i = 0; i < INTEGRAL_LINE_TABULAR_SIZE; i++) {
884 if(integralLine->at(i).trajectory.size() > 0) {
885 intervalSize = static_cast<ttk::SimplexId>(
886 integralLine->at(i).trajectory.size());
887 outputVertexNumber += intervalSize;
888 outputCellNumber += intervalSize - 1;
889 if(triangulation->getVertexRank(integralLine->at(i).trajectory.at(0))
890 != ttk::MPIrank_) {
891 intervalSize--;
892 }
893 if(integralLine->at(i).trajectory.size() > 1) {
894 realCellNumber += intervalSize - 1;
895 if(triangulation->getVertexRank(
896 integralLine->at(i).trajectory.back())
897 != ttk::MPIrank_) {
898 intervalSize--;
899 }
900 }
901 realVertexNumber += intervalSize;
902 } else {
903 break;
904 }
905 }
906 integralLine++;
907 }
908 }
909
910 ttk::SimplexId vertIndex;
911 ttk::SimplexId cellIndex;
912
913 // Perform exclusive prefix sum to find local offset for vertices and
914 // cells
915 MPI_Exscan(&realVertexNumber, &vertIndex, 1, ttk::getMPIType(vertIndex),
916 MPI_SUM, ttk::MPIcomm_);
917 MPI_Exscan(&realCellNumber, &cellIndex, 1, ttk::getMPIType(cellIndex),
918 MPI_SUM, ttk::MPIcomm_);
919
920 // Rank 0 received garbage values, it is replaced by the correct offset
921 // (always 0)
922 if(ttk::MPIrank_ == 0) {
923 vertIndex = 0;
924 cellIndex = 0;
925 }
926
927 // Generate Global identifers and package ghost data for the process that owns
928 // the ghost data
929 ttk::SimplexId startCellId = 0;
930 ttk::SimplexId startVertexId = 0;
931 globalVertexId.resize(outputVertexNumber);
932 globalCellId.resize(outputCellNumber);
933 std::vector<std::vector<ttk::intgl::GhostElementsToSort>> unmatchedGhosts(
934 neighborNumber_);
935 std::vector<std::vector<ttk::SimplexId>> unmatchedGhostsVertexLocalId(
936 neighborNumber_);
937 std::vector<std::vector<ttk::SimplexId>> unmatchedGhostsEdgeLocalId(
938 neighborNumber_);
939 for(int thread = 0; thread < threadNumber_; thread++) {
940 std::list<std::array<ttk::intgl::IntegralLine,
941 INTEGRAL_LINE_TABULAR_SIZE>>::const_iterator
942 integralLine
943 = integralLines[thread].list_.begin();
944 while(integralLine != integralLines[thread].list_.end()) {
945 for(int i = 0; i < INTEGRAL_LINE_TABULAR_SIZE; i++) {
946 if(integralLine->at(i).trajectory.size() > 0) {
947 if(triangulation->getVertexRank(integralLine->at(i).trajectory.at(0))
948 != ttk::MPIrank_) {
949 globalVertexId.at(startVertexId) = -1;
950 } else {
951 globalVertexId.at(startVertexId) = vertIndex;
952 vertIndex++;
953 }
954 startVertexId++;
955 if(integralLine->at(i).trajectory.size() > 1) {
956 if(triangulation->getVertexRank(
957 integralLine->at(i).trajectory.at(0))
958 != ttk::MPIrank_) {
959 globalCellId.at(startCellId) = -1;
960 unmatchedGhosts
961 .at(neighborsToId_[triangulation->getVertexRank(
962 integralLine->at(i).trajectory.at(0))])
963 .push_back(ttk::intgl::GhostElementsToSort{
964 vertIndex, integralLine->at(i).localVertexIdentifier.at(0),
965 integralLine->at(i).seedIdentifier,
966 integralLine->at(i).forkIdentifier, -1, startVertexId - 1,
967 startCellId});
968 } else {
969 globalCellId.at(startCellId) = cellIndex;
970 cellIndex++;
971 }
972 startCellId++;
973 if(integralLine->at(i).trajectory.size() > 2) {
974 std::iota(globalCellId.begin() + startCellId,
975 globalCellId.begin() + startCellId
976 + integralLine->at(i).trajectory.size() - 2,
977 cellIndex);
978 std::iota(globalVertexId.begin() + startVertexId,
979 globalVertexId.begin() + startVertexId
980 + integralLine->at(i).trajectory.size() - 2,
981 vertIndex);
982 startCellId += integralLine->at(i).trajectory.size() - 2;
983 startVertexId += integralLine->at(i).trajectory.size() - 2;
984 vertIndex += integralLine->at(i).trajectory.size() - 2;
985 cellIndex += integralLine->at(i).trajectory.size() - 2;
986 }
987 if(triangulation->getVertexRank(
988 integralLine->at(i).trajectory.back())
989 != ttk::MPIrank_) {
990 globalVertexId.at(startVertexId) = -1;
991 unmatchedGhosts
992 .at(neighborsToId_[triangulation->getVertexRank(
993 integralLine->at(i).trajectory.back())])
994 .push_back(ttk::intgl::GhostElementsToSort{
995 globalVertexId.at(startVertexId - 1),
996 integralLine->at(i).localVertexIdentifier.at(
997 integralLine->at(i).localVertexIdentifier.size() - 2),
998 integralLine->at(i).seedIdentifier,
999 integralLine->at(i).forkIdentifier,
1000 globalCellId.at(startCellId - 1), startVertexId,
1001 startCellId - 1});
1002 } else {
1003 globalVertexId.at(startVertexId) = vertIndex;
1004 vertIndex++;
1005 }
1006 startVertexId++;
1007 }
1008 } else {
1009 break;
1010 }
1011 }
1012 integralLine++;
1013 }
1014 }
1015 // unmatchedGhosts contains the ghost data for each process.
1016 // For two processes i and j that are neighbors, their respective
1017 // unmatchedGhosts[k_i] and unmatchedGhosts[k_j] have the same
1018 // number of elements and each element represents a segment of
1019 // integral line. This means that if unmatchedGhosts[k_i] and
1020 // unmatchedGhosts[k_j] are sorted using the same comparator,
1021 // element l of unmatchedGhosts[k_i] and element l of unmatchedGhosts[k_j]
1022 // represent the same segment of integral line. Therefore, once
1023 // unmatchedGhosts vectors are sorted for each process, all that
1024 // is left to do is exchange the global identifiers of vertices
1025 // and edges, as the receiving process already knows to which local
1026 // vertices in its domain the ghosts corresponds to.
1027 if(ttk::MPIsize_ > 1) {
1028 for(int i = 0; i < neighborNumber_; i++) {
1029 TTK_PSORT(threadNumber_, unmatchedGhosts.at(i).begin(),
1030 unmatchedGhosts.at(i).end());
1031 }
1032 this->exchangeGhosts(unmatchedGhosts, globalVertexId, globalCellId);
1033 }
1034
1035 return 0;
1036}
1037
1038inline int ttk::IntegralLines::exchangeGhosts(
1039 const std::vector<std::vector<ttk::intgl::GhostElementsToSort>>
1040 &unmatchedGhosts,
1041 std::vector<ttk::SimplexId> &globalVertexId,
1042 std::vector<ttk::SimplexId> &globalCellId) {
1043 ttk::SimplexId id{0};
1044 std::vector<std::vector<ttk::SimplexId>> globalIdsToSend(neighborNumber_);
1045 std::vector<std::vector<ttk::SimplexId>> globalIdsToReceive(neighborNumber_);
1046 for(int i = 0; i < neighborNumber_; i++) {
1047 globalIdsToSend[i].resize(2 * unmatchedGhosts[i].size());
1048 globalIdsToReceive[i].resize(2 * unmatchedGhosts[i].size());
1049 }
1050 // The sending buffer is prepared
1051#ifdef TTK_ENABLE_OPENMP4
1052#pragma omp parallel for
1053#endif
1054 for(int i = 0; i < neighborNumber_; i++) {
1055 for(size_t j = 0; j < unmatchedGhosts[i].size(); j++) {
1056 globalIdsToSend[i][2 * j] = unmatchedGhosts[i][j].ownedGlobalId;
1057 globalIdsToSend[i][2 * j + 1] = unmatchedGhosts[i][j].globalEdgeId;
1058 }
1059 }
1060 // Each process sends and receives global identifiers of ghosts
1061 // to its neighbors
1062 for(int neighbor : neighbors_) {
1063 MPI_Sendrecv(globalIdsToSend[neighborsToId_[neighbor]].data(),
1064 globalIdsToSend[neighborsToId_[neighbor]].size(),
1065 ttk::getMPIType(id), neighbor, ttk::MPIrank_,
1066 globalIdsToReceive[neighborsToId_[neighbor]].data(),
1067 globalIdsToSend[neighborsToId_[neighbor]].size(),
1068 ttk::getMPIType(id), neighbor, neighbor, ttk::MPIcomm_,
1069 MPI_STATUS_IGNORE);
1070 }
1071 // The global identifiers of ghosts are inserted in the globalVertexId
1072 // and globalCellId vectors
1073#ifdef TTK_ENABLE_OPENMP4
1074#pragma omp parallel for
1075#endif
1076 for(int i = 0; i < neighborNumber_; i++) {
1077 for(size_t j = 0; j < unmatchedGhosts[i].size(); j++) {
1078 globalVertexId.at(unmatchedGhosts[i][j].ghostVertexLocalId)
1079 = globalIdsToReceive[i][2 * j];
1080 if(globalCellId.at(unmatchedGhosts[i][j].ghostEdgeLocalId) == -1) {
1081 globalCellId.at(unmatchedGhosts[i][j].ghostEdgeLocalId)
1082 = globalIdsToReceive[i][j * 2 + 1];
1083 }
1084 }
1085 }
1086 return 0;
1087}
1088
1089#endif // TTK_ENABLE_MPI
ttkNotUsed
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition BaseClass.h:47
SimplexId
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
INTEGRAL_LINE_TABULAR_SIZE
#define INTEGRAL_LINE_TABULAR_SIZE
Definition IntegralLines.h:36
TTK_PSORT
#define TTK_PSORT(NTHREADS,...)
Parallel sort macro.
Definition OpenMP.h:46
ttk::AbstractTriangulation
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
Definition AbstractTriangulation.h:60
ttk::ArrayLinkedList
This class describes a dynamic size data structure for thread safe computation. It is a linked list o...
Definition ArrayLinkedList.h:21
ttk::IntegralLines::vertexIdentifierScalarField_
std::vector< ttk::SimplexId > * vertexIdentifierScalarField_
Definition IntegralLines.h:378
ttk::IntegralLines::scalarFieldCriticalPoints_
ttk::ScalarFieldCriticalPoints scalarFieldCriticalPoints_
Definition IntegralLines.h:382
ttk::IntegralLines::buildScalarFieldCriticalPoints
void buildScalarFieldCriticalPoints()
Definition IntegralLines.h:367
ttk::IntegralLines::setInputOffsets
void setInputOffsets(const SimplexId *const data)
Definition IntegralLines.h:347
ttk::IntegralLines::chunkSize_
ttk::SimplexId chunkSize_
Definition IntegralLines.h:374
ttk::IntegralLines::inputOffsets_
const ttk::SimplexId * inputOffsets_
Definition IntegralLines.h:377
ttk::IntegralLines::createTask
void createTask(const triangulationType *triangulation, std::vector< ttk::intgl::IntegralLine * > &chunkIntegralLine, const ttk::SimplexId *offsets, int nbElement) const
Create an OpenMP task that contains the computation of nbElement integral lines.
Definition IntegralLines.h:682
ttk::IntegralLines::setDirection
void setDirection(int direction)
Definition IntegralLines.h:187
ttk::IntegralLines::direction_
ttk::SimplexId direction_
Definition IntegralLines.h:375
ttk::IntegralLines::setInputScalarField
void setInputScalarField(void *data)
Definition IntegralLines.h:335
ttk::IntegralLines::setChunkSize
void setChunkSize(int size)
Definition IntegralLines.h:363
ttk::IntegralLines::computeIntegralLine
void computeIntegralLine(const triangulationType *triangulation, ttk::intgl::IntegralLine *integralLine, const ttk::SimplexId *offsets) const
Computes the integral line starting at the vertex of global id seedIdentifier.
Definition IntegralLines.h:509
ttk::IntegralLines::findNextVertex
void findNextVertex(ttk::SimplexId &vnext, ttk::SimplexId &fnext, std::vector< ttk::SimplexId > &component, const SimplexId *offsets) const
Finds the vertex of highest or lowest offsets (depending on the direction of the integral line) in th...
Definition IntegralLines.h:488
ttk::IntegralLines::preconditionTriangulation
int preconditionTriangulation(ttk::AbstractTriangulation *triangulation)
Definition IntegralLines.h:326
ttk::IntegralLines::setVertexNumber
void setVertexNumber(const SimplexId &vertexNumber)
Definition IntegralLines.h:179
ttk::IntegralLines::seedNumber_
ttk::SimplexId seedNumber_
Definition IntegralLines.h:373
ttk::IntegralLines::execute
int execute(triangulationType *triangulation)
Definition IntegralLines.h:701
ttk::IntegralLines::setOutputIntegralLines
void setOutputIntegralLines(std::vector< ttk::ArrayLinkedList< ttk::intgl::IntegralLine, INTEGRAL_LINE_TABULAR_SIZE > > *integralLines)
Definition IntegralLines.h:356
ttk::IntegralLines::setVertexIdentifierScalarField
void setVertexIdentifierScalarField(std::vector< SimplexId > *const data)
Definition IntegralLines.h:352
ttk::IntegralLines::outputIntegralLines_
std::vector< ttk::ArrayLinkedList< ttk::intgl::IntegralLine, INTEGRAL_LINE_TABULAR_SIZE > > * outputIntegralLines_
Definition IntegralLines.h:381
ttk::IntegralLines::~IntegralLines
~IntegralLines() override
ttk::IntegralLines::vertexNumber_
ttk::SimplexId vertexNumber_
Definition IntegralLines.h:372
ttk::IntegralLines::prepareForTask
void prepareForTask(const triangulationType *triangulation, std::vector< ttk::intgl::IntegralLine * > &chunkIntegralLine, int startingIndex, int nbElement, std::vector< SimplexId > *seeds) const
Initializes the three attributes of an integral line: the global id of its seed, its trajectory,...
ttk::IntegralLines::setSeedNumber
void setSeedNumber(const SimplexId &seedNumber)
Definition IntegralLines.h:183
ttk::ScalarFieldCriticalPoints
TTK processing package for the computation of critical points in PL scalar fields defined on PL manif...
Definition ScalarFieldCriticalPoints.h:58
ttk::Timer::getElapsedTime
double getElapsedTime()
Definition Timer.h:15
ttk::Geometry::distance
T distance(const T *p0, const T *p1, const int &dimension=3)
Definition Geometry.cpp:362
ttk::rpd::operator<
bool operator<(const FiltratedEdge &e1, const FiltratedEdge &e2)
Definition RipsPersistenceDiagramUtils.h:58
ttk
TTK base package defining the standard types.
Definition AbstractTriangulation.h:51
ttk::MPIsize_
COMMON_EXPORTS int MPIsize_
Definition BaseClass.cpp:10
ttk::SimplexId
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
ttk::MPIrank_
COMMON_EXPORTS int MPIrank_
Definition BaseClass.cpp:9
end
T end(std::pair< T, T > &p)
Definition ripser.cpp:503
ttk::intgl::IntegralLine
Struct containing the data of an integral line. trajectories: vector of identifiers of each vertex th...
Definition IntegralLines.h:54
ttk::intgl::IntegralLine::localVertexIdentifier
std::vector< ttk::SimplexId > localVertexIdentifier
Definition IntegralLines.h:57
ttk::intgl::IntegralLine::forkIdentifier
ttk::SimplexId forkIdentifier
Definition IntegralLines.h:59
ttk::intgl::IntegralLine::distanceFromSeed
std::vector< double > distanceFromSeed
Definition IntegralLines.h:56
ttk::intgl::IntegralLine::seedIdentifier
ttk::SimplexId seedIdentifier
Definition IntegralLines.h:58
ttk::intgl::IntegralLine::trajectory
std::vector< ttk::SimplexId > trajectory
Definition IntegralLines.h:55
printMsg
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/| (_) |"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)
Generated on Fri Apr 10 2026 10:02:07 for TTK by doxygen 1.13.2