Tensor.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _TENSOR_H
#define _TENSOR_H
 
#include <vector>
#include <iostream>
 
class Vec3;
 
 
class Tensor
{
protected:
  using t_ind = unsigned short int; 
 
  const t_ind       n; 
  std::vector<Real> v; 
 
  virtual t_ind index(t_ind i, t_ind j) const { return i-1 + n*(j-1); }
 
private:
  void define3Dtransform(const Vec3& v1, const Vec3& v2, const Vec3& v3);
 
public:
  explicit Tensor(const t_ind nsd, bool identity = false);
  explicit Tensor(const Vec3& vn, bool vnIsX = false);
  Tensor(const Vec3& t1, const Vec3& t2,
         bool t1isZ = false, bool t2isXZ = false);
  Tensor(const Vec3& v1, const Vec3& v2, const Vec3& v3);
  Tensor(Real a1, Real a2, Real a3);
  Tensor(const Tensor& T, bool transpose = false);
  Tensor(const std::vector<Real>& a, bool transpose = false);
 
  virtual ~Tensor() {}
 
  void zero() { std::fill(v.begin(),v.end(),Real(0)); }
 
  void diag(Real value = Real(1));
  void diag(const Vec3& diagonal);
 
  operator const std::vector<Real>&() const { return v; }
  operator std::vector<Real>&() { return v; }
 
  const Real* ptr() const { return &v.front(); }
 
  const Real& operator()(t_ind i, t_ind j) const { return v[this->index(i,j)]; }
  Real& operator()(t_ind i, t_ind j) { return v[this->index(i,j)]; }
  Vec3 operator[](t_ind i) const;
 
  Tensor& operator=(const Tensor& T);
  Tensor& operator=(const std::vector<Real>& val);
  Tensor& operator=(Real val);
 
  Tensor& operator+=(const Tensor& T);
  Tensor& operator+=(Real val);
 
  Tensor& operator-=(const Tensor& T);
  Tensor& operator-=(Real val);
 
  Tensor& operator*=(Real val);
  Tensor& operator*=(const Tensor& B) { return this->postMult(B); }
  Tensor& postMult(const Tensor& B);
  Tensor& preMult(const Tensor& A);
  Tensor& rotate(Real alpha, t_ind axis);
 
  Tensor& outerProd(const Vec3& a, const Vec3& b);
 
  Real innerProd(const Tensor& T) const;
 
  t_ind dim() const { return n; }
 
  size_t size() const { return v.size(); }
 
  virtual bool symmetric() const { return false; }
 
  bool equal(const Tensor& T, Real tol = Real(1.0e-6)) const;
 
  bool isZero(Real tol = Real(1.0e-6)) const;
 
  virtual Tensor& transpose();
  virtual Tensor& symmetrize();
  Tensor& shift(short int idx = 1);
 
  virtual Real trace() const;
 
  virtual Real det() const;
 
  virtual Real inverse(Real tol = Real(0));
 
  Vec3 rotVec() const;
 
  virtual std::ostream& print(std::ostream& os, int prec = 0) const;
 
  // Global operators
 
  friend Vec3 operator*(const Tensor& T, const Vec3& v);
  friend Vec3 operator*(const Vec3& v, const Tensor& T);
  friend Tensor operator*(const Tensor& A, const Tensor& B);
  friend Tensor operator*(Real a, const Tensor& T);
 
  friend std::ostream& operator<<(std::ostream& os, const Tensor& T)
  {
    return T.print(os);
  }
};
 
 
class SymmTensor : public Tensor
{
  bool redim(const t_ind nsd, bool with33 = false);
 
protected:
  t_ind index(t_ind i, t_ind j) const override
  {
    if (i == j)
      return i-1; // diagonal term
    else if (n == 2)
      return v.size()-1; // off-diagonal term (2D)
 
    if (i == j+1 || i+2 == j) std::swap(i,j);
    return i+2; // upper triangular term (3D)
  }
 
public:
  explicit SymmTensor(const t_ind nsd, bool with33 = false);
  SymmTensor(const std::vector<Real>& vec);
  SymmTensor(const SymmTensor& T) : Tensor(0) { this->copy(T); }
 
  SymmTensor& operator=(Real val);
 
  void copy(const SymmTensor& T);
 
  bool symmetric() const override { return true; }
 
  Tensor& transpose() override { return *this; }
  Tensor& symmetrize() override { return *this; }
 
  Real trace() const override;
 
  Real det() const override;
 
  Real inverse(Real tol = Real(0)) override;
 
  SymmTensor& transform(const Tensor& T);
 
  SymmTensor& rightCauchyGreen(const Tensor& F);
 
  SymmTensor& outerProd(const Vec3& u);
 
  Real L2norm(bool doSqrt = true) const;
  Real vonMises(bool doSqrt = true) const;
  bool principal(Vec3& p, bool sorted = true) const;
  bool principal(Vec3& p, Tensor& pdir) const;
  bool principal(Vec3& p, Vec3* pdir, int ndir = 0) const;
  bool principal(Vec3& p, SymmTensor* M) const;
 
  std::ostream& print(std::ostream& os, int prec = 0) const override;
 
  // Global operators
 
  friend SymmTensor operator+(const SymmTensor& T, Real a);
  friend SymmTensor operator-(const SymmTensor& T, Real a);
 
  friend SymmTensor operator*(Real a, const SymmTensor& T);
  friend SymmTensor operator*(const SymmTensor& A, const SymmTensor& B);
};
 
 
inline SymmTensor operator+(const SymmTensor& A, const SymmTensor& B)
{
  SymmTensor C(A); C += B; return C;
}
 
inline SymmTensor operator-(const SymmTensor& A, const SymmTensor& B)
{
  SymmTensor C(A); C -= B; return C;
}
 
inline Real ddot(const SymmTensor& A, const SymmTensor& B)
{
  return A.innerProd(B);
}
 
 
class LocalSystem
{
protected:
  LocalSystem() {}
 
public:
  virtual ~LocalSystem() {}
 
  virtual const Tensor& getTmat(const Vec3& X) const = 0;
 
  static int patch; 
};
 
#endif