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