elastic.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_ELASTIC_H
#define MECHSYS_ELASTIC_H

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

#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;

class Elastic : public EquilibModel
{
public:
    // Constructor
    Elastic(Array<REAL> const & Prms, Array<REAL> const & IniData);
    String Name() const { return "Elastic"; }
    void TgStiffness(Tensors::Tensor4 & D) const;
    void BackupState();
    void StressUpdate(Tensors::Tensor2 const & DEps, Tensors::Tensor2 & DSig);
    void Actualize   (Tensors::Tensor2 const & DSig, Tensors::Tensor2 & DEps) {}
    void RestoreState();
    int  nInternalStateValues() const { return 0; };
    void InternalStateValues(Array<REAL> & IntStateVals)  const {}
    void InternalStateNames(Array<String> & IntStateNames)  const {}
private:
    // Data
    REAL _E;       // Young
    REAL _nu;      // Poisson
    REAL _kappa;   // Edometric recompression inclination
    bool _non_lin; // Lin/non-Lin model
    REAL _v_ini;   // Initial specific volume
    // Methods
    int _D_times_dEps(Tensor2 const & Eps,           // In: (actual) strain tensor
                      Tensor2 const & Sig,           // In: (actual) stress tensor
                      REAL    const & dummy1,        // In: dummy
                      Tensor2 const & dEps,          // In: Strain increment
                      Tensor2       & dSig,          // Out: Stress increment
                      REAL          & dummy2) const; // Out: dummy
    // 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)
                      REAL    const & dummy1 ,        // In: dummy
                      REAL    const & dummy2) const;  // In: dummy

}; // class Elastic




inline Elastic::Elastic(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
        : _E  (Prms[0]),
          _nu (Prms[1])
{
    // Check number of parameters (Lin/non-Lin model)
    if (Prms.size()==3) { _non_lin=true; _kappa=Prms[2]; }
    else _non_lin=false;

    // Check number of initial data
    if (IniData.size()!=4) // Layered analysis
        throw new Fatal(_("Elastic::Constructor: Initial data not sufficiet (4) for a layered analysis. { Sv' Sh' vini OCR }\n"));

    // Read initial data
    REAL sig_vert = IniData[0]; // Vertical stress component
    REAL sig_horz = IniData[1]; // Horizontal stress component
           _v_ini = IniData[2]; // Initial specific volume
           //_OCR = IniData[3]; // Over consolidation rate (not used here)

    // Fill stress and strain tensors
    _sig = sig_horz, sig_horz, sig_vert, 0.0, 0.0, 0.0; // Mandel representation
    _eps = 0.0,0.0,0.0,0.0,0.0,0.0;
} // }}}

inline void Elastic::TgStiffness(Tensors::Tensor4 & De) const // {{{
{
    // Calculate Young modulus
    REAL E;
    if (_non_lin) E=(_sig(0)+_sig(1)+_sig(2)) * (1.0-2.0*_nu) * _v_ini / _kappa;
    else          E=_E;

    // Calculate elastic tangent tensor
    Tensors::AddScaled(     E/(1.0+_nu)                   ,  Tensors::IIsym ,
                        _nu*E/( (1.0+_nu)*(1.0-2.0*_nu) ) , Tensors::IdyI  , De ); // Z = a*X + b*Y
} // }}}

inline void Elastic::BackupState() // {{{
{
    _sig_bkp = _sig;
    _eps_bkp = _eps;
} // }}}

inline void Elastic::StressUpdate(Tensors::Tensor2 const & DEps, Tensors::Tensor2 & DSig) // {{{
{
    // Check if this model is linear or non-linear
    if (_non_lin)
    {
        // Allocate a Time Integration object
        AutoStepME<Tensor2, Tensor2, REAL, Elastic > TI(EquilibModel::_isc);

        // Save a temporary copy of _sig
        DSig = _sig;

        // Update stress (_sig) and strain (_eps) state
        REAL dummy=0;
        int num_sub_steps = TI.Solve(this,                    // Address of this instance of Elastic class
                                     &Elastic::_D_times_dEps, // Pointer to a function that calculates dSig = D:dEps
                                     &Elastic::_local_error , // Pointer to a function that calculates local error on dSig
                                     _eps, _sig, dummy, DEps);

        //cout << TI.Results() << endl;

        // Check num_sub_steps
        if (num_sub_steps==-1)
            throw new Fatal(_("Elastic::StressUpdate: Number of max sub-steps (%d) not sufficient for AutoStepME"), EquilibModel::_isc.ME_maxSS());

        // Calculate the total increment of stress
        DSig = _sig - DSig;
    }
    else
    {
        // Just a forward-Euler step
        Tensors::Tensor4 De;
        TgStiffness(De);
        Tensors::Dot(De,DEps, DSig);
    }
} // }}}

inline int Elastic::_D_times_dEps(Tensor2 const & Eps,           // In: (actual) strain tensor {{{
                                  Tensor2 const & Sig,           // In: (actual) stress tensor
                                  REAL    const & dummy1,        // In: dummy
                                  Tensor2 const & dEps,          // In: Strain increment
                                  Tensor2       & dSig,          // Out: Stress increment
                                  REAL          & dummy2) const  // Out: dummy
{
    Tensors::Tensor4 De;
    TgStiffness(De);
    Tensors::Dot(De,dEps, dSig);
    return 0; // returns: -1 if failed, 0 to continue, 1 to stop
} // }}}

inline REAL Elastic::_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)
                                  REAL    const & dummy1 ,        // In: dummy
                                  REAL    const & dummy2) const   // In: dummy
{
    return Tensors::Norm(Ey)/Tensors::Norm(y_high);
} // }}}

inline void Elastic::RestoreState() // {{{
{
    _sig = _sig_bkp;
    _eps = _eps_bkp;
} // }}}




// Allocate a new Elastic model
EquilibModel * ElasticMaker(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
{
    return new Elastic(Prms, IniData);
} // }}}

// Register an Elastic model into ModelFactory array map
int ElasticRegister() // {{{
{
    EquilibModelFactory["Elastic"] = ElasticMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __Elastic_dummy_int = ElasticRegister();

#endif // MECHSYS_ELASTIC_H

// vim:fdm=marker

---