| 67 | } |
| 68 | |
| 69 | int |
| 70 | SecantLineSearch::search(double s0, |
| 71 | double s1, |
| 72 | LinearSOE &theSOE, |
| 73 | IncrementalIntegrator &theIntegrator) |
| 74 | { |
| 75 | double r0 = 0.0; |
| 76 | |
| 77 | if ( s0 != 0.0 ) |
| 78 | r0 = fabs( s1 / s0 ); |
| 79 | |
| 80 | if (r0 <= tolerance ) |
| 81 | return 0; // Line Search Not Required Residual Decrease Less Than Tolerance |
| 82 | |
| 83 | if (s1 == s0) |
| 84 | return 0; // Secant will have a divide-by-zero if continue |
| 85 | |
| 86 | // set some variables |
| 87 | double eta = 1.0; |
| 88 | double s = s1; |
| 89 | double etaJ = 1.0; |
| 90 | double etaJm1 = 0.0; |
| 91 | double sJ = s1; |
| 92 | double sJm1 = s0; |
| 93 | double r = r0; |
| 94 | |
| 95 | const Vector &dU = theSOE.getX(); |
| 96 | |
| 97 | if (printFlag == 0) { |
| 98 | opserr << "Secant Line Search - initial: " |
| 99 | << " eta(0) : " << eta << " , Ratio |s/s0| = " << r0 << endln; |
| 100 | } |
| 101 | |
| 102 | // perform the secant iterations: |
| 103 | // |
| 104 | // eta(j+1) = eta(j) - s(j) * (eta(j-1)-eta(j)) |
| 105 | // ------------------------ |
| 106 | // s(j-1) - s(j) |
| 107 | |
| 108 | int count = 0; //initial value of iteration counter |
| 109 | while ( r > tolerance && count < maxIter ) { |
| 110 | |
| 111 | count++; |
| 112 | |
| 113 | eta = etaJ - sJ * (etaJm1-etaJ) / (sJm1 - sJ); |
| 114 | |
| 115 | //-- want to put limits on eta and stop solution blowing up |
| 116 | if (eta > maxEta) eta = maxEta; |
| 117 | if (r > r0 ) eta = 1.0; |
| 118 | if (eta < minEta) eta = minEta; |
| 119 | |
| 120 | //update the incremental difference in response and determine new unbalance |
| 121 | if (eta == etaJ) |
| 122 | break; // no change in response |
| 123 | |
| 124 | *x = dU; |
| 125 | *x *= eta-etaJ; |
| 126 |
no test coverage detected