Page MenuHome
Differential D3850 Diff 12284 intern/mantaflow/intern/manta_pp/omp/util/vector4d.h

Changeset View

Standalone View

View Options

intern/mantaflow/intern/manta_pp/omp/util/vector4d.h

  • This file was added.
/******************************************************************************
*
* MantaFlow fluid solver framework
* Copyright 2011 Tobias Pfaff, Nils Thuerey
*
* This program is free software, distributed under the terms of the
* GNU General Public License (GPL)
* http://www.gnu.org/licenses
*
* 4D vector class
*
******************************************************************************/
#ifndef _VECTOR4D_H
#define _VECTOR4D_H
#include "vectorbase.h"
namespace Manta
{
//! Basic inlined vector class
template<class S>
class Vector4D
{
public:
//! Constructor
inline Vector4D() : x(0),y(0),z(0),t(0) {}
//! Copy-Constructor
inline Vector4D ( const Vector4D<S> &v ) : x(v.x), y(v.y), z(v.z),t(v.t) {}
//! Copy-Constructor
inline Vector4D ( const float * v) : x((S)v[0]), y((S)v[1]), z((S)v[2]), t((S)v[3]) {}
//! Copy-Constructor
inline Vector4D ( const double * v) : x((S)v[0]), y((S)v[1]), z((S)v[2]), t((S)v[3]) {}
//! Construct a vector from one S
inline Vector4D ( S v) : x(v), y(v), z(v), t(v) {}
//! Construct a vector from three Ss
inline Vector4D ( S vx, S vy, S vz, S vw) : x(vx), y(vy), z(vz), t(vw) {}
// Operators
//! Assignment operator
inline const Vector4D<S>& operator= ( const Vector4D<S>& v ) {
x = v.x;
y = v.y;
z = v.z;
t = v.t;
return *this;
}
//! Assignment operator
inline const Vector4D<S>& operator= ( S s ) {
x = y = z = t = s;
return *this;
}
//! Assign and add operator
inline const Vector4D<S>& operator+= ( const Vector4D<S>& v ) {
x += v.x;
y += v.y;
z += v.z;
t += v.t;
return *this;
}
//! Assign and add operator
inline const Vector4D<S>& operator+= ( S s ) {
x += s;
y += s;
z += s;
t += s;
return *this;
}
//! Assign and sub operator
inline const Vector4D<S>& operator-= ( const Vector4D<S>& v ) {
x -= v.x;
y -= v.y;
z -= v.z;
t -= v.t;
return *this;
}
//! Assign and sub operator
inline const Vector4D<S>& operator-= ( S s ) {
x -= s;
y -= s;
z -= s;
t -= s;
return *this;
}
//! Assign and mult operator
inline const Vector4D<S>& operator*= ( const Vector4D<S>& v ) {
x *= v.x;
y *= v.y;
z *= v.z;
t *= v.t;
return *this;
}
//! Assign and mult operator
inline const Vector4D<S>& operator*= ( S s ) {
x *= s;
y *= s;
z *= s;
t *= s;
return *this;
}
//! Assign and div operator
inline const Vector4D<S>& operator/= ( const Vector4D<S>& v ) {
x /= v.x;
y /= v.y;
z /= v.z;
t /= v.t;
return *this;
}
//! Assign and div operator
inline const Vector4D<S>& operator/= ( S s ) {
x /= s;
y /= s;
z /= s;
t /= s;
return *this;
}
//! Negation operator
inline Vector4D<S> operator- () const {
return Vector4D<S> (-x, -y, -z, -t);
}
//! Get smallest component
//inline S min() const { return ( x<y ) ? ( ( x<z ) ? x:z ) : ( ( y<z ) ? y:z ); }
//! Get biggest component
//inline S max() const { return ( x>y ) ? ( ( x>z ) ? x:z ) : ( ( y>z ) ? y:z ); }
//! Test if all components are zero
inline bool empty() {
return x==0 && y==0 && z==0 && t==0;
}
//! access operator
inline S& operator[] ( unsigned int i ) {
return value[i];
}
//! constant access operator
inline const S& operator[] ( unsigned int i ) const {
return value[i];
}
//! debug output vector to a string
std::string toString() const;
//! test if nans are present
bool isValid() const;
//! actual values
union {
S value[4];
struct {
S x;
S y;
S z;
S t;
};
struct {
S X;
S Y;
S Z;
S T;
};
};
// zero element
static const Vector4D<S> Zero, Invalid;
protected:
};
//************************************************************************
// Additional operators
//************************************************************************
//! Addition operator
template<class S>
inline Vector4D<S> operator+ ( const Vector4D<S> &v1, const Vector4D<S> &v2 ) {
return Vector4D<S> ( v1.x+v2.x, v1.y+v2.y, v1.z+v2.z, v1.t+v2.t );
}
//! Addition operator
template<class S, class S2>
inline Vector4D<S> operator+ ( const Vector4D<S>& v, S2 s ) {
return Vector4D<S> ( v.x+s, v.y+s, v.z+s, v.t+s );
}
//! Addition operator
template<class S, class S2>
inline Vector4D<S> operator+ ( S2 s, const Vector4D<S>& v ) {
return Vector4D<S> ( v.x+s, v.y+s, v.z+s, v.t+s );
}
//! Subtraction operator
template<class S>
inline Vector4D<S> operator- ( const Vector4D<S> &v1, const Vector4D<S> &v2 ) {
return Vector4D<S> ( v1.x-v2.x, v1.y-v2.y, v1.z-v2.z, v1.t-v2.t );
}
//! Subtraction operator
template<class S, class S2>
inline Vector4D<S> operator- ( const Vector4D<S>& v, S2 s ) {
return Vector4D<S> ( v.x-s, v.y-s, v.z-s, v.t-s );
}
//! Subtraction operator
template<class S, class S2>
inline Vector4D<S> operator- ( S2 s, const Vector4D<S>& v ) {
return Vector4D<S> ( s-v.x, s-v.y, s-v.z, s-v.t );
}
//! Multiplication operator
template<class S>
inline Vector4D<S> operator* ( const Vector4D<S> &v1, const Vector4D<S> &v2 ) {
return Vector4D<S> ( v1.x*v2.x, v1.y*v2.y, v1.z*v2.z, v1.t*v2.t );
}
//! Multiplication operator
template<class S, class S2>
inline Vector4D<S> operator* ( const Vector4D<S>& v, S2 s ) {
return Vector4D<S> ( v.x*s, v.y*s, v.z*s , v.t*s );
}
//! Multiplication operator
template<class S, class S2>
inline Vector4D<S> operator* ( S2 s, const Vector4D<S>& v ) {
return Vector4D<S> ( s*v.x, s*v.y, s*v.z, s*v.t );
}
//! Division operator
template<class S>
inline Vector4D<S> operator/ ( const Vector4D<S> &v1, const Vector4D<S> &v2 ) {
return Vector4D<S> ( v1.x/v2.x, v1.y/v2.y, v1.z/v2.z, v1.t/v2.t );
}
//! Division operator
template<class S, class S2>
inline Vector4D<S> operator/ ( const Vector4D<S>& v, S2 s ) {
return Vector4D<S> ( v.x/s, v.y/s, v.z/s, v.t/s );
}
//! Division operator
template<class S, class S2>
inline Vector4D<S> operator/ ( S2 s, const Vector4D<S>& v ) {
return Vector4D<S> ( s/v.x, s/v.y, s/v.z, s/v.t );
}
//! Comparison operator
template<class S>
inline bool operator== (const Vector4D<S>& s1, const Vector4D<S>& s2) {
return s1.x == s2.x && s1.y == s2.y && s1.z == s2.z && s1.t == s2.t;
}
//! Comparison operator
template<class S>
inline bool operator!= (const Vector4D<S>& s1, const Vector4D<S>& s2) {
return s1.x != s2.x || s1.y != s2.y || s1.z != s2.z || s1.t != s2.t;
}
//************************************************************************
// External functions
//************************************************************************
//! Dot product
template<class S>
inline S dot ( const Vector4D<S> &t, const Vector4D<S> &v ) {
return t.x*v.x + t.y*v.y + t.z*v.z + t.t*v.t;
}
//! Cross product
/*template<class S>
inline Vector4D<S> cross ( const Vector4D<S> &t, const Vector4D<S> &v ) {
NYI Vector4D<S> cp (
( ( t.y*v.z ) - ( t.z*v.y ) ),
( ( t.z*v.x ) - ( t.x*v.z ) ),
( ( t.x*v.y ) - ( t.y*v.x ) ) );
return cp;
}*/
//! Compute the magnitude (length) of the vector
template<class S>
inline S norm ( const Vector4D<S>& v ) {
S l = v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
return ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON ) ? 1. : sqrt ( l );
}
//! Compute squared magnitude
template<class S>
inline S normSquare ( const Vector4D<S>& v ) {
return v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
}
//! Returns a normalized vector
template<class S>
inline Vector4D<S> getNormalized ( const Vector4D<S>& v ) {
S l = v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON )
return v; /* normalized "enough"... */
else if ( l > VECTOR_EPSILON*VECTOR_EPSILON )
{
S fac = 1./sqrt ( l );
return Vector4D<S> ( v.x*fac, v.y*fac, v.z*fac , v.t*fac );
}
else
return Vector4D<S> ( ( S ) 0 );
}
//! Compute the norm of the vector and normalize it.
/*! \return The value of the norm */
template<class S>
inline S normalize ( Vector4D<S> &v ) {
S norm;
S l = v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON ) {
norm = 1.;
} else if ( l > VECTOR_EPSILON*VECTOR_EPSILON ) {
norm = sqrt ( l );
v *= 1./norm;
} else {
v = Vector4D<S>::Zero;
norm = 0.;
}
return ( S ) norm;
}
//! Outputs the object in human readable form as string
template<class S> std::string Vector4D<S>::toString() const {
char buf[256];
snprintf ( buf,256,"[%+4.6f,%+4.6f,%+4.6f,%+4.6f]", ( double ) ( *this ) [0], ( double ) ( *this ) [1], ( double ) ( *this ) [2] , ( double ) ( *this ) [3] );
// for debugging, optionally increase precision:
//snprintf ( buf,256,"[%+4.16f,%+4.16f,%+4.16f,%+4.16f]", ( double ) ( *this ) [0], ( double ) ( *this ) [1], ( double ) ( *this ) [2], ( double ) ( *this ) [3] );
return std::string ( buf );
}
template<> std::string Vector4D<int>::toString() const;
//! Outputs the object in human readable form to stream
template<class S>
std::ostream& operator<< ( std::ostream& os, const Vector4D<S>& i ) {
os << i.toString();
return os;
}
//! Reads the contents of the object from a stream
template<class S>
std::istream& operator>> ( std::istream& is, Vector4D<S>& i ) {
char c;
char dummy[4];
is >> c >> i[0] >> dummy >> i[1] >> dummy >> i[2] >> dummy >> i[3] >> c;
return is;
}
/**************************************************************************/
// Define default vector alias
/**************************************************************************/
//! 3D vector class of type Real (typically float)
typedef Vector4D<Real> Vec4;
//! 3D vector class of type int
typedef Vector4D<int> Vec4i;
//! convert to Real Vector
template<class T> inline Vec4 toVec4 ( T v ) {
return Vec4 ( v[0],v[1],v[2],v[3] );
}
template<class T> inline Vec4i toVec4i ( T v ) {
return Vec4i( v[0],v[1],v[2],v[3] );
}
/**************************************************************************/
// Specializations for common math functions
/**************************************************************************/
template<> inline Vec4 clamp<Vec4>(const Vec4& a, const Vec4& b, const Vec4& c) {
return Vec4 ( clamp(a.x, b.x, c.x),
clamp(a.y, b.y, c.y),
clamp(a.z, b.z, c.z),
clamp(a.t, b.t, c.t) );
}
template<> inline Vec4 safeDivide<Vec4>(const Vec4 &a, const Vec4& b) {
return Vec4(safeDivide(a.x,b.x), safeDivide(a.y,b.y), safeDivide(a.z,b.z), safeDivide(a.t,b.t));
}
template<> inline Vec4 nmod<Vec4>(const Vec4& a, const Vec4& b) {
return Vec4(nmod(a.x,b.x),nmod(a.y,b.y),nmod(a.z,b.z) ,nmod(a.t,b.t) );
}
/**************************************************************************/
// 4d interpolation (note only 4d here, 2d/3d interpolations are in interpol.h)
/**************************************************************************/
#define BUILD_INDEX_4D \
Real px=pos.x-0.5f, py=pos.y-0.5f, pz=pos.z-0.5f, pt=pos.t-0.5f; \
int xi = (int)px; \
int yi = (int)py; \
int zi = (int)pz; \
int ti = (int)pt; \
Real s1 = px-(Real)xi, s0 = 1.-s1; \
Real t1 = py-(Real)yi, t0 = 1.-t1; \
Real f1 = pz-(Real)zi, f0 = 1.-f1; \
Real g1 = pt-(Real)ti, g0 = 1.-g1; \
/* clamp to border */ \
if (px < 0.) { xi = 0; s0 = 1.0; s1 = 0.0; } \
if (py < 0.) { yi = 0; t0 = 1.0; t1 = 0.0; } \
if (pz < 0.) { zi = 0; f0 = 1.0; f1 = 0.0; } \
if (pt < 0.) { ti = 0; g0 = 1.0; g1 = 0.0; } \
if (xi >= size.x-1) { xi = size.x-2; s0 = 0.0; s1 = 1.0; } \
if (yi >= size.y-1) { yi = size.y-2; t0 = 0.0; t1 = 1.0; } \
if (zi >= size.z-1) { zi = size.z-2; f0 = 0.0; f1 = 1.0; } \
if (ti >= size.t-1) { ti = size.t-2; g0 = 0.0; g1 = 1.0; } \
const int sX = 1; \
const int sY = size.x;
static inline void checkIndexInterpol4d(const Vec4i& size, int idx) {
if (idx<0 || idx > size.x * size.y * size.z * size.t) {
std::ostringstream s;
s << "Grid interpol4d dim " << size << " : index " << idx << " out of bound ";
errMsg(s.str());
}
}
template <class T>
inline T interpol4d(const T* data, const Vec4i& size, const IndexInt sZ, const IndexInt sT, const Vec4& pos) {
BUILD_INDEX_4D
IndexInt idx = (IndexInt)xi + sY * (IndexInt)yi + sZ * (IndexInt)zi + sT * (IndexInt)ti;
DEBUG_ONLY(checkIndexInterpol4d(size,idx));
DEBUG_ONLY(checkIndexInterpol4d(size,idx+sX+sY+sZ+sT));
return( ((data[idx] *t0 + data[idx+sY] *t1) * s0
+ (data[idx+sX] *t0 + data[idx+sX+sY] *t1) * s1) * f0
+((data[idx+sZ] *t0 + data[idx+sY+sZ] *t1) * s0
+ (data[idx+sX+sZ]*t0 + data[idx+sX+sY+sZ]*t1) * s1) * f1 ) * g0
+
( ((data[idx+sT] *t0 + data[idx+sT+sY] *t1) * s0
+ (data[idx+sT+sX] *t0 + data[idx+sT+sX+sY] *t1) * s1) * f0
+((data[idx+sT+sZ] *t0 + data[idx+sT+sY+sZ] *t1) * s0
+ (data[idx+sT+sX+sZ]*t0 + data[idx+sT+sX+sY+sZ]*t1) * s1) * f1 ) * g1 ;
}
}; // namespace
#endif