el-pl-test.cppGo to the documentation of this file.00001 00002 //################################################################################ 00003 //# COPY-YES (C): :-)) # 00004 //# PROJECT: Object Oriented Finite Element Program # 00005 //# PURPOSE: # 00006 //# CLASS: # 00007 //# # 00008 //# VERSION: # 00009 //# LANGUAGE: C++.ver >= 2.0 # 00010 //# TARGET OS: DOS || UNIX || . . . # 00011 //# DESIGNER(S): Boris Jeremic # 00012 //# PROGRAMMER(S): Boris Jeremic # 00013 //# # 00014 //# # 00015 //# DATE: May 2004 # 00016 //# UPDATE HISTORY: # 00017 //# # 00018 //# # 00019 //################################################################################ 00020 //*/ 00021 #include "time.h" 00022 #include "limits.h" 00023 #include "math.h" 00024 #include "float.h" 00025 // nDarray tools 00026 #include "BJtensor.h" 00027 #include "stresst.h" 00028 #include "straint.h" 00029 #include "BJvector.h" 00030 #include "BJmatrix.h" 00031 //#include "iostream.h" 00032 00033 #include <Tensor.h> 00034 00035 00036 int main(void) 00037 { 00038 ::printf("\n\n------------- MACHINE (CPU) DEPENDENT THINGS --------------\n\n"); 00039 00040 // defining machine epsilon for different built in data types supported by C++ 00041 float float_eps = f_macheps(); // or use float.h values FLT_EPSILON 00042 double double_eps = d_macheps(); // or use float.h values DBL_EPSILON 00043 //long double long_double_eps = ld_macheps(); // or use float.h values LDBL_EPSILON 00044 00045 ::printf("\n float macheps = %.20e \n",float_eps); 00046 ::printf(" double macheps = %.20e \n",double_eps); 00047 //::printf(" long double macheps = %.20Le \n\n\n",long_double_eps); 00048 00049 00050 // Usual tolerance is defined as square root of machine epsilon 00051 float sqrt_f_eps = sqrt(f_macheps()); 00052 double sqrt_d_eps = sqrt(d_macheps()); 00053 //long double sqrt_ld_eps = sqrt(ld_macheps()); 00054 00055 ::printf("\n sqrt float macheps = %.20e \n",sqrt_f_eps); 00056 ::printf(" sqrt double macheps = %.20e \n",sqrt_d_eps); 00057 //::printf(" sqrt long double macheps = %.20Le \n",sqrt_ld_eps); 00058 00059 00060 double one = 1.0; 00061 double onepmc = 1.0+double_eps; 00062 ::printf("\n one = %.32e \n",one); 00063 ::printf("onepmc = %.32e \n",onepmc); 00064 if (one == onepmc) 00065 { 00066 printf("1.0 == 1.0+double_eps -->> OK \n"); 00067 } 00068 else 00069 { 00070 printf("1.0 == 1.0+double_eps -->> NOT OK \n"); 00071 } 00072 00073 00074 00075 double one2 = 10e20; 00076 double onepmc2 = 10e20*(1.0+double_eps); 00077 ::printf("\n one2 = %.32e \n",one2); 00078 ::printf("onepmc2 = %.32e \n",onepmc2); 00079 if (one2 == onepmc2) 00080 { 00081 printf("10e20 == 10e20*(1.0+double_eps) -->> OK \n"); 00082 } 00083 else 00084 { 00085 printf("10e20 == 10e20*(1.0+double_eps) -->> NOT OK \n"); 00086 } 00087 00088 00089 00090 double one3 = 1.0; 00091 double onepmc3 = 1.0+double_eps*10000.0; 00092 ::printf("\n one3 = %.32e \n",one3); 00093 ::printf("onepmc3 = %.32e\n",onepmc3); 00094 if (one3 == onepmc3) 00095 { 00096 printf("1.0 == 1.0+double_eps*10.0 -->> OK \n"); 00097 } 00098 else 00099 { 00100 printf("1.0 == 1.0+double_eps*10.0 -->> NOT OK \n"); 00101 } 00102 00103 00104 00105 ::printf("\n\n----------------- TENSOR --------------------\n\n"); 00106 00107 // Tensor t0(); 00108 // Tensor *ta = new Tensor; 00109 00110 00111 //.................................................................... 00112 // There are several built in tensor types. One of them is: 00113 // Levi-Civita permutation tensor 00114 tensor e("e",3,def_dim_3); 00115 e.print("e","\nLevi-Civita permutation tensor e"); 00116 00117 // The other one is: 00118 // Kronecker delta tensor 00119 tensor I2("I", 2, def_dim_2); 00120 I2.print("I2","\ntensor I2 (2nd order unit tensor: Kronecker Delta -> d_ij)"); 00121 00122 00123 00124 tensor Ipmnq = I2("pm") * I2("nq"); 00125 tensor Imnpq = I2("mn") * I2("pq"); 00126 00127 00128 00130 // tensor X = (1.0/3.0) * Imnpq ;// - Apqmn * (1.0/3.0); 00131 // X.print(); 00132 00133 // some of the 4. order tensor used in conituum mechanics: 00134 // the most general representation of the fourth order isotropic tensor 00135 // includes the following fourth orther unit isotropic tensors: 00136 tensor I_ijkl( 4, def_dim_4, 0.0 ); 00137 I_ijkl = I2("ij")*I2("kl"); 00138 // I_ijkl.print("ijkl","\ntensor I_ijkl = d_ij d_kl)"); 00139 tensor I_ikjl( 4, def_dim_4, 0.0 ); 00140 // I_ikjl = I2("ij")*I2("kl"); 00141 I_ikjl = I_ijkl.transpose0110(); 00142 // I_ikjl.print("ikjl","\ntensor I_ikjl) "); 00143 tensor I_ikjl_inv_2 = I_ikjl.inverse_2(); 00144 // I_ikjl_inv_2.print("I_ikjl","\ntensor I_ikjl_inverse_2)"); 00145 if ( I_ikjl == I_ikjl_inv_2 ) ::printf("\n\n I_ikjl == I_ikjl_inv_2 !\n"); 00146 else ::printf("\a\n\n I_ikjl != I_ikjl_inv_2 !\n"); 00147 00148 tensor I_iljk( 4, def_dim_4, 0.0 ); 00149 // I_ikjl = I2("ij")*I2("kl"); 00150 //.. I_iljk = I_ikjl.transpose0101(); 00151 I_iljk = I_ijkl.transpose0111(); 00152 // I_iljk.print("iljk","\ntensor I_iljk) "); 00153 // To check out this three fourth-order isotropic tensor please see 00154 // W.Michael Lai, David Rubin, Erhard Krempl 00155 // " Introduction to Continuum Mechanics" 00156 // QA808.2 00157 // ISBN 0-08-022699-X 00158 00159 // tensor additions ans scalar multiplications 00160 // symmetric part of fourth oder unit isotropic tensor 00161 tensor I4s = (I_ikjl+I_iljk)*0.5; 00162 // tensor I4s = 0.5*(I_ikjl+I_iljk); 00163 I4s.print("I4s","\n symmetric tensor I4s = (I_ikjl+I_iljk)*0.5 "); 00164 // skew-symmetric part of fourth oder unit isotropic tensor 00165 tensor I4sk = (I_ikjl-I_iljk)*0.5; 00166 I4sk.print("I4sk","\n skew-symmetric tensor I4sk = (I_ikjl-I_iljk)*0.5 "); 00167 00168 //.................................................................... 00169 // Linear Isotropic Elasticity Tensor 00170 // building elasticity tensor 00171 double Ey = 20000; 00172 double nu = 0.2; 00173 // building stiffness tensor in one command line 00174 tensor Eel = (I_ijkl*((Ey*nu*2.0)/(2.0*(1.0+nu)*(1-2.0*nu)))) + 00175 ( (I_ikjl+I_iljk)*(Ey/(2.0*(1.0+nu)))); 00176 // building compliance tensor in one command line 00177 tensor Del= (I_ijkl * (-nu/Ey)) + ( (I_ikjl+I_iljk) * ((1.0+nu)/(2.0*Ey))); 00178 00179 00180 double PIo3 = PI/3.0; 00181 //::printf("Pi/3.0 = %f \n",PIo3); 00182 //------------------------------------------- 00183 double p_start = 5000.000; 00184 double q_start = 2000.0; 00185 double theta_start = PIo3; 00186 //..// double theta_start = 0.4; 00187 00188 stresstensor stress = stress.pqtheta2stress(p_start, q_start, theta_start); 00189 stress.report("stress\n"); 00190 00191 00192 straintensor strain = D2("ijkl") * stress("kl"); 00193 strain.report("strain = D2(\"ijkl\") * stress(\"kl\") \n"); 00194 00195 00196 00197 00198 00199 00200 ::printf("\n\n\n\n\n"); 00201 ::printf("\n\n----------------- Derivatives of Yield function --------\n\n"); 00202 00203 00204 00205 static double stressvalues[] = { 1.0, 2.0, 3.0, 00206 2.0, 2.0, 5.0, 00207 3.0, 5.0, 9.0 }; 00208 00209 stresstensor start_sigma(stressvalues); 00210 start_sigma.report("\n\n starting stress tensor start_sigma )"); 00211 00212 00213 00214 static double strainvalues[] = { 0.1, 0.0, 0.0, 00215 0.0, 0.1, 0.0, 00216 0.0, 0.0, 0.1 }; 00217 straintensor epsilon(strainvalues); 00218 epsilon.report("\n\n strain tensor epsilon (2nd-order with values assignment)"); 00219 00220 00221 stresstensor stress_increment = E("ijpq") * strain_incr("pq"); 00222 stress_increment.null_indices(); 00223 stress_increment.report("\n\n stress increment"); 00224 00225 double f_start = 0.0; 00226 double f_pred = 0.0; 00227 00228 stresstensor elastic_predictor_stress = start_stress + stress_increment; 00229 elastic_predictor_stress.report("\n\n elastic predictor stress"); 00230 00231 00232 // Yield surface evaluation for von Mises 00233 // double k = 10.0: //cohesion 00234 // double k2 = k * k; 00235 // stresstensor s = elastic_predictor_stress.deviator(); 00236 // stresstensor temp1 = s("ij") * s("ij"); 00237 // double temp = temp1.trace(); 00238 // temp = temp * 3.0 / 2.0; 00239 // 00240 // double f = temp - k2; 00241 00242 00243 // f_start = 00244 //::printf("\n############## f_start = %.10e ",f_start); 00245 00246 // f_pred = 00247 //::printf("############## f_pred = %.10e\n\n",f_pred); 00248 00249 00250 stresstensor dFods( 2, def_dim_2, 0.0); 00251 // stresstensor s; // deviator 00252 tensor temp1( 2, def_dim_2, 0.0); 00253 tensor temp2( 2, def_dim_2, 0.0); 00254 double lower = 0.0; 00255 tensor temp3( 2, def_dim_2, 0.0); 00256 00257 double Delta_lambda = 0.0; 00258 00259 00260 stresstensor plastic_stress; 00261 straintensor plastic_strain; 00262 stresstensor elastic_plastic_stress; 00263 // ::printf("\n...... felpred = %lf\n",felpred); 00264 00265 if ( f_pred >= 0 ) 00266 { 00267 00268 // dFods = 00269 00270 // tensor upperE1 = Eel("pqkl")*dQods("kl"); 00271 // upperE1.null_indices(); 00272 // tensor upperE2 = dFods("ij")*Eel("ijmn"); 00273 // upperE2.null_indices(); 00274 // 00275 // tensor upperE = upperE1("pq") * upperE2("mn"); 00276 // upperE.null_indices(); 00277 // 00278 // 00279 // temp2 = upperE2("ij")*dQods("ij"); // L_ij * E_ijkl * R_kl 00280 // temp2.null_indices(); 00281 // lower = temp2.trace(); 00282 // 00283 // tensor Epl = upperE*(1./lower); 00284 // 00285 // tensor Eep = Eel - Epl; 00286 00287 00288 } 00289 00290 ::printf("\n\n\n"); 00291 00292 00293 00294 // -------------- 00295 00296 exit(1); 00297 // return 1; 00298 }
Generated on Mon Oct 23 15:05:24 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents