NewmarkExplicit.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.2 $
00022 // $Date: 2005/12/21 00:32:57 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/integrator/NewmarkExplicit.cpp,v $
00024 
00025 
00026 // File: ~/analysis/integrator/NewmarkExplicit.cpp
00027 // 
00028 // Written: Andreas Schellenberg (andreas.schellenberg@gmx.net)
00029 // Created: 02/05
00030 // Revision: A
00031 //
00032 // Description: This file contains the implementation of the NewmarkExplicit class.
00033 //
00034 // What: "@(#) NewmarkExplicit.cpp, revA"
00035 
00036 #include <NewmarkExplicit.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 
00048 NewmarkExplicit::NewmarkExplicit()
00049     : TransientIntegrator(INTEGRATOR_TAGS_NewmarkExplicit),
00050     gamma(0), 
00051     alphaM(0.0), betaK(0.0), betaKi(0.0), betaKc(0.0),
00052     c2(0.0), c3(0.0),
00053     Ut(0), Utdot(0), Utdotdot(0), U(0), Udot(0), Udotdot(0)
00054 {
00055     
00056 }
00057 
00058 
00059 NewmarkExplicit::NewmarkExplicit(double _gamma)
00060     : TransientIntegrator(INTEGRATOR_TAGS_NewmarkExplicit),
00061     gamma(_gamma), 
00062     alphaM(0.0), betaK(0.0), betaKi(0.0), betaKc(0.0),
00063     c2(0.0), c3(0.0),
00064     Ut(0), Utdot(0), Utdotdot(0), U(0), Udot(0), Udotdot(0)
00065 {
00066     
00067 }
00068 
00069 
00070 NewmarkExplicit::NewmarkExplicit(double _gamma, 
00071     double _alphaM, double _betaK, double _betaKi , double _betaKc)
00072     : TransientIntegrator(INTEGRATOR_TAGS_NewmarkExplicit),
00073     gamma(_gamma), 
00074     alphaM(_alphaM), betaK(_betaK), betaKi(_betaKi), betaKc(_betaKc),
00075     c2(0.0), c3(0.0),
00076     Ut(0), Utdot(0), Utdotdot(0), U(0), Udot(0), Udotdot(0)
00077 {
00078     
00079 }
00080 
00081 
00082 NewmarkExplicit::~NewmarkExplicit()
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 }
00098 
00099 
00100 int NewmarkExplicit::newStep(double deltaT)
00101 {
00102     updateCount = 0;
00103 
00104     if (gamma == 0)  {
00105         opserr << "NewmarkExplicit::newStep() - error in variable\n";
00106         opserr << "gamma = " << gamma << endln;
00107         return -1;
00108     }
00109     
00110     if (deltaT <= 0.0)  {
00111         opserr << "NewmarkExplicit::newStep() - error in variable\n";
00112         opserr << "dT = " << deltaT << endln;
00113         return -2;
00114     }
00115     
00116     // get a pointer to the AnalysisModel
00117     AnalysisModel *theModel = this->getAnalysisModelPtr();
00118 
00119     // set the constants
00120     c2 = gamma*deltaT;
00121     c3 = 1.0;    
00122     
00123     if (U == 0)  {
00124         opserr << "NewmarkExplicit::newStep() - domainChange() failed or hasn't been called\n";
00125         return -3;      
00126     }
00127     
00128     // set response at t to be that at t+deltaT of previous step
00129     (*Ut) = *U;
00130     (*Utdot) = *Udot;
00131     (*Utdotdot) = *Udotdot;
00132 
00133     // determine new displacements and velocities at time t+deltaT
00134     U->addVector(1.0, *Utdot, deltaT);
00135     double a1 = 0.5*deltaT*deltaT;
00136     U->addVector(1.0, *Utdotdot, a1);
00137     
00138     double a2 = deltaT*(1.0 - gamma);
00139     Udot->addVector(1.0, *Utdotdot, a2);
00140     
00141     // set the trial response quantities for the elements
00142     theModel->setDisp(*U);
00143     theModel->setVel(*Udot);
00144     
00145     // increment the time and apply the load
00146     double time = theModel->getCurrentDomainTime();
00147     time += deltaT;
00148     if (theModel->updateDomain(time, deltaT) < 0)  {
00149         opserr << "NewmarkExplicit::newStep() - failed to update the domain\n";
00150         return -4;
00151     }
00152     
00153     // determine the accelerations at t+deltaT
00154     Udotdot->Zero();
00155         
00156     // set the trial response quantities for the nodes
00157     theModel->setAccel(*Udotdot);
00158     
00159     return 0;
00160 }
00161 
00162 
00163 int NewmarkExplicit::revertToLastStep()
00164 {
00165     // set response at t+deltaT to be that at t .. for next step
00166     if (U != 0)  {
00167         (*U) = *Ut;
00168         (*Udot) = *Utdot;
00169         (*Udotdot) = *Utdotdot;
00170     }
00171 
00172     return 0;
00173 }
00174 
00175 
00176 int NewmarkExplicit::formEleTangent(FE_Element *theEle)
00177 {
00178     theEle->zeroTangent();
00179 
00180     theEle->addCtoTang(c2);
00181     theEle->addMtoTang(c3);
00182     
00183     return 0;
00184 }    
00185 
00186 
00187 int NewmarkExplicit::formNodTangent(DOF_Group *theDof)
00188 {
00189     theDof->zeroTangent();
00190 
00191     theDof->addCtoTang(c2);
00192     theDof->addMtoTang(c3);
00193     
00194     return 0;
00195 }
00196 
00197 
00198 int NewmarkExplicit::domainChanged()
00199 {
00200     AnalysisModel *myModel = this->getAnalysisModelPtr();
00201     LinearSOE *theLinSOE = this->getLinearSOEPtr();
00202     const Vector &x = theLinSOE->getX();
00203     int size = x.Size();
00204     
00205     // if damping factors exist set them in the element & node of the domain
00206     if (alphaM != 0.0 || betaK != 0.0 || betaKi != 0.0 || betaKc != 0.0)
00207         myModel->setRayleighDampingFactors(alphaM, betaK, betaKi, betaKc);
00208     
00209     // create the new Vector objects
00210     if (U == 0 || U->Size() != size)  {
00211         
00212         // delete the old
00213         if (Ut != 0)
00214             delete Ut;
00215         if (Utdot != 0)
00216             delete Utdot;
00217         if (Utdotdot != 0)
00218             delete Utdotdot;
00219         if (U != 0)
00220             delete U;
00221         if (Udot != 0)
00222             delete Udot;
00223         if (Udotdot != 0)
00224             delete Udotdot;
00225         
00226         // create the new
00227         Ut = new Vector(size);
00228         Utdot = new Vector(size);
00229         Utdotdot = new Vector(size);
00230         U = new Vector(size);
00231         Udot = new Vector(size);
00232         Udotdot = new Vector(size);
00233         
00234         // check we obtained the new
00235         if (Ut == 0 || Ut->Size() != size ||
00236             Utdot == 0 || Utdot->Size() != size ||
00237             Utdotdot == 0 || Utdotdot->Size() != size ||
00238             U == 0 || U->Size() != size ||
00239             Udot == 0 || Udot->Size() != size ||
00240             Udotdot == 0 || Udotdot->Size() != size)  {
00241             
00242             opserr << "NewmarkExplicit::domainChanged - ran out of memory\n";
00243             
00244             // delete the old
00245             if (Ut != 0)
00246                 delete Ut;
00247             if (Utdot != 0)
00248                 delete Utdot;
00249             if (Utdotdot != 0)
00250                 delete Utdotdot;
00251             if (U != 0)
00252                 delete U;
00253             if (Udot != 0)
00254                 delete Udot;
00255             if (Udotdot != 0)
00256                 delete Udotdot;
00257             
00258             Ut = 0; Utdot = 0; Utdotdot = 0;
00259             U = 0; Udot = 0; Udotdot = 0;
00260 
00261             return -1;
00262         }
00263     }        
00264     
00265     // now go through and populate U, Udot and Udotdot by iterating through
00266     // the DOF_Groups and getting the last committed velocity and acceleration
00267     DOF_GrpIter &theDOFs = myModel->getDOFs();
00268     DOF_Group *dofPtr;
00269     
00270     while ((dofPtr = theDOFs()) != 0)  {
00271         const ID &id = dofPtr->getID();
00272         int idSize = id.Size();
00273         
00274         int i;
00275         const Vector &disp = dofPtr->getCommittedDisp();        
00276         for (i=0; i < idSize; i++)  {
00277             int loc = id(i);
00278             if (loc >= 0)  {
00279                 (*U)(loc) = disp(i);            
00280             }
00281         }
00282         
00283         const Vector &vel = dofPtr->getCommittedVel();
00284         for (i=0; i < idSize; i++)  {
00285             int loc = id(i);
00286             if (loc >= 0)  {
00287                 (*Udot)(loc) = vel(i);
00288             }
00289         }
00290         
00291         const Vector &accel = dofPtr->getCommittedAccel();      
00292         for (i=0; i < idSize; i++)  {
00293             int loc = id(i);
00294             if (loc >= 0)  {
00295                 (*Udotdot)(loc) = accel(i);
00296             }
00297         }        
00298     }  
00299     
00300     return 0;
00301 }
00302 
00303 
00304 int NewmarkExplicit::update(const Vector &aiPlusOne)
00305 {
00306     updateCount++;
00307     if (updateCount > 1)  {
00308         opserr << "WARNING NewmarkExplicit::update() - called more than once -";
00309         opserr << " NewmarkExplicit integration scheme requires a LINEAR solution algorithm\n";
00310         return -1;
00311     }
00312     
00313     AnalysisModel *theModel = this->getAnalysisModelPtr();
00314     if (theModel == 0)  {
00315         opserr << "WARNING NewmarkExplicit::update() - no AnalysisModel set\n";
00316         return -1;
00317     }   
00318     
00319     // check domainChanged() has been called, i.e. Ut will not be zero
00320     if (Ut == 0)  {
00321         opserr << "WARNING NewmarkExplicit::update() - domainChange() failed or not called\n";
00322         return -2;
00323     }   
00324     
00325     // check aiPlusOne is of correct size
00326     if (aiPlusOne.Size() != U->Size())  {
00327         opserr << "WARNING NewmarkExplicit::update() - Vectors of incompatible size ";
00328         opserr << " expecting " << U->Size() << " obtained " << aiPlusOne.Size() << endln;
00329         return -3;
00330     }
00331     
00332     //  determine the response at t+deltaT
00333     Udot->addVector(1.0, aiPlusOne, c2);
00334 
00335     (*Udotdot) = aiPlusOne;
00336     
00337     // update the response at the DOFs
00338     theModel->setVel(*Udot);
00339     theModel->setAccel(*Udotdot);
00340     
00341     return 0;
00342 }    
00343 
00344 
00345 int NewmarkExplicit::sendSelf(int cTag, Channel &theChannel)
00346 {
00347     Vector data(5);
00348     data(0) = gamma;
00349     data(1) = alphaM;
00350     data(2) = betaK;
00351     data(3) = betaKi;
00352     data(4) = betaKc;
00353     
00354     if (theChannel.sendVector(this->getDbTag(), cTag, data) < 0)  {
00355         opserr << "WARNING NewmarkExplicit::sendSelf() - could not send data\n";
00356         return -1;
00357     }
00358 
00359     return 0;
00360 }
00361 
00362 
00363 int NewmarkExplicit::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
00364 {
00365     Vector data(5);
00366     if (theChannel.recvVector(this->getDbTag(), cTag, data) < 0)  {
00367         opserr << "WARNING NewmarkExplicit::recvSelf() - could not receive data\n";
00368         gamma = 0.5; 
00369         return -1;
00370     }
00371     
00372     gamma  = data(0);
00373     alphaM = data(1);
00374     betaK  = data(2);
00375     betaKi = data(3);
00376     betaKc = data(4);
00377     
00378     return 0;
00379 }
00380 
00381 
00382 void NewmarkExplicit::Print(OPS_Stream &s, int flag)
00383 {
00384     AnalysisModel *theModel = this->getAnalysisModelPtr();
00385     if (theModel != 0)  {
00386         double currentTime = theModel->getCurrentDomainTime();
00387         s << "\t NewmarkExplicit - currentTime: " << currentTime << endln;
00388         s << "  gamma: " << gamma << endln;
00389         s << "  c2: " << c2 << " c3: " << c3 << endln;
00390         s << "  Rayleigh Damping - alphaM: " << alphaM;
00391         s << "  betaK: " << betaK << "  betaKi: " << betaKi << endln;       
00392     } else 
00393         s << "\t NewmarkExplicit - no associated AnalysisModel\n";
00394 }
00395 
00396 
00397 
00398 
00399 

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