Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected]
25 : :
26 : : ***************************************************************** */
27 : : #ifndef MESQUITE_VECTOR3D_HPP
28 : : #define MESQUITE_VECTOR3D_HPP
29 : :
30 : : #include "Mesquite.hpp"
31 : :
32 : : #include <iosfwd>
33 : : #include <cassert>
34 : : #include <cstring>
35 : : #include <vector>
36 : :
37 : : /*! \file Vector3D.hpp
38 : : \brief Vector object with exactly 3 dimensions.
39 : :
40 : : This is as fast as a C array.
41 : :
42 : : \author Darryl Melander
43 : : \author Thomas Leurent
44 : : */
45 : : namespace MBMesquite
46 : : {
47 : : class Matrix3D;
48 : : class MsqError;
49 : :
50 : : /*!
51 : : \class Vector3D
52 : : \brief Vector3D is the object that effeciently stores information about
53 : : about three-deminsional vectors. It is also the parent class of
54 : : MsqVertex. */
55 : : class MESQUITE_EXPORT Vector3D
56 : : {
57 : : public:
58 : : // Constructors
59 : : Vector3D();
60 : : Vector3D( double xyz );
61 : : Vector3D( double x, double y, double z );
62 : : Vector3D( const double xyz[3] );
63 : : Vector3D( const Vector3D& to_copy );
64 : :
65 : : // *** virtual destructor *** Do not use for Vector3D, we need to keep
66 : : // the deallocation of those objects very fast
67 : :
68 : : // Functions to get the coordinates
69 : : double x() const;
70 : : double y() const;
71 : : double z() const;
72 : : void get_coordinates( double& x, double& y, double& z ) const;
73 : : void get_coordinates( double xyz[3] ) const;
74 : : const double& operator[]( size_t index ) const; // 0-based
75 : :
76 : : // Functions to set the coordinates.
77 : : void x( const double x );
78 : : void y( const double y );
79 : : void z( const double z );
80 : : void set( const double x, const double y, const double z );
81 : : void set( const double xyz[3] );
82 : : void set( const Vector3D& to_copy );
83 : : // Subscripts on non-consts both get and set coords
84 : : double& operator[]( size_t index ); // 0-based
85 : : Vector3D& operator=( const Vector3D& to_copy );
86 : : Vector3D& operator=( const double& to_copy );
87 : :
88 : : // Functions to modify existing coordinates
89 : : Vector3D operator-() const; //- unary negation.
90 : : Vector3D& operator*=( const double scalar );
91 : : Vector3D& operator/=( const double scalar );
92 : : Vector3D& operator*=( const Vector3D& rhs ); //- cross product
93 : : Vector3D& operator+=( const Vector3D& rhs );
94 : : Vector3D& operator-=( const Vector3D& rhs );
95 : :
96 : : // Binary operators (like a+b).
97 : : friend const Vector3D operator+( const Vector3D& lhs, const Vector3D& rhs );
98 : : friend const Vector3D operator-( const Vector3D& lhs, const Vector3D& rhs );
99 : : friend const Vector3D operator*( const Vector3D& lhs,
100 : : const double scalar ); //!< lhs * scalar
101 : : friend const Vector3D operator*( const double scalar,
102 : : const Vector3D& rhs ); //!< scalar * rhs
103 : : friend const Vector3D operator/( const Vector3D& lhs,
104 : : const double scalar ); //- lhs / scalar
105 : : friend double operator%( const Vector3D& v1,
106 : : const Vector3D& v2 ); //!< dot product
107 : : friend double inner( const Vector3D v1[], const Vector3D v2[],
108 : : int n ); //!< dot product for array
109 : : friend double operator%( const double scalar,
110 : : const Vector3D& v2 ); //!< scalar * sum_i v2[i]
111 : : friend double operator%( const Vector3D& v1,
112 : : const double scalar ); //!< scalar * sum_i v1[i]
113 : : friend const Vector3D operator*( const Vector3D& v1,
114 : : const Vector3D& v2 ); //!< cross product
115 : :
116 : : //! \f$ v = A*x \f$
117 : : friend void eqAx( Vector3D& v, const Matrix3D& A, const Vector3D& x );
118 : : //! \f$ v += A*x \f$
119 : : friend void plusEqAx( Vector3D& v, const Matrix3D& A, const Vector3D& x );
120 : : //! \f$ v += A^T*x \f$
121 : : friend void plusEqTransAx( Vector3D& v, const Matrix3D& A, const Vector3D& x );
122 : : friend void eqTransAx( Vector3D& v, const Matrix3D& A, const Vector3D& x );
123 : :
124 : : // Comparison functions
125 : : friend bool operator==( const Vector3D& lhs, const Vector3D& rhs );
126 : : friend bool operator!=( const Vector3D& lhs, const Vector3D& rhs );
127 : : static double distance_between( const Vector3D& p1, const Vector3D& p2 );
128 : : int within_tolerance_box( const Vector3D& compare_to, double tolerance ) const;
129 : : //- Compare two Vector3Ds to see if they are spatially equal.
130 : : // Return TRUE if difference in x, y, and z are all within tolerance.
131 : : // Essentially checks to see if 'this' lies within a box centered on
132 : : // 'compare_to' with sides of length ('tolerance' * 2).
133 : :
134 : : // Length functions
135 : : inline double length_squared() const;
136 : : inline double length() const;
137 : : friend double length( const Vector3D* v, int n ); //!< L2 norm for an array of Vector3Ds
138 : : friend double Linf( const Vector3D* v, int n ); //!< L inf norm for array of Vector3Ds
139 : :
140 : : inline void set_length( const double new_length );
141 : : inline void normalize();
142 : 190150 : Vector3D operator~() const
143 : : {
144 : 190150 : return *this * ( 1.0 / length() );
145 : : }
146 : :
147 : : // Utility functions. All angle functions work in radians.
148 : : static double interior_angle( const Vector3D& a, const Vector3D& b, MsqError& err );
149 : : //- Interior angle: acos((a%b)/(|a||b|))
150 : : static Vector3D interpolate( const double param, const Vector3D& p1, const Vector3D& p2 );
151 : : //- Interpolate between two points. Returns (1-param)*v1 + param*v2.
152 : :
153 : 66075488 : const double* to_array() const
154 : : {
155 : 66075488 : return mCoords;
156 : : }
157 : :
158 : 5258815 : double* to_array()
159 : : {
160 : 5258815 : return mCoords;
161 : : }
162 : :
163 : : protected:
164 : : double mCoords[3];
165 : : };
166 : :
167 : : // Constructors
168 : 102806541 : inline Vector3D::Vector3D()
169 : : {
170 : 102806541 : mCoords[0] = 0;
171 : 102806541 : mCoords[1] = 0;
172 : 102806541 : mCoords[2] = 0;
173 : 102806541 : }
174 : 231271 : inline Vector3D::Vector3D( double p_x )
175 : : {
176 : 231271 : mCoords[0] = p_x;
177 : 231271 : mCoords[1] = p_x;
178 : 231271 : mCoords[2] = p_x;
179 : 231271 : }
180 : 132995990 : inline Vector3D::Vector3D( double p_x, double p_y, double p_z )
181 : : {
182 : 132995990 : mCoords[0] = p_x;
183 : 132995990 : mCoords[1] = p_y;
184 : 132995990 : mCoords[2] = p_z;
185 : 132995990 : }
186 : 30151119 : inline Vector3D::Vector3D( const double xyz[3] )
187 : : {
188 : 30151119 : std::memcpy( mCoords, xyz, 3 * sizeof( double ) );
189 : 30151119 : }
190 : 83461808 : inline Vector3D::Vector3D( const Vector3D& to_copy )
191 : : {
192 : 83461808 : std::memcpy( mCoords, to_copy.mCoords, 3 * sizeof( double ) );
193 : 83461808 : }
194 : :
195 : : // Functions to get coordinates
196 : 148658489 : inline double Vector3D::x() const
197 : : {
198 : 148658489 : return mCoords[0];
199 : : }
200 : 148381895 : inline double Vector3D::y() const
201 : : {
202 : 148381895 : return mCoords[1];
203 : : }
204 : 148381895 : inline double Vector3D::z() const
205 : : {
206 : 148381895 : return mCoords[2];
207 : : }
208 : 3877976 : inline void Vector3D::get_coordinates( double& p_x, double& p_y, double& p_z ) const
209 : : {
210 : 3877976 : p_x = mCoords[0];
211 : 3877976 : p_y = mCoords[1];
212 : 3877976 : p_z = mCoords[2];
213 : 3877976 : }
214 : 11311 : inline void Vector3D::get_coordinates( double xyz[3] ) const
215 : : {
216 : 11311 : std::memcpy( xyz, mCoords, 3 * sizeof( double ) );
217 : 11311 : }
218 : 69919362 : inline const double& Vector3D::operator[]( size_t index ) const
219 : : {
220 : 69919362 : return mCoords[index];
221 : : }
222 : :
223 : : // Functions to set coordinates
224 : 0 : inline void Vector3D::x( const double p_x )
225 : : {
226 : 0 : mCoords[0] = p_x;
227 : 0 : }
228 : 0 : inline void Vector3D::y( const double p_y )
229 : : {
230 : 0 : mCoords[1] = p_y;
231 : 0 : }
232 : 0 : inline void Vector3D::z( const double p_z )
233 : : {
234 : 0 : mCoords[2] = p_z;
235 : 0 : }
236 : 7688 : inline void Vector3D::set( const double p_x, const double p_y, const double p_z )
237 : : {
238 : 7688 : mCoords[0] = p_x;
239 : 7688 : mCoords[1] = p_y;
240 : 7688 : mCoords[2] = p_z;
241 : 7688 : }
242 : 19301556 : inline void Vector3D::set( const double xyz[3] )
243 : : {
244 : 19301556 : std::memcpy( mCoords, xyz, 3 * sizeof( double ) );
245 : 19301556 : }
246 : 111162 : inline void Vector3D::set( const Vector3D& to_copy )
247 : : {
248 : 111162 : std::memcpy( mCoords, to_copy.mCoords, 3 * sizeof( double ) );
249 : 111162 : }
250 : 50692204 : inline double& Vector3D::operator[]( size_t index )
251 : : {
252 : 50692204 : return mCoords[index];
253 : : }
254 : :
255 : 148746396 : inline Vector3D& Vector3D::operator=( const Vector3D& to_copy )
256 : : {
257 : 148746396 : mCoords[0] = to_copy.mCoords[0];
258 : 148746396 : mCoords[1] = to_copy.mCoords[1];
259 : 148746396 : mCoords[2] = to_copy.mCoords[2];
260 : : // memcpy(mCoords, to_copy.mCoords, 3*sizeof(double));
261 : 148746396 : return *this;
262 : : }
263 : :
264 : 613677 : inline Vector3D& Vector3D::operator=( const double& to_copy )
265 : : {
266 : 613677 : mCoords[0] = to_copy;
267 : 613677 : mCoords[1] = to_copy;
268 : 613677 : mCoords[2] = to_copy;
269 : 613677 : return *this;
270 : : }
271 : :
272 : : // Functions that modify existing coordinates
273 : 385679 : inline Vector3D Vector3D::operator-() const
274 : : {
275 : 385679 : return Vector3D( -mCoords[0], -mCoords[1], -mCoords[2] );
276 : : }
277 : 25031753 : inline Vector3D& Vector3D::operator*=( const double scalar )
278 : : {
279 : 25031753 : mCoords[0] *= scalar;
280 : 25031753 : mCoords[1] *= scalar;
281 : 25031753 : mCoords[2] *= scalar;
282 : 25031753 : return *this;
283 : : }
284 : : //! divides each Vector3D entry by the given scalar.
285 : 480713 : inline Vector3D& Vector3D::operator/=( const double scalar )
286 : : {
287 : 480713 : mCoords[0] /= scalar;
288 : 480713 : mCoords[1] /= scalar;
289 : 480713 : mCoords[2] /= scalar;
290 : 480713 : return *this;
291 : : }
292 : : inline Vector3D& Vector3D::operator*=( const Vector3D& rhs )
293 : : {
294 : : double new_coords[3] = { mCoords[1] * rhs.mCoords[2] - mCoords[2] * rhs.mCoords[1],
295 : : mCoords[2] * rhs.mCoords[0] - mCoords[0] * rhs.mCoords[2],
296 : : mCoords[0] * rhs.mCoords[1] - mCoords[1] * rhs.mCoords[0] };
297 : : std::memcpy( mCoords, new_coords, 3 * sizeof( double ) );
298 : : return *this;
299 : : }
300 : 15878468 : inline Vector3D& Vector3D::operator+=( const Vector3D& rhs )
301 : : {
302 : 15878468 : mCoords[0] += rhs.mCoords[0];
303 : 15878468 : mCoords[1] += rhs.mCoords[1];
304 : 15878468 : mCoords[2] += rhs.mCoords[2];
305 : 15878468 : return *this;
306 : : }
307 : 51336415 : inline Vector3D& Vector3D::operator-=( const Vector3D& rhs )
308 : : {
309 : 51336415 : mCoords[0] -= rhs.mCoords[0];
310 : 51336415 : mCoords[1] -= rhs.mCoords[1];
311 : 51336415 : mCoords[2] -= rhs.mCoords[2];
312 : 51336415 : return *this;
313 : : }
314 : :
315 : : // Binary operators
316 : 8410423 : inline const Vector3D operator+( const Vector3D& lhs, const Vector3D& rhs )
317 : : {
318 : 8410423 : return Vector3D( lhs.x() + rhs.x(), lhs.y() + rhs.y(), lhs.z() + rhs.z() );
319 : : }
320 : 53230915 : inline const Vector3D operator-( const Vector3D& lhs, const Vector3D& rhs )
321 : : {
322 : 53230915 : return Vector3D( lhs.x() - rhs.x(), lhs.y() - rhs.y(), lhs.z() - rhs.z() );
323 : : }
324 : 1840998 : inline const Vector3D operator*( const Vector3D& lhs, const double scalar )
325 : : {
326 : 1840998 : return Vector3D( lhs.x() * scalar, lhs.y() * scalar, lhs.z() * scalar );
327 : : }
328 : 15951040 : inline const Vector3D operator*( const double scalar, const Vector3D& rhs )
329 : : {
330 : 15951040 : return Vector3D( rhs.x() * scalar, rhs.y() * scalar, rhs.z() * scalar );
331 : : }
332 : 7307181 : inline const Vector3D operator/( const Vector3D& lhs, const double scalar )
333 : : {
334 [ - + ]: 7307181 : assert( scalar != 0 );
335 : 7307181 : return Vector3D( lhs.x() / scalar, lhs.y() / scalar, lhs.z() / scalar );
336 : : }
337 : 78517273 : inline double operator%( const Vector3D& lhs,
338 : : const Vector3D& rhs ) // Dot Product
339 : : {
340 : 78517273 : return ( lhs.mCoords[0] * rhs.mCoords[0] + lhs.mCoords[1] * rhs.mCoords[1] + lhs.mCoords[2] * rhs.mCoords[2] );
341 : : }
342 : :
343 : : /*! Dot product for arrays of Vector3Ds. see also operator% .*/
344 : 6517 : inline double inner( const Vector3D lhs[], const Vector3D rhs[], int n )
345 : : {
346 : : int i;
347 : 6517 : double dot = 0;
348 [ + + ]: 4743012 : for( i = 0; i < n; ++i )
349 : 9472990 : dot += lhs[i].mCoords[0] * rhs[i].mCoords[0] + lhs[i].mCoords[1] * rhs[i].mCoords[1] +
350 : 9472990 : lhs[i].mCoords[2] * rhs[i].mCoords[2];
351 : 6517 : return dot;
352 : : }
353 : : /*! Dot product for arrays of Vector3Ds. see also operator% .*/
354 : 319 : inline double inner( const std::vector< Vector3D >& lhs, const std::vector< Vector3D >& rhs )
355 : : {
356 : 319 : double dot = 0;
357 [ - + ]: 319 : assert( lhs.size() == rhs.size() );
358 [ + + ]: 34136 : for( size_t i = 0; i < lhs.size(); ++i )
359 : 33817 : dot = lhs[i] % rhs[i];
360 : 319 : return dot;
361 : : }
362 : :
363 : : inline double operator%( const double scalar,
364 : : const Vector3D& rhs ) // Dot Product
365 : : {
366 : : return ( scalar * ( rhs.mCoords[0] + rhs.mCoords[1] + rhs.mCoords[2] ) );
367 : : }
368 : : inline double operator%( const Vector3D& lhs,
369 : : const double scalar ) // Dot Product
370 : : {
371 : : return ( scalar * ( lhs.mCoords[0] + lhs.mCoords[1] + lhs.mCoords[2] ) );
372 : : }
373 : 25327592 : inline const Vector3D operator*( const Vector3D& lhs,
374 : : const Vector3D& rhs ) // Cross Product
375 : : {
376 : 25327592 : return Vector3D( lhs.mCoords[1] * rhs.mCoords[2] - lhs.mCoords[2] * rhs.mCoords[1],
377 : 25327592 : lhs.mCoords[2] * rhs.mCoords[0] - lhs.mCoords[0] * rhs.mCoords[2],
378 : 25327592 : lhs.mCoords[0] * rhs.mCoords[1] - lhs.mCoords[1] * rhs.mCoords[0] );
379 : : }
380 : :
381 : : // output operator
382 : : MESQUITE_EXPORT std::ostream& operator<<( std::ostream& s, const MBMesquite::Vector3D& v );
383 : :
384 : 223 : inline double Vector3D::distance_between( const Vector3D& p1, const Vector3D& p2 )
385 : : {
386 [ + - ]: 223 : Vector3D v = p2 - p1;
387 [ + - ]: 223 : return v.length();
388 : : }
389 : 39 : inline int Vector3D::within_tolerance_box( const Vector3D& compare_to, double tolerance ) const
390 : : {
391 [ + - ]: 39 : return ( ( std::fabs( this->mCoords[0] - compare_to.mCoords[0] ) < tolerance ) &&
392 [ + - ][ + - ]: 78 : ( std::fabs( this->mCoords[1] - compare_to.mCoords[1] ) < tolerance ) &&
393 : 39 : ( std::fabs( this->mCoords[2] - compare_to.mCoords[2] ) < tolerance ) );
394 : : }
395 : :
396 : : // Length functions
397 : 4429605 : inline double Vector3D::length_squared() const
398 : : {
399 : 4429605 : return ( mCoords[0] * mCoords[0] + mCoords[1] * mCoords[1] + mCoords[2] * mCoords[2] );
400 : : }
401 : 3359297 : inline double Vector3D::length() const
402 : : {
403 : 3359297 : return std::sqrt( mCoords[0] * mCoords[0] + mCoords[1] * mCoords[1] + mCoords[2] * mCoords[2] );
404 : : }
405 : :
406 : : inline double inner_product( const Vector3D* v1, const Vector3D* v2, size_t n )
407 : : {
408 : : double result = 0.0;
409 : : const Vector3D* const end = v1 + n;
410 : : while( v1 < end )
411 : : {
412 : : result += *v1 % *v2;
413 : : ++v1;
414 : : ++v2;
415 : : }
416 : : return result;
417 : : }
418 : :
419 : 3357 : inline double length_squared( const Vector3D* v, int n )
420 : : {
421 : 3357 : double sum = 0.0;
422 [ + + ]: 2155143 : for( int i = 0; i < n; ++i )
423 : 2151786 : sum += v[i].length_squared();
424 : 3357 : return sum;
425 : : }
426 : 96054 : inline double length_squared( const std::vector< Vector3D >& v )
427 : : {
428 : 96054 : double sum = 0.0;
429 [ + + ]: 200885 : for( size_t i = 0; i < v.size(); ++i )
430 : 104831 : sum += v[i].length_squared();
431 : 96054 : return sum;
432 : : }
433 : :
434 : 2779 : inline double length( const Vector3D* v, int n ) // norm for an array of Vector3Ds
435 : : {
436 : 2779 : return std::sqrt( length_squared( v, n ) );
437 : : }
438 : : inline double length( const std::vector< Vector3D >& v )
439 : : {
440 : : return std::sqrt( length_squared( v ) );
441 : : }
442 : :
443 : 0 : inline double Linf( const Vector3D* v, int n ) // max entry for an array of Vector3Ds
444 : : {
445 : 0 : double max = 0;
446 : : // loop over the length of the array
447 [ # # ]: 0 : for( int i = 0; i < n; ++i )
448 : : {
449 [ # # ]: 0 : if( max < std::fabs( v[i][0] ) ) max = std::fabs( v[i][0] );
450 [ # # ]: 0 : if( max < std::fabs( v[i][1] ) ) max = std::fabs( v[i][1] );
451 [ # # ]: 0 : if( max < std::fabs( v[i][2] ) ) max = std::fabs( v[i][2] );
452 : : }
453 : : // return the value of the largest entry in the array
454 : 0 : return max;
455 : : }
456 : :
457 : 0 : inline double Linf( const std::vector< Vector3D >& v ) // max entry for an array of Vector3Ds
458 : : {
459 : 0 : double max = 0;
460 : : // loop over the length of the array
461 [ # # ]: 0 : for( size_t i = 0; i < v.size(); ++i )
462 : : {
463 [ # # ]: 0 : if( max < std::fabs( v[i][0] ) ) max = std::fabs( v[i][0] );
464 [ # # ]: 0 : if( max < std::fabs( v[i][1] ) ) max = std::fabs( v[i][1] );
465 [ # # ]: 0 : if( max < std::fabs( v[i][2] ) ) max = std::fabs( v[i][2] );
466 : : }
467 : : // return the value of the largest entry in the array
468 : 0 : return max;
469 : : }
470 : :
471 : 175114 : inline void Vector3D::set_length( const double new_length )
472 : : {
473 : 175114 : double factor = new_length / length();
474 : 175114 : *this *= factor;
475 : 175114 : }
476 : 175114 : inline void Vector3D::normalize()
477 : : {
478 : 175114 : set_length( 1.0 );
479 : 175114 : }
480 : :
481 : : // Utility functions.
482 : : inline Vector3D Vector3D::interpolate( const double param, const Vector3D& p1, const Vector3D& p2 )
483 : : {
484 : : return ( 1 - param ) * p1 + param * p2;
485 : : }
486 : :
487 : : inline bool operator==( const Vector3D& v1, const Vector3D& v2 )
488 : : {
489 : : return v1.mCoords[0] == v2.mCoords[0] && v1.mCoords[1] == v2.mCoords[1] && v1.mCoords[2] == v2.mCoords[2];
490 : : }
491 : :
492 : : inline bool operator!=( const Vector3D& v1, const Vector3D& v2 )
493 : : {
494 : : return v1.mCoords[0] != v2.mCoords[0] || v1.mCoords[1] != v2.mCoords[1] || v1.mCoords[2] != v2.mCoords[2];
495 : : }
496 : :
497 : : } // namespace MBMesquite
498 : :
499 : : #endif
|
