00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #ifndef SYM_TENZOR_H
00012 #define SYM_TENZOR_H
00013
00014
00015 #include "Message/Message.h"
00016 #include "Utility/PAssert.h"
00017 #include "PETE/IpplExpressions.h"
00018 #include "AppTypes/TSVMeta.h"
00019 #include "AppTypes/Tenzor.h"
00020
00021 #ifdef IPPL_USE_STANDARD_HEADERS
00022 #include <iostream>
00023 using namespace std;
00024 #else
00025 #include <iostream.h>
00026 #endif
00027
00028
00030
00031
00032
00034
00035
00036
00037
00038
00039
00040
00041 template<class T, unsigned D>
00042 class SymTenzor
00043 {
00044 public:
00045
00046 typedef T Element_t;
00047 enum { ElemDim = 2 };
00048 enum { Size = D*(D+1)/2 };
00049
00050
00051 SymTenzor() {
00052 TSV_MetaAssignScalar< SymTenzor<T,D>,T,OpAssign>::apply(*this,T(0));
00053 }
00054
00055
00056 class DontInitialize {};
00057 SymTenzor(DontInitialize) {}
00058
00059
00060 SymTenzor(const T& x00) {
00061 TSV_MetaAssignScalar< SymTenzor<T,D>,T,OpAssign>::apply(*this,x00);
00062 }
00063
00064
00065 SymTenzor(const T& x00, const T& x10, const T& x11) {
00066 PInsist(D==2, "Number of arguments does not match SymTenzor dimension!!");
00067 X[0] = x00; X[1] = x10; X[2] = x11;
00068 }
00069
00070 SymTenzor(const T& x00, const T& x10, const T& x11, const T& x20,
00071 const T& x21, const T& x22) {
00072 PInsist(D==3, "Number of arguments does not match SymTenzor dimension!!");
00073 X[0]= x00; X[1]= x10; X[2]= x11;
00074 X[3]= x20; X[4]= x21; X[5]= x22;
00075 }
00076
00077
00078 SymTenzor(const SymTenzor<T,D> &rhs) {
00079 TSV_MetaAssign< SymTenzor<T,D> , SymTenzor<T,D> ,OpAssign > ::
00080 apply(*this,rhs);
00081 }
00082
00083
00084
00085 SymTenzor(const Tenzor<T,D>& t) {
00086 for (int i=0; i<D; ++i) {
00087 (*this)(i,i) = t(i,i);
00088 for (int j=i+1; j<D; ++j)
00089 (*this)(i,j) = (t(i,j)+t(j,i))*0.5;
00090 }
00091 }
00092
00093
00094 ~SymTenzor() { };
00095
00096
00097 const SymTenzor<T,D>& operator= (const SymTenzor<T,D> &rhs) {
00098 TSV_MetaAssign< SymTenzor<T,D> , SymTenzor<T,D> ,OpAssign > ::
00099 apply(*this,rhs);
00100 return *this;
00101 }
00102 template<class T1>
00103 const SymTenzor<T,D>& operator= (const SymTenzor<T1,D> &rhs) {
00104 TSV_MetaAssign< SymTenzor<T,D> , SymTenzor<T1,D> ,OpAssign > ::
00105 apply(*this,rhs);
00106 return *this;
00107 }
00108 const SymTenzor<T,D>& operator= (const T& rhs) {
00109 TSV_MetaAssignScalar< SymTenzor<T,D> , T ,OpAssign > :: apply(*this,rhs);
00110 return *this;
00111 }
00112 const SymTenzor<T,D>& operator= (const Tenzor<T,D> &rhs) {
00113 for (int i=0; i<D; ++i) {
00114 (*this)(i,i) = rhs(i,i);
00115 for (int j=i+1; j<D; ++j)
00116 (*this)(i,j) = (rhs(i,j)+rhs(j,i))*0.5;
00117 }
00118 return *this;
00119 }
00120
00121
00122 template<class T1>
00123 SymTenzor<T,D>& operator+=(const SymTenzor<T1,D> &rhs)
00124 {
00125 TSV_MetaAssign< SymTenzor<T,D> , SymTenzor<T1,D> , OpAddAssign > ::
00126 apply(*this,rhs);
00127 return *this;
00128 }
00129 SymTenzor<T,D>& operator+=(const T& rhs)
00130 {
00131 TSV_MetaAssignScalar< SymTenzor<T,D> , T , OpAddAssign > ::
00132 apply(*this,rhs);
00133 return *this;
00134 }
00135
00136 template<class T1>
00137 SymTenzor<T,D>& operator-=(const SymTenzor<T1,D> &rhs)
00138 {
00139 TSV_MetaAssign< SymTenzor<T,D> , SymTenzor<T1,D> , OpSubtractAssign > ::
00140 apply(*this,rhs);
00141 return *this;
00142 }
00143 SymTenzor<T,D>& operator-=(const T& rhs)
00144 {
00145 TSV_MetaAssignScalar< SymTenzor<T,D> , T , OpSubtractAssign > ::
00146 apply(*this,rhs);
00147 return *this;
00148 }
00149
00150 template<class T1>
00151 SymTenzor<T,D>& operator*=(const SymTenzor<T1,D> &rhs)
00152 {
00153 TSV_MetaAssign< SymTenzor<T,D> , SymTenzor<T1,D> , OpMultipplyAssign > ::
00154 apply(*this,rhs);
00155 return *this;
00156 }
00157 SymTenzor<T,D>& operator*=(const T& rhs)
00158 {
00159 TSV_MetaAssignScalar< SymTenzor<T,D> , T , OpMultipplyAssign > ::
00160 apply(*this,rhs);
00161 return *this;
00162 }
00163
00164 template<class T1>
00165 SymTenzor<T,D>& operator/=(const SymTenzor<T1,D> &rhs)
00166 {
00167 TSV_MetaAssign< SymTenzor<T,D> , SymTenzor<T1,D> , OpDivideAssign > ::
00168 apply(*this,rhs);
00169 return *this;
00170 }
00171 SymTenzor<T,D>& operator/=(const T& rhs)
00172 {
00173 TSV_MetaAssignScalar< SymTenzor<T,D> , T , OpDivideAssign > ::
00174 apply(*this,rhs);
00175 return *this;
00176 }
00177
00178
00179
00180 void diagonal(const T& rhs) {
00181 for ( int i = 0 ; i < D ; i++ ) {
00182 X[((i+1)*i/2) + i] = rhs;
00183 }
00184 }
00185
00186 int len(void) const { return Size; }
00187 int size(void) const { return sizeof(*this); }
00188 int get_Size(void) const { return Size; }
00189
00190
00191
00192 Element_t operator()(unsigned int i, unsigned int j) const {
00193 int lo = i < j ? i : j;
00194 int hi = i > j ? i : j;
00195 return X[((hi+1)*hi/2) + lo];
00196 }
00197
00198 Element_t& operator()(unsigned int i, unsigned int j) {
00199 int lo = i < j ? i : j;
00200 int hi = i > j ? i : j;
00201 return X[((hi+1)*hi/2) + lo];
00202 }
00203
00204 Element_t& operator()(pair<int,int> a) {
00205 int i = a.first;
00206 int j = a.second;
00207 int lo = i < j ? i : j;
00208 int hi = i > j ? i : j;
00209 return X[((hi+1)*hi/2) + lo];
00210 }
00211
00212 Element_t operator()( pair<int,int> a) const {
00213 int i = a.first;
00214 int j = a.second;
00215 int lo = i < j ? i : j;
00216 int hi = i > j ? i : j;
00217 return X[((hi+1)*hi/2) + lo];
00218 }
00219
00220 Element_t HL(unsigned int hi, unsigned int lo) const {
00221 PAssert( hi >= lo );
00222 PAssert( hi<D );
00223 PAssert( lo<D );
00224 return X[hi*(hi+1)/2 + lo];
00225 }
00226 Element_t& HL(unsigned int hi, unsigned int lo) {
00227 PAssert( hi >= lo );
00228 PAssert( hi<D );
00229 PAssert( lo<D );
00230 return X[hi*(hi+1)/2 + lo];
00231 }
00232
00233 Element_t& operator[](unsigned int i) {
00234 PAssert (i < Size);
00235 return X[i];
00236 }
00237
00238 Element_t operator[](unsigned int i) const {
00239 PAssert (i < Size);
00240 return X[i];
00241 }
00242
00243
00244
00245
00246 Element_t& operator()(unsigned int i) {
00247 PAssert (i < Size);
00248 return X[i];
00249 }
00250
00251 Element_t operator()(unsigned int i) const {
00252 PAssert (i < Size);
00253 return X[i];
00254 }
00255
00256
00257
00258
00259 bool operator==(const SymTenzor<T,D>& that) const {
00260 return TSV_MetaCompareArrays<T,T,D*(D+1)/2>::apply(X,that.X);
00261 }
00262 bool operator!=(const SymTenzor<T,D>& that) const {
00263 return !(*this == that);
00264 }
00265
00266
00267
00268 Message& putMessage(Message& m) const {
00269 m.setCopy(true);
00270 ::putMessage(m, X, X + ((D*(D + 1)/2)));
00271 return m;
00272 }
00273
00274 Message& getMessage(Message& m) {
00275 ::getMessage(m, X, X + ((D*(D + 1)/2)));
00276 return m;
00277 }
00278
00279 private:
00280
00281
00282 T X[Size];
00283
00284 };
00285
00286
00288
00289
00290
00292
00293 template <class T, unsigned D>
00294 T trace(const SymTenzor<T,D> &rhs) {
00295 T result = 0.0;
00296 for (int i = 0 ; i < D ; i++ )
00297 result += rhs(i,i);
00298 return result;
00299 }
00300
00301 template <class T, unsigned D>
00302 SymTenzor<T,D> transpose(const SymTenzor<T,D> &rhs) {
00303 return rhs;
00304 }
00305
00306
00307 template <class T, unsigned D>
00308 T det(const SymTenzor<T,D>& rhs) {
00309 PInsist(D<3, "Tenzor det() function not implemented for D>3!");
00310 return T(-999999.999999);
00311 }
00312
00313 template <class T>
00314 T det(const SymTenzor<T,3>& rhs) {
00315 T result;
00316 result =
00317 rhs(0,0)*(rhs(1,1)*rhs(2,2) - rhs(1,2)*rhs(2,1)) +
00318 rhs(0,1)*(rhs(1,2)*rhs(2,0) - rhs(1,0)*rhs(2,2)) +
00319 rhs(0,2)*(rhs(1,0)*rhs(2,1) - rhs(1,1)*rhs(2,0));
00320 return result;
00321 }
00322
00323 template <class T>
00324 T det(const SymTenzor<T,2>& rhs) {
00325 T result;
00326 result = rhs(0,0)*rhs(1,1) - rhs(0,1)*rhs(1,0);
00327 return result;
00328 }
00329
00330 template <class T>
00331 T det(const SymTenzor<T,1>& rhs) {
00332 T result = rhs(0,0);
00333 return result;
00334 }
00335
00336 // cofactors() - pow(-1, i+j)*M(i,j), where M(i,j) is a minor of the tensor.
00337 // See, for example, Arfken, Mathematical Methods for Physicists, 2nd Edition,
00338 // p. 157 (the section where the determinant of a matrix is defined).
00339
00340 // Only implement for 1D, 2D, 3D:
00341
00342 template <class T, unsigned D>
00343 Tenzor<T,D> cofactors(const SymTenzor<T,D>& rhs) {
00344 PInsist(D<4, "SymTenzor cofactors() function not implemented for D>3!");
00345 return Tenzor<T,D>(-999999.999999);
00346 }
00347
00348 template <class T>
00349 Tenzor<T,3> cofactors(const SymTenzor<T,3>& rhs) {
00350 // check for typename keyword bug
00351 #if (defined(IPPL_SGI_CC_721_TYPENAME_BUG) || defined(__MWERKS__))
00352 Tenzor<T,3> result = Tenzor<T,3>::DontInitialize();
00353 #else
00354 Tenzor<T,3> result = typename Tenzor<T,3>::DontInitialize();
00355 #endif
00356 result(0,0) = rhs(1,1)*rhs(2,2) - rhs(1,2)*rhs(2,1);
00357 result(1,0) = rhs(0,2)*rhs(2,1) - rhs(0,1)*rhs(2,2);
00358 result(2,0) = rhs(0,1)*rhs(1,2) - rhs(1,1)*rhs(0,2);
00359 result(0,1) = rhs(2,0)*rhs(1,2) - rhs(1,0)*rhs(2,2);
00360 result(1,1) = rhs(0,0)*rhs(2,2) - rhs(0,2)*rhs(2,0);
00361 result(2,1) = rhs(1,0)*rhs(0,2) - rhs(0,0)*rhs(1,2);
00362 result(0,2) = rhs(1,0)*rhs(2,1) - rhs(2,0)*rhs(1,1);
00363 result(1,2) = rhs(0,1)*rhs(2,0) - rhs(0,0)*rhs(2,1);
00364 result(2,2) = rhs(0,0)*rhs(1,1) - rhs(1,0)*rhs(0,1);
00365 return result;
00366 }
00367
00368 template <class T>
00369 Tenzor<T,2> cofactors(const SymTenzor<T,2>& rhs) {
00370 // check for typename keyword bug
00371 #if (defined(IPPL_SGI_CC_721_TYPENAME_BUG) || defined(__MWERKS__))
00372 Tenzor<T,2> result = Tenzor<T,2>::DontInitialize();
00373 #else
00374 Tenzor<T,2> result = typename Tenzor<T,2>::DontInitialize();
00375 #endif
00376 result(0,0) = rhs(1,1);
00377 result(1,0) = -rhs(0,1);
00378 result(0,1) = -rhs(1,0);
00379 result(1,1) = rhs(0,0);
00380 return result;
00381 }
00382
00383 // For D=1, cofactor is the unit tensor, because det = single tensor element
00384 // value:
00385 template <class T>
00386 Tenzor<T,1> cofactors(const SymTenzor<T,1>& rhs) {
00387 Tenzor<T,1> result = Tenzor<T,1>(1);
00388 return result;
00389 }
00390
00391
00393 //
00394 // Unary Operators
00395 //
00397
00398 //----------------------------------------------------------------------
00399 // unary operator-
00400 template<class T, unsigned D>
00401 inline SymTenzor<T,D> operator-(const SymTenzor<T,D> &op)
00402 {
00403 return TSV_MetaUnary< SymTenzor<T,D> , OpUnaryMinus > :: apply(op);
00404 }
00405
00406 //----------------------------------------------------------------------
00407 // unary operator+
00408 template<class T, unsigned D>
00409 inline const SymTenzor<T,D> &operator+(const SymTenzor<T,D> &op)
00410 {
00411 return op;
00412 }
00413
00415 //
00416 // Binary Operators
00417 //
00419
00420 //
00421 // Elementwise operators.
00422 //
00423
00424 TSV_ELEMENTWISE_OPERATOR(SymTenzor,operator+,OpAdd)
00425 TSV_ELEMENTWISE_OPERATOR(SymTenzor,operator-,OpSubtract)
00426 TSV_ELEMENTWISE_OPERATOR(SymTenzor,operator*,OpMultipply)
00427 TSV_ELEMENTWISE_OPERATOR(SymTenzor,operator/,OpDivide)
00428 TSV_ELEMENTWISE_OPERATOR(SymTenzor,Min,FnMin)
00429 TSV_ELEMENTWISE_OPERATOR(SymTenzor,Max,FnMax)
00430
00431 TSV_ELEMENTWISE_OPERATOR2(SymTenzor,Tenzor,operator+,OpAdd)
00432 TSV_ELEMENTWISE_OPERATOR2(Tenzor,SymTenzor,operator+,OpAdd)
00433 TSV_ELEMENTWISE_OPERATOR2(SymTenzor,Tenzor,operator-,OpSubtract)
00434 TSV_ELEMENTWISE_OPERATOR2(Tenzor,SymTenzor,operator-,OpSubtract)
00435 TSV_ELEMENTWISE_OPERATOR2(SymTenzor,Tenzor,operator*,OpMultipply)
00436 TSV_ELEMENTWISE_OPERATOR2(Tenzor,SymTenzor,operator*,OpMultipply)
00437 TSV_ELEMENTWISE_OPERATOR2(SymTenzor,Tenzor,operator/,OpDivide)
00438 TSV_ELEMENTWISE_OPERATOR2(Tenzor,SymTenzor,operator/,OpDivide)
00439
00440 //----------------------------------------------------------------------
00441 // dot products.
00442 //----------------------------------------------------------------------
00443
00444 template < class T1, class T2, unsigned D >
00445 inline Tenzor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D>
00446 dot(const SymTenzor<T1,D> &lhs, const SymTenzor<T2,D> &rhs)
00447 {
00448 return TSV_MetaDot< SymTenzor<T1,D> , SymTenzor<T2,D> > :: apply(lhs,rhs);
00449 }
00450
00451 template < class T1, class T2, unsigned D >
00452 inline Tenzor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D>
00453 dot(const SymTenzor<T1,D> &lhs, const Tenzor<T2,D> &rhs)
00454 {
00455 return TSV_MetaDot< SymTenzor<T1,D> , Tenzor<T2,D> > :: apply(lhs,rhs);
00456 }
00457
00458 template < class T1, class T2, unsigned D >
00459 inline Tenzor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D>
00460 dot(const Tenzor<T1,D> &lhs, const SymTenzor<T2,D> &rhs)
00461 {
00462 return TSV_MetaDot< Tenzor<T1,D> , SymTenzor<T2,D> > :: apply(lhs,rhs);
00463 }
00464
00465 template < class T1, class T2, unsigned D >
00466 inline Vektor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D>
00467 dot(const Vektor<T1,D> &lhs, const SymTenzor<T2,D> &rhs)
00468 {
00469 return TSV_MetaDot< Vektor<T1,D> , SymTenzor<T2,D> > :: apply(lhs,rhs);
00470 }
00471
00472 template < class T1, class T2, unsigned D >
00473 inline Vektor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D>
00474 dot(const SymTenzor<T1,D> &lhs, const Vektor<T2,D> &rhs)
00475 {
00476 return TSV_MetaDot< SymTenzor<T1,D> , Vektor<T2,D> > :: apply(lhs,rhs);
00477 }
00478
00479 //----------------------------------------------------------------------
00480 // double dot products.
00481 //----------------------------------------------------------------------
00482
00483 template < class T1, class T2, unsigned D >
00484 inline typename PETEBinaryReturn<T1,T2,OpMultipply>::type
00485 dotdot(const SymTenzor<T1,D> &lhs, const SymTenzor<T2,D> &rhs)
00486 {
00487 return TSV_MetaDotDot< SymTenzor<T1,D> , SymTenzor<T2,D> > :: apply(lhs,rhs);
00488 }
00489
00490 template < class T1, class T2, unsigned D >
00491 inline typename PETEBinaryReturn<T1,T2,OpMultipply>::type
00492 dotdot(const SymTenzor<T1,D> &lhs, const Tenzor<T2,D> &rhs)
00493 {
00494 return TSV_MetaDotDot< SymTenzor<T1,D> , Tenzor<T2,D> > :: apply(lhs,rhs);
00495 }
00496
00497 template < class T1, class T2, unsigned D >
00498 inline typename PETEBinaryReturn<T1,T2,OpMultipply>::type
00499 dotdot(const Tenzor<T1,D> &lhs, const SymTenzor<T2,D> &rhs)
00500 {
00501 return TSV_MetaDotDot< Tenzor<T1,D> , SymTenzor<T2,D> > :: apply(lhs,rhs);
00502 }
00503
00504 //----------------------------------------------------------------------
00505 // I/O
00506 template<class T, unsigned D>
00507 ostream& operator<<(ostream& out, const SymTenzor<T,D>& rhs) {
00508 if (D >= 1) {
00509 for (int i=0; i<D; i++) {
00510 out << "(";
00511 for (int j=0; j<D-1; j++) {
00512 out << rhs(i,j) << " , ";
00513 }
00514 out << rhs(i,D-1) << ")";
00515 // I removed this. --TJW: if (i < D - 1) out << endl;
00516 }
00517 } else {
00518 out << "( " << rhs(0,0) << " )";
00519 }
00520 return out;
00521 }
00522
00523 #endif // SYM_TENZOR_H
00524
00525 /***************************************************************************
00526 * $RCSfile: SymTenzor.h,v $ $Author: adelmann $
00527 * $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:24 $
00528 * IPPL_VERSION_ID: $Id: SymTenzor.h,v 1.1.1.1 2003/01/23 07:40:24 adelmann Exp $
00529 ***************************************************************************/
00530