symtensor.h

Go to the documentation of this file.

#pragma once
#include "axes.h"
#include <cassert>
#include <limits>
#include <algorithm>
 
class SymTensorInfo;
 
PUSHWARNING
VSWARNING(26495) // member init (not initialized on purpose)
VSWARNING(4702) // Weird unreachable code warnings in release
class BASE_EXPORT SymTensor {
public:
#ifdef DEBUG
 SymTensor() { d11_ = d22_ = d33_ = d12_ = d13_ = d23_ = initVal<T>(); }
#else
 SymTensor() { } // NOLINT
#endif
 SymTensor(const SymTensor &s) { operator=(s); }
 explicit SymTensor(const double &i11,const double &i22=0.0,const double &i33=0.0,
 const double &i12=0.0,const double &i13=0.0,const double &i23=0.0)
 : d11_(i11), d22_(i22), d33_(i33), d12_(i12), d13_(i13), d23_(i23) { }
 const SymTensor &operator=(const SymTensor &s);
 
 bool operator==(const SymTensor &s) const;
 bool operator<(const SymTensor &s) const;
 
 const double &s11() const { return d11_; }
 const double &s22() const { return d22_; } 
 const double &s33() const { return d33_; } 
 const double &s12() const { return d12_; } 
 const double &s13() const { return d13_; } 
 const double &s23() const { return d23_; } 
 const double &s21() const { return d12_; } 
 const double &s31() const { return d13_; } 
 const double &s32() const { return d23_; } 
 double &rs11() { return d11_; }
 double &rs22() { return d22_; } 
 double &rs33() { return d33_; } 
 double &rs12() { return d12_; } 
 double &rs13() { return d13_; } 
 double &rs23() { return d23_; } 
 double &rs21() { return d12_; } 
 double &rs31() { return d13_; } 
 double &rs32() { return d23_; } 
 inline const double &operator[](unsigned int i) const;
 inline double &operator[](unsigned int i);
 
 const double &operator()(unsigned int i,unsigned int j) const;
 double &operator()(unsigned int i,unsigned int j);
 
 DVect3 getEigenInfo(SymTensorInfo *si=0) const;
 double getTrace() const { return d11_+d22_+d33_; }
 SymTensor getDeviatoric() const { double p=getTrace()/3.0; return SymTensor(d11_-p, d22_-p, d33_-p, d12_, d13_, d23_); }
 double getI1() const { return getTrace(); }
 double getI2() const { return d11_*d22_ + d22_*d33_ + d11_*d33_ - d12_*d12_ - d23_*d23_ - d13_*d13_; }
 double getI3() const { return getDeterminate(); }
 double getJ2() const;
 double getJ2(SymTensor *dev, double *I1=nullptr);
 //Keep only getJ2(), Von Mises and Octahedral stress/strain are calculated from J2.
 //double getVonMises() const { return sqrt(3.0*getJ2()); }
 //double getVonMisesStrain() const { return sqrt(4.0*getJ2()/9.0); }
 //double getOctahedralShear() const { return sqrt(2.0*getJ2()/3.0); }
 //double getOctahedralShearStrain() const { return sqrt(8.0*getJ2()/3.0); }
 double getJ3() const;
 double getLode(double *I1 = nullptr, double *J2 = nullptr, double *J3 = nullptr);
 double getDeterminate() const;
 double getNorm2() const { return d11_*d11_ + d22_*d22_ + d33_*d33_ + 2.0*(d12_*d12_ + d13_*d13_ + d23_*d23_); }
 double getTotalMeasure() const;
 SymTensor toGlobal(const Axes3D &a) const;
 inline DVect3 operator*(const DVect3 &input) const;
 inline DVect2 operator*(const DVect2 &input) const;
 inline SymTensor operator*(const double &mul) const;
 inline const SymTensor &operator*=(const double &mul);
 inline SymTensor operator/(const double &mul) const;
 inline const SymTensor &operator/=(const double &mul);
 SymTensor mul(const double &d) const { SymTensor ret(d11_*d,d22_*d,d33_*d,d12_*d,d13_*d,d23_*d); return ret; }
 const SymTensor &operator+=(const SymTensor &s) { d11_+=s.d11_; d22_+=s.d22_; d33_+=s.d33_; d12_+=s.d12_; d13_+=s.d13_; d23_+=s.d23_; return *this; }
 const SymTensor &operator-=(const SymTensor &s) { d11_-=s.d11_; d22_-=s.d22_; d33_-=s.d33_; d12_-=s.d12_; d13_-=s.d13_; d23_-=s.d23_; return *this; }
 inline SymTensor operator+(const SymTensor &s) const;
 inline SymTensor operator-(const SymTensor &s) const;
 static SymTensor fromPrincipal(const DVect3 &prin,const Axes3D &axes);
 static SymTensor fromForceNormal(const DVect3 &normal, const DVect3 &force);
 static inline uint32 doubleToSingleComponent(uint32 dof1,uint32 dof2);
 bool isDiagonal(double tol = std::numeric_limits<double>::epsilon()*1000.0) const { return (abs(d12_) > tol || abs(d13_) > tol || abs(d23_) > tol) ? false : true; }
 inline bool isIsotropic(double tol = std::numeric_limits<double>::epsilon()*1000.0) const;
 inline void adjustTrace(const double &newTrace);
 inline void incrementDiagonal(const double &increment) { d11_ += increment; d22_ += increment; d33_ += increment; }
 inline void rotate(const DVect3 &rot);
 inline double maxAbs() const { return std::max(std::abs(d11_),std::max(std::abs(d22_),std::max(std::abs(d33_),std::max(std::abs(d12_),std::max(std::abs(d13_),std::abs(d23_)))))); }
 inline bool zero() const { return (d11_==0.0) and (d22_==0.0) and (d33_==0.0) and (d12_==0.0) and (d13_==0.0) and (d23_==0.0); }
private:
 double d11_;
 double d22_;
 double d33_;
 double d12_;
 double d13_;
 double d23_;
};
POPWARNING
 
