SecantLineSearch.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/SecantLineSearch.cpp,v $ 00024 00025 // Written: fmk 00026 // Created: 11/01 00027 // 00028 // What: "@(#)SecantLineSearch.h, revA" 00029 00030 #include <SecantLineSearch.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 SecantLineSearch::SecantLineSearch(double tol, int mIter, double mnEta, double mxEta, int pFlag) 00039 :LineSearch(LINESEARCH_TAGS_SecantLineSearch), 00040 x(0), tolerance(tol), maxIter(mIter), minEta(mnEta), maxEta(mxEta), printFlag(pFlag) 00041 { 00042 00043 } 00044 00045 SecantLineSearch::~SecantLineSearch() 00046 { 00047 if (x != 0) 00048 delete x; 00049 } 00050 00051 00052 int 00053 SecantLineSearch::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 SecantLineSearch::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; // Secant will have a divide-by-zero if continue 00084 00085 // set some variables 00086 double eta = 1.0; 00087 double s = s1; 00088 double etaJ = 1.0; 00089 double etaJm1 = 0.0; 00090 double sJ = s1; 00091 double sJm1 = s0; 00092 double r = r0; 00093 00094 const Vector &dU = theSOE.getX(); 00095 00096 if (printFlag == 0) { 00097 opserr << "Secant Line Search - initial: " 00098 << " eta(0) : " << eta << " , Ratio |s/s0| = " << r0 << endln; 00099 } 00100 00101 // perform the secant iterations: 00102 // 00103 // eta(j+1) = eta(j) - s(j) * (eta(j-1)-eta(j)) 00104 // ------------------------ 00105 // s(j-1) - s(j) 00106 00107 int count = 0; //intial value of iteration counter 00108 while ( r > tolerance && count < maxIter ) { 00109 00110 count++; 00111 00112 eta = etaJ - sJ * (etaJm1-etaJ) / (sJm1 - sJ); 00113 00114 //-- want to put limits on eta and stop solution blowing up 00115 if (eta > maxEta) eta = maxEta; 00116 if (r > r0 ) eta = 1.0; 00117 if (eta < minEta) eta = minEta; 00118 00119 //update the incremental difference in response and determine new unbalance 00120 *x = dU; 00121 *x *= eta-etaJ; 00122 00123 if (theIntegrator.update(*x) < 0) { 00124 opserr << "WARNING SecantLineSearch::search() -"; 00125 opserr << "the Integrator failed in update()\n"; 00126 return -1; 00127 } 00128 00129 if (theIntegrator.formUnbalance() < 0) { 00130 opserr << "WARNING SecantLineSearch::search() -"; 00131 opserr << "the Integrator failed in formUnbalance()\n"; 00132 return -2; 00133 } 00134 00135 //new residual 00136 const Vector &ResidJ = theSOE.getB(); 00137 00138 //new value of s 00139 s = dU ^ ResidJ; 00140 00141 //new value of r 00142 r = fabs( s / s0 ); 00143 00144 if (printFlag == 0) { 00145 opserr << "Secant Line Search - iteration: " << count 00146 << " , eta(j) : " << eta << " , Ratio |sj/s0| = " << r << endln; 00147 } 00148 00149 if (etaJ == eta) 00150 count = maxIter; 00151 00152 // set variables for next iteration 00153 etaJm1 = etaJ; 00154 etaJ = eta; 00155 sJm1 = sJ; 00156 sJ = s; 00157 00158 if (sJm1 == sJ) 00159 count = maxIter; 00160 00161 } //end while 00162 00163 // set X in the SOE for the revised dU, needed for convergence tests 00164 *x = dU; 00165 *x *= eta; 00166 theSOE.setX(*x); 00167 00168 return 0; 00169 } 00170 00171 00172 int 00173 SecantLineSearch::sendSelf(int cTag, Channel &theChannel) 00174 { 00175 return 0; 00176 } 00177 00178 int 00179 SecantLineSearch::recvSelf(int cTag, 00180 Channel &theChannel, 00181 FEM_ObjectBroker &theBroker) 00182 { 00183 return 0; 00184 } 00185 00186 00187 void 00188 SecantLineSearch::Print(OPS_Stream &s, int flag) 00189 { 00190 if (flag == 0) { 00191 s << "SecantLineSearch :: Line Search Tolerance = " << tolerance << endln; 00192 s << " max num Iterations = " << maxIter << endln; 00193 s << " max value on eta = " << maxEta << endln; 00194 } 00195 } 00196 00197 00198 00199 00200 00201 00202 00203 00204
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