MechSys: models/barcelona.h Source File

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_BARCELONA_H
00023 #define MECHSYS_BARCELONA_H
00024 
00025 // STL
00026 #include <iostream>
00027 #include <blitz/tinyvec-et.h>
00028 
00029 // MechSys
00030 #include "models/coupled/coupledmodel.h"
00031 #include "tensors/tensors.h"
00032 #include "tensors/operators.h"
00033 #include "tensors/functions.h"
00034 #include "numerical/autostepme.h"
00035 #include "util/array.h"
00036 #include "util/util.h"
00037 
00038 using Tensors::I;
00039 using Tensors::Tensor2;
00040 using Tensors::Tensor4;
00041 using Tensors::Stress_p_q;
00042 using Tensors::Strain_Ev_Ed;
00043 
00045 class Barcelona : public CoupledModel
00046 {
00047 public:
00048 
00050     typedef blitz::TinyVector<REAL,2> IntVars; // {z0,z1} == {p0star,s0}
00051     typedef blitz::TinyVector<REAL,2> HardMod; // {HH0,HH1}
00052     typedef blitz::TinyVector<REAL,2> T_dFdz;  // {y0,y1} derivative of f w.r.t. intvars
00053     typedef blitz::TinyVector<REAL,7> AugVec;  // {Sig/Eps,Suc} Augmented tensor with suction
00054 
00055     // Constructor
00056     Barcelona(Array<REAL> const & Prms, Array<REAL> const & IniData);
00057 
00058     // Methods
00059     String Name            () const { return "Barcelona"; }
00060     void   TgStiffness     (Tensor4 & D, Tensor2 & d) const {}
00061     REAL   Sr              () const { return 0; }
00062     REAL   phi             () const { return 0; }
00063     REAL   n_times_dSr_dPp () const { return 0; }
00064     void   TgPermeability  (Tensor2 & k) const {} // TODO: write this method 
00065     void   FlowVelocity    (Tensor1 const & Grad, Tensor1 & Vel) const {};
00066     void   FlowUpdate      (REAL const & DPp) {};
00067     void   Actualize       (Tensor2 const & DSig, REAL const & DPp, Tensor2 & DEps, REAL & DnSr);
00068     void   StressUpdate    (Tensor2 const & DEps, REAL const & DPp, Tensor2 & DSig, REAL & DnSr) {}
00069     void   BackupState     (){}
00070     void   RestoreState    (){}
00071 
00072     // Methods
00073     Tensor2 const & Sig                  () const { return _sig;       }
00074     Tensor2 const & Eps                  () const { return _eps;       }
00075     REAL            Pp                   () const { return 0.0 - _suc; }
00076     REAL            GammaW               () const { return 0.0;        } // TODO: write this method
00077     int             nInternalStateValues () const { return 3;          }
00078     void            InternalStateValues  (Array<REAL>   & IntStateVals ) const;
00079     void            InternalStateNames   (Array<String> & IntStateNames) const;
00080     
00081 private:
00082     // State vars
00083     REAL    _v;   
00084     Tensor2 _sig; 
00085     REAL    _suc; 
00086     Tensor2 _eps; 
00087     IntVars _z;   
00088     
00089     // Parameters
00090     REAL _lam0;
00091     REAL _kap;
00092     REAL _G;
00093     REAL _phics;
00094     REAL _r;
00095     REAL _bet;
00096     REAL _k;
00097     REAL _lams;
00098     REAL _kaps;
00099     REAL _patm;
00100     REAL _pref;
00101     REAL _A;
00102     REAL _B;
00103 
00104     // Derived constants
00105     REAL _Mcs;
00106     REAL _alp;
00107 
00108     mutable REAL _EDMAX;
00109             REAL _MAX_ERR_F;
00110 
00111     // Private methods
00112     REAL _lam        (REAL const & Suc) const;
00113     REAL _psi        (REAL const & Suc) const;
00114     REAL _p0         (REAL const & Suc, REAL const & z0) const;
00115     REAL _yield_func (Tensor2 const & Sig, IntVars const & z, REAL const & Suc) const;
00116 
00117     void _calc_Ce_ce(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & Ce, Tensor2 & ce) const;
00118 
00119     void _calc_dF_dSig_dF_dSuc(Tensor2 const & Sig, IntVars const & z, REAL const & Suc, Tensor2 & V, REAL & S) const;
00120 
00121     void _calc_gradients_and_hardening(REAL    const & v  ,        // Specific volume
00122                                        Tensor2 const & Sig,        // Stress
00123                                        IntVars const & z  ,        // Internal variables
00124                                        REAL    const & Suc,        // Suction
00125                                        Tensor2       & r  ,        // Plastic flow direction
00126                                        Tensor2       & V  ,        // dF_dSig
00127                                        T_dFdz        & y  ,        // dF_dz
00128                                        REAL          & S  ,        // dF_dSuc
00129                                        HardMod       & HH ,        // Hardening
00130                                        REAL          & hp ) const; // Plastic coeficient
00131 
00132     int _tg_deps_dz_given_dsig_and_dsuc(AugVec  const & SigSuc ,
00133                                         Tensor2 const & Eps    ,
00134                                         IntVars const & z      ,
00135                                         AugVec  const & dSigSuc,
00136                                         Tensor2       & dEps   ,
00137                                         IntVars       & dz     ) const;
00138 
00139     REAL _local_error(Tensor2 const & Ey, Tensor2 const & y_high,
00140                       IntVars const & Ez, IntVars const & z_high) const;
00141 
00142     // Find intersection Newton-Raphson. Start from the actual state {_sig,_suc}. Returns Alpha.
00143     REAL _find_intersection_NR(Tensor2 const & dSig_tr,
00144                                REAL    const & dSuc_tr,
00145                                REAL    const & F      ,
00146                                REAL    const & F_tr   ,
00147                                Tensor2       & Sig_int,
00148                                REAL          & Suc_int) const;
00149 
00150     // Find intersection Brent's method
00151     REAL _find_intersection(Tensor2 const & dSig_tr,
00152                             REAL    const & dSuc_tr,
00153                             REAL    const & F      ,
00154                             REAL    const & F_tr   ,
00155                             Tensor2       & Sig_int,
00156                             REAL          & Suc_int) const;
00157 
00158 }; // class Barcelona
00159 
00160 
00162 
00163 
00164 inline Barcelona::Barcelona(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
00165 {
00166 
00167     // ##################################################################### Setup Parameters
00168 
00169     /* example.mat
00170        #------------------------------------------------------------------------------------------
00171        #             0     1     2       3     4     5    6     7      8      9     10    11    12
00172        #          lam0   kap     G   phics     r   bet    k  lams   kaps   patm   pref     A     B
00173        Barcelona  0.20  0.02  10.0  25.377  0.75  12.5  0.6  0.08  0.008  0.101  0.101  1000  1000  
00174     */
00175 
00176     // Check number of parameters
00177     const size_t n_prms = 13;
00178     if (Prms.size()!=n_prms)
00179         throw new Fatal(_("Barcelona::Barcelona: The number of parameters (%d) is incorrect. It must be equal to %d"), Prms.size(), n_prms);
00180 
00181     // Parameters
00182     _lam0  = Prms[0];
00183     _kap   = Prms[1];
00184     _G     = Prms[2];
00185     _phics = Prms[3];
00186     _r     = Prms[4];
00187     _bet   = Prms[5];
00188     _k     = Prms[6];
00189     _lams  = Prms[7];
00190     _kaps  = Prms[8];
00191     _patm  = Prms[9];
00192     _pref  = Prms[10];
00193     _A     = Prms[11];
00194     _B     = Prms[12];
00195 
00196     // Derived constants
00197     REAL sin_phi = sin(Util::ToRad(_phics));
00198     _Mcs = 6.0*sin_phi/(3.0-sin_phi);
00199     _alp = (_Mcs*(_Mcs-9.0)*(_Mcs-3.0)) / (9.0*(6.0-_Mcs)*(1.0-_kap/_lam0));
00200     
00201     // ##################################################################### Setup Initial State
00202 
00203     // Check number of initial data
00204     if (IniData.size()!=7) // Layered analysis
00205         throw new Fatal(_("Barcelona::Barcelona: The number of initial data is not sufficiet (it must be equal to 5). { SigX SigY SigZ vini OCR z1 Pp }"));
00206 
00207     // Read initial data
00208     REAL sig_x = IniData[0]; // X stress component
00209     REAL sig_y = IniData[1]; // Y stress component
00210     REAL sig_z = IniData[2]; // Z stress component
00211             _v = IniData[3]; // Initial specific volume
00212     REAL   OCR = IniData[4]; // Over Consolidation Ratio
00213     REAL    z1 = IniData[5]; // Initial "position" of surface for suction
00214     REAL    Pp = IniData[6]; // Initial pore-pressure
00215 
00216     // Fill stress and strain tensors
00217     REAL sq2=Util::Sq2();
00218     _sig = sig_x, sig_y, sig_z, 0.0*sq2, 0.0*sq2, 0.0*sq2; // Mandel representation
00219     _eps = 0.0,0.0,0.0,0.0*sq2,0.0*sq2,0.0*sq2;
00220 
00221     // Calculate initial state
00222     _suc = 0.0 - Pp; // Initial suction; suc = ua - uw
00223 
00224     // Calculate internal variables
00225     REAL p,q;
00226     Stress_p_q(_sig,p,q);
00227     _z(0) = OCR*p; // Initial "size" of Cam-clay elipse
00228     _z(1) = z1;    // Initial "position" of surface for suction
00229 
00230     // Check consistency
00231     const REAL MAX_F=1.0e-5;
00232     REAL F = _yield_func(_sig, _z, _suc);
00233     if (F>MAX_F) throw new Fatal(_("Barcelona::Barcelona: Stress point is outside of yield surface (F_ini = %f)"),F);
00234 
00235     // Constants
00236     _MAX_ERR_F = (0.5/100.0)*_pref;
00237 
00238 } // }}}
00239 
00240 inline void Barcelona::Actualize(Tensor2 const & DSig, REAL const & DPp, Tensor2 & DEps, REAL & DnSr) // {{{
00241 {
00242     // Trial state
00243     Tensor2 dsig_tr = DSig;
00244     Tensor2  sig_tr = _sig + dsig_tr;
00245     REAL    dsuc_tr = 0.0 - DPp;
00246     REAL     suc_tr = _suc + dsuc_tr;
00247 
00248     // Verify if suc_tr crosses SI (Suction-Increase) surface
00249     if (suc_tr>_z(1)) // Actualize model considering that DSig is null
00250     {
00251         if (_z(1)<0.0) throw new Fatal("Barcelona::Actualize: ERROR: _z(1)<0.0");
00252 
00253         // Check if DSig is null
00254         const REAL DSIGTOL = 1.0e-7;
00255         if (Tensors::Norm(dsig_tr)>DSIGTOL)
00256             throw new Fatal("Barcelona::Actualize: When suc_tr (%f) crosses SI surface, Norm(dsig_tr) (%f) must be null",suc_tr,Tensors::Norm(dsig_tr));
00257 
00258         // Intersection
00259         REAL suc_int = _z(1);
00260 
00261         // Elastic and Elastoplastic parts of suction
00262         REAL  dsuc_el = suc_int - _suc;
00263         REAL dsuc_epl = dsuc_tr - dsuc_el;
00264 
00265         // Actualize (Elastic)
00266         REAL deve = _kaps*dsuc_el/(_v*(_z(1)+_patm));
00267              DEps = (deve/3.0)*I;
00268             _v   += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v
00269             _suc += dsuc_el;
00270             _eps += DEps;
00271 
00272         // Forward-Euler - Actualize (Elastoplastic)
00273         int n_div=10;
00274         REAL dsuc = dsuc_epl/n_div;
00275         for (int i=0; i<n_div; ++i)
00276         {
00277             // FE Increments
00278             IntVars dz;
00279             REAL deve =        _kaps *dsuc/(_v*(_z(1)+_patm));
00280             REAL devp = (_lams-_kaps)*dsuc/(_v*(_z(1)+_patm));
00281             REAL  dev = deve + devp;
00282 
00283             // Actualize
00284             dz(0) = _v* _z(0)       *devp/(_lam0-_kap );
00285             dz(1) = _v*(_z(1)+_patm)*devp/(_lams-_kaps);
00286              DEps = (dev/3.0)*I;
00287             _v   += -(DEps(0)+DEps(1)+DEps(2))*_v;
00288             _suc += dsuc;
00289             _eps += DEps;
00290             _z   += dz;
00291         }
00292 
00293         return;
00294     }
00295 
00296     // Check loading condition for suction
00297     enum Situation {EL, EL_EPL, EPL}; // EL:     Elastic-Loading
00298                                       // EL_EPL: Elastic-Loading + ElastoPlastic-Loading
00299                                       // EPL:    ElastoPlastic-Loading
00300     REAL F    = _yield_func(_sig   ,_z,_suc   );
00301     REAL F_tr = _yield_func( sig_tr,_z, suc_tr);
00302     Situation sit;
00303     if (F<0.0) // current state is inside the surface
00304     {
00305         if (F_tr<=0.0) sit=EL;
00306         else           sit=EL_EPL;
00307     }
00308     else if (F>=0.0 && F<=_MAX_ERR_F) // current state near/on the surface
00309     {
00310         if (F_tr<=F) sit=EL;
00311         else
00312         {
00313             // Gradients
00314             Tensor2 df_dsig;
00315             REAL    df_dsuc;
00316             _calc_dF_dSig_dF_dSuc(_sig,_z,_suc, df_dsig,df_dsuc);
00317             REAL scalar = Tensors::Dot(df_dsig,dsig_tr) + df_dsuc*dsuc_tr;
00318             if (scalar>=0.0) sit=EPL;
00319             else             sit=EL_EPL;
00320         }
00321     }
00322     else // current state is outside the surface
00323         throw new Fatal("Barcelona::Actualize: F>0.0");
00324 
00325     // Actualize
00326     if (sit==EL)
00327     {
00328         // Ce
00329         Tensor4 Ce;
00330         Tensor2 ce;
00331         _calc_Ce_ce(_v,_sig,_suc, Ce,ce);
00332 
00333         // DEps
00334         Tensors::Dot(Ce,dsig_tr, DEps);
00335         DEps += ce * dsuc_tr; // dEps = Ce:dSig + ce*dSuc
00336         
00337         // Actualize (Elastic)
00338         _v   += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v
00339         _sig += dsig_tr;
00340         _suc += dsuc_tr;
00341         _eps += DEps;
00342         // _z remains unchanged
00343     }
00344     else // EL_EPL or EPL
00345     {
00346         if (sit==EL_EPL)
00347         {
00348             // Find intersection
00349             Tensor2 sig_int;
00350             REAL    suc_int;
00351 
00352             //std::cout << "F = " << F << ", F_tr = " << F_tr << std::endl;
00353 
00354             _find_intersection(dsig_tr,dsuc_tr,F,F_tr,sig_int,suc_int);
00355             //_find_intersection_NR(dsig_tr,dsuc_tr,F,F_tr,sig_int,suc_int);
00356 
00357             if (suc_int>_z(1)) throw new Fatal("Barcelona::Actualize: For DSig not null, SI (suction) surface must not be intersected");
00358 
00359             // Elastic and Elastoplastic parts of stress
00360             Tensor2  dsig_el;  dsig_el = sig_int - _sig;
00361             Tensor2 dsig_epl; dsig_epl = dsig_tr - dsig_el;
00362 
00363             // Elastic and Elastoplastic parts of suction
00364             REAL  dsuc_el = suc_int - _suc;
00365             REAL dsuc_epl = dsuc_tr - dsuc_el;
00366 
00367             // Ce
00368             Tensor4 Ce;
00369             Tensor2 ce;
00370             _calc_Ce_ce(_v,_sig,_suc, Ce,ce);
00371 
00372             // DEps
00373             Tensors::Dot(Ce,dsig_el, DEps);
00374             DEps += ce * dsuc_el; // dEps = Ce:dSig + ce*dSuc
00375             
00376             // Actualize (Elastic)
00377             _v   += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v
00378             _sig += dsig_el;
00379             _suc += dsuc_el;
00380             _eps += DEps;
00381             // _z remains unchanged
00382 
00383             // Continue with elastoplastic part
00384             dsig_tr = dsig_epl;
00385             dsuc_tr = dsuc_epl;
00386         }
00387 
00388 #ifdef USE_FORWARDEULER
00389         // Forward-Euler - Actualize (Elastoplastic)
00390         int n_div=10;
00391         Tensor2 dsig;     dsig    = dsig_tr/n_div;
00392         REAL    dsuc;     dsuc    = dsuc_tr/n_div;
00393         AugVec  dsigsuc;  dsigsuc = dsig(0),dsig(1),dsig(2),dsig(3),dsig(4),dsig(5),dsuc;
00394         for (int i=0; i<n_div; ++i)
00395         {
00396             IntVars dz;
00397             AugVec sigsuc;  sigsuc = _sig(0),_sig(1),_sig(2),_sig(3),_sig(4),_sig(5),_suc;
00398             _tg_deps_dz_given_dsig_and_dsuc(sigsuc, _eps, _z, dsigsuc, DEps, dz);
00399             _v   += -(DEps(0)+DEps(1)+DEps(2))*_v;
00400             _sig += dsig;
00401             _suc += dsuc;
00402             _eps += DEps;
00403             _z   += dz;
00404         }
00405 #else
00406         // Augmented vector
00407         AugVec  sigsuc;   sigsuc = _sig   (0),_sig   (1),_sig   (2),_sig   (3),_sig   (4),_sig   (5),_suc;
00408         AugVec dsigsuc;  dsigsuc = dsig_tr(0),dsig_tr(1),dsig_tr(2),dsig_tr(3),dsig_tr(4),dsig_tr(5),dsuc_tr;
00409 
00410         // TimeInteg
00411         AutoStepME<AugVec,Tensor2,IntVars,Barcelona> TI(CoupledModel::_isc);
00412 
00413         // Initial Eps
00414         Tensor2 eps_ini;  eps_ini = _eps;
00415 
00416         // Update stress/suction (_sig,_suc), internal variables (_z) and strain (_eps) state
00417         int nss = TI.Solve(this,                                        // Address of this instance of GenericEP class
00418                            &Barcelona::_tg_deps_dz_given_dsig_and_dsuc, // Pointer to a function that calculates {dEps,dz}=[[Cep cep];[bk Pk]]:{dSig,dSuc}
00419                            &Barcelona::_local_error,                    // Pointer to a function that calculates local error on dEps
00420                            sigsuc, _eps, _z, dsigsuc);
00421         // Check num_sub_steps
00422         if (nss==-1)
00423             throw new Fatal(_("Barcelona::Actualize: (Loading) Number of max sub-steps (%d) was not sufficient for AutoStepME."), CoupledModel::_isc.ME_maxSS());
00424 
00425         // Calculate the total increment of strain
00426         DEps = _eps - eps_ini;
00427 
00428         // Update internal state
00429         _v += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v
00430 
00431         // Set stress/suction state
00432         _sig = sigsuc(0),sigsuc(1),sigsuc(2),sigsuc(3),sigsuc(4),sigsuc(5);
00433         _suc = sigsuc(6);
00434 #endif
00435     }
00436 
00437 } // }}}
00438 
00439 inline REAL Barcelona::_lam(REAL const & Suc) const // {{{
00440 {
00441     return _lam0*((1.0-_r)*exp(-_bet*Suc) + _r);
00442 } // }}}
00443 
00444 inline REAL Barcelona::_psi(REAL const & Suc) const // {{{
00445 {
00446     return (_lam0-_kap)/(_lam(Suc)-_kap);
00447 } // }}}
00448 
00449 inline REAL Barcelona::_p0(REAL const & Suc, REAL const & z0) const // {{{
00450 {
00451     return _pref*pow(z0/_pref, _psi(Suc));
00452 } // }}}
00453 
00454 inline REAL Barcelona::_yield_func(Tensor2 const & Sig, IntVars const & z, REAL const & Suc) const // {{{
00455 {
00456     REAL p,q;
00457     Stress_p_q(Sig,p,q);
00458     REAL g = - _Mcs*_Mcs * (p+_k*Suc) * (_p0(Suc,z(0))-p);
00459     REAL F = q*q + g;
00460     if (fabs(F)<=_MAX_ERR_F) F=0.0;
00461     return F;
00462 } // }}}
00463 
00464 inline void Barcelona::_calc_Ce_ce(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & Ce, Tensor2 & ce) const // {{{
00465 {
00466     // Invariants
00467     REAL p = (Sig(0)+Sig(1)+Sig(2))/3.0;
00468 
00469     // Calculate Bulk modulae
00470     REAL K  =  p         *v/_kap;
00471     REAL Ks = (Suc+_patm)*v/_kaps;
00472 
00473     // Calc Ce
00474     Tensors::AddScaled(1.0/(2.0*_G), Tensors::Psd, 1.0/(9.0*K), Tensors::IdyI,  Ce);
00475 
00476     // Calc ce
00477     ce = (1.0/(3.0*Ks)) * I;
00478 
00479 } // }}}
00480 
00481 inline void Barcelona::_calc_dF_dSig_dF_dSuc(Tensor2 const & Sig, IntVars const & z, REAL const & Suc, Tensor2 & V, REAL & S) const // {{{
00482 {
00483     // Invariants
00484     Tensor2 devS;
00485     REAL p = (Sig(0)+Sig(1)+Sig(2))/3.0;
00486       devS = Sig - p*I;
00487 
00488     // Suction-related variables
00489     REAL lam = _lam(Suc);
00490     REAL psi = _psi(Suc);
00491     REAL  p0 = _p0(Suc,z(0));
00492     REAL  ps = _k*Suc;
00493 
00494     // Gradients
00495     REAL       MM = _Mcs*_Mcs;
00496     REAL dp0_dSuc = p0*log(z(0)/_pref)*psi*_bet*_lam0*(1.0-_r)*exp(-_bet*Suc)/(lam-_kap);
00497                 V = (MM*(2.0*p - p0 + ps)/3.0)*I + 3.0*devS; // dF_dSig
00498                 S = MM*((p-p0)*_k - (p+ps)*dp0_dSuc);        // dF_dSuc
00499 } // }}}
00500 
00501 inline void Barcelona::_calc_gradients_and_hardening(REAL    const & v, // {{{
00502                                                      Tensor2 const & Sig,
00503                                                      IntVars const & z,
00504                                                      REAL    const & Suc,
00505                                                      Tensor2       & r,
00506                                                      Tensor2       & V,
00507                                                      T_dFdz        & y, 
00508                                                      REAL          & S,
00509                                                      HardMod       & HH,
00510                                                      REAL          & hp) const
00511 {
00512     // Invariants
00513     Tensor2 devS;
00514     REAL p = (Sig(0)+Sig(1)+Sig(2))/3.0;
00515       devS = Sig - p*I;
00516 
00517     // Suction-related variables
00518     REAL lam = _lam(Suc);
00519     REAL psi = _psi(Suc);
00520     REAL  p0 = _p0(Suc,z(0));
00521     REAL  ps = _k*Suc;
00522 
00523     // Gradients
00524     REAL       MM = _Mcs*_Mcs;
00525     REAL      tmp = MM*(2.0*p - p0 + ps)/3.0;
00526     REAL dp0_dSuc = p0*log(z(0)/_pref)*psi*_bet*_lam0*(1.0-_r)*exp(-_bet*Suc)/(lam-_kap);
00527                 r = tmp*I + (3.0*_alp)*devS;          // flow rule
00528                 V = tmp*I +  3.0      *devS;          // dF_dSig
00529                 S = MM*((p-p0)*_k - (p+ps)*dp0_dSuc); // dF_dSuc
00530              y(0) = -MM*psi*p0*(p+ps)/z(0);
00531              y(1) = 0.0;
00532 
00533     // Hardening
00534     REAL  tr_r = r(0)+r(1)+r(2);
00535     REAL  chi0 = (_lam0 - _kap )/v;
00536     REAL  chis = (_lams - _kaps)/v;
00537          HH(0) =  z(0)       *tr_r/chi0;
00538          HH(1) = (z(1)+_patm)*tr_r/chis;
00539 
00540     // Plastic coeficient
00541     hp = -y(0)*HH(0);
00542     
00543 } // }}}
00544 
00545 inline int Barcelona::_tg_deps_dz_given_dsig_and_dsuc(AugVec  const & SigSuc, // {{{
00546                                                       Tensor2 const & Eps,
00547                                                       IntVars const & z,
00548                                                       AugVec  const & dSigSuc,
00549                                                       Tensor2       & dEps,
00550                                                       IntVars       & dz) const
00551 {
00552     // Stop if excessive deviatoric strains
00553     REAL ev,ed;
00554     Tensors::Strain_Ev_Ed(Eps, ev,ed);
00555     if (ed>_EDMAX) return 1; // -1 => failed;  0 => continue;  1 => stop
00556 
00557     // Input values
00558     Tensor2  sig;   sig =  SigSuc(0), SigSuc(1), SigSuc(2), SigSuc(3), SigSuc(4), SigSuc(5);
00559     Tensor2 dsig;  dsig = dSigSuc(0),dSigSuc(1),dSigSuc(2),dSigSuc(3),dSigSuc(4),dSigSuc(5);
00560     REAL     suc;   suc =  SigSuc(6);
00561     REAL    dsuc;  dsuc = dSigSuc(6);
00562     
00563     // Gradients and Hardening
00564     REAL    v=_v; // (constant) current specific volume
00565     Tensor2 r;    // plastic flow direction
00566     Tensor2 V;    // dF_dSig
00567     T_dFdz  y;    // dF_dz
00568     REAL    S;    // dF_dSuc
00569     HardMod HH;   // Hardening modulae
00570     REAL    hp;   // Plastic coeficient
00571     _calc_gradients_and_hardening(v, sig, z, suc, r, V, y, S, HH, hp);
00572 
00573     // Ce and ce
00574     Tensor4 Cep; // Cep <- Ce
00575     Tensor2 ce;
00576     _calc_Ce_ce(v,sig,suc, Cep,ce); // Cep <- Ce
00577 
00578     // Cep
00579     Tensors::Ger(1.0/hp,r,V, Cep); // Cep <- (1/hp)*(r dyadic V) + Ce
00580 
00581     // cep
00582     Tensor2 cep;
00583     cep = (S/hp)*r + ce;
00584 
00585     // dEps
00586     Tensors::Dot(Cep,dsig, dEps);
00587     dEps += cep*dsuc;            // dEps = Cep:dsig + cep*dsuc
00588 
00589     // dLam and dz
00590     REAL dLam = (Tensors::Dot(V,dsig) + S*dsuc)/hp; // dLam = (V:dsig + S*dsuc)/hp
00591            dz = dLam*HH;
00592 
00593     return 0;
00594 } // }}}
00595 
00596 inline REAL Barcelona::_local_error(Tensor2 const & Ey, Tensor2 const & y_high, // {{{
00597                             IntVars const & Ez, IntVars const & z_high) const
00598 {
00599     REAL Err_y = Tensors::Norm(Ey)/Tensors::Norm(y_high);
00600     REAL Err_z = sqrt(blitz::dot(Ez,Ez))/sqrt(blitz::dot(z_high,z_high));
00601     return (Err_y>Err_z ? Err_y : Err_z);
00602 } // }}}
00603 
00604 inline REAL Barcelona::_find_intersection_NR(Tensor2 const & dSig_tr, // {{{
00605                                              REAL    const & dSuc_tr,
00606                                              REAL    const & F      ,
00607                                              REAL    const & F_tr   ,
00608                                              Tensor2       & Sig_int,
00609                                              REAL          & Suc_int) const
00610 {
00611     // Constants
00612     const REAL  ITOL = 1.0e-12;
00613     const int  MAXIT = 10;
00614 
00615     // Find intersection
00616     REAL alpha = -F/(F_tr-F);
00617     REAL F_int;
00618     REAL F_err;
00619     for (int it=1; it<=MAXIT; ++it)
00620     {
00621         // Initial intersection
00622         Sig_int = _sig + alpha*dSig_tr;
00623         Suc_int = _suc + alpha*dSuc_tr;
00624 
00625         // Yield function
00626         F_int = _yield_func(Sig_int,_z,Suc_int);
00627 
00628         // Error
00629         F_err = fabs(F_int)/fabs(_z(0));
00630 
00631         // Check convergence
00632         if (F_err<ITOL) return alpha;
00633 
00634         // Gradients and Hardening
00635         Tensor2 V;
00636         REAL    S;
00637         _calc_dF_dSig_dF_dSuc(Sig_int,_z,Suc_int, V,S);
00638 
00639         // Next coeficient
00640         alpha += -F_int/(Tensors::Dot(V,dSig_tr) + S*dSuc_tr);
00641 
00642         // Check
00643         if (alpha<0.0 || alpha>1.0) throw new Fatal("Barcelona::_find_intersection_NR: Alpha must be inside [0,1].");
00644     }
00645 
00646     // Did not converge
00647     throw new Fatal("Barcelona::_find_intersection_NR: Intersection not found.");
00648 
00649 } // }}}
00650 
00651 inline REAL Barcelona::_find_intersection(Tensor2 const & dSig_tr, // {{{
00652                                           REAL    const & dSuc_tr,
00653                                           REAL    const & F      ,
00654                                           REAL    const & F_tr   ,
00655                                           Tensor2       & Sig_int,
00656                                           REAL          & Suc_int) const
00657 {
00658     // Vars
00659     REAL alpha;
00660     REAL F_int;
00661     
00662     // Find alpha for intersection and the intersection (SIG_int) itself
00663 
00664     //Using Brent’s method, find the root of a function func known to lie between x1 and x2. The
00665     //root, returned as zbrent, will be refined until its accuracy is tol.
00666 
00667     const int  ITMAX=20;   // Maximum allowed number of iterations.
00668     const REAL EPS=3.0e-8; // Machine REALing-point precision.
00669     const REAL tol=1.0e-5;
00670 
00671     REAL a=0.0; // x1
00672     REAL b=1.0; // x2
00673     REAL c=1.0; // x2
00674     REAL d,e,min1,min2;
00675     REAL fa=F;    //(*func)(a);
00676     REAL fb=F_tr; //(*func)(b);
00677     REAL fc,p,q,r,s,tol1,xm;
00678     if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
00679         throw new Fatal("Barcelona::_find_intersection: Root must be bracketed in zbrent");
00680     fc=fb;
00681     int iter;
00682     for (iter=1;iter<=ITMAX;iter++)
00683     {
00684         if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0))
00685         {
00686             c=a;   //Rename a, b, c and adjust bounding interval
00687             fc=fa; // d.
00688             e=d=b-a;
00689         }
00690         if (fabs(fc) < fabs(fb))
00691         {
00692             a=b;
00693             b=c;
00694             c=a;
00695             fa=fb;
00696             fb=fc;
00697             fc=fa;
00698         }
00699         tol1=2.0*EPS*fabs(b)+0.5*tol; // Convergence check.
00700         xm=0.5*(c-b);
00701         if (fabs(xm) <= tol1 || fb == 0.0)
00702         {
00703             //std::cout << "iter = " << iter << std::endl;
00704             //return b; 
00705               alpha = b;
00706             Sig_int = _sig + alpha*dSig_tr;
00707             Suc_int = _suc + alpha*dSuc_tr;
00708             break;
00709         }
00710         if (fabs(e) >= tol1 && fabs(fa) > fabs(fb))
00711         {
00712             s=fb/fa; // Attempt inverse quadratic interpolation.
00713             if (a == c)
00714             {
00715                 p=2.0*xm*s;
00716                 q=1.0-s;
00717             }
00718             else
00719             {
00720                 q=fa/fc;
00721                 r=fb/fc;
00722                 p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
00723                 q=(q-1.0)*(r-1.0)*(s-1.0);
00724             }
00725             if (p > 0.0) q = -q; // Check whether in bounds.
00726             p=fabs(p);
00727             min1=3.0*xm*q-fabs(tol1*q);
00728             min2=fabs(e*q);
00729             if (2.0*p < (min1 < min2 ? min1 : min2))
00730             {
00731                 e=d; // Accept interpolation.
00732                 d=p/q;
00733             }
00734             else
00735             {
00736                 d=xm; // Interpolation failed, use bisection.
00737                 e=d;
00738             }
00739         }
00740         else
00741         { // Bounds decreasing too slowly, use bisection.
00742             d=xm;
00743             e=d;
00744         }
00745         a=b; // Move last best guess to a.
00746         fa=fb;
00747         if (fabs(d) > tol1) // Evaluate new trial root.
00748             b += d;
00749         else
00750             b += Util::Sign(tol1,xm);
00751 
00752         // Calculate function value at intersection
00753         //fb=(*func)(b);
00754         // Trial Intersection
00755           alpha = b;
00756         Sig_int = _sig + alpha*dSig_tr;
00757         Suc_int = _suc + alpha*dSuc_tr;
00758           F_int = _yield_func(Sig_int, _z, Suc_int);
00759              fb = F_int;
00760     }
00761     if (iter>ITMAX)
00762         throw new Fatal("Barcelona::_find_intersection: Maximum number of iterations exceeded in zbrent");
00763 
00764     return alpha;
00765 
00766 } // }}}
00767 
00768 inline void Barcelona::InternalStateValues(Array<REAL> & IntStateVals) const // {{{
00769 {
00770     IntStateVals.resize(3); // 3 = 2 nIntVars plus 1 for the specific volume component
00771     IntStateVals[0] = _z(0);
00772     IntStateVals[1] = _z(1);
00773     IntStateVals[2] = _v;
00774 } // }}}
00775 
00776 inline void Barcelona::InternalStateNames(Array<String> & IntStateNames) const // {{{
00777 {
00778     IntStateNames.resize(3); // 3 = 2 nIntVars plus 1 for the specific volume component
00779     IntStateNames[0] = "z0";
00780     IntStateNames[1] = "z1";
00781     IntStateNames[2] = "v";
00782 } // }}}
00783 
00784 #endif // MECHSYS_BARCELONA_H
00785 
00786 // vim:fdm=marker