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 
11 #include <vtkm/cont/DataSet.h>
12 #include <vtkm/cont/Initialize.h>
13 #include <vtkm/filter/Streamline.h>
14 #include <vtkm/io/VTKDataSetReader.h>
15 #include <vtkm/io/VTKDataSetWriter.h>
16 
17 // Example computing streamlines.
18 // An example vector field is available in the vtk-m data directory: magField.vtk
19 // Example usage:
20 //   this will advect 200 particles 50 steps using a step size of 0.01
21 //
22 // Particle_Advection <path-to-data-dir>/magField.vtk vec 200 50 0.01 output.vtk
23 //
24 
main(int argc,char ** argv)25 int main(int argc, char** argv)
26 {
27   auto opts = vtkm::cont::InitializeOptions::DefaultAnyDevice;
28   auto config = vtkm::cont::Initialize(argc, argv, opts);
29 
30   if (argc < 8)
31   {
32     std::cerr << "Usage: " << argv[0]
33               << "dataFile varName numSeeds numSteps stepSize outputFile [options]" << std::endl;
34     std::cerr << "where options are: " << std::endl << config.Usage << std::endl;
35     return -1;
36   }
37 
38   std::string dataFile = argv[1];
39   std::string varName = argv[2];
40   vtkm::Id numSeeds = std::stoi(argv[3]);
41   vtkm::Id numSteps = std::stoi(argv[4]);
42   vtkm::FloatDefault stepSize = std::stof(argv[5]);
43   std::string outputFile = argv[6];
44 
45   vtkm::cont::DataSet ds;
46 
47   if (dataFile.find(".vtk") != std::string::npos)
48   {
49     vtkm::io::VTKDataSetReader rdr(dataFile);
50     ds = rdr.ReadDataSet();
51   }
52   else
53   {
54     std::cerr << "Unsupported data file: " << dataFile << std::endl;
55     return -1;
56   }
57 
58   //create seeds randomly placed withing the bounding box of the data.
59   vtkm::Bounds bounds = ds.GetCoordinateSystem().GetBounds();
60   std::vector<vtkm::Particle> seeds;
61 
62   for (vtkm::Id i = 0; i < numSeeds; i++)
63   {
64     vtkm::Particle p;
65     vtkm::FloatDefault rx = (vtkm::FloatDefault)rand() / (vtkm::FloatDefault)RAND_MAX;
66     vtkm::FloatDefault ry = (vtkm::FloatDefault)rand() / (vtkm::FloatDefault)RAND_MAX;
67     vtkm::FloatDefault rz = (vtkm::FloatDefault)rand() / (vtkm::FloatDefault)RAND_MAX;
68     p.Pos[0] = static_cast<vtkm::FloatDefault>(bounds.X.Min + rx * bounds.X.Length());
69     p.Pos[1] = static_cast<vtkm::FloatDefault>(bounds.Y.Min + ry * bounds.Y.Length());
70     p.Pos[2] = static_cast<vtkm::FloatDefault>(bounds.Z.Min + rz * bounds.Z.Length());
71     p.ID = i;
72     seeds.push_back(p);
73   }
74   auto seedArray = vtkm::cont::make_ArrayHandle(seeds, vtkm::CopyFlag::Off);
75 
76   //compute streamlines
77   vtkm::filter::Streamline streamline;
78 
79   streamline.SetStepSize(stepSize);
80   streamline.SetNumberOfSteps(numSteps);
81   streamline.SetSeeds(seedArray);
82 
83   streamline.SetActiveField(varName);
84   auto output = streamline.Execute(ds);
85 
86   vtkm::io::VTKDataSetWriter wrt(outputFile);
87   wrt.WriteDataSet(output);
88 
89   return 0;
90 }
91