diff --git a/data/baseline/filter/amrArrays2D.png b/data/baseline/filter/amrArrays2D.png new file mode 100644 index 000000000..bc32f1d96 --- /dev/null +++ b/data/baseline/filter/amrArrays2D.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0994e1747beaddba1ab6e97a48e17956e64a12ebfdffac1c2ae8734f3578956f +size 20716 diff --git a/data/baseline/filter/amrArrays3D.png b/data/baseline/filter/amrArrays3D.png new file mode 100644 index 000000000..bca21202b --- /dev/null +++ b/data/baseline/filter/amrArrays3D.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6a610fe2de049c6a2402f7b19ea5a03f2e637e90d62a6df38832e62b89a691bb +size 25953 diff --git a/vtkm/Bounds.h b/vtkm/Bounds.h index b72547e8d..cebb602ae 100644 --- a/vtkm/Bounds.h +++ b/vtkm/Bounds.h @@ -102,6 +102,42 @@ struct Bounds return (this->X.Contains(point[0]) && this->Y.Contains(point[1]) && this->Z.Contains(point[2])); } + /// \b Returns the volume of the bounds. + /// + /// \c Volume computes the product of the lengths of the ranges in each dimension. If the bounds + /// are empty, 0 is returned. + /// + VTKM_EXEC_CONT + vtkm::Float64 Volume() const + { + if (this->IsNonEmpty()) + { + return (this->X.Length() * this->Y.Length() * this->Z.Length()); + } + else + { + return 0.0; + } + } + + /// \b Returns the area of the bounds in the X-Y-plane. + /// + /// \c Area computes the product of the lengths of the ranges in dimensions X and Y. If the bounds + /// are empty, 0 is returned. + /// + VTKM_EXEC_CONT + vtkm::Float64 Area() const + { + if (this->IsNonEmpty()) + { + return (this->X.Length() * this->Y.Length()); + } + else + { + return 0.0; + } + } + /// \b Returns the center of the range. /// /// \c Center computes the point at the middle of the bounds. If the bounds @@ -152,6 +188,16 @@ struct Bounds return unionBounds; } + /// \b Return the intersection of this and another range. + /// + VTKM_EXEC_CONT + vtkm::Bounds Intersection(const vtkm::Bounds& otherBounds) const + { + return vtkm::Bounds(this->X.Intersection(otherBounds.X), + this->Y.Intersection(otherBounds.Y), + this->Z.Intersection(otherBounds.Z)); + } + /// \b Operator for union /// VTKM_EXEC_CONT diff --git a/vtkm/CellClassification.h b/vtkm/CellClassification.h index ccd414695..6ab306c6e 100644 --- a/vtkm/CellClassification.h +++ b/vtkm/CellClassification.h @@ -19,7 +19,7 @@ enum CellClassification : vtkm::UInt8 GHOST = 1 << 0, //Ghost cell INVALID = 1 << 1, //Cell is invalid UNUSED0 = 1 << 2, - UNUSED1 = 1 << 3, + BLANKED = 1 << 3, //Blanked cell in AMR UNUSED3 = 1 << 4, UNUSED4 = 1 << 5, UNUSED5 = 1 << 6, diff --git a/vtkm/Range.h b/vtkm/Range.h index 0be2c3557..27f45e9dc 100644 --- a/vtkm/Range.h +++ b/vtkm/Range.h @@ -153,6 +153,14 @@ struct Range return unionRange; } + /// \b Return the intersection of this and another range. + /// + VTKM_EXEC_CONT + vtkm::Range Intersection(const vtkm::Range& otherRange) const + { + return vtkm::Range(vtkm::Max(this->Min, otherRange.Min), vtkm::Min(this->Max, otherRange.Max)); + } + /// \b Operator for union /// VTKM_EXEC_CONT diff --git a/vtkm/filter/AmrArrays.h b/vtkm/filter/AmrArrays.h new file mode 100644 index 000000000..b706953e1 --- /dev/null +++ b/vtkm/filter/AmrArrays.h @@ -0,0 +1,86 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +#ifndef vtk_m_filter_AmrArrays_h +#define vtk_m_filter_AmrArrays_h + +#include + +namespace vtkm +{ +namespace filter +{ + +class AmrArrays : public vtkm::filter::FilterDataSet +{ +public: + using SupportedTypes = vtkm::List; + + VTKM_CONT + AmrArrays(); + + template + vtkm::cont::PartitionedDataSet PrepareForExecution( + const vtkm::cont::PartitionedDataSet& input, + const vtkm::filter::PolicyBase&); + + template + VTKM_CONT bool MapFieldOntoOutput(vtkm::cont::DataSet& result, + const vtkm::cont::Field& field, + const vtkm::filter::PolicyBase&) + { + result.AddField(field); + return true; + } + +private: + /// the list of ids contains all amrIds of the level above/below that have an overlap + VTKM_CONT + void GenerateParentChildInformation(); + + /// the corresponding template function based on the dimension of this dataset + VTKM_CONT + template + void ComputeGenerateParentChildInformation(); + + /// generate the vtkGhostType array based on the overlap analogously to vtk + /// blanked cells: 8 normal cells: 0 + VTKM_CONT + void GenerateGhostType(); + + /// the corresponding template function based on the dimension of this dataset + VTKM_CONT + template + void ComputeGenerateGhostType(); + + /// Add helper arrays like in ParaView + VTKM_CONT + void GenerateIndexArrays(); + + /// the input partitioned dataset + vtkm::cont::PartitionedDataSet AmrDataSet; + + /// per level + /// contains the index where the PartitionIds start for each level + std::vector> PartitionIds; + + /// per partitionId + /// contains all PartitonIds of the level above that have an overlap + std::vector> ParentsIdsVector; + + /// per partitionId + /// contains all PartitonIds of the level below that have an overlap + std::vector> ChildrenIdsVector; +}; +} +} // namespace vtkm::filter + +#include + +#endif //vtk_m_filter_AmrArrays_h diff --git a/vtkm/filter/AmrArrays.hxx b/vtkm/filter/AmrArrays.hxx new file mode 100644 index 000000000..667dd86e2 --- /dev/null +++ b/vtkm/filter/AmrArrays.hxx @@ -0,0 +1,307 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +#ifndef vtk_m_filter_AmrArrays_hxx +#define vtk_m_filter_AmrArrays_hxx + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace vtkm +{ +namespace worklet +{ +template +struct GenerateGhostTypeWorklet : vtkm::worklet::WorkletVisitCellsWithPoints +{ + using ControlSignature = void(CellSetIn cellSet, + FieldInPoint pointArray, + FieldInOutCell ghostArray); + using ExecutionSignature = void(PointCount, _2, _3); + using InputDomain = _1; + + GenerateGhostTypeWorklet(vtkm::Bounds boundsChild) + : BoundsChild(boundsChild) + { + } + + template + VTKM_EXEC void operator()(vtkm::IdComponent numPoints, + const pointArrayType pointArray, + cellArrayType& ghostArray) const + { + vtkm::Bounds boundsCell = vtkm::Bounds(); + for (vtkm::IdComponent pointId = 0; pointId < numPoints; pointId++) + { + boundsCell.Include(pointArray[pointId]); + } + vtkm::Bounds boundsIntersection = boundsCell.Intersection(BoundsChild); + if ((Dim == 2 && boundsIntersection.Area() > 0.5 * boundsCell.Area()) || + (Dim == 3 && boundsIntersection.Volume() > 0.5 * boundsCell.Volume())) + { + // std::cout<AmrDataSet); + if (bounds.Z.Max - bounds.Z.Min < vtkm::Epsilon()) + { + ComputeGenerateParentChildInformation<2>(); + } + else + { + ComputeGenerateParentChildInformation<3>(); + } +} + +template void AmrArrays::ComputeGenerateParentChildInformation<2>(); + +template void AmrArrays::ComputeGenerateParentChildInformation<3>(); + +VTKM_CONT +template +void AmrArrays::ComputeGenerateParentChildInformation() +{ + // read out spacings in decreasing order to infer levels + std::set> spacings; + for (vtkm::Id p = 0; p < this->AmrDataSet.GetNumberOfPartitions(); p++) + { + vtkm::cont::ArrayHandleUniformPointCoordinates uniformCoords = + this->AmrDataSet.GetPartition(p) + .GetCoordinateSystem() + .GetData() + .AsArrayHandle(); + spacings.insert(uniformCoords.GetSpacing()[0]); + } + std::set>::iterator itr; + // for (itr = spacings.begin(); itr != spacings.end(); itr++) + // { + // std::cout << *itr << "\n"; + // } + + /// contains the partitionIds of each level and blockId + this->PartitionIds.resize(spacings.size()); + for (vtkm::Id p = 0; p < this->AmrDataSet.GetNumberOfPartitions(); p++) + { + vtkm::cont::ArrayHandleUniformPointCoordinates uniformCoords = + this->AmrDataSet.GetPartition(p) + .GetCoordinateSystem() + .GetData() + .AsArrayHandle(); + int index = -1; + for (itr = spacings.begin(); itr != spacings.end(); itr++) + { + index++; + if (*itr == uniformCoords.GetSpacing()[0]) + { + break; + } + } + this->PartitionIds.at(index).push_back(p); + // std::cout <ParentsIdsVector.resize(this->AmrDataSet.GetNumberOfPartitions()); + this->ChildrenIdsVector.resize(this->AmrDataSet.GetNumberOfPartitions()); + for (unsigned int l = 0; l < this->PartitionIds.size() - 1; l++) + { + for (unsigned int bParent = 0; bParent < this->PartitionIds.at(l).size(); bParent++) + { + // std::cout << std::endl << "level " << l << " block " << bParent << std::endl; + vtkm::Bounds boundsParent = + this->AmrDataSet.GetPartition(this->PartitionIds.at(l).at(bParent)) + .GetCoordinateSystem() + .GetBounds(); + + // compute size of a cell to compare overlap against + auto coords = this->AmrDataSet.GetPartition(this->PartitionIds.at(l).at(bParent)) + .GetCoordinateSystem() + .GetDataAsMultiplexer(); + vtkm::cont::CellSetStructured cellset; + vtkm::Id ptids[8]; + this->AmrDataSet.GetPartition(this->PartitionIds.at(l).at(bParent)) + .GetCellSet() + .CopyTo(cellset); + cellset.GetCellPointIds(0, ptids); + vtkm::Bounds boundsCell = vtkm::Bounds(); + for (vtkm::IdComponent pointId = 0; pointId < cellset.GetNumberOfPointsInCell(0); pointId++) + { + boundsCell.Include(coords.ReadPortal().Get(ptids[pointId])); + } + + // see if there is overlap of at least one half of a cell + for (unsigned int bChild = 0; bChild < this->PartitionIds.at(l + 1).size(); bChild++) + { + vtkm::Bounds boundsChild = + this->AmrDataSet.GetPartition(this->PartitionIds.at(l + 1).at(bChild)) + .GetCoordinateSystem() + .GetBounds(); + vtkm::Bounds boundsIntersection = boundsParent.Intersection(boundsChild); + if ((Dim == 2 && boundsIntersection.Area() > 0.5 * boundsCell.Area()) || + (Dim == 3 && boundsIntersection.Volume() >= 0.5 * boundsCell.Volume())) + { + this->ParentsIdsVector.at(this->PartitionIds.at(l + 1).at(bChild)) + .push_back(this->PartitionIds.at(l).at(bParent)); + this->ChildrenIdsVector.at(this->PartitionIds.at(l).at(bParent)) + .push_back(this->PartitionIds.at(l + 1).at(bChild)); + // std::cout << " overlaps with level " << l + 1 << " block " << bChild << " " + // << boundsParent << " " << boundsChild << " " << boundsIntersection << " " + // << boundsIntersection.Area() << " " << boundsIntersection.Volume() << std::endl; + } + // else + // { + // std::cout << " does not overlap with level " << l + 1 << " block " << bChild << " " + // << boundsParent << " " << boundsChild << " " << boundsIntersection << " " + // << boundsIntersection.Area() << " " << boundsIntersection.Volume() << std::endl; + // } + } + } + } +} + +VTKM_CONT +void AmrArrays::GenerateGhostType() +{ + vtkm::Bounds bounds = vtkm::cont::BoundsCompute(this->AmrDataSet); + if (bounds.Z.Max - bounds.Z.Min < vtkm::Epsilon()) + { + ComputeGenerateGhostType<2>(); + } + else + { + ComputeGenerateGhostType<3>(); + } +} + +template void AmrArrays::ComputeGenerateGhostType<2>(); + +template void AmrArrays::ComputeGenerateGhostType<3>(); + +VTKM_CONT +template +void AmrArrays::ComputeGenerateGhostType() +{ + for (unsigned int l = 0; l < this->PartitionIds.size(); l++) + { + for (unsigned int bParent = 0; bParent < this->PartitionIds.at(l).size(); bParent++) + { + // std::cout<ChildrenIdsVector.at(this->PartitionIds.at(l).at(bParent)).size()<<" children"<AmrDataSet.GetPartition(this->PartitionIds.at(l).at(bParent)); + vtkm::cont::CellSetStructured cellset; + partition.GetCellSet().CopyTo(cellset); + vtkm::cont::ArrayHandle ghostField; + if (!partition.HasField("vtkGhostType", vtkm::cont::Field::Association::CELL_SET)) + { + vtkm::cont::ArrayCopy( + vtkm::cont::ArrayHandleConstant(0, partition.GetNumberOfCells()), + ghostField); + } + // leave field unchanged if it is the highest level + if (l < this->PartitionIds.size() - 1) + { + auto pointField = partition.GetCoordinateSystem().GetDataAsMultiplexer(); + + for (unsigned int bChild = 0; + bChild < this->ChildrenIdsVector.at(this->PartitionIds.at(l).at(bParent)).size(); + bChild++) + { + vtkm::Bounds boundsChild = + this->AmrDataSet + .GetPartition( + this->ChildrenIdsVector.at(this->PartitionIds.at(l).at(bParent)).at(bChild)) + .GetCoordinateSystem() + .GetBounds(); + // std::cout<<" is (partly) contained in level "<ChildrenIdsVector.at(this->PartitionIds.at(l).at(bParent)).at(bChild)<<" with bounds "<{ boundsChild }, + cellset, + pointField, + ghostField); + } + } + partition.AddCellField("vtkGhostType", ghostField); + this->AmrDataSet.ReplacePartition(this->PartitionIds.at(l).at(bParent), partition); + } + } +} + +// Add helper arrays like in ParaView +VTKM_CONT +void AmrArrays::GenerateIndexArrays() +{ + for (unsigned int l = 0; l < this->PartitionIds.size(); l++) + { + for (unsigned int b = 0; b < this->PartitionIds.at(l).size(); b++) + { + vtkm::cont::DataSet partition = this->AmrDataSet.GetPartition(this->PartitionIds.at(l).at(b)); + + vtkm::cont::ArrayHandle fieldAmrLevel; + vtkm::cont::ArrayCopy( + vtkm::cont::ArrayHandleConstant(l, partition.GetNumberOfCells()), fieldAmrLevel); + partition.AddCellField("vtkAmrLevel", fieldAmrLevel); + + vtkm::cont::ArrayHandle fieldBlockId; + vtkm::cont::ArrayCopy( + vtkm::cont::ArrayHandleConstant(b, partition.GetNumberOfCells()), fieldBlockId); + partition.AddCellField("vtkAmrIndex", fieldBlockId); + + vtkm::cont::ArrayHandle fieldPartitionIndex; + vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant( + this->PartitionIds.at(l).at(b), partition.GetNumberOfCells()), + fieldPartitionIndex); + partition.AddCellField("vtkCompositeIndex", fieldPartitionIndex); + + this->AmrDataSet.ReplacePartition(this->PartitionIds.at(l).at(b), partition); + } + } +} + +template +vtkm::cont::PartitionedDataSet AmrArrays::PrepareForExecution( + const vtkm::cont::PartitionedDataSet& input, + const vtkm::filter::PolicyBase&) +{ + this->AmrDataSet = input; + this->GenerateParentChildInformation(); + this->GenerateGhostType(); + this->GenerateIndexArrays(); + return this->AmrDataSet; +} + +} +} +#endif diff --git a/vtkm/filter/CMakeLists.txt b/vtkm/filter/CMakeLists.txt index bb78c89a3..cdb6727a4 100644 --- a/vtkm/filter/CMakeLists.txt +++ b/vtkm/filter/CMakeLists.txt @@ -77,6 +77,7 @@ set(common_sources_device ) set(extra_headers + AmrArrays.h CellSetConnectivity.h ClipWithField.h ClipWithImplicitFunction.h @@ -134,6 +135,7 @@ set(extra_headers ) set(extra_header_template_sources + AmrArrays.hxx CellSetConnectivity.hxx ClipWithField.hxx ClipWithImplicitFunction.hxx diff --git a/vtkm/filter/testing/CMakeLists.txt b/vtkm/filter/testing/CMakeLists.txt index 3195fd64d..a5453da38 100644 --- a/vtkm/filter/testing/CMakeLists.txt +++ b/vtkm/filter/testing/CMakeLists.txt @@ -98,6 +98,7 @@ if (VTKm_ENABLE_RENDERING) list(APPEND libraries vtkm_rendering) list(APPEND unit_tests + RegressionTestAmrArrays.cxx RegressionTestContourFilter.cxx RegressionTestPointTransform.cxx RegressionTestSliceFilter.cxx diff --git a/vtkm/filter/testing/RegressionTestAmrArrays.cxx b/vtkm/filter/testing/RegressionTestAmrArrays.cxx new file mode 100644 index 000000000..1c3a1821e --- /dev/null +++ b/vtkm/filter/testing/RegressionTestAmrArrays.cxx @@ -0,0 +1,78 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace +{ + +void TestAmrArraysExecute(int dim, int numberOfLevels, int cellsPerDimension) +{ + std::cout << "Generate Image for AMR" << std::endl; + + using M = vtkm::rendering::MapperRayTracer; + using C = vtkm::rendering::CanvasRayTracer; + using V3 = vtkm::rendering::View3D; + + // Generate AMR + vtkm::source::Amr source(dim, cellsPerDimension, numberOfLevels); + vtkm::cont::PartitionedDataSet amrDataSet = source.Execute(); + // std::cout << "amr " << std::endl; + // amrDataSet.PrintSummary(std::cout); + + // Remove blanked cells + vtkm::filter::Threshold threshold; + threshold.SetLowerThreshold(0); + threshold.SetUpperThreshold(1); + threshold.SetActiveField("vtkGhostType"); + vtkm::cont::PartitionedDataSet derivedDataSet = threshold.Execute(amrDataSet); + // std::cout << "derived " << std::endl; + // derivedDataSet.PrintSummary(std::cout); + + // Extract surface for efficient 3D pipeline + vtkm::filter::ExternalFaces surface; + surface.SetFieldsToPass("RTDataCells"); + derivedDataSet = surface.Execute(derivedDataSet); + + // Merge dataset + vtkm::cont::DataSet result = vtkm::cont::MergePartitionedDataSet(derivedDataSet); + // std::cout << "merged " << std::endl; + // result.PrintSummary(std::cout); + + vtkm::rendering::testing::RenderAndRegressionTest(result, + "RTDataCells", + vtkm::cont::ColorTable("inferno"), + "filter/amrArrays" + + std::to_string(dim) + "D.png", + false); +} + +void TestAmrArrays() +{ + int numberOfLevels = 5; + int cellsPerDimension = 6; + TestAmrArraysExecute(2, numberOfLevels, cellsPerDimension); + TestAmrArraysExecute(3, numberOfLevels, cellsPerDimension); +} +} // namespace + +int RegressionTestAmrArrays(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(TestAmrArrays, argc, argv); +} diff --git a/vtkm/source/Amr.cxx b/vtkm/source/Amr.cxx new file mode 100644 index 000000000..a0a1c437e --- /dev/null +++ b/vtkm/source/Amr.cxx @@ -0,0 +1,111 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ + + +#include +#include +#include +#include +#include + + +namespace vtkm +{ +namespace source +{ + +Amr::Amr(vtkm::IdComponent dimension, + vtkm::IdComponent cellsPerDimension, + vtkm::IdComponent numberOfLevels) + : Dimension(dimension) + , CellsPerDimension(cellsPerDimension) + , NumberOfLevels(numberOfLevels) +{ +} + +Amr::~Amr() = default; + +template +vtkm::cont::DataSet Amr::GenerateDataSet(unsigned int level, unsigned int amrIndex) const +{ + vtkm::Id3 extent = { vtkm::Id(this->CellsPerDimension / 2) }; + vtkm::Id3 dimensions = { this->CellsPerDimension + 1 }; + vtkm::Vec3f origin = { float(1. / pow(2, level) * amrIndex) }; + vtkm::Vec3f spacing = { float(1. / this->CellsPerDimension / pow(2, level)) }; + vtkm::Vec3f center = 0.5f - (origin + spacing * extent); + vtkm::Vec3f frequency = { 60.f, 30.f, 40.f }; + frequency = frequency * this->CellsPerDimension; + vtkm::FloatDefault deviation = 0.5f / this->CellsPerDimension; + + if (Dim == 2) + { + extent[2] = 0; + dimensions[2] = 1; + origin[2] = 0; + spacing[2] = 1; + center[2] = 0; + } + + vtkm::source::Wavelet waveletSource(-extent, extent); + waveletSource.SetOrigin(origin); + waveletSource.SetSpacing(spacing); + waveletSource.SetCenter(center); + waveletSource.SetFrequency(frequency); + waveletSource.SetStandardDeviation(deviation); + vtkm::cont::DataSet wavelet = waveletSource.Execute(); + + vtkm::filter::CellAverage cellAverage; + cellAverage.SetActiveField("RTData", vtkm::cont::Field::Association::POINTS); + cellAverage.SetOutputFieldName("RTDataCells"); + return cellAverage.Execute(wavelet); +} + +vtkm::cont::PartitionedDataSet Amr::Execute() const +{ + assert(this->CellsPerDimension > 1); + assert(this->CellsPerDimension % 2 == 0); + + // Generate AMR + std::vector> blocksPerLevel(this->NumberOfLevels); + unsigned int counter = 0; + for (unsigned int l = 0; l < blocksPerLevel.size(); l++) + { + for (unsigned int b = 0; b < pow(2, l); b++) + { + blocksPerLevel.at(l).push_back(counter++); + } + } + vtkm::cont::PartitionedDataSet amrDataSet; + + // Fill AMR with data from the wavelet + for (unsigned int l = 0; l < blocksPerLevel.size(); l++) + { + for (unsigned int b = 0; b < blocksPerLevel.at(l).size(); b++) + { + if (this->Dimension == 2) + { + amrDataSet.AppendPartition(this->GenerateDataSet<2>(l, b)); + } + else if (this->Dimension == 3) + { + amrDataSet.AppendPartition(this->GenerateDataSet<3>(l, b)); + } + } + } + + // Generate helper arrays + vtkm::filter::AmrArrays amrArrays; + amrDataSet = amrArrays.Execute(amrDataSet); + + return amrDataSet; +} + +} // namespace source +} // namespace vtkm diff --git a/vtkm/source/Amr.h b/vtkm/source/Amr.h new file mode 100644 index 000000000..bd451f84a --- /dev/null +++ b/vtkm/source/Amr.h @@ -0,0 +1,81 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ + +#ifndef vtk_m_source_Amr_h +#define vtk_m_source_Amr_h + +#include +#include + +namespace vtkm +{ +namespace source +{ +/** + * @brief The Amr source creates a dataset similar to VTK's + * vtkRTAnalyticSource. + * + * This class generates a predictable structured dataset with a smooth yet + * interesting set of scalars, which is useful for testing and benchmarking. + * + * The Execute method creates a complete structured dataset that have a + * point field names 'scalars' + * + * The scalars are computed as: + * + * ``` + * MaxVal * Gauss + MagX * sin(FrqX*x) + MagY * sin(FrqY*y) + MagZ * cos(FrqZ*z) + * ``` + * + * The dataset properties are determined by: + * - `Minimum/MaximumExtent`: The logical point extents of the dataset. + * - `Spacing`: The distance between points of the dataset. + * - `Center`: The center of the dataset. + * + * The scalar functions is control via: + * - `Center`: The center of a Gaussian contribution to the scalars. + * - `StandardDeviation`: The unscaled width of a Gaussian contribution. + * - `MaximumValue`: Upper limit of the scalar range. + * - `Frequency`: The Frq[XYZ] parameters of the periodic contributions. + * - `Magnitude`: The Mag[XYZ] parameters of the periodic contributions. + * + * By default, the following parameters are used: + * - `Extents`: { -10, -10, -10 } `-->` { 10, 10, 10 } + * - `Spacing`: { 1, 1, 1 } + * - `Center`: { 0, 0, 0 } + * - `StandardDeviation`: 0.5 + * - `MaximumValue`: 255 + * - `Frequency`: { 60, 30, 40 } + * - `Magnitude`: { 10, 18, 5 } + */ +class VTKM_SOURCE_EXPORT Amr +{ +public: + VTKM_CONT + Amr(vtkm::IdComponent dimension = 2, + vtkm::IdComponent cellsPerDimension = 6, + vtkm::IdComponent numberOfLevels = 4); + VTKM_CONT + ~Amr(); + + vtkm::cont::PartitionedDataSet Execute() const; + +private: + template + vtkm::cont::DataSet GenerateDataSet(unsigned int level, unsigned int amrIndex) const; + + vtkm::IdComponent Dimension; + vtkm::IdComponent CellsPerDimension; + vtkm::IdComponent NumberOfLevels; +}; +} //namespace source +} //namespace vtkm + +#endif //vtk_m_source_Amr_h diff --git a/vtkm/source/CMakeLists.txt b/vtkm/source/CMakeLists.txt index 03de54190..2a1c61461 100644 --- a/vtkm/source/CMakeLists.txt +++ b/vtkm/source/CMakeLists.txt @@ -9,6 +9,7 @@ ##============================================================================ set(headers + Amr.h Oscillator.h Source.h Tangle.h @@ -17,6 +18,7 @@ set(headers ) set(device_sources + Amr.cxx Oscillator.cxx Source.cxx Tangle.cxx @@ -25,11 +27,12 @@ set(device_sources ) vtkm_library(NAME vtkm_source + SOURCES ${sources} DEVICE_SOURCES ${device_sources} HEADERS ${headers} ) -target_link_libraries(vtkm_source PUBLIC vtkm_cont) +target_link_libraries(vtkm_source PUBLIC vtkm_cont vtkm_filter) #----------------------------------------------------------------------------- if (VTKm_ENABLE_TESTING) diff --git a/vtkm/source/Wavelet.cxx b/vtkm/source/Wavelet.cxx index 3eb1b4b06..88e09fa34 100644 --- a/vtkm/source/Wavelet.cxx +++ b/vtkm/source/Wavelet.cxx @@ -98,6 +98,7 @@ struct WaveletField : public vtkm::worklet::WorkletVisitPointsWithCells Wavelet::Wavelet(vtkm::Id3 minExtent, vtkm::Id3 maxExtent) : Center{ minExtent - ((minExtent - maxExtent) / 2) } + , Origin{ minExtent } , Spacing{ 1. } , Frequency{ 60., 30., 40. } , Magnitude{ 10., 18., 5. } @@ -108,17 +109,16 @@ Wavelet::Wavelet(vtkm::Id3 minExtent, vtkm::Id3 maxExtent) { } -vtkm::cont::DataSet Wavelet::Execute() const +template +vtkm::cont::DataSet Wavelet::GenerateDataSet(vtkm::cont::CoordinateSystem coords) const { - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - // Create points: - const vtkm::Id3 dims{ this->MaximumExtent - this->MinimumExtent + vtkm::Id3{ 1 } }; - const vtkm::Vec3f origin{ this->MinimumExtent }; - vtkm::cont::CoordinateSystem coords{ "coordinates", dims, origin, this->Spacing }; - // And cells: - vtkm::cont::CellSetStructured<3> cellSet; + vtkm::Vec dims; + for (unsigned int d = 0; d < Dim; d++) + { + dims[d] = this->MaximumExtent[d] - this->MinimumExtent[d] + 1; + } + vtkm::cont::CellSetStructured cellSet; cellSet.SetPointDimensions(dims); // Compile the dataset: @@ -133,7 +133,27 @@ vtkm::cont::DataSet Wavelet::Execute() const return dataSet; } -vtkm::cont::Field Wavelet::GeneratePointField(const vtkm::cont::CellSetStructured<3>& cellset, +vtkm::cont::DataSet Wavelet::Execute() const +{ + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + // Create points: + const vtkm::Id3 dims{ this->MaximumExtent - this->MinimumExtent + vtkm::Id3{ 1 } }; + vtkm::cont::CoordinateSystem coords{ "coordinates", dims, this->Origin, this->Spacing }; + + // Compile the dataset: + if (this->MaximumExtent[2] - this->MinimumExtent[2] < vtkm::Epsilon()) + { + return this->GenerateDataSet<2>(coords); + } + else + { + return this->GenerateDataSet<3>(coords); + } +} + +template +vtkm::cont::Field Wavelet::GeneratePointField(const vtkm::cont::CellSetStructured& cellset, const std::string& name) const { const vtkm::Id3 dims{ this->MaximumExtent - this->MinimumExtent + vtkm::Id3{ 1 } }; diff --git a/vtkm/source/Wavelet.h b/vtkm/source/Wavelet.h index 8bd30cb73..001d3b2af 100644 --- a/vtkm/source/Wavelet.h +++ b/vtkm/source/Wavelet.h @@ -53,6 +53,8 @@ namespace source * - `MaximumValue`: 255 * - `Frequency`: { 60, 30, 40 } * - `Magnitude`: { 10, 18, 5 } + * + * If the extent has zero length in the z-direction, a 2D dataset is generated. */ class VTKM_SOURCE_EXPORT Wavelet final : public vtkm::source::Source { @@ -62,6 +64,8 @@ public: VTKM_CONT void SetCenter(const vtkm::Vec& center) { this->Center = center; } + VTKM_CONT void SetOrigin(const vtkm::Vec& origin) { this->Origin = origin; } + VTKM_CONT void SetSpacing(const vtkm::Vec& spacing) { this->Spacing = spacing; } VTKM_CONT void SetFrequency(const vtkm::Vec& frequency) @@ -94,10 +98,15 @@ public: vtkm::cont::DataSet Execute() const; private: - vtkm::cont::Field GeneratePointField(const vtkm::cont::CellSetStructured<3>& cellset, + template + vtkm::cont::Field GeneratePointField(const vtkm::cont::CellSetStructured& cellset, const std::string& name) const; + template + vtkm::cont::DataSet GenerateDataSet(vtkm::cont::CoordinateSystem coords) const; + vtkm::Vec3f Center; + vtkm::Vec3f Origin; vtkm::Vec3f Spacing; vtkm::Vec3f Frequency; vtkm::Vec3f Magnitude;