ElTawil2D.cppGo to the documentation of this file.00001 // ElTawil2D.cpp: implementation of the YieldSurfaceBC class. 00002 // 00004 00005 #include "ElTawil2D.h" 00006 #include <math.h> 00007 00008 #define ELTAWIL2D_CLASS_TAG -1 00009 00011 // Construction/Destruction 00013 00014 ElTawil2D::ElTawil2D(int tag, double xbal, double ybal, double ypos, double yneg, 00015 YS_Evolution &model, double cz_, double ty_) 00016 :YieldSurface_BC2D(tag, ELTAWIL2D_CLASS_TAG, 0, 0, model), 00017 xBal(xbal), yBal(ybal), yPosCap(ypos), yNegCap(yneg), cz(cz_), 00018 yPosCap_orig(ypos), yNegCap_orig(yneg), ty(ty_), qy(0.005) //0.01 00019 { 00020 capY = yPosCap; 00021 // capY = yPosCap - yBal; 00022 00023 // set to origin 00024 yPosCap -= yBal; 00025 yNegCap -= yBal; 00026 00027 // set translation 00028 double transY = yBal/capY; 00029 Vector t(2); 00030 t(0) = 0; 00031 //t(1) = yBal/capY; 00032 t(1) = transY; 00033 hModel->setInitTranslation(t); 00034 00035 // opserr << "Translation Y = " << transY << endln; 00036 00037 //capX = xBal*(1 - pow( fabs((transY*capY)/yNegCap) , ty)); 00038 capX = xBal; 00039 00040 // bad: 00041 capX_orig = capX; 00042 capY_orig = capY; 00043 capXdim = capX; 00044 capYdim = capY; 00045 00046 } 00047 00048 ElTawil2D::~ElTawil2D() 00049 { 00050 00051 } 00052 00054 // YS specific methods 00056 void ElTawil2D::setExtent() 00057 { 00058 // Extent along the axis 00059 // xPos = 1; 00060 // xNeg = -1; 00061 00062 xPos = xBal/capX; 00063 xNeg = -xPos; 00064 00065 yPos = yPosCap/capY - qy; //0.02 0.98; 00066 yNeg = yNegCap/capY + qy; //-0.98; 00067 00068 ytPos = yPos - 0.005; // 0.01 // 0.03 // 0.95 00069 ytNeg = yNeg + 0.005; 00070 00071 double yValPos = ytPos*capY; 00072 double yValNeg = ytNeg*capY; 00073 00074 xtPos = xBal*(1 - pow((yValPos/yPosCap) , cz)); 00075 xtNeg = xBal*(1 - pow( fabs(yValNeg/yNegCap) , ty)); 00076 00077 xtPos = xtPos/capX; 00078 xtNeg = xtNeg/capX; 00079 /* 00080 opserr << "ytPos = " << ytPos << ", ytNeg = " << ytNeg << ", xtPos = " << xtPos 00081 << ", xtNeg = " << xtNeg << endln; 00082 */ 00083 } 00084 00085 00086 void ElTawil2D::getGradient(double &gx, double &gy, double x, double y) 00087 { 00088 // check if the point is on the surface 00089 double drift = getDrift(x, y); 00090 double loc = forceLocation(drift); 00091 double capx = capXdim; 00092 double capy = capYdim; 00093 00094 if(loc != 0) 00095 { 00096 opserr << "ERROR - ElTawil2D::getGradient(double &gx, double &gy, double x, double y)\n"; 00097 opserr << "Force point not on yield surface, drift = " << drift << " loc = " << loc <<"\n"; 00098 // opserr << "\a"; 00099 gx = 1.0; 00100 gy = 1.0; 00101 } 00102 else 00103 { 00104 double a = 10.277;//3.043;//4.29293; //by compatibility, check err in gradient 00105 // double yt = 0.95; 00106 00107 if(y > ytPos) 00108 { 00109 gx = 2*a*x/capx; 00110 gy = 1; 00111 } 00112 else if(y < ytNeg) 00113 { 00114 gx = 2*a*x/capx; 00115 gy = -1; //-1*y/capy; <-- why did i write this?// -1 ? 00116 } 00117 else 00118 { 00119 // double xVal = x*capx; 00120 double yVal = fabs(y*capy); 00121 // double yVal = y*capy; //!! 00122 00123 gx = 1/xBal; 00124 if(x < 0) 00125 gx = -1*gx; 00126 00127 if(y < 0) 00128 gy = -1*(1/pow( fabs(yNegCap), ty))*ty*(pow(yVal,ty-1)); 00129 else 00130 gy = (1/pow(yPosCap, cz))*cz*(pow(yVal,cz-1)); 00131 } 00132 } 00133 00134 // opserr << "gx = " << gx << ", gy = " << gy << ", capy = " << capy << "\n"; 00135 // opserr << "\a"; 00136 } 00137 00138 double ElTawil2D::getSurfaceDrift(double x, double y) 00139 { 00140 double phi; 00141 double a = 5;//10.277; //3.043;//4.29293; --> effects convergence 00142 00143 00144 if(y > ytPos && fabs(x) < fabs(y*xtPos/ytPos) ) 00145 { 00146 phi = a*x*x + y + qy; 00147 } 00148 else if(y < ytNeg && fabs(x) < fabs(y*xtNeg/ytNeg)) 00149 { 00150 phi = a*x*x - y + qy; 00151 } 00152 else 00153 { 00154 double xVal = x*capX; 00155 double yVal = y*capY; 00156 00157 if(y < 0) 00158 phi = fabs(xVal/xBal) + pow(fabs(yVal/yNegCap), ty); 00159 else 00160 phi = fabs(xVal/xBal) + pow(yVal/yPosCap, cz); 00161 } 00162 00163 double drift = phi - 1; 00164 00165 // opserr << "Eltawil2D - x = " << x << ", y = " << y << ", drift = " << drift << endln; 00166 00167 return drift; 00168 } 00169 00170 // Always call the base class method first 00171 void ElTawil2D::customizeInterpolate(double &xi, double &yi, double &xj, double &yj) 00172 { 00173 this->YieldSurface_BC2D::customizeInterpolate(xi, yi, xj, yj); 00174 00175 if(yj >= ytPos && fabs(xj) <= fabs(yj*xtPos/ytPos) ) 00176 { 00177 xi = 0; 00178 yi = 0; 00179 } 00180 else if(yj <= ytNeg && fabs(xj) <= fabs(yj*xtNeg/ytNeg)) 00181 { 00182 xi = 0; 00183 yi = 0; 00184 } 00185 else 00186 ; // values are okay 00187 00188 00189 } 00190 00191 00192 YieldSurface_BC *ElTawil2D::getCopy(void) 00193 { 00194 ElTawil2D *theCopy = new ElTawil2D(this->getTag(), xBal, yBal, yPosCap_orig, yNegCap_orig, *hModel, 00195 cz, ty); 00196 //later copy all the state variables 00197 return theCopy; 00198 } 00199 00200 int ElTawil2D::displaySelf(Renderer &theViewer, int displayMode, float fact) 00201 { 00202 this->YieldSurface_BC2D::displaySelf(theViewer, displayMode, fact); 00203 Vector pOld(3), pCurr(3); 00204 Vector rgb(3); 00205 rgb(0) = 0.1; rgb(1) = 0.5; rgb(2) = 0.5; 00206 if(displayMode == this->SurfOnly) 00207 { 00208 rgb(0) = 0.7; rgb(1) = 0.7; rgb(2) = 1.0; 00209 } 00210 00211 00212 // double incr = ((yPosCap/capY) - (yNegCap/capY)) / 10; 00213 double incr = fabs(0.33333333*yNegCap/capY); 00214 if(fact < 1) incr = fact; 00215 00216 double xOld = 0; 00217 double yOld = yNegCap/capY; 00218 // hModel->toDeformedCoord(xOld, yOld); 00219 00220 double err = 1e-4; 00221 for(double yc = yNegCap/capY; yc <= yPosCap/capY + err; yc = yc+incr) 00222 { 00223 double y = yc; 00224 double yVal = y*capY; 00225 double xVal; 00226 00227 if(y < 0) 00228 { 00229 xVal = xBal*(1 - pow( fabs(yVal/yNegCap) , ty)); 00230 } 00231 else 00232 { 00233 xVal = xBal*(1 - pow( (yVal/yPosCap) , cz)); 00234 } 00235 00236 double x = xVal/capX; 00237 00238 if(displayMode==100) 00239 opserr << "(undeformed) x = " << x << ", y = " << y; 00240 00241 double x1 = x; 00242 double y1 = y; 00243 double x2 = -x; 00244 double y2 = y; 00245 00246 double x1Old = xOld; 00247 double y1Old = yOld; 00248 double x2Old = -xOld; 00249 double y2Old = yOld; 00250 00251 00252 hModel->toDeformedCoord(x1, y1); 00253 hModel->toDeformedCoord(x1Old, y1Old); 00254 hModel->toDeformedCoord(x2, y2); 00255 hModel->toDeformedCoord(x2Old, y2Old); 00256 00257 // if(displayMode==100) 00258 // opserr << " (deformed) x = " << x << ", y = " << y << endln; 00259 00260 pCurr(0) = x1; 00261 pCurr(1) = y1; 00262 pOld(0) = x1Old; 00263 pOld(1) = y1Old; 00264 00265 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00266 00267 pCurr(0) = x2; 00268 pCurr(1) = y2; 00269 pOld(0) = x2Old; 00270 pOld(1) = y2Old; 00271 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00272 00273 xOld = x; 00274 yOld = y; 00275 } 00276 00277 00278 // displayForcePoint(theViewer, displayMode, fact); 00279 return 0; 00280 } 00281 00282 void ElTawil2D::Print(OPS_Stream &s, int flag) 00283 { 00284 s << "\nYield Surface No: " << this->getTag() << " type: ElTawil2D\n"; 00285 this->YieldSurface_BC::Print(s, flag); 00286 } 00287 00288 00289 /* 00290 // return 0; 00291 00292 Vector pOld(3), pCurr(3); 00293 Vector rgb(3); 00294 rgb(0) = 0.1; rgb(1) = 0.5; rgb(2) = 0.5; 00295 double incr = 0.2;// 0.02, incr = 0.0001; 00296 double x1, y1, xOld, yOld, x2, y2; 00297 double t; 00298 00299 if(fact < 1) incr = fact; 00300 xOld = 0; 00301 //yOld = 1; 00302 00303 t = interpolate(0, 0, 0, 1); 00304 yOld = t*1; 00305 00306 double x, y = 10; 00307 double err = 1e-5; 00308 00309 for(double xc = 0; xc <= xPos+err; xc = xc+incr) // xc <= xPos + err 00310 { 00311 if(xc > xPos) xc = int(xPos); 00312 double yc = sqrt(xPos - xc*xc); // eqn of a circle 2.1 -> 1 00313 00314 t = interpolate(0, 0, xc, yc); //radial return 00315 x = t*xc; 00316 y = t*yc; 00317 00318 if(fact >=1 && x >= 0.64) incr = 0.05; 00319 // if(fact >=1 && x >= 0.84) incr = 0.02; 00320 // if(fact >=1 && x >= 0.94) incr = 0.01; 00321 00322 //if(x < 0.06 || x > 0.9) 00323 { 00324 00326 x1 = x; 00327 y1 = y; 00328 hModel->toDeformedCoord(x1, y1); 00329 if(displayMode==100) 00330 { 00331 opserr << " x = " << x << ", y = " << y << " "; 00332 opserr << " x1 = " << x1 << ", y1 = " << y1 << "\n"; 00333 } 00334 pCurr(0) = x1; 00335 pCurr(1) = y1; 00336 00337 x2 = xOld; 00338 y2 = yOld; 00339 hModel->toDeformedCoord(x2, y2); 00340 pOld(0) = x2; 00341 pOld(1) = y2; 00342 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00343 //opserr << "pOld = " << pOld << " pCurr = " << pCurr << "\n"; 00344 00346 x1 = -1*x; 00347 y1 = y; 00348 hModel->toDeformedCoord(x1, y1); 00349 pCurr(0) = x1; 00350 pCurr(1) = y1; 00351 00352 x2 = -1*xOld; 00353 y2 = yOld; 00354 hModel->toDeformedCoord(x2, y2); 00355 pOld(0) = x2; 00356 pOld(1) = y2; 00357 00358 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00359 00361 t = interpolate(0, 0, xc, -yc); //radial return 00362 x = t*xc; 00363 y = -1*t*yc; 00364 00365 x1 = x; 00366 y1 = y; 00367 hModel->toDeformedCoord(x1, y1); 00368 pCurr(0) = x1; 00369 pCurr(1) = y1; 00370 00371 x2 = xOld; 00372 y2 = -1*yOld; 00373 hModel->toDeformedCoord(x2, y2); 00374 pOld(0) = x2; 00375 pOld(1) = y2; 00376 00377 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00378 00380 x1 = -1*x; 00381 y1 = -1*y; 00382 hModel->toDeformedCoord(x1, y1); 00383 pCurr(0) = x1; 00384 pCurr(1) = y1; 00385 00386 x2 = -1*xOld; 00387 y2 = -1*yOld; 00388 hModel->toDeformedCoord(x2, y2); 00389 pOld(0) = x2; 00390 pOld(1) = y2; 00391 00392 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00393 00394 00395 xOld = x; 00396 yOld = y; 00397 }//x > 0 00398 00399 }//while 00400 */ 00401
Generated on Mon Oct 23 15:05:23 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents