FORMAnalysis.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 2001, 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 ** Reliability module developed by: ** 00020 ** Terje Haukaas (haukaas@ce.berkeley.edu) ** 00021 ** Armen Der Kiureghian (adk@ce.berkeley.edu) ** 00022 ** ** 00023 ** ****************************************************************** */ 00024 00025 // $Revision: 1.7 $ 00026 // $Date: 2003/10/27 23:45:41 $ 00027 // $Source: /usr/local/cvs/OpenSees/SRC/reliability/analysis/analysis/FORMAnalysis.cpp,v $ 00028 00029 00030 // 00031 // Written by Terje Haukaas (haukaas@ce.berkeley.edu) 00032 // 00033 00034 #include <FORMAnalysis.h> 00035 #include <FindDesignPointAlgorithm.h> 00036 #include <ReliabilityDomain.h> 00037 #include <ReliabilityAnalysis.h> 00038 #include <Vector.h> 00039 #include <Matrix.h> 00040 #include <MatrixOperations.h> 00041 #include <NormalRV.h> 00042 #include <RandomVariable.h> 00043 #include <math.h> 00044 #include <ProbabilityTransformation.h> 00045 00046 #include <fstream> 00047 #include <iomanip> 00048 #include <iostream> 00049 using std::ifstream; 00050 using std::ios; 00051 using std::setw; 00052 using std::setprecision; 00053 using std::setiosflags; 00054 00055 00056 FORMAnalysis::FORMAnalysis(ReliabilityDomain *passedReliabilityDomain, 00057 FindDesignPointAlgorithm *passedFindDesignPointAlgorithm, 00058 ProbabilityTransformation *passedProbabilityTransformation, 00059 TCL_Char *passedFileName, 00060 int p_relSensTag) 00061 :ReliabilityAnalysis() 00062 { 00063 theReliabilityDomain = passedReliabilityDomain; 00064 theFindDesignPointAlgorithm = passedFindDesignPointAlgorithm; 00065 theProbabilityTransformation = passedProbabilityTransformation; 00066 fileName = new char[256]; 00067 strcpy(fileName,passedFileName); 00068 relSensTag = p_relSensTag; 00069 } 00070 00071 00072 FORMAnalysis::~FORMAnalysis() 00073 { 00074 if (fileName != 0) 00075 delete [] fileName; 00076 } 00077 00078 00079 00080 int 00081 FORMAnalysis::analyze(void) 00082 { 00083 00084 // Alert the user that the FORM analysis has started 00085 opserr << "FORM Analysis is running ... " << endln; 00086 00087 00088 // Declare variables used in this method 00089 Vector xStar; 00090 Vector uStar; 00091 Vector alpha; 00092 Vector gamma; 00093 double stdv,mean; 00094 int i; 00095 double Go, Glast; 00096 Vector uSecondLast; 00097 Vector alphaSecondLast; 00098 Vector lastSearchDirection; 00099 double beta; 00100 double pf1; 00101 int numberOfEvaluations; 00102 int lsf; 00103 int numRV = theReliabilityDomain->getNumberOfRandomVariables(); 00104 int numLsf = theReliabilityDomain->getNumberOfLimitStateFunctions(); 00105 RandomVariable *aRandomVariable; 00106 LimitStateFunction *theLimitStateFunction; 00107 NormalRV *aStdNormRV=0; 00108 aStdNormRV = new NormalRV(1,0.0,1.0,0.0); 00109 Vector delta(numRV); 00110 Vector eta(numRV); 00111 Vector kappa(numRV); 00112 00113 00114 // Check if computer ran out of memory 00115 if (aStdNormRV==0) { 00116 opserr << "FORMAnalysis::analyze() - out of memory while instantiating internal objects." << endln; 00117 return -1; 00118 } 00119 00120 00121 // Open output file 00122 ofstream outputFile( fileName, ios::out ); 00123 00124 00125 // Loop over number of limit-state functions and perform FORM analysis 00126 for (lsf=1; lsf<=numLsf; lsf++ ) { 00127 00128 00129 // Inform the user which limit-state function is being evaluated 00130 opserr << "Limit-state function number: " << lsf << endln; 00131 00132 00133 // Set tag of "active" limit-state function 00134 theReliabilityDomain->setTagOfActiveLimitStateFunction(lsf); 00135 00136 00137 // Get the limit-state function pointer 00138 theLimitStateFunction = 0; 00139 lsf = theReliabilityDomain->getTagOfActiveLimitStateFunction(); 00140 theLimitStateFunction = theReliabilityDomain->getLimitStateFunctionPtr(lsf); 00141 if (theLimitStateFunction == 0) { 00142 opserr << "FORMAnalysis::analyze() - could not find" << endln 00143 << " limit-state function with tag #" << lsf << "." << endln; 00144 return -1; 00145 } 00146 00147 00148 // Find the design point 00149 if (theFindDesignPointAlgorithm->findDesignPoint(theReliabilityDomain) < 0){ 00150 opserr << "FORMAnalysis::analyze() - failed while finding the" << endln 00151 << " design point for limit-state function number " << lsf << "." << endln; 00152 outputFile << "#######################################################################" << endln; 00153 outputFile << "# FORM ANALYSIS RESULTS, LIMIT-STATE FUNCTION NUMBER " 00154 <<setiosflags(ios::left)<<setprecision(1)<<setw(4)<<lsf <<" #" << endln; 00155 outputFile << "# #" << endln; 00156 outputFile << "# No convergence! #" << endln; 00157 outputFile << "# #" << endln; 00158 outputFile << "#######################################################################" << endln << endln << endln; 00159 } 00160 else { 00161 00162 // Get results from the "find desingn point algorithm" 00163 xStar = theFindDesignPointAlgorithm->get_x(); 00164 uStar = theFindDesignPointAlgorithm->get_u(); 00165 alpha = theFindDesignPointAlgorithm->get_alpha(); 00166 gamma = theFindDesignPointAlgorithm->get_gamma(); 00167 i = theFindDesignPointAlgorithm->getNumberOfSteps(); 00168 Go = theFindDesignPointAlgorithm->getFirstGFunValue(); 00169 Glast = theFindDesignPointAlgorithm->getLastGFunValue(); 00170 uSecondLast = theFindDesignPointAlgorithm->getSecondLast_u(); 00171 alphaSecondLast = theFindDesignPointAlgorithm->getSecondLast_alpha(); 00172 lastSearchDirection = theFindDesignPointAlgorithm->getLastSearchDirection(); 00173 numberOfEvaluations = theFindDesignPointAlgorithm->getNumberOfEvaluations(); 00174 00175 00176 // Postprocessing 00177 beta = alpha ^ uStar; 00178 pf1 = 1.0 - aStdNormRV->getCDFvalue(beta); 00179 00180 00181 // Reliability sensitivity analysis wrt. mean/stdv 00182 if (relSensTag == 1) { 00183 Vector DuStarDmean; 00184 Vector DuStarDstdv; 00185 double dBetaDmean; 00186 double dBetaDstdv; 00187 00188 for ( int j=1; j<=numRV; j++ ) 00189 { 00190 DuStarDmean = theProbabilityTransformation->meanSensitivityOf_x_to_u(xStar,j); 00191 DuStarDstdv = theProbabilityTransformation->stdvSensitivityOf_x_to_u(xStar,j); 00192 dBetaDmean = alpha^DuStarDmean; 00193 dBetaDstdv = alpha^DuStarDstdv; 00194 aRandomVariable = theReliabilityDomain->getRandomVariablePtr(j); 00195 stdv = aRandomVariable->getStdv(); 00196 mean = aRandomVariable->getMean(); 00197 delta(j-1) = stdv * dBetaDmean; 00198 eta(j-1) = stdv * dBetaDstdv; 00199 // (Kappa is the sensitivity wrt. the coefficient of variation) 00200 if (mean != 0.0) { 00201 kappa(j-1) = -dBetaDmean*stdv/((stdv/mean)*(stdv/mean))+dBetaDstdv*mean; 00202 } 00203 else { 00204 kappa(j-1) = 0.0; 00205 } 00206 } 00207 // Don't scale them: 00208 // delta = delta * (1.0/delta.Norm()); 00209 // eta = eta * (1.0/eta.Norm()); 00210 // kappa = kappa * (1.0/kappa.Norm()); 00211 } 00212 00213 00214 // Store key results in the limit-state functions 00215 theLimitStateFunction->FORMReliabilityIndexBeta = beta; 00216 theLimitStateFunction->FORMProbabilityOfFailure_pf1 = pf1; 00217 theLimitStateFunction->designPoint_x_inOriginalSpace = xStar; 00218 theLimitStateFunction->designPoint_u_inStdNormalSpace = uStar; 00219 theLimitStateFunction->normalizedNegativeGradientVectorAlpha= alpha; 00220 theLimitStateFunction->importanceVectorGamma = gamma; 00221 theLimitStateFunction->numberOfStepsToFindDesignPointAlgorithm = i; 00222 theLimitStateFunction->GFunValueAtStartPt = Go; 00223 theLimitStateFunction->GFunValueAtEndPt = Glast; 00224 theLimitStateFunction->secondLast_u = uSecondLast; 00225 theLimitStateFunction->secondLastAlpha = alphaSecondLast; 00226 theLimitStateFunction->lastSearchDirection = lastSearchDirection; 00227 00228 00229 // Print FORM results to the output file 00230 outputFile << "#######################################################################" << endln; 00231 outputFile << "# FORM ANALYSIS RESULTS, LIMIT-STATE FUNCTION NUMBER " 00232 <<setiosflags(ios::left)<<setprecision(1)<<setw(4)<<lsf <<" #" << endln; 00233 outputFile << "# #" << endln; 00234 outputFile << "# Limit-state function value at start point: ......... " 00235 <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<Go 00236 << " #" << endln; 00237 outputFile << "# Limit-state function value at end point: ........... " 00238 <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<Glast 00239 << " #" << endln; 00240 outputFile << "# Number of steps: ................................... " 00241 <<setiosflags(ios::left)<<setw(12)<<i 00242 << " #" << endln; 00243 outputFile << "# Number of g-function evaluations: .................. " 00244 <<setiosflags(ios::left)<<setw(12)<<numberOfEvaluations 00245 << " #" << endln; 00246 outputFile << "# Reliability index beta: ............................ " 00247 <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<beta 00248 << " #" << endln; 00249 outputFile << "# FO approx. probability of failure, pf1: ............ " 00250 <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<pf1 00251 << " #" << endln; 00252 outputFile << "# #" << endln; 00253 outputFile << "# rv# x* u* alpha gamma delta eta #" << endln; 00254 outputFile.setf( ios::scientific, ios::floatfield ); 00255 for (int i=0; i<xStar.Size(); i++) { 00256 outputFile << "# " <<setw(3)<<(i+1)<<" "; 00257 outputFile.setf(ios::scientific, ios::floatfield); 00258 if (xStar(i)<0.0) { outputFile << "-"; } 00259 else { outputFile << " "; } 00260 outputFile <<setprecision(3)<<setw(11)<<fabs(xStar(i)); 00261 if (uStar(i)<0.0) { outputFile << "-"; } 00262 else { outputFile << " "; } 00263 outputFile <<setprecision(3)<<setw(11)<<fabs(uStar(i)); 00264 outputFile.unsetf( ios::scientific ); 00265 outputFile.setf(ios::fixed, ios::floatfield); 00266 00267 if (alpha(i)<0.0) { outputFile << "-"; } 00268 else { outputFile << " "; } 00269 outputFile<<setprecision(5)<<setw(8)<<fabs(alpha(i)); 00270 00271 if (gamma(i)<0.0) { outputFile << "-"; } 00272 else { outputFile << " "; } 00273 outputFile<<setprecision(5)<<setw(8)<<fabs(gamma(i)); 00274 00275 if (relSensTag == 1) { 00276 if (delta(i)<0.0) { outputFile << "-"; } 00277 else { outputFile << " "; } 00278 outputFile<<setprecision(5)<<setw(8)<<fabs(delta(i)); 00279 00280 if (eta(i)<0.0) { outputFile << "-"; } 00281 else { outputFile << " "; } 00282 outputFile<<setprecision(5)<<setw(8)<<fabs(eta(i)); 00283 00284 // Printing reliability sensitivity wrt. coefficient of variation 00285 // aRandomVariable = theReliabilityDomain->getRandomVariablePtr(i+1); 00286 // mean = aRandomVariable->getMean(); 00287 // if (mean==0.0) { outputFile << " - "; } 00288 // else { 00289 // if (kappa(i)<0.0) { outputFile << "-"; } 00290 // else { outputFile << " "; } 00291 // outputFile<<setprecision(5)<<setw(8)<<fabs(kappa(i)); 00292 // } 00293 } 00294 else { 00295 outputFile << " - - "; 00296 } 00297 00298 outputFile<<" #" << endln; 00299 } 00300 outputFile << "# #" << endln; 00301 outputFile << "#######################################################################" << endln << endln << endln; 00302 00303 00304 // Inform the user that we're done with this limit-state function 00305 opserr << "Done analyzing limit-state function " << lsf << ", beta=" << beta << endln; 00306 } 00307 00308 } 00309 00310 00311 // Clean up 00312 outputFile.close(); 00313 delete aStdNormRV; 00314 00315 // Print summary of results to screen (more here!!!) 00316 opserr << "FORMAnalysis completed." << endln; 00317 00318 return 0; 00319 } 00320
Generated on Mon Oct 23 15:05:25 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents