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

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