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