 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
105 void VerdictVector::xy_to_rtheta()
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
116 void VerdictVector::rtheta_to_xy()
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
150 void VerdictVector::reflect_about_xaxis( double, double )
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
188 double VerdictVector::vector_angle_quick( const VerdictVector& vec1, const VerdictVector& vec2 )
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 {
220 angle += TWO_VERDICT_PI;
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 {
338 angle += TWO_VERDICT_PI;
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
354 void VerdictVector::orthogonal_vectors( VerdictVector& vector2, VerdictVector& vector3 )
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] ) {}
1.9.1.