Tenzor.h 13.3 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4
// -*- C++ -*-
/***************************************************************************
 *
 * The IPPL Framework
kraus's avatar
kraus committed
5
 *
gsell's avatar
gsell committed
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
 *
 * Visit http://people.web.psi.ch/adelmann/ for more details
 *
 ***************************************************************************/

#ifndef	TENZOR_H
#define	TENZOR_H

// include files
#include "Utility/PAssert.h"
#include "Message/Message.h"
#include "PETE/IpplExpressions.h"
#include "AppTypes/TSVMeta.h"

#include <iostream>

// forward declarations
template <class T, unsigned D> class SymTenzor;
template <class T, unsigned D> class AntiSymTenzor;


//////////////////////////////////////////////////////////////////////
//
// Definition of class Tenzor.
//
//////////////////////////////////////////////////////////////////////

template<class T, unsigned D>
class Tenzor
{
public:

  typedef T Element_t;
  enum { ElemDim = 2 };
  enum { Size = D*D };

kraus's avatar
kraus committed
42
  // Default Constructor
gsell's avatar
gsell committed
43 44 45 46 47 48 49 50
  Tenzor() {
    TSV_MetaAssignScalar<Tenzor<T,D>,T,OpAssign>::apply(*this,T(0));
  }

  // A noninitializing ctor.
  class DontInitialize {};
  Tenzor(DontInitialize) {}

kraus's avatar
kraus committed
51
  // Copy Constructor
gsell's avatar
gsell committed
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
  Tenzor(const Tenzor<T,D> &rhs) {
    TSV_MetaAssign< Tenzor<T,D> , Tenzor<T,D> ,OpAssign >::apply(*this,rhs);
  }

  // constructor from a single T
  Tenzor(const T& x00) {
    TSV_MetaAssignScalar< Tenzor<T,D> , T ,OpAssign >::apply(*this,x00);
  }

  // constructors for fixed dimension
  Tenzor(const T& x00, const T& x10, const T& x01, const T& x11) {
    PInsist(D==2, "Number of arguments does not match Tenzor dimension!!");
    X[0] = x00;
    X[1] = x10;
    X[2] = x01;
    X[3] = x11;
  }
  Tenzor(const T& x00, const T& x10, const T& x20, const T& x01, const T& x11,
         const T& x21, const T& x02, const T& x12, const T& x22) {
    PInsist(D==3, "Number of arguments does not match Tenzor dimension!!");
    X[0] = x00;
    X[1] = x10;
    X[2] = x20;
    X[3] = x01;
    X[4] = x11;
    X[5] = x21;
    X[6] = x02;
    X[7] = x12;
    X[8] = x22;
  }

  // constructor from SymTenzor
  Tenzor(const SymTenzor<T,D>&);
kraus's avatar
kraus committed
85

gsell's avatar
gsell committed
86 87
  // constructor from AntiSymTenzor
  Tenzor(const AntiSymTenzor<T,D>&);
kraus's avatar
kraus committed
88

gsell's avatar
gsell committed
89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
  // destructor
  ~Tenzor() { };

  // assignment operators
  const Tenzor<T,D>& operator= (const Tenzor<T,D> &rhs) {
    TSV_MetaAssign< Tenzor<T,D> , Tenzor<T,D> ,OpAssign> :: apply(*this,rhs);
    return *this;
  }
  template<class T1>
  const Tenzor<T,D>& operator= (const Tenzor<T1,D> &rhs) {
    TSV_MetaAssign< Tenzor<T,D> , Tenzor<T1,D> ,OpAssign> :: apply(*this,rhs);
    return *this;
  }
  const Tenzor<T,D>& operator= (const T& rhs) {
    TSV_MetaAssignScalar< Tenzor<T,D> , T ,OpAssign> :: apply(*this,rhs);
    return *this;
  }

  // accumulation operators
  template<class T1>
  Tenzor<T,D>& operator+=(const Tenzor<T1,D> &rhs)
  {
kraus's avatar
kraus committed
111
    TSV_MetaAssign< Tenzor<T,D> , Tenzor<T1,D> , OpAddAssign > ::
gsell's avatar
gsell committed
112 113 114 115 116
      apply(*this,rhs);
    return *this;
  }
  Tenzor<T,D>& operator+=(const T& rhs)
  {
kraus's avatar
kraus committed
117
    TSV_MetaAssignScalar< Tenzor<T,D> , T , OpAddAssign > ::
gsell's avatar
gsell committed
118 119 120 121 122 123 124
      apply(*this,rhs);
    return *this;
  }

