CurvaturesBySearchAlgorithm.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.4 $ 00026 // $Date: 2003/03/04 00:38:55 $ 00027 // $Source: /usr/local/cvs/OpenSees/SRC/reliability/analysis/curvature/CurvaturesBySearchAlgorithm.cpp,v $ 00028 00029 00030 // 00031 // Written by Terje Haukaas (haukaas@ce.berkeley.edu) 00032 // 00033 00034 #include <CurvaturesBySearchAlgorithm.h> 00035 #include <FindCurvatures.h> 00036 #include <LimitStateFunction.h> 00037 #include <FindDesignPointAlgorithm.h> 00038 #include <RandomVariable.h> 00039 #include <Vector.h> 00040 00041 00042 CurvaturesBySearchAlgorithm::CurvaturesBySearchAlgorithm(int passedNumberOfCurvatures, 00043 FindDesignPointAlgorithm *passedFindDesignPointAlgorithm) 00044 :FindCurvatures(), curvatures(passedNumberOfCurvatures) 00045 { 00046 numberOfCurvatures = passedNumberOfCurvatures; 00047 theFindDesignPointAlgorithm = passedFindDesignPointAlgorithm; 00048 } 00049 00050 CurvaturesBySearchAlgorithm::~CurvaturesBySearchAlgorithm() 00051 { 00052 } 00053 00054 00055 int 00056 CurvaturesBySearchAlgorithm::computeCurvatures(ReliabilityDomain *theReliabilityDomain) 00057 { 00058 00059 // Initial declarations 00060 Vector last_u; 00061 Vector secondLast_u; 00062 Vector lastAlpha; 00063 Vector secondLastAlpha; 00064 Vector lastDirection; 00065 Vector uLastMinus_u = last_u - secondLast_u; 00066 double signumProduct; 00067 double alphaProduct; 00068 double sumSquared; 00069 double norm_uLastMinus_u; 00070 00071 // "Download" limit-state function from reliability domain 00072 int lsf = theReliabilityDomain->getTagOfActiveLimitStateFunction(); 00073 LimitStateFunction *theLimitStateFunction = 00074 theReliabilityDomain->getLimitStateFunctionPtr(lsf); 00075 00076 // The design point in the original space 00077 Vector x_star = theLimitStateFunction->designPoint_x_inOriginalSpace; 00078 00079 // Number of random variables 00080 int nrv = x_star.Size(); 00081 00082 // Declare vector to store all principal axes 00083 Vector principalAxes( nrv * numberOfCurvatures ); 00084 00085 // Start point of new searches: Perturb x^star by 10% of each standard deviation 00086 Vector startPoint(nrv); 00087 RandomVariable *aRandomVariable; 00088 int numberOfRandomVariables = theReliabilityDomain->getNumberOfRandomVariables(); 00089 int i; 00090 for (i=0; i<nrv; i++) { 00091 aRandomVariable = theReliabilityDomain->getRandomVariablePtr(i+1); 00092 startPoint(i) = aRandomVariable->getStartValue(); 00093 } 00094 00095 // Compute curvatures by repeated searches 00096 for (i=0; i<numberOfCurvatures; i++) { 00097 00098 // Get hold of 'u', 'alpha' and search direction at the two last steps 00099 last_u = theLimitStateFunction->designPoint_u_inStdNormalSpace; 00100 secondLast_u = theLimitStateFunction->secondLast_u; 00101 lastAlpha = theLimitStateFunction->normalizedNegativeGradientVectorAlpha; 00102 secondLastAlpha = theLimitStateFunction->secondLastAlpha; 00103 lastDirection = theLimitStateFunction->lastSearchDirection; 00104 00105 // Compute curvature according to Der Kiureghian & De Stefano (1992), Eq.26: 00106 00107 // Initial computations 00108 uLastMinus_u = last_u - secondLast_u; 00109 signumProduct = secondLastAlpha ^ uLastMinus_u; 00110 alphaProduct = secondLastAlpha ^ lastAlpha; 00111 sumSquared = 0.0; 00112 00113 // Compute norm of the difference vector 00114 for ( int k=0; k<last_u.Size(); k++ ) { 00115 sumSquared += uLastMinus_u(k)*uLastMinus_u(k); 00116 } 00117 norm_uLastMinus_u = sqrt(sumSquared); 00118 00119 // Check sign and compute curvature 00120 if (fabs(signumProduct)==(signumProduct)) { 00121 curvatures(i) = acos(alphaProduct) / norm_uLastMinus_u; 00122 } 00123 else { 00124 curvatures(i) = -acos(alphaProduct) / norm_uLastMinus_u; 00125 } 00126 00127 // Run a new search in the subspace orthogonal to previous principal directions 00128 if ( i!=(numberOfCurvatures-1) ) { 00129 00130 // Store all previous principal axes in a vector 00131 for (int j=0; j<nrv; j++ ) { 00132 principalAxes((i*nrv)+j) = lastDirection(j); 00133 } 00134 00135 // To be completed... 00136 } 00137 } 00138 00139 return 0; 00140 } 00141 00142 00143 00144 Vector 00145 CurvaturesBySearchAlgorithm::getCurvatures() 00146 { 00147 return curvatures; 00148 } 00149 00150 00151 |