Tensor.h

#ifndef UTIL_TENSOR_H
#define UTIL_TENSOR_H

/*
* Util Package - C++ Utilities for Scientific Computation
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/


#include <util/space/Dimension.h>
#include <util/global.h>

#ifdef UTIL_MPI
#include <util/mpi/MpiTraits.h>
#endif

#include <iostream>
#include <cmath>

namespace Util
{

   class Vector;

   class Tensor
   {

   public:


      
      Tensor();

      Tensor(const Tensor& t);

      explicit Tensor(double scalar);

      explicit Tensor(const double a[][Dimension]);

      

      
      Tensor& operator=(const Tensor& t);

      Tensor& operator=(const double a[][Dimension]);


      
      void operator+=(const Tensor& dt);

      void operator-=(const Tensor& dt);

      void operator*=(double s);

      void operator/=(double s);


      
      const double& operator ()(int i, int j) const;

      double& operator () (int i, int j);


      
      Tensor& zero();

      Tensor& identity();

      Tensor& setRow(int i, const Vector& r);

      Tensor& setColumn(int i, const Vector& r);

      Tensor& add(const Tensor& t1, const Tensor& t2);

      Tensor& subtract(const Tensor& t1, const Tensor& t2);

      Tensor& multiply(const Tensor& t, double s);

      Tensor& divide(const Tensor& t, double s);

      Tensor& transpose(const Tensor& t);

      Tensor& transpose();

      Tensor& symmetrize(const Tensor& t);

      Tensor& symmetrize();

      Tensor& dyad(const Vector& v1, const Vector& v2);



      double trace() const; 

      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);


     
      static const Tensor Zero;

      static const Tensor Identity;

      static void initStatic();

      #ifdef UTIL_MPI

      static void commitMpiType();
      #endif


   private:

      static const int Width = 13;

      static const int Precision = 5;

      double elem_[DimensionSq];

   //friends:

      friend bool operator == (const Tensor& t1, const Tensor& t2);

      friend bool operator == (const Tensor& t1, const double t2[][Dimension]);

      friend std::istream& operator >> (std::istream& in, Tensor &tensor);

      friend std::ostream& operator << (std::ostream& out, const Tensor &tensor);

   };

   bool operator == (const Tensor& t1, const Tensor& t2);

   bool operator == (const Tensor& t1, const double t2[][Dimension]);

   bool operator == (const double t1[][Dimension], const Tensor& t2);

   bool operator != (const Tensor& t1, const Tensor& t2);

   bool operator != (const Tensor& t1, const double t2[][Dimension]);

   bool operator != (const double t1[][Dimension], const Tensor& t2);

   std::istream& operator >> (std::istream& in, Tensor &tensor);

   std::ostream& operator << (std::ostream& out, const Tensor &tensor);

   #ifdef UTIL_MPI

   template <>
   class MpiTraits<Tensor>
   {  
   public:  
      static MPI::Datatype type;   
      static bool hasType;         
   };
   #endif

   // Inline methods

   /*
   * Default constructor
   */
   inline
   Tensor::Tensor()
   {}

   /*
   * Copy constructor
   */
   inline
   Tensor::Tensor(const Tensor& t)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] = t.elem_[i];
      }
   }

   /*
   * Constructor, initialize all elements to a scalar value.
   */
   inline
   Tensor::Tensor(double scalar)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] = scalar;
      }
   }

   /*
   * Construct Tensor from C double [][Dimension] array.
   */
   inline
   Tensor::Tensor(const double a[][Dimension])
   {
      for (int i = 0; i < Dimension; ++i) {
         for (int j = 0; j < Dimension; ++j) {
            elem_[i*Dimension + j] = a[i][j];
         }
      }
   }

   /*
   * Set all elements of this tensor to zero.
   */
   inline
   Tensor& Tensor::zero()
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] = 0.0;
      }
      return *this;
   }

   /*
   * Set this to the identity (unity) tensor.
   */
   inline
   Tensor& Tensor::identity()
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] = 0.0;
      }
      for (int i = 0; i < Dimension; ++i) {
         elem_[i*Dimension + i] = 1.0;
      }
      return *this;
   }

   /*
   * Assignment.
   */
   inline
   Tensor& Tensor::operator=(const Tensor& t)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] = t.elem_[i];
      }
      return *this;
   }

   /*
   * Assignment from C double [][Dimension] array.
   */
   inline
   Tensor& Tensor::operator=(const double a[][Dimension])
   {
      int i, j;
      for (i = 0; i < Dimension; ++i) {
         for (j = 0; j < Dimension; ++j) {
            elem_[i*Dimension + j] = a[i][j];
         }
      }
      return *this;
   }

   /*
   * Add tensor dt to this tensor.
   */
   inline
   void Tensor::operator+=(const Tensor& dt)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] += dt.elem_[i];
      }
   }

   /*
   * Subtract tensor dt from this tensor.
   */
   inline
   void Tensor::operator-=(const Tensor& dt)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] -= dt.elem_[i]; 
      }
   }

   /*
   * Multiply this tensor by scalar s.
   */
   inline
   void Tensor::operator*=(double s)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] *= s; 
      }
   }

   /*
   * Divide this tensor by scalar s.
   */
   inline
   void Tensor::operator/=(double s)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] /= s; 
      }
   }

   /*
   * Return one element by value.
   */
   inline
   const double& Tensor::operator()(int i, int j) const
   {
      assert(i >=0);  
      assert(i < Dimension);  
      assert(j >=0);  
      assert(j < Dimension);  
      return elem_[i*Dimension + j]; 
   }

   /*
   * Return a reference to one element of the tensor.
   */
   inline
   double& Tensor::operator()(int i, int j)
   {  
      assert(i >=0);  
      assert(i < Dimension);  
      assert(j >=0);  
      assert(j < Dimension);  
      return elem_[i*Dimension + j]; 
   }

   /*
   * Add tensors t1 and t2.
   *
   * Upon return, *this = t1 + t2.
   */
   inline
   Tensor& Tensor::add(const Tensor& t1, const Tensor& t2)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] = t1.elem_[i] + t2.elem_[i];
      }
      return *this;
   }

   /*
   * Subtract tensor t2 from t1.
   *
   * Upon return, *this = t1 - t2.
   */
   inline
   Tensor& Tensor::subtract(const Tensor& t1, const Tensor& t2)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] = t1.elem_[i] - t2.elem_[i];
      }
      return *this;
   }

   /*
   * Multiply a tensor t by a scalar s.
   *
   * Upon return, *this = t*s.
   */
   inline
   Tensor& Tensor::multiply(const Tensor& t, double s)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] = t.elem_[i]*s;
      }
      return *this;
   }

   /*
   * Divide tensor t by scalar s.
   *
   * Upon return, *this = t/s;
   */
   inline
   Tensor& Tensor::divide(const Tensor& t, double s)
   {
      for (int i = 0; i < DimensionSq; ++i) {
         elem_[i] = t.elem_[i]/s;
      }
      return *this;
   }

   /*
   * Return trace of a tensor.
   */
   inline
   double Tensor::trace() const
   {
      double trace = 0.0;
      for (int i = 0; i < Dimension; ++i) {
         trace += (*this)(i, i);
      }
      return trace;
   }

   /*
   * Compute the transpose of a tensor.
   *
   * Upon return, this tensor is the transpose of t.
   */
   inline
   Tensor& Tensor::transpose(const Tensor& t)
   {
      int i, j;
      for (i = 0; i < Dimension; ++i) {
         (*this)(i, i) = t(i, i);
      }
      for (i = 1; i < Dimension; ++i) {
         for (j = 0; j < i; ++j) {
            (*this)(i, j) = t(j, i);  
            (*this)(j, i) = t(i, j);
         }
      }
      return *this;
   }

   /*
   * Transpose this tensor.
   *
   * Upon return, this tensor is transposed.
   */
   inline
   Tensor& Tensor::transpose()
   {
      double save;
      int i, j;
      for (i = 1; i < Dimension; ++i) {
         for (j = 0; j < i; ++j) {
            save = (*this)(i, j);
            (*this)(i, j) = (*this)(j, i);  
            (*this)(j, i) = save;
         }
      }
      return *this;
   }

   /*
   * Compute symmetric part of a tensor.
   *
   * Upon return, this is symmetric part of t:
   * *this = 0.5*(t + t.transpose());
   */
   inline
   Tensor& Tensor::symmetrize(const Tensor& t)
   {
      double ave;
      int i, j;
      for (i = 0; i < Dimension; ++i) {
         (*this)(i, i) = t(i, i);
      }
      for (i = 1; i < Dimension; ++i) {
         for (j = 0; j < i; ++j) {
            ave = 0.5*( t(i, j) + t(j, i) );
            (*this)(i, j) = ave;  
            (*this)(j, i) = ave;
         }
      }
      return *this;
   }

   /*
   * Symmetrize this tensor.
   *
   * Upon return, t is symmetrized:
   */
   inline
   Tensor& Tensor::symmetrize()
   {
      double ave;
      int i, j;
      for (i = 1; i < Dimension; ++i) {
         for (j = 0; j < i; ++j) {
            ave = 0.5*( (*this)(i, j) + (*this)(j, i) );
            (*this)(i, j) = ave;  
            (*this)(j, i) = ave;
         }
      }
      return *this;
   }

   /*
   * Serialize to/from an archive.
   */
   template <class Archive>
   inline 
   void Tensor::serialize(Archive& ar, const unsigned int version)
   { 
      for (int i = 0; i < DimensionSq; ++i) {
         ar & elem_[i];
      }
   }

}
 
#include <util/space/Vector.h>
namespace Util
{

   /*
   * Set row i of the tensor to elements of vector r.
   */
   inline
   Tensor& Tensor::setRow(int i, const Vector& r)
   {
      for (int j = 0; j < Dimension; j++) {
         elem_[i*Dimension+j] = r[j];
      }
      return *this;
   }

   /*
   * Set column i of the tensor to elements of vector r.
   */
   inline
   Tensor& Tensor::setColumn(int j, const Vector& r)
   {
      for (int i = 0; i < Dimension; i++) {
         elem_[i*Dimension+j] = r[i];
      }
      return *this;
   }

   /*
   * Multiply two vectors to create a dyadic tensor.
   *
   * Upon return, *this = t*s.
   */
   inline
   Tensor& Tensor::dyad(const Vector& v1, const Vector& v2)
   {
      int i, j;
      for (i = 0; i < Dimension; ++i) {
         for (j = 0; j < Dimension; ++j) {
            elem_[i*Dimension + j] = v1[i]*v2[j];
         }
      }
      return *this;
   }

} 
#endif