mantaflow  0.10
A framework for fluid simulation
vectorbase.h
Go to the documentation of this file.
1 
2 /******************************************************************************
3  *
4  * MantaFlow fluid solver framework
5  * Copyright 2011-2016 Tobias Pfaff, Nils Thuerey
6  *
7  * This program is free software, distributed under the terms of the
8  * GNU General Public License (GPL)
9  * http://www.gnu.org/licenses
10  *
11  * Basic vector class
12  *
13  ******************************************************************************/
14 
15 #ifndef _VECTORBASE_H
16 #define _VECTORBASE_H
17 
18 // get rid of windos min/max defines
19 #if defined(WIN32) || defined(_WIN32)
20 # define NOMINMAX
21 #endif
22 
23 #include <stdio.h>
24 #include <string>
25 #include <limits>
26 #include <iostream>
27 #include "general.h"
28 
29 // if min/max are still around...
30 #if defined(WIN32) || defined(_WIN32)
31 # undef min
32 # undef max
33 #endif
34 
35 // redefine usage of some windows functions
36 #if defined(WIN32) || defined(_WIN32)
37 # ifndef snprintf
38 # define snprintf _snprintf
39 # endif
40 #endif
41 
42 // use which fp-precision? 1=float, 2=double
43 #ifndef FLOATINGPOINT_PRECISION
44 # define FLOATINGPOINT_PRECISION 1
45 #endif
46 
47 // VECTOR_EPSILON is the minimal vector length
48 // In order to be able to discriminate floating point values near zero, and
49 // to be sure not to fail a comparison because of roundoff errors, use this
50 // value as a threshold.
51 #if FLOATINGPOINT_PRECISION==1
52  typedef float Real;
53 # define VECTOR_EPSILON (1e-6f)
54 #else
55  typedef double Real;
56 # define VECTOR_EPSILON (1e-10)
57 #endif
58 
59 #ifndef M_PI
60 # define M_PI 3.1415926536
61 #endif
62 #ifndef M_E
63 # define M_E 2.7182818284
64 #endif
65 
66 namespace Manta
67 {
68 
70 template<class S>
71 class Vector3D
72 {
73 public:
75  inline Vector3D() : x(0),y(0),z(0) {}
76 
78  inline Vector3D ( const Vector3D<S> &v ) : x(v.x), y(v.y), z(v.z) {}
79 
81  inline Vector3D ( const float * v) : x((S)v[0]), y((S)v[1]), z((S)v[2]) {}
82 
84  inline Vector3D ( const double * v) : x((S)v[0]), y((S)v[1]), z((S)v[2]) {}
85 
87  inline Vector3D ( S v) : x(v), y(v), z(v) {}
88 
90  inline Vector3D ( S vx, S vy, S vz) : x(vx), y(vy), z(vz) {}
91 
92  // Operators
93 
95  inline const Vector3D<S>& operator= ( const Vector3D<S>& v ) {
96  x = v.x;
97  y = v.y;
98  z = v.z;
99  return *this;
100  }
102  inline const Vector3D<S>& operator= ( S s ) {
103  x = y = z = s;
104  return *this;
105  }
107  inline const Vector3D<S>& operator+= ( const Vector3D<S>& v ) {
108  x += v.x;
109  y += v.y;
110  z += v.z;
111  return *this;
112  }
114  inline const Vector3D<S>& operator+= ( S s ) {
115  x += s;
116  y += s;
117  z += s;
118  return *this;
119  }
121  inline const Vector3D<S>& operator-= ( const Vector3D<S>& v ) {
122  x -= v.x;
123  y -= v.y;
124  z -= v.z;
125  return *this;
126  }
128  inline const Vector3D<S>& operator-= ( S s ) {
129  x -= s;
130  y -= s;
131  z -= s;
132  return *this;
133  }
135  inline const Vector3D<S>& operator*= ( const Vector3D<S>& v ) {
136  x *= v.x;
137  y *= v.y;
138  z *= v.z;
139  return *this;
140  }
142  inline const Vector3D<S>& operator*= ( S s ) {
143  x *= s;
144  y *= s;
145  z *= s;
146  return *this;
147  }
149  inline const Vector3D<S>& operator/= ( const Vector3D<S>& v ) {
150  x /= v.x;
151  y /= v.y;
152  z /= v.z;
153  return *this;
154  }
156  inline const Vector3D<S>& operator/= ( S s ) {
157  x /= s;
158  y /= s;
159  z /= s;
160  return *this;
161  }
163  inline Vector3D<S> operator- () const {
164  return Vector3D<S> (-x, -y, -z);
165  }
166 
168  inline S min() const {
169  return ( x<y ) ? ( ( x<z ) ? x:z ) : ( ( y<z ) ? y:z );
170  }
172  inline S max() const {
173  return ( x>y ) ? ( ( x>z ) ? x:z ) : ( ( y>z ) ? y:z );
174  }
175 
177  inline bool empty() {
178  return x==0 && y==0 && z==0;
179  }
180 
182  inline S& operator[] ( unsigned int i ) {
183  return value[i];
184  }
186  inline const S& operator[] ( unsigned int i ) const {
187  return value[i];
188  }
189 
191  std::string toString() const;
192 
194  bool isValid() const;
195 
197  union {
198  S value[3];
199  struct {
200  S x;
201  S y;
202  S z;
203  };
204  struct {
205  S X;
206  S Y;
207  S Z;
208  };
209  };
210 
212  static const Vector3D<S> Zero, Invalid;
213 
215  inline Vector3D ( S vx, S vy, S vz, S vDummy) : x(vx), y(vy), z(vz) {}
216 
217 protected:
218 
219 };
220 
221 //************************************************************************
222 // Additional operators
223 //************************************************************************
224 
226 template<class S>
227 inline Vector3D<S> operator+ ( const Vector3D<S> &v1, const Vector3D<S> &v2 ) {
228  return Vector3D<S> ( v1.x+v2.x, v1.y+v2.y, v1.z+v2.z );
229 }
231 template<class S, class S2>
232 inline Vector3D<S> operator+ ( const Vector3D<S>& v, S2 s ) {
233  return Vector3D<S> ( v.x+s, v.y+s, v.z+s );
234 }
236 template<class S, class S2>
237 inline Vector3D<S> operator+ ( S2 s, const Vector3D<S>& v ) {
238  return Vector3D<S> ( v.x+s, v.y+s, v.z+s );
239 }
240 
242 template<class S>
243 inline Vector3D<S> operator- ( const Vector3D<S> &v1, const Vector3D<S> &v2 ) {
244  return Vector3D<S> ( v1.x-v2.x, v1.y-v2.y, v1.z-v2.z );
245 }
247 template<class S, class S2>
248 inline Vector3D<S> operator- ( const Vector3D<S>& v, S2 s ) {
249  return Vector3D<S> ( v.x-s, v.y-s, v.z-s );
250 }
252 template<class S, class S2>
253 inline Vector3D<S> operator- ( S2 s, const Vector3D<S>& v ) {
254  return Vector3D<S> ( s-v.x, s-v.y, s-v.z );
255 }
256 
258 template<class S>
259 inline Vector3D<S> operator* ( const Vector3D<S> &v1, const Vector3D<S> &v2 ) {
260  return Vector3D<S> ( v1.x*v2.x, v1.y*v2.y, v1.z*v2.z );
261 }
263 template<class S, class S2>
264 inline Vector3D<S> operator* ( const Vector3D<S>& v, S2 s ) {
265  return Vector3D<S> ( v.x*s, v.y*s, v.z*s );
266 }
268 template<class S, class S2>
269 inline Vector3D<S> operator* ( S2 s, const Vector3D<S>& v ) {
270  return Vector3D<S> ( s*v.x, s*v.y, s*v.z );
271 }
272 
274 template<class S>
275 inline Vector3D<S> operator/ ( const Vector3D<S> &v1, const Vector3D<S> &v2 ) {
276  return Vector3D<S> ( v1.x/v2.x, v1.y/v2.y, v1.z/v2.z );
277 }
279 template<class S, class S2>
280 inline Vector3D<S> operator/ ( const Vector3D<S>& v, S2 s ) {
281  return Vector3D<S> ( v.x/s, v.y/s, v.z/s );
282 }
284 template<class S, class S2>
285 inline Vector3D<S> operator/ ( S2 s, const Vector3D<S>& v ) {
286  return Vector3D<S> ( s/v.x, s/v.y, s/v.z );
287 }
288 
290 template<class S>
291 inline bool operator== (const Vector3D<S>& s1, const Vector3D<S>& s2) {
292  return s1.x == s2.x && s1.y == s2.y && s1.z == s2.z;
293 }
294 
296 template<class S>
297 inline bool operator!= (const Vector3D<S>& s1, const Vector3D<S>& s2) {
298  return s1.x != s2.x || s1.y != s2.y || s1.z != s2.z;
299 }
300 
301 //************************************************************************
302 // External functions
303 //************************************************************************
304 
306 template<class S>
307 inline Vector3D<S> vmin (const Vector3D<S>& s1, const Vector3D<S>& s2) {
308  return Vector3D<S>(std::min(s1.x,s2.x), std::min(s1.y,s2.y), std::min(s1.z,s2.z));
309 }
310 
312 template<class S, class S2>
313 inline Vector3D<S> vmin (const Vector3D<S>& s1, S2 s2) {
314  return Vector3D<S>(std::min(s1.x,s2), std::min(s1.y,s2), std::min(s1.z,s2));
315 }
316 
318 template<class S1, class S>
319 inline Vector3D<S> vmin (S1 s1, const Vector3D<S>& s2) {
320  return Vector3D<S>(std::min(s1,s2.x), std::min(s1,s2.y), std::min(s1,s2.z));
321 }
322 
324 template<class S>
325 inline Vector3D<S> vmax (const Vector3D<S>& s1, const Vector3D<S>& s2) {
326  return Vector3D<S>(std::max(s1.x,s2.x), std::max(s1.y,s2.y), std::max(s1.z,s2.z));
327 }
328 
330 template<class S, class S2>
331 inline Vector3D<S> vmax (const Vector3D<S>& s1, S2 s2) {
332  return Vector3D<S>(std::max(s1.x,s2), std::max(s1.y,s2), std::max(s1.z,s2));
333 }
334 
336 template<class S1, class S>
337 inline Vector3D<S> vmax (S1 s1, const Vector3D<S>& s2) {
338  return Vector3D<S>(std::max(s1,s2.x), std::max(s1,s2.y), std::max(s1,s2.z));
339 }
340 
342 template<class S>
343 inline S dot ( const Vector3D<S> &t, const Vector3D<S> &v ) {
344  return t.x*v.x + t.y*v.y + t.z*v.z;
345 }
346 
348 template<class S>
349 inline Vector3D<S> cross ( const Vector3D<S> &t, const Vector3D<S> &v ) {
350  Vector3D<S> cp (
351  ( ( t.y*v.z ) - ( t.z*v.y ) ),
352  ( ( t.z*v.x ) - ( t.x*v.z ) ),
353  ( ( t.x*v.y ) - ( t.y*v.x ) ) );
354  return cp;
355 }
356 
358 
363 template<class S>
364 inline const Vector3D<S>& projectNormalTo ( const Vector3D<S>& v, const Vector3D<S> &n) {
365  S sprod = dot (v, n);
366  return v - n * dot(v, n);
367 }
368 
371 template<class S>
372 inline S norm ( const Vector3D<S>& v ) {
373  S l = v.x*v.x + v.y*v.y + v.z*v.z;
374  if ( l <= VECTOR_EPSILON*VECTOR_EPSILON ) return(0.);
375  return ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON ) ? 1. : sqrt ( l );
376 }
377 
379 template<class S>
380 inline S normSquare ( const Vector3D<S>& v ) {
381  return v.x*v.x + v.y*v.y + v.z*v.z;
382 }
383 
385 inline Real norm(const Real v) { return fabs(v); }
386 inline Real normSquare(const Real v) { return square(v); }
387 inline Real norm(const int v) { return abs(v); }
388 inline Real normSquare(const int v) { return square(v); }
389 
391 template<class S>
392 inline Vector3D<S> getNormalized ( const Vector3D<S>& v ) {
393  S l = v.x*v.x + v.y*v.y + v.z*v.z;
394  if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON )
395  return v; /* normalized "enough"... */
396  else if ( l > VECTOR_EPSILON*VECTOR_EPSILON )
397  {
398  S fac = 1./sqrt ( l );
399  return Vector3D<S> ( v.x*fac, v.y*fac, v.z*fac );
400  }
401  else
402  return Vector3D<S> ( ( S ) 0 );
403 }
404 
406 
407 template<class S>
408 inline S normalize ( Vector3D<S> &v ) {
409  S norm;
410  S l = v.x*v.x + v.y*v.y + v.z*v.z;
411  if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON ) {
412  norm = 1.;
413  } else if ( l > VECTOR_EPSILON*VECTOR_EPSILON ) {
414  norm = sqrt ( l );
415  v *= 1./norm;
416  } else {
417  v = Vector3D<S>::Zero;
418  norm = 0.;
419  }
420  return ( S ) norm;
421 }
422 
424 
427 template<class S>
428 Vector3D<S> getOrthogonalVector(const Vector3D<S>& v) {
429  // Determine the component with max. absolute value
430  int maxIndex= ( fabs ( v.x ) > fabs ( v.y ) ) ? 0 : 1;
431  maxIndex= ( fabs ( v[maxIndex] ) > fabs ( v.z ) ) ? maxIndex : 2;
432 
433  // Choose another axis than the one with max. component and project
434  // orthogonal to self
435  Vector3D<S> o ( 0.0 );
436  o[ ( maxIndex+1 ) %3]= 1;
437 
438  Vector3D<S> c = cross(v, o);
439  normalize(c);
440  return c;
441 }
442 
444 
449 template<class S>
450 inline void vecToAngle ( const Vector3D<S>& v, S& phi, S& theta )
451 {
452  if ( fabs ( v.y ) < VECTOR_EPSILON )
453  theta = M_PI/2;
454  else if ( fabs ( v.x ) < VECTOR_EPSILON && fabs ( v.z ) < VECTOR_EPSILON )
455  theta = ( v.y>=0 ) ? 0:M_PI;
456  else
457  theta = atan ( sqrt ( v.x*v.x+v.z*v.z ) /v.y );
458  if ( theta<0 ) theta+=M_PI;
459 
460  if ( fabs ( v.x ) < VECTOR_EPSILON )
461  phi = M_PI/2;
462  else
463  phi = atan ( v.z/v.x );
464  if ( phi<0 ) phi+=M_PI;
465  if ( fabs ( v.z ) < VECTOR_EPSILON )
466  phi = ( v.x>=0 ) ? 0 : M_PI;
467  else if ( v.z < 0 )
468  phi += M_PI;
469 }
470 
472 
479 template<class S>
480 inline Vector3D<S> reflectVector ( const Vector3D<S>& t, const Vector3D<S>& n ) {
481  Vector3D<S> nn= ( dot ( t, n ) > 0.0 ) ? ( n*-1.0 ) : n;
482  return ( t - nn * ( 2.0 * dot ( nn, t ) ) );
483 }
484 
486 
493 template<class S>
494 inline Vector3D<S> refractVector ( const Vector3D<S> &t, const Vector3D<S> &normal, S nt, S nair, int &refRefl ) {
495  // from Glassner's book, section 5.2 (Heckberts method)
496  S eta = nair / nt;
497  S n = -dot ( t, normal );
498  S tt = 1.0 + eta*eta* ( n*n-1.0 );
499  if ( tt<0.0 ) {
500  // we have total reflection!
501  refRefl = 1;
502  } else {
503  // normal reflection
504  tt = eta*n - sqrt ( tt );
505  return ( t*eta + normal*tt );
506  }
507  return t;
508 }
509 
511 template<class S> std::string Vector3D<S>::toString() const {
512  char buf[256];
513  snprintf ( buf,256,"[%+4.6f,%+4.6f,%+4.6f]", ( double ) ( *this ) [0], ( double ) ( *this ) [1], ( double ) ( *this ) [2] );
514  // for debugging, optionally increase precision:
515  //snprintf ( buf,256,"[%+4.16f,%+4.16f,%+4.16f]", ( double ) ( *this ) [0], ( double ) ( *this ) [1], ( double ) ( *this ) [2] );
516  return std::string ( buf );
517 }
518 
519 template<> std::string Vector3D<int>::toString() const;
520 
521 
523 
524 template<class S>
525 std::ostream& operator<< ( std::ostream& os, const Vector3D<S>& i ) {
526  os << i.toString();
527  return os;
528 }
529 
531 
532 template<class S>
533 std::istream& operator>> ( std::istream& is, Vector3D<S>& i ) {
534  char c;
535  char dummy[3];
536  is >> c >> i[0] >> dummy >> i[1] >> dummy >> i[2] >> c;
537  return is;
538 }
539 
540 /**************************************************************************/
541 // Define default vector alias
542 /**************************************************************************/
543 
546 
549 
551 template<class T> inline Vec3 toVec3 ( T v ) {
552  return Vec3 ( v[0],v[1],v[2] );
553 }
554 
556 template<class T> inline Vec3i toVec3i ( T v ) {
557  return Vec3i ( ( int ) v[0], ( int ) v[1], ( int ) v[2] );
558 }
559 
561 template<class T> inline Vec3i toVec3i ( T v0, T v1, T v2 ) {
562  return Vec3i ( ( int ) v0, ( int ) v1, ( int ) v2 );
563 }
564 
566 template<class T> inline Vec3i toVec3iRound ( T v ) {
567  return Vec3i ( ( int ) round ( v[0] ), ( int ) round ( v[1] ), ( int ) round ( v[2] ) );
568 }
569 
571 template<class T> inline Vec3i toVec3iChecked ( T v ) {
572  Vec3i ret;
573  for (size_t i=0; i<3; i++) {
574  Real a = v[i];
575  if (fabs(a-floor(a+0.5)) > 1e-5)
576  errMsg("argument is not an int, cannot convert");
577  ret[i] = (int) (a+0.5);
578  }
579  return ret;
580 }
581 
583 template<class T> inline Vector3D<double> toVec3d ( T v ) {
584  return Vector3D<double> ( v[0], v[1], v[2] );
585 }
586 
588 template<class T> inline Vector3D<float> toVec3f ( T v ) {
589  return Vector3D<float> ( v[0], v[1], v[2] );
590 }
591 
592 
593 /**************************************************************************/
594 // Specializations for common math functions
595 /**************************************************************************/
596 
597 template<> inline Vec3 clamp<Vec3>(const Vec3& a, const Vec3& b, const Vec3& c) {
598  return Vec3 ( clamp(a.x, b.x, c.x),
599  clamp(a.y, b.y, c.y),
600  clamp(a.z, b.z, c.z) );
601 }
602 template<> inline Vec3 safeDivide<Vec3>(const Vec3 &a, const Vec3& b) {
603  return Vec3(safeDivide(a.x,b.x), safeDivide(a.y,b.y), safeDivide(a.z,b.z));
604 }
605 template<> inline Vec3 nmod<Vec3>(const Vec3& a, const Vec3& b) {
606  return Vec3(nmod(a.x,b.x),nmod(a.y,b.y),nmod(a.z,b.z));
607 }
608 
609 }; // namespace
610 
611 
612 #endif
613 
614 
Vector3D()
Constructor.
Definition: vectorbase.h:75
Definition: commonkernels.h:22
Vector3D(const Vector3D< S > &v)
Copy-Constructor.
Definition: vectorbase.h:78
const Vector3D< S > & operator=(const Vector3D< S > &v)
Assignment operator.
Definition: vectorbase.h:95
std::string toString() const
debug output vector to a string
S & operator[](unsigned int i)
access operator
Definition: vectorbase.h:182
Basic inlined vector class.
Definition: vectorbase.h:71
static const Vector3D< S > Zero
zero element
Definition: vectorbase.h:212
Vector3D(S vx, S vy, S vz)
Construct a vector from three Ss.
Definition: vectorbase.h:90
Vector3D< S > operator-() const
Negation operator.
Definition: vectorbase.h:163
Vector3D(S vx, S vy, S vz, S vDummy)
For compatibility with 4d vectors (discards 4th comp)
Definition: vectorbase.h:215
Vector3D(const double *v)
Copy-Constructor.
Definition: vectorbase.h:84
const Vector3D< S > & operator-=(const Vector3D< S > &v)
Assign and sub operator.
Definition: vectorbase.h:121
const Vector3D< S > & operator*=(const Vector3D< S > &v)
Assign and mult operator.
Definition: vectorbase.h:135
S max() const
Get biggest component.
Definition: vectorbase.h:172
bool isValid() const
test if nans are present
Vector3D(const float *v)
Copy-Constructor.
Definition: vectorbase.h:81
const Vector3D< S > & operator/=(const Vector3D< S > &v)
Assign and div operator.
Definition: vectorbase.h:149
S min() const
Get smallest component.
Definition: vectorbase.h:168
Vector3D(S v)
Construct a vector from one S.
Definition: vectorbase.h:87
bool empty()
Test if all components are zero.
Definition: vectorbase.h:177
const Vector3D< S > & operator+=(const Vector3D< S > &v)
Assign and add operator.
Definition: vectorbase.h:107