 |
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.hpp,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.hpp contains declarations of vector operations
18 *
19 * This file is part of VERDICT
20 *
21 */
22
23 #ifndef VERDICTVECTOR_HPP
24 #define VERDICTVECTOR_HPP
25
26 #include "moab/verdict.h"
27 #include <cassert>
28 #include <cmath>
29
30 class VerdictVector;
31 typedef void ( VerdictVector::*transform_function )( double gamma, double gamma2 );
32 // a pointer to some function that transforms the point,
33 // taking a double parameter. e.g. blow_out, rotate, and scale_angle
34
35 class VerdictVector
36 {
37 public:
38 //- Heading: Constructors and Destructor
39 VerdictVector(); //- Default constructor.
40
41 VerdictVector( const double x, const double y, const double z );
42 //- Constructor: create vector from three components
43
44 VerdictVector( const double xyz[3] );
45 //- Constructor: create vector from tuple
46
47 VerdictVector( const VerdictVector& tail, const VerdictVector& head );
48 //- Constructor for a VerdictVector starting at tail and pointing
49 //- to head.
50
51 VerdictVector( const VerdictVector& copy_from ); //- Copy Constructor
52
53 //- Heading: Set and Inquire Functions
54 void set( const double xv, const double yv, const double zv );
55 //- Change vector components to {x}, {y}, and {z}
56
57 void set( const double xyz[3] );
58 //- Change vector components to xyz[0], xyz[1], xyz[2]
59
60 void set( const VerdictVector& tail, const VerdictVector& head );
61 //- Change vector to go from tail to head.
62
63 void set( const VerdictVector& to_copy );
64 //- Same as operator=(const VerdictVector&)
65
66 double x() const; //- Return x component of vector
67
68 double y() const; //- Return y component of vector
69
70 double z() const; //- Return z component of vector
71
72 void get_xyz( double& x, double& y, double& z ); //- Get x, y, z components
73 void get_xyz( double xyz[3] ); //- Get xyz tuple
74
75 double& r(); //- Return r component of vector, if (r,theta) format
76
77 double& theta(); //- Return theta component of vector, if (r,theta) format
78
79 void x( const double xv ); //- Set x component of vector
80
81 void y( const double yv ); //- Set y component of vector
82
83 void z( const double zv ); //- Set z component of vector
84
85 void r( const double xv ); //- Set r component of vector, if (r,theta) format
86
87 void theta( const double yv ); //- Set theta component of vector, if (r,theta) format
88
89 void xy_to_rtheta();
90 //- convert from cartesian to polar coordinates, just 2d for now
91 //- theta is in [0,2 PI)
92
93 void rtheta_to_xy();
94 //- convert from polar to cartesian coordinates, just 2d for now
95
96 void scale_angle( double gamma, double );
97 //- tranform_function.
98 //- transform (x,y) to (r,theta) to (r,gamma*theta) to (x',y')
99 //- plus some additional scaling so long chords won't cross short ones
100
101 void blow_out( double gamma, double gamma2 = 0.0 );
102 //- transform_function
103 //- blow points radially away from the origin,
104
105 void rotate( double angle, double );
106 //- transform function.
107 //- transform (x,y) to (r,theta) to (r,theta+angle) to (x',y')
108
109 void reflect_about_xaxis( double dummy, double );
110 //- dummy argument to make it a transform function
111
112 double normalize();
113 //- Normalize (set magnitude equal to 1) vector - return the magnitude
114
115 VerdictVector& length( const double new_length );
116 //- Change length of vector to {new_length}. Can be used to move a
117 //- location a specified distance from the origin in the current
118 //- orientation.
119
120 double length() const;
121 //- Calculate the length of the vector.
122 //- Use {length_squared()} if only comparing lengths, not adding.
123
124 double distance_between( const VerdictVector& test_vector );
125 //- Calculate the distance from the head of one vector
126 // to the head of the test_vector.
127
128 double length_squared() const;
129 //- Calculate the squared length of the vector.
130 //- Faster than {length()} since it eliminates the square root if
131 //- only comparing other lengths.
132
133 double interior_angle( const VerdictVector& otherVector );
134 //- Calculate the interior angle: acos((a%b)/(|a||b|))
135 //- Returns angle in degrees.
136
137 double vector_angle_quick( const VerdictVector& vec1, const VerdictVector& vec2 );
138 //- Calculate the interior angle between the projections of
139 //- {vec1} and {vec2} onto the plane defined by the {this} vector.
140 //- The angle returned is the right-handed angle around the {this}
141 //- vector from {vec1} to {vec2}. Angle is in radians.
142
143 double vector_angle( const VerdictVector& vector1, const VerdictVector& vector2 ) const;
144 //- Compute the angle between the projections of {vector1} and {vector2}
145 //- onto the plane defined by *this. The angle is the
146 //- right-hand angle, in radians, about *this from {vector1} to {vector2}.
147
148 void perpendicular_z();
149 //- Transform this vector to a perpendicular one, leaving
150 //- z-component alone. Rotates clockwise about the z-axis by pi/2.
151
152 void print_me();
153 //- Prints out the coordinates of this vector.
154
155 void orthogonal_vectors( VerdictVector& vector2, VerdictVector& vector3 );
156 //- Finds 2 (arbitrary) vectors that are orthogonal to this one
157
158 void next_point( const VerdictVector& direction, double distance, VerdictVector& out_point );
159 //- Finds the next point in space based on *this* point (starting point),
160 //- a direction and the distance to extend in the direction. The
161 //- direction vector need not be a unit vector. The out_point can be
162 //- "*this" (i.e., modify point in place).
163
164 bool within_tolerance( const VerdictVector& vectorPtr2, double tolerance ) const;
165 //- Compare two vectors to see if they are spatially equal. They
166 //- compare if x, y, and z are all within tolerance.
167
168 //- Heading: Operator Overloads *****************************
169 VerdictVector& operator+=( const VerdictVector& vec );
170 //- Compound Assignment: addition: {this = this + vec}
171
172 VerdictVector& operator-=( const VerdictVector& vec );
173 //- Compound Assignment: subtraction: {this = this - vec}
174
175 VerdictVector& operator*=( const VerdictVector& vec );
176 //- Compound Assignment: cross product: {this = this * vec},
177 //- non-commutative
178
179 VerdictVector& operator*=( const double scalar );
180 //- Compound Assignment: multiplication: {this = this * scalar}
181
182 VerdictVector& operator/=( const double scalar );
183 //- Compound Assignment: division: {this = this / scalar}
184
185 VerdictVector operator-() const;
186 //- unary negation.
187
188 friend VerdictVector operator~( const VerdictVector& vec );
189 //- normalize. Returns a new vector which is a copy of {vec},
190 //- scaled such that {|vec|=1}. Uses overloaded bitwise NOT operator.
191
192 friend VerdictVector operator+( const VerdictVector& v1, const VerdictVector& v2 );
193 //- vector addition
194
195 friend VerdictVector operator-( const VerdictVector& v1, const VerdictVector& v2 );
196 //- vector subtraction
197
198 friend VerdictVector operator*( const VerdictVector& v1, const VerdictVector& v2 );
199 //- vector cross product, non-commutative
200
201 friend VerdictVector operator*( const VerdictVector& v1, const double sclr );
202 //- vector * scalar
203
204 friend VerdictVector operator*( const double sclr, const VerdictVector& v1 );
205 //- scalar * vector
206
207 friend double operator%( const VerdictVector& v1, const VerdictVector& v2 );
208 //- dot product
209
210 friend VerdictVector operator/( const VerdictVector& v1, const double sclr );
211 //- vector / scalar
212
213 friend int operator==( const VerdictVector& v1, const VerdictVector& v2 );
214 //- Equality operator
215
216 friend int operator!=( const VerdictVector& v1, const VerdictVector& v2 );
217 //- Inequality operator
218
219 friend VerdictVector interpolate( const double param, const VerdictVector& v1, const VerdictVector& v2 );
220 //- Interpolate between two vectors. Returns (1-param)*v1 + param*v2
221
222 VerdictVector& operator=( const VerdictVector& from );
223
224 private:
225 double xVal; //- x component of vector.
226 double yVal; //- y component of vector.
227 double zVal; //- z component of vector.
228 };
229
230 VerdictVector vectorRotate( const double angle, const VerdictVector& normalAxis, const VerdictVector& referenceAxis );
231 //- A new coordinate system is created with the xy plane corresponding
232 //- to the plane normal to {normalAxis}, and the x axis corresponding to
233 //- the projection of {referenceAxis} onto the normal plane. The normal
234 //- plane is the tangent plane at the root point. A unit vector is
235 //- constructed along the local x axis and then rotated by the given
236 //- ccw angle to form the new point. The new point, then is a unit
237 //- distance from the global origin in the tangent plane.
238 //- {angle} is in radians.
239
240 inline double VerdictVector::x() const
241 {
242 return xVal;
243 }
244 inline double VerdictVector::y() const
245 {
246 return yVal;
247 }
248 inline double VerdictVector::z() const
249 {
250 return zVal;
251 }
252 inline void VerdictVector::get_xyz( double xyz[3] )
253 {
254 xyz[0] = xVal;
255 xyz[1] = yVal;
256 xyz[2] = zVal;
257 }
258 inline void VerdictVector::get_xyz( double& xv, double& yv, double& zv )
259 {
260 xv = xVal;
261 yv = yVal;
262 zv = zVal;
263 }
264 inline double& VerdictVector::r()
265 {
266 return xVal;
267 }
268 inline double& VerdictVector::theta()
269 {
270 return yVal;
271 }
272 inline void VerdictVector::x( const double xv )
273 {
274 xVal = xv;
275 }
276 inline void VerdictVector::y( const double yv )
277 {
278 yVal = yv;
279 }
280 inline void VerdictVector::z( const double zv )
281 {
282 zVal = zv;
283 }
284 inline void VerdictVector::r( const double xv )
285 {
286 xVal = xv;
287 }
288 inline void VerdictVector::theta( const double yv )
289 {
290 yVal = yv;
291 }
292 inline VerdictVector& VerdictVector::operator+=( const VerdictVector& vector )
293 {
294 xVal += vector.x();
295 yVal += vector.y();
296 zVal += vector.z();
297 return *this;
298 }
299
300 inline VerdictVector& VerdictVector::operator-=( const VerdictVector& vector )
301 {
302 xVal -= vector.x();
303 yVal -= vector.y();
304 zVal -= vector.z();
305 return *this;
306 }
307
308 inline VerdictVector& VerdictVector::operator*=( const VerdictVector& vector )
309 {
310 double xcross, ycross, zcross;
311 xcross = yVal * vector.z() - zVal * vector.y();
312 ycross = zVal * vector.x() - xVal * vector.z();
313 zcross = xVal * vector.y() - yVal * vector.x();
314 xVal = xcross;
315 yVal = ycross;
316 zVal = zcross;
317 return *this;
318 }
319
320 inline VerdictVector::VerdictVector( const VerdictVector& copy_from )
321 : xVal( copy_from.xVal ), yVal( copy_from.yVal ), zVal( copy_from.zVal )
322 {
323 }
324
325 inline VerdictVector::VerdictVector() : xVal( 0 ), yVal( 0 ), zVal( 0 ) {}
326
327 inline VerdictVector::VerdictVector( const VerdictVector& tail, const VerdictVector& head )
328 : xVal( head.xVal - tail.xVal ), yVal( head.yVal - tail.yVal ), zVal( head.zVal - tail.zVal )
329 {
330 }
331
332 inline VerdictVector::VerdictVector( const double xIn, const double yIn, const double zIn )
333 : xVal( xIn ), yVal( yIn ), zVal( zIn )
334 {
335 }
336
337 // This sets the vector to be perpendicular to it's current direction.
338 // NOTE:
339 // This is a 2D function. It only works in the XY Plane.
340 inline void VerdictVector::perpendicular_z()
341 {
342 double temp = x();
343 x( y() );
344 y( -temp );
345 }
346
347 inline void VerdictVector::set( const double xv, const double yv, const double zv )
348 {
349 xVal = xv;
350 yVal = yv;
351 zVal = zv;
352 }
353
354 inline void VerdictVector::set( const double xyz[3] )
355 {
356 xVal = xyz[0];
357 yVal = xyz[1];
358 zVal = xyz[2];
359 }
360
361 inline void VerdictVector::set( const VerdictVector& tail, const VerdictVector& head )
362 {
363 xVal = head.xVal - tail.xVal;
364 yVal = head.yVal - tail.yVal;
365 zVal = head.zVal - tail.zVal;
366 }
367
368 inline VerdictVector& VerdictVector::operator=( const VerdictVector& from )
369 {
370 xVal = from.xVal;
371 yVal = from.yVal;
372 zVal = from.zVal;
373 return *this;
374 }
375
376 inline void VerdictVector::set( const VerdictVector& to_copy )
377 {
378 *this = to_copy;
379 }
380
381 // Scale all values by scalar.
382 inline VerdictVector& VerdictVector::operator*=( const double scalar )
383 {
384 xVal *= scalar;
385 yVal *= scalar;
386 zVal *= scalar;
387 return *this;
388 }
389
390 // Scales all values by 1/scalar
391 inline VerdictVector& VerdictVector::operator/=( const double scalar )
392 {
393 assert( scalar != 0 );
394 xVal /= scalar;
395 yVal /= scalar;
396 zVal /= scalar;
397 return *this;
398 }
399
400 // Returns the normalized 'this'.
401 inline VerdictVector operator~( const VerdictVector& vec )
402 {
403 double mag = sqrt( vec.xVal * vec.xVal + vec.yVal * vec.yVal + vec.zVal * vec.zVal );
404
405 VerdictVector temp = vec;
406 if( mag != 0.0 )
407 {
408 temp /= mag;
409 }
410 return temp;
411 }
412
413 // Unary minus. Negates all values in vector.
414 inline VerdictVector VerdictVector::operator-() const
415 {
416 return VerdictVector( -xVal, -yVal, -zVal );
417 }
418
419 inline VerdictVector operator+( const VerdictVector& vector1, const VerdictVector& vector2 )
420 {
421 double xv = vector1.x() + vector2.x();
422 double yv = vector1.y() + vector2.y();
423 double zv = vector1.z() + vector2.z();
424 return VerdictVector( xv, yv, zv );
425 // return VerdictVector(vector1) += vector2;
426 }
427
428 inline VerdictVector operator-( const VerdictVector& vector1, const VerdictVector& vector2 )
429 {
430 double xv = vector1.x() - vector2.x();
431 double yv = vector1.y() - vector2.y();
432 double zv = vector1.z() - vector2.z();
433 return VerdictVector( xv, yv, zv );
434 // return VerdictVector(vector1) -= vector2;
435 }
436
437 // Cross products.
438 // vector1 cross vector2
439 inline VerdictVector operator*( const VerdictVector& vector1, const VerdictVector& vector2 )
440 {
441 return VerdictVector( vector1 ) *= vector2;
442 }
443
444 // Returns a scaled vector.
445 inline VerdictVector operator*( const VerdictVector& vector1, const double scalar )
446 {
447 return VerdictVector( vector1 ) *= scalar;
448 }
449
450 // Returns a scaled vector
451 inline VerdictVector operator*( const double scalar, const VerdictVector& vector1 )
452 {
453 return VerdictVector( vector1 ) *= scalar;
454 }
455
456 // Returns a vector scaled by 1/scalar
457 inline VerdictVector operator/( const VerdictVector& vector1, const double scalar )
458 {
459 return VerdictVector( vector1 ) /= scalar;
460 }
461
462 inline int operator==( const VerdictVector& v1, const VerdictVector& v2 )
463 {
464 return ( v1.xVal == v2.xVal && v1.yVal == v2.yVal && v1.zVal == v2.zVal );
465 }
466
467 inline int operator!=( const VerdictVector& v1, const VerdictVector& v2 )
468 {
469 return ( v1.xVal != v2.xVal || v1.yVal != v2.yVal || v1.zVal != v2.zVal );
470 }
471
472 inline double VerdictVector::length_squared() const
473 {
474 return ( xVal * xVal + yVal * yVal + zVal * zVal );
475 }
476
477 inline double VerdictVector::length() const
478 {
479 return ( sqrt( xVal * xVal + yVal * yVal + zVal * zVal ) );
480 }
481
482 inline double VerdictVector::normalize()
483 {
484 double mag = length();
485 if( mag != 0 )
486 {
487 xVal = xVal / mag;
488 yVal = yVal / mag;
489 zVal = zVal / mag;
490 }
491 return mag;
492 }
493
494 // Dot Product.
495 inline double operator%( const VerdictVector& vector1, const VerdictVector& vector2 )
496 {
497 return ( vector1.x() * vector2.x() + vector1.y() * vector2.y() + vector1.z() * vector2.z() );
498 }
499
500 #endif
1.9.1.