Streamlines from vtkm

src/tests/vtkh/CMakeLists.txt

@@ -93,11 +93,12 @@ set(MPI_TESTS t_vtk-h_smoke_par
               t_vtk-h_statistics_par
               t_vtk-h_marching_cubes_par
               t_vtk-h_multi_render_par
+              t_vtk-h_particle_advection_par
               t_vtk-h_scalar_renderer_par
               t_vtk-h_raytracer_par
               t_vtk-h_sampling_par
               t_vtk-h_volume_renderer_par
-              
+
               )
 
 # Contour tree depends on the VTK-m build type.

src/tests/vtkh/t_vtk-h_particle_advection_par.cpp

@@ -0,0 +1,145 @@
+//-----------------------------------------------------------------------------
+///
+/// file: t_vtk-h_particle_advection_par.cpp
+///
+//-----------------------------------------------------------------------------
+
+#include "gtest/gtest.h"
+
+#include <vtkh/vtkh.hpp>
+#include <vtkh/DataSet.hpp>
+#include <vtkh/filters/ParticleAdvection.hpp>
+#include <vtkh/filters/Streamline.hpp>
+#include <vtkm/io/writer/VTKDataSetWriter.h>
+#include <vtkm/cont/DataSet.h>
+#include <vtkm/cont/DataSetFieldAdd.h>
+#include <vtkm/cont/CellSetSingleType.h>
+#include "t_test_utils.hpp"
+#include <iostream>
+#include <mpi.h>
+
+void checkValidity(vtkh::DataSet *data, const int maxSteps, bool isSL)
+{
+  int numDomains = data->GetNumberOfDomains();
+
+  //Check all domains
+  for(int i = 0; i < numDomains; i++)
+  {
+    auto currentDomain = data->GetDomain(i);
+    auto cs = currentDomain.GetCellSet();
+    if (isSL)
+    {
+      auto cellSet = cs.Cast<vtkm::cont::CellSetExplicit<>>();
+      //Ensure that streamlines took <= to the max number of steps
+      for(int j = 0; j < cellSet.GetNumberOfCells(); j++)
+      {
+        EXPECT_LE(cellSet.GetNumberOfPointsInCell(j), maxSteps);
+      }
+    }
+    else
+    {
+      if (!cs.IsType<vtkm::cont::CellSetSingleType<>>())
+        EXPECT_TRUE(false);
+    }
+  }
+}
+
+void writeDataSet(vtkh::DataSet *data, std::string fName, int rank)
+{
+  int numDomains = data->GetNumberOfDomains();
+  std::cerr << "num domains " << numDomains << std::endl;
+  for(int i = 0; i < numDomains; i++)
+  {
+    char fileNm[128];
+    sprintf(fileNm, "%s.rank%d.domain%d.vtk", fName.c_str(), rank, i);
+    vtkm::io::writer::VTKDataSetWriter write(fileNm);
+    write.WriteDataSet(data->GetDomain(i));
+  }
+}
+
+static inline vtkm::FloatDefault
+rand01()
+{
+  return (vtkm::FloatDefault)rand() / (RAND_MAX+1.0f);
+}
+
+static inline vtkm::FloatDefault
+randRange(const vtkm::FloatDefault &a, const vtkm::FloatDefault &b)
+{
+    return a + (b-a)*rand01();
+}
+
+template <typename FilterType>
+vtkh::DataSet *
+RunFilter(vtkh::DataSet& input,
+          const std::string& fieldName,
+          const std::vector<vtkm::Particle>& seeds,
+          int maxAdvSteps,
+          double stepSize)
+{
+  FilterType filter;
+
+  filter.SetInput(&input);
+  filter.SetField(fieldName);
+  filter.SetNumberOfSteps(maxAdvSteps);
+  filter.SetStepSize(stepSize);
+  filter.SetSeeds(seeds);
+  filter.Update();
+
+  return filter.GetOutput();
+}
+
+//----------------------------------------------------------------------------
+TEST(vtkh_particle_advection, vtkh_serial_particle_advection)
+{
+  MPI_Init(NULL, NULL);
+  int comm_size, rank;
+  MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  vtkh::SetMPICommHandle(MPI_Comm_c2f(MPI_COMM_WORLD));
+
+  std::cout << "Running parallel Particle Advection, vtkh - with " << comm_size << " ranks" << std::endl;
+
+  vtkh::DataSet data_set;
+  const int base_size = 32;
+  const int blocks_per_rank = 1;
+  const int maxAdvSteps = 1000;
+  const int num_blocks = comm_size * blocks_per_rank;
+
+  for(int i = 0; i < blocks_per_rank; ++i)
+  {
+    int domain_id = rank * blocks_per_rank + i;
+    data_set.AddDomain(CreateTestDataRectilinear(domain_id, num_blocks, base_size), domain_id);
+  }
+
+  std::vector<vtkm::Particle> seeds;
+
+  vtkm::Bounds bounds = data_set.GetGlobalBounds();
+  std::cout<<"Bounds= "<<bounds<<std::endl;
+
+  for (int i = 0; i < 100; i++)
+  {
+    vtkm::Particle p;
+    p.Pos = vtkm::Vec3f(randRange(bounds.X.Min, bounds.X.Max),
+                        randRange(bounds.Y.Min, bounds.Y.Max),
+                        randRange(bounds.Z.Min, bounds.Z.Max));
+    p.ID = static_cast<vtkm::Id>(i);
+
+    seeds.push_back(p);
+  }
+
+  vtkh::DataSet *outPA=NULL, *outSL=NULL;
+
+  outPA = RunFilter<vtkh::ParticleAdvection>(data_set, "vector_data_Float64", seeds, maxAdvSteps, 0.1);
+  outPA->PrintSummary(std::cout);
+  checkValidity(outPA, maxAdvSteps+1, false);
+
+  outSL = RunFilter<vtkh::Streamline>(data_set, "vector_data_Float64", seeds, maxAdvSteps, 0.1);
+  outSL->PrintSummary(std::cout);
+  checkValidity(outSL, maxAdvSteps+1, true);
+
+  writeDataSet(outSL, "advection_SeedsRandomWhole", rank);
+
+  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Finalize();
+}

