OpenVDB  2.1.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 
601 
602 // --- Operators defined in the Range of a given map
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 // Partial template specialization of Gradient
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 
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 
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 
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 
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 
824 
825 // Computes |Grad[Phi]| using upwinding
826 template<typename MapType, BiasedGradientScheme GradScheme>
828 {
831 
832 
833  // random access version
834  template<typename Accessor>
835  static typename Accessor::ValueType
836  result(const MapType& map, const Accessor& grid, const Coord& ijk)
837  {
838  typedef typename Accessor::ValueType ValueType;
839  typedef math::Vec3<ValueType> Vec3Type;
840 
841  Vec3Type up = Gradient<MapType, FD>::result(map, grid, ijk);
842  Vec3Type down = Gradient<MapType, BD>::result(map, grid, ijk);
843  return math::GudonovsNormSqrd(grid.getValue(ijk)>0, down, up);
844  }
845 
846  // stencil access version
847  template<typename StencilT>
848  static typename StencilT::ValueType
849  result(const MapType& map, const StencilT& stencil)
850  {
851  typedef typename StencilT::ValueType ValueType;
852  typedef math::Vec3<ValueType> Vec3Type;
853 
854  Vec3Type up = Gradient<MapType, FD>::result(map, stencil);
855  Vec3Type down = Gradient<MapType, BD>::result(map, stencil);
856  return math::GudonovsNormSqrd(stencil.template getValue<0, 0, 0>()>0, down, up);
857  }
858 };
859 
861 template<BiasedGradientScheme GradScheme>
863 {
864  // random access version
865  template<typename Accessor>
866  static typename Accessor::ValueType
867  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
868  {
869  typedef typename Accessor::ValueType ValueType;
870 
871  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
872  return invdxdx * ISGradientNormSqrd<GradScheme>::result(grid, ijk);
873  }
874 
875  // stencil access version
876  template<typename StencilT>
877  static typename StencilT::ValueType
878  result(const UniformScaleMap& map, const StencilT& stencil)
879  {
880  typedef typename StencilT::ValueType ValueType;
881 
882  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
883  return invdxdx * ISGradientNormSqrd<GradScheme>::result(stencil);
884  }
885 };
886 
888 template<BiasedGradientScheme GradScheme>
890 {
891  // random access version
892  template<typename Accessor>
893  static typename Accessor::ValueType
894  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
895  {
896  typedef typename Accessor::ValueType ValueType;
897 
898  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
899  return invdxdx * ISGradientNormSqrd<GradScheme>::result(grid, ijk);
900  }
901 
902  // stencil access version
903  template<typename StencilT>
904  static typename StencilT::ValueType
905  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
906  {
907  typedef typename StencilT::ValueType ValueType;
908 
909  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
910  return invdxdx * ISGradientNormSqrd<GradScheme>::result(stencil);
911  }
912 };
913 
914 
916 template<typename MapType, DScheme DiffScheme>
920 {
921  // random access version
922  template<typename Accessor> static typename Accessor::ValueType::value_type
923  result(const MapType& map, const Accessor& grid, const Coord& ijk)
924  {
925  typedef typename Accessor::ValueType Vec3Type;
926  typedef typename Accessor::ValueType::value_type ValueType;
927 
928  ValueType div(0);
929  for (int i=0; i < 3; i++) {
930  Vec3d vec( D1Vec<DiffScheme>::inX(grid, ijk, i),
931  D1Vec<DiffScheme>::inY(grid, ijk, i),
932  D1Vec<DiffScheme>::inZ(grid, ijk, i) );
933  div += ValueType(map.applyIJT(vec, ijk.asVec3d())[i]);
934  }
935  return div;
936  }
937 
938  // stencil access version
939  template<typename StencilT> static typename StencilT::ValueType::value_type
940  result(const MapType& map, const StencilT& stencil)
941  {
942  typedef typename StencilT::ValueType Vec3Type;
943  typedef typename StencilT::ValueType::value_type ValueType;
944 
945  ValueType div(0);
946  for (int i=0; i < 3; i++) {
947  Vec3d vec( D1Vec<DiffScheme>::inX(stencil, i),
948  D1Vec<DiffScheme>::inY(stencil, i),
949  D1Vec<DiffScheme>::inZ(stencil, i) );
950  div += ValueType(map.applyIJT(vec, stencil.getCenterCoord().asVec3d())[i]);
951  }
952  return div;
953  }
954 };
955 
958 template<DScheme DiffScheme>
959 struct Divergence<TranslationMap, DiffScheme>
960 {
961  // random access version
962  template<typename Accessor> static typename Accessor::ValueType::value_type
963  result(const TranslationMap&, const Accessor& grid, const Coord& ijk)
964  {
965  typedef typename Accessor::ValueType Vec3Type;
966  typedef typename Accessor::ValueType::value_type ValueType;
967 
968  ValueType div(0);
969  div =ISDivergence<DiffScheme>::result(grid, ijk);
970  return div;
971  }
972 
973  // stencil access version
974  template<typename StencilT> static typename StencilT::ValueType::value_type
975  result(const TranslationMap&, const StencilT& stencil)
976  {
977  typedef typename StencilT::ValueType Vec3Type;
978  typedef typename StencilT::ValueType::value_type ValueType;
979 
980  ValueType div(0);
981  div =ISDivergence<DiffScheme>::result(stencil);
982  return div;
983  }
984 };
985 
988 template<DScheme DiffScheme>
989 struct Divergence<UniformScaleMap, DiffScheme>
990 {
991  // random access version
992  template<typename Accessor> static typename Accessor::ValueType::value_type
993  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
994  {
995  typedef typename Accessor::ValueType Vec3Type;
996  typedef typename Accessor::ValueType::value_type ValueType;
997 
998  ValueType div(0);
999 
1000  div =ISDivergence<DiffScheme>::result(grid, ijk);
1001  ValueType invdx = ValueType(map.getInvScale()[0]);
1002  return div * invdx;
1003  }
1004 
1005  // stencil access version
1006  template<typename StencilT> static typename StencilT::ValueType::value_type
1007  result(const UniformScaleMap& map, const StencilT& stencil)
1008  {
1009  typedef typename StencilT::ValueType Vec3Type;
1010  typedef typename StencilT::ValueType::value_type ValueType;
1011 
1012  ValueType div(0);
1013 
1014  div =ISDivergence<DiffScheme>::result(stencil);
1015  ValueType invdx = ValueType(map.getInvScale()[0]);
1016  return div * invdx;
1017  }
1018 };
1019 
1022 template<DScheme DiffScheme>
1024 {
1025  // random access version
1026  template<typename Accessor> static typename Accessor::ValueType::value_type
1027  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1028  {
1029  typedef typename Accessor::ValueType Vec3Type;
1030  typedef typename Accessor::ValueType::value_type ValueType;
1031 
1032  ValueType div(0);
1033 
1034  div =ISDivergence<DiffScheme>::result(grid, ijk);
1035  ValueType invdx = ValueType(map.getInvScale()[0]);
1036  return div * invdx;
1037  }
1038 
1039  // stencil access version
1040  template<typename StencilT> static typename StencilT::ValueType::value_type
1041  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1042  {
1043  typedef typename StencilT::ValueType Vec3Type;
1044  typedef typename StencilT::ValueType::value_type ValueType;
1045 
1046  ValueType div(0);
1047 
1048  div =ISDivergence<DiffScheme>::result(stencil);
1049  ValueType invdx = ValueType(map.getInvScale()[0]);
1050  return div * invdx;
1051  }
1052 };
1053 
1056 template<>
1058 {
1059  // random access version
1060  template<typename Accessor> static typename Accessor::ValueType::value_type
1061  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1062  {
1063  typedef typename Accessor::ValueType Vec3Type;
1064  typedef typename Accessor::ValueType::value_type ValueType;
1065 
1066  ValueType div(0);
1067  div =ISDivergence<CD_2NDT>::result(grid, ijk);
1068  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1069  return div * inv2dx;
1070  }
1071 
1072  // stencil access version
1073  template<typename StencilT> static typename StencilT::ValueType::value_type
1074  result(const UniformScaleMap& map, const StencilT& stencil)
1075  {
1076  typedef typename StencilT::ValueType Vec3Type;
1077  typedef typename StencilT::ValueType::value_type ValueType;
1078 
1079  ValueType div(0);
1080  div =ISDivergence<CD_2NDT>::result(stencil);
1081  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1082  return div * inv2dx;
1083  }
1084 };
1085 
1088 template<>
1090 {
1091  // random access version
1092  template<typename Accessor> static typename Accessor::ValueType::value_type
1093  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1094  {
1095  typedef typename Accessor::ValueType Vec3Type;
1096  typedef typename Accessor::ValueType::value_type ValueType;
1097 
1098  ValueType div(0);
1099 
1100  div =ISDivergence<CD_2NDT>::result(grid, ijk);
1101  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1102  return div * inv2dx;
1103  }
1104 
1105  // stencil access version
1106  template<typename StencilT> static typename StencilT::ValueType::value_type
1107  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1108  {
1109  typedef typename StencilT::ValueType Vec3Type;
1110  typedef typename StencilT::ValueType::value_type ValueType;
1111 
1112  ValueType div(0);
1113 
1114  div =ISDivergence<CD_2NDT>::result(stencil);
1115  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1116  return div * inv2dx;
1117  }
1118 };
1119 
1122 template<DScheme DiffScheme>
1123 struct Divergence<ScaleMap, DiffScheme>
1124 {
1125  // random access version
1126  template<typename Accessor> static typename Accessor::ValueType::value_type
1127  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
1128  {
1129  typedef typename Accessor::ValueType Vec3Type;
1130  typedef typename Accessor::ValueType::value_type ValueType;
1131 
1132  ValueType div = ValueType(
1133  D1Vec<DiffScheme>::inX(grid, ijk, 0) * (map.getInvScale()[0]) +
1134  D1Vec<DiffScheme>::inY(grid, ijk, 1) * (map.getInvScale()[1]) +
1135  D1Vec<DiffScheme>::inZ(grid, ijk, 2) * (map.getInvScale()[2]));
1136  return div;
1137  }
1138 
1139  // stencil access version
1140  template<typename StencilT> static typename StencilT::ValueType::value_type
1141  result(const ScaleMap& map, const StencilT& stencil)
1142  {
1143  typedef typename StencilT::ValueType Vec3Type;
1144  typedef typename StencilT::ValueType::value_type ValueType;
1145 
1146  ValueType div(0);
1147  div = ValueType(
1148  D1Vec<DiffScheme>::inX(stencil, 0) * (map.getInvScale()[0]) +
1149  D1Vec<DiffScheme>::inY(stencil, 1) * (map.getInvScale()[1]) +
1150  D1Vec<DiffScheme>::inZ(stencil, 2) * (map.getInvScale()[2]) );
1151  return div;
1152  }
1153 };
1154 
1157 template<DScheme DiffScheme>
1158 struct Divergence<ScaleTranslateMap, DiffScheme>
1159 {
1160  // random access version
1161  template<typename Accessor> static typename Accessor::ValueType::value_type
1162  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1163  {
1164  typedef typename Accessor::ValueType Vec3Type;
1165  typedef typename Accessor::ValueType::value_type ValueType;
1166 
1167  ValueType div = ValueType(
1168  D1Vec<DiffScheme>::inX(grid, ijk, 0) * (map.getInvScale()[0]) +
1169  D1Vec<DiffScheme>::inY(grid, ijk, 1) * (map.getInvScale()[1]) +
1170  D1Vec<DiffScheme>::inZ(grid, ijk, 2) * (map.getInvScale()[2]));
1171  return div;
1172  }
1173 
1174  // stencil access version
1175  template<typename StencilT> static typename StencilT::ValueType::value_type
1176  result(const ScaleTranslateMap& map, const StencilT& stencil)
1177  {
1178  typedef typename StencilT::ValueType Vec3Type;
1179  typedef typename StencilT::ValueType::value_type ValueType;
1180 
1181  ValueType div(0);
1182  div = ValueType(
1183  D1Vec<DiffScheme>::inX(stencil, 0) * (map.getInvScale()[0]) +
1184  D1Vec<DiffScheme>::inY(stencil, 1) * (map.getInvScale()[1]) +
1185  D1Vec<DiffScheme>::inZ(stencil, 2) * (map.getInvScale()[2]) );
1186  return div;
1187  }
1188 };
1189 
1192 template<>
1194 {
1195  // random access version
1196  template<typename Accessor> static typename Accessor::ValueType::value_type
1197  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
1198  {
1199  typedef typename Accessor::ValueType Vec3Type;
1200  typedef typename Accessor::ValueType::value_type ValueType;
1201 
1202  ValueType div = ValueType(
1203  D1Vec<CD_2NDT>::inX(grid, ijk, 0) * (map.getInvTwiceScale()[0]) +
1204  D1Vec<CD_2NDT>::inY(grid, ijk, 1) * (map.getInvTwiceScale()[1]) +
1205  D1Vec<CD_2NDT>::inZ(grid, ijk, 2) * (map.getInvTwiceScale()[2]) );
1206  return div;
1207  }
1208 
1209  // stencil access version
1210  template<typename StencilT> static typename StencilT::ValueType::value_type
1211  result(const ScaleMap& map, const StencilT& stencil)
1212  {
1213  typedef typename StencilT::ValueType Vec3Type;
1214  typedef typename StencilT::ValueType::value_type ValueType;
1215 
1216  ValueType div = ValueType(
1217  D1Vec<CD_2NDT>::inX(stencil, 0) * (map.getInvTwiceScale()[0]) +
1218  D1Vec<CD_2NDT>::inY(stencil, 1) * (map.getInvTwiceScale()[1]) +
1219  D1Vec<CD_2NDT>::inZ(stencil, 2) * (map.getInvTwiceScale()[2]) );
1220  return div;
1221  }
1222 };
1223 
1226 template<>
1228 {
1229  // random access version
1230  template<typename Accessor> static typename Accessor::ValueType::value_type
1231  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1232  {
1233  typedef typename Accessor::ValueType Vec3Type;
1234  typedef typename Accessor::ValueType::value_type ValueType;
1235 
1236  ValueType div = ValueType(
1237  D1Vec<CD_2NDT>::inX(grid, ijk, 0) * (map.getInvTwiceScale()[0]) +
1238  D1Vec<CD_2NDT>::inY(grid, ijk, 1) * (map.getInvTwiceScale()[1]) +
1239  D1Vec<CD_2NDT>::inZ(grid, ijk, 2) * (map.getInvTwiceScale()[2]) );
1240  return div;
1241  }
1242 
1243  // stencil access version
1244  template<typename StencilT> static typename StencilT::ValueType::value_type
1245  result(const ScaleTranslateMap& map, const StencilT& stencil)
1246  {
1247  typedef typename StencilT::ValueType Vec3Type;
1248  typedef typename StencilT::ValueType::value_type ValueType;
1249 
1250  ValueType div = ValueType(
1251  D1Vec<CD_2NDT>::inX(stencil, 0) * (map.getInvTwiceScale()[0]) +
1252  D1Vec<CD_2NDT>::inY(stencil, 1) * (map.getInvTwiceScale()[1]) +
1253  D1Vec<CD_2NDT>::inZ(stencil, 2) * (map.getInvTwiceScale()[2]) );
1254  return div;
1255  }
1256 };
1258 
1259 
1261 template<typename MapType, DScheme DiffScheme>
1264 struct Curl
1265 {
1266  // random access version
1267  template<typename Accessor> static typename Accessor::ValueType
1268  result(const MapType& map, const Accessor& grid, const Coord& ijk)
1269  {
1270  typedef typename Accessor::ValueType Vec3Type;
1271  Vec3Type mat[3];
1272  for (int i = 0; i < 3; i++) {
1273  Vec3d vec(
1274  D1Vec<DiffScheme>::inX(grid, ijk, i),
1275  D1Vec<DiffScheme>::inY(grid, ijk, i),
1276  D1Vec<DiffScheme>::inZ(grid, ijk, i));
1277  // dF_i/dx_j (x_1 = x, x_2 = y, x_3 = z)
1278  mat[i] = Vec3Type(map.applyIJT(vec, ijk.asVec3d()));
1279  }
1280  return Vec3Type(mat[2][1] - mat[1][2], // dF_3/dx_2 - dF_2/dx_3
1281  mat[0][2] - mat[2][0], // dF_1/dx_3 - dF_3/dx_1
1282  mat[1][0] - mat[0][1]); // dF_2/dx_1 - dF_1/dx_2
1283  }
1284 
1285  // stencil access version
1286  template<typename StencilT> static typename StencilT::ValueType
1287  result(const MapType& map, const StencilT& stencil)
1288  {
1289  typedef typename StencilT::ValueType Vec3Type;
1290  Vec3Type mat[3];
1291  for (int i = 0; i < 3; i++) {
1292  Vec3d vec(
1293  D1Vec<DiffScheme>::inX(stencil, i),
1294  D1Vec<DiffScheme>::inY(stencil, i),
1295  D1Vec<DiffScheme>::inZ(stencil, i));
1296  // dF_i/dx_j (x_1 = x, x_2 = y, x_3 = z)
1297  mat[i] = Vec3Type(map.applyIJT(vec, stencil.getCenterCoord().asVec3d()));
1298  }
1299  return Vec3Type(mat[2][1] - mat[1][2], // dF_3/dx_2 - dF_2/dx_3
1300  mat[0][2] - mat[2][0], // dF_1/dx_3 - dF_3/dx_1
1301  mat[1][0] - mat[0][1]); // dF_2/dx_1 - dF_1/dx_2
1302  }
1303 };
1304 
1306 template<DScheme DiffScheme>
1307 struct Curl<UniformScaleMap, DiffScheme>
1308 {
1309  // random access version
1310  template<typename Accessor> static typename Accessor::ValueType
1311  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1312  {
1313  typedef typename Accessor::ValueType Vec3Type;
1314  typedef typename Vec3Type::value_type ValueType;
1315  return ISCurl<DiffScheme>::result(grid, ijk) * ValueType(map.getInvScale()[0]);
1316  }
1317 
1318  // Stencil access version
1319  template<typename StencilT> static typename StencilT::ValueType
1320  result(const UniformScaleMap& map, const StencilT& stencil)
1321  {
1322  typedef typename StencilT::ValueType Vec3Type;
1323  typedef typename Vec3Type::value_type ValueType;
1324  return ISCurl<DiffScheme>::result(stencil) * ValueType(map.getInvScale()[0]);
1325  }
1326 };
1327 
1329 template<DScheme DiffScheme>
1330 struct Curl<UniformScaleTranslateMap, DiffScheme>
1331 {
1332  // random access version
1333  template<typename Accessor> static typename Accessor::ValueType
1334  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1335  {
1336  typedef typename Accessor::ValueType Vec3Type;
1337  typedef typename Vec3Type::value_type ValueType;
1338 
1339  return ISCurl<DiffScheme>::result(grid, ijk) * ValueType(map.getInvScale()[0]);
1340  }
1341 
1342  // stencil access version
1343  template<typename StencilT> static typename StencilT::ValueType
1344  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1345  {
1346  typedef typename StencilT::ValueType Vec3Type;
1347  typedef typename Vec3Type::value_type ValueType;
1348 
1349  return ISCurl<DiffScheme>::result(stencil) * ValueType(map.getInvScale()[0]);
1350  }
1351 };
1352 
1354 template<>
1356 {
1357  // random access version
1358  template<typename Accessor> static typename Accessor::ValueType
1359  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1360  {
1361  typedef typename Accessor::ValueType Vec3Type;
1362  typedef typename Vec3Type::value_type ValueType;
1363 
1364  return ISCurl<CD_2NDT>::result(grid, ijk) * ValueType(map.getInvTwiceScale()[0]);
1365  }
1366 
1367  // stencil access version
1368  template<typename StencilT> static typename StencilT::ValueType
1369  result(const UniformScaleMap& map, const StencilT& stencil)
1370  {
1371  typedef typename StencilT::ValueType Vec3Type;
1372  typedef typename Vec3Type::value_type ValueType;
1373 
1374  return ISCurl<CD_2NDT>::result(stencil) * ValueType(map.getInvTwiceScale()[0]);
1375  }
1376 };
1377 
1379 template<>
1381 {
1382  // random access version
1383  template<typename Accessor> static typename Accessor::ValueType
1384  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1385  {
1386  typedef typename Accessor::ValueType Vec3Type;
1387  typedef typename Vec3Type::value_type ValueType;
1388 
1389  return ISCurl<CD_2NDT>::result(grid, ijk) * ValueType(map.getInvTwiceScale()[0]);
1390  }
1391 
1392  // stencil access version
1393  template<typename StencilT> static typename StencilT::ValueType
1394  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1395  {
1396  typedef typename StencilT::ValueType Vec3Type;
1397  typedef typename Vec3Type::value_type ValueType;
1398 
1399  return ISCurl<CD_2NDT>::result(stencil) * ValueType(map.getInvTwiceScale()[0]);
1400  }
1401 };
1403 
1404 
1406 template<typename MapType, DDScheme DiffScheme>
1410 {
1411  // random access version
1412  template<typename Accessor>
1413  static typename Accessor::ValueType result(const MapType& map,
1414  const Accessor& grid, const Coord& ijk)
1415  {
1416  typedef typename Accessor::ValueType ValueType;
1417  // all the second derivatives in index space
1418  ValueType iddx = D2<DiffScheme>::inX(grid, ijk);
1419  ValueType iddy = D2<DiffScheme>::inY(grid, ijk);
1420  ValueType iddz = D2<DiffScheme>::inZ(grid, ijk);
1421 
1422  ValueType iddxy = D2<DiffScheme>::inXandY(grid, ijk);
1423  ValueType iddyz = D2<DiffScheme>::inYandZ(grid, ijk);
1424  ValueType iddxz = D2<DiffScheme>::inXandZ(grid, ijk);
1425 
1426  // second derivatives in index space
1427  Mat3d d2_is(iddx, iddxy, iddxz,
1428  iddxy, iddy, iddyz,
1429  iddxz, iddyz, iddz);
1430 
1431  Mat3d d2_rs; // to hold the second derivative matrix in range space
1433  d2_rs = map.applyIJC(d2_is);
1434  } else {
1435  // compute the first derivatives with 2nd order accuracy.
1436  Vec3d d1_is(D1<CD_2ND>::inX(grid, ijk),
1437  D1<CD_2ND>::inY(grid, ijk),
1438  D1<CD_2ND>::inZ(grid, ijk) );
1439 
1440  d2_rs = map.applyIJC(d2_is, d1_is, ijk.asVec3d());
1441  }
1442 
1443  // the trace of the second derivative (range space) matrix is laplacian
1444  return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2));
1445  }
1446 
1447  // stencil access version
1448  template<typename StencilT>
1449  static typename StencilT::ValueType result(const MapType& map, const StencilT& stencil)
1450  {
1451  typedef typename StencilT::ValueType ValueType;
1452  // all the second derivatives in index space
1453  ValueType iddx = D2<DiffScheme>::inX(stencil);
1454  ValueType iddy = D2<DiffScheme>::inY(stencil);
1455  ValueType iddz = D2<DiffScheme>::inZ(stencil);
1456 
1457  ValueType iddxy = D2<DiffScheme>::inXandY(stencil);
1458  ValueType iddyz = D2<DiffScheme>::inYandZ(stencil);
1459  ValueType iddxz = D2<DiffScheme>::inXandZ(stencil);
1460 
1461  // second derivatives in index space
1462  Mat3d d2_is(iddx, iddxy, iddxz,
1463  iddxy, iddy, iddyz,
1464  iddxz, iddyz, iddz);
1465 
1466  Mat3d d2_rs; // to hold the second derivative matrix in range space
1468  d2_rs = map.applyIJC(d2_is);
1469  } else {
1470  // compute the first derivatives with 2nd order accuracy.
1471  Vec3d d1_is(D1<CD_2ND>::inX(stencil),
1472  D1<CD_2ND>::inY(stencil),
1473  D1<CD_2ND>::inZ(stencil) );
1474 
1475  d2_rs = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d());
1476  }
1477 
1478  // the trace of the second derivative (range space) matrix is laplacian
1479  return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2));
1480  }
1481 };
1482 
1483 
1484 template<DDScheme DiffScheme>
1485 struct Laplacian<TranslationMap, DiffScheme>
1486 {
1487  // random access version
1488  template<typename Accessor>
1489  static typename Accessor::ValueType result(const TranslationMap&,
1490  const Accessor& grid, const Coord& ijk)
1491  {
1492  return ISLaplacian<DiffScheme>::result(grid, ijk);
1493  }
1494 
1495  // stencil access version
1496  template<typename StencilT>
1497  static typename StencilT::ValueType result(const TranslationMap&, const StencilT& stencil)
1498  {
1499  return ISLaplacian<DiffScheme>::result(stencil);
1500  }
1501 };
1502 
1503 
1504 // The Laplacian is invariant to rotation or reflection.
1505 template<DDScheme DiffScheme>
1506 struct Laplacian<UnitaryMap, DiffScheme>
1507 {
1508  // random access version
1509  template<typename Accessor>
1510  static typename Accessor::ValueType result(const UnitaryMap&,
1511  const Accessor& grid, const Coord& ijk)
1512  {
1513  return ISLaplacian<DiffScheme>::result(grid, ijk);
1514  }
1515 
1516  // stencil access version
1517  template<typename StencilT>
1518  static typename StencilT::ValueType result(const UnitaryMap&, const StencilT& stencil)
1519  {
1520  return ISLaplacian<DiffScheme>::result(stencil);
1521  }
1522 };
1523 
1524 
1525 template<DDScheme DiffScheme>
1526 struct Laplacian<UniformScaleMap, DiffScheme>
1527 {
1528  // random access version
1529  template<typename Accessor> static typename Accessor::ValueType
1530  result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk)
1531  {
1532  typedef typename Accessor::ValueType ValueType;
1533  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1534  return ISLaplacian<DiffScheme>::result(grid, ijk) * invdxdx;
1535  }
1536 
1537  // stencil access version
1538  template<typename StencilT> static typename StencilT::ValueType
1539  result(const UniformScaleMap& map, const StencilT& stencil)
1540  {
1541  typedef typename StencilT::ValueType ValueType;
1542  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1543  return ISLaplacian<DiffScheme>::result(stencil) * invdxdx;
1544  }
1545 };
1546 
1547 
1548 template<DDScheme DiffScheme>
1550 {
1551  // random access version
1552  template<typename Accessor> static typename Accessor::ValueType
1553  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1554  {
1555  typedef typename Accessor::ValueType ValueType;
1556  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1557  return ISLaplacian<DiffScheme>::result(grid, ijk) * invdxdx;
1558  }
1559 
1560  // stencil access version
1561  template<typename StencilT> static typename StencilT::ValueType
1562  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
1563  {
1564  typedef typename StencilT::ValueType ValueType;
1565  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1566  return ISLaplacian<DiffScheme>::result(stencil) * invdxdx;
1567  }
1568 };
1569 
1570 
1571 template<DDScheme DiffScheme>
1572 struct Laplacian<ScaleMap, DiffScheme>
1573 {
1574  // random access version
1575  template<typename Accessor> static typename Accessor::ValueType
1576  result(const ScaleMap& map, const Accessor& grid, const Coord& ijk)
1577  {
1578  typedef typename Accessor::ValueType ValueType;
1579 
1580  // compute the second derivatives in index space
1581  ValueType iddx = D2<DiffScheme>::inX(grid, ijk);
1582  ValueType iddy = D2<DiffScheme>::inY(grid, ijk);
1583  ValueType iddz = D2<DiffScheme>::inZ(grid, ijk);
1584  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1585  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1586  return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1587  }
1588 
1589  // stencil access version
1590  template<typename StencilT> static typename StencilT::ValueType
1591  result(const ScaleMap& map, const StencilT& stencil)
1592  {
1593  typedef typename StencilT::ValueType ValueType;
1594 
1595  // compute the second derivatives in index space
1596  ValueType iddx = D2<DiffScheme>::inX(stencil);
1597  ValueType iddy = D2<DiffScheme>::inY(stencil);
1598  ValueType iddz = D2<DiffScheme>::inZ(stencil);
1599  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1600  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1601  return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1602  }
1603 };
1604 
1605 
1606 template<DDScheme DiffScheme>
1607 struct Laplacian<ScaleTranslateMap, DiffScheme>
1608 {
1609  // random access version
1610  template<typename Accessor> static typename Accessor::ValueType
1611  result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
1612  {
1613  typedef typename Accessor::ValueType ValueType;
1614  // compute the second derivatives in index space
1615  ValueType iddx = D2<DiffScheme>::inX(grid, ijk);
1616  ValueType iddy = D2<DiffScheme>::inY(grid, ijk);
1617  ValueType iddz = D2<DiffScheme>::inZ(grid, ijk);
1618  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1619  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1620  return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1621  }
1622 
1623  // stencil access version
1624  template<typename StencilT> static typename StencilT::ValueType
1625  result(const ScaleTranslateMap& map, const StencilT& stencil)
1626  {
1627  typedef typename StencilT::ValueType ValueType;
1628  // compute the second derivatives in index space
1629  ValueType iddx = D2<DiffScheme>::inX(stencil);
1630  ValueType iddy = D2<DiffScheme>::inY(stencil);
1631  ValueType iddz = D2<DiffScheme>::inZ(stencil);
1632  const Vec3d& invScaleSqr = map.getInvScaleSqr();
1633  // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum
1634  return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1635  }
1636 };
1637 
1638 
1642 template<typename MapType, DScheme DiffScheme>
1643 struct CPT
1644 {
1645  // random access version
1646  template<typename Accessor> static math::Vec3<typename Accessor::ValueType>
1647  result(const MapType& map, const Accessor& grid, const Coord& ijk)
1648  {
1649  typedef typename Accessor::ValueType ValueType;
1650  typedef Vec3<ValueType> Vec3Type;
1651 
1652  // current distance
1653  ValueType d = grid.getValue(ijk);
1654  // compute gradient in physical space where it is a unit normal
1655  // since the grid holds a distance level set.
1656  Vec3d vectorFromSurface(d*Gradient<MapType,DiffScheme>::result(map, grid, ijk));
1658  Vec3d result = ijk.asVec3d() - map.applyInverseMap(vectorFromSurface);
1659  return Vec3Type(result);
1660  } else {
1661  Vec3d location = map.applyMap(ijk.asVec3d());
1662  Vec3d result = map.applyInverseMap(location - vectorFromSurface);
1663  return Vec3Type(result);
1664  }
1665  }
1666 
1667  // stencil access version
1668  template<typename StencilT> static math::Vec3<typename StencilT::ValueType>
1669  result(const MapType& map, const StencilT& stencil)
1670  {
1671  typedef typename StencilT::ValueType ValueType;
1672  typedef Vec3<ValueType> Vec3Type;
1673 
1674  // current distance
1675  ValueType d = stencil.template getValue<0, 0, 0>();
1676  // compute gradient in physical space where it is a unit normal
1677  // since the grid holds a distance level set.
1678  Vec3d vectorFromSurface(d*Gradient<MapType, DiffScheme>::result(map, stencil));
1680  Vec3d result = stencil.getCenterCoord().asVec3d()
1681  - map.applyInverseMap(vectorFromSurface);
1682  return Vec3Type(result);
1683  } else {
1684  Vec3d location = map.applyMap(stencil.getCenterCoord().asVec3d());
1685  Vec3d result = map.applyInverseMap(location - vectorFromSurface);
1686  return Vec3Type(result);
1687  }
1688  }
1689 };
1690 
1691 
1695 template<typename MapType, DScheme DiffScheme>
1697 {
1698  // random access version
1699  template<typename Accessor> static Vec3<typename Accessor::ValueType>
1700  result(const MapType& map, const Accessor& grid, const Coord& ijk)
1701  {
1702  typedef typename Accessor::ValueType ValueType;
1703  typedef Vec3<ValueType> Vec3Type;
1704  // current distance
1705  ValueType d = grid.getValue(ijk);
1706  // compute gradient in physical space where it is a unit normal
1707  // since the grid holds a distance level set.
1708  Vec3Type vectorFromSurface =
1709  d*Gradient<MapType,DiffScheme>::result(map, grid, ijk);
1710  Vec3d result = map.applyMap(ijk.asVec3d()) - vectorFromSurface;
1711 
1712  return Vec3Type(result);
1713  }
1714 
1715  // stencil access version
1716  template<typename StencilT> static Vec3<typename StencilT::ValueType>
1717  result(const MapType& map, const StencilT& stencil)
1718  {
1719  typedef typename StencilT::ValueType ValueType;
1720  typedef Vec3<ValueType> Vec3Type;
1721  // current distance
1722  ValueType d = stencil.template getValue<0, 0, 0>();
1723  // compute gradient in physical space where it is a unit normal
1724  // since the grid holds a distance level set.
1725  Vec3Type vectorFromSurface =
1726  d*Gradient<MapType, DiffScheme>::result(map, stencil);
1727  Vec3d result = map.applyMap(stencil.getCenterCoord().asVec3d()) - vectorFromSurface;
1728 
1729  return Vec3Type(result);
1730  }
1731 };
1732 
1733 
1737 template<typename MapType, DDScheme DiffScheme2, DScheme DiffScheme1>
1739 {
1740  // random access version
1741  template<typename Accessor>
1742  static void compute(const MapType& map, const Accessor& grid, const Coord& ijk,
1743  double& alpha, double& beta)
1744  {
1745  typedef typename Accessor::ValueType ValueType;
1746  typedef Vec3<ValueType> Vec3Type;
1747 
1748  // compute the gradient in index space
1749  Vec3d d1_is(D1<DiffScheme1>::inX(grid, ijk),
1750  D1<DiffScheme1>::inY(grid, ijk),
1751  D1<DiffScheme1>::inZ(grid, ijk) );
1752 
1753  // all the second derivatives in index space
1754  ValueType iddx = D2<DiffScheme2>::inX(grid, ijk);
1755  ValueType iddy = D2<DiffScheme2>::inY(grid, ijk);
1756  ValueType iddz = D2<DiffScheme2>::inZ(grid, ijk);
1757 
1758  ValueType iddxy = D2<DiffScheme2>::inXandY(grid, ijk);
1759  ValueType iddyz = D2<DiffScheme2>::inYandZ(grid, ijk);
1760  ValueType iddxz = D2<DiffScheme2>::inXandZ(grid, ijk);
1761 
1762  // second derivatives in index space
1763  Mat3d d2_is(iddx, iddxy, iddxz,
1764  iddxy, iddy, iddyz,
1765  iddxz, iddyz, iddz);
1766 
1767  // convert to range space
1768  Mat3d d2_rs;
1769  Vec3d d1_rs;
1771  d2_rs = map.applyIJC(d2_is);
1772  d1_rs = map.applyIJT(d1_is);
1773  } else {
1774  d2_rs = map.applyIJC(d2_is, d1_is, ijk.asVec3d());
1775  d1_rs = map.applyIJT(d1_is, ijk.asVec3d());
1776  }
1777 
1778  // assemble the mean curvature
1779  double Dx2 = d1_rs(0)*d1_rs(0);
1780  double Dy2 = d1_rs(1)*d1_rs(1);
1781  double Dz2 = d1_rs(2)*d1_rs(2);
1782 
1783  // for return
1784  alpha = (Dx2*(d2_rs(1,1)+d2_rs(2,2))+Dy2*(d2_rs(0,0)+d2_rs(2,2))
1785  +Dz2*(d2_rs(0,0)+d2_rs(1,1))
1786  -2*(d1_rs(0)*(d1_rs(1)*d2_rs(0,1)+d1_rs(2)*d2_rs(0,2))
1787  +d1_rs(1)*d1_rs(2)*d2_rs(1,2)));
1788  beta = std::sqrt(Dx2 + Dy2 + Dz2); // * 1/dx
1789  }
1790 
1791  template<typename Accessor>
1792  static typename Accessor::ValueType result(const MapType& map,
1793  const Accessor& grid, const Coord& ijk)
1794  {
1795  typedef typename Accessor::ValueType ValueType;
1796  double alpha, beta;
1797  compute(map, grid, ijk, alpha, beta);
1798 
1799  return ValueType(alpha/(2. *math::Pow3(beta)));
1800  }
1801 
1802  template<typename Accessor>
1803  static typename Accessor::ValueType normGrad(const MapType& map,
1804  const Accessor& grid, const Coord& ijk)
1805  {
1806  typedef typename Accessor::ValueType ValueType;
1807  double alpha, beta;
1808  compute(map, grid, ijk, alpha, beta);
1809 
1810  return ValueType(alpha/(2. *math::Pow2(beta)));
1811  }
1812 
1813  // stencil access version
1814  template<typename StencilT>
1815  static void compute(const MapType& map, const StencilT& stencil, double& alpha, double& beta)
1816  {
1817  typedef typename StencilT::ValueType ValueType;
1818  typedef Vec3<ValueType> Vec3Type;
1819 
1820  // compute the gradient in index space
1821  Vec3d d1_is(D1<DiffScheme1>::inX(stencil),
1822  D1<DiffScheme1>::inY(stencil),
1823  D1<DiffScheme1>::inZ(stencil) );
1824 
1825  // all the second derivatives in index space
1826  ValueType iddx = D2<DiffScheme2>::inX(stencil);
1827  ValueType iddy = D2<DiffScheme2>::inY(stencil);
1828  ValueType iddz = D2<DiffScheme2>::inZ(stencil);
1829 
1830  ValueType iddxy = D2<DiffScheme2>::inXandY(stencil);
1831  ValueType iddyz = D2<DiffScheme2>::inYandZ(stencil);
1832  ValueType iddxz = D2<DiffScheme2>::inXandZ(stencil);
1833 
1834  // second derivatives in index space
1835  Mat3d d2_is(iddx, iddxy, iddxz,
1836  iddxy, iddy, iddyz,
1837  iddxz, iddyz, iddz);
1838 
1839  // convert to range space
1840  Mat3d d2_rs;
1841  Vec3d d1_rs;
1843  d2_rs = map.applyIJC(d2_is);
1844  d1_rs = map.applyIJT(d1_is);
1845  } else {
1846  d2_rs = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d());
1847  d1_rs = map.applyIJT(d1_is, stencil.getCenterCoord().asVec3d());
1848  }
1849 
1850  // assemble the mean curvature
1851  double Dx2 = d1_rs(0)*d1_rs(0);
1852  double Dy2 = d1_rs(1)*d1_rs(1);
1853  double Dz2 = d1_rs(2)*d1_rs(2);
1854 
1855  // for return
1856  alpha = (Dx2*(d2_rs(1,1)+d2_rs(2,2))+Dy2*(d2_rs(0,0)+d2_rs(2,2))
1857  +Dz2*(d2_rs(0,0)+d2_rs(1,1))
1858  -2*(d1_rs(0)*(d1_rs(1)*d2_rs(0,1)+d1_rs(2)*d2_rs(0,2))
1859  +d1_rs(1)*d1_rs(2)*d2_rs(1,2)));
1860  beta = std::sqrt(Dx2 + Dy2 + Dz2); // * 1/dx
1861  }
1862 
1863  template<typename StencilT>
1864  static typename StencilT::ValueType
1865  result(const MapType& map, const StencilT stencil)
1866  {
1867  typedef typename StencilT::ValueType ValueType;
1868  double alpha, beta;
1869  compute(map, stencil, alpha, beta);
1870  return ValueType(alpha/(2*math::Pow3(beta)));
1871  }
1872 
1873  template<typename StencilT>
1874  static typename StencilT::ValueType normGrad(const MapType& map, const StencilT stencil)
1875  {
1876  typedef typename StencilT::ValueType ValueType;
1877  double alpha, beta;
1878  compute(map, stencil, alpha, beta);
1879 
1880  return ValueType(alpha/(2*math::Pow2(beta)));
1881  }
1882 };
1883 
1884 
1885 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1886 struct MeanCurvature<TranslationMap, DiffScheme2, DiffScheme1>
1887 {
1888  // random access version
1889  template<typename Accessor>
1890  static typename Accessor::ValueType result(const TranslationMap&,
1891  const Accessor& grid, const Coord& ijk)
1892  {
1893  typedef typename Accessor::ValueType ValueType;
1894 
1895  ValueType alpha, beta;
1896  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
1897 
1898  return ValueType(alpha /(2*math::Pow3(beta)));
1899  }
1900 
1901  template<typename Accessor>
1902  static typename Accessor::ValueType normGrad(const TranslationMap&,
1903  const Accessor& grid, const Coord& ijk)
1904  {
1905  typedef typename Accessor::ValueType ValueType;
1906 
1907  ValueType alpha, beta;
1908  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
1909 
1910  return ValueType(alpha/(2*math::Pow2(beta)));
1911  }
1912 
1913  // stencil access version
1914  template<typename StencilT>
1915  static typename StencilT::ValueType result(const TranslationMap&, const StencilT& stencil)
1916  {
1917  typedef typename StencilT::ValueType ValueType;
1918 
1919  ValueType alpha, beta;
1921 
1922  return ValueType(alpha /(2*math::Pow3(beta)));
1923  }
1924 
1925  template<typename StencilT>
1926  static typename StencilT::ValueType normGrad(const TranslationMap&, const StencilT& stencil)
1927  {
1928  typedef typename StencilT::ValueType ValueType;
1929 
1930  ValueType alpha, beta;
1932 
1933  return ValueType(alpha/(2*math::Pow2(beta)));
1934  }
1935 };
1936 
1937 
1938 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1939 struct MeanCurvature<UniformScaleMap, DiffScheme2, DiffScheme1>
1940 {
1941  // random access version
1942  template<typename Accessor>
1943  static typename Accessor::ValueType result(const UniformScaleMap& map,
1944  const Accessor& grid, const Coord& ijk)
1945  {
1946  typedef typename Accessor::ValueType ValueType;
1947 
1948  ValueType alpha, beta;
1949  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
1950  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1951 
1952  return ValueType(alpha*inv2dx/math::Pow3(beta));
1953  }
1954 
1955  template<typename Accessor>
1956  static typename Accessor::ValueType normGrad(const UniformScaleMap& map,
1957  const Accessor& grid, const Coord& ijk)
1958  {
1959  typedef typename Accessor::ValueType ValueType;
1960 
1961  ValueType alpha, beta;
1962  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
1963  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1964 
1965  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
1966  }
1967 
1968  // stencil access version
1969  template<typename StencilT>
1970  static typename StencilT::ValueType result(const UniformScaleMap& map, const StencilT& stencil)
1971  {
1972  typedef typename StencilT::ValueType ValueType;
1973 
1974  ValueType alpha, beta;
1976  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
1977 
1978  return ValueType(alpha*inv2dx/math::Pow3(beta));
1979  }
1980 
1981  template<typename StencilT>
1982  static typename StencilT::ValueType normGrad(const UniformScaleMap& map, const StencilT& stencil)
1983  {
1984  typedef typename StencilT::ValueType ValueType;
1985 
1986  ValueType alpha, beta;
1988  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
1989 
1990  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
1991  }
1992 };
1993 
1994 
1995 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1996 struct MeanCurvature<UniformScaleTranslateMap, DiffScheme2, DiffScheme1>
1997 {
1998  // random access version
1999  template<typename Accessor> static typename Accessor::ValueType
2000  result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
2001  {
2002  typedef typename Accessor::ValueType ValueType;
2003 
2004  ValueType alpha, beta;
2005  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
2006  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
2007 
2008  return ValueType(alpha*inv2dx/math::Pow3(beta));
2009  }
2010 
2011  template<typename Accessor> static typename Accessor::ValueType
2012  normGrad(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk)
2013  {
2014  typedef typename Accessor::ValueType ValueType;
2015 
2016  ValueType alpha, beta;
2017  ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta);
2018  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
2019 
2020  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
2021  }
2022 
2023  // stencil access version
2024  template<typename StencilT> static typename StencilT::ValueType
2025  result(const UniformScaleTranslateMap& map, const StencilT& stencil)
2026  {
2027  typedef typename StencilT::ValueType ValueType;
2028 
2029  ValueType alpha, beta;
2031  ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]);
2032 
2033  return ValueType(alpha*inv2dx/math::Pow3(beta));
2034  }
2035 
2036  template<typename StencilT> static typename StencilT::ValueType
2037  normGrad(const UniformScaleTranslateMap& map, const StencilT& stencil)
2038  {
2039  typedef typename StencilT::ValueType ValueType;
2040 
2041  ValueType alpha, beta;
2043  ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]);
2044 
2045  return ValueType(alpha*invdxdx/(2*math::Pow2(beta)));
2046  }
2047 };
2048 
2049 
2055 {
2056 public:
2057  template<typename GridType>
2058  GenericMap(const GridType& g): mMap(g.transform().baseMap()) {}
2059 
2060  GenericMap(const Transform& t): mMap(t.baseMap()) {}
2061  GenericMap(MapBase::Ptr map): mMap(boost::const_pointer_cast<const MapBase>(map)) {}
2064 
2065  Vec3d applyMap(const Vec3d& in) const { return mMap->applyMap(in); }
2066  Vec3d applyInverseMap(const Vec3d& in) const { return mMap->applyInverseMap(in); }
2067 
2068  Vec3d applyIJT(const Vec3d& in) const { return mMap->applyIJT(in); }
2069  Vec3d applyIJT(const Vec3d& in, const Vec3d& pos) const { return mMap->applyIJT(in, pos); }
2070  Mat3d applyIJC(const Mat3d& m) const { return mMap->applyIJC(m); }
2071  Mat3d applyIJC(const Mat3d& m, const Vec3d& v, const Vec3d& pos) const
2072  { return mMap->applyIJC(m,v,pos); }
2073 
2074  double determinant() const { return mMap->determinant(); }
2075  double determinant(const Vec3d& in) const { return mMap->determinant(in); }
2076 
2077  Vec3d voxelSize() const { return mMap->voxelSize(); }
2078  Vec3d voxelSize(const Vec3d&v) const { return mMap->voxelSize(v); }
2079 
2080 private:
2082 };
2083 
2084 } // end math namespace
2085 } // namespace OPENVDB_VERSION_NAME
2086 } // end openvdb namespace
2087 
2088 #endif // OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED
2089 
2090 // Copyright (c) 2012-2013 DreamWorks Animation LLC
2091 // All rights reserved. This software is distributed under the
2092 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
A specialized Affine transform that scales along the principal axis the scaling is uniform in the thr...
Definition: Maps.h:926
Real GudonovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)
Definition: FiniteDifference.h:351
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1394
boost::shared_ptr< const MapBase > ConstPtr
Definition: Maps.h:162
static Accessor::ValueType result(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1890
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1530
Definition: FiniteDifference.h:196
static double result(const MapT &map, const AccessorType &grid, const Coord &ijk)
Definition: Operators.h:94
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:404
static StencilT::ValueType normGrad(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:1926
A specialized Affine transform that uniformaly scales along the principal axis and then translates th...
Definition: Maps.h:1476
math::Vec3< ValueType > Vec3Type
Definition: Operators.h:111
static Accessor::ValueType normGrad(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1956
static StencilT::ValueType normGrad(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1982
const Vec3d & getInvScaleSqr() const
Return the square of the scale. Used to optimize some finite difference calculations.
Definition: Maps.h:1341
static StencilT::ValueType::value_type result(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:975
const MapT & mMap
Definition: GridOperators.h:384
Definition: FiniteDifference.h:68
static StencilT::ValueType result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:1449
static StencilT::ValueType normGrad(const MapType &map, const StencilT stencil)
Definition: Operators.h:1874
static Accessor::ValueType::value_type result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:923
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1311
static Accessor::ValueType::value_type result(const ScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1197
static StencilT::ValueType::value_type result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1107
static Accessor::ValueType normGrad(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:2012
static Accessor::ValueType normGrad(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1902
static internal::ReturnValue< Accessor >::Vec3Type result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:664
T::ValueType ValueType
Definition: Operators.h:110
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:905
GenericMap(MapBase::ConstPtr map)
Definition: Operators.h:2062
Compute the mean curvature.
Definition: Operators.h:1738
static StencilT::ValueType::value_type result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1041
static Vec3< typename StencilT::ValueType > result(const StencilT &stencil, const Vec3Bias &V)
Definition: Operators.h:238
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:438
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1943
static internal::ReturnValue< StencilT >::Vec3Type result(const ScaleMap &map, const StencilT &stencil)
Definition: Operators.h:742
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:380
static internal::ReturnValue< Accessor >::Vec3Type result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:696
Vec3d applyIJT(const Vec3d &in, const Vec3d &pos) const
Definition: Operators.h:2069
Definition: FiniteDifference.h:62
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1413
double determinant() const
Definition: Operators.h:2074
static StencilT::ValueType::value_type result(const ScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1176
static Accessor::ValueType result(const ScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1576
~GenericMap()
Definition: Operators.h:2063
Vec3d applyIJT(const Vec3d &in) const
Definition: Operators.h:2068
const Vec3d & getInvTwiceScale() const
Return 1/(2 scale). Used to optimize some finite difference calculations.
Definition: Maps.h:1343
static Accessor::ValueType result(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1489
Type Pow2(Type x)
Return .
Definition: Math.h:467
static StencilT::ValueType::value_type result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1007
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:419
Vec3d applyMap(const Vec3d &in) const
Definition: Operators.h:2065
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:65
Definition: FiniteDifference.h:179
const Vec3d & getInvScaleSqr() const
Return the square of the scale. Used to optimize some finite difference calculations.
Definition: Maps.h:813
static internal::ReturnValue< StencilT >::Vec3Type result(const ScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:776
static Accessor::ValueType::value_type result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:993
Mat3d applyIJC(const Mat3d &m) const
Definition: Operators.h:2070
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:894
static Accessor::ValueType result(const ScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1611
Compute the Laplacian at a given location in a grid using finite differencing of various orders...
Definition: Operators.h:1409
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1369
Vec3d voxelSize() const
Definition: Operators.h:2077
DScheme
Different discrete schemes used in the first derivatives.
Definition: FiniteDifference.h:59
Definition: FiniteDifference.h:198
const MapType map
Definition: Operators.h:72
static StencilT::ValueType result(const ScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1625
Adapter for vector-valued world-space operators to return the vector magnitude.
Definition: Operators.h:92
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:523
static math::Vec3< typename StencilT::ValueType > result(const MapType &map, const StencilT &stencil, const Vec3< typename StencilT::ValueType > &V)
Definition: Operators.h:810
static StencilT::ValueType result(const UnitaryMap &, const StencilT &stencil)
Definition: Operators.h:1518
static Vec3< typename Accessor::ValueType > result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:126
Definition: FiniteDifference.h:70
static void result(const StencilT &stencil, typename StencilT::ValueType &alpha, typename StencilT::ValueType &beta)
Definition: Operators.h:573
static Accessor::ValueType::value_type result(const ScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1162
Compute the divergence of a vector-valued grid using differencing of various orders, the result defined with respect to the range-space of the map.
Definition: Operators.h:919
MapAdapter(const MapType &m)
Definition: Operators.h:62
static Vec3< typename Accessor::ValueType > result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1700
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1334
Definition: FiniteDifference.h:71
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:457
ResultType result(const AccessorType &grid, const Coord &ijk)
Definition: Operators.h:66
Vec3d asVec3d() const
Definition: Coord.h:136
Definition: FiniteDifference.h:1474
#define OPENVDB_VERSION_NAME
Definition: version.h:45
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:510
static void compute(const MapType &map, const Accessor &grid, const Coord &ijk, double &alpha, double &beta)
Definition: Operators.h:1742
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:47
static Accessor::ValueType::value_type result(const ScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1127
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1553
static Accessor::ValueType::value_type result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:484
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1970
Vec3d voxelSize(const Vec3d &v) const
Definition: Operators.h:2078
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:390
Biased Gradient Operators, using upwinding defined by the Vec3Bias input.
Definition: Operators.h:219
static internal::ReturnValue< Accessor >::Vec3Type result(const ScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:762
Definition: FiniteDifference.h:181
static Accessor::ValueType::value_type result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1027
Definition: FiniteDifference.h:67
static double result(const AccessorType &grid, const Coord &ijk)
Definition: Operators.h:80
static StencilT::ValueType::value_type result(const StencilT &stencil)
Definition: Operators.h:493
Definition: Operators.h:51
static internal::ReturnValue< Accessor >::Vec3Type result(const ScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:728
Compute the closest-point transform to a level set.
Definition: Operators.h:1696
static void compute(const MapType &map, const StencilT &stencil, double &alpha, double &beta)
Definition: Operators.h:1815
static math::Vec3< typename Accessor::ValueType > result(const MapType &map, const Accessor &grid, const Coord &ijk, const Vec3< typename Accessor::ValueType > &V)
Definition: Operators.h:798
Definition: FiniteDifference.h:439
Definition: FiniteDifference.h:65
static StencilT::ValueType result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:849
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:867
static StencilT::ValueType result(const ScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1591
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1792
Vec3d applyInverseMap(const Vec3d &in) const
Definition: Operators.h:2066
static internal::ReturnValue< StencilT >::Vec3Type result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:709
Compute the closest-point transform to a level set.
Definition: Operators.h:1643
Divergence operator defined in index space using various first derivative schemes.
Definition: Operators.h:480
ResultType result(const StencilType &stencil)
Definition: Operators.h:70
A specialized linear transform that performs a unitary maping i.e. rotation and or reflection...
Definition: Maps.h:1619
Definition: FiniteDifference.h:72
static math::Vec3< typename Accessor::ValueType > result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1647
static StencilT::ValueType::value_type result(const ScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1245
Definition: Operators.h:48
static internal::ReturnValue< StencilT >::Vec3Type result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:625
Compute the curl of a vector-valued grid using differencing of various orders in the space defined by...
Definition: Operators.h:1264
NineteenPointStencil< GridType > StencilType
Definition: Operators.h:200
Map traits.
Definition: Maps.h:79
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:878
double determinant(const Vec3d &in) const
Definition: Operators.h:2075
A specialized Affine transform that scales along the principal axis the scaling need not be uniform i...
Definition: Maps.h:1173
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:2025
static StencilT::ValueType result(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:1915
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1384
Definition: FiniteDifference.h:73
static internal::ReturnValue< StencilT >::Vec3Type result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:677
static internal::ReturnValue< StencilT >::Vec3Type result(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:650
Compute the mean curvature in index space.
Definition: Operators.h:540
static Accessor::ValueType::value_type result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1093
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1344
static StencilT::ValueType::value_type result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1074
static Accessor::ValueType result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1359
static StencilT::ValueType::value_type result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:940
static StencilT::ValueType result(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:1562
A wrapper that holds a MapBase::ConstPtr and exposes a reduced set of functionality needed by the mat...
Definition: Operators.h:2054
static Vec3< typename StencilT::ValueType > result(const StencilT &stencil)
Definition: Operators.h:137
static Accessor::ValueType::value_type result(const ScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1231
Curl operator defined in index space using various first derivative schemes.
Definition: Operators.h:506
NineteenPointStencil< GridType > StencilType
Definition: Operators.h:190
Definition: FiniteDifference.h:180
ThirteenPointStencil< GridType > StencilType
Definition: Operators.h:180
Definition: Operators.h:152
Definition: FiniteDifference.h:1720
static Vec3< typename Accessor::ValueType > result(const Accessor &grid, const Coord &ijk, const Vec3Bias &V)
Definition: Operators.h:226
Definition: FiniteDifference.h:74
static Accessor::ValueType result(const UnitaryMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1510
ValueType WENO5(const ValueType &v1, const ValueType &v2, const ValueType &v3, const ValueType &v4, const ValueType &v5, float scale2=0.01)
implimentation of nonimally fith-order finite-difference WENO. This function returns the numerical fl...
Definition: FiniteDifference.h:331
static Accessor::ValueType normGrad(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1803
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:836
GenericMap(const GridType &g)
Definition: Operators.h:2058
static math::Vec3< typename StencilT::ValueType > result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:1669
static StencilT::ValueType result(const TranslationMap &, const StencilT &stencil)
Definition: Operators.h:1497
A specialized Affine transform that scales along the principal axis the scaling need not be uniform i...
Definition: Maps.h:684
const Vec3d & getInvScale() const
Return 1/(scale)
Definition: Maps.h:1345
Biased gradient operators, defined with respect to the range-space of the map.
Definition: Operators.h:794
NineteenPointStencil< GridType > StencilType
Definition: Operators.h:210
Definition: FiniteDifference.h:66
Type Pow3(Type x)
Return .
Definition: Math.h:471
Gradient operators defined in index space of various orders.
Definition: Operators.h:122
Definition: FiniteDifference.h:195
const Vec3d & getInvTwiceScale() const
Return 1/(2 scale). Used to optimize some finite difference calculations.
Definition: Maps.h:815
static Vec3< typename StencilT::ValueType > result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:1717
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:260
static StencilT::ValueType normGrad(const UniformScaleTranslateMap &map, const StencilT &stencil)
Definition: Operators.h:2037
static double result(const MapT &map, const StencilType &stencil)
Definition: Operators.h:99
Adapter to associate a map with a world-space operator, giving it the same call signature as an index...
Definition: Operators.h:61
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition: Coord.h:102
Mat3d applyIJC(const Mat3d &m, const Vec3d &v, const Vec3d &pos) const
Definition: Operators.h:2071
GenericMap(MapBase::Ptr map)
Definition: Operators.h:2061
const Vec3d & getInvScale() const
Return 1/(scale)
Definition: Maps.h:817
static StencilT::ValueType result(const MapType &map, const StencilT &stencil)
Definition: Operators.h:1287
static double result(const StencilType &stencil)
Definition: Operators.h:85
static StencilT::ValueType::value_type result(const ScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1211
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
static internal::ReturnValue< Accessor >::Vec3Type result(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:642
static void result(const Accessor &grid, const Coord &ijk, typename Accessor::ValueType &alpha, typename Accessor::ValueType &beta)
Definition: Operators.h:544
static Accessor::ValueType::value_type result(const UniformScaleMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1061
Definition: Operators.h:827
static StencilT::ValueType result(const StencilT &stencil)
Definition: Operators.h:273
Abstract base class for maps.
Definition: Maps.h:158
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:1268
boost::shared_ptr< MapBase > Ptr
Definition: Maps.h:161
static internal::ReturnValue< Accessor >::Vec3Type result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:614
static StencilT::ValueType result(const MapType &map, const StencilT stencil)
Definition: Operators.h:1865
Laplacian defined in index space, using various center-difference stencils.
Definition: Operators.h:363
Adapter for vector-valued index-space operators to return the vector magnitude.
Definition: Operators.h:78
Center difference gradient operators, defined with respect to the range-space of the map...
Definition: Operators.h:609
static Accessor::ValueType result(const UniformScaleTranslateMap &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:2000
SevenPointStencil< GridType > StencilType
Definition: Operators.h:158
static StencilT::ValueType::value_type result(const ScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1141
Definition: FiniteDifference.h:197
Definition: FiniteDifference.h:69
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1320
Definition: FiniteDifference.h:194
SevenPointStencil< GridType > StencilType
Definition: Operators.h:169
static StencilT::ValueType result(const UniformScaleMap &map, const StencilT &stencil)
Definition: Operators.h:1539
GenericMap(const Transform &t)
Definition: Operators.h:2060
static Accessor::ValueType::value_type result(const TranslationMap &, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:963
A specialized linear transform that performs a translation.
Definition: Maps.h:998