VTK-m  2.0
worklet/LagrangianStructures.h
Go to the documentation of this file.
1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10 
11 #ifndef vtk_m_filter_flow_worklet_LagrangianStructures_h
12 #define vtk_m_filter_flow_worklet_LagrangianStructures_h
13 
14 #include <vtkm/Matrix.h>
15 #include <vtkm/Types.h>
18 
21 
22 namespace vtkm
23 {
24 namespace worklet
25 {
26 namespace flow
27 {
28 
29 template <vtkm::IdComponent dimensions>
31 
32 template <>
34 {
35 public:
37 
38  VTKM_CONT
40  : EndTime(endTime)
41  , GridData(cellSet)
42  {
43  }
44 
45  using ControlSignature = void(WholeArrayIn, WholeArrayIn, FieldOut);
46 
47  using ExecutionSignature = void(WorkIndex, _1, _2, _3);
48 
49  template <typename PointArray>
50  VTKM_EXEC void operator()(const vtkm::Id index,
51  const PointArray& input,
52  const PointArray& output,
53  Scalar& outputField) const
54  {
55  const vtkm::Vec<vtkm::Id, 6> neighborIndices = this->GridData.GetNeighborIndices(index);
56 
57  // Calculate Stretching / Squeezing
58  auto xin1 = input.Get(neighborIndices[0]);
59  auto xin2 = input.Get(neighborIndices[1]);
60  auto yin1 = input.Get(neighborIndices[2]);
61  auto yin2 = input.Get(neighborIndices[3]);
62 
63  Scalar xDiff = 1.0f / (xin2[0] - xin1[0]);
64  Scalar yDiff = 1.0f / (yin2[1] - yin1[1]);
65 
66  auto xout1 = output.Get(neighborIndices[0]);
67  auto xout2 = output.Get(neighborIndices[1]);
68  auto yout1 = output.Get(neighborIndices[2]);
69  auto yout2 = output.Get(neighborIndices[3]);
70 
71  // Total X gradient w.r.t X, Y
72  Scalar f1x = (xout2[0] - xout1[0]) * xDiff;
73  Scalar f1y = (yout2[0] - yout1[0]) * yDiff;
74 
75  // Total Y gradient w.r.t X, Y
76  Scalar f2x = (xout2[1] - xout1[1]) * xDiff;
77  Scalar f2y = (yout2[1] - yout1[1]) * yDiff;
78 
80  vtkm::MatrixSetRow(jacobian, 0, vtkm::Vec<Scalar, 2>(f1x, f1y));
81  vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec<Scalar, 2>(f2x, f2y));
82 
83  vtkm::filter::flow::internal::ComputeLeftCauchyGreenTensor(jacobian);
84 
85  vtkm::Vec<Scalar, 2> eigenValues;
86  vtkm::filter::flow::internal::Jacobi(jacobian, eigenValues);
87 
88  Scalar delta = eigenValues[0];
89  // Check if we need to clamp these values
90  // Also provide options.
91  // 1. FTLE
92  // 2. FLLE
93  // 3. Eigen Values (Min/Max)
94  //Scalar delta = trace + sqrtr;
95  // Given endTime is in units where start time is 0,
96  // else do endTime-startTime
97  // return value for computation
98  outputField = vtkm::Log(delta) / (static_cast<Scalar>(2.0f) * EndTime);
99  }
100 
101 public:
102  // To calculate FTLE field
104  // To assist in calculation of indices
105  vtkm::filter::flow::internal::GridMetaData GridData;
106 };
107 
108 template <>
110 {
111 public:
113 
114  VTKM_CONT
116  : EndTime(endTime)
117  , GridData(cellSet)
118  {
119  }
120 
121  using ControlSignature = void(WholeArrayIn, WholeArrayIn, FieldOut);
122 
123  using ExecutionSignature = void(WorkIndex, _1, _2, _3);
124 
125  /*
126  * Point position arrays are the input and the output positions of the particle advection.
127  */
128  template <typename PointArray>
129  VTKM_EXEC void operator()(const vtkm::Id index,
130  const PointArray& input,
131  const PointArray& output,
132  Scalar& outputField) const
133  {
134  const vtkm::Vec<vtkm::Id, 6> neighborIndices = this->GridData.GetNeighborIndices(index);
135 
136  auto xin1 = input.Get(neighborIndices[0]);
137  auto xin2 = input.Get(neighborIndices[1]);
138  auto yin1 = input.Get(neighborIndices[2]);
139  auto yin2 = input.Get(neighborIndices[3]);
140  auto zin1 = input.Get(neighborIndices[4]);
141  auto zin2 = input.Get(neighborIndices[5]);
142 
143  Scalar xDiff = 1.0f / (xin2[0] - xin1[0]);
144  Scalar yDiff = 1.0f / (yin2[1] - yin1[1]);
145  Scalar zDiff = 1.0f / (zin2[2] - zin1[2]);
146 
147  auto xout1 = output.Get(neighborIndices[0]);
148  auto xout2 = output.Get(neighborIndices[1]);
149  auto yout1 = output.Get(neighborIndices[2]);
150  auto yout2 = output.Get(neighborIndices[3]);
151  auto zout1 = output.Get(neighborIndices[4]);
152  auto zout2 = output.Get(neighborIndices[5]);
153 
154  // Total X gradient w.r.t X, Y, Z
155  Scalar f1x = (xout2[0] - xout1[0]) * xDiff;
156  Scalar f1y = (yout2[0] - yout1[0]) * yDiff;
157  Scalar f1z = (zout2[0] - zout1[0]) * zDiff;
158 
159  // Total Y gradient w.r.t X, Y, Z
160  Scalar f2x = (xout2[1] - xout1[1]) * xDiff;
161  Scalar f2y = (yout2[1] - yout1[1]) * yDiff;
162  Scalar f2z = (zout2[1] - zout1[1]) * zDiff;
163 
164  // Total Z gradient w.r.t X, Y, Z
165  Scalar f3x = (xout2[2] - xout1[2]) * xDiff;
166  Scalar f3y = (yout2[2] - yout1[2]) * yDiff;
167  Scalar f3z = (zout2[2] - zout1[2]) * zDiff;
168 
170  vtkm::MatrixSetRow(jacobian, 0, vtkm::Vec<Scalar, 3>(f1x, f1y, f1z));
171  vtkm::MatrixSetRow(jacobian, 1, vtkm::Vec<Scalar, 3>(f2x, f2y, f2z));
172  vtkm::MatrixSetRow(jacobian, 2, vtkm::Vec<Scalar, 3>(f3x, f3y, f3z));
173 
174  vtkm::filter::flow::internal::ComputeLeftCauchyGreenTensor(jacobian);
175 
176  vtkm::Vec<Scalar, 3> eigenValues;
177  vtkm::filter::flow::internal::Jacobi(jacobian, eigenValues);
178 
179  Scalar delta = eigenValues[0];
180  // Given endTime is in units where start time is 0. else do endTime-startTime
181  // return value for ftle computation
182  outputField = vtkm::Log(delta) / (static_cast<Scalar>(2.0f) * EndTime);
183  }
184 
185 public:
186  // To calculate FTLE field
188  // To assist in calculation of indices
189  vtkm::filter::flow::internal::GridMetaData GridData;
190 };
191 
192 }
193 }
194 } //vtkm::worklet::flow
195 
196 #endif //vtk_m_filter_flow_worklet_LagrangianStructures_h
vtkm::worklet::flow::LagrangianStructures< 3 >::LagrangianStructures
VTKM_CONT LagrangianStructures(Scalar endTime, vtkm::cont::UnknownCellSet cellSet)
Definition: worklet/LagrangianStructures.h:115
vtkm::worklet::flow::LagrangianStructures< 2 >::operator()
VTKM_EXEC void operator()(const vtkm::Id index, const PointArray &input, const PointArray &output, Scalar &outputField) const
Definition: worklet/LagrangianStructures.h:50
VTKM_EXEC
#define VTKM_EXEC
Definition: ExportMacros.h:51
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
LagrangianStructureHelpers.h
Types.h
WorkletMapField.h
vtkm::worklet::WorkletMapField::FieldOut
A control signature tag for output fields.
Definition: WorkletMapField.h:60
vtkm::worklet::flow::LagrangianStructures
Definition: worklet/LagrangianStructures.h:30
vtkm::worklet::flow::LagrangianStructures< 3 >::Scalar
vtkm::FloatDefault Scalar
Definition: worklet/LagrangianStructures.h:112
vtkm::MatrixSetRow
VTKM_EXEC_CONT void MatrixSetRow(vtkm::Matrix< T, NumRow, NumCol > &matrix, vtkm::IdComponent rowIndex, const vtkm::Vec< T, NumCol > &rowValues)
Convenience function for setting a row of a matrix.
Definition: Matrix.h:130
Matrix.h
vtkm::cont::UnknownCellSet
A CellSet of an unknown type.
Definition: UnknownCellSet.h:48
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
DispatcherMapField.h
vtkm::worklet::flow::LagrangianStructures< 3 >::EndTime
Scalar EndTime
Definition: worklet/LagrangianStructures.h:187
vtkm::worklet::flow::LagrangianStructures< 3 >::operator()
VTKM_EXEC void operator()(const vtkm::Id index, const PointArray &input, const PointArray &output, Scalar &outputField) const
Definition: worklet/LagrangianStructures.h:129
vtkm::worklet::flow::LagrangianStructures< 2 >::GridData
vtkm::filter::flow::internal::GridMetaData GridData
Definition: worklet/LagrangianStructures.h:105
vtkm::worklet::flow::LagrangianStructures< 3 >::ControlSignature
void(WholeArrayIn, WholeArrayIn, FieldOut) ControlSignature
Definition: worklet/LagrangianStructures.h:121
vtkm::worklet::flow::LagrangianStructures< 2 >::ExecutionSignature
void(WorkIndex, _1, _2, _3) ExecutionSignature
Definition: worklet/LagrangianStructures.h:47
vtkm::worklet::flow::LagrangianStructures< 2 >::Scalar
vtkm::FloatDefault Scalar
Definition: worklet/LagrangianStructures.h:36
VTKM_CONT
#define VTKM_CONT
Definition: ExportMacros.h:57
vtkm::worklet::flow::LagrangianStructures< 3 >::ExecutionSignature
void(WorkIndex, _1, _2, _3) ExecutionSignature
Definition: worklet/LagrangianStructures.h:123
vtkm::Vec
A short fixed-length array.
Definition: Types.h:767
vtkm::Matrix
Basic Matrix type.
Definition: Matrix.h:33
vtkm::FloatDefault
vtkm::Float32 FloatDefault
The floating point type to use when no other precision is specified.
Definition: Types.h:198
vtkm::worklet::flow::LagrangianStructures< 3 >::GridData
vtkm::filter::flow::internal::GridMetaData GridData
Definition: worklet/LagrangianStructures.h:189
vtkm::worklet::flow::LagrangianStructures< 2 >::EndTime
Scalar EndTime
Definition: worklet/LagrangianStructures.h:103
vtkm::Log
VTKM_EXEC_CONT vtkm::Float32 Log(vtkm::Float32 x)
Computes the natural logarithm of x.
Definition: Math.h:1455
vtkm::worklet::flow::LagrangianStructures< 2 >::ControlSignature
void(WholeArrayIn, WholeArrayIn, FieldOut) ControlSignature
Definition: worklet/LagrangianStructures.h:45
vtkm::worklet::WorkletMapField
Base class for worklets that do a simple mapping of field arrays.
Definition: WorkletMapField.h:38
vtkm::exec::arg::WorkIndex
The ExecutionSignature tag to use to get the work index.
Definition: WorkIndex.h:39
GridMetaData.h
vtkm::worklet::flow::LagrangianStructures< 2 >::LagrangianStructures
VTKM_CONT LagrangianStructures(Scalar endTime, vtkm::cont::UnknownCellSet cellSet)
Definition: worklet/LagrangianStructures.h:39