ArcLength1.cpp

Go 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.4 $
00022 // $Date: 2003/02/14 23:00:46 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/integrator/ArcLength1.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/analysis/integrator/ArcLength1.C
00027 // 
00028 // Written: fmk 
00029 // Created: 07/98
00030 // Revision: A
00031 //
00032 // Description: This file contains the class definition for ArcLength1.
00033 // ArcLength1 is an algorithmic class for perfroming a static analysis
00034 // using the arc length scheme, that is within a load step the follwing
00035 // constraint is enforced: dU^TdU + alpha^2*dLambda^2 = ArcLength1^2
00036 // where dU is change in nodal displacements for step, dLambda is
00037 // change in applied load and ArcLength1 is a control parameter.
00038 //
00039 // What: "@(#) ArcLength1.C, revA"
00040 
00041 
00042 #include <ArcLength1.h>
00043 #include <AnalysisModel.h>
00044 #include <LinearSOE.h>
00045 #include <Vector.h>
00046 #include <Channel.h>
00047 #include <math.h>
00048 
00049 ArcLength1::ArcLength1(double arcLength, double alpha)
00050 :StaticIntegrator(INTEGRATOR_TAGS_ArcLength1),
00051  arcLength2(arcLength*arcLength), alpha2(alpha*alpha),
00052  deltaUhat(0), deltaUbar(0), deltaU(0), deltaUstep(0), 
00053  phat(0), deltaLambdaStep(0.0), currentLambda(0.0), 
00054  signLastDeltaLambdaStep(1)
00055 {
00056 
00057 }
00058 
00059 ArcLength1::~ArcLength1()
00060 {
00061     // delete any vector object created
00062     if (deltaUhat != 0)
00063         delete deltaUhat;
00064     if (deltaU != 0)
00065         delete deltaU;
00066     if (deltaUstep != 0)
00067         delete deltaUstep;
00068     if (deltaUbar != 0)
00069         delete deltaUbar;
00070         if (phat != 0)
00071         delete phat;
00072 }
00073 
00074 int
00075 ArcLength1::newStep(void)
00076 {
00077     // get pointers to AnalysisModel and LinearSOE
00078     AnalysisModel *theModel = this->getAnalysisModelPtr();
00079     LinearSOE *theLinSOE = this->getLinearSOEPtr();    
00080     if (theModel == 0 || theLinSOE == 0) {
00081         opserr << "WARNING ArcLength1::newStep() ";
00082         opserr << "No AnalysisModel or LinearSOE has been set\n";
00083         return -1;
00084     }
00085 
00086     // get the current load factor
00087     currentLambda = theModel->getCurrentDomainTime();
00088 
00089     if (deltaLambdaStep < 0)
00090         signLastDeltaLambdaStep = -1;
00091     else
00092         signLastDeltaLambdaStep = +1;
00093 
00094     // determine dUhat
00095     this->formTangent();
00096     theLinSOE->setB(*phat);
00097     theLinSOE->solve();
00098     (*deltaUhat) = theLinSOE->getX();
00099     Vector &dUhat = *deltaUhat;
00100     
00101     // determine delta lambda(1) == dlambda
00102     double dLambda = sqrt(arcLength2/((dUhat^dUhat)+alpha2));
00103     dLambda *= signLastDeltaLambdaStep; // base sign of load change
00104                                         // on what was happening last step
00105     deltaLambdaStep = dLambda;
00106     currentLambda += dLambda;
00107 
00108     // determine delta U(1) == dU
00109     (*deltaU) = dUhat;
00110     (*deltaU) *= dLambda;
00111     (*deltaUstep) = (*deltaU);
00112 
00113     // update model with delta lambda and delta U
00114     theModel->incrDisp(*deltaU);    
00115     theModel->applyLoadDomain(currentLambda);    
00116     theModel->updateDomain();
00117 
00118     return 0;
00119 }
00120 
00121 int
00122 ArcLength1::update(const Vector &dU)
00123 {
00124     AnalysisModel *theModel = this->getAnalysisModelPtr();
00125     LinearSOE *theLinSOE = this->getLinearSOEPtr();    
00126     if (theModel == 0 || theLinSOE == 0) {
00127         opserr << "WARNING ArcLength1::update() ";
00128         opserr << "No AnalysisModel or LinearSOE has been set\n";
00129         return -1;
00130     }
00131 
00132     (*deltaUbar) = dU; // have to do this as the SOE is gonna change
00133 
00134     // determine dUhat    
00135     theLinSOE->setB(*phat);
00136     theLinSOE->solve();
00137     (*deltaUhat) = theLinSOE->getX();    
00138 
00139     // determine delta lambda(i)
00140     double a = (*deltaUstep)^(*deltaUbar);
00141     double b = (*deltaUstep)^(*deltaUhat) + alpha2*deltaLambdaStep;
00142     if (b == 0) {
00143       opserr << "ArcLength1::update() - zero denominator,";
00144       opserr << " alpha was set to 0.0 and zero reference load\n";
00145       return -1;
00146     }
00147     double dLambda = -a/b;
00148 
00149     // determine delta U(i)
00150     (*deltaU) = (*deltaUbar);    
00151     deltaU->addVector(1.0, *deltaUhat,dLambda);
00152     
00153     // update dU and dlambda
00154     (*deltaUstep) += *deltaU;
00155     deltaLambdaStep += dLambda;
00156     currentLambda += dLambda;
00157 
00158     // update the model
00159     theModel->incrDisp(*deltaU);    
00160     theModel->applyLoadDomain(currentLambda);    
00161     theModel->updateDomain();
00162     
00163     // set the X soln in linearSOE to be deltaU for convergence Test
00164     theLinSOE->setX(*deltaU);
00165 
00166     return 0;
00167 }
00168 
00169 
00170 
00171 int 
00172 ArcLength1::domainChanged(void)
00173 {
00174     // we first create the Vectors needed
00175     AnalysisModel *theModel = this->getAnalysisModelPtr();
00176     LinearSOE *theLinSOE = this->getLinearSOEPtr();    
00177     if (theModel == 0 || theLinSOE == 0) {
00178         opserr << "WARNING ArcLength1::update() ";
00179         opserr << "No AnalysisModel or LinearSOE has been set\n";
00180         return -1;
00181     }    
00182     int size = theModel->getNumEqn(); // ask model in case N+1 space
00183 
00184     if (deltaUhat == 0 || deltaUhat->Size() != size) { // create new Vector
00185         if (deltaUhat != 0)
00186             delete deltaUhat;   // delete the old
00187         deltaUhat = new Vector(size);
00188         if (deltaUhat == 0 || deltaUhat->Size() != size) { // check got it
00189             opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00190             opserr << " deltaUhat Vector of size " << size << endln;
00191             exit(-1);
00192         }
00193     }
00194 
00195     if (deltaUbar == 0 || deltaUbar->Size() != size) { // create new Vector
00196         if (deltaUbar != 0)
00197             delete deltaUbar;   // delete the old
00198         deltaUbar = new Vector(size);
00199         if (deltaUbar == 0 || deltaUbar->Size() != size) { // check got it
00200             opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00201             opserr << " deltaUbar Vector of size " << size << endln;
00202             exit(-1);
00203         }
00204     }
00205 
00206     
00207     if (deltaU == 0 || deltaU->Size() != size) { // create new Vector
00208         if (deltaU != 0)
00209             delete deltaU;   // delete the old
00210         deltaU = new Vector(size);
00211         if (deltaU == 0 || deltaU->Size() != size) { // check got it
00212             opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00213             opserr << " deltaU Vector of size " << size << endln;
00214             exit(-1);
00215         }
00216     }
00217 
00218     if (deltaUstep == 0 || deltaUstep->Size() != size) { 
00219         if (deltaUstep != 0)
00220             delete deltaUstep;  
00221         deltaUstep = new Vector(size);
00222         if (deltaUstep == 0 || deltaUstep->Size() != size) { 
00223             opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00224             opserr << " deltaUstep Vector of size " << size << endln;
00225             exit(-1);
00226         }
00227     }
00228 
00229     if (phat == 0 || phat->Size() != size) { 
00230         if (phat != 0)
00231             delete phat;  
00232         phat = new Vector(size);
00233         if (phat == 0 || phat->Size() != size) { 
00234             opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00235             opserr << " phat Vector of size " << size << endln;
00236             exit(-1);
00237         }
00238     }    
00239 
00240     // now we have to determine phat
00241     // do this by incrementing lambda by 1, applying load
00242     // and getting phat from unbalance.
00243     currentLambda = theModel->getCurrentDomainTime();
00244     currentLambda += 1.0;
00245     theModel->applyLoadDomain(currentLambda);    
00246     this->formUnbalance(); // NOTE: this assumes unbalance at last was 0
00247     (*phat) = theLinSOE->getB();
00248     currentLambda -= 1.0;
00249     theModel->setCurrentDomainTime(currentLambda);    
00250     
00251     return 0;
00252 }
00253 
00254 int
00255 ArcLength1::sendSelf(int cTag,
00256                     Channel &theChannel)
00257 {
00258   Vector data(5);
00259   data(0) = arcLength2;
00260   data(1) = alpha2;
00261   data(2) = deltaLambdaStep;
00262   data(3) = currentLambda;
00263   data(4)  = signLastDeltaLambdaStep;
00264 
00265   if (theChannel.sendVector(this->getDbTag(), cTag, data) < 0) {
00266       opserr << "ArcLength1::sendSelf() - failed to send the data\n";
00267       return -1;
00268   }
00269   return 0;
00270 }
00271 
00272 
00273 int
00274 ArcLength1::recvSelf(int cTag,
00275                     Channel &theChannel, FEM_ObjectBroker &theBroker)
00276 {
00277   Vector data(5);
00278   if (theChannel.recvVector(this->getDbTag(), cTag, data) < 0) {
00279       opserr << "ArcLength1::sendSelf() - failed to send the data\n";
00280       return -1;
00281   }      
00282 
00283   // set the data
00284   arcLength2 = data(0);
00285   alpha2 = data(1);
00286   deltaLambdaStep = data(2);
00287   currentLambda = data(3);
00288   signLastDeltaLambdaStep = data(4);
00289   return 0;
00290 }
00291 
00292 void
00293 ArcLength1::Print(OPS_Stream &s, int flag)
00294 {
00295     AnalysisModel *theModel = this->getAnalysisModelPtr();
00296     if (theModel != 0) {
00297         double cLambda = theModel->getCurrentDomainTime();
00298         s << "\t ArcLength1 - currentLambda: " << cLambda;
00299         s << "  ArcLength1: " << sqrt(arcLength2) <<  "  alpha: ";
00300         s << sqrt(alpha2) << endln;
00301     } else 
00302         s << "\t ArcLength1 - no associated AnalysisModel\n";
00303 }
00304 
00305 
00306 
00307 
00308 
00309 
00310 
00311 

Generated on Mon Oct 23 15:04:58 2006 for OpenSees by doxygen 1.5.0