tijtensor.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_TIJTENSOR_H
#define MECHSYS_TIJTENSOR_H

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

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

using Tensors::Tensor2;
using Tensors::CharInvs;
using Tensors::Stress_tn_ts;
using Tensors::Sqrt;
using Tensors::Inv;
using Tensors::Mult;
using Tensors::I;

namespace Tensors
{

class TijTensor
{
    static REAL ZERO;
public:

    TijTensor(REAL TsMin) : _ts_min(TsMin) {}

    void Derivs(const Tensor2 & Sig,       // In: Stress tensor
                REAL          & tn,        // Out: Mean stress invariant
                REAL          & ts,        // Out: Deviatoric stress invariant
                REAL            Is[3],     // Out: Characterist invariants of stress (Sig)
                REAL          & sqrt_rz,   // Out: Sqrt of I3/I2
                Tensor2       & Tau,       // Out: Tensor Sqrt of Sig
                Tensor2       & inv_Tau,   // Out: Inverse of tensor Tau
                Tensor2       & a,         // Out: "normal-to-SMP" tensor == Derivative of tn w.r.t tij tensor
                Tensor2       & t,         // Out: Modified tij tensor
                Tensor2       & dtn_ds,    // Out: Derivative of tn w.r.t Sig
                Tensor2       & dts_ds,    // Out: Derivative of ts w.r.t Sig
                Tensor2       & dts_dt     // Out: Derivative of ts w.r.t tij tensor
    );
private:
    REAL _ts_min; 
}; // class TijTensor

REAL TijTensor::ZERO = 1.0e-10;




inline void TijTensor::Derivs(const Tensor2 & Sig,     // {{{
                              REAL          & tn, 
                              REAL          & ts, 
                              REAL            Is[3],
                              REAL          & sqrt_rz, 
                              Tensor2       & Tau,
                              Tensor2       & inv_Tau,
                              Tensor2       & a,
                              Tensor2       & t,
                              Tensor2       & dtn_ds,
                              Tensor2       & dts_ds,
                              Tensor2       & dts_dt)
{
    // ---------------------------------------------------------------------- tij calculation
    
    // Characteristic invariants:  I1=Is[0], I2=Is[1], I3=Is[2]
    CharInvs(Sig, Is);
    if (fabs(Is[1])<ZERO)
        throw new Fatal(_("TijTensor::Derivs: assert(I2>0.0) failed (I2=%.6f). Sig=<%.6f; %.6f; %.6f; %.6f; %.6f; %.6f>"), Is[1],
                Sig(0), Sig(1), Sig(2), Sig(3)/sqrt(2.0), Sig(4)/sqrt(2.0), Sig(5)/sqrt(2.0));

    // Square root of I3/I2
    sqrt_rz = Is[2]/Is[1];
    if (fabs(sqrt_rz)<ZERO)
        throw new Fatal(_("TijTensor::Derivs: assert(I3/I2>0) failed (I3/I2=%.6f). Sig=<%.6f; %.6f; %.6f; %.6f; %.6f; %.6f>"), sqrt_rz,
                Sig(0), Sig(1), Sig(2), Sig(3)/sqrt(2.0), Sig(4)/sqrt(2.0), Sig(5)/sqrt(2.0));
    sqrt_rz = sqrt(sqrt_rz); // sqrt_rz = sqrt(I3/I2)

    // Professr Nakai's stress invariants
    if (!Stress_tn_ts(Sig, tn,ts))
        throw new Fatal(_("TijTensor::Derivs: tn and ts calculation failed. Sig=<%.6f; %.6f; %.6f; %.6f; %.6f; %.6f>"),
                Sig(0), Sig(1), Sig(2), Sig(3)/sqrt(2.0), Sig(4)/sqrt(2.0), Sig(5)/sqrt(2.0));

    // Tau = square root of Sig
    if (!Sqrt(Sig, Tau))
        throw new Fatal(_("TijTensor::Derivs: Tau = Sqrt(Sig) failed. Sig=<%.6f; %.6f; %.6f; %.6f; %.6f; %.6f>"),
                Sig(0), Sig(1), Sig(2), Sig(3)/sqrt(2.0), Sig(4)/sqrt(2.0), Sig(5)/sqrt(2.0));

    // inverse of Tau
    if (!Inv(Tau, inv_Tau))
        throw new Fatal(_("TijTensor::Derivs: inverse of Tau failed. Sig=<%.6f; %.6f; %.6f; %.6f; %.6f; %.6f>"),
                Sig(0), Sig(1), Sig(2), Sig(3)/sqrt(2.0), Sig(4)/sqrt(2.0), Sig(5)/sqrt(2.0));

    // a tensor
    a = inv_Tau * sqrt_rz;

    // tij tensor
    Mult(a,Sig, t);
    
    // -------------------------------------------------------------- derivatives calculation

    Tensor2 SigSig;
    Mult(Sig,Sig, SigSig);
    
    Tensor2 dI1_ds; dI1_ds = I*Is[0] - Sig;
    Tensor2 dI2_ds; dI2_ds = SigSig  - Sig*Is[0] + I*Is[1];
    
    dtn_ds = dI2_ds*(3.0/Is[1]) - dI1_ds*(tn/Is[1]);

    // Check ts
    if (fabs(ts)<_ts_min)
    {
        dts_ds = 0.0; // In fact, it is undetermined
        dts_dt = 0.0; // In fact, it is undetermined
    }
    else
    {
        REAL alp = tn/(6.0*ts);
        REAL bet = tn*(6.0*tn-Is[0])/(6.0*ts*Is[1]);
        REAL gam = (Is[0]-6.0*tn)/(2.0*ts*Is[1]);
        dts_ds = I*alp + dI1_ds*bet + dI2_ds*gam;
        dts_dt = t*(1.0/ts) - a*(tn/ts);
    }
        
} // }}}

}; // namespace Tensors

#endif // MECHSYS_TIJTENSOR_H

// vim:fdm=marker

---