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