1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10 // Copyright (c) 2018, The Regents of the University of California, through
11 // Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
12 // from the U.S. Dept. of Energy). All rights reserved.
13 //
14 // Redistribution and use in source and binary forms, with or without modification,
15 // are permitted provided that the following conditions are met:
16 //
17 // (1) Redistributions of source code must retain the above copyright notice, this
18 // list of conditions and the following disclaimer.
19 //
20 // (2) Redistributions in binary form must reproduce the above copyright notice,
21 // this list of conditions and the following disclaimer in the documentation
22 // and/or other materials provided with the distribution.
23 //
24 // (3) Neither the name of the University of California, Lawrence Berkeley National
25 // Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
26 // used to endorse or promote products derived from this software without
27 // specific prior written permission.
28 //
29 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
30 // ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
31 // WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
32 // IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
33 // INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
34 // BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
35 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
36 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
37 // OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
38 // OF THE POSSIBILITY OF SUCH DAMAGE.
39 //
40 //=============================================================================
41 //
42 // This code is an extension of the algorithm presented in the paper:
43 // Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
44 // Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
45 // Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
46 // (LDAV), October 2016, Baltimore, Maryland.
47 //
48 // The PPP2 algorithm and software were jointly developed by
49 // Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
50 // Oliver Ruebel (LBNL)
51 //==============================================================================
52
53 #ifndef _vtk_m_filter_testing_TestingContourTreeUniformDistributedFilter_h_
54 #define _vtk_m_filter_testing_TestingContourTreeUniformDistributedFilter_h_
55
56 #include <vtkm/Types.h>
57 #include <vtkm/cont/DataSet.h>
58 #include <vtkm/cont/DataSetBuilderUniform.h>
59 #include <vtkm/cont/PartitionedDataSet.h>
60 #include <vtkm/cont/Serialization.h>
61 #include <vtkm/cont/testing/MakeTestDataSet.h>
62 #include <vtkm/cont/testing/Testing.h>
63 #include <vtkm/filter/ContourTreeUniformDistributed.h>
64 #include <vtkm/filter/MapFieldPermutation.h>
65 #include <vtkm/io/ErrorIO.h>
66 #include <vtkm/io/VTKDataSetReader.h>
67 #include <vtkm/worklet/contourtree_augmented/Types.h>
68 #include <vtkm/worklet/contourtree_distributed/TreeCompiler.h>
69
70 namespace vtkm
71 {
72 namespace filter
73 {
74 namespace testing
75 {
76 namespace contourtree_uniform_distributed
77 {
78 // numberOf Blocks must be a power of 2
ComputeNumberOfBlocksPerAxis(vtkm::Id3 globalSize,vtkm::Id numberOfBlocks)79 inline vtkm::Id3 ComputeNumberOfBlocksPerAxis(vtkm::Id3 globalSize, vtkm::Id numberOfBlocks)
80 {
81 // DEBUG: std::cout << "GlobalSize: " << globalSize << " numberOfBlocks:" << numberOfBlocks << " -> ";
82 // Inefficient way to compute log2 of numberOfBlocks, i.e., number of total splits
83 vtkm::Id numSplits = 0;
84 vtkm::Id currNumberOfBlock = numberOfBlocks;
85 bool isPowerOfTwo = true;
86 while (currNumberOfBlock > 1)
87 {
88 if (currNumberOfBlock % 2 != 0)
89 {
90 isPowerOfTwo = false;
91 break;
92 }
93 currNumberOfBlock /= 2;
94 ++numSplits;
95 }
96
97 if (isPowerOfTwo)
98 {
99 vtkm::Id3 splitsPerAxis{ 0, 0, 0 };
100 while (numSplits > 0)
101 {
102 // Find split axis as axis with largest extent
103 vtkm::IdComponent splitAxis = 0;
104 for (vtkm::IdComponent d = 1; d < 3; ++d)
105 if (globalSize[d] > globalSize[splitAxis])
106 splitAxis = d;
107 // Split in half along that axis
108 // DEBUG: std::cout << splitAxis << " " << globalSize << std::endl;
109 VTKM_ASSERT(globalSize[splitAxis] > 1);
110 ++splitsPerAxis[splitAxis];
111 globalSize[splitAxis] /= 2;
112 --numSplits;
113 }
114 // DEBUG: std::cout << "splitsPerAxis: " << splitsPerAxis;
115 vtkm::Id3 blocksPerAxis;
116 for (vtkm::IdComponent d = 0; d < 3; ++d)
117 blocksPerAxis[d] = vtkm::Id{ 1 } << splitsPerAxis[d];
118 // DEBUG: std::cout << " blocksPerAxis: " << blocksPerAxis << std::endl;
119 return blocksPerAxis;
120 }
121 else
122 {
123 std::cout << "numberOfBlocks is not a power of two. Splitting along longest axis." << std::endl;
124 vtkm::IdComponent splitAxis = 0;
125 for (vtkm::IdComponent d = 1; d < 3; ++d)
126 if (globalSize[d] > globalSize[splitAxis])
127 splitAxis = d;
128 vtkm::Id3 blocksPerAxis{ 1, 1, 1 };
129 blocksPerAxis[splitAxis] = numberOfBlocks;
130 // DEBUG: std::cout << " blocksPerAxis: " << blocksPerAxis << std::endl;
131 return blocksPerAxis;
132 }
133 }
134
ComputeBlockExtents(vtkm::Id3 globalSize,vtkm::Id3 blocksPerAxis,vtkm::Id blockNo)135 inline std::tuple<vtkm::Id3, vtkm::Id3, vtkm::Id3> ComputeBlockExtents(vtkm::Id3 globalSize,
136 vtkm::Id3 blocksPerAxis,
137 vtkm::Id blockNo)
138 {
139 // DEBUG: std::cout << "ComputeBlockExtents("<<globalSize <<", " << blocksPerAxis << ", " << blockNo << ")" << std::endl;
140 // DEBUG: std::cout << "Block " << blockNo;
141
142 vtkm::Id3 blockIndex, blockOrigin, blockSize;
143 for (vtkm::IdComponent d = 0; d < 3; ++d)
144 {
145 blockIndex[d] = blockNo % blocksPerAxis[d];
146 blockNo /= blocksPerAxis[d];
147
148 float dx = float(globalSize[d] - 1) / float(blocksPerAxis[d]);
149 blockOrigin[d] = vtkm::Id(blockIndex[d] * dx);
150 vtkm::Id maxIdx =
151 blockIndex[d] < blocksPerAxis[d] - 1 ? vtkm::Id((blockIndex[d] + 1) * dx) : globalSize[d] - 1;
152 blockSize[d] = maxIdx - blockOrigin[d] + 1;
153 // DEBUG: std::cout << " " << blockIndex[d] << dx << " " << blockOrigin[d] << " " << maxIdx << " " << blockSize[d] << "; ";
154 }
155 // DEBUG: std::cout << " -> " << blockIndex << " " << blockOrigin << " " << blockSize << std::endl;
156 return std::make_tuple(blockIndex, blockOrigin, blockSize);
157 }
158
CreateSubDataSet(const vtkm::cont::DataSet & ds,vtkm::Id3 blockOrigin,vtkm::Id3 blockSize,const std::string & fieldName)159 inline vtkm::cont::DataSet CreateSubDataSet(const vtkm::cont::DataSet& ds,
160 vtkm::Id3 blockOrigin,
161 vtkm::Id3 blockSize,
162 const std::string& fieldName)
163 {
164 vtkm::Id3 globalSize;
165 ds.GetCellSet().CastAndCall(vtkm::worklet::contourtree_augmented::GetPointDimensions(),
166 globalSize);
167 const vtkm::Id nOutValues = blockSize[0] * blockSize[1] * blockSize[2];
168
169 const auto inDataArrayHandle = ds.GetPointField(fieldName).GetData();
170
171 vtkm::cont::ArrayHandle<vtkm::Id> copyIdsArray;
172 copyIdsArray.Allocate(nOutValues);
173 auto copyIdsPortal = copyIdsArray.WritePortal();
174
175 vtkm::Id3 outArrIdx;
176 for (outArrIdx[2] = 0; outArrIdx[2] < blockSize[2]; ++outArrIdx[2])
177 for (outArrIdx[1] = 0; outArrIdx[1] < blockSize[1]; ++outArrIdx[1])
178 for (outArrIdx[0] = 0; outArrIdx[0] < blockSize[0]; ++outArrIdx[0])
179 {
180 vtkm::Id3 inArrIdx = outArrIdx + blockOrigin;
181 vtkm::Id inIdx = (inArrIdx[2] * globalSize[1] + inArrIdx[1]) * globalSize[0] + inArrIdx[0];
182 vtkm::Id outIdx =
183 (outArrIdx[2] * blockSize[1] + outArrIdx[1]) * blockSize[0] + outArrIdx[0];
184 VTKM_ASSERT(inIdx >= 0 && inIdx < inDataArrayHandle.GetNumberOfValues());
185 VTKM_ASSERT(outIdx >= 0 && outIdx < nOutValues);
186 copyIdsPortal.Set(outIdx, inIdx);
187 }
188 // DEBUG: std::cout << copyIdsPortal.GetNumberOfValues() << std::endl;
189
190 vtkm::cont::Field permutedField;
191 bool success =
192 vtkm::filter::MapFieldPermutation(ds.GetPointField(fieldName), copyIdsArray, permutedField);
193 if (!success)
194 throw vtkm::cont::ErrorBadType("Field copy failed (probably due to invalid type)");
195
196 vtkm::cont::DataSetBuilderUniform dsb;
197 if (globalSize[2] <= 1) // 2D Data Set
198 {
199 vtkm::Id2 dimensions{ blockSize[0], blockSize[1] };
200 vtkm::cont::DataSet dataSet = dsb.Create(dimensions);
201 dataSet.AddField(permutedField);
202 return dataSet;
203 }
204 else
205 {
206 vtkm::cont::DataSet dataSet = dsb.Create(blockSize);
207 dataSet.AddField(permutedField);
208 return dataSet;
209 }
210 }
211
ReadGroundTruthContourTree(std::string filename)212 inline std::vector<vtkm::worklet::contourtree_distributed::Edge> ReadGroundTruthContourTree(
213 std::string filename)
214 {
215 std::ifstream ct_file(filename);
216 if (!ct_file.is_open())
217 {
218 throw vtkm::io::ErrorIO("Unable to open data file: " + filename);
219 }
220 vtkm::Id val1, val2;
221 std::vector<vtkm::worklet::contourtree_distributed::Edge> result;
222 while (ct_file >> val1 >> val2)
223 {
224 result.push_back(vtkm::worklet::contourtree_distributed::Edge(val1, val2));
225 }
226 std::sort(result.begin(), result.end());
227 return result;
228 }
229
230 inline vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed(
231 const vtkm::cont::DataSet& ds,
232 std::string fieldName,
233 bool useMarchingCubes,
234 int numberOfBlocks,
235 int rank = 0,
236 int numberOfRanks = 1)
237 {
238 // Get dimensions of data set
239 vtkm::Id3 globalSize;
240 ds.GetCellSet().CastAndCall(vtkm::worklet::contourtree_augmented::GetPointDimensions(),
241 globalSize);
242
243 // Determine split
244 vtkm::Id3 blocksPerAxis = ComputeNumberOfBlocksPerAxis(globalSize, numberOfBlocks);
245 vtkm::Id blocksPerRank = numberOfBlocks / numberOfRanks;
246 vtkm::Id numRanksWithExtraBlock = numberOfBlocks % numberOfRanks;
247 vtkm::Id blocksOnThisRank, startBlockNo;
248 if (rank < numRanksWithExtraBlock)
249 {
250 blocksOnThisRank = blocksPerRank + 1;
251 startBlockNo = (blocksPerRank + 1) * rank;
252 }
253 else
254 {
255 blocksOnThisRank = blocksPerRank;
256 startBlockNo = numRanksWithExtraBlock * (blocksPerRank + 1) +
257 (rank - numRanksWithExtraBlock) * blocksPerRank;
258 }
259
260 // Created partitioned (split) data set
261 vtkm::cont::PartitionedDataSet pds;
262 vtkm::cont::ArrayHandle<vtkm::Id3> localBlockIndices;
263 vtkm::cont::ArrayHandle<vtkm::Id3> localBlockOrigins;
264 vtkm::cont::ArrayHandle<vtkm::Id3> localBlockSizes;
265 localBlockIndices.Allocate(blocksOnThisRank);
266 localBlockOrigins.Allocate(blocksOnThisRank);
267 localBlockSizes.Allocate(blocksOnThisRank);
268 auto localBlockIndicesPortal = localBlockIndices.WritePortal();
269 auto localBlockOriginsPortal = localBlockOrigins.WritePortal();
270 auto localBlockSizesPortal = localBlockSizes.WritePortal();
271
272 for (vtkm::Id blockNo = 0; blockNo < blocksOnThisRank; ++blockNo)
273 {
274 vtkm::Id3 blockOrigin, blockSize, blockIndex;
275 std::tie(blockIndex, blockOrigin, blockSize) =
276 ComputeBlockExtents(globalSize, blocksPerAxis, startBlockNo + blockNo);
277 pds.AppendPartition(CreateSubDataSet(ds, blockOrigin, blockSize, fieldName));
278 localBlockOriginsPortal.Set(blockNo, blockOrigin);
279 localBlockSizesPortal.Set(blockNo, blockSize);
280 localBlockIndicesPortal.Set(blockNo, blockIndex);
281 }
282
283 // Run the contour tree analysis
284 vtkm::filter::ContourTreeUniformDistributed filter(blocksPerAxis,
285 globalSize,
286 localBlockIndices,
287 localBlockOrigins,
288 localBlockSizes,
289 // Freudenthal: Only use boundary extrema
290 // MC: use all points on boundary
291 // TODO/FIXME: Figure out why MC does not
292 // work when only using boundary extrema
293 !useMarchingCubes,
294 useMarchingCubes,
295 false,
296 vtkm::cont::LogLevel::UserVerboseLast,
297 vtkm::cont::LogLevel::UserVerboseLast);
298 filter.SetActiveField(fieldName);
299 auto result = filter.Execute(pds);
300
301 if (numberOfRanks == 1)
302 {
303 // Serial or only one parallel rank -> Result is already
304 // everything we need
305 return result;
306 }
307 else
308 {
309 // Mutiple ranks -> Some assembly required. Collect data
310 // on rank 0, all other ranks return empty data sets
311 using FieldTypeList =
312 vtkm::ListAppend<vtkm::filter::ContourTreeUniformDistributed::SupportedTypes,
313 vtkm::List<vtkm::Id>>;
314 using DataSetWrapper =
315 vtkm::cont::SerializableDataSet<FieldTypeList, vtkm::cont::CellSetListStructured>;
316
317 // Communicate results to rank 0
318 auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
319 vtkmdiy::Master master(comm, 1);
320 struct EmptyBlock
321 {
322 }; // Dummy block structure, since we need block data for DIY
323 master.add(comm.rank(), new EmptyBlock, new vtkmdiy::Link);
324 // .. Send data to rank 0
325 master.foreach ([result, filter](void*, const vtkmdiy::Master::ProxyWithLink& p) {
326 vtkmdiy::BlockID root{ 0, 0 }; // Rank 0
327 p.enqueue(root, result.GetNumberOfPartitions());
328 for (const vtkm::cont::DataSet& curr_ds : result)
329 {
330 auto curr_sds = DataSetWrapper(curr_ds);
331 p.enqueue(root, curr_sds);
332 }
333 });
334 // Exchange data, i.e., send to rank 0 (pass "true" to exchange data between
335 // *all* blocks, not just neighbors)
336 master.exchange(true);
337
338 if (comm.rank() == 0)
339 {
340 // Receive data on rank zero and return combined results
341 vtkm::cont::PartitionedDataSet combined_result;
342 master.foreach ([&combined_result, filter, numberOfRanks](
343 void*, const vtkmdiy::Master::ProxyWithLink& p) {
344 for (int receiveFromRank = 0; receiveFromRank < numberOfRanks; ++receiveFromRank)
345 {
346 vtkm::Id numberOfDataSetsToReceive;
347 p.dequeue({ receiveFromRank, receiveFromRank }, numberOfDataSetsToReceive);
348 for (vtkm::Id currReceiveDataSetNo = 0; currReceiveDataSetNo < numberOfDataSetsToReceive;
349 ++currReceiveDataSetNo)
350 {
351 auto sds = vtkm::filter::MakeSerializableDataSet(filter);
352 p.dequeue({ receiveFromRank, receiveFromRank }, sds);
353 combined_result.AppendPartition(sds.DataSet);
354 }
355 }
356 });
357 return combined_result; // Return combined result on rank 0
358 }
359 else
360 {
361 // Return an empty data set on all other ranks
362 return vtkm::cont::PartitionedDataSet{};
363 }
364 }
365 }
366
367 inline void TestContourTreeUniformDistributed8x9(int nBlocks, int rank = 0, int size = 1)
368 {
369 if (rank == 0)
370 {
371 std::cout << "Testing ContourTreeUniformDistributed on 2D 8x9 data set divided into " << nBlocks
372 << " blocks." << std::endl;
373 }
374 vtkm::cont::DataSet in_ds = vtkm::cont::testing::MakeTestDataSet().Make2DUniformDataSet3();
375 vtkm::cont::PartitionedDataSet result =
376 RunContourTreeDUniformDistributed(in_ds, "pointvar", false, nBlocks, rank, size);
377
378 if (vtkm::cont::EnvironmentTracker::GetCommunicator().rank() == 0)
379 {
380 vtkm::worklet::contourtree_distributed::TreeCompiler treeCompiler;
381 for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
382 {
383 treeCompiler.AddHierarchicalTree(result.GetPartition(ds_no));
384 }
385 treeCompiler.ComputeSuperarcs();
386
387 // Print the contour tree we computed
388 std::cout << "Computed Contour Tree" << std::endl;
389 treeCompiler.PrintSuperarcs();
390
391 // Print the expected contour tree
392 std::cout << "Expected Contour Tree" << std::endl;
393 std::cout << " 10 20" << std::endl;
394 std::cout << " 20 34" << std::endl;
395 std::cout << " 20 38" << std::endl;
396 std::cout << " 20 61" << std::endl;
397 std::cout << " 23 34" << std::endl;
398 std::cout << " 24 34" << std::endl;
399 std::cout << " 50 61" << std::endl;
400 std::cout << " 61 71" << std::endl;
401
402 using Edge = vtkm::worklet::contourtree_distributed::Edge;
403 VTKM_TEST_ASSERT(test_equal(treeCompiler.superarcs.size(), 8),
404 "Wrong result for ContourTreeUniformDistributed filter");
405 VTKM_TEST_ASSERT(treeCompiler.superarcs[0] == Edge{ 10, 20 },
406 "Wrong result for ContourTreeUniformDistributed filter");
407 VTKM_TEST_ASSERT(treeCompiler.superarcs[1] == Edge{ 20, 34 },
408 "Wrong result for ContourTreeUniformDistributed filter");
409 VTKM_TEST_ASSERT(treeCompiler.superarcs[2] == Edge{ 20, 38 },
410 "Wrong result for ContourTreeUniformDistributed filter");
411 VTKM_TEST_ASSERT(treeCompiler.superarcs[3] == Edge{ 20, 61 },
412 "Wrong result for ContourTreeUniformDistributed filter");
413 VTKM_TEST_ASSERT(treeCompiler.superarcs[4] == Edge{ 23, 34 },
414 "Wrong result for ContourTreeUniformDistributed filter");
415 VTKM_TEST_ASSERT(treeCompiler.superarcs[5] == Edge{ 24, 34 },
416 "Wrong result for ContourTreeUniformDistributed filter");
417 VTKM_TEST_ASSERT(treeCompiler.superarcs[6] == Edge{ 50, 61 },
418 "Wrong result for ContourTreeUniformDistributed filter");
419 VTKM_TEST_ASSERT(treeCompiler.superarcs[7] == Edge{ 61, 71 },
420 "Wrong result for ContourTreeUniformDistributed filter");
421 }
422 }
423
424 inline void TestContourTreeUniformDistributed5x6x7(int nBlocks,
425 bool marchingCubes,
426 int rank = 0,
427 int size = 1)
428 {
429 if (rank == 0)
430 {
431 std::cout << "Testing ContourTreeUniformDistributed with "
432 << (marchingCubes ? "marching cubes" : "Freudenthal")
433 << " mesh connectivity on 3D 5x6x7 data set divided into " << nBlocks << " blocks."
434 << std::endl;
435 }
436
437 vtkm::cont::DataSet in_ds = vtkm::cont::testing::MakeTestDataSet().Make3DUniformDataSet4();
438 vtkm::cont::PartitionedDataSet result =
439 RunContourTreeDUniformDistributed(in_ds, "pointvar", marchingCubes, nBlocks, rank, size);
440
441 if (rank == 0)
442 {
443 vtkm::worklet::contourtree_distributed::TreeCompiler treeCompiler;
444 for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
445 {
446 treeCompiler.AddHierarchicalTree(result.GetPartition(ds_no));
447 }
448 treeCompiler.ComputeSuperarcs();
449
450 // Print the contour tree we computed
451 std::cout << "Computed Contour Tree" << std::endl;
452 treeCompiler.PrintSuperarcs();
453
454 // Print the expected contour tree
455 using Edge = vtkm::worklet::contourtree_distributed::Edge;
456 std::cout << "Expected Contour Tree" << std::endl;
457 if (!marchingCubes)
458 {
459 std::cout << " 0 112" << std::endl;
460 std::cout << " 71 72" << std::endl;
461 std::cout << " 72 78" << std::endl;
462 std::cout << " 72 101" << std::endl;
463 std::cout << " 101 112" << std::endl;
464 std::cout << " 101 132" << std::endl;
465 std::cout << " 107 112" << std::endl;
466 std::cout << " 131 132" << std::endl;
467 std::cout << " 132 138" << std::endl;
468
469 VTKM_TEST_ASSERT(test_equal(treeCompiler.superarcs.size(), 9),
470 "Wrong result for ContourTreeUniformDistributed filter");
471 VTKM_TEST_ASSERT(treeCompiler.superarcs[0] == Edge{ 0, 112 },
472 "Wrong result for ContourTreeUniformDistributed filter");
473 VTKM_TEST_ASSERT(treeCompiler.superarcs[1] == Edge{ 71, 72 },
474 "Wrong result for ContourTreeUniformDistributed filter");
475 VTKM_TEST_ASSERT(treeCompiler.superarcs[2] == Edge{ 72, 78 },
476 "Wrong result for ContourTreeUniformDistributed filter");
477 VTKM_TEST_ASSERT(treeCompiler.superarcs[3] == Edge{ 72, 101 },
478 "Wrong result for ContourTreeUniformDistributed filter");
479 VTKM_TEST_ASSERT(treeCompiler.superarcs[4] == Edge{ 101, 112 },
480 "Wrong result for ContourTreeUniformDistributed filter");
481 VTKM_TEST_ASSERT(treeCompiler.superarcs[5] == Edge{ 101, 132 },
482 "Wrong result for ContourTreeUniformDistributed filter");
483 VTKM_TEST_ASSERT(treeCompiler.superarcs[6] == Edge{ 107, 112 },
484 "Wrong result for ContourTreeUniformDistributed filter");
485 VTKM_TEST_ASSERT(treeCompiler.superarcs[7] == Edge{ 131, 132 },
486 "Wrong result for ContourTreeUniformDistributed filter");
487 VTKM_TEST_ASSERT(treeCompiler.superarcs[8] == Edge{ 132, 138 },
488 "Wrong result for ContourTreeUniformDistributed filter");
489 }
490 else
491 {
492 std::cout << " 0 203" << std::endl;
493 std::cout << " 71 72" << std::endl;
494 std::cout << " 72 78" << std::endl;
495 std::cout << " 72 101" << std::endl;
496 std::cout << " 101 112" << std::endl;
497 std::cout << " 101 132" << std::endl;
498 std::cout << " 107 112" << std::endl;
499 std::cout << " 112 203" << std::endl;
500 std::cout << " 131 132" << std::endl;
501 std::cout << " 132 138" << std::endl;
502 std::cout << " 203 209" << std::endl;
503
504 VTKM_TEST_ASSERT(test_equal(treeCompiler.superarcs.size(), 11),
505 "Wrong result for ContourTreeUniformDistributed filter");
506 VTKM_TEST_ASSERT(treeCompiler.superarcs[0] == Edge{ 0, 203 },
507 "Wrong result for ContourTreeUniformDistributed filter");
508 VTKM_TEST_ASSERT(treeCompiler.superarcs[1] == Edge{ 71, 72 },
509 "Wrong result for ContourTreeUniformDistributed filter");
510 VTKM_TEST_ASSERT(treeCompiler.superarcs[2] == Edge{ 72, 78 },
511 "Wrong result for ContourTreeUniformDistributed filter");
512 VTKM_TEST_ASSERT(treeCompiler.superarcs[3] == Edge{ 72, 101 },
513 "Wrong result for ContourTreeUniformDistributed filter");
514 VTKM_TEST_ASSERT(treeCompiler.superarcs[4] == Edge{ 101, 112 },
515 "Wrong result for ContourTreeUniformDistributed filter");
516 VTKM_TEST_ASSERT(treeCompiler.superarcs[5] == Edge{ 101, 132 },
517 "Wrong result for ContourTreeUniformDistributed filter");
518 VTKM_TEST_ASSERT(treeCompiler.superarcs[6] == Edge{ 107, 112 },
519 "Wrong result for ContourTreeUniformDistributed filter");
520 VTKM_TEST_ASSERT(treeCompiler.superarcs[7] == Edge{ 112, 203 },
521 "Wrong result for ContourTreeUniformDistributed filter");
522 VTKM_TEST_ASSERT(treeCompiler.superarcs[8] == Edge{ 131, 132 },
523 "Wrong result for ContourTreeUniformDistributed filter");
524 VTKM_TEST_ASSERT(treeCompiler.superarcs[9] == Edge{ 132, 138 },
525 "Wrong result for ContourTreeUniformDistributed filter");
526 VTKM_TEST_ASSERT(treeCompiler.superarcs[10] == Edge{ 203, 209 },
527 "Wrong result for ContourTreeUniformDistributed filter");
528 }
529 }
530 }
531
532 inline void TestContourTreeFile(std::string ds_filename,
533 std::string fieldName,
534 std::string gtct_filename,
535 int nBlocks,
536 bool marchingCubes = false,
537 int rank = 0,
538 int size = 1)
539 {
540 if (rank == 0)
541 {
542 std::cout << "Testing ContourTreeUniformDistributed with "
543 << (marchingCubes ? "marching cubes" : "Freudenthal") << " mesh connectivity on \""
544 << ds_filename << "\" divided into " << nBlocks << " blocks." << std::endl;
545 }
546
547 vtkm::io::VTKDataSetReader reader(ds_filename);
548 vtkm::cont::DataSet ds;
549 try
550 {
551 ds = reader.ReadDataSet();
552 }
catch(vtkm::io::ErrorIO & e)553 catch (vtkm::io::ErrorIO& e)
554 {
555 std::string message("Error reading: ");
556 message += ds_filename;
557 message += ", ";
558 message += e.GetMessage();
559
560 VTKM_TEST_FAIL(message.c_str());
561 }
562
563 vtkm::cont::PartitionedDataSet result =
564 RunContourTreeDUniformDistributed(ds, fieldName, marchingCubes, nBlocks, rank, size);
565
566 if (rank == 0)
567 {
568 vtkm::worklet::contourtree_distributed::TreeCompiler treeCompiler;
569 for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
570 {
571 treeCompiler.AddHierarchicalTree(result.GetPartition(ds_no));
572 }
573 treeCompiler.ComputeSuperarcs();
574
575 std::vector<vtkm::worklet::contourtree_distributed::Edge> groundTruthSuperarcs =
576 ReadGroundTruthContourTree(gtct_filename);
577 if (groundTruthSuperarcs.size() < 50)
578 {
579 std::cout << "Computed Contour Tree" << std::endl;
580 treeCompiler.PrintSuperarcs();
581
582 // Print the expected contour tree
583 std::cout << "Expected Contour Tree" << std::endl;
584 vtkm::worklet::contourtree_distributed::TreeCompiler::PrintSuperarcArray(
585 groundTruthSuperarcs);
586 }
587 else
588 {
589 std::cout << "Not printing computed and expected contour tree due to size." << std::endl;
590 }
591
592 VTKM_TEST_ASSERT(treeCompiler.superarcs == groundTruthSuperarcs,
593 "Test failed for data set " + ds_filename);
594 }
595 }
596
597 }
598 }
599 }
600 } // vtkm::filter::testing::contourtree_uniform_distributed
601
602 #endif
603