src/vtkh/filters/CMakeLists.txt

@@ -18,13 +18,15 @@ set(vtkh_filters_headers
   NoOp.hpp
   Lagrangian.hpp
   MarchingCubes.hpp
+  ParticleAdvection.hpp
   PointAverage.hpp
   PointTransform.hpp
   Recenter.hpp
   Tetrahedralize.hpp
   Threshold.hpp
   Triangulate.hpp
   Statistics.hpp
+  Streamline.hpp
   Slice.hpp
   VectorComponent.hpp
   VectorMagnitude.hpp
@@ -59,6 +61,7 @@ set(vtkh_filters_sources
   NoOp.cpp
   Lagrangian.cpp
   MarchingCubes.cpp
+  ParticleAdvection.cpp
   PointAverage.cpp
   PointTransform.cpp
   Recenter.cpp
@@ -67,6 +70,7 @@ set(vtkh_filters_sources
   Triangulate.cpp
   Slice.cpp
   Statistics.cpp
+  Streamline.cpp
   VectorComponent.cpp
   VectorMagnitude.cpp
   communication/MemStream.cpp
@@ -76,7 +80,6 @@ if(ENABLE_FILTER_CONTOUR_TREE)
   list(APPEND vtkh_filters_sources ContourTree.cpp)
 endif()
 
-
 set(vtkh_comm_filters_sources
   communication/MemStream.cpp
   )

src/vtkh/filters/ParticleAdvection.cpp

@@ -0,0 +1,76 @@
+#include <iostream>
+#include <vtkh/filters/ParticleAdvection.hpp>
+#include <vtkm/filter/ParticleAdvection.h>
+#include <vtkh/vtkh.hpp>
+#include <vtkh/Error.hpp>
+
+namespace vtkh
+{
+
+ParticleAdvection::ParticleAdvection()
+{
+}
+
+ParticleAdvection::~ParticleAdvection()
+{
+
+}
+
+void ParticleAdvection::PreExecute()
+{
+  Filter::PreExecute();
+  Filter::CheckForRequiredField(m_field_name);
+}
+
+void ParticleAdvection::PostExecute()
+{
+  Filter::PostExecute();
+}
+
+void ParticleAdvection::DoExecute()
+{
+  this->m_output = new DataSet();
+#ifndef VTKH_BYPASS_VTKM_BIH
+  const int num_domains = this->m_input->GetNumberOfDomains();
+
+  vtkm::cont::PartitionedDataSet inputs;
+  for(int i = 0; i < num_domains; ++i)
+  {
+    vtkm::Id domain_id;
+    vtkm::cont::DataSet dom;
+    this->m_input->GetDomain(i, dom, domain_id);
+    if(dom.HasField(m_field_name))
+    {
+      using vectorField_d = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 3>>;
+      using vectorField_f = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>>;
+      auto field = dom.GetField(m_field_name).GetData();
+      if(!field.IsType<vectorField_d>() && !field.IsType<vectorField_f>())
+      {
+        throw Error("Vector field type does not match <vtkm::Vec<vtkm::Float32,3>> or <vtkm::Vec<vtkm::Float64,3>>");
+      }
+    }
+    else
+    {
+      throw Error("Domain does not contain specified vector field for ParticleAdvection analysis.");
+    }
+
+    inputs.AppendPartition(dom);
+  }
+
+  vtkm::filter::ParticleAdvection particleAdvectionFilter;
+  auto seedsAH = vtkm::cont::make_ArrayHandle(m_seeds, vtkm::CopyFlag::Off);
+
+  particleAdvectionFilter.SetStepSize(m_step_size);
+  particleAdvectionFilter.SetActiveField(m_field_name);
+  particleAdvectionFilter.SetSeeds(seedsAH);
+  particleAdvectionFilter.SetNumberOfSteps(m_num_steps);
+  auto out = particleAdvectionFilter.Execute(inputs);
+
+  for (vtkm::Id i = 0; i < out.GetNumberOfPartitions(); i++)
+  {
+    this->m_output->AddDomain(out.GetPartition(i), i);
+  }
+#endif
+}
+
+} //  namespace vtkh

src/vtkh/filters/ParticleAdvection.hpp

@@ -0,0 +1,37 @@
+#ifndef VTK_H_PARTICLE_ADVECTION_HPP
+#define VTK_H_PARTICLE_ADVECTION_HPP
+
+#include <vtkh/vtkh_exports.h>
+#include <vtkh/vtkh.hpp>
+#include <vtkh/filters/Filter.hpp>
+#include <vtkh/DataSet.hpp>
+
+#include <vtkm/Particle.h>
+
+namespace vtkh
+{
+
+class VTKH_API ParticleAdvection : public Filter
+{
+public:
+  ParticleAdvection();
+  virtual ~ParticleAdvection();
+  std::string GetName() const override { return "vtkh::ParticleAdvection";}
+  void SetField(const std::string &field_name) {  m_field_name = field_name; }
+  void SetStepSize(const double &step_size) {   m_step_size = step_size; }
+  void SetSeeds(const std::vector<vtkm::Particle>& seeds) { m_seeds = seeds; }
+  void SetNumberOfSteps(int numSteps) { m_num_steps = numSteps; }
+
+protected:
+  void PreExecute() override;
+  void PostExecute() override;
+  void DoExecute() override;
+
+  std::string m_field_name;
+  double m_step_size;
+  int m_num_steps;
+  std::vector<vtkm::Particle> m_seeds;
+};
+
+} //namespace vtkh
+#endif

src/vtkh/filters/Streamline.cpp

@@ -0,0 +1,76 @@
+#include <iostream>
+#include <vtkh/filters/Streamline.hpp>
+#include <vtkm/filter/Streamline.h>
+#include <vtkh/vtkh.hpp>
+#include <vtkh/Error.hpp>
+
+namespace vtkh
+{
+
+Streamline::Streamline()
+{
+}
+
+Streamline::~Streamline()
+{
+
+}
+
+void Streamline::PreExecute()
+{
+  Filter::PreExecute();
+  Filter::CheckForRequiredField(m_field_name);
+}
+
+void Streamline::PostExecute()
+{
+  Filter::PostExecute();
+}
+
+void Streamline::DoExecute()
+{
+  this->m_output = new DataSet();
+#ifndef VTKH_BYPASS_VTKM_BIH
+  const int num_domains = this->m_input->GetNumberOfDomains();
+
+  vtkm::cont::PartitionedDataSet inputs;
+  for(int i = 0; i < num_domains; ++i)
+  {
+    vtkm::Id domain_id;
+    vtkm::cont::DataSet dom;
+    this->m_input->GetDomain(i, dom, domain_id);
+    if(dom.HasField(m_field_name))
+    {
+      using vectorField_d = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 3>>;
+      using vectorField_f = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>>;
+      auto field = dom.GetField(m_field_name).GetData();
+      if(!field.IsType<vectorField_d>() && !field.IsType<vectorField_f>())
+      {
+        throw Error("Vector field type does not match <vtkm::Vec<vtkm::Float32,3>> or <vtkm::Vec<vtkm::Float64,3>>");
+      }
+    }
+    else
+    {
+      throw Error("Domain does not contain specified vector field for Streamline analysis.");
+    }
+
+    inputs.AppendPartition(dom);
+  }
+
+  vtkm::filter::Streamline streamlineFilter;
+  auto seedsAH = vtkm::cont::make_ArrayHandle(m_seeds, vtkm::CopyFlag::Off);
+
+  streamlineFilter.SetStepSize(m_step_size);
+  streamlineFilter.SetActiveField(m_field_name);
+  streamlineFilter.SetSeeds(seedsAH);
+  streamlineFilter.SetNumberOfSteps(m_num_steps);
+  auto out = streamlineFilter.Execute(inputs);
+
+  for (vtkm::Id i = 0; i < out.GetNumberOfPartitions(); i++)
+  {
+    this->m_output->AddDomain(out.GetPartition(i), i);
+  }
+#endif
+}
+
+} //  namespace vtkh