Newmark1.cpp

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.8 $
00022 // $Date: 2003/02/14 23:00:49 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/integrator/Newmark1.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/analysis/integrator/Newmark1.C
00027 // 
00028 // Written: fmk 
00029 // Created: 11/98
00030 // Revision: A
00031 //
00032 // Description: This file contains the implementation of the Newmark1 class.
00033 //
00034 // What: "@(#) Newmark1.C, revA"
00035 
00036 #include <Newmark1.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 Newmark1::Newmark1()
00048 :TransientIntegrator(INTEGRATOR_TAGS_Newmark1),
00049  gamma(0), beta(0), 
00050  alphaM(0.0), betaK(0.0), betaKi(0.0),
00051  c1(0.0), c2(0.0), c3(0.0), c4(0.0),
00052  Up(0), Updot(0),  U(0), Udot(0), Udotdot(0)
00053 {
00054     
00055 }
00056 
00057 Newmark1::Newmark1(double theGamma, double theBeta, bool dispFlag)
00058 :TransientIntegrator(INTEGRATOR_TAGS_Newmark1),
00059  gamma(theGamma), beta(theBeta), 
00060  alphaM(0.0), betaK(0.0), betaKi(0.0), betaKc(0.0),
00061  c1(0.0), c2(0.0), c3(0.0), c4(0.0),
00062  Up(0), Updot(0),  U(0), Udot(0), Udotdot(0)
00063 {
00064 
00065 }
00066 
00067 Newmark1::Newmark1(double theGamma, double theBeta, 
00068                    double alpham, double betak, double betaki , double betakc)
00069 :TransientIntegrator(INTEGRATOR_TAGS_Newmark1),
00070  gamma(theGamma), beta(theBeta), 
00071  alphaM(alpham), betaK(betak), betaKi(betaki), betaKc(betakc),
00072  c1(0.0), c2(0.0), c3(0.0), c4(0.0),
00073  Up(0), Updot(0),  U(0), Udot(0), Udotdot(0)
00074 {
00075 
00076 }
00077 
00078 Newmark1::~Newmark1()
00079 {
00080   // clean up the memory created
00081   if (Up != 0)
00082     delete Up;
00083   if (Updot != 0)
00084     delete Updot;
00085   if (U != 0)
00086     delete U;
00087   if (Udot != 0)
00088     delete Udot;
00089   if (Udotdot != 0)
00090     delete Udotdot;
00091 }
00092 
00093 
00094 int
00095 Newmark1::initialize(void)
00096 {
00097   return 0;
00098 }
00099 
00100 int
00101 Newmark1::newStep(double deltaT)
00102 {
00103 
00104   if (beta == 0 || gamma == 0 ) {
00105     opserr << "Newton::newStep() - error in variable\n";
00106     opserr << "gamma = " << gamma << " beta= " << beta << endln;
00107     return -1;
00108   }
00109 
00110   if (deltaT <= 0.0) {
00111     opserr << "Newmark1::newStep() - error in variable\n";
00112     opserr << "dT = " << deltaT << endln;
00113     return -2;  
00114   }
00115 
00116   // set the constants
00117   c1 = 1.0;
00118   c2 = gamma/(beta*deltaT);
00119   c3 = 1.0/(beta*deltaT*deltaT);
00120   c4 = gamma*deltaT;
00121 
00122 
00123   // set the new trial response quantities
00124   AnalysisModel *theModel = this->getAnalysisModelPtr();
00125 
00126   if (U == 0) {
00127     opserr << "Newton::newStep() - domainChange() failed or hasn't been called\n";
00128     return -3;  
00129   }
00130 
00131   // determine predicted quantities at time t + delta t
00132   U->addVector(1.0, *Udot, deltaT);
00133   double a1 = deltaT * deltaT * (.5 - beta);
00134   U->addVector(1.0, *Udotdot, a1);
00135 
00136   double a2 = deltaT * (1 - gamma);
00137   Udot->addVector(1.0, *Udotdot,a2);
00138 
00139   Udotdot->Zero();
00140 
00141   *Up = *U;
00142   *Updot = *Udot;
00143 
00144   theModel->setResponse(*U,*Udot,*Udotdot);        
00145 
00146   // increment the time and apply the load
00147   double time = theModel->getCurrentDomainTime();
00148   time +=deltaT;
00149 
00150   if (theModel->updateDomain(time, deltaT) < 0) {
00151     opserr << "Newmark1::newStep() - failed to update the domain\n";
00152     return -4;
00153   }
00154   
00155   return 0;
00156 }
00157 
00158 int
00159 Newmark1::revertToLastStep() {
00160   return this->domainChanged();
00161 }
00162 
00163 int
00164 Newmark1::formEleTangent(FE_Element *theEle)
00165 {
00166   theEle->zeroTangent();
00167   if (statusFlag == CURRENT_TANGENT) {
00168     theEle->addKtToTang(c1);
00169     theEle->addCtoTang(c2);
00170     theEle->addMtoTang(c3);
00171   } else if (statusFlag == INITIAL_TANGENT) {
00172     theEle->addKiToTang(c1);
00173     theEle->addCtoTang(c2);
00174     theEle->addMtoTang(c3);
00175   }
00176 
00177   return 0;
00178 }    
00179 
00180 
00181 int
00182 Newmark1::formNodTangent(DOF_Group *theDof)
00183 {
00184   theDof->zeroTangent();
00185   theDof->addMtoTang(c3);
00186   theDof->addCtoTang(c2);      
00187 
00188   return(0);
00189 }    
00190 
00191 
00192 
00193 int 
00194 Newmark1::domainChanged()
00195 {
00196   AnalysisModel *myModel = this->getAnalysisModelPtr();
00197   LinearSOE *theLinSOE = this->getLinearSOEPtr();
00198   const Vector &x = theLinSOE->getX();
00199   int size = x.Size();
00200 
00201   // if damping factors exist set them in the ele & node of the domain
00202   if (alphaM != 0.0 || betaK != 0.0 || betaKi != 0.0 || betaKc != 0.0)
00203     myModel->setRayleighDampingFactors(alphaM, betaK, betaKi, betaKc);
00204   
00205   // create the new Vector objects
00206   if (U == 0 || U ->Size() != size) {
00207 
00208     // delete the old
00209     if (Up != 0)
00210       delete Up;
00211     if (Updot != 0)
00212       delete Updot;
00213     if (U != 0)
00214       delete U;
00215     if (Udot != 0)
00216       delete Udot;
00217     if (Udotdot != 0)
00218       delete Udotdot;
00219     
00220     // create the new
00221     Up = new Vector(size);
00222     Updot = new Vector(size);
00223     U = new Vector(size);
00224     Udot = new Vector(size);
00225     Udotdot = new Vector(size);
00226 
00227     // cheack we obtained the new
00228     if (Up == 0 || Up->Size() != size ||
00229         Updot == 0 || Updot->Size() != size ||
00230         U == 0 || U->Size() != size ||
00231         Udot == 0 || Udot->Size() != size ||
00232         Udotdot == 0 || Udotdot->Size() != size) {
00233       
00234       opserr << "Newmark1::domainChanged - ran out of memory\n";
00235 
00236       // delete the old
00237       if (Up != 0)
00238         delete Up;
00239       if (Updot != 0)
00240         delete Updot;
00241       if (U != 0)
00242         delete U;
00243       if (Udot != 0)
00244         delete Udot;
00245       if (Udotdot != 0)
00246         delete Udotdot;
00247 
00248       Up = 0; Updot = 0; 
00249       U = 0; Udot = 0; Udotdot = 0;
00250       return -1;
00251     }
00252   }        
00253     
00254   // now go through and populate U, Udot and Udotdot by iterating through
00255   // the DOF_Groups and getting the last committed velocity and accel
00256 
00257   DOF_GrpIter &theDOFs = myModel->getDOFs();
00258   DOF_Group *dofPtr;
00259     
00260   while ((dofPtr = theDOFs()) != 0) {
00261     const ID &id = dofPtr->getID();
00262     int idSize = id.Size();
00263 
00264 
00265         int i;
00266     const Vector &disp = dofPtr->getCommittedDisp();    
00267     for (i=0; i < idSize; i++) {
00268       int loc = id(i);
00269       if (loc >= 0) {
00270         (*U)(loc) = disp(i);            
00271       }
00272     }
00273 
00274     const Vector &vel = dofPtr->getCommittedVel();
00275     for (i=0; i < idSize; i++) {
00276       int loc = id(i);
00277       if (loc >= 0) {
00278         (*Udot)(loc) = vel(i);
00279       }
00280     }
00281 
00282     const Vector &accel = dofPtr->getCommittedAccel();  
00283     for (i=0; i < idSize; i++) {
00284       int loc = id(i);
00285       if (loc >= 0) {
00286         (*Udotdot)(loc) = accel(i);
00287       }
00288     }
00289 
00301   }    
00302   return 0;
00303 }
00304 
00305 
00306 int
00307 Newmark1::update(const Vector &deltaU)
00308 {
00309   AnalysisModel *theModel = this->getAnalysisModelPtr();
00310   if (theModel == 0) {
00311     opserr << "WARNING Newmark1::update() - no AnalysisModel set\n";
00312     return -1;
00313   }     
00314 
00315   // check domainChanged() has been called, i.e. Ut will not be zero
00316   if (U == 0) {
00317     opserr << "WARNING Newmark1::update() - domainChange() failed or not called\n";
00318     return -2;
00319   }     
00320 
00321   // check deltaU is of correct size
00322   if (deltaU.Size() != U->Size()) {
00323     opserr << "WARNING Newmark1::update() - Vectors of incompatable size ";
00324     opserr << " expecting " << U->Size() << " obtained " << deltaU.Size() << endln;
00325     return -3;
00326   }
00327     
00328   //  determine the response at t+delta t
00329   (*U) += deltaU;
00330 
00331   (*Udotdot) = (*U);
00332   (*Udotdot) -= (*Up);
00333   (*Udotdot) *= c3;
00334 
00335   (*Udot) = (*Updot);
00336   Udot->addVector(1.0, *Udotdot,c4);
00337   
00338   // update the responses at the DOFs
00339   theModel->setResponse(*U,*Udot,*Udotdot);        
00340   if (theModel->updateDomain() < 0) {
00341     opserr << "Newmark1::newStep() - failed to update the domain\n";
00342     return -4;
00343   }
00344                     
00345   return 0;
00346 }    
00347 
00348 int
00349 Newmark1::sendSelf(int cTag, Channel &theChannel)
00350 {
00351     Vector data(7);
00352     data(0) = gamma;
00353     data(1) = beta;
00354     data(2) = 1.0;      
00355     data(3) = alphaM;
00356     data(4) = betaK;
00357     data(5) = betaKi;
00358     data(6) = betaKc;
00359     
00360     if (theChannel.sendVector(this->getDbTag(), cTag, data) < 0) {
00361         opserr << "WARNING Newmark1::sendSelf() - could not send data\n";
00362         return -1;
00363     }   
00364     return 0;
00365 }
00366 
00367 int
00368 Newmark1::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
00369 {
00370     Vector data(7);
00371     if (theChannel.recvVector(this->getDbTag(), cTag, data) < 0) {
00372         opserr << "WARNING Newmark1::recvSelf() - could not receive data\n";
00373         gamma = 0.5; beta = 0.25; 
00374         return -1;
00375     }
00376     
00377     gamma = data(0);
00378     beta = data(1);
00379     alphaM = data(3);
00380     betaK = data(4);
00381     betaKi = data(5);
00382     betaKc = data(6);
00383       
00384     return 0;
00385     
00386 }
00387 
00388 void
00389 Newmark1::Print(OPS_Stream &s, int flag)
00390 {
00391     AnalysisModel *theModel = this->getAnalysisModelPtr();
00392     if (theModel != 0) {
00393         double currentTime = theModel->getCurrentDomainTime();
00394         s << "\t Newmark1 - currentTime: " << currentTime;
00395         s << "  gamma: " << gamma << "  beta: " << beta << endln;
00396         s << " c1: " << c1 << " c2: " << c2 << " c3: " << c3 << endln;
00397         s << "  Rayleigh Damping - alphaM: " << alphaM;
00398         s << "  betaK: " << betaK << "  betaKi: " << betaKi << endln;       
00399     } else 
00400         s << "\t Newmark1 - no associated AnalysisModel\n";
00401 }
00402 
00403 
00404 
00405
Generated on Mon Oct 23 15:04:59 2006 for OpenSees by 1.5.0
opensees-support @ berkeley.edu
Supported by the National Science Foundation