RealTensor2.H

Go to the documentation of this file.

 
#ifndef _REALTENSOR2_H_
#define _REALTENSOR2_H_
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <vector>
#include "Types.H"
#include "Vec3.H"
#include "gsl/gsl_eigen.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_vector.h"
 
namespace criptic {
 
  class RealTensor2 {
 
  public:
 
    // Storage
    
    Real t[9];            
    RealTensor2() {
      t[0] = t[1] = t[2] = t[3] = t[4] = t[5] = t[6] = t[7] = t[8] = 0.0;
    }
    
    RealTensor2(const Real& xx, const Real& xy, const Real& xz,
        const Real& yx, const Real& yy, const Real& yz,
        const Real& zx, const Real& zy, const Real& zz)
    {
      t[0] = xx; t[1] = xy; t[2] = xz;
      t[3] = yx; t[4] = yy; t[5] = yz;
      t[6] = zx; t[7] = zy; t[8] = zz;
    }
 
    RealTensor2(const RealVec &a, const RealVec &b) {
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++)
      t[3*i+j] = a[i] * b[j];
    }
 
    // Element access
    
    Real& operator () (int i, int j) {
      assert(i >= 0 && i < 3 && j >= 0 && j < 3);
      return t[3*i+j];
    }
    const Real& operator () (int i, int j) const {
      assert(i >= 0 && i < 3 && j >= 0 && j < 3);
      return t[3*i+j];
    }
    Real& operator [] (int i) {
      assert(i >= 0 && i < 9);
      return t[i];
    }
    const Real& operator [] (int i) const {
      assert(i >= 0 && i < 9);
      return t[i];
    }
    
    // Unary operators
 
    RealTensor2 operator-() const {
      RealTensor2 rt( -t[0], -t[1], -t[2],
              -t[3], -t[4], -t[5],
              -t[6], -t[7], -t[8] );
      return rt;
    }
    
    // Unary operators with scalars
 
    RealTensor2& operator+=(const Real a) {
      for (int i=0; i<9; i++) t[i] += a;
      return *this;
    }
    RealTensor2& operator-=(const Real a) {
      for (int i=0; i<9; i++) t[i] -= a;
      return *this;
    }
    RealTensor2& operator*=(const Real a) {
      for (int i=0; i<9; i++) t[i] *= a;
      return *this;
    }
    RealTensor2& operator/=(const Real a) {
      for (int i=0; i<9; i++) t[i] /= a;
      return *this;
    }
 
    // Unary operators with tensors
    
    RealTensor2& operator+=(const RealTensor2& a) {
      for (int i=0; i<9; i++) t[i] += a[i];
      return *this;
    }
    RealTensor2& operator-=(const RealTensor2& a) {
      for (int i=0; i<9; i++) t[i] -= a[i];
      return *this;
    }
    RealTensor2& operator*=(const RealTensor2& a) {
      for (int i=0; i<9; i++) t[i] *= a[i];
      return *this;
    }
    RealTensor2& operator/=(const RealTensor2& a) {
      for (int i=0; i<9; i++) t[i] /= a[i];
      return *this;
    }
    
    // Binary operations with scalars
 
    RealTensor2& operator=(const Real a) {
      for (int i=0; i<9; i++) t[i] = a;
      return *this;
    }
    
    RealTensor2 operator+(const Real a) const {
      RealTensor2 rt(*this);
      for (int i=0; i<9; i++) rt.t[i] += a;
      return rt;
    }
    RealTensor2 operator-(const Real a) const {
      RealTensor2 rt(*this);
      for (int i=0; i<9; i++) rt.t[i] -= a;
      return rt;
    }
    RealTensor2 operator*(const Real a) const {
      RealTensor2 rt(*this);
      for (int i=0; i<9; i++) rt.t[i] *= a;
      return rt;
    }
    RealTensor2 operator/(const Real a) const {
      RealTensor2 rt(*this);
      for (int i=0; i<9; i++) rt.t[i] /= a;
      return rt;
    }
 
    // Matrix-vector multiplication
 
