ContourTreeMesh.h 39.7 KB
Newer Older
Sujin Philip's avatar
Sujin Philip committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
//============================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
//  This software is distributed WITHOUT ANY WARRANTY; without even
//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
//  PURPOSE.  See the above copyright notice for more information.
//
//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
//  Copyright 2014 UT-Battelle, LLC.
//  Copyright 2014 Los Alamos National Security.
//
//  Under the terms of Contract DE-NA0003525 with NTESS,
//  the U.S. Government retains certain rights in this software.
//
//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
//  Laboratory (LANL), the U.S. Government retains certain rights in
//  this software.
//============================================================================
// Copyright (c) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy).  All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// (1) Redistributions of source code must retain the above copyright notice, this
//     list of conditions and the following disclaimer.
//
// (2) Redistributions in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
// (3) Neither the name of the University of California, Lawrence Berkeley National
//     Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
//     used to endorse or promote products derived from this software without
//     specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
//=============================================================================
//
//  This code is an extension of the algorithm presented in the paper:
//  Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
//  Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
//  Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
//  (LDAV), October 2016, Baltimore, Maryland.
//
//  The PPP2 algorithm and software were jointly developed by
//  Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
//  Oliver Ruebel (LBNL)
//==============================================================================

63
64
#ifndef vtk_m_worklet_contourtree_augmented_mesh_dem_contour_tree_mesh_h
#define vtk_m_worklet_contourtree_augmented_mesh_dem_contour_tree_mesh_h
Sujin Philip's avatar
Sujin Philip committed
65
66
67
68
69

#include <cstdlib>
#include <fstream>
#include <iostream>
#include <vtkm/Types.h>
70
#include <vtkm/cont/Algorithm.h>
Sujin Philip's avatar
Sujin Philip committed
71
72
73
74
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/cont/ArrayPortalToIterators.h>
#include <vtkm/cont/ArrayRangeCompute.h>
#include <vtkm/cont/EnvironmentTracker.h>
75
#include <vtkm/io/ErrorIO.h>
Sujin Philip's avatar
Sujin Philip committed
76
77
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/contourtree_augmented/ArrayTransforms.h>
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#include <vtkm/worklet/contourtree_augmented/data_set_mesh/IdRelabeler.h> // This is needed only as an unused default argument.
#include <vtkm/worklet/contourtree_augmented/meshtypes/MeshStructureContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ArcComparator.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedOtherStartIndexNNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedSimulatedSimplicityIndexComparator.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CombinedVectorDifferentFromNext.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/CompressNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ComputeMaxNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/FindStartIndexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/InitToCombinedSortOrderArraysWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/MergeCombinedOtherStartIndexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/ReplaceArcNumWithToVertexWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/SubtractAssignWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/contourtreemesh/UpdateCombinedNeighboursWorklet.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/ComputeMeshBoundaryContourTreeMesh.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundaryContourTreeMesh.h>
94

Sujin Philip's avatar
Sujin Philip committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h> // TODO remove should not be needed

#include <vtkm/cont/ExecutionObjectBase.h>

#include <vtkm/cont/Invoker.h>

#include <algorithm>

VTKM_THIRDPARTY_PRE_INCLUDE
#include <vtkm/thirdparty/diy/Configure.h>
#include <vtkm/thirdparty/diy/diy.h>
VTKM_THIRDPARTY_POST_INCLUDE
// clang-format on

namespace contourtree_mesh_inc_ns =
  vtkm::worklet::contourtree_augmented::mesh_dem_contourtree_mesh_inc;

namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{

template <typename FieldType>
class ContourTreeMesh : public vtkm::cont::ExecutionObjectBase
{ // class ContourTreeMesh
public:
  //Mesh dependent helper functions
124
  void SetPrepareForExecutionBehavior(bool getMax);
Sujin Philip's avatar
Sujin Philip committed
125

126
127
  contourtree_mesh_inc_ns::MeshStructureContourTreeMesh PrepareForExecution(
    vtkm::cont::DeviceAdapterId,
128
    vtkm::cont::Token& token) const;
Sujin Philip's avatar
Sujin Philip committed
129
130
131

  ContourTreeMesh() {}

132
  // Constructor
Sujin Philip's avatar
Sujin Philip committed
133
134
135
136
137
  ContourTreeMesh(const IdArrayType& arcs,
                  const IdArrayType& inSortOrder,
                  const vtkm::cont::ArrayHandle<FieldType>& values,
                  const IdArrayType& inGlobalMeshIndex);

138
139
140
141
142
143
144
  // Constructor
  ContourTreeMesh(const IdArrayType& nodes,
                  const IdArrayType& arcs,
                  const IdArrayType& inSortOrder,
                  const vtkm::cont::ArrayHandle<FieldType>& values,
                  const IdArrayType& inGlobalMeshIndex);

145
  //  Construct a ContourTreeMesh from nodes/arcs and another ContourTreeMesh (instead of a DataSetMesh)
Sujin Philip's avatar
Sujin Philip committed
146
147
148
149
150
151
152
153
154
155
156
157
158
159
  //     nodes/arcs: From the contour tree
  //     ContourTreeMesh: the contour tree mesh used to compute the contour tree described by nodes/arcs
  ContourTreeMesh(const IdArrayType& nodes,
                  const IdArrayType& arcs,
                  const ContourTreeMesh<FieldType>& mesh);

  // Initalize contour tree mesh from mesh and arcs. For fully augmented contour tree with all
  // mesh vertices as nodes. Same as using { 0, 1, ..., nodes.size()-1 } as nodes for the
  // ContourTreeMeshh(nodes, arcsm mesh) constructor above
  ContourTreeMesh(const IdArrayType& arcs, const ContourTreeMesh<FieldType>& mesh);

  // Load contour tree mesh from file
  ContourTreeMesh(const char* filename)
  {
160
161
    Load(filename);
    this->NumVertices = this->SortedValues.GetNumberOfValues();
Sujin Philip's avatar
Sujin Philip committed
162
163
  }

164
  vtkm::Id GetNumberOfVertices() const { return this->NumVertices; }
165

Sujin Philip's avatar
Sujin Philip committed
166
  // Combine two ContourTreeMeshes
167
  void MergeWith(ContourTreeMesh<FieldType>& other);
Sujin Philip's avatar
Sujin Philip committed
168

169
170
171
  // Save/Load the mesh helpers
  void Save(const char* filename) const;
  void Load(const char* filename);
Sujin Philip's avatar
Sujin Philip committed
172
173
174
175
176
177
178
179
180
181
182
183
184

  // Empty placeholder function to ensure compliance of this class with the interface
  // the other mesh classes. This is a no-op here since this class is initalized
  // from a known contour tree so sort is already done
  template <typename T, typename StorageType>
  void SortData(const vtkm::cont::ArrayHandle<T, StorageType>& values) const
  {
    (void)values; // Do nothink but avoid unsused param warning
  }

  // Public fields
  static const int MAX_OUTDEGREE = 20;

185
  vtkm::Id NumVertices;
186
187
  vtkm::cont::ArrayHandleIndex SortOrder;
  vtkm::cont::ArrayHandleIndex SortIndices;
188
189
  vtkm::cont::ArrayHandle<FieldType> SortedValues;
  IdArrayType GlobalMeshIndex;
Sujin Philip's avatar
Sujin Philip committed
190
191
192
193
194
195
  // neighbours stores for each vertex the indices of its neighbours. For each vertex
  // the indices are sorted by value, i.e, the first neighbour has the lowest and
  // the last neighbour the highest value for the vertex. In the array we just
  // concatinate the list of neighbours from all vertices, i.e., we first
  // have the list of neighbours of the first vertex, then the second vertex and so on, i.e.:
  // [ n_1_1, n_1_2, n_2_1, n_2_2, n_2_3, etc.]
196
  IdArrayType Neighbours;
Sujin Philip's avatar
Sujin Philip committed
197
198
  // firstNeighour gives us for each vertex an index into the neighours array indicating
  // the index where the list of neighbours for the vertex begins
199
  IdArrayType FirstNeighbour;
Sujin Philip's avatar
Sujin Philip committed
200
  // the maximum number of neighbours of a vertex
201
  vtkm::Id MaxNeighbours;
Sujin Philip's avatar
Sujin Philip committed
202

203
204
205
  // Print Contents
  void PrintContent(std::ostream& outStream = std::cout) const;

Sujin Philip's avatar
Sujin Philip committed
206
  // Debug print routine
207
  void DebugPrint(const char* message, const char* fileName, long lineNum) const;
Sujin Philip's avatar
Sujin Philip committed
208

209
  // Get boundary execution object
210
  MeshBoundaryContourTreeMeshExec GetMeshBoundaryExecutionObject(vtkm::Id3 globalSize,
211
212
213
                                                                 vtkm::Id3 minIdx,
                                                                 vtkm::Id3 maxIdx) const;

214
215
216
  void GetBoundaryVertices(IdArrayType& boundaryVertexArray,                    // output
                           IdArrayType& boundarySortIndexArray,                 // output
                           MeshBoundaryContourTreeMeshExec* meshBoundaryExecObj //input
217
  ) const;
218

219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
  /// copies the global IDs for a set of sort IDs
  /// notice that the sort ID is the same as the mesh ID for the ContourTreeMesh class.
  /// To reduce memory usage we here use a fancy array handle rather than copy data
  /// as is needed for the DataSetMesh types.
  /// We here return a fancy array handle to convert values on-the-fly without requiring additional memory
  /// @param[in] sortIds Array with sort Ids to be converted from local to global Ids
  /// @param[in] localToGlobalIdRelabeler This parameter is here only for
  ///            consistency with the DataSetMesh types but is not
  ///            used here and as such can simply be set to nullptr
  inline vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType> GetGlobalIdsFromSortIndices(
    const IdArrayType& sortIds,
    const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler =
      nullptr) const
  {                                 // GetGlobalIDsFromSortIndices()
    (void)localToGlobalIdRelabeler; // avoid compiler warning
    return vtkm::cont::make_ArrayHandlePermutation(sortIds, this->GlobalMeshIndex);
  } // GetGlobalIDsFromSortIndices()

  /// copies the global IDs for a set of mesh IDs
  /// notice that the sort ID is the same as the mesh ID for the ContourTreeMesh class.
  /// To reduce memory usage we here use a fancy array handle rather than copy data
  /// as is needed for the DataSetMesh types.
241
242
  /// MeshIdArrayType must be an array if Ids. Usually this is a vtkm::worklet::contourtree_augmented::IdArrayType
  /// but in some cases it may also be a fancy array to avoid memory allocation
243
244
245
246
247
  /// We here return a fancy array handle to convert values on-the-fly without requiring additional memory
  /// @param[in] meshIds Array with mesh Ids to be converted from local to global Ids
  /// @param[in] localToGlobalIdRelabeler This parameter is here only for
  ///            consistency with the DataSetMesh types but is not
  ///            used here and as such can simply be set to nullptr
248
249
250
251
252
  template <typename MeshIdArrayType>
  inline vtkm::cont::ArrayHandlePermutation<MeshIdArrayType, IdArrayType>
  GetGlobalIdsFromMeshIndices(const MeshIdArrayType& meshIds,
                              const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler*
                                localToGlobalIdRelabeler = nullptr) const
253
254
255
256
257
  {                                 // GetGlobalIDsFromMeshIndices()
    (void)localToGlobalIdRelabeler; // avoid compiler warning
    return vtkm::cont::make_ArrayHandlePermutation(meshIds, this->GlobalMeshIndex);
  } // GetGlobalIDsFromMeshIndices()

Sujin Philip's avatar
Sujin Philip committed
258
259
260
261
262
263
264
265
266
267
268
private:
  vtkm::cont::Invoker Invoke;

  bool mGetMax; // Define the behavior for the PrepareForExecution function

  // Private init and helper functions
  void InitialiseNeighboursFromArcs(const IdArrayType& arcs);
  void ComputeNNeighboursVector(IdArrayType& nNeighbours) const;
  void ComputeMaxNeighbours();

  // Private helper functions for saving data vectors
269
  // Internal helper function to save 1D index array to file
Sujin Philip's avatar
Sujin Philip committed
270
  template <typename ValueType>
271
  void SaveVector(std::ostream& os, const vtkm::cont::ArrayHandle<ValueType>& vec) const;
Sujin Philip's avatar
Sujin Philip committed
272

273
274
  // Internal helper function to Load 1D index array from file
  template <typename ValueType>
275
  void LoadVector(std::istream& is, vtkm::cont::ArrayHandle<ValueType>& vec);
Sujin Philip's avatar
Sujin Philip committed
276
277
278

}; // ContourTreeMesh

279
280
// print content
template <typename FieldType>
281
inline void ContourTreeMesh<FieldType>::PrintContent(std::ostream& outStream /*= std::cout*/) const
282
283
284
285
286
287
288
289
290
{ // PrintContent()
  PrintHeader(this->NumVertices, outStream);
  //PrintIndices("SortOrder", SortOrder, outStream);
  PrintValues("SortedValues", SortedValues, -1, outStream);
  PrintIndices("GlobalMeshIndex", GlobalMeshIndex, -1, outStream);
  PrintIndices("Neighbours", Neighbours, -1, outStream);
  PrintIndices("FirstNeighbour", FirstNeighbour, -1, outStream);
  outStream << "MaxNeighbours=" << MaxNeighbours << std::endl;
  outStream << "mGetMax=" << mGetMax << std::endl;
291
} // PrintContent()
Sujin Philip's avatar
Sujin Philip committed
292
293
294

// debug routine
template <typename FieldType>
295
296
297
inline void ContourTreeMesh<FieldType>::DebugPrint(const char* message,
                                                   const char* fileName,
                                                   long lineNum) const
Sujin Philip's avatar
Sujin Philip committed
298
299
300
301
302
303
304
305
306
{ // DebugPrint()
  std::cout << "---------------------------" << std::endl;
  std::cout << std::setw(30) << std::left << fileName << ":" << std::right << std::setw(4)
            << lineNum << std::endl;
  std::cout << std::left << std::string(message) << std::endl;
  std::cout << "Contour Tree Mesh Contains:     " << std::endl;
  std::cout << "---------------------------" << std::endl;
  std::cout << std::endl;

307
  PrintContent(std::cout);
Sujin Philip's avatar
Sujin Philip committed
308
309
310
311
312
313
314
315
} // DebugPrint()

// create the contour tree mesh from contour tree data
template <typename FieldType>
ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& arcs,
                                            const IdArrayType& inSortOrder,
                                            const vtkm::cont::ArrayHandle<FieldType>& values,
                                            const IdArrayType& inGlobalMeshIndex)
316
  : SortOrder()
317
318
319
320
  , SortedValues()
  , GlobalMeshIndex(inGlobalMeshIndex)
  , Neighbours()
  , FirstNeighbour()
Sujin Philip's avatar
Sujin Philip committed
321
{
322
323
324
325
  this->NumVertices = inSortOrder.GetNumberOfValues();
  // Initalize the SortedIndices as a smart array handle
  this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
  this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
326
  // values permuted by SortOrder to sort the values
327
  auto permutedValues = vtkm::cont::make_ArrayHandlePermutation(inSortOrder, values);
Sujin Philip's avatar
Sujin Philip committed
328
  // TODO check if we actually need to make this copy here. we could just store the permutedValues array to save memory
329
  vtkm::cont::Algorithm::Copy(permutedValues, this->SortedValues);
Sujin Philip's avatar
Sujin Philip committed
330
  this->InitialiseNeighboursFromArcs(arcs);
331
332
333
334
335
336
337
338
#ifdef DEBUG_PRINT
  // Print the contents fo this for debugging
  DebugPrint("ContourTreeMesh Initialized", __FILE__, __LINE__);
#endif
}


template <typename FieldType>
339
340
341
342
343
inline ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& nodes,
                                                   const IdArrayType& arcs,
                                                   const IdArrayType& inSortOrder,
                                                   const vtkm::cont::ArrayHandle<FieldType>& values,
                                                   const IdArrayType& inGlobalMeshIndex)
344
345
346
  : GlobalMeshIndex(inGlobalMeshIndex)
  , Neighbours()
  , FirstNeighbour()
347
{
348
349
  // Initialize the SortedValues array the values permutted by the SortOrder permutted by the nodes, i.e.,
  // this->SortedValues[v] = values[inSortOrder[nodes[v]]];
350
351
352
  vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType> permutedSortOrder(nodes,
                                                                                 inSortOrder);
  auto permutedValues = vtkm::cont::make_ArrayHandlePermutation(permutedSortOrder, values);
353
  vtkm::cont::Algorithm::Copy(permutedValues, this->SortedValues);
354
  // Initalize the SortedIndices as a smart array handle
355
  this->NumVertices = this->SortedValues.GetNumberOfValues();
356
357
  this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
  this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
358
359
360
361
362
  this->InitialiseNeighboursFromArcs(arcs);
#ifdef DEBUG_PRINT
  // Print the contents fo this for debugging
  DebugPrint("ContourTreeMesh Initialized", __FILE__, __LINE__);
#endif
Sujin Philip's avatar
Sujin Philip committed
363
364
365
}

template <typename FieldType>
366
367
inline ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& arcs,
                                                   const ContourTreeMesh<FieldType>& mesh)
368
  : SortedValues(mesh.SortedValues)
369
370
371
  , GlobalMeshIndex(mesh.GlobalMeshIndex)
  , Neighbours()
  , FirstNeighbour()
Sujin Philip's avatar
Sujin Philip committed
372
{
373
  // Initalize the SortedIndices as a smart array handle
374
  this->NumVertices = this->SortedValues.GetNumberOfValues();
375
376
  this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
  this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
Sujin Philip's avatar
Sujin Philip committed
377
  this->InitialiseNeighboursFromArcs(arcs);
378
379
380
381
#ifdef DEBUG_PRINT
  // Print the contents fo this for debugging
  DebugPrint("ContourTreeMesh Initialized", __FILE__, __LINE__);
#endif
Sujin Philip's avatar
Sujin Philip committed
382
383
384
385
}


template <typename FieldType>
386
387
388
inline ContourTreeMesh<FieldType>::ContourTreeMesh(const IdArrayType& nodes,
                                                   const IdArrayType& arcs,
                                                   const ContourTreeMesh<FieldType>& mesh)
389
  : Neighbours()
390
  , FirstNeighbour()
Sujin Philip's avatar
Sujin Philip committed
391
{
392
  // Initatlize the global mesh index with the GlobalMeshIndex permutted by the nodes
Sujin Philip's avatar
Sujin Philip committed
393
  vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType> permutedGlobalMeshIndex(
394
395
396
397
398
    nodes, mesh.GlobalMeshIndex);
  vtkm::cont::Algorithm::Copy(permutedGlobalMeshIndex, this->GlobalMeshIndex);
  // Initialize the SortedValues array with the SortedValues permutted by the nodes
  auto permutedSortedValues = vtkm::cont::make_ArrayHandlePermutation(nodes, mesh.SortedValues);
  vtkm::cont::Algorithm::Copy(permutedSortedValues, this->SortedValues);
Sujin Philip's avatar
Sujin Philip committed
399
  // Initialize the neighbours from the arcs
400
  this->NumVertices = this->SortedValues.GetNumberOfValues();
401
402
  this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
  this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
Sujin Philip's avatar
Sujin Philip committed
403
  this->InitialiseNeighboursFromArcs(arcs);
404
405
406
407
#ifdef DEBUG_PRINT
  // Print the contents fo this for debugging
  DebugPrint("ContourTreeMesh Initialized", __FILE__, __LINE__);
#endif
Sujin Philip's avatar
Sujin Philip committed
408
409
410
411
412
}


// Initalize the contour tree from the arcs array and sort order
template <typename FieldType>
413
inline void ContourTreeMesh<FieldType>::InitialiseNeighboursFromArcs(const IdArrayType& arcs)
Sujin Philip's avatar
Sujin Philip committed
414
415
416
417
{
  // Find target indices for valid arcs in neighbours array ...
  IdArrayType arcTargetIndex;
  arcTargetIndex.Allocate(arcs.GetNumberOfValues());
418
  OneIfArcValid oneIfArcValidFunctor;
Sujin Philip's avatar
Sujin Philip committed
419
  auto oneIfArcValidArrayHandle =
420
    vtkm::cont::ArrayHandleTransform<IdArrayType, OneIfArcValid>(arcs, oneIfArcValidFunctor);
Sujin Philip's avatar
Sujin Philip committed
421
  vtkm::cont::Algorithm::ScanExclusive(oneIfArcValidArrayHandle, arcTargetIndex);
422
423
  vtkm::Id nValidArcs = arcTargetIndex.ReadPortal().Get(arcTargetIndex.GetNumberOfValues() - 1) +
    oneIfArcValidFunctor(arcs.ReadPortal().Get(arcs.GetNumberOfValues() - 1));
Sujin Philip's avatar
Sujin Philip committed
424
425

  // ... and compress array
426
427
  this->Neighbours.ReleaseResources();
  this->Neighbours.Allocate(2 * nValidArcs);
Sujin Philip's avatar
Sujin Philip committed
428
429

  contourtree_mesh_inc_ns::CompressNeighboursWorklet compressNeighboursWorklet;
430
  this->Invoke(compressNeighboursWorklet, arcs, arcTargetIndex, this->Neighbours);
Sujin Philip's avatar
Sujin Philip committed
431
432
433

  // Sort arcs so that all arcs from the same vertex are adjacent. All arcs are in neighbours array based on
  // sort index of their 'from' vertex (and then within a run sorted by sort index of their 'to' vertex).
434
  vtkm::cont::Algorithm::Sort(this->Neighbours, contourtree_mesh_inc_ns::ArcComparator(arcs));
Sujin Philip's avatar
Sujin Philip committed
435
436

  // Find start index for each vertex into neighbours array
437
  this->FirstNeighbour.Allocate(this->NumVertices);
Sujin Philip's avatar
Sujin Philip committed
438
439
440

  contourtree_mesh_inc_ns::FindStartIndexWorklet findStartIndexWorklet;
  this->Invoke(findStartIndexWorklet,
441
               this->Neighbours,
Sujin Philip's avatar
Sujin Philip committed
442
               arcs,
443
               this->FirstNeighbour // output
444
  );
Sujin Philip's avatar
Sujin Philip committed
445
446
447
448
449


  // Replace arc number with 'to' vertex in neighbours array
  contourtree_mesh_inc_ns::ReplaceArcNumWithToVertexWorklet replaceArcNumWithToVertexWorklet;
  this->Invoke(replaceArcNumWithToVertexWorklet,
450
               this->Neighbours, // input/output
Sujin Philip's avatar
Sujin Philip committed
451
               arcs              // input
452
  );
Sujin Philip's avatar
Sujin Philip committed
453
454
455
456
457
458
459

  // Compute maximum number of neighbours
  this->ComputeMaxNeighbours();

#ifdef DEBUG_PRINT
  std::cout << std::setw(30) << std::left << __FILE__ << ":" << std::right << std::setw(4)
            << __LINE__ << std::endl;
460
461
  auto firstNeighbourPortal = this->FirstNeighbour.ReadPortal();
  auto neighboursPortal = this->Neighbours.ReadPortal();
462
  for (vtkm::Id vtx = 0; vtx < FirstNeighbour.GetNumberOfValues(); ++vtx)
Sujin Philip's avatar
Sujin Philip committed
463
464
465
  {
    std::cout << vtx << ": ";
    vtkm::Id neighboursBeginIndex = firstNeighbourPortal.Get(vtx);
466
    vtkm::Id neighboursEndIndex = (vtx < this->NumVertices - 1) ? firstNeighbourPortal.Get(vtx + 1)
467
                                                                : Neighbours.GetNumberOfValues();
Sujin Philip's avatar
Sujin Philip committed
468
469
470
471
472
473
474

    for (vtkm::Id ni = neighboursBeginIndex; ni < neighboursEndIndex; ++ni)
    {
      std::cout << neighboursPortal.Get(ni) << " ";
    }
    std::cout << std::endl;
  }
475
  std::cout << "Max neighbours: " << this->MaxNeighbours << std::endl;
Sujin Philip's avatar
Sujin Philip committed
476
477
478
479
#endif
}

template <typename FieldType>
480
inline void ContourTreeMesh<FieldType>::ComputeNNeighboursVector(IdArrayType& nNeighbours) const
Sujin Philip's avatar
Sujin Philip committed
481
{
482
  nNeighbours.Allocate(this->FirstNeighbour.GetNumberOfValues()); // same as this->NumVertices
Sujin Philip's avatar
Sujin Philip committed
483
  contourtree_mesh_inc_ns::ComputeMaxNeighboursWorklet computeMaxNeighboursWorklet(
484
485
    this->Neighbours.GetNumberOfValues());
  this->Invoke(computeMaxNeighboursWorklet, this->FirstNeighbour, nNeighbours);
Sujin Philip's avatar
Sujin Philip committed
486
487
488
}

template <typename FieldType>
489
inline void ContourTreeMesh<FieldType>::ComputeMaxNeighbours()
Sujin Philip's avatar
Sujin Philip committed
490
491
492
493
494
{
  // Compute maximum number of neighbours
  IdArrayType nNeighbours;
  this->ComputeNNeighboursVector(nNeighbours);
  vtkm::cont::ArrayHandle<vtkm::Range> rangeArray = vtkm::cont::ArrayRangeCompute(nNeighbours);
495
  this->MaxNeighbours = static_cast<vtkm::Id>(rangeArray.ReadPortal().Get(0).Max);
Sujin Philip's avatar
Sujin Philip committed
496
497
498
499
}

// Define the behavior for the execution object generate by the PrepareForExecution function
template <typename FieldType>
500
inline void ContourTreeMesh<FieldType>::SetPrepareForExecutionBehavior(bool getMax)
Sujin Philip's avatar
Sujin Philip committed
501
502
503
504
505
506
{
  this->mGetMax = getMax;
}

// Get VTKM execution object that represents the structure of the mesh and provides the mesh helper functions on the device
template <typename FieldType>
507
508
509
contourtree_mesh_inc_ns::MeshStructureContourTreeMesh inline ContourTreeMesh<
  FieldType>::PrepareForExecution(vtkm::cont::DeviceAdapterId device,
                                  vtkm::cont::Token& token) const
Sujin Philip's avatar
Sujin Philip committed
510
{
511
512
  return contourtree_mesh_inc_ns::MeshStructureContourTreeMesh(
    this->Neighbours, this->FirstNeighbour, this->MaxNeighbours, this->mGetMax, device, token);
Sujin Philip's avatar
Sujin Philip committed
513
514
515
516
517
518
519
520
521
}

struct NotNoSuchElement
{
  VTKM_EXEC_CONT bool operator()(vtkm::Id x) const { return x != NO_SUCH_ELEMENT; }
};

// Combine two ContourTreeMeshes
template <typename FieldType>
522
inline void ContourTreeMesh<FieldType>::MergeWith(ContourTreeMesh<FieldType>& other)
523
{ // Merge With
Sujin Philip's avatar
Sujin Philip committed
524
525
526
527
528
529
#ifdef DEBUG_PRINT
  this->DebugPrint("THIS ContourTreeMesh", __FILE__, __LINE__);
  other.DebugPrint("OTHER ContourTreeMesh", __FILE__, __LINE__);
#endif

  // Create combined sort order
530
531
  // TODO This vector could potentially be implemented purely as a smart array handle to reduce memory usage
  IdArrayType overallSortOrder;
532
  overallSortOrder.Allocate(this->NumVertices + other.NumVertices);
Sujin Philip's avatar
Sujin Philip committed
533
534

  { // Create a new scope so that the following two vectors get deleted when leaving the scope
535
    auto thisIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices); // A regular index array
536
    MarkOther markOtherFunctor;
Sujin Philip's avatar
Sujin Philip committed
537
    auto otherIndices = vtkm::cont::make_ArrayHandleTransform(
538
      vtkm::cont::ArrayHandleIndex(other.NumVertices), markOtherFunctor);
539
540
541
542
543
544
545
546
547
    contourtree_mesh_inc_ns::CombinedSimulatedSimplicityIndexComparator<FieldType>
      cssicFunctorExecObj(
        this->GlobalMeshIndex, other.GlobalMeshIndex, this->SortedValues, other.SortedValues);
    // TODO FIXME We here need to force the arrays for the comparator onto the CPU by using DeviceAdapterTagSerial
    //            Instead we should implement the merge of the arrays on the device and not use std::merge
    vtkm::cont::Token tempToken;
    auto cssicFunctor =
      cssicFunctorExecObj.PrepareForExecution(vtkm::cont::DeviceAdapterTagSerial(), tempToken);
    // Merge the arrays
548
549
550
551
552
    std::merge(vtkm::cont::ArrayPortalToIteratorBegin(thisIndices.ReadPortal()),
               vtkm::cont::ArrayPortalToIteratorEnd(thisIndices.ReadPortal()),
               vtkm::cont::ArrayPortalToIteratorBegin(otherIndices.ReadPortal()),
               vtkm::cont::ArrayPortalToIteratorEnd(otherIndices.ReadPortal()),
               vtkm::cont::ArrayPortalToIteratorBegin(overallSortOrder.WritePortal()),
Sujin Philip's avatar
Sujin Philip committed
553
554
555
556
557
               cssicFunctor);
  }

#ifdef DEBUG_PRINT
  std::cout << "OverallSortOrder.size  " << overallSortOrder.GetNumberOfValues() << std::endl;
558
  PrintIndices("overallSortOrder", overallSortOrder);
Sujin Philip's avatar
Sujin Philip committed
559
560
561
562
563
564
  std::cout << std::endl;
#endif

  IdArrayType overallSortIndex;
  overallSortIndex.Allocate(overallSortOrder.GetNumberOfValues());
  {
565
566
567
568
569
570
571
572
    // Array decorator with functor returning 0, 1 for each element depending
    // on whethern the current value is different from the next
    auto differentFromNextArr = vtkm::cont::make_ArrayHandleDecorator(
      overallSortIndex.GetNumberOfValues() - 1,
      mesh_dem_contourtree_mesh_inc::CombinedVectorDifferentFromNextDecoratorImpl{},
      overallSortOrder,
      this->GlobalMeshIndex,
      other.GlobalMeshIndex);
Sujin Philip's avatar
Sujin Philip committed
573
574

    // Compute the exclusive scan of our transformed combined vector
575
    overallSortIndex.WritePortal().Set(0, 0);
Sujin Philip's avatar
Sujin Philip committed
576
    IdArrayType tempArr;
577
578

    vtkm::cont::Algorithm::ScanInclusive(differentFromNextArr, tempArr);
Sujin Philip's avatar
Sujin Philip committed
579
    vtkm::cont::Algorithm::CopySubRange(
580
      tempArr, 0, tempArr.GetNumberOfValues(), overallSortIndex, 1);
Sujin Philip's avatar
Sujin Philip committed
581
  }
