genericep.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_GENERICEP_H
#define MECHSYS_GENERICEP_H

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

#define Z0_MINUS_Z1_TOL  1.0e-7
#define CP_ZERO          1.0e-5

#include "models/equilibmodel.h"
#include "tensors/tensors.h"
#include "tensors/operators.h"
#include "tensors/functions.h"
#include "numerical/autostepme.h"
#include "util/string.h"

using Tensors::Tensor2;      // Second order tensor
using Tensors::Reduce;       // Used in: Phi <- V:De:r - y0*HH0  and  dLam <- V:De:dEps/Phi
using Tensors::GerX;         // Used in: Dep <- (-1/Phi)(De:r)dy(V:De) + De
using Tensors::Dot;          // Dot product between 4th order tensors and 2nd order tensors (Dep:dEps, Cep:dSig)
using Tensors::Norm;         // (Second Euclidian) norm of a second order tensor
using blitz::dot;            // Blitz++: Inner product between TinyVectors
using Tensors::Ger;          // Used in: Cep <- (1/cp)*(r dyadic V) + Ce

template<unsigned int nIntVars>
class GenericEP : public EquilibModel
{
public:
    typedef blitz::TinyVector<REAL,nIntVars> IntVars;
    typedef blitz::TinyVector<REAL,nIntVars> HardMod;
    typedef blitz::TinyVector<REAL,nIntVars> T_dfdz;

    // Constructor
    GenericEP() : _num_dLam(1.0) {}

    // Destructor
    virtual ~GenericEP() {}

    // Methods
    void TgStiffness  (Tensors::Tensor4 & D) const;
    void Actualize    (Tensors::Tensor2 const & DSig, Tensors::Tensor2 & DEps);
    void StressUpdate (Tensors::Tensor2 const & DEps, Tensors::Tensor2 & DSig);
    void BackupState  ();
    void RestoreState ();

    // Can be overriden
    virtual int  nAdditionalInternalStateValues () const { return 0; }
    virtual void AdditionalInternalStateValues  (Array<REAL>   & IntStateVals)  const {};
    virtual void AdditionalInternalStateNames   (Array<String> & IntStateNames) const {};

    // Access methods
    int  nInternalStateValues () const { return nIntVars+1+nAdditionalInternalStateValues(); }; // +1 => specific volume
    void InternalStateValues  (Array<REAL>   & IntStateVals ) const;
    void InternalStateNames   (Array<String> & IntStateNames) const;

protected:
    // Data
    REAL    _v; // State variable: Specific volume
    IntVars _z; // State variable: Internal quantities
    REAL    _v_bkp;
    IntVars _z_bkp;

    // Methods
    virtual void _calc_grads_hards(REAL    const & v  ,           // In:  Specific volume
                                   Tensor2 const & Sig,           // In:  Stress tensor
                                   IntVars const & z  ,           // In:  Internal variables
                                   Tensor2       & r  ,           // Out: Plastic strain direction
                                   Tensor2       & V  ,           // Out: Derivative of f ("loading" surface) w.r.t Sig (df_dSig)
                                   T_dfdz        & y  ,           // Out: Derivative of f ("loading" surface) w.r.t z
                                   HardMod       & HH ,           // Out: Hardening modulae
                                   REAL          & cp ) const =0; // Out: Plastic coeficient

    virtual void _calc_De(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & De) const =0;
    virtual void _calc_Ce(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & Ce) const =0;
    virtual void _unloading_update_int_vars() =0;
    virtual void _calc_dz(Tensors::Tensor2 const & dSig, Tensors::Tensor2 const & dEps, REAL const & dLam, HardMod const & HH, IntVars & dz) const =0;

private:
    // Data
    REAL _num_dLam; // numerator of dLam
    REAL _bkp_num_dLam;

    // Methods
    int _Dep_times_dEps(Tensor2 const & Eps,       // In: (actual) strain tensor
                        Tensor2 const & Sig,       // In: (actual) stress tensor
                        IntVars const & z,         // In: (actual) internal variables
                        Tensor2 const & dEps,      // In: Strain increment
                        Tensor2       & dSig,      // Out: Stress increment
                        IntVars       & dz) const; // Out: Internal variables increment
    // returns: -1 if failed, 0 to continue, 1 to stop
    
    int _Cep_times_dSig(Tensor2 const & Sig,       // In: (actual) stress tensor
                        Tensor2 const & Eps,       // In: (actual) strain tensor
                        IntVars const & z,         // In: (actual) internal variables
                        Tensor2 const & dSig,      // In: Stress increment
                        Tensor2       & dEps,      // Out: Strain increment
                        IntVars       & dz) const; // Out: Internal variables increment
    // returns: -1 if failed, 0 to continue, 1 to stop
    
