50 const triangulationType *triangulation,
55 const std::vector<ttk::SimplexId> &globalVertexId,
56 const std::vector<ttk::SimplexId> &globalCellId,
58 vtkUnstructuredGrid *output) {
59 if(input ==
nullptr || output ==
nullptr
60 || input->GetPointData() ==
nullptr) {
61 this->
printErr(
"Null pointers in getTrajectories parameters");
65 vtkNew<vtkUnstructuredGrid> ug{};
66 vtkNew<vtkPoints> pts{};
67 vtkNew<vtkDoubleArray> dist{};
68 vtkNew<vtkIdTypeArray> identifier{};
70 vtkNew<vtkIdTypeArray> vtkEdgeIdentifiers{};
71 vtkNew<vtkIdTypeArray> vtkVertexGlobalIdArray{};
72 vtkNew<vtkIntArray> vtkVertexRankArray{};
73 vtkNew<vtkIntArray> vtkEdgeRankArray{};
75 vtkNew<vtkIdTypeArray> vtkForkIdentifiers{};
76 vtkNew<vtkUnsignedCharArray> outputMaskField{};
78 outputMaskField->SetNumberOfComponents(1);
81 dist->SetNumberOfComponents(1);
82 dist->SetName(
"DistanceFromSeed");
83 identifier->SetNumberOfComponents(1);
84 identifier->SetName(
"SeedIdentifier");
85 vtkForkIdentifiers->SetNumberOfComponents(1);
86 vtkForkIdentifiers->SetName(
"ForkIdentifiers");
89 vtkVertexGlobalIdArray->SetNumberOfComponents(1);
90 vtkVertexGlobalIdArray->SetName(
"GlobalPointIds");
91 vtkEdgeIdentifiers->SetNumberOfComponents(1);
92 vtkEdgeIdentifiers->SetName(
"GlobalCellIds");
93 vtkEdgeRankArray->SetNumberOfComponents(1);
94 vtkEdgeRankArray->SetName(
"RankArray");
95 vtkVertexRankArray->SetNumberOfComponents(1);
96 vtkVertexRankArray->SetName(
"RankArray");
98 const auto numberOfArrays = input->GetPointData()->GetNumberOfArrays();
100 std::vector<vtkDataArray *> scalarArrays{};
101 scalarArrays.reserve(numberOfArrays);
102 for(
int k = 0; k < numberOfArrays; ++k) {
103 const auto a = input->GetPointData()->GetArray(k);
104 if(a->GetNumberOfComponents() == 1) {
106 scalarArrays.push_back(a);
110 std::vector<vtkSmartPointer<vtkDataArray>> inputScalars(scalarArrays.size());
111 for(
size_t k = 0; k < scalarArrays.size(); ++k) {
114 inputScalars[k]->SetNumberOfComponents(1);
115 inputScalars[k]->SetName(scalarArrays[k]->GetName());
117 std::array<float, 3> p;
118 std::array<vtkIdType, 2> ids;
124 auto integralLine = integralLines[thread].list_.begin();
125 while(integralLine != integralLines[thread].list_.end()) {
127 if(integralLine->at(i).trajectory.size() > 0) {
129 triangulation->getVertexPoint(vertex, p[0], p[1], p[2]);
130 ids[0] = pts->InsertNextPoint(p.data());
132 dist->InsertNextTuple1(integralLine->at(i).distanceFromSeed.at(0));
134 outputMaskField->InsertNextTuple1(0);
135 vtkVertexGlobalIdArray->InsertNextTuple1(
136 globalVertexId.at(vertexCounter));
138 vtkVertexRankArray->InsertNextTuple1(
139 triangulation->getVertexRank(vertex));
141 outputMaskField->InsertNextTuple1(0);
143 identifier->InsertNextTuple1(integralLine->at(i).seedIdentifier);
144 vtkForkIdentifiers->InsertNextTuple1(
145 integralLine->at(i).forkIdentifier);
147 for(
size_t k = 0; k < scalarArrays.size(); ++k) {
148 inputScalars[k]->InsertNextTuple1(
149 scalarArrays[k]->GetTuple1(vertex));
151 for(
size_t j = 1; j < integralLine->at(i).trajectory.size(); ++j) {
152 vertex = integralLine->at(i).trajectory.at(j);
154 vtkVertexGlobalIdArray->InsertNextTuple1(
155 globalVertexId.at(vertexCounter));
157 vtkEdgeIdentifiers->InsertNextTuple1(globalCellId.at(edgeCounter));
159 vtkEdgeRankArray->InsertNextTuple1(triangulation->getVertexRank(
160 integralLine->at(i).trajectory.at(j - 1)));
161 vtkVertexRankArray->InsertNextTuple1(
162 triangulation->getVertexRank(vertex));
164 outputMaskField->InsertNextTuple1(1);
165 vtkForkIdentifiers->InsertNextTuple1(
166 integralLine->at(i).forkIdentifier);
167 triangulation->getVertexPoint(vertex, p[0], p[1], p[2]);
168 ids[1] = pts->InsertNextPoint(p.data());
170 dist->InsertNextTuple1(integralLine->at(i).distanceFromSeed.at(j));
171 identifier->InsertNextTuple1(integralLine->at(i).seedIdentifier);
173 for(
unsigned int k = 0; k < scalarArrays.size(); ++k)
174 inputScalars[k]->InsertNextTuple1(
175 scalarArrays[k]->GetTuple1(vertex));
176 ug->InsertNextCell(VTK_LINE, 2, ids.data());
180 outputMaskField->SetTuple1(
181 outputMaskField->GetNumberOfTuples() - 1, 0);
191 ug->GetPointData()->AddArray(dist);
192 ug->GetPointData()->AddArray(identifier);
193 ug->GetPointData()->AddArray(vtkForkIdentifiers);
194 ug->GetPointData()->AddArray(outputMaskField);
196 for(
unsigned int k = 0; k < scalarArrays.size(); ++k) {
197 ug->GetPointData()->AddArray(inputScalars[k]);
200 ug->GetPointData()->AddArray(vtkVertexRankArray);
201 ug->GetCellData()->AddArray(vtkEdgeRankArray);
202 ug->GetCellData()->SetGlobalIds(vtkEdgeIdentifiers);
203 ug->GetPointData()->SetGlobalIds(vtkVertexGlobalIdArray);
205 output->ShallowCopy(ug);
211 vtkInformationVector **inputVector,
212 vtkInformationVector *outputVector) {
214 vtkDataSet *domain = vtkDataSet::GetData(inputVector[0], 0);
215 vtkPointSet *seeds = vtkPointSet::GetData(inputVector[1], 0);
216 vtkUnstructuredGrid *output = vtkUnstructuredGrid::GetData(outputVector, 0);
219 vtkDataArray *inputScalars = this->GetInputArrayToProcess(0, domain);
221 int const keepGoing = checkEmptyMPIInput<ttk::Triangulation>(triangulation);
226 domain, 0, triangulation,
false, 1, ForceInputOffsetScalarField);
228 const ttk::SimplexId numberOfPointsInDomain = domain->GetNumberOfPoints();
230 int numberOfPointsInSeeds = seeds->GetNumberOfPoints();
231#ifndef TTK_ENABLE_KAMIKAZE
239#ifdef TTK_ENABLE_MPI_TIME
245 std::vector<ttk::SimplexId> inputIdentifiers{};
247 MPI_Reduce(&numberOfPointsInSeeds, &totalSeeds, 1, MPI_INTEGER, MPI_SUM, 0,
252 isDistributed = numberOfPointsInSeeds != totalSeeds;
254 MPI_Bcast(&isDistributed, 1, MPI_INTEGER, 0, ttk::MPIcomm_);
255 MPI_Bcast(&totalSeeds, 1, MPI_INTEGER, 0, ttk::MPIcomm_);
258 this->setGlobalElementCounter(totalSeeds);
259 vtkDataArray *globalSeedsId;
261 globalSeedsId = seeds->GetPointData()->GetArray(
"GlobalPointIds");
263 globalSeedsId = vtkDataArray::CreateDataArray(VTK_ID_TYPE);
267 globalSeedsId->SetNumberOfComponents(1);
268 globalSeedsId->SetNumberOfTuples(totalSeeds);
271 MPI_Bcast(ttkUtils::GetPointer<ttk::LongSimplexId>(globalSeedsId),
272 totalSeeds, ttk::getMPIType(
id), 0, ttk::MPIcomm_);
274 for(
int i = 0; i < totalSeeds; i++) {
275 localId = triangulation->getVertexLocalId(globalSeedsId->GetTuple1(i));
279 inputIdentifiers.push_back(localId);
282 numberOfPointsInSeeds = inputIdentifiers.size();
284 std::vector<ttk::SimplexId> idSpareStorage{};
290 for(
int i = 0; i < numberOfPointsInSeeds; i++) {
291 localId = triangulation->getVertexLocalId(inputIdentifierGlobalId[i]);
293 inputIdentifiers.push_back(localId);
296 numberOfPointsInSeeds = inputIdentifiers.size();
297 MPI_Allreduce(&numberOfPointsInSeeds, &totalSeeds, 1, MPI_INTEGER,
298 MPI_SUM, ttk::MPIcomm_);
299 this->setGlobalElementCounter(totalSeeds);
302 this->setGlobalElementCounter(numberOfPointsInSeeds);
303 inputIdentifiers.resize(numberOfPointsInSeeds);
304 totalSeeds = numberOfPointsInSeeds;
305 std::vector<ttk::SimplexId> idSpareStorage{};
310 for(
int i = 0; i < numberOfPointsInSeeds; i++) {
311 inputIdentifiers.at(i)
312 = triangulation->getVertexLocalId(inputIdentifierGlobalId[i]);
316 std::vector<ttk::SimplexId> idSpareStorage{};
320 std::unordered_set<ttk::SimplexId> isSeed;
322 isSeed.insert(identifiers[k]);
324 std::vector<ttk::SimplexId> inputIdentifiers(isSeed.begin(), isSeed.end());
325#ifndef TTK_ENABLE_KAMIKAZE
326 totalSeeds = inputIdentifiers.size();
341 this->
setInputOffsets(ttkUtils::GetPointer<ttk::SimplexId>(inputOffsets));
346 std::max(std::max(std::min(1000, (
int)numberOfPointsInSeeds),
350 std::vector<std::vector<std::vector<ttk::intgl::ElementToBeSent>>> toSend(
352 this->setNeighbors(triangulation->getNeighborRanks());
354 toSend.resize(this->neighborNumber_);
355 for(
int i = 0; i < this->neighborNumber_; i++) {
358 toSend[i][j].reserve((
int)numberOfPointsInSeeds * 0.005
359 / this->threadNumber_);
363 this->setToSend(&toSend);
364 this->createMessageType();
366#ifdef TTK_ENABLE_MPI_TIME
370 +
" MPI processes lasted :" + std::to_string(elapsedTime));
373#ifndef TTK_ENABLE_KAMIKAZE
376 this->
printErr(
"wrong scalar field.");
380 if(inputOffsets->GetDataType() != VTK_INT
381 and inputOffsets->GetDataType() != VTK_ID_TYPE) {
382 this->
printErr(
"input offset field type not supported.");
387 if(numberOfPointsInDomain <= 0) {
388 this->
printErr(
"domain has no points.");
392 if(totalSeeds <= 0) {
393 this->
printErr(
"seeds have no points.");
397#ifdef TTK_ENABLE_MPI_TIME
402 (status = this->execute<TTK_TT>(
403 static_cast<TTK_TT *
>(triangulation->
getData()))));
405#ifndef TTK_ENABLE_KAMIKAZE
408 this->
printErr(
"IntegralLines.execute() error code : "
409 + std::to_string(status));
414 std::vector<ttk::SimplexId> globalVertexId;
415 std::vector<ttk::SimplexId> globalCellId;
417 (getGlobalIdentifiers<TTK_TT>(
418 globalVertexId, globalCellId, integralLines,
419 static_cast<TTK_TT *
>(triangulation->
getData()))));
420#ifdef TTK_ENABLE_MPI_TIME
424 +
" MPI processes lasted :" + std::to_string(elapsedTime));
431 (getTrajectories<TTK_TT>(
432 domain,
static_cast<TTK_TT *
>(triangulation->
getData()),
433 integralLines, globalVertexId, globalCellId, output)));
436 (getTrajectories<TTK_TT>(
437 domain,
static_cast<TTK_TT *
>(triangulation->
getData()),
438 integralLines, output)));
441 return (
int)(status == 0);