MOAB: Mesh Oriented datABase  (version 5.4.1)
VerdictVector.cpp
Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Module:    $RCSfile: VerdictVector.cpp,v $
00004 
00005   Copyright (c) 2006 Sandia Corporation.
00006   All rights reserved.
00007   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
00008 
00009      This software is distributed WITHOUT ANY WARRANTY; without even
00010      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00011      PURPOSE.  See the above copyright notice for more information.
00012 
00013 =========================================================================*/
00014 
00015 /*
00016  *
00017  * VerdictVector.cpp contains implementation of Vector operations
00018  *
00019  * This file is part of VERDICT
00020  *
00021  */
00022 
00023 #define VERDICT_EXPORTS
00024 
00025 #include "moab/verdict.h"
00026 #include <cmath>
00027 #include "VerdictVector.hpp"
00028 #include <cfloat>
00029 
00030 #if defined( __BORLANDC__ )
00031 #pragma warn - 8004 /* "assigned a value that is never used" */
00032 #endif
00033 
00034 const double TWO_VERDICT_PI = 2.0 * VERDICT_PI;
00035 
00036 VerdictVector& VerdictVector::length( const double new_length )
00037 {
00038     double len = this->length();
00039     xVal *= new_length / len;
00040     yVal *= new_length / len;
00041     zVal *= new_length / len;
00042     return *this;
00043 }
00044 
00045 double VerdictVector::distance_between( const VerdictVector& test_vector )
00046 {
00047     double xv = xVal - test_vector.x();
00048     double yv = yVal - test_vector.y();
00049     double zv = zVal - test_vector.z();
00050 
00051     return ( sqrt( xv * xv + yv * yv + zv * zv ) );
00052 }
00053 
00054 /*
00055 void VerdictVector::print_me()
00056 {
00057   PRINT_INFO("X: %f\n",xVal);
00058   PRINT_INFO("Y: %f\n",yVal);
00059   PRINT_INFO("Z: %f\n",zVal);
00060 
00061 }
00062 */
00063 
00064 double VerdictVector::interior_angle( const VerdictVector& otherVector )
00065 {
00066     double cosAngle = 0., angleRad = 0., len1, len2 = 0.;
00067 
00068     if( ( ( len1 = this->length() ) > 0 ) && ( ( len2 = otherVector.length() ) > 0 ) )
00069         cosAngle = ( *this % otherVector ) / ( len1 * len2 );
00070     else
00071     {
00072         assert( len1 > 0 );
00073         assert( len2 > 0 );
00074     }
00075 
00076     if( ( cosAngle > 1.0 ) && ( cosAngle < 1.0001 ) )
00077     {
00078         cosAngle = 1.0;
00079         angleRad = acos( cosAngle );
00080     }
00081     else if( cosAngle < -1.0 && cosAngle > -1.0001 )
00082     {
00083         cosAngle = -1.0;
00084         angleRad = acos( cosAngle );
00085     }
00086     else if( cosAngle >= -1.0 && cosAngle <= 1.0 )
00087         angleRad = acos( cosAngle );
00088     else
00089     {
00090         assert( cosAngle < 1.0001 && cosAngle > -1.0001 );
00091     }
00092 
00093     return ( ( angleRad * 180. ) / VERDICT_PI );
00094 }
00095 
00096 // Interpolate between two vectors.
00097 // Returns (1-param)*v1 + param*v2
00098 VerdictVector interpolate( const double param, const VerdictVector& v1, const VerdictVector& v2 )
00099 {
00100     VerdictVector temp = ( 1.0 - param ) * v1;
00101     temp += param * v2;
00102     return temp;
00103 }
00104 
00105 void VerdictVector::xy_to_rtheta()
00106 {
00107     // careful about overwriting
00108     double r_     = length();
00109     double theta_ = atan2( y(), x() );
00110     if( theta_ < 0.0 ) theta_ += TWO_VERDICT_PI;
00111 
00112     r( r_ );
00113     theta( theta_ );
00114 }
00115 
00116 void VerdictVector::rtheta_to_xy()
00117 {
00118     // careful about overwriting
00119     double x_ = r() * cos( theta() );
00120     double y_ = r() * sin( theta() );
00121 
00122     x( x_ );
00123     y( y_ );
00124 }
00125 
00126 void VerdictVector::rotate( double angle, double )
00127 {
00128     xy_to_rtheta();
00129     theta() += angle;
00130     rtheta_to_xy();
00131 }
00132 
00133 void VerdictVector::blow_out( double gamma, double rmin )
00134 {
00135     // if gamma == 1, then
00136     // map on a circle : r'^2 = sqrt( 1 - (1-r)^2 )
00137     // if gamma ==0, then map back to itself
00138     // in between, linearly interpolate
00139     xy_to_rtheta();
00140     //  r() = sqrt( (2. - r()) * r() ) * gamma  + r() * (1-gamma);
00141     assert( gamma > 0.0 );
00142     // the following limits should really be roundoff-based
00143     if( r() > rmin * 1.001 && r() < 1.001 )
00144     {
00145         r() = rmin + pow( r(), gamma ) * ( 1.0 - rmin );
00146     }
00147     rtheta_to_xy();
00148 }
00149 
00150 void VerdictVector::reflect_about_xaxis( double, double )
00151 {
00152     yVal = -yVal;
00153 }
00154 
00155 void VerdictVector::scale_angle( double gamma, double )
00156 {
00157     const double r_factor     = 0.3;
00158     const double theta_factor = 0.6;
00159 
00160     xy_to_rtheta();
00161 
00162     // if neary 2pi, treat as zero
00163     // some near zero stuff strays due to roundoff
00164     if( theta() > TWO_VERDICT_PI - 0.02 ) theta() = 0;
00165     // the above screws up on big sheets - need to overhaul at the sheet level
00166 
00167     if( gamma < 1 )
00168     {
00169         // squeeze together points of short radius so that
00170         // long chords won't cross them
00171         theta() += ( VERDICT_PI - theta() ) * ( 1 - gamma ) * theta_factor * ( 1 - r() );
00172 
00173         // push away from center of circle, again so long chords won't cross
00174         r( ( r_factor + r() ) / ( 1 + r_factor ) );
00175 
00176         // scale angle by gamma
00177         theta() *= gamma;
00178     }
00179     else
00180     {
00181         // scale angle by gamma, making sure points nearly 2pi are treated as zero
00182         double new_theta = theta() * gamma;
00183         if( new_theta < 2.5 * VERDICT_PI || r() < 0.2 ) theta( new_theta );
00184     }
00185     rtheta_to_xy();
00186 }
00187 
00188 double VerdictVector::vector_angle_quick( const VerdictVector& vec1, const VerdictVector& vec2 )
00189 {
00190     //- compute the angle between two vectors in the plane defined by this vector
00191     // build yAxis and xAxis such that xAxis is the projection of
00192     // vec1 onto the normal plane of this vector
00193 
00194     // NOTE: vec1 and vec2 are Vectors from the vertex of the angle along
00195     //       the two sides of the angle.
00196     //       The angle returned is the right-handed angle around this vector
00197     //       from vec1 to vec2.
00198 
00199     // NOTE: vector_angle_quick gives exactly the same answer as vector_angle below
00200     //       providing this vector is normalized.  It does so with two fewer
00201     //       cross-product evaluations and two fewer vector normalizations.
00202     //       This can be a substantial time savings if the function is called
00203     //       a significant number of times (e.g Hexer) ... (jrh 11/28/94)
00204     // NOTE: vector_angle() is much more robust. Do not use vector_angle_quick()
00205     //       unless you are very sure of the safety of your input vectors.
00206 
00207     VerdictVector ry = ( *this ) * vec1;
00208     VerdictVector rx = ry * ( *this );
00209 
00210     double xv = vec2 % rx;
00211     double yv = vec2 % ry;
00212 
00213     double angle;
00214     assert( xv != 0.0 || yv != 0.0 );
00215 
00216     angle = atan2( yv, xv );
00217 
00218     if( angle < 0.0 )
00219     {
00220         angle += TWO_VERDICT_PI;
00221     }
00222     return angle;
00223 }
00224 
00225 VerdictVector vectorRotate( const double angle, const VerdictVector& normalAxis, const VerdictVector& referenceAxis )
00226 {
00227     // A new coordinate system is created with the xy plane corresponding
00228     // to the plane normal to the normal axis, and the x axis corresponding to
00229     // the projection of the reference axis onto the normal plane.  The normal
00230     // plane is the tangent plane at the root point.  A unit vector is
00231     // constructed along the local x axis and then rotated by the given
00232     // ccw angle to form the new point.  The new point, then is a unit
00233     // distance from the global origin in the tangent plane.
00234 
00235     double x, y;
00236 
00237     // project a unit distance from root along reference axis
00238 
00239     VerdictVector yAxis = normalAxis * referenceAxis;
00240     VerdictVector xAxis = yAxis * normalAxis;
00241     yAxis.normalize();
00242     xAxis.normalize();
00243 
00244     x = cos( angle );
00245     y = sin( angle );
00246 
00247     xAxis *= x;
00248     yAxis *= y;
00249     return VerdictVector( xAxis + yAxis );
00250 }
00251 
00252 double VerdictVector::vector_angle( const VerdictVector& vector1, const VerdictVector& vector2 ) const
00253 {
00254     // This routine does not assume that any of the input vectors are of unit
00255     // length. This routine does not normalize the input vectors.
00256     // Special cases:
00257     //     If the normal vector is zero length:
00258     //         If a new one can be computed from vectors 1 & 2:
00259     //             the normal is replaced with the vector cross product
00260     //         else the two vectors are colinear and zero or 2PI is returned.
00261     //     If the normal is colinear with either (or both) vectors
00262     //         a new one is computed with the cross products
00263     //         (and checked again).
00264 
00265     // Check for zero length normal vector
00266     VerdictVector normal = *this;
00267     double normal_lensq  = normal.length_squared();
00268     double len_tol       = 0.0000001;
00269     if( normal_lensq <= len_tol )
00270     {
00271         // null normal - make it the normal to the plane defined by vector1
00272         // and vector2. If still null, the vectors are colinear so check
00273         // for zero or 180 angle.
00274         normal       = vector1 * vector2;
00275         normal_lensq = normal.length_squared();
00276         if( normal_lensq <= len_tol )
00277         {
00278             double cosine = vector1 % vector2;
00279             if( cosine > 0.0 )
00280                 return 0.0;
00281             else
00282                 return VERDICT_PI;
00283         }
00284     }
00285 
00286     // Trap for normal vector colinear to one of the other vectors. If so,
00287     // use a normal defined by the two vectors.
00288     double dot_tol = 0.985;
00289     double dot     = vector1 % normal;
00290     if( dot * dot >= vector1.length_squared() * normal_lensq * dot_tol )
00291     {
00292         normal       = vector1 * vector2;
00293         normal_lensq = normal.length_squared();
00294 
00295         // Still problems if all three vectors were colinear
00296         if( normal_lensq <= len_tol )
00297         {
00298             double cosine = vector1 % vector2;
00299             if( cosine >= 0.0 )
00300                 return 0.0;
00301             else
00302                 return VERDICT_PI;
00303         }
00304     }
00305     else
00306     {
00307         // The normal and vector1 are not colinear, now check for vector2
00308         dot = vector2 % normal;
00309         if( dot * dot >= vector2.length_squared() * normal_lensq * dot_tol )
00310         {
00311             normal = vector1 * vector2;
00312         }
00313     }
00314 
00315     // Assume a plane such that the normal vector is the plane's normal.
00316     // Create yAxis perpendicular to both the normal and vector1. yAxis is
00317     // now in the plane. Create xAxis as the perpendicular to both yAxis and
00318     // the normal. xAxis is in the plane and is the projection of vector1
00319     // into the plane.
00320 
00321     normal.normalize();
00322     VerdictVector yAxis = normal;
00323     yAxis *= vector1;
00324     double yv = vector2 % yAxis;
00325     //  yAxis memory slot will now be used for xAxis
00326     yAxis *= normal;
00327     double xv = vector2 % yAxis;
00328 
00329     //  assert(x != 0.0 || y != 0.0);
00330     if( xv == 0.0 && yv == 0.0 )
00331     {
00332         return 0.0;
00333     }
00334     double angle = atan2( yv, xv );
00335 
00336     if( angle < 0.0 )
00337     {
00338         angle += TWO_VERDICT_PI;
00339     }
00340     return angle;
00341 }
00342 
00343 bool VerdictVector::within_tolerance( const VerdictVector& vectorPtr2, double tolerance ) const
00344 {
00345     if( ( fabs( this->x() - vectorPtr2.x() ) < tolerance ) && ( fabs( this->y() - vectorPtr2.y() ) < tolerance ) &&
00346         ( fabs( this->z() - vectorPtr2.z() ) < tolerance ) )
00347     {
00348         return true;
00349     }
00350 
00351     return false;
00352 }
00353 
00354 void VerdictVector::orthogonal_vectors( VerdictVector& vector2, VerdictVector& vector3 )
00355 {
00356     double xv[3];
00357     unsigned short i    = 0;
00358     unsigned short imin = 0;
00359     double rmin         = 1.0E20;
00360     unsigned short iperm1[3];
00361     unsigned short iperm2[3];
00362     unsigned short cont_flag = 1;
00363     double vec1[3], vec2[3];
00364     double rmag;
00365 
00366     // Copy the input vector and normalize it
00367     VerdictVector vector1 = *this;
00368     vector1.normalize();
00369 
00370     // Initialize perm flags
00371     iperm1[0] = 1;
00372     iperm1[1] = 2;
00373     iperm1[2] = 0;
00374     iperm2[0] = 2;
00375     iperm2[1] = 0;
00376     iperm2[2] = 1;
00377 
00378     // Get into the array format we can work with
00379     vector1.get_xyz( vec1 );
00380 
00381     while( i < 3 && cont_flag )
00382     {
00383         if( fabs( vec1[i] ) < 1e-6 )
00384         {
00385             vec2[i]         = 1.0;
00386             vec2[iperm1[i]] = 0.0;
00387             vec2[iperm2[i]] = 0.0;
00388             cont_flag       = 0;
00389         }
00390 
00391         if( fabs( vec1[i] ) < rmin )
00392         {
00393             imin = i;
00394             rmin = fabs( vec1[i] );
00395         }
00396         ++i;
00397     }
00398 
00399     if( cont_flag )
00400     {
00401         xv[imin]         = 1.0;
00402         xv[iperm1[imin]] = 0.0;
00403         xv[iperm2[imin]] = 0.0;
00404 
00405         // Determine cross product
00406         vec2[0] = vec1[1] * xv[2] - vec1[2] * xv[1];
00407         vec2[1] = vec1[2] * xv[0] - vec1[0] * xv[2];
00408         vec2[2] = vec1[0] * xv[1] - vec1[1] * xv[0];
00409 
00410         // Unitize
00411         rmag = sqrt( vec2[0] * vec2[0] + vec2[1] * vec2[1] + vec2[2] * vec2[2] );
00412         vec2[0] /= rmag;
00413         vec2[1] /= rmag;
00414         vec2[2] /= rmag;
00415     }
00416 
00417     // Copy 1st orthogonal vector into VerdictVector vector2
00418     vector2.set( vec2 );
00419 
00420     // Cross vectors to determine last orthogonal vector
00421     vector3 = vector1 * vector2;
00422 }
00423 
00424 //- Find next point from this point using a direction and distance
00425 void VerdictVector::next_point( const VerdictVector& direction, double distance, VerdictVector& out_point )
00426 {
00427     VerdictVector my_direction = direction;
00428     my_direction.normalize();
00429 
00430     // Determine next point in space
00431     out_point.x( xVal + ( distance * my_direction.x() ) );
00432     out_point.y( yVal + ( distance * my_direction.y() ) );
00433     out_point.z( zVal + ( distance * my_direction.z() ) );
00434 
00435     return;
00436 }
00437 
00438 VerdictVector::VerdictVector( const double xyz[3] ) : xVal( xyz[0] ), yVal( xyz[1] ), zVal( xyz[2] ) {}
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines