models

subbar.h

00001 /*************************************************************************************
00002  * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
00003  * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
00004  * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
00005  *                                                                                   *
00006  * This file is part of MechSys.                                                     *
00007  *                                                                                   *
00008  * MechSys is free software; you can redistribute it and/or modify it under the      *
00009  * terms of the GNU General Public License as published by the Free Software         *
00010  * Foundation; either version 2 of the License, or (at your option) any later        *
00011  * version.                                                                          *
00012  *                                                                                   *
00013  * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
00014  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
00015  * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
00016  *                                                                                   *
00017  * You should have received a copy of the GNU General Public License along with      *
00018  * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
00019  * Fifth Floor, Boston, MA 02110-1301, USA                                           *
00020  *************************************************************************************/
00021 
00022 #ifndef MECHSYS_SUBBAR_H
00023 #define MECHSYS_SUBBAR_H
00024 
00025 #ifdef HAVE_CONFIG_H
00026   #include "config.h"
00027 #else
00028   #ifndef REAL
00029     #define REAL double
00030   #endif
00031 #endif
00032 
00033 //#define SUBBAR_INIT1 1
00034 
00035 // STL
00036 #include <iostream>
00037 
00038 // Blitz++
00039 #include <blitz/tinyvec-et.h>
00040 
00041 // MechSys
00042 #include "models/coupled/coupledmodel.h"
00043 #include "tensors/tensors.h"
00044 #include "tensors/operators.h"
00045 #include "tensors/functions.h"
00046 #include "numerical/integschemesctes.h"
00047 #include "numerical/autostepme.h"
00048 #include "util/string.h"
00049 #include "util/array.h"
00050 #include "util/util.h"
00051 
00052 using Tensors::I;
00053 using Tensors::Tensor2;
00054 using Tensors::Tensor4;
00055 using Tensors::Stress_p_q;
00056 using Tensors::Strain_Ev_Ed;
00057 using Tensors::Stress_p_q_S_sin3th;
00058 using Tensors::Mult;
00059 
00061 class SubBar : public CoupledModel
00062 {
00063     static REAL Q_ZERO; // min value of q (used to calc derivatives)
00064 public:
00065 
00067     typedef blitz::TinyVector<REAL,4> IntVars; // {z0,z1,z2,z3} == {p0star_subload, s0_subload, p0star, s0}
00068     typedef blitz::TinyVector<REAL,4> HardMod; // {HH0,HH1,HH2,HH3}
00069     typedef blitz::TinyVector<REAL,4> T_dFdz;  // {y0,y1,y2,y3} derivative of f w.r.t. intvars
00070     typedef blitz::TinyVector<REAL,7> AugVec;  // {Sig/Eps,Suc} Augmented tensor with suction
00071 
00072     // Constructor
00073     SubBar(Array<REAL> const & Prms, Array<REAL> const & IniData);
00074 
00075     // Methods
00076     String Name            () const { return "SubBar"; }
00077     void   TgStiffness     (Tensor4 & D, Tensor2 & d) const;
00078     REAL   Sr              () const; // Degree of saturation
00079     REAL   phi             () const; // Sig_net = Sig_tot - phi*Suc*I
00080     REAL   n_times_dSr_dPp () const; // n:porosity times dSr_dPp (derivative of degree of sat. w.r.t. porepressure)
00081     void   TgPermeability  (Tensor2 & k) const;
00082     void   FlowVelocity    (Tensor1 const & Grad, Tensor1 & Vel) const;
00083     void   Actualize       (Tensor2 const & DSig, REAL const & DPp, Tensor2 & DEps, REAL & DnSr);
00084     void   StressUpdate    (Tensor2 const & DEps, REAL const & DPp, Tensor2 & DSig, REAL & DnSr);
00085     void   FlowUpdate      (REAL const & DPp) {}
00086     void   BackupState     ();
00087     void   RestoreState    ();
00088 
00089     // Methods
00090     REAL YieldFunc          (Tensor2 const & Sig, IntVars const & z, REAL const & Suc)                        const;
00091     REAL YieldFunc_Normal   (Tensor2 const & Sig, IntVars const & z, REAL const & Suc)                        const;
00092     void Calc_dFdSig_dFdSuc (Tensor2 const & Sig, IntVars const & z, REAL const & Suc, Tensor2 & V, REAL & S) const;
00093     REAL Calc_p0            (IntVars const & z, REAL const & Suc) const; // p0(z0,Suc)
00094     REAL Calc_C             (IntVars const & z, REAL const & Suc) const; // C (z1,Suc)
00095 
00096     // Derived Access methods
00097     Tensor2 const & Sig() const { return  _sig;   }
00098     Tensor2 const & Eps() const { return  _eps;   }
00099     REAL            Pp () const { return  _ppr;   }
00100     REAL        GammaW () const { return  _gamaW; } 
00101     int  nInternalStateValues() const { return 5; }
00102     void InternalStateValues (Array<REAL>   & IntStateVals ) const;
00103     void InternalStateNames  (Array<String> & IntStateNames) const;
00104 
00105     // Methods to help plot
00106     REAL GetM(REAL Theta) const { REAL sin3th=sin(3.0*Theta); return _calc_M(sin3th); }
00107 
00108 private:
00109     // State vars
00110     REAL    _v;   
00111     Tensor2 _sig; 
00112     REAL    _suc; 
00113     Tensor2 _eps; 
00114     IntVars _z;   
00115     REAL    _ppr; 
00116     
00117     // Backup state vars
00118     REAL    _bkp_v;   
00119     Tensor2 _bkp_sig; 
00120     REAL    _bkp_suc; 
00121     Tensor2 _bkp_eps; 
00122     IntVars _bkp_z;   
00123     REAL    _bkp_ppr; 
00124 
00125     // Parameters (Mechanic)
00126     REAL _lam0;
00127     REAL _kap;
00128     REAL _phics;
00129     REAL _G;
00130     REAL _r;
00131     REAL _bet;
00132     REAL _k;
00133     REAL _lams;
00134     REAL _kaps;
00135     REAL _patm;
00136     REAL _pref;
00137     REAL _BB;
00138 
00139     // Parameters (Hydraulic)
00140     REAL _a;
00141     REAL _b;
00142     REAL _c;
00143     REAL _m;
00144     REAL _ks;
00145     REAL _gamaW;
00146 
00147     // Parameters (Subloading)
00148     REAL _c0;
00149     REAL _cs;
00150 
00151     // Derived constants
00152     REAL _Mcs;
00153     REAL _wcs;
00154     REAL _alp;
00155 
00156     // Auxiliar variables
00157     REAL _num_dLam; // numerator of dLam
00158     REAL _bkp_num_dLam;
00159 
00160     // Constants
00161     REAL _MAX_ERR_F; // Max error on F (Yield Function)
00162     REAL _MIN_HP;    // Min value of hp (plastic coeficient)
00163 
00164     // Private methods
00165     REAL _lam        (REAL const & Suc) const;
00166     REAL _psi        (REAL const & Suc) const;
00167 
00168     void _calc_Ce_ce(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & Ce, Tensor2 & ce) const;
00169     void _calc_De_de(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & De, Tensor2 & de) const;
00170 
00171     void _calc_f0_grads(Tensor2 const & Sig     ,
00172                         IntVars const & z       ,
00173                         REAL    const & Suc     ,
00174                         Tensor2       & df0_dSig,
00175                         REAL          & df0_dSuc,
00176                         REAL          & df0_dz0 ) const;
00177 
00178     void _calc_gradients_and_hardening(REAL    const & v  ,        // Specific volume
00179                                        Tensor2 const & Sig,        // Stress
00180                                        IntVars const & z  ,        // Internal variables
00181                                        REAL    const & Suc,        // Suction
00182                                        Tensor2       & r  ,        // Plastic flow direction
00183                                        Tensor2       & V  ,        // dF_dSig
00184                                        T_dFdz        & y  ,        // dF_dz
00185                                        REAL          & S  ,        // dF_dSuc
00186                                        HardMod       & HH ,        // Hardening
00187                                        REAL          & hp ) const; // Plastic coeficient
00188 
00189     int _tg_deps_dz_given_dsig_and_dsuc(AugVec  const & SigSuc ,
00190                                         Tensor2 const & Eps    ,
00191                                         IntVars const & z      ,
00192                                         AugVec  const & dSigSuc,
00193                                         Tensor2       & dEps   ,
00194                                         IntVars       & dz     ) const;
00195 
00196     int _tg_dsig_dz_given_deps_and_dsuc(AugVec  const & EpsSuc ,
00197                                         Tensor2 const & Sig    ,
00198                                         IntVars const & z      ,
00199                                         AugVec  const & dEpsSuc,
00200                                         Tensor2       & dSig   ,
00201                                         IntVars       & dz     ) const;
00202 
00203     REAL _local_error(Tensor2 const & Ey, Tensor2 const & y_high,
00204                       IntVars const & Ez, IntVars const & z_high) const;
00205 
00207     void _force_consistency();
00208 
00209     // Local
00210     REAL _calc_M(REAL const & sin3th) const;
00211 
00212 }; // class SubBar
00213 
00214 REAL SubBar::Q_ZERO = 1.0e-7;
00215 
00216 
00218 
00219 
00220 inline SubBar::SubBar(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
00221     : _num_dLam(1.0)
00222 {
00223 
00224     // ##################################################################### Setup Parameters
00225 
00226     /* example.mat
00227        #------------------------------------------------------------------------------------------------------------------------------
00228        #          0     1       2     3     4     5    6     7      8      9     10   11   12   13   14   15     16     17    18    19
00229        #       lam0   kap   phics     G     r   bet    k  lams   kaps   patm   pref    B    a    b    c    m     ks  gamaW    c0    cs
00230        SubBar  0.20  0.02  25.377  10.0  0.75  12.5  0.6  0.08  0.008  0.101  0.101 1000  0.1  0.5  1.0  2.0  10e-8   0.01  5000  1000
00231     */
00232 
00233     // Check number of parameters
00234     const size_t n_prms = 20;
00235     if (Prms.size()!=n_prms)
00236         throw new Fatal(_("SubBar::SubBar: The number of parameters (%d) is incorrect. It must be equal to %d"), Prms.size(), n_prms);
00237 
00238     // Parameters (mechanical)
00239     _lam0  = Prms[0];
00240     _kap   = Prms[1];
00241     _phics = Prms[2];
00242     _G     = Prms[3];
00243     _r     = Prms[4];
00244     _bet   = Prms[5];
00245     _k     = Prms[6];
00246     _lams  = Prms[7];
00247     _kaps  = Prms[8];
00248     _patm  = Prms[9];
00249     _pref  = Prms[10];
00250     _BB    = Prms[11];
00251 
00252     // Parameters (hydraulic)
00253     _a     = Prms[12];
00254     _b     = Prms[13];
00255     _c     = Prms[14];
00256     _m     = Prms[15];
00257     _ks    = Prms[16];
00258     _gamaW = Prms[17];
00259 
00260     // Parameters (subloading)
00261     _c0 = Prms[18];
00262     _cs = Prms[19];
00263 
00264     // Derived constants
00265     REAL sin_phi = sin(Util::ToRad(_phics));
00266     _Mcs = 6.0*sin_phi/(3.0-sin_phi);
00267     _wcs = pow((3.0-sin_phi)/(3.0+sin_phi),4.0);
00268     _alp = (_Mcs*(_Mcs-9.0)*(_Mcs-3.0)) / (9.0*(6.0-_Mcs)*(1.0-_kap/_lam0));
00269     
00270     // ##################################################################### Setup Initial State
00271 
00272     // Check number of initial data
00273     if (IniData.size()!=7) // Layered analysis
00274         throw new Fatal(_("SubBar::SubBar: The number of initial data is not sufficiet (it must be equal to 5). { SigX SigY SigZ vini OCR OSI Pp }"));
00275 
00276     // Read initial data
00277     REAL sig_x = IniData[0]; // X stress component
00278     REAL sig_y = IniData[1]; // Y stress component
00279     REAL sig_z = IniData[2]; // Z stress component
00280             _v = IniData[3]; // Initial specific volume
00281     REAL   OCR = IniData[4]; // Over Consolidation Ratio
00282     REAL   OSI = IniData[5]; // Over Suction Increment = Suc_{max} - Suc_{actual}
00283     REAL    Pp = IniData[6]; // Initial pore-pressure
00284 
00285     // Fill stress and strain tensors
00286     REAL sq2=Util::Sq2();
00287     _sig = sig_x, sig_y, sig_z, 0.0*sq2, 0.0*sq2, 0.0*sq2; // Mandel representation
00288     _eps = 0.0,0.0,0.0,0.0*sq2,0.0*sq2,0.0*sq2;
00289 
00290     // Calculate initial state
00291     _ppr = Pp;       // Initial pore-pressure
00292     _suc = 0.0 - Pp; // Initial suction; suc = ua - uw
00293     if (_suc<0.0) _suc=0.0;
00294 
00295     // Calculate internal variables
00296     //REAL p = (_sig(0)+_sig(1)+_sig(2))/3.0;
00297     _z    = 0;
00298 #ifdef SUBBAR_INIT1
00299     _z(3) = _suc+OSI;
00300     _z(1) = _suc;
00301     _force_consistency(); // Pass subloading surface on stress point (find z(0) and z(1))
00302     _z(2) = OCR*_z(0);    // Initial "size" of (normal) Cam-clay elipse
00303 #else
00304     _force_consistency(); // Pass subloading surface on stress point (find z(0) and z(1))
00305     _z(2) = OCR*_z(0);    // Initial "size" of (normal) Cam-clay elipse
00306     _z(3) = OSI+_z(1);    // Initial size of (normal) SI cap plane
00307 #endif
00308 
00309     // Check if subloading surface is inside normal surface
00310     if (_z(0)>_z(2)) throw new Fatal(_("SubBar::SubBar: Subloading surface (z0=%f) must be smaller than normal surface (z2=%f)"),_z(0),_z(2));
00311     if (_z(1)>_z(3)) throw new Fatal(_("SubBar::SubBar: Subloading surface (z1=%f) must be smaller than normal surface (z3=%f)"),_z(1),_z(3));
00312 
00313     // Constants
00314     _MAX_ERR_F = (1.0/100.0)*_pref*_pref; // Max error on F (Yield Function)
00315     _MIN_HP    = (1.0/100.0)*_pref*_pref; // Min value of hp
00316 
00317 #ifndef SUBBAR_INIT1
00318     // Check consistency
00319     REAL F = YieldFunc(_sig, _z, _suc);
00320     if (fabs(F)>_MAX_ERR_F) throw new Fatal(_("SubBar::SubBar: Stress point is outside yield surface (F_ini = %f)"),F);
00321     //if (fabs(F)>_MAX_ERR_F) std::cout << _("SubBar::SubBar: Stress point is outside yield surface (F_ini = ") << F << ")" << std::endl;
00322 #endif
00323 
00324 } // }}}
00325 
00326 inline void SubBar::TgStiffness(Tensor4 & D, Tensor2 & d) const // {{{
00327 {
00328     if (_num_dLam<0) // Elastic tangent stiffness
00329     {
00330         _calc_De_de(_v,_sig,_suc, D,d); // D<-De, d<-de
00331     }
00332     else // Elastoplastic tangent stiffness
00333     {
00334         // Gradients and Hardening
00335         Tensor2 r;    // plastic flow direction
00336         Tensor2 V;    // dF_dSig
00337         T_dFdz  y;    // dF_dz
00338         REAL    S;    // dF_dSuc
00339         HardMod HH;   // Hardening modulae
00340         REAL    hp;   // Plastic coeficient
00341         _calc_gradients_and_hardening(_v, _sig, _z, _suc, r, V, y, S, HH, hp);
00342 
00343         // De and de
00344         Tensor4 De;
00345         Tensor2 de;
00346         _calc_De_de(_v,_sig,_suc, De,de);
00347 
00348         // Phi and Dep
00349         REAL Phi = Tensors::Reduce(V,De,r)+hp;  // Phi = V:De:r - yk*HHk
00350         Tensors::GerX(-1.0/Phi,De,r,V,De,De,D); // D <- Dep = (-1/Phi) * (De:r)dy(V:De) + De
00351 
00352         // dep
00353         Tensors::Dot(De,r, d);                   // d <- De:r
00354         d = de - ((Tensors::Dot(V,de)+S)/Phi)*d; // d <- dep = de - ((V:de)/Phi)*De:r
00355     }
00356 } // }}}
00357 
00358 inline REAL SubBar::Sr() const // {{{
00359 {
00360     if (_suc <= 0.0) // if saturated
00361         return 1.0;
00362     else
00363         return 1.0/pow(1.0+pow(_suc/_a,_b),_c); // Degree of saturation
00364 } // }}}
00365 
00366 inline REAL SubBar::phi() const // {{{
00367 {
00368     return Sr();
00369 } // }}}
00370 
00371 inline REAL SubBar::n_times_dSr_dPp() const // {{{
00372 {
00373     if (_suc <= 0.0) // if saturated
00374         return 0.0;
00375     else
00376     {
00377         REAL dSr_dPp = -_b*_c*pow(_suc/_a,_b-1.0)/(_a*pow(1.0+pow(_suc/_a,_b),_c+1.0));
00378         REAL       n = (_v-1.0)/_v;
00379         return n*dSr_dPp;
00380     }
00381 } // }}}
00382 
00383 inline void SubBar::TgPermeability(Tensor2 & k) const // {{{
00384 {
00385     k = (_ks*pow(Sr(),_m))*I;
00386 } // }}}
00387 
00388 inline void SubBar::FlowVelocity(Tensor1 const & Grad, Tensor1 & Vel) const // {{{
00389 {
00390     Tensor2 K;
00391     TgPermeability(K); 
00392     Tensors::Dot(-K,Grad, Vel);
00393 } // }}}
00394 
00395 inline void SubBar::Actualize(Tensor2 const & DSig, REAL const & DPp, Tensor2 & DEps, REAL & DnSr) // {{{
00396 {
00397     // Initial _eps and nSr
00398     Tensor2 eps_ini = _eps;
00399     REAL    n_ini   = (_v-1.0)/_v;
00400     REAL    Sr_ini  = Sr();
00401     REAL    nSr_ini = n_ini * Sr_ini;
00402 
00403     // Increment of suction
00404     REAL DSuc = 0.0 - DPp;
00405     if ((_suc+DSuc)<0.0) DSuc=0.0-_suc;
00406 
00407     // Gradients and Hardening
00408     Tensor2 r;    // plastic flow direction
00409     Tensor2 V;    // dF_dSig
00410     T_dFdz  y;    // dF_dz
00411     REAL    S;    // dF_dSuc
00412     HardMod HH;   // Hardening modulae
00413     REAL    hp;   // Plastic coeficient
00414     _calc_gradients_and_hardening(_v, _sig, _z, _suc, r, V, y, S, HH, hp);
00415 
00416     // Plastic (Lagrange) multiplier
00417     REAL dLam = (Tensors::Dot(V,DSig) + S*DSuc)/hp; // dLam = (V:dsig + S*dsuc)/hp
00418 
00419     //std::cout << "dLam=" << dLam << std::endl;
00420 
00421     // Actualize
00422     if (dLam<0.0) // Unloading
00423     {
00424         //if (CoupledModel::_isc.Type()==IntegSchemesCtes::FE)
00425         //{
00426             //std::cout << "Actualize:: Unloading:: ForwardEuler (FE): n_div=" << _isc.FE_ndiv() << ", DSuc=" << DSuc << std::endl;
00427 
00428             // Forward-Euler - Actualize (Unloading)
00429             int n_div=_isc.FE_ndiv();
00430             Tensor2 dsig = DSig/n_div;
00431             REAL    dsuc = DSuc/n_div;
00432             for (int i=0; i<n_div; ++i)
00433             {
00434                 // Ce, ce and DEps
00435                 Tensor4 Ce; Tensor2 ce; _calc_Ce_ce(_v,_sig,_suc, Ce,ce);
00436                 Tensors::Dot(Ce,dsig, DEps); DEps += ce * dsuc; // dEps = Ce:dSig + ce*dSuc
00437                 
00438                 // Actualize (Elastic)
00439                 _sig += dsig;
00440                 _suc += dsuc;
00441                 _eps += DEps;
00442 
00443                 _z(1) += dsuc;
00444                 //if (_z(3)<_z(1)) _z(3)=_z(1); // TODO: Check this
00445                 _force_consistency();
00446             }
00447         //}
00448         //throw new Fatal(_("Actualize:: Unloading:: ModifiedEuler (ME): not yet . . ."));
00449 
00450         /*
00451         // Check consistency
00452         REAL F = YieldFunc(_sig, _z, _suc);
00453         if (fabs(F)>_MAX_ERR_F)
00454         {
00455             REAL p,q;
00456             Stress_p_q(_sig,p,q);
00457             REAL f0 = q*q - _Mcs*_Mcs * (p+_k*_suc) * (Calc_p0(_z,_suc)-p);
00458             REAL f1 = _suc - _z(1);
00459             String str; str.Printf(_("SubBar::Actualize: Stress point is not on yield surface (f0=%f, f1=%f, F=%f)"),f0,f1,F);
00460             std::cout << "" << str << "" << std::endl;
00461         }
00462         */
00463 
00464     }
00465     else // Loading
00466     {
00467         if (CoupledModel::_isc.Type()==IntegSchemesCtes::FE)
00468         {
00469             //std::cout << "Actualize:: Loading:: ForwardEuler (FE): n_div=" << _isc.FE_ndiv() << std::endl;
00470 
00471             // Forward-Euler - Actualize (Loading)
00472             int n_div=_isc.FE_ndiv();
00473             Tensor2 dsig;     dsig    = DSig/n_div;
00474             REAL    dsuc;     dsuc    = DSuc/n_div;
00475             AugVec  dsigsuc;  dsigsuc = dsig(0),dsig(1),dsig(2),dsig(3),dsig(4),dsig(5),dsuc;
00476             for (int i=0; i<n_div; ++i)
00477             {
00478                 IntVars dz;
00479                 AugVec sigsuc;  sigsuc = _sig(0),_sig(1),_sig(2),_sig(3),_sig(4),_sig(5),_suc;
00480                 _tg_deps_dz_given_dsig_and_dsuc(sigsuc, _eps, _z, dsigsuc, DEps, dz);
00481                 _sig += dsig;
00482                 _suc += dsuc;
00483                 _eps += DEps;
00484                 _z   += dz;
00485             }
00486         }
00487         else // suppose it is AutoME
00488         {
00489             // Augmented vector
00490             AugVec  sigsuc;   sigsuc = _sig(0),_sig(1),_sig(2),_sig(3),_sig(4),_sig(5),_suc;
00491             AugVec dsigsuc;  dsigsuc = DSig(0),DSig(1),DSig(2),DSig(3),DSig(4),DSig(5),DSuc;
00492 
00493             // TimeInteg
00494             AutoStepME<AugVec,Tensor2,IntVars,SubBar> TI(CoupledModel::_isc);
00495 
00496             // Update stress/suction (_sig,_suc), internal variables (_z) and strain (_eps) state
00497             int nss = TI.Solve(this,                                        // Address of this instance of GenericEP class
00498                                &SubBar::_tg_deps_dz_given_dsig_and_dsuc, // Pointer to a function that calculates {dEps,dz}=[[Cep cep];[bk Pk]]:{dSig,dSuc}
00499                                &SubBar::_local_error,                    // Pointer to a function that calculates local error on dEps
00500                                sigsuc, _eps, _z, dsigsuc);
00501             // Check num_sub_steps
00502             if (nss==-1)
00503                 throw new Fatal(_("SubBar::Actualize: (Loading) Number of max sub-steps (%d) was not sufficient for AutoStepME."), CoupledModel::_isc.ME_maxSS());
00504 
00505             // Set stress/suction state
00506             _sig = sigsuc(0),sigsuc(1),sigsuc(2),sigsuc(3),sigsuc(4),sigsuc(5);
00507             _suc = sigsuc(6);
00508         }
00509     }
00510 
00511     // Calculate the total (secant) increment of strain
00512     DEps = _eps - eps_ini;
00513 
00514     // Update internal state (secant)
00515     _v += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v
00516 
00517     // Final pore-pressure
00518     _ppr += DPp;
00519 
00520     // Final nSr
00521     REAL n_fin   = (_v-1.0)/_v;
00522     REAL Sr_fin  = Sr();
00523     REAL nSr_fin = n_fin * Sr_fin;
00524     DnSr = nSr_fin - nSr_ini; 
00525 
00526 } // }}}
00527 
00528 inline void SubBar::StressUpdate(Tensor2 const & DEps, REAL const & DPp, Tensor2 & DSig, REAL & DnSr) // {{{
00529 {
00530     // Initial _sig and nSr
00531     Tensor2 sig_ini = _sig;
00532     REAL    n_ini   = (_v-1.0)/_v;
00533     REAL    Sr_ini  = Sr();
00534     REAL    nSr_ini = n_ini * Sr_ini;
00535 
00536     // Increment of suction
00537     REAL DSuc = 0.0 - DPp;
00538     if ((_suc+DSuc)<0.0) DSuc=0.0-_suc;
00539 
00540     // Gradients and Hardening
00541     Tensor2 r;    // plastic flow direction
00542     Tensor2 V;    // dF_dSig
00543     T_dFdz  y;    // dF_dz
00544     REAL    S;    // dF_dSuc
00545     HardMod HH;   // Hardening modulae
00546     REAL    hp;   // Plastic coeficient
00547     _calc_gradients_and_hardening(_v, _sig, _z, _suc, r, V, y, S, HH, hp);
00548 
00549     // De and de
00550     Tensor4 De;
00551     Tensor2 de;
00552     _calc_De_de(_v,_sig,_suc, De,de);
00553 
00554     // Numerator of Plastic (Lagrange) multiplier
00555     _num_dLam = (Tensors::Reduce(V,De,DEps) + Tensors::Dot(V,de)*DSuc + S*DSuc);
00556 
00557     // Update
00558     if (_num_dLam<0.0) // Unloading
00559     {
00560         //if (CoupledModel::_isc.Type()==IntegSchemesCtes::FE)
00561         //{
00562             //std::cout << "StressUpdate:: Unloading:: ForwardEuler (FE): n_div=" << _isc.FE_ndiv() << ", DSuc=" << DSuc << std::endl;
00563 
00564             // Forward-Euler - StressUpdate (Unloading)
00565             int n_div=_isc.FE_ndiv();
00566             Tensor2 deps = DEps/n_div;
00567             REAL    dsuc = DSuc/n_div;
00568             for (int i=0; i<n_div; ++i)
00569             {
00570                 // De, de and DSig
00571                 _calc_De_de(_v,_sig,_suc, De,de);
00572                 Tensors::Dot(De,deps, DSig); DSig += de * dsuc; // dSig = De:dEps + de*dSuc
00573                 
00574                 // Update (Elastic)
00575                 _sig += DSig;
00576                 _suc += dsuc;
00577                 _eps += deps;
00578 
00579                 _z(1) += dsuc;
00580                 //if (_z(3)<_z(1)) _z(3)=_z(1); // TODO: Check this
00581                 _force_consistency();
00582             }
00583         //}
00584     }
00585     else // Loading
00586     {
00587         if (CoupledModel::_isc.Type()==IntegSchemesCtes::FE)
00588         {
00589             // Forward-Euler - Actualize (Loading)
00590             int n_div = _isc.FE_ndiv();
00591             Tensor2 deps = DEps/n_div;
00592             REAL    dsuc = DSuc/n_div;
00593             AugVec depssuc;  depssuc = deps(0),deps(1),deps(2),deps(3),deps(4),deps(5),dsuc;
00594             for (int i=0; i<n_div; ++i)
00595             {
00596                 IntVars dz;
00597                 AugVec epssuc;  epssuc = _eps(0),_eps(1),_eps(2),_eps(3),_eps(4),_eps(5),_suc;
00598                 _tg_dsig_dz_given_deps_and_dsuc(epssuc, _sig, _z, depssuc, DSig, dz);
00599                 _sig += DSig;
00600                 _suc += dsuc;
00601                 _eps += deps;
00602                 _z   += dz;
00603             }
00604         }
00605         else
00606         {
00607             // Augmented vector
00608             AugVec  epssuc;   epssuc = _eps(0),_eps(1),_eps(2),_eps(3),_eps(4),_eps(5),_suc;
00609             AugVec depssuc;  depssuc = DEps(0),DEps(1),DEps(2),DEps(3),DEps(4),DEps(5),DSuc;
00610 
00611             // TimeInteg
00612             AutoStepME<AugVec,Tensor2,IntVars,SubBar> TI(CoupledModel::_isc);
00613 
00614             // Update stress/suction (_sig,_suc), internal variables (_z) and strain (_eps) state
00615             int nss = TI.Solve(this,                                         // Address of this instance of GenericEP class
00616                                &SubBar::_tg_dsig_dz_given_deps_and_dsuc, // Pointer to a function that calculates {dSig,dz}=[[Dep dep];[Bk Qk]]:{dEps,dSuc}
00617                                &SubBar::_local_error,                    // Pointer to a function that calculates local error on dEps
00618                                epssuc, _sig, _z, depssuc);
00619             // Check num_sub_steps
00620             if (nss==-1)
00621                 throw new Fatal(_("SubBar::StressUpdate: (Loading) Number of max sub-steps (%d) was not sufficient for AutoStepME."), CoupledModel::_isc.ME_maxSS());
00622 
00623             // Set strain/suction state
00624             _eps = epssuc(0),epssuc(1),epssuc(2),epssuc(3),epssuc(4),epssuc(5);
00625             _suc = epssuc(6);
00626         }
00627     }
00628 
00629     // Calculate the total (secant) increment of stress
00630     DSig = _sig - sig_ini;
00631 
00632     // Update internal state (secant)
00633     _v += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v
00634 
00635     // Final pore-pressure
00636     _ppr += DPp;
00637 
00638     // Final nSr
00639     REAL n_fin   = (_v-1.0)/_v;
00640     REAL Sr_fin  = Sr();
00641     REAL nSr_fin = n_fin * Sr_fin;
00642     DnSr = nSr_fin - nSr_ini; 
00643 
00644 } // }}}
00645 
00646 inline void SubBar::BackupState() // {{{
00647 {
00648     _bkp_v   = _v;
00649     _bkp_sig = _sig;
00650     _bkp_suc = _suc;
00651     _bkp_eps = _eps;
00652     _bkp_z   = _z;
00653     _bkp_ppr = _ppr;
00654     _bkp_num_dLam = _num_dLam;
00655 } // }}}
00656 
00657 inline void SubBar::RestoreState() // {{{
00658 {
00659       _v = _bkp_v;
00660     _sig = _bkp_sig;
00661     _suc = _bkp_suc;
00662     _eps = _bkp_eps;
00663       _z = _bkp_z;
00664     _ppr = _bkp_ppr;
00665     _num_dLam = _bkp_num_dLam;
00666 } // }}}
00667 
00668 inline REAL SubBar::_lam(REAL const & Suc) const // {{{
00669 {
00670     return _lam0*((1.0-_r)*exp(-_bet*Suc) + _r);
00671 } // }}}
00672 
00673 inline REAL SubBar::_psi(REAL const & Suc) const // {{{
00674 {
00675     return (_lam0-_kap)/(_lam(Suc)-_kap);
00676 } // }}}
00677 
00678 inline REAL SubBar::Calc_p0(IntVars const & z, REAL const & Suc) const // {{{
00679 {
00680     return _pref*pow(z(0)/_pref, _psi(Suc));
00681 } // }}}
00682 
00683 inline REAL SubBar::Calc_C(IntVars const & z, REAL const & Suc) const // {{{
00684 {
00685     return (_pref*_pref) * ( exp(_BB*(Suc-z(1))/_pref) - exp(-_BB*z(1)/_pref) );
00686 } // }}}
00687 
00688 inline REAL SubBar::YieldFunc(Tensor2 const & Sig, IntVars const & z, REAL const & Suc) const // {{{
00689 {
00690     REAL     p,q,sin3th;
00691     Tensor2  devS;
00692     Stress_p_q_S_sin3th(Sig, p,q,devS,sin3th);
00693     REAL M = _calc_M(sin3th);
00694     REAL g = - M*M * (p+_k*Suc) * (Calc_p0(z,Suc)-p);
00695     REAL F = q*q + g + Calc_C(z,Suc);
00696     //if (fabs(F)<=_MAX_ERR_F) F=0.0;
00697     return F;
00698 } // }}}
00699 
00700 inline REAL SubBar::YieldFunc_Normal(Tensor2 const & Sig, IntVars const & z, REAL const & Suc) const // {{{
00701 {
00702     IntVars z_tmp; z_tmp=z(2),z(3),0,0;
00703     REAL     p,q,sin3th;
00704     Tensor2  devS;
00705     Stress_p_q_S_sin3th(Sig, p,q,devS,sin3th);
00706     REAL M = _calc_M(sin3th);
00707     REAL g = - M*M * (p+_k*Suc) * (Calc_p0(z_tmp,Suc)-p);
00708     REAL F = q*q + g + Calc_C(z_tmp,Suc);
00709     //if (fabs(F)<=_MAX_ERR_F) F=0.0;
00710     return F;
00711 } // }}}
00712 
00713 inline void SubBar::_calc_Ce_ce(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & Ce, Tensor2 & ce) const // {{{
00714 {
00715     // Invariants
00716     REAL p = (Sig(0)+Sig(1)+Sig(2))/3.0;
00717 
00718     // Calculate Bulk modulae
00719     REAL K  =  p         *v/_kap;
00720     REAL Ks = (Suc+_patm)*v/_kaps;
00721 
00722     // Calc Ce
00723     Tensors::AddScaled(1.0/(2.0*_G), Tensors::Psd, 1.0/(9.0*K), Tensors::IdyI,  Ce);
00724 
00725     // Calc ce
00726     ce = (1.0/(3.0*Ks)) * I;
00727 
00728 } // }}}
00729 
00730 inline void SubBar::_calc_De_de(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & De, Tensor2 & de) const // {{{
00731 {
00732     // Invariants
00733     REAL p = (Sig(0)+Sig(1)+Sig(2))/3.0;
00734 
00735     // Calculate Bulk modulae
00736     REAL K  =  p         *v/_kap;
00737     REAL Ks = (Suc+_patm)*v/_kaps;
00738 
00739     // Calc De
00740     Tensors::AddScaled(2.0*_G,Tensors::Psd, K,Tensors::IdyI, De);
00741 
00742     // Calc de
00743     de = (-K/Ks)*I;
00744 
00745 } // }}}
00746 
00747 inline void SubBar::Calc_dFdSig_dFdSuc(Tensor2 const & Sig, IntVars const & z, REAL const & Suc, Tensor2 & V, REAL & S) const // {{{
00748 {
00749     // Invariants
00750     REAL     p,q,sin3th;
00751     Tensor2  devS;
00752     Stress_p_q_S_sin3th(Sig, p,q,devS,sin3th);
00753 
00754     // Suction-related variables
00755     REAL lam = _lam(Suc);
00756     REAL psi = _psi(Suc);
00757     REAL  p0 = Calc_p0(z,Suc);
00758     REAL  ps = _k*Suc;
00759 
00760     // Gradients
00761     REAL        M = _calc_M(sin3th);
00762     REAL       MM = M*M;
00763     REAL      tmp = MM*(2.0*p-p0+ps)/3.0;
00764     REAL dp0_dSuc = p0*log(z(0)/_pref)*psi*_bet*_lam0*(1.0-_r)*exp(-_bet*Suc)/(lam-_kap);
00765     REAL  dC_dSuc = _BB*_pref*exp(_BB*(Suc-z(1))/_pref);
00766     REAL  dg_dSuc = MM*((p-p0)*_k - (p+ps)*dp0_dSuc);
00767                 V = tmp*I + 3.0*devS;  // dF_dSig
00768                 S = dg_dSuc + dC_dSuc; // dF_dSuc
00769     if (q>Q_ZERO)
00770     {
00771         REAL  ss3th = pow(sin3th,2.0);
00772         REAL cos3th = (ss3th>1.0 ? 0.0 : sqrt(1.0-ss3th));
00773         if (cos3th>Q_ZERO) // if q is Nan of Inf, prs may become a lot slower
00774         {
00775             REAL                dM_dth = (0.75*M*(1.0-_wcs)*cos3th) / (1.0+_wcs+(_wcs-1.0)*sin3th);
00776             Tensor2 SS;         Mult(devS,devS,SS);
00777             Tensor2 dev_SS;     dev_SS = SS - (I * (SS(0)+SS(1)+SS(2))/3.0);
00778             Tensor2 dth_dsig; dth_dsig = (1.5/(q*q*cos3th)) * (dev_SS*(3.0/q) - devS*sin3th);
00779                                     V += dth_dsig*(2.0*M*(p-p0)*(p+ps)*dM_dth);
00780         }
00781     }
00782 
00783 } // }}}
00784 
00785 inline void SubBar::_calc_f0_grads(Tensor2 const & Sig     , // {{{
00786                                    IntVars const & z       ,
00787                                    REAL    const & Suc     ,
00788                                    Tensor2       & df0_dSig,
00789                                    REAL          & df0_dSuc, 
00790                                    REAL          & df0_dz0 ) const
00791 {
00792     // Invariants
00793     REAL     p,q,sin3th;
00794     Tensor2  devS;
00795     Stress_p_q_S_sin3th(Sig, p,q,devS,sin3th);
00796 
00797     // Suction-related variables
00798     REAL lam = _lam(Suc);
00799     REAL psi = _psi(Suc);
00800     REAL  p0 = Calc_p0(z,Suc);
00801     REAL  ps = _k*Suc;
00802 
00803     // Gradients
00804     REAL        M = _calc_M(sin3th);
00805     REAL       MM = M*M;
00806     REAL dp0_dSuc = p0*log(z(0)/_pref)*psi*_bet*_lam0*(1.0-_r)*exp(-_bet*Suc)/(lam-_kap);
00807          df0_dSuc = MM*((p-p0)*_k - (p+ps)*dp0_dSuc);
00808          df0_dSig = (MM*(2.0*p-p0+ps)/3.0)*I + 3.0*devS;
00809           df0_dz0 = -MM*psi*p0*(p+ps)/z(0);
00810     
00811 } // }}}
00812 
00813 inline void SubBar::_calc_gradients_and_hardening(REAL    const & v, // {{{
00814                                                   Tensor2 const & Sig,
00815                                                   IntVars const & z,
00816                                                   REAL    const & Suc,
00817                                                   Tensor2       & r,
00818                                                   Tensor2       & V,
00819                                                   T_dFdz        & y, 
00820                                                   REAL          & S,
00821                                                   HardMod       & HH,
00822                                                   REAL          & hp) const
00823 {
00824     // Invariants
00825     REAL     p,q,sin3th;
00826     Tensor2  devS;
00827     Stress_p_q_S_sin3th(Sig, p,q,devS,sin3th);
00828 
00829     // Suction-related variables
00830     REAL lam = _lam(Suc);
00831     REAL psi = _psi(Suc);
00832     REAL  p0 = Calc_p0(z,Suc);
00833     REAL  ps = _k*Suc;
00834 
00835     //std::cout << "p0=" << p0 << std::endl;
00836     
00837     // Gradients
00838     REAL        M = _calc_M(sin3th);
00839     REAL       MM = M*M;
00840     REAL      tmp = MM*(2.0*p-p0+ps)/3.0;
00841     REAL dp0_dSuc = p0*log(z(0)/_pref)*psi*_bet*_lam0*(1.0-_r)*exp(-_bet*Suc)/(lam-_kap);
00842     REAL  dC_dSuc = _BB*_pref*exp(_BB*(Suc-z(1))/_pref);
00843     REAL  dg_dSuc = MM*((p-p0)*_k - (p+ps)*dp0_dSuc);
00844                 r = tmp*I + (3.0*_alp)*devS;  // flow rule
00845                 V = tmp*I +  3.0      *devS;  // dF_dSig
00846                 S = dg_dSuc + dC_dSuc;        // dF_dSuc
00847     if (q>Q_ZERO)
00848     {
00849         REAL  ss3th = pow(sin3th,2.0);
00850         REAL cos3th = (ss3th>1.0 ? 0.0 : sqrt(1.0-ss3th));
00851         if (cos3th>Q_ZERO) // if q is Nan of Inf, prs may become a lot slower
00852         {
00853             REAL                dM_dth = (0.75*M*(1.0-_wcs)*cos3th) / (1.0+_wcs+(_wcs-1.0)*sin3th);
00854             Tensor2 SS;         Mult(devS,devS,SS);
00855             Tensor2 dev_SS;     dev_SS = SS - (I * (SS(0)+SS(1)+SS(2))/3.0);
00856             Tensor2 dth_dsig; dth_dsig = (1.5/(q*q*cos3th)) * (dev_SS*(3.0/q) - devS*sin3th);
00857                                     V += dth_dsig*(2.0*M*(p-p0)*(p+ps)*dM_dth);
00858                                     r += dth_dsig*(2.0*M*(p-p0)*(p+ps)*dM_dth);
00859         }
00860     }
00861     y(0) = -MM*psi*p0*(p+ps)/z(0);
00862     y(1) = -_BB*Calc_C(z,Suc)/_pref;
00863     y(2) = 0.0;
00864     y(3) = 0.0;
00865 
00866     // Hardening
00867     REAL  tr_r = r(0)+r(1)+r(2);
00868     REAL  chi0 = (_lam0 - _kap )/v;
00869     REAL  chis = (_lams - _kaps)/v;
00870     REAL     L = _c0*pow(chi0*log( z(2)       / z(0)       ),2.0)/p;
00871     REAL    Ls = _cs*pow(chis*log((z(3)+_patm)/(z(1)+_patm)),2.0)/(Suc+_patm);
00872          HH(0) =  z(0)       *(tr_r+L )/chi0;
00873          HH(1) = (z(1)+_patm)*(tr_r+Ls)/chis;
00874          HH(2) =  z(2)       * tr_r    /chi0;
00875          HH(3) = (z(3)+_patm)* tr_r    /chis;
00876 
00877     // Plastic coeficient
00878     hp = -(y(0)*HH(0)+y(1)*HH(1));
00879 
00880     //std::cout << "z=" << z << ", L=" << L << ", Ls=" << Ls << ", hp=" << hp << std::endl;
00881     
00882 } // }}}
00883 
00884 inline int SubBar::_tg_deps_dz_given_dsig_and_dsuc(AugVec  const & SigSuc, // {{{
00885                                                    Tensor2 const & Eps,
00886                                                    IntVars const & z,
00887                                                    AugVec  const & dSigSuc,
00888                                                    Tensor2       & dEps,
00889                                                    IntVars       & dz) const
00890 {
00891     // Stop if excessive deviatoric strains
00892     REAL ev,ed;
00893     Tensors::Strain_Ev_Ed(Eps, ev,ed);
00894     if (ed>CoupledModel::_isc.EdMax()) // return 1; // -1 => failed;  0 => continue;  1 => stop
00895         throw new Warning(_("SubBar::_tg_deps_dz_given_dsig_and_dsuc: Max devaitoric strain Ed=%f greater than EdMax=%f"),ed,_isc.EdMax());
00896 
00897     // Input values
00898     Tensor2  sig;   sig =  SigSuc(0), SigSuc(1), SigSuc(2), SigSuc(3), SigSuc(4), SigSuc(5);
00899     Tensor2 dsig;  dsig = dSigSuc(0),dSigSuc(1),dSigSuc(2),dSigSuc(3),dSigSuc(4),dSigSuc(5);
00900     REAL     suc;   suc =  SigSuc(6);
00901     REAL    dsuc;  dsuc = dSigSuc(6);
00902     
00903     // Gradients and Hardening
00904     REAL    v=_v; // (constant) current specific volume
00905     Tensor2 r;    // plastic flow direction
00906     Tensor2 V;    // dF_dSig
00907     T_dFdz  y;    // dF_dz
00908     REAL    S;    // dF_dSuc
00909     HardMod HH;   // Hardening modulae
00910     REAL    hp;   // Plastic coeficient
00911     _calc_gradients_and_hardening(v, sig, z, suc, r, V, y, S, HH, hp);
00912 
00913     // Try to check for softening initial point
00914     if (fabs(hp/(sig(0)+sig(1)+sig(2)))<=_MIN_HP)
00915         throw new Warning(_("SubBar::_tg_deps_dz_given_dsig_and_dsuc: hp=%f (zero). Maybe it is the start of softening behaviour"),hp);
00916 
00917     // Ce and ce
00918     Tensor4 Cep; // Cep <- Ce
00919     Tensor2 ce;
00920     _calc_Ce_ce(v,sig,suc, Cep,ce); // Cep <- Ce
00921 
00922     // Cep
00923     Tensors::Ger(1.0/hp,r,V, Cep); // Cep <- (1/hp)*(r dyadic V) + Ce
00924 
00925     // cep
00926     Tensor2 cep;
00927     cep = (S/hp)*r + ce;
00928 
00929     // dEps
00930     Tensors::Dot(Cep,dsig, dEps);
00931     dEps += cep*dsuc;            // dEps = Cep:dsig + cep*dsuc
00932 
00933     // dLam and dz
00934     REAL dLam = (Tensors::Dot(V,dsig) + S*dsuc)/hp; // dLam = (V:dsig + S*dsuc)/hp
00935            dz = dLam*HH;
00936 
00937     return 0;
00938 
00939 } // }}}
00940 
00941 inline int SubBar::_tg_dsig_dz_given_deps_and_dsuc(AugVec  const & EpsSuc, // {{{
00942                                                    Tensor2 const & Sig,
00943                                                    IntVars const & z,
00944                                                    AugVec  const & dEpsSuc,
00945                                                    Tensor2       & dSig,
00946                                                    IntVars       & dz) const
00947 {
00948     // Input values
00949     Tensor2  eps;   eps =  EpsSuc(0), EpsSuc(1), EpsSuc(2), EpsSuc(3), EpsSuc(4), EpsSuc(5);
00950     Tensor2 deps;  deps = dEpsSuc(0),dEpsSuc(1),dEpsSuc(2),dEpsSuc(3),dEpsSuc(4),dEpsSuc(5);
00951     REAL     suc;   suc =  EpsSuc(6);
00952     REAL    dsuc;  dsuc = dEpsSuc(6);
00953     
00954     // Gradients and Hardening
00955     REAL    v=_v; // (constant) current specific volume
00956     Tensor2 r;    // plastic flow direction
00957     Tensor2 V;    // dF_dSig
00958     T_dFdz  y;    // dF_dz
00959     REAL    S;    // dF_dSuc
00960     HardMod HH;   // Hardening modulae
00961     REAL    hp;   // Plastic coeficient
00962     _calc_gradients_and_hardening(v, Sig, z, suc, r, V, y, S, HH, hp);
00963 
00964     // De and de
00965     Tensor4 De;
00966     Tensor2 de;
00967     _calc_De_de(v,Sig,suc, De,de);
00968 
00969     // Phi and Dep
00970     Tensor4 Dep;
00971     REAL Phi = Tensors::Reduce(V,De,r)+hp;    // Phi = V:De:r - yk*HHk
00972     Tensors::GerX(-1.0/Phi,De,r,V,De,De,Dep); // Dep = (-1/Phi) * (De:r)dy(V:De) + De
00973 
00974     // dep
00975     Tensor2 dep;
00976     Tensors::Dot(De,r, dep);                     // dep <- De:r
00977     dep = de - ((Tensors::Dot(V,de)+S)/Phi)*dep; // dep <- de - ((V:de+S)/Phi)*De:r
00978 
00979     // dSig
00980     Tensors::Dot(Dep,deps, dSig); // dSig <- Dep:deps
00981     dSig += dep*dsuc;             // dSig <- Dep:deps + dep*dsuc
00982 
00983     // dLam and dz
00984     REAL dLam = (Tensors::Reduce(V,De,deps) + Tensors::Dot(V,de)*dsuc + S*dsuc)/Phi; // dLam = (V:De:deps + V:de*dsuc + S*dsuc)/Phi
00985            dz = dLam*HH;
00986 
00987     return 0;
00988 
00989 } // }}}
00990 
00991 inline REAL SubBar::_local_error(Tensor2 const & Ey, Tensor2 const & y_high, // {{{
00992                                      IntVars const & Ez, IntVars const & z_high) const
00993 {
00994     REAL Err_y = Tensors::Norm(Ey)/Tensors::Norm(y_high);
00995     REAL Err_z = sqrt(blitz::dot(Ez,Ez))/sqrt(blitz::dot(z_high,z_high));
00996     return (Err_y>Err_z ? Err_y : Err_z);
00997 } // }}}
00998 
00999 inline void SubBar::_force_consistency() // {{{
01000 {
01001     // Invariants
01002     REAL     p,q,sin3th;
01003     Tensor2  devS;
01004     Stress_p_q_S_sin3th(_sig, p,q,devS,sin3th);
01005 
01006     // Auxiliar variables
01007     REAL   M = _calc_M(sin3th);
01008     REAL  MM = M*M;
01009     REAL psi = _psi(_suc);
01010     REAL  p0 = p + (q*q)/(MM*(p+_k*_suc));
01011 
01012     // Internal values
01013     _z(0) = _pref*pow(p0/_pref,1.0/psi);
01014 
01015 #ifndef SUBBAR_INIT1
01016     // Find z(1) (Newton-Raphson)
01017           int          k = 0;
01018     const int      maxIt = 20;
01019     const REAL resid_TOL = 1.0e-7;
01020                    _z(1) = _suc; // Trial
01021           REAL         g = -MM*(p+_k*_suc)*(p0-p);
01022           REAL    C_star = -q*q - g;
01023     while (k<=maxIt)
01024     {
01025         REAL    Ck = Calc_C(_z, _suc);
01026         REAL resid = C_star - Ck;
01027         if (fabs(resid)<=resid_TOL)
01028         {
01029             //std::cout << "SubBar::_force_consistency: k=" << k << ", resid=" << resid << ", z=" << _z << std::endl;
01030             return;
01031         }
01032         REAL dresid_dz1 = _BB*Ck/_pref;
01033         _z(1) = _z(1) - resid/dresid_dz1;
01034         k++;
01035     }
01036     throw new Fatal (_("SubBar::_force_consistency Newton-Raphson failed"));
01037 #endif
01038 
01039 } // }}}
01040 
01041 inline void SubBar::InternalStateValues(Array<REAL> & IntStateVals) const // {{{
01042 {
01043     IntStateVals.resize(5); // 5 = 4 nIntVars plus 1 for the specific volume component
01044     IntStateVals[0] = _z(0);
01045     IntStateVals[1] = _z(1);
01046     IntStateVals[2] = _z(2);
01047     IntStateVals[3] = _z(3);
01048     IntStateVals[4] = _v;
01049 } // }}}
01050 
01051 inline void SubBar::InternalStateNames(Array<String> & IntStateNames) const // {{{
01052 {
01053     IntStateNames.resize(5); // 5 = 4 nIntVars plus 1 for the specific volume component
01054     IntStateNames[0] = "z0";
01055     IntStateNames[1] = "z1";
01056     IntStateNames[2] = "z2";
01057     IntStateNames[3] = "z3";
01058     IntStateNames[4] = "v";
01059 } // }}}
01060 
01061 inline REAL SubBar::_calc_M(REAL const & sin3th) const // {{{
01062 {
01063     return _Mcs*pow( 2.0*_wcs/(1.0+_wcs+(_wcs-1.0)*sin3th) ,0.25);
01064 } // }}}
01065 
01066 
01068 
01069 
01070 // Allocate a new SubBar model
01071 CoupledModel * SubBarMaker(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
01072 {
01073     return new SubBar(Prms, IniData);
01074 } // }}}
01075 
01076 // Register an SubBar model into ModelFactory array map
01077 int SubBarRegister() // {{{
01078 {
01079     CoupledModelFactory["SubBar"] = SubBarMaker;
01080     return 0;
01081 } // }}}
01082 
01083 // Execute the autoregistration
01084 int __SubBar_dummy_int = SubBarRegister();
01085 
01086 #endif // MECHSYS_SUBBAR_H
01087 
01088 // vim:fdm=marker
Generated on Wed Jan 24 15:56:27 2007 for MechSys by  doxygen 1.4.7