linearelastic.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_LINEARELASTIC_H
#define MECHSYS_LINEARELASTIC_H

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

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

class LinearElastic : public EquilibModel
{
public:
    // Constants
    static const size_t NPRMS = 2;
    static const size_t NIDAT = 5;

    // Constructor
    LinearElastic(Array<REAL> const & Prms, Array<REAL> const & IniData);

    // Derived Methods
    String Name         () const { return "LinearElastic"; }
    void   TgStiffness  (Tensors::Tensor4 & D) const;
    void   StressUpdate (Tensors::Tensor2 const & DEps, Tensors::Tensor2 & DSig);
    void   Actualize    (Tensors::Tensor2 const & DSig, Tensors::Tensor2 & DEps);
    void   BackupState  ();
    void   RestoreState ();

    // Access Methods
    int  nInternalStateValues ()                              const { return 1;                                        }
    void InternalStateValues  (Array<REAL>   & IntStateVals ) const { IntStateVals .resize(1); IntStateVals [0] = _v;  }
    void InternalStateNames   (Array<String> & IntStateNames) const { IntStateNames.resize(1); IntStateNames[0] = "v"; }

private:
    // Parameters
    REAL _E;  // Young
    REAL _nu; // Poisson
    REAL _v;  // Specific volume
    REAL _v_bkp;

}; // class LinearElastic




inline LinearElastic::LinearElastic(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
{

    // ##################################################################### Setup Parameters

    /* example.mat 
     * #----------------------------
     * #              Young  Poisson 
     * LinearElastic  20000     0.25 
     */
    
    // Check number of parameters 
    if (Prms.size()!=NPRMS) 
        throw new Fatal(_("LinearElastic::Constructor: The number of parameters (%d) is incorrect. It must be equal to %d \n"));

    // Fill the parameters
    _E  = Prms[0];
    _nu = Prms[1];

    // Check number of initial data
    if (IniData.size()!=NIDAT) 
        throw new Fatal(_("LinearElastic::Constructor: The number of initial data is not sufficient (it must be equal to %d). { Sx' Sy' Sz' vini(not-used) OCR(not-used) }"), NIDAT);

    // Fill stress and strain tensors
    _sig = IniData[0], IniData[1], IniData[2], 0.0, 0.0, 0.0; // Mandel representation
    _eps = 0.0,0.0,0.0,0.0,0.0,0.0;

    // Specific volume
    _v = IniData[3];

} // }}}

inline void LinearElastic::TgStiffness(Tensors::Tensor4 & De) const // {{{
{
    // 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 LinearElastic::StressUpdate(Tensors::Tensor2 const & DEps, Tensors::Tensor2 & DSig) // {{{
{
    // Just a forward-Euler step
    Tensors::Tensor4 De;
    TgStiffness(De);
    Tensors::Dot(De,DEps, DSig);

    // Update stress and strains
    _sig += DSig;
    _eps += DEps;

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

} // }}}

inline void LinearElastic::Actualize(Tensors::Tensor2 const & DSig, Tensors::Tensor2 & DEps) // {{{
{
    // Calculate elastic tangent tensor
    Tensors::Tensor4 Ce;
    Tensors::AddScaled( (1.0+_nu)/_E , Tensors::IIsym ,
                            - _nu/_E , Tensors::IdyI  , Ce ); // Z = a*X + b*Y

    // Just a forward-Euler step
    Tensors::Dot(Ce,DSig, DEps);

    // Update stress and strains
    _sig += DSig;
    _eps += DEps;

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

} // }}}

inline void LinearElastic::BackupState() // {{{
{
    _sig_bkp = _sig;
    _eps_bkp = _eps;
    _v_bkp   = _v;
} // }}}

inline void LinearElastic::RestoreState() // {{{
{
    _sig = _sig_bkp;
    _eps = _eps_bkp;
    _v   = _v_bkp;
} // }}}




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

// Register an LinearElastic model into ModelFactory array map
int LinearElasticRegister() // {{{
{
    EquilibModelFactory["LinearElastic"] = LinearElasticMaker;
    return 0;
} // }}}

// Execute the autoregistration
int __LinearElastic_dummy_int = LinearElasticRegister();

#endif // MECHSYS_LINEARELASTIC_H

// vim:fdm=marker

---