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