HG-MD
1
|
00001 #ifndef MATRIX_H 00002 #define MATRIX_H 00003 00004 #include <cmath> 00005 #include <sstream> 00006 #include "Vector.h" 00007 00010 class MatrixSymmetric3D 00011 { 00012 public: 00013 00014 Mdouble XX, XY, XZ, YY, YZ, ZZ; 00015 00016 inline MatrixSymmetric3D(void){} 00017 00018 inline MatrixSymmetric3D (const Mdouble xx, const Mdouble xy, const Mdouble xz, const Mdouble yy, const Mdouble yz, const Mdouble zz) 00019 { 00020 XX = xx; XY = xy; XZ = xz; 00021 YY = yy; YZ = yz; 00022 ZZ = zz; 00023 } 00024 00025 inline void set_zero () 00026 { 00027 XX = XY = XZ = YY = YZ = ZZ = 0.0; 00028 } 00029 00030 inline Mdouble trace () const 00031 { 00032 return (XX+YY+ZZ)/3; 00033 } 00034 00035 inline MatrixSymmetric3D operator + (const MatrixSymmetric3D &A) const 00036 { 00037 return MatrixSymmetric3D(XX + A.XX, XY + A.XY, XZ + A.XZ, 00038 YY + A.YY, YZ + A.YZ, ZZ + A.ZZ); 00039 } 00040 00041 inline MatrixSymmetric3D operator - (const MatrixSymmetric3D &A) const 00042 { 00043 return MatrixSymmetric3D(XX - A.XX, XY - A.XY, XZ - A.XZ, 00044 YY - A.YY, YZ - A.YZ, ZZ - A.ZZ); 00045 } 00046 00047 inline MatrixSymmetric3D operator + (const Mdouble A) const 00048 { 00049 return MatrixSymmetric3D(XX + A, XY + A, XZ + A, 00050 YY + A, YZ + A, ZZ + A); 00051 } 00052 00053 inline MatrixSymmetric3D operator - (const Mdouble A) const 00054 { 00055 return MatrixSymmetric3D(XX - A, XY - A, XZ - A, 00056 YY - A, YZ - A, ZZ - A); 00057 } 00058 00059 inline MatrixSymmetric3D operator * (const Mdouble A) const 00060 { 00061 return MatrixSymmetric3D(XX * A, XY * A, XZ * A, 00062 YY * A, YZ * A, ZZ * A); 00063 } 00064 00065 inline MatrixSymmetric3D operator / (const Mdouble A) const 00066 { 00067 return MatrixSymmetric3D(XX / A, XY / A, XZ / A, 00068 YY / A, YZ / A, ZZ / A); 00069 } 00070 00071 friend inline std::ostream& operator<<(std::ostream& os, const MatrixSymmetric3D &A) 00072 { 00073 os << A.XX << ' ' << A.XY << ' ' << A.XZ << " " << A.YY << ' ' << A.YZ << " " << A.ZZ; 00074 return os; 00075 } 00076 00077 friend inline std::istream& operator>>(std::istream& is, MatrixSymmetric3D &A) 00078 { 00079 is >> A.XX >> A.XY >> A.XZ >> A.YY >> A.YZ >> A.ZZ; 00080 return is; 00081 } 00082 00083 inline MatrixSymmetric3D& operator+=(const MatrixSymmetric3D &A) 00084 { 00085 XX += A.XX; 00086 XY += A.XY; 00087 XZ += A.XZ; 00088 YY += A.YY; 00089 YZ += A.YZ; 00090 ZZ += A.ZZ; 00091 return *this; 00092 } 00093 00094 inline MatrixSymmetric3D& operator-=(const MatrixSymmetric3D &A) 00095 { 00096 XX -= A.XX; 00097 XY -= A.XY; 00098 XZ -= A.XZ; 00099 YY -= A.YY; 00100 YZ -= A.YZ; 00101 ZZ -= A.ZZ; 00102 return *this; 00103 } 00104 00105 inline MatrixSymmetric3D& operator/=(const Mdouble a) 00106 { 00107 XX /= a; 00108 XY /= a; 00109 XZ /= a; 00110 YY /= a; 00111 YZ /= a; 00112 ZZ /= a; 00113 return *this; 00114 } 00115 00116 // Pointwise square 00117 friend inline MatrixSymmetric3D square(const MatrixSymmetric3D &A) 00118 { 00119 return MatrixSymmetric3D(sqr(A.XX), sqr(A.XY), sqr(A.XZ), sqr(A.YY), sqr(A.YZ), sqr(A.ZZ)); 00120 } 00121 00122 // Pointwise square root 00123 friend inline MatrixSymmetric3D sqrt(const MatrixSymmetric3D &A) 00124 { 00125 return MatrixSymmetric3D(sqrt(A.XX), sqrt(A.XY), sqrt(A.XZ), sqrt(A.YY), sqrt(A.YZ), sqrt(A.ZZ)); 00126 } 00127 00128 //calculates A*A' 00129 //friend inline MatrixSymmetric3D Dyadic(Vec3D A) { 00130 // return MatrixSymmetric3D(A.X*A.X, A.X*A.Y, A.X*A.Z, A.Y*A.Y, A.Y*A.Z, A.Z*A.Z); 00131 //} 00132 }; 00133 00136 00137 MatrixSymmetric3D SymmetrizedDyadic(Vec3D A, Vec3D B) { 00138 return MatrixSymmetric3D(A.X*B.X, 0.5*(A.X*B.Y+B.X*A.Y), 0.5*(A.X*B.Z+B.X*A.Z), A.Y*B.Y, 0.5*(A.Y*B.Z+B.Y*A.Z), A.Z*B.Z); 00139 } 00140 00141 00144 class Matrix3D 00145 { 00146 public: 00147 00148 Mdouble XX, XY, XZ, YX, YY, YZ, ZX, ZY, ZZ; 00149 00150 inline Matrix3D(void){} 00151 00152 inline Matrix3D (const Mdouble xx, const Mdouble xy, const Mdouble xz, const Mdouble yx, const Mdouble yy, const Mdouble yz, const Mdouble zx, const Mdouble zy, const Mdouble zz) 00153 { 00154 XX = xx; XY = xy; XZ = xz; 00155 YX = yx; YY = yy; YZ = yz; 00156 ZX = zx; ZY = zy; ZZ = zz; 00157 } 00158 00159 inline void set_zero () 00160 { 00161 XX = XY = XZ = YX = YY = YZ = ZX = ZY = ZZ = 0.0; 00162 } 00163 00164 inline Mdouble trace () const 00165 { 00166 return (XX+YY+ZZ)/3; 00167 } 00168 00169 inline Matrix3D operator + (const Matrix3D &A) const 00170 { 00171 return Matrix3D(XX + A.XX, XY + A.XY, XZ + A.XZ, 00172 YX + A.YX, YY + A.YY, YZ + A.YZ, 00173 ZX + A.ZX, ZY + A.ZY, ZZ + A.ZZ); 00174 } 00175 00176 inline Matrix3D operator - (const Matrix3D &A) const 00177 { 00178 return Matrix3D(XX - A.XX, XY - A.XY, XZ - A.XZ, 00179 YX - A.YX, YY - A.YY, YZ - A.YZ, 00180 ZX - A.ZX, ZY - A.ZY, ZZ - A.ZZ); 00181 } 00182 00183 inline Matrix3D operator + (const Mdouble A) const 00184 { 00185 return Matrix3D(XX + A, XY + A, XZ + A, 00186 YX + A, YY + A, YZ + A, 00187 ZX + A, ZY + A, ZZ + A); 00188 } 00189 00190 inline Matrix3D operator - (const Mdouble A) const 00191 { 00192 return Matrix3D(XX - A, XY - A, XZ - A, 00193 YX - A, YY - A, YZ - A, 00194 ZX - A, ZY - A, ZZ - A); 00195 } 00196 00197 inline Matrix3D operator * (const Mdouble A) const 00198 { 00199 return Matrix3D(XX * A, XY * A, XZ * A, 00200 YX * A, YY * A, YZ * A, 00201 ZX * A, ZY * A, ZZ * A); 00202 } 00203 00204 inline Vec3D operator * (const Vec3D A) const 00205 { 00206 return Vec3D(XX * A.X + XY * A.Y + XZ * A.Z, 00207 YX * A.X + YY * A.Y + YZ * A.Z, 00208 ZX * A.X + ZY * A.Y + ZZ * A.Z ); 00209 } 00210 00211 inline Matrix3D operator / (const Mdouble A) const 00212 { 00213 return Matrix3D(XX / A, XY / A, XZ / A, 00214 YX / A, YY / A, YZ / A, 00215 ZX / A, ZY / A, ZZ / A); 00216 } 00217 00218 friend inline std::ostream& operator<<(std::ostream& os, const Matrix3D &A) 00219 { 00220 os << A.XX << ' ' << A.XY << ' ' << A.XZ << ' ' 00221 << A.YX << ' ' << A.YY << ' ' << A.YZ << ' ' 00222 << A.ZX << ' ' << A.ZY << ' ' << A.ZZ; 00223 return os; 00224 } 00225 00226 friend inline std::istream& operator>>(std::istream& is, Matrix3D &A) 00227 { 00228 is >> A.XX >> A.XY >> A.XZ >> A.YX >> A.YY >> A.YZ >> A.ZX >> A.ZY >> A.ZZ; 00229 return is; 00230 } 00231 00232 inline Matrix3D& operator+=(const Matrix3D &A) 00233 { 00234 XX += A.XX; 00235 XY += A.XY; 00236 XZ += A.XZ; 00237 YX += A.YX; 00238 YY += A.YY; 00239 YZ += A.YZ; 00240 ZX += A.ZX; 00241 ZY += A.ZY; 00242 ZZ += A.ZZ; 00243 return *this; 00244 } 00245 00246 inline Matrix3D& operator-=(const Matrix3D &A) 00247 { 00248 XX -= A.XX; 00249 XY -= A.XY; 00250 XZ -= A.XZ; 00251 YX -= A.YX; 00252 YY -= A.YY; 00253 YZ -= A.YZ; 00254 ZX -= A.ZX; 00255 ZY -= A.ZY; 00256 ZZ -= A.ZZ; 00257 return *this; 00258 } 00259 00260 inline Matrix3D& operator/=(const Mdouble A) 00261 { 00262 XX /= A; 00263 XY /= A; 00264 XZ /= A; 00265 YX /= A; 00266 YY /= A; 00267 YZ /= A; 00268 ZX /= A; 00269 ZY /= A; 00270 ZZ /= A; 00271 return *this; 00272 } 00273 00274 // Pointwise square 00275 friend inline Matrix3D square(const Matrix3D &A) 00276 { 00277 return Matrix3D(sqr(A.XX), sqr(A.XY), sqr(A.XZ), 00278 sqr(A.YX), sqr(A.YY), sqr(A.YZ), 00279 sqr(A.ZX), sqr(A.ZY), sqr(A.ZZ)); 00280 } 00281 00282 // Pointwise square root 00283 friend inline Matrix3D sqrt(const Matrix3D &A) 00284 { 00285 return Matrix3D(sqrt(A.XX), sqrt(A.XY), sqrt(A.XZ), 00286 sqrt(A.YX), sqrt(A.YY), sqrt(A.YZ), 00287 sqrt(A.ZX), sqrt(A.ZY), sqrt(A.ZZ)); 00288 } 00289 00290 }; 00291 00293 Matrix3D Dyadic(Vec3D A, Vec3D B) { 00294 return Matrix3D(A.X*B.X, A.X*B.Y, A.X*B.Z, 00295 A.Y*B.X, A.Y*B.Y, A.Y*B.Z, 00296 A.Z*B.X, A.Z*B.Y, A.Z*B.Z); 00297 } 00298 00299 00300 #endif