Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'change_lagrangian_output_structuredgrid'

Browse Source 
e68b283e8 Delete commented, unused class variables
4d358e5b2 Remove coordinates from class member variable list
e6647efaf Fix condition to populate coordinate variables
2e2a847f5 Change Lagrangian output format to structured grid

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Kenneth Moreland <kmorel@sandia.gov>
Merge-request: !1945
This commit is contained in:
Matt Larsen 2020-02-12 21:20:07 +00:00 committed by Kitware Robot
commit 2123a56ab3
3 changed files with 79 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
vtkm/filter/Lagrangian.h
View File

@ -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<Float64>& xC,
std::vector<Float64>& yC,
std::vector<Float64>& zC);
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(

136
vtkm/filter/Lagrangian.hxx
View File

@ -18,6 +18,8 @@
#include <vtkm/cont/ArrayPortalToIterators.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/DataSetBuilderRectilinear.h>
#include <vtkm/cont/DataSetFieldAdd.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/io/writer/VTKDataSetWriter.h>
@ -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 <typename DisplacementType>
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<Float64>& xC,
std::vector<Float64>& yC,
std::vector<Float64>& 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<vtkm::FloatDefault>(x * x_spacing);
xC.push_back(bounds.X.Min + xi);
}
for (int y = 0; y < this->SeedRes[1]; y++)
{
vtkm::FloatDefault yi = static_cast<vtkm::FloatDefault>(y * y_spacing);
yC.push_back(bounds.Y.Min + yi);
}
for (int z = 0; z < this->SeedRes[2]; z++)
{
vtkm::FloatDefault zi = static_cast<vtkm::FloatDefault>(z * z_spacing);
zC.push_back(bounds.Z.Min + zi);
}
}
//-----------------------------------------------------------------------------
template <typename T, typename StorageType, typename DerivedPolicy>
@ -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<vtkm::Float32>(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<vtkm::Id> connectivity;
std::vector<vtkm::Vec<T, 3>> pointCoordinates;
std::vector<vtkm::UInt8> shapes;
std::vector<vtkm::IdComponent> numIndices;
/* Steps to create a structured dataset */
vtkm::cont::ArrayHandle<vtkm::Vec3f> BasisParticlesDisplacement;
BasisParticlesDisplacement.Allocate(this->SeedRes[0] * this->SeedRes[1] * this->SeedRes[2]);
DisplacementCalculation displacement;
this->Invoke(displacement, particles, BasisParticlesOriginal, BasisParticlesDisplacement);
std::vector<Float64> 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<vtkm::FloatDefault>(start_point.Pos[0]),
static_cast<vtkm::FloatDefault>(start_point.Pos[1]),
static_cast<vtkm::FloatDefault>(start_point.Pos[2])));
pointCoordinates.push_back(vtkm::Vec3f(static_cast<vtkm::FloatDefault>(end_point[0]),
static_cast<vtkm::FloatDefault>(end_point[1]),
static_cast<vtkm::FloatDefault>(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);
}

8
vtkm/filter/testing/UnitTestLagrangianFilter.cxx
View File

@ -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.");
}
}
}