My title

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