From 755137a82200a35220b61cf0f7125c7c6cde31ee Mon Sep 17 00:00:00 2001 From: NAThompson Date: Sat, 13 Jun 2020 10:58:45 -0400 Subject: [PATCH] Particle advection tests with file. --- data/data/rectilinear/fishtank.vtk | 3 + .../rectilinear/fishtank_double_ascii.vtk | 3 + .../fishtank_double_big_endian.vtk | 3 + data/data/rectilinear/fusion.vtk | 3 + .../contour_tree_augmented/ContourTreeApp.cxx | 7 +- .../testing/UnitTestStreamlineFilter.cxx | 91 +++++++ vtkm/io/VTKRectilinearGridReader.cxx | 67 ++++- vtkm/io/testing/UnitTestVTKDataSetReader.cxx | 256 ++++++++++++++++++ vtkm/worklet/testing/CMakeLists.txt | 2 +- .../testing/UnitTestParticleAdvection.cxx | 133 +++++++++ 10 files changed, 549 insertions(+), 19 deletions(-) create mode 100644 data/data/rectilinear/fishtank.vtk create mode 100644 data/data/rectilinear/fishtank_double_ascii.vtk create mode 100644 data/data/rectilinear/fishtank_double_big_endian.vtk create mode 100644 data/data/rectilinear/fusion.vtk diff --git a/data/data/rectilinear/fishtank.vtk b/data/data/rectilinear/fishtank.vtk new file mode 100644 index 000000000..e1ab1c481 --- /dev/null +++ b/data/data/rectilinear/fishtank.vtk @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ef3dfd79f0c8d18780d0749014d71c0226134041283d33de0bcd994e343dd421 +size 2001070 diff --git a/data/data/rectilinear/fishtank_double_ascii.vtk b/data/data/rectilinear/fishtank_double_ascii.vtk new file mode 100644 index 000000000..fce0d5b9c --- /dev/null +++ b/data/data/rectilinear/fishtank_double_ascii.vtk @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2bb3d36ea5ecef5e7ef1057d0dddebbc590424915083091ead3dac2928000524 +size 2904465 diff --git a/data/data/rectilinear/fishtank_double_big_endian.vtk b/data/data/rectilinear/fishtank_double_big_endian.vtk new file mode 100644 index 000000000..7ecfbc133 --- /dev/null +++ b/data/data/rectilinear/fishtank_double_big_endian.vtk @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bffad7dae3dd6ef018ad7a9e109464ced0f3b9bc15cf1fb5d555f6d0d00b621f +size 3001624 diff --git a/data/data/rectilinear/fusion.vtk b/data/data/rectilinear/fusion.vtk new file mode 100644 index 000000000..5c5212da1 --- /dev/null +++ b/data/data/rectilinear/fusion.vtk @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2cbdf56fd5445ddc5b6bc05507b8825fb8d74fe1ccce894bde03e5ff2ecf5fb6 +size 525141 diff --git a/examples/contour_tree_augmented/ContourTreeApp.cxx b/examples/contour_tree_augmented/ContourTreeApp.cxx index 5674e719b..658700993 100644 --- a/examples/contour_tree_augmented/ContourTreeApp.cxx +++ b/examples/contour_tree_augmented/ContourTreeApp.cxx @@ -64,7 +64,6 @@ #include #include #include -#include #include #include #include @@ -459,17 +458,17 @@ int main(int argc, char* argv[]) // TODO All we should need to do to implement BOV support is to copy the values // in the values vector and copy the dimensions in the dims vector vtkm::Id nRows, nCols, nSlices; - vtkm::worklet::contourtree_augmented::GetRowsColsSlices temp; + vtkm::filter::GetRowsColsSlices temp; temp(inDataSet.GetCellSet(), nRows, nCols, nSlices); dims[0] = nRows; dims[1] = nCols; dims[2] = nSlices; auto tempField = inDataSet.GetField("values").GetData(); - values.resize(static_cast(tempField.GetNumberOfValues())); + values.resize(tempField.GetNumberOfValues()); auto tempFieldHandle = tempField.AsVirtual().ReadPortal(); for (vtkm::Id i = 0; i < tempField.GetNumberOfValues(); i++) { - values[static_cast(i)] = static_cast(tempFieldHandle.Get(i)); + values[i] = static_cast(tempFieldHandle.Get(i)); } VTKM_LOG_S(vtkm::cont::LogLevel::Error, "BOV reader not yet support in MPI mode by this example"); diff --git a/vtkm/filter/testing/UnitTestStreamlineFilter.cxx b/vtkm/filter/testing/UnitTestStreamlineFilter.cxx index 8422e55ba..dd6340fe8 100644 --- a/vtkm/filter/testing/UnitTestStreamlineFilter.cxx +++ b/vtkm/filter/testing/UnitTestStreamlineFilter.cxx @@ -12,6 +12,8 @@ #include #include #include +#include +#include namespace { @@ -102,10 +104,99 @@ void TestPathline() VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == 3, "Wrong number of cells"); } +void TestStreamlineFile(const std::string& fname, + const std::vector& pts, + vtkm::FloatDefault stepSize, + vtkm::Id maxSteps, + const std::vector& endPts) +{ + vtkm::io::VTKDataSetReader reader(fname); + vtkm::cont::DataSet ds; + try + { + ds = reader.ReadDataSet(); + } + catch (vtkm::io::ErrorIO& e) + { + std::string message("Error reading: "); + message += fname; + message += ", "; + message += e.GetMessage(); + + VTKM_TEST_FAIL(message.c_str()); + } + vtkm::Id numPoints = static_cast(pts.size()); + + std::vector seeds; + for (vtkm::Id i = 0; i < numPoints; i++) + seeds.push_back(vtkm::Particle(pts[static_cast(i)], i)); + auto seedArray = vtkm::cont::make_ArrayHandle(seeds); + + vtkm::filter::Streamline streamline; + streamline.SetStepSize(stepSize); + streamline.SetNumberOfSteps(maxSteps); + streamline.SetSeeds(seedArray); + + streamline.SetActiveField("vec"); + auto output = streamline.Execute(ds); + + auto coords = output.GetCoordinateSystem().GetData(); + vtkm::cont::DynamicCellSet dcells = output.GetCellSet(); + VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == numPoints, "Wrong number of cells"); + VTKM_TEST_ASSERT(dcells.IsType>(), "Wrong cell type"); + + auto cells = dcells.Cast>(); + auto cPortal = coords.ReadPortal(); + const vtkm::FloatDefault eps = 1e-3; + + for (vtkm::Id i = 0; i < numPoints; i++) + { + vtkm::Id numPts = cells.GetNumberOfPointsInCell(i); + std::vector ids(static_cast(numPts)); + cells.GetCellPointIds(i, ids.data()); + + vtkm::Vec3f e = endPts[static_cast(i)]; + vtkm::Vec3f pt = cPortal.Get(ids[ids.size() - 1]); + if (vtkm::Magnitude(pt - e) > eps) + { + VTKM_LOG_S(vtkm::cont::LogLevel::Error, + "Expected magnitude of error <=" << eps << ", but got error " + << vtkm::Magnitude(pt - e)); + VTKM_TEST_ASSERT(false, "Streamline end point is wrong."); + } + } +} + void TestStreamlineFilters() { TestStreamline(); TestPathline(); + + std::string basePath = vtkm::cont::testing::Testing::GetTestDataBasePath(); + + //Fusion test. + std::vector fusionPts, fusionEndPts; + fusionPts.push_back(vtkm::Vec3f(0.8f, 0.6f, 0.6f)); + fusionPts.push_back(vtkm::Vec3f(0.8f, 0.8f, 0.6f)); + fusionPts.push_back(vtkm::Vec3f(0.8f, 0.8f, 0.3f)); + fusionEndPts.push_back(vtkm::Vec3f(0.5335789918f, 0.87112802267f, 0.6723330020f)); + fusionEndPts.push_back(vtkm::Vec3f(0.5601879954f, 0.91389900446f, 0.43989110522f)); + fusionEndPts.push_back(vtkm::Vec3f(0.7004770041f, 0.63193398714f, 0.64524400234f)); + vtkm::FloatDefault fusionStep = 0.005f; + std::string fusionFile = basePath + "/rectilinear/fusion.vtk"; + TestStreamlineFile(fusionFile, fusionPts, fusionStep, 1000, fusionEndPts); + + //Fishtank test. + std::vector fishPts, fishEndPts; + fishPts.push_back(vtkm::Vec3f(0.75f, 0.5f, 0.01f)); + fishPts.push_back(vtkm::Vec3f(0.4f, 0.2f, 0.7f)); + fishPts.push_back(vtkm::Vec3f(0.5f, 0.3f, 0.8f)); + fishEndPts.push_back(vtkm::Vec3f(0.7734669447f, 0.4870159328f, 0.8979591727f)); + fishEndPts.push_back(vtkm::Vec3f(0.7257543206f, 0.1277695596f, 0.7468645573f)); + fishEndPts.push_back(vtkm::Vec3f(0.8347796798f, 0.1276152730f, 0.4985143244f)); + vtkm::FloatDefault fishStep = 0.001f; + std::string fishFile = basePath + "/rectilinear/fishtank.vtk"; + TestStreamlineFile(fishFile, fishPts, fishStep, 100, fishEndPts); } } diff --git a/vtkm/io/VTKRectilinearGridReader.cxx b/vtkm/io/VTKRectilinearGridReader.cxx index 9dcf3a5ca..81c038954 100644 --- a/vtkm/io/VTKRectilinearGridReader.cxx +++ b/vtkm/io/VTKRectilinearGridReader.cxx @@ -46,27 +46,34 @@ void VTKRectilinearGridReader::Read() this->DataFile->Stream >> dim[0] >> dim[1] >> dim[2] >> std::ws; //Read the points. - std::string dataType; + std::string fileStorageDataType; std::size_t numPoints[3]; vtkm::cont::VariantArrayHandle X, Y, Z; - // Always read coordinates as vtkm::FloatDefault - std::string readDataType = vtkm::io::internal::DataTypeName::Name(); - - this->DataFile->Stream >> tag >> numPoints[0] >> dataType >> std::ws; + this->DataFile->Stream >> tag >> numPoints[0] >> fileStorageDataType >> std::ws; if (tag != "X_COORDINATES") throw vtkm::io::ErrorIO("X_COORDINATES tag not found"); - X = this->DoReadArrayVariant(vtkm::cont::Field::Association::ANY, readDataType, numPoints[0], 1); - this->DataFile->Stream >> tag >> numPoints[1] >> dataType >> std::ws; + // In binary mode, we must read the data as they are stored in the file. + // In text mode we can parse as FloatDefault no matter the precision of the storage. + X = this->DoReadArrayVariant( + vtkm::cont::Field::Association::ANY, fileStorageDataType, numPoints[0], 1); + + + this->DataFile->Stream >> tag >> numPoints[1] >> fileStorageDataType >> std::ws; if (tag != "Y_COORDINATES") throw vtkm::io::ErrorIO("Y_COORDINATES tag not found"); - Y = this->DoReadArrayVariant(vtkm::cont::Field::Association::ANY, readDataType, numPoints[1], 1); - this->DataFile->Stream >> tag >> numPoints[2] >> dataType >> std::ws; + Y = this->DoReadArrayVariant( + vtkm::cont::Field::Association::ANY, fileStorageDataType, numPoints[1], 1); + + this->DataFile->Stream >> tag >> numPoints[2] >> fileStorageDataType >> std::ws; if (tag != "Z_COORDINATES") throw vtkm::io::ErrorIO("Z_COORDINATES tag not found"); - Z = this->DoReadArrayVariant(vtkm::cont::Field::Association::ANY, readDataType, numPoints[2], 1); + + Z = this->DoReadArrayVariant( + vtkm::cont::Field::Association::ANY, fileStorageDataType, numPoints[2], 1); + if (dim != vtkm::Id3(static_cast(numPoints[0]), static_cast(numPoints[1]), @@ -78,14 +85,46 @@ void VTKRectilinearGridReader::Read() vtkm::cont::ArrayHandle> coords; + // We need to store all coordinate arrays as FloatDefault. vtkm::cont::ArrayHandle Xc, Yc, Zc; - X.CopyTo(Xc); - Y.CopyTo(Yc); - Z.CopyTo(Zc); + // But the VariantArrayHandle has type fileStorageDataType. + // If the fileStorageDataType is the same as the storage type, no problem: + if (fileStorageDataType == vtkm::io::internal::DataTypeName::Name()) + { + X.CopyTo(Xc); + Y.CopyTo(Yc); + Z.CopyTo(Zc); + } + else + { + // Two cases if the data in the file differs from FloatDefault: + if (fileStorageDataType == "float") + { + vtkm::cont::ArrayHandle Xcf, Ycf, Zcf; + X.CopyTo(Xcf); + Y.CopyTo(Ycf); + Z.CopyTo(Zcf); + vtkm::cont::ArrayCopy(Xcf, Xc); + vtkm::cont::ArrayCopy(Ycf, Yc); + vtkm::cont::ArrayCopy(Zcf, Zc); + } + else + { + vtkm::cont::ArrayHandle Xcd, Ycd, Zcd; + X.CopyTo(Xcd); + Y.CopyTo(Ycd); + Z.CopyTo(Zcd); + vtkm::cont::ArrayCopy(Xcd, Xc); + vtkm::cont::ArrayCopy(Ycd, Yc); + vtkm::cont::ArrayCopy(Zcd, Zc); + } + } + // As a postscript to this somewhat branchy code, I thought that X.CopyTo(Xc) should be able to cast to the value_type of Xc. + // But that would break the ability to make the cast lazy, and hence if you change it you induce performance bugs. + coords = vtkm::cont::make_ArrayHandleCartesianProduct(Xc, Yc, Zc); vtkm::cont::CoordinateSystem coordSys("coordinates", coords); this->DataSet.AddCoordinateSystem(coordSys); - this->DataSet.SetCellSet(internal::CreateCellSetStructured(dim)); // Read points and cell attributes diff --git a/vtkm/io/testing/UnitTestVTKDataSetReader.cxx b/vtkm/io/testing/UnitTestVTKDataSetReader.cxx index 185503a4b..06747a884 100644 --- a/vtkm/io/testing/UnitTestVTKDataSetReader.cxx +++ b/vtkm/io/testing/UnitTestVTKDataSetReader.cxx @@ -578,6 +578,254 @@ void TestReadingStructuredGridBin() "Incorrect cellset type"); } +void TestReadingFishTank() +{ + std::string fishtank = + vtkm::cont::testing::Testing::GetTestDataBasePath() + "/rectilinear/fishtank.vtk"; + vtkm::cont::DataSet ds = readVTKDataSet(fishtank.c_str()); + + // This is information you can glean by running 'strings' on fishtank.vtk: + VTKM_TEST_ASSERT(ds.GetCellSet().IsType>(), + "Incorrect cellset type"); + VTKM_TEST_ASSERT(ds.GetNumberOfPoints() == 50 * 50 * 50, "Incorrect number of points"); + VTKM_TEST_ASSERT(ds.GetCellSet().GetNumberOfPoints() == 50 * 50 * 50, + "Incorrect number of points (from cell set)"); + VTKM_TEST_ASSERT(ds.GetNumberOfFields() == 2, "Incorrect number of fields"); + VTKM_TEST_ASSERT(ds.HasField("vec"), "The vtk file has a field 'vec', but the dataset does not."); + VTKM_TEST_ASSERT(ds.HasField("vec_magnitude"), + "The vtk file has a field 'vec_magnitude', but the dataset does not."); + + // I believe the coordinate system is implicitly given by the first element of X_COORDINATES: + VTKM_TEST_ASSERT(ds.GetNumberOfCoordinateSystems() == 1, + "Need one and only one coordinate system."); + // In order to get the data from the coordinate system, I used the following workflow: + // First, I deleted all ascii header lines just past 'X_COORDINATES 50 float'. + // Once this is done, I can get the binary data from + // $ od -tfF --endian=big fishtank_copy.vtk + // The result is: + // 0 0.020408163 ... 0.9591837 0.97959185 1 + // So monotone increasing, bound [0,1]. + const vtkm::cont::CoordinateSystem& coordinateSystem = ds.GetCoordinateSystem(); + vtkm::Vec ranges = coordinateSystem.GetRange(); + vtkm::Range xRange = ranges[0]; + VTKM_TEST_ASSERT(xRange.Min == 0); + VTKM_TEST_ASSERT(xRange.Max == 1); + // Do the same past 'Y_COORDINATES 50 float'. + // You get exactly the same as the x data. + vtkm::Range yRange = ranges[1]; + VTKM_TEST_ASSERT(yRange.Min == 0); + VTKM_TEST_ASSERT(yRange.Max == 1); + // And finally, do it past 'Z_COORDINATES 50 float': + vtkm::Range zRange = ranges[2]; + VTKM_TEST_ASSERT(zRange.Min == 0); + VTKM_TEST_ASSERT(zRange.Max == 1); + + // Now delete the text up to LOOKUP TABLE default. + // I see: + // 0 0 0 0 3.5267966 . . . + // This is a vector magnitude, so all values must be >= 0. + // A cursory glance shows that 124.95 is a large value, so we can sanity check the data with the bounds + // [0, ~130]. + // And if we open the file in Paraview, we can observe the bounds [0, 156.905]. + const vtkm::cont::Field& vec_magnitude = ds.GetField("vec_magnitude"); + VTKM_TEST_ASSERT(vec_magnitude.GetName() == "vec_magnitude"); + VTKM_TEST_ASSERT(vec_magnitude.IsFieldPoint()); + + vtkm::Range mag_range; + vec_magnitude.GetRange(&mag_range); + VTKM_TEST_ASSERT(mag_range.Min == 0); + VTKM_TEST_ASSERT(mag_range.Max <= 156.906); + + // This info was gleaned from the Paraview Information panel: + const vtkm::cont::Field& vec = ds.GetField("vec"); + VTKM_TEST_ASSERT(vec.GetName() == "vec"); + VTKM_TEST_ASSERT(vec.IsFieldPoint()); + // Bounds from Information panel: + // [-65.3147, 86.267], [-88.0325, 78.7217], [-67.0969, 156.867] + const vtkm::cont::ArrayHandle& vecRanges = vec.GetRange(); + VTKM_TEST_ASSERT(vecRanges.GetNumberOfValues() == 3); + auto vecRangesReadPortal = vecRanges.ReadPortal(); + + auto xVecRange = vecRangesReadPortal.Get(0); + VTKM_TEST_ASSERT(xVecRange.Min >= -65.3148 && xVecRange.Min <= -65.3146); + VTKM_TEST_ASSERT(xVecRange.Max >= 86.26 && xVecRange.Min <= 86.268); + + auto yVecRange = vecRangesReadPortal.Get(1); + VTKM_TEST_ASSERT(yVecRange.Min >= -88.0326 && yVecRange.Min <= -88.0324); + VTKM_TEST_ASSERT(yVecRange.Max >= 78.721); + VTKM_TEST_ASSERT(yVecRange.Max <= 78.7218); + + auto zVecRange = vecRangesReadPortal.Get(2); + VTKM_TEST_ASSERT(zVecRange.Min >= -67.097 && zVecRange.Min <= -67.096); + VTKM_TEST_ASSERT(zVecRange.Max >= 156.866 && zVecRange.Max <= 156.868); +} + +void TestReadingDoublePrecisionFishTank() +{ + std::string fishtank = vtkm::cont::testing::Testing::GetTestDataBasePath() + + "/rectilinear/fishtank_double_big_endian.vtk"; + vtkm::cont::DataSet ds = readVTKDataSet(fishtank.c_str()); + + // This is information you can glean by running 'strings' on fishtank.vtk: + VTKM_TEST_ASSERT(ds.GetCellSet().IsType>(), + "Incorrect cellset type"); + VTKM_TEST_ASSERT(ds.GetNumberOfPoints() == 50 * 50 * 50, "Incorrect number of points"); + VTKM_TEST_ASSERT(ds.GetCellSet().GetNumberOfPoints() == 50 * 50 * 50, + "Incorrect number of points (from cell set)"); + + VTKM_TEST_ASSERT(ds.HasField("vec"), "The vtk file has a field 'vec', but the dataset does not."); + + + VTKM_TEST_ASSERT(ds.GetNumberOfCoordinateSystems() == 1, + "fishtank has one and only one coordinate system."); + // See the single precision version for info: + const vtkm::cont::CoordinateSystem& coordinateSystem = ds.GetCoordinateSystem(); + vtkm::Vec ranges = coordinateSystem.GetRange(); + vtkm::Range xRange = ranges[0]; + VTKM_TEST_ASSERT(xRange.Min == 0); + VTKM_TEST_ASSERT(xRange.Max == 1); + vtkm::Range yRange = ranges[1]; + VTKM_TEST_ASSERT(yRange.Min == 0); + VTKM_TEST_ASSERT(yRange.Max == 1); + vtkm::Range zRange = ranges[2]; + VTKM_TEST_ASSERT(zRange.Min == 0); + VTKM_TEST_ASSERT(zRange.Max == 1); + + // This info was gleaned from the Paraview Information panel: + const vtkm::cont::Field& vec = ds.GetField("vec"); + VTKM_TEST_ASSERT(vec.GetName() == "vec"); + VTKM_TEST_ASSERT(vec.IsFieldPoint()); + // Bounds from Information panel: + // [-65.3147, 86.267], [-88.0325, 78.7217], [-67.0969, 156.867] + const vtkm::cont::ArrayHandle& vecRanges = vec.GetRange(); + VTKM_TEST_ASSERT(vecRanges.GetNumberOfValues() == 3); + auto vecRangesReadPortal = vecRanges.ReadPortal(); + + auto xVecRange = vecRangesReadPortal.Get(0); + VTKM_TEST_ASSERT(xVecRange.Min >= -65.3148 && xVecRange.Min <= -65.3146); + VTKM_TEST_ASSERT(xVecRange.Max >= 86.26 && xVecRange.Min <= 86.268); + + auto yVecRange = vecRangesReadPortal.Get(1); + VTKM_TEST_ASSERT(yVecRange.Min >= -88.0326 && yVecRange.Min <= -88.0324); + VTKM_TEST_ASSERT(yVecRange.Max >= 78.721); + VTKM_TEST_ASSERT(yVecRange.Max <= 78.7218); + + auto zVecRange = vecRangesReadPortal.Get(2); + VTKM_TEST_ASSERT(zVecRange.Min >= -67.097 && zVecRange.Min <= -67.096); + VTKM_TEST_ASSERT(zVecRange.Max >= 156.866 && zVecRange.Max <= 156.868); +} + +void TestReadingASCIIFishTank() +{ + std::string fishtank = + vtkm::cont::testing::Testing::GetTestDataBasePath() + "/rectilinear/fishtank_double_ascii.vtk"; + vtkm::cont::DataSet ds = readVTKDataSet(fishtank.c_str()); + VTKM_TEST_ASSERT(ds.GetCellSet().IsType>(), + "Incorrect cellset type"); + VTKM_TEST_ASSERT(ds.GetNumberOfPoints() == 50 * 50 * 50, "Incorrect number of points"); + VTKM_TEST_ASSERT(ds.GetCellSet().GetNumberOfPoints() == 50 * 50 * 50, + "Incorrect number of points (from cell set)"); + VTKM_TEST_ASSERT(ds.HasField("vec"), "The vtk file has a field 'vec', but the dataset does not."); + VTKM_TEST_ASSERT(ds.GetNumberOfCoordinateSystems() == 1, + "fishtank has one and only one coordinate system."); + const vtkm::cont::CoordinateSystem& coordinateSystem = ds.GetCoordinateSystem(); + vtkm::Vec ranges = coordinateSystem.GetRange(); + vtkm::Range xRange = ranges[0]; + VTKM_TEST_ASSERT(xRange.Min == 0); + VTKM_TEST_ASSERT(xRange.Max == 1); + vtkm::Range yRange = ranges[1]; + VTKM_TEST_ASSERT(yRange.Min == 0); + VTKM_TEST_ASSERT(yRange.Max == 1); + vtkm::Range zRange = ranges[2]; + VTKM_TEST_ASSERT(zRange.Min == 0); + VTKM_TEST_ASSERT(zRange.Max == 1); + + const vtkm::cont::Field& vec = ds.GetField("vec"); + VTKM_TEST_ASSERT(vec.GetName() == "vec"); + VTKM_TEST_ASSERT(vec.IsFieldPoint()); + // Bounds from Paraview information panel: + // [-65.3147, 86.267], [-88.0325, 78.7217], [-67.0969, 156.867] + const vtkm::cont::ArrayHandle& vecRanges = vec.GetRange(); + VTKM_TEST_ASSERT(vecRanges.GetNumberOfValues() == 3); + auto vecRangesReadPortal = vecRanges.ReadPortal(); + auto xVecRange = vecRangesReadPortal.Get(0); + VTKM_TEST_ASSERT(xVecRange.Min >= -65.3148 && xVecRange.Min <= -65.3146); + VTKM_TEST_ASSERT(xVecRange.Max >= 86.26 && xVecRange.Min <= 86.268); + + auto yVecRange = vecRangesReadPortal.Get(1); + VTKM_TEST_ASSERT(yVecRange.Min >= -88.0326 && yVecRange.Min <= -88.0324); + VTKM_TEST_ASSERT(yVecRange.Max >= 78.721); + VTKM_TEST_ASSERT(yVecRange.Max <= 78.7218); + + auto zVecRange = vecRangesReadPortal.Get(2); + VTKM_TEST_ASSERT(zVecRange.Min >= -67.097 && zVecRange.Min <= -67.096); + VTKM_TEST_ASSERT(zVecRange.Max >= 156.866 && zVecRange.Max <= 156.868); +} + +void TestReadingFusion() +{ + std::string fusion = + vtkm::cont::testing::Testing::GetTestDataBasePath() + "/rectilinear/fusion.vtk"; + vtkm::cont::DataSet ds = readVTKDataSet(fusion.c_str()); + + VTKM_TEST_ASSERT(ds.GetCellSet().IsType>(), + "Incorrect cellset type"); + VTKM_TEST_ASSERT(ds.GetNumberOfPoints() == 32 * 32 * 32, "Incorrect number of points"); + VTKM_TEST_ASSERT(ds.GetCellSet().GetNumberOfPoints() == 32 * 32 * 32, + "Incorrect number of points (from cell set)"); + VTKM_TEST_ASSERT(ds.HasField("vec_magnitude"), + "The vtk file has a field 'vec_magnitude', but the dataset does not."); + VTKM_TEST_ASSERT(ds.HasField("vec"), "The vtk file has a field 'vec', but the dataset does not."); + VTKM_TEST_ASSERT(ds.GetNumberOfCoordinateSystems() == 1, + "The vtk file has a field 'vec', but the dataset does not."); + + // Taken from Paraview + clicking Data Axes Grid: + const vtkm::cont::CoordinateSystem& coordinateSystem = ds.GetCoordinateSystem(); + vtkm::Vec ranges = coordinateSystem.GetRange(); + vtkm::Range xRange = ranges[0]; + VTKM_TEST_ASSERT(xRange.Min == 0); + VTKM_TEST_ASSERT(xRange.Max == 1); + vtkm::Range yRange = ranges[1]; + VTKM_TEST_ASSERT(yRange.Min == 0); + VTKM_TEST_ASSERT(yRange.Max == 1); + vtkm::Range zRange = ranges[2]; + VTKM_TEST_ASSERT(zRange.Min == 0); + VTKM_TEST_ASSERT(zRange.Max == 1); + + // Paraview Information Panel of this file: + // vec_magnitude [0, 3.73778] + vtkm::cont::Field vec_magnitude = ds.GetField("vec_magnitude"); + VTKM_TEST_ASSERT(vec_magnitude.GetName() == "vec_magnitude"); + VTKM_TEST_ASSERT(vec_magnitude.IsFieldPoint()); + + vtkm::Range mag_range; + vec_magnitude.GetRange(&mag_range); + VTKM_TEST_ASSERT(mag_range.Min == 0); + VTKM_TEST_ASSERT(mag_range.Max <= 3.73779); + VTKM_TEST_ASSERT(mag_range.Max >= 3.73777); + + vtkm::cont::Field vec = ds.GetField("vec"); + VTKM_TEST_ASSERT(vec.GetName() == "vec"); + VTKM_TEST_ASSERT(vec.IsFieldPoint()); + const vtkm::cont::ArrayHandle& vecRanges = vec.GetRange(); + VTKM_TEST_ASSERT(vecRanges.GetNumberOfValues() == 3); + auto vecRangesReadPortal = vecRanges.ReadPortal(); + + // vec float [-3.41054, 3.40824], [-3.41018, 3.41036], [-0.689022, 0.480726] + auto xVecRange = vecRangesReadPortal.Get(0); + VTKM_TEST_ASSERT(test_equal(xVecRange.Min, -3.41054)); + VTKM_TEST_ASSERT(test_equal(xVecRange.Max, 3.40824)); + + auto yVecRange = vecRangesReadPortal.Get(1); + + VTKM_TEST_ASSERT(test_equal(yVecRange.Min, -3.41018)); + VTKM_TEST_ASSERT(test_equal(yVecRange.Max, 3.41036)); + + auto zVecRange = vecRangesReadPortal.Get(2); + VTKM_TEST_ASSERT(test_equal(zVecRange.Min, -0.689022)); + VTKM_TEST_ASSERT(test_equal(zVecRange.Max, 0.480726)); +} + void TestReadingVTKDataSet() { std::cout << "Test reading VTK Polydata file in ASCII" << std::endl; @@ -607,6 +855,14 @@ void TestReadingVTKDataSet() TestReadingStructuredGridASCII(); std::cout << "Test reading VTK StructuredGrid file in BINARY" << std::endl; TestReadingStructuredGridBin(); + std::cout << "Test reading float precision fishtank" << std::endl; + TestReadingFishTank(); + std::cout << "Test reading double precision fishtank" << std::endl; + TestReadingDoublePrecisionFishTank(); + std::cout << "Test ASCII fishtank" << std::endl; + TestReadingASCIIFishTank(); + std::cout << "Test reading fusion" << std::endl; + TestReadingFusion(); } int UnitTestVTKDataSetReader(int argc, char* argv[]) diff --git a/vtkm/worklet/testing/CMakeLists.txt b/vtkm/worklet/testing/CMakeLists.txt index b94137b72..f00f86a7e 100644 --- a/vtkm/worklet/testing/CMakeLists.txt +++ b/vtkm/worklet/testing/CMakeLists.txt @@ -89,7 +89,7 @@ set(unit_tests vtkm_unit_tests( SOURCES ${unit_tests} - LIBRARIES vtkm_source vtkm_worklet vtkm_filter + LIBRARIES vtkm_source vtkm_worklet vtkm_filter vtkm_io ALL_BACKENDS USE_VTKM_JOB_POOL ) diff --git a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx index e52a2e02e..381a25ec8 100644 --- a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx +++ b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -847,6 +848,111 @@ void TestWorkletsBasic() } } +template +void ValidateResult(const ResultType& res, + vtkm::Id maxSteps, + const std::vector& endPts) +{ + const vtkm::FloatDefault eps = 1e-3; + vtkm::Id numPts = static_cast(endPts.size()); + + VTKM_TEST_ASSERT(res.Particles.GetNumberOfValues() == numPts, + "Wrong number of points in particle advection result."); + + auto portal = res.Particles.ReadPortal(); + bool fail = false; + for (vtkm::Id i = 0; i < 3; i++) + { + vtkm::Vec3f p = portal.Get(i).Pos; + vtkm::Vec3f e = endPts[static_cast(i)]; + if (vtkm::Magnitude(p - e) > eps) + { + std::cout << std::setprecision(15) << "P_" << i << " p " << p << " e " << e << " diff " + << vtkm::Magnitude(p - e) << " eps= " << eps << std::endl; + fail = true; + } + + //VTKM_TEST_ASSERT(vtkm::Magnitude(p - e) <= eps, "Particle advection point is wrong"); + VTKM_TEST_ASSERT(portal.Get(i).NumSteps == maxSteps, "Particle advection NumSteps is wrong"); + VTKM_TEST_ASSERT(portal.Get(i).Status.CheckOk(), "Particle advection Status is wrong"); + VTKM_TEST_ASSERT(portal.Get(i).Status.CheckTerminate(), + "Particle advection particle did not terminate"); + } + VTKM_TEST_ASSERT(fail == false, "Particle advection point is wrong"); +} + + +void TestParticleAdvectionFile(const std::string& fname, + const std::vector& pts, + vtkm::FloatDefault stepSize, + vtkm::Id maxSteps, + const std::vector& endPts) +{ + VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Testing particle advection on file " << fname); + vtkm::io::VTKDataSetReader reader(fname); + vtkm::cont::DataSet ds; + try + { + ds = reader.ReadDataSet(); + } + catch (vtkm::io::ErrorIO& e) + { + std::string message("Error reading: "); + message += fname; + message += ", "; + message += e.GetMessage(); + + VTKM_TEST_FAIL(message.c_str()); + } + + using FieldHandle = vtkm::cont::ArrayHandle; + using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator; + using RK4Type = vtkm::worklet::particleadvection::RK4Integrator; + + VTKM_TEST_ASSERT(ds.HasField("vec")); + vtkm::cont::Field& field = ds.GetField("vec"); + auto fieldData = field.GetData(); + + if (!fieldData.IsType()) + { + VTKM_LOG_S(vtkm::cont::LogLevel::Error, + "The field data is of type " + << vtkm::cont::TypeToString() + << ", but we expect type vtkm::cont::ArrayHandle"); + VTKM_TEST_FAIL("No field with correct type found."); + } + + + FieldHandle fieldArray = fieldData.Cast(); + GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), fieldArray); + RK4Type rk4(eval, stepSize); + + for (int i = 0; i < 2; i++) + { + std::vector seeds; + for (size_t j = 0; j < pts.size(); j++) + seeds.push_back(vtkm::Particle(pts[j], static_cast(j))); + auto seedArray = vtkm::cont::make_ArrayHandle(seeds); + + if (i == 0) + { + vtkm::worklet::ParticleAdvection pa; + vtkm::worklet::ParticleAdvectionResult res; + + res = pa.Run(rk4, seedArray, maxSteps); + ValidateResult(res, maxSteps, endPts); + } + else if (i == 1) + { + vtkm::worklet::Streamline s; + vtkm::worklet::StreamlineResult res; + + res = s.Run(rk4, seedArray, maxSteps); + ValidateResult(res, maxSteps, endPts); + } + } +} + void TestParticleAdvection() { TestIntegrators(); @@ -854,6 +960,33 @@ void TestParticleAdvection() TestParticleStatus(); TestWorkletsBasic(); TestParticleWorkletsWithDataSetTypes(); + + std::string basePath = vtkm::cont::testing::Testing::GetTestDataBasePath(); + + //Fusion test. + std::vector fusionPts, fusionEndPts; + fusionPts.push_back(vtkm::Vec3f(0.8f, 0.6f, 0.6f)); + fusionPts.push_back(vtkm::Vec3f(0.8f, 0.8f, 0.6f)); + fusionPts.push_back(vtkm::Vec3f(0.8f, 0.8f, 0.3f)); + fusionEndPts.push_back(vtkm::Vec3f(0.5335789918f, 0.87112802267f, 0.6723330020f)); + fusionEndPts.push_back(vtkm::Vec3f(0.5601879954f, 0.91389900446f, 0.43989110522f)); + fusionEndPts.push_back(vtkm::Vec3f(0.7004770041f, 0.63193398714f, 0.64524400234f)); + vtkm::FloatDefault fusionStep = 0.005f; + std::string fusionFile = basePath + "/rectilinear/fusion.vtk"; + TestParticleAdvectionFile(fusionFile, fusionPts, fusionStep, 1000, fusionEndPts); + + //Fishtank test. + std::vector fishPts, fishEndPts; + fishPts.push_back(vtkm::Vec3f(0.75f, 0.5f, 0.01f)); + fishPts.push_back(vtkm::Vec3f(0.4f, 0.2f, 0.7f)); + fishPts.push_back(vtkm::Vec3f(0.5f, 0.3f, 0.8f)); + fishEndPts.push_back(vtkm::Vec3f(0.7734669447f, 0.4870159328f, 0.8979591727f)); + fishEndPts.push_back(vtkm::Vec3f(0.7257543206f, 0.1277695596f, 0.7468645573f)); + fishEndPts.push_back(vtkm::Vec3f(0.8347796798f, 0.1276152730f, 0.4985143244f)); + + vtkm::FloatDefault fishStep = 0.001f; + std::string fishFile = basePath + "/rectilinear/fishtank.vtk"; + TestParticleAdvectionFile(fishFile, fishPts, fishStep, 100, fishEndPts); } int UnitTestParticleAdvection(int argc, char* argv[])