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