  template<class T1>
  Tenzor<T,D>& operator-=(const Tenzor<T1,D> &rhs)
  {
kraus's avatar
kraus committed
125
    TSV_MetaAssign< Tenzor<T,D> , Tenzor<T1,D> , OpSubtractAssign > ::
gsell's avatar
gsell committed
126 127 128 129 130
      apply(*this,rhs);
    return *this;
  }
  Tenzor<T,D>& operator-=(const T& rhs)
  {
kraus's avatar
kraus committed
131
    TSV_MetaAssignScalar< Tenzor<T,D> , T , OpSubtractAssign > ::
gsell's avatar
gsell committed
132 133 134 135 136 137 138
      apply(*this,rhs);
    return *this;
  }

  template<class T1>
  Tenzor<T,D>& operator*=(const Tenzor<T1,D> &rhs)
  {
kraus's avatar
kraus committed
139
    TSV_MetaAssign< Tenzor<T,D> , Tenzor<T1,D> , OpMultipplyAssign > ::
gsell's avatar
gsell committed
140 141 142 143 144
      apply(*this,rhs);
    return *this;
  }
  Tenzor<T,D>& operator*=(const T& rhs)
  {
kraus's avatar
kraus committed
145
    TSV_MetaAssignScalar< Tenzor<T,D> , T , OpMultipplyAssign > ::
gsell's avatar
gsell committed
146 147 148 149 150 151 152
      apply(*this,rhs);
    return *this;
  }

  template<class T1>
  Tenzor<T,D>& operator/=(const Tenzor<T1,D> &rhs)
  {
kraus's avatar
kraus committed
153
    TSV_MetaAssign< Tenzor<T,D> , Tenzor<T1,D> , OpDivideAssign > ::
gsell's avatar
gsell committed
154 155 156 157 158
      apply(*this,rhs);
    return *this;
  }
  Tenzor<T,D>& operator/=(const T& rhs)
  {
kraus's avatar
kraus committed
159
    TSV_MetaAssignScalar< Tenzor<T,D> , T , OpDivideAssign > ::
gsell's avatar
gsell committed
160 161 162 163 164 165 166
      apply(*this,rhs);
    return *this;
  }

  // Methods

  void diagonal(const T& rhs) {
kraus's avatar
kraus committed
167
    for (unsigned int i = 0 ; i < D ; i++ )
gsell's avatar
gsell committed
168 169 170 171 172 173 174 175 176 177
      (*this)(i,i) = rhs;
  }

  int len(void)  const { return Size; }
  int size(void) const { return sizeof(*this); }

  // Operators

  Element_t &operator[]( unsigned int i )
  {
178
    PAssert_LT(i, Size);
gsell's avatar
gsell committed
179 180 181 182 183
    return X[i];
  }

  Element_t operator[]( unsigned int i ) const
  {
184
    PAssert_LT(i, Size);
gsell's avatar
gsell committed
185 186 187 188 189 190
    return X[i];
  }

  //TJW: add these 12/16/97 to help with NegReflectAndZeroFace BC:
  // These are the same as operator[] but with () instead:

kraus's avatar
kraus committed
191
  Element_t& operator()(unsigned int i) {
gsell's avatar
gsell committed
192 193 194 195
    PAssert (i < Size);
    return X[i];
  }

kraus's avatar
kraus committed
196
  Element_t operator()(unsigned int i) const {
gsell's avatar
gsell committed
197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263
    PAssert (i < Size);
    return X[i];
  }
  //TJW.

  Element_t operator()( unsigned int i,  unsigned int j ) const {
    PAssert ( (i<D) && (j<D) );
    return X[i*D+j];
  }

  Element_t& operator()( unsigned int i, unsigned int j ) {
    PAssert ( (i<D) && (j<D) );
    return X[i*D+j];
  }

