Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
VerdictVector.cpp
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Module: $RCSfile: VerdictVector.cpp,v $
4 
5  Copyright (c) 2006 Sandia Corporation.
6  All rights reserved.
7  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
8 
9  This software is distributed WITHOUT ANY WARRANTY; without even
10  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11  PURPOSE. See the above copyright notice for more information.
12 
13 =========================================================================*/
14 
15 /*
16  *
17  * VerdictVector.cpp contains implementation of Vector operations
18  *
19  * This file is part of VERDICT
20  *
21  */
22 
23 #define VERDICT_EXPORTS
24 
25 #include "moab/verdict.h"
26 #include <cmath>
27 #include "VerdictVector.hpp"
28 #include <cfloat>
29 
30 #if defined( __BORLANDC__ )
31 #pragma warn - 8004 /* "assigned a value that is never used" */
32 #endif
33 
34 const double TWO_VERDICT_PI = 2.0 * VERDICT_PI;
35 
36 VerdictVector& VerdictVector::length( const double new_length )
37 {
38  double len = this->length();
39  xVal *= new_length / len;
40  yVal *= new_length / len;
41  zVal *= new_length / len;
42  return *this;
43 }
44 
45 double VerdictVector::distance_between( const VerdictVector& test_vector )
46 {
47  double xv = xVal - test_vector.x();
48  double yv = yVal - test_vector.y();
49  double zv = zVal - test_vector.z();
50 
51  return ( sqrt( xv * xv + yv * yv + zv * zv ) );
52 }
53 
54 /*
55 void VerdictVector::print_me()
56 {
57  PRINT_INFO("X: %f\n",xVal);
58  PRINT_INFO("Y: %f\n",yVal);
59  PRINT_INFO("Z: %f\n",zVal);
60 
61 }
62 */
63 
64 double VerdictVector::interior_angle( const VerdictVector& otherVector )
65 {
66  double cosAngle = 0., angleRad = 0., len1, len2 = 0.;
67 
68  if( ( ( len1 = this->length() ) > 0 ) && ( ( len2 = otherVector.length() ) > 0 ) )
69  cosAngle = ( *this % otherVector ) / ( len1 * len2 );
70  else
71  {
72  assert( len1 > 0 );
73  assert( len2 > 0 );
74  }
75 
76  if( ( cosAngle > 1.0 ) && ( cosAngle < 1.0001 ) )
77  {
78  cosAngle = 1.0;
79  angleRad = acos( cosAngle );
80  }
81  else if( cosAngle < -1.0 && cosAngle > -1.0001 )
82  {
83  cosAngle = -1.0;
84  angleRad = acos( cosAngle );
85  }
86  else if( cosAngle >= -1.0 && cosAngle <= 1.0 )
87  angleRad = acos( cosAngle );
88  else
89  {
90  assert( cosAngle < 1.0001 && cosAngle > -1.0001 );
91  }
92 
93  return ( ( angleRad * 180. ) / VERDICT_PI );
94 }
95 
96 // Interpolate between two vectors.
97 // Returns (1-param)*v1 + param*v2
98 VerdictVector interpolate( const double param, const VerdictVector& v1, const VerdictVector& v2 )
99 {
100  VerdictVector temp = ( 1.0 - param ) * v1;
101  temp += param * v2;
102  return temp;
103 }
104 
106 {
107  // careful about overwriting
108  double r_ = length();
109  double theta_ = atan2( y(), x() );
110  if( theta_ < 0.0 ) theta_ += TWO_VERDICT_PI;
111 
112  r( r_ );
113  theta( theta_ );
114 }
115 
117 {
118  // careful about overwriting
119  double x_ = r() * cos( theta() );
120  double y_ = r() * sin( theta() );
121 
122  x( x_ );
123  y( y_ );
124 }
125 
126 void VerdictVector::rotate( double angle, double )
127 {
128  xy_to_rtheta();
129  theta() += angle;
130  rtheta_to_xy();
131 }
132 
133 void VerdictVector::blow_out( double gamma, double rmin )
134 {
135  // if gamma == 1, then
136  // map on a circle : r'^2 = sqrt( 1 - (1-r)^2 )
137  // if gamma ==0, then map back to itself
138  // in between, linearly interpolate
139  xy_to_rtheta();
140  // r() = sqrt( (2. - r()) * r() ) * gamma + r() * (1-gamma);
141  assert( gamma > 0.0 );
142  // the following limits should really be roundoff-based
143  if( r() > rmin * 1.001 && r() < 1.001 )
144  {
145  r() = rmin + pow( r(), gamma ) * ( 1.0 - rmin );
146  }
147  rtheta_to_xy();
148 }
149 
151 {
152  yVal = -yVal;
153 }
154 
155 void VerdictVector::scale_angle( double gamma, double )
156 {
157  const double r_factor = 0.3;
158  const double theta_factor = 0.6;
159 
160  xy_to_rtheta();
161 
162  // if neary 2pi, treat as zero
163  // some near zero stuff strays due to roundoff
164  if( theta() > TWO_VERDICT_PI - 0.02 ) theta() = 0;
165  // the above screws up on big sheets - need to overhaul at the sheet level
166 
167  if( gamma < 1 )
168  {
169  // squeeze together points of short radius so that
170  // long chords won't cross them
171  theta() += ( VERDICT_PI - theta() ) * ( 1 - gamma ) * theta_factor * ( 1 - r() );
172 
173  // push away from center of circle, again so long chords won't cross
174  r( ( r_factor + r() ) / ( 1 + r_factor ) );
175 
176  // scale angle by gamma
177  theta() *= gamma;
178  }
179  else
180  {
181  // scale angle by gamma, making sure points nearly 2pi are treated as zero
182  double new_theta = theta() * gamma;
183  if( new_theta < 2.5 * VERDICT_PI || r() < 0.2 ) theta( new_theta );
184  }
185  rtheta_to_xy();
186 }
187 
189 {
190  //- compute the angle between two vectors in the plane defined by this vector
191  // build yAxis and xAxis such that xAxis is the projection of
192  // vec1 onto the normal plane of this vector
193 
194  // NOTE: vec1 and vec2 are Vectors from the vertex of the angle along
195  // the two sides of the angle.
196  // The angle returned is the right-handed angle around this vector
197  // from vec1 to vec2.
198 
199  // NOTE: vector_angle_quick gives exactly the same answer as vector_angle below
200  // providing this vector is normalized. It does so with two fewer
201  // cross-product evaluations and two fewer vector normalizations.
202  // This can be a substantial time savings if the function is called
203  // a significant number of times (e.g Hexer) ... (jrh 11/28/94)
204  // NOTE: vector_angle() is much more robust. Do not use vector_angle_quick()
205  // unless you are very sure of the safety of your input vectors.
206 
207  VerdictVector ry = ( *this ) * vec1;
208  VerdictVector rx = ry * ( *this );
209 
210  double xv = vec2 % rx;
211  double yv = vec2 % ry;
212 
213  double angle;
214  assert( xv != 0.0 || yv != 0.0 );
215 
216  angle = atan2( yv, xv );
217 
218  if( angle < 0.0 )
219  {
221  }
222  return angle;
223 }
224 
225 VerdictVector vectorRotate( const double angle, const VerdictVector& normalAxis, const VerdictVector& referenceAxis )
226 {
227  // A new coordinate system is created with the xy plane corresponding
228  // to the plane normal to the normal axis, and the x axis corresponding to
229  // the projection of the reference axis onto the normal plane. The normal
230  // plane is the tangent plane at the root point. A unit vector is
231  // constructed along the local x axis and then rotated by the given
232  // ccw angle to form the new point. The new point, then is a unit
233  // distance from the global origin in the tangent plane.
234 
235  double x, y;
236 
237  // project a unit distance from root along reference axis
238 
239  VerdictVector yAxis = normalAxis * referenceAxis;
240  VerdictVector xAxis = yAxis * normalAxis;
241  yAxis.normalize();
242  xAxis.normalize();
243 
244  x = cos( angle );
245  y = sin( angle );
246 
247  xAxis *= x;
248  yAxis *= y;
249  return VerdictVector( xAxis + yAxis );
250 }
251 
252 double VerdictVector::vector_angle( const VerdictVector& vector1, const VerdictVector& vector2 ) const
253 {
254  // This routine does not assume that any of the input vectors are of unit
255  // length. This routine does not normalize the input vectors.
256  // Special cases:
257  // If the normal vector is zero length:
258  // If a new one can be computed from vectors 1 & 2:
259  // the normal is replaced with the vector cross product
260  // else the two vectors are colinear and zero or 2PI is returned.
261  // If the normal is colinear with either (or both) vectors
262  // a new one is computed with the cross products
263  // (and checked again).
264 
265  // Check for zero length normal vector
266  VerdictVector normal = *this;
267  double normal_lensq = normal.length_squared();
268  double len_tol = 0.0000001;
269  if( normal_lensq <= len_tol )
270  {
271  // null normal - make it the normal to the plane defined by vector1
272  // and vector2. If still null, the vectors are colinear so check
273  // for zero or 180 angle.
274  normal = vector1 * vector2;
275  normal_lensq = normal.length_squared();
276  if( normal_lensq <= len_tol )
277  {
278  double cosine = vector1 % vector2;
279  if( cosine > 0.0 )
280  return 0.0;
281  else
282  return VERDICT_PI;
283  }
284  }
285 
286  // Trap for normal vector colinear to one of the other vectors. If so,
287  // use a normal defined by the two vectors.
288  double dot_tol = 0.985;
289  double dot = vector1 % normal;
290  if( dot * dot >= vector1.length_squared() * normal_lensq * dot_tol )
291  {
292  normal = vector1 * vector2;
293  normal_lensq = normal.length_squared();
294 
295  // Still problems if all three vectors were colinear
296  if( normal_lensq <= len_tol )
297  {
298  double cosine = vector1 % vector2;
299  if( cosine >= 0.0 )
300  return 0.0;
301  else
302  return VERDICT_PI;
303  }
304  }
305  else
306  {
307  // The normal and vector1 are not colinear, now check for vector2
308  dot = vector2 % normal;
309  if( dot * dot >= vector2.length_squared() * normal_lensq * dot_tol )
310  {
311  normal = vector1 * vector2;
312  }
313  }
314 
315  // Assume a plane such that the normal vector is the plane's normal.
316  // Create yAxis perpendicular to both the normal and vector1. yAxis is
317  // now in the plane. Create xAxis as the perpendicular to both yAxis and
318  // the normal. xAxis is in the plane and is the projection of vector1
319  // into the plane.
320 
321  normal.normalize();
322  VerdictVector yAxis = normal;
323  yAxis *= vector1;
324  double yv = vector2 % yAxis;
325  // yAxis memory slot will now be used for xAxis
326  yAxis *= normal;
327  double xv = vector2 % yAxis;
328 
329  // assert(x != 0.0 || y != 0.0);
330  if( xv == 0.0 && yv == 0.0 )
331  {
332  return 0.0;
333  }
334  double angle = atan2( yv, xv );
335 
336  if( angle < 0.0 )
337  {
339  }
340  return angle;
341 }
342 
343 bool VerdictVector::within_tolerance( const VerdictVector& vectorPtr2, double tolerance ) const
344 {
345  if( ( fabs( this->x() - vectorPtr2.x() ) < tolerance ) && ( fabs( this->y() - vectorPtr2.y() ) < tolerance ) &&
346  ( fabs( this->z() - vectorPtr2.z() ) < tolerance ) )
347  {
348  return true;
349  }
350 
351  return false;
352 }
353 
355 {
356  double xv[3];
357  unsigned short i = 0;
358  unsigned short imin = 0;
359  double rmin = 1.0E20;
360  unsigned short iperm1[3];
361  unsigned short iperm2[3];
362  unsigned short cont_flag = 1;
363  double vec1[3], vec2[3];
364  double rmag;
365 
366  // Copy the input vector and normalize it
367  VerdictVector vector1 = *this;
368  vector1.normalize();
369 
370  // Initialize perm flags
371  iperm1[0] = 1;
372  iperm1[1] = 2;
373  iperm1[2] = 0;
374  iperm2[0] = 2;
375  iperm2[1] = 0;
376  iperm2[2] = 1;
377 
378  // Get into the array format we can work with
379  vector1.get_xyz( vec1 );
380 
381  while( i < 3 && cont_flag )
382  {
383  if( fabs( vec1[i] ) < 1e-6 )
384  {
385  vec2[i] = 1.0;
386  vec2[iperm1[i]] = 0.0;
387  vec2[iperm2[i]] = 0.0;
388  cont_flag = 0;
389  }
390 
391  if( fabs( vec1[i] ) < rmin )
392  {
393  imin = i;
394  rmin = fabs( vec1[i] );
395  }
396  ++i;
397  }
398 
399  if( cont_flag )
400  {
401  xv[imin] = 1.0;
402  xv[iperm1[imin]] = 0.0;
403  xv[iperm2[imin]] = 0.0;
404 
405  // Determine cross product
406  vec2[0] = vec1[1] * xv[2] - vec1[2] * xv[1];
407  vec2[1] = vec1[2] * xv[0] - vec1[0] * xv[2];
408  vec2[2] = vec1[0] * xv[1] - vec1[1] * xv[0];
409 
410  // Unitize
411  rmag = sqrt( vec2[0] * vec2[0] + vec2[1] * vec2[1] + vec2[2] * vec2[2] );
412  vec2[0] /= rmag;
413  vec2[1] /= rmag;
414  vec2[2] /= rmag;
415  }
416 
417  // Copy 1st orthogonal vector into VerdictVector vector2
418  vector2.set( vec2 );
419 
420  // Cross vectors to determine last orthogonal vector
421  vector3 = vector1 * vector2;
422 }
423 
424 //- Find next point from this point using a direction and distance
425 void VerdictVector::next_point( const VerdictVector& direction, double distance, VerdictVector& out_point )
426 {
427  VerdictVector my_direction = direction;
428  my_direction.normalize();
429 
430  // Determine next point in space
431  out_point.x( xVal + ( distance * my_direction.x() ) );
432  out_point.y( yVal + ( distance * my_direction.y() ) );
433  out_point.z( zVal + ( distance * my_direction.z() ) );
434 
435  return;
436 }
437 
438 VerdictVector::VerdictVector( const double xyz[3] ) : xVal( xyz[0] ), yVal( xyz[1] ), zVal( xyz[2] ) {}