    RealVec operator*(const RealVec& v) const {
      RealVec r;
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) r[i] += (*this)(i,j) * v[j];
      return r;
    }
    
    // Binary operations with tensors
    
    RealTensor2 operator+(const RealTensor2& a) const {
      RealTensor2 rt(*this);
      for (int i=0; i<9; i++) rt.t[i] += a[i];
      return rt;
    }
    RealTensor2 operator-(const RealTensor2& a) const {
      RealTensor2 rt(*this);
      for (int i=0; i<9; i++) rt.t[i] -= a[i];
      return rt;
    }
    RealTensor2 operator*(const RealTensor2& a) const {
      RealTensor2 rt(*this);
      for (int i=0; i<9; i++) rt.t[i] *= a[i];
      return rt;
    }
    RealTensor2 operator/(const RealTensor2& a) const {
      RealTensor2 rt(*this);
      for (int i=0; i<9; i++) rt.t[i] /= a[i];
      return rt;
    }
    
    // Binary operations with vectors and tensors
 
    RealVec contract1(const RealVec& v) const {
      RealVec rv;
      for (int j=0; j<3; j++)
    for (int i=0; i<3; i++) rv(j) += v(i) * (*this)(i,j);
      return rv;
    }
 
    RealVec contract2(const RealVec& v) const {
      return (*this) * v;
    }
 
    Real contract(const RealVec& a, const RealVec& b) const {
      Real c = 0.0;
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++)
      c += a[i] * b[j] * (*this)(i,j);
      return c;
    }
 
    Real contract(const RealTensor2& a) const {
      Real c = 0.0;
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) c += (*this)(i,j) * a(i,j);
      return c;
    }
 
    Real divTangent(const RealVec& v) const {
 
      // This routine makes use of the relation
      // d/dx_i eT_i = d/dx_i (v_i / |v|)
      //    = (1 / |v|) [ (dv_i / dx_i) - (v_i v_j dv_j / dx_i) / |v| ]
      Real divT = 0.0;
      Real vmag = v.mag();
      for (int i=0; i<3; i++) {
    for (int j=0; j<3; j++) {
      divT += -v[i] * v[j] * (*this)(j,i) / vmag;
    }
    divT += (*this)(i,i);
      }
      divT /= vmag;
      return divT;
    }
    
    // Unary operations
 
    Real trace() const {
      Real tr = 0.0;
      for (int i=0; i<3; i++) tr += (*this)(i,i);
      return tr;
    }
 
    RealTensor2 T() const {
      RealTensor2 rt;
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) rt(i,j) = (*this)(j,i);
      return rt;
    }
 
    void transpose() {
      for (int i=1; i<3; i++)
    for (int j=0; j<i; j++)
      std::swap( (*this)(i,j), (*this)(j,i) );
    }
 
    Real mag() {
      Real m2 = 0.0;
      for (int i=0; i<9; i++) m2 += t[i]*t[i];
      return std::sqrt(m2);
    }
 
    Real det() const {
      return (*this)(0,0) * ( (*this)(1,1) * (*this)(2,2) -
                  (*this)(1,2) * (*this)(2,1) ) +
    (*this)(0,1) * ( (*this)(1,2) * (*this)(2,0) -
             (*this)(1,0) * (*this)(2,2) ) +
    (*this)(0,2) * ( (*this)(1,0) * (*this)(2,1) -
             (*this)(1,1) * (*this)(2,0) );
    }
 
    RealVec eVal() const {
 
      // Allocate workspace
      gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(3);
      gsl_matrix *T = gsl_matrix_alloc(3,3);
      gsl_matrix *Q = gsl_matrix_alloc(3,3);
      gsl_vector *D = gsl_vector_alloc(3);
      
      // Create a GSL matrix corresponding to our tensor      
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) gsl_matrix_set(T, i, j, (*this)(i,j));
      
      // Get eigenvalues and eigenvectors
      gsl_eigen_symmv(T, D, Q, w);
 
      // Pack the results into our vector
      RealVec eV;
      for (int i=0; i<3; i++) {
    eV[i] = gsl_vector_get(D, i);
      }
 
      // Free workspace
      gsl_matrix_free(T);
      gsl_matrix_free(Q);
      gsl_vector_free(D);
      gsl_eigen_symmv_free(w);
 
      // Return
      return eV;
    }
    
    RealTensor2 eVec() const {
 
      // Allocate workspace
      gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(3);
      gsl_matrix *T = gsl_matrix_alloc(3,3);
      gsl_matrix *Q = gsl_matrix_alloc(3,3);
      gsl_vector *D = gsl_vector_alloc(3);
      
      // Create a GSL matrix corresponding to our tensor      
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) gsl_matrix_set(T, i, j, (*this)(i,j));
      
      // Get eigenvalues and eigenvectors
      gsl_eigen_symmv(T, D, Q, w);
 
      // Pack the result into a new tensor
      RealTensor2 ev;
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) ev(i,j) = gsl_matrix_get(Q, i, j);
 
      // Free workspace
      gsl_matrix_free(T);
      gsl_matrix_free(Q);
      gsl_vector_free(D);
      gsl_eigen_symmv_free(w);
 
      // Return
      return ev;
    }
      
    void eigenDecompose(RealVec& eVal, RealTensor2& eVec) const {
 
      // Allocate workspace
      gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(3);
      gsl_matrix *T = gsl_matrix_alloc(3,3);
      gsl_matrix *Q = gsl_matrix_alloc(3,3);
      gsl_vector *D = gsl_vector_alloc(3);
      
      // Create a GSL matrix corresponding to our tensor      
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) gsl_matrix_set(T, i, j, (*this)(i,j));
      
      // Get eigenvalues and eigenvectors
      gsl_eigen_symmv(T, D, Q, w);
 
      // Pack the results into our vector and tensor structures
      for (int i=0; i<3; i++) {
    eVal[i] = gsl_vector_get(D, i);
    for (int j=0; j<3; j++) eVec(i,j) = gsl_matrix_get(Q, i, j);
      }
 
      // Free workspace
      gsl_matrix_free(T);
      gsl_matrix_free(Q);
      gsl_vector_free(D);
      gsl_eigen_symmv_free(w);
    }
    
    RealTensor2 inv() const {
      // This method works by finding the eigendecomposition of the
      // matrix, as T = Q D Q^T, where Q is an orthogonal matrix whose
      // columns are the eigenvalues of T, and D is a diagonal matrix
      // whose diagonal elements are the corresponding
      // eigenvalues. The inverse is then trivial to compute: since
      // I = T T^-1 = (Q D Q^T) (Q D^-1 Q^T), it immeidately follows
      // that T^-1 = Q D^-1 Q^T.
 
      // Allocate workspace
      gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(3);
      gsl_matrix *T = gsl_matrix_alloc(3,3);
      gsl_matrix *Q = gsl_matrix_alloc(3,3);
      gsl_vector *D = gsl_vector_alloc(3);
      
      // Create a GSL matrix corresponding to our tensor      
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) gsl_matrix_set(T, i, j, (*this)(i,j));
      
      // Get eigenvalues and eigenvectors
      gsl_eigen_symmv(T, D, Q, w);
      
      // Get D^-1
      for (int i=0; i<3; i++)
    gsl_vector_set(D, i, 1/(gsl_vector_get(D, i)));
 
      // Compute final tensor; note that we take advantage of the fact
      // that D^(-1) is diagonal to simplify as follows:
      // (Q D^(-1) Q^T)_ij = Q_in (D^(-1))_nm Q^T_mj
      //                    = Q_in Q_jm (D^(-1))_nm
      //                    = Q_in Q_jm (D^(-1))_nm delta_nm
      //                    = Q_in Q_jn (D^(-1))_nn
      // Also note that, since the final tensor is also symmetric, we
      // can just compute the diagonal and lower triangular elements,
      // then fill in the rest by symmetry.
      RealTensor2 Tinv;
      for (int i=0; i<3; i++) {
    for (int j=0; j<=i; j++) {
      for (int n=0; n<3; n++) {
        Tinv(i,j) += gsl_vector_get(D, n) *
          gsl_matrix_get(Q, i, n) *
          gsl_matrix_get(Q, j, n);
      }
    }
      }
      for (int i=0; i<3; i++) {
    for (int j=i+1; j<3; j++) {
      Tinv(i,j) = Tinv(j,i);
    }
      }
 
      // Free workspace
      gsl_matrix_free(T);
      gsl_matrix_free(Q);
      gsl_vector_free(D);
      gsl_eigen_symmv_free(w);
 
      // Return
      return Tinv;      
    }
    
    RealTensor2 sqrt() const {
      // Note: this method computes the square root by finding the
      // eigenvalue and eigenvectors of the tensor; for a real,
      // symmetric matrix T, we can uniquely decomposite it as T = Q D
      // Q^T, where Q is an orthogonal matrix whose columns are the
      // eigenvalues of T, and D is a diagonal matrix whose diagonal
      // elements are the corresponding eigenvalues. Then it
      // immediately follows that T = (Q D^(1/2) Q^T) (Q D^(1/2) Q^T),
      // since Q^T = Q^-1 for an orthogonal matrix, and therefore that
      // Q D^(1/2) Q^T = T^(1/2).
      
      // Allocate workspace
      gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(3);
      gsl_matrix *T = gsl_matrix_alloc(3,3);
      gsl_matrix *Q = gsl_matrix_alloc(3,3);
      gsl_vector *D = gsl_vector_alloc(3);
      
      // Create a GSL matrix corresponding to our tensor      
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) gsl_matrix_set(T, i, j, (*this)(i,j));
      
      // Get eigenvalues and eigenvectors
      gsl_eigen_symmv(T, D, Q, w);
      
      // Take square root of the eigenvalues to form D^(1/2)
      for (int i=0; i<3; i++)
    gsl_vector_set(D, i, std::sqrt(gsl_vector_get(D, i)));
 
      // Compute final tensor; note that we take advantage of the fact
      // that D^(1/2) is diagonal to simplify as follows:
      // (Q D^(1/2) Q^T)_ij = Q_in (D^(1/2))_nm Q^T_mj
      //                    = Q_in Q_jm (D^(1/2))_nm
      //                    = Q_in Q_jm (D^(1/2))_nm delta_nm
      //                    = Q_in Q_jn (D^(1/2))_nn
      // Also note that, since the final tensor is also symmetric, we
      // can just compute the diagonal and lower triangular elements,
      // then fill in the rest by symmetry.
      RealTensor2 Tsqrt;
      for (int i=0; i<3; i++) {
    for (int j=0; j<3; j++) {
      for (int n=0; n<3; n++) {
        Tsqrt(i,j) += gsl_vector_get(D, n) *
          gsl_matrix_get(Q, i, n) *
          gsl_matrix_get(Q, j, n);
      }
    }
      }
      for (int i=0; i<3; i++) {
    for (int j=i+1; j<3; j++) {
      Tsqrt(i,j) = Tsqrt(j,i);
    }
      }
 
      // Free workspace
      gsl_matrix_free(T);
      gsl_matrix_free(Q);
      gsl_vector_free(D);
      gsl_eigen_symmv_free(w);
 
      // Return
      return Tsqrt;
    }
 
    RealTensor2 invSqrt() const {
      // Note: this method computes the inverse of the square root by
      // eigendecomposition, using a slight modification of the
      // algorithm used to compute the square root. Given that T = Q D
      // Q^T is the eigendecomposition of T, and further as shown
      // above that T^(1/2) = Q D^(1/2) Q^T, it further follows that I
      // = (Q D^(1/2) Q^T) (Q D^(-1/2) Q^T) = T^(1/2) (Q D^(-1/2) Q^T),
      // and therefore that (Q D^(-1/2) Q^T) = (T^(1/2))^(-1).
      
      // Allocate workspace
      gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(3);
      gsl_matrix *T = gsl_matrix_alloc(3,3);
      gsl_matrix *Q = gsl_matrix_alloc(3,3);
      gsl_vector *D = gsl_vector_alloc(3);
      
      // Create a GSL matrix corresponding to our tensor      
      for (int i=0; i<3; i++)
    for (int j=0; j<3; j++) gsl_matrix_set(T, i, j, (*this)(i,j));
      
      // Get eigenvalues and eigenvectors
      gsl_eigen_symmv(T, D, Q, w);
      
      // Compute D^(-1/2)
      for (int i=0; i<3; i++)
    gsl_vector_set(D, i, 1/std::sqrt(gsl_vector_get(D, i)));
 
      // Compute final tensor; note that we take advantage of the fact
      // that D^(1/2) is diagonal to simplify as follows:
      // (Q D^(1/2) Q^T)_ij = Q_in (D^(-1/2))_nm Q^T_mj
      //                    = Q_in Q_jm (D^(-1/2))_nm
      //                    = Q_in Q_jm (D^(-1/2))_nm delta_nm
      //                    = Q_in Q_jn (D^(-1/2))_nn
      // Also note that, since the final tensor is also symmetric, we
      // can just compute the diagonal and lower triangular elements,
      // then fill in the rest by symmetry.
      RealTensor2 TInvSqrt;
      for (int i=0; i<3; i++) {
    for (int j=0; j<=i; j++) {
      for (int n=0; n<3; n++) {
        TInvSqrt(i,j) += gsl_vector_get(D, n) *
          gsl_matrix_get(Q, i, n) *
          gsl_matrix_get(Q, j, n);
      }
    }
      }
      for (int i=0; i<3; i++) {
    for (int j=i+1; j<3; j++) {
      TInvSqrt(i,j) = TInvSqrt(j,i);
    }
      }
 
      // Free workspace
      gsl_matrix_free(T);
      gsl_matrix_free(Q);
      gsl_vector_free(D);
      gsl_eigen_symmv_free(w);
 
      // Return
      return TInvSqrt;
    }
    
  };
 
  // Define the zero and identity tensors
  static const criptic::RealTensor2
  zeroTensor(0,0,0,0,0,0,0,0,0); 
  static const criptic::RealTensor2
  identityTensor(1,0,0,0,1,0,0,0,1); 
  // Overload the basic arithmetic operators for Real numbers to allow
  // them to operate on RealTensor2's, so for example 0.5 * t works as
  // well as t * 0.5.
  inline criptic::RealTensor2 operator+(const criptic::Real a,
                    const criptic::RealTensor2& t) {
    RealTensor2 rt( a + t[0], a + t[1], a + t[2],
            a + t[3], a + t[4], a + t[5],
            a + t[6], a + t[7], a + t[8] );
    return rt;
  }
  inline criptic::RealTensor2 operator-(const criptic::Real a,
                    const criptic::RealTensor2& t) {
    RealTensor2 rt( a - t[0], a - t[1], a - t[2],
            a - t[3], a - t[4], a - t[5],
            a - t[6], a - t[7], a - t[8] );
    return rt;
  }
  inline criptic::RealTensor2 operator*(const criptic::Real a,
                    const criptic::RealTensor2& t) {
    RealTensor2 rt( a * t[0], a * t[1], a * t[2],
            a * t[3], a * t[4], a * t[5],
            a * t[6], a * t[7], a * t[8] );
    return rt;
  }
  inline criptic::RealTensor2 operator/(const criptic::Real a,
                    const criptic::RealTensor2& t) {
    RealTensor2 rt( a / t[0], a / t[1], a / t[2],
            a / t[3], a / t[4], a / t[5],
            a / t[6], a / t[7], a / t[8] );
    return rt;
  }
 
  inline criptic::RealTensor2 outer(const criptic::RealVec& v1,
                    const criptic::RealVec& v2) {
    RealTensor2 rt(v1, v2);
    return rt;
  }
 
  typedef struct {
    RealVec eT;          
    RealVec eN;          
    RealVec eB;          
  } TNBBasis;
 
  inline TNBBasis TNBVectors(const RealVec& v,
                 const RealTensor2& vGrad) {
 
    // Basis object to return
    TNBBasis b;
      
    // First set tangent vector
    Real vmag = v.mag();
    b.eT = v/vmag;
 
    // Next set normal vector, using the relation
    // eN_i \propto eT_j (d eT_i / dx_j)
    //      \propto eT_j ( |v| dv_i / dx_j - v_i v_k dv_k / dx_j )
    RealVec vDotGradv;   // (vDotGradv)_j = v_k dv_k / dx_j
    for (int j=0; j<3; j++) {
      for (int k=0; k<3; k++) vDotGradv[j] += v[k] * vGrad[3*k + j];
    }
    Real norm = 0.0;
    for (int i=0; i<3; i++) {
      b.eN[i] = 0.0;
      for (int j=0; j<3; j++)
    b.eN[i] += b.eT[j] * (vmag*vmag * vGrad[3*i+j] - v[i] * vDotGradv[j]);
      norm += b.eN[i] * b.eN[i];
    }
 
    // Handle special case where field is straight so grad(v) = 0, in
    // which case eN is degenerate, and the calculation we have just
    // done will yield a norm of zero
    if (norm > 0.0) {
      // Non-degerate case
      b.eN.normalize();
    } else {
      // Degenerate case; break the degeneracy by choosing eN to be in
      // the xy plane, unless eT = \hat{z}, in which case we just
      // choose eN = \hat{x}
      Real xynorm = std::sqrt( b.eT[0]*b.eT[0] + b.eT[1]*b.eT[1] );
      if (xynorm == 0.0) {
    b.eN[0] = 1.0;
    b.eN[1] = b.eN[2] = 0.0;
      } else {
    b.eN[0] = b.eT[1] / xynorm;
    b.eN[1] = -b.eT[0] / xynorm;
      }
    }
    
    // Get binormal vector
    b.eB = b.eT.cross(b.eN);
 
    // Return
    return b;
  }
}
 
#endif
// _REALTENSOR2_H_