diff --git a/vtkm/filter/Lagrangian.h b/vtkm/filter/Lagrangian.h index eebfaea40..06d1819bf 100644 --- a/vtkm/filter/Lagrangian.h +++ b/vtkm/filter/Lagrangian.h @@ -63,10 +63,13 @@ public: void UpdateSeedResolution(vtkm::cont::DataSet input); VTKM_CONT - void WriteDataSet(vtkm::Id cycle, const std::string& filename, vtkm::cont::DataSet dataset); + void InitializeSeedPositions(const vtkm::cont::DataSet& input); VTKM_CONT - void InitializeUniformSeeds(const vtkm::cont::DataSet& input); + void InitializeCoordinates(const vtkm::cont::DataSet& input, + std::vector& xC, + std::vector& yC, + std::vector& zC); template VTKM_CONT vtkm::cont::DataSet DoExecute( diff --git a/vtkm/filter/Lagrangian.hxx b/vtkm/filter/Lagrangian.hxx index fde238a93..d32d29f58 100644 --- a/vtkm/filter/Lagrangian.hxx +++ b/vtkm/filter/Lagrangian.hxx @@ -18,6 +18,8 @@ #include #include #include +#include +#include #include #include #include @@ -56,9 +58,7 @@ public: vtkm::Id steps = end_point.NumSteps; if (steps > 0 && res == 1) { - if (end_point.Pos[0] >= bounds.X.Min && end_point.Pos[0] <= bounds.X.Max && - end_point.Pos[1] >= bounds.Y.Min && end_point.Pos[1] <= bounds.Y.Max && - end_point.Pos[2] >= bounds.Z.Min && end_point.Pos[2] <= bounds.Z.Max) + if (bounds.Contains(end_point.Pos)) { res = 1; } @@ -76,6 +76,24 @@ public: private: vtkm::Bounds bounds; }; + +class DisplacementCalculation : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldIn end_point, FieldIn start_point, FieldInOut output); + using ExecutionSignature = void(_1, _2, _3); + using InputDomain = _1; + + template + VTKM_EXEC void operator()(const vtkm::Particle& end_point, + const vtkm::Particle& start_point, + DisplacementType& res) const + { + res[0] = end_point.Pos[0] - start_point.Pos[0]; + res[1] = end_point.Pos[1] - start_point.Pos[1]; + res[2] = end_point.Pos[2] - start_point.Pos[2]; + } +}; } namespace vtkm @@ -100,17 +118,6 @@ inline VTKM_CONT Lagrangian::Lagrangian() { } -//----------------------------------------------------------------------------- -inline void Lagrangian::WriteDataSet(vtkm::Id cycle, - const std::string& filename, - vtkm::cont::DataSet dataset) -{ - std::stringstream file_stream; - file_stream << filename << cycle << ".vtk"; - vtkm::io::writer::VTKDataSetWriter writer(file_stream.str()); - writer.WriteDataSet(dataset); -} - //----------------------------------------------------------------------------- inline void Lagrangian::UpdateSeedResolution(const vtkm::cont::DataSet input) { @@ -156,7 +163,7 @@ inline void Lagrangian::UpdateSeedResolution(const vtkm::cont::DataSet input) //----------------------------------------------------------------------------- -inline void Lagrangian::InitializeUniformSeeds(const vtkm::cont::DataSet& input) +inline void Lagrangian::InitializeSeedPositions(const vtkm::cont::DataSet& input) { vtkm::Bounds bounds = input.GetCoordinateSystem().GetBounds(); @@ -200,6 +207,39 @@ inline void Lagrangian::InitializeUniformSeeds(const vtkm::cont::DataSet& input) } } +//----------------------------------------------------------------------------- +inline void Lagrangian::InitializeCoordinates(const vtkm::cont::DataSet& input, + std::vector& xC, + std::vector& yC, + std::vector& zC) +{ + vtkm::Bounds bounds = input.GetCoordinateSystem().GetBounds(); + + vtkm::Float64 x_spacing = 0.0, y_spacing = 0.0, z_spacing = 0.0; + if (this->SeedRes[0] > 1) + x_spacing = (double)(bounds.X.Max - bounds.X.Min) / (double)(this->SeedRes[0] - 1); + if (this->SeedRes[1] > 1) + y_spacing = (double)(bounds.Y.Max - bounds.Y.Min) / (double)(this->SeedRes[1] - 1); + if (this->SeedRes[2] > 1) + z_spacing = (double)(bounds.Z.Max - bounds.Z.Min) / (double)(this->SeedRes[2] - 1); + // Divide by zero handling for 2D data set. How is this handled + + for (int x = 0; x < this->SeedRes[0]; x++) + { + vtkm::FloatDefault xi = static_cast(x * x_spacing); + xC.push_back(bounds.X.Min + xi); + } + for (int y = 0; y < this->SeedRes[1]; y++) + { + vtkm::FloatDefault yi = static_cast(y * y_spacing); + yC.push_back(bounds.Y.Min + yi); + } + for (int z = 0; z < this->SeedRes[2]; z++) + { + vtkm::FloatDefault zi = static_cast(z * z_spacing); + zC.push_back(bounds.Z.Min + zi); + } +} //----------------------------------------------------------------------------- template @@ -212,7 +252,7 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute( if (cycle == 0) { - InitializeUniformSeeds(input); + InitializeSeedPositions(input); BasisParticlesOriginal.Allocate(this->SeedRes[0] * this->SeedRes[1] * this->SeedRes[2]); vtkm::cont::ArrayCopy(BasisParticles, BasisParticlesOriginal); } @@ -246,66 +286,28 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute( RK4Type rk4(gridEval, static_cast(this->stepSize)); res = particleadvection.Run(rk4, basisParticleArray, 1); // Taking a single step - auto particles = res.Particles; - auto particlePortal = particles.GetPortalControl(); - - auto start_position = BasisParticlesOriginal.GetPortalControl(); - auto portal_validity = BasisParticlesValidity.GetPortalControl(); vtkm::cont::DataSet outputData; - vtkm::cont::DataSetBuilderExplicit dataSetBuilder; + vtkm::cont::DataSetBuilderRectilinear dataSetBuilder; if (cycle % this->writeFrequency == 0) { - int connectivity_index = 0; - std::vector connectivity; - std::vector> pointCoordinates; - std::vector shapes; - std::vector numIndices; + /* Steps to create a structured dataset */ + vtkm::cont::ArrayHandle BasisParticlesDisplacement; + BasisParticlesDisplacement.Allocate(this->SeedRes[0] * this->SeedRes[1] * this->SeedRes[2]); + DisplacementCalculation displacement; + this->Invoke(displacement, particles, BasisParticlesOriginal, BasisParticlesDisplacement); + std::vector xC, yC, zC; + InitializeCoordinates(input, xC, yC, zC); + outputData = dataSetBuilder.Create(xC, yC, zC); + vtkm::cont::DataSetFieldAdd dataSetFieldAdd; + dataSetFieldAdd.AddPointField(outputData, "valid", BasisParticlesValidity); + dataSetFieldAdd.AddPointField(outputData, "displacement", BasisParticlesDisplacement); - for (vtkm::Id index = 0; index < particlePortal.GetNumberOfValues(); index++) - { - auto start_point = start_position.Get(index); - auto end_point = particlePortal.Get(index).Pos; - auto steps = particlePortal.Get(index).NumSteps; - - if (steps > 0 && portal_validity.Get(index) == 1) - { - if (bounds.Contains(end_point)) - { - connectivity.push_back(connectivity_index); - connectivity.push_back(connectivity_index + 1); - connectivity_index += 2; - pointCoordinates.push_back( - vtkm::Vec3f(static_cast(start_point.Pos[0]), - static_cast(start_point.Pos[1]), - static_cast(start_point.Pos[2]))); - pointCoordinates.push_back(vtkm::Vec3f(static_cast(end_point[0]), - static_cast(end_point[1]), - static_cast(end_point[2]))); - shapes.push_back(vtkm::CELL_SHAPE_LINE); - numIndices.push_back(2); - } - else - { - portal_validity.Set(index, 0); - } - } - else - { - portal_validity.Set(index, 0); - } - } - - outputData = dataSetBuilder.Create(pointCoordinates, shapes, numIndices, connectivity); - std::stringstream file_path; - file_path << "output/basisflows_" << this->rank << "_"; - auto f_path = file_path.str(); - WriteDataSet(cycle, f_path, outputData); if (this->resetParticles) { - InitializeUniformSeeds(input); + InitializeSeedPositions(input); BasisParticlesOriginal.Allocate(this->SeedRes[0] * this->SeedRes[1] * this->SeedRes[2]); vtkm::cont::ArrayCopy(BasisParticles, BasisParticlesOriginal); } diff --git a/vtkm/filter/testing/UnitTestLagrangianFilter.cxx b/vtkm/filter/testing/UnitTestLagrangianFilter.cxx index 15be1315e..79657e883 100644 --- a/vtkm/filter/testing/UnitTestLagrangianFilter.cxx +++ b/vtkm/filter/testing/UnitTestLagrangianFilter.cxx @@ -80,15 +80,17 @@ void TestLagrangianFilterMultiStepInterval() { VTKM_TEST_ASSERT(extractedBasisFlows.GetNumberOfCoordinateSystems() == 1, "Wrong number of coordinate systems in the output dataset."); - VTKM_TEST_ASSERT(extractedBasisFlows.GetNumberOfCells() == 3375, + VTKM_TEST_ASSERT(extractedBasisFlows.GetNumberOfPoints() == 4096, "Wrong number of basis flows extracted."); + VTKM_TEST_ASSERT(extractedBasisFlows.GetNumberOfFields() == 2, "Wrong number of fields."); } else { - VTKM_TEST_ASSERT(extractedBasisFlows.GetNumberOfCells() == 0, - "Output dataset should have no cells."); + VTKM_TEST_ASSERT(extractedBasisFlows.GetNumberOfPoints() == 0, + "Output dataset should have no points."); VTKM_TEST_ASSERT(extractedBasisFlows.GetNumberOfCoordinateSystems() == 0, "Wrong number of coordinate systems in the output dataset."); + VTKM_TEST_ASSERT(extractedBasisFlows.GetNumberOfFields() == 0, "Wrong number of fields."); } } }