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.3 $
00022 // $Date: 2001/03/29 05:23:32 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/integrator/HHT.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/analysis/integrator/HHT.C
00027 // 
00028 // Written: fmk 
00029 // Created: 11/98
00030 // Revision: A
00031 //
00032 // Description: This file contains the implementation of the HHT class.
00033 //
00034 // What: "@(#) HHT.C, revA"
00035 
00036 #include <HHT.h>
00037 #include <FE_Element.h>
00038 #include <LinearSOE.h>
00039 #include <AnalysisModel.h>
00040 #include <Vector.h>
00041 #include <DOF_Group.h>
00042 #include <DOF_GrpIter.h>
00043 #include <AnalysisModel.h>
00044 #include <Channel.h>
00045 #include <FEM_ObjectBroker.h>
00046 
00047 HHT::HHT()
00048 :TransientIntegrator(INTEGRATOR_TAGS_HHT),
00049  alpha(0.5), gamma(1.0), beta(0), 
00050  rayleighDamping(false), alphaM(0.0), betaK(0.0), betaKi(0.0), betaKc(0.0),
00051  c1(0.0), c2(0.0), c3(0.0), 
00052  Ut(0), Utdot(0), Utdotdot(0),  U(0), Udot(0), Udotdot(0),
00053  Ualpha(0),Udotalpha(0)
00054 {
00055     
00056 }
00057 
00058 HHT::HHT(double _alpha)
00059 :TransientIntegrator(INTEGRATOR_TAGS_HHT),
00060  alpha(_alpha), gamma(1.5-_alpha), beta((2-_alpha)*(2-_alpha)*0.25),
00061  rayleighDamping(false), alphaM(0.0), betaK(0.0), betaKi(0.0), betaKc(0.0),
00062  c1(0.0), c2(0.0), c3(0.0), 
00063  Ut(0), Utdot(0), Utdotdot(0),  U(0), Udot(0), Udotdot(0),
00064  Ualpha(0),Udotalpha(0)
00065 {
00066 
00067 }
00068 
00069 HHT::HHT(double _alpha, double alpham, double betak, double betaki, double betakc)
00070 :TransientIntegrator(INTEGRATOR_TAGS_HHT),
00071  alpha(_alpha), gamma(1.5-_alpha), beta((2-_alpha)*(2-_alpha)*0.25),  
00072  rayleighDamping(true), 
00073  alphaM(alpham), betaK(betak), betaKi(betaki), betaKc(betakc),
00074  c1(0.0), c2(0.0), c3(0.0), 
00075  Ut(0), Utdot(0), Utdotdot(0),  U(0), Udot(0), Udotdot(0),
00076  Ualpha(0),Udotalpha(0)
00077 {
00078     if (alpham == 0.0 && betak == 0.0 && betaki == 0.0 && betakc == 0.0)
00079  rayleighDamping = false;
00080 }
00081 
00082 HHT::~HHT()
00083 {
00084   // clean up the memory created
00085   if (Ut != 0)
00086     delete Ut;
00087   if (Utdot != 0)
00088     delete Utdot;
00089   if (Utdotdot != 0)
00090     delete Utdotdot;
00091   if (U != 0)
00092     delete U;
00093   if (Udot != 0)
00094     delete Udot;
00095   if (Udotdot != 0)
00096     delete Udotdot;
00097   if (Ualpha != 0)
00098     delete Ualpha;
00099   if (Udotalpha != 0)
00100     delete Udotalpha;
00101 }
00102 
00103 
00104 int
00105 HHT::formEleResidual(FE_Element *theEle)
00106 {
00107   theEle->zeroResidual();
00108   if (rayleighDamping == false) {
00109       theEle->addRIncInertiaToResidual();
00110   } else {
00111       theEle->addRIncInertiaToResidual();
00112       theEle->addKtForce(*Udot,-betaK);
00113       theEle->addKcForce(*Udot, -betaKc);
00114       theEle->addKiForce(*Udot, -betaKi);
00115       theEle->addM_Force(*Udot,-alphaM);
00116   }    
00117   return 0;
00118 }    
00119 
00120 int
00121 HHT::formNodUnbalance(DOF_Group *theDof)
00122 {
00123   theDof->zeroUnbalance();
00124   if (rayleighDamping == false) 
00125       theDof->addPIncInertiaToUnbalance();
00126   else {
00127       theDof->addPIncInertiaToUnbalance();
00128       theDof->addM_Force(*Udot,-alphaM);
00129   }
00130   return 0;
00131 }    
00132 
00133 
00134 int
00135 HHT::initialize(void)
00136 {
00137   return 0;
00138 }
00139 
00140 
00141 
00142 int
00143 HHT::newStep(double deltaT)
00144 {
00145 
00146   if (beta == 0 || gamma == 0 ) {
00147     cerr << "HHT::newStep() - error in variable\n";
00148     cerr << "gamma = " << gamma << " beta= " << beta << endl;
00149     return -1;
00150   }
00151     
00152   if (deltaT <= 0.0) {
00153     cerr << "HHT::newStep() - error in variable\n";
00154     cerr << "dT = " << deltaT << endl;
00155     return -2; 
00156   }
00157   c1 = 1.0;
00158   c2 = gamma/(beta*deltaT);
00159   c3 = 1.0/(beta*deltaT*deltaT);
00160 
00161   AnalysisModel *theModel = this->getAnalysisModelPtr();
00162 
00163   
00164   // if alphaKc is specified .. 
00165   //    loop over FE_Elements getting them to do a setKc()
00166   if (rayleighDamping == true && betaKc != 0.0) {
00167     FE_Element *elePtr;
00168     FE_EleIter &theEles = theModel->getFEs();    
00169     while((elePtr = theEles()) != 0)     
00170       elePtr->setKc();
00171   }
00172 
00173   if (U == 0) {
00174     cerr << "HHT::newStep() - domainChange() failed or hasn't been called\n";
00175     return -3; 
00176   }
00177 
00178   // set response at t to be that at t+delta t of previous step
00179   (*Ut) = *U;        
00180   (*Utdot) = *Udot;  
00181   (*Utdotdot) = *Udotdot;  
00182     
00183   // set new velocity and accelerations at t + delta t
00184   double a1 = (1.0 - gamma/beta); 
00185 
00186   double a2 = (deltaT)*(1.0 - 0.5*gamma/beta);
00187   // (*Udot) *= a1;  
00188   Udot->addVector(a1,*Utdotdot,a2);
00189 
00190   double a3 = -1.0/(beta*deltaT);
00191   double a4 = 1 - 0.5/beta;
00192   //  (*Udotdot) *= a4;  
00193    Udotdot->addVector(a4,*Utdot,a3);
00194 
00195   (*Ualpha) = *Ut;
00196   (*Udotalpha) = *Utdot;
00197   //  (*Udotalpha) *= (1 - alpha);
00198   Udotalpha->addVector((1-alpha),*Udot, alpha);
00199 
00200   // set the new trial response quantities
00201   theModel->setResponse(*Ualpha,*Udotalpha,*Udotdot);        
00202 
00203   // increment the time and apply the load
00204   double time = theModel->getCurrentDomainTime();
00205   time +=deltaT;
00206   theModel->applyLoadDomain(time);
00207 
00208   theModel->updateDomain();
00209 
00210   return 0;
00211 }
00212 
00213 
00214 int
00215 HHT::revertToLastStep()
00216 {
00217   // set response at t+delta t to be that at t .. for next newStep
00218   if (U != 0) {
00219     (*U) = *Ut;        
00220     (*Udot) = *Utdot;  
00221     (*Udotdot) = *Utdotdot;  
00222   }
00223   return 0;
00224 }
00225 
00226 
00227 
00228 int
00229 HHT::formEleTangent(FE_Element *theEle)
00230 {
00231   theEle->zeroTangent();
00232   if (statusFlag == CURRENT_TANGENT) {
00233     if (rayleighDamping == false) {
00234       theEle->addKtToTang(c1);
00235       theEle->addCtoTang(c2);
00236       theEle->addMtoTang(c3);
00237     } else {
00238       theEle->addKtToTang(c1 + c2*betaK);
00239       theEle->addMtoTang(c3 + c2*alphaM);
00240       theEle->addKiToTang(c2*betaKi);
00241       theEle->addKcToTang(c2*betaKc);
00242     }    
00243   } else if (statusFlag == INITIAL_TANGENT) {
00244     if (rayleighDamping == false) {
00245       theEle->addKiToTang(c1);
00246       theEle->addCtoTang(c2);
00247       theEle->addMtoTang(c3);
00248     } else {
00249       theEle->addKtToTang(c2*betaK);
00250       theEle->addMtoTang(c3 + c2*alphaM);
00251       theEle->addKiToTang(c1 + c2*betaKi);
00252       theEle->addKcToTang(c2*betaKc);
00253     }    
00254   } else if (statusFlag == CURRENT_SECANT) {
00255     if (rayleighDamping == false) {
00256       theEle->addKsToTang(c1);
00257       theEle->addCtoTang(c2);
00258       theEle->addMtoTang(c3);
00259     } else {
00260       theEle->addKsToTang(c1 + c2*betaK);
00261       theEle->addMtoTang(c3 + c2*alphaM);
00262       theEle->addKiToTang(c2*betaKi);
00263       theEle->addKcToTang(c2*betaKc);
00264     }    
00265   }
00266 
00267   return 0;
00268 }    
00269 
00270 
00271 
00272 int
00273 HHT::formNodTangent(DOF_Group *theDof)
00274 {
00275   theDof->zeroTangent();
00276   if (rayleighDamping == false) 
00277       theDof->addMtoTang(c3);
00278   else
00279       theDof->addMtoTang(c3 + c2*alphaM);        
00280   
00281   return(0);
00282 }    
00283 
00284 
00285 
00286 int 
00287 HHT::domainChanged()
00288 {
00289   AnalysisModel *myModel = this->getAnalysisModelPtr();
00290   LinearSOE *theLinSOE = this->getLinearSOEPtr();
00291   const Vector &x = theLinSOE->getX();
00292   int size = x.Size();
00293 
00294   // create the new Vector objects
00295   if (Ut == 0 || Ut->Size() != size) {
00296 
00297     // delete the old
00298     if (Ut != 0)
00299       delete Ut;
00300     if (Utdot != 0)
00301       delete Utdot;
00302     if (Utdotdot != 0)
00303       delete Utdotdot;
00304     if (U != 0)
00305       delete U;
00306     if (Udot != 0)
00307       delete Udot;
00308     if (Udotdot != 0)
00309       delete Udotdot;
00310     if (Ualpha != 0)
00311       delete Ualpha;
00312     if (Udotalpha != 0)
00313       delete Udotalpha;
00314     
00315     // create the new
00316     Ut = new Vector(size);
00317     Utdot = new Vector(size);
00318     Utdotdot = new Vector(size);
00319     U = new Vector(size);
00320     Udot = new Vector(size);
00321     Udotdot = new Vector(size);
00322     Ualpha = new Vector(size);
00323     Udotalpha = new Vector(size);
00324 
00325     // check we obtained the new
00326     if (Ut == 0 || Ut->Size() != size ||
00327  Utdot == 0 || Utdot->Size() != size ||
00328  Utdotdot == 0 || Utdotdot->Size() != size ||
00329  U == 0 || U->Size() != size ||
00330  Udot == 0 || Udot->Size() != size ||
00331  Udotdot == 0 || Udotdot->Size() != size ||
00332  Ualpha == 0 || Ualpha->Size() != size ||
00333  Udotalpha == 0 || Udotalpha->Size() != size) {
00334   
00335       cerr << "HHT::domainChanged - ran out of memory\n";
00336 
00337       // delete the old
00338       if (Ut != 0)
00339  delete Ut;
00340       if (Utdot != 0)
00341  delete Utdot;
00342       if (Utdotdot != 0)
00343  delete Utdotdot;
00344       if (U != 0)
00345  delete U;
00346       if (Udot != 0)
00347  delete Udot;
00348       if (Udotdot != 0)
00349  delete Udotdot;
00350     if (Ualpha != 0)
00351       delete Ualpha;
00352     if (Udotalpha != 0)
00353       delete Udotalpha;
00354 
00355       Ut = 0; Utdot = 0; Utdotdot = 0;
00356       U = 0; Udot = 0; Udotdot = 0; Udotalpha=0; Ualpha =0;
00357       return -1;
00358     }
00359   }        
00360     
00361   // now go through and populate U, Udot and Udotdot by iterating through
00362   // the DOF_Groups and getting the last committed velocity and accel
00363 
00364   DOF_GrpIter &theDOFs = myModel->getDOFs();
00365   DOF_Group *dofPtr;
00366 
00367   while ((dofPtr = theDOFs()) != 0) {
00368     const ID &id = dofPtr->getID();
00369     int idSize = id.Size();
00370 
00371 
00372  int i;
00373     const Vector &disp = dofPtr->getCommittedDisp(); 
00374     for (i=0; i < idSize; i++) {
00375       int loc = id(i);
00376       if (loc >= 0) {
00377   (*U)(loc) = disp(i);  
00378       }
00379     }
00380 
00381     const Vector &vel = dofPtr->getCommittedVel();
00382     for (i=0; i < idSize; i++) {
00383       int loc = id(i);
00384       if (loc >= 0) {
00385   (*Udot)(loc) = vel(i);
00386       }
00387     }
00388 
00389     const Vector &accel = dofPtr->getCommittedAccel(); 
00390     for (i=0; i < idSize; i++) {
00391       int loc = id(i);
00392       if (loc >= 0) {
00393   (*Udotdot)(loc) = accel(i);
00394       }
00395     }
00407   }    
00408 
00409   return 0;
00410 }
00411 
00412 
00413 int
00414 HHT::update(const Vector &deltaU)
00415 {
00416   AnalysisModel *theModel = this->getAnalysisModelPtr();
00417   if (theModel == 0) {
00418     cerr << "WARNING HHT::update() - no AnalysisModel set\n";
00419     return -1;
00420   } 
00421 
00422   // check domainChanged() has been called, i.e. Ut will not be zero
00423   if (Ut == 0) {
00424     cerr << "WARNING HHT::update() - domainChange() failed or not called\n";
00425     return -2;
00426   } 
00427 
00428   // check deltaU is of correct size
00429   if (deltaU.Size() != U->Size()) {
00430     cerr << "WARNING HHT::update() - Vectors of incompatable size ";
00431     cerr << " expecting " << U->Size() << " obtained " << deltaU.Size() << endl;
00432     return -3;
00433   }
00434     
00435   //  determine the response at t+delta t
00436   (*U) += deltaU;
00437   Udot->addVector(1.0,deltaU,c2);
00438   Udotdot->addVector(1.0,deltaU,c3);
00439   Ualpha->addVector(1.0,deltaU, alpha);
00440   Udotalpha->addVector(1.0,deltaU, c2*alpha);
00441   
00442   // update the responses at the DOFs
00443   theModel->setResponse(*Ualpha,*Udotalpha,*Udotdot);        
00444   theModel->updateDomain();
00445       
00446   return 0;
00447 }    
00448 
00449 int
00450 HHT::commit(void)
00451 {
00452   AnalysisModel *theModel = this->getAnalysisModelPtr();
00453   if (theModel == 0) {
00454     cerr << "WARNING HHT::commit() - no AnalysisModel set\n";
00455     return -1;
00456   }   
00457 
00458   // update the responses at the DOFs
00459   theModel->setResponse(*U,*Udot,*Udotdot);        
00460   theModel->updateDomain();
00461 
00462   return theModel->commitDomain();
00463 }
00464 
00465 int
00466 HHT::sendSelf(int cTag, Channel &theChannel)
00467 {
00468     Vector data(8);
00469     data(0) = alpha;
00470     data(1) = beta;
00471     data(2) = gamma;
00472     if (rayleighDamping == true) {
00473  data(3) = 1.0; 
00474  data(4) = alphaM;
00475  data(5) = betaK;
00476  data(6) = betaKi;
00477  data(7) = betaKc;
00478     } else
00479  data(3) = 0.0; 
00480     
00481     if (theChannel.sendVector(this->getDbTag(), cTag, data) < 0) {
00482  cerr << "WARNING HHT::sendSelf() - could not send data\n";
00483  return -1;
00484     } 
00485     return 0;
00486 }
00487 
00488 int
00489 HHT::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
00490 {
00491     Vector data(8);
00492     if (theChannel.recvVector(this->getDbTag(), cTag, data) < 0) {
00493  cerr << "WARNING HHT::recvSelf() - could not receive data\n";
00494  return -1;
00495     }
00496 
00497     alpha = data(0);
00498     beta = data(1);
00499     gamma = data(2);
00500     
00501     if (data(3) == 1.0) {
00502  rayleighDamping = true;
00503  alphaM = data(4);
00504  betaK = data(5);
00505  betaKi = data(6);
00506  betaKc = data(7);
00507     } else
00508  rayleighDamping = false; 
00509     
00510     return 0;
00511 }
00512 
00513 void
00514 HHT::Print(ostream &s, int flag)
00515 {
00516     AnalysisModel *theModel = this->getAnalysisModelPtr();
00517     if (theModel != 0) {
00518  double currentTime = theModel->getCurrentDomainTime();
00519  s << "\t HHT - currentTime: " << currentTime << " alpha: ";
00520  s << alpha << " gamma: " << gamma << "  beta: " << beta << endl;
00521  if (rayleighDamping == true) {
00522      s << "  Rayleigh Damping - alphaM: " << alphaM;
00523      s << "  betaK: " << betaK << endl;     
00524  }
00525     } else 
00526  s << "\t HHT - no associated AnalysisModel\n";
00527 }
00528