    int _De_times_dEps(Tensor2 const & Eps,           // In: (actual) strain tensor
                       Tensor2 const & Sig,           // In: (actual) stress tensor
                       float   const & dummy1,        // In: 
                       Tensor2 const & dEps,          // In: Strain increment
                       Tensor2       & dSig,          // Out: Stress increment
                       float         & dummy2) const; // Out: 
    // returns: -1 if failed, 0 to continue, 1 to stop
    
    int _Ce_times_dSig(Tensor2 const & Sig,           // In: (actual) stress tensor
                       Tensor2 const & Eps,           // In: (actual) strain tensor
                       float   const & dummy1,        // In:
                       Tensor2 const & dSig,          // In: stress increment
                       Tensor2       & dEps,          // Out: strain increment
                       float         & dummy2) const; // Out:
    // returns: -1 if failed, 0 to continue, 1 to stop

    REAL _local_error(Tensor2 const & Ey    ,         // In: Stress or Strain error between low and high order evaluation (FE and ME)
                      Tensor2 const & y_high,         // In: Stress or Strain high order evaluation - "correct" (ME)
                      IntVars const & Ez    ,         // In: Error on internal variables
                      IntVars const & z_high) const;  // In: Higher evaluation of internal variables

    REAL _local_error(Tensor2 const & Ey    ,         // In: Stress or Strain error between low and high order evaluation (FE and ME)
                      Tensor2 const & y_high,         // In: Stress or Strain high order evaluation - "correct" (ME)
                      float   const & dummy1,         // In: dummy
                      float   const & dummy2) const;  // In: dummy
}; // class GenericEP




template<unsigned int nIntVars>
inline void GenericEP<nIntVars>::TgStiffness(Tensors::Tensor4 & D) const // {{{
{
    if (_num_dLam<0) // Elastic tangent stiffness
    {
        _calc_De(_v, _sig, _z, true, D); // D <- De
    }
    else // Elastoplastic tangent stiffness
    {
        // Gradients and Hardening
        Tensor2 r,V;    T_dfdz y;    HardMod HH;    REAL cp;
        _calc_grads_hards(_v, _sig, _z, r, V, y, HH, cp);

        // Dep and Phi
        Tensors::Tensor4 De;
        _calc_De(_v, _sig, _z, false, De);
        REAL Phi = Reduce(V,De,r)+cp;   // Phi = V:De:r - y0*HH0
        GerX(-1.0/Phi,De,r,V,De,De, D); // D <- Dep = (-1/Phi) * (De:r)dy(V:De) + De
    }
} // }}}

template<unsigned int nIntVars>
inline void GenericEP<nIntVars>::Actualize(Tensors::Tensor2 const & DSig, Tensors::Tensor2 & DEps) // {{{
{
    // Initial _eps
    Tensor2 eps_ini = _eps;

    // Gradients and Hardening
    Tensor2 r,V;    T_dfdz y;    HardMod HH;    REAL cp;
    _calc_grads_hards(_v, _sig, _z, r, V, y, HH, cp);

    // Plastic (Lagrange) multiplier
    REAL dLam = Dot(V,DSig)/cp; // dLam = V:DSig/cp

    // Check unloading/loading condition
    if (dLam<0.0) // Unloading
    {
        if (EquilibModel::_isc.Type()==IntegSchemesCtes::FE)
        {
            // Forward-Euler - Actualize (Elastic)
            int    n_div = _isc.FE_ndiv();
            Tensor2 dsig = DSig/n_div;
            float dummy;
            for (int i=0; i<n_div; ++i)
            {
                _Ce_times_dSig(_sig, _eps, dummy, dsig, DEps, dummy);
                _sig += dsig;
                _eps += DEps;
            }
        }
        else
        {
            // Allocate a Time Integration object
            AutoStepME<Tensor2, Tensor2, float, GenericEP<nIntVars> > TI(EquilibModel::_isc);

            // Update stress (_sig) and strain (_eps) state
            float dummy;
            int nss = TI.Solve(this,                                  // Address of this instance of GenericEP class
                               &GenericEP<nIntVars>::_Ce_times_dSig , // Pointer to a function that calculates dEps = C:dSig
                               &GenericEP<nIntVars>::_local_error   , // Pointer to a function that calculates local error on dEps
                               _sig, _eps, dummy, DSig);
            // Check num_sub_steps
            if (nss==-1)
                throw new Fatal(_("GenericEP::Actualize: (Unloading) Number of max sub-steps (%d) was not sufficient for AutoStepME."), EquilibModel::_isc.ME_maxSS());
        }
        // Update internal state
        _unloading_update_int_vars();
    }
    else // Loading
    {
        if (EquilibModel::_isc.Type()==IntegSchemesCtes::FE)
        {
            // Forward-Euler - Actualize (Elastoplastic)
            int    n_div = _isc.FE_ndiv();
            Tensor2 dsig = DSig/n_div;
            for (int i=0; i<n_div; ++i)
            {
                IntVars dz;
                _Cep_times_dSig(_sig, _eps, _z, dsig, DEps, dz);
                _sig += dsig;
                _eps += DEps;
                _z   += dz;
            }
        }
        else
        {
            // Allocate a Time Integration object
            AutoStepME<Tensor2, Tensor2, IntVars, GenericEP<nIntVars> > TI(EquilibModel::_isc);

            // Update stress (_sig), internal variables (_z) and strain (_eps) state
            int nss = TI.Solve(this,                                  // Address of this instance of GenericEP class
                               &GenericEP<nIntVars>::_Cep_times_dSig, // Pointer to a function that calculates dEps = C:dSig
                               &GenericEP<nIntVars>::_local_error   , // Pointer to a function that calculates local error on dEps
                               _sig, _eps, _z, DSig);
            // Check num_sub_steps
            if (nss==-1)
                throw new Fatal(_("GenericEP::Actualize: (Loading) Number of max sub-steps (%d) was not sufficient for AutoStepME."), EquilibModel::_isc.ME_maxSS());
        }
    }

    // Calculate the total (secant) increment of strain
    DEps = _eps - eps_ini;

    // Update internal state (secant)
    _v += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v

} // }}}

