diff options
Diffstat (limited to 'extern/mantaflow/helper/util/vectorbase.h')
| -rw-r--r-- | extern/mantaflow/helper/util/vectorbase.h | 679 |
1 files changed, 679 insertions, 0 deletions
diff --git a/extern/mantaflow/helper/util/vectorbase.h b/extern/mantaflow/helper/util/vectorbase.h new file mode 100644 index 00000000000..a3135431eb3 --- /dev/null +++ b/extern/mantaflow/helper/util/vectorbase.h @@ -0,0 +1,679 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011-2016 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Basic vector class + * + ******************************************************************************/ + +#ifndef _VECTORBASE_H +#define _VECTORBASE_H + +// get rid of windos min/max defines +#if defined(WIN32) || defined(_WIN32) +# define NOMINMAX +#endif + +#include <stdio.h> +#include <string> +#include <limits> +#include <iostream> +#include "general.h" + +// if min/max are still around... +#if defined(WIN32) || defined(_WIN32) +# undef min +# undef max +#endif + +// redefine usage of some windows functions +#if defined(WIN32) || defined(_WIN32) +# ifndef snprintf +# define snprintf _snprintf +# endif +#endif + +// use which fp-precision? 1=float, 2=double +#ifndef FLOATINGPOINT_PRECISION +# define FLOATINGPOINT_PRECISION 1 +#endif + +// VECTOR_EPSILON is the minimal vector length +// In order to be able to discriminate floating point values near zero, and +// to be sure not to fail a comparison because of roundoff errors, use this +// value as a threshold. +#if FLOATINGPOINT_PRECISION == 1 +typedef float Real; +# define VECTOR_EPSILON (1e-6f) +#else +typedef double Real; +# define VECTOR_EPSILON (1e-10) +#endif + +#ifndef M_PI +# define M_PI 3.1415926536 +#endif +#ifndef M_E +# define M_E 2.7182818284 +#endif + +namespace Manta { + +//! Basic inlined vector class +template<class S> class Vector3D { + public: + //! Constructor + inline Vector3D() : x(0), y(0), z(0) + { + } + + //! Copy-Constructor + inline Vector3D(const Vector3D<S> &v) : x(v.x), y(v.y), z(v.z) + { + } + + //! Copy-Constructor + inline Vector3D(const int *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]) + { + } + + //! Copy-Constructor + inline Vector3D(const float *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]) + { + } + + //! Copy-Constructor + inline Vector3D(const double *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]) + { + } + + //! Construct a vector from one S + inline Vector3D(S v) : x(v), y(v), z(v) + { + } + + //! Construct a vector from three Ss + inline Vector3D(S vx, S vy, S vz) : x(vx), y(vy), z(vz) + { + } + + // Operators + + //! Assignment operator + inline const Vector3D<S> &operator=(const Vector3D<S> &v) + { + x = v.x; + y = v.y; + z = v.z; + return *this; + } + //! Assignment operator + inline const Vector3D<S> &operator=(S s) + { + x = y = z = s; + return *this; + } + //! Assign and add operator + inline const Vector3D<S> &operator+=(const Vector3D<S> &v) + { + x += v.x; + y += v.y; + z += v.z; + return *this; + } + //! Assign and add operator + inline const Vector3D<S> &operator+=(S s) + { + x += s; + y += s; + z += s; + return *this; + } + //! Assign and sub operator + inline const Vector3D<S> &operator-=(const Vector3D<S> &v) + { + x -= v.x; + y -= v.y; + z -= v.z; + return *this; + } + //! Assign and sub operator + inline const Vector3D<S> &operator-=(S s) + { + x -= s; + y -= s; + z -= s; + return *this; + } + //! Assign and mult operator + inline const Vector3D<S> &operator*=(const Vector3D<S> &v) + { + x *= v.x; + y *= v.y; + z *= v.z; + return *this; + } + //! Assign and mult operator + inline const Vector3D<S> &operator*=(S s) + { + x *= s; + y *= s; + z *= s; + return *this; + } + //! Assign and div operator + inline const Vector3D<S> &operator/=(const Vector3D<S> &v) + { + x /= v.x; + y /= v.y; + z /= v.z; + return *this; + } + //! Assign and div operator + inline const Vector3D<S> &operator/=(S s) + { + x /= s; + y /= s; + z /= s; + return *this; + } + //! Negation operator + inline Vector3D<S> operator-() const + { + return Vector3D<S>(-x, -y, -z); + } + + //! 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; + } + + //! 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[3]; + struct { + S x; + S y; + S z; + }; + struct { + S X; + S Y; + S Z; + }; + }; + + //! zero element + static const Vector3D<S> Zero, Invalid; + + //! For compatibility with 4d vectors (discards 4th comp) + inline Vector3D(S vx, S vy, S vz, S vDummy) : x(vx), y(vy), z(vz) + { + } + + protected: +}; + +//! helper to check whether float/double value is non-zero +inline bool notZero(Real f) +{ + if (std::abs(f) > VECTOR_EPSILON) + return true; + return false; +} + +//************************************************************************ +// Additional operators +//************************************************************************ + +//! Addition operator +template<class S> inline Vector3D<S> operator+(const Vector3D<S> &v1, const Vector3D<S> &v2) +{ + return Vector3D<S>(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z); +} +//! Addition operator +template<class S, class S2> inline Vector3D<S> operator+(const Vector3D<S> &v, S2 s) +{ + return Vector3D<S>(v.x + s, v.y + s, v.z + s); +} +//! Addition operator +template<class S, class S2> inline Vector3D<S> operator+(S2 s, const Vector3D<S> &v) +{ + return Vector3D<S>(v.x + s, v.y + s, v.z + s); +} + +//! Subtraction operator +template<class S> inline Vector3D<S> operator-(const Vector3D<S> &v1, const Vector3D<S> &v2) +{ + return Vector3D<S>(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z); +} +//! Subtraction operator +template<class S, class S2> inline Vector3D<S> operator-(const Vector3D<S> &v, S2 s) +{ + return Vector3D<S>(v.x - s, v.y - s, v.z - s); +} +//! Subtraction operator +template<class S, class S2> inline Vector3D<S> operator-(S2 s, const Vector3D<S> &v) +{ + return Vector3D<S>(s - v.x, s - v.y, s - v.z); +} + +//! Multiplication operator +template<class S> inline Vector3D<S> operator*(const Vector3D<S> &v1, const Vector3D<S> &v2) +{ + return Vector3D<S>(v1.x * v2.x, v1.y * v2.y, v1.z * v2.z); +} +//! Multiplication operator +template<class S, class S2> inline Vector3D<S> operator*(const Vector3D<S> &v, S2 s) +{ + return Vector3D<S>(v.x * s, v.y * s, v.z * s); +} +//! Multiplication operator +template<class S, class S2> inline Vector3D<S> operator*(S2 s, const Vector3D<S> &v) +{ + return Vector3D<S>(s * v.x, s * v.y, s * v.z); +} + +//! Division operator +template<class S> inline Vector3D<S> operator/(const Vector3D<S> &v1, const Vector3D<S> &v2) +{ + return Vector3D<S>(v1.x / v2.x, v1.y / v2.y, v1.z / v2.z); +} +//! Division operator +template<class S, class S2> inline Vector3D<S> operator/(const Vector3D<S> &v, S2 s) +{ + return Vector3D<S>(v.x / s, v.y / s, v.z / s); +} +//! Division operator +template<class S, class S2> inline Vector3D<S> operator/(S2 s, const Vector3D<S> &v) +{ + return Vector3D<S>(s / v.x, s / v.y, s / v.z); +} + +//! Comparison operator +template<class S> inline bool operator==(const Vector3D<S> &s1, const Vector3D<S> &s2) +{ + return s1.x == s2.x && s1.y == s2.y && s1.z == s2.z; +} + +//! Comparison operator +template<class S> inline bool operator!=(const Vector3D<S> &s1, const Vector3D<S> &s2) +{ + return s1.x != s2.x || s1.y != s2.y || s1.z != s2.z; +} + +//************************************************************************ +// External functions +//************************************************************************ + +//! Min operator +template<class S> inline Vector3D<S> vmin(const Vector3D<S> &s1, const Vector3D<S> &s2) +{ + return Vector3D<S>(std::min(s1.x, s2.x), std::min(s1.y, s2.y), std::min(s1.z, s2.z)); +} + +//! Min operator +template<class S, class S2> inline Vector3D<S> vmin(const Vector3D<S> &s1, S2 s2) +{ + return Vector3D<S>(std::min(s1.x, s2), std::min(s1.y, s2), std::min(s1.z, s2)); +} + +//! Min operator +template<class S1, class S> inline Vector3D<S> vmin(S1 s1, const Vector3D<S> &s2) +{ + return Vector3D<S>(std::min(s1, s2.x), std::min(s1, s2.y), std::min(s1, s2.z)); +} + +//! Max operator +template<class S> inline Vector3D<S> vmax(const Vector3D<S> &s1, const Vector3D<S> &s2) +{ + return Vector3D<S>(std::max(s1.x, s2.x), std::max(s1.y, s2.y), std::max(s1.z, s2.z)); +} + +//! Max operator +template<class S, class S2> inline Vector3D<S> vmax(const Vector3D<S> &s1, S2 s2) +{ + return Vector3D<S>(std::max(s1.x, s2), std::max(s1.y, s2), std::max(s1.z, s2)); +} + +//! Max operator +template<class S1, class S> inline Vector3D<S> vmax(S1 s1, const Vector3D<S> &s2) +{ + return Vector3D<S>(std::max(s1, s2.x), std::max(s1, s2.y), std::max(s1, s2.z)); +} + +//! Dot product +template<class S> inline S dot(const Vector3D<S> &t, const Vector3D<S> &v) +{ + return t.x * v.x + t.y * v.y + t.z * v.z; +} + +//! Cross product +template<class S> inline Vector3D<S> cross(const Vector3D<S> &t, const Vector3D<S> &v) +{ + Vector3D<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; +} + +//! Project a vector into a plane, defined by its normal +/*! Projects a vector into a plane normal to the given vector, which must + have unit length. Self is modified. + \param v The vector to project + \param n The plane normal + \return The projected vector */ +template<class S> +inline const Vector3D<S> &projectNormalTo(const Vector3D<S> &v, const Vector3D<S> &n) +{ + S sprod = dot(v, n); + return v - n * dot(v, n); +} + +//! Compute the magnitude (length) of the vector +//! (clamps to 0 and 1 with VECTOR_EPSILON) +template<class S> inline S norm(const Vector3D<S> &v) +{ + S l = v.x * v.x + v.y * v.y + v.z * v.z; + if (l <= VECTOR_EPSILON * VECTOR_EPSILON) + return (0.); + return (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) ? 1. : sqrt(l); +} + +//! Compute squared magnitude +template<class S> inline S normSquare(const Vector3D<S> &v) +{ + return v.x * v.x + v.y * v.y + v.z * v.z; +} + +//! compatibility, allow use of int, Real and Vec inputs with norm/normSquare +inline Real norm(const Real v) +{ + return fabs(v); +} +inline Real normSquare(const Real v) +{ + return square(v); +} +inline Real norm(const int v) +{ + return abs(v); +} +inline Real normSquare(const int v) +{ + return square(v); +} + +//! Returns a normalized vector +template<class S> inline Vector3D<S> getNormalized(const Vector3D<S> &v) +{ + S l = v.x * v.x + v.y * v.y + v.z * v.z; + 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 Vector3D<S>(v.x * fac, v.y * fac, v.z * fac); + } + else + return Vector3D<S>((S)0); +} + +//! Compute the norm of the vector and normalize it. +/*! \return The value of the norm */ +template<class S> inline S normalize(Vector3D<S> &v) +{ + S norm; + S l = v.x * v.x + v.y * v.y + v.z * v.z; + 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 = Vector3D<S>::Zero; + norm = 0.; + } + return (S)norm; +} + +//! Obtain an orthogonal vector +/*! Compute a vector that is orthonormal to the given vector. + * Nothing else can be assumed for the direction of the new vector. + * \return The orthonormal vector */ +template<class S> Vector3D<S> getOrthogonalVector(const Vector3D<S> &v) +{ + // Determine the component with max. absolute value + int maxIndex = (fabs(v.x) > fabs(v.y)) ? 0 : 1; + maxIndex = (fabs(v[maxIndex]) > fabs(v.z)) ? maxIndex : 2; + + // Choose another axis than the one with max. component and project + // orthogonal to self + Vector3D<S> o(0.0); + o[(maxIndex + 1) % 3] = 1; + + Vector3D<S> c = cross(v, o); + normalize(c); + return c; +} + +//! Convert vector to polar coordinates +/*! Stable vector to angle conversion + *\param v vector to convert + \param phi unique angle [0,2PI] + \param theta unique angle [0,PI] + */ +template<class S> inline void vecToAngle(const Vector3D<S> &v, S &phi, S &theta) +{ + if (fabs(v.y) < VECTOR_EPSILON) + theta = M_PI / 2; + else if (fabs(v.x) < VECTOR_EPSILON && fabs(v.z) < VECTOR_EPSILON) + theta = (v.y >= 0) ? 0 : M_PI; + else + theta = atan(sqrt(v.x * v.x + v.z * v.z) / v.y); + if (theta < 0) + theta += M_PI; + + if (fabs(v.x) < VECTOR_EPSILON) + phi = M_PI / 2; + else + phi = atan(v.z / v.x); + if (phi < 0) + phi += M_PI; + if (fabs(v.z) < VECTOR_EPSILON) + phi = (v.x >= 0) ? 0 : M_PI; + else if (v.z < 0) + phi += M_PI; +} + +//! Compute vector reflected at a surface +/*! Compute a vector, that is self (as an incoming vector) + * reflected at a surface with a distinct normal vector. + * Note that the normal is reversed, if the scalar product with it is positive. + \param t The incoming vector + \param n The surface normal + \return The new reflected vector + */ +template<class S> inline Vector3D<S> reflectVector(const Vector3D<S> &t, const Vector3D<S> &n) +{ + Vector3D<S> nn = (dot(t, n) > 0.0) ? (n * -1.0) : n; + return (t - nn * (2.0 * dot(nn, t))); +} + +//! Compute vector refracted at a surface +/*! \param t The incoming vector + * \param n The surface normal + * \param nt The "inside" refraction index + * \param nair The "outside" refraction index + * \param refRefl Set to 1 on total reflection + * \return The refracted vector + */ +template<class S> +inline Vector3D<S> refractVector( + const Vector3D<S> &t, const Vector3D<S> &normal, S nt, S nair, int &refRefl) +{ + // from Glassner's book, section 5.2 (Heckberts method) + S eta = nair / nt; + S n = -dot(t, normal); + S tt = 1.0 + eta * eta * (n * n - 1.0); + if (tt < 0.0) { + // we have total reflection! + refRefl = 1; + } + else { + // normal reflection + tt = eta * n - sqrt(tt); + return (t * eta + normal * tt); + } + return t; +} + +//! Outputs the object in human readable form as string +template<class S> std::string Vector3D<S>::toString() const +{ + char buf[256]; + snprintf(buf, + 256, + "[%+4.6f,%+4.6f,%+4.6f]", + (double)(*this)[0], + (double)(*this)[1], + (double)(*this)[2]); + // for debugging, optionally increase precision: + // snprintf ( buf,256,"[%+4.16f,%+4.16f,%+4.16f]", ( double ) ( *this ) [0], ( double ) ( *this ) + // [1], ( double ) ( *this ) [2] ); + return std::string(buf); +} + +template<> std::string Vector3D<int>::toString() const; + +//! Outputs the object in human readable form to stream +/*! Output format [x,y,z] */ +template<class S> std::ostream &operator<<(std::ostream &os, const Vector3D<S> &i) +{ + os << i.toString(); + return os; +} + +//! Reads the contents of the object from a stream +/*! Input format [x,y,z] */ +template<class S> std::istream &operator>>(std::istream &is, Vector3D<S> &i) +{ + char c; + char dummy[3]; + is >> c >> i[0] >> dummy >> i[1] >> dummy >> i[2] >> c; + return is; +} + +/**************************************************************************/ +// Define default vector alias +/**************************************************************************/ + +//! 3D vector class of type Real (typically float) +typedef Vector3D<Real> Vec3; + +//! 3D vector class of type int +typedef Vector3D<int> Vec3i; + +//! convert to Real Vector +template<class T> inline Vec3 toVec3(T v) +{ + return Vec3(v[0], v[1], v[2]); +} + +//! convert to int Vector +template<class T> inline Vec3i toVec3i(T v) +{ + return Vec3i((int)v[0], (int)v[1], (int)v[2]); +} + +//! convert to int Vector +template<class T> inline Vec3i toVec3i(T v0, T v1, T v2) +{ + return Vec3i((int)v0, (int)v1, (int)v2); +} + +//! round, and convert to int Vector +template<class T> inline Vec3i toVec3iRound(T v) +{ + return Vec3i((int)round(v[0]), (int)round(v[1]), (int)round(v[2])); +} + +//! convert to int Vector if values are close enough to an int +template<class T> inline Vec3i toVec3iChecked(T v) +{ + Vec3i ret; + for (size_t i = 0; i < 3; i++) { + Real a = v[i]; + if (fabs(a - floor(a + 0.5)) > 1e-5) + errMsg("argument is not an int, cannot convert"); + ret[i] = (int)(a + 0.5); + } + return ret; +} + +//! convert to double Vector +template<class T> inline Vector3D<double> toVec3d(T v) +{ + return Vector3D<double>(v[0], v[1], v[2]); +} + +//! convert to float Vector +template<class T> inline Vector3D<float> toVec3f(T v) +{ + return Vector3D<float>(v[0], v[1], v[2]); +} + +/**************************************************************************/ +// Specializations for common math functions +/**************************************************************************/ + +template<> inline Vec3 clamp<Vec3>(const Vec3 &a, const Vec3 &b, const Vec3 &c) +{ + return Vec3(clamp(a.x, b.x, c.x), clamp(a.y, b.y, c.y), clamp(a.z, b.z, c.z)); +} +template<> inline Vec3 safeDivide<Vec3>(const Vec3 &a, const Vec3 &b) +{ + return Vec3(safeDivide(a.x, b.x), safeDivide(a.y, b.y), safeDivide(a.z, b.z)); +} +template<> inline Vec3 nmod<Vec3>(const Vec3 &a, const Vec3 &b) +{ + return Vec3(nmod(a.x, b.x), nmod(a.y, b.y), nmod(a.z, b.z)); +} + +}; // namespace Manta + +#endif |