ESyS-Particle  2.3.2
vec3.hpp
Go to the documentation of this file.
1 
2 // //
3 // Copyright (c) 2003-2014 by The University of Queensland //
4 // Centre for Geoscience Computing //
5 // http://earth.uq.edu.au/centre-geoscience-computing //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 
14 #ifndef __VEC3_HPP
15 #define __VEC3_HPP
16 
17 #include "Foundation/Matrix3.h"
18 
19 //the error...
20 VEC3_INLINE VecErr::VecErr(const string& m):MError(m)
21 {
22  message.insert(0,"Vec3 ");
23 }
24 
25 // constructors
27 {
28  data[0]=0;
29  data[1]=0;
30  data[2]=0;
31 }
32 
34 {
35  data[0]=s;
36  data[1]=s;
37  data[2]=s;
38 }
39 
40 VEC3_INLINE Vec3::Vec3(double a,double b,double c)
41 {
42  data[0]=a;
43  data[1]=b;
44  data[2]=c;
45 }
46 
48 {
49  data[0]=rhs.data[0];
50  data[1]=rhs.data[1];
51  data[2]=rhs.data[2];
52 }
53 
54 // operators
55 
57 {
58  data[0]=rhs.data[0];
59  data[1]=rhs.data[1];
60  data[2]=rhs.data[2];
61  return *this;
62 }
63 
65 {
66  data[0]=s;
67  data[1]=s;
68  data[2]=s;
69  return *this;
70 }
71 
73 {
74  data[0]-=rhs.data[0];
75  data[1]-=rhs.data[1];
76  data[2]-=rhs.data[2];
77  return *this;
78 }
79 
81 {
82  data[0]+=rhs.data[0];
83  data[1]+=rhs.data[1];
84  data[2]+=rhs.data[2];
85  return *this;
86 }
87 
89 {
90  return Vec3(data[0]+rhs.data[0], data[1]+rhs.data[1], data[2]+rhs.data[2]);
91 }
92 
94 {
95  return Vec3(data[0]-rhs.data[0], data[1]-rhs.data[1], data[2]-rhs.data[2]);
96 }
97 
99 {
100  return Vec3( -data[0],-data[1],-data[2] );
101 }
102 
104 {
105  const double x = m(0,0)*data[0] + m(1,0)*data[1] + m(2,0)*data[2];
106  const double y = m(0,1)*data[0] + m(1,1)*data[1] + m(2,1)*data[2];
107  const double z = m(0,2)*data[0] + m(1,2)*data[1] + m(2,2)*data[2];
108 
109  return Vec3(x,y,z);
110 }
111 
112 VEC3_INLINE double Vec3::operator*(const Vec3& rhs) const
113 {
114  return data[0]*rhs.data[0]+data[1]*rhs.data[1]+data[2]*rhs.data[2];
115 }
116 
118 {
119  return Vec3(data[0]*s,data[1]*s,data[2]*s) ;
120 }
121 
123 {
124  data[0]*=rhs;
125  data[1]*=rhs;
126  data[2]*=rhs;
127  return *this;
128 }
129 
131 {
132  data[0] /= c;
133  data[1] /= c;
134  data[2] /= c;
135 
136  return *this;
137 }
138 
140 {
141  return Vec3(data[0]/s,data[1]/s,data[2]/s) ;
142 }
143 
145 {
146  return Vec3(data[0]+s, data[1]+s, data[2]+s) ;
147 }
148 
150 {
151  return Vec3(data[0]-s, data[1]-s, data[2]-s);
152 }
153 
154 VEC3_INLINE Vec3 Vec3::rotate(const Vec3 &axis, const Vec3 &axisPt) const
155 {
156  const double phi = axis.norm();
157  if (phi > 0.0)
158  {
159  const Vec3 r = *this - axisPt;
160  const Vec3 n = axis/phi;
161  const double cosPhi = cos(phi);
162  const Vec3 rotatedR =
163  r*cosPhi + n*((dot(n, r))*(1-cosPhi)) + cross(r, n)*sin(phi);
164  return rotatedR + axisPt;
165  }
166  return *this;
167 }
168 
170 {
171  data[0] += s;
172  data[1] += s;
173  data[2] += s;
174  return *this;
175 }
176 
178 {
179  data[0] -= s;
180  data[1] -= s;
181  data[2] -= s;
182  return *this;
183 }
184 
185 // vector product
186 // 9 Flops ( 6 mult, 3 sub )
187 VEC3_INLINE Vec3 cross(const Vec3& lhs,const Vec3& rhs)
188 {
189  return Vec3(lhs.data[1]*rhs.data[2]-lhs.data[2]*rhs.data[1],
190  lhs.data[2]*rhs.data[0]-lhs.data[0]*rhs.data[2],
191  lhs.data[0]*rhs.data[1]-lhs.data[1]*rhs.data[0]);
192 }
193 
194 // dot product
195 
196 VEC3_INLINE double dot(const Vec3& v1, const Vec3& v2)
197 {
198  return v1.data[0] * v2.data[0] +
199  v1.data[1] * v2.data[1] +
200  v1.data[2] * v2.data[2];
201 }
202 
203 VEC3_INLINE Vec3 operator*(double f,const Vec3& rhs)
204 {
205  return Vec3(f*rhs.data[0], f*rhs.data[1], f*rhs.data[2]);
206 }
207 
208 
209 // euclidian norm
210 // 6 Flops ( 3 mult, 2 add, 1 sqrt )
211 VEC3_INLINE double Vec3::norm() const
212 {
213  return sqrt(data[0]*data[0]+data[1]*data[1]+data[2]*data[2]);
214 }
215 
216 // square of the euclidian norm
217 // 5 Flops ( 3 mult, 2 add)
218 VEC3_INLINE double Vec3::norm2() const
219 {
220  return data[0]*data[0]+data[1]*data[1]+data[2]*data[2];
221 }
222 
223 // returns unit vector in direction of the original vector
224 // 9 Flops ( 3 mult, 2 add, 3 div, 1 sqrt )
226 {
227  return (*this)/norm();
228 }
229 
230 // per element min/max
231 VEC3_INLINE Vec3 cmax(const Vec3& v1,const Vec3& v2)
232 {
233  Vec3 res;
234  res.data[0]=v1.data[0]>v2.data[0] ? v1.data[0] : v2.data[0];
235  res.data[1]=v1.data[1]>v2.data[1] ? v1.data[1] : v2.data[1];
236  res.data[2]=v1.data[2]>v2.data[2] ? v1.data[2] : v2.data[2];
237  return res;
238 }
239 
240 VEC3_INLINE Vec3 cmin(const Vec3& v1,const Vec3& v2)
241 {
242  Vec3 res;
243  res.data[0]=v1.data[0]<v2.data[0] ? v1.data[0] : v2.data[0];
244  res.data[1]=v1.data[1]<v2.data[1] ? v1.data[1] : v2.data[1];
245  res.data[2]=v1.data[2]<v2.data[2] ? v1.data[2] : v2.data[2];
246  return res;
247 }
248 
249 // save version, throws exception if norm()==0
251 {
252  double n=norm();
253  if(n==0) throw VecErr("norm() of data[2]ero-vector");
254  Vec3 res(data[0],data[1],data[2]);
255  return res/n;
256 }
257 
258 VEC3_INLINE double Vec3::max() const
259 {
260  double m = ( data[0]>data[1] ? data[0] : data[1] );
261 
262  return ( m>data[2] ? m : data[2] );
263 
264 }
265 
266 VEC3_INLINE double Vec3::min() const
267 {
268  double m = ( data[0]<data[1] ? data[0] : data[1] );
269 
270  return ( m<data[2] ? m : data[2] );
271 }
272 
273 VEC3_INLINE bool Vec3::operator==(const Vec3& V) const
274 {
275  return((data[0]==V.data[0])&&(data[1]==V.data[1])&&(data[2]==V.data[2]));
276 }
277 
278 VEC3_INLINE bool Vec3::operator!=(const Vec3& V) const
279 {
280  return((data[0]!=V.data[0])||(data[1]!=V.data[1])||(data[2]!=V.data[2]));
281 }
282 
283 // per component min/max
284 
285 VEC3_INLINE Vec3 comp_max(const Vec3& V1,const Vec3& V2)
286 {
287  double x=(V1.X() > V2.X()) ? V1.X() : V2.X();
288  double y=(V1.Y() > V2.Y()) ? V1.Y() : V2.Y();
289  double z=(V1.Z() > V2.Z()) ? V1.Z() : V2.Z();
290 
291  return Vec3(x,y,z);
292 }
293 
294 VEC3_INLINE Vec3 comp_min(const Vec3& V1,const Vec3& V2)
295 {
296  double x=(V1.X() < V2.X()) ? V1.X() : V2.X();
297  double y=(V1.Y() < V2.Y()) ? V1.Y() : V2.Y();
298  double z=(V1.Z() < V2.Z()) ? V1.Z() : V2.Z();
299 
300  return Vec3(x,y,z);
301 }
302 
303 // in/output
304 
305 VEC3_INLINE std::ostream& operator << (std::ostream& ostr,const Vec3& V)
306 {
307  const char delimiter = ' ';
308  ostr
309  << V.data[0] << delimiter
310  << V.data[1] << delimiter
311  << V.data[2];
312 
313  return ostr;
314 }
315 
316 VEC3_INLINE std::istream& operator >> (std::istream& istr,Vec3& V)
317 {
318  istr
319  >> V.data[0]
320  >> V.data[1]
321  >> V.data[2];
322 
323  return istr;
324 }
325 
326 #endif // __VEC3_HPP