template<unsigned int nIntVars>
inline void GenericEP<nIntVars>::StressUpdate(Tensors::Tensor2 const & DEps, Tensors::Tensor2 & DSig) // {{{
{
    // Initial _sig
    Tensor2 sig_ini = _sig;

    // Gradients and Hardening
    Tensor2 r,V;    T_dfdz y;    HardMod HH;    REAL cp;
    _calc_grads_hards(_v, _sig, _z, r, V, y, HH, cp);

    // Calculate numerator of dLam
    Tensors::Tensor4 De; 
    if (_num_dLam<0) _calc_De(_v, _sig, _z, true,  De);
    else             _calc_De(_v, _sig, _z, false, De);
    _num_dLam = Reduce(V,De,DEps);

    // Check unloading/loading condition
    if (_num_dLam<0.0) // Unloading
    {
        if (EquilibModel::_isc.Type()==IntegSchemesCtes::FE)
        {
            // Forward-Euler - Update (Elastic)
            int    n_div = _isc.FE_ndiv();
            Tensor2 deps = DEps/n_div;
            float dummy;
            for (int i=0; i<n_div; ++i)
            {
                _De_times_dEps(_eps, _sig, dummy, deps, DSig, dummy);
                _sig += DSig;
                _eps += deps;
            }
        }
        else
        {
            // Allocate a Time Integration object
            AutoStepME<Tensor2, Tensor2, float, GenericEP<nIntVars> > TI(EquilibModel::_isc);

            // Update stress (_sig) and strain (_eps) state
            float dummy;
            int nss = TI.Solve(this,                                  // Address of this instance of GenericEP class
                               &GenericEP<nIntVars>::_De_times_dEps , // Pointer to a function that calculates dEps = C:dSig
                               &GenericEP<nIntVars>::_local_error   , // Pointer to a function that calculates local error on dEps
                               _eps, _sig, dummy, DEps);
            // Check num_sub_steps
            if (nss==-1)
                throw new Fatal(_("GenericEP::Actualize: (Unloading) Number of max sub-steps (%d) was not sufficient for AutoStepME."), EquilibModel::_isc.ME_maxSS());
        }
        // Update internal state
        _unloading_update_int_vars();
    }
    else // Loading
    {
        if (EquilibModel::_isc.Type()==IntegSchemesCtes::FE)
        {
            // Forward-Euler - Update (Elastoplastic)
            int    n_div = _isc.FE_ndiv();
            Tensor2 deps = DEps/n_div;
            for (int i=0; i<n_div; ++i)
            {
                IntVars dz;
                _Dep_times_dEps(_eps, _sig, _z, deps, DSig, dz);
                _sig += DSig;
                _eps += deps;
                _z   += dz;
            }
        }
        else
        {
            // Allocate a Time Integration object
            AutoStepME<Tensor2, Tensor2, IntVars, GenericEP<nIntVars> > TI(EquilibModel::_isc);

            // Update stress (_sig), internal variables (_z) and strain (_eps) state
            int nss = TI.Solve(this,                                  // Address of this instance of GenericEP class
                               &GenericEP<nIntVars>::_Dep_times_dEps, // Pointer to a function that calculates dSig = D:dEps
                               &GenericEP<nIntVars>::_local_error ,   // Pointer to a function that calculates local error on dSig
                               _eps, _sig, _z, DEps);
            // Check num_sub_steps
            if (nss==-1)
                throw new Fatal(_("GenericEP::StressUpdate: (Loading) Number of max sub-steps (%d) was not sufficient for AutoStepME."), EquilibModel::_isc.ME_maxSS());
        }
    }

    // Calculate the total (secant) increment of stress
    DSig = _sig - sig_ini;

    // Update internal state
    _v += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v

} // }}}