  Element_t operator()(const std::pair<int,int> i) const {
    PAssert ( (i.first>=0) && (i.second>=0) && (i.first<D) && (i.second<D) );
    return (*this)(i.first,i.second);
  }

  Element_t& operator()(const std::pair<int,int> i) {
    PAssert ( (i.first>=0) && (i.second>=0) && (i.first<D) && (i.second<D) );
    return (*this)(i.first,i.second);
  }


  //----------------------------------------------------------------------
  // Comparison operators.
  bool operator==(const Tenzor<T,D>& that) const {
    return TSV_MetaCompareArrays<T,T,D*D>::apply(X,that.X);
  }
  bool operator!=(const Tenzor<T,D>& that) const {
    return !(*this == that);
  }

  //----------------------------------------------------------------------
  // parallel communication
  Message& putMessage(Message& m) const {
    m.setCopy(true);
    const T *p = X;
    ::putMessage(m, p, p + D*D);
    return m;
  }

  Message& getMessage(Message& m) {
    T *p = X;
    ::getMessage(m, p, p + D*D);
    return m;
  }

private:

  // The elements themselves.
  T X[Size];

};


//////////////////////////////////////////////////////////////////////
//
// Free functions
//
//////////////////////////////////////////////////////////////////////

// trace() - sum of diagonal elements.

template <class T, unsigned D>
264
inline T trace(const Tenzor<T,D>& rhs) {
gsell's avatar
gsell committed
265
  T result = 0.0;
kraus's avatar
kraus committed
266
  for (unsigned int i = 0 ; i < D ; i++ )
gsell's avatar
gsell committed
267 268 269 270 271 272 273
    result += rhs(i,i);
  return result;
}

// transpose(). transpose(i,j) = input(j,i).

template <class T, unsigned D>
274
inline Tenzor<T,D> transpose(const Tenzor<T,D>& rhs) {
gsell's avatar
gsell committed
275 276
  Tenzor<T,D> result = typename Tenzor<T,D>::DontInitialize();

kraus's avatar
kraus committed
277 278
  for (unsigned int j = 0 ; j < D ; j++ )
    for (unsigned int i = 0 ; i < D ; i++ )
gsell's avatar
gsell committed
279 280 281 282 283 284 285
      result(i,j) = rhs(j,i);
  return result;
}

// Determinant: only implement for 1D, 2D, 3D:

template <class T, unsigned D>
286
inline T det(const Tenzor<T,D>& rhs) {
gsell's avatar
gsell committed
287 288 289 290 291
  PInsist(D<4, "Tenzor det() function not implemented for D>3!");
  return T(-999999.999999);
}

template <class T>
292
inline T det(const Tenzor<T,3>& rhs) {
gsell's avatar
gsell committed
293
  T result;
kraus's avatar
kraus committed
294
  result =
gsell's avatar
gsell committed
295 296 297 298 299 300 301
    rhs(0,0)*(rhs(1,1)*rhs(2,2) - rhs(1,2)*rhs(2,1)) +
    rhs(0,1)*(rhs(1,2)*rhs(2,0) - rhs(1,0)*rhs(2,2)) +
    rhs(0,2)*(rhs(1,0)*rhs(2,1) - rhs(1,1)*rhs(2,0));
  return result;
}

template <class T>
302
inline T det(const Tenzor<T,2>& rhs) {
gsell's avatar
gsell committed
303 304 305 306 307 308
  T result;
  result = rhs(0,0)*rhs(1,1) - rhs(0,1)*rhs(1,0);
  return result;
}

template <class T>
309
inline T det(const Tenzor<T,1>& rhs) {
gsell's avatar
gsell committed
310 311 312 313 314 315 316 317 318 319 320
  T result = rhs(0,0);
  return result;
}

// cofactors() - pow(-1, i+j)*M(i,j), where M(i,j) is a minor of the tensor.
// See, for example, Arfken, Mathematical Methods for Physicists, 2nd Edition,
// p. 157 (the section where the determinant of a matrix is defined).

// Only implement for 1D, 2D, 3D:

template <class T, unsigned D>
321
inline Tenzor<T,D> cofactors(const Tenzor<T,D>& rhs) {
gsell's avatar
gsell committed
322 323 324 325 326
  PInsist(D<4, "Tenzor cofactors() function not implemented for D>3!");
  return Tenzor<T,D>(-999999.999999);
}

template <class T>
327
inline Tenzor<T,3> cofactors(const Tenzor<T,3>& rhs) {
gsell's avatar
gsell committed
328 329 330 331 332 333 334 335 336 337 338 339 340 341 342
  Tenzor<T,3> result = typename Tenzor<T,3>::DontInitialize();

  result(0,0) = rhs(1,1)*rhs(2,2) - rhs(1,2)*rhs(2,1);
  result(1,0) = rhs(0,2)*rhs(2,1) - rhs(0,1)*rhs(2,2);
  result(2,0) = rhs(0,1)*rhs(1,2) - rhs(1,1)*rhs(0,2);
  result(0,1) = rhs(2,0)*rhs(1,2) - rhs(1,0)*rhs(2,2);
  result(1,1) = rhs(0,0)*rhs(2,2) - rhs(0,2)*rhs(2,0);
  result(2,1) = rhs(1,0)*rhs(0,2) - rhs(0,0)*rhs(1,2);
  result(0,2) = rhs(1,0)*rhs(2,1) - rhs(2,0)*rhs(1,1);
  result(1,2) = rhs(0,1)*rhs(2,0) - rhs(0,0)*rhs(2,1);
  result(2,2) = rhs(0,0)*rhs(1,1) - rhs(1,0)*rhs(0,1);
  return result;
}

template <class T>
343
inline Tenzor<T,2> cofactors(const Tenzor<T,2>& rhs) {
gsell's avatar
gsell committed
344 345 346 347 348 349 350 351 352 353 354 355
  Tenzor<T,2> result = typename Tenzor<T,2>::DontInitialize();

  result(0,0) =  rhs(1,1);
  result(1,0) = -rhs(0,1);
  result(0,1) = -rhs(1,0);
  result(1,1) =  rhs(0,0);
  return result;
}

// For D=1, cofactor is the unit tensor, because det = single tensor element
// value:
template <class T>
356
inline Tenzor<T,1> cofactors(const Tenzor<T,1>& rhs) {
gsell's avatar
gsell committed
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406
  Tenzor<T,1> result = Tenzor<T,1>(1);
  return result;
}


//////////////////////////////////////////////////////////////////////
//
// Unary Operators
//
//////////////////////////////////////////////////////////////////////

//----------------------------------------------------------------------
// unary operator-
template<class T, unsigned D>
inline Tenzor<T,D> operator-(const Tenzor<T,D> &op)
{
  return TSV_MetaUnary< Tenzor<T,D> , OpUnaryMinus > :: apply(op);
}

//----------------------------------------------------------------------
// unary operator+
template<class T, unsigned D>
inline const Tenzor<T,D> &operator+(const Tenzor<T,D> &op)
{
  return op;
}

//////////////////////////////////////////////////////////////////////
//
// Binary Operators
//
//////////////////////////////////////////////////////////////////////

//
// Elementwise operators.
//

TSV_ELEMENTWISE_OPERATOR(Tenzor,operator+,OpAdd)
TSV_ELEMENTWISE_OPERATOR(Tenzor,operator-,OpSubtract)
TSV_ELEMENTWISE_OPERATOR(Tenzor,operator*,OpMultipply)
TSV_ELEMENTWISE_OPERATOR(Tenzor,operator/,OpDivide)
TSV_ELEMENTWISE_OPERATOR(Tenzor,Min,FnMin)
TSV_ELEMENTWISE_OPERATOR(Tenzor,Max,FnMax)

//----------------------------------------------------------------------
// dot products.
//----------------------------------------------------------------------

template < class T1, class T2, unsigned D >
inline Tenzor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D>
kraus's avatar
kraus committed
407
dot(const Tenzor<T1,D> &lhs, const Tenzor<T2,D> &rhs)
gsell's avatar
gsell committed
408 409 410 411 412 413
{
  return TSV_MetaDot< Tenzor<T1,D> , Tenzor<T2,D> > :: apply(lhs,rhs);
}