582
  vtkm::Id numVerticesCombined =
583
    overallSortIndex.ReadPortal().Get(overallSortIndex.GetNumberOfValues() - 1) + 1;
584

Sujin Philip's avatar
Sujin Philip committed
585
586
#ifdef DEBUG_PRINT
  std::cout << "OverallSortIndex.size  " << overallSortIndex.GetNumberOfValues() << std::endl;
587
  PrintIndices("overallSortIndex", overallSortIndex);
588
  std::cout << "numVerticesCombined: " << numVerticesCombined << std::endl;
Sujin Philip's avatar
Sujin Philip committed
589
590
591
592
593
  std::cout << std::endl;
#endif

  // thisToCombinedSortOrder and otherToCombinedSortOrder
  IdArrayType thisToCombinedSortOrder;
594
  thisToCombinedSortOrder.Allocate(this->FirstNeighbour.GetNumberOfValues());
Sujin Philip's avatar
Sujin Philip committed
595
  IdArrayType otherToCombinedSortOrder;
596
  otherToCombinedSortOrder.Allocate(other.FirstNeighbour.GetNumberOfValues());
Sujin Philip's avatar
Sujin Philip committed
597
598
599
600
601
602
603
604
605
  contourtree_mesh_inc_ns::InitToCombinedSortOrderArraysWorklet
    initToCombinedSortOrderArraysWorklet;
  this->Invoke(initToCombinedSortOrderArraysWorklet,
               overallSortIndex,
               overallSortOrder,
               thisToCombinedSortOrder,
               otherToCombinedSortOrder);

#ifdef DEBUG_PRINT
606
607
  PrintIndices("thisToCombinedSortOrder", thisToCombinedSortOrder);
  PrintIndices("otherToCombinedSortOrder", otherToCombinedSortOrder);
Sujin Philip's avatar
Sujin Philip committed
608
609
610
#endif

  IdArrayType combinedNNeighbours;
611
  vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, numVerticesCombined),
Sujin Philip's avatar
Sujin Philip committed
612
613
614
615
616
617
618
619
620
621
                              combinedNNeighbours);
  { // New scope so that array gets deleted when leaving scope
    IdArrayType nNeighbours;
    this->ComputeNNeighboursVector(nNeighbours);
    auto permutedCombinedNNeighbours =
      vtkm::cont::make_ArrayHandlePermutation(thisToCombinedSortOrder, combinedNNeighbours);
    vtkm::cont::Algorithm::Copy(nNeighbours, permutedCombinedNNeighbours);
  }

  IdArrayType combinedOtherStartIndex;
622
  vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, numVerticesCombined),
Sujin Philip's avatar
Sujin Philip committed
623
624
625
626
627
628
629
630
631
632
633
                              combinedOtherStartIndex);
  { // New scope so that array gets deleted when leaving scope
    IdArrayType nNeighbours;
    other.ComputeNNeighboursVector(nNeighbours);
    contourtree_mesh_inc_ns::CombinedOtherStartIndexNNeighboursWorklet
      combinedOtherStartIndexNNeighboursWorklet;
    this->Invoke(combinedOtherStartIndexNNeighboursWorklet,
                 nNeighbours,              // input
                 otherToCombinedSortOrder, // input
                 combinedNNeighbours,      // input/output
                 combinedOtherStartIndex   // input/output
634
    );
Sujin Philip's avatar
Sujin Philip committed
635
636
637
  }

#ifdef DEBUG_PRINT
638
639
  PrintIndices("combinedNNeighbours", combinedNNeighbours);
  PrintIndices("combinedOtherStartIndex", combinedOtherStartIndex);
Sujin Philip's avatar
Sujin Philip committed
640
641
642
#endif

  IdArrayType combinedFirstNeighbour;
643
  combinedFirstNeighbour.Allocate(numVerticesCombined);
Sujin Philip's avatar
Sujin Philip committed
644
  vtkm::cont::Algorithm::ScanExclusive(combinedNNeighbours, combinedFirstNeighbour);
645
646
647
  vtkm::Id nCombinedNeighbours =
    combinedFirstNeighbour.ReadPortal().Get(combinedFirstNeighbour.GetNumberOfValues() - 1) +
    combinedNNeighbours.ReadPortal().Get(combinedNNeighbours.GetNumberOfValues() - 1);
Sujin Philip's avatar
Sujin Philip committed
648
649
650
651
652
653
654
655
656

  IdArrayType combinedNeighbours;
  combinedNeighbours.Allocate(nCombinedNeighbours);

  // Update combined neighbours
  contourtree_mesh_inc_ns::UpdateCombinedNeighboursWorklet updateCombinedNeighboursWorklet;
  // Updata neighbours from this
  this->Invoke(
    updateCombinedNeighboursWorklet,
657
658
    this->FirstNeighbour,
    this->Neighbours,
Sujin Philip's avatar
Sujin Philip committed
659
660
661
662
    thisToCombinedSortOrder,
    combinedFirstNeighbour,
    vtkm::cont::ArrayHandleConstant<vtkm::Id>(
      0,
663
      numVerticesCombined), // Constant 0 array. Just needed so we can use the same worklet for both cases
Sujin Philip's avatar
Sujin Philip committed
664
665
666
    combinedNeighbours);
  // Update neighbours from other
  this->Invoke(updateCombinedNeighboursWorklet,
667
668
               other.FirstNeighbour,
               other.Neighbours,
Sujin Philip's avatar
Sujin Philip committed
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
               otherToCombinedSortOrder,
               combinedFirstNeighbour,
               combinedOtherStartIndex,
               combinedNeighbours);

  // TODO VTKM -Version MergedCombinedOtherStartIndex. Replace 1r block with the 1s block. Need to check for Segfault in contourtree_mesh_inc_ns::MergeCombinedOtherStartIndexWorklet. This workloat also still uses a number of stl algorithms that should be replaced with VTKm code (which is porbably also why the worklet fails).
  /* // 1s--start
    contourtree_mesh_inc_ns::MergeCombinedOtherStartIndexWorklet<DeviceAdapter> mergeCombinedOtherStartIndexWorklet;
    vtkm::worklet::DispatcherMapField< contourtree_mesh_inc_ns::MergeCombinedOtherStartIndexWorklet<DeviceAdapter>> mergeCombinedOtherStartIndexDispatcher(mergeCombinedOtherStartIndexWorklet);
    this->Invoke(mergeCombinedOtherStartIndexWorklet,
       combinedOtherStartIndex, // (input/output)
       combinedNeighbours,      // (input/output)
       combinedFirstNeighbour   // (input)
       );
    // 1s--end
    */

  // TODO Fix the MergedCombinedOtherStartIndex worklet and remove //1r block below
  // 1r--start