template<unsigned int nIntVars>
inline void GenericEP<nIntVars>::BackupState() // {{{
{
    _sig_bkp = _sig;
    _eps_bkp = _eps;
      _v_bkp = _v;
      _z_bkp = _z;
    _bkp_num_dLam = _num_dLam;
} // }}}

template<unsigned int nIntVars>
inline void GenericEP<nIntVars>::RestoreState() // {{{
{
    _sig = _sig_bkp;
    _eps = _eps_bkp;
      _v =   _v_bkp;
      _z =   _z_bkp;
    _num_dLam = _bkp_num_dLam;
} // }}}

template<unsigned int nIntVars>
inline int GenericEP<nIntVars>::_Dep_times_dEps(Tensor2 const & Eps,       // In: (actual) strain tensor {{{
                                                Tensor2 const & Sig,       // In: (actual) stress tensor
                                                IntVars const & z,         // In: (actual) internal variables
                                                Tensor2 const & dEps,      // In: Strain increment
                                                Tensor2       & dSig,      // Out: Stress increment
                                                IntVars       & dz) const  // Out: Internal variables increment
{
    // Gradients and Hardening
    REAL v=_v; // using (constant) current specific volume
    Tensor2 r,V;    T_dfdz y;    HardMod HH;    REAL cp;
    _calc_grads_hards(v, Sig, z, r, V, y, HH, cp);

    // Dep and Phi
    Tensors::Tensor4 De,Dep;
    _calc_De(v, Sig, z, false, De);
    REAL Phi = Reduce(V,De,r)+cp;    // Phi = V:De:r - y0*HH0
    GerX(-1.0/Phi,De,r,V,De,De,Dep); // Dep = (-1/Phi) * (De:r)dy(V:De) + De

    // dSig and dLam
    Dot(Dep, dEps, dSig);              // dSig = Dep:dEps
    REAL dLam = Reduce(V,De,dEps)/Phi; // dLam = V:De:dEps/Phi

    // dz
    _calc_dz(dSig,dEps,dLam,HH, dz);

    return 0;
} // }}}

template<unsigned int nIntVars>
inline int GenericEP<nIntVars>::_Cep_times_dSig(Tensor2 const & Sig,       // In: (actual) stress tensor {{{
                                                Tensor2 const & Eps,       // In: (actual) strain tensor
                                                IntVars const & z,         // In: (actual) internal variables
                                                Tensor2 const & dSig,      // In: Stress increment
                                                Tensor2       & dEps,      // Out: Strain increment
                                                IntVars       & dz) const  // Out: Internal variables increment
{
    // Gradients and Hardening
    REAL v=_v; // using (constant) current specific volume
    Tensor2 r,V;    T_dfdz y;    HardMod HH;    REAL cp;
    _calc_grads_hards(v, Sig, z, r, V, y, HH, cp);

    // Try to check for softening initial point
    if (fabs(cp/(_sig(0)+_sig(1)+_sig(2)))<=CP_ZERO)
        throw new Warning(_("GenericEP<%s>::_Cep_times_dSig: hp=%f. Maybe it is the start of softening behaviour"),this->Name().GetSTL().c_str(),cp);

    // Stop if excessive deviatoric strains
    REAL Ev,Ed; Tensors::Strain_Ev_Ed(Eps, Ev, Ed);
    if (Ed>EquilibModel::_isc.EdMax())
        throw new Warning(_("GenericEP<%s>::_Cep_times_dSig: Ed=%f (%) greater than EdMax=%f (%)"),this->Name().GetSTL().c_str(),Ed*100.0,EquilibModel::_isc.EdMax()*100);
    //std::cout << v << Sig << z << r << V << y << HH << cp << std::endl;

    // Cep
    Tensors::Tensor4 Cep;
    _calc_Ce(v, Sig, z, false, Cep); // Cep <- Ce
    Ger(1.0/cp,r,V,Cep);          // Cep <- (1/cp)*(r dyadic V) + Ce

    // dEps and dLam
    Dot(Cep, dSig, dEps);       // dEps <- Cep:dSig
    REAL dLam = Dot(V,dSig)/cp; // dLam <- V:dSig/cp

    // dz
    _calc_dz(dSig,dEps,dLam,HH, dz);

    return 0;
} // }}}