template < class T1, class T2, unsigned D >
inline Vektor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D>
kraus's avatar
kraus committed
414
dot(const Vektor<T1,D> &lhs, const Tenzor<T2,D> &rhs)
gsell's avatar
gsell committed
415 416 417 418 419 420
{
  return TSV_MetaDot< Vektor<T1,D> , Tenzor<T2,D> > :: apply(lhs,rhs);
}

template < class T1, class T2, unsigned D >
inline Vektor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D>
kraus's avatar
kraus committed
421
dot(const Tenzor<T1,D> &lhs, const Vektor<T2,D> &rhs)
gsell's avatar
gsell committed
422 423 424 425 426 427 428 429 430 431
{
  return TSV_MetaDot< Tenzor<T1,D> , Vektor<T2,D> > :: apply(lhs,rhs);
}

//----------------------------------------------------------------------
// double dot product.
//----------------------------------------------------------------------

template < class T1, class T2, unsigned D >
inline typename PETEBinaryReturn<T1,T2,OpMultipply>::type
kraus's avatar
kraus committed
432
dotdot(const Tenzor<T1,D> &lhs, const Tenzor<T2,D> &rhs)
gsell's avatar
gsell committed
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456
{
  return TSV_MetaDotDot< Tenzor<T1,D> , Tenzor<T2,D> > :: apply(lhs,rhs);
}

//----------------------------------------------------------------------
// Outer product.
//----------------------------------------------------------------------

template<class T1, class T2, unsigned int D>
inline Tenzor<typename PETEBinaryReturn<T1,T2,OpMultipply>::type,D >
outerProduct(const Vektor<T1,D>& v1, const Vektor<T2,D>& v2)
{
  typedef typename PETEBinaryReturn<T1,T2,OpMultipply>::type T0;
  Tenzor<T0,D> ret = typename Tenzor<T0,D>::DontInitialize();

  for (unsigned int i=0; i<D; ++i)
    for (unsigned int j=0; j<D; ++j)
      ret(i,j) = v1[i]*v2[j];
  return ret;
}

//----------------------------------------------------------------------
// I/O
template<class T, unsigned D>
457
inline std::ostream& operator<<(std::ostream& out, const Tenzor<T,D>& rhs)
gsell's avatar
gsell committed
458 459
{
  if (D >= 1) {
kraus's avatar
kraus committed
460
    for (unsigned int i=0; i<D; i++) {
gsell's avatar
gsell committed
461
      out << "(";
kraus's avatar
kraus committed
462
      for (unsigned int j=0; j<D-1; j++) {
gsell's avatar
gsell committed
463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480
	out << rhs(i,j) << " , ";
      }
      out << rhs(i,D-1) << ")";
      // I removed this. --TJW: if (i < D - 1) out << endl;
    }
  } else {
    out << "( " << rhs(0,0) << " )";
  }
  return out;
}

// include header files for SymTenzor and AntiSymTenzor in order
// to define constructors for Tenzor using these types
#include "AppTypes/SymTenzor.h"
#include "AppTypes/AntiSymTenzor.h"

template <class T, unsigned D>
Tenzor<T,D>::Tenzor(const SymTenzor<T,D>& rhs) {
kraus's avatar
kraus committed
481 482
  for (unsigned int i=0; i<D; ++i)
    for (unsigned int j=0; j<D; ++j)
gsell's avatar
gsell committed
483 484 485 486 487
      (*this)(i,j) = rhs(i,j);
}

template <class T, unsigned D>
Tenzor<T,D>::Tenzor(const AntiSymTenzor<T,D>& rhs) {
kraus's avatar
kraus committed
488 489
  for (unsigned int i=0; i<D; ++i)
    for (unsigned int j=0; j<D; ++j)
gsell's avatar
gsell committed
490 491 492 493 494 495 496 497 498
      (*this)(i,j) = rhs(i,j);
}


#endif // TENZOR_H

/***************************************************************************
 * $RCSfile: Tenzor.h,v $   $Author: adelmann $
 * $Revision: 1.1.1.1 $   $Date: 2003/01/23 07:40:24 $
kraus's avatar
kraus committed
499
 * IPPL_VERSION_ID: $Id: Tenzor.h,v 1.1.1.1 2003/01/23 07:40:24 adelmann Exp $
gsell's avatar
gsell committed
500
 ***************************************************************************/