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