TTK
Loading...
Searching...
No Matches
ttkIntegralLines.cpp
Go to the documentation of this file.
1#include "Triangulation.h"
2#include <string>
3#include <ttkIntegralLines.h>
4#include <ttkMacros.h>
5#include <ttkUtils.h>
6
7#include <ArrayLinkedList.h>
8#include <vtkCellData.h>
9#include <vtkDataArray.h>
10#include <vtkDataObject.h>
11#include <vtkDataSet.h>
12#include <vtkDoubleArray.h>
13#include <vtkInformation.h>
14#include <vtkObjectFactory.h>
15#include <vtkPointData.h>
16#include <vtkPointSet.h>
17#include <vtkUnstructuredGrid.h>
18
19#include <array>
20
22
24 this->SetNumberOfInputPorts(2);
25 this->SetNumberOfOutputPorts(1);
26}
27
29
30int ttkIntegralLines::FillInputPortInformation(int port, vtkInformation *info) {
31 if(port == 0)
32 info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkDataSet");
33 if(port == 1)
34 info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPointSet");
35
36 return 1;
37}
38
40 vtkInformation *info) {
41 if(port == 0)
42 info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid");
43
44 return 1;
45}
46
47template <typename triangulationType>
49 vtkDataSet *input,
50 const triangulationType *triangulation,
51 const std::vector<
53 &integralLines,
54#ifdef TTK_ENABLE_MPI
55 const std::vector<ttk::SimplexId> &globalVertexId,
56 const std::vector<ttk::SimplexId> &globalCellId,
57#endif
58 vtkUnstructuredGrid *output) {
59 if(input == nullptr || output == nullptr
60 || input->GetPointData() == nullptr) {
61 this->printErr("Null pointers in getTrajectories parameters");
62 return 0;
63 }
64
65 vtkNew<vtkUnstructuredGrid> ug{};
66 vtkNew<vtkPoints> pts{};
67 vtkNew<vtkDoubleArray> dist{};
68 vtkNew<vtkIdTypeArray> identifier{};
69#ifdef TTK_ENABLE_MPI
70 vtkNew<vtkIdTypeArray> vtkEdgeIdentifiers{};
71 vtkNew<vtkIdTypeArray> vtkVertexGlobalIdArray{};
72 vtkNew<vtkIntArray> vtkVertexRankArray{};
73 vtkNew<vtkIntArray> vtkEdgeRankArray{};
74#endif
75 vtkNew<vtkIdTypeArray> vtkForkIdentifiers{};
76 vtkNew<vtkUnsignedCharArray> outputMaskField{};
77
78 outputMaskField->SetNumberOfComponents(1);
79 outputMaskField->SetName(ttk::MaskScalarFieldName);
80
81 dist->SetNumberOfComponents(1);
82 dist->SetName("DistanceFromSeed");
83 identifier->SetNumberOfComponents(1);
84 identifier->SetName("SeedIdentifier");
85 vtkForkIdentifiers->SetNumberOfComponents(1);
86 vtkForkIdentifiers->SetName("ForkIdentifiers");
87
88#ifdef TTK_ENABLE_MPI
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");
97#endif
98 const auto numberOfArrays = input->GetPointData()->GetNumberOfArrays();
99
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) {
105 // only keep scalar arrays
106 scalarArrays.push_back(a);
107 }
108 }
109
110 std::vector<vtkSmartPointer<vtkDataArray>> inputScalars(scalarArrays.size());
111 for(size_t k = 0; k < scalarArrays.size(); ++k) {
112 inputScalars[k]
113 = vtkSmartPointer<vtkDataArray>::Take(scalarArrays[k]->NewInstance());
114 inputScalars[k]->SetNumberOfComponents(1);
115 inputScalars[k]->SetName(scalarArrays[k]->GetName());
116 }
117 std::array<float, 3> p;
118 std::array<vtkIdType, 2> ids;
119#ifdef TTK_ENABLE_MPI
120 ttk::SimplexId vertexCounter = 0;
121 ttk::SimplexId edgeCounter = 0;
122#endif
123 for(int thread = 0; thread < threadNumber_; thread++) {
124 auto integralLine = integralLines[thread].list_.begin();
125 while(integralLine != integralLines[thread].list_.end()) {
126 for(int i = 0; i < INTEGRAL_LINE_TABULAR_SIZE; i++) {
127 if(integralLine->at(i).trajectory.size() > 0) {
128 ttk::SimplexId vertex = integralLine->at(i).trajectory.at(0);
129 triangulation->getVertexPoint(vertex, p[0], p[1], p[2]);
130 ids[0] = pts->InsertNextPoint(p.data());
131 // distanceScalars
132 dist->InsertNextTuple1(integralLine->at(i).distanceFromSeed.at(0));
133#ifdef TTK_ENABLE_MPI
134 outputMaskField->InsertNextTuple1(0);
135 vtkVertexGlobalIdArray->InsertNextTuple1(
136 globalVertexId.at(vertexCounter));
137 vertexCounter++;
138 vtkVertexRankArray->InsertNextTuple1(
139 triangulation->getVertexRank(vertex));
140#else
141 outputMaskField->InsertNextTuple1(0);
142#endif
143 identifier->InsertNextTuple1(integralLine->at(i).seedIdentifier);
144 vtkForkIdentifiers->InsertNextTuple1(
145 integralLine->at(i).forkIdentifier);
146 // inputScalars
147 for(size_t k = 0; k < scalarArrays.size(); ++k) {
148 inputScalars[k]->InsertNextTuple1(
149 scalarArrays[k]->GetTuple1(vertex));
150 }
151 for(size_t j = 1; j < integralLine->at(i).trajectory.size(); ++j) {
152 vertex = integralLine->at(i).trajectory.at(j);
153#ifdef TTK_ENABLE_MPI
154 vtkVertexGlobalIdArray->InsertNextTuple1(
155 globalVertexId.at(vertexCounter));
156 vertexCounter++;
157 vtkEdgeIdentifiers->InsertNextTuple1(globalCellId.at(edgeCounter));
158 edgeCounter++;
159 vtkEdgeRankArray->InsertNextTuple1(triangulation->getVertexRank(
160 integralLine->at(i).trajectory.at(j - 1)));
161 vtkVertexRankArray->InsertNextTuple1(
162 triangulation->getVertexRank(vertex));
163#endif
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());
169 // distanceScalars
170 dist->InsertNextTuple1(integralLine->at(i).distanceFromSeed.at(j));
171 identifier->InsertNextTuple1(integralLine->at(i).seedIdentifier);
172 // inputScalars
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());
177 // iteration
178 ids[0] = ids[1];
179 }
180 outputMaskField->SetTuple1(
181 outputMaskField->GetNumberOfTuples() - 1, 0);
182 } else {
183 break;
184 }
185 }
186 integralLine++;
187 }
188 }
189
190 ug->SetPoints(pts);
191 ug->GetPointData()->AddArray(dist);
192 ug->GetPointData()->AddArray(identifier);
193 ug->GetPointData()->AddArray(vtkForkIdentifiers);
194 ug->GetPointData()->AddArray(outputMaskField);
195
196 for(unsigned int k = 0; k < scalarArrays.size(); ++k) {
197 ug->GetPointData()->AddArray(inputScalars[k]);
198 }
199#ifdef TTK_ENABLE_MPI
200 ug->GetPointData()->AddArray(vtkVertexRankArray);
201 ug->GetCellData()->AddArray(vtkEdgeRankArray);
202 ug->GetCellData()->SetGlobalIds(vtkEdgeIdentifiers);
203 ug->GetPointData()->SetGlobalIds(vtkVertexGlobalIdArray);
204#endif
205 output->ShallowCopy(ug);
206
207 return 1;
208}
209
210int ttkIntegralLines::RequestData(vtkInformation *ttkNotUsed(request),
211 vtkInformationVector **inputVector,
212 vtkInformationVector *outputVector) {
213
214 vtkDataSet *domain = vtkDataSet::GetData(inputVector[0], 0);
215 vtkPointSet *seeds = vtkPointSet::GetData(inputVector[1], 0);
216 vtkUnstructuredGrid *output = vtkUnstructuredGrid::GetData(outputVector, 0);
217
218 ttk::Triangulation *triangulation = ttkAlgorithm::GetTriangulation(domain);
219 vtkDataArray *inputScalars = this->GetInputArrayToProcess(0, domain);
220
221 int const keepGoing = checkEmptyMPIInput<ttk::Triangulation>(triangulation);
222 if(keepGoing < 2) {
223 return keepGoing;
224 }
225 vtkDataArray *inputOffsets = this->GetOrderArray(
226 domain, 0, triangulation, false, 1, ForceInputOffsetScalarField);
227
228 const ttk::SimplexId numberOfPointsInDomain = domain->GetNumberOfPoints();
229 this->setVertexNumber(numberOfPointsInDomain);
230 int numberOfPointsInSeeds = seeds->GetNumberOfPoints();
231#ifndef TTK_ENABLE_KAMIKAZE
232 int totalSeeds;
233#else
234#ifdef TTK_ENABLE_MPI
235 int totalSeeds;
236#endif
237#endif
238
239#ifdef TTK_ENABLE_MPI_TIME
240 ttk::Timer t_mpi;
241 ttk::startMPITimer(t_mpi, ttk::MPIrank_, ttk::MPIsize_);
242#endif
243#ifdef TTK_ENABLE_MPI
244 // Necessary when using MPI
245 std::vector<ttk::SimplexId> inputIdentifiers{};
246 if(ttk::MPIsize_ > 1) {
247 MPI_Reduce(&numberOfPointsInSeeds, &totalSeeds, 1, MPI_INTEGER, MPI_SUM, 0,
248 ttk::MPIcomm_);
249 int isDistributed;
250
251 if(ttk::MPIrank_ == 0) {
252 isDistributed = numberOfPointsInSeeds != totalSeeds;
253 }
254 MPI_Bcast(&isDistributed, 1, MPI_INTEGER, 0, ttk::MPIcomm_);
255 MPI_Bcast(&totalSeeds, 1, MPI_INTEGER, 0, ttk::MPIcomm_);
256
257 if(!isDistributed) {
258 this->setGlobalElementCounter(totalSeeds);
259 vtkDataArray *globalSeedsId;
260 if(ttk::MPIrank_ == 0) {
261 globalSeedsId = seeds->GetPointData()->GetArray("GlobalPointIds");
262 } else {
263 globalSeedsId = vtkDataArray::CreateDataArray(VTK_ID_TYPE);
264 }
265
266 if(ttk::MPIrank_ != 0) {
267 globalSeedsId->SetNumberOfComponents(1);
268 globalSeedsId->SetNumberOfTuples(totalSeeds);
269 }
270 ttk::LongSimplexId id = 0;
271 MPI_Bcast(ttkUtils::GetPointer<ttk::LongSimplexId>(globalSeedsId),
272 totalSeeds, ttk::getMPIType(id), 0, ttk::MPIcomm_);
273 ttk::SimplexId localId = -1;
274 for(int i = 0; i < totalSeeds; i++) {
275 localId = triangulation->getVertexLocalId(globalSeedsId->GetTuple1(i));
276
277 if(localId != -1
278 && triangulation->getVertexRank(localId) == ttk::MPIrank_) {
279 inputIdentifiers.push_back(localId);
280 }
281 }
282 numberOfPointsInSeeds = inputIdentifiers.size();
283 } else {
284 std::vector<ttk::SimplexId> idSpareStorage{};
285 ttk::SimplexId *inputIdentifierGlobalId;
286 inputIdentifierGlobalId = this->GetIdentifierArrayPtr(
287 ForceInputVertexScalarField, 2, ttk::VertexScalarFieldName, seeds,
288 idSpareStorage);
289 ttk::SimplexId localId = 0;
290 for(int i = 0; i < numberOfPointsInSeeds; i++) {
291 localId = triangulation->getVertexLocalId(inputIdentifierGlobalId[i]);
292 if(localId != -1) {
293 inputIdentifiers.push_back(localId);
294 }
295 }
296 numberOfPointsInSeeds = inputIdentifiers.size();
297 MPI_Allreduce(&numberOfPointsInSeeds, &totalSeeds, 1, MPI_INTEGER,
298 MPI_SUM, ttk::MPIcomm_);
299 this->setGlobalElementCounter(totalSeeds);
300 }
301 } else {
302 this->setGlobalElementCounter(numberOfPointsInSeeds);
303 inputIdentifiers.resize(numberOfPointsInSeeds);
304 totalSeeds = numberOfPointsInSeeds;
305 std::vector<ttk::SimplexId> idSpareStorage{};
306 ttk::SimplexId *inputIdentifierGlobalId;
307 inputIdentifierGlobalId = this->GetIdentifierArrayPtr(
308 ForceInputVertexScalarField, 2, ttk::VertexScalarFieldName, seeds,
309 idSpareStorage);
310 for(int i = 0; i < numberOfPointsInSeeds; i++) {
311 inputIdentifiers.at(i)
312 = triangulation->getVertexLocalId(inputIdentifierGlobalId[i]);
313 }
314 }
315#else
316 std::vector<ttk::SimplexId> idSpareStorage{};
317 ttk::SimplexId *identifiers = this->GetIdentifierArrayPtr(
318 ForceInputVertexScalarField, 2, ttk::VertexScalarFieldName, seeds,
319 idSpareStorage);
320 std::unordered_set<ttk::SimplexId> isSeed;
321 for(ttk::SimplexId k = 0; k < numberOfPointsInSeeds; ++k) {
322 isSeed.insert(identifiers[k]);
323 }
324 std::vector<ttk::SimplexId> inputIdentifiers(isSeed.begin(), isSeed.end());
325#ifndef TTK_ENABLE_KAMIKAZE
326 totalSeeds = inputIdentifiers.size();
327#endif
328 isSeed.clear();
329#endif
330
331 std::vector<
333 integralLines(
336
337 this->setVertexNumber(numberOfPointsInDomain);
338 this->setSeedNumber(numberOfPointsInSeeds);
339 this->setDirection(Direction);
340 this->setInputScalarField(inputScalars->GetVoidPointer(0));
341 this->setInputOffsets(ttkUtils::GetPointer<ttk::SimplexId>(inputOffsets));
342 this->setVertexIdentifierScalarField(&inputIdentifiers);
343 this->setOutputIntegralLines(&integralLines);
344 this->preconditionTriangulation(triangulation);
345 this->setChunkSize(
346 std::max(std::max(std::min(1000, (int)numberOfPointsInSeeds),
347 (int)numberOfPointsInSeeds / (threadNumber_ * 100)),
348 1));
349#ifdef TTK_ENABLE_MPI
350 std::vector<std::vector<std::vector<ttk::intgl::ElementToBeSent>>> toSend(
352 this->setNeighbors(triangulation->getNeighborRanks());
353 if(ttk::MPIsize_ > 1) {
354 toSend.resize(this->neighborNumber_);
355 for(int i = 0; i < this->neighborNumber_; i++) {
356 toSend[i].resize(this->threadNumber_);
357 for(int j = 0; j < this->threadNumber_; j++) {
358 toSend[i][j].reserve((int)numberOfPointsInSeeds * 0.005
359 / this->threadNumber_);
360 }
361 }
362 }
363 this->setToSend(&toSend);
364 this->createMessageType();
365#endif
366#ifdef TTK_ENABLE_MPI_TIME
367 double elapsedTime = ttk::endMPITimer(t_mpi, ttk::MPIrank_, ttk::MPIsize_);
368 if(ttk::MPIrank_ == 0) {
369 printMsg("Preparation performed using " + std::to_string(ttk::MPIsize_)
370 + " MPI processes lasted :" + std::to_string(elapsedTime));
371 }
372#endif
373#ifndef TTK_ENABLE_KAMIKAZE
374 // field problem
375 if(!inputScalars) {
376 this->printErr("wrong scalar field.");
377 return -1;
378 }
379 // field problem
380 if(inputOffsets->GetDataType() != VTK_INT
381 and inputOffsets->GetDataType() != VTK_ID_TYPE) {
382 this->printErr("input offset field type not supported.");
383 return -1;
384 }
385
386 // no points.
387 if(numberOfPointsInDomain <= 0) {
388 this->printErr("domain has no points.");
389 return -1;
390 }
391 // no points.
392 if(totalSeeds <= 0) {
393 this->printErr("seeds have no points.");
394 return -1;
395 }
396#endif
397#ifdef TTK_ENABLE_MPI_TIME
398 ttk::startMPITimer(t_mpi, ttk::MPIrank_, ttk::MPIsize_);
399#endif
400 int status = 0;
401 ttkTemplateMacro(triangulation->getType(),
402 (status = this->execute<TTK_TT>(
403 static_cast<TTK_TT *>(triangulation->getData()))));
404
405#ifndef TTK_ENABLE_KAMIKAZE
406 // something wrong in baseCode
407 if(status != 0) {
408 this->printErr("IntegralLines.execute() error code : "
409 + std::to_string(status));
410 return 0;
411 }
412#endif
413#ifdef TTK_ENABLE_MPI
414 std::vector<ttk::SimplexId> globalVertexId;
415 std::vector<ttk::SimplexId> globalCellId;
416 ttkTemplateMacro(triangulation->getType(),
417 (getGlobalIdentifiers<TTK_TT>(
418 globalVertexId, globalCellId, integralLines,
419 static_cast<TTK_TT *>(triangulation->getData()))));
420#ifdef TTK_ENABLE_MPI_TIME
421 elapsedTime = ttk::endMPITimer(t_mpi, ttk::MPIrank_, ttk::MPIsize_);
422 if(ttk::MPIrank_ == 0) {
423 printMsg("Computation performed using " + std::to_string(ttk::MPIsize_)
424 + " MPI processes lasted :" + std::to_string(elapsedTime));
425 }
426#endif
427#endif
428 // make the vtk trajectories
429#ifdef TTK_ENABLE_MPI
430 ttkTemplateMacro(triangulation->getType(),
431 (getTrajectories<TTK_TT>(
432 domain, static_cast<TTK_TT *>(triangulation->getData()),
433 integralLines, globalVertexId, globalCellId, output)));
434#else
435 ttkTemplateMacro(triangulation->getType(),
436 (getTrajectories<TTK_TT>(
437 domain, static_cast<TTK_TT *>(triangulation->getData()),
438 integralLines, output)));
439#endif
440
441 return (int)(status == 0);
442}
#define ttkTemplateMacro(triangulationType, call)
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition BaseClass.h:47
#define INTEGRAL_LINE_TABULAR_SIZE
ttk::SimplexId * GetIdentifierArrayPtr(const bool &enforceArrayIndex, const int &arrayIndex, const std::string &arrayName, vtkDataSet *const inputData, std::vector< ttk::SimplexId > &spareStorage, const int inputPort=0, const bool printErr=true)
ttk::Triangulation * GetTriangulation(vtkDataSet *dataSet)
vtkDataArray * GetOrderArray(vtkDataSet *const inputData, const int scalarArrayIdx, ttk::Triangulation *triangulation, const bool getGlobalOrder=false, const int orderArrayIdx=0, const bool enforceOrderArrayIdx=false)
TTK VTK-filter for the computation of edge-based integral lines of the gradient of an input scalar fi...
int getTrajectories(vtkDataSet *input, const triangulationType *triangulation, const std::vector< ttk::ArrayLinkedList< ttk::intgl::IntegralLine, INTEGRAL_LINE_TABULAR_SIZE > > &integralLines, vtkUnstructuredGrid *output)
int FillOutputPortInformation(int port, vtkInformation *info) override
int RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) override
int FillInputPortInformation(int port, vtkInformation *info) override
~ttkIntegralLines() override
This class describes a dynamic size data structure for thread safe computation. It is a linked list o...
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:149
void setInputOffsets(const SimplexId *const data)
void setDirection(int direction)
void setInputScalarField(void *data)
void setChunkSize(int size)
int preconditionTriangulation(ttk::AbstractTriangulation *triangulation)
void setVertexNumber(const SimplexId &vertexNumber)
void setOutputIntegralLines(std::vector< ttk::ArrayLinkedList< ttk::intgl::IntegralLine, INTEGRAL_LINE_TABULAR_SIZE > > *integralLines)
void setVertexIdentifierScalarField(std::vector< SimplexId > *const data)
void setSeedNumber(const SimplexId &seedNumber)
Triangulation is a class that provides time and memory efficient traversal methods on triangulations ...
AbstractTriangulation * getData()
Triangulation::Type getType() const
COMMON_EXPORTS int MPIsize_
Definition BaseClass.cpp:10
COMMON_EXPORTS int MPIrank_
Definition BaseClass.cpp:9
const char VertexScalarFieldName[]
default name for vertex scalar field
Definition DataTypes.h:35
long long int LongSimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:15
const char MaskScalarFieldName[]
default name for mask scalar field
Definition DataTypes.h:32
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
Struct containing the data of an integral line. trajectories: vector of identifiers of each vertex th...
vtkStandardNewMacro(ttkIntegralLines)
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)