functions.h

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys is free software; you can redistribute it and/or modify it under the      *
 * terms of the GNU General Public License as published by the Free Software         *
 * Foundation; either version 2 of the License, or (at your option) any later        *
 * version.                                                                          *
 *                                                                                   *
 * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
 * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
 *                                                                                   *
 * You should have received a copy of the GNU General Public License along with      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

#ifndef MECHSYS_TENSORS_FUNCTIONS_H
#define MECHSYS_TENSORS_FUNCTIONS_H

#ifdef HAVE_CONFIG_H
  #include "config.h"
#else
  #ifndef REAL
    #define REAL double
  #endif
#endif

#include <cassert> // for assert()
#include <cmath>   // for sqrt()

#include "tensors/tensors.h"
#include "tensors/jacobirot.h"
#include "util/exception.h"

namespace Tensors
{



void Mult(Tensor2 const & S,  
          Tensor2 const & T,  
          Tensor2       & R   
);


REAL Norm(Tensor2 const & T 
);

REAL Det(Tensor2 const & T);


bool Inv(Tensor2 const & T,  
         Tensor2       & R   
);


bool Eigenvals(Tensor2 const & T,    
               REAL            L[3]  
);

// {{{
// }}}
bool Eigenvp(Tensor2 const & T,    
             REAL            L[3], 
             Tensor2         P[3]  
);


bool Sqrt(Tensor2 const & T, 
          Tensor2       & R  
);

void CharInvs(Tensor2 const & T,   
              REAL            I[3] 
);


void Strain_Ev_Ed(Tensor2 const & Eps, 
                  REAL          & Ev,  
                  REAL          & Ed   
);


void Stress_p_q(Tensor2 const & Sig, 
                REAL          & p,   
                REAL          & q    
);


inline REAL Stress_q(Tensor2 const & Sig 
);


inline REAL Sin3ThDev(Tensor2 const & S);

void Stress_p_q_S_sin3th(Tensor2 const & Sig   ,
                         REAL          & p     ,
                         REAL          & q     ,
                         Tensor2       & S     , 
                         REAL          & sin3th  
);


bool Stress_tn_ts(Tensor2 const & Sig, 
                  REAL          & tn,  
                  REAL          & ts   
);


void Stress_P_Q(Tensor2 const & Sig, 
                REAL          & P,   
                REAL          & Q    
);

void Stress_P_Q_S_sin3th(Tensor2 const & Sig   , 
                         REAL          & P     , 
                         REAL          & Q     , 
                         Tensor2       & S     , 
                         REAL          & sin3th
);

inline REAL Sin3Th(REAL SI, REAL SII, REAL SIII);

inline void Hid2Sig(REAL const &    p, REAL const &    q, REAL const &   th,
                    REAL       & Sig1, REAL       & Sig2, REAL       & Sig3);

inline void Hid2Sig(REAL const *  P, REAL const *   Q, REAL const *    T,
                    REAL       * SI, REAL       * SII, REAL       * SIII, int Size); // S# must be already allocated !!!

inline void Hid2Sig_(REAL const & SX, REAL const &  SY, REAL const &    p,
                     REAL       & SI, REAL       & SII, REAL       & SIII);
 // end of group TensorFunctions




// Tensor multiplication:  R = S*T
inline void Mult(Tensor2 const & S, Tensor2 const & T, Tensor2 & R)  // {{{
{
    REAL sq2 = sqrt(2.0);

    R(0) = S(0)*T(0)     + S(3)*T(3)/2.0 + S(5)*T(5)/2.0;
    R(1) = S(3)*T(3)/2.0 + S(1)*T(1)     + S(4)*T(4)/2.0;
    R(2) = S(5)*T(5)/2.0 + S(4)*T(4)/2.0 + S(2)*T(2);

    R(3) = S(0)*T(3)     + S(3)*T(1)     + S(5)*T(4)/sq2;
    R(4) = S(3)*T(5)/sq2 + S(1)*T(4)     + S(4)*T(2);
    R(5) = S(0)*T(5)     + S(3)*T(4)/sq2 + S(5)*T(2);

#ifndef NDEBUG
    REAL A = S(3)*T(0)     + S(1)*T(3)     + S(4)*T(5)/sq2;
    REAL B = S(5)*T(3)/sq2 + S(4)*T(1)     + S(2)*T(4);
    REAL C = S(5)*T(0)     + S(4)*T(3)/sq2 + S(2)*T(5);
    //std::cout << _15_8 R(3)-A << " | " << _15_8 R(4)-B << " | " << _15_8 R(5)-C << std::endl;
    const REAL ZERO=1.0e-10;
    if (!(fabs(R(3)-A)<ZERO && fabs(R(4)-B)<ZERO && fabs(R(5)-C)<ZERO))
        throw new Fatal("Tensors::Mult: Result of matrix multiplication is NOT symmetric. assert(fabs(R(3)-A)<ZERO & fabs(R(4)-B)<ZERO & fabs(R(5)-C)<ZERO) failed)");
#endif
} // }}}

inline REAL Norm(Tensor2 const & T)  // {{{
{
    return sqrt(T(0)*T(0) + T(1)*T(1) + T(2)*T(2) + T(3)*T(3) + T(4)*T(4) + T(5)*T(5));
} // }}}
    
// Determinant of T
inline REAL Det(Tensor2 const & T)  // {{{
{
    return T(0)*T(1)*T(2) + T(3)*T(4)*T(5)/sqrt(2.0) - T(0)*T(4)*T(4)/2.0 - T(1)*T(5)*T(5)/2.0 - T(2)*T(3)*T(3)/2.0;
} // }}}

// Inverse of a tensor T:  R = T^(-1)
inline bool Inv(Tensor2 const & T, Tensor2 & R)  // {{{
{
    REAL sq2 = sqrt(2.0);
    REAL det = Det(T);
    
    const REAL ZERO=1.0e-10;
    if (fabs(det)<ZERO) return false;

    R(0) =     ( T(1)*T(2)     - T(4)*T(4)/2.0 )/det;
    R(1) =     ( T(2)*T(0)     - T(5)*T(5)/2.0 )/det;
    R(2) =     ( T(0)*T(1)     - T(3)*T(3)/2.0 )/det;
    R(3) = sq2*( T(4)*T(5)/2.0 - T(2)*T(3)/sq2 )/det;
    R(4) = sq2*( T(5)*T(3)/2.0 - T(0)*T(4)/sq2 )/det;
    R(5) = sq2*( T(3)*T(4)/2.0 - T(1)*T(5)/sq2 )/det;

    return true;
} // }}}

// Eigenvalues (L) of a tensor T
inline bool Eigenvals(Tensor2 const & T, REAL L[3])  // {{{
{
    // Calculate eigenvalues and eigenvectors
    if (JacobiRot(T,L)==-1) return false;

    return true;
} // }}}

// Eigenvalues (L) and Eigenprojectors (P) of a tensor T
inline bool Eigenvp(Tensor2 const & T, REAL L[3], Tensor2 P[3])  // {{{
{
    REAL V0[3]; // eigenvector
    REAL V1[3]; // eigenvector
    REAL V2[3]; // eigenvector

    // Calculate eigenvalues and eigenvectors
    if (JacobiRot(T,V0,V1,V2,L)==-1) return false;

    REAL sq2 = sqrt(2.0);

    // Calculate eigenprojectors (Pi = Vi dyad Vi)
    // Eigenprojector 0
    P[0](0) = V0[0] * V0[0];
    P[0](1) = V0[1] * V0[1];
    P[0](2) = V0[2] * V0[2];
    P[0](3) = V0[0] * V0[1] * sq2;
    P[0](4) = V0[1] * V0[2] * sq2;
    P[0](5) = V0[0] * V0[2] * sq2;
    // Eigenprojector 1
    P[1](0) = V1[0] * V1[0];
    P[1](1) = V1[1] * V1[1];
    P[1](2) = V1[2] * V1[2];
    P[1](3) = V1[0] * V1[1] * sq2;
    P[1](4) = V1[1] * V1[2] * sq2;
    P[1](5) = V1[0] * V1[2] * sq2;
    // Eigenprojector 2
    P[2](0) = V2[0] * V2[0];
    P[2](1) = V2[1] * V2[1];
    P[2](2) = V2[2] * V2[2];
    P[2](3) = V2[0] * V2[1] * sq2;
    P[2](4) = V2[1] * V2[2] * sq2;
    P[2](5) = V2[0] * V2[2] * sq2;

    return true;
} // }}}

// Square root of a tensor T:  R = T^(0.5)
inline bool Sqrt(Tensor2 const & T, Tensor2 & R)  // {{{
{
    REAL    L[3]; // Eigenvalues
    Tensor2 P[3]; // Eigenprojectors

    // Calculate eigenvalues and eigenprojectors
    if (!Eigenvp(T,L,P)) return false;
    
    // Calculate R = sqrt(T)
    R(0)=0.0;
    R(1)=0.0;
    R(2)=0.0;
    R(3)=0.0;
    R(4)=0.0;
    R(5)=0.0;
    for (int i=0; i<3; ++i)
        R = sqrt(L[i])*P[i] + R; // R <- sqrt(L[i])*P[i] + R

    return true;
} // }}}

// Characteristic invariants of a symmetric second order tensor
inline void CharInvs(Tensor2 const & T, REAL I[3])  // {{{
{
    I[0] = T(0)+T(1)+T(2);
    I[1] = T(0)*T(1) + T(1)*T(2) + T(2)*T(0) - (T(3)*T(3) + T(4)*T(4) + T(5)*T(5))/2.0;
    I[2] = T(0)*T(1)*T(2) + T(3)*T(4)*T(5)/sqrt(2.0) - (T(0)*T(4)*T(4) + T(1)*T(5)*T(5) + T(2)*T(3)*T(3))/2.0;
} // }}}

inline void Strain_Ev_Ed (Tensor2 const & Eps, REAL & Ev, REAL & Ed)  // {{{
{
    Ev = Eps(0)+Eps(1)+Eps(2);
    Ed = sqrt(2.0*(   (Eps(0)-Eps(1))*(Eps(0)-Eps(1))
                    + (Eps(1)-Eps(2))*(Eps(1)-Eps(2))
                    + (Eps(2)-Eps(0))*(Eps(2)-Eps(0)) + 3.0*(Eps(3)*Eps(3) + Eps(4)*Eps(4) + Eps(5)*Eps(5)) ))/3.0;
} // }}}

inline void Stress_p_q(Tensor2 const & Sig, REAL & p, REAL & q)  // {{{
{
    p = (Sig(0)+Sig(1)+Sig(2))/3.0;
    q = sqrt((   (Sig(0)-Sig(1))*(Sig(0)-Sig(1))
               + (Sig(1)-Sig(2))*(Sig(1)-Sig(2))
               + (Sig(2)-Sig(0))*(Sig(2)-Sig(0)) + 3.0*(Sig(3)*Sig(3) + Sig(4)*Sig(4) + Sig(5)*Sig(5)) )/2.0);
} // }}}

inline REAL Stress_q(Tensor2 const & Sig) // {{{
{
    return sqrt((   (Sig(0)-Sig(1))*(Sig(0)-Sig(1))
                  + (Sig(1)-Sig(2))*(Sig(1)-Sig(2))
                  + (Sig(2)-Sig(0))*(Sig(2)-Sig(0)) + 3.0*(Sig(3)*Sig(3) + Sig(4)*Sig(4) + Sig(5)*Sig(5)) )/2.0);
} // }}}

inline REAL Sin3ThDev(Tensor2 const & S) // {{{
{
    const REAL ZZERO = 1.0e-5;
    REAL normS = Norm(S);
    if (normS<=ZZERO) return -1.0;
    else
    {
        REAL     L = 9.0*sqrt(2.0)*Det(S)/(sqrt(3.0)*pow(normS,3.0));
        if      (L>= 1.0) return  1.0;
        else if (L<=-1.0) return -1.0;
        else              return  L;
    }
} // }}}

inline void Stress_p_q_S_sin3th(Tensor2 const & Sig, REAL & p,  REAL & q, Tensor2 & S, REAL & sin3th) // {{{
{
    Stress_p_q(Sig,p,q);
    S = Sig - I * p;
    const REAL    ZZERO  = 1.0e-7;
    if (q<=ZZERO) sin3th = -1.0; // 3th==-90 => th=30
    else          sin3th = 27.0*Det(S)/(2.0*pow(q,3.0));
    if (sin3th>= 1.0) sin3th =  1.0;
    if (sin3th<=-1.0) sin3th = -1.0;
} // }}}

inline bool Stress_tn_ts(Tensor2 const & Sig, REAL & tn, REAL & ts)  // {{{
{
    REAL I[3];
    CharInvs(Sig,I);
    const REAL ZERO=1.0e-10;
    if (fabs(I[1])<=ZERO) return false; // Second invariant equal zero

    REAL s = I[0]*I[1]*I[2]-9.0*I[2]*I[2];
    if (s<0.0) return false; // in ts, avoid sqrt(s) NaN

    tn = 3.0*I[2]/I[1];
    if (tn<ZERO) return false; // tn must be >= 0

    ts = sqrt(s)/I[1];
    if (ts!=ts) return false;  // NaN values

    return true;
} // }}}

inline void Stress_P_Q(Tensor2 const & Sig, REAL & P,  REAL & Q) // {{{
{
    P = (Sig(0)+Sig(1)+Sig(2))/sqrt(3.0);
    Q = sqrt( ((Sig(0)-Sig(1))*(Sig(0)-Sig(1))
             + (Sig(1)-Sig(2))*(Sig(1)-Sig(2))
             + (Sig(2)-Sig(0))*(Sig(2)-Sig(0)))/3.0 + Sig(3)*Sig(3) + Sig(4)*Sig(4) + Sig(5)*Sig(5) );
} // }}}

inline void Stress_P_Q_S_sin3th(Tensor2 const & Sig, REAL & P, REAL & Q, Tensor2 & S, REAL & sin3th) // {{{
{
    Stress_P_Q(Sig,P,Q);
    S = Sig - I * (P/sqrt(3.0));
    const REAL    ZZERO  = 1.0e-7;
    if (Q<=ZZERO) sin3th = -1.0; // 3th==-90 => th=30
    else          sin3th = 27.0*Det(S)/(2.0*pow(sqrt(3.0)*Q/sqrt(2.0),3.0));
    if (sin3th>= 1.0) sin3th =  1.0;
    if (sin3th<=-1.0) sin3th = -1.0;
} // }}}

inline REAL Sin3Th(REAL SI, REAL SII, REAL SIII)  // {{{
{
    /*       \->T   | S1  
     *        \  A<-|   
     *         \    |     
     *     -30  \   |         
     *       ',  \  |
     *         '. \ |
     *           '.\|
     *            .' -.                
     *          .'     '.
     *        .'         '.
     *   S2 .'             '. S3    */

    REAL       p = (SI+SII+SIII)/3.0;
    REAL devS[3] = {SI-p, SII-p, SIII-p};
    REAL    det  = devS[0]*devS[1]*devS[2];
    REAL    norm = sqrt(devS[0]*devS[0] + devS[1]*devS[1] + devS[2]*devS[2]);
    if (norm<=1.0e-5) return -1.0;
    else
    {
        REAL     L = 9.0*sqrt(2.0)*det/(sqrt(3.0)*pow(norm,3.0));
        if      (L>= 1.0) return  1.0;
        else if (L<=-1.0) return -1.0;
        else              return  L;
    }
} // }}}

inline void Hid2Sig(REAL const &    p, REAL const &    q, REAL const &   th, // {{{
                    REAL       & Sig1, REAL       & Sig2, REAL       & Sig3)
{
    /*  Converts hidrostatic coordinates to sigma (principal) coord (T in radians)

                   SI|      r = (Sx^2+Sy^2)^0.5
                     |Sy    z = Sz
                 \ T |
                  \  |
                   \ |
                    \|
                   .' -.--------> Sx
                 .'     '.
               .'         '.
         SII .'             '. SIII      */

    // Constants
    REAL sq2    = sqrt(2.0);
    REAL sq6    = sqrt(6.0);
    REAL sq2by3 = sq2/sqrt(3.0);
    REAL     r  = q*sq2by3;
    REAL     Sx = -r*sin(th);
    REAL     Sy =  r*cos(th);
           Sig1 = p + 2.0*Sy/sq6;
           Sig2 = p -     Sy/sq6 - Sx/sq2;
           Sig3 = p -     Sy/sq6 + Sx/sq2;
} // }}}

inline void Hid2Sig(REAL const *  P, REAL const *   Q, REAL const *    T, // {{{
                    REAL       * SI, REAL       * SII, REAL       * SIII, int Size)
{
    /*  Converts hidrostatic coordinates to sigma (principal) coord (T in radians)

                   SI|      r = (Sx^2+Sy^2)^0.5
                     |Sy    z = Sz
                 \ T |
                  \  |
                   \ |
                    \|
                   .' -.--------> Sx
                 .'     '.
               .'         '.
         SII .'             '. SIII      */

    // Constants
    REAL sq2    = sqrt(2.0);
    REAL sq6    = sqrt(6.0);
    REAL sq2by3 = sq2/sqrt(3.0);
    for (int i=0; i<Size; ++i)
    {
        REAL r  = Q[i]*sq2by3;
        REAL Sx = -r*sin(T[i]);
        REAL Sy =  r*cos(T[i]);
        SI  [i] = P[i] + 2.0*Sy/sq6;
        SII [i] = P[i] -     Sy/sq6 - Sx/sq2;
        SIII[i] = P[i] -     Sy/sq6 + Sx/sq2;
    }
} // }}}

inline void Hid2Sig_(REAL const & SX, REAL const &  SY, REAL const &    p, // {{{
                     REAL       & SI, REAL       & SII, REAL       & SIII)
{
    REAL sq2  = sqrt(2.0);
    REAL sq6  = sqrt(6.0);
         SI   = p + 2.0*SY/sq6;
         SII  = p -     SY/sq6 - SX/sq2;
         SIII = p -     SY/sq6 + SX/sq2;
} // }}}

}; // namespace Tensors

#endif // MECHSYS_TENSORS_FUNCTIONS_H

// vim:fdm=marker

---