Lagrangian posthoc interpolation

src/tests/vtkh/CMakeLists.txt

@@ -94,6 +94,7 @@ set(MPI_TESTS t_vtk-h_smoke_par
               t_vtk-h_no_op_par
               t_vtk-h_histogram_par
               t_vtk-h_statistics_par
+              t_vtk-h_lagrangian_interpolation
               t_vtk-h_marching_cubes_par
               t_vtk-h_multi_render_par
               t_vtk-h_particle_advection_par

src/tests/vtkh/t_vtk-h_lagrangian_interpolation.cpp

@@ -0,0 +1,180 @@
+//-----------------------------------------------------------------------------
+///
+/// file: t_vtk-h_dataset.cpp
+///
+//-----------------------------------------------------------------------------
+
+#include "gtest/gtest.h"
+
+#include <vtkh/vtkh.hpp>
+#include <vtkh/DataSet.hpp>
+#include <vtkh/filters/Lagrangian.hpp>
+#include <vtkh/filters/LagrangianInterpolation.hpp>
+#include <vtkh/rendering/LineRenderer.hpp>
+#include <vtkh/rendering/Scene.hpp>
+#include "t_test_utils.hpp"
+#include <vtkm/cont/DataSet.h>
+#include <vtkm/cont/DataSetBuilderUniform.h>
+#include <vtkm/cont/DataSetFieldAdd.h>
+#include <vtkm/io/writer/VTKDataSetWriter.h> 
+#include <cstring>
+#include <sstream>
+#include <fstream>
+#include <string.h>
+#include <iostream>
+#include <mpi.h>
+
+bool remove_directory(const std::string &path)
+{
+  return ( remove(path.c_str()) == 0 ); 
+}
+
+bool create_directory(const std::string &path)
+{
+  return (mkdir(path.c_str(),S_IRWXU | S_IRWXG) == 0);
+}
+
+vtkm::cont::DataSet MakeTestUniformDataSet(vtkm::Id time)
+{
+  vtkm::Float64 xmin, xmax, ymin, ymax, zmin, zmax;
+  xmin = 0.0;
+  ymin = 0.0;
+  zmin = 0.0;
+
+  xmax = 10.0;
+  ymax = 10.0;
+  zmax = 10.0;
+
+  const vtkm::Id3 DIMS(16, 16, 16);
+
+  vtkm::cont::DataSetBuilderUniform dsb;
+
+  vtkm::Float64 xdiff = (xmax - xmin) / (static_cast<vtkm::Float64>(DIMS[0] - 1));
+  vtkm::Float64 ydiff = (ymax - ymin) / (static_cast<vtkm::Float64>(DIMS[1] - 1));
+  vtkm::Float64 zdiff = (zmax - zmin) / (static_cast<vtkm::Float64>(DIMS[2] - 1));
+
+  vtkm::Vec<vtkm::Float64, 3> ORIGIN(0, 0, 0);
+  vtkm::Vec<vtkm::Float64, 3> SPACING(xdiff, ydiff, zdiff);
+
+  vtkm::cont::DataSet dataset = dsb.Create(DIMS, ORIGIN, SPACING);
+  vtkm::cont::DataSetFieldAdd dsf;
+
+  vtkm::Id numPoints = DIMS[0] * DIMS[1] * DIMS[2];
+
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 3>> velocityField;
+  velocityField.Allocate(numPoints);
+
+  vtkm::Id count = 0;
+  for (vtkm::Id i = 0; i < DIMS[0]; i++)
+  {
+    for (vtkm::Id j = 0; j < DIMS[1]; j++)
+    {
+      for (vtkm::Id k = 0; k < DIMS[2]; k++)
+      {
+        velocityField.WritePortal().Set(count, vtkm::Vec<vtkm::Float64, 3>(0.01, 0.0, 0.0));
+        count++;
+      }
+    }
+  }
+  dsf.AddPointField(dataset, "velocity", velocityField);
+  return dataset;
+}
+
+void render_output(vtkh::DataSet *data, std::string file_name)
+{
+  data->AddConstantPointField(1.f,"lines");
+
+  vtkm::Bounds bounds = data->GetGlobalBounds();
+
+  vtkm::rendering::Camera camera;
+  camera.ResetToBounds(bounds);
+  float bg_color[4] = { 0.f, 0.f, 0.f, 1.f};
+  vtkh::Render render = vtkh::MakeRender(512,
+                                         512,
+                                         camera,
+                                         *data,
+                                         file_name,
+                                         bg_color);
+
+  vtkh::Scene scene;
+  scene.AddRender(render);
+
+  vtkh::LineRenderer tracer;
+  tracer.SetRadius(.1f);
+  tracer.SetInput(data);
+  tracer.SetField("lines");
+
+  scene.AddRenderer(&tracer);
+  scene.Render();
+}
+
+//----------------------------------------------------------------------------
+TEST(vtkh_lagrangian_interpolation, vtkh_serial_lagrangian_interpolation)
+{
+
+  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));
+
+  vtkh::Lagrangian lagrangian;
+  lagrangian.SetField("velocity");
+  lagrangian.SetStepSize(0.1);
+  lagrangian.SetWriteFrequency(10);
+  lagrangian.SetCustomSeedResolution(1);
+  lagrangian.SetSeedResolutionInX(1);
+  lagrangian.SetSeedResolutionInY(1);
+  lagrangian.SetSeedResolutionInZ(1);
+
+   std::cout << "Running Lagrangian filter test - vtkh" << std::endl;
+
+  for(vtkm::Id time = 1; time <= 10; ++time)
+  {
+    vtkh::DataSet data_set;
+    data_set.AddDomain(MakeTestUniformDataSet(time),0);
+    lagrangian.SetInput(&data_set);
+    lagrangian.Update();
+    vtkh::DataSet *extracted_basis = lagrangian.GetOutput();
+//    extracted_basis->PrintSummary(std::cout);
+    if(time == 10)  
+    {
+      create_directory("output");      
+      vtkm::cont::DataSet ds = extracted_basis->GetDomain(0);
+
+      std::stringstream file_path;
+      file_path << "output/Lagrangian_Structured_0_10.vtk";
+      vtkm::io::writer::VTKDataSetWriter writer(file_path.str());
+      writer.WriteDataSet(ds); 
+
+      std::ofstream seedfile;
+      seedfile.open("Seed.txt"); 
+      seedfile << "5 5 5\n";
+      seedfile.close();
+      vtkh::LagrangianInterpolation interp;
+      interp.SetField("displacement");
+      interp.SetInput(extracted_basis);
+      interp.SetRadius(0.1);
+      interp.SetNumSeeds(1);
+      interp.SetInterval(10);
+      interp.SetStartCycle(10);
+      interp.SetEndCycle(10);
+      interp.SetSeedPath("Seed.txt");
+      interp.SetBasisPath("output/");
+      interp.SetOutputPath("output/Pathlines_");  
+      interp.Update();
+    }
+  }
+
+
+  remove_directory("output/Lagrangian_Structured_0_10.vtk");
+  remove_directory("output/Pathlines_0_10.vtk");
+  remove_directory("output");
+  remove_directory("Seed.txt");
+
+  MPI_Finalize();
+
+}
+

src/vtkh/filters/CMakeLists.txt

@@ -18,6 +18,7 @@ set(vtkh_filters_headers
   IsoVolume.hpp
   NoOp.hpp
   Lagrangian.hpp
+  LagrangianInterpolation.hpp
   MarchingCubes.hpp
   Particle.hpp
   ParticleAdvection.hpp
@@ -56,6 +57,7 @@ set(vtkh_filters_sources
   IsoVolume.cpp
   NoOp.cpp
   Lagrangian.cpp
+  LagrangianInterpolation.cpp
   MarchingCubes.cpp
   ParticleAdvection.cpp
   PointAverage.cpp