688
689
690
  auto combinedOtherStartIndexPortal = combinedOtherStartIndex.WritePortal();
  auto combinedFirstNeighbourPortal = combinedFirstNeighbour.ReadPortal();
  auto combinedNeighboursPortal = combinedNeighbours.WritePortal();
Sujin Philip's avatar
Sujin Philip committed
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
  std::vector<vtkm::Id> tempCombinedNeighours((std::size_t)combinedNeighbours.GetNumberOfValues());
  for (vtkm::Id vtx = 0; vtx < combinedNeighbours.GetNumberOfValues(); ++vtx)
  {
    tempCombinedNeighours[(std::size_t)vtx] = combinedNeighboursPortal.Get(vtx);
  }
  for (vtkm::Id vtx = 0; vtx < combinedFirstNeighbour.GetNumberOfValues(); ++vtx)
  {
    if (combinedOtherStartIndexPortal.Get(vtx)) // Needs merge
    {
      auto neighboursBegin = tempCombinedNeighours.begin() + combinedFirstNeighbourPortal.Get(vtx);
      auto neighboursEnd = (vtx < combinedFirstNeighbour.GetNumberOfValues() - 1)
        ? tempCombinedNeighours.begin() + combinedFirstNeighbourPortal.Get(vtx + 1)
        : tempCombinedNeighours.end();
      std::inplace_merge(
        neighboursBegin, neighboursBegin + combinedOtherStartIndexPortal.Get(vtx), neighboursEnd);
      auto it = std::unique(neighboursBegin, neighboursEnd);
707
      combinedOtherStartIndexPortal.Set(vtx, static_cast<vtkm::Id>(neighboursEnd - it));
Sujin Philip's avatar
Sujin Philip committed
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
      while (it != neighboursEnd)
        *(it++) = NO_SUCH_ELEMENT;
    }
  }
  // copy the values back to VTKm
  for (vtkm::Id vtx = 0; vtx < combinedNeighbours.GetNumberOfValues(); ++vtx)
  {
    combinedNeighboursPortal.Set(vtx, tempCombinedNeighours[(std::size_t)vtx]);
  }
  // 1r--end

  IdArrayType combinedFirstNeighbourShift;
  combinedFirstNeighbourShift.Allocate(combinedFirstNeighbour.GetNumberOfValues());
  vtkm::cont::Algorithm::ScanExclusive(combinedOtherStartIndex, combinedFirstNeighbourShift);

  {
    IdArrayType tempCombinedNeighbours;
    vtkm::cont::Algorithm::CopyIf(
      combinedNeighbours, combinedNeighbours, tempCombinedNeighbours, NotNoSuchElement());
    combinedNeighbours = tempCombinedNeighbours; // Swap the two arrays
  }

  // Adjust firstNeigbour indices by deleted elements
  { // make sure variables created are deleted after the context
    contourtree_mesh_inc_ns::SubtractAssignWorklet subAssignWorklet;
    this->Invoke(subAssignWorklet, combinedFirstNeighbour, combinedFirstNeighbourShift);
  }

  // Compute combined global mesh index arrays
  IdArrayType combinedGlobalMeshIndex;
  combinedGlobalMeshIndex.Allocate(combinedFirstNeighbour.GetNumberOfValues());
  { // make sure arrays used for copy go out of scope
    auto permutedCombinedGlobalMeshIndex =
      vtkm::cont::make_ArrayHandlePermutation(thisToCombinedSortOrder, combinedGlobalMeshIndex);
742
    vtkm::cont::Algorithm::Copy(GlobalMeshIndex, permutedCombinedGlobalMeshIndex);
Sujin Philip's avatar
Sujin Philip committed
743
744
745
746
  }
  { // make sure arrays used for copy go out of scope
    auto permutedCombinedGlobalMeshIndex =
      vtkm::cont::make_ArrayHandlePermutation(otherToCombinedSortOrder, combinedGlobalMeshIndex);
747
    vtkm::cont::Algorithm::Copy(other.GlobalMeshIndex, permutedCombinedGlobalMeshIndex);
Sujin Philip's avatar
Sujin Philip committed
748
749
750
751
752
753
754
755
  }

  // Compute combined sorted values
  vtkm::cont::ArrayHandle<FieldType> combinedSortedValues;
  combinedSortedValues.Allocate(combinedFirstNeighbour.GetNumberOfValues());
  { // make sure arrays used for copy go out of scope
    auto permutedCombinedSortedValues =
      vtkm::cont::make_ArrayHandlePermutation(thisToCombinedSortOrder, combinedSortedValues);
756
    vtkm::cont::Algorithm::Copy(SortedValues, permutedCombinedSortedValues);
Sujin Philip's avatar
Sujin Philip committed
757
758
759
760
  }
  { // make sure arrays used for copy go out of scope
    auto permutedCombinedSortedValues =
      vtkm::cont::make_ArrayHandlePermutation(otherToCombinedSortOrder, combinedSortedValues);
761
    vtkm::cont::Algorithm::Copy(other.SortedValues, permutedCombinedSortedValues);
Sujin Philip's avatar
Sujin Philip committed
762
763
764
  }

  // Swap in combined version. VTKM ArrayHandles are smart so we can just swap in the new for the old