class SymTensorInfo {
public:
 friend class SymTensor;
 BASE_EXPORT SymTensorInfo(): type_(Type::ThreeDCube) { }
 BASE_EXPORT SymTensorInfo(const SymTensorInfo &si);
 BASE_EXPORT const SymTensorInfo &operator=(const SymTensorInfo &si);
 BASE_EXPORT Axes3D getAxes() const;
 BASE_EXPORT SymTensor resolve(const DVect3 &prin) const;
private:
 enum class Type { ThreeDCube, ThreeDJacobi, ZMax, ZMid, ZMin };
 Type type_;
 Axes3D axes_;
};
 
inline const SymTensor &SymTensor::operator=(const SymTensor &s) {
 d11_ = s.d11_;
 d22_ = s.d22_;
 d33_ = s.d33_;
 d12_ = s.d12_;
 d13_ = s.d13_;
 d23_ = s.d23_;
 return *this;
}
 
inline bool SymTensor::operator==(const SymTensor &s) const {
 return (d11_ == s.d11_ && d22_ == s.d22_ &&
 d33_ == s.d33_ && d12_ == s.d12_ && d13_ == s.d13_ &&
 d23_ == s.d23_);
}
 
inline const double &SymTensor::operator[](unsigned int i) const {
 assert(i<6);
 switch (i) {
 case 0: return d11_;
 case 1: return d22_;
 case 2: return d33_;
 case 3: return d12_;
 case 4: return d13_;
 case 5: return d23_;
 }
 return d11_;
}
 
inline double &SymTensor::operator[](unsigned int i) {
 assert(i<6);
 switch (i) {
 case 0: return d11_;
 case 1: return d22_;
 case 2: return d33_;
 case 3: return d12_;
 case 4: return d13_;
 case 5: return d23_;
 }
 return d11_;
}
 
inline double SymTensor::getTotalMeasure() const {
 double I1 = getTrace();
 double J2 = getJ2();
 return sqrt(I1*I1/3.0 + 2.0*J2);
}
 
inline DVect3 SymTensor::operator*(const DVect3 &normal) const {
 return DVect3(normal.x()*s11() + normal.y()*s12() + normal.z()*s13(),
 normal.x()*s21() + normal.y()*s22() + normal.z()*s23(),
 normal.x()*s31() + normal.y()*s32() + normal.z()*s33());
}
 
inline DVect2 SymTensor::operator*(const DVect2 &normal) const {
 return DVect2(normal.x()*s11() + normal.y()*s12(),
 normal.x()*s21() + normal.y()*s22());
}
 
inline SymTensor SymTensor::operator*(const double &mul) const {
 SymTensor ret;
 ret.d11_ = d11_ * mul;
 ret.d22_ = d22_ * mul;
 ret.d33_ = d33_ * mul;
 ret.d12_ = d12_ * mul;
 ret.d13_ = d13_ * mul;
 ret.d23_ = d23_ * mul;
 return ret;
}
 
