diff --git a/vtkm/filter/Lagrangian.hxx b/vtkm/filter/Lagrangian.hxx index 8d464a84a..b9b129cf0 100644 --- a/vtkm/filter/Lagrangian.hxx +++ b/vtkm/filter/Lagrangian.hxx @@ -25,9 +25,9 @@ #include #include #include -#include #include #include +#include #include #include @@ -279,12 +279,14 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute( using FieldType = vtkm::worklet::particleadvection::VelocityField; using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator; using RK4Type = vtkm::worklet::particleadvection::RK4Integrator; + using Stepper = vtkm::worklet::particleadvection::Stepper; + vtkm::worklet::ParticleAdvection particleadvection; vtkm::worklet::ParticleAdvectionResult res; FieldType velocities(field); GridEvalType gridEval(coords, cells, velocities); - RK4Type rk4(gridEval, static_cast(this->stepSize)); + Stepper rk4(gridEval, static_cast(this->stepSize)); res = particleadvection.Run(rk4, basisParticleArray, 1); // Taking a single step auto particles = res.Particles; diff --git a/vtkm/filter/LagrangianStructures.h b/vtkm/filter/LagrangianStructures.h index 7ecd3a6a9..9d1c2293c 100644 --- a/vtkm/filter/LagrangianStructures.h +++ b/vtkm/filter/LagrangianStructures.h @@ -17,7 +17,7 @@ #include #include -#include +#include namespace vtkm { diff --git a/vtkm/filter/LagrangianStructures.hxx b/vtkm/filter/LagrangianStructures.hxx index 80b234552..962f4a30b 100644 --- a/vtkm/filter/LagrangianStructures.hxx +++ b/vtkm/filter/LagrangianStructures.hxx @@ -16,9 +16,9 @@ #include #include #include -#include #include #include +#include #include @@ -86,7 +86,8 @@ inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute( using FieldHandle = vtkm::cont::ArrayHandle, StorageType>; using FieldType = vtkm::worklet::particleadvection::VelocityField; using GridEvaluator = vtkm::worklet::particleadvection::GridEvaluator; - using Integrator = vtkm::worklet::particleadvection::RK4Integrator; + using IntegratorType = vtkm::worklet::particleadvection::RK4Integrator; + using Stepper = vtkm::worklet::particleadvection::Stepper; vtkm::FloatDefault stepSize = this->GetStepSize(); vtkm::Id numberOfSteps = this->GetNumberOfSteps(); @@ -138,7 +139,7 @@ inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute( FieldType velocities(field); GridEvaluator evaluator(input.GetCoordinateSystem(), input.GetCellSet(), velocities); - Integrator integrator(evaluator, stepSize); + Stepper integrator(evaluator, stepSize); vtkm::worklet::ParticleAdvection particles; vtkm::worklet::ParticleAdvectionResult advectionResult; vtkm::cont::ArrayHandle advectionPoints; diff --git a/vtkm/filter/Pathline.h b/vtkm/filter/Pathline.h index 8800b42c0..baed55be8 100644 --- a/vtkm/filter/Pathline.h +++ b/vtkm/filter/Pathline.h @@ -14,7 +14,7 @@ #include #include #include -#include +#include namespace vtkm { diff --git a/vtkm/filter/StreamSurface.h b/vtkm/filter/StreamSurface.h index 3c9b8ccf0..bfb2dcdc4 100644 --- a/vtkm/filter/StreamSurface.h +++ b/vtkm/filter/StreamSurface.h @@ -15,7 +15,7 @@ #include #include #include -#include +#include namespace vtkm { diff --git a/vtkm/filter/StreamSurface.hxx b/vtkm/filter/StreamSurface.hxx index 75f4e8ac4..d4512e581 100644 --- a/vtkm/filter/StreamSurface.hxx +++ b/vtkm/filter/StreamSurface.hxx @@ -21,6 +21,7 @@ #include #include #include +#include namespace vtkm { @@ -57,11 +58,12 @@ inline VTKM_CONT vtkm::cont::DataSet StreamSurface::DoExecute( using FieldType = vtkm::worklet::particleadvection::VelocityField; using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator; using RK4Type = vtkm::worklet::particleadvection::RK4Integrator; + using Stepper = vtkm::worklet::particleadvection::Stepper; //compute streamlines FieldType velocities(field); GridEvalType eval(coords, cells, velocities); - RK4Type rk4(eval, this->StepSize); + Stepper rk4(eval, this->StepSize); vtkm::worklet::Streamline streamline; diff --git a/vtkm/filter/particleadvection/DataSetIntegrator.h b/vtkm/filter/particleadvection/DataSetIntegrator.h index 4db0a4326..4069d5701 100644 --- a/vtkm/filter/particleadvection/DataSetIntegrator.h +++ b/vtkm/filter/particleadvection/DataSetIntegrator.h @@ -15,6 +15,7 @@ #include #include #include +#include #include #include @@ -51,17 +52,18 @@ public: { auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off); auto seedArray = vtkm::cont::make_ArrayHandle(v, copyFlag); - RK4Type rk4(*this->Eval, stepSize); + Stepper rk4(*this->Eval, stepSize); this->DoAdvect(seedArray, rk4, maxSteps, result); } protected: using RK4Type = vtkm::worklet::particleadvection::RK4Integrator; + using Stepper = vtkm::worklet::particleadvection::Stepper; using FieldHandleType = vtkm::cont::ArrayHandle; template inline void DoAdvect(vtkm::cont::ArrayHandle& seeds, - const RK4Type& rk4, + const Stepper& rk4, vtkm::Id maxSteps, ResultType& result) const; diff --git a/vtkm/filter/particleadvection/DataSetIntegrator.hxx b/vtkm/filter/particleadvection/DataSetIntegrator.hxx index fce782b23..620c2fdbe 100644 --- a/vtkm/filter/particleadvection/DataSetIntegrator.hxx +++ b/vtkm/filter/particleadvection/DataSetIntegrator.hxx @@ -23,7 +23,10 @@ using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator< using TemporalGridEvalType = vtkm::worklet::particleadvection::TemporalGridEvaluator< vtkm::worklet::particleadvection::VelocityField>>; using RK4Type = vtkm::worklet::particleadvection::RK4Integrator; +using Stepper = vtkm::worklet::particleadvection::Stepper; using TemporalRK4Type = vtkm::worklet::particleadvection::RK4Integrator; +using TemporalStepper = + vtkm::worklet::particleadvection::Stepper; //----- // Specialization for ParticleAdvection worklet @@ -31,7 +34,7 @@ template <> template <> inline void DataSetIntegratorBase::DoAdvect( vtkm::cont::ArrayHandle& seeds, - const RK4Type& rk4, + const Stepper& rk4, vtkm::Id maxSteps, vtkm::worklet::ParticleAdvectionResult& result) const { @@ -45,7 +48,7 @@ template <> template <> inline void DataSetIntegratorBase::DoAdvect( vtkm::cont::ArrayHandle& seeds, - const RK4Type& rk4, + const Stepper& rk4, vtkm::Id maxSteps, vtkm::worklet::StreamlineResult& result) const { @@ -59,7 +62,7 @@ template <> template <> inline void DataSetIntegratorBase::DoAdvect( vtkm::cont::ArrayHandle& seeds, - const TemporalRK4Type& rk4, + const TemporalStepper& rk4, vtkm::Id maxSteps, vtkm::worklet::StreamlineResult& result) const { diff --git a/vtkm/worklet/particleadvection/CMakeLists.txt b/vtkm/worklet/particleadvection/CMakeLists.txt index d78d4935e..9afc70707 100644 --- a/vtkm/worklet/particleadvection/CMakeLists.txt +++ b/vtkm/worklet/particleadvection/CMakeLists.txt @@ -14,7 +14,7 @@ set(headers Field.h GridEvaluators.h GridEvaluatorStatus.h - IntegratorBase.h + Stepper.h IntegratorStatus.h Particles.h ParticleAdvectionWorklets.h diff --git a/vtkm/worklet/particleadvection/CellInterpolationHelper.h b/vtkm/worklet/particleadvection/CellInterpolationHelper.h index 930edf32d..4426a6abc 100644 --- a/vtkm/worklet/particleadvection/CellInterpolationHelper.h +++ b/vtkm/worklet/particleadvection/CellInterpolationHelper.h @@ -31,158 +31,144 @@ namespace vtkm namespace exec { -class CellInterpolationHelper : public vtkm::VirtualObjectBase +class CellInterpolationHelper { -public: - VTKM_EXEC_CONT virtual ~CellInterpolationHelper() noexcept override - { - // This must not be defaulted, since defaulted virtual destructors are - // troublesome with CUDA __host__ __device__ markup. - } - - VTKM_EXEC - virtual void GetCellInfo(const vtkm::Id& cellId, - vtkm::UInt8& cellShape, - vtkm::IdComponent& numVerts, - vtkm::VecVariable& indices) const = 0; -}; - -class StructuredCellInterpolationHelper : public vtkm::exec::CellInterpolationHelper -{ -public: - StructuredCellInterpolationHelper() = default; - - VTKM_CONT - StructuredCellInterpolationHelper(vtkm::Id3 cellDims, vtkm::Id3 pointDims, bool is3D) - : CellDims(cellDims) - , PointDims(pointDims) - , Is3D(is3D) - { - } - - VTKM_EXEC - void GetCellInfo(const vtkm::Id& cellId, - vtkm::UInt8& cellShape, - vtkm::IdComponent& numVerts, - vtkm::VecVariable& indices) const override - { - vtkm::Id3 logicalCellId; - logicalCellId[0] = cellId % this->CellDims[0]; - logicalCellId[1] = (cellId / this->CellDims[0]) % this->CellDims[1]; - if (this->Is3D) - { - logicalCellId[2] = cellId / (this->CellDims[0] * this->CellDims[1]); - indices.Append((logicalCellId[2] * this->PointDims[1] + logicalCellId[1]) * - this->PointDims[0] + - logicalCellId[0]); - indices.Append(indices[0] + 1); - indices.Append(indices[1] + this->PointDims[0]); - indices.Append(indices[2] - 1); - indices.Append(indices[0] + this->PointDims[0] * this->PointDims[1]); - indices.Append(indices[4] + 1); - indices.Append(indices[5] + this->PointDims[0]); - indices.Append(indices[6] - 1); - cellShape = static_cast(vtkm::CELL_SHAPE_HEXAHEDRON); - numVerts = 8; - } - else - { - indices.Append(logicalCellId[1] * this->PointDims[0] + logicalCellId[0]); - indices.Append(indices[0] + 1); - indices.Append(indices[1] + this->PointDims[0]); - indices.Append(indices[2] - 1); - cellShape = static_cast(vtkm::CELL_SHAPE_QUAD); - numVerts = 4; - } - } - private: - vtkm::Id3 CellDims; - vtkm::Id3 PointDims; - bool Is3D = true; -}; - -class SingleCellTypeInterpolationHelper : public vtkm::exec::CellInterpolationHelper -{ - using ConnType = vtkm::cont::ArrayHandle; - using ConnPortalType = typename ConnType::ReadPortalType; - -public: - SingleCellTypeInterpolationHelper() = default; - - VTKM_CONT - SingleCellTypeInterpolationHelper(vtkm::UInt8 cellShape, - vtkm::IdComponent pointsPerCell, - const ConnType& connectivity, - vtkm::cont::DeviceAdapterId device, - vtkm::cont::Token& token) - : CellShape(cellShape) - , PointsPerCell(pointsPerCell) - , Connectivity(connectivity.PrepareForInput(device, token)) - { - } - - VTKM_EXEC - void GetCellInfo(const vtkm::Id& cellId, - vtkm::UInt8& cellShape, - vtkm::IdComponent& numVerts, - vtkm::VecVariable& indices) const override - { - cellShape = CellShape; - numVerts = PointsPerCell; - vtkm::Id n = static_cast(PointsPerCell); - vtkm::Id offset = cellId * n; - - for (vtkm::Id i = 0; i < n; i++) - indices.Append(Connectivity.Get(offset + i)); - } - -private: - vtkm::UInt8 CellShape; - vtkm::IdComponent PointsPerCell; - ConnPortalType Connectivity; -}; - -class ExplicitCellInterpolationHelper : public vtkm::exec::CellInterpolationHelper -{ using ShapeType = vtkm::cont::ArrayHandle; using OffsetType = vtkm::cont::ArrayHandle; using ConnType = vtkm::cont::ArrayHandle; - using ShapePortalType = typename ShapeType::ReadPortalType; using OffsetPortalType = typename OffsetType::ReadPortalType; using ConnPortalType = typename ConnType::ReadPortalType; public: - ExplicitCellInterpolationHelper() = default; + enum class HelperType + { + STRUCTURED, + EXPSINGLE, + EXPLICIT + }; VTKM_CONT - ExplicitCellInterpolationHelper(const ShapeType& shape, - const OffsetType& offset, - const ConnType& connectivity, - vtkm::cont::DeviceAdapterId device, - vtkm::cont::Token& token) + CellInterpolationHelper() = default; + + VTKM_CONT + CellInterpolationHelper(const vtkm::Id3& cellDims, const vtkm::Id3& pointDims, bool is3D) + : CellDims(cellDims) + , PointDims(pointDims) + , Is3D(is3D) + { + this->Type = HelperType::STRUCTURED; + } + + CellInterpolationHelper(const vtkm::UInt8 cellShape, + const vtkm::IdComponent pointsPerCell, + const ConnType& connectivity, + vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) + : CellShape(cellShape) + , PointsPerCell(pointsPerCell) + , Connectivity(connectivity.PrepareForInput(device, token)) + { + this->Type = HelperType::EXPSINGLE; + } + + + VTKM_CONT + CellInterpolationHelper(const ShapeType& shape, + const OffsetType& offset, + const ConnType& connectivity, + vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) : Shape(shape.PrepareForInput(device, token)) , Offset(offset.PrepareForInput(device, token)) , Connectivity(connectivity.PrepareForInput(device, token)) { + this->Type = HelperType::EXPLICIT; } VTKM_EXEC void GetCellInfo(const vtkm::Id& cellId, vtkm::UInt8& cellShape, vtkm::IdComponent& numVerts, - vtkm::VecVariable& indices) const override + vtkm::VecVariable& indices) const { - cellShape = this->Shape.Get(cellId); - const vtkm::Id offset = this->Offset.Get(cellId); - numVerts = static_cast(this->Offset.Get(cellId + 1) - offset); + switch (this->Type) + { + case HelperType::STRUCTURED: + { + vtkm::Id3 logicalCellId; + logicalCellId[0] = cellId % this->CellDims[0]; + logicalCellId[1] = (cellId / this->CellDims[0]) % this->CellDims[1]; + if (this->Is3D) + { + logicalCellId[2] = cellId / (this->CellDims[0] * this->CellDims[1]); + indices.Append((logicalCellId[2] * this->PointDims[1] + logicalCellId[1]) * + this->PointDims[0] + + logicalCellId[0]); + indices.Append(indices[0] + 1); + indices.Append(indices[1] + this->PointDims[0]); + indices.Append(indices[2] - 1); + indices.Append(indices[0] + this->PointDims[0] * this->PointDims[1]); + indices.Append(indices[4] + 1); + indices.Append(indices[5] + this->PointDims[0]); + indices.Append(indices[6] - 1); + cellShape = static_cast(vtkm::CELL_SHAPE_HEXAHEDRON); + numVerts = 8; + } + else + { + indices.Append(logicalCellId[1] * this->PointDims[0] + logicalCellId[0]); + indices.Append(indices[0] + 1); + indices.Append(indices[1] + this->PointDims[0]); + indices.Append(indices[2] - 1); + cellShape = static_cast(vtkm::CELL_SHAPE_QUAD); + numVerts = 4; + } + } + break; - for (vtkm::IdComponent i = 0; i < numVerts; i++) - indices.Append(this->Connectivity.Get(offset + i)); + case HelperType::EXPSINGLE: + { + cellShape = this->CellShape; + numVerts = this->PointsPerCell; + vtkm::Id n = static_cast(PointsPerCell); + vtkm::Id offset = cellId * n; + for (vtkm::Id i = 0; i < n; i++) + indices.Append(Connectivity.Get(offset + i)); + } + break; + + case HelperType::EXPLICIT: + { + cellShape = this->Shape.Get(cellId); + const vtkm::Id offset = this->Offset.Get(cellId); + numVerts = static_cast(this->Offset.Get(cellId + 1) - offset); + for (vtkm::IdComponent i = 0; i < numVerts; i++) + indices.Append(this->Connectivity.Get(offset + i)); + } + break; + + default: + { + // Code path not expected to execute in correct cases + // Supress unused variable warning + cellShape = vtkm::UInt8(0); + numVerts = vtkm::IdComponent(0); + } + } } private: + HelperType Type; + // variables for structured type + vtkm::Id3 CellDims; + vtkm::Id3 PointDims; + bool Is3D = true; + // variables for single explicit type + vtkm::UInt8 CellShape; + vtkm::IdComponent PointsPerCell; + // variables for explicit type ShapePortalType Shape; OffsetPortalType Offset; ConnPortalType Connectivity; @@ -198,26 +184,19 @@ namespace cont class CellInterpolationHelper : public vtkm::cont::ExecutionObjectBase { -public: - using HandleType = vtkm::cont::VirtualObjectHandle; - - virtual ~CellInterpolationHelper() = default; - - VTKM_CONT virtual const vtkm::exec::CellInterpolationHelper* PrepareForExecution( - vtkm::cont::DeviceAdapterId device, - vtkm::cont::Token& token) const = 0; -}; - -class StructuredCellInterpolationHelper : public vtkm::cont::CellInterpolationHelper -{ -public: +private: + using ExecutionType = vtkm::exec::CellInterpolationHelper; using Structured2DType = vtkm::cont::CellSetStructured<2>; using Structured3DType = vtkm::cont::CellSetStructured<3>; + using SingleExplicitType = vtkm::cont::CellSetSingleType<>; + using ExplicitType = vtkm::cont::CellSetExplicit<>; - StructuredCellInterpolationHelper() = default; +public: + VTKM_CONT + CellInterpolationHelper() = default; VTKM_CONT - StructuredCellInterpolationHelper(const vtkm::cont::DynamicCellSet& cellSet) + CellInterpolationHelper(const vtkm::cont::DynamicCellSet& cellSet) { if (cellSet.IsSameType(Structured2DType())) { @@ -228,6 +207,7 @@ public: cellSet.Cast().GetSchedulingRange(vtkm::TopologyElementTagPoint()); this->CellDims = vtkm::Id3(cellDims[0], cellDims[1], 0); this->PointDims = vtkm::Id3(pointDims[0], pointDims[1], 1); + this->Type = vtkm::exec::CellInterpolationHelper::HelperType::STRUCTURED; } else if (cellSet.IsSameType(Structured3DType())) { @@ -236,113 +216,22 @@ public: cellSet.Cast().GetSchedulingRange(vtkm::TopologyElementTagCell()); this->PointDims = cellSet.Cast().GetSchedulingRange(vtkm::TopologyElementTagPoint()); + this->Type = vtkm::exec::CellInterpolationHelper::HelperType::STRUCTURED; } - else - throw vtkm::cont::ErrorBadType("Cell set is not of type CellSetStructured"); - } - - VTKM_CONT - const vtkm::exec::CellInterpolationHelper* PrepareForExecution( - vtkm::cont::DeviceAdapterId deviceId, - vtkm::cont::Token& token) const override - { - auto& tracker = vtkm::cont::GetRuntimeDeviceTracker(); - const bool valid = tracker.CanRunOn(deviceId); - if (!valid) - { - throwFailedRuntimeDeviceTransfer("StructuredCellInterpolationHelper", deviceId); - } - - using ExecutionType = vtkm::exec::StructuredCellInterpolationHelper; - ExecutionType* execObject = new ExecutionType(this->CellDims, this->PointDims, this->Is3D); - this->ExecHandle.Reset(execObject); - - return this->ExecHandle.PrepareForExecution(deviceId, token); - } - -private: - vtkm::Id3 CellDims; - vtkm::Id3 PointDims; - bool Is3D = true; - mutable HandleType ExecHandle; -}; - -class SingleCellTypeInterpolationHelper : public vtkm::cont::CellInterpolationHelper -{ -public: - using SingleExplicitType = vtkm::cont::CellSetSingleType<>; - - SingleCellTypeInterpolationHelper() = default; - - VTKM_CONT - SingleCellTypeInterpolationHelper(const vtkm::cont::DynamicCellSet& cellSet) - { - if (cellSet.IsSameType(SingleExplicitType())) + else if (cellSet.IsSameType(SingleExplicitType())) { SingleExplicitType CellSet = cellSet.Cast(); - const auto cellShapes = CellSet.GetShapesArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint()); const auto numIndices = CellSet.GetNumIndicesArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint()); - CellShape = vtkm::cont::ArrayGetValue(0, cellShapes); PointsPerCell = vtkm::cont::ArrayGetValue(0, numIndices); Connectivity = CellSet.GetConnectivityArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint()); + this->Type = vtkm::exec::CellInterpolationHelper::HelperType::EXPSINGLE; } - else - throw vtkm::cont::ErrorBadType("Cell set is not of type CellSetSingleType"); - } - - struct SingleCellTypeFunctor - { - VTKM_CONT bool operator()(vtkm::cont::DeviceAdapterId device, - const vtkm::cont::SingleCellTypeInterpolationHelper& contInterpolator, - HandleType& execInterpolator, - vtkm::cont::Token& token) const - { - using ExecutionType = vtkm::exec::SingleCellTypeInterpolationHelper; - ExecutionType* execObject = new ExecutionType(contInterpolator.CellShape, - contInterpolator.PointsPerCell, - contInterpolator.Connectivity, - device, - token); - execInterpolator.Reset(execObject); - return true; - } - }; - - VTKM_CONT - const vtkm::exec::CellInterpolationHelper* PrepareForExecution( - vtkm::cont::DeviceAdapterId deviceId, - vtkm::cont::Token& token) const override - { - const bool success = vtkm::cont::TryExecuteOnDevice( - deviceId, SingleCellTypeFunctor(), *this, this->ExecHandle, token); - if (!success) - { - throwFailedRuntimeDeviceTransfer("SingleCellTypeInterpolationHelper", deviceId); - } - return this->ExecHandle.PrepareForExecution(deviceId, token); - } - -private: - vtkm::UInt8 CellShape; - vtkm::IdComponent PointsPerCell; - vtkm::cont::ArrayHandle Connectivity; - mutable HandleType ExecHandle; -}; - -class ExplicitCellInterpolationHelper : public vtkm::cont::CellInterpolationHelper -{ -public: - ExplicitCellInterpolationHelper() = default; - - VTKM_CONT - ExplicitCellInterpolationHelper(const vtkm::cont::DynamicCellSet& cellSet) - { - if (cellSet.IsSameType(vtkm::cont::CellSetExplicit<>())) + else if (cellSet.IsSameType(ExplicitType())) { vtkm::cont::CellSetExplicit<> CellSet = cellSet.Cast>(); Shape = @@ -351,28 +240,44 @@ public: CellSet.GetOffsetsArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint()); Connectivity = CellSet.GetConnectivityArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint()); + this->Type = vtkm::exec::CellInterpolationHelper::HelperType::EXPLICIT; } else - throw vtkm::cont::ErrorBadType("Cell set is not of type CellSetExplicit"); + throw vtkm::cont::ErrorInternal("Unsupported cellset type"); } VTKM_CONT - const vtkm::exec::CellInterpolationHelper* PrepareForExecution( - vtkm::cont::DeviceAdapterId device, - vtkm::cont::Token& token) const override + const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) const { - using ExecutionType = vtkm::exec::ExplicitCellInterpolationHelper; - ExecutionType* execObject = - new ExecutionType(this->Shape, this->Offset, this->Connectivity, device, token); - this->ExecHandle.Reset(execObject); - return this->ExecHandle.PrepareForExecution(device, token); + switch (this->Type) + { + case ExecutionType::HelperType::STRUCTURED: + return ExecutionType(this->CellDims, this->PointDims, this->Is3D); + + case ExecutionType::HelperType::EXPSINGLE: + return ExecutionType( + this->CellShape, this->PointsPerCell, this->Connectivity, device, token); + + case ExecutionType::HelperType::EXPLICIT: + return ExecutionType(this->Shape, this->Offset, this->Connectivity, device, token); + } + throw vtkm::cont::ErrorInternal("Undefined case for building cell interpolation helper"); } private: + // Variables required for strucutred grids + vtkm::Id3 CellDims; + vtkm::Id3 PointDims; + bool Is3D = true; + // variables for single explicit type + vtkm::UInt8 CellShape; + vtkm::IdComponent PointsPerCell; + // Variables required for unstructured grids vtkm::cont::ArrayHandle Shape; vtkm::cont::ArrayHandle Offset; vtkm::cont::ArrayHandle Connectivity; - mutable HandleType ExecHandle; + ExecutionType::HelperType Type; }; } //namespace cont diff --git a/vtkm/worklet/particleadvection/EulerIntegrator.h b/vtkm/worklet/particleadvection/EulerIntegrator.h index 26990643d..2c909156a 100644 --- a/vtkm/worklet/particleadvection/EulerIntegrator.h +++ b/vtkm/worklet/particleadvection/EulerIntegrator.h @@ -13,8 +13,6 @@ #ifndef vtk_m_worklet_particleadvection_EulerIntegrator_h #define vtk_m_worklet_particleadvection_EulerIntegrator_h -#include - namespace vtkm { namespace worklet @@ -22,73 +20,56 @@ namespace worklet namespace particleadvection { -template -class EulerIntegrator : public IntegratorBase +template +class ExecEulerIntegrator { +public: + VTKM_EXEC_CONT + ExecEulerIntegrator(const EvaluatorType& evaluator) + : Evaluator(evaluator) + { + } + + template + VTKM_EXEC IntegratorStatus CheckStep(Particle& particle, + vtkm::FloatDefault stepLength, + vtkm::Vec3f& velocity) const + { + auto time = particle.Time; + auto inpos = particle.Pos; + vtkm::VecVariable vectors; + GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors); + if (status.CheckOk()) + velocity = particle.Velocity(vectors, stepLength); + return IntegratorStatus(status); + } + +private: + EvaluatorType Evaluator; +}; + +template +class EulerIntegrator +{ +private: + EvaluatorType Evaluator; + public: EulerIntegrator() = default; VTKM_CONT - EulerIntegrator(const FieldEvaluateType& evaluator, const vtkm::FloatDefault stepLength) - : IntegratorBase(stepLength) - , Evaluator(evaluator) + EulerIntegrator(const EvaluatorType& evaluator) + : Evaluator(evaluator) { } - template - class ExecObject - : public IntegratorBase::ExecObjectBaseImpl< - vtkm::cont::internal::ExecutionObjectType, - typename EulerIntegrator::template ExecObject> + VTKM_CONT auto PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) const + -> ExecEulerIntegratorEvaluator.PrepareForExecution(device, token))> { - VTKM_IS_DEVICE_ADAPTER_TAG(Device); - - using FieldEvaluateExecType = - vtkm::cont::internal::ExecutionObjectType; - using Superclass = - IntegratorBase::ExecObjectBaseImpl>; - - public: - VTKM_EXEC_CONT - ExecObject(const FieldEvaluateExecType& evaluator, - vtkm::FloatDefault stepLength, - vtkm::FloatDefault tolerance) - : Superclass(evaluator, stepLength, tolerance) - { - } - - VTKM_EXEC - IntegratorStatus CheckStep(vtkm::Particle* particle, - vtkm::FloatDefault stepLength, - vtkm::Vec3f& velocity) const - { - auto time = particle->Time; - auto inpos = particle->Pos; - vtkm::VecVariable vectors; - GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors); - if (status.CheckOk()) - velocity = particle->Velocity(vectors, stepLength); - return IntegratorStatus(status); - } - }; - -private: - FieldEvaluateType Evaluator; - -protected: - VTKM_CONT virtual void PrepareForExecutionImpl( - vtkm::cont::DeviceAdapterId device, - vtkm::cont::VirtualObjectHandle& execObjectHandle, - vtkm::cont::Token& token) const override - { - vtkm::cont::TryExecuteOnDevice(device, - detail::IntegratorPrepareForExecutionFunctor(), - execObjectHandle, - this->Evaluator, - this->StepLength, - this->Tolerance, - token); + auto evaluator = this->Evaluator.PrepareForExecution(device, token); + using ExecEvaluatorType = decltype(evaluator); + return ExecEulerIntegrator(evaluator); } }; //EulerIntegrator diff --git a/vtkm/worklet/particleadvection/Field.h b/vtkm/worklet/particleadvection/Field.h index d776ede20..d2f94699f 100644 --- a/vtkm/worklet/particleadvection/Field.h +++ b/vtkm/worklet/particleadvection/Field.h @@ -26,22 +26,8 @@ namespace worklet namespace particleadvection { -class ExecutionField : public vtkm::VirtualObjectBase -{ -public: - VTKM_EXEC_CONT - virtual ~ExecutionField() noexcept override {} - - VTKM_EXEC - virtual void GetValue(const vtkm::VecVariable& indices, - const vtkm::Id vertices, - const vtkm::Vec3f& parametric, - const vtkm::UInt8 cellShape, - vtkm::VecVariable& value) const = 0; -}; - template -class ExecutionVelocityField : public vtkm::worklet::particleadvection::ExecutionField +class ExecutionVelocityField { public: using FieldPortalType = typename FieldArrayType::ReadPortalType; @@ -73,7 +59,7 @@ private: }; template -class ExecutionElectroMagneticField : public vtkm::worklet::particleadvection::ExecutionField +class ExecutionElectroMagneticField { public: using FieldPortalType = typename FieldArrayType::ReadPortalType; @@ -112,22 +98,15 @@ private: FieldPortalType MagneticValues; }; -class Field : public vtkm::cont::ExecutionObjectBase +template +class VelocityField : public vtkm::cont::ExecutionObjectBase { public: - using HandleType = vtkm::cont::VirtualObjectHandle; - - virtual ~Field() = default; + using ExecutionType = ExecutionVelocityField; VTKM_CONT - virtual const ExecutionField* PrepareForExecution(vtkm::cont::DeviceAdapterId deviceId, - vtkm::cont::Token& token) const = 0; -}; + VelocityField() = default; -template -class VelocityField : public vtkm::worklet::particleadvection::Field -{ -public: VTKM_CONT VelocityField(const FieldArrayType& fieldValues) : FieldValues(fieldValues) @@ -135,24 +114,25 @@ public: } VTKM_CONT - const ExecutionField* PrepareForExecution(vtkm::cont::DeviceAdapterId device, - vtkm::cont::Token& token) const override + const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) const { - using ExecutionType = ExecutionVelocityField; - ExecutionType* execObject = new ExecutionType(this->FieldValues, device, token); - this->ExecHandle.Reset(execObject); - return this->ExecHandle.PrepareForExecution(device, token); + return ExecutionType(this->FieldValues, device, token); } private: FieldArrayType FieldValues; - mutable HandleType ExecHandle; }; template -class ElectroMagneticField : public vtkm::worklet::particleadvection::Field +class ElectroMagneticField : public vtkm::cont::ExecutionObjectBase { public: + using ExecutionType = ExecutionElectroMagneticField; + + VTKM_CONT + ElectroMagneticField() = default; + VTKM_CONT ElectroMagneticField(const FieldArrayType& electricField, const FieldArrayType& magneticField) : ElectricField(electricField) @@ -161,20 +141,15 @@ public: } VTKM_CONT - const ExecutionField* PrepareForExecution(vtkm::cont::DeviceAdapterId device, - vtkm::cont::Token& token) const override + const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) const { - using ExecutionType = ExecutionElectroMagneticField; - ExecutionType* execObject = - new ExecutionType(this->ElectricField, this->MagneticField, device, token); - this->ExecHandle.Reset(execObject); - return this->ExecHandle.PrepareForExecution(device, token); + return ExecutionType(this->ElectricField, this->MagneticField, device, token); } private: FieldArrayType ElectricField; FieldArrayType MagneticField; - mutable HandleType ExecHandle; }; } // namespace particleadvection diff --git a/vtkm/worklet/particleadvection/GridEvaluatorStatus.h b/vtkm/worklet/particleadvection/GridEvaluatorStatus.h index c3b238944..44182cec1 100644 --- a/vtkm/worklet/particleadvection/GridEvaluatorStatus.h +++ b/vtkm/worklet/particleadvection/GridEvaluatorStatus.h @@ -22,9 +22,6 @@ #include #include -#include -#include - namespace vtkm { namespace worklet diff --git a/vtkm/worklet/particleadvection/GridEvaluators.h b/vtkm/worklet/particleadvection/GridEvaluators.h index 177c95758..e886eb5e4 100644 --- a/vtkm/worklet/particleadvection/GridEvaluators.h +++ b/vtkm/worklet/particleadvection/GridEvaluators.h @@ -26,7 +26,6 @@ #include #include #include -#include namespace vtkm { @@ -46,7 +45,7 @@ public: VTKM_CONT ExecutionGridEvaluator(const vtkm::cont::CellLocatorGeneral& locator, - std::shared_ptr interpolationHelper, + const vtkm::cont::CellInterpolationHelper interpolationHelper, const vtkm::Bounds& bounds, const FieldType& field, const GhostCellArrayType& ghostCells, @@ -56,7 +55,7 @@ public: , Field(field.PrepareForExecution(device, token)) , GhostCells(ghostCells.PrepareForInput(device, token)) , HaveGhostCells(ghostCells.GetNumberOfValues() > 0) - , InterpolationHelper(interpolationHelper->PrepareForExecution(device, token)) + , InterpolationHelper(interpolationHelper.PrepareForExecution(device, token)) , Locator(locator.PrepareForExecution(device, token)) { } @@ -125,9 +124,8 @@ public: vtkm::IdComponent nVerts; vtkm::VecVariable ptIndices; vtkm::VecVariable fieldValues; - InterpolationHelper->GetCellInfo(cellId, cellShape, nVerts, ptIndices); - - this->Field->GetValue(ptIndices, nVerts, parametric, cellShape, out); + this->InterpolationHelper.GetCellInfo(cellId, cellShape, nVerts, ptIndices); + this->Field.GetValue(ptIndices, nVerts, parametric, cellShape, out); status.SetOk(); } @@ -146,10 +144,10 @@ private: using GhostCellPortal = typename vtkm::cont::ArrayHandle::ReadPortalType; vtkm::Bounds Bounds; - const vtkm::worklet::particleadvection::ExecutionField* Field; + typename FieldType::ExecutionType Field; GhostCellPortal GhostCells; bool HaveGhostCells; - const vtkm::exec::CellInterpolationHelper* InterpolationHelper; + vtkm::exec::CellInterpolationHelper InterpolationHelper; typename vtkm::cont::CellLocatorGeneral::ExecObjType Locator; }; @@ -217,32 +215,13 @@ private: this->Locator.SetCoordinates(coordinates); this->Locator.SetCellSet(cellset); this->Locator.Update(); - if (cellset.IsSameType(Structured2DType()) || cellset.IsSameType(Structured3DType())) - { - vtkm::cont::StructuredCellInterpolationHelper interpolationHelper(cellset); - this->InterpolationHelper = - std::make_shared(interpolationHelper); - } - else if (cellset.IsSameType(vtkm::cont::CellSetSingleType<>())) - { - vtkm::cont::SingleCellTypeInterpolationHelper interpolationHelper(cellset); - this->InterpolationHelper = - std::make_shared(interpolationHelper); - } - else if (cellset.IsSameType(vtkm::cont::CellSetExplicit<>())) - { - vtkm::cont::ExplicitCellInterpolationHelper interpolationHelper(cellset); - this->InterpolationHelper = - std::make_shared(interpolationHelper); - } - else - throw vtkm::cont::ErrorInternal("Unsupported cellset type."); + this->InterpolationHelper = vtkm::cont::CellInterpolationHelper(cellset); } vtkm::Bounds Bounds; FieldType Field; GhostCellArrayType GhostCellArray; - std::shared_ptr InterpolationHelper; + vtkm::cont::CellInterpolationHelper InterpolationHelper; vtkm::cont::CellLocatorGeneral Locator; }; diff --git a/vtkm/worklet/particleadvection/IntegratorBase.h b/vtkm/worklet/particleadvection/IntegratorBase.h deleted file mode 100644 index 608e4d52d..000000000 --- a/vtkm/worklet/particleadvection/IntegratorBase.h +++ /dev/null @@ -1,252 +0,0 @@ -//============================================================================= -// -// 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_worklet_particleadvection_IntegratorBase_h -#define vtk_m_worklet_particleadvection_IntegratorBase_h - -#include - -#include -#include -#include -#include - -#include -#include - -#include -#include -#include - -namespace vtkm -{ -namespace worklet -{ -namespace particleadvection -{ - -class IntegratorBase : public vtkm::cont::ExecutionObjectBase -{ -protected: - VTKM_CONT - IntegratorBase() = default; - - VTKM_CONT - IntegratorBase(vtkm::FloatDefault stepLength) - : StepLength(stepLength) - { - } - -public: - class ExecObject : public vtkm::VirtualObjectBase - { - protected: - VTKM_EXEC_CONT - ExecObject(const vtkm::FloatDefault stepLength, vtkm::FloatDefault tolerance) - : StepLength(stepLength) - , Tolerance(tolerance) - { - } - - public: - VTKM_EXEC - virtual IntegratorStatus Step(vtkm::Particle* inpos, - vtkm::FloatDefault& time, - vtkm::Vec3f& outpos) const = 0; - - VTKM_EXEC - virtual IntegratorStatus SmallStep(vtkm::Particle* inpos, - vtkm::FloatDefault& time, - vtkm::Vec3f& outpos) const = 0; - - protected: - vtkm::FloatDefault StepLength = 1.0f; - vtkm::FloatDefault Tolerance = 0.001f; - }; - - template - VTKM_CONT const ExecObject* PrepareForExecution(Device, vtkm::cont::Token& token) const - { - this->PrepareForExecutionImpl( - Device(), - const_cast&>(this->ExecObjectHandle), - token); - return this->ExecObjectHandle.PrepareForExecution(Device(), token); - } - -private: - vtkm::cont::VirtualObjectHandle ExecObjectHandle; - -protected: - vtkm::FloatDefault StepLength; - vtkm::FloatDefault Tolerance = - std::numeric_limits::epsilon() * static_cast(100.0f); - - VTKM_CONT virtual void PrepareForExecutionImpl( - vtkm::cont::DeviceAdapterId device, - vtkm::cont::VirtualObjectHandle& execObjectHandle, - vtkm::cont::Token& token) const = 0; - - template - class ExecObjectBaseImpl : public ExecObject - { - protected: - VTKM_EXEC_CONT - ExecObjectBaseImpl(const FieldEvaluateType& evaluator, - vtkm::FloatDefault stepLength, - vtkm::FloatDefault tolerance) - : ExecObject(stepLength, tolerance) - , Evaluator(evaluator) - { - } - - public: - VTKM_EXEC - IntegratorStatus Step(vtkm::Particle* particle, - vtkm::FloatDefault& time, - vtkm::Vec3f& outpos) const override - { - vtkm::Vec3f velocity(0, 0, 0); - auto status = this->CheckStep(particle, this->StepLength, velocity); - if (status.CheckOk()) - { - outpos = particle->Pos + this->StepLength * velocity; - time += this->StepLength; - } - else - outpos = particle->Pos; - - return status; - } - - VTKM_EXEC - IntegratorStatus SmallStep(vtkm::Particle* particle, - vtkm::FloatDefault& time, - vtkm::Vec3f& outpos) const override - { - //Stepping by this->StepLength goes beyond the bounds of the dataset. - //We need to take an Euler step that goes outside of the dataset. - //Use a binary search to find the largest step INSIDE the dataset. - //Binary search uses a shrinking bracket of inside / outside, so when - //we terminate, the outside value is the stepsize that will nudge - //the particle outside the dataset. - - //The binary search will be between {0, this->StepLength} - vtkm::FloatDefault stepRange[2] = { 0, this->StepLength }; - - vtkm::Vec3f currPos(particle->Pos); - vtkm::Vec3f currVelocity(0, 0, 0); - vtkm::VecVariable currValue, tmp; - auto evalStatus = this->Evaluator.Evaluate(currPos, particle->Time, currValue); - if (evalStatus.CheckFail()) - return IntegratorStatus(evalStatus); - - const vtkm::FloatDefault eps = vtkm::Epsilon() * 10; - vtkm::FloatDefault div = 1; - while ((stepRange[1] - stepRange[0]) > eps) - { - //Try a step midway between stepRange[0] and stepRange[1] - div *= 2; - vtkm::FloatDefault currStep = stepRange[0] + (this->StepLength / div); - - //See if we can step by currStep - IntegratorStatus status = this->CheckStep(particle, currStep, currVelocity); - - if (status.CheckOk()) //Integration step succedded. - { - //See if this point is in/out. - auto newPos = particle->Pos + currStep * currVelocity; - evalStatus = this->Evaluator.Evaluate(newPos, particle->Time + currStep, tmp); - if (evalStatus.CheckOk()) - { - //Point still in. Update currPos and set range to {currStep, stepRange[1]} - currPos = newPos; - stepRange[0] = currStep; - } - else - { - //The step succedded, but the next point is outside. - //Step too long. Set range to: {stepRange[0], currStep} and continue. - stepRange[1] = currStep; - } - } - else - { - //Step too long. Set range to: {stepRange[0], stepCurr} and continue. - stepRange[1] = currStep; - } - } - - evalStatus = this->Evaluator.Evaluate(currPos, particle->Time + stepRange[0], currValue); - //The eval at Time + stepRange[0] better be *inside* - VTKM_ASSERT(evalStatus.CheckOk() && !evalStatus.CheckSpatialBounds()); - if (evalStatus.CheckFail() || evalStatus.CheckSpatialBounds()) - return IntegratorStatus(evalStatus); - - //Update the position and time. - outpos = currPos + stepRange[1] * particle->Velocity(currValue, stepRange[1]); - time += stepRange[1]; - - //Get the evaluation status for the point that is *just* outside of the data. - evalStatus = this->Evaluator.Evaluate(outpos, time, currValue); - - //The eval should fail, and the point should be outside either spatially or temporally. - VTKM_ASSERT(evalStatus.CheckFail() && - (evalStatus.CheckSpatialBounds() || evalStatus.CheckTemporalBounds())); - - IntegratorStatus status(evalStatus); - status.SetOk(); //status is ok. - - return status; - } - - VTKM_EXEC - IntegratorStatus CheckStep(vtkm::Particle* particle, - vtkm::FloatDefault stepLength, - vtkm::Vec3f& velocity) const - { - return static_cast(this)->CheckStep(particle, stepLength, velocity); - } - - protected: - FieldEvaluateType Evaluator; - }; -}; - -namespace detail -{ - -template