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