OpenVDB  1.2.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  typedef typename Accessor::ValueType ValueType;
294  typedef math::Vec3<ValueType> Vec3Type;
295 
296  // SSE optimized
297  const simd::Float4
298  v1(grid.getValue(ijk.offsetBy(-2, 0, 0)) - grid.getValue(ijk.offsetBy(-3, 0, 0)),
299  grid.getValue(ijk.offsetBy( 0,-2, 0)) - grid.getValue(ijk.offsetBy( 0,-3, 0)),
300  grid.getValue(ijk.offsetBy( 0, 0,-2)) - grid.getValue(ijk.offsetBy( 0, 0,-3)), 0),
301  v2(grid.getValue(ijk.offsetBy(-1, 0, 0)) - grid.getValue(ijk.offsetBy(-2, 0, 0)),
302  grid.getValue(ijk.offsetBy( 0,-1, 0)) - grid.getValue(ijk.offsetBy( 0,-2, 0)),
303  grid.getValue(ijk.offsetBy( 0, 0,-1)) - grid.getValue(ijk.offsetBy( 0, 0,-2)), 0),
304  v3(grid.getValue(ijk ) - grid.getValue(ijk.offsetBy(-1, 0, 0)),
305  grid.getValue(ijk ) - grid.getValue(ijk.offsetBy( 0,-1, 0)),
306  grid.getValue(ijk ) - grid.getValue(ijk.offsetBy( 0, 0,-1)), 0),
307  v4(grid.getValue(ijk.offsetBy( 1, 0, 0)) - grid.getValue(ijk ),
308  grid.getValue(ijk.offsetBy( 0, 1, 0)) - grid.getValue(ijk ),
309  grid.getValue(ijk.offsetBy( 0, 0, 1)) - grid.getValue(ijk ), 0),
310  v5(grid.getValue(ijk.offsetBy( 2, 0, 0)) - grid.getValue(ijk.offsetBy( 1, 0, 0)),
311  grid.getValue(ijk.offsetBy( 0, 2, 0)) - grid.getValue(ijk.offsetBy( 0, 1, 0)),
312  grid.getValue(ijk.offsetBy( 0, 0, 2)) - grid.getValue(ijk.offsetBy( 0, 0, 1)), 0),
313  v6(grid.getValue(ijk.offsetBy( 3, 0, 0)) - grid.getValue(ijk.offsetBy( 2, 0, 0)),
314  grid.getValue(ijk.offsetBy( 0, 3, 0)) - grid.getValue(ijk.offsetBy( 0, 2, 0)),
315  grid.getValue(ijk.offsetBy( 0, 0, 3)) - grid.getValue(ijk.offsetBy( 0, 0, 2)), 0),
316  down = math::WENO5(v1, v2, v3, v4, v5),
317  up = math::WENO5(v6, v5, v4, v3, v2);
318 
319  return math::GudonovsNormSqrd(grid.getValue(ijk)>0, down, up);
320  }
321 
322  // stencil access version
323  template<typename StencilT>
324  static typename StencilT::ValueType
325  result(const StencilT& s)
326  {
327  typedef typename StencilT::ValueType ValueType;
328  typedef math::Vec3<ValueType> Vec3Type;
329 
330  // SSE optimized
331  const simd::Float4
332  v1(s.template getValue<-2, 0, 0>() - s.template getValue<-3, 0, 0>(),
333  s.template getValue< 0,-2, 0>() - s.template getValue< 0,-3, 0>(),
334  s.template getValue< 0, 0,-2>() - s.template getValue< 0, 0,-3>(), 0),
335  v2(s.template getValue<-1, 0, 0>() - s.template getValue<-2, 0, 0>(),
336  s.template getValue< 0,-1, 0>() - s.template getValue< 0,-2, 0>(),
337  s.template getValue< 0, 0,-1>() - s.template getValue< 0, 0,-2>(), 0),
338  v3(s.template getValue< 0, 0, 0>() - s.template getValue<-1, 0, 0>(),
339  s.template getValue< 0, 0, 0>() - s.template getValue< 0,-1, 0>(),
340  s.template getValue< 0, 0, 0>() - s.template getValue< 0, 0,-1>(), 0),
341  v4(s.template getValue< 1, 0, 0>() - s.template getValue< 0, 0, 0>(),
342  s.template getValue< 0, 1, 0>() - s.template getValue< 0, 0, 0>(),
343  s.template getValue< 0, 0, 1>() - s.template getValue< 0, 0, 0>(), 0),
344  v5(s.template getValue< 2, 0, 0>() - s.template getValue< 1, 0, 0>(),
345  s.template getValue< 0, 2, 0>() - s.template getValue< 0, 1, 0>(),
346  s.template getValue< 0, 0, 2>() - s.template getValue< 0, 0, 1>(), 0),
347  v6(s.template getValue< 3, 0, 0>() - s.template getValue< 2, 0, 0>(),
348  s.template getValue< 0, 3, 0>() - s.template getValue< 0, 2, 0>(),
349  s.template getValue< 0, 0, 3>() - s.template getValue< 0, 0, 2>(), 0),
350  down = math::WENO5(v1, v2, v3, v4, v5),
351  up = math::WENO5(v6, v5, v4, v3, v2);
352 
353  return math::GudonovsNormSqrd(s.template getValue<0, 0, 0>()>0, down, up);
354  }
355 };
356 #endif //DWA_OPENVDB // for SIMD - note will do the computations in float
357 
358 
359 
361 template<DDScheme DiffScheme>
364 {
365  // random access version
366  template<typename Accessor>
367  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk);
368 
369  // stencil access version
370  template<typename StencilT>
371  static typename StencilT::ValueType result(const StencilT& stencil);
372 };
373 
374 
375 template<>
377 {
378  // random access version
379  template<typename Accessor>
380  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk)
381  {
382  return grid.getValue(ijk.offsetBy(1,0,0)) + grid.getValue(ijk.offsetBy(-1, 0, 0)) +
383  grid.getValue(ijk.offsetBy(0,1,0)) + grid.getValue(ijk.offsetBy(0, -1, 0)) +
384  grid.getValue(ijk.offsetBy(0,0,1)) + grid.getValue(ijk.offsetBy(0, 0,-1))
385  - 6*grid.getValue(ijk);
386  }
387 
388  // stencil access version
389  template<typename StencilT>
390  static typename StencilT::ValueType result(const StencilT& stencil)
391  {
392  return stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() +
393  stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() +
394  stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>()
395  - 6*stencil.template getValue< 0, 0, 0>();
396  }
397 };
398 
399 template<>
401 {
402  // random access version
403  template<typename Accessor>
404  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk)
405  {
406  return (-1./12.)*(
407  grid.getValue(ijk.offsetBy(2,0,0)) + grid.getValue(ijk.offsetBy(-2, 0, 0)) +
408  grid.getValue(ijk.offsetBy(0,2,0)) + grid.getValue(ijk.offsetBy( 0,-2, 0)) +
409  grid.getValue(ijk.offsetBy(0,0,2)) + grid.getValue(ijk.offsetBy( 0, 0,-2)) )
410  + (4./3.)*(
411  grid.getValue(ijk.offsetBy(1,0,0)) + grid.getValue(ijk.offsetBy(-1, 0, 0)) +
412  grid.getValue(ijk.offsetBy(0,1,0)) + grid.getValue(ijk.offsetBy( 0,-1, 0)) +
413  grid.getValue(ijk.offsetBy(0,0,1)) + grid.getValue(ijk.offsetBy( 0, 0,-1)) )
414  - 7.5*grid.getValue(ijk);
415  }
416 
417  // stencil access version
418  template<typename StencilT>
419  static typename StencilT::ValueType result(const StencilT& stencil)
420  {
421  return (-1./12.)*(
422  stencil.template getValue< 2, 0, 0>() + stencil.template getValue<-2, 0, 0>() +
423  stencil.template getValue< 0, 2, 0>() + stencil.template getValue< 0,-2, 0>() +
424  stencil.template getValue< 0, 0, 2>() + stencil.template getValue< 0, 0,-2>() )
425  + (4./3.)*(
426  stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() +
427  stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() +
428  stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>() )
429  - 7.5*stencil.template getValue< 0, 0, 0>();
430  }
431 };
432 
433 template<>
435 {
436  // random access version
437  template<typename Accessor>
438  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk)
439  {
440  return (1./90.)*(
441  grid.getValue(ijk.offsetBy(3,0,0)) + grid.getValue(ijk.offsetBy(-3, 0, 0)) +
442  grid.getValue(ijk.offsetBy(0,3,0)) + grid.getValue(ijk.offsetBy( 0,-3, 0)) +
443  grid.getValue(ijk.offsetBy(0,0,3)) + grid.getValue(ijk.offsetBy( 0, 0,-3)) )
444  - (3./20.)*(
445  grid.getValue(ijk.offsetBy(2,0,0)) + grid.getValue(ijk.offsetBy(-2, 0, 0)) +
446  grid.getValue(ijk.offsetBy(0,2,0)) + grid.getValue(ijk.offsetBy( 0,-2, 0)) +
447  grid.getValue(ijk.offsetBy(0,0,2)) + grid.getValue(ijk.offsetBy( 0, 0,-2)) )
448  + 1.5 *(
449  grid.getValue(ijk.offsetBy(1,0,0)) + grid.getValue(ijk.offsetBy(-1, 0, 0)) +
450  grid.getValue(ijk.offsetBy(0,1,0)) + grid.getValue(ijk.offsetBy( 0,-1, 0)) +
451  grid.getValue(ijk.offsetBy(0,0,1)) + grid.getValue(ijk.offsetBy( 0, 0,-1)) )
452  - (3*49/18.)*grid.getValue(ijk);
453  }
454 
455  // stencil access version
456  template<typename StencilT>
457  static typename StencilT::ValueType result(const StencilT& stencil)
458  {
459  return (1./90.)*(
460  stencil.template getValue< 3, 0, 0>() + stencil.template getValue<-3, 0, 0>() +
461  stencil.template getValue< 0, 3, 0>() + stencil.template getValue< 0,-3, 0>() +
462  stencil.template getValue< 0, 0, 3>() + stencil.template getValue< 0, 0,-3>() )
463  - (3./20.)*(
464  stencil.template getValue< 2, 0, 0>() + stencil.template getValue<-2, 0, 0>() +
465  stencil.template getValue< 0, 2, 0>() + stencil.template getValue< 0,-2, 0>() +
466  stencil.template getValue< 0, 0, 2>() + stencil.template getValue< 0, 0,-2>() )
467  + 1.5 *(
468  stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() +
469  stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() +
470  stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>() )
471  - (3*49/18.)*stencil.template getValue< 0, 0, 0>();
472  }
473 };
475 
476 
478 template<DScheme DiffScheme>
481 {
482  // random access version
483  template<typename Accessor> static typename Accessor::ValueType::value_type
484  result(const Accessor& grid, const Coord& ijk)
485  {
486  return D1Vec<DiffScheme>::inX(grid, ijk, 0) +
487  D1Vec<DiffScheme>::inY(grid, ijk, 1) +
488  D1Vec<DiffScheme>::inZ(grid, ijk, 2);
489  }
490 
491  // stencil access version
492  template<typename StencilT> static typename StencilT::ValueType::value_type
493  result(const StencilT& stencil)
494  {
495  return D1Vec<DiffScheme>::inX(stencil, 0) +
496  D1Vec<DiffScheme>::inY(stencil, 1) +
497  D1Vec<DiffScheme>::inZ(stencil, 2);
498  }
499 };
501 
502 
504 template<DScheme DiffScheme>
506 struct ISCurl
507 {
508  // random access version
509  template<typename Accessor>
510  static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk)
511  {
512  typedef typename Accessor::ValueType Vec3Type;
513  return Vec3Type( D1Vec<DiffScheme>::inY(grid, ijk, 2) - //dw/dy - dv/dz
514  D1Vec<DiffScheme>::inZ(grid, ijk, 1),
515  D1Vec<DiffScheme>::inZ(grid, ijk, 0) - //du/dz - dw/dx
516  D1Vec<DiffScheme>::inX(grid, ijk, 2),
517  D1Vec<DiffScheme>::inX(grid, ijk, 1) - //dv/dx - du/dy
518  D1Vec<DiffScheme>::inY(grid, ijk, 0) );
519  }
520 
521  // stencil access version
522  template<typename StencilT>
523  static typename StencilT::ValueType result(const StencilT& stencil)
524  {
525  typedef typename StencilT::ValueType Vec3Type;
526  return Vec3Type( D1Vec<DiffScheme>::inY(stencil, 2) - //dw/dy - dv/dz
527  D1Vec<DiffScheme>::inZ(stencil, 1),
528  D1Vec<DiffScheme>::inZ(stencil, 0) - //du/dz - dw/dx
529  D1Vec<DiffScheme>::inX(stencil, 2),
530  D1Vec<DiffScheme>::inX(stencil, 1) - //dv/dx - du/dy
531  D1Vec<DiffScheme>::inY(stencil, 0) );
532  }
533 };
535 
536 
538 template<DDScheme DiffScheme2, DScheme DiffScheme1>
541 {
542  // random access version
543  template<typename Accessor>
544  static void result(const Accessor& grid, const Coord& ijk,
545  typename Accessor::ValueType& alpha,
546  typename Accessor::ValueType& beta)
547  {
548  typedef typename Accessor::ValueType ValueType;
549 
550  ValueType Dx = D1<DiffScheme1>::inX(grid, ijk);
551  ValueType Dy = D1<DiffScheme1>::inY(grid, ijk);
552  ValueType Dz = D1<DiffScheme1>::inZ(grid, ijk);
553 
554  ValueType Dx2 = Dx*Dx;
555  ValueType Dy2 = Dy*Dy;
556  ValueType Dz2 = Dz*Dz;
557 
558  ValueType Dxx = D2<DiffScheme2>::inX(grid, ijk);
559  ValueType Dyy = D2<DiffScheme2>::inY(grid, ijk);
560  ValueType Dzz = D2<DiffScheme2>::inZ(grid, ijk);
561 
562  ValueType Dxy = D2<DiffScheme2>::inXandY(grid, ijk);
563  ValueType Dyz = D2<DiffScheme2>::inYandZ(grid, ijk);
564  ValueType Dxz = D2<DiffScheme2>::inXandZ(grid, ijk);
565 
566  // for return
567  alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz));
568  beta = ValueType(std::sqrt(double(Dx2 + Dy2 + Dz2))); // * 1/dx
569  }
570 
571  // stencil access version
572  template<typename StencilT>
573  static void result(const StencilT& stencil,
574  typename StencilT::ValueType& alpha,
575  typename StencilT::ValueType& beta)
576  {
577  typedef typename StencilT::ValueType ValueType;
578  ValueType Dx = D1<DiffScheme1>::inX(stencil);
579  ValueType Dy = D1<DiffScheme1>::inY(stencil);
580  ValueType Dz = D1<DiffScheme1>::inZ(stencil);
581 
582  ValueType Dx2 = Dx*Dx;
583  ValueType Dy2 = Dy*Dy;
584  ValueType Dz2 = Dz*Dz;
585 
586  ValueType Dxx = D2<DiffScheme2>::inX(stencil);
587  ValueType Dyy = D2<DiffScheme2>::inY(stencil);
588  ValueType Dzz = D2<DiffScheme2>::inZ(stencil);
589 
590  ValueType Dxy = D2<DiffScheme2>::inXandY(stencil);
591  ValueType Dyz = D2<DiffScheme2>::inYandZ(stencil);
592  ValueType Dxz = D2<DiffScheme2>::inXandZ(stencil);
593 
594  // for return
595  alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz));
596  beta = ValueType(std::sqrt(double(Dx2 + Dy2 + Dz2))); // * 1/dx
597  }
598 };
599 
600 
601 // --- Operators defined in the Range of a given map
602 
603 
605 template<typename MapType, DScheme DiffScheme>
609 struct Gradient
610 {
611  // random access version
612  template<typename Accessor>
614  result(const MapType& map, const Accessor& grid, const Coord& ijk)
615  {
616  typedef typename internal::ReturnValue<Accessor>::Vec3Type Vec3Type;
617 
618  Vec3d iGradient( ISGradient<DiffScheme>::result(grid, ijk) );
619  return Vec3Type(map.applyIJT(iGradient, ijk.asVec3d()));
620  }
621 
622  // stencil access version
623  template<typename StencilT>
625  result(const MapType& map, const StencilT& stencil)
626  {
627  typedef typename internal::ReturnValue<StencilT>::Vec3Type Vec3Type;
628 
629  Vec3d iGradient( ISGradient<DiffScheme>::result(stencil) );
630  return Vec3Type(map.applyIJT(iGradient, stencil.getCenterCoord().asVec3d()));
631  }
632 };
633 
634 
635 // translation, any order
636 template<DScheme DiffScheme>
637 struct Gradient<TranslationMap, DiffScheme>
638 {
639  // random access version
640  template<typename Accessor>
642  result(const TranslationMap&, const Accessor& grid, const Coord& ijk)
643  {
644  return ISGradient<DiffScheme>::result(grid, ijk);
645  }
646 
647  // stencil access version
648  template<typename StencilT>
650  result(const TranslationMap&, const StencilT& stencil)
651  {
652  return ISGradient<DiffScheme>::result(stencil);
653  }
654 };
655 
656 
657 // uniform scale, 2nd order
658 template<>
660 {
661  // random access version
662  template<typename Accessor>
664  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
665  {
666  typedef typename internal::ReturnValue<Accessor>::ValueType ValueType;
667  typedef typename internal::ReturnValue<Accessor>::Vec3Type Vec3Type;
668 
669  Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) );
670  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
671  return iGradient * inv2dx;
672  }
673 
674  // stencil access version
675  template<typename StencilT>
677  result(const UniformScaleMap& map, const StencilT& stencil)
678  {
679  typedef typename internal::ReturnValue<StencilT>::ValueType ValueType;
680  typedef typename internal::ReturnValue<StencilT>::Vec3Type Vec3Type;
681 
682  Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) );
683  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
684  return iGradient * inv2dx;
685  }
686 };
687 
688 
689 // uniform scale translate, 2nd order
690 template<>
692 {
693  // random access version
694  template<typename Accessor>
696  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
697  {
698  typedef typename internal::ReturnValue<Accessor>::ValueType ValueType;
699  typedef typename internal::ReturnValue<Accessor>::Vec3Type Vec3Type;
700 
701  Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) );
702  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
703  return iGradient * inv2dx;
704  }
705 
706  // stencil access version
707  template<typename StencilT>
709  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
710  {
711  typedef typename internal::ReturnValue<StencilT>::ValueType ValueType;
712  typedef typename internal::ReturnValue<StencilT>::Vec3Type Vec3Type;
713 
714  Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) );
715  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
716  return iGradient * inv2dx;
717  }
718 };
719 
720 
721 // scale, 2nd order
722 template<>
724 {
725  // random access version
726  template<typename Accessor>
728  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
729  {
730  typedef typename internal::ReturnValue<Accessor>::ValueType ValueType;
731  typedef typename internal::ReturnValue<Accessor>::Vec3Type Vec3Type;
732 
733  Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) );
734  return Vec3Type(ValueType(iGradient[0] * map.getInvTwiceScale()[0]),
735  ValueType(iGradient[1] * map.getInvTwiceScale()[1]),
736  ValueType(iGradient[2] * map.getInvTwiceScale()[2]) );
737  }
738 
739  // stencil access version
740  template<typename StencilT>
742  result(const ScaleMap& map, const StencilT& stencil)
743  {
744  typedef typename internal::ReturnValue<StencilT>::ValueType ValueType;
745  typedef typename internal::ReturnValue<StencilT>::Vec3Type Vec3Type;
746 
747  Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) );
748  return Vec3Type(ValueType(iGradient[0] * map.getInvTwiceScale()[0]),
749  ValueType(iGradient[1] * map.getInvTwiceScale()[1]),
750  ValueType(iGradient[2] * map.getInvTwiceScale()[2]) );
751  }
752 };
753 
754 
755 // scale translate, 2nd order
756 template<>
758 {
759  // random access version
760  template<typename Accessor>
762  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
763  {
764  typedef typename internal::ReturnValue<Accessor>::ValueType ValueType;
765  typedef typename internal::ReturnValue<Accessor>::Vec3Type Vec3Type;
766 
767  Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) );
768  return Vec3Type(ValueType(iGradient[0] * map.getInvTwiceScale()[0]),
769  ValueType(iGradient[1] * map.getInvTwiceScale()[1]),
770  ValueType(iGradient[2] * map.getInvTwiceScale()[2]) );
771  }
772 
773  // Stencil access version
774  template<typename StencilT>
776  result(const ScaleTranslateMap& map, const StencilT& stencil)
777  {
778  typedef typename internal::ReturnValue<StencilT>::ValueType ValueType;
779  typedef typename internal::ReturnValue<StencilT>::Vec3Type Vec3Type;
780 
781  Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) );
782  return Vec3Type(ValueType(iGradient[0] * map.getInvTwiceScale()[0]),
783  ValueType(iGradient[1] * map.getInvTwiceScale()[1]),
784  ValueType(iGradient[2] * map.getInvTwiceScale()[2]) );
785  }
786 };
788 
789 
791 template<typename MapType, BiasedGradientScheme GradScheme>
795 {
796  // random access version
797  template<typename Accessor> static math::Vec3<typename Accessor::ValueType>
798  result(const MapType& map, const Accessor& grid, const Coord& ijk,
800  {
801  typedef typename Accessor::ValueType ValueType;
802  typedef math::Vec3<ValueType> Vec3Type;
803 
804  Vec3d iGradient( ISGradientBiased<GradScheme, Vec3Type>::result(grid, ijk, V) );
805  return Vec3Type(map.applyIJT(iGradient, ijk.asVec3d()));
806  }
807 
808  // stencil access version
809  template<typename StencilT> static math::Vec3<typename StencilT::ValueType>
810  result(const MapType& map, const StencilT& stencil,
812  {
813  typedef typename StencilT::ValueType ValueType;
814  typedef math::Vec3<ValueType> Vec3Type;
815 
816  Vec3d iGradient( ISGradientBiased<GradScheme, Vec3Type>::result(stencil, V) );
817  return Vec3Type(map.applyIJT(iGradient, stencil.getCenterCoord().asVec3d()));
818  }
819 };
821 
822 
823 template<typename MapType, BiasedGradientScheme GradScheme>
825 {
828 
829 
830  // random access version
831  template<typename Accessor>
832  static typename Accessor::ValueType
833  result(const MapType& map, const Accessor& grid, const Coord& ijk)
834  {
835  typedef typename Accessor::ValueType ValueType;
836  typedef math::Vec3<ValueType> Vec3Type;
837 
838  Vec3Type up = Gradient<MapType, FD>::result(map, grid, ijk);
839  Vec3Type down = Gradient<MapType, BD>::result(map, grid, ijk);
840  return math::GudonovsNormSqrd(grid.getValue(ijk)>0, down, up);
841  }
842 
843  // stencil access version
844  template<typename StencilT>
845  static typename StencilT::ValueType
846  result(const MapType& map, const StencilT& stencil)
847  {
848  typedef typename StencilT::ValueType ValueType;
849  typedef math::Vec3<ValueType> Vec3Type;
850 
851  Vec3Type up = Gradient<MapType, FD>::result(map, stencil);
852  Vec3Type down = Gradient<MapType, BD>::result(map, stencil);
853  return math::GudonovsNormSqrd(stencil.template getValue<0, 0, 0>()>0, down, up);
854  }
855 };
856 
857 
858 template<BiasedGradientScheme GradScheme>
860 {
861  // random access version
862  template<typename Accessor>
863  static typename Accessor::ValueType
864  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
865  {
866  typedef typename Accessor::ValueType ValueType;
867 
868  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
869  return invdxdx * ISGradientNormSqrd<GradScheme>::result(grid, ijk);
870  }
871 
872  // stencil access version
873  template<typename StencilT>
874  static typename StencilT::ValueType
875  result(const UniformScaleMap& map, const StencilT& stencil)
876  {
877  typedef typename StencilT::ValueType ValueType;
878 
879  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
880  return invdxdx * ISGradientNormSqrd<GradScheme>::result(stencil);
881  }
882 };
883 
884 template<BiasedGradientScheme GradScheme>
886 {
887  // random access version
888  template<typename Accessor>
889  static typename Accessor::ValueType
890  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
891  {
892  typedef typename Accessor::ValueType ValueType;
893 
894  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
895  return invdxdx * ISGradientNormSqrd<GradScheme>::result(grid, ijk);
896  }
897 
898  // stencil access version
899  template<typename StencilT>
900  static typename StencilT::ValueType
901  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
902  {
903  typedef typename StencilT::ValueType ValueType;
904 
905  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
906  return invdxdx * ISGradientNormSqrd<GradScheme>::result(stencil);
907  }
908 };
909 
910 
912 template<typename MapType, DScheme DiffScheme>
916 {
917  // random access version
918  template<typename Accessor> static typename Accessor::ValueType::value_type
919  result(const MapType& map, const Accessor& grid, const Coord& ijk)
920  {
921  typedef typename Accessor::ValueType Vec3Type;
922  typedef typename Accessor::ValueType::value_type ValueType;
923 
924  ValueType div(0);
925  for (int i=0; i < 3; i++) {
926  Vec3d vec( D1Vec<DiffScheme>::inX(grid, ijk, i),
927  D1Vec<DiffScheme>::inY(grid, ijk, i),
928  D1Vec<DiffScheme>::inZ(grid, ijk, i) );
929  div += ValueType(map.applyIJT(vec, ijk.asVec3d())[i]);
930  }
931  return div;
932  }
933 
934  // stencil access version
935  template<typename StencilT> static typename StencilT::ValueType::value_type
936  result(const MapType& map, const StencilT& stencil)
937  {
938  typedef typename StencilT::ValueType Vec3Type;
939  typedef typename StencilT::ValueType::value_type ValueType;
940 
941  ValueType div(0);
942  for (int i=0; i < 3; i++) {
943  Vec3d vec( D1Vec<DiffScheme>::inX(stencil, i),
944  D1Vec<DiffScheme>::inY(stencil, i),
945  D1Vec<DiffScheme>::inZ(stencil, i) );
946  div += ValueType(map.applyIJT(vec, stencil.getCenterCoord().asVec3d())[i]);
947  }
948  return div;
949  }
950 };
951 
952 
953 // translation, any scheme
954 template<DScheme DiffScheme>
955 struct Divergence<TranslationMap, DiffScheme>
956 {
957  // random access version
958  template<typename Accessor> static typename Accessor::ValueType::value_type
959  result(const TranslationMap&, const Accessor& grid, const Coord& ijk)
960  {
961  typedef typename Accessor::ValueType Vec3Type;
962  typedef typename Accessor::ValueType::value_type ValueType;
963 
964  ValueType div(0);
965  div =ISDivergence<DiffScheme>::result(grid, ijk);
966  return div;
967  }
968 
969  // stencil access version
970  template<typename StencilT> static typename StencilT::ValueType::value_type
971  result(const TranslationMap&, const StencilT& stencil)
972  {
973  typedef typename StencilT::ValueType Vec3Type;
974  typedef typename StencilT::ValueType::value_type ValueType;
975 
976  ValueType div(0);
977  div =ISDivergence<DiffScheme>::result(stencil);
978  return div;
979  }
980 };
981 
982 
983 // uniform scale, any scheme
984 template<DScheme DiffScheme>
985 struct Divergence<UniformScaleMap, DiffScheme>
986 {
987  // random access version
988  template<typename Accessor> static typename Accessor::ValueType::value_type
989  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
990  {
991  typedef typename Accessor::ValueType Vec3Type;
992  typedef typename Accessor::ValueType::value_type ValueType;
993 
994  ValueType div(0);
995 
996  div =ISDivergence<DiffScheme>::result(grid, ijk);
997  ValueType invdx = ValueType(map.getInvScale()[0]);
998  return div * invdx;
999  }
1000 
1001  // stencil access version
1002  template<typename StencilT> static typename StencilT::ValueType::value_type
1003  result(const UniformScaleMap& map, const StencilT& stencil)
1004  {
1005  typedef typename StencilT::ValueType Vec3Type;
1006  typedef typename StencilT::ValueType::value_type ValueType;
1007 
1008  ValueType div(0);
1009 
1010  div =ISDivergence<DiffScheme>::result(stencil);
1011  ValueType invdx = ValueType(map.getInvScale()[0]);
1012  return div * invdx;
1013  }
1014 };
1015 
1016 
1017 // uniform scale and translation, any scheme
1018 template<DScheme DiffScheme>
1020 {
1021  // random access version
1022  template<typename Accessor> static typename Accessor::ValueType::value_type
1023  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1024  {
1025  typedef typename Accessor::ValueType Vec3Type;
1026  typedef typename Accessor::ValueType::value_type ValueType;
1027 
1028  ValueType div(0);
1029 
1030  div =ISDivergence<DiffScheme>::result(grid, ijk);
1031  ValueType invdx = ValueType(map.getInvScale()[0]);
1032  return div * invdx;
1033  }
1034 
1035  // stencil access version
1036  template<typename StencilT> static typename StencilT::ValueType::value_type
1037  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1038  {
1039  typedef typename StencilT::ValueType Vec3Type;
1040  typedef typename StencilT::ValueType::value_type ValueType;
1041 
1042  ValueType div(0);
1043 
1044  div =ISDivergence<DiffScheme>::result(stencil);
1045  ValueType invdx = ValueType(map.getInvScale()[0]);
1046  return div * invdx;
1047  }
1048 };
1049 
1050 
1051 // uniform scale 2nd order
1052 template<>
1054 {
1055  // random access version
1056  template<typename Accessor> static typename Accessor::ValueType::value_type
1057  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1058  {
1059  typedef typename Accessor::ValueType Vec3Type;
1060  typedef typename Accessor::ValueType::value_type ValueType;
1061 
1062  ValueType div(0);
1063  div =ISDivergence<CD_2NDT>::result(grid, ijk);
1064  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1065  return div * inv2dx;
1066  }
1067 
1068  // stencil access version
1069  template<typename StencilT> static typename StencilT::ValueType::value_type
1070  result(const UniformScaleMap& map, const StencilT& stencil)
1071  {
1072  typedef typename StencilT::ValueType Vec3Type;
1073  typedef typename StencilT::ValueType::value_type ValueType;
1074 
1075  ValueType div(0);
1076  div =ISDivergence<CD_2NDT>::result(stencil);
1077  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1078  return div * inv2dx;
1079  }
1080 };
1081 
1082 
1083 // uniform scale translate 2nd order
1084 template<>
1086 {
1087  // random access version
1088  template<typename Accessor> static typename Accessor::ValueType::value_type
1089  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1090  {
1091  typedef typename Accessor::ValueType Vec3Type;
1092  typedef typename Accessor::ValueType::value_type ValueType;
1093 
1094  ValueType div(0);
1095 
1096  div =ISDivergence<CD_2NDT>::result(grid, ijk);
1097  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1098  return div * inv2dx;
1099  }
1100 
1101  // stencil access version
1102  template<typename StencilT> static typename StencilT::ValueType::value_type
1103  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1104  {
1105  typedef typename StencilT::ValueType Vec3Type;
1106  typedef typename StencilT::ValueType::value_type ValueType;
1107 
1108  ValueType div(0);
1109 
1110  div =ISDivergence<CD_2NDT>::result(stencil);
1111  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1112  return div * inv2dx;
1113  }
1114 };
1115 
1116 
1117 // scale, any scheme
1118 template<DScheme DiffScheme>
1119 struct Divergence<ScaleMap, DiffScheme>
1120 {
1121  // random access version
1122  template<typename Accessor> static typename Accessor::ValueType::value_type
1123  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
1124  {
1125  typedef typename Accessor::ValueType Vec3Type;
1126  typedef typename Accessor::ValueType::value_type ValueType;
1127 
1128  ValueType div = ValueType(
1129  D1Vec<DiffScheme>::inX(grid, ijk, 0) * (map.getInvScale()[0]) +
1130  D1Vec<DiffScheme>::inY(grid, ijk, 1) * (map.getInvScale()[1]) +
1131  D1Vec<DiffScheme>::inZ(grid, ijk, 2) * (map.getInvScale()[2]));
1132  return div;
1133  }
1134 
1135  // stencil access version
1136  template<typename StencilT> static typename StencilT::ValueType::value_type
1137  result(const ScaleMap& map, const StencilT& stencil)
1138  {
1139  typedef typename StencilT::ValueType Vec3Type;
1140  typedef typename StencilT::ValueType::value_type ValueType;
1141 
1142  ValueType div(0);
1143  div = ValueType(
1144  D1Vec<DiffScheme>::inX(stencil, 0) * (map.getInvScale()[0]) +
1145  D1Vec<DiffScheme>::inY(stencil, 1) * (map.getInvScale()[1]) +
1146  D1Vec<DiffScheme>::inZ(stencil, 2) * (map.getInvScale()[2]) );
1147  return div;
1148  }
1149 };
1150 
1151 
1152 // scale translate, any scheme
1153 template<DScheme DiffScheme>
1154 struct Divergence<ScaleTranslateMap, DiffScheme>
1155 {
1156  // random access version
1157  template<typename Accessor> static typename Accessor::ValueType::value_type
1158  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1159  {
1160  typedef typename Accessor::ValueType Vec3Type;
1161  typedef typename Accessor::ValueType::value_type ValueType;
1162 
1163  ValueType div = ValueType(
1164  D1Vec<DiffScheme>::inX(grid, ijk, 0) * (map.getInvScale()[0]) +
1165  D1Vec<DiffScheme>::inY(grid, ijk, 1) * (map.getInvScale()[1]) +
1166  D1Vec<DiffScheme>::inZ(grid, ijk, 2) * (map.getInvScale()[2]));
1167  return div;
1168  }
1169 
1170  // stencil access version
1171  template<typename StencilT> static typename StencilT::ValueType::value_type
1172  result(const ScaleTranslateMap& map, const StencilT& stencil)
1173  {
1174  typedef typename StencilT::ValueType Vec3Type;
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 
1186 
1187 // scale 2nd order
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 Vec3Type;
1196  typedef typename Accessor::ValueType::value_type ValueType;
1197 
1198  ValueType div = ValueType(
1199  D1Vec<CD_2NDT>::inX(grid, ijk, 0) * (map.getInvTwiceScale()[0]) +
1200  D1Vec<CD_2NDT>::inY(grid, ijk, 1) * (map.getInvTwiceScale()[1]) +
1201  D1Vec<CD_2NDT>::inZ(grid, ijk, 2) * (map.getInvTwiceScale()[2]) );
1202  return div;
1203  }
1204 
1205  // stencil access version
1206  template<typename StencilT> static typename StencilT::ValueType::value_type
1207  result(const ScaleMap& map, const StencilT& stencil)
1208  {
1209  typedef typename StencilT::ValueType Vec3Type;
1210  typedef typename StencilT::ValueType::value_type ValueType;
1211 
1212  ValueType div = ValueType(
1213  D1Vec<CD_2NDT>::inX(stencil, 0) * (map.getInvTwiceScale()[0]) +
1214  D1Vec<CD_2NDT>::inY(stencil, 1) * (map.getInvTwiceScale()[1]) +
1215  D1Vec<CD_2NDT>::inZ(stencil, 2) * (map.getInvTwiceScale()[2]) );
1216  return div;
1217  }
1218 };
1219 
1220 
1221 // scale and translate, 2nd order
1222 template<>
1224 {
1225  // random access version
1226  template<typename Accessor> static typename Accessor::ValueType::value_type
1227  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1228  {
1229  typedef typename Accessor::ValueType Vec3Type;
1230  typedef typename Accessor::ValueType::value_type ValueType;
1231 
1232  ValueType div = ValueType(
1233  D1Vec<CD_2NDT>::inX(grid, ijk, 0) * (map.getInvTwiceScale()[0]) +
1234  D1Vec<CD_2NDT>::inY(grid, ijk, 1) * (map.getInvTwiceScale()[1]) +
1235  D1Vec<CD_2NDT>::inZ(grid, ijk, 2) * (map.getInvTwiceScale()[2]) );
1236  return div;
1237  }
1238 
1239  // stencil access version
1240  template<typename StencilT> static typename StencilT::ValueType::value_type
1241  result(const ScaleTranslateMap& map, const StencilT& stencil)
1242  {
1243  typedef typename StencilT::ValueType Vec3Type;
1244  typedef typename StencilT::ValueType::value_type ValueType;
1245 
1246  ValueType div = ValueType(
1247  D1Vec<CD_2NDT>::inX(stencil, 0) * (map.getInvTwiceScale()[0]) +
1248  D1Vec<CD_2NDT>::inY(stencil, 1) * (map.getInvTwiceScale()[1]) +
1249  D1Vec<CD_2NDT>::inZ(stencil, 2) * (map.getInvTwiceScale()[2]) );
1250  return div;
1251  }
1252 };
1254 
1255 
1257 template<typename MapType, DScheme DiffScheme>
1260 struct Curl
1261 {
1262  // random access version
1263  template<typename Accessor> static typename Accessor::ValueType
1264  result(const MapType& map, const Accessor& grid, const Coord& ijk)
1265  {
1266  typedef typename Accessor::ValueType Vec3Type;
1267  Vec3Type mat[3];
1268  for (int i = 0; i < 3; i++) {
1269  Vec3d vec(
1270  D1Vec<DiffScheme>::inX(grid, ijk, i),
1271  D1Vec<DiffScheme>::inY(grid, ijk, i),
1272  D1Vec<DiffScheme>::inZ(grid, ijk, i));
1273  // dF_i/dx_j (x_1 = x, x_2 = y, x_3 = z)
1274  mat[i] = Vec3Type(map.applyIJT(vec, ijk.asVec3d()));
1275  }
1276  return Vec3Type(mat[2][1] - mat[1][2], // dF_3/dx_2 - dF_2/dx_3
1277  mat[0][2] - mat[2][0], // dF_1/dx_3 - dF_3/dx_1
1278  mat[1][0] - mat[0][1]); // dF_2/dx_1 - dF_1/dx_2
1279  }
1280 
1281  // stencil access version
1282  template<typename StencilT> static typename StencilT::ValueType
1283  result(const MapType& map, const StencilT& stencil)
1284  {
1285  typedef typename StencilT::ValueType Vec3Type;
1286  Vec3Type mat[3];
1287  for (int i = 0; i < 3; i++) {
1288  Vec3d vec(
1289  D1Vec<DiffScheme>::inX(stencil, i),
1290  D1Vec<DiffScheme>::inY(stencil, i),
1291  D1Vec<DiffScheme>::inZ(stencil, i));
1292  // dF_i/dx_j (x_1 = x, x_2 = y, x_3 = z)
1293  mat[i] = Vec3Type(map.applyIJT(vec, stencil.getCenterCoord().asVec3d()));
1294  }
1295  return Vec3Type(mat[2][1] - mat[1][2], // dF_3/dx_2 - dF_2/dx_3
1296  mat[0][2] - mat[2][0], // dF_1/dx_3 - dF_3/dx_1
1297  mat[1][0] - mat[0][1]); // dF_2/dx_1 - dF_1/dx_2
1298  }
1299 };
1300 
1301 
1302 template<DScheme DiffScheme>
1303 struct Curl<UniformScaleMap, DiffScheme>
1304 {
1305  // random access version
1306  template<typename Accessor> static typename Accessor::ValueType
1307  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1308  {
1309  typedef typename Accessor::ValueType Vec3Type;
1310  typedef typename Vec3Type::value_type ValueType;
1311  return ISCurl<DiffScheme>::result(grid, ijk) * ValueType(map.getInvScale()[0]);
1312  }
1313 
1314  // Stencil access version
1315  template<typename StencilT> static typename StencilT::ValueType
1316  result(const UniformScaleMap& map, const StencilT& stencil)
1317  {
1318  typedef typename StencilT::ValueType Vec3Type;
1319  typedef typename Vec3Type::value_type ValueType;
1320  return ISCurl<DiffScheme>::result(stencil) * ValueType(map.getInvScale()[0]);
1321  }
1322 };
1323 
1324 
1325 template<DScheme DiffScheme>
1326 struct Curl<UniformScaleTranslateMap, DiffScheme>
1327 {
1328  // random access version
1329  template<typename Accessor> static typename Accessor::ValueType
1330  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1331  {
1332  typedef typename Accessor::ValueType Vec3Type;
1333  typedef typename Vec3Type::value_type ValueType;
1334 
1335  return ISCurl<DiffScheme>::result(grid, ijk) * ValueType(map.getInvScale()[0]);
1336  }
1337 
1338  // stencil access version
1339  template<typename StencilT> static typename StencilT::ValueType
1340  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1341  {
1342  typedef typename StencilT::ValueType Vec3Type;
1343  typedef typename Vec3Type::value_type ValueType;
1344 
1345  return ISCurl<DiffScheme>::result(stencil) * ValueType(map.getInvScale()[0]);
1346  }
1347 };
1348 
1349 template<>
1351 {
1352  // random access version
1353  template<typename Accessor> static typename Accessor::ValueType
1354  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1355  {
1356  typedef typename Accessor::ValueType Vec3Type;
1357  typedef typename Vec3Type::value_type ValueType;
1358 
1359  return ISCurl<CD_2NDT>::result(grid, ijk) * ValueType(map.getInvTwiceScale()[0]);
1360  }
1361 
1362  // stencil access version
1363  template<typename StencilT> static typename StencilT::ValueType
1364  result(const UniformScaleMap& map, const StencilT& stencil)
1365  {
1366  typedef typename StencilT::ValueType Vec3Type;
1367  typedef typename Vec3Type::value_type ValueType;
1368 
1369  return ISCurl<CD_2NDT>::result(stencil) * ValueType(map.getInvTwiceScale()[0]);
1370  }
1371 };
1372 
1373 template<>
1375 {
1376  // random access version
1377  template<typename Accessor> static typename Accessor::ValueType
1378  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1379  {
1380  typedef typename Accessor::ValueType Vec3Type;
1381  typedef typename Vec3Type::value_type ValueType;
1382 
1383  return ISCurl<CD_2NDT>::result(grid, ijk) * ValueType(map.getInvTwiceScale()[0]);
1384  }
1385 
1386  // stencil access version
1387  template<typename StencilT> static typename StencilT::ValueType
1388  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1389  {
1390  typedef typename StencilT::ValueType Vec3Type;
1391  typedef typename Vec3Type::value_type ValueType;
1392 
1393  return ISCurl<CD_2NDT>::result(stencil) * ValueType(map.getInvTwiceScale()[0]);
1394  }
1395 };
1397 
1398 
1400 template<typename MapType, DDScheme DiffScheme>
1404 {
1405  // random access version
1406  template<typename Accessor>
1407  static typename Accessor::ValueType result(const MapType& map,
1408  const Accessor& grid, const Coord& ijk)
1409  {
1410  typedef typename Accessor::ValueType ValueType;
1411  // all the second derivatives in index space
1412  ValueType iddx = D2<DiffScheme>::inX(grid, ijk);
1413  ValueType iddy = D2<DiffScheme>::inY(grid, ijk);
1414  ValueType iddz = D2<DiffScheme>::inZ(grid, ijk);
1415 
1416  ValueType iddxy = D2<DiffScheme>::inXandY(grid, ijk);
1417  ValueType iddyz = D2<DiffScheme>::inYandZ(grid, ijk);
1418  ValueType iddxz = D2<DiffScheme>::inXandZ(grid, ijk);
1419 
1420  // second derivatives in index space
1421  Mat3d d2_is(iddx, iddxy, iddxz,
1422  iddxy, iddy, iddyz,
1423  iddxz, iddyz, iddz);
1424 
1425  Mat3d d2_rs; // to hold the second derivative matrix in range space
1427  d2_rs = map.applyIJC(d2_is);
1428  } else {
1429  // compute the first derivatives with 2nd order accuracy.
1430  Vec3d d1_is(D1<CD_2ND>::inX(grid, ijk),
1431  D1<CD_2ND>::inY(grid, ijk),
1432  D1<CD_2ND>::inZ(grid, ijk) );
1433 
1434  d2_rs = map.applyIJC(d2_is, d1_is, ijk.asVec3d());
1435  }
1436 
1437  // the trace of the second derivative (range space) matrix is laplacian
1438  return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2));
1439  }
1440 
1441  // stencil access version
1442  template<typename StencilT>
1443  static typename StencilT::ValueType result(const MapType& map, const StencilT& stencil)
1444  {
1445  typedef typename StencilT::ValueType ValueType;
1446  // all the second derivatives in index space
1447  ValueType iddx = D2<DiffScheme>::inX(stencil);
1448  ValueType iddy = D2<DiffScheme>::inY(stencil);
1449  ValueType iddz = D2<DiffScheme>::inZ(stencil);
1450 
1451  ValueType iddxy = D2<DiffScheme>::inXandY(stencil);
1452  ValueType iddyz = D2<DiffScheme>::inYandZ(stencil);
1453  ValueType iddxz = D2<DiffScheme>::inXandZ(stencil);
1454 
1455  // second derivatives in index space
1456  Mat3d d2_is(iddx, iddxy, iddxz,
1457  iddxy, iddy, iddyz,
1458  iddxz, iddyz, iddz);
1459 
1460  Mat3d d2_rs; // to hold the second derivative matrix in range space
1462  d2_rs = map.applyIJC(d2_is);
1463  } else {
1464  // compute the first derivatives with 2nd order accuracy.
1465  Vec3d d1_is(D1<CD_2ND>::inX(stencil),
1466  D1<CD_2ND>::inY(stencil),
1467  D1<CD_2ND>::inZ(stencil) );
1468 
1469  d2_rs = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d());
1470  }
1471 
1472  // the trace of the second derivative (range space) matrix is laplacian
1473  return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2));
1474  }
1475 };
1476 
1477 
1478 template<DDScheme DiffScheme>
1479 struct Laplacian<TranslationMap, DiffScheme>
1480 {
1481  // random access version
1482  template<typename Accessor>
1483  static typename Accessor::ValueType result(const TranslationMap&,
1484  const Accessor& grid, const Coord& ijk)
1485  {
1486  return ISLaplacian<DiffScheme>::result(grid, ijk);
1487  }
1488 
1489  // stencil access version
1490  template<typename StencilT>
1491  static typename StencilT::ValueType result(const TranslationMap&, const StencilT& stencil)
1492  {
1493  return ISLaplacian<DiffScheme>::result(stencil);
1494  }
1495 };
1496 
1497 
1498 // The Laplacian is invariant to rotation or reflection.
1499 template<DDScheme DiffScheme>
1500 struct Laplacian<UnitaryMap, DiffScheme>
1501 {
1502  // random access version
1503  template<typename Accessor>
1504  static typename Accessor::ValueType result(const UnitaryMap&,
1505  const Accessor& grid, const Coord& ijk)
1506  {
1507  return ISLaplacian<DiffScheme>::result(grid, ijk);
1508  }
1509 
1510  // stencil access version
1511  template<typename StencilT>
1512  static typename StencilT::ValueType result(const UnitaryMap&, const StencilT& stencil)
1513  {
1514  return ISLaplacian<DiffScheme>::result(stencil);
1515  }
1516 };
1517 
1518 
1519 template<DDScheme DiffScheme>
1520 struct Laplacian<UniformScaleMap, DiffScheme>
1521 {
1522  // random access version
1523  template<typename Accessor> static typename Accessor::ValueType
1524  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1525  {
1526  typedef typename Accessor::ValueType ValueType;
1527  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1528  return ISLaplacian<DiffScheme>::result(grid, ijk) * invdxdx;
1529  }
1530 
1531  // stencil access version
1532  template<typename StencilT> static typename StencilT::ValueType
1533  result(const UniformScaleMap& map, const StencilT& stencil)
1534  {
1535  typedef typename StencilT::ValueType ValueType;
1536  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1537  return ISLaplacian<DiffScheme>::result(stencil) * invdxdx;
1538  }
1539 };
1540 
1541 
1542 template<DDScheme DiffScheme>
1544 {
1545  // random access version
1546  template<typename Accessor> static typename Accessor::ValueType
1547  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1548  {
1549  typedef typename Accessor::ValueType ValueType;
1550  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1551  return ISLaplacian<DiffScheme>::result(grid, ijk) * invdxdx;
1552  }
1553 
1554  // stencil access version
1555  template<typename StencilT> static typename StencilT::ValueType
1556  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1557  {
1558  typedef typename StencilT::ValueType ValueType;
1559  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1560  return ISLaplacian<DiffScheme>::result(stencil) * invdxdx;
1561  }
1562 };
1563 
1564 
1565 template<DDScheme DiffScheme>
1566 struct Laplacian<ScaleMap, DiffScheme>
1567 {
1568  // random access version
1569  template<typename Accessor> static typename Accessor::ValueType
1570  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
1571  {
1572  typedef typename Accessor::ValueType ValueType;
1573 
1574  // compute the second derivatives in index space
1575  ValueType iddx = D2<DiffScheme>::inX(grid, ijk);
1576  ValueType iddy = D2<DiffScheme>::inY(grid, ijk);
1577  ValueType iddz = D2<DiffScheme>::inZ(grid, ijk);
1578  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1579  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1580  return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1581  }
1582 
1583  // stencil access version
1584  template<typename StencilT> static typename StencilT::ValueType
1585  result(const ScaleMap& map, const StencilT& stencil)
1586  {
1587  typedef typename StencilT::ValueType ValueType;
1588 
1589  // compute the second derivatives in index space
1590  ValueType iddx = D2<DiffScheme>::inX(stencil);
1591  ValueType iddy = D2<DiffScheme>::inY(stencil);
1592  ValueType iddz = D2<DiffScheme>::inZ(stencil);
1593  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1594  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1595  return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1596  }
1597 };
1598 
1599 
1600 template<DDScheme DiffScheme>
1601 struct Laplacian<ScaleTranslateMap, DiffScheme>
1602 {
1603  // random access version
1604  template<typename Accessor> static typename Accessor::ValueType
1605  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1606  {
1607  typedef typename Accessor::ValueType ValueType;
1608  // compute the second derivatives in index space
1609  ValueType iddx = D2<DiffScheme>::inX(grid, ijk);
1610  ValueType iddy = D2<DiffScheme>::inY(grid, ijk);
1611  ValueType iddz = D2<DiffScheme>::inZ(grid, ijk);
1612  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1613  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1614  return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1615  }
1616 
1617  // stencil access version
1618  template<typename StencilT> static typename StencilT::ValueType
1619  result(const ScaleTranslateMap& map, const StencilT& stencil)
1620  {
1621  typedef typename StencilT::ValueType ValueType;
1622  // compute the second derivatives in index space
1623  ValueType iddx = D2<DiffScheme>::inX(stencil);
1624  ValueType iddy = D2<DiffScheme>::inY(stencil);
1625  ValueType iddz = D2<DiffScheme>::inZ(stencil);
1626  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1627  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1628  return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1629  }
1630 };
1631 
1632 
1636 template<typename MapType, DScheme DiffScheme>
1637 struct CPT
1638 {
1639  // random access version
1640  template<typename Accessor> static math::Vec3<typename Accessor::ValueType>
1641  result(const MapType& map, const Accessor& grid, const Coord& ijk)
1642  {
1643  typedef typename Accessor::ValueType ValueType;
1644  typedef Vec3<ValueType> Vec3Type;
1645 
1646  // current distance
1647  ValueType d = grid.getValue(ijk);
1648  // compute gradient in physical space where it is a unit normal
1649  // since the grid holds a distance level set.
1650  Vec3d vectorFromSurface(d*Gradient<MapType,DiffScheme>::result(map, grid, ijk));
1652  Vec3d result = ijk.asVec3d() - map.applyInverseMap(vectorFromSurface);
1653  return Vec3Type(result);
1654  } else {
1655  Vec3d location = map.applyMap(ijk.asVec3d());
1656  Vec3d result = map.applyInverseMap(location - vectorFromSurface);
1657  return Vec3Type(result);
1658  }
1659  }
1660 
1661  // stencil access version
1662  template<typename StencilT> static math::Vec3<typename StencilT::ValueType>
1663  result(const MapType& map, const StencilT& stencil)
1664  {
1665  typedef typename StencilT::ValueType ValueType;
1666  typedef Vec3<ValueType> Vec3Type;
1667 
1668  // current distance
1669  ValueType d = stencil.template getValue<0, 0, 0>();
1670  // compute gradient in physical space where it is a unit normal
1671  // since the grid holds a distance level set.
1672  Vec3d vectorFromSurface(d*Gradient<MapType, DiffScheme>::result(map, stencil));
1674  Vec3d result = stencil.getCenterCoord().asVec3d()
1675  - map.applyInverseMap(vectorFromSurface);
1676  return Vec3Type(result);
1677  } else {
1678  Vec3d location = map.applyMap(stencil.getCenterCoord().asVec3d());
1679  Vec3d result = map.applyInverseMap(location - vectorFromSurface);
1680  return Vec3Type(result);
1681  }
1682  }
1683 };
1684 
1685 
1689 template<typename MapType, DScheme DiffScheme>
1691 {
1692  // random access version
1693  template<typename Accessor> static Vec3<typename Accessor::ValueType>
1694  result(const MapType& map, const Accessor& grid, const Coord& ijk)
1695  {
1696  typedef typename Accessor::ValueType ValueType;
1697  typedef Vec3<ValueType> Vec3Type;
1698  // current distance
1699  ValueType d = grid.getValue(ijk);
1700  // compute gradient in physical space where it is a unit normal
1701  // since the grid holds a distance level set.
1702  Vec3Type vectorFromSurface =
1703  d*Gradient<MapType,DiffScheme>::result(map, grid, ijk);
1704  Vec3d result = map.applyMap(ijk.asVec3d()) - vectorFromSurface;
1705 
1706  return Vec3Type(result);
1707  }
1708 
1709  // stencil access version
1710  template<typename StencilT> static Vec3<typename StencilT::ValueType>
1711  result(const MapType& map, const StencilT& stencil)
1712  {
1713  typedef typename StencilT::ValueType ValueType;
1714  typedef Vec3<ValueType> Vec3Type;
1715  // current distance
1716  ValueType d = stencil.template getValue<0, 0, 0>();
1717  // compute gradient in physical space where it is a unit normal
1718  // since the grid holds a distance level set.
1719  Vec3Type vectorFromSurface =
1720  d*Gradient<MapType, DiffScheme>::result(map, stencil);
1721  Vec3d result = map.applyMap(stencil.getCenterCoord().asVec3d()) - vectorFromSurface;
1722 
1723  return Vec3Type(result);
1724  }
1725 };
1726 
1727 
1731 template<typename MapType, DDScheme DiffScheme2, DScheme DiffScheme1>
1733 {
1734  // random access version
1735  template<typename Accessor>
1736  static void compute(const MapType& map, const Accessor& grid, const Coord& ijk,
1737  double& alpha, double& beta)
1738  {
1739  typedef typename Accessor::ValueType ValueType;
1740  typedef Vec3<ValueType> Vec3Type;
1741 
1742  // compute the gradient in index space
1743  Vec3d d1_is(D1<DiffScheme1>::inX(grid, ijk),
1744  D1<DiffScheme1>::inY(grid, ijk),
1745  D1<DiffScheme1>::inZ(grid, ijk) );
1746 
1747  // all the second derivatives in index space
1748  ValueType iddx = D2<DiffScheme2>::inX(grid, ijk);
1749  ValueType iddy = D2<DiffScheme2>::inY(grid, ijk);
1750  ValueType iddz = D2<DiffScheme2>::inZ(grid, ijk);
1751 
1752  ValueType iddxy = D2<DiffScheme2>::inXandY(grid, ijk);
1753  ValueType iddyz = D2<DiffScheme2>::inYandZ(grid, ijk);
1754  ValueType iddxz = D2<DiffScheme2>::inXandZ(grid, ijk);
1755 
1756  // second derivatives in index space
1757  Mat3d d2_is(iddx, iddxy, iddxz,
1758  iddxy, iddy, iddyz,
1759  iddxz, iddyz, iddz);
1760 
1761  // convert to range space
1762  Mat3d d2_rs;
1763  Vec3d d1_rs;
1765  d2_rs = map.applyIJC(d2_is);
1766  d1_rs = map.applyIJT(d1_is);
1767  } else {
1768  d2_rs = map.applyIJC(d2_is, d1_is, ijk.asVec3d());
1769  d1_rs = map.applyIJT(d1_is, ijk.asVec3d());
1770  }
1771 
1772  // assemble the mean curvature
1773  double Dx2 = d1_rs(0)*d1_rs(0);
1774  double Dy2 = d1_rs(1)*d1_rs(1);
1775  double Dz2 = d1_rs(2)*d1_rs(2);
1776 
1777  // for return
1778  alpha = (Dx2*(d2_rs(1,1)+d2_rs(2,2))+Dy2*(d2_rs(0,0)+d2_rs(2,2))
1779  +Dz2*(d2_rs(0,0)+d2_rs(1,1))
1780  -2*(d1_rs(0)*(d1_rs(1)*d2_rs(0,1)+d1_rs(2)*d2_rs(0,2))
1781  +d1_rs(1)*d1_rs(2)*d2_rs(1,2)));
1782  beta = std::sqrt(Dx2 + Dy2 + Dz2); // * 1/dx
1783  }
1784 
1785  template<typename Accessor>
1786  static typename Accessor::ValueType result(const MapType& map,
1787  const Accessor& grid, const Coord& ijk)
1788  {
1789  typedef typename Accessor::ValueType ValueType;
1790  double alpha, beta;
1791  compute(map, grid, ijk, alpha, beta);
1792 
1793  return ValueType(alpha/(2. *math::Pow3(beta)));
1794  }
1795 
1796  template<typename Accessor>
1797  static typename Accessor::ValueType normGrad(const MapType& map,
1798  const Accessor& grid, const Coord& ijk)
1799  {
1800  typedef typename Accessor::ValueType ValueType;
1801  double alpha, beta;
1802  compute(map, grid, ijk, alpha, beta);
1803 
1804  return ValueType(alpha/(2. *math::Pow2(beta)));
1805  }
1806 
1807  // stencil access version
1808  template<typename StencilT>
1809  static void compute(const MapType& map, const StencilT& stencil, double& alpha, double& beta)
1810  {
1811  typedef typename StencilT::ValueType ValueType;
1812  typedef Vec3<ValueType> Vec3Type;
1813 
1814  // compute the gradient in index space
1815  Vec3d d1_is(D1<DiffScheme1>::inX(stencil),
1816  D1<DiffScheme1>::inY(stencil),
1817  D1<DiffScheme1>::inZ(stencil) );
1818 
1819  // all the second derivatives in index space
1820  ValueType iddx = D2<DiffScheme2>::inX(stencil);
1821  ValueType iddy = D2<DiffScheme2>::inY(stencil);
1822  ValueType iddz = D2<DiffScheme2>::inZ(stencil);
1823 
1824  ValueType iddxy = D2<DiffScheme2>::inXandY(stencil);
1825  ValueType iddyz = D2<DiffScheme2>::inYandZ(stencil);
1826  ValueType iddxz = D2<DiffScheme2>::inXandZ(stencil);
1827 
1828  // second derivatives in index space
1829  Mat3d d2_is(iddx, iddxy, iddxz,
1830  iddxy, iddy, iddyz,
1831  iddxz, iddyz, iddz);
1832 
1833  // convert to range space
1834  Mat3d d2_rs;
1835  Vec3d d1_rs;
1837  d2_rs = map.applyIJC(d2_is);
1838  d1_rs = map.applyIJT(d1_is);
1839  } else {
1840  d2_rs = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d());
1841  d1_rs = map.applyIJT(d1_is, stencil.getCenterCoord().asVec3d());
1842  }
1843 
1844  // assemble the mean curvature
1845  double Dx2 = d1_rs(0)*d1_rs(0);
1846  double Dy2 = d1_rs(1)*d1_rs(1);
1847  double Dz2 = d1_rs(2)*d1_rs(2);
1848 
1849  // for return
1850  alpha = (Dx2*(d2_rs(1,1)+d2_rs(2,2))+Dy2*(d2_rs(0,0)+d2_rs(2,2))
1851  +Dz2*(d2_rs(0,0)+d2_rs(1,1))
1852  -2*(d1_rs(0)*(d1_rs(1)*d2_rs(0,1)+d1_rs(2)*d2_rs(0,2))
1853  +d1_rs(1)*d1_rs(2)*d2_rs(1,2)));
1854  beta = std::sqrt(Dx2 + Dy2 + Dz2); // * 1/dx
1855  }
1856 
1857  template<typename StencilT>
1858  static typename StencilT::ValueType
1859  result(const MapType& map, const StencilT stencil)
1860  {
1861  typedef typename StencilT::ValueType ValueType;
1862  double alpha, beta;
1863  compute(map, stencil, alpha, beta);
1864  return ValueType(alpha/(2*math::Pow3(beta)));
1865  }
1866 
1867  template<typename StencilT>
1868  static typename StencilT::ValueType normGrad(const MapType& map, const StencilT stencil)
1869  {
1870  typedef typename StencilT::ValueType ValueType;
1871  double alpha, beta;
1872  compute(map, stencil, alpha, beta);
1873 
1874  return ValueType(alpha/(2*math::Pow2(beta)));
1875  }
1876 };
1877 
1878 
1879 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1880 struct MeanCurvature<TranslationMap, DiffScheme2, DiffScheme1>
1881 {
1882  // random access version
1883  template<typename Accessor>
1884  static typename Accessor::ValueType result(const TranslationMap&,
1885  const Accessor& grid, const Coord& ijk)
1886  {
1887  typedef typename Accessor::ValueType ValueType;
1888 
1889  ValueType alpha, beta;
1890  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
1891 
1892  return ValueType(alpha /(2*math::Pow3(beta)));
1893  }
1894 
1895  template<typename Accessor>
1896  static typename Accessor::ValueType normGrad(const TranslationMap&,
1897  const Accessor& grid, const Coord& ijk)
1898  {
1899  typedef typename Accessor::ValueType ValueType;
1900 
1901  ValueType alpha, beta;
1902  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
1903 
1904  return ValueType(alpha/(2*math::Pow2(beta)));
1905  }
1906 
1907  // stencil access version
1908  template<typename StencilT>
1909  static typename StencilT::ValueType result(const TranslationMap&, const StencilT& stencil)
1910  {
1911  typedef typename StencilT::ValueType ValueType;
1912 
1913  ValueType alpha, beta;
1915 
1916  return ValueType(alpha /(2*math::Pow3(beta)));
1917  }
1918 
1919  template<typename StencilT>
1920  static typename StencilT::ValueType normGrad(const TranslationMap&, const StencilT& stencil)
1921  {
1922  typedef typename StencilT::ValueType ValueType;
1923 
1924  ValueType alpha, beta;
1926 
1927  return ValueType(alpha/(2*math::Pow2(beta)));
1928  }
1929 };
1930 
1931 
1932 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1933 struct MeanCurvature<UniformScaleMap, DiffScheme2, DiffScheme1>
1934 {
1935  // random access version
1936  template<typename Accessor>
1937  static typename Accessor::ValueType result(const UniformScaleMap& map,
1938  const Accessor& grid, const Coord& ijk)
1939  {
1940  typedef typename Accessor::ValueType ValueType;
1941 
1942  ValueType alpha, beta;
1943  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
1944  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1945 
1946  return ValueType(alpha*inv2dx/math::Pow3(beta));
1947  }
1948 
1949  template<typename Accessor>
1950  static typename Accessor::ValueType normGrad(const UniformScaleMap& map,
1951  const Accessor& grid, const Coord& ijk)
1952  {
1953  typedef typename Accessor::ValueType ValueType;
1954 
1955  ValueType alpha, beta;
1956  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
1957  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1958 
1959  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
1960  }
1961 
1962  // stencil access version
1963  template<typename StencilT>
1964  static typename StencilT::ValueType result(const UniformScaleMap& map, const StencilT& stencil)
1965  {
1966  typedef typename StencilT::ValueType ValueType;
1967 
1968  ValueType alpha, beta;
1970  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1971 
1972  return ValueType(alpha*inv2dx/math::Pow3(beta));
1973  }
1974 
1975  template<typename StencilT>
1976  static typename StencilT::ValueType normGrad(const UniformScaleMap& map, const StencilT& stencil)
1977  {
1978  typedef typename StencilT::ValueType ValueType;
1979 
1980  ValueType alpha, beta;
1982  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1983 
1984  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
1985  }
1986 };
1987 
1988 
1989 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1990 struct MeanCurvature<UniformScaleTranslateMap, DiffScheme2, DiffScheme1>
1991 {
1992  // random access version
1993  template<typename Accessor> static typename Accessor::ValueType
1994  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1995  {
1996  typedef typename Accessor::ValueType ValueType;
1997 
1998  ValueType alpha, beta;
1999  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
2000  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
2001 
2002  return ValueType(alpha*inv2dx/math::Pow3(beta));
2003  }
2004 
2005  template<typename Accessor> static typename Accessor::ValueType
2006  normGrad(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
2007  {
2008  typedef typename Accessor::ValueType ValueType;
2009 
2010  ValueType alpha, beta;
2011  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
2012  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
2013 
2014  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
2015  }
2016 
2017  // stencil access version
2018  template<typename StencilT> static typename StencilT::ValueType
2019  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
2020  {
2021  typedef typename StencilT::ValueType ValueType;
2022 
2023  ValueType alpha, beta;
2025  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
2026 
2027  return ValueType(alpha*inv2dx/math::Pow3(beta));
2028  }
2029 
2030  template<typename StencilT> static typename StencilT::ValueType
2031  normGrad(const UniformScaleTranslateMap& map, const StencilT& stencil)
2032  {
2033  typedef typename StencilT::ValueType ValueType;
2034 
2035  ValueType alpha, beta;
2037  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
2038 
2039  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
2040  }
2041 };
2042 
2043 
2049 {
2050 public:
2051  template<typename GridType>
2052  GenericMap(const GridType& g): mMap(g.transform().baseMap()) {}
2053 
2054  GenericMap(const Transform& t): mMap(t.baseMap()) {}
2055  GenericMap(MapBase::Ptr map): mMap(boost::const_pointer_cast<const MapBase>(map)) {}
2058 
2059  Vec3d applyMap(const Vec3d& in) const { return mMap->applyMap(in); }
2060  Vec3d applyInverseMap(const Vec3d& in) const { return mMap->applyInverseMap(in); }
2061 
2062  Vec3d applyIJT(const Vec3d& in) const { return mMap->applyIJT(in); }
2063  Vec3d applyIJT(const Vec3d& in, const Vec3d& pos) const { return mMap->applyIJT(in, pos); }
2064  Mat3d applyIJC(const Mat3d& m) const { return mMap->applyIJC(m); }
2065  Mat3d applyIJC(const Mat3d& m, const Vec3d& v, const Vec3d& pos) const
2066  { return mMap->applyIJC(m,v,pos); }
2067 
2068  double determinant() const { return mMap->determinant(); }
2069  double determinant(const Vec3d& in) const { return mMap->determinant(in); }
2070 
2071  Vec3d voxelSize() const { return mMap->voxelSize(); }
2072  Vec3d voxelSize(const Vec3d&v) const { return mMap->voxelSize(v); }
2073 
2074 private:
2076 };
2077 
2078 } // end math namespace
2079 } // namespace OPENVDB_VERSION_NAME
2080 } // end openvdb namespace
2081 
2082 #endif // OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED
2083 
2084 // Copyright (c) 2012-2013 DreamWorks Animation LLC
2085 // All rights reserved. This software is distributed under the
2086 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )