1 //============================================================================ 2 // Copyright (c) Kitware, Inc. 3 // All rights reserved. 4 // See LICENSE.txt for details. 5 // This software is distributed WITHOUT ANY WARRANTY; without even 6 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 7 // PURPOSE. See the above copyright notice for more information. 8 // 9 // Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS). 10 // Copyright 2014 UT-Battelle, LLC. 11 // Copyright 2014 Los Alamos National Security. 12 // 13 // Under the terms of Contract DE-NA0003525 with NTESS, 14 // the U.S. Government retains certain rights in this software. 15 // 16 // Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National 17 // Laboratory (LANL), the U.S. Government retains certain rights in 18 // this software. 19 //============================================================================ 20 21 #ifndef vtk_m_worklet_Gradient_h 22 #define vtk_m_worklet_Gradient_h 23 24 #include <vtkm/worklet/DispatcherMapTopology.h> 25 #include <vtkm/worklet/DispatcherPointNeighborhood.h> 26 27 #include <vtkm/worklet/gradient/CellGradient.h> 28 #include <vtkm/worklet/gradient/Divergence.h> 29 #include <vtkm/worklet/gradient/GradientOutput.h> 30 #include <vtkm/worklet/gradient/PointGradient.h> 31 #include <vtkm/worklet/gradient/QCriterion.h> 32 #include <vtkm/worklet/gradient/StructuredPointGradient.h> 33 #include <vtkm/worklet/gradient/Transpose.h> 34 #include <vtkm/worklet/gradient/Vorticity.h> 35 36 namespace vtkm 37 { 38 namespace worklet 39 { 40 41 template <typename T> 42 struct GradientOutputFields; 43 44 namespace gradient 45 { 46 47 //----------------------------------------------------------------------------- 48 template <typename CoordinateSystem, typename T, typename S, typename Device> 49 struct DeducedPointGrad 50 { DeducedPointGradDeducedPointGrad51 DeducedPointGrad(const CoordinateSystem& coords, 52 const vtkm::cont::ArrayHandle<T, S>& field, 53 GradientOutputFields<T>* result) 54 : Points(&coords) 55 , Field(&field) 56 , Result(result) 57 { 58 } 59 60 template <typename CellSetType> operatorDeducedPointGrad61 void operator()(const CellSetType& cellset) const 62 { 63 vtkm::worklet::DispatcherMapTopology<PointGradient<T>> dispatcher; 64 dispatcher.SetDevice(Device()); 65 dispatcher.Invoke(cellset, //topology to iterate on a per point basis 66 cellset, //whole cellset in 67 *this->Points, 68 *this->Field, 69 *this->Result); 70 } 71 operatorDeducedPointGrad72 void operator()(const vtkm::cont::CellSetStructured<3>& cellset) const 73 { 74 vtkm::worklet::DispatcherPointNeighborhood<StructuredPointGradient<T>> dispatcher; 75 dispatcher.SetDevice(Device()); 76 dispatcher.Invoke(cellset, //topology to iterate on a per point basis 77 *this->Points, 78 *this->Field, 79 *this->Result); 80 } 81 82 template <typename PermIterType> operatorDeducedPointGrad83 void operator()(const vtkm::cont::CellSetPermutation<vtkm::cont::CellSetStructured<3>, 84 PermIterType>& cellset) const 85 { 86 vtkm::worklet::DispatcherPointNeighborhood<StructuredPointGradient<T>> dispatcher; 87 dispatcher.SetDevice(Device()); 88 dispatcher.Invoke(cellset, //topology to iterate on a per point basis 89 *this->Points, 90 *this->Field, 91 *this->Result); 92 } 93 operatorDeducedPointGrad94 void operator()(const vtkm::cont::CellSetStructured<2>& cellset) const 95 { 96 vtkm::worklet::DispatcherPointNeighborhood<StructuredPointGradient<T>> dispatcher; 97 dispatcher.SetDevice(Device()); 98 dispatcher.Invoke(cellset, //topology to iterate on a per point basis 99 *this->Points, 100 *this->Field, 101 *this->Result); 102 } 103 104 template <typename PermIterType> operatorDeducedPointGrad105 void operator()(const vtkm::cont::CellSetPermutation<vtkm::cont::CellSetStructured<2>, 106 PermIterType>& cellset) const 107 { 108 vtkm::worklet::DispatcherPointNeighborhood<StructuredPointGradient<T>> dispatcher; 109 dispatcher.SetDevice(Device()); 110 dispatcher.Invoke(cellset, //topology to iterate on a per point basis 111 *this->Points, 112 *this->Field, 113 *this->Result); 114 } 115 116 117 const CoordinateSystem* const Points; 118 const vtkm::cont::ArrayHandle<T, S>* const Field; 119 GradientOutputFields<T>* Result; 120 121 private: 122 void operator=(const DeducedPointGrad<CoordinateSystem, T, S, Device>&) = delete; 123 }; 124 125 } //namespace gradient 126 127 template <typename T> 128 struct GradientOutputFields : public vtkm::cont::ExecutionObjectBase 129 { 130 131 using ValueType = T; 132 using BaseTType = typename vtkm::BaseComponent<T>::Type; 133 134 template <typename DeviceAdapter> 135 struct ExecutionTypes 136 { 137 using Portal = vtkm::exec::GradientOutput<T>; 138 }; 139 GradientOutputFieldsGradientOutputFields140 GradientOutputFields() 141 : Gradient() 142 , Divergence() 143 , Vorticity() 144 , QCriterion() 145 , StoreGradient(true) 146 , ComputeDivergence(false) 147 , ComputeVorticity(false) 148 , ComputeQCriterion(false) 149 { 150 } 151 GradientOutputFieldsGradientOutputFields152 GradientOutputFields(bool store, bool divergence, bool vorticity, bool qc) 153 : Gradient() 154 , Divergence() 155 , Vorticity() 156 , QCriterion() 157 , StoreGradient(store) 158 , ComputeDivergence(divergence) 159 , ComputeVorticity(vorticity) 160 , ComputeQCriterion(qc) 161 { 162 } 163 164 /// Add divergence field to the output data. 165 /// The input array must have 3 components in order to compute this. 166 /// The default is off. SetComputeDivergenceGradientOutputFields167 void SetComputeDivergence(bool enable) { ComputeDivergence = enable; } GetComputeDivergenceGradientOutputFields168 bool GetComputeDivergence() const { return ComputeDivergence; } 169 170 /// Add voriticity/curl field to the output data. 171 /// The input array must have 3 components in order to compute this. 172 /// The default is off. SetComputeVorticityGradientOutputFields173 void SetComputeVorticity(bool enable) { ComputeVorticity = enable; } GetComputeVorticityGradientOutputFields174 bool GetComputeVorticity() const { return ComputeVorticity; } 175 176 /// Add Q-criterion field to the output data. 177 /// The input array must have 3 components in order to compute this. 178 /// The default is off. SetComputeQCriterionGradientOutputFields179 void SetComputeQCriterion(bool enable) { ComputeQCriterion = enable; } GetComputeQCriterionGradientOutputFields180 bool GetComputeQCriterion() const { return ComputeQCriterion; } 181 182 /// Add gradient field to the output data. 183 /// The input array must have 3 components in order to disable this. 184 /// The default is on. SetComputeGradientGradientOutputFields185 void SetComputeGradient(bool enable) { StoreGradient = enable; } GetComputeGradientGradientOutputFields186 bool GetComputeGradient() const { return StoreGradient; } 187 188 //todo fix this for scalar PrepareForOutputGradientOutputFields189 vtkm::exec::GradientOutput<T> PrepareForOutput(vtkm::Id size) 190 { 191 vtkm::exec::GradientOutput<T> portal(this->StoreGradient, 192 this->ComputeDivergence, 193 this->ComputeVorticity, 194 this->ComputeQCriterion, 195 this->Gradient, 196 this->Divergence, 197 this->Vorticity, 198 this->QCriterion, 199 size); 200 return portal; 201 } 202 203 vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> Gradient; 204 vtkm::cont::ArrayHandle<BaseTType> Divergence; 205 vtkm::cont::ArrayHandle<vtkm::Vec<BaseTType, 3>> Vorticity; 206 vtkm::cont::ArrayHandle<BaseTType> QCriterion; 207 208 private: 209 bool StoreGradient; 210 bool ComputeDivergence; 211 bool ComputeVorticity; 212 bool ComputeQCriterion; 213 }; 214 class PointGradient 215 { 216 public: 217 template <typename CellSetType, 218 typename CoordinateSystem, 219 typename T, 220 typename S, 221 typename DeviceAdapter> Run(const CellSetType & cells,const CoordinateSystem & coords,const vtkm::cont::ArrayHandle<T,S> & field,DeviceAdapter device)222 vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> Run(const CellSetType& cells, 223 const CoordinateSystem& coords, 224 const vtkm::cont::ArrayHandle<T, S>& field, 225 DeviceAdapter device) 226 { 227 vtkm::worklet::GradientOutputFields<T> extraOutput(true, false, false, false); 228 return this->Run(cells, coords, field, extraOutput, device); 229 } 230 231 template <typename CellSetType, 232 typename CoordinateSystem, 233 typename T, 234 typename S, 235 typename DeviceAdapter> Run(const CellSetType & cells,const CoordinateSystem & coords,const vtkm::cont::ArrayHandle<T,S> & field,GradientOutputFields<T> & extraOutput,DeviceAdapter)236 vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> Run(const CellSetType& cells, 237 const CoordinateSystem& coords, 238 const vtkm::cont::ArrayHandle<T, S>& field, 239 GradientOutputFields<T>& extraOutput, 240 DeviceAdapter) 241 { 242 //we are using cast and call here as we pass the cells twice to the invoke 243 //and want the type resolved once before hand instead of twice 244 //by the dispatcher ( that will cost more in time and binary size ) 245 gradient::DeducedPointGrad<CoordinateSystem, T, S, DeviceAdapter> func( 246 coords, field, &extraOutput); 247 vtkm::cont::CastAndCall(cells, func); 248 return extraOutput.Gradient; 249 } 250 }; 251 252 class CellGradient 253 { 254 public: 255 template <typename CellSetType, 256 typename CoordinateSystem, 257 typename T, 258 typename S, 259 typename DeviceAdapter> Run(const CellSetType & cells,const CoordinateSystem & coords,const vtkm::cont::ArrayHandle<T,S> & field,DeviceAdapter device)260 vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> Run(const CellSetType& cells, 261 const CoordinateSystem& coords, 262 const vtkm::cont::ArrayHandle<T, S>& field, 263 DeviceAdapter device) 264 { 265 vtkm::worklet::GradientOutputFields<T> extra(true, false, false, false); 266 return this->Run(cells, coords, field, extra, device); 267 } 268 269 template <typename CellSetType, 270 typename CoordinateSystem, 271 typename T, 272 typename S, 273 typename DeviceAdapter> Run(const CellSetType & cells,const CoordinateSystem & coords,const vtkm::cont::ArrayHandle<T,S> & field,GradientOutputFields<T> & extraOutput,DeviceAdapter)274 vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> Run(const CellSetType& cells, 275 const CoordinateSystem& coords, 276 const vtkm::cont::ArrayHandle<T, S>& field, 277 GradientOutputFields<T>& extraOutput, 278 DeviceAdapter) 279 { 280 using DispatcherType = 281 vtkm::worklet::DispatcherMapTopology<vtkm::worklet::gradient::CellGradient<T>>; 282 283 vtkm::worklet::gradient::CellGradient<T> worklet; 284 DispatcherType dispatcher(worklet); 285 dispatcher.SetDevice(DeviceAdapter()); 286 287 dispatcher.Invoke(cells, coords, field, extraOutput); 288 return extraOutput.Gradient; 289 } 290 }; 291 } 292 } // namespace vtkm::worklet 293 #endif 294