ESyS-Particle  2.3.2
Matrix3.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 #ifndef __MATRIX3_HPP
14 #define __MATRIX3_HPP
15 
16 #include "Foundation/Matrix3.h"
17 #include "Foundation/vec3.h"
18 
19 
21 {
22  for(int i=0;i<3;i++)
23  {
24  for(int j=0;j<3;j++)
25  {
26  m[i][j]=0;
27  }
28  }
29 }
30 
31 MATRIX3_INLINE Matrix3::Matrix3(const Vec3& V1, const Vec3& V2,const Vec3& V3)
32 {
33  m[0][0]=V1.data[0];
34  m[0][1]=V2.data[0];
35  m[0][2]=V3.data[0];
36  m[1][0]=V1.data[1];
37  m[1][1]=V2.data[1];
38  m[1][2]=V3.data[1];
39  m[2][0]=V1.data[2];
40  m[2][1]=V2.data[2];
41  m[2][2]=V3.data[2];
42 
43 }
44 
45 MATRIX3_INLINE Matrix3::Matrix3(const double a[3][3])
46 {
47  m[0][0]=a[0][0];
48  m[0][1]=a[0][1];
49  m[0][2]=a[0][2];
50  m[1][0]=a[1][0];
51  m[1][1]=a[1][1];
52  m[1][2]=a[1][2];
53  m[2][0]=a[2][0];
54  m[2][1]=a[2][1];
55  m[2][2]=a[2][2];
56 }
57 
59 {
60  for(int i=0;i<3;i++)
61  {
62  for(int j=0;j<3;j++)
63  {
64  m[i][j]=rhs.m[i][j];
65  }
66  }
67 }
68 
70 {
71 }
72 
77 {
78  return m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])+
79  m[0][1]*(m[1][2]*m[2][0]-m[1][0]*m[2][2])+
80  m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
81 }
82 
83 
85 {
86  Matrix3 res=*this;
87 
88  res.invert();
89 
90  return res;
91 }
92 
94 {
95  double h;
96 
97  h=m[1][0];
98  m[1][0]=m[0][1];
99  m[0][1]=h;
100  h=m[2][0];
101  m[2][0]=m[0][2];
102  m[0][2]=h;
103  h=m[2][1];
104  m[2][1]=m[1][2];
105  m[1][2]=h;
106 }
107 
109 {
110  Matrix3 res;
111 
112  res.m[0][0]=m[0][0];
113  res.m[0][1]=m[1][0];
114  res.m[0][2]=m[2][0];
115  res.m[1][0]=m[0][1];
116  res.m[1][1]=m[1][1];
117  res.m[1][2]=m[2][1];
118  res.m[2][0]=m[0][2];
119  res.m[2][1]=m[1][2];
120  res.m[2][2]=m[2][2];
121 
122  return res;
123 }
124 
125 // 15 Flops ( 9 mult,6 add)
127 {
128  double x=m[0][0]*V.data[0]+m[0][1]*V.data[1]+m[0][2]*V.data[2];
129  double y=m[1][0]*V.data[0]+m[1][1]*V.data[1]+m[1][2]*V.data[2];
130  double z=m[2][0]*V.data[0]+m[2][1]*V.data[1]+m[2][2]*V.data[2];
131 
132  return Vec3(x,y,z);
133 }
134 
135 // 9 Flops
137 {
138  Matrix3 res;
139 
140  res.m[0][0]=m[0][0]*d;
141  res.m[0][1]=m[0][1]*d;
142  res.m[0][2]=m[0][2]*d;
143  res.m[1][0]=m[1][0]*d;
144  res.m[1][1]=m[1][1]*d;
145  res.m[1][2]=m[1][2]*d;
146  res.m[2][0]=m[2][0]*d;
147  res.m[2][1]=m[2][1]*d;
148  res.m[2][2]=m[2][2]*d;
149 
150  return res;
151 }
152 
153 // 9 Flops
155 {
156  Matrix3 res;
157 
158  res.m[0][0]=m[0][0]/d;
159  res.m[0][1]=m[0][1]/d;
160  res.m[0][2]=m[0][2]/d;
161  res.m[1][0]=m[1][0]/d;
162  res.m[1][1]=m[1][1]/d;
163  res.m[1][2]=m[1][2]/d;
164  res.m[2][0]=m[2][0]/d;
165  res.m[2][1]=m[2][1]/d;
166  res.m[2][2]=m[2][2]/d;
167 
168  return res;
169 }
170 
172 {
173  for(int i = 0; i < 3; i++) {
174  for(int j = 0; j < 3; j++) {
175  m[i][j] = other.m[i][j];
176  }
177  }
178  return *this;
179 }
180 
182 {
183  for(int i = 0; i < 3; i++) {
184  for(int j = 0; j < 3; j++) {
185  if (m[i][j] != other.m[i][j]) {
186  return false;
187  }
188  }
189  }
190  return true;
191 }
192 
193 MATRIX3_INLINE std::ostream &operator<<(std::ostream &oStream, const Matrix3 &m)
194 {
195  oStream << m(0,0);
196  for(int i = 1; i < 9; i++) {
197  oStream << " " << m(i/3, i%3);
198  }
199  return oStream;
200 }
201 
202 // 45 Flops (27mult, 18add)
204 {
205  Matrix3 res;
206 
207  for(int i=0;i<3;i++){
208  for(int j=0;j<3;j++){
209  res.m[i][j]=m[i][0]*rhs.m[0][j]+m[i][1]*rhs.m[1][j]+m[i][2]*rhs.m[2][j];
210  }
211  }
212  return res;
213 }
214 
215 // 9 Flops
217 {
218  for(int i=0;i<3;i++){
219  for(int j=0;j<3;j++){
220  m[i][j]=m[i][j]+rhs.m[i][j];
221  }
222  }
223  return *this;
224 }
225 
230 {
231  Matrix3 res;
232 
233  for(int i=0;i<3;i++){
234  for(int j=0;j<3;j++){
235  res.m[i][j]=m[i][j]+rhs.m[i][j];
236  }
237  }
238 
239  return res;
240 }
241 
246 {
247  Matrix3 res;
248 
249  for(int i=0;i<3;i++) {
250  for(int j=0;j<3;j++) {
251  res.m[i][j]=m[i][j]-rhs.m[i][j];
252  }
253  }
254 
255  return res;
256 }
257 
262 {
263  return m[0][0]+m[1][1]+m[2][2];
264 }
265 
270 {
271  double res=0.0;
272 
273  for(int i=0;i<3;i++){
274  for(int j=0;j<3;j++){
275  res+=m[i][j]*m[i][j];
276  }
277  }
278  return res;
279 }
280 
281 // matrix from vector
283 {
284  Matrix3 res;
285 
286  res.m[0][1]=-V.Z();
287  res.m[0][2]=V.Y();
288  res.m[1][0]=V.Z();
289  res.m[1][2]=-V.X();
290  res.m[2][0]=-V.Y();
291  res.m[2][1]=V.X();
292 
293  return res;
294 }
295 
296 // unit matrix
298 {
299  Matrix3 res;
300 
301  res.m[0][0]=1.0;
302  res.m[1][1]=1.0;
303  res.m[2][2]=1.0;
304 
305  return res;
306 }
307 
312 {
313  Matrix3 res;
314 
315  for(int i=0;i<3;i++){
316  for(int j=0;j<3;j++){
317  res.m[i][j]=d*M.m[i][j];
318  }
319  }
320 
321  return res;
322 }
323 
324 
325 
326 #endif // __MATRIX3_HPP