ESyS-Particle
4.0.1
|
00001 00002 // // 00003 // Copyright (c) 2003-2011 by The University of Queensland // 00004 // Earth Systems Science Computational Centre (ESSCC) // 00005 // http://www.uq.edu.au/esscc // 00006 // // 00007 // Primary Business: Brisbane, Queensland, Australia // 00008 // Licensed under the Open Software License version 3.0 // 00009 // http://www.opensource.org/licenses/osl-3.0.php // 00010 // // 00012 00013 #ifndef __MATRIX3_HPP 00014 #define __MATRIX3_HPP 00015 00016 #include "Foundation/Matrix3.h" 00017 #include "Foundation/vec3.h" 00018 00019 00020 MATRIX3_INLINE Matrix3::Matrix3() 00021 { 00022 for(int i=0;i<3;i++) 00023 { 00024 for(int j=0;j<3;j++) 00025 { 00026 m[i][j]=0; 00027 } 00028 } 00029 } 00030 00031 MATRIX3_INLINE Matrix3::Matrix3(const Vec3& V1, const Vec3& V2,const Vec3& V3) 00032 { 00033 m[0][0]=V1.data[0]; 00034 m[0][1]=V2.data[0]; 00035 m[0][2]=V3.data[0]; 00036 m[1][0]=V1.data[1]; 00037 m[1][1]=V2.data[1]; 00038 m[1][2]=V3.data[1]; 00039 m[2][0]=V1.data[2]; 00040 m[2][1]=V2.data[2]; 00041 m[2][2]=V3.data[2]; 00042 00043 } 00044 00045 MATRIX3_INLINE Matrix3::Matrix3(const double a[3][3]) 00046 { 00047 m[0][0]=a[0][0]; 00048 m[0][1]=a[0][1]; 00049 m[0][2]=a[0][2]; 00050 m[1][0]=a[1][0]; 00051 m[1][1]=a[1][1]; 00052 m[1][2]=a[1][2]; 00053 m[2][0]=a[2][0]; 00054 m[2][1]=a[2][1]; 00055 m[2][2]=a[2][2]; 00056 } 00057 00058 MATRIX3_INLINE Matrix3::Matrix3(const Matrix3& rhs) 00059 { 00060 for(int i=0;i<3;i++) 00061 { 00062 for(int j=0;j<3;j++) 00063 { 00064 m[i][j]=rhs.m[i][j]; 00065 } 00066 } 00067 } 00068 00069 MATRIX3_INLINE Matrix3::~Matrix3() 00070 { 00071 } 00072 00076 MATRIX3_INLINE double Matrix3::det() 00077 { 00078 return m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])+ 00079 m[0][1]*(m[1][2]*m[2][0]-m[1][0]*m[2][2])+ 00080 m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]); 00081 } 00082 00083 00084 MATRIX3_INLINE Matrix3 Matrix3::inv() 00085 { 00086 Matrix3 res=*this; 00087 00088 res.invert(); 00089 00090 return res; 00091 } 00092 00093 MATRIX3_INLINE void Matrix3::transpose() 00094 { 00095 double h; 00096 00097 h=m[1][0]; 00098 m[1][0]=m[0][1]; 00099 m[0][1]=h; 00100 h=m[2][0]; 00101 m[2][0]=m[0][2]; 00102 m[0][2]=h; 00103 h=m[2][1]; 00104 m[2][1]=m[1][2]; 00105 m[1][2]=h; 00106 } 00107 00108 MATRIX3_INLINE Matrix3 Matrix3::trans() const 00109 { 00110 Matrix3 res; 00111 00112 res.m[0][0]=m[0][0]; 00113 res.m[0][1]=m[1][0]; 00114 res.m[0][2]=m[2][0]; 00115 res.m[1][0]=m[0][1]; 00116 res.m[1][1]=m[1][1]; 00117 res.m[1][2]=m[2][1]; 00118 res.m[2][0]=m[0][2]; 00119 res.m[2][1]=m[1][2]; 00120 res.m[2][2]=m[2][2]; 00121 00122 return res; 00123 } 00124 00125 // 15 Flops ( 9 mult,6 add) 00126 MATRIX3_INLINE Vec3 Matrix3::operator *(const Vec3& V) const 00127 { 00128 double x=m[0][0]*V.data[0]+m[0][1]*V.data[1]+m[0][2]*V.data[2]; 00129 double y=m[1][0]*V.data[0]+m[1][1]*V.data[1]+m[1][2]*V.data[2]; 00130 double z=m[2][0]*V.data[0]+m[2][1]*V.data[1]+m[2][2]*V.data[2]; 00131 00132 return Vec3(x,y,z); 00133 } 00134 00135 // 9 Flops 00136 MATRIX3_INLINE Matrix3 Matrix3::operator *(double d) const 00137 { 00138 Matrix3 res; 00139 00140 res.m[0][0]=m[0][0]*d; 00141 res.m[0][1]=m[0][1]*d; 00142 res.m[0][2]=m[0][2]*d; 00143 res.m[1][0]=m[1][0]*d; 00144 res.m[1][1]=m[1][1]*d; 00145 res.m[1][2]=m[1][2]*d; 00146 res.m[2][0]=m[2][0]*d; 00147 res.m[2][1]=m[2][1]*d; 00148 res.m[2][2]=m[2][2]*d; 00149 00150 return res; 00151 } 00152 00153 // 9 Flops 00154 MATRIX3_INLINE Matrix3 Matrix3::operator /(double d) const 00155 { 00156 Matrix3 res; 00157 00158 res.m[0][0]=m[0][0]/d; 00159 res.m[0][1]=m[0][1]/d; 00160 res.m[0][2]=m[0][2]/d; 00161 res.m[1][0]=m[1][0]/d; 00162 res.m[1][1]=m[1][1]/d; 00163 res.m[1][2]=m[1][2]/d; 00164 res.m[2][0]=m[2][0]/d; 00165 res.m[2][1]=m[2][1]/d; 00166 res.m[2][2]=m[2][2]/d; 00167 00168 return res; 00169 } 00170 00171 MATRIX3_INLINE Matrix3 &Matrix3::operator=(const Matrix3 &other) 00172 { 00173 for(int i = 0; i < 3; i++) { 00174 for(int j = 0; j < 3; j++) { 00175 m[i][j] = other.m[i][j]; 00176 } 00177 } 00178 return *this; 00179 } 00180 00181 MATRIX3_INLINE bool Matrix3::operator==(const Matrix3 &other) const 00182 { 00183 for(int i = 0; i < 3; i++) { 00184 for(int j = 0; j < 3; j++) { 00185 if (m[i][j] != other.m[i][j]) { 00186 return false; 00187 } 00188 } 00189 } 00190 return true; 00191 } 00192 00193 MATRIX3_INLINE std::ostream &operator<<(std::ostream &oStream, const Matrix3 &m) 00194 { 00195 oStream << m(0,0); 00196 for(int i = 1; i < 9; i++) { 00197 oStream << " " << m(i/3, i%3); 00198 } 00199 return oStream; 00200 } 00201 00202 // 45 Flops (27mult, 18add) 00203 MATRIX3_INLINE Matrix3 Matrix3::operator *(const Matrix3& rhs) const 00204 { 00205 Matrix3 res; 00206 00207 for(int i=0;i<3;i++){ 00208 for(int j=0;j<3;j++){ 00209 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]; 00210 } 00211 } 00212 return res; 00213 } 00214 00215 // 9 Flops 00216 MATRIX3_INLINE Matrix3& Matrix3::operator+=(const Matrix3& rhs) 00217 { 00218 for(int i=0;i<3;i++){ 00219 for(int j=0;j<3;j++){ 00220 m[i][j]=m[i][j]+rhs.m[i][j]; 00221 } 00222 } 00223 return *this; 00224 } 00225 00229 MATRIX3_INLINE Matrix3 Matrix3::operator +(const Matrix3& rhs) const 00230 { 00231 Matrix3 res; 00232 00233 for(int i=0;i<3;i++){ 00234 for(int j=0;j<3;j++){ 00235 res.m[i][j]=m[i][j]+rhs.m[i][j]; 00236 } 00237 } 00238 00239 return res; 00240 } 00241 00245 MATRIX3_INLINE Matrix3 Matrix3::operator-(const Matrix3& rhs) const 00246 { 00247 Matrix3 res; 00248 00249 for(int i=0;i<3;i++) { 00250 for(int j=0;j<3;j++) { 00251 res.m[i][j]=m[i][j]-rhs.m[i][j]; 00252 } 00253 } 00254 00255 return res; 00256 } 00257 00261 MATRIX3_INLINE double Matrix3::trace() const 00262 { 00263 return m[0][0]+m[1][1]+m[2][2]; 00264 } 00265 00269 MATRIX3_INLINE double Matrix3::norm() const 00270 { 00271 double res=0.0; 00272 00273 for(int i=0;i<3;i++){ 00274 for(int j=0;j<3;j++){ 00275 res+=m[i][j]*m[i][j]; 00276 } 00277 } 00278 return res; 00279 } 00280 00281 // matrix from vector 00282 MATRIX3_INLINE Matrix3 star(const Vec3& V) 00283 { 00284 Matrix3 res; 00285 00286 res.m[0][1]=-V.Z(); 00287 res.m[0][2]=V.Y(); 00288 res.m[1][0]=V.Z(); 00289 res.m[1][2]=-V.X(); 00290 res.m[2][0]=-V.Y(); 00291 res.m[2][1]=V.X(); 00292 00293 return res; 00294 } 00295 00296 // unit matrix 00297 MATRIX3_INLINE Matrix3 Matrix3::Unit() 00298 { 00299 Matrix3 res; 00300 00301 res.m[0][0]=1.0; 00302 res.m[1][1]=1.0; 00303 res.m[2][2]=1.0; 00304 00305 return res; 00306 } 00307 00311 MATRIX3_INLINE Matrix3 operator*(double d,const Matrix3& M) 00312 { 00313 Matrix3 res; 00314 00315 for(int i=0;i<3;i++){ 00316 for(int j=0;j<3;j++){ 00317 res.m[i][j]=d*M.m[i][j]; 00318 } 00319 } 00320 00321 return res; 00322 } 00323 00324 00325 00326 #endif // __MATRIX3_HPP