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

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