RMC01.cppGo to the documentation of this file.00001 //================================================================================ 00002 //# COPY LEFT and RIGHT: # 00003 //# Commercial use of this program without express permission of the # 00004 //# University of California, is strictly encouraged. Copyright and Copyleft # 00005 //# are covered by the following clause: # 00006 //# # 00007 //# Woody's license: # 00008 //# ``This source code is Copyrighted in U.S., by the The Regents of the # 00009 //# University of California, for an indefinite period, and anybody caught # 00010 //# using it without our permission, will be mighty good friends of ourn, # 00011 //# cause we don't give a darn. Hack it. Compile it. Debug it. Run it. Yodel # 00012 //# it. Enjoy it. We wrote it, that's all we wanted to do.'' bj # 00013 //# # 00014 //# # 00015 //# # 00016 //# PROJECT: Object Oriented Finite Element Program # 00017 //# PURPOSE: Rounded Mohr Coulomb base functions # 00018 //# CLASS: # 00019 //# # 00020 //# VERSION: # 00021 //# LANGUAGE: C++ # 00022 //# TARGET OS: DOS || UNIX || . . . # 00023 //# DESIGNER(S): Boris Jeremic jeremic@ucdavis.edu # 00024 //# Zhao Cheng, # 00025 //# Zhaohui Yang # 00026 //# PROGRAMMER(S): Zhao Cheng, Boris Jeremic # 00027 //# # 00028 //# # 00029 //# DATE: 12 Feb. 2003 # 00030 //# UPDATE HISTORY: # 00031 //# # 00032 //# # 00033 //# # 00034 //# # 00035 //# # 00036 //################################################################################ 00037 //*/ 00038 00039 #include "RMC01.h" 00040 00041 //############################################################################# 00042 double g_0( double theta, double e ) 00043 { 00044 double g = 0.0; 00045 00046 double ct = cos(theta); 00047 double e_ = e; 00048 00049 double temp3 = 2.0 * e_ - 1.0; 00050 double temp4 = 4.0 * ( 1.0 - e_ * e_ ) * ct * ct; 00051 00052 double upper = temp4 + temp3 * temp3; 00053 double lower = 00054 2.0 * ( 1.0 - e_ *e_ ) * ct + 00055 temp3 *sqrt ( temp4 + 5.0 * e_ * e_ - 4. * e_ ); 00056 00057 g = upper/lower; 00058 return g; 00059 } 00060 00061 //############################################################################# 00062 double g_prime( double theta, double e ) 00063 { 00064 double g_prime = 0.0; 00065 00066 double ct = cos(theta); 00067 double st = sin(theta); 00068 double e_ = e; 00069 00070 double temp2 = 1.0 - e_ * e_; 00071 double temp3 = 2.0 * e_ - 1.0; 00072 double temp4 = 4.0 * temp2 * ct * ct; 00073 00074 double N = temp4 + temp3*temp3; 00075 00076 double tempsqrt = temp4 + 5.0*e_*e_ - 4.0 * e_ ; 00077 // double EPS = sqrt(d_macheps()); 00078 // this is because it might be close 00079 // to zero ( -1e-19 ) numericaly 00080 // but sqrt() does not accept it 00081 // if ( tempsqrt < 0.0 || fabs(tempsqrt) < EPS ) 00082 // { 00083 //::fprintf(stdout,"tempsqrt < 0.0 || fabs(tempsqrt) < EPS in "); 00084 //::fprintf(stdout,"double Material_Model::g_Willam_Warnke_prime( double theta, double e ) 00085 //const \a\a\n"); 00086 //::fprintf(stderr,"tempsqrt < 0.0 || fabs(tempsqrt) < EPS in "); 00087 //::fprintf(stderr,"double Material_Model::g_Willam_Warnke_prime( double theta, double e ) 00088 //const \a\a\n"); 00089 // ::exit(1); 00090 // } 00091 00092 00093 double temp1 = sqrt( tempsqrt ); 00094 double D = 2.0*temp2*ct + temp3*temp1; 00095 00096 double N_prime = -8.0*temp2*ct*st; 00097 00098 double temp5 = 2.0*temp2*st; 00099 double temp6 = 4.0*temp3*temp2*ct*st; 00100 double D_prime = - temp5 - temp6/temp1; 00101 00102 g_prime = - N*D_prime/(D*D) + N_prime/D; 00103 00104 return g_prime; 00105 } 00106 00107 //############################################################################# 00108 double g_second( double theta, double e ) 00109 { 00110 double g_second = 0.0; 00111 00112 double ct = cos(theta); 00113 double st = sin(theta); 00114 double e_ = e; 00115 00116 00117 double temp2 = 1.0-e_*e_; 00118 double temp3 = 2.0*e_-1.0; 00119 double temp4 = 4.0*temp2*ct*ct; 00120 00121 double N = temp4 + temp3*temp3; 00122 00123 double tempsqrt = temp4+5.0*e_*e_-4.0*e_ ; 00124 // double EPS = sqrt(d_macheps()); 00125 // this is because it might be close 00126 // to zero ( -1e-19 ) numericaly 00127 // but sqrt() does not accept it 00128 // if ( tempsqrt < 0.0 || fabs(tempsqrt) < EPS ) 00129 // { 00130 //::fprintf(stdout,"tempsqrt < 0.0 || fabs(tempsqrt) < EPS in "); 00131 //::fprintf(stdout," double Material_Model::g_Willam_Warnke_second( double theta, double e ) 00132 //const \a\a\n"); 00133 //::fprintf(stderr,"tempsqrt < 0.0 || fabs(tempsqrt) < EPS in "); 00134 //::fprintf(stderr," double Material_Model::g_Willam_Warnke_second( double theta, double e ) 00135 //const \a\a\n"); 00136 // ::exit(1); 00137 // } 00138 00139 double temp1 = sqrt(tempsqrt); 00140 double D = 2.0*temp2*ct + temp3*temp1; 00141 00142 double N_prime = -8.0*temp2*ct*st; 00143 00144 double temp5 = 2.0*temp2*st; 00145 double temp6 = 4.0*temp3*temp2*ct*st; 00146 double D_prime = - temp5 - temp6/temp1; 00147 00148 double N_second = - 8.0*temp2*ct*ct + 8.0*temp2*st*st; 00149 00150 double temp7 = 4.0*temp3*temp2*ct*ct; 00151 double temp8 = 16.0*temp3*temp2*temp2*ct*ct*st*st; 00152 double temp9 = 4.0*temp3*temp2*st*st; 00153 double temp10 = 2.0*temp2*ct; 00154 double D_second = - temp10 - 00155 temp7/temp1 - 00156 temp8/(temp1*temp1*temp1) + 00157 temp9/temp1; 00158 00159 g_second = 2.0 * N * D_prime * D_prime / (D*D*D) - 00160 2.0 * D_prime * N_prime / (D*D) - 00161 N * D_second / (D*D) + 00162 N_second / D; 00163 00164 return g_second; 00165 }
Generated on Mon Oct 23 15:05:16 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents