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:
    SymTensor() { } 
    SymTensor(const SymTensor &s)=default;
    explicit SymTensor(double i11,double i22=0.0,double i33=0.0,double i12=0.0,double i13=0.0,double i23=0.0)
                      : d11_(i11), d22_(i22), d33_(i33), d12_(i12), d13_(i13), d23_(i23) { }
    SymTensor &operator=(const SymTensor &s)=default;
 
    bool operator==(const SymTensor &s) const;
    bool operator<(const SymTensor &s) const;
 
    double s11() const { return d11_; }
    double s22() const { return d22_; } 
    double s33() const { return d33_; } 
    double s12() const { return d12_; } 
    double s13() const { return d13_; } 
    double s23() const { return d23_; } 
    double s21() const { return d12_; } 
    double s31() const { return d13_; } 
    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 double operator[](unsigned int i) const;
    inline double &operator[](unsigned int i);
 
    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_ = 0.0;
    double d22_ = 0.0;
    double d33_ = 0.0;
    double d12_ = 0.0;
    double d13_ = 0.0;
    double d23_ = 0.0;
};
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 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 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