inline const SymTensor &SymTensor::operator*=(const double &mul) {
 d11_ *= mul;
 d22_ *= mul;
 d33_ *= mul;
 d12_ *= mul;
 d13_ *= mul;
 d23_ *= mul;
 return *this;
}
 
inline SymTensor SymTensor::operator/(const double &mul) const {
 SymTensor ret;
 ret.d11_ = d11_ / mul;
 ret.d22_ = d22_ / mul;
 ret.d33_ = d33_ / mul;
 ret.d12_ = d12_ / mul;
 ret.d13_ = d13_ / mul;
 ret.d23_ = d23_ / mul;
 return ret;
}
 
inline const SymTensor &SymTensor::operator/=(const double &mul) {
 d11_ /= mul;
 d22_ /= mul;
 d33_ /= mul;
 d12_ /= mul;
 d13_ /= mul;
 d23_ /= mul;
 return *this;
}
 
inline SymTensor SymTensor::operator+(const SymTensor &s) const {
 SymTensor ret;
 ret.d11_ = d11_ + s.d11_;
 ret.d22_ = d22_ + s.d22_;
 ret.d33_ = d33_ + s.d33_;
 ret.d12_ = d12_ + s.d12_;
 ret.d13_ = d13_ + s.d13_;
 ret.d23_ = d23_ + s.d23_;
 return ret;
}
 
inline SymTensor SymTensor::operator-(const SymTensor &s) const {
 SymTensor ret;
 ret.d11_ = d11_ - s.d11_;
 ret.d22_ = d22_ - s.d22_;
 ret.d33_ = d33_ - s.d33_;
 ret.d12_ = d12_ - s.d12_;
 ret.d13_ = d13_ - s.d13_;
 ret.d23_ = d23_ - s.d23_;
 return ret;
}
 
inline bool SymTensor::isIsotropic(double tol) const {
 double dtol = std::max(std::max(std::abs(s11()),std::abs(s22())),std::abs(s33())) * tol;
 if (std::abs(s11()-s22()) > dtol) return false;
 if (std::abs(s11()-s33()) > dtol) return false;
 if (std::abs(s12()) > dtol) return false;
 if (std::abs(s13()) > dtol) return false;
 if (std::abs(s23()) > dtol) return false;
 return true;
}
 
inline void SymTensor::adjustTrace(const double &newTrace) {
 static constexpr double d1d3 = 1.0 / 3.0;
 /* --- strain differences --- */
 double dx = s11() - s22();
 double dy = s22() - s33();
 double dz = s33() - s11();
 rs11() = (newTrace + dx - dz) * d1d3;
 rs22() = (newTrace + dy - dx) * d1d3;
 rs33() = (newTrace + dz - dy) * d1d3;
}
 
inline void SymTensor::rotate(const DVect3 &rot) {
 SymTensor copy(*this);
 d11_ += 2.0*( copy.s12()*rot.x() + copy.s13()*rot.y());
 d22_ += 2.0*(-copy.s12()*rot.x() + copy.s23()*rot.z());
 d33_ += -2.0*( copy.s13()*rot.y() + copy.s23()*rot.z());
 d12_ += ((copy.s22() - copy.s11())*rot.x() + copy.s23()*rot.x() + copy.s13()*rot.z());
 d13_ += ( copy.s23()*rot.x() + (copy.s33() - copy.s11())*rot.x() - copy.s12()*rot.z());
 d23_ += ( -copy.s13()*rot.x() - copy.s12()*rot.x() + (copy.s33() - copy.s22())*rot.z());
}
 
// 11 0
// 22 1
// 33 2
// 12 3
// 13 4
// 23 5
uint32 SymTensor::doubleToSingleComponent(uint32 dof1,uint32 dof2) {
 dof1 = std::clamp<uint32>(dof1,0,2);
 dof2 = std::clamp<uint32>(dof2,0,2);
 switch (dof1) {
 case 0:
 switch (dof2) {
 case 0: return 0;
 case 1: return 3;
 case 2: return 4;
 }
 break;
 case 1:
 switch (dof2) {
 case 0: return 3;
 case 1: return 1;
 case 2: return 5;
 }
 break;
 case 2:
 switch (dof2) {
 case 0: return 4;
 case 1: return 5;
 case 2: return 3;
 }
 break;
 }
 return 0;
}
 
namespace base {
 BASE_EXPORT string ts(const SymTensor &s, int width=0, char notation = '\0', int precision = -1, char fill = ' ');
}
 
// EoF