RegulaFalsiLineSearch.cppGo to the documentation of this file.00001 /* ****************************************************************** ** 00002 ** OpenSees - Open System for Earthquake Engineering Simulation ** 00003 ** Pacific Earthquake Engineering Research Center ** 00004 ** ** 00005 ** ** 00006 ** (C) Copyright 1999, The Regents of the University of California ** 00007 ** All Rights Reserved. ** 00008 ** ** 00009 ** Commercial use of this program without express permission of the ** 00010 ** University of California, Berkeley, is strictly prohibited. See ** 00011 ** file 'COPYRIGHT' in main directory for information on usage and ** 00012 ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** 00013 ** ** 00014 ** Developed by: ** 00015 ** Frank McKenna (fmckenna@ce.berkeley.edu) ** 00016 ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** 00017 ** Filip C. Filippou (filippou@ce.berkeley.edu) ** 00018 ** ** 00019 ** ****************************************************************** */ 00020 00021 // $Revision: 1.3 $ 00022 // $Date: 2003/04/02 22:02:33 $ 00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/algorithm/equiSolnAlgo/RegulaFalsiLineSearch.cpp,v $ 00024 00025 // Written: fmk 00026 // Created: 11/01 00027 // 00028 // What: "@(#)RegulaFalsiLineSearch.h, revA" 00029 00030 #include <RegulaFalsiLineSearch.h> 00031 #include <IncrementalIntegrator.h> 00032 #include <LinearSOE.h> 00033 #include <Channel.h> 00034 #include <FEM_ObjectBroker.h> 00035 #include <Vector.h> 00036 #include <math.h> 00037 00038 RegulaFalsiLineSearch::RegulaFalsiLineSearch(double tol, int mIter, double mnEta, double mxEta, int pFlag) 00039 :LineSearch(LINESEARCH_TAGS_RegulaFalsiLineSearch), 00040 x(0), tolerance(tol), maxIter(mIter), minEta(mnEta), maxEta(mxEta), printFlag(pFlag) 00041 { 00042 00043 } 00044 00045 RegulaFalsiLineSearch::~RegulaFalsiLineSearch() 00046 { 00047 if (x != 0) 00048 delete x; 00049 } 00050 00051 00052 int 00053 RegulaFalsiLineSearch::newStep(LinearSOE &theSOE) 00054 { 00055 const Vector &dU = theSOE.getX(); 00056 00057 if (x == 0) 00058 x = new Vector(dU); 00059 00060 if (x->Size() != dU.Size()) { 00061 delete x; 00062 x = new Vector(dU); 00063 } 00064 00065 return 0; 00066 } 00067 00068 int 00069 RegulaFalsiLineSearch::search(double s0, 00070 double s1, 00071 LinearSOE &theSOE, 00072 IncrementalIntegrator &theIntegrator) 00073 { 00074 double r0 = 0.0; 00075 00076 if ( s0 != 0.0 ) 00077 r0 = fabs( s1 / s0 ); 00078 00079 if (r0 <= tolerance ) 00080 return 0; // Line Search Not Required Residual Decrease Less Than Tolerance 00081 00082 if (s1 == s0) 00083 return 0; // RegulaFalsi will have a divide-by-zero error if continue 00084 00085 // set some variables 00086 double eta = 1.0; 00087 double s = s1; 00088 double etaU = 1.0; 00089 double etaL = 0.0; 00090 double sU = s1; 00091 double sL = s0; 00092 double r = r0; 00093 double etaJ = 1.0; 00094 00095 00096 const Vector &dU = theSOE.getX(); 00097 00098 if (printFlag == 0) { 00099 opserr << "RegulaFalsi Line Search - initial: " 00100 << " eta(0) : " << eta << " , Ratio |s/s0| = " << r0 << endln; 00101 } 00102 00103 // perform the secant iterations: 00104 // 00105 // eta(j+1) = eta(u) - s(u) * (eta(l) -eta(u)) 00106 // ------------------------ 00107 // s(l) - s(u) 00108 00109 int count = 0; //intial value of iteration counter 00110 while ( r > tolerance && count < maxIter ) { 00111 00112 count++; 00113 00114 eta = etaU - sU * (etaL-etaU) / (sL - sU); 00115 00116 00117 //-- want to put limits on eta(i) 00118 if (eta > maxEta) eta = maxEta; 00119 if ( r > r0 ) eta = 1.0; 00120 if (eta < minEta) eta = minEta; 00121 00122 00123 //update the incremental difference in response and determine new unbalance 00124 *x = dU; 00125 *x *= eta-etaJ; 00126 00127 if (theIntegrator.update(*x) < 0) { 00128 opserr << "WARNING RegulaFalsiLineSearch::search() -"; 00129 opserr << "the Integrator failed in update()\n"; 00130 return -1; 00131 } 00132 00133 if (theIntegrator.formUnbalance() < 0) { 00134 opserr << "WARNING RegulaFalsiLineSearch::search() -"; 00135 opserr << "the Integrator failed in formUnbalance()\n"; 00136 return -2; 00137 } 00138 00139 //new residual 00140 const Vector &ResidJ = theSOE.getB(); 00141 00142 //new value of s 00143 s = dU ^ ResidJ; 00144 00145 //new value of r 00146 r = fabs( s / s0 ); 00147 00148 00149 if (printFlag == 0) { 00150 opserr << "RegulaFalsi Line Search - iteration: " << count 00151 << " , eta(j) : " << eta << " , Ratio |sj/s0| = " << r << endln; 00152 } 00153 00154 00155 if (etaJ == eta) 00156 count = maxIter; 00157 00158 // set variables for next iteration 00159 etaJ = eta; 00160 00161 if (s*sU < 0.0) { 00162 etaL = eta; 00163 sL = s; 00164 } else if (s*sU == 0.0) 00165 count = maxIter; 00166 else { 00167 etaU = eta; 00168 sU = s; 00169 } 00170 00171 if (sL == sU) 00172 count = maxIter; 00173 00174 } //end while 00175 00176 // set X in the SOE for the revised dU, needed for convergence tests 00177 *x = dU; 00178 *x *= eta; 00179 theSOE.setX(*x); 00180 00181 return 0; 00182 } 00183 00184 00185 int 00186 RegulaFalsiLineSearch::sendSelf(int cTag, Channel &theChannel) 00187 { 00188 return 0; 00189 } 00190 00191 int 00192 RegulaFalsiLineSearch::recvSelf(int cTag, 00193 Channel &theChannel, 00194 FEM_ObjectBroker &theBroker) 00195 { 00196 return 0; 00197 } 00198 00199 00200 void 00201 RegulaFalsiLineSearch::Print(OPS_Stream &s, int flag) 00202 { 00203 if (flag == 0) { 00204 s << "RegulaFalsiLineSearch :: Line Search Tolerance = " << tolerance << endln; 00205 s << " max num Iterations = " << maxIter << endln; 00206 s << " max value on eta = " << maxEta << endln; 00207 } 00208 } 00209 00210 00211 00212 00213 00214 00215 00216 00217
Generated on Mon Oct 23 15:04:57 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents