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