ElTawil2DUnSym.cppGo to the documentation of this file.00001 // ElTawil2DUnSym.cpp: implementation of the YieldSurfaceBC class. 00002 // 00004 00005 #include "ElTawil2DUnSym.h" 00006 #include <math.h> 00007 00008 #define ELTAWIL2DUNSYM_CLASS_TAG -1 00009 00011 // Construction/Destruction 00013 00014 ElTawil2DUnSym::ElTawil2DUnSym 00015 (int tag, double xPosbal, double yPosbal, 00016 double xNegbal, double yNegbal, 00017 double ypos, double yneg, 00018 YS_Evolution &model, 00019 double cz_pos, double ty_pos, 00020 double cz_neg, double ty_neg) 00021 00022 :YieldSurface_BC2D(tag, ELTAWIL2DUNSYM_CLASS_TAG, 0, 0, model), 00023 xPosBal(xPosbal), yPosBal(yPosbal), 00024 xNegBal(xNegbal), yNegBal(yNegbal), 00025 yPosCap(ypos), yNegCap(yneg), 00026 yPosCap_orig(ypos), yNegCap_orig(yneg), 00027 czPos(cz_pos), tyPos(ty_pos), czNeg(cz_neg), tyNeg(ty_neg), qy(0.005) 00028 { 00029 00030 // This tight coupling with base-class can be a diaster to debug! 00031 if(yPosBal < 0 || yNegBal < 0) 00032 opserr << "WARNING - ElTawil2DUnSym() - yBalance < 0" << endln; 00033 00034 //capY = yPosCap; // why did I change this? 00035 yBal = yPosBal; 00036 00037 if(yNegBal < yBal) 00038 yBal = yNegBal; 00039 00040 // opserr << "yBal= " << yBal << ", yPosBal= " << yPosBal 00041 // << ", yNegBal= " << yNegBal << endl << endln; 00042 00043 // capY = yPosCap - yBal;//!! 00044 capY = yPosCap; 00045 00046 // set to origin 00047 yPosCap -= yBal; 00048 yNegCap -= yBal; 00049 00050 // larger of the two will align with the axis 00051 // the other one will be above 00052 yPosBal -= yBal; 00053 yNegBal -= yBal; 00054 /* 00055 opserr << "yBal= " << yBal << ", yPosBal= " << yPosBal 00056 << ", yNegBal= " << yNegBal << endln; 00057 */ 00058 00059 // set translation 00060 double transY = yBal/capY; 00061 Vector t(2); 00062 t(0) = 0; 00063 //t(1) = yBal/capY; 00064 t(1) = transY; 00065 00066 hModel->setInitTranslation(t); 00067 00068 capX = xPosBal; 00069 if( fabs(xNegBal) > capX) 00070 capX = fabs(xNegBal); 00071 00072 // bad: 00073 capX_orig = capX; 00074 capY_orig = capY; 00075 capXdim = capX; 00076 capYdim = capY; 00077 00078 } 00079 00080 ElTawil2DUnSym::~ElTawil2DUnSym() 00081 { 00082 00083 } 00084 00086 // YS specific methods 00088 void ElTawil2DUnSym::setExtent() 00089 { 00090 // Extent along the axis 00091 // xPos = 1; 00092 // xNeg = -1; 00093 00094 xPos = min_(xPosBal/capX, 00095 (xPosBal*(1 - pow( fabs(yPosBal/(yNegCap - yPosBal)) , tyPos)))/capX); 00096 // xNeg = xNegBal/capX; 00097 xNeg = max_(xNegBal/capX, 00098 (xNegBal*(1 - pow( fabs(yNegBal/(yNegCap - yNegBal)) , tyNeg)))/capX); 00099 00100 yPos = yPosCap/capY - qy; //0.02 0.98; 00101 yNeg = yNegCap/capY + qy; //-0.98; 00102 00103 ytPos = yPos - 0.005; // 0.01 // 0.03 // 0.95 00104 ytNeg = yNeg + 0.005; 00105 00106 double yVal1 = ytPos*capY - yPosBal; 00107 double yVal4 = ytNeg*capY - yPosBal; 00108 double yVal2 = ytPos*capY - yNegBal; 00109 double yVal3 = ytNeg*capY - yNegBal; 00110 00111 /* 00112 char c = ' '; 00113 00114 opserr << "yVals = " << yVal1 << c << yVal2 << c 00115 << yVal3 << c << yVal4 << endln; 00116 00117 */ 00118 double xt1 = xPosBal*(1 - pow((yVal1/(yPosCap - yPosBal)) , czPos)); 00119 double xt4 = xPosBal*(1 - pow( fabs(yVal4/(yNegCap - yPosBal)) , tyPos)); 00120 double xt2 = xNegBal*(1 - pow((yVal2/(yPosCap - yNegBal)) , czNeg)); 00121 double xt3 = xNegBal*(1 - pow( fabs(yVal3/(yNegCap - yNegBal)) , tyNeg)); 00122 00123 00124 xt1 = xt1/capX; 00125 xt2 = xt2/capX; 00126 xt3 = xt3/capX; 00127 xt4 = xt4/capX; 00128 00129 /* 00130 opserr << "xPos = " << xPos << ", xNeg = " << xNeg << endl 00131 << "ytPos = " << ytPos << ", ytNeg = " << ytNeg << endl 00132 << "xt1 = " << xt1 << ", xt2 = " << xt2 << endl 00133 << "xt3 = " << xt3 << ", xt4 = " << xt4 << endln; 00134 */ 00135 } 00136 00137 00138 void ElTawil2DUnSym::getGradient(double &gx, double &gy, double x, double y) 00139 { 00140 // check if the point is on the surface 00141 double drift = getDrift(x, y); 00142 double loc = forceLocation(drift); 00143 double capx = capXdim; 00144 double capy = capYdim; 00145 00146 // opserr << "ElTawil2DUnSym::getGradient drift:" << this->YieldSurface_BC2D::getDrift(x, y) << endln; 00147 00148 if(loc != 0) 00149 { 00150 opserr << "ERROR - ElTawil2D::getGradient(double &gx, double &gy, double x, double y)\n"; 00151 opserr << "Force point not on yield surface, drift = " << drift << " loc = " << loc <<"\n"; 00152 opserr << "\a"; 00153 } 00154 else 00155 { 00156 double a = 10.277;//3.043;//4.29293; //by compatibility, check err in gradient 00157 // double yt = 0.95; 00158 00159 if(y > ytPos) 00160 { 00161 gx = 2*a*x/capx; 00162 gy = 1; 00163 } 00164 else if(y < ytNeg) 00165 { 00166 gx = 2*a*x/capx; 00167 gy = -1; //*y/capy; 00168 } 00169 else 00170 { 00171 double xVal = x*capx; 00172 double yVal = y*capy; 00173 00174 if(xVal >= 0 && yVal >= yPosBal) // quadrant 1 00175 { 00176 gx = 1/xPosBal; 00177 gy = (1/pow( (yPosCap-yPosBal), czPos))*czPos*(pow(yVal-yPosBal,czPos-1)); 00178 } 00179 else if(xVal >=0 && yVal < yPosBal) // quadrant 1 or 4 00180 { 00181 gx = 1/xPosBal; 00182 gy = -1*(1/pow( fabs(yNegCap-yPosBal), tyPos))*tyPos*(pow(fabs(yVal-yPosBal),tyPos-1)); 00183 } 00184 else if(xVal< 0 && yVal >= yNegBal) // quadrant 2 00185 { 00186 gx = 1/xNegBal; // xNegBal should be < 0 00187 gy = (1/pow( (yPosCap-yNegBal), czNeg))*czNeg*(pow(yVal-yNegBal,czNeg-1)); 00188 } 00189 else if(xVal<0 && yVal < yNegBal) // quadrant 2 or 3 00190 { 00191 gx = 1/xNegBal; 00192 gy = -1*(1/pow( fabs(yNegCap-yNegBal), tyNeg))*tyNeg*(pow(fabs(yVal-yNegBal),tyNeg-1)); 00193 } 00194 else 00195 { 00196 opserr << "Eltawil2DUnsym - condition not possible" << endln; 00197 opserr << "\a"; 00198 } 00199 00200 /* gx = 1/xBal; 00201 if(x < 0) 00202 gx = -1*gx; 00203 00204 if(y < 0) 00205 gy = -1*(1/pow( fabs(yNegCap), ty))*ty*(pow(yVal,ty-1)); 00206 else 00207 gy = (1/pow(yPosCap, cz))*cz*(pow(yVal,cz-1)); 00208 */ 00209 } 00210 } 00211 00212 // opserr << "gx = " << gx << "gy = " << gy << "\n"; 00213 // opserr << "\a"; 00214 } 00215 00216 double ElTawil2DUnSym::getSurfaceDrift(double x, double y) 00217 { 00218 double phi; 00219 double a = 5;//10.277; //3.043;//4.29293; --> effects convergence 00220 // why isnt a = 10.77 here or 5 in getGradient 00221 double capx = capX; 00222 double capy = capY; 00223 00224 // ignore the small difference between xt1, xt2 00225 // and xt4, xt3 for determining whether the 00226 // quadratic should be used 00227 if(y > ytPos && fabs(x) < fabs(y*xt1/ytPos) ) 00228 { 00229 phi = a*x*x + y + qy; 00230 } 00231 else if(y < ytNeg && fabs(x) < fabs(y*xt4/ytNeg)) 00232 { 00233 phi = a*x*x - y + qy; 00234 } 00235 else 00236 { 00237 double xVal = x*capx; 00238 double yVal = y*capy; 00239 00240 if(xVal >=0 && yVal >= yPosBal) // quad 1 00241 { 00242 phi = fabs(xVal/xPosBal) + pow((yVal - yPosBal)/(yPosCap - yPosBal), czPos); 00243 } 00244 else if(xVal >=0 && yVal < yPosBal) // quad 1 or 4 00245 { 00246 phi = fabs(xVal/xPosBal) + pow(fabs((yVal - yPosBal)/(yNegCap - yPosBal)), tyPos); 00247 } 00248 else if(xVal < 0 && yVal >= yNegBal)// quad 2 00249 { 00250 phi = fabs(xVal/xNegBal) + pow((yVal - yNegBal)/(yPosCap - yNegBal), czNeg); 00251 } 00252 else if(xVal < 0 && yVal < yNegBal) // quad 2 or 3 00253 { 00254 00255 phi = fabs(xVal/xNegBal) + pow(fabs((yVal - yNegBal)/(yNegCap - yNegBal)), tyNeg); 00256 } 00257 else 00258 { 00259 opserr << "ElTawil2DUnSym::getSurfaceDrift(..) - cond not possible\n"; 00260 opserr << "x=" << x << ", y=" << y << ", capx=" << capx << ", capy=" << capy << endln; 00261 opserr << "xVal = " << xVal << ", yVal = " << yVal << endln; 00262 opserr << "\a"; 00263 } 00264 /* 00265 if(y < 0) 00266 phi = fabs(xVal/xBal) + pow(fabs(yVal/yNegCap), ty); 00267 else 00268 phi = fabs(xVal/xBal) + pow(yVal/yPosCap, cz); 00269 */ 00270 } 00271 00272 double drift = phi - 1; 00273 return drift; 00274 } 00275 00276 // Always call the base class method first 00277 void ElTawil2DUnSym::customizeInterpolate(double &xi, double &yi, double &xj, double &yj) 00278 { 00279 this->YieldSurface_BC2D::customizeInterpolate(xi, yi, xj, yj); 00280 00281 // again, ignore the small difference between xt1, xt2 00282 // and xt4, xt3 00283 00284 if(yj >= ytPos && fabs(xj) <= fabs(yj*xt1/ytPos) ) 00285 { 00286 xi = 0; 00287 yi = 0; 00288 } 00289 else if(yj <= ytNeg && fabs(xj) <= fabs(yj*xt4/ytNeg)) 00290 { 00291 xi = 0; 00292 yi = 0; 00293 } 00294 else 00295 ; // values are okay 00296 00297 00298 } 00299 00300 00301 YieldSurface_BC *ElTawil2DUnSym::getCopy(void) 00302 { 00303 // ElTawil2D *theCopy = new ElTawil2D(this->getTag(), xBal, yBal, yPosCap_orig, yNegCap_orig, *hModel, 00304 // cz, ty); 00305 //later copy all the state variables 00306 00307 00308 ElTawil2DUnSym *theCopy = new ElTawil2DUnSym(this->getTag(), xPosBal, yPosBal+yBal, 00309 xNegBal, yNegBal+yBal, yPosCap+yBal, yNegCap+yBal, 00310 *hModel, czPos, tyPos, czNeg, tyNeg); 00311 00312 return theCopy; 00313 } 00314 00315 int ElTawil2DUnSym::displaySelf(Renderer &theViewer, int displayMode, float fact) 00316 { 00317 this->YieldSurface_BC2D::displaySelf(theViewer, displayMode, fact); 00318 Vector pOld(3), pCurr(3); 00319 Vector rgb(3); 00320 rgb(0) = 0.1; rgb(1) = 0.5; rgb(2) = 0.5; 00321 // rgb(0) = 0.0; rgb(1) = 0.0; rgb(2) = 0.0; 00322 00323 if(displayMode == this->SurfOnly) 00324 { 00325 rgb(0) = 0.7; rgb(1) = 0.7; rgb(2) = 1.0; 00326 //opserr << "Displaying self" << endln; 00327 } 00328 00329 00330 // double incr = 0.101; 00331 double incr = fabs(0.33333333*yNegCap/capY); 00332 // double incr = fabs(0.1*yNegCap/capY); 00333 00334 if(fact < 1) incr = fact; 00335 00336 double xOld = 0; 00337 double yOld = yNegCap/capY; 00338 hModel->toDeformedCoord(xOld, yOld); 00339 00340 double err = 0.01; 00341 // quad 1 and 4 00342 double yc; 00343 for(yc = yNegCap/capY; yc <= yPosCap/capY + err; yc = yc+incr) 00344 //int iter = 0; 00345 //while(1) 00346 { 00347 //opserr << "incr = " << incr; 00348 /*switch (iter) 00349 { 00350 case 0:{yc = yNegCap/capY; break;} 00351 case 1:{yc = 0.333*yNegCap/capY; break;} 00352 case 2:{yc = 0.666*yNegCap/capY; break;} 00353 case 2:{yc = 0.0; break;} 00354 case 3:{yc = yPosBal/capY; break;} 00355 case 4:{yc = 0.5*yPosCap/capY; break;} 00356 case 5:{yc = yPosCap/capY; break;} 00357 default: break; 00358 } 00359 if(iter == 6) 00360 break; 00361 iter++; 00362 */ 00363 00364 double y = yc; 00365 if(y > yPosCap/capY) 00366 y = yPosCap/capY; 00367 00368 double yVal = y*capY; 00369 double xVal; 00370 00371 if(yVal < yPosBal) 00372 { 00373 xVal = xPosBal*(1 - pow( fabs((yVal- yPosBal)/(yNegCap - yPosBal)) , tyPos)); 00374 // xVal = 0; 00375 } 00376 else 00377 { 00378 xVal = xPosBal*(1 - pow( ((yVal - yPosBal)/(yPosCap - yPosBal)) , czPos)); 00379 } 00380 00381 double x = xVal/capX; 00382 00383 if(displayMode==100) 00384 opserr << "(undeformed) x = " << x << ", y = " << y; 00385 00386 hModel->toDeformedCoord(x, y); 00387 if(displayMode==100) 00388 opserr << " (deformed) x = " << x << ", y = " << y << endln; 00389 00390 pCurr(0) = x; 00391 pCurr(1) = y; 00392 pOld(0) = xOld; 00393 pOld(1) = yOld; 00394 00395 // x >= 0 00396 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00397 xOld = x; 00398 yOld = y; 00399 } 00400 00401 xOld = 0; 00402 yOld = yNegCap/capY; 00403 hModel->toDeformedCoord(xOld, yOld); 00404 00405 // opserr << "Plotting negative side\n"; 00406 // xNegBal itself is < 0 00407 // quad 2 and 3 00408 //iter = 0; 00409 for(yc = yNegCap/capY; yc <= yPosCap/capY + err; yc = yc+incr) 00410 //while(1) 00411 { 00412 /*switch (iter) 00413 { 00414 case 0:{yc = yNegCap/capY; break;} 00415 case 1:{yc = 0.5*yNegCap/capY; break;} 00416 case 2:{yc = 0.0; break;} 00417 case 3:{yc = yNegBal/capY; break;} 00418 case 4:{yc = 0.50*yPosCap/capY; break;} 00419 case 5:{yc = yPosCap/capY; break;} 00420 default: break; 00421 } 00422 if(iter == 6) 00423 break; 00424 iter++; 00425 */ 00426 double y = yc; 00427 if(y > yPosCap/capY) 00428 y = yPosCap/capY; 00429 00430 double yVal = y*capY; 00431 double xVal; 00432 00433 if(yVal < yNegBal) 00434 { 00435 xVal = xNegBal*(1 - pow( fabs((yVal-yNegBal)/(yNegCap-yNegBal)) , tyNeg)); 00436 } 00437 else 00438 { 00439 xVal = xNegBal*(1 - pow( ((yVal-yNegBal)/(yPosCap - yNegBal)) , czNeg)); 00440 } 00441 00442 double x = xVal/capX; 00443 00444 if(displayMode==100) 00445 opserr << "(undeformed) x = " << x << ", y = " << y; 00446 00447 hModel->toDeformedCoord(x, y); 00448 if(displayMode==100) 00449 opserr << " (deformed) x = " << x << ", y = " << y << endln; 00450 00451 pCurr(0) = x; 00452 pCurr(1) = y; 00453 pOld(0) = xOld; 00454 pOld(1) = yOld; 00455 00456 // x >= 0 00457 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00458 xOld = x; 00459 yOld = y; 00460 } 00461 00462 //displayForcePoint(theViewer, displayMode, fact); 00463 return 0; 00464 } 00465 00466 void ElTawil2DUnSym::Print(OPS_Stream &s, int flag) 00467 { 00468 s << "\nYield Surface No: " << this->getTag() << " type: ElTawil2DUnSym\n"; 00469 } 00470 00471 00472 /* 00473 // return 0; 00474 00475 Vector pOld(3), pCurr(3); 00476 Vector rgb(3); 00477 rgb(0) = 0.1; rgb(1) = 0.5; rgb(2) = 0.5; 00478 double incr = 0.2;// 0.02, incr = 0.0001; 00479 double x1, y1, xOld, yOld, x2, y2; 00480 double t; 00481 00482 if(fact < 1) incr = fact; 00483 xOld = 0; 00484 //yOld = 1; 00485 00486 t = interpolate(0, 0, 0, 1); 00487 yOld = t*1; 00488 00489 double x, y = 10; 00490 double err = 1e-5; 00491 00492 for(double xc = 0; xc <= xPos+err; xc = xc+incr) // xc <= xPos + err 00493 { 00494 if(xc > xPos) xc = int(xPos); 00495 double yc = sqrt(xPos - xc*xc); // eqn of a circle 2.1 -> 1 00496 00497 t = interpolate(0, 0, xc, yc); //radial return 00498 x = t*xc; 00499 y = t*yc; 00500 00501 if(fact >=1 && x >= 0.64) incr = 0.05; 00502 // if(fact >=1 && x >= 0.84) incr = 0.02; 00503 // if(fact >=1 && x >= 0.94) incr = 0.01; 00504 00505 //if(x < 0.06 || x > 0.9) 00506 { 00507 00509 x1 = x; 00510 y1 = y; 00511 hModel->toDeformedCoord(x1, y1); 00512 if(displayMode==100) 00513 { 00514 opserr << " x = " << x << ", y = " << y << " "; 00515 opserr << " x1 = " << x1 << ", y1 = " << y1 << "\n"; 00516 } 00517 pCurr(0) = x1; 00518 pCurr(1) = y1; 00519 00520 x2 = xOld; 00521 y2 = yOld; 00522 hModel->toDeformedCoord(x2, y2); 00523 pOld(0) = x2; 00524 pOld(1) = y2; 00525 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00526 //opserr << "pOld = " << pOld << " pCurr = " << pCurr << "\n"; 00527 00529 x1 = -1*x; 00530 y1 = y; 00531 hModel->toDeformedCoord(x1, y1); 00532 pCurr(0) = x1; 00533 pCurr(1) = y1; 00534 00535 x2 = -1*xOld; 00536 y2 = yOld; 00537 hModel->toDeformedCoord(x2, y2); 00538 pOld(0) = x2; 00539 pOld(1) = y2; 00540 00541 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00542 00544 t = interpolate(0, 0, xc, -yc); //radial return 00545 x = t*xc; 00546 y = -1*t*yc; 00547 00548 x1 = x; 00549 y1 = y; 00550 hModel->toDeformedCoord(x1, y1); 00551 pCurr(0) = x1; 00552 pCurr(1) = y1; 00553 00554 x2 = xOld; 00555 y2 = -1*yOld; 00556 hModel->toDeformedCoord(x2, y2); 00557 pOld(0) = x2; 00558 pOld(1) = y2; 00559 00560 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00561 00563 x1 = -1*x; 00564 y1 = -1*y; 00565 hModel->toDeformedCoord(x1, y1); 00566 pCurr(0) = x1; 00567 pCurr(1) = y1; 00568 00569 x2 = -1*xOld; 00570 y2 = -1*yOld; 00571 hModel->toDeformedCoord(x2, y2); 00572 pOld(0) = x2; 00573 pOld(1) = y2; 00574 00575 theViewer.drawLine(pOld, pCurr, rgb, rgb); 00576 00577 00578 xOld = x; 00579 yOld = y; 00580 }//x > 0 00581 00582 }//while 00583 */ 00584
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