template<unsigned int nIntVars>
inline int GenericEP<nIntVars>::_De_times_dEps(Tensor2 const & Eps,           // In: (actual) strain tensor {{{
                                               Tensor2 const & Sig,           // In: (actual) stress tensor
                                               float   const & dummy1,        // In:
                                               Tensor2 const & dEps,          // In: Strain increment
                                               Tensor2       & dSig,          // Out: Stress increment
                                               float         & dummy2) const  // Out:
{
    Tensors::Tensor4 De;
    _calc_De(_v, Sig, _z, true, De);
    Dot(De,dEps,dSig);
    return 0; // returns: -1 if failed, 0 to continue, 1 to stop
} // }}}

template<unsigned int nIntVars>
inline int GenericEP<nIntVars>::_Ce_times_dSig(Tensor2 const & Sig, // {{{
                                               Tensor2 const & Eps,
                                               float   const & dummy1,
                                               Tensor2 const & dSig,
                                               Tensor2       & dEps,
                                               float         & dummy2) const
{
    Tensors::Tensor4 Ce;
    _calc_Ce(_v, Sig, _z, true, Ce);
    Dot(Ce,dSig,dEps);
    return 0; // returns: -1 if failed, 0 to continue, 1 to stop
} // }}}

template<unsigned int nIntVars>
inline REAL GenericEP<nIntVars>::_local_error(Tensor2 const & Ey    ,         // In: Stress or Strain error between low and high order evaluation (FE and ME) {{{
                                              Tensor2 const & y_high,         // In: Stress or Strain high order evaluation - "correct" (ME)
                                              IntVars const & Ez    ,         // In: Error on internal variables
                                              IntVars const & z_high) const   // In: Higher evaluation of internal variables
{
    REAL Err_y = Norm(Ey)/Norm(y_high);
    REAL Err_z = sqrt(dot(Ez,Ez))/sqrt(dot(z_high,z_high));
    return (Err_y>Err_z ? Err_y : Err_z);
} // }}}

template<unsigned int nIntVars>
inline REAL GenericEP<nIntVars>::_local_error(Tensor2 const & Ey    ,         // In: Stress or Strain error between low and high order evaluation (FE and ME) {{{
                                              Tensor2 const & y_high,         // In: Stress or Strain high order evaluation - "correct" (ME)
                                              float   const & dummy1,         // In:
                                              float   const & dummy2) const   // In:
{
    return Norm(Ey)/Norm(y_high);
} // }}}

template<unsigned int nIntVars>
inline void GenericEP<nIntVars>::InternalStateValues(Array<REAL> & IntStateVals) const // {{{
{
    int nIntStateVals = nIntVars+1+nAdditionalInternalStateValues(); // +1 => specific volume component
    IntStateVals.resize(nIntStateVals);
    for (unsigned int i_intvar=0; i_intvar<nIntVars; ++i_intvar)
        IntStateVals[i_intvar] = _z(i_intvar);
        IntStateVals[nIntVars] = _v;
    if (nAdditionalInternalStateValues()>0)
    {
        Array<REAL> values;
        AdditionalInternalStateValues(values);
        for (int i=0; i<nAdditionalInternalStateValues(); ++i)
            IntStateVals[nIntVars+1+i] = values[i];
    }
} // }}}

template<unsigned int nIntVars>
inline void GenericEP<nIntVars>::InternalStateNames(Array<String> & IntStateNames) const // {{{
{
    int nIntStateNames = nIntVars+1+nAdditionalInternalStateValues(); // +1 => specific volume component
    IntStateNames.resize(nIntStateNames);
    for (unsigned int i_intvar=0; i_intvar<nIntVars; ++i_intvar)
    {
        std::ostringstream str_z; str_z << "z" << i_intvar;
        IntStateNames[i_intvar] = str_z.str();
    }
    IntStateNames[nIntVars] = "v";
    if (nAdditionalInternalStateValues()>0)
    {
        Array<String> names;
        AdditionalInternalStateNames(names);
        for (int i=0; i<nAdditionalInternalStateValues(); ++i)
            IntStateNames[nIntVars+1+i] = names[i];
    }
} // }}}

#endif // MECHSYS_GENERICEP_H

// vim:fdm=marker

---