765
766
767
768
769
  this->SortedValues = combinedSortedValues;
  this->GlobalMeshIndex = combinedGlobalMeshIndex;
  this->Neighbours = combinedNeighbours;
  this->FirstNeighbour = combinedFirstNeighbour;
  this->NumVertices = SortedValues.GetNumberOfValues();
770
771
  this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
  this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
Sujin Philip's avatar
Sujin Philip committed
772
773
774
775
776
777
778
779

  // Re-compute maximum number of neigbours
  ComputeMaxNeighbours();

#ifdef DEBUG_PRINT
  // Print the contents fo this for debugging
  DebugPrint("ContourTreeMeshes merged", __FILE__, __LINE__);
#endif
780
} // Merge With
Sujin Philip's avatar
Sujin Philip committed
781
782
783


template <typename FieldType>
784
inline void ContourTreeMesh<FieldType>::Save(const char* filename) const
Sujin Philip's avatar
Sujin Philip committed
785
786
{
  std::ofstream os(filename);
787
788
789
790
  SaveVector(os, this->SortedValues);
  SaveVector(os, this->GlobalMeshIndex);
  SaveVector(os, this->Neighbours);
  SaveVector(os, this->FirstNeighbour);
Sujin Philip's avatar
Sujin Philip committed
791
792
793
}

template <typename FieldType>
794
inline void ContourTreeMesh<FieldType>::Load(const char* filename)
Sujin Philip's avatar
Sujin Philip committed
795
796
{
  std::ifstream is(filename);
797
798
799
800
  if (!is.is_open())
  {
    throw vtkm::io::ErrorIO(std::string("Unable to open file: ") + std::string(filename));
  }
801
802
803
804
  LoadVector(is, this->SortedValues);
  LoadVector(is, this->GlobalMeshIndex);
  LoadVector(is, this->Neighbours);
  LoadVector(is, this->FirstNeighbour);
Sujin Philip's avatar
Sujin Philip committed
805
  this->ComputeMaxNeighbours();
806
807
808
  this->NumVertices = this->SortedValues.GetNumberOfValues();
  this->SortOrder = vtkm::cont::ArrayHandleIndex(this->NumVertices);
  this->SortIndices = vtkm::cont::ArrayHandleIndex(this->NumVertices);
Sujin Philip's avatar
Sujin Philip committed
809
810
811
812
}

template <typename FieldType>
template <typename ValueType>
813
814
815
inline void ContourTreeMesh<FieldType>::SaveVector(
  std::ostream& os,
  const vtkm::cont::ArrayHandle<ValueType>& vec) const
Sujin Philip's avatar
Sujin Philip committed
816
817
{
  vtkm::Id numVals = vec.GetNumberOfValues();
818
  //os.write(rXeinterpret_cast<const char*>(&numVals), sizeof(ValueType));
819
  os << numVals << ": ";
820
  auto vecPortal = vec.ReadPortal();
Sujin Philip's avatar
Sujin Philip committed
821
  for (vtkm::Id i = 0; i < numVals; ++i)
822
823
824
    os << vecPortal.Get(i) << " ";
  //os.write(reinterpret_cast<const char*>(vecPortal.Get(i)), sizeof(ValueType));
  os << std::endl;
Sujin Philip's avatar
Sujin Philip committed
825
826
827
828
}

template <typename FieldType>
template <typename ValueType>
829
inline void ContourTreeMesh<FieldType>::LoadVector(std::istream& is,
830
                                                   vtkm::cont::ArrayHandle<ValueType>& vec)
Sujin Philip's avatar
Sujin Philip committed
831
832
{
  vtkm::Id numVals;
833
834
835
836
837
838
839
  is >> numVals;
  char colon = is.get();
  if (colon != ':')
  {
    throw vtkm::io::ErrorIO("Error parsing file");
  }

Sujin Philip's avatar
Sujin Philip committed
840
  vec.Allocate(numVals);
841
  auto vecPortal = vec.WritePortal();
842
  ValueType val;
Sujin Philip's avatar
Sujin Philip committed
843
844
  for (vtkm::Id i = 0; i < numVals; ++i)
  {
845
    is >> val;
Sujin Philip's avatar
Sujin Philip committed
846
847
848
849
    vecPortal.Set(i, val);
  }
}

850
template <typename FieldType>
851
inline MeshBoundaryContourTreeMeshExec ContourTreeMesh<FieldType>::GetMeshBoundaryExecutionObject(
852
  vtkm::Id3 globalSize,
853
854
855
  vtkm::Id3 minIdx,
  vtkm::Id3 maxIdx) const
{
856
  return MeshBoundaryContourTreeMeshExec(this->GlobalMeshIndex, globalSize, minIdx, maxIdx);
857
858
}

859
template <typename FieldType>
860
inline void ContourTreeMesh<FieldType>::GetBoundaryVertices(
861
862
863
  IdArrayType& boundaryVertexArray,                    // output
  IdArrayType& boundarySortIndexArray,                 // output
  MeshBoundaryContourTreeMeshExec* meshBoundaryExecObj //input
864
) const
865
866
867
868
869
870
871
872
873
874
{
  // start by generating a temporary array of indices
  auto indexArray = vtkm::cont::ArrayHandleIndex(this->GlobalMeshIndex.GetNumberOfValues());
  // compute the boolean array indicating which values lie on the boundary
  vtkm::cont::ArrayHandle<bool> isOnBoundary;
  ComputeMeshBoundaryContourTreeMesh computeMeshBoundaryContourTreeMeshWorklet;
  this->Invoke(computeMeshBoundaryContourTreeMeshWorklet,
               indexArray,           // input
               *meshBoundaryExecObj, // input
               isOnBoundary          // outut
875
  );
876
877
878
879
880
881
882

  // we will conditionally copy the boundary vertices' indices, capturing the end iterator to compute the # of boundary vertices
  vtkm::cont::Algorithm::CopyIf(indexArray, isOnBoundary, boundaryVertexArray);
  // duplicate these into the index array, since the BRACT uses indices into the underlying mesh anyway
  vtkm::cont::Algorithm::Copy(boundaryVertexArray, boundarySortIndexArray);
}

Sujin Philip's avatar
Sujin Philip committed
883
884
885
886
887
} // namespace contourtree_augmented
} // worklet
} // vtkm

#endif