OpenSees » SRC » material » nD » NewTemplate3Dep
FT3Dep_test.cpp

00001 
00002 //   COPYLEFT (C): Woody's viral GPL-like license (by BJ):
00003 //                 ``This    source  code is Copyrighted in
00004 //                 U.S.,  for  an  indefinite  period,  and anybody
00005 //                 caught  using it without our permission, will be
00006 //                 mighty good friends of ourn, cause we don't give
00007 //                 a  darn.  Hack it. Compile it. Debug it. Run it.
00008 //                 Yodel  it.  Enjoy it. We wrote it, that's all we
00009 //                 wanted to do.''
00010 //
00011 //
00012 // COPYRIGHT (C):     :-))
00013 // PROJECT:           Object Oriented Finite Element Program
00014 // FILE:              
00015 // CLASS:             
00016 // MEMBER FUNCTIONS:
00017 //
00018 // MEMBER VARIABLES
00019 //
00020 // PURPOSE:           
00021 //
00022 // RETURN:
00023 // VERSION:
00024 // LANGUAGE:          C++
00025 // TARGET OS:         
00026 // DESIGNER:          Zhao Cheng, Boris Jeremic
00027 // PROGRAMMER:        Zhao Cheng, 
00028 // DATE:              Fall 2005
00029 // UPDATE HISTORY:    
00030 //
00032 //
00033 
00034 
00035 #include <string.h>
00036 #include <math.h>
00037 #include <stdio.h>
00038 #include <stdlib.h>
00039 
00040 #include "NewTemplate3Dep.h"
00041 
00042 #include "ElasticState.h"
00043 #include "Isotropic_Elastic.h"
00044 #include "elnp_Elastic.h"
00045 #include "DM04_Elastic.h"
00046 
00047 #include "YieldFunction.h"
00048 #include "DP_YF.h"
00049 #include "VM_YF.h"
00050 #include "CC_YF.h"
00051 #include "DM04_YF.h"
00052 
00053 #include "PlasticFlow.h"
00054 #include "DP_PF.h"
00055 #include "VM_PF.h"
00056 #include "CC_PF.h"
00057 #include "DM04_PF.h"
00058 
00059 #include "ScalarEvolution.h"
00060 #include "Linear_Eeq.h"
00061 #include "CC_Ev.h"
00062 
00063 #include "TensorEvolution.h"
00064 #include "Linear_Eij.h"
00065 #include "AF_Eij.h"
00066 #include "DM04_alpha_Eij.h"
00067 #include "DM04_z_Eij.h"
00068 
00069 #include <G3Globals.h>
00070 #include <OPS_Globals.h>
00071 #include <ConsoleErrorHandler.h>
00072 #include <OPS_Stream.h>
00073 
00074 ErrorHandler *g3ErrorHandler =0;
00075 double        ops_Dt = 0;
00076 Domain       *ops_TheActiveDomain = 0;
00077 Element      *ops_TheActiveElement = 0;
00078 
00079 OPS_Stream *opserrPtr;
00080 
00081 int main() {
00082 
00083 opserr << "\n\n\n*** EP Finite Deformations: T E S T ***\n\n" << "\n";
00084 
00085 ofstream outStress ("Results.txt");
00086 
00087 
00088 
00089 double rho = 0.0;
00090 double e0 = 0.735;
00091 double G0 = 125;
00092 double v = 0.05;
00093 double Pat = 100.0;
00094 double kc = 1.0;
00095 double M_cal = 1.25;
00096 double c = 0.712;
00097 double lambda_c = 0.019;
00098 double xi = 0.7;
00099 double er = 0.934;
00100 double m = 0.01;
00101 double h0 = 7.05;
00102 double ch = 0.968;
00103 double nb = 1.1;
00104 double A0 = 0.704;
00105 double nd = 3.5;
00106 double z_max = 4.0;
00107 double cz = 600.0;
00108 stresstensor zeroT1;
00109 stresstensor zeroT2;
00110 stresstensor initStress;
00111 double p = 1000.0;
00112 initStress.val(1,1) = -p; initStress.val(2,2) = -p; initStress.val(3,3) = -p;
00113 //----------------1    2   3  4  5    6   7      8  9         10  11 12 13   14  15  16  17  18    19
00114 double MC[19] = {rho, e0, G0, v, Pat, kc, M_cal, c, lambda_c, xi, er, m, h0, ch, nb, A0, nd, z_max, cz};
00115 stresstensor TS[2] = {zeroT1, zeroT2};
00116 MaterialParameter matpar(MC,19, TS,2);
00117 DM04_Elastic le(3, 4, 5, 6, 2, initStress);
00118 DM04_YF dpy(0,12, 2,1);
00119 DM04_PF dpf(0, 2, 0,11, 0,9, 0,10, 0,5, 0,12, 0,7, 0,8, 0,16, 0,17, 2,1, 2,2 );
00120 DM04_alpha_Eij Eij1(2, 11, 9, 10, 5, 12, 7, 8, 15, 13, 14, 3, 1, 2 );
00121 DM04_z_Eij Eij2(12, 19, 18, 1, 2 );
00122 TensorEvolution *TE[2] = {&Eij1, &Eij2};
00123 
00124 NewTemplate3Dep FTEP(1, &matpar, &le, &dpy, &dpf, TE);
00125 
00126 stresstensor thisE;
00127 stresstensor tStress;
00128 straintensor tStrain;
00129 
00130 double E11 = 0.0;    double E12 = 0.0;    double E13 = 0.0;
00131 double E21 = 0.0;    double E22 = 0.0;    double E23 = 0.0;
00132 double E31 = 0.0;    double E32 = 0.0;    double E33 = 0.0;
00133 
00134 int Num_cyc = 30;
00135 int I_cyc = 0;
00136 
00137 //double ee = 0.0;
00138 double pp = 0.0;
00139 double qq = 0.0;
00140 
00141 double  q_cut = 300.0;
00142 //double ep_cut = 0.015;
00143 
00144 double d_g = 0.00005;
00145  
00146  while ( I_cyc < Num_cyc) {
00147     
00148     //undrained tri-axial compression 
00149     E11 += (- d_g);    E22 += (d_g*0.5);    E33 += (d_g*0.5);
00150     
00151     thisE.val(1,1) = E11; thisE.val(1,2) = E12; thisE.val(1,3) = E13;
00152     thisE.val(2,1) = E21; thisE.val(2,2) = E22; thisE.val(2,3) = E23;
00153     thisE.val(3,1) = E31; thisE.val(3,2) = E32; thisE.val(3,3) = E33;
00154 
00155     FTEP.setTrialStrain(thisE);
00156 
00157     tStress = FTEP.getStressTensor();
00158     tStrain = FTEP.getStrainTensor();
00159     
00160     pp  = tStress.p_hydrostatic();
00161     qq = - tStress.val(1,1) + tStress.val(2,2);
00162 
00163     outStress << -tStrain.val(1,1)  <<  "  "  <<  pp << "  " << qq << endln;
00164 
00165     //if (fabs(E11) >= ep_cut) {
00166     if (fabs(qq) >= q_cut) {
00167         I_cyc++;
00168         d_g *= (-1.0);
00169     }
00170 
00171  }
00172 
00173   return 0;
00174 
00175 }
00176 
00177 
00178 /*
00179 double rho = 0.0;
00180 double e0 = 0.735;
00181 double G0 = 125;
00182 double v = 0.25;
00183 double Pat = 100.0;
00184 double kc = 0.1;
00185 double M_cal = 1.25;
00186 double c = 0.712;
00187 double lambda_c = 0.019;
00188 double xi = 0.7;
00189 double er = 0.934;
00190 double m = 0.01;
00191 double h0 = 7.05;
00192 double ch = 0.968;
00193 double nb = 1.1;
00194 double A0 = 0.704;
00195 double nd = 3.5;
00196 double z_max = 4.0;
00197 double cz = 600.0;
00198 stresstensor zeroT1;
00199 stresstensor zeroT2;
00200 stresstensor initStress;
00201 double p = 100;
00202 initStress.val(1,1) = -p; initStress.val(2,2) = -p; initStress.val(3,3) = -p;
00203 //----------------1    2   3  4  5    6   7      8  9         10  11 12 13   14  15  16  17  18    19
00204 double MC[19] = {rho, e0, G0, v, Pat, kc, M_cal, c, lambda_c, xi, er, m, h0, ch, nb, A0, nd, z_max, cz};
00205 stresstensor TS[2] = {zeroT1, zeroT2};
00206 MaterialParameter inpar(MC,19, TS,2);
00207 DM04_Elastic le(3, 4, 5, 6, 2, initStress);
00208 DM04_YF dpy(0,12, 2,1);
00209 DM04_PF dpf(0,2, 0,11, 0,9, 0,10, 0,5, 0,12, 0,7, 0,8, 0,16, 0,17, 2,1, 2,2 );
00210 DM04_alpha_Eij Eij1(2, 11, 9, 10, 5, 12, 7, 8, 15, 13, 14, 3, 1, 2 );
00211 DM04_z_Eij Eij2(12, 19, 18, 1, 2 );
00212 TensorEvolution *TE[2] = {&Eij1, &Eij2};
00213 
00214 NewTemplate3Dep FTEP(1, &inpar, &le, &dpy, &dpf, TE);
00215 
00216 stresstensor thisE;
00217 stresstensor tStress;
00218 straintensor tStrain;
00219 
00220 double E11 = 0.0;    double E12 = 0.0;    double E13 = 0.0;
00221 double E21 = 0.0;    double E22 = 0.0;    double E23 = 0.0;
00222 double E31 = 0.0;    double E32 = 0.0;    double E33 = 0.0;
00223 
00224   //double ee = 0.0;
00225   double pp = 0.0;
00226   double qq = 0.0;
00227 
00228 int Num_cyc = 20;
00229 int I_cyc = 0;
00230 
00231 //double  q_cut = 5.0;
00232 double ep_cut = 0.0005;  
00233   
00234 double d_g = 0.000001;
00235     
00236 //for (int i = 0; i < 600; i++) {
00237 while ( I_cyc < Num_cyc) {
00238  
00239     //undrained tri-axial compression 
00240     E11 += (- d_g);
00241     
00242     thisE.val(1,1) = E11; thisE.val(1,2) = E12; thisE.val(1,3) = E13;
00243     thisE.val(2,1) = E21; thisE.val(2,2) = E22; thisE.val(2,3) = E23;
00244     thisE.val(3,1) = E31; thisE.val(3,2) = E32; thisE.val(3,3) = E33;
00245 
00246     FTEP.setTrialStrain(thisE);
00247 
00248     tStress = FTEP.getStressTensor();
00249     tStrain = FTEP.getStrainTensor();
00250     
00251     //ee = FTEP.getVoid();    
00252     pp  = tStress.p_hydrostatic();
00253     //qq  = tStress.q_deviatoric();
00254     qq = - tStress.val(1,1) + tStress.val(2,2);    
00255 
00256     outStress << -tStrain.val(1,1)  <<  "  "  <<  pp << "  " << qq << endlnn;
00257     
00258     if (fabs(E11) >= ep_cut) {
00259     //if (fabs(qq) >= q_cut) {
00260         I_cyc++;
00261         d_g *= (-1.0);
00262     }    
00263     
00264 
00265 }
00266     
00267   return 0;
00268 
00269 }
00270 */
00271 
00272 
00273 
00274 /*
00275  // Test for Von_Mises
00276  double rho = 0.0;
00277  double E = 1.0e4;
00278  double v = 0.25;
00279  double k = 4.0;
00280  double H = 0.0;
00281  double Hij = 500.0;
00282  double ha = 5000.0;
00283  double Cr = 1000.0;
00284  stresstensor zero;
00285  double MC[7] = {rho, E, v, H, Hij, ha, Cr};
00286  double IS[1] = {k};
00287  stresstensor TS[1] = {zero};
00288  MaterialParameter matpar(MC,7, IS,1, TS,1);
00289  Isotropic_Elastic le(2, 3);
00290  VM_YF dpy(1,1, 2,1);
00291  VM_PF dpf(2,1);
00292  Linear_Eeq Eep(4);
00293  //Linear_Eij Eij(5);
00294  AF_Eij Eij(6, 7, 1);
00295  ScalarEvolution *SE = {&Eep};
00296  TensorEvolution *TE = {&Eij};
00297 
00298  NewTemplate3Dep FTEP(1, &matpar, &le, &dpy, &dpf, &SE, &TE);
00299 
00300  int Num_step = 200;
00301  double d_g = 0.00005;
00302 
00303  stresstensor thisE;
00304  stresstensor tStress;
00305  straintensor tStrain;
00306  
00307  double pp = 0.0;
00308  double qq = 0.0;
00309 
00310  tensor aStress(2, def_dim_2, 0.0);
00311 
00312  for ( int i = 0; i <= Num_step; i++ )
00313  {
00314 
00315  //Uniaxial Loading
00316  double E11 = d_g*float(i);    double E12 = 0.0;    double E13 = 0.0;
00317  double E21 = 0.0;    double E22 = -d_g*float(i);    double E23 = 0.0;
00318  double E31 = 0.0;    double E32 = 0.0;    double E33 = -d_g*float(i);
00319  //
00320 
00321  thisE.val(1,1) = E11; thisE.val(1,2) = E12; thisE.val(1,3) = E13;
00322  thisE.val(2,1) = E21; thisE.val(2,2) = E22; thisE.val(2,3) = E23;
00323  thisE.val(3,1) = E31; thisE.val(3,2) = E32; thisE.val(3,3) = E33;
00324  // thisE.print("F","F:");
00325 
00326  FTEP.setTrialStrain(thisE);
00327  
00328  tStress = FTEP.getStressTensor();
00329  tStrain = FTEP.getStrainTensor();
00330     
00331  pp  = tStress.p_hydrostatic();
00332  qq = tStress.val(1,1) - tStress.val(2,2);
00333 
00334  outStress << -tStrain.val(1,1)  <<  "  "  <<  pp << "  " << qq << endlnn;
00335 
00336  }
00337 
00338  return 0;
00339  }
00340 */
00341 
00342 
00343 /*
00344 
00345  // Test for Cam Clay 
00346 
00347  double rho = 0.0;
00348  double e0 = 0.8;
00349  double M = 0.8;
00350  double lambda = 0.15;
00351  double kappa = 0.05;
00352  double v = 0.25;
00353  double Kc = 200.0;
00354  double p0 = 200.0;
00355  stresstensor initStress;
00356  initStress.val(1,1) = -p0/1.1; initStress.val(2,2) = -p0/1.1; initStress.val(3,3) = -p0/1.1;
00357 
00358  stresstensor zero;
00359  double MC[7] = {rho, e0, M, lambda, kappa, v, Kc};
00360  double IS[1] = {p0};
00361  MaterialParameter inpar(MC,7, IS,1);
00362  elnp_Elastic le(5, 6, 7, 2, initStress);
00363  CC_YF dpy(0,3, 1,1);
00364  CC_PF dpf(0,3, 1,1);
00365  CC_Ev Ev(4, 5, 2, 1);
00366  ScalarEvolution *SE = {&Ev};
00367                                                                             
00368  NewTemplate3Dep FTEP(1, &inpar, &le, &dpy, &dpf, &SE);
00369 
00370  int Num_step = 1000;
00371  double d_g = -0.0001;
00372 
00373  stresstensor thisE;
00374  stresstensor tStress;
00375  straintensor tStrain;
00376  
00377  double pp = 0.0;
00378  double qq = 0.0;
00379 
00380  tensor aStress(2, def_dim_2, 0.0);
00381 
00382  for ( int i = 0; i <= Num_step; i++ )
00383  {
00384 
00385  //Uniaxial Loading
00386  double E11 = d_g*float(i);    double E12 = 0.0;    double E13 = 0.0;
00387  double E21 = 0.0;    double E22 = -d_g*float(i)*0.5;    double E23 = 0.0;
00388  double E31 = 0.0;    double E32 = 0.0;    double E33 = -d_g*float(i)*0.5;
00389  //
00390 
00391  thisE.val(1,1) = E11; thisE.val(1,2) = E12; thisE.val(1,3) = E13;
00392  thisE.val(2,1) = E21; thisE.val(2,2) = E22; thisE.val(2,3) = E23;
00393  thisE.val(3,1) = E31; thisE.val(3,2) = E32; thisE.val(3,3) = E33;
00394 
00395  FTEP.setTrialStrain(thisE);
00396  
00397  tStress = FTEP.getStressTensor();
00398  tStrain = FTEP.getStrainTensor();
00399     
00400  pp  = tStress.p_hydrostatic();
00401  qq = -tStress.val(1,1) + tStress.val(2,2);
00402 
00403  outStress << -tStrain.val(1,1)  <<  "  "  <<  pp << "  " << qq << endlnn;
00404 
00405  }
00406 
00407  return 0;
00408  }
00409 */
00410 
00411 
00412 //  //*********************   Simple Shear ******************************
00413 //  double E11 = -0.001;    double E12 = 0.0;    double E13 = d_g*float(i);
00414 //  double E21 = 0.0;    double E22 = -0.001;    double E23 = 0.0;
00415 //  double E31 = d_g*float(i);    double E32 = 0.0;    double E33 = -0.001;
00416 //  //*******************************************************************
00417 
00418  
00419 // /***********************************************************************************/
00420 // /*   Uniaxial                                                                      */
00421 // double E11 = d_g*float(i);    double E12 = 0.0;    double E13 = 0.0;
00422 // double E21 = 0.0;    double E22 = d_g*float(i)*(-v);    double E23 = 0.0;
00423 // double E31 = 0.0;    double E32 = 0.0;    double E33 = d_g*float(i)*(-v); 
00424 
00425 // /***********************************************************************************/
00426 // /*   Triaxial Compression                                                          */
00427 // double E11 = 0.0;    double E12 = 0.0;    double E13 = 0.0;
00428 // double E21 = 0.0;    double E22 = 0.0;    double E23 = 0.0;
00429 // double E31 = 0.0;    double E32 = 0.0;    double E33 = -d_g*float(i);
00430 
00432 //double E11 = d_g*float(i);    double E12 = 0.0;    double E13 = 0.0;
00433 //double E21 = 0.0;    double E22 = d_g*float(i)*(-v);    double E23 = 0.0;
00434 //double E31 = 0.0;    double E32 = 0.0;    double E33 = d_g*float(i)*(-v);
Generated on Mon Oct 23 15:05:15 2006 for OpenSees by 1.5.0
opensees-support @ berkeley.edu
Supported by the National Science Foundation