 |
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 #ifndef MB_CART_VECT_HPP
2 #define MB_CART_VECT_HPP
3
4 #include <cmath>
5 #include <iosfwd>
6 #include <cfloat>
7
8 namespace moab
9 {
10
11 /**
12 * \brief Cartesian Vector
13 */
14 class CartVect
15 {
16 private:
17 double d[3];
18
19 public:
20 inline CartVect() {}
21 /**Initialze all three values to same scalar (typically zero)*/
22 explicit inline CartVect( double v )
23 {
24 d[0] = d[1] = d[2] = v;
25 }
26 inline CartVect( double i, double j, double k )
27 {
28 d[0] = i;
29 d[1] = j;
30 d[2] = k;
31 }
32 /**Initialze from array*/
33 explicit inline CartVect( const double a[3] )
34 {
35 d[0] = a[0];
36 d[1] = a[1];
37 d[2] = a[2];
38 }
39 inline CartVect& operator=( const double v[3] )
40 {
41 d[0] = v[0];
42 d[1] = v[1];
43 d[2] = v[2];
44 return *this;
45 }
46
47 inline double& operator[]( unsigned i )
48 {
49 return d[i];
50 }
51 inline double operator[]( unsigned i ) const
52 {
53 return d[i];
54 }
55
56 inline CartVect& operator+=( const CartVect& v )
57 {
58 d[0] += v.d[0];
59 d[1] += v.d[1];
60 d[2] += v.d[2];
61 return *this;
62 }
63 inline CartVect& operator-=( const CartVect& v )
64 {
65 d[0] -= v.d[0];
66 d[1] -= v.d[1];
67 d[2] -= v.d[2];
68 return *this;
69 }
70 /** Assign cross product to this */
71 inline CartVect& operator*=( const CartVect& v );
72
73 inline CartVect& operator+=( double s )
74 {
75 d[0] += s;
76 d[1] += s;
77 d[2] += s;
78 return *this;
79 }
80 inline CartVect& operator-=( double s )
81 {
82 d[0] -= s;
83 d[1] -= s;
84 d[2] -= s;
85 return *this;
86 }
87 inline CartVect& operator*=( double s )
88 {
89 d[0] *= s;
90 d[1] *= s;
91 d[2] *= s;
92 return *this;
93 }
94 inline CartVect& operator/=( double s )
95 {
96 d[0] /= s;
97 d[1] /= s;
98 d[2] /= s;
99 return *this;
100 }
101 inline bool operator==( const CartVect& v ) const
102 {
103 return d[0] == v[0] && d[1] == v[1] && d[2] == v[2];
104 }
105 inline bool operator==( double val ) const
106 {
107 return d[0] == val && d[1] == val && d[2] == val;
108 }
109
110 inline double length() const; //!< vector length
111
112 inline double length_squared() const;
113
114 inline void normalize(); //!< make unit length, or 0 if length < DBL_MIN
115
116 inline void flip(); //!< flip direction
117
118 /** per-element scalar multiply (this[0] *= v[0], this[1] *= v[1], ...) */
119 inline void scale( const CartVect& v )
120 {
121 d[0] *= v.d[0];
122 d[1] *= v.d[1];
123 d[2] *= v.d[2];
124 }
125
126 // get pointer to array of three doubles
127 inline double* array()
128 {
129 return d;
130 }
131 inline const double* array() const
132 {
133 return d;
134 }
135
136 /** initialize double array from this */
137 inline void get( double v[3] ) const
138 {
139 v[0] = d[0];
140 v[1] = d[1];
141 v[2] = d[2];
142 }
143
144 /** initialize float array from this */
145 inline void get( float v[3] ) const
146 {
147 v[0] = static_cast< float >( d[0] );
148 v[1] = static_cast< float >( d[1] );
149 v[2] = static_cast< float >( d[2] );
150 }
151 };
152
153 inline CartVect operator+( const CartVect& u, const CartVect& v )
154 {
155 return CartVect( u[0] + v[0], u[1] + v[1], u[2] + v[2] );
156 }
157
158 inline CartVect operator-( const CartVect& u, const CartVect& v )
159 {
160 return CartVect( u[0] - v[0], u[1] - v[1], u[2] - v[2] );
161 }
162
163 /** cross product */
164 inline CartVect operator*( const CartVect& u, const CartVect& v )
165 {
166 return CartVect( u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0] );
167 }
168
169 //! Dot product
170 inline double operator%( const CartVect& u, const CartVect& v )
171 {
172 return u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
173 }
174
175 inline CartVect& CartVect::operator*=( const CartVect& v )
176 {
177 return *this = *this * v;
178 }
179
180 inline double CartVect::length() const
181 {
182 return std::sqrt( *this % *this );
183 }
184
185 inline double CartVect::length_squared() const
186 {
187 return d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
188 }
189
190 inline void CartVect::normalize()
191 {
192 double tmp = length();
193 if( tmp < DBL_MIN )
194 d[0] = d[1] = d[2] = 0;
195 else
196 *this /= tmp;
197 }
198
199 inline void CartVect::flip()
200 {
201 d[0] = -d[0];
202 d[1] = -d[1];
203 d[2] = -d[2];
204 }
205
206 //! Interior angle between two vectors
207 inline double angle( const CartVect& u, const CartVect& v )
208 {
209 double tmp = ( u % v ) / std::sqrt( ( u % u ) * ( v % v ) );
210 if( tmp > 1. ) tmp = 1.;
211 if( tmp < -1. ) tmp = -1.;
212 return std::acos( tmp );
213 }
214
215 //! Interior angle between two vectors
216 inline double angle_robust( CartVect u, CartVect v )
217 {
218 u.normalize();
219 v.normalize();
220
221 double tmp = ( u % v ) ;
222 if( tmp - 1. >= - 1.e-12 )
223 {
224 double dist = (u - v).length();
225 return dist; // approximate sin(x) with x for very small dist
226 }
227
228 if( tmp < -1. ) tmp = -1.;
229 return std::acos( tmp );
230 }
231
232 inline CartVect operator-( const CartVect& v )
233 {
234 return CartVect( -v[0], -v[1], -v[2] );
235 }
236 inline CartVect operator+( const CartVect& v, double s )
237 {
238 return CartVect( v[0] + s, v[1] + s, v[2] + s );
239 }
240 inline CartVect operator-( const CartVect& v, double s )
241 {
242 return CartVect( v[0] - s, v[1] - s, v[2] - s );
243 }
244 inline CartVect operator*( const CartVect& v, double s )
245 {
246 return CartVect( v[0] * s, v[1] * s, v[2] * s );
247 }
248 inline CartVect operator/( const CartVect& v, double s )
249 {
250 return CartVect( v[0] / s, v[1] / s, v[2] / s );
251 }
252 inline CartVect operator+( double s, const CartVect& v )
253 {
254 return CartVect( v[0] + s, v[1] + s, v[2] + s );
255 }
256 inline CartVect operator-( double s, const CartVect& v )
257 {
258 return CartVect( v[0] - s, v[1] - s, v[2] - s );
259 }
260 inline CartVect operator*( double s, const CartVect& v )
261 {
262 return CartVect( v[0] * s, v[1] * s, v[2] * s );
263 }
264
265 //! Get unit vector in same direction
266 inline CartVect unit( const CartVect& v )
267 {
268 const double len = v.length();
269 return CartVect( v[0] / len, v[1] / len, v[2] / len );
270 }
271
272 std::ostream& operator<<( std::ostream& s, const CartVect& v );
273
274 } // namespace moab
275
276 #endif
1.9.1.