models

barcelonax.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_BARCELONAX_H
00023 #define MECHSYS_BARCELONAX_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 // STL
00034 #include <iostream>
00035 
00036 // Blitz++
00037 #include <blitz/tinyvec-et.h>
00038 
00039 // MechSys
00040 #include "models/coupled/coupledmodel.h"
00041 #include "tensors/tensors.h"
00042 #include "tensors/operators.h"
00043 #include "tensors/functions.h"
00044 #include "numerical/integschemesctes.h"
00045 #include "numerical/autostepme.h"
00046 #include "util/string.h"
00047 #include "util/array.h"
00048 #include "util/util.h"
00049 
00050 using Tensors::I;
00051 using Tensors::Tensor2;
00052 using Tensors::Tensor4;
00053 using Tensors::Stress_p_q_S_sin3th; // Mean and deviatoric Cambridge stress invariants
00054 using Tensors::Strain_Ev_Ed;
00055 using Tensors::Mult;
00056 
00058 class BarcelonaX : public CoupledModel
00059 {
00060     static REAL Q_ZERO; // min value of q (used to calc derivatives)
00061 public:
00062 
00064     typedef blitz::TinyVector<REAL,2> IntVars; // {z0,z1} == {p0star,s0}
00065     typedef blitz::TinyVector<REAL,2> HardMod; // {HH0,HH1}
00066     typedef blitz::TinyVector<REAL,2> T_dFdz;  // {y0,y1} derivative of f w.r.t. intvars
00067     typedef blitz::TinyVector<REAL,7> AugVec;  // {Sig/Eps,Suc} Augmented tensor with suction
00068 
00069     // Constructor
00070     BarcelonaX(Array<REAL> const & Prms, Array<REAL> const & IniData);
00071 
00072     // Methods
00073     String Name            () const { return "BarcelonaX"; }
00074     void   TgStiffness     (Tensor4 & D, Tensor2 & d) const;
00075     REAL   Sr              () const; // Degree of saturation
00076     REAL   phi             () const; // Sig_net = Sig_tot - phi*Suc*I
00077     REAL   n_times_dSr_dPp () const; // n:porosity times dSr_dPp (derivative of degree of sat. w.r.t. porepressure)
00078     void   TgPermeability  (Tensor2 & k) const;
00079     void   FlowVelocity    (Tensor1 const & Grad, Tensor1 & Vel) const;
00080     void   Actualize       (Tensor2 const & DSig, REAL const & DPp, Tensor2 & DEps, REAL & DnSr);
00081     void   StressUpdate    (Tensor2 const & DEps, REAL const & DPp, Tensor2 & DSig, REAL & DnSr);
00082     void   FlowUpdate      (REAL const & DPp) {}
00083     void   BackupState     ();
00084     void   RestoreState    ();
00085 
00086     // Methods
00087     REAL YieldFunc          (Tensor2 const & Sig, IntVars const & z, REAL const & Suc)                        const;
00088     void Calc_dFdSig_dFdSuc (Tensor2 const & Sig, IntVars const & z, REAL const & Suc, Tensor2 & V, REAL & S) const;
00089     REAL Calc_p0            (REAL const & Suc, REAL const & z0) const;
00090 
00091     // Derived Access methods
00092     Tensor2 const & Sig() const { return       _sig; };
00093     Tensor2 const & Eps() const { return       _eps; };
00094     REAL            Pp () const { return       _ppr; };
00095     REAL        GammaW () const { return       _gamaW; } 
00096     int  nInternalStateValues() const { return 3; }
00097     void InternalStateValues (Array<REAL>   & IntStateVals ) const;
00098     void InternalStateNames  (Array<String> & IntStateNames) const;
00099 
00100     // Methods to help plot
00101     REAL GetM(REAL Theta) const { REAL sin3th=sin(3.0*Theta); return _calc_M(sin3th); }
00102 
00103 private:
00104     // State vars
00105     REAL    _v;   
00106     Tensor2 _sig; 
00107     REAL    _suc; 
00108     Tensor2 _eps; 
00109     IntVars _z;   
00110     REAL    _ppr; 
00111     
00112     // Backup state vars
00113     REAL    _bkp_v;   
00114     Tensor2 _bkp_sig; 
00115     REAL    _bkp_suc; 
00116     Tensor2 _bkp_eps; 
00117     IntVars _bkp_z;   
00118     REAL    _bkp_ppr; 
00119 
00120     // Parameters (Mechanic)
00121     REAL _lam0;
00122     REAL _kap;
00123     REAL _phics;
00124     REAL _G;
00125     REAL _r;
00126     REAL _bet;
00127     REAL _k;
00128     REAL _lams;
00129     REAL _kaps;
00130     REAL _patm;
00131     REAL _pref;
00132     REAL _BB;
00133 
00134     // Parameters (Hydraulic)
00135     REAL _a;
00136     REAL _b;
00137     REAL _c;
00138     REAL _m;
00139     REAL _ks;
00140     REAL _gamaW;
00141 
00142     // Derived constants
00143     REAL _Mcs;
00144     REAL _wcs;
00145     REAL _alp;
00146 
00147     // Auxiliar variables
00148     REAL _num_dLam; // numerator of dLam
00149     REAL _bkp_num_dLam;
00150 
00151     // Constants
00152     REAL _MAX_ERR_F; // Max error on F (Yield Function)
00153     REAL _MIN_HP;    // Min value of hp (plastic coeficient)
00154 
00155     // Private methods
00156     REAL _lam        (REAL const & Suc) const;
00157     REAL _psi        (REAL const & Suc) const;
00158 
00159     void _calc_Ce_ce(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & Ce, Tensor2 & ce) const;
00160     void _calc_De_de(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & De, Tensor2 & de) const;
00161 
00162     void _calc_gradients_and_hardening(REAL    const & v  ,        // Specific volume
00163                                        Tensor2 const & Sig,        // Stress
00164                                        IntVars const & z  ,        // Internal variables
00165                                        REAL    const & Suc,        // Suction
00166                                        Tensor2       & r  ,        // Plastic flow direction
00167                                        Tensor2       & V  ,        // dF_dSig
00168                                        T_dFdz        & y  ,        // dF_dz
00169                                        REAL          & S  ,        // dF_dSuc
00170                                        HardMod       & HH ,        // Hardening
00171                                        REAL          & hp ) const; // Plastic coeficient
00172 
00173     int _tg_deps_dz_given_dsig_and_dsuc(AugVec  const & SigSuc ,
00174                                         Tensor2 const & Eps    ,
00175                                         IntVars const & z      ,
00176                                         AugVec  const & dSigSuc,
00177                                         Tensor2       & dEps   ,
00178                                         IntVars       & dz     ) const;
00179 
00180     int _tg_dsig_dz_given_deps_and_dsuc(AugVec  const & EpsSuc ,
00181                                         Tensor2 const & Sig    ,
00182                                         IntVars const & z      ,
00183                                         AugVec  const & dEpsSuc,
00184                                         Tensor2       & dSig   ,
00185                                         IntVars       & dz     ) const;
00186 
00187     REAL _local_error(Tensor2 const & Ey, Tensor2 const & y_high,
00188                       IntVars const & Ez, IntVars const & z_high) const;
00189 
00190     // Find intersection Newton-Raphson. Start from the actual state {_sig,_suc}. Returns Alpha.
00191     REAL _find_intersection_NR(Tensor2 const & dSig_tr,
00192                                REAL    const & dSuc_tr,
00193                                REAL    const & F      ,
00194                                REAL    const & F_tr   ,
00195                                Tensor2       & Sig_int,
00196                                REAL          & Suc_int) const;
00197 
00198     // Find intersection Brent's method
00199     REAL _find_intersection(Tensor2 const & dSig_tr,
00200                             REAL    const & dSuc_tr,
00201                             REAL    const & F      ,
00202                             REAL    const & F_tr   ,
00203                             Tensor2       & Sig_int,
00204                             REAL          & Suc_int) const;
00205     // Local
00206     REAL _calc_M(REAL const & sin3th) const;
00207 
00208 }; // class BarcelonaX
00209 
00210 REAL BarcelonaX::Q_ZERO = 1.0e-7;
00211 
00212 
00214 
00215 
00216 inline BarcelonaX::BarcelonaX(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
00217     : _num_dLam(-1.0/*elastic*/) // Deve ser -1 pois a matrix Dep NAO eh definida no inicio
00218 {
00219 
00220     // ##################################################################### Setup Parameters
00221 
00222     /* example.mat
00223        #-----------------------------------------------------------------------------------------------------------------------
00224        #              0     1       2     3     4     5    6     7      8      9     10    11   12   13   14   15     16     17
00225        #           lam0   kap   phics     G     r   bet    k  lams   kaps   patm   pref     B    a    b    c    m     ks  gamaW
00226        BarcelonaX  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
00227     */
00228 
00229     // Check number of parameters
00230     const size_t n_prms = 18;
00231     if (Prms.size()!=n_prms)
00232         throw new Fatal(_("BarcelonaX::BarcelonaX: The number of parameters (%d) is incorrect. It must be equal to %d"), Prms.size(), n_prms);
00233 
00234     // Parameters (mechanical)
00235     _lam0  = Prms[0];
00236     _kap   = Prms[1];
00237     _phics = Prms[2];
00238     _G     = Prms[3];
00239     _r     = Prms[4];
00240     _bet   = Prms[5];
00241     _k     = Prms[6];
00242     _lams  = Prms[7];
00243     _kaps  = Prms[8];
00244     _patm  = Prms[9];
00245     _pref  = Prms[10];
00246     _BB    = Prms[11];
00247 
00248     // Parameters (hydraulic)
00249     _a     = Prms[12];
00250     _b     = Prms[13];
00251     _c     = Prms[14];
00252     _m     = Prms[15];
00253     _ks    = Prms[16];
00254     _gamaW = Prms[17];
00255 
00256     // Derived constants
00257     REAL sin_phi = sin(Util::ToRad(_phics));
00258     _Mcs = 6.0*sin_phi/(3.0-sin_phi);
00259     _wcs = pow((3.0-sin_phi)/(3.0+sin_phi),4.0);
00260     _alp = (_Mcs*(_Mcs-9.0)*(_Mcs-3.0)) / (9.0*(6.0-_Mcs)*(1.0-_kap/_lam0));
00261     
00262     // ##################################################################### Setup Initial State
00263 
00264     // Check number of initial data
00265     if (IniData.size()!=7) // Layered analysis
00266         throw new Fatal(_("BarcelonaX::BarcelonaX: The number of initial data is not sufficiet (it must be equal to 5). { SigX SigY SigZ vini OCR OSI Pp }"));
00267 
00268     // Read initial data
00269     REAL sig_x = IniData[0]; // X stress component
00270     REAL sig_y = IniData[1]; // Y stress component
00271     REAL sig_z = IniData[2]; // Z stress component
00272             _v = IniData[3]; // Initial specific volume
00273     REAL   OCR = IniData[4]; // Over Consolidation Ratio
00274     REAL   OSI = IniData[5]; // Over Suction Increment = Suc_{max} - Suc_{actual}
00275     REAL    Pp = IniData[6]; // Initial pore-pressure
00276 
00277     // Fill stress and strain tensors
00278     REAL sq2=Util::Sq2();
00279     _sig = sig_x, sig_y, sig_z, 0.0*sq2, 0.0*sq2, 0.0*sq2; // Mandel representation
00280     _eps = 0.0,0.0,0.0,0.0*sq2,0.0*sq2,0.0*sq2;
00281 
00282     // Calculate initial state
00283     _ppr = Pp;       // Initial pore-pressure
00284     _suc = 0.0 - Pp; // Initial suction; suc = ua - uw
00285     if (_suc<0.0) _suc=0.0;
00286 
00287     // Calculate internal variables
00288     REAL p,q;
00289     Tensors::Stress_p_q(_sig,p,q);
00290     _z(0) = OCR*p;    // Initial "size" of Cam-clay elipse
00291     _z(1) = OSI+_suc; // Initial "position" of surface for suction
00292 
00293     // Check consistency
00294     //const REAL MAX_F=1.0e-5;
00295     //REAL F = YieldFunc(_sig, _z, _suc);
00296     //if (F>MAX_F) throw new Fatal(_("BarcelonaX::BarcelonaX: Stress point is outside of yield surface (F_ini = %f)"),F);
00297 
00298     // Constants
00299     _MAX_ERR_F = (0.5/100.0)*_pref*_pref; // Max error on F (Yield Function)
00300     _MIN_HP    = (5.0/100.0)*_pref*_pref; // Min value of hp
00301 
00302 } // }}}
00303 
00304 inline void BarcelonaX::TgStiffness(Tensor4 & D, Tensor2 & d) const // {{{
00305 {
00306     if (_num_dLam<0) // Elastic tangent stiffness
00307     {
00308         _calc_De_de(_v,_sig,_suc, D,d); // D<-De, d<-de
00309     }
00310     else // Elastoplastic tangent stiffness
00311     {
00312         // Gradients and Hardening
00313         REAL    v=_v; // (constant) current specific volume
00314         Tensor2 r;    // plastic flow direction
00315         Tensor2 V;    // dF_dSig
00316         T_dFdz  y;    // dF_dz
00317         REAL    S;    // dF_dSuc
00318         HardMod HH;   // Hardening modulae
00319         REAL    hp;   // Plastic coeficient
00320         _calc_gradients_and_hardening(v, _sig, _z, _suc, r, V, y, S, HH, hp);
00321 
00322         // De and de
00323         Tensor4 De;
00324         Tensor2 de;
00325         _calc_De_de(v,_sig,_suc, De,de);
00326 
00327         // Phi and Dep
00328         REAL Phi = Tensors::Reduce(V,De,r)+hp;  // Phi = V:De:r - yk*HHk
00329         Tensors::GerX(-1.0/Phi,De,r,V,De,De,D); // D <- Dep = (-1/Phi) * (De:r)dy(V:De) + De
00330 
00331         // dep
00332         Tensors::Dot(De,r, d); // d <- De:r
00333         if (_suc > 0.0) // if unsaturated 
00334             d = de - ((Tensors::Dot(V,de)+S)/Phi)*d; // d <- dep = de - ((V:de)/Phi)*De:r
00335     }
00336 } // }}}
00337 
00338 inline REAL BarcelonaX::Sr() const // {{{
00339 {
00340     if (_suc <= 0.0) // if saturated
00341         return 1.0;
00342     else
00343         return 1.0/pow(1.0+pow(_suc/_a,_b),_c); // Degree of saturation
00344 } // }}}
00345 
00346 inline REAL BarcelonaX::phi() const // {{{
00347 {
00348     return Sr();
00349 } // }}}
00350 
00351 inline REAL BarcelonaX::n_times_dSr_dPp() const // {{{
00352 {
00353     if (_suc <= 0.0) // if saturated
00354         return 0.0;
00355     else
00356     {
00357         REAL dSr_dsuc = -_b*_c*pow(_suc/_a,_b-1.0)/(_a*pow(1.0+pow(_suc/_a,_b),_c+1.0));
00358         REAL dsuc_dPp = -1.0;
00359         REAL  dSr_dPp = dSr_dsuc*dsuc_dPp;
00360         REAL        n = (_v-1.0)/_v;
00361         return n*dSr_dPp;
00362     }
00363 } // }}}
00364 
00365 inline void BarcelonaX::TgPermeability(Tensor2 & k) const // {{{
00366 {
00367     k = (_ks*pow(Sr(),_m))*I;
00368 } // }}}
00369 
00370 inline void BarcelonaX::FlowVelocity(Tensor1 const & Grad, Tensor1 & Vel) const // {{{
00371 {
00372     Tensor2 K;
00373     TgPermeability(K); 
00374     Tensors::Dot(-K,Grad, Vel);
00375 } // }}}
00376 
00377 inline void BarcelonaX::Actualize(Tensor2 const & DSig, REAL const & DPp, Tensor2 & DEps, REAL & DnSr) // {{{
00378 {
00379     // Initial _eps and nSr
00380     Tensor2 eps_ini = _eps;
00381     REAL    n_ini   = (_v-1.0)/_v;
00382     REAL    Sr_ini  = Sr();
00383     REAL    nSr_ini = n_ini * Sr_ini;
00384 
00385     // Trial state
00386     Tensor2 dsig_tr = DSig;
00387     Tensor2  sig_tr = _sig + dsig_tr;
00388     REAL    dsuc_tr = 0.0 - DPp;
00389     REAL     suc_tr = _suc + dsuc_tr;
00390     if (suc_tr<0.0) { dsuc_tr=0.0-_suc; suc_tr=0.0; }
00391 
00392     // Check loading condition for suction
00393     enum Situation {EL, EL_EPL, EPL}; // EL:     Elastic-Loading
00394                                       // EL_EPL: Elastic-Loading + ElastoPlastic-Loading
00395                                       // EPL:    ElastoPlastic-Loading
00396     REAL F    = YieldFunc(_sig   ,_z,_suc   );
00397     REAL F_tr = YieldFunc( sig_tr,_z, suc_tr);
00398     Situation sit;
00399     if (F<0.0) // current state is inside the surface
00400     {
00401         if (F_tr<=0.0) sit=EL;
00402         else           sit=EL_EPL;
00403     }
00404     else //if (F>=0.0 && F<=_MAX_ERR_F) // current state near/on the surface
00405     {
00406         if (F_tr<=F) sit=EL;
00407         else
00408         {
00409             // Gradients
00410             Tensor2 df_dsig;
00411             REAL    df_dsuc;
00412             Calc_dFdSig_dFdSuc(_sig,_z,_suc, df_dsig,df_dsuc);
00413             REAL scalar = Tensors::Dot(df_dsig,dsig_tr) + df_dsuc*dsuc_tr;
00414             if (scalar>=0.0) sit=EPL;
00415             else             sit=EL_EPL;
00416         }
00417     }
00418     //else // current state is outside the surface
00419         //throw new Fatal("BarcelonaX::Actualize: F>0.0");
00420 
00421     // Actualize
00422     if (sit==EL) // Elastic
00423     {
00424         // Forward-Euler - Actualize (Elastoplastic)
00425         int n_div=_isc.FE_ndiv();
00426         Tensor2 dsig;  dsig = dsig_tr/n_div;
00427         REAL    dsuc;  dsuc = dsuc_tr/n_div;
00428         for (int i=0; i<n_div; ++i)
00429         {
00430             // Ce
00431             Tensor4 Ce;
00432             Tensor2 ce;
00433             _calc_Ce_ce(_v,_sig,_suc, Ce,ce);
00434 
00435             // DEps
00436             Tensors::Dot(Ce,dsig, DEps);
00437             DEps += ce * dsuc; // dEps = Ce:dSig + ce*dSuc
00438             
00439             // Actualize (Elastic)
00440             _sig += dsig;
00441             _suc += dsuc;
00442             _eps += DEps;
00443             // _z remains unchanged
00444         }
00445     }
00446     else // EL_EPL or EPL : Elastoplastic
00447     {
00448         if (sit==EL_EPL)
00449         {
00450             // Find intersection
00451             Tensor2 sig_int;
00452             REAL    suc_int;
00453 
00454             //std::cout << "(Actualize): F = " << F << ", F_tr = " << F_tr << std::endl;
00455 
00456             _find_intersection(dsig_tr,dsuc_tr,F,F_tr,sig_int,suc_int);
00457             //_find_intersection_NR(dsig_tr,dsuc_tr,F,F_tr,sig_int,suc_int);
00458 
00459             // Elastic and Elastoplastic parts of stress
00460             Tensor2  dsig_el;  dsig_el = sig_int - _sig;
00461             Tensor2 dsig_epl; dsig_epl = dsig_tr - dsig_el;
00462 
00463             // Elastic and Elastoplastic parts of suction
00464             REAL  dsuc_el = suc_int - _suc;
00465             REAL dsuc_epl = dsuc_tr - dsuc_el;
00466 
00467             // Ce
00468             Tensor4 Ce;
00469             Tensor2 ce;
00470             _calc_Ce_ce(_v,_sig,_suc, Ce,ce);
00471 
00472             // DEps
00473             Tensors::Dot(Ce,dsig_el, DEps);
00474             DEps += ce * dsuc_el; // dEps = Ce:dSig + ce*dSuc
00475             
00476             // Actualize (Elastic)
00477             _sig += dsig_el;
00478             _suc += dsuc_el;
00479             _eps += DEps;
00480             // _z remains unchanged
00481 
00482             // Continue with elastoplastic part
00483             dsig_tr = dsig_epl;
00484             dsuc_tr = dsuc_epl;
00485         }
00486 
00487         if (CoupledModel::_isc.Type()==IntegSchemesCtes::FE)
00488         {
00489             // Forward-Euler - Actualize (Elastoplastic)
00490             int n_div=_isc.FE_ndiv();
00491             Tensor2 dsig;     dsig    = dsig_tr/n_div;
00492             REAL    dsuc;     dsuc    = dsuc_tr/n_div;
00493             AugVec  dsigsuc;  dsigsuc = dsig(0),dsig(1),dsig(2),dsig(3),dsig(4),dsig(5),dsuc;
00494             for (int i=0; i<n_div; ++i)
00495             {
00496                 IntVars dz;
00497                 AugVec sigsuc;  sigsuc = _sig(0),_sig(1),_sig(2),_sig(3),_sig(4),_sig(5),_suc;
00498                 _tg_deps_dz_given_dsig_and_dsuc(sigsuc, _eps, _z, dsigsuc, DEps, dz);
00499                 _sig += dsig;
00500                 _suc += dsuc;
00501                 _eps += DEps;
00502                 _z   += dz;
00503             }
00504         }
00505         else // suppose it is AutoME
00506         {
00507             // Augmented vector
00508             AugVec  sigsuc;   sigsuc = _sig   (0),_sig   (1),_sig   (2),_sig   (3),_sig   (4),_sig   (5),_suc;
00509             AugVec dsigsuc;  dsigsuc = dsig_tr(0),dsig_tr(1),dsig_tr(2),dsig_tr(3),dsig_tr(4),dsig_tr(5),dsuc_tr;
00510 
00511             // TimeInteg
00512             AutoStepME<AugVec,Tensor2,IntVars,BarcelonaX> TI(CoupledModel::_isc);
00513 
00514             // Update stress/suction (_sig,_suc), internal variables (_z) and strain (_eps) state
00515             int nss = TI.Solve(this,                                        // Address of this instance of GenericEP class
00516                                &BarcelonaX::_tg_deps_dz_given_dsig_and_dsuc, // Pointer to a function that calculates {dEps,dz}=[[Cep cep];[bk Pk]]:{dSig,dSuc}
00517                                &BarcelonaX::_local_error,                    // Pointer to a function that calculates local error on dEps
00518                                sigsuc, _eps, _z, dsigsuc);
00519             // Check num_sub_steps
00520             if (nss==-1)
00521                 throw new Fatal(_("BarcelonaX::Actualize: (Loading) Number of max sub-steps (%d) was not sufficient for AutoStepME."), CoupledModel::_isc.ME_maxSS());
00522 
00523             // Set stress/suction state
00524             _sig = sigsuc(0),sigsuc(1),sigsuc(2),sigsuc(3),sigsuc(4),sigsuc(5);
00525             _suc = sigsuc(6);
00526         }
00527     }
00528 
00529     // Calculate the total (secant) increment of strain
00530     DEps = _eps - eps_ini;
00531 
00532     // Update internal state (secant)
00533     _v += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v
00534 
00535     // Final pore-pressure
00536     _ppr += DPp;
00537 
00538     // Final nSr
00539     REAL n_fin   = (_v-1.0)/_v;
00540     REAL Sr_fin  = Sr();
00541     REAL nSr_fin = n_fin * Sr_fin;
00542     DnSr = nSr_fin - nSr_ini; 
00543 
00544 } // }}}
00545 
00546 inline void BarcelonaX::StressUpdate(Tensor2 const & DEps, REAL const & DPp, Tensor2 & DSig, REAL & DnSr) // {{{
00547 {
00548     // Initial _sig and nSr
00549     Tensor2 sig_ini = _sig;
00550     REAL    Sr_ini  = Sr();
00551 
00552     // Initial n_times_dSr_dPp
00553     REAL n_times_dSr_dPp_ini = n_times_dSr_dPp();
00554 
00555     // De and de
00556     Tensor4 De;
00557     Tensor2 de;
00558     _calc_De_de(_v,_sig,_suc, De,de);
00559 
00560     // Trial state
00561     Tensor2 dsig_tr; Tensors::Dot(De,DEps, dsig_tr); dsig_tr += de*(-DPp);
00562     Tensor2  sig_tr = _sig + dsig_tr;
00563     REAL    dsuc_tr = 0.0 - DPp;
00564     REAL     suc_tr = _suc + dsuc_tr;
00565     if (suc_tr<0.0) { dsuc_tr=0.0-_suc; suc_tr=0.0; }
00566 
00567     // Check loading condition for suction
00568     enum Situation {EL, EL_EPL, EPL}; // EL:     Elastic-Loading
00569                                       // EL_EPL: Elastic-Loading + ElastoPlastic-Loading
00570                                       // EPL:    ElastoPlastic-Loading
00571     REAL F    = YieldFunc(_sig   ,_z,_suc   );
00572     REAL F_tr = YieldFunc( sig_tr,_z, suc_tr);
00573     Situation sit;
00574     if (F<0.0) // current state is inside the surface
00575     {
00576         if (F_tr<=0.0) sit=EL;
00577         else           sit=EL_EPL;
00578     }
00579     else //if (F>=0.0 && F<=_MAX_ERR_F) // current state near/on the surface
00580     {
00581         if (F_tr<=F) sit=EL;
00582         else
00583         {
00584             // Calculate the numerator of dLam
00585             Tensor2 V; // df_dsig
00586             REAL    S; // df_dsuc
00587             Calc_dFdSig_dFdSuc(_sig,_z,_suc, V,S);
00588             _num_dLam = (Tensors::Dot(V,dsig_tr) + S*dsuc_tr); // dLam = (V:De:deps + V:de*dsuc + S*dsuc)/Phi
00589             if (_num_dLam>=0.0) sit=EPL;
00590             else                sit=EL_EPL;
00591         }
00592     }
00593     //else // current state is outside the surface
00594         //throw new Fatal("BarcelonaX::Actualize: F>0.0");
00595 
00596     // Update
00597     if (sit==EL)
00598     {
00599         // Set the numerator of dLam so that it may be used in TgStiffness
00600         _num_dLam = -1.0;
00601 
00602         // DSig
00603         DSig = dsig_tr;
00604         
00605         // Update (Elastic)
00606         _sig += DSig;
00607         _suc += dsuc_tr;
00608         _eps += DEps;
00609         // _z remains unchanged
00610     }
00611     else // EL_EPL or EPL
00612     {
00613         // Set the numerator of dLam so that it may be used in TgStiffness
00614         _num_dLam = +1.0;
00615 
00616         // Auxiliar tensor
00617         Tensor2 deps = DEps;
00618         REAL    dsuc = dsuc_tr;
00619 
00620         // Find intersection
00621         if (sit==EL_EPL)
00622         {
00623             //std::cout << "(StressUpdate): F = " << F << ", F_tr = " << F_tr << std::endl;
00624 
00625             // Find intersection
00626             Tensor2 sig_int;
00627             REAL    suc_int;
00628             _find_intersection(dsig_tr,dsuc_tr,F,F_tr,sig_int,suc_int);
00629 
00630             // Ce
00631             Tensor4 Ce;
00632             Tensor2 ce;
00633             _calc_Ce_ce(_v,_sig,_suc, Ce,ce);
00634 
00635             // Elastic parts of stress, strain and suction
00636             Tensor2 dsig_el;  dsig_el = sig_int - _sig;
00637             REAL              dsuc_el = suc_int - _suc;
00638             Tensor2 deps_el;  Tensors::Dot(Ce,dsig_el, deps_el);
00639             deps_el += ce*dsuc_el; // deps_el = Ce:dsig_el + ce*dsuc_el
00640 
00641             // Actualize (Elastic)
00642             _sig += dsig_el;
00643             _suc += dsuc_el;
00644             _eps += deps_el;
00645             // _z remains unchanged
00646 
00647             // Continue with elastoplastic part
00648             deps = DEps    - deps_el;
00649             dsuc = dsuc_tr - dsuc_el;
00650         }
00651 
00652         if (CoupledModel::_isc.Type()==IntegSchemesCtes::FE)
00653         {
00654             // Forward-Euler - Actualize (Elastoplastic)
00655             int n_div=_isc.FE_ndiv();
00656             deps = deps/n_div;
00657             dsuc = dsuc/n_div;
00658             AugVec depssuc;  depssuc = deps(0),deps(1),deps(2),deps(3),deps(4),deps(5),dsuc;
00659             for (int i=0; i<n_div; ++i)
00660             {
00661                 IntVars dz;
00662                 AugVec epssuc;  epssuc = _eps(0),_eps(1),_eps(2),_eps(3),_eps(4),_eps(5),_suc;
00663                 _tg_dsig_dz_given_deps_and_dsuc(epssuc, _sig, _z, depssuc, DSig, dz);
00664                 _sig += DSig;
00665                 _suc += dsuc;
00666                 _eps += deps;
00667                 _z   += dz;
00668             }
00669         }
00670         else
00671         {
00672             // Augmented vector
00673             AugVec  epssuc;   epssuc = _eps(0),_eps(1),_eps(2),_eps(3),_eps(4),_eps(5),_suc;
00674             AugVec depssuc;  depssuc = deps(0),deps(1),deps(2),deps(3),deps(4),deps(5),dsuc;
00675 
00676             // TimeInteg
00677             AutoStepME<AugVec,Tensor2,IntVars,BarcelonaX> TI(CoupledModel::_isc);
00678 
00679             // Update stress/suction (_sig,_suc), internal variables (_z) and strain (_eps) state
00680             int nss = TI.Solve(this,                                         // Address of this instance of GenericEP class
00681                                &BarcelonaX::_tg_dsig_dz_given_deps_and_dsuc, // Pointer to a function that calculates {dSig,dz}=[[Dep dep];[Bk Qk]]:{dEps,dSuc}
00682                                &BarcelonaX::_local_error,                    // Pointer to a function that calculates local error on dEps
00683                                epssuc, _sig, _z, depssuc);
00684             // Check num_sub_steps
00685             if (nss==-1)
00686                 throw new Fatal(_("BarcelonaX::StressUpdate: (Loading) Number of max sub-steps (%d) was not sufficient for AutoStepME."), CoupledModel::_isc.ME_maxSS());
00687 
00688             // Set strain/suction state
00689             _eps = epssuc(0),epssuc(1),epssuc(2),epssuc(3),epssuc(4),epssuc(5);
00690             _suc = epssuc(6);
00691         }
00692     }
00693 
00694     // Calculate the total (secant) increment of stress
00695     DSig = _sig - sig_ini;
00696 
00697     // Update internal state (secant)
00698     _v += -(DEps(0)+DEps(1)+DEps(2))*_v; // dv=-dEv*v
00699 
00700     // Final pore-pressure
00701     _ppr += DPp;
00702 
00703     // Final nSr
00704     DnSr = -(DEps(0)+DEps(1)+DEps(2))*Sr_ini + n_times_dSr_dPp_ini*DPp; // DnSr = dn*Sr + n*dSr_dPp*DPp ;  dn = -dEv
00705 
00706 } // }}}
00707 
00708 inline void BarcelonaX::BackupState() // {{{
00709 {
00710     _bkp_v   = _v;
00711     _bkp_sig = _sig;
00712     _bkp_suc = _suc;
00713     _bkp_eps = _eps;
00714     _bkp_z   = _z;
00715     _bkp_ppr = _ppr;
00716     _bkp_num_dLam = _num_dLam;
00717 } // }}}
00718 
00719 inline void BarcelonaX::RestoreState() // {{{
00720 {
00721       _v = _bkp_v;
00722     _sig = _bkp_sig;
00723     _suc = _bkp_suc;
00724     _eps = _bkp_eps;
00725       _z = _bkp_z;
00726     _ppr = _bkp_ppr;
00727     _num_dLam = _bkp_num_dLam;
00728 } // }}}
00729 
00730 inline REAL BarcelonaX::_lam(REAL const & Suc) const // {{{
00731 {
00732     return _lam0*((1.0-_r)*exp(-_bet*Suc) + _r);
00733 } // }}}
00734 
00735 inline REAL BarcelonaX::_psi(REAL const & Suc) const // {{{
00736 {
00737     return (_lam0-_kap)/(_lam(Suc)-_kap);
00738 } // }}}
00739 
00740 inline REAL BarcelonaX::Calc_p0(REAL const & Suc, REAL const & z0) const // {{{
00741 {
00742     return _pref*pow(z0/_pref, _psi(Suc));
00743 } // }}}
00744 
00745 inline REAL BarcelonaX::YieldFunc(Tensor2 const & Sig, IntVars const & z, REAL const & Suc) const // {{{
00746 {
00747     REAL     p,q,sin3th;
00748     Tensor2  devS;
00749     Stress_p_q_S_sin3th(Sig, p,q,devS,sin3th);
00750     REAL M = _calc_M(sin3th);
00751     REAL g = - M*M * (p+_k*Suc) * (Calc_p0(Suc,z(0))-p);
00752     REAL C = _pref*_pref*(exp(_BB*(Suc-z(1))/_pref) - exp(-_BB*z(1)/_pref));
00753     REAL F = q*q + g + C;
00754     if (fabs(F)<=_MAX_ERR_F) F=0.0;
00755     return F;
00756 } // }}}
00757 
00758 inline void BarcelonaX::_calc_Ce_ce(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & Ce, Tensor2 & ce) const // {{{
00759 {
00760     // Invariants
00761     REAL p = (Sig(0)+Sig(1)+Sig(2))/3.0;
00762 
00763     // Calculate Bulk modulae
00764     REAL K  =  p         *v/_kap;
00765     REAL Ks = (Suc+_patm)*v/_kaps;
00766 
00767     // Calc Ce
00768     Tensors::AddScaled(1.0/(2.0*_G), Tensors::Psd, 1.0/(9.0*K), Tensors::IdyI,  Ce);
00769 
00770     // Calc ce
00771     ce = (1.0/(3.0*Ks)) * I;
00772 
00773 } // }}}
00774 
00775 inline void BarcelonaX::_calc_De_de(REAL const & v, Tensor2 const & Sig, REAL const & Suc, Tensor4 & De, Tensor2 & de) const // {{{
00776 {
00777     // Invariants
00778     REAL p = (Sig(0)+Sig(1)+Sig(2))/3.0;
00779 
00780     // Calculate Bulk modulae
00781     REAL K  =  p         *v/_kap;
00782     REAL Ks = (Suc+_patm)*v/_kaps;
00783 
00784     // Calc De
00785     Tensors::AddScaled(2.0*_G,Tensors::Psd, K,Tensors::IdyI, De);
00786 
00787     // Calc de
00788     if (_suc <= 0.0) // if saturated
00789         de = 0.0;
00790     else
00791         de = (-K/Ks)*I;
00792 } // }}}
00793 
00794 inline void BarcelonaX::Calc_dFdSig_dFdSuc(Tensor2 const & Sig, IntVars const & z, REAL const & Suc, Tensor2 & V, REAL & S) const // {{{
00795 {
00796     // Invariants
00797     REAL     p,q,sin3th;
00798     Tensor2  devS;
00799     Stress_p_q_S_sin3th(Sig, p,q,devS,sin3th);
00800 
00801     // Suction-related variables
00802     REAL lam = _lam(Suc);
00803     REAL psi = _psi(Suc);
00804     REAL  p0 = Calc_p0(Suc,z(0));
00805     REAL  ps = _k*Suc;
00806 
00807     // Gradients
00808     REAL        M = _calc_M(sin3th);
00809     REAL       MM = M*M;
00810     REAL      tmp = MM*(2.0*p-p0+ps)/3.0;
00811     REAL dp0_dSuc = p0*log(z(0)/_pref)*psi*_bet*_lam0*(1.0-_r)*exp(-_bet*Suc)/(lam-_kap);
00812     REAL  dC_dSuc = _BB*_pref*exp(_BB*(Suc-z(1))/_pref);
00813     REAL  dg_dSuc = MM*((p-p0)*_k - (p+ps)*dp0_dSuc);
00814                 V = tmp*I + 3.0*devS;  // dF_dSig
00815                 S = dg_dSuc + dC_dSuc; // dF_dSuc
00816     if (q>Q_ZERO)
00817     {
00818         REAL  ss3th = pow(sin3th,2.0);
00819         REAL cos3th = (ss3th>1.0 ? 0.0 : sqrt(1.0-ss3th));
00820         if (cos3th>Q_ZERO) // if q is Nan of Inf, prs may become a lot slower
00821         {
00822             REAL                dM_dth = (0.75*M*(1.0-_wcs)*cos3th) / (1.0+_wcs+(_wcs-1.0)*sin3th);
00823             Tensor2 SS;         Mult(devS,devS,SS);
00824             Tensor2 dev_SS;     dev_SS = SS - (I * (SS(0)+SS(1)+SS(2))/3.0);
00825             Tensor2 dth_dsig; dth_dsig = (1.5/(q*q*cos3th)) * (dev_SS*(3.0/q) - devS*sin3th);
00826                                     V += dth_dsig*(2.0*M*(p-p0)*(p+ps)*dM_dth);
00827         }
00828     }
00829 
00830 } // }}}
00831 
00832 inline void BarcelonaX::_calc_gradients_and_hardening(REAL    const & v, // {{{
00833                                                       Tensor2 const & Sig,
00834                                                       IntVars const & z,
00835                                                       REAL    const & Suc,
00836                                                       Tensor2       & r,
00837                                                       Tensor2       & V,
00838                                                       T_dFdz        & y, 
00839                                                       REAL          & S,
00840                                                       HardMod       & HH,
00841                                                       REAL          & hp) const
00842 {
00843     // Invariants
00844     REAL     p,q,sin3th;
00845     Tensor2  devS;
00846     Stress_p_q_S_sin3th(Sig, p,q,devS,sin3th);
00847 
00848     // Suction-related variables
00849     REAL lam = _lam(Suc);
00850     REAL psi = _psi(Suc);
00851     REAL  p0 = Calc_p0(Suc,z(0));
00852     REAL  ps = _k*Suc;
00853 
00854     // Gradients
00855     REAL        C = _pref*_pref*(exp(_BB*(Suc-z(1))/_pref) - exp(-_BB*z(1)/_pref));
00856     REAL        M = _calc_M(sin3th);
00857     REAL       MM = M*M;
00858     REAL      tmp = MM*(2.0*p-p0+ps)/3.0;
00859     REAL dp0_dSuc = p0*log(z(0)/_pref)*psi*_bet*_lam0*(1.0-_r)*exp(-_bet*Suc)/(lam-_kap);
00860     REAL  dC_dSuc = _BB*_pref*exp(_BB*(Suc-z(1))/_pref);
00861     REAL  dg_dSuc = MM*((p-p0)*_k - (p+ps)*dp0_dSuc);
00862                 r = tmp*I + (3.0*_alp)*devS; // flow rule
00863                 V = tmp*I +  3.0      *devS; // dF_dSig
00864                 S = dg_dSuc + dC_dSuc;       // dF_dSuc
00865     if (q>Q_ZERO)
00866     {
00867         REAL  ss3th = pow(sin3th,2.0);
00868         REAL cos3th = (ss3th>1.0 ? 0.0 : sqrt(1.0-ss3th));
00869         if (cos3th>Q_ZERO) // if q is Nan of Inf, prs may become a lot slower
00870         {
00871             REAL                dM_dth = (0.75*M*(1.0-_wcs)*cos3th) / (1.0+_wcs+(_wcs-1.0)*sin3th);
00872             Tensor2 SS;         Mult(devS,devS,SS);
00873             Tensor2 dev_SS;     dev_SS = SS - (I * (SS(0)+SS(1)+SS(2))/3.0);
00874             Tensor2 dth_dsig; dth_dsig = (1.5/(q*q*cos3th)) * (dev_SS*(3.0/q) - devS*sin3th);
00875                                     V += dth_dsig*(2.0*M*(p-p0)*(p+ps)*dM_dth);
00876                                     r += dth_dsig*(2.0*M*(p-p0)*(p+ps)*dM_dth);
00877         }
00878     }
00879     y(0) = -MM*psi*p0*(p+ps)/z(0);
00880     y(1) = -_BB*C/_pref;
00881 
00882     // Hardening
00883     REAL  tr_r = r(0)+r(1)+r(2);
00884     REAL  chi0 = (_lam0 - _kap )/v;
00885     REAL  chis = (_lams - _kaps)/v;
00886          HH(0) =  z(0)       *tr_r/chi0;
00887          HH(1) = (z(1)+_patm)*tr_r/chis;
00888 
00889     // Plastic coeficient
00890     hp = -(y(0)*HH(0)+y(1)*HH(1));
00891     
00892 } // }}}
00893 
00894 inline int BarcelonaX::_tg_deps_dz_given_dsig_and_dsuc(AugVec  const & SigSuc, // {{{
00895                                                        Tensor2 const & Eps,
00896                                                        IntVars const & z,
00897                                                        AugVec  const & dSigSuc,
00898                                                        Tensor2       & dEps,
00899                                                        IntVars       & dz) const
00900 {
00901     // Stop if excessive deviatoric strains
00902     REAL ev,ed;
00903     Tensors::Strain_Ev_Ed(Eps, ev,ed);
00904     if (ed>CoupledModel::_isc.EdMax()) // return 1; // -1 => failed;  0 => continue;  1 => stop
00905         throw new Warning(_("BarcelonaX::_tg_deps_dz_given_dsig_and_dsuc: Max devaitoric strain Ed=%f greater than EdMax=%f"),ed,_isc.EdMax());
00906 
00907     // Input values
00908     Tensor2  sig;   sig =  SigSuc(0), SigSuc(1), SigSuc(2), SigSuc(3), SigSuc(4), SigSuc(5);
00909     Tensor2 dsig;  dsig = dSigSuc(0),dSigSuc(1),dSigSuc(2),dSigSuc(3),dSigSuc(4),dSigSuc(5);
00910     REAL     suc;   suc =  SigSuc(6);
00911     REAL    dsuc;  dsuc = dSigSuc(6);
00912     
00913     // Gradients and Hardening
00914     REAL    v=_v; // (constant) current specific volume
00915     Tensor2 r;    // plastic flow direction
00916     Tensor2 V;    // dF_dSig
00917     T_dFdz  y;    // dF_dz
00918     REAL    S;    // dF_dSuc
00919     HardMod HH;   // Hardening modulae
00920     REAL    hp;   // Plastic coeficient
00921     _calc_gradients_and_hardening(v, sig, z, suc, r, V, y, S, HH, hp);
00922 
00923     // Try to check for softening initial point
00924     if (fabs(hp/(sig(0)+sig(1)+sig(2)))<=_MIN_HP)
00925         throw new Warning(_("BarcelonaX::_tg_deps_dz_given_dsig_and_dsuc: hp=%f (zero). Maybe it is the start of softening behaviour"),hp);
00926 
00927     // Ce and ce
00928     Tensor4 Cep; // Cep <- Ce
00929     Tensor2 ce;
00930     _calc_Ce_ce(v,sig,suc, Cep,ce); // Cep <- Ce
00931 
00932     // Cep
00933     Tensors::Ger(1.0/hp,r,V, Cep); // Cep <- (1/hp)*(r dyadic V) + Ce
00934 
00935     // cep
00936     Tensor2 cep;
00937     cep = (S/hp)*r + ce;
00938 
00939     // dEps
00940     Tensors::Dot(Cep,dsig, dEps);
00941     dEps += cep*dsuc;            // dEps = Cep:dsig + cep*dsuc
00942 
00943     // dLam and dz
00944     REAL dLam = (Tensors::Dot(V,dsig) + S*dsuc)/hp; // dLam = (V:dsig + S*dsuc)/hp
00945            dz = dLam*HH;
00946 
00947     return 0;
00948 } // }}}
00949 
00950 inline int BarcelonaX::_tg_dsig_dz_given_deps_and_dsuc(AugVec  const & EpsSuc, // {{{
00951                                                        Tensor2 const & Sig,
00952                                                        IntVars const & z,
00953                                                        AugVec  const & dEpsSuc,
00954                                                        Tensor2       & dSig,
00955                                                        IntVars       & dz) const
00956 {
00957     // Input values
00958     Tensor2  eps;   eps =  EpsSuc(0), EpsSuc(1), EpsSuc(2), EpsSuc(3), EpsSuc(4), EpsSuc(5);
00959     Tensor2 deps;  deps = dEpsSuc(0),dEpsSuc(1),dEpsSuc(2),dEpsSuc(3),dEpsSuc(4),dEpsSuc(5);
00960     REAL     suc;   suc =  EpsSuc(6);
00961     REAL    dsuc;  dsuc = dEpsSuc(6);
00962     
00963     // Gradients and Hardening
00964     REAL    v=_v; // (constant) current specific volume
00965     Tensor2 r;    // plastic flow direction
00966     Tensor2 V;    // dF_dSig
00967     T_dFdz  y;    // dF_dz
00968     REAL    S;    // dF_dSuc
00969     HardMod HH;   // Hardening modulae
00970     REAL    hp;   // Plastic coeficient
00971     _calc_gradients_and_hardening(v, Sig, z, suc, r, V, y, S, HH, hp);
00972 
00973     // De and de
00974     Tensor4 De;
00975     Tensor2 de;
00976     _calc_De_de(v,Sig,suc, De,de);
00977 
00978     // Phi and Dep
00979     Tensor4 Dep;
00980     REAL Phi = Tensors::Reduce(V,De,r)+hp;    // Phi = V:De:r - yk*HHk
00981     Tensors::GerX(-1.0/Phi,De,r,V,De,De,Dep); // Dep = (-1/Phi) * (De:r)dy(V:De) + De
00982 
00983     // dep
00984     Tensor2 dep;
00985     Tensors::Dot(De,r, dep); // dep <- De:r
00986     dep = de - ((Tensors::Dot(V,de)+S)/Phi)*dep; // dep <- de - ((V:de)/Phi)*De:r
00987 
00988     // dSig
00989     Tensors::Dot(Dep,deps, dSig);
00990     dSig += dep*dsuc; // dSig = Dep:deps + dep*dsuc
00991 
00992     // dLam and dz
00993     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
00994            dz = dLam*HH;
00995 
00996     return 0;
00997 } // }}}
00998 
00999 inline REAL BarcelonaX::_local_error(Tensor2 const & Ey, Tensor2 const & y_high, // {{{
01000                                      IntVars const & Ez, IntVars const & z_high) const
01001 {
01002     REAL Err_y = Tensors::Norm(Ey)/Tensors::Norm(y_high);
01003     REAL Err_z = sqrt(blitz::dot(Ez,Ez))/sqrt(blitz::dot(z_high,z_high));
01004     return (Err_y>Err_z ? Err_y : Err_z);
01005 } // }}}
01006 
01007 inline REAL BarcelonaX::_find_intersection_NR(Tensor2 const & dSig_tr, // {{{
01008                                               REAL    const & dSuc_tr,
01009                                               REAL    const & F      ,
01010                                               REAL    const & F_tr   ,
01011                                               Tensor2       & Sig_int,
01012                                               REAL          & Suc_int) const
01013 {
01014     // Constants
01015     const REAL  ITOL = 1.0e-12;
01016     const int  MAXIT = 10;
01017 
01018     // Find intersection
01019     REAL alpha = -F/(F_tr-F);
01020     REAL F_int;
01021     REAL F_err;
01022     for (int it=1; it<=MAXIT; ++it)
01023     {
01024         // Initial intersection
01025         Sig_int = _sig + alpha*dSig_tr;
01026         Suc_int = _suc + alpha*dSuc_tr;
01027 
01028         // Yield function
01029         F_int = YieldFunc(Sig_int,_z,Suc_int);
01030 
01031         // Error
01032         F_err = fabs(F_int)/fabs(_z(0));
01033 
01034         // Check convergence
01035         if (F_err<ITOL) return alpha;
01036 
01037         // Gradients and Hardening
01038         Tensor2 V;
01039         REAL    S;
01040         Calc_dFdSig_dFdSuc(Sig_int,_z,Suc_int, V,S);
01041 
01042         // Next coeficient
01043         alpha += -F_int/(Tensors::Dot(V,dSig_tr) + S*dSuc_tr);
01044 
01045         // Check
01046         if (alpha<0.0 || alpha>1.0) throw new Fatal("BarcelonaX::_find_intersection_NR: Alpha must be inside [0,1].");
01047     }
01048 
01049     // Did not converge
01050     throw new Fatal("BarcelonaX::_find_intersection_NR: Intersection not found.");
01051 
01052 } // }}}
01053 
01054 inline REAL BarcelonaX::_find_intersection(Tensor2 const & dSig_tr, // {{{
01055                                            REAL    const & dSuc_tr,
01056                                            REAL    const & F      ,
01057                                            REAL    const & F_tr   ,
01058                                            Tensor2       & Sig_int,
01059                                            REAL          & Suc_int) const
01060 {
01061     // Vars
01062     REAL alpha;
01063     REAL F_int;
01064     
01065     // Find alpha for intersection and the intersection (SIG_int) itself
01066 
01067     //Using Brent’s method, find the root of a function func known to lie between x1 and x2. The
01068     //root, returned as zbrent, will be refined until its accuracy is tol.
01069 
01070     const int  ITMAX=40;   // Maximum allowed number of iterations.
01071     const REAL EPS=3.0e-8; // Machine REALing-point precision.
01072     const REAL tol=1.0e-5;
01073 
01074     REAL a=0.0; // x1
01075     REAL b=1.0; // x2
01076     REAL c=1.0; // x2
01077     REAL d=0.0,e=0.0,min1=0.0,min2=0.0;
01078     REAL fa=F;    //(*func)(a);
01079     REAL fb=F_tr; //(*func)(b);
01080     REAL fc,p,q,r,s,tol1,xm;
01081     if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
01082         throw new Fatal("BarcelonaX::_find_intersection: Root must be bracketed in zbrent");
01083     fc=fb;
01084     int iter;
01085     for (iter=1;iter<=ITMAX;iter++)
01086     {
01087         if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0))
01088         {
01089             c=a;   //Rename a, b, c and adjust bounding interval
01090             fc=fa; // d.
01091             e=d=b-a;
01092         }
01093         if (fabs(fc) < fabs(fb))
01094         {
01095             a=b;
01096             b=c;
01097             c=a;
01098             fa=fb;
01099             fb=fc;
01100             fc=fa;
01101         }
01102         tol1=2.0*EPS*fabs(b)+0.5*tol; // Convergence check.
01103         xm=0.5*(c-b);
01104         if (fabs(xm) <= tol1 || fb == 0.0)
01105         {
01106               alpha = b;
01107             Sig_int = _sig + alpha*dSig_tr;
01108             Suc_int = _suc + alpha*dSuc_tr;
01109             break;
01110         }
01111         if (fabs(e) >= tol1 && fabs(fa) > fabs(fb))
01112         {
01113             s=fb/fa; // Attempt inverse quadratic interpolation.
01114             if (a == c)
01115             {
01116                 p=2.0*xm*s;
01117                 q=1.0-s;
01118             }
01119             else
01120             {
01121                 q=fa/fc;
01122                 r=fb/fc;
01123                 p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
01124                 q=(q-1.0)*(r-1.0)*(s-1.0);
01125             }
01126             if (p > 0.0) q = -q; // Check whether in bounds.
01127             p=fabs(p);
01128             min1=3.0*xm*q-fabs(tol1*q);
01129             min2=fabs(e*q);
01130             if (2.0*p < (min1 < min2 ? min1 : min2))
01131             {
01132                 e=d; // Accept interpolation.
01133                 d=p/q;
01134             }
01135             else
01136             {
01137                 d=xm; // Interpolation failed, use bisection.
01138                 e=d;
01139             }
01140         }
01141         else
01142         { // Bounds decreasing too slowly, use bisection.
01143             d=xm;
01144             e=d;
01145         }
01146         a=b; // Move last best guess to a.
01147         fa=fb;
01148         if (fabs(d) > tol1) // Evaluate new trial root.
01149             b += d;
01150         else
01151             b += Util::Sign(tol1,xm);
01152 
01153         // Calculate function value at intersection
01154         //fb=(*func)(b);
01155         // Trial Intersection
01156           alpha = b;
01157         Sig_int = _sig + alpha*dSig_tr;
01158         Suc_int = _suc + alpha*dSuc_tr;
01159           F_int = YieldFunc(Sig_int, _z, Suc_int);
01160              fb = F_int;
01161     }
01162     if (iter>ITMAX)
01163         throw new Fatal("BarcelonaX::_find_intersection: Maximum number of iterations exceeded in zbrent");
01164 
01165     return alpha;
01166 
01167 } // }}}
01168 
01169 inline void BarcelonaX::InternalStateValues(Array<REAL> & IntStateVals) const // {{{
01170 {
01171     IntStateVals.resize(3); // 3 = 2 nIntVars plus 1 for the specific volume component
01172     IntStateVals[0] = _z(0);
01173     IntStateVals[1] = _z(1);
01174     IntStateVals[2] = _v;
01175 } // }}}
01176 
01177 inline void BarcelonaX::InternalStateNames(Array<String> & IntStateNames) const // {{{
01178 {
01179     IntStateNames.resize(3); // 3 = 2 nIntVars plus 1 for the specific volume component
01180     IntStateNames[0] = "z0";
01181     IntStateNames[1] = "z1";
01182     IntStateNames[2] = "v";
01183 } // }}}
01184 
01185 inline REAL BarcelonaX::_calc_M(REAL const & sin3th) const // {{{
01186 {
01187     return _Mcs*pow( 2.0*_wcs/(1.0+_wcs+(_wcs-1.0)*sin3th) ,0.25);
01188 } // }}}
01189 
01190 
01192 
01193 
01194 // Allocate a new BarcelonaX model
01195 CoupledModel * BarcelonaXMaker(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
01196 {
01197     return new BarcelonaX(Prms, IniData);
01198 } // }}}
01199 
01200 // Register an BarcelonaX model into ModelFactory array map
01201 int BarcelonaXRegister() // {{{
01202 {
01203     CoupledModelFactory["BarcelonaX"] = BarcelonaXMaker;
01204     return 0;
01205 } // }}}
01206 
01207 // Execute the autoregistration
01208 int __BarcelonaX_dummy_int = BarcelonaXRegister();
01209 
01210 #endif // MECHSYS_BARCELONAX_H
01211 
01212 // vim:fdm=marker
Generated on Wed Jan 24 15:56:27 2007 for MechSys by  doxygen 1.4.7