constmod

failcrit.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_FAILCRIT_H
00023 #define MECHSYS_FAILCRIT_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 #include <vector>
00034 
00035 // Tensors
00036 #include "tensors/tensor1.h"
00037 #include "tensors/tensor2.h"
00038 #include "tensors/tensoperators.h"
00039 
00040 // VTKWrap
00041 #include "vtkwrap/structgrid.h"
00042 #include "vtkwrap/sgridisosurf.h"
00043 #include "vtkwrap/colors.h"
00044 #include "vtkwrap/vtkwin.h"
00045 
00046 // MechSys
00047 #include "util/string.h"
00048 #include "util/util.h"
00049 #include "tensors/functions.h"
00050 
00051 using Tensors::Sin3Th;
00052 using Tensors::Hid2Sig;
00053 using Util::ToRad;
00054 using Util::Sq2;
00055 using Util::Sq3;
00056 using Util::Sq6;
00057 
00058 class FailCrit
00059 {
00060     static REAL BIGNUM;
00061 public:
00062     // Enum
00063     enum Type {FailCrit_TR,  // Tresca
00064                FailCrit_VM,  // von Mises
00065                FailCrit_MC,  // Mohr-Coulomb
00066                FailCrit_LD,  // Lade-Duncan
00067                FailCrit_DP,  // Drucker-Prager
00068                FailCrit_MN,  // Matsuoka-Nakai
00069                FailCrit_SM,  // SMP (Novo)
00070                FailCrit_AR,  // Argyris-Sheng et al.
00071                FailCrit_AN}; // General Anisotropic
00072 
00073     // Anisotropic criteria n unit normal to the mobilized plane
00074     enum nType {FailCrit_Oct,  // octaedral plane
00075                 FailCrit_SMP}; // SMP (Spatially Mobilized Plane)
00076 
00077     // non-linear envelop
00078     enum EnvType {FailCrit_LIN, FailCrit_NONLIN_1};
00079 
00080     // Struct
00081     struct Prms
00082     {
00083         REAL Phi;  // Shear angle (radians)
00084         REAL c;    // Cohesion
00085         REAL Su;   // Undrained strength
00086         REAL pu;   // p for Su calculation; equal to NONLIN_1's Pref
00087         REAL ax;   // Bedding plane x component
00088         REAL ay;   // Bedding plane y component
00089         REAL az;   // Bedding plane z component
00090         REAL Alp;  // Anisotropic Alpha coeficient
00091         REAL R;    // Anisotropic R (shear angle tangent)
00092         REAL C;    // Anisotropic C (cohesion)
00093         REAL Psi;  // Envelop nonlinearity coeficient
00094     }; // struct Prms
00095 
00096     // Constructor & Destructor
00097      FailCrit(Type const & type, nType const & n_type=FailCrit_SMP, Prms const * pParameters=NULL);
00098     ~FailCrit();
00099 
00100     // Config Methods
00101     FailCrit & SetPhi     (REAL           PhiDeg ) { _prms.Phi = ToRad(PhiDeg); _calc_constants(); return (*this); }
00102     FailCrit & Setc       (REAL           c      ) { _prms.c   = c;             _calc_constants(); return (*this); }
00103     FailCrit & SetSu      (REAL           Su     ) { _prms.Su  = Su;            _calc_constants(); return (*this); }
00104     FailCrit & Setpu      (REAL           pu     ) { _prms.pu  = pu;            _calc_constants(); return (*this); }
00105     FailCrit & Setax      (REAL           ax     ) { _prms.ax  = ax;            _calc_constants(); return (*this); }
00106     FailCrit & Setay      (REAL           ay     ) { _prms.ay  = ay;            _calc_constants(); return (*this); }
00107     FailCrit & Setaz      (REAL           az     ) { _prms.az  = az;            _calc_constants(); return (*this); }
00108     FailCrit & SetAlp     (REAL           Alp    ) { _prms.Alp = Alp;           _calc_constants(); return (*this); }
00109     FailCrit & SetR       (REAL           R      ) { _prms.R   = R;             _calc_constants(); return (*this); }
00110     FailCrit & SetC       (REAL           C      ) { _prms.C   = C;             _calc_constants(); return (*this); }
00111     FailCrit & SetPsi     (REAL           Psi    ) { _prms.Psi = Psi;                              return (*this); }
00112     FailCrit & SetnType   (nType          n_type ) { _n_type   = n_type;                           return (*this); }
00113     FailCrit & SetEnvType (EnvType        envType) { _env_type = envType;                          return (*this); }
00114     FailCrit & SetMinp    (REAL           Minp   ) { _min_p    = Minp;                             return (*this); }
00115     FailCrit & SetMaxp    (REAL           Maxp   ) { _max_p    = Maxp;                             return (*this); }
00116     FailCrit & SetColor   (String const & Color  ) { _iso_clr  = Color;                            return (*this); }
00117     FailCrit & SetOpac    (REAL           Opacity) { _iso_opac = Opacity;                          return (*this); }
00118 
00119     // Methods
00120     REAL         Func       (REAL SI, REAL SII, REAL SIII);
00121     REAL         AR_Calc_M  (REAL sin3th) const;
00122     String       Name       () const;
00123     Type         GetType    () const { return _type;     }
00124     nType        GetnType   () const { return _n_type;   }
00125     EnvType      GetEnvType () const { return _env_type; }
00126     REAL         Eta        () const { return _eta;      }
00127     REAL         Minp       () const { return _min_p;    }
00128     REAL         Maxp       () const { return _max_p;    }
00129     REAL         kDP        () const { return _k_DP;     }
00130     REAL         kSM        () const { return _k_SM;     }
00131     REAL         AR_w       () const { return _AR_w;     }
00132     Prms const & GetPrms    () const { return _prms;     }
00133 
00134     // VTK
00135     vtkActor * GenIsoSurf    (int nPoints);
00136     void       AddActorsTo   (VTKWin & Win);
00137     void       DelActorsFrom (VTKWin & Win);
00138 
00139 private:
00140     // Data
00141     Type           _type;
00142     nType          _n_type;
00143     EnvType        _env_type;
00144     Prms           _prms;
00145     REAL           _eta;
00146     REAL           _sinPhi;
00147     REAL           _k_DP;
00148     REAL           _k_LD;
00149     REAL           _k_MN;
00150     REAL           _k_SM;
00151     REAL           _AR_Mmax;
00152     REAL           _AR_w;
00153     REAL           _min_p;
00154     REAL           _max_p;
00155     // VTK
00156     SGridIsoSurf * _iso_surf;
00157     String         _iso_clr;
00158     REAL           _iso_opac;
00159 
00160     // Methods
00161     void _calc_constants();
00162 }; // class FailCrit
00163 
00164 REAL FailCrit::BIGNUM = 1.0e+10;
00165 
00166 
00168 
00169 
00170 inline FailCrit::FailCrit(Type const & type, nType const & n_type, Prms const * pParameters) // {{{
00171     : _type     (type),
00172       _n_type   (n_type),
00173       _env_type (FailCrit_LIN),
00174       _iso_surf (NULL),
00175       _iso_clr  (String("blue_light")),
00176       _iso_opac (1.0)
00177 {
00178     // Parameters
00179     if (pParameters==NULL)
00180     {
00181         _prms.Phi = ToRad(26.0);
00182         _prms.c   = 0.0;
00183         _prms.Su  = 0.0;
00184         _prms.pu  = 1.0;
00185         _prms.ax  = 0.0;
00186         _prms.ay  = 0.0;
00187         _prms.az  = 1.0;
00188         _prms.Alp = 0.1;
00189         _prms.R   = 0.37;
00190         _prms.C   = 0.0;
00191         _prms.Psi = 0.5;
00192     }
00193     else _prms = (*pParameters);
00194 
00195     // Min and Max p for generation of the isosurface
00196     _min_p = 1.0e-5;
00197     _max_p = _prms.pu;
00198 
00199     // Constants
00200     _calc_constants();
00201 
00202 } // }}}
00203 
00204 inline FailCrit::~FailCrit() // {{{
00205 {
00206     if (_iso_surf!=NULL) delete _iso_surf;
00207 } // }}}
00208 
00209 inline REAL FailCrit::Func(REAL SI, REAL SII, REAL SIII) // {{{
00210 {
00211     REAL ff;
00212     switch (_type)
00213     {
00214         case FailCrit_TR:
00215         {
00216             std::vector<REAL> s;
00217             s.push_back(fabs(SI   - SII ));
00218             s.push_back(fabs(SII  - SIII));
00219             s.push_back(fabs(SIII - SI  ));
00220             REAL s_max = *max_element(s.begin(),s.end());
00221             ff = s_max/2.0 - _prms.Su;
00222             break;
00223         }
00224         case FailCrit_VM:
00225         {
00226             REAL q  = sqrt(pow(SI-SII,2.0)+pow(SII-SIII,2.0)+pow(SIII-SI,2.0))/Sq2();
00227                  ff = q - 2.0*_prms.Su;
00228             break;
00229         }
00230         case FailCrit_MC:
00231         {
00232             REAL s[] = {SI,SII,SIII};
00233             Util::Sort(s,3); // ascending
00234             ff = (s[2]-s[0]) - (s[2]+s[0])*_sinPhi;
00235             break;
00236         }
00237         case FailCrit_LD:
00238         {
00239             if (SI<=0.0 || SII<=0.0 || SIII<=0.0) return BIGNUM;
00240             REAL I1 = SI+SII+SIII;
00241             REAL I3 = SI*SII*SIII;
00242                  ff = pow(I1,3.0)/I3 - _k_LD;
00243             break;
00244         }
00245         case FailCrit_DP:
00246         {
00247             REAL P  = (SI+SII+SIII)/Sq3();
00248             REAL Q  = sqrt(SI*SI + SII*SII + SIII*SIII - P*P);
00249                  ff = Q/P - _k_DP;
00250             break;
00251         }
00252         case FailCrit_MN:
00253         {
00254             if (SI<=0.0 || SII<=0.0 || SIII<=0.0) return BIGNUM;
00255             REAL I1 = SI+SII+SIII;
00256             REAL I2 = SI*SII + SII*SIII + SIII*SI;
00257             REAL I3 = SI*SII*SIII;
00258                  ff = I1*I2/I3 - _k_MN;
00259             break;
00260         }
00261         case FailCrit_SM:
00262         {
00263             if (SI<=0.0 || SII<=0.0 || SIII<=0.0) return BIGNUM;
00264             REAL I2 = SI*SII + SII*SIII + SIII*SI;
00265             REAL I3 = SI*SII*SIII;
00266             REAL P  = sqrt(I3/I2)*(sqrt(SI)+sqrt(SII)+sqrt(SIII));
00267             REAL Q  = sqrt(SI*SI + SII*SII + SIII*SIII - P*P);
00268             if (_env_type==FailCrit_NONLIN_1) ff = Q/(P*exp(-_prms.Psi*(P-_prms.pu))) - _k_SM;
00269             else                              ff = Q/P - _k_SM;
00270             break;
00271         }
00272         case FailCrit_AR:
00273         {
00274             REAL P  = (SI+SII+SIII)/Sq3();
00275             REAL Q  = sqrt(SI*SI + SII*SII + SIII*SIII - P*P);
00276             REAL mu = AR_Calc_M(Sin3Th(SI,SII,SIII));
00277                  ff = Q/P - mu;
00278             break;
00279         }
00280         case FailCrit_AN:
00281         {
00282             if (SI<=0.0 || SII<=0.0 || SIII<=0.0) return BIGNUM;
00283 
00284             REAL  a0 = _prms.az;
00285             REAL  a1 = _prms.ax;
00286             REAL  a2 = _prms.ay;
00287             REAL alp = _prms.Alp;
00288             REAL   R = _prms.R;
00289             REAL   C = _prms.C;
00290 
00291             using std::vector;
00292 
00293             // First and Second order 3D Tensors
00294             typedef TensorsLib::Tensor1<REAL,3> T1; // <<< from Tensors C++ library
00295             typedef TensorsLib::Tensor2<REAL,3> T2; // <<< from Tensors C++ library
00296 
00297             // Declarations
00298             T2           s(0.0); // Stress tensor
00299             T2           ss;     // Stress tensor squared     
00300             vector<REAL> L;      // Eigenvalues               
00301             vector<T2>   P;      // Eigenprojectors           
00302             vector<REAL> Iv;     // Invariants                
00303             T1           n;      // Normal to Mobilized Plane             
00304             T2           Pnn;    // Mobilized Plane projector             
00305             T1           a;      // Normal to bedding planes  
00306             T1           b;      // "Aniso" oblique direction 
00307             T1           bet;    // "Normalized" direction    
00308             T2           Pbn;    // Oblique projector         
00309             REAL         k;      // Temporary = (s % Pnn)
00310             REAL         sig;    // Normal (oblique) stress on n plane
00311             REAL         tau;    // Shear (oblique) stress on n plane
00312 
00313             // Initialization
00314             s[0][0]=SI; s[1][1]=SII; s[2][2]=SIII;
00315                a[0]=a0;    a[1]=a1;     a[2]=a2;
00316             ss = s*s;
00317 
00318             // Get eigen{values,projectors} of stress tensor
00319             s.GetEigen(L,P);
00320 
00321             // Characteristic invariants of stress tensor
00322             Iv.resize(3);
00323             Iv[0] = s.Trace();
00324             Iv[1] = (Iv[0]*Iv[0] - ss.Trace())/2.0;
00325             Iv[2] = s.Det();
00326 
00327             // Normal to Mobilized Plane
00328             switch (_n_type)
00329             {
00330             case FailCrit_Oct:
00331                 for (int i=0; i<3; ++i) n[i]=1.0/Sq3();
00332                 break;
00333             case FailCrit_SMP:
00334                 for (int i=0; i<3; ++i) n[i]=sqrt(Iv[2]/Iv[1])/sqrt(L[i]);
00335                 break;
00336             default:
00337                 throw "FailCrit::Func: Invalid unit normal nType";
00338             }
00339 
00340             // Mobilized Plane projector
00341             Pnn = (n & n);
00342 
00343             // Oblique direction of projection
00344             b   = a*alp + n;
00345             bet = b/(b*n);   // (b*n)   == scalar product
00346             Pbn = (bet & n); // (bet&n) == dyadic
00347 
00348             // "Anisotropic" stress measures on "mobilized" plane
00349             k   = (s % Pnn); // (s%Pnn) == REAL dot
00350             sig = k*bet.Norm(); 
00351             tau = sqrt((ss%Pnn) - 2.0*k*(s%Pbn) + k*k*(bet*bet));
00352 
00353             // Criteria
00354             ff = tau - R*sig - C;
00355 
00356             break;
00357         }
00358         default:
00359             throw "FailCrit::Func: Invalid FailCrit Type";
00360     }
00361     return ff;
00362 } // }}}
00363 
00364 inline REAL FailCrit::AR_Calc_M(REAL sin3th) const // {{{
00365 {
00366     return _AR_Mmax*pow(2.0*_AR_w/(1.0+_AR_w-(1.0-_AR_w)*sin3th), 0.25);
00367 } // }}}
00368 
00369 inline String FailCrit::Name() const // {{{
00370 {
00371     switch (_type)
00372     {
00373         case FailCrit_TR: { return String("Tresca")              ; }
00374         case FailCrit_VM: { return String("von Mises")           ; }
00375         case FailCrit_MC: { return String("Mohr-Coulomb")        ; }
00376         case FailCrit_LD: { return String("Lade-Duncan")         ; }
00377         case FailCrit_DP: { return String("Drucker-Prager")      ; }
00378         case FailCrit_MN: { return String("Matsuoka-Nakai")      ; }
00379         case FailCrit_SM:
00380         {
00381             if (_env_type==FailCrit_NONLIN_1) return String("SMP (Novo) Non-Linear"); 
00382             else                              return String("SMP (Novo)"); 
00383         }
00384         case FailCrit_AR: { return String("Argyris-Sheng et al."); }
00385         case FailCrit_AN: { return String("General Anisotropic") ; }
00386         default:
00387             throw "FailCrit::Name: Invalid FailCrit Type";
00388     }
00389 } // }}}
00390 
00391 inline vtkActor * FailCrit::GenIsoSurf(int nPoints) // {{{
00392 {
00393     // Delete old objects
00394     if (_iso_surf!=NULL) delete _iso_surf;
00395 
00396     // Create a temporary grid around hydrostatic axis
00397     REAL minP=_min_p; REAL maxP=_max_p;
00398     REAL minQ=0.0;    REAL maxQ=_max_p*2.0;
00399     REAL minT=0.0;    REAL maxT=ToRad(360.0);
00400     MeshGrid mg(minP,maxP,nPoints, minQ,maxQ,nPoints, minT,maxT,nPoints);
00401     REAL const * P = mg.X();
00402     REAL const * Q = mg.Y();
00403     REAL const * T = mg.Z();
00404     int       Size = mg.Length();
00405 
00406     // Rotate (convert) P,Q,T values into SI,SII,SIII values
00407     REAL * SI   = new REAL [Size];
00408     REAL * SII  = new REAL [Size];
00409     REAL * SIII = new REAL [Size];
00410     REAL * FF   = new REAL [Size];
00411     Hid2Sig(P,Q,T, SI,SII,SIII, Size);
00412     
00413     // Calculate failure function
00414     for (int i=0; i<Size; ++i) FF[i]=Func(SI[i], SII[i], SIII[i]);
00415 
00416     // Create a VTK structured grid
00417     StructGrid sg(SII,nPoints, SIII,nPoints, SI,nPoints, FF);
00418 
00419     // LookupTable for colors
00420     vtkLookupTable * lt = vtkLookupTable::New();
00421     lt->SetNumberOfColors(1);
00422     lt->Build();
00423     lt->SetTableValue(0,CLR[_iso_clr.GetSTL()].C);
00424 
00425     // Generate the Surfaces
00426     _iso_surf = new SGridIsoSurf (sg.GetGrid(), /*F=0*/0.0, lt);
00427     _iso_surf->GetActor()->GetProperty()->SetOpacity (_iso_opac);
00428 
00429     // Clean up
00430     delete [] SI;
00431     delete [] SII;
00432     delete [] SIII;
00433     delete [] FF;
00434     lt->Delete();
00435 
00436     // Return
00437     return _iso_surf->GetActor();
00438 
00439 } // }}}
00440 
00441 inline void FailCrit::AddActorsTo(VTKWin & Win) // {{{
00442 {
00443     if (_iso_surf!=NULL) Win.AddActor(_iso_surf->GetActor());
00444 } // }}}
00445 
00446 inline void FailCrit::DelActorsFrom(VTKWin & Win) // {{{
00447 {
00448     if (_iso_surf!=NULL) Win.DelActor(_iso_surf->GetActor());
00449 } // }}}
00450 
00451 inline void FailCrit::_calc_constants() // {{{
00452 {
00453     // Mohr-Coulomb sin(Phi)
00454     _sinPhi = sin(_prms.Phi);
00455 
00456     // Tresca and von Mises
00457     if (_prms.Su==0.0) _prms.Su = 3.0*_sinPhi*_prms.pu/(3.0-_sinPhi);
00458 
00459     // Eta constant (Mohr-Coulomb equivalence)
00460     REAL N = 2.0*Sq2()*_sinPhi/(3.0-_sinPhi);
00461     _eta = N;
00462 
00463     // Drucker-Prager
00464     _k_DP = N;
00465 
00466     // Lade-Duncan
00467     _k_LD = 54.0/(Sq2()*N*N*N-3.0*N*N+2.0);
00468 
00469     // Matsuoka-Nakai
00470     _k_MN = (9.0*Sq2()*N+18.0)/(2.0-2.0*N*N+Sq2()*N);
00471 
00472     // SMP (Novo)
00473     REAL k1 = sqrt((2.0+Sq2()*N-2.0*N*N)/(3.0*Sq3()*(Sq2()*N+2.0)));
00474     REAL k2 = sqrt((2.0*N+Sq2())/Sq6());
00475     REAL k3 = sqrt((Sq2()-N)/Sq6());
00476       _k_SM = sqrt((N*N+1.0)/(pow(k2+2.0*k3,2.0)*k1*k1) - 1.0);
00477 
00478     // Argyris-Sheng et al.
00479     _AR_Mmax = N;
00480     _AR_w    = pow((3.0-_sinPhi)/(3.0+_sinPhi),4.0);
00481 
00482 } // }}}
00483 
00484 #endif // MECHSYS_FAILCRIT_H
00485 
00486 // vim:fdm=marker
Generated on Wed Jan 24 15:56:25 2007 for MechSys by  doxygen 1.4.7