cgma
CubitVector.hpp
Go to the documentation of this file.
00001 //- Class: CubitVector
00002 //-
00003 //- Description: This file defines the CubitVector class which is a
00004 //- standard three-dimensional vector. All relevant arithmetic
00005 //- operators are overloaded so CubitVectors can be used similar to
00006 //- built-in types.
00007 //-
00008 //- Owner: Greg Sjaardema
00009 //- Checked by: Randy Lober, January 94
00010 //- Version: $Id: 
00011 
00012 #ifndef CUBITVECTOR_HPP
00013 #define CUBITVECTOR_HPP
00014 
00015 #include "CubitDefines.h"
00016 #include "CubitVectorStruct.h"
00017 #include "CGMUtilConfigure.h"
00018 
00019 class CubitVector;
00020 typedef void ( CubitVector::*transform_function )( double gamma,
00021                                                    double gamma2);
00022 // a pointer to some function that transforms the point,
00023 // taking a double parameter.  e.g. blow_out, rotate, and scale_angle
00024 
00025 class CUBIT_UTIL_EXPORT CubitVector
00026 {
00027 public:
00028   
00029     //- Heading: Constructors and Destructor
00030   CubitVector();  //- Default constructor.
00031   
00032   explicit CubitVector(const double x, const double y, const double z);
00033     //- Constructor: create vector from three components
00034   
00035   explicit CubitVector( const double xyz[3] );
00036     //- Constructor: create vector from tuple
00037 
00038   explicit CubitVector (const CubitVector& tail, const CubitVector& head);
00039     //- Constructor for a CubitVector starting at tail and pointing
00040     //- to head.
00041   
00042   CubitVector(const CubitVector& copy_from);  //- Copy Constructor
00043   
00044   explicit CubitVector(const CubitVectorStruct& from);
00045   
00046     //- Heading: Set and Inquire Functions
00047   void set(const double x, const double y, const double z);
00048     //- Change vector components to {x}, {y}, and {z}
00049   
00050   void set( const double xyz[3] );
00051     //- Change vector components to xyz[0], xyz[1], xyz[2]
00052 
00053   void set(const CubitVector& tail, const CubitVector& head);
00054     //- Change vector to go from tail to head.
00055   
00056   void set(const CubitVector& to_copy);
00057     //- Same as i.=(const CubitVector&)
00058   
00059   double x() const; //- Return x component of vector
00060   
00061   double y() const; //- Return y component of vector
00062   
00063   double z() const; //- Return z component of vector
00064   
00065   double& x();
00066   double& y();
00067   double& z();
00068   
00069   void get_xyz( double &x, double &y, double &z ) const; //- Get x, y, z components
00070   void get_xyz( double xyz[3] ) const; //- Get xyz tuple
00071   
00072   double &r(); //- Return r component of vector, if (r,theta) format
00073   
00074   double &theta();  //- Return theta component of vector, if (r,theta) format
00075   
00076   void x( const double x ); //- Set x component of vector
00077   
00078   void y( const double y ); //- Set y component of vector
00079   
00080   void z( const double z ); //- Set z component of vector
00081   
00082   void r( const double x ); //- Set r component of vector, if (r,theta) format
00083   
00084   void theta( const double y ); //- Set theta component of vector, if (r,theta) format
00085   
00086   void xy_to_rtheta();
00087     //- convert from cartesian to polar coordinates, just 2d for now
00088     //- theta is in [0,2 PI)
00089   
00090   void rtheta_to_xy();
00091     //- convert from  polar to cartesian coordinates, just 2d for now
00092   
00093   void scale_angle(double gamma, double );
00094     //- tranform_function.
00095     //- transform  (x,y) to (r,theta) to (r,gamma*theta) to (x',y')
00096     //- plus some additional scaling so long chords won't cross short ones
00097   
00098   void blow_out(double gamma, double gamma2 = 0.0);
00099     //- transform_function
00100     //- blow points radially away from the origin, 
00101   
00102   void rotate(double angle, double );
00103     //- transform function.
00104     //- transform  (x,y) to (r,theta) to (r,theta+angle) to (x',y')
00105 
00106   void project_to_plane( const CubitVector &planenormal );
00107     //- project this vector onto a plane.
00108   
00109   void project_to_line_segment( const CubitVector &pt0, const CubitVector &pt1 );
00110     //- project this vector onto a line segment.
00111   
00112   void reflect_about_xaxis(double dummy, double );
00113     //- dummy argument to make it a transform function
00114   
00115   double normalize();
00116     //- Normalize (set magnitude equal to 1) vector - return the magnitude
00117 
00118   CubitVector& length(const double new_length);
00119     //- Change length of vector to {new_length}. Can be used to move a
00120     //- location a specified distance from the origin in the current
00121     //- orientation.
00122   
00123   double length() const;
00124     //- Calculate the length of the vector.
00125     //- Use {length_squared()} if only comparing lengths, not adding.
00126   
00127   double distance_between(const CubitVector& test_vector) const;
00128     //- Calculate the distance from the head of one vector
00129     //  to the head of the test_vector.
00130   
00131   double distance_between_squared(const CubitVector& test_vector) const;
00132     //- Calculate the distance squared from the head of one vector
00133     //  to the head of the test_vector.
00134 
00135   double distance_from_infinite_line(const CubitVector& point_on_line,
00136                                      const CubitVector& line_direction) const;
00137     //- Calculate the minimum distance between the head of this vector and a 
00138     // line of infinite length.
00139   
00140   double distance_from_infinite_line_squared(const CubitVector& point_on_line,
00141                                      const CubitVector& line_direction) const;
00142     //- Calculate the square of the minimum distance between the head of this 
00143     //  vector and a line of infinite length.
00144   
00145   double length_squared() const;
00146     //- Calculate the squared length of the vector.
00147     //- Faster than {length()} since it eliminates the square root if
00148     //- only comparing other lengths.
00149   
00150   double interior_angle(const CubitVector &otherVector) const;
00151     //- Calculate the interior angle: acos((a%b)/(|a||b|))
00152     //- Returns angle in degrees.
00153   
00154   static bool colinear( const CubitVector &p0,
00155                         const CubitVector &p1,
00156                         const CubitVector &p2 );
00157     //- determine if 3 points are colinear.
00158 
00159   static CubitVector normal( const CubitVector &p0,
00160                              const CubitVector &p1,
00161                              const CubitVector &p2 );
00162     //- normal given 3 positions.
00163 
00164   static bool barycentric_coordinates( const CubitVector &v1,
00165                                        const CubitVector &v2,
00166                                        const CubitVector &v3,
00167                                        const CubitVector &point,
00168                                        double &coord_A,
00169                                        double &coord_B,
00170                                        double &coord_C );
00171 
00172   double vector_angle_quick(const CubitVector& vec1, const CubitVector& vec2);
00173     //- Calculate the interior angle between the projections of
00174     //- {vec1} and {vec2} onto the plane defined by the {this} vector.
00175     //- The angle returned is the right-handed angle around the {this}
00176     //- vector from {vec1} to {vec2}. Angle is in radians.
00177   
00178   double vector_angle( const CubitVector &vector1,
00179                        const CubitVector &vector2) const;
00180     //- Compute the angle between the projections of {vector1} and {vector2}
00181     //- onto the plane defined by *this. The angle is the
00182     //- right-hand angle, in radians, about *this from {vector1} to {vector2}.
00183   
00184   void perpendicular_z();
00185     //- Transform this vector to a perpendicular one, leaving
00186     //- z-component alone. Rotates clockwise about the z-axis by pi/2.
00187 
00188   void print_me();
00189     //- Prints out the coordinates of this vector.
00190   
00191   void orthogonal_vectors( CubitVector &vector2, CubitVector &vector3 ) const;
00192     //- Finds 2 (arbitrary) vectors that are orthogonal to this one
00193   
00194   void next_point( const CubitVector &direction, double distance, 
00195                    CubitVector& out_point ) const;
00196     //- Finds the next point in space based on *this* point (starting point), 
00197     //- a direction and the distance to extend in the direction. The
00198     //- direction vector need not be a unit vector.  The out_point can be
00199     //- "*this" (i.e., modify point in place).
00200   
00201   CubitBoolean about_equal( const CubitVector &w,
00202     const double relative_tolerance = 1.0e-6,
00203     const double absolute_tolerance = 1.0e-6 ) const;
00204     // Return true if vectors are equal within either relative tolerance
00205     // or absolute tolerance.
00206     //
00207     // More specifically:
00208     // Return true if the magnitude of the difference is less than
00209     // 1. absolute_tolerance, 
00210     // OR 
00211     // 2. relative_tolerance times the magnitude of the vectors
00212     //
00213     // E.g.
00214     // if v = <1, 1.0e-7, 0>, w = <1, -1.0e-7, 0> and relative_tol = 1.0e-6,
00215     // then return true.
00216 
00217   CubitBoolean within_tolerance(const CubitVector &vectorPtr2,
00218                                 double tolerance) const;
00219     //- Compare two vectors to see if they are spatially equal.  
00220     // Return TRUE if difference in x, y, and z are all within tolerance.
00221   
00222   CubitBoolean within_scaled_tolerance(const CubitVector &v2, double tol) const;
00223   // Return true if vectors are within_tolerance() or if, for each coord, 
00224   // the ratio of the two vector's coords differs no more than tol from 1.0.
00225   // Can be used to see if two vectors are equal up to a certain number of
00226   // significant digits.  For example, passing in a tolerance of 1e-6 will
00227   // return true if the two vectors are equal with 6 significant digits.
00228 
00229     //- Heading: Operator Overloads  *****************************
00230   CubitVector&  operator+=(const CubitVector &vec);
00231     //- Compound Assignment: addition: {this = this + vec}
00232   
00233   CubitVector& operator-=(const CubitVector &vec);
00234     //- Compound Assignment: subtraction: {this = this - vec}
00235   
00236   CubitVector& operator*=(const CubitVector &vec);
00237     //- Compound Assignment: cross product: {this = this * vec},
00238     //- non-commutative
00239   
00240   CubitVector& operator*=(const double scalar);
00241     //- Compound Assignment: multiplication: {this = this * scalar}
00242   
00243   CubitVector& operator/=(const double scalar);
00244     //- Compound Assignment: division: {this = this / scalar}
00245   
00246   CubitVector operator-() const;
00247     //- unary negation.
00248 
00249   double operator[](int i) const;
00250     //- return the ith value of the vector (x, y, z)
00251   
00252   friend CubitVector operator~(const CubitVector &vec);
00253     //- normalize. Returns a new vector which is a copy of {vec},
00254     //- scaled such that {|vec|=1}. Uses overloaded bitwise NOT operator.
00255   
00256   friend CubitVector operator+(const CubitVector &v1, 
00257                                const CubitVector &v2);
00258     //- vector addition
00259   
00260   friend CubitVector operator-(const CubitVector &v1, 
00261                                const CubitVector &v2);
00262     //- vector subtraction
00263   
00264   friend CubitVector operator*(const CubitVector &v1, 
00265                                const CubitVector &v2);
00266     //- vector cross product, non-commutative
00267   
00268   friend CubitVector operator*(const CubitVector &v1, const double sclr);
00269     //- vector * scalar
00270   
00271   friend CubitVector operator*(const double sclr, const CubitVector &v1);
00272     //- scalar * vector
00273   
00274   friend double operator%(const CubitVector &v1, const CubitVector &v2);
00275     //- dot product
00276   
00277   friend CubitVector operator/(const CubitVector &v1, const double sclr);
00278     //- vector / scalar
00279   
00280   friend int operator==(const CubitVector &v1, const CubitVector &v2);
00281     //- Equality operator
00282   
00283   friend int operator!=(const CubitVector &v1, const CubitVector &v2);
00284     //- Inequality operator
00285   
00286   friend CubitVector interpolate(const double param, const CubitVector &v1,
00287                                  const CubitVector &v2);
00288     //- Interpolate between two vectors. Returns (1-param)*v1 + param*v2
00289 
00290   CubitVector &operator=(const CubitVectorStruct &from);
00291 
00292   operator CubitVectorStruct() 
00293     {
00294       CubitVectorStruct to = {xVal, yVal, zVal};
00295       return to;
00296     }
00297 
00298   CubitVector &operator=(const CubitVector& from);
00299   
00300 private:
00301   
00302   double xVal;  //- x component of vector.
00303   double yVal;  //- y component of vector.
00304   double zVal;  //- z component of vector.
00305 
00306   static bool attempt_barycentric_coordinates_adjustment( const CubitVector &v1,
00307                                                           const CubitVector &v2,
00308                                                           const CubitVector &v3,
00309                                                           const CubitVector &point,
00310                                                           double &coord_A,
00311                                                           double &coord_B,
00312                                                           double &coord_C );
00313   static bool private_barycentric_coordinates( bool adjust_on_fail,
00314                                                const CubitVector &v1,
00315                                                const CubitVector &v2,
00316                                                const CubitVector &v3,
00317                                                const CubitVector &point,
00318                                                double &coord_A,
00319                                                double &coord_B,
00320                                                double &coord_C );
00321 
00322 };
00323 
00324 CUBIT_UTIL_EXPORT CubitVector vectorRotate(const double angle, 
00325                          const CubitVector &normalAxis,
00326                          const CubitVector &referenceAxis);
00327   //- A new coordinate system is created with the xy plane corresponding
00328   //- to the plane normal to {normalAxis}, and the x axis corresponding to
00329   //- the projection of {referenceAxis} onto the normal plane.  The normal
00330   //- plane is the tangent plane at the root point.  A unit vector is
00331   //- constructed along the local x axis and then rotated by the given
00332   //- ccw angle to form the new point.  The new point, then is a unit
00333   //- distance from the global origin in the tangent plane.
00334   //- {angle} is in radians.
00335 
00336 inline double CubitVector::x() const
00337 { return xVal; }
00338 inline double CubitVector::y() const
00339 { return yVal; }
00340 inline double CubitVector::z() const
00341 { return zVal; }
00342 inline double& CubitVector::x()
00343 { return xVal; }
00344 inline double& CubitVector::y()
00345 { return yVal; }
00346 inline double& CubitVector::z()
00347 { return zVal; }
00348 inline void CubitVector::get_xyz(double xyz[3]) const
00349 {
00350   xyz[0] = xVal;
00351   xyz[1] = yVal;
00352   xyz[2] = zVal;
00353 }
00354 inline void CubitVector::get_xyz(double &xOut, double &yOut, double &zOut) const
00355 {
00356   xOut = xVal;
00357   yOut = yVal;
00358   zOut = zVal;
00359 }
00360 inline double &CubitVector::r()
00361 { return xVal; }
00362 inline double &CubitVector::theta()
00363 { return yVal; }
00364 inline void CubitVector::x( const double xIn )
00365 { xVal = xIn; }
00366 inline void CubitVector::y( const double yIn )
00367 { yVal = yIn; }
00368 inline void CubitVector::z( const double zIn )
00369 { zVal = zIn; }
00370 inline void CubitVector::r( const double xIn )
00371 { xVal = xIn; }
00372 inline void CubitVector::theta( const double yIn )
00373 { yVal = yIn; }
00374 inline CubitVector& CubitVector::operator+=(const CubitVector &v)
00375 {
00376   xVal += v.xVal;
00377   yVal += v.yVal;
00378   zVal += v.zVal;
00379   return *this;
00380 }
00381 
00382 inline CubitVector& CubitVector::operator-=(const CubitVector &v)
00383 {
00384   xVal -= v.xVal;
00385   yVal -= v.yVal;
00386   zVal -= v.zVal;
00387   return *this;
00388 }
00389 
00390 inline CubitVector& CubitVector::operator*=(const CubitVector &v)
00391 {
00392   double xcross, ycross, zcross;
00393   xcross = yVal * v.zVal - zVal * v.yVal;
00394   ycross = zVal * v.xVal - xVal * v.zVal;
00395   zcross = xVal * v.yVal - yVal * v.xVal;
00396   xVal = xcross;
00397   yVal = ycross;
00398   zVal = zcross;
00399   return *this;
00400 }
00401 
00402 inline CubitVector::CubitVector(const CubitVector& copy_from)
00403     : xVal(copy_from.xVal), yVal(copy_from.yVal), zVal(copy_from.zVal)
00404 {}
00405 
00406 inline CubitVector::CubitVector()
00407     : xVal(0), yVal(0), zVal(0)
00408 {}
00409 
00410 inline CubitVector::CubitVector (const CubitVector& tail,
00411                                  const CubitVector& head)
00412     : xVal(head.xVal - tail.xVal),
00413       yVal(head.yVal - tail.yVal),
00414       zVal(head.zVal - tail.zVal)
00415 {}
00416 
00417 inline CubitVector::CubitVector(const double xIn,
00418                                 const double yIn,
00419                                 const double zIn)
00420     : xVal(xIn), yVal(yIn), zVal(zIn)
00421 {}
00422 
00423 inline CubitVector::CubitVector(const double xyz[3]) 
00424     : xVal(xyz[0]), yVal(xyz[1]), zVal(xyz[2])
00425 {}
00426 
00427 // This sets the vector to be perpendicular to it's current direction.
00428 // NOTE:
00429 //      This is a 2D function.  It only works in the XY Plane.
00430 inline void CubitVector::perpendicular_z()
00431 {
00432   double temp = xVal;
00433   x( yVal );
00434   y( -temp );
00435 }
00436 
00437 inline void CubitVector::set(const double xIn,
00438                              const double yIn,
00439                              const double zIn)
00440 {
00441   xVal = xIn;
00442   yVal = yIn;
00443   zVal = zIn;
00444 }
00445 
00446 inline void CubitVector::set(const double xyz[3])
00447 {
00448   xVal = xyz[0];
00449   yVal = xyz[1];
00450   zVal = xyz[2];
00451 }
00452 
00453 inline void CubitVector::set(const CubitVector& tail,
00454                              const CubitVector& head)
00455 {
00456   xVal = head.xVal - tail.xVal;
00457   yVal = head.yVal - tail.yVal;
00458   zVal = head.zVal - tail.zVal;
00459 }
00460 
00461 inline CubitVector& CubitVector::operator=(const CubitVector &from)  
00462 {
00463   xVal = from.xVal; 
00464   yVal = from.yVal; 
00465   zVal = from.zVal; 
00466   return *this;
00467 }
00468 
00469 inline void CubitVector::set(const CubitVector& to_copy)
00470 {
00471   *this = to_copy;
00472 }
00473 
00474 // Scale all values by scalar.
00475 inline CubitVector& CubitVector::operator*=(const double scalar)
00476 {
00477   xVal *= scalar;
00478   yVal *= scalar;
00479   zVal *= scalar;
00480   return *this;
00481 }
00482 
00483 // Scales all values by 1/scalar
00484 inline CubitVector& CubitVector::operator/=(const double scalar)
00485 {
00486   if(scalar == 0)
00487     throw ("Cannot divide by zero.");
00488   //assert (scalar != 0);
00489   xVal /= scalar;
00490   yVal /= scalar;
00491   zVal /= scalar;
00492   return *this;
00493 }
00494 
00495 inline CubitVector& CubitVector::operator=(const CubitVectorStruct &from)  
00496 {
00497   xVal = from.xVal; 
00498   yVal = from.yVal; 
00499   zVal = from.zVal; 
00500   return *this;
00501 }
00502 
00503 inline CubitVector::CubitVector(const CubitVectorStruct &from)  
00504 {
00505   xVal = from.xVal; 
00506   yVal = from.yVal; 
00507   zVal = from.zVal; 
00508 }
00509 
00510 // Returns the normalized 'this'.
00511 inline CubitVector operator~(const CubitVector &vec)
00512 {
00513   double mag = sqrt(vec.xVal*vec.xVal +
00514                     vec.yVal*vec.yVal +
00515                     vec.zVal*vec.zVal);
00516   
00517   CubitVector temp = vec;
00518   if (mag != 0.0)
00519   {
00520     temp /= mag;
00521   }
00522   return temp;
00523 }
00524 
00525 // Unary minus.  Negates all values in vector.
00526 inline CubitVector CubitVector::operator-() const
00527 {
00528   return CubitVector(-xVal, -yVal, -zVal);
00529 }
00530 
00531 inline double CubitVector::operator[](int i) const
00532 {
00533   if(i < 0 || i > 2)
00534     throw ("Index Out of Bounds");
00535   //assert(i > -1 && i < 3);
00536   if      (i == 0) return xVal;
00537   else if (i == 1) return yVal;
00538   else             return zVal;
00539 }
00540 
00541 inline CubitVector operator+(const CubitVector &vector1,
00542                       const CubitVector &vector2)
00543 {
00544   double xv = vector1.xVal + vector2.xVal;
00545   double yv = vector1.yVal + vector2.yVal;
00546   double zv = vector1.zVal + vector2.zVal;
00547   return CubitVector(xv,yv,zv);
00548 //  return CubitVector(vector1) += vector2;
00549 }
00550 
00551 inline CubitVector operator-(const CubitVector &vector1,
00552                       const CubitVector &vector2)
00553 {
00554   double xv = vector1.xVal - vector2.xVal;
00555   double yv = vector1.yVal - vector2.yVal;
00556   double zv = vector1.zVal - vector2.zVal;
00557   return CubitVector(xv,yv,zv);
00558 //  return CubitVector(vector1) -= vector2;
00559 }
00560 
00561 // Cross products.
00562 // vector1 cross vector2
00563 inline CubitVector operator*(const CubitVector &vector1,
00564                       const CubitVector &vector2)
00565 {
00566   return CubitVector(vector1) *= vector2;
00567 }
00568 
00569 // Returns a scaled vector.
00570 inline CubitVector operator*(const CubitVector &vector1,
00571                       const double scalar)
00572 {
00573   return CubitVector(vector1) *= scalar;
00574 }
00575 
00576 // Returns a scaled vector
00577 inline CubitVector operator*(const double scalar,
00578                              const CubitVector &vector1)
00579 {
00580   return CubitVector(vector1) *= scalar;
00581 }
00582 
00583 // Returns a vector scaled by 1/scalar
00584 inline CubitVector operator/(const CubitVector &vector1,
00585                              const double scalar)
00586 {
00587   return CubitVector(vector1) /= scalar;
00588 }
00589 
00590 inline int operator==(const CubitVector &v1, const CubitVector &v2)
00591 {
00592   return (v1.xVal == v2.xVal && v1.yVal == v2.yVal && v1.zVal == v2.zVal);
00593 }
00594 
00595 inline int operator!=(const CubitVector &v1, const CubitVector &v2)
00596 {
00597   return (v1.xVal != v2.xVal || v1.yVal != v2.yVal || v1.zVal != v2.zVal);
00598 }
00599 
00600 inline double CubitVector::length_squared() const
00601 {
00602   return( xVal*xVal + yVal*yVal + zVal*zVal );
00603 }
00604 
00605 inline double CubitVector::length() const
00606 {
00607   return( sqrt(xVal*xVal + yVal*yVal + zVal*zVal) );
00608 }
00609 
00610 inline double CubitVector::normalize()
00611 {
00612   double mag = length();
00613   if (mag != 0)
00614   {
00615     xVal = xVal / mag;
00616     yVal = yVal / mag;
00617     zVal = zVal / mag;
00618   }
00619   return mag;
00620 }
00621 
00622 
00623 // Dot Product.
00624 inline double operator%(const CubitVector &vector1,
00625                         const CubitVector &vector2)
00626 {
00627   return( vector1.xVal * vector2.xVal +
00628           vector1.yVal * vector2.yVal +
00629           vector1.zVal * vector2.zVal );
00630 }
00631 
00632 // Interpolate between two vectors.
00633 // Returns (1-param)*v1 + param*v2
00634 inline CubitVector interpolate(const double param, const CubitVector &v1,
00635                                const CubitVector &v2)
00636 {
00637   CubitVector temp = (1.0 - param) * v1;
00638   temp += param * v2;
00639   return temp;
00640 }
00641 
00642 inline CubitVector CubitVector::normal( const CubitVector &p0,
00643                                         const CubitVector &p1,
00644                                         const CubitVector &p2 )
00645 {
00646    CubitVector edge0( p0, p1 );
00647    CubitVector edge1( p0, p2 );
00648    
00649    return edge0 * edge1; // not normalized.
00650 }
00651 
00652 #endif
00653 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines