SoilFootingSection2d.cppGo to the documentation of this file.00001 // Date: 12/08/2004 00002 // 00003 // File name: SoilFootingSection2d.cpp 00004 // 00005 // Coded by: Sivapalan Gajan <sgajan@ucdavis.edu> 00006 // 00007 // Description: This file contains the members and methods for 00008 // SoilFootingSection2d class. 00009 00010 00011 #include <math.h> 00012 #include <stdlib.h> 00013 00014 #include <SoilFootingSection2d.h> 00015 #include <Matrix.h> 00016 #include <Vector.h> 00017 #include <ID.h> 00018 #include <FEM_ObjectBroker.h> 00019 #include <MatrixUtil.h> 00020 00021 #include <classTags.h> 00022 00023 ID SoilFootingSection2d::code(3); 00024 00025 00026 00027 // default constructor 00028 00029 SoilFootingSection2d::SoilFootingSection2d(void) 00030 :SectionForceDeformation(0, SEC_TAG_SoilFooting2d), 00031 e(3), s(3),eCommit(3), sCommit(3), deModel(3), ks(3,3), ksE(3,3), ini_size(3) 00032 { 00033 code(0) = SECTION_RESPONSE_P; // vertical load 00034 code(1) = SECTION_RESPONSE_VY; // shear 00035 code(2) = SECTION_RESPONSE_MZ; // moment 00036 } 00037 00038 00039 // constructor 00040 00041 SoilFootingSection2d::SoilFootingSection2d 00042 (int tag, double fs, double vult, double l, double kv, double kh, double rv, double deltaL) 00043 :SectionForceDeformation(tag, SEC_TAG_SoilFooting2d), 00044 e(3), s(3), eCommit(3), sCommit(3), deModel(3), ks(3,3), ksE(3,3), ini_size(3), 00045 FS(fs), Vult(vult), L(l), Kv(kv), Kh(kh), Rv(rv), dL(deltaL) 00046 { 00047 if (FS <= 1.0) 00048 { 00049 opserr <<"SoilFootingSection:: Invalid input for FS\n" 00050 <<"FS should satisfy: FS > 1.0\n"; 00051 } 00052 00053 code(0) = SECTION_RESPONSE_P; // axial load 00054 code(1) = SECTION_RESPONSE_VY; // shear 00055 code(2) = SECTION_RESPONSE_MZ; // moment 00056 00057 V = Vult / FS; 00058 qult = Vult/L; 00059 ecc = 0.0; 00060 dTheta = 0.0; 00061 dThetaPrev = 0.0; 00062 dVtemp = 0.0; 00063 tolerance = 0.0001 * (1.0/FS); 00064 tolerance = 0.01; 00065 Kt = Kv*pow(L, 3.0)/12.0; 00066 00067 incr = 0; 00068 noNodes = (int)(ceil(L/dL))+1; 00069 00070 c1 = 0; 00071 c1T = 0; 00072 c2 = noNodes; 00073 c2T = noNodes; 00074 00075 c1Commit = c1; 00076 c1TCommit = c1T; 00077 c2Commit = c2; 00078 c2TCommit = c2T; 00079 00080 eccCommit = ecc; 00081 00082 hPrev = -10.0; 00083 hCurr = -10.0; 00084 dVtemp = 0.0; 00085 Mmaxpast = 0.0; 00086 Mult = 0.0; 00087 Melastic = 0.0; 00088 00089 00090 initializeBoundingSurface(); 00091 initializeInternalVariables(); 00092 00093 00094 isOver = 0; 00095 isdV = 0; 00096 isElastic = 1; 00097 00098 } 00099 00100 00101 00102 // destructor 00103 00104 00105 SoilFootingSection2d::~SoilFootingSection2d(void) 00106 { 00107 if (foot != 0) 00108 for (int i = 0; i <= noNodes; i++) 00109 delete [] foot[i]; 00110 if (soilMin != 0) 00111 for (int i = 0; i <= noNodes; i++) 00112 delete [] soilMin[i]; 00113 if (soilMax != 0) 00114 for (int i = 0; i <= noNodes; i++) 00115 delete [] soilMax[i]; 00116 if (pressure != 0) 00117 for (int i = 0; i <= noNodes; i++) 00118 delete [] pressure[i]; 00119 if (pressMax != 0) 00120 for (int i = 0; i <= noNodes; i++) 00121 delete [] pressMax[i]; 00122 00123 00124 } 00125 00126 00127 // initializing the bounding surface 00128 00129 void 00130 SoilFootingSection2d::initializeBoundingSurface (void) 00131 { 00132 a = 0.32; 00133 b = 0.37; 00134 ccc = 0.25; 00135 d = 0.55; 00136 eee = 0.8; 00137 f = 0.8; 00138 00139 Fv = V/Vult; 00140 A = a * pow(Fv, ccc) * pow(1-Fv, d); 00141 B = b * pow(Fv, eee) * pow(1-Fv, f); 00142 00143 beta = (A * h) / (pow(A*A*h*h + B*B*L*L, 0.5)); 00144 00145 if (beta < 0.0) 00146 beta *= (-1.0); 00147 00148 alpha = Fv / (1 - beta*(1-Fv)); 00149 pult = alpha; 00150 00151 qult *= alpha; 00152 00153 00154 } 00155 00156 00157 // initializing internal variables 00158 00159 void 00160 SoilFootingSection2d::initializeInternalVariables (void) 00161 { 00162 00163 00164 foot = new double* [noNodes+1]; 00165 soilMin = new double* [noNodes+1]; 00166 soilMax = new double* [noNodes+1]; 00167 pressure = new double* [noNodes+1]; 00168 pressMax = new double* [noNodes+1]; 00169 00170 for (int j = 0; j <= noNodes; j++) 00171 { 00172 foot[j] = new double [ini_size]; 00173 soilMin[j] = new double [ini_size]; 00174 soilMax[j] = new double [ini_size]; 00175 pressure[j] = new double [ini_size]; 00176 pressMax[j] = new double [ini_size]; 00177 } 00178 00179 00180 for (int i = 0; i <= noNodes; i++) 00181 for (int j = 0; j < ini_size; j++) 00182 { 00183 foot[i][j] = V/Kv; 00184 soilMin[i][j] = V/Kv; 00185 soilMax[i][j] = V/Kv; 00186 pressure[i][j] = 1.0/FS; 00187 pressMax[i][j] = 1.0/FS; 00188 } 00189 00190 00191 e.Zero(); 00192 s.Zero(); 00193 eCommit.Zero(); 00194 sCommit.Zero(); 00195 ks.Zero(); 00196 ksE.Zero(); 00197 00198 ks(0,0) = Kv; 00199 ks(1,1) = Kh; 00200 ks(2,2) = Kt; 00201 ksE = ks; 00202 00203 dTh = 0.0; 00204 dThP = 0.0; 00205 00206 Mlimit = V*L/6.0; 00207 thetaPlus = Mlimit / (Kv * pow(L, 2.0) / 12.0); 00208 thetaMinus = -1.0 * Mlimit / (Kv * pow(L, 2.0) / 12.0); 00209 00210 thetaRange = 2 * thetaPlus; 00211 00212 00213 00214 } 00215 00216 00217 void 00218 SoilFootingSection2d::tempFunction (void) 00219 { 00220 00221 } 00222 00223 00224 00225 // updating committed forces and displacements and internal variables 00226 00227 int 00228 SoilFootingSection2d::commitState(void) 00229 { 00230 00231 00232 incr++; 00233 00234 00235 // cout<<"d/din = "<<Mmaxpast<<endl; 00236 00237 00238 if (fabs(s(2)) > Mmaxpast) 00239 Mmaxpast = fabs(s(2)); 00240 00241 if (Mmaxpast > Melastic) 00242 isElastic = 0; 00243 00244 00245 thetaPlusPrev = thetaPlus; 00246 thetaMinusPrev = thetaMinus; 00247 00248 00249 double e_2 = e(2); 00250 00251 00252 if (e(2) > thetaPlus) 00253 { 00254 thetaPlus = e_2; 00255 thetaMinus = thetaPlus - thetaRange; 00256 } 00257 00258 if (e(2) < thetaMinus) 00259 { 00260 thetaMinus = e_2; 00261 thetaPlus = thetaMinus + thetaRange; 00262 } 00263 00264 00265 HPrevCommit = sCommit(1); 00266 MPrevCommit = sCommit(2); 00267 00268 00269 eCommit = e; 00270 sCommit = s; 00271 ksE = ks; 00272 dThetaPrev = dTheta; 00273 00274 c1Commit = c1; 00275 c1TCommit = c1T; 00276 c2Commit = c2; 00277 c2TCommit = c2T; 00278 00279 eccCommit = ecc; 00280 00281 hPrev = hCurr; 00282 00283 for (int i = 0; i <= noNodes; i++) 00284 for (int j = 2; j > 0; j--) 00285 { 00286 foot[i][j] = foot[i][j-1]; 00287 soilMin[i][j] = soilMin[i][j-1]; 00288 soilMax[i][j] = soilMax[i][j-1]; 00289 pressure[i][j] = pressure[i][j-1]; 00290 pressMax[i][j] = pressMax[i][j-1]; 00291 } 00292 00293 tolerance = 0.0000000000001 * (1/FS); 00294 00295 00296 isOver = 1; 00297 isdV = 0; 00298 return 0; 00299 } 00300 00301 00302 // going back to the last committed values 00303 00304 int 00305 SoilFootingSection2d::revertToLastCommit(void) 00306 { 00307 00308 thetaPlus = thetaPlusPrev; 00309 thetaMinus = thetaMinusPrev; 00310 00311 e = eCommit; 00312 s = sCommit; 00313 ks = ksE; 00314 dTheta = dThetaPrev; 00315 00316 c1 = c1Commit; 00317 c1T = c1TCommit; 00318 c2 = c2Commit; 00319 c2T = c2TCommit; 00320 ecc = eccCommit; 00321 00322 hCurr = hPrev; 00323 00324 for (int i = 0; i <= noNodes; i++) 00325 { 00326 foot[i][1] = foot[i][2]; 00327 soilMin[i][1] = soilMin[i][2]; 00328 soilMax[i][1] = soilMax[i][2]; 00329 pressure[i][1] = pressure[i][2]; 00330 pressMax[i][1] = pressMax[i][2]; 00331 } 00332 00333 return 0; 00334 } 00335 00336 00337 // going back to the initial conditions 00338 00339 int 00340 SoilFootingSection2d::revertToStart(void) 00341 { 00342 00343 eCommit.Zero(); 00344 sCommit.Zero(); 00345 00346 c1 = 0; 00347 c1T = 0; 00348 c2 = noNodes; 00349 c2T = noNodes; 00350 00351 c1Commit = c1; 00352 c1TCommit = c1T; 00353 c2Commit = c2; 00354 c2TCommit = c2T; 00355 eccCommit = ecc; 00356 00357 dTheta = 0.0; 00358 dThetaPrev = 0.0; 00359 00360 for (int i = 0; i <= noNodes; i++) 00361 for (int j = 0; j < ini_size; j++) 00362 { 00363 foot[i][j] = V/Kv; 00364 soilMin[i][j] = V/Kv; 00365 soilMax[i][j] = V/Kv; 00366 pressure[i][j] = 1/FS; 00367 pressMax[i][j] = 1/FS; 00368 } 00369 00370 00371 return 0; 00372 } 00373 00374 00375 00376 00377 SectionForceDeformation* SoilFootingSection2d::getCopy() 00378 { 00379 00380 SoilFootingSection2d *theCopy = 00381 new SoilFootingSection2d (this->getTag(), FS, Vult, L, Kv, Kh, Rv, dL); 00382 00383 return theCopy; 00384 } 00385 00386 00387 // most important method in the class !! 00388 // this function is called at the beginning of and at the middle of iterations 00389 00390 int 00391 SoilFootingSection2d::setTrialSectionDeformation (const Vector &def) 00392 { 00393 00394 int temp = 0; 00395 Vector de(3), ds(3); 00396 double epsilon = pow(10.0, -20.0); 00397 00398 00399 00400 e = def; 00401 de = e - eCommit; 00402 00403 00404 00405 if (fabs(de(0)) < epsilon) de(0) = 0.0; 00406 if (fabs(de(1)) < epsilon) de(1) = 0.0; 00407 if (fabs(de(2)) < epsilon) de(2) = 0.0; 00408 00409 00410 00411 deModel.Zero(); 00412 00413 dThP = dTh; 00414 dTh = de(2); 00415 00416 if (de(0) == 0.0 && de(1) == 0.0 && de(2) == 0.0) 00417 { 00418 00419 // give out the previous ks rather than ks_elastic 00420 00421 } 00422 else 00423 temp = applyLoading(de); 00424 00425 00426 ds = ks * deModel; 00427 00428 if (fabs(ds(0)) < epsilon) ds(0) = 0.0; 00429 if (fabs(ds(1)) < epsilon) ds(1) = 0.0; 00430 if (fabs(ds(2)) < epsilon) ds(2) = 0.0; 00431 00432 00433 s = sCommit + ds; 00434 00435 return 0; 00436 } 00437 00438 00439 00440 00441 // applyLoading method 00442 00443 int 00444 SoilFootingSection2d::applyLoading(Vector de) 00445 { 00446 00447 int c, nn = 0, switch1; 00448 double area1, area2, area1Prev = 0.0; 00449 double s_recover = Rv; 00450 double q_recover = 1.0 / L; 00451 double LcOverL; 00452 00453 00454 00455 soilFree = 0.0; 00456 00457 double ds, du, dTheta1, dss, theta; 00458 double dVt, dHt, dMt, dHt1; 00459 double Vinit, dMcal; 00460 double epsilon = pow(10.0, -20.0); 00461 char tempKey; 00462 double detKs; 00463 double expo = 0; 00464 double n_load, n_unload; 00465 double e_2; 00466 00467 if (fabs(de(0)) < epsilon) de(0) = 0.0; 00468 if (fabs(de(1)) < epsilon) de(1) = 0.0; 00469 if (fabs(de(2)) < epsilon) de(2) = 0.0; 00470 00471 00472 du = de(1); 00473 dTheta = de(2); 00474 dTheta1 = de(2); 00475 theta = eCommit(2); 00476 dss = 0.0; 00477 00478 00479 c1 = c1Commit; 00480 c1T = c1TCommit; 00481 c2 = c2Commit; 00482 c2T = c2TCommit; 00483 ecc = eccCommit; 00484 00485 00486 00487 00488 double *footTemp = new double[noNodes+1]; 00489 double *soilMinTemp = new double[noNodes+1]; 00490 double *soilMaxTemp = new double[noNodes+1]; 00491 double *pressureTemp = new double[noNodes+1]; 00492 double *pressMaxTemp = new double[noNodes+1]; 00493 00494 double *ddH = new double[200]; 00495 00496 // for shear sliding model 00497 00498 double deltaIn, delta; 00499 double a11, b11, c11, gr; 00500 double Fh1, Fh2, Fm1, Fm2; 00501 double ptFh1, ptFh2, ptFm1, ptFm2; 00502 double dist1, dist2, factor; 00503 double hNew, hNew1, uH, uM, dVt1; 00504 int ii; 00505 00506 00507 Vector tempLoad(3); 00508 00509 // extracting the forces - 00510 // should be close enough to the actual applied external loads 00511 00512 tempLoad = ks * de; 00513 00514 dVt = tempLoad(0); 00515 dHt = tempLoad(1); 00516 dMt = tempLoad(2); 00517 00518 00519 if (fabs(dVt) < epsilon) dVt = 0.0; 00520 if (fabs(dHt) < epsilon) dHt = 0.0; 00521 if (fabs(dMt) < epsilon) dMt = 0.0; 00522 00523 00524 // set-up the state of switch depending on external loading 00525 00526 switch1 = -1; 00527 00528 if (dVt == 0.0) 00529 if (dHt == 0.0) 00530 switch1 = (dMt == 0) ? 0 : 1; 00531 else 00532 switch1 = (dMt == 0) ? 2 : 3; 00533 else 00534 if (dHt == 0.0) 00535 switch1 = (dMt == 0) ? 4 : 5; 00536 else 00537 switch1 = (dMt == 0) ? 6 : 7; 00538 00539 00540 00541 if (isOver == 0) 00542 switch1 = 7; 00543 00544 00545 ds = 0.0; 00546 deModel.Zero(); 00547 ks.Zero(); 00548 00549 00550 for (int i = 0; i <= noNodes; i++) 00551 { 00552 footTemp[i] = foot[i][1]; 00553 soilMinTemp[i] = soilMin[i][1]; 00554 soilMaxTemp[i] = soilMax[i][1]; 00555 pressureTemp[i] = pressure[i][1]; 00556 pressMaxTemp[i] = pressMax[i][1]; 00557 } 00558 00559 for (int j = 0; j < 190; j++) 00560 ddH[j] = 0.0; 00561 00562 00563 Vinit = sCommit(0); 00564 00565 00566 00567 00568 if (switch1 == 0) 00569 { 00570 // ALL ZERO ! 00571 00572 ks(0,0) = epsilon; 00573 ks(1,1) = epsilon; 00574 ks(2,2) = epsilon; 00575 00576 ks(0,0) = Kv; 00577 ks(1,1) = Kh; 00578 ks(2,2) = Kt; 00579 00580 00581 ks(0,0) = 1.0; 00582 ks(1,1) = 1.0; 00583 ks(2,2) = 1.0; 00584 00585 00586 deModel(0) = de(0); 00587 deModel(1) = de(1); 00588 deModel(2) = de(2); 00589 00590 return 0; 00591 } 00592 00593 else if (switch1 == 4) 00594 { 00595 // ONLY dVt 00596 00597 deModel(0) = dVt / Kv; 00598 00599 ks(0,0) = Kv; 00600 ks(1,1) = Kh; 00601 ks(2,2) = Kt; 00602 00603 return 0; 00604 } 00605 00606 else 00607 { 00608 // APPLY dVt in any case 00609 00610 00611 if (dVt == 0.0) 00612 if (dHt == 0.0) 00613 switch1 = (dMt == 0) ? 0 : 1; 00614 else 00615 switch1 = (dMt == 0) ? 2 : 3; 00616 00617 if (isOver == 0) 00618 switch1 = 7; 00619 00620 00621 00622 if (isOver != 0) 00623 isOver++; 00624 00625 00626 dVt = 0.0; 00627 00628 // do { 00629 00630 dVt1 = dVt; 00631 Vinit = sCommit(0)+dVt1; 00632 00633 00634 if (sCommit(1) != 0.0) 00635 hCurr = sCommit(2)/sCommit(1); 00636 // else if (dHt != 0.0) 00637 // hCurr = dMt / dHt; 00638 else 00639 hCurr = hPrev; 00640 00641 00642 00643 00644 if (hCurr < 0.01 && hCurr >= 0.0) 00645 hCurr = 0.01; 00646 if (hCurr > -0.01 && hCurr <= 0.0) 00647 hCurr = -0.01; 00648 00649 00650 00651 FS = Vult/Vinit; 00652 Fv = Vinit/Vult; 00653 00654 00655 LcOverL = 1.0 / (FS * alpha); 00656 s_recover = Rv * (1.0 - LcOverL); 00657 00658 00659 00660 00661 A = a * pow(Fv, ccc) * pow(1-Fv, d); 00662 B = b * pow(Fv, eee) * pow(1-Fv, f); 00663 00664 beta = (A * hCurr) / (pow(A*A*hCurr*hCurr + B*B*L*L, 0.5)); 00665 00666 if (beta < 0.0) 00667 beta *= (-1.0); 00668 00669 // if no shear 00670 if (switch1 == 1 || switch1 == 5) 00671 beta = 1.0; 00672 00673 alpha = Fv / (1 - beta*(1-Fv)); 00674 pult = alpha; 00675 00676 Mult = (Vinit*L/2.0) * beta * (1.0 - Fv); 00677 00678 Melastic = Mult / (3.0 * (1.0 - Fv)); 00679 00680 00681 // apply M-Theta model 00682 00683 00684 00685 if (dTheta > 0.0) 00686 { 00687 c = c1; 00688 00689 00690 if (c > noNodes) 00691 c = noNodes; 00692 00693 00694 do // check for Vinit 00695 { 00696 int i; 00697 for (i = 0; i <= noNodes; i++) 00698 { 00699 ds = (i-c)*(L/noNodes)*tan(dTheta); 00700 00701 foot[i][0] = footTemp[i] + ds; 00702 00703 if (foot[i][0] >= soilMaxTemp[i]) 00704 soilMax[i][0] = foot[i][0]; 00705 else 00706 soilMax[i][0] = soilMaxTemp[i]; 00707 00708 if (foot[i][0] >= soilMinTemp[i]) 00709 soilMin[i][0] = foot[i][0]; 00710 else 00711 soilMin[i][0] = soilMinTemp[i]; 00712 } 00713 00714 c1 = c2 = -1; 00715 00716 for (i = 0; i <= noNodes; i++) 00717 if (foot[i][0] >= soilMax[i][0]) 00718 { 00719 c1 = i; 00720 break; 00721 } 00722 00723 for (i = noNodes; i >= 0; i--) 00724 if (foot[i][0] >= soilMax[i][0]) 00725 { 00726 c2 = i; 00727 break; 00728 } 00729 00730 00731 00732 if ((c1 == -1) && (c2 == -1)) 00733 { 00734 00735 if ((c >= 0) && (c <= noNodes)) 00736 c1 = c2 = c; 00737 else 00738 c1 = c2 = c1T; 00739 } 00740 00741 00742 for (i = 0; i <= c1; i++) 00743 { 00744 soilMin[i][0] = soilMax[i][0] 00745 - (soilMax[i][0] - soilFree)*s_recover; 00746 if (foot[i][0] >= soilMin[i][0]) 00747 soilMin[i][0] = foot[i][0]; 00748 } 00749 00750 for (i = 0; i <= noNodes; i++) 00751 if (foot[i][0] >= soilMin[i][0]) 00752 { 00753 c1T = i; 00754 break; 00755 } 00756 00757 for (i = noNodes; i >= 0; i--) 00758 if (foot[i][0] >= soilMin[i][0]) 00759 { 00760 c2T = i; 00761 break; 00762 } 00763 00764 if (c1 == 0) 00765 c1T = 0; 00766 if (c2 == noNodes) 00767 c2T = noNodes; 00768 00769 if (c1T != 0) 00770 for (int i = 0; i <= c1T; i++) 00771 pressure[i][0] = 0.0; 00772 00773 for (i = c1; i <= c2; i++) 00774 pressure[i][0] = pressureTemp[i] + 00775 (soilMax[i][0]-soilMin[i][1])*q_recover* Kv/ Vult; 00776 00777 00778 expo = ((10.0 - 0.5)/(0.66*noNodes))*c1T + 0.5; 00779 00780 n_load = 0.5; 00781 n_unload = 10.0; 00782 00783 e_2 = eCommit(2); 00784 n_load = 9.5 * pow((thetaPlus-e_2)/thetaRange, 1.0) + 0.5; 00785 n_unload = 9.5 * pow((e_2-thetaMinus)/thetaRange, 1.0) + 0.5; 00786 00787 00788 00789 00790 if (c1 != c1T) 00791 for (int i = c1T; i <= c1; i++) 00792 pressure[i][0] = pow((double)(i-c1T)/(c1-c1T), n_unload) 00793 * pressure[c1][0]; 00794 00795 00796 if (c2 != c2T) 00797 for (int i = c2; i <= c2T; i++) 00798 pressure[i][0] = pow((double)(c2T-i)/(c2T-c2), n_load) 00799 * pressure[c2][0]; 00800 00801 for (i = c2T+1; i <= noNodes; i++) 00802 pressure[i][0] = 0.0; 00803 00804 int a; 00805 00806 for (a = 0; a <= noNodes; a++) 00807 { 00808 if (pressure[a][0] > pult) 00809 pressure[a][0] = pult; 00810 if (pressure[a][0] < 0.0) 00811 pressure[a][0] = 0.0; 00812 } 00813 00814 area1 = 0.0; 00815 for (a = 0; a <= noNodes; a++) 00816 area1 += pressure[a][0]; 00817 area1 *= (L/noNodes); 00818 00819 area2 = 0.0; 00820 for (a = 0; a <= noNodes; a++) 00821 area2 += (pressure[a][0] * (L/noNodes) 00822 * (a - noNodes/2)*(L/noNodes)); 00823 00824 for (a = 0; a <= noNodes; a++) 00825 if (pressure[a][0] > pressMax[a][1]) 00826 pressMax[a][0] = pressure[a][0]; 00827 else 00828 pressMax[a][0] = pressMax[a][1]; 00829 00830 ecc = area2/area1; 00831 area1 *= cos(theta); 00832 00833 if ((area1/L) > (Vinit/Vult)) 00834 c += 1; 00835 else 00836 nn = -2; 00837 00838 00839 if (c >= noNodes) 00840 nn = -2; 00841 if (c <= 0) 00842 nn = -2; 00843 00844 00845 nn++; 00846 if ((area1/L-Vinit/Vult < 0.0) && (area1Prev/L-Vinit/Vult > 0.0)) 00847 { 00848 tempKey = 'y'; 00849 if (tempKey == 'y' || tempKey == 'Y') 00850 { 00851 00852 tolerance *= 2.0; 00853 } 00854 else 00855 { 00856 } 00857 } 00858 00859 area1Prev = area1; 00860 00861 } while ((fabs(area1/L - Vinit/Vult) > tolerance) && (nn != -1)); 00862 00863 } 00864 00865 else if (dTheta < 0.0) 00866 { 00867 00868 00869 c = c2; 00870 00871 00872 if (c < 0) 00873 c = 0; 00874 00875 00876 do // check for Vinit 00877 { int i; 00878 for (i = 0; i <= noNodes; i++) 00879 { 00880 ds = (i-c)*(L/noNodes)*tan(dTheta); 00881 00882 00883 foot[i][0] = footTemp[i] + ds; 00884 00885 if (foot[i][0] >= soilMaxTemp[i]) 00886 soilMax[i][0] = foot[i][0]; 00887 else 00888 soilMax[i][0] = soilMaxTemp[i]; 00889 00890 if (foot[i][0] >= soilMinTemp[i]) 00891 soilMin[i][0] = foot[i][0]; 00892 else 00893 soilMin[i][0] = soilMinTemp[i]; 00894 } 00895 00896 c1 = c2 = -1; 00897 00898 for (i = 0; i <= noNodes; i++) 00899 if (foot[i][0] >= soilMax[i][0]) 00900 { 00901 c1 = i; 00902 break; 00903 } 00904 00905 for (i = noNodes; i >= 0; i--) 00906 if (foot[i][0] >= soilMax[i][0]) 00907 { 00908 c2 = i; 00909 break; 00910 } 00911 00912 if ((c1 == -1) && (c2 == -1)) 00913 { 00914 00915 if ((c >= 0) && (c <= noNodes)) 00916 c1 = c2 = c; 00917 else 00918 c1 = c2 = c2T; 00919 } 00920 00921 00922 for (i = c2; i <= noNodes; i++) 00923 { 00924 soilMin[i][0] = soilMax[i][0] 00925 - (soilMax[i][0] - soilFree)*s_recover; 00926 if (foot[i][0] >= soilMin[i][0]) 00927 soilMin[i][0] = foot[i][0]; 00928 } 00929 00930 for (i = 0; i <= noNodes; i++) 00931 if (foot[i][0] >= soilMin[i][0]) 00932 { 00933 c1T = i; 00934 break; 00935 } 00936 00937 for (i = noNodes; i >= 0; i--) 00938 if (foot[i][0] >= soilMin[i][0]) 00939 { 00940 c2T = i; 00941 break; 00942 } 00943 00944 if (c1 <= 0) 00945 c1T = 0; 00946 if (c2 >= noNodes) 00947 c2T = noNodes; 00948 00949 00950 if (c1T != 0) 00951 for (int i = 0; i <= c1T; i++) 00952 pressure[i][0] = 0.0; 00953 00954 for (i = c1; i <= c2; i++) 00955 pressure[i][0] = pressureTemp[i] + 00956 (soilMax[i][0]-soilMin[i][1])*q_recover* Kv/ Vult; 00957 00958 00959 expo = ((10.0 - 0.5)/(0.66*noNodes))*(noNodes - c2T) + 0.5; 00960 00961 00962 n_load = 0.5; 00963 n_unload = 10.0; 00964 00965 e_2 = eCommit(2); 00966 n_unload = 9.5 * pow((thetaPlus-e_2)/thetaRange, 1.0) + 0.5; 00967 n_load = 9.5 * pow((e_2-thetaMinus)/thetaRange, 1.0) + 0.5; 00968 00969 00970 00971 if (c2 != c2T) 00972 for (int i = c2; i <= c2T; i++) 00973 pressure[i][0] = pow((double)(c2T-i)/(c2T-c2), n_unload) 00974 * pressure[c2][0]; 00975 00976 if (c1 != c1T) 00977 for (int i = c1T; i <= c1; i++) 00978 pressure[i][0] = pow((double)(i-c1T)/(c1-c1T), n_load) 00979 * pressure[c1][0]; 00980 00981 for (i = c2T+1; i <= noNodes; i++) 00982 pressure[i][0] = 0.0; 00983 int a; 00984 for (a = 0; a <= noNodes; a++) 00985 { 00986 if (pressure[a][0] > pult) 00987 pressure[a][0] = pult; 00988 if (pressure[a][0] < 0.0) 00989 pressure[a][0] = 0.0; 00990 } 00991 00992 area1 = 0.0; 00993 for (a = 0; a <= noNodes; a++) 00994 area1 += pressure[a][0]; 00995 area1 *= (L/noNodes); 00996 00997 area2 = 0.0; 00998 for (a = 0; a <= noNodes; a++) 00999 area2 += (pressure[a][0] * (L/noNodes) 01000 * (a - noNodes/2)*(L/noNodes)); 01001 01002 01003 for (a = 0; a <= noNodes; a++) 01004 if (pressure[a][0] > pressMax[a][1]) 01005 pressMax[a][0] = pressure[a][0]; 01006 else 01007 pressMax[a][0] = pressMax[a][1]; 01008 01009 ecc = area2/area1; 01010 area1 *= cos(theta); 01011 01012 if ((area1/L) > (Vinit/Vult)) 01013 c -= 1; 01014 else 01015 nn = -2; 01016 01017 01018 if (c <= 0) 01019 nn = -2; 01020 if (c >= noNodes) 01021 nn = -2; 01022 01023 01024 01025 01026 nn++; 01027 if ((area1/L-Vinit/Vult < 0.0) && (area1Prev/L-Vinit/Vult > 0.0)) 01028 { 01029 tempKey = 'y'; 01030 01031 if (tempKey == 'y' || tempKey == 'Y') 01032 { 01033 01034 tolerance *= 2.0; 01035 } 01036 else 01037 { 01038 } 01039 } 01040 01041 area1Prev = area1; 01042 01043 01044 } while ((fabs(area1/L - Vinit/Vult) > tolerance) && (nn != -1)); 01045 01046 } 01047 01048 else 01049 { 01050 dMcal = 0.0; 01051 } 01052 01053 // Update Moment 01054 01055 M = Vinit * ecc; 01056 dMcal = M - sCommit(2); 01057 01058 if (dTheta == 0.0) 01059 dMcal = 0.0; 01060 01061 01062 dMt = dMcal; 01063 ds = foot[noNodes/2][0] - foot[noNodes/2][1]; 01064 01065 01066 /* 01067 01068 if (((dTheta >= 0.0) && (eCommit(2)+dTheta <= thetaPlus)) || 01069 ((dTheta < 0.0) && (eCommit(2)+dTheta >= thetaMinus))) 01070 ds = 0.0; 01071 01072 */ 01073 01074 01075 if (isOver == 0) 01076 ds = 0.0; 01077 01078 dVt = Kv * (de(0) - ds); 01079 // dVt = 1.0 * (de(0) - ds); 01080 01081 01082 if (sCommit(0)+dVt > Vult) 01083 dVt = 0.99*Vult - sCommit(0); 01084 if (sCommit(0)+dVt < 0.0) 01085 dVt = -1.0*sCommit(0) + 0.01*Vult; 01086 01087 01088 // } while ((fabs(dVt - dVt1) > 1.0) && (isOver != 0)); 01089 01090 01091 01093 01094 Fm2 = (sCommit(2)+dMt)/Vult/L; 01095 01096 if (Fm2 > B) 01097 { 01098 dMt = 0.0; 01099 // exit(0); 01100 } 01101 01102 01103 for (int i = 0; i <= 180; i++) 01104 { 01105 01106 01107 Fm2 = (sCommit(2)+dMt)/Vult/L; 01108 01109 01110 if (A*A*B*B - Fm2*Fm2*A*A > 0.0) 01111 { 01112 hNew = fabs(pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)*Vult); 01113 dHt = -1.0 * hNew + (2.0*hNew/180.0)*i - sCommit(1); 01114 } 01115 01116 else 01117 { 01118 /* 01119 if (isOver != 0) 01120 { 01121 hNew = sCommit(1); 01122 cout <<"TOP\n"; 01123 cout <<"A, B, Fm2, dMt = "<<A<<" "<<B<<" "<<Fm2<<" "<<dMt<<endl; 01124 cout <<"A*A*B*B - Fm2*Fm2*A*A = "<<A*A*B*B - Fm2*Fm2*A*A<<endl; 01125 exit(0); 01126 } 01127 */ 01128 } 01129 01130 // dHt = -1.0 * hNew + (2.0*hNew/180.0)*i - sCommit(1); 01131 01132 Fh1 = sCommit(1)/Vult; 01133 Fm1 = sCommit(2)/Vult/L; 01134 Fh2 = (sCommit(1)+dHt)/Vult; 01135 Fm2 = (sCommit(2)+dMt)/Vult/L; 01136 01137 01138 if (Fh1 != Fh2) 01139 { 01140 gr = (Fm2 - Fm1)/(Fh2 - Fh1); 01141 01142 a11 = B*B + A*A*gr*gr; 01143 b11 = 2.0*gr*A*A * (Fm1 - gr*Fh1); 01144 c11 = A*A * (Fh1*Fh1*gr*gr + Fm1*Fm1 - 2.0*gr*Fh1*Fm1 - B*B); 01145 01146 factor = b11*b11 - 4.0*a11*c11; 01147 01148 if (factor < 0.0) 01149 factor = 0.0; 01150 01151 01152 if (factor >= 0.0) 01153 { 01154 ptFh1 = (-b11 + pow(factor, 0.5)) / 2.0 / a11; 01155 ptFh2 = (-b11 - pow(factor, 0.5)) / 2.0 / a11; 01156 ptFm1 = gr*(ptFh1 - Fh1) + Fm1; 01157 ptFm2 = gr*(ptFh2 - Fh1) + Fm1; 01158 } 01159 else 01160 { 01161 01162 } 01163 } 01164 else 01165 { 01166 01167 ptFh1 = fabs(Fh1); 01168 ptFh2 = -1.0 * fabs(Fh1); 01169 if ((A*A*B*B-B*B*Fh1*Fh1)/(A*A) < 0.0) 01170 ptFm1 = ptFm2 = 0.0; 01171 else 01172 { 01173 ptFm1 = pow((A*A*B*B-B*B*Fh1*Fh1)/(A*A), 0.5); 01174 ptFm2 = -1.0 * pow((A*A*B*B-B*B*Fh1*Fh1)/(A*A), 0.5); 01175 } 01176 } 01177 01178 01179 if (fabs(ptFh1) > 0.25 || fabs(ptFh2) > 0.25 || 01180 fabs(ptFm1) > 0.15 || fabs(ptFm2) > 0.15 ) 01181 { 01182 } 01183 01184 01186 // One way to make sure that d <= din is to make Fh1 <= ptFh1 .... and so on 01188 01189 if (ptFh1 > ptFh2) 01190 { 01191 if (Fh1 > ptFh1) 01192 Fh1 = ptFh1; 01193 if (Fh2 > ptFh1) 01194 Fh2 = ptFh1; 01195 if (Fh1 < ptFh2) 01196 Fh1 = ptFh2; 01197 if (Fh2 < ptFh2) 01198 Fh2 = ptFh2; 01199 } 01200 else 01201 { 01202 if (Fh1 < ptFh1) 01203 Fh1 = ptFh1; 01204 if (Fh2 < ptFh1) 01205 Fh2 = ptFh1; 01206 if (Fh1 > ptFh2) 01207 Fh1 = ptFh2; 01208 if (Fh2 > ptFh2) 01209 Fh2 = ptFh2; 01210 } 01211 01212 if (ptFm1 > ptFm2) 01213 { 01214 if (Fm1 > ptFm1) 01215 Fm1 = ptFm1; 01216 if (Fm2 > ptFm1) 01217 Fm2 = ptFm1; 01218 if (Fm1 < ptFm2) 01219 Fm1 = ptFm2; 01220 if (Fm2 < ptFm2) 01221 Fm2 = ptFm2; 01222 } 01223 else 01224 { 01225 if (Fm1 < ptFm1) 01226 Fm1 = ptFm1; 01227 if (Fm2 < ptFm1) 01228 Fm2 = ptFm1; 01229 if (Fm1 > ptFm2) 01230 Fm1 = ptFm2; 01231 if (Fm2 > ptFm2) 01232 Fm2 = ptFm2; 01233 } 01234 01235 01236 deltaIn = pow(pow(ptFm2-ptFm1, 2.0) + pow(ptFh2-ptFh1, 2.0), 0.5); 01237 01238 deltaIn = fabs(2*A); 01239 01240 dist1 = pow(pow(ptFm1-Fm1, 2.0) + pow(ptFh1-Fh1, 2.0), 0.5); 01241 dist2 = pow(pow(ptFm1-Fm2, 2.0) + pow(ptFh1-Fh2, 2.0), 0.5); 01242 01243 01244 if (dist2 < dist1) 01245 delta = dist2; 01246 else 01247 delta = pow(pow(ptFm2-Fm2, 2.0) + pow(ptFh2-Fh2, 2.0), 0.5); 01248 01249 if (dist1 == dist2) 01250 delta = 0.0; 01251 01252 if (deltaIn - delta <= 0.0) 01253 { 01254 // exit (0); 01255 01256 if (deltaIn != delta) 01257 { 01258 } 01259 01260 deltaIn = delta + 0.0001; 01261 } 01262 01263 /* 01264 if (dist2 < dist1) 01265 hNew = ptFm1*L/ptFh1; 01266 else 01267 hNew = ptFm2*L/ptFh2; 01268 */ 01269 01270 01271 if (delta/deltaIn > 0.05) 01272 { 01273 uM = (2.0 * L * (dTheta) * (L/hCurr)*pow(B/A, 2.0))/((delta/(deltaIn-delta))/0.05); 01274 uM = (2.0 * L * (dTheta) * (L/hCurr)*pow(B/A, 2.0)) * pow(1.0 - delta/deltaIn, 2.0); 01275 // uM = 0.0; 01276 uH = de(1) - uM; 01277 dHt1 = Kh * uH * pow(delta/(deltaIn), 1.0); 01278 dHt1 = Kh * uH; 01279 } 01280 else 01281 { 01282 uM = (2.0 * L * (dTheta) * (L/hCurr)*pow(B/A, 2.0)); 01283 uM = (2.0 * L * (dTheta) * (L/hCurr)*pow(B/A, 2.0)) * pow(1.0 - delta/deltaIn, 2.0); 01284 // uM = 0.0; 01285 uH = de(1) - uM; 01286 dHt1 = Kh * uH * pow(delta/(deltaIn), 1.0); 01287 dHt1 = Kh * uH; 01288 } 01289 01290 01291 /* 01292 if (delta/deltaIn > 0.02) 01293 { 01294 uH = dHt / (Kh * pow(delta/(deltaIn-delta), 1.0)); 01295 ddH[i] = fabs(de(1) - uH); 01296 } 01297 else 01298 { 01299 uM = (2.0 * L * (dTheta) * (L/hCurr)*pow(B/A, 2.0)); 01300 ddH[i] = fabs(de(1) - uM); 01301 } 01302 01303 */ 01304 01305 01306 ddH[i] = fabs(dHt - dHt1); 01307 01308 01309 01310 } // for loop for finding dHt ends here 01311 01312 01313 dHt1 = 100000000000.0; 01314 01315 for (int j = 0; j <= 180; j++) 01316 { 01317 if (ddH[j] <= dHt1) 01318 { 01319 dHt1 = ddH[j]; 01320 ii = j; 01321 } 01322 } 01323 01324 01325 01326 01328 01329 01330 if (A*A*B*B - Fm2*Fm2*A*A > 0.0) 01331 { 01332 hNew = fabs(pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)*Vult); 01333 dHt = -1.0 * hNew + (2.0*hNew/180.0)*ii - sCommit(1); 01334 } 01335 01336 else 01337 { 01338 /* 01339 if (isOver != 0) 01340 { 01341 hNew = sCommit(1); 01342 cout <<"BOT\n"; 01343 exit(0); 01344 } 01345 */ 01346 } 01347 01348 // dHt = -1.0 * hNew + (2.0*hNew/180.0)*ii - sCommit(1); 01349 01350 if (sCommit(1) != 0.0) 01351 hCurr = sCommit(2)/sCommit(1); 01352 // else if (dHt != 0.0) 01353 // hCurr = dMt / dHt; 01354 else 01355 hCurr = hPrev; 01356 01357 01358 01359 if (hCurr < 0.01 && hCurr >= 0.0) 01360 hCurr = 0.01; 01361 if (hCurr > -0.01 && hCurr <= 0.0) 01362 hCurr = -0.01; 01363 01364 01365 01366 hNew = 0.0; 01367 01368 int gate = 1; 01369 dHt = dMt / ((hCurr+hPrev)/2.0); 01370 dHt1 = dMt / ((hCurr+hPrev)/2.0); 01371 01372 01373 do { 01374 01375 // dHt = (dHt + dHt1)/2.0; 01376 // dHt = dHt/2.0; 01377 // dHt = dHt1; 01378 01379 01380 hNew = dHt - dHt1; 01381 01382 /* 01383 if (dHt1 > dHt) 01384 dHt += 10.0; 01385 else if (dHt1 < dHt) 01386 dHt -= 10.0; 01387 */ 01388 01389 dHt = (dHt + dHt1)/2.0; 01390 01391 dHt = dMt / hCurr; 01392 01393 01394 // cout <<"dHt assumed = "<<dHt<<endl; 01395 01396 01397 Fh1 = sCommit(1)/Vult; 01398 Fm1 = sCommit(2)/Vult/L; 01399 Fh2 = (sCommit(1)+dHt)/Vult; 01400 Fm2 = (sCommit(2)+dMt)/Vult/L; 01401 01402 01403 dHt1 = fabs(pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)); 01404 if (Fh2 > dHt1) 01405 dHt = dHt1*Vult - sCommit(1); 01406 if (Fh2 < -1.0*dHt1) 01407 dHt = -1.0*dHt1*Vult - sCommit(1); 01408 01409 01410 01411 01412 if (Fh2*Fh2/A/A + Fm2*Fm2/B/B - 1.0 > 0.0) 01413 { 01414 01415 /* 01416 // if (dHt > 0.0) 01417 if (Fh2 > 0.0) 01418 dHt = fabs(pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)*Vult)-sCommit(1); 01419 // if (dHt < 0.0) 01420 else if (Fh2 < 0.0) 01421 dHt = -1.0 * fabs(pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)*Vult)-sCommit(1); 01422 01423 01424 01425 if (dMt > 0.0) 01426 dMt = pow(pow(A*B*hCurr, 2.0)/(A*A*hCurr*hCurr+B*B*L*L), 0.5)*Vult*L - sCommit(2); 01427 else if (dMt < 0.0) 01428 dMt = -1.0*pow(pow(A*B*hCurr, 2.0)/(A*A*hCurr*hCurr+B*B*L*L), 0.5)*Vult*L - sCommit(2); 01429 dHt = dMt/hCurr; 01430 01431 01432 cout<<"In top\n"; 01433 cout <<"sCommit(1) and (2) = "<<sCommit(1)<<" "<<sCommit(2)<<endl; 01434 cout <<"dHt and dMt = "<<dHt<<" "<<dMt<<endl; 01435 cout <<"sCommit(0) and dVt = "<<sCommit(0)<<" "<<dVt<<endl; 01436 cout <<"Fh2 A Fm2 B = "<<Fh2<<" "<<A<<" "<<Fm2<<" "<<B<<endl; 01437 exit(0); 01438 */ 01439 01440 } 01441 01442 01443 Fh2 = (sCommit(1)+dHt)/Vult; 01444 Fm2 = (sCommit(2)+dMt)/Vult/L; 01445 01446 01447 01448 if (Fh1 != Fh2) 01449 { 01450 gr = (Fm2 - Fm1)/(Fh2 - Fh1); 01451 01452 a11 = B*B + A*A*gr*gr; 01453 b11 = 2.0*gr*A*A * (Fm1 - gr*Fh1); 01454 c11 = A*A * (Fh1*Fh1*gr*gr + Fm1*Fm1 - 2.0*gr*Fh1*Fm1 - B*B); 01455 01456 factor = b11*b11 - 4.0*a11*c11; 01457 01458 if (factor < 0.0) 01459 factor = 0.0; 01460 01461 01462 if (factor >= 0.0) 01463 { 01464 ptFh1 = (-b11 + pow(factor, 0.5)) / 2.0 / a11; 01465 ptFh2 = (-b11 - pow(factor, 0.5)) / 2.0 / a11; 01466 ptFm1 = gr*(ptFh1 - Fh1) + Fm1; 01467 ptFm2 = gr*(ptFh2 - Fh1) + Fm1; 01468 } 01469 else 01470 { 01471 01472 } 01473 } 01474 else 01475 { 01476 01477 ptFh1 = fabs(Fh1); 01478 ptFh2 = -1.0 * fabs(Fh1); 01479 if ((A*A*B*B-B*B*Fh1*Fh1)/(A*A) < 0.0) 01480 ptFm1 = ptFm2 = 0.0; 01481 else 01482 { 01483 ptFm1 = pow((A*A*B*B-B*B*Fh1*Fh1)/(A*A), 0.5); 01484 ptFm2 = -1.0 * pow((A*A*B*B-B*B*Fh1*Fh1)/(A*A), 0.5); 01485 } 01486 } 01487 01488 01489 if (fabs(ptFh1) > 0.25 || fabs(ptFh2) > 0.25 || 01490 fabs(ptFm1) > 0.15 || fabs(ptFm2) > 0.15 ) 01491 { 01492 } 01493 01494 01496 // One way to make sure that d <= din is to make Fh1 <= ptFh1 .... and so on 01498 01499 if (ptFh1 > ptFh2) 01500 { 01501 if (Fh1 > ptFh1) 01502 Fh1 = ptFh1; 01503 if (Fh2 > ptFh1) 01504 Fh2 = ptFh1; 01505 if (Fh1 < ptFh2) 01506 Fh1 = ptFh2; 01507 if (Fh2 < ptFh2) 01508 Fh2 = ptFh2; 01509 } 01510 else 01511 { 01512 if (Fh1 < ptFh1) 01513 Fh1 = ptFh1; 01514 if (Fh2 < ptFh1) 01515 Fh2 = ptFh1; 01516 if (Fh1 > ptFh2) 01517 Fh1 = ptFh2; 01518 if (Fh2 > ptFh2) 01519 Fh2 = ptFh2; 01520 } 01521 01522 if (ptFm1 > ptFm2) 01523 { 01524 if (Fm1 > ptFm1) 01525 Fm1 = ptFm1; 01526 if (Fm2 > ptFm1) 01527 Fm2 = ptFm1; 01528 if (Fm1 < ptFm2) 01529 Fm1 = ptFm2; 01530 if (Fm2 < ptFm2) 01531 Fm2 = ptFm2; 01532 } 01533 else 01534 { 01535 if (Fm1 < ptFm1) 01536 Fm1 = ptFm1; 01537 if (Fm2 < ptFm1) 01538 Fm2 = ptFm1; 01539 if (Fm1 > ptFm2) 01540 Fm1 = ptFm2; 01541 if (Fm2 > ptFm2) 01542 Fm2 = ptFm2; 01543 } 01544 01545 01546 01547 deltaIn = pow(pow(ptFm2-ptFm1, 2.0) + pow(ptFh2-ptFh1, 2.0), 0.5); 01548 deltaIn = fabs(2.0*A); 01549 // deltaIn = 0.184; 01550 01551 01552 dist1 = pow(pow(ptFm1-Fm1, 2.0) + pow(ptFh1-Fh1, 2.0), 0.5); 01553 dist2 = pow(pow(ptFm1-Fm2, 2.0) + pow(ptFh1-Fh2, 2.0), 0.5); 01554 01555 01556 if (dist2 < dist1) 01557 { 01558 // delta = dist2; 01559 // delta = (dist1+dist2)/2.0; 01560 delta = dist1; 01561 01562 delta = pow(pow(ptFh1-Fh1, 2.0), 0.5); 01563 01564 01565 } 01566 else 01567 { 01568 // delta = pow(pow(ptFm2-Fm2, 2.0) + pow(ptFh2-Fh2, 2.0), 0.5); 01569 dist1 = pow(pow(ptFm2-Fm1, 2.0) + pow(ptFh2-Fh1, 2.0), 0.5); 01570 dist2 = pow(pow(ptFm2-Fm2, 2.0) + pow(ptFh2-Fh2, 2.0), 0.5); 01571 // delta = (dist1+dist2)/2.0; 01572 delta = dist1; 01573 01574 delta = pow(pow(ptFh2-Fh1, 2.0), 0.5); 01575 01576 01577 } 01578 01579 // if (dist1 == dist2) 01580 // delta = 0.0; 01581 // delta = deltaIn; 01582 01583 01584 01585 /* 01586 01587 Mult = Vult*L*(Fv/2.0)*(hCurr/L)*(1.0 - Fv) / pow(pow(B/A, 2.0) + pow(hCurr/L, 2.0), 0.5); 01588 deltaIn = fabs(2.0 * Mult); 01589 01590 if (dMt >= 0.0) 01591 delta = fabs(-1.0*Mult - (sCommit(2)+dMt)); 01592 else 01593 delta = fabs(Mult - (sCommit(2)+dMt)); 01594 01595 */ 01596 01597 01598 if (deltaIn - delta <= 0.0) 01599 { 01600 01601 // exit(0); 01602 01603 if (deltaIn != delta) 01604 { 01605 } 01606 01607 deltaIn = delta + 0.0001; 01608 } 01609 01610 /* 01611 if (incr > 910) 01612 { 01613 cout <<"incr = "<<incr<<" "<<delta<<" "<<deltaIn<<endl; 01614 // cout <<"delta = "<<delta<<endl; 01615 // cout <<"deltaIn = "<<deltaIn<<endl; 01616 } 01617 */ 01618 01619 01620 /* 01621 01622 if (delta/deltaIn > 0.05) 01623 { 01624 uM = (2.0 * L * (dTheta) * (L/hCurr)*pow(B/A, 2.0))/((delta/(deltaIn-delta))/0.05); 01625 uM = (2.0 * L * (dTheta) * (L/hCurr)*pow(B/A, 2.0)) * pow(1.0 - delta/deltaIn, 10.0); 01626 // uM = 0.0; 01627 uH = de(1) - uM; 01628 dHt1 = Kh * uH; 01629 // cout <<delta/deltaIn<<endl; 01630 // cout <<"TOP\n"; 01631 } 01632 01633 else 01634 { 01635 uM = (2.0 * L * (dTheta) * (L/hCurr)*pow(B/A, 2.0)); 01636 uM = (2.0 * L * (dTheta) * (L/hCurr)*pow(B/A, 2.0)) * pow(1.0 - delta/deltaIn, 10.0); 01637 // uM = 0.0; 01638 uH = de(1) - uM; 01639 dHt1 = Kh * uH; 01640 // cout <<delta/deltaIn<<endl; 01641 // cout <<"BOT\n"; 01642 } 01643 01644 */ 01645 01646 01647 // if (delta/deltaIn < 0.1) 01648 // delta = 0.1*deltaIn; 01649 01650 // if (delta/deltaIn > 0.5) 01651 // delta = deltaIn; 01652 01653 01654 Mmaxpast = delta/deltaIn; 01655 01656 if (delta/deltaIn > 1.0) 01657 uM = 0.0; 01658 else 01659 uM = (L * (dTheta-dMt/Kt) * (L/hCurr)*pow(B/A, 2.0)) 01660 // * (pow(1.0 - delta/deltaIn, 4.0)); 01661 * (pow(1.0 - delta/deltaIn, 10.0)); 01662 01663 01664 01665 // uM = 0.0; 01666 01667 uH = de(1) - uM; 01668 // dHt1 = Kh * uH * pow((delta+0.0001)/(deltaIn - delta), 1.0); 01669 // dHt1 = Kh * uH * pow(0.0 + delta/deltaIn, 0.2); 01670 // dHt1 = Kh * uH * pow(epsilon + delta/deltaIn, 2.0); 01671 // dHt1 = Kh * uH * pow(delta/deltaIn, 2.0); 01672 // dHt1 = Kh * uH * pow(delta/(deltaIn-delta), 2.0); 01673 dHt1 = Kh * uH * (delta/deltaIn); 01674 dHt1 = Kh * uH; 01675 dHt1 = Kh* ((double)(c2-c1)/noNodes) * uH; 01676 01677 01678 // if (uH == 0.0) exit(0); 01679 01680 01681 // cout <<"dMt and dHt = "<<dMt<<" "<<dHt<<endl; 01682 // cout <<"delta and deltaIn = "<<delta<<" "<<deltaIn<<endl; 01683 01684 01685 // if (dHt1 == 0.0) 01686 // dHt1 = de(1); 01687 01688 01689 Fh2 = (sCommit(1)+dHt1)/Vult; 01690 Fm2 = (sCommit(2)+dMt)/Vult/L; 01691 01692 // if (Fh2*Fh2/A/A + Fm2*Fm2/B/B - 1.0 >= 0.0) 01693 if (Fh2*Fh2/A/A + Fm2*Fm2/B/B - 1.0 > 0.01) 01694 { 01695 // gate = 1; 01696 /* 01697 cout <<"dHt1 = "<<dHt1<<endl; 01698 cout <<"dHt = "<<dHt<<endl; 01699 cout <<"dMt = "<<dMt<<endl; 01700 cout <<"Fh2 = "<<Fh2<<endl; 01701 cout <<"Fm2 = "<<Fm2<<endl; 01702 cout <<"A = "<<A<<endl; 01703 cout <<"B = "<<B<<endl; 01704 cout <<"Fh2*Fh2/A/A + Fm2*Fm2/B/B = "<<Fh2*Fh2/A/A + Fm2*Fm2/B/B<<endl; 01705 */ 01706 01707 // dHt1 = dMt / ((hCurr+hPrev)/2.0); 01708 // dHt1 = uH; 01709 01710 /* 01711 cout <<"delta = "<<delta<<endl; 01712 cout <<"deltaIn = "<<deltaIn<<endl; 01713 cout <<"dist1 = "<<dist1<<endl; 01714 cout <<"dist2 = "<<dist2<<endl; 01715 cout <<"dHt1 = "<<dHt1<<endl; 01716 cout <<"dHt = "<<dHt<<endl; 01717 cout <<"dMt = "<<dMt<<endl; 01718 cout <<"hCurr = "<<hCurr<<endl; 01719 cout <<"hPrev = "<<hPrev<<endl; 01720 */ 01721 01722 01723 /* 01724 if (delta/deltaIn < 0.05) 01725 dHt1 = uH; 01726 else 01727 dHt1 = dHt; 01728 */ 01729 01730 /* 01731 if (-1.0 * pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5) < Fh2) 01732 dHt1 = -1.0 * pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)*Vult - sCommit(1); 01733 else if (1.0 * pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5) > Fh2) 01734 dHt1 = pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)*Vult - sCommit(1); 01735 else 01736 exit (0); 01737 */ 01738 01739 01740 // dMt = 0.0; 01741 // dHt1 = 0.0; 01742 01743 01744 01745 // gate = 0; 01746 } 01747 else 01748 gate = 101; 01749 01750 01751 // cout <<"gate = "<<gate<<endl; 01752 01753 01754 /* 01755 if ((fabs(dHt1) < epsilon) || (fabs(dHt) < epsilon)) 01756 { 01757 dHt1 = dMt / ((hCurr+hPrev)/2.0); 01758 gate = 0; 01759 } 01760 */ 01761 01762 /* 01763 cout <<"dHt = "<<dHt<<endl; 01764 cout <<"dHt1 = "<<dHt1<<endl; 01765 cout <<"delta = "<<delta<<endl; 01766 */ 01767 01768 01769 /* 01770 if (dHt1 > dHt) 01771 dHt += 1.0; 01772 else if (dHt1 < dHt) 01773 dHt -= 1.0; 01774 */ 01775 01776 01777 01778 if ((dHt - dHt1 > 0.0) && (hNew < 0.0)) 01779 gate = 0; 01780 01781 01782 01783 01784 /* 01785 Fh2 = (sCommit(1)+dHt1)/Vult; 01786 Fm2 = (sCommit(2)+dMt)/Vult/L; 01787 01788 if (Fh2*Fh2/A/A + Fm2*Fm2/B/B - 1.0 > 0.0) 01789 { 01790 dHt1 = 0.0; 01791 } 01792 01793 */ 01794 01795 01796 /* 01797 Fm2 = (sCommit(2)+dMt)/Vult/L; 01798 dHt = pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)*Vult; 01799 01800 if ((sCommit(1)+dHt1 > dHt) && (dHt1 > 0.0)) 01801 // dHt1 = dHt - sCommit(1); 01802 // dHt1 = dMt/hCurr; 01803 dHt1 = 0.0; 01804 if ((sCommit(1)+dHt1 < -1.0*dHt) && (dHt1 < 0.0)) 01805 // dHt1 = -1.0*dHt - sCommit(1); 01806 // dHt1 = dMt/hCurr; 01807 dHt1 = 0.0; 01808 */ 01809 01810 01811 // if (fabs(dHt) > 100000) gate = 0; 01812 01813 01814 // } while ((fabs(dHt-dHt1) > 100.0) && (gate)); 01815 } while (0); 01816 01817 // if (dHt1 != 0.0) 01818 dHt = dHt1; 01819 01820 01821 01822 01823 Fh2 = (sCommit(1)+dHt)/Vult; 01824 Fm2 = (sCommit(2)+dMt)/Vult/L; 01825 01826 01827 01828 01829 dHt1 = fabs(pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)); 01830 if (Fh2 > dHt1) 01831 dHt = dHt1*Vult - sCommit(1); 01832 if (Fh2 < -1.0*dHt1) 01833 dHt = -1.0*dHt1*Vult - sCommit(1); 01834 01835 01836 /* 01837 dHt1 = fabs(pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)); 01838 if (Fh2 > dHt1) 01839 dHt = (Fv/2.0)*(1.0-Fv)/pow(B*B/A/A+hCurr*hCurr/L/L, 0.5)*Vult - sCommit(1); 01840 if (Fh2 < -1.0*dHt1) 01841 dHt = -1.0*(Fv/2.0)*(1.0-Fv)/pow(B*B/A/A+hCurr*hCurr/L/L, 0.5)*Vult - sCommit(1); 01842 */ 01843 01844 01845 01846 01847 if (du == 0.0) 01848 dHt = 0.0; 01849 01850 01851 01852 01853 01854 01855 if (Fh2*Fh2/A/A + Fm2*Fm2/B/B - 1.0 > 0.0) 01856 { 01857 01858 /* 01859 01860 // if (dHt > 0.0) 01861 if (Fh2 > 0.0) 01862 dHt = fabs(pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)*Vult)-sCommit(1); 01863 // if (dHt < 0.0) 01864 else if (Fh2 < 0.0) 01865 dHt = -1.0 * fabs(pow((A*A*B*B - Fm2*Fm2*A*A)/(B*B), 0.5)*Vult)-sCommit(1); 01866 01867 01868 01869 01870 if (dMt > 0.0) 01871 dMt = pow(pow(A*B*hCurr, 2.0)/(A*A*hCurr*hCurr+B*B*L*L), 0.5)*Vult*L - sCommit(2); 01872 else if (dMt < 0.0) 01873 dMt = -1.0*pow(pow(A*B*hCurr, 2.0)/(A*A*hCurr*hCurr+B*B*L*L), 0.5)*Vult*L - sCommit(2); 01874 dHt = dMt/hCurr; 01875 01876 01877 cout<<"In bottom\n"; 01878 cout <<"sCommit(1) and (2) = "<<sCommit(1)<<" "<<sCommit(2)<<endl; 01879 cout <<"dHt and dMt = "<<dHt<<" "<<dMt<<endl; 01880 cout <<"sCommit(0) and dVt = "<<sCommit(0)<<" "<<dVt<<endl; 01881 // exit(0); 01882 01883 */ 01884 01885 } 01886 01887 01888 01889 01890 01892 01893 01894 01895 /* 01896 if (isOver == 0) 01897 ds = 0.0; 01898 01899 dVt = Kv * (de(0) - ds); 01900 01901 01902 if (sCommit(0)+dVt > Vult) 01903 dVt = 0.99*Vult - sCommit(0); 01904 if (sCommit(0)+dVt < 0.0) 01905 dVt = -1.0*sCommit(0) + 0.01*Vult; 01906 */ 01907 01908 01909 if (isOver == 0) 01910 01911 for (int i = 0; i <= noNodes; i++) 01912 { 01913 01914 foot[i][0] = de(0); 01915 soilMin[i][0] = de(0); 01916 soilMax[i][0] = de(0); 01917 pressure[i][0] = dVt/Vult; 01918 pressMax[i][0] = dVt/Vult; 01919 } 01920 01921 } 01922 01923 01924 01925 01926 01927 // setup the stiffness matrix 01928 01929 deModel = de; 01930 ks.Zero(); 01931 01932 01933 if (de(0)-ds != 0.0) 01934 ks(0,0) = dVt/(de(0)-ds); 01935 else 01936 ks(0,0) = Kv; 01937 01938 // if (ks(0,0) == 0.0) 01939 // ks(0,0) = 1.0; 01940 01941 if ((dTheta != 0.0) && (de(0)-ds != 0.0)) 01942 ks(0,2) = -1.0 * (ds/dTheta) * (dVt/(de(0)-ds)); 01943 else if (dTheta != 0.0) 01944 ks(0,2) = -1.0 * Kv * (de(0)/dTheta); 01945 else 01946 ks(0,2) = 0.0; 01947 01948 01949 if (isOver != 0) 01950 { 01951 if (uH != 0.0) 01952 ks(1,1) = dHt/uH; 01953 else 01954 { 01955 ks(1,1) = Kh; 01956 } 01957 } 01958 else 01959 { 01960 if (de(1) != 0.0) 01961 ks(1,1) = dHt/de(1); 01962 else 01963 { 01964 ks(1,1) = Kh; 01965 } 01966 } 01967 01968 01969 01970 if ((dTheta != 0.0) && (uH != 0.0)) 01971 ks(1,2) = -1.0 * (uM/uH) * (dHt/dTheta); 01972 else if (dTheta != 0.0) 01973 ks(1,2) = -1.0 * Kh * (de(1)/dTheta); 01974 else 01975 ks(1,2) = 0.0; 01976 01977 01978 01979 if (de(2) != 0.0) 01980 ks(2,2) = dMt/de(2); 01981 else 01982 ks(2,2) = Kt; 01983 01984 01985 01986 01987 // make sure det|ks| is non-zero 01988 01989 detKs = ks(0,0) * (ks(1,1)*ks(2,2) - ks(1,2)*ks(2,1)) 01990 + ks(0,1) * (ks(2,0)*ks(1,2) - ks(1,0)*ks(2,2)) 01991 + ks(0,2) * (ks(1,0)*ks(2,1) - ks(2,0)*ks(1,1)); 01992 01993 01994 01995 01996 01997 delete [] footTemp; 01998 delete [] soilMinTemp; 01999 delete [] soilMaxTemp; 02000 delete [] pressureTemp; 02001 delete [] pressMaxTemp; 02002 02003 delete [] ddH; 02004 02005 return 0; 02006 } 02007 02008 02009 const Vector & 02010 SoilFootingSection2d::getSectionDeformation (void) 02011 { 02012 return e; 02013 } 02014 02015 const Vector & 02016 SoilFootingSection2d::getStressResultant (void) 02017 { 02018 return s; 02019 } 02020 02021 const Matrix & 02022 SoilFootingSection2d::getSectionTangent(void) 02023 { 02024 return ks; 02025 } 02026 02027 const Matrix & 02028 SoilFootingSection2d::getSectionFlexibility (void) 02029 { 02030 static Matrix fs(3,3); 02031 // invert2by2Matrix(ks, fs); 02032 02033 return fs; 02034 } 02035 02036 const ID& 02037 SoilFootingSection2d::getType () 02038 { 02039 return code; 02040 } 02041 02042 int 02043 SoilFootingSection2d::getOrder () const 02044 { 02045 return 3; 02046 } 02047 02048 02049 const Matrix& 02050 SoilFootingSection2d::getInitialTangent() 02051 { 02052 static Matrix IniTan(3,3); 02053 02054 return IniTan; 02055 } 02056 02057 int 02058 SoilFootingSection2d::sendSelf(int commitTag, Channel &theChannel) 02059 { 02060 return -1; 02061 } 02062 02063 int 02064 SoilFootingSection2d::recvSelf(int commitTag, Channel &theChannel, 02065 FEM_ObjectBroker &theBroker) 02066 { 02067 return -1; 02068 } 02069 02070 void 02071 SoilFootingSection2d::Print(OPS_Stream &s, int flag) 02072 { 02073 s << "YieldSurfaceSection2d, tag: " << this->getTag() << endln; 02074 s << "\tSection Force:" << sCommit; 02075 s << "\tSection Defom:" << eCommit; 02076 } 02077
Generated on Mon Oct 23 15:05:19 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents