MOAB: Mesh Oriented datABase  (version 5.2.1)
Vector3D.hpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 #ifndef MESQUITE_VECTOR3D_HPP
00028 #define MESQUITE_VECTOR3D_HPP
00029 
00030 #include "Mesquite.hpp"
00031 
00032 #include <iosfwd>
00033 #include <cassert>
00034 #include <cstring>
00035 #include <vector>
00036 
00037 /*! \file Vector3D.hpp
00038   \brief Vector object with exactly 3 dimensions.
00039 
00040   This is as fast as a C array.
00041 
00042   \author Darryl Melander
00043   \author Thomas Leurent
00044 */
00045 namespace MBMesquite
00046 {
00047 class Matrix3D;
00048 class MsqError;
00049 
00050 /*!
00051    \class Vector3D
00052    \brief Vector3D is the object that effeciently stores information about
00053    about three-deminsional vectors.  It is also the parent class of
00054    MsqVertex.      */
00055 class MESQUITE_EXPORT Vector3D
00056 {
00057   public:
00058     // Constructors
00059     Vector3D();
00060     Vector3D( double xyz );
00061     Vector3D( double x, double y, double z );
00062     Vector3D( const double xyz[3] );
00063     Vector3D( const Vector3D& to_copy );
00064 
00065     // *** virtual destructor *** Do not use for Vector3D, we need to keep
00066     // the deallocation of those objects very fast
00067 
00068     // Functions to get the coordinates
00069     double x() const;
00070     double y() const;
00071     double z() const;
00072     void get_coordinates( double& x, double& y, double& z ) const;
00073     void get_coordinates( double xyz[3] ) const;
00074     const double& operator[]( size_t index ) const;  // 0-based
00075 
00076     // Functions to set the coordinates.
00077     void x( const double x );
00078     void y( const double y );
00079     void z( const double z );
00080     void set( const double x, const double y, const double z );
00081     void set( const double xyz[3] );
00082     void set( const Vector3D& to_copy );
00083     // Subscripts on non-consts both get and set coords
00084     double& operator[]( size_t index );  // 0-based
00085     Vector3D& operator=( const Vector3D& to_copy );
00086     Vector3D& operator=( const double& to_copy );
00087 
00088     // Functions to modify existing coordinates
00089     Vector3D operator-() const;  //- unary negation.
00090     Vector3D& operator*=( const double scalar );
00091     Vector3D& operator/=( const double scalar );
00092     Vector3D& operator*=( const Vector3D& rhs );  //- cross product
00093     Vector3D& operator+=( const Vector3D& rhs );
00094     Vector3D& operator-=( const Vector3D& rhs );
00095 
00096     // Binary operators (like a+b).
00097     friend const Vector3D operator+( const Vector3D& lhs, const Vector3D& rhs );
00098     friend const Vector3D operator-( const Vector3D& lhs, const Vector3D& rhs );
00099     friend const Vector3D operator*( const Vector3D& lhs,
00100                                      const double scalar );  //!< lhs * scalar
00101     friend const Vector3D operator*( const double scalar,
00102                                      const Vector3D& rhs );  //!< scalar * rhs
00103     friend const Vector3D operator/( const Vector3D& lhs,
00104                                      const double scalar );  //- lhs / scalar
00105     friend double operator%( const Vector3D& v1,
00106                              const Vector3D& v2 );  //!< dot product
00107     friend double inner( const Vector3D v1[], const Vector3D v2[],
00108                          int n );  //!< dot product for array
00109     friend double operator%( const double scalar,
00110                              const Vector3D& v2 );  //!< scalar * sum_i v2[i]
00111     friend double operator%( const Vector3D& v1,
00112                              const double scalar );  //!< scalar * sum_i v1[i]
00113     friend const Vector3D operator*( const Vector3D& v1,
00114                                      const Vector3D& v2 );  //!< cross product
00115 
00116     //! \f$ v = A*x \f$
00117     friend void eqAx( Vector3D& v, const Matrix3D& A, const Vector3D& x );
00118     //! \f$ v += A*x \f$
00119     friend void plusEqAx( Vector3D& v, const Matrix3D& A, const Vector3D& x );
00120     //! \f$ v += A^T*x \f$
00121     friend void plusEqTransAx( Vector3D& v, const Matrix3D& A, const Vector3D& x );
00122     friend void eqTransAx( Vector3D& v, const Matrix3D& A, const Vector3D& x );
00123 
00124     // Comparison functions
00125     friend bool operator==( const Vector3D& lhs, const Vector3D& rhs );
00126     friend bool operator!=( const Vector3D& lhs, const Vector3D& rhs );
00127     static double distance_between( const Vector3D& p1, const Vector3D& p2 );
00128     int within_tolerance_box( const Vector3D& compare_to, double tolerance ) const;
00129     //- Compare two Vector3Ds to see if they are spatially equal.
00130     // Return TRUE if difference in x, y, and z are all within tolerance.
00131     // Essentially checks to see if 'this' lies within a box centered on
00132     // 'compare_to' with sides of length ('tolerance' * 2).
00133 
00134     // Length functions
00135     inline double length_squared() const;
00136     inline double length() const;
00137     friend double length( const Vector3D* v, int n );  //!< L2 norm for an array of Vector3Ds
00138     friend double Linf( const Vector3D* v, int n );    //!< L inf norm for array of Vector3Ds
00139 
00140     inline void set_length( const double new_length );
00141     inline void normalize();
00142     Vector3D operator~() const
00143     {
00144         return *this * ( 1.0 / length() );
00145     }
00146 
00147     // Utility functions.  All angle functions work in radians.
00148     static double interior_angle( const Vector3D& a, const Vector3D& b, MsqError& err );
00149     //- Interior angle: acos((a%b)/(|a||b|))
00150     static Vector3D interpolate( const double param, const Vector3D& p1, const Vector3D& p2 );
00151     //- Interpolate between two points. Returns (1-param)*v1 + param*v2.
00152 
00153     const double* to_array() const
00154     {
00155         return mCoords;
00156     }
00157 
00158     double* to_array()
00159     {
00160         return mCoords;
00161     }
00162 
00163   protected:
00164     double mCoords[3];
00165 };
00166 
00167 // Constructors
00168 inline Vector3D::Vector3D()
00169 {
00170     mCoords[0] = 0;
00171     mCoords[1] = 0;
00172     mCoords[2] = 0;
00173 }
00174 inline Vector3D::Vector3D( double p_x )
00175 {
00176     mCoords[0] = p_x;
00177     mCoords[1] = p_x;
00178     mCoords[2] = p_x;
00179 }
00180 inline Vector3D::Vector3D( double p_x, double p_y, double p_z )
00181 {
00182     mCoords[0] = p_x;
00183     mCoords[1] = p_y;
00184     mCoords[2] = p_z;
00185 }
00186 inline Vector3D::Vector3D( const double xyz[3] )
00187 {
00188     std::memcpy( mCoords, xyz, 3 * sizeof( double ) );
00189 }
00190 inline Vector3D::Vector3D( const Vector3D& to_copy )
00191 {
00192     std::memcpy( mCoords, to_copy.mCoords, 3 * sizeof( double ) );
00193 }
00194 
00195 // Functions to get coordinates
00196 inline double Vector3D::x() const
00197 {
00198     return mCoords[0];
00199 }
00200 inline double Vector3D::y() const
00201 {
00202     return mCoords[1];
00203 }
00204 inline double Vector3D::z() const
00205 {
00206     return mCoords[2];
00207 }
00208 inline void Vector3D::get_coordinates( double& p_x, double& p_y, double& p_z ) const
00209 {
00210     p_x = mCoords[0];
00211     p_y = mCoords[1];
00212     p_z = mCoords[2];
00213 }
00214 inline void Vector3D::get_coordinates( double xyz[3] ) const
00215 {
00216     std::memcpy( xyz, mCoords, 3 * sizeof( double ) );
00217 }
00218 inline const double& Vector3D::operator[]( size_t index ) const
00219 {
00220     return mCoords[index];
00221 }
00222 
00223 // Functions to set coordinates
00224 inline void Vector3D::x( const double p_x )
00225 {
00226     mCoords[0] = p_x;
00227 }
00228 inline void Vector3D::y( const double p_y )
00229 {
00230     mCoords[1] = p_y;
00231 }
00232 inline void Vector3D::z( const double p_z )
00233 {
00234     mCoords[2] = p_z;
00235 }
00236 inline void Vector3D::set( const double p_x, const double p_y, const double p_z )
00237 {
00238     mCoords[0] = p_x;
00239     mCoords[1] = p_y;
00240     mCoords[2] = p_z;
00241 }
00242 inline void Vector3D::set( const double xyz[3] )
00243 {
00244     std::memcpy( mCoords, xyz, 3 * sizeof( double ) );
00245 }
00246 inline void Vector3D::set( const Vector3D& to_copy )
00247 {
00248     std::memcpy( mCoords, to_copy.mCoords, 3 * sizeof( double ) );
00249 }
00250 inline double& Vector3D::operator[]( size_t index )
00251 {
00252     return mCoords[index];
00253 }
00254 
00255 inline Vector3D& Vector3D::operator=( const Vector3D& to_copy )
00256 {
00257     mCoords[0] = to_copy.mCoords[0];
00258     mCoords[1] = to_copy.mCoords[1];
00259     mCoords[2] = to_copy.mCoords[2];
00260     //    memcpy(mCoords, to_copy.mCoords, 3*sizeof(double));
00261     return *this;
00262 }
00263 
00264 inline Vector3D& Vector3D::operator=( const double& to_copy )
00265 {
00266     mCoords[0] = to_copy;
00267     mCoords[1] = to_copy;
00268     mCoords[2] = to_copy;
00269     return *this;
00270 }
00271 
00272 // Functions that modify existing coordinates
00273 inline Vector3D Vector3D::operator-() const
00274 {
00275     return Vector3D( -mCoords[0], -mCoords[1], -mCoords[2] );
00276 }
00277 inline Vector3D& Vector3D::operator*=( const double scalar )
00278 {
00279     mCoords[0] *= scalar;
00280     mCoords[1] *= scalar;
00281     mCoords[2] *= scalar;
00282     return *this;
00283 }
00284 //! divides each Vector3D entry by the given scalar.
00285 inline Vector3D& Vector3D::operator/=( const double scalar )
00286 {
00287     mCoords[0] /= scalar;
00288     mCoords[1] /= scalar;
00289     mCoords[2] /= scalar;
00290     return *this;
00291 }
00292 inline Vector3D& Vector3D::operator*=( const Vector3D& rhs )
00293 {
00294     double new_coords[3] = { mCoords[1] * rhs.mCoords[2] - mCoords[2] * rhs.mCoords[1],
00295                              mCoords[2] * rhs.mCoords[0] - mCoords[0] * rhs.mCoords[2],
00296                              mCoords[0] * rhs.mCoords[1] - mCoords[1] * rhs.mCoords[0] };
00297     std::memcpy( mCoords, new_coords, 3 * sizeof( double ) );
00298     return *this;
00299 }
00300 inline Vector3D& Vector3D::operator+=( const Vector3D& rhs )
00301 {
00302     mCoords[0] += rhs.mCoords[0];
00303     mCoords[1] += rhs.mCoords[1];
00304     mCoords[2] += rhs.mCoords[2];
00305     return *this;
00306 }
00307 inline Vector3D& Vector3D::operator-=( const Vector3D& rhs )
00308 {
00309     mCoords[0] -= rhs.mCoords[0];
00310     mCoords[1] -= rhs.mCoords[1];
00311     mCoords[2] -= rhs.mCoords[2];
00312     return *this;
00313 }
00314 
00315 // Binary operators
00316 inline const Vector3D operator+( const Vector3D& lhs, const Vector3D& rhs )
00317 {
00318     return Vector3D( lhs.x() + rhs.x(), lhs.y() + rhs.y(), lhs.z() + rhs.z() );
00319 }
00320 inline const Vector3D operator-( const Vector3D& lhs, const Vector3D& rhs )
00321 {
00322     return Vector3D( lhs.x() - rhs.x(), lhs.y() - rhs.y(), lhs.z() - rhs.z() );
00323 }
00324 inline const Vector3D operator*( const Vector3D& lhs, const double scalar )
00325 {
00326     return Vector3D( lhs.x() * scalar, lhs.y() * scalar, lhs.z() * scalar );
00327 }
00328 inline const Vector3D operator*( const double scalar, const Vector3D& rhs )
00329 {
00330     return Vector3D( rhs.x() * scalar, rhs.y() * scalar, rhs.z() * scalar );
00331 }
00332 inline const Vector3D operator/( const Vector3D& lhs, const double scalar )
00333 {
00334     assert( scalar != 0 );
00335     return Vector3D( lhs.x() / scalar, lhs.y() / scalar, lhs.z() / scalar );
00336 }
00337 inline double operator%( const Vector3D& lhs,
00338                          const Vector3D& rhs )  // Dot Product
00339 {
00340     return ( lhs.mCoords[0] * rhs.mCoords[0] + lhs.mCoords[1] * rhs.mCoords[1] + lhs.mCoords[2] * rhs.mCoords[2] );
00341 }
00342 
00343 /*! Dot product for arrays of Vector3Ds. see also operator% .*/
00344 inline double inner( const Vector3D lhs[], const Vector3D rhs[], int n )
00345 {
00346     int i;
00347     double dot = 0;
00348     for( i = 0; i < n; ++i )
00349         dot += lhs[i].mCoords[0] * rhs[i].mCoords[0] + lhs[i].mCoords[1] * rhs[i].mCoords[1] +
00350                lhs[i].mCoords[2] * rhs[i].mCoords[2];
00351     return dot;
00352 }
00353 /*! Dot product for arrays of Vector3Ds. see also operator% .*/
00354 inline double inner( const std::vector< Vector3D >& lhs, const std::vector< Vector3D >& rhs )
00355 {
00356     double dot = 0;
00357     assert( lhs.size() == rhs.size() );
00358     for( size_t i = 0; i < lhs.size(); ++i )
00359         dot = lhs[i] % rhs[i];
00360     return dot;
00361 }
00362 
00363 inline double operator%( const double scalar,
00364                          const Vector3D& rhs )  // Dot Product
00365 {
00366     return ( scalar * ( rhs.mCoords[0] + rhs.mCoords[1] + rhs.mCoords[2] ) );
00367 }
00368 inline double operator%( const Vector3D& lhs,
00369                          const double scalar )  // Dot Product
00370 {
00371     return ( scalar * ( lhs.mCoords[0] + lhs.mCoords[1] + lhs.mCoords[2] ) );
00372 }
00373 inline const Vector3D operator*( const Vector3D& lhs,
00374                                  const Vector3D& rhs )  // Cross Product
00375 {
00376     return Vector3D( lhs.mCoords[1] * rhs.mCoords[2] - lhs.mCoords[2] * rhs.mCoords[1],
00377                      lhs.mCoords[2] * rhs.mCoords[0] - lhs.mCoords[0] * rhs.mCoords[2],
00378                      lhs.mCoords[0] * rhs.mCoords[1] - lhs.mCoords[1] * rhs.mCoords[0] );
00379 }
00380 
00381 // output operator
00382 MESQUITE_EXPORT std::ostream& operator<<( std::ostream& s, const MBMesquite::Vector3D& v );
00383 
00384 inline double Vector3D::distance_between( const Vector3D& p1, const Vector3D& p2 )
00385 {
00386     Vector3D v = p2 - p1;
00387     return v.length();
00388 }
00389 inline int Vector3D::within_tolerance_box( const Vector3D& compare_to, double tolerance ) const
00390 {
00391     return ( ( std::fabs( this->mCoords[0] - compare_to.mCoords[0] ) < tolerance ) &&
00392              ( std::fabs( this->mCoords[1] - compare_to.mCoords[1] ) < tolerance ) &&
00393              ( std::fabs( this->mCoords[2] - compare_to.mCoords[2] ) < tolerance ) );
00394 }
00395 
00396 // Length functions
00397 inline double Vector3D::length_squared() const
00398 {
00399     return ( mCoords[0] * mCoords[0] + mCoords[1] * mCoords[1] + mCoords[2] * mCoords[2] );
00400 }
00401 inline double Vector3D::length() const
00402 {
00403     return std::sqrt( mCoords[0] * mCoords[0] + mCoords[1] * mCoords[1] + mCoords[2] * mCoords[2] );
00404 }
00405 
00406 inline double inner_product( const Vector3D* v1, const Vector3D* v2, size_t n )
00407 {
00408     double result             = 0.0;
00409     const Vector3D* const end = v1 + n;
00410     while( v1 < end )
00411     {
00412         result += *v1 % *v2;
00413         ++v1;
00414         ++v2;
00415     }
00416     return result;
00417 }
00418 
00419 inline double length_squared( const Vector3D* v, int n )
00420 {
00421     double sum = 0.0;
00422     for( int i = 0; i < n; ++i )
00423         sum += v[i].length_squared();
00424     return sum;
00425 }
00426 inline double length_squared( const std::vector< Vector3D >& v )
00427 {
00428     double sum = 0.0;
00429     for( size_t i = 0; i < v.size(); ++i )
00430         sum += v[i].length_squared();
00431     return sum;
00432 }
00433 
00434 inline double length( const Vector3D* v, int n )  // norm for an array of Vector3Ds
00435 {
00436     return std::sqrt( length_squared( v, n ) );
00437 }
00438 inline double length( const std::vector< Vector3D >& v )
00439 {
00440     return std::sqrt( length_squared( v ) );
00441 }
00442 
00443 inline double Linf( const Vector3D* v, int n )  // max entry for an array of Vector3Ds
00444 {
00445     double max = 0;
00446     // loop over the length of the array
00447     for( int i = 0; i < n; ++i )
00448     {
00449         if( max < std::fabs( v[i][0] ) ) max = std::fabs( v[i][0] );
00450         if( max < std::fabs( v[i][1] ) ) max = std::fabs( v[i][1] );
00451         if( max < std::fabs( v[i][2] ) ) max = std::fabs( v[i][2] );
00452     }
00453     // return the value of the largest entry in the array
00454     return max;
00455 }
00456 
00457 inline double Linf( const std::vector< Vector3D >& v )  // max entry for an array of Vector3Ds
00458 {
00459     double max = 0;
00460     // loop over the length of the array
00461     for( size_t i = 0; i < v.size(); ++i )
00462     {
00463         if( max < std::fabs( v[i][0] ) ) max = std::fabs( v[i][0] );
00464         if( max < std::fabs( v[i][1] ) ) max = std::fabs( v[i][1] );
00465         if( max < std::fabs( v[i][2] ) ) max = std::fabs( v[i][2] );
00466     }
00467     // return the value of the largest entry in the array
00468     return max;
00469 }
00470 
00471 inline void Vector3D::set_length( const double new_length )
00472 {
00473     double factor = new_length / length();
00474     *this *= factor;
00475 }
00476 inline void Vector3D::normalize()
00477 {
00478     set_length( 1.0 );
00479 }
00480 
00481 // Utility functions.
00482 inline Vector3D Vector3D::interpolate( const double param, const Vector3D& p1, const Vector3D& p2 )
00483 {
00484     return ( 1 - param ) * p1 + param * p2;
00485 }
00486 
00487 inline bool operator==( const Vector3D& v1, const Vector3D& v2 )
00488 {
00489     return v1.mCoords[0] == v2.mCoords[0] && v1.mCoords[1] == v2.mCoords[1] && v1.mCoords[2] == v2.mCoords[2];
00490 }
00491 
00492 inline bool operator!=( const Vector3D& v1, const Vector3D& v2 )
00493 {
00494     return v1.mCoords[0] != v2.mCoords[0] || v1.mCoords[1] != v2.mCoords[1] || v1.mCoords[2] != v2.mCoords[2];
00495 }
00496 
00497 }  // namespace MBMesquite
00498 
00499 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines