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