OpenVDB  7.0.0
Operators.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
5 
6 #ifndef OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED
7 #define OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED
8 
9 #include "FiniteDifference.h"
10 #include "Stencils.h"
11 #include "Maps.h"
12 #include "Transform.h"
13 #include <cmath> // for std::sqrt()
14 
15 
16 namespace openvdb {
18 namespace OPENVDB_VERSION_NAME {
19 namespace math {
20 
21 // Simple tools to help determine when type conversions are needed
22 template<typename Vec3T> struct is_vec3d { static const bool value = false; };
23 template<> struct is_vec3d<Vec3d> { static const bool value = true; };
24 
25 template<typename T> struct is_double { static const bool value = false; };
26 template<> struct is_double<double> { static const bool value = true; };
27 
28 
34 template<typename MapType, typename OpType, typename ResultType>
35 struct MapAdapter {
36  MapAdapter(const MapType& m): map(m) {}
37 
38  template<typename AccessorType>
39  inline ResultType
40  result(const AccessorType& grid, const Coord& ijk) { return OpType::result(map, grid, ijk); }
41 
42  template<typename StencilType>
43  inline ResultType
44  result(const StencilType& stencil) { return OpType::result(map, stencil); }
45 
46  const MapType map;
47 };
48 
49 
51 template<typename OpType>
52 struct ISOpMagnitude {
53  template<typename AccessorType>
54  static inline double result(const AccessorType& grid, const Coord& ijk) {
55  return double(OpType::result(grid, ijk).length());
56  }
57 
58  template<typename StencilType>
59  static inline double result(const StencilType& stencil) {
60  return double(OpType::result(stencil).length());
61  }
62 };
63 
65 template<typename OpType, typename MapT>
66 struct OpMagnitude {
67  template<typename AccessorType>
68  static inline double result(const MapT& map, const AccessorType& grid, const Coord& ijk) {
69  return double(OpType::result(map, grid, ijk).length());
70  }
71 
72  template<typename StencilType>
73  static inline double result(const MapT& map, const StencilType& stencil) {
74  return double(OpType::result(map, stencil).length());
75  }
76 };
77 
78 
79 namespace internal {
80 
81 // This additional layer is necessary for Visual C++ to compile.
82 template<typename T>
83 struct ReturnValue {
84  using ValueType = typename T::ValueType;
86 };
87 
88 } // namespace internal
89 
90 // ---- Operators defined in index space
91 
92 
94 template<DScheme DiffScheme>
96 struct ISGradient
97 {
98  // random access version
99  template<typename Accessor> static Vec3<typename Accessor::ValueType>
100  result(const Accessor& grid, const Coord& ijk)
101  {
102  using ValueType = typename Accessor::ValueType;
103  using Vec3Type = Vec3<ValueType>;
104  return Vec3Type( D1<DiffScheme>::inX(grid, ijk),
105  D1<DiffScheme>::inY(grid, ijk),
106  D1<DiffScheme>::inZ(grid, ijk) );
107  }
108 
109  // stencil access version
110  template<typename StencilT> static Vec3<typename StencilT::ValueType>
111  result(const StencilT& stencil)
112  {
113  using ValueType = typename StencilT::ValueType;
114  using Vec3Type = Vec3<ValueType>;
115  return Vec3Type( D1<DiffScheme>::inX(stencil),
116  D1<DiffScheme>::inY(stencil),
117  D1<DiffScheme>::inZ(stencil) );
118  }
119 };
121 
125 template<BiasedGradientScheme bgs>
126 struct BIAS_SCHEME {
127  static const DScheme FD = FD_1ST;
128  static const DScheme BD = BD_1ST;
129 
130  template<typename GridType, bool IsSafe = true>
131  struct ISStencil {
133  };
134 };
135 
136 template<> struct BIAS_SCHEME<FIRST_BIAS>
137 {
138  static const DScheme FD = FD_1ST;
139  static const DScheme BD = BD_1ST;
140 
141  template<typename GridType, bool IsSafe = true>
142  struct ISStencil {
144  };
145 };
146 
147 template<> struct BIAS_SCHEME<SECOND_BIAS>
148 {
149  static const DScheme FD = FD_2ND;
150  static const DScheme BD = BD_2ND;
151 
152  template<typename GridType, bool IsSafe = true>
153  struct ISStencil {
155  };
156 };
157 template<> struct BIAS_SCHEME<THIRD_BIAS>
158 {
159  static const DScheme FD = FD_3RD;
160  static const DScheme BD = BD_3RD;
161 
162  template<typename GridType, bool IsSafe = true>
163  struct ISStencil {
165  };
166 };
167 template<> struct BIAS_SCHEME<WENO5_BIAS>
168 {
169  static const DScheme FD = FD_WENO5;
170  static const DScheme BD = BD_WENO5;
171 
172  template<typename GridType, bool IsSafe = true>
173  struct ISStencil {
175  };
176 };
177 template<> struct BIAS_SCHEME<HJWENO5_BIAS>
178 {
179  static const DScheme FD = FD_HJWENO5;
180  static const DScheme BD = BD_HJWENO5;
181 
182  template<typename GridType, bool IsSafe = true>
183  struct ISStencil {
185  };
186 };
187 
188 
190 
192 template<BiasedGradientScheme GradScheme, typename Vec3Bias>
194 {
197 
198  // random access version
199  template<typename Accessor>
201  result(const Accessor& grid, const Coord& ijk, const Vec3Bias& V)
202  {
203  using ValueType = typename Accessor::ValueType;
204  using Vec3Type = Vec3<ValueType>;
205 
206  return Vec3Type(V[0]<0 ? D1<FD>::inX(grid,ijk) : D1<BD>::inX(grid,ijk),
207  V[1]<0 ? D1<FD>::inY(grid,ijk) : D1<BD>::inY(grid,ijk),
208  V[2]<0 ? D1<FD>::inZ(grid,ijk) : D1<BD>::inZ(grid,ijk) );
209  }
210 
211  // stencil access version
212  template<typename StencilT>
214  result(const StencilT& stencil, const Vec3Bias& V)
215  {
216  using ValueType = typename StencilT::ValueType;
217  using Vec3Type = Vec3<ValueType>;
218 
219  return Vec3Type(V[0]<0 ? D1<FD>::inX(stencil) : D1<BD>::inX(stencil),
220  V[1]<0 ? D1<FD>::inY(stencil) : D1<BD>::inY(stencil),
221  V[2]<0 ? D1<FD>::inZ(stencil) : D1<BD>::inZ(stencil) );
222  }
223 };
224 
225 
226 template<BiasedGradientScheme GradScheme>
228 {
231 
232 
233  // random access version
234  template<typename Accessor>
235  static typename Accessor::ValueType
236  result(const Accessor& grid, const Coord& ijk)
237  {
238  using ValueType = typename Accessor::ValueType;
239  using Vec3Type = math::Vec3<ValueType>;
240 
241  Vec3Type up = ISGradient<FD>::result(grid, ijk);
242  Vec3Type down = ISGradient<BD>::result(grid, ijk);
243  return math::GodunovsNormSqrd(grid.getValue(ijk)>0, down, up);
244  }
245 
246  // stencil access version
247  template<typename StencilT>
248  static typename StencilT::ValueType
249  result(const StencilT& stencil)
250  {
251  using ValueType = typename StencilT::ValueType;
252  using Vec3Type = math::Vec3<ValueType>;
253 
254  Vec3Type up = ISGradient<FD>::result(stencil);
255  Vec3Type down = ISGradient<BD>::result(stencil);
256  return math::GodunovsNormSqrd(stencil.template getValue<0, 0, 0>()>0, down, up);
257  }
258 };
259 
260 #ifdef DWA_OPENVDB // for SIMD - note will do the computations in float
261 template<>
262 struct ISGradientNormSqrd<HJWENO5_BIAS>
263 {
264  // random access version
265  template<typename Accessor>
266  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk)
267  {
268  struct GetValue
269  {
270  const Accessor& acc;
271  GetValue(const Accessor& acc_): acc(acc_) {}
272  // Return the grid value at ijk converted to simd::Float4::value_type (= float).
273  inline simd::Float4::value_type operator()(const Coord& ijk_) {
274  return static_cast<simd::Float4::value_type>(acc.getValue(ijk_));
275  }
276  }
277  valueAt(grid);
278 
279  // SSE optimized
280  const simd::Float4
281  v1(valueAt(ijk.offsetBy(-2, 0, 0)) - valueAt(ijk.offsetBy(-3, 0, 0)),
282  valueAt(ijk.offsetBy( 0,-2, 0)) - valueAt(ijk.offsetBy( 0,-3, 0)),
283  valueAt(ijk.offsetBy( 0, 0,-2)) - valueAt(ijk.offsetBy( 0, 0,-3)), 0),
284  v2(valueAt(ijk.offsetBy(-1, 0, 0)) - valueAt(ijk.offsetBy(-2, 0, 0)),
285  valueAt(ijk.offsetBy( 0,-1, 0)) - valueAt(ijk.offsetBy( 0,-2, 0)),
286  valueAt(ijk.offsetBy( 0, 0,-1)) - valueAt(ijk.offsetBy( 0, 0,-2)), 0),
287  v3(valueAt(ijk ) - valueAt(ijk.offsetBy(-1, 0, 0)),
288  valueAt(ijk ) - valueAt(ijk.offsetBy( 0,-1, 0)),
289  valueAt(ijk ) - valueAt(ijk.offsetBy( 0, 0,-1)), 0),
290  v4(valueAt(ijk.offsetBy( 1, 0, 0)) - valueAt(ijk ),
291  valueAt(ijk.offsetBy( 0, 1, 0)) - valueAt(ijk ),
292  valueAt(ijk.offsetBy( 0, 0, 1)) - valueAt(ijk ), 0),
293  v5(valueAt(ijk.offsetBy( 2, 0, 0)) - valueAt(ijk.offsetBy( 1, 0, 0)),
294  valueAt(ijk.offsetBy( 0, 2, 0)) - valueAt(ijk.offsetBy( 0, 1, 0)),
295  valueAt(ijk.offsetBy( 0, 0, 2)) - valueAt(ijk.offsetBy( 0, 0, 1)), 0),
296  v6(valueAt(ijk.offsetBy( 3, 0, 0)) - valueAt(ijk.offsetBy( 2, 0, 0)),
297  valueAt(ijk.offsetBy( 0, 3, 0)) - valueAt(ijk.offsetBy( 0, 2, 0)),
298  valueAt(ijk.offsetBy( 0, 0, 3)) - valueAt(ijk.offsetBy( 0, 0, 2)), 0),
299  down = math::WENO5(v1, v2, v3, v4, v5),
300  up = math::WENO5(v6, v5, v4, v3, v2);
301 
302  return math::GodunovsNormSqrd(grid.getValue(ijk)>0, down, up);
303  }
304 
305  // stencil access version
306  template<typename StencilT>
307  static typename StencilT::ValueType result(const StencilT& s)
308  {
309  using F4Val = simd::Float4::value_type;
310 
311  // SSE optimized
312  const simd::Float4
313  v1(F4Val(s.template getValue<-2, 0, 0>()) - F4Val(s.template getValue<-3, 0, 0>()),
314  F4Val(s.template getValue< 0,-2, 0>()) - F4Val(s.template getValue< 0,-3, 0>()),
315  F4Val(s.template getValue< 0, 0,-2>()) - F4Val(s.template getValue< 0, 0,-3>()), 0),
316  v2(F4Val(s.template getValue<-1, 0, 0>()) - F4Val(s.template getValue<-2, 0, 0>()),
317  F4Val(s.template getValue< 0,-1, 0>()) - F4Val(s.template getValue< 0,-2, 0>()),
318  F4Val(s.template getValue< 0, 0,-1>()) - F4Val(s.template getValue< 0, 0,-2>()), 0),
319  v3(F4Val(s.template getValue< 0, 0, 0>()) - F4Val(s.template getValue<-1, 0, 0>()),
320  F4Val(s.template getValue< 0, 0, 0>()) - F4Val(s.template getValue< 0,-1, 0>()),
321  F4Val(s.template getValue< 0, 0, 0>()) - F4Val(s.template getValue< 0, 0,-1>()), 0),
322  v4(F4Val(s.template getValue< 1, 0, 0>()) - F4Val(s.template getValue< 0, 0, 0>()),
323  F4Val(s.template getValue< 0, 1, 0>()) - F4Val(s.template getValue< 0, 0, 0>()),
324  F4Val(s.template getValue< 0, 0, 1>()) - F4Val(s.template getValue< 0, 0, 0>()), 0),
325  v5(F4Val(s.template getValue< 2, 0, 0>()) - F4Val(s.template getValue< 1, 0, 0>()),
326  F4Val(s.template getValue< 0, 2, 0>()) - F4Val(s.template getValue< 0, 1, 0>()),
327  F4Val(s.template getValue< 0, 0, 2>()) - F4Val(s.template getValue< 0, 0, 1>()), 0),
328  v6(F4Val(s.template getValue< 3, 0, 0>()) - F4Val(s.template getValue< 2, 0, 0>()),
329  F4Val(s.template getValue< 0, 3, 0>()) - F4Val(s.template getValue< 0, 2, 0>()),
330  F4Val(s.template getValue< 0, 0, 3>()) - F4Val(s.template getValue< 0, 0, 2>()), 0),
331  down = math::WENO5(v1, v2, v3, v4, v5),
332  up = math::WENO5(v6, v5, v4, v3, v2);
333 
334  return math::GodunovsNormSqrd(s.template getValue<0, 0, 0>()>0, down, up);
335  }
336 };
337 #endif //DWA_OPENVDB // for SIMD - note will do the computations in float
338 
339 
340 
342 template<DDScheme DiffScheme>
345 {
346  // random access version
347  template<typename Accessor>
348  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk);
349 
350  // stencil access version
351  template<typename StencilT>
352  static typename StencilT::ValueType result(const StencilT& stencil);
353 };
354 
355 
356 template<>
358 {
359  // random access version
360  template<typename Accessor>
361  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk)
362  {
363  return grid.getValue(ijk.offsetBy(1,0,0)) + grid.getValue(ijk.offsetBy(-1, 0, 0)) +
364  grid.getValue(ijk.offsetBy(0,1,0)) + grid.getValue(ijk.offsetBy(0, -1, 0)) +
365  grid.getValue(ijk.offsetBy(0,0,1)) + grid.getValue(ijk.offsetBy(0, 0,-1))
366  - 6*grid.getValue(ijk);
367  }
368 
369  // stencil access version
370  template<typename StencilT>
371  static typename StencilT::ValueType result(const StencilT& stencil)
372  {
373  return stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() +
374  stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() +
375  stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>()
376  - 6*stencil.template getValue< 0, 0, 0>();
377  }
378 };
379 
380 template<>
382 {
383  // random access version
384  template<typename Accessor>
385  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk)
386  {
387  using ValueT = typename Accessor::ValueType;
388  return static_cast<ValueT>(
389  (-1./12.)*(
390  grid.getValue(ijk.offsetBy(2,0,0)) + grid.getValue(ijk.offsetBy(-2, 0, 0)) +
391  grid.getValue(ijk.offsetBy(0,2,0)) + grid.getValue(ijk.offsetBy( 0,-2, 0)) +
392  grid.getValue(ijk.offsetBy(0,0,2)) + grid.getValue(ijk.offsetBy( 0, 0,-2)) )
393  + (4./3.)*(
394  grid.getValue(ijk.offsetBy(1,0,0)) + grid.getValue(ijk.offsetBy(-1, 0, 0)) +
395  grid.getValue(ijk.offsetBy(0,1,0)) + grid.getValue(ijk.offsetBy( 0,-1, 0)) +
396  grid.getValue(ijk.offsetBy(0,0,1)) + grid.getValue(ijk.offsetBy( 0, 0,-1)) )
397  - 7.5*grid.getValue(ijk));
398  }
399 
400  // stencil access version
401  template<typename StencilT>
402  static typename StencilT::ValueType result(const StencilT& stencil)
403  {
404  using ValueT = typename StencilT::ValueType;
405  return static_cast<ValueT>(
406  (-1./12.)*(
407  stencil.template getValue< 2, 0, 0>() + stencil.template getValue<-2, 0, 0>() +
408  stencil.template getValue< 0, 2, 0>() + stencil.template getValue< 0,-2, 0>() +
409  stencil.template getValue< 0, 0, 2>() + stencil.template getValue< 0, 0,-2>() )
410  + (4./3.)*(
411  stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() +
412  stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() +
413  stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>() )
414  - 7.5*stencil.template getValue< 0, 0, 0>());
415  }
416 };
417 
418 template<>
420 {
421  // random access version
422  template<typename Accessor>
423  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk)
424  {
425  using ValueT = typename Accessor::ValueType;
426  return static_cast<ValueT>(
427  (1./90.)*(
428  grid.getValue(ijk.offsetBy(3,0,0)) + grid.getValue(ijk.offsetBy(-3, 0, 0)) +
429  grid.getValue(ijk.offsetBy(0,3,0)) + grid.getValue(ijk.offsetBy( 0,-3, 0)) +
430  grid.getValue(ijk.offsetBy(0,0,3)) + grid.getValue(ijk.offsetBy( 0, 0,-3)) )
431  - (3./20.)*(
432  grid.getValue(ijk.offsetBy(2,0,0)) + grid.getValue(ijk.offsetBy(-2, 0, 0)) +
433  grid.getValue(ijk.offsetBy(0,2,0)) + grid.getValue(ijk.offsetBy( 0,-2, 0)) +
434  grid.getValue(ijk.offsetBy(0,0,2)) + grid.getValue(ijk.offsetBy( 0, 0,-2)) )
435  + 1.5 *(
436  grid.getValue(ijk.offsetBy(1,0,0)) + grid.getValue(ijk.offsetBy(-1, 0, 0)) +
437  grid.getValue(ijk.offsetBy(0,1,0)) + grid.getValue(ijk.offsetBy( 0,-1, 0)) +
438  grid.getValue(ijk.offsetBy(0,0,1)) + grid.getValue(ijk.offsetBy( 0, 0,-1)) )
439  - (3*49/18.)*grid.getValue(ijk));
440  }
441 
442  // stencil access version
443  template<typename StencilT>
444  static typename StencilT::ValueType result(const StencilT& stencil)
445  {
446  using ValueT = typename StencilT::ValueType;
447  return static_cast<ValueT>(
448  (1./90.)*(
449  stencil.template getValue< 3, 0, 0>() + stencil.template getValue<-3, 0, 0>() +
450  stencil.template getValue< 0, 3, 0>() + stencil.template getValue< 0,-3, 0>() +
451  stencil.template getValue< 0, 0, 3>() + stencil.template getValue< 0, 0,-3>() )
452  - (3./20.)*(
453  stencil.template getValue< 2, 0, 0>() + stencil.template getValue<-2, 0, 0>() +
454  stencil.template getValue< 0, 2, 0>() + stencil.template getValue< 0,-2, 0>() +
455  stencil.template getValue< 0, 0, 2>() + stencil.template getValue< 0, 0,-2>() )
456  + 1.5 *(
457  stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() +
458  stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() +
459  stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>() )
460  - (3*49/18.)*stencil.template getValue< 0, 0, 0>());
461  }
462 };
464 
465 
467 template<DScheme DiffScheme>
470 {
471  // random access version
472  template<typename Accessor> static typename Accessor::ValueType::value_type
473  result(const Accessor& grid, const Coord& ijk)
474  {
475  return D1Vec<DiffScheme>::inX(grid, ijk, 0) +
476  D1Vec<DiffScheme>::inY(grid, ijk, 1) +
477  D1Vec<DiffScheme>::inZ(grid, ijk, 2);
478  }
479 
480  // stencil access version
481  template<typename StencilT> static typename StencilT::ValueType::value_type
482  result(const StencilT& stencil)
483  {
484  return D1Vec<DiffScheme>::inX(stencil, 0) +
485  D1Vec<DiffScheme>::inY(stencil, 1) +
486  D1Vec<DiffScheme>::inZ(stencil, 2);
487  }
488 };
490 
491 
493 template<DScheme DiffScheme>
495 struct ISCurl
496 {
497  // random access version
498  template<typename Accessor>
499  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk)
500  {
501  using Vec3Type = typename Accessor::ValueType;
502  return Vec3Type( D1Vec<DiffScheme>::inY(grid, ijk, 2) - //dw/dy - dv/dz
503  D1Vec<DiffScheme>::inZ(grid, ijk, 1),
504  D1Vec<DiffScheme>::inZ(grid, ijk, 0) - //du/dz - dw/dx
505  D1Vec<DiffScheme>::inX(grid, ijk, 2),
506  D1Vec<DiffScheme>::inX(grid, ijk, 1) - //dv/dx - du/dy
507  D1Vec<DiffScheme>::inY(grid, ijk, 0) );
508  }
509 
510  // stencil access version
511  template<typename StencilT>
512  static typename StencilT::ValueType result(const StencilT& stencil)
513  {
514  using Vec3Type = typename StencilT::ValueType;
515  return Vec3Type( D1Vec<DiffScheme>::inY(stencil, 2) - //dw/dy - dv/dz
516  D1Vec<DiffScheme>::inZ(stencil, 1),
517  D1Vec<DiffScheme>::inZ(stencil, 0) - //du/dz - dw/dx
518  D1Vec<DiffScheme>::inX(stencil, 2),
519  D1Vec<DiffScheme>::inX(stencil, 1) - //dv/dx - du/dy
520  D1Vec<DiffScheme>::inY(stencil, 0) );
521  }
522 };
524 
525 
527 template<DDScheme DiffScheme2, DScheme DiffScheme1>
530 {
535  template<typename Accessor>
536  static bool result(const Accessor& grid, const Coord& ijk,
537  typename Accessor::ValueType& alpha,
538  typename Accessor::ValueType& beta)
539  {
540  using ValueType = typename Accessor::ValueType;
541 
542  const ValueType Dx = D1<DiffScheme1>::inX(grid, ijk);
543  const ValueType Dy = D1<DiffScheme1>::inY(grid, ijk);
544  const ValueType Dz = D1<DiffScheme1>::inZ(grid, ijk);
545 
546  const ValueType Dx2 = Dx*Dx;
547  const ValueType Dy2 = Dy*Dy;
548  const ValueType Dz2 = Dz*Dz;
549  const ValueType normGrad = Dx2 + Dy2 + Dz2;
550  if (normGrad <= math::Tolerance<ValueType>::value()) {
551  alpha = beta = 0;
552  return false;
553  }
554 
555  const ValueType Dxx = D2<DiffScheme2>::inX(grid, ijk);
556  const ValueType Dyy = D2<DiffScheme2>::inY(grid, ijk);
557  const ValueType Dzz = D2<DiffScheme2>::inZ(grid, ijk);
558 
559  const ValueType Dxy = D2<DiffScheme2>::inXandY(grid, ijk);
560  const ValueType Dyz = D2<DiffScheme2>::inYandZ(grid, ijk);
561  const ValueType Dxz = D2<DiffScheme2>::inXandZ(grid, ijk);
562 
563  // for return
564  alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz));
565  beta = ValueType(std::sqrt(double(normGrad))); // * 1/dx
566  return true;
567  }
568 
573  template<typename StencilT>
574  static bool result(const StencilT& stencil,
575  typename StencilT::ValueType& alpha,
576  typename StencilT::ValueType& beta)
577  {
578  using ValueType = typename StencilT::ValueType;
579  const ValueType Dx = D1<DiffScheme1>::inX(stencil);
580  const ValueType Dy = D1<DiffScheme1>::inY(stencil);
581  const ValueType Dz = D1<DiffScheme1>::inZ(stencil);
582 
583  const ValueType Dx2 = Dx*Dx;
584  const ValueType Dy2 = Dy*Dy;
585  const ValueType Dz2 = Dz*Dz;
586  const ValueType normGrad = Dx2 + Dy2 + Dz2;
587  if (normGrad <= math::Tolerance<ValueType>::value()) {
588  alpha = beta = 0;
589  return false;
590  }
591 
592  const ValueType Dxx = D2<DiffScheme2>::inX(stencil);
593  const ValueType Dyy = D2<DiffScheme2>::inY(stencil);
594  const ValueType Dzz = D2<DiffScheme2>::inZ(stencil);
595 
596  const ValueType Dxy = D2<DiffScheme2>::inXandY(stencil);
597  const ValueType Dyz = D2<DiffScheme2>::inYandZ(stencil);
598  const ValueType Dxz = D2<DiffScheme2>::inXandZ(stencil);
599 
600  // for return
601  alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz));
602  beta = ValueType(std::sqrt(double(normGrad))); // * 1/dx
603  return true;
604  }
605 };
606 
608 
609 // --- Operators defined in the Range of a given map
610 
612 template<typename MapType, DScheme DiffScheme>
616 struct Gradient
617 {
618  // random access version
619  template<typename Accessor>
621  result(const MapType& map, const Accessor& grid, const Coord& ijk)
622  {
623  using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type;
624 
625  Vec3d iGradient( ISGradient<DiffScheme>::result(grid, ijk) );
626  return Vec3Type(map.applyIJT(iGradient, ijk.asVec3d()));
627  }
628 
629  // stencil access version
630  template<typename StencilT>
632  result(const MapType& map, const StencilT& stencil)
633  {
634  using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type;
635 
636  Vec3d iGradient( ISGradient<DiffScheme>::result(stencil) );
637  return Vec3Type(map.applyIJT(iGradient, stencil.getCenterCoord().asVec3d()));
638  }
639 };
640 
641 // Partial template specialization of Gradient
642 // translation, any order
643 template<DScheme DiffScheme>
644 struct Gradient<TranslationMap, DiffScheme>
645 {
646  // random access version
647  template<typename Accessor>
649  result(const TranslationMap&, const Accessor& grid, const Coord& ijk)
650  {
651  return ISGradient<DiffScheme>::result(grid, ijk);
652  }
653 
654  // stencil access version
655  template<typename StencilT>
657  result(const TranslationMap&, const StencilT& stencil)
658  {
659  return ISGradient<DiffScheme>::result(stencil);
660  }
661 };
662 
665 template<>
667 {
668  // random access version
669  template<typename Accessor>
671  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
672  {
673  using ValueType = typename internal::ReturnValue<Accessor>::ValueType;
674  using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type;
675 
676  Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) );
677  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
678  return iGradient * inv2dx;
679  }
680 
681  // stencil access version
682  template<typename StencilT>
684  result(const UniformScaleMap& map, const StencilT& stencil)
685  {
686  using ValueType = typename internal::ReturnValue<StencilT>::ValueType;
687  using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type;
688 
689  Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) );
690  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
691  return iGradient * inv2dx;
692  }
693 };
694 
697 template<>
699 {
700  // random access version
701  template<typename Accessor>
703  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
704  {
705  using ValueType = typename internal::ReturnValue<Accessor>::ValueType;
706  using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type;
707 
708  Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) );
709  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
710  return iGradient * inv2dx;
711  }
712 
713  // stencil access version
714  template<typename StencilT>
716  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
717  {
718  using ValueType = typename internal::ReturnValue<StencilT>::ValueType;
719  using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type;
720 
721  Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) );
722  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
723  return iGradient * inv2dx;
724  }
725 };
726 
729 template<>
731 {
732  // random access version
733  template<typename Accessor>
735  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
736  {
737  using ValueType = typename internal::ReturnValue<Accessor>::ValueType;
738  using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type;
739 
740  Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) );
742  const auto gradient0 = iGradient[0] * map.getInvTwiceScale()[0];
743  const auto gradient1 = iGradient[1] * map.getInvTwiceScale()[1];
744  const auto gradient2 = iGradient[2] * map.getInvTwiceScale()[2];
746  return Vec3Type(ValueType(gradient0),
747  ValueType(gradient1),
748  ValueType(gradient2));
749  }
750 
751  // stencil access version
752  template<typename StencilT>
754  result(const ScaleMap& map, const StencilT& stencil)
755  {
756  using ValueType = typename internal::ReturnValue<StencilT>::ValueType;
757  using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type;
758 
759  Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) );
761  const auto gradient0 = iGradient[0] * map.getInvTwiceScale()[0];
762  const auto gradient1 = iGradient[1] * map.getInvTwiceScale()[1];
763  const auto gradient2 = iGradient[2] * map.getInvTwiceScale()[2];
765  return Vec3Type(ValueType(gradient0),
766  ValueType(gradient1),
767  ValueType(gradient2));
768  }
769 };
770 
773 template<>
775 {
776  // random access version
777  template<typename Accessor>
779  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
780  {
781  using ValueType = typename internal::ReturnValue<Accessor>::ValueType;
782  using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type;
783 
784  Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) );
786  const auto gradient0 = iGradient[0] * map.getInvTwiceScale()[0];
787  const auto gradient1 = iGradient[1] * map.getInvTwiceScale()[1];
788  const auto gradient2 = iGradient[2] * map.getInvTwiceScale()[2];
790  return Vec3Type(ValueType(gradient0),
791  ValueType(gradient1),
792  ValueType(gradient2));
793  }
794 
795  // Stencil access version
796  template<typename StencilT>
798  result(const ScaleTranslateMap& map, const StencilT& stencil)
799  {
800  using ValueType = typename internal::ReturnValue<StencilT>::ValueType;
801  using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type;
802 
803  Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) );
805  const auto gradient0 = iGradient[0] * map.getInvTwiceScale()[0];
806  const auto gradient1 = iGradient[1] * map.getInvTwiceScale()[1];
807  const auto gradient2 = iGradient[2] * map.getInvTwiceScale()[2];
809  return Vec3Type(ValueType(gradient0),
810  ValueType(gradient1),
811  ValueType(gradient2));
812  }
813 };
815 
816 
818 template<typename MapType, BiasedGradientScheme GradScheme>
822 {
823  // random access version
824  template<typename Accessor> static math::Vec3<typename Accessor::ValueType>
825  result(const MapType& map, const Accessor& grid, const Coord& ijk,
827  {
828  using ValueType = typename Accessor::ValueType;
829  using Vec3Type = math::Vec3<ValueType>;
830 
831  Vec3d iGradient( ISGradientBiased<GradScheme, Vec3Type>::result(grid, ijk, V) );
832  return Vec3Type(map.applyIJT(iGradient, ijk.asVec3d()));
833  }
834 
835  // stencil access version
836  template<typename StencilT> static math::Vec3<typename StencilT::ValueType>
837  result(const MapType& map, const StencilT& stencil,
839  {
840  using ValueType = typename StencilT::ValueType;
841  using Vec3Type = math::Vec3<ValueType>;
842 
843  Vec3d iGradient( ISGradientBiased<GradScheme, Vec3Type>::result(stencil, V) );
844  return Vec3Type(map.applyIJT(iGradient, stencil.getCenterCoord().asVec3d()));
845  }
846 };
848 
849 
851 
852 // Computes |Grad[Phi]| using upwinding
853 template<typename MapType, BiasedGradientScheme GradScheme>
855 {
858 
859 
860  // random access version
861  template<typename Accessor>
862  static typename Accessor::ValueType
863  result(const MapType& map, const Accessor& grid, const Coord& ijk)
864  {
865  using ValueType = typename Accessor::ValueType;
866  using Vec3Type = math::Vec3<ValueType>;
867 
868  Vec3Type up = Gradient<MapType, FD>::result(map, grid, ijk);
869  Vec3Type down = Gradient<MapType, BD>::result(map, grid, ijk);
870  return math::GodunovsNormSqrd(grid.getValue(ijk)>0, down, up);
871  }
872 
873  // stencil access version
874  template<typename StencilT>
875  static typename StencilT::ValueType
876  result(const MapType& map, const StencilT& stencil)
877  {
878  using ValueType = typename StencilT::ValueType;
879  using Vec3Type = math::Vec3<ValueType>;
880 
881  Vec3Type up = Gradient<MapType, FD>::result(map, stencil);
882  Vec3Type down = Gradient<MapType, BD>::result(map, stencil);
883  return math::GodunovsNormSqrd(stencil.template getValue<0, 0, 0>()>0, down, up);
884  }
885 };
886 
888 template<BiasedGradientScheme GradScheme>
890 {
891  // random access version
892  template<typename Accessor>
893  static typename Accessor::ValueType
894  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
895  {
896  using ValueType = typename Accessor::ValueType;
897 
898  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
899  return invdxdx * ISGradientNormSqrd<GradScheme>::result(grid, ijk);
900  }
901 
902  // stencil access version
903  template<typename StencilT>
904  static typename StencilT::ValueType
905  result(const UniformScaleMap& map, const StencilT& stencil)
906  {
907  using ValueType = typename StencilT::ValueType;
908 
909  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
910  return invdxdx * ISGradientNormSqrd<GradScheme>::result(stencil);
911  }
912 };
913 
915 template<BiasedGradientScheme GradScheme>
917 {
918  // random access version
919  template<typename Accessor>
920  static typename Accessor::ValueType
921  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
922  {
923  using ValueType = typename Accessor::ValueType;
924 
925  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
926  return invdxdx * ISGradientNormSqrd<GradScheme>::result(grid, ijk);
927  }
928 
929  // stencil access version
930  template<typename StencilT>
931  static typename StencilT::ValueType
932  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
933  {
934  using ValueType = typename StencilT::ValueType;
935 
936  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
937  return invdxdx * ISGradientNormSqrd<GradScheme>::result(stencil);
938  }
939 };
940 
941 
943 template<typename MapType, DScheme DiffScheme>
947 {
948  // random access version
949  template<typename Accessor> static typename Accessor::ValueType::value_type
950  result(const MapType& map, const Accessor& grid, const Coord& ijk)
951  {
952  using ValueType = typename Accessor::ValueType::value_type;
953 
954  ValueType div(0);
955  for (int i=0; i < 3; i++) {
956  Vec3d vec( D1Vec<DiffScheme>::inX(grid, ijk, i),
957  D1Vec<DiffScheme>::inY(grid, ijk, i),
958  D1Vec<DiffScheme>::inZ(grid, ijk, i) );
959  div += ValueType(map.applyIJT(vec, ijk.asVec3d())[i]);
960  }
961  return div;
962  }
963 
964  // stencil access version
965  template<typename StencilT> static typename StencilT::ValueType::value_type
966  result(const MapType& map, const StencilT& stencil)
967  {
968  using ValueType = typename StencilT::ValueType::value_type;
969 
970  ValueType div(0);
971  for (int i=0; i < 3; i++) {
972  Vec3d vec( D1Vec<DiffScheme>::inX(stencil, i),
973  D1Vec<DiffScheme>::inY(stencil, i),
974  D1Vec<DiffScheme>::inZ(stencil, i) );
975  div += ValueType(map.applyIJT(vec, stencil.getCenterCoord().asVec3d())[i]);
976  }
977  return div;
978  }
979 };
980 
983 template<DScheme DiffScheme>
984 struct Divergence<TranslationMap, DiffScheme>
985 {
986  // random access version
987  template<typename Accessor> static typename Accessor::ValueType::value_type
988  result(const TranslationMap&, const Accessor& grid, const Coord& ijk)
989  {
990  using ValueType = typename Accessor::ValueType::value_type;
991 
992  ValueType div(0);
993  div =ISDivergence<DiffScheme>::result(grid, ijk);
994  return div;
995  }
996 
997  // stencil access version
998  template<typename StencilT> static typename StencilT::ValueType::value_type
999  result(const TranslationMap&, const StencilT& stencil)
1000  {
1001  using ValueType = typename StencilT::ValueType::value_type;
1002 
1003  ValueType div(0);
1004  div =ISDivergence<DiffScheme>::result(stencil);
1005  return div;
1006  }
1007 };
1008 
1011 template<DScheme DiffScheme>
1012 struct Divergence<UniformScaleMap, DiffScheme>
1013 {
1014  // random access version
1015  template<typename Accessor> static typename Accessor::ValueType::value_type
1016  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1017  {
1018  using ValueType = typename Accessor::ValueType::value_type;
1019 
1020  ValueType div(0);
1021 
1022  div =ISDivergence<DiffScheme>::result(grid, ijk);
1023  ValueType invdx = ValueType(map.getInvScale()[0]);
1024  return div * invdx;
1025  }
1026 
1027  // stencil access version
1028  template<typename StencilT> static typename StencilT::ValueType::value_type
1029  result(const UniformScaleMap& map, const StencilT& stencil)
1030  {
1031  using ValueType = typename StencilT::ValueType::value_type;
1032 
1033  ValueType div(0);
1034 
1035  div =ISDivergence<DiffScheme>::result(stencil);
1036  ValueType invdx = ValueType(map.getInvScale()[0]);
1037  return div * invdx;
1038  }
1039 };
1040 
1043 template<DScheme DiffScheme>
1045 {
1046  // random access version
1047  template<typename Accessor> static typename Accessor::ValueType::value_type
1048  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1049  {
1050  using ValueType = typename Accessor::ValueType::value_type;
1051 
1052  ValueType div(0);
1053 
1054  div =ISDivergence<DiffScheme>::result(grid, ijk);
1055  ValueType invdx = ValueType(map.getInvScale()[0]);
1056  return div * invdx;
1057  }
1058 
1059  // stencil access version
1060  template<typename StencilT> static typename StencilT::ValueType::value_type
1061  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1062  {
1063  using ValueType = typename StencilT::ValueType::value_type;
1064 
1065  ValueType div(0);
1066 
1067  div =ISDivergence<DiffScheme>::result(stencil);
1068  ValueType invdx = ValueType(map.getInvScale()[0]);
1069  return div * invdx;
1070  }
1071 };
1072 
1075 template<>
1077 {
1078  // random access version
1079  template<typename Accessor> static typename Accessor::ValueType::value_type
1080  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1081  {
1082  using ValueType = typename Accessor::ValueType::value_type;
1083 
1084  ValueType div(0);
1085  div =ISDivergence<CD_2NDT>::result(grid, ijk);
1086  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1087  return div * inv2dx;
1088  }
1089 
1090  // stencil access version
1091  template<typename StencilT> static typename StencilT::ValueType::value_type
1092  result(const UniformScaleMap& map, const StencilT& stencil)
1093  {
1094  using ValueType = typename StencilT::ValueType::value_type;
1095 
1096  ValueType div(0);
1097  div =ISDivergence<CD_2NDT>::result(stencil);
1098  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1099  return div * inv2dx;
1100  }
1101 };
1102 
1105 template<>
1107 {
1108  // random access version
1109  template<typename Accessor> static typename Accessor::ValueType::value_type
1110  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1111  {
1112  using ValueType = typename Accessor::ValueType::value_type;
1113 
1114  ValueType div(0);
1115 
1116  div =ISDivergence<CD_2NDT>::result(grid, ijk);
1117  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1118  return div * inv2dx;
1119  }
1120 
1121  // stencil access version
1122  template<typename StencilT> static typename StencilT::ValueType::value_type
1123  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1124  {
1125  using ValueType = typename StencilT::ValueType::value_type;
1126 
1127  ValueType div(0);
1128 
1129  div =ISDivergence<CD_2NDT>::result(stencil);
1130  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1131  return div * inv2dx;
1132  }
1133 };
1134 
1137 template<DScheme DiffScheme>
1138 struct Divergence<ScaleMap, DiffScheme>
1139 {
1140  // random access version
1141  template<typename Accessor> static typename Accessor::ValueType::value_type
1142  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
1143  {
1144  using ValueType = typename Accessor::ValueType::value_type;
1145 
1146  ValueType div = ValueType(
1147  D1Vec<DiffScheme>::inX(grid, ijk, 0) * (map.getInvScale()[0]) +
1148  D1Vec<DiffScheme>::inY(grid, ijk, 1) * (map.getInvScale()[1]) +
1149  D1Vec<DiffScheme>::inZ(grid, ijk, 2) * (map.getInvScale()[2]));
1150  return div;
1151  }
1152 
1153  // stencil access version
1154  template<typename StencilT> static typename StencilT::ValueType::value_type
1155  result(const ScaleMap& map, const StencilT& stencil)
1156  {
1157  using ValueType = typename StencilT::ValueType::value_type;
1158 
1159  ValueType div(0);
1160  div = ValueType(
1161  D1Vec<DiffScheme>::inX(stencil, 0) * (map.getInvScale()[0]) +
1162  D1Vec<DiffScheme>::inY(stencil, 1) * (map.getInvScale()[1]) +
1163  D1Vec<DiffScheme>::inZ(stencil, 2) * (map.getInvScale()[2]) );
1164  return div;
1165  }
1166 };
1167 
1170 template<DScheme DiffScheme>
1171 struct Divergence<ScaleTranslateMap, DiffScheme>
1172 {
1173  // random access version
1174  template<typename Accessor> static typename Accessor::ValueType::value_type
1175  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1176  {
1177  using ValueType = typename Accessor::ValueType::value_type;
1178 
1179  ValueType div = ValueType(
1180  D1Vec<DiffScheme>::inX(grid, ijk, 0) * (map.getInvScale()[0]) +
1181  D1Vec<DiffScheme>::inY(grid, ijk, 1) * (map.getInvScale()[1]) +
1182  D1Vec<DiffScheme>::inZ(grid, ijk, 2) * (map.getInvScale()[2]));
1183  return div;
1184  }
1185 
1186  // stencil access version
1187  template<typename StencilT> static typename StencilT::ValueType::value_type
1188  result(const ScaleTranslateMap& map, const StencilT& stencil)
1189  {
1190  using ValueType = typename StencilT::ValueType::value_type;
1191 
1192  ValueType div(0);
1193  div = ValueType(
1194  D1Vec<DiffScheme>::inX(stencil, 0) * (map.getInvScale()[0]) +
1195  D1Vec<DiffScheme>::inY(stencil, 1) * (map.getInvScale()[1]) +
1196  D1Vec<DiffScheme>::inZ(stencil, 2) * (map.getInvScale()[2]) );
1197  return div;
1198  }
1199 };
1200 
1203 template<>
1205 {
1206  // random access version
1207  template<typename Accessor> static typename Accessor::ValueType::value_type
1208  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
1209  {
1210  using ValueType = typename Accessor::ValueType::value_type;
1211 
1212  ValueType div = ValueType(
1213  D1Vec<CD_2NDT>::inX(grid, ijk, 0) * (map.getInvTwiceScale()[0]) +
1214  D1Vec<CD_2NDT>::inY(grid, ijk, 1) * (map.getInvTwiceScale()[1]) +
1215  D1Vec<CD_2NDT>::inZ(grid, ijk, 2) * (map.getInvTwiceScale()[2]) );
1216  return div;
1217  }
1218 
1219  // stencil access version
1220  template<typename StencilT> static typename StencilT::ValueType::value_type
1221  result(const ScaleMap& map, const StencilT& stencil)
1222  {
1223  using ValueType = typename StencilT::ValueType::value_type;
1224 
1225  ValueType div = ValueType(
1226  D1Vec<CD_2NDT>::inX(stencil, 0) * (map.getInvTwiceScale()[0]) +
1227  D1Vec<CD_2NDT>::inY(stencil, 1) * (map.getInvTwiceScale()[1]) +
1228  D1Vec<CD_2NDT>::inZ(stencil, 2) * (map.getInvTwiceScale()[2]) );
1229  return div;
1230  }
1231 };
1232 
1235 template<>
1237 {
1238  // random access version
1239  template<typename Accessor> static typename Accessor::ValueType::value_type
1240  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1241  {
1242  using ValueType = typename Accessor::ValueType::value_type;
1243 
1244  ValueType div = ValueType(
1245  D1Vec<CD_2NDT>::inX(grid, ijk, 0) * (map.getInvTwiceScale()[0]) +
1246  D1Vec<CD_2NDT>::inY(grid, ijk, 1) * (map.getInvTwiceScale()[1]) +
1247  D1Vec<CD_2NDT>::inZ(grid, ijk, 2) * (map.getInvTwiceScale()[2]) );
1248  return div;
1249  }
1250 
1251  // stencil access version
1252  template<typename StencilT> static typename StencilT::ValueType::value_type
1253  result(const ScaleTranslateMap& map, const StencilT& stencil)
1254  {
1255  using ValueType = typename StencilT::ValueType::value_type;
1256 
1257  ValueType div = ValueType(
1258  D1Vec<CD_2NDT>::inX(stencil, 0) * (map.getInvTwiceScale()[0]) +
1259  D1Vec<CD_2NDT>::inY(stencil, 1) * (map.getInvTwiceScale()[1]) +
1260  D1Vec<CD_2NDT>::inZ(stencil, 2) * (map.getInvTwiceScale()[2]) );
1261  return div;
1262  }
1263 };
1265 
1266 
1268 template<typename MapType, DScheme DiffScheme>
1271 struct Curl
1272 {
1273  // random access version
1274  template<typename Accessor> static typename Accessor::ValueType
1275  result(const MapType& map, const Accessor& grid, const Coord& ijk)
1276  {
1277  using Vec3Type = typename Accessor::ValueType;
1278  Vec3Type mat[3];
1279  for (int i = 0; i < 3; i++) {
1280  Vec3d vec(
1281  D1Vec<DiffScheme>::inX(grid, ijk, i),
1282  D1Vec<DiffScheme>::inY(grid, ijk, i),
1283  D1Vec<DiffScheme>::inZ(grid, ijk, i));
1284  // dF_i/dx_j (x_1 = x, x_2 = y, x_3 = z)
1285  mat[i] = Vec3Type(map.applyIJT(vec, ijk.asVec3d()));
1286  }
1287  return Vec3Type(mat[2][1] - mat[1][2], // dF_3/dx_2 - dF_2/dx_3
1288  mat[0][2] - mat[2][0], // dF_1/dx_3 - dF_3/dx_1
1289  mat[1][0] - mat[0][1]); // dF_2/dx_1 - dF_1/dx_2
1290  }
1291 
1292  // stencil access version
1293  template<typename StencilT> static typename StencilT::ValueType
1294  result(const MapType& map, const StencilT& stencil)
1295  {
1296  using Vec3Type = typename StencilT::ValueType;
1297  Vec3Type mat[3];
1298  for (int i = 0; i < 3; i++) {
1299  Vec3d vec(
1300  D1Vec<DiffScheme>::inX(stencil, i),
1301  D1Vec<DiffScheme>::inY(stencil, i),
1302  D1Vec<DiffScheme>::inZ(stencil, i));
1303  // dF_i/dx_j (x_1 = x, x_2 = y, x_3 = z)
1304  mat[i] = Vec3Type(map.applyIJT(vec, stencil.getCenterCoord().asVec3d()));
1305  }
1306  return Vec3Type(mat[2][1] - mat[1][2], // dF_3/dx_2 - dF_2/dx_3
1307  mat[0][2] - mat[2][0], // dF_1/dx_3 - dF_3/dx_1
1308  mat[1][0] - mat[0][1]); // dF_2/dx_1 - dF_1/dx_2
1309  }
1310 };
1311 
1313 template<DScheme DiffScheme>
1314 struct Curl<UniformScaleMap, DiffScheme>
1315 {
1316  // random access version
1317  template<typename Accessor> static typename Accessor::ValueType
1318  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1319  {
1320  using Vec3Type = typename Accessor::ValueType;
1321  using ValueType = typename Vec3Type::value_type;
1322  return ISCurl<DiffScheme>::result(grid, ijk) * ValueType(map.getInvScale()[0]);
1323  }
1324 
1325  // Stencil access version
1326  template<typename StencilT> static typename StencilT::ValueType
1327  result(const UniformScaleMap& map, const StencilT& stencil)
1328  {
1329  using Vec3Type = typename StencilT::ValueType;
1330  using ValueType = typename Vec3Type::value_type;
1331  return ISCurl<DiffScheme>::result(stencil) * ValueType(map.getInvScale()[0]);
1332  }
1333 };
1334 
1336 template<DScheme DiffScheme>
1337 struct Curl<UniformScaleTranslateMap, DiffScheme>
1338 {
1339  // random access version
1340  template<typename Accessor> static typename Accessor::ValueType
1341  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1342  {
1343  using Vec3Type = typename Accessor::ValueType;
1344  using ValueType = typename Vec3Type::value_type;
1345 
1346  return ISCurl<DiffScheme>::result(grid, ijk) * ValueType(map.getInvScale()[0]);
1347  }
1348 
1349  // stencil access version
1350  template<typename StencilT> static typename StencilT::ValueType
1351  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1352  {
1353  using Vec3Type = typename StencilT::ValueType;
1354  using ValueType = typename Vec3Type::value_type;
1355 
1356  return ISCurl<DiffScheme>::result(stencil) * ValueType(map.getInvScale()[0]);
1357  }
1358 };
1359 
1361 template<>
1363 {
1364  // random access version
1365  template<typename Accessor> static typename Accessor::ValueType
1366  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1367  {
1368  using Vec3Type = typename Accessor::ValueType;
1369  using ValueType = typename Vec3Type::value_type;
1370 
1371  return ISCurl<CD_2NDT>::result(grid, ijk) * ValueType(map.getInvTwiceScale()[0]);
1372  }
1373 
1374  // stencil access version
1375  template<typename StencilT> static typename StencilT::ValueType
1376  result(const UniformScaleMap& map, const StencilT& stencil)
1377  {
1378  using Vec3Type = typename StencilT::ValueType;
1379  using ValueType = typename Vec3Type::value_type;
1380 
1381  return ISCurl<CD_2NDT>::result(stencil) * ValueType(map.getInvTwiceScale()[0]);
1382  }
1383 };
1384 
1386 template<>
1388 {
1389  // random access version
1390  template<typename Accessor> static typename Accessor::ValueType
1391  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1392  {
1393  using Vec3Type = typename Accessor::ValueType;
1394  using ValueType = typename Vec3Type::value_type;
1395 
1396  return ISCurl<CD_2NDT>::result(grid, ijk) * ValueType(map.getInvTwiceScale()[0]);
1397  }
1398 
1399  // stencil access version
1400  template<typename StencilT> static typename StencilT::ValueType
1401  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1402  {
1403  using Vec3Type = typename StencilT::ValueType;
1404  using ValueType = typename Vec3Type::value_type;
1405 
1406  return ISCurl<CD_2NDT>::result(stencil) * ValueType(map.getInvTwiceScale()[0]);
1407  }
1408 };
1410 
1411 
1413 template<typename MapType, DDScheme DiffScheme>
1417 {
1418  // random access version
1419  template<typename Accessor>
1420  static typename Accessor::ValueType result(const MapType& map,
1421  const Accessor& grid, const Coord& ijk)
1422  {
1423  using ValueType = typename Accessor::ValueType;
1424  // all the second derivatives in index space
1425  ValueType iddx = D2<DiffScheme>::inX(grid, ijk);
1426  ValueType iddy = D2<DiffScheme>::inY(grid, ijk);
1427  ValueType iddz = D2<DiffScheme>::inZ(grid, ijk);
1428 
1429  ValueType iddxy = D2<DiffScheme>::inXandY(grid, ijk);
1430  ValueType iddyz = D2<DiffScheme>::inYandZ(grid, ijk);
1431  ValueType iddxz = D2<DiffScheme>::inXandZ(grid, ijk);
1432 
1433  // second derivatives in index space
1434  Mat3d d2_is(iddx, iddxy, iddxz,
1435  iddxy, iddy, iddyz,
1436  iddxz, iddyz, iddz);
1437 
1438  Mat3d d2_rs; // to hold the second derivative matrix in range space
1440  d2_rs = map.applyIJC(d2_is);
1441  } else {
1442  // compute the first derivatives with 2nd order accuracy.
1443  Vec3d d1_is(static_cast<double>(D1<CD_2ND>::inX(grid, ijk)),
1444  static_cast<double>(D1<CD_2ND>::inY(grid, ijk)),
1445  static_cast<double>(D1<CD_2ND>::inZ(grid, ijk)));
1446 
1447  d2_rs = map.applyIJC(d2_is, d1_is, ijk.asVec3d());
1448  }
1449 
1450  // the trace of the second derivative (range space) matrix is laplacian
1451  return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2));
1452  }
1453 
1454  // stencil access version
1455  template<typename StencilT>
1456  static typename StencilT::ValueType result(const MapType& map, const StencilT& stencil)
1457  {
1458  using ValueType = typename StencilT::ValueType;
1459  // all the second derivatives in index space
1460  ValueType iddx = D2<DiffScheme>::inX(stencil);
1461  ValueType iddy = D2<DiffScheme>::inY(stencil);
1462  ValueType iddz = D2<DiffScheme>::inZ(stencil);
1463 
1464  ValueType iddxy = D2<DiffScheme>::inXandY(stencil);
1465  ValueType iddyz = D2<DiffScheme>::inYandZ(stencil);
1466  ValueType iddxz = D2<DiffScheme>::inXandZ(stencil);
1467 
1468  // second derivatives in index space
1469  Mat3d d2_is(iddx, iddxy, iddxz,
1470  iddxy, iddy, iddyz,
1471  iddxz, iddyz, iddz);
1472 
1473  Mat3d d2_rs; // to hold the second derivative matrix in range space
1475  d2_rs = map.applyIJC(d2_is);
1476  } else {
1477  // compute the first derivatives with 2nd order accuracy.
1478  Vec3d d1_is(D1<CD_2ND>::inX(stencil),
1479  D1<CD_2ND>::inY(stencil),
1480  D1<CD_2ND>::inZ(stencil) );
1481 
1482  d2_rs = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d());
1483  }
1484 
1485  // the trace of the second derivative (range space) matrix is laplacian
1486  return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2));
1487  }
1488 };
1489 
1490 
1491 template<DDScheme DiffScheme>
1492 struct Laplacian<TranslationMap, DiffScheme>
1493 {
1494  // random access version
1495  template<typename Accessor>
1496  static typename Accessor::ValueType result(const TranslationMap&,
1497  const Accessor& grid, const Coord& ijk)
1498  {
1499  return ISLaplacian<DiffScheme>::result(grid, ijk);
1500  }
1501 
1502  // stencil access version
1503  template<typename StencilT>
1504  static typename StencilT::ValueType result(const TranslationMap&, const StencilT& stencil)
1505  {
1506  return ISLaplacian<DiffScheme>::result(stencil);
1507  }
1508 };
1509 
1510 
1511 // The Laplacian is invariant to rotation or reflection.
1512 template<DDScheme DiffScheme>
1513 struct Laplacian<UnitaryMap, DiffScheme>
1514 {
1515  // random access version
1516  template<typename Accessor>
1517  static typename Accessor::ValueType result(const UnitaryMap&,
1518  const Accessor& grid, const Coord& ijk)
1519  {
1520  return ISLaplacian<DiffScheme>::result(grid, ijk);
1521  }
1522 
1523  // stencil access version
1524  template<typename StencilT>
1525  static typename StencilT::ValueType result(const UnitaryMap&, const StencilT& stencil)
1526  {
1527  return ISLaplacian<DiffScheme>::result(stencil);
1528  }
1529 };
1530 
1531 
1532 template<DDScheme DiffScheme>
1533 struct Laplacian<UniformScaleMap, DiffScheme>
1534 {
1535  // random access version
1536  template<typename Accessor> static typename Accessor::ValueType
1537  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1538  {
1539  using ValueType = typename Accessor::ValueType;
1540  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1541  return ISLaplacian<DiffScheme>::result(grid, ijk) * invdxdx;
1542  }
1543 
1544  // stencil access version
1545  template<typename StencilT> static typename StencilT::ValueType
1546  result(const UniformScaleMap& map, const StencilT& stencil)
1547  {
1548  using ValueType = typename StencilT::ValueType;
1549  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1550  return ISLaplacian<DiffScheme>::result(stencil) * invdxdx;
1551  }
1552 };
1553 
1554 
1555 template<DDScheme DiffScheme>
1557 {
1558  // random access version
1559  template<typename Accessor> static typename Accessor::ValueType
1560  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1561  {
1562  using ValueType = typename Accessor::ValueType;
1563  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1564  return ISLaplacian<DiffScheme>::result(grid, ijk) * invdxdx;
1565  }
1566 
1567  // stencil access version
1568  template<typename StencilT> static typename StencilT::ValueType
1569  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1570  {
1571  using ValueType = typename StencilT::ValueType;
1572  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1573  return ISLaplacian<DiffScheme>::result(stencil) * invdxdx;
1574  }
1575 };
1576 
1577 
1578 template<DDScheme DiffScheme>
1579 struct Laplacian<ScaleMap, DiffScheme>
1580 {
1581  // random access version
1582  template<typename Accessor> static typename Accessor::ValueType
1583  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
1584  {
1585  using ValueType = typename Accessor::ValueType;
1586 
1587  // compute the second derivatives in index space
1588  ValueType iddx = D2<DiffScheme>::inX(grid, ijk);
1589  ValueType iddy = D2<DiffScheme>::inY(grid, ijk);
1590  ValueType iddz = D2<DiffScheme>::inZ(grid, ijk);
1591  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1593  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1594  const ValueType value = iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2];
1596  return value;
1597  }
1598 
1599  // stencil access version
1600  template<typename StencilT> static typename StencilT::ValueType
1601  result(const ScaleMap& map, const StencilT& stencil)
1602  {
1603  using ValueType = typename StencilT::ValueType;
1604 
1605  // compute the second derivatives in index space
1606  ValueType iddx = D2<DiffScheme>::inX(stencil);
1607  ValueType iddy = D2<DiffScheme>::inY(stencil);
1608  ValueType iddz = D2<DiffScheme>::inZ(stencil);
1609  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1611  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1612  const ValueType value = iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2];
1614  return value;
1615  }
1616 };
1617 
1618 
1619 template<DDScheme DiffScheme>
1620 struct Laplacian<ScaleTranslateMap, DiffScheme>
1621 {
1622  // random access version
1623  template<typename Accessor> static typename Accessor::ValueType
1624  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1625  {
1626  using ValueType = typename Accessor::ValueType;
1627  // compute the second derivatives in index space
1628  ValueType iddx = D2<DiffScheme>::inX(grid, ijk);
1629  ValueType iddy = D2<DiffScheme>::inY(grid, ijk);
1630  ValueType iddz = D2<DiffScheme>::inZ(grid, ijk);
1631  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1633  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1634  const ValueType value = iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2];
1636  return value;
1637  }
1638 
1639  // stencil access version
1640  template<typename StencilT> static typename StencilT::ValueType
1641  result(const ScaleTranslateMap& map, const StencilT& stencil)
1642  {
1643  using ValueType = typename StencilT::ValueType;
1644  // compute the second derivatives in index space
1645  ValueType iddx = D2<DiffScheme>::inX(stencil);
1646  ValueType iddy = D2<DiffScheme>::inY(stencil);
1647  ValueType iddz = D2<DiffScheme>::inZ(stencil);
1648  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1650  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1651  const ValueType value = iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2];
1653  return value;
1654  }
1655 };
1656 
1657 
1661 template<typename MapType, DScheme DiffScheme>
1662 struct CPT
1663 {
1664  // random access version
1665  template<typename Accessor> static math::Vec3<typename Accessor::ValueType>
1666  result(const MapType& map, const Accessor& grid, const Coord& ijk)
1667  {
1668  using ValueType = typename Accessor::ValueType;
1669  using Vec3Type = Vec3<ValueType>;
1670 
1671  // current distance
1672  ValueType d = grid.getValue(ijk);
1673  // compute gradient in physical space where it is a unit normal
1674  // since the grid holds a distance level set.
1675  Vec3d vectorFromSurface(d*Gradient<MapType,DiffScheme>::result(map, grid, ijk));
1677  Vec3d result = ijk.asVec3d() - map.applyInverseMap(vectorFromSurface);
1678  return Vec3Type(result);
1679  } else {
1680  Vec3d location = map.applyMap(ijk.asVec3d());
1681  Vec3d result = map.applyInverseMap(location - vectorFromSurface);
1682  return Vec3Type(result);
1683  }
1684  }
1685 
1686  // stencil access version
1687  template<typename StencilT> static math::Vec3<typename StencilT::ValueType>
1688  result(const MapType& map, const StencilT& stencil)
1689  {
1690  using ValueType = typename StencilT::ValueType;
1691  using Vec3Type = Vec3<ValueType>;
1692 
1693  // current distance
1694  ValueType d = stencil.template getValue<0, 0, 0>();
1695  // compute gradient in physical space where it is a unit normal
1696  // since the grid holds a distance level set.
1697  Vec3d vectorFromSurface(d*Gradient<MapType, DiffScheme>::result(map, stencil));
1699  Vec3d result = stencil.getCenterCoord().asVec3d()
1700  - map.applyInverseMap(vectorFromSurface);
1701  return Vec3Type(result);
1702  } else {
1703  Vec3d location = map.applyMap(stencil.getCenterCoord().asVec3d());
1704  Vec3d result = map.applyInverseMap(location - vectorFromSurface);
1705  return Vec3Type(result);
1706  }
1707  }
1708 };
1709 
1710 
1714 template<typename MapType, DScheme DiffScheme>
1716 {
1717  // random access version
1718  template<typename Accessor> static Vec3<typename Accessor::ValueType>
1719  result(const MapType& map, const Accessor& grid, const Coord& ijk)
1720  {
1721  using ValueType = typename Accessor::ValueType;
1722  using Vec3Type = Vec3<ValueType>;
1723  // current distance
1724  ValueType d = grid.getValue(ijk);
1725  // compute gradient in physical space where it is a unit normal
1726  // since the grid holds a distance level set.
1727  Vec3Type vectorFromSurface =
1728  d*Gradient<MapType,DiffScheme>::result(map, grid, ijk);
1729  Vec3d result = map.applyMap(ijk.asVec3d()) - vectorFromSurface;
1730 
1731  return Vec3Type(result);
1732  }
1733 
1734  // stencil access version
1735  template<typename StencilT> static Vec3<typename StencilT::ValueType>
1736  result(const MapType& map, const StencilT& stencil)
1737  {
1738  using ValueType = typename StencilT::ValueType;
1739  using Vec3Type = Vec3<ValueType>;
1740  // current distance
1741  ValueType d = stencil.template getValue<0, 0, 0>();
1742  // compute gradient in physical space where it is a unit normal
1743  // since the grid holds a distance level set.
1744  Vec3Type vectorFromSurface =
1745  d*Gradient<MapType, DiffScheme>::result(map, stencil);
1746  Vec3d result = map.applyMap(stencil.getCenterCoord().asVec3d()) - vectorFromSurface;
1747 
1748  return Vec3Type(result);
1749  }
1750 };
1751 
1752 
1757 template<typename MapType, DDScheme DiffScheme2, DScheme DiffScheme1>
1759 {
1764  template<typename Accessor>
1765  static bool compute(const MapType& map, const Accessor& grid, const Coord& ijk,
1766  double& alpha, double& beta)
1767  {
1768  using ValueType = typename Accessor::ValueType;
1769 
1770  // compute the gradient in index and world space
1771  Vec3d d1_is(static_cast<double>(D1<DiffScheme1>::inX(grid, ijk)),
1772  static_cast<double>(D1<DiffScheme1>::inY(grid, ijk)),
1773  static_cast<double>(D1<DiffScheme1>::inZ(grid, ijk))), d1_ws;
1774  if (is_linear<MapType>::value) {//resolved at compiletime
1775  d1_ws = map.applyIJT(d1_is);
1776  } else {
1777  d1_ws = map.applyIJT(d1_is, ijk.asVec3d());
1778  }
1779  const double Dx2 = d1_ws(0)*d1_ws(0);
1780  const double Dy2 = d1_ws(1)*d1_ws(1);
1781  const double Dz2 = d1_ws(2)*d1_ws(2);
1782  const double normGrad = Dx2 + Dy2 + Dz2;
1783  if (normGrad <= math::Tolerance<double>::value()) {
1784  alpha = beta = 0;
1785  return false;
1786  }
1787 
1788  // all the second derivatives in index space
1789  ValueType iddx = D2<DiffScheme2>::inX(grid, ijk);
1790  ValueType iddy = D2<DiffScheme2>::inY(grid, ijk);
1791  ValueType iddz = D2<DiffScheme2>::inZ(grid, ijk);
1792 
1793  ValueType iddxy = D2<DiffScheme2>::inXandY(grid, ijk);
1794  ValueType iddyz = D2<DiffScheme2>::inYandZ(grid, ijk);
1795  ValueType iddxz = D2<DiffScheme2>::inXandZ(grid, ijk);
1796 
1797  // second derivatives in index space
1798  Mat3d d2_is(iddx, iddxy, iddxz,
1799  iddxy, iddy, iddyz,
1800  iddxz, iddyz, iddz);
1801 
1802  // convert second derivatives to world space
1803  Mat3d d2_ws;
1804  if (is_linear<MapType>::value) {//resolved at compiletime
1805  d2_ws = map.applyIJC(d2_is);
1806  } else {
1807  d2_ws = map.applyIJC(d2_is, d1_is, ijk.asVec3d());
1808  }
1809 
1810  // assemble the nominator and denominator for mean curvature
1811  alpha = (Dx2*(d2_ws(1,1)+d2_ws(2,2))+Dy2*(d2_ws(0,0)+d2_ws(2,2))
1812  +Dz2*(d2_ws(0,0)+d2_ws(1,1))
1813  -2*(d1_ws(0)*(d1_ws(1)*d2_ws(0,1)+d1_ws(2)*d2_ws(0,2))
1814  +d1_ws(1)*d1_ws(2)*d2_ws(1,2)));
1815  beta = std::sqrt(normGrad); // * 1/dx
1816  return true;
1817  }
1818 
1819  template<typename Accessor>
1820  static typename Accessor::ValueType result(const MapType& map,
1821  const Accessor& grid, const Coord& ijk)
1822  {
1823  using ValueType = typename Accessor::ValueType;
1824  double alpha, beta;
1825  return compute(map, grid, ijk, alpha, beta) ?
1826  ValueType(alpha/(2. *math::Pow3(beta))) : 0;
1827  }
1828 
1829  template<typename Accessor>
1830  static typename Accessor::ValueType normGrad(const MapType& map,
1831  const Accessor& grid, const Coord& ijk)
1832  {
1833  using ValueType = typename Accessor::ValueType;
1834  double alpha, beta;
1835  return compute(map, grid, ijk, alpha, beta) ?
1836  ValueType(alpha/(2. *math::Pow2(beta))) : 0;
1837  }
1838 
1843  template<typename StencilT>
1844  static bool compute(const MapType& map, const StencilT& stencil,
1845  double& alpha, double& beta)
1846  {
1847  using ValueType = typename StencilT::ValueType;
1848 
1849  // compute the gradient in index and world space
1850  Vec3d d1_is(D1<DiffScheme1>::inX(stencil),
1851  D1<DiffScheme1>::inY(stencil),
1852  D1<DiffScheme1>::inZ(stencil) ), d1_ws;
1853  if (is_linear<MapType>::value) {//resolved at compiletime
1854  d1_ws = map.applyIJT(d1_is);
1855  } else {
1856  d1_ws = map.applyIJT(d1_is, stencil.getCenterCoord().asVec3d());
1857  }
1858  const double Dx2 = d1_ws(0)*d1_ws(0);
1859  const double Dy2 = d1_ws(1)*d1_ws(1);
1860  const double Dz2 = d1_ws(2)*d1_ws(2);
1861  const double normGrad = Dx2 + Dy2 + Dz2;
1862  if (normGrad <= math::Tolerance<double>::value()) {
1863  alpha = beta = 0;
1864  return false;
1865  }
1866 
1867  // all the second derivatives in index space
1868  ValueType iddx = D2<DiffScheme2>::inX(stencil);
1869  ValueType iddy = D2<DiffScheme2>::inY(stencil);
1870  ValueType iddz = D2<DiffScheme2>::inZ(stencil);
1871 
1872  ValueType iddxy = D2<DiffScheme2>::inXandY(stencil);
1873  ValueType iddyz = D2<DiffScheme2>::inYandZ(stencil);
1874  ValueType iddxz = D2<DiffScheme2>::inXandZ(stencil);
1875 
1876  // second derivatives in index space
1877  Mat3d d2_is(iddx, iddxy, iddxz,
1878  iddxy, iddy, iddyz,
1879  iddxz, iddyz, iddz);
1880 
1881  // convert second derivatives to world space
1882  Mat3d d2_ws;
1883  if (is_linear<MapType>::value) {//resolved at compiletime
1884  d2_ws = map.applyIJC(d2_is);
1885  } else {
1886  d2_ws = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d());
1887  }
1888 
1889  // for return
1890  alpha = (Dx2*(d2_ws(1,1)+d2_ws(2,2))+Dy2*(d2_ws(0,0)+d2_ws(2,2))
1891  +Dz2*(d2_ws(0,0)+d2_ws(1,1))
1892  -2*(d1_ws(0)*(d1_ws(1)*d2_ws(0,1)+d1_ws(2)*d2_ws(0,2))
1893  +d1_ws(1)*d1_ws(2)*d2_ws(1,2)));
1894  beta = std::sqrt(normGrad); // * 1/dx
1895  return true;
1896  }
1897 
1898  template<typename StencilT>
1899  static typename StencilT::ValueType
1900  result(const MapType& map, const StencilT stencil)
1901  {
1902  using ValueType = typename StencilT::ValueType;
1903  double alpha, beta;
1904  return compute(map, stencil, alpha, beta) ?
1905  ValueType(alpha/(2*math::Pow3(beta))) : 0;
1906  }
1907 
1908  template<typename StencilT>
1909  static typename StencilT::ValueType normGrad(const MapType& map, const StencilT stencil)
1910  {
1911  using ValueType = typename StencilT::ValueType;
1912  double alpha, beta;
1913  return compute(map, stencil, alpha, beta) ?
1914  ValueType(alpha/(2*math::Pow2(beta))) : 0;
1915  }
1916 };
1917 
1918 
1919 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1920 struct MeanCurvature<TranslationMap, DiffScheme2, DiffScheme1>
1921 {
1922  // random access version
1923  template<typename Accessor>
1924  static typename Accessor::ValueType result(const TranslationMap&,
1925  const Accessor& grid, const Coord& ijk)
1926  {
1927  using ValueType = typename Accessor::ValueType;
1928 
1929  ValueType alpha, beta;
1930  return ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta) ?
1931  ValueType(alpha /(2*math::Pow3(beta))) : 0;
1932  }
1933 
1934  template<typename Accessor>
1935  static typename Accessor::ValueType normGrad(const TranslationMap&,
1936  const Accessor& grid, const Coord& ijk)
1937  {
1938  using ValueType = typename Accessor::ValueType;
1939 
1940  ValueType alpha, beta;
1941  return ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta) ?
1942  ValueType(alpha/(2*math::Pow2(beta))) : 0;
1943  }
1944 
1945  // stencil access version
1946  template<typename StencilT>
1947  static typename StencilT::ValueType result(const TranslationMap&, const StencilT& stencil)
1948  {
1949  using ValueType = typename StencilT::ValueType;
1950 
1951  ValueType alpha, beta;
1952  return ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta) ?
1953  ValueType(alpha /(2*math::Pow3(beta))) : 0;
1954  }
1955 
1956  template<typename StencilT>
1957  static typename StencilT::ValueType normGrad(const TranslationMap&, const StencilT& stencil)
1958  {
1959  using ValueType = typename StencilT::ValueType;
1960 
1961  ValueType alpha, beta;
1962  return ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta) ?
1963  ValueType(alpha/(2*math::Pow2(beta))) : 0;
1964  }
1965 };
1966 
1967 
1968 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1969 struct MeanCurvature<UniformScaleMap, DiffScheme2, DiffScheme1>
1970 {
1971  // random access version
1972  template<typename Accessor>
1973  static typename Accessor::ValueType result(const UniformScaleMap& map,
1974  const Accessor& grid, const Coord& ijk)
1975  {
1976  using ValueType = typename Accessor::ValueType;
1977 
1978  ValueType alpha, beta;
1979  if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta)) {
1980  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1981  return ValueType(alpha*inv2dx/math::Pow3(beta));
1982  }
1983  return 0;
1984  }
1985 
1986  template<typename Accessor>
1987  static typename Accessor::ValueType normGrad(const UniformScaleMap& map,
1988  const Accessor& grid, const Coord& ijk)
1989  {
1990  using ValueType = typename Accessor::ValueType;
1991 
1992  ValueType alpha, beta;
1993  if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta)) {
1994  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1995  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
1996  }
1997  return 0;
1998  }
1999 
2000  // stencil access version
2001  template<typename StencilT>
2002  static typename StencilT::ValueType result(const UniformScaleMap& map, const StencilT& stencil)
2003  {
2004  using ValueType = typename StencilT::ValueType;
2005 
2006  ValueType alpha, beta;
2007  if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta)) {
2008  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
2009  return ValueType(alpha*inv2dx/math::Pow3(beta));
2010  }
2011  return 0;
2012  }
2013 
2014  template<typename StencilT>
2015  static typename StencilT::ValueType normGrad(const UniformScaleMap& map, const StencilT& stencil)
2016  {
2017  using ValueType = typename StencilT::ValueType;
2018 
2019  ValueType alpha, beta;
2020  if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta)) {
2021  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
2022  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
2023  }
2024  return 0;
2025  }
2026 };
2027 
2028 
2029 template<DDScheme DiffScheme2, DScheme DiffScheme1>
2030 struct MeanCurvature<UniformScaleTranslateMap, DiffScheme2, DiffScheme1>
2031 {
2032  // random access version
2033  template<typename Accessor> static typename Accessor::ValueType
2034  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
2035  {
2036  using ValueType = typename Accessor::ValueType;
2037 
2038  ValueType alpha, beta;
2039  if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta)) {
2040  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
2041  return ValueType(alpha*inv2dx/math::Pow3(beta));
2042  }
2043  return 0;
2044  }
2045 
2046  template<typename Accessor> static typename Accessor::ValueType
2047  normGrad(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
2048  {
2049  using ValueType = typename Accessor::ValueType;
2050 
2051  ValueType alpha, beta;
2052  if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta)) {
2053  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
2054  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
2055  }
2056  return 0;
2057  }
2058 
2059  // stencil access version
2060  template<typename StencilT> static typename StencilT::ValueType
2061  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
2062  {
2063  using ValueType = typename StencilT::ValueType;
2064 
2065  ValueType alpha, beta;
2066  if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta)) {
2067  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
2068  return ValueType(alpha*inv2dx/math::Pow3(beta));
2069  }
2070  return 0;
2071  }
2072 
2073  template<typename StencilT> static typename StencilT::ValueType
2074  normGrad(const UniformScaleTranslateMap& map, const StencilT& stencil)
2075  {
2076  using ValueType = typename StencilT::ValueType;
2077 
2078  ValueType alpha, beta;
2079  if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta)) {
2080  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
2081  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
2082  }
2083  return 0;
2084  }
2085 };
2086 
2087 
2093 {
2094 public:
2095  template<typename GridType>
2096  GenericMap(const GridType& g): mMap(g.transform().baseMap()) {}
2097 
2098  GenericMap(const Transform& t): mMap(t.baseMap()) {}
2099  GenericMap(MapBase::Ptr map): mMap(ConstPtrCast<const MapBase>(map)) {}
2100  GenericMap(MapBase::ConstPtr map): mMap(map) {}
2102 
2103  Vec3d applyMap(const Vec3d& in) const { return mMap->applyMap(in); }
2104  Vec3d applyInverseMap(const Vec3d& in) const { return mMap->applyInverseMap(in); }
2105 
2106  Vec3d applyIJT(const Vec3d& in) const { return mMap->applyIJT(in); }
2107  Vec3d applyIJT(const Vec3d& in, const Vec3d& pos) const { return mMap->applyIJT(in, pos); }
2108  Mat3d applyIJC(const Mat3d& m) const { return mMap->applyIJC(m); }
2109  Mat3d applyIJC(const Mat3d& m, const Vec3d& v, const Vec3d& pos) const
2110  { return mMap->applyIJC(m,v,pos); }
2111 
2112  double determinant() const { return mMap->determinant(); }
2113  double determinant(const Vec3d& in) const { return mMap->determinant(in); }
2114 
2115  Vec3d voxelSize() const { return mMap->voxelSize(); }
2116  Vec3d voxelSize(const Vec3d&v) const { return mMap->voxelSize(v); }
2117 
2118 private:
2119  MapBase::ConstPtr mMap;
2120 };
2121 
2122 } // end math namespace
2123 } // namespace OPENVDB_VERSION_NAME
2124 } // end openvdb namespace
2125 
2126 #endif // OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED
static Accessor::ValueType normGrad(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1830
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
Bracket code with OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN/_END, to inhibit warnings about type conve...
Definition: Platform.h:196
static StencilT::ValueType::value_type result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:966
static Accessor::ValueType result(const ScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1624
Definition: FiniteDifference.h:47
static Vec3< typename StencilT::ValueType > result(const StencilT &stencil, const Vec3Bias &V)
Definition: Operators.h:214
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1546
GenericMap(const Transform &t)
Definition: Operators.h:2098
Compute the closest-point transform to a level set.
Definition: Operators.h:1662
static Vec3< typename Accessor::ValueType > result(const Accessor &grid, const Coord &ijk, const Vec3Bias &V)
Definition: Operators.h:201
static StencilT::ValueType::value_type result(const ScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1253
static internal::ReturnValue< StencilT >::Vec3Type result(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:657
ResultType result(const StencilType &stencil)
Definition: Operators.h:44
static double result(const MapT &map, const AccessorType &grid, const Coord &ijk)
Definition: Operators.h:68
static bool compute(const MapType &map, const Accessor &grid, const Coord &ijk, double &alpha, double &beta)
Random access version.
Definition: Operators.h:1765
static internal::ReturnValue< StencilT >::Vec3Type result(const ScaleMap &map, const StencilT &stencil)
Definition: Operators.h:754
~GenericMap()
Definition: Operators.h:2101
static math::Vec3< typename StencilT::ValueType > result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:1688
static Accessor::ValueType::value_type result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1048
static double result(const StencilType &stencil)
Definition: Operators.h:59
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1973
Vec3d applyInverseMap(const Vec3d &in) const
Definition: Operators.h:2104
Definition: FiniteDifference.h:42
static StencilT::ValueType::value_type result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1029
Map traits.
Definition: Maps.h:55
double determinant() const
Definition: Operators.h:2112
A specialized Affine transform that scales along the principal axis the scaling need not be uniform i...
Definition: Maps.h:1155
static StencilT::ValueType result(const ScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1601
static Accessor::ValueType::value_type result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:473
Gradient operators defined in index space of various orders.
Definition: Operators.h:96
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:512
Adapter for vector-valued world-space operators to return the vector magnitude.
Definition: Operators.h:66
Definition: FiniteDifference.h:35
Definition: FiniteDifference.h:38
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:371
Laplacian defined in index space, using various center-difference stencils.
Definition: Operators.h:344
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1275
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:905
Abstract base class for maps.
Definition: Maps.h:134
Adapter to associate a map with a world-space operator, giving it the same call signature as an index...
Definition: Operators.h:35
Definition: FiniteDifference.h:1484
const Vec3d & getInvTwiceScale() const
Return 1/(2 scale). Used to optimize some finite difference calculations.
Definition: Maps.h:1330
typename T::ValueType ValueType
Definition: Operators.h:84
SharedPtr< MapBase > Ptr
Definition: Maps.h:137
Definition: FiniteDifference.h:40
static Accessor::ValueType::value_type result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1016
static Accessor::ValueType result(const ScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1583
Type Pow2(Type x)
Return x2.
Definition: Math.h:495
static internal::ReturnValue< Accessor >::Vec3Type result(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:649
static StencilT::ValueType result(const ScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1641
Definition: FiniteDifference.h:39
Divergence operator defined in index space using various first derivative schemes.
Definition: Operators.h:469
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:423
Compute the closest-point transform to a level set.
Definition: Operators.h:1715
static bool compute(const MapType &map, const StencilT &stencil, double &alpha, double &beta)
Stencil access version.
Definition: Operators.h:1844
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
Compute the mean curvature.
Definition: Operators.h:1758
static Accessor::ValueType::value_type result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1110
static StencilT::ValueType normGrad(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:2074
static StencilT::ValueType::value_type result(const StencilT &stencil)
Definition: Operators.h:482
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition: Coord.h:92
Compute the Laplacian at a given location in a grid using finite differencing of various orders...
Definition: Operators.h:1416
static Accessor::ValueType result(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1924
Definition: FiniteDifference.h:154
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_END
Definition: Platform.h:197
static Accessor::ValueType::value_type result(const ScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1208
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:444
static Accessor::ValueType::value_type result(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:988
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:361
DScheme
Different discrete schemes used in the first derivatives.
Definition: FiniteDifference.h:32
Definition: FiniteDifference.h:168
Biased gradient operators, defined with respect to the range-space of the map.
Definition: Operators.h:821
Compute the mean curvature in index space.
Definition: Operators.h:529
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1401
const Vec3d & getInvScaleSqr() const
Return the square of the scale. Used to optimize some finite difference calculations.
Definition: Maps.h:1328
Vec3d voxelSize(const Vec3d &v) const
Definition: Operators.h:2116
Vec3d applyIJT(const Vec3d &in, const Vec3d &pos) const
Definition: Operators.h:2107
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:2061
A specialized Affine transform that scales along the principal axis the scaling need not be uniform i...
Definition: Maps.h:656
static Vec3< typename StencilT::ValueType > result(const StencilT &stencil)
Definition: Operators.h:111
SharedPtr< T > ConstPtrCast(const SharedPtr< U > &ptr)
Return a new shared pointer that points to the same object as the given pointer but with possibly dif...
Definition: Types.h:103
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1327
Definition: FiniteDifference.h:169
Vec3d voxelSize() const
Definition: Operators.h:2115
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
Curl operator defined in index space using various first derivative schemes.
Definition: Operators.h:495
static Vec3< typename Accessor::ValueType > result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:100
Definition: Operators.h:854
A specialized linear transform that performs a unitary maping i.e. rotation and or reflection...
Definition: Maps.h:1605
GenericMap(MapBase::ConstPtr map)
Definition: Operators.h:2100
static math::Vec3< typename StencilT::ValueType > result(const MapType &map, const StencilT &stencil, const Vec3< typename StencilT::ValueType > &V)
Definition: Operators.h:837
Vec3d asVec3d() const
Definition: Coord.h:144
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1569
static internal::ReturnValue< StencilT >::Vec3Type result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:684
static double result(const MapT &map, const StencilType &stencil)
Definition: Operators.h:73
static Accessor::ValueType::value_type result(const ScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1142
static internal::ReturnValue< Accessor >::Vec3Type result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:703
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:402
Definition: FiniteDifference.h:153
Definition: FiniteDifference.h:1762
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:2002
static StencilT::ValueType::value_type result(const ScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1155
const Vec3d & getInvScale() const
Return 1/(scale)
Definition: Maps.h:1332
Definition: Stencils.h:240
MapAdapter(const MapType &m)
Definition: Operators.h:36
Definition: Operators.h:126
Definition: FiniteDifference.h:43
static internal::ReturnValue< Accessor >::Vec3Type result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:621
Type Pow3(Type x)
Return x3.
Definition: Math.h:499
Definition: FiniteDifference.h:41
static internal::ReturnValue< Accessor >::Vec3Type result(const ScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:779
SharedPtr< const MapBase > ConstPtr
Definition: Maps.h:138
static StencilT::ValueType result(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:1504
Definition: Exceptions.h:13
static StencilT::ValueType::value_type result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1061
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:863
static Vec3< typename StencilT::ValueType > result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:1736
const Vec3d & getInvTwiceScale() const
Return 1/(2 scale). Used to optimize some finite difference calculations.
Definition: Maps.h:795
static StencilT::ValueType result(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:1947
Compute the curl of a vector-valued grid using differencing of various orders in the space defined by...
Definition: Operators.h:1271
A specialized Affine transform that uniformaly scales along the principal axis and then translates th...
Definition: Maps.h:1463
static StencilT::ValueType normGrad(const MapType &map, const StencilT stencil)
Definition: Operators.h:1909
GenericMap(const GridType &g)
Definition: Operators.h:2096
static Accessor::ValueType normGrad(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1987
static StencilT::ValueType result(const UnitaryMap &, const StencilT &stencil)
Definition: Operators.h:1525
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1318
static StencilT::ValueType::value_type result(const ScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1221
Biased Gradient Operators, using upwinding defined by the Vec3Bias input.
Definition: Operators.h:193
Definition: FiniteDifference.h:152
Vec3d applyMap(const Vec3d &in) const
Definition: Operators.h:2103
static Vec3< typename Accessor::ValueType > result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1719
static internal::ReturnValue< Accessor >::Vec3Type result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:671
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1537
Definition: FiniteDifference.h:170
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:249
static Accessor::ValueType result(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1496
Definition: FiniteDifference.h:171
static StencilT::ValueType result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:876
Center difference gradient operators, defined with respect to the range-space of the map...
Definition: Operators.h:616
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:385
Definition: Transform.h:39
A specialized linear transform that performs a translation.
Definition: Maps.h:972
Definition: Operators.h:22
static Accessor::ValueType result(const UnitaryMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1517
Tolerance for floating-point comparison.
Definition: Math.h:90
static math::Vec3< typename Accessor::ValueType > result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1666
static Accessor::ValueType::value_type result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1080
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1420
Definition: FiniteDifference.h:416
static Accessor::ValueType::value_type result(const ScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1175
Real GodunovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)
Definition: FiniteDifference.h:326
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1376
static double result(const AccessorType &grid, const Coord &ijk)
Definition: Operators.h:54
static math::Vec3< typename Accessor::ValueType > result(const MapType &map, const Accessor &grid, const Coord &ijk, const Vec3< typename Accessor::ValueType > &V)
Definition: Operators.h:825
const Vec3d & getInvScale() const
Return 1/(scale)
Definition: Maps.h:797
static internal::ReturnValue< StencilT >::Vec3Type result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:716
A specialized Affine transform that scales along the principal axis the scaling is uniform in the thr...
Definition: Maps.h:900
double determinant(const Vec3d &in) const
Definition: Operators.h:2113
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1820
static StencilT::ValueType::value_type result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1123
static StencilT::ValueType::value_type result(const ScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1188
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1341
Definition: Operators.h:227
static bool result(const StencilT &stencil, typename StencilT::ValueType &alpha, typename StencilT::ValueType &beta)
Stencil access version.
Definition: Operators.h:574
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:236
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:921
static Accessor::ValueType normGrad(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:2047
static Accessor::ValueType::value_type result(const ScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1240
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
static StencilT::ValueType normGrad(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:2015
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:894
Definition: FiniteDifference.h:167
Mat3d applyIJC(const Mat3d &m) const
Definition: Operators.h:2108
static StencilT::ValueType result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:1294
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1560
ResultType result(const AccessorType &grid, const Coord &ijk)
Definition: Operators.h:40
static StencilT::ValueType result(const MapType &map, const StencilT stencil)
Definition: Operators.h:1900
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1391
Definition: FiniteDifference.h:46
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1351
static internal::ReturnValue< Accessor >::Vec3Type result(const ScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:735
A wrapper that holds a MapBase::ConstPtr and exposes a reduced set of functionality needed by the mat...
Definition: Operators.h:2092
static StencilT::ValueType::value_type result(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:999
GenericMap(MapBase::Ptr map)
Definition: Operators.h:2099
static Accessor::ValueType::value_type result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:950
static internal::ReturnValue< StencilT >::Vec3Type result(const ScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:798
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:932
Adapter for vector-valued index-space operators to return the vector magnitude.
Definition: Operators.h:52
static bool result(const Accessor &grid, const Coord &ijk, typename Accessor::ValueType &alpha, typename Accessor::ValueType &beta)
Random access version.
Definition: Operators.h:536
static StencilT::ValueType::value_type result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1092
Definition: FiniteDifference.h:45
Vec3d applyIJT(const Vec3d &in) const
Definition: Operators.h:2106
static internal::ReturnValue< StencilT >::Vec3Type result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:632
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:499
static StencilT::ValueType result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:1456
Mat3d applyIJC(const Mat3d &m, const Vec3d &v, const Vec3d &pos) const
Definition: Operators.h:2109
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1366
Compute the divergence of a vector-valued grid using differencing of various orders, the result defined with respect to the range-space of the map.
Definition: Operators.h:946
const MapType map
Definition: Operators.h:46
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:2034
static Accessor::ValueType normGrad(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1935
static StencilT::ValueType normGrad(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:1957
Definition: FiniteDifference.h:44
const Vec3d & getInvScaleSqr() const
Return the square of the scale. Used to optimize some finite difference calculations.
Definition: Maps.h:793
ValueType WENO5(const ValueType &v1, const ValueType &v2, const ValueType &v3, const ValueType &v4, const ValueType &v5, float scale2=0.01f)
Implementation of nominally fifth-order finite-difference WENO.
Definition: FiniteDifference.h:304
Definition: Operators.h:25