OpenSees » SRC » analysis » integrator

Collocation.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/Collocation.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 Collocation.
00030 //
00031 // What: "@(#) Collocation.cpp, revA"
00032 
00033 #include <Collocation.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 #include <math.h>
00044 
00045 
00046 Collocation::Collocation()
00047     : TransientIntegrator(INTEGRATOR_TAGS_Collocation),
00048     theta(1.0), gamma(0.0), beta(0.0), deltaT(0.0),
00049     alphaM(0.0), betaK(0.0), betaKi(0.0), betaKc(0.0),
00050     c1(0.0), c2(0.0), c3(0.0), 
00051     Ut(0), Utdot(0), Utdotdot(0), U(0), Udot(0), Udotdot(0)
00052 {
00053     
00054 }
00055 
00056 
00057 Collocation::Collocation(double _theta)
00058     : TransientIntegrator(INTEGRATOR_TAGS_Collocation),
00059     theta(_theta), gamma(0.5), beta(0.0), deltaT(0.0),
00060     alphaM(0.0), betaK(0.0), betaKi(0.0), betaKc(0.0),
00061     c1(0.0), c2(0.0), c3(0.0), 
00062     Ut(0), Utdot(0), Utdotdot(0), U(0), Udot(0), Udotdot(0)
00063 {
00064     beta = -6.018722044382699e+002 * pow(theta,9) +
00065             6.618777151634235e+003 * pow(theta,8) +
00066            -3.231561059595987e+004 * pow(theta,7) +
00067             9.195359004558867e+004 * pow(theta,6) +
00068            -1.680788908312227e+005 * pow(theta,5) +
00069             2.047005794710718e+005 * pow(theta,4) +
00070            -1.661421563528177e+005 * pow(theta,3) +
00071             8.667950092619179e+004 * pow(theta,2) +
00072            -2.638652989051994e+004 * theta +
00073             3.572862280471971e+003;
00074 }
00075 
00076 
00077 Collocation::Collocation(double _theta,
00078     double _alphaM, double _betaK, double _betaKi, double _betaKc)
00079     : TransientIntegrator(INTEGRATOR_TAGS_Collocation),
00080     theta(_theta), gamma(0.5), beta(0.0), deltaT(0.0),
00081     alphaM(_alphaM), betaK(_betaK), betaKi(_betaKi), betaKc(_betaKc),
00082     c1(0.0), c2(0.0), c3(0.0), 
00083     Ut(0), Utdot(0), Utdotdot(0), U(0), Udot(0), Udotdot(0)
00084 {
00085     beta = -6.018722044382699e+002 * pow(theta,9) +
00086             6.618777151634235e+003 * pow(theta,8) +
00087            -3.231561059595987e+004 * pow(theta,7) +
00088             9.195359004558867e+004 * pow(theta,6) +
00089            -1.680788908312227e+005 * pow(theta,5) +
00090             2.047005794710718e+005 * pow(theta,4) +
00091            -1.661421563528177e+005 * pow(theta,3) +
00092             8.667950092619179e+004 * pow(theta,2) +
00093            -2.638652989051994e+004 * theta +
00094             3.572862280471971e+003;
00095 }
00096 
00097 
00098 Collocation::Collocation(double _theta, double _beta, double _gamma)
00099     : TransientIntegrator(INTEGRATOR_TAGS_Collocation),
00100     theta(_theta), gamma(_gamma), beta(_beta), deltaT(0.0),
00101     alphaM(0.0), betaK(0.0), betaKi(0.0), betaKc(0.0),
00102     c1(0.0), c2(0.0), c3(0.0), 
00103     Ut(0), Utdot(0), Utdotdot(0), U(0), Udot(0), Udotdot(0)
00104 {
00105     
00106 }
00107 
00108 
00109 Collocation::Collocation(double _theta, double _beta, double _gamma,
00110     double _alphaM, double _betaK, double _betaKi, double _betaKc)
00111     : TransientIntegrator(INTEGRATOR_TAGS_Collocation),
00112     theta(_theta), gamma(_gamma), beta(_beta), deltaT(0.0),
00113     alphaM(_alphaM), betaK(_betaK), betaKi(_betaKi), betaKc(_betaKc),
00114     c1(0.0), c2(0.0), c3(0.0), 
00115     Ut(0), Utdot(0), Utdotdot(0), U(0), Udot(0), Udotdot(0)
00116 {
00117     
00118 }
00119 
00120 
00121 Collocation::~Collocation()
00122 {
00123     // clean up the memory created
00124     if (Ut != 0)
00125         delete Ut;
00126     if (Utdot != 0)
00127         delete Utdot;
00128     if (Utdotdot != 0)
00129         delete Utdotdot;
00130     if (U != 0)
00131         delete U;
00132     if (Udot != 0)
00133         delete Udot;
00134     if (Udotdot != 0)
00135         delete Udotdot;
00136 }
00137 
00138 
00139 int Collocation::newStep(double _deltaT)
00140 {
00141     deltaT = _deltaT;
00142     if (theta <= 0.0 )  {
00143         opserr << "Collocation::newStep() - error in variable\n";
00144         opserr << "theta: " << theta << " <= 0.0\n";
00145         return -1;
00146     }
00147     
00148     if (deltaT <= 0.0)  {
00149         opserr << "Collocation::newStep() - error in variable\n";
00150         opserr << "dT = " << deltaT << endln;
00151         return -2;
00152     }
00153     
00154     // get a pointer to the AnalysisModel
00155     AnalysisModel *theModel = this->getAnalysisModelPtr();
00156 
00157     // set the constants
00158     c1 = 1.0;
00159     c2 = gamma/(beta*theta*deltaT);
00160     c3 = 1.0/(beta*theta*theta*deltaT*deltaT);
00161     
00162     if (U == 0)  {
00163         opserr << "Collocation::newStep() - domainChange() failed or hasn't been called\n"; 
00164         return -3;      
00165     }
00166     
00167     // set response at t to be that at t+deltaT of previous step
00168     (*Ut) = *U;
00169     (*Utdot) = *Udot;
00170     (*Utdotdot) = *Udotdot;
00171 
00172     // increment the time to t+theta*deltaT and apply the load
00173     double time = theModel->getCurrentDomainTime();
00174     time += theta*deltaT;
00175 //    theModel->applyLoadDomain(time);
00176     if (theModel->updateDomain(time, deltaT) < 0)  {
00177         opserr << "Collocation::newStep() - failed to update the domain\n";
00178         return -4;
00179     }
00180     
00181     // determine new velocities and accelerations at t+theta*deltaT
00182     double a1 = (1.0 - gamma/beta); 
00183     double a2 = theta*deltaT*(1.0 - 0.5*gamma/beta);
00184     Udot->addVector(a1, *Utdotdot, a2);
00185       
00186     double a3 = -1.0/(beta*theta*deltaT);
00187     double a4 = 1.0 - 0.5/beta;
00188     Udotdot->addVector(a4, *Utdot, a3);
00189     
00190     // set the trial response quantities for the nodes
00191     theModel->setVel(*Udot);
00192     theModel->setAccel(*Udotdot);
00193     
00194     return 0;
00195 }
00196 
00197 
00198 int Collocation::revertToLastStep()
00199 {
00200     // set response at t+deltaT to be that at t .. for next step
00201     if (U != 0)  {
00202         (*U) = *Ut;
00203         (*Udot) = *Utdot;
00204         (*Udotdot) = *Utdotdot;
00205     }
00206 
00207     return 0;
00208 }
00209 
00210 
00211 int Collocation::formEleTangent(FE_Element *theEle)
00212 {
00213     theEle->zeroTangent();
00214     if (statusFlag == CURRENT_TANGENT)  {
00215         theEle->addKtToTang(c1);
00216         theEle->addCtoTang(c2);
00217         theEle->addMtoTang(c3);
00218     } else if (statusFlag == INITIAL_TANGENT)  {
00219         theEle->addKiToTang(c1);
00220         theEle->addCtoTang(c2);
00221         theEle->addMtoTang(c3);
00222     }
00223     
00224     return 0;
00225 }
00226 
00227 
00228 int Collocation::formNodTangent(DOF_Group *theDof)
00229 {
00230     theDof->zeroTangent();
00231 
00232     theDof->addCtoTang(c2);
00233     theDof->addMtoTang(c3);
00234     
00235     return 0;
00236 }    
00237 
00238 
00239 int Collocation::domainChanged()
00240 {
00241     AnalysisModel *myModel = this->getAnalysisModelPtr();
00242     LinearSOE *theLinSOE = this->getLinearSOEPtr();
00243     const Vector &x = theLinSOE->getX();
00244     int size = x.Size();
00245     
00246     // if damping factors exist set them in the ele & node of the domain
00247     if (alphaM != 0.0 || betaK != 0.0 || betaKi != 0.0 || betaKc != 0.0)
00248         myModel->setRayleighDampingFactors(alphaM, betaK, betaKi, betaKc);
00249     
00250     // create the new Vector objects
00251     if (Ut == 0 || Ut->Size() != size)  {
00252         
00253         // delete the old
00254         if (Ut != 0)
00255             delete Ut;
00256         if (Utdot != 0)
00257             delete Utdot;
00258         if (Utdotdot != 0)
00259             delete Utdotdot;
00260         if (U != 0)
00261             delete U;
00262         if (Udot != 0)
00263             delete Udot;
00264         if (Udotdot != 0)
00265             delete Udotdot;
00266         
00267         // create the new
00268         Ut = new Vector(size);
00269         Utdot = new Vector(size);
00270         Utdotdot = new Vector(size);
00271         U = new Vector(size);
00272         Udot = new Vector(size);
00273         Udotdot = new Vector(size);
00274         
00275         // check we obtained the new
00276         if (Ut == 0 || Ut->Size() != size ||
00277             Utdot == 0 || Utdot->Size() != size ||
00278             Utdotdot == 0 || Utdotdot->Size() != size ||
00279             U == 0 || U->Size() != size ||
00280             Udot == 0 || Udot->Size() != size ||
00281             Udotdot == 0 || Udotdot->Size() != size)  {
00282             
00283             opserr << "Collocation::domainChanged - ran out of memory\n";
00284             
00285             // delete the old
00286             if (Ut != 0)
00287                 delete Ut;
00288             if (Utdot != 0)
00289                 delete Utdot;
00290             if (Utdotdot != 0)
00291                 delete Utdotdot;
00292             if (U != 0)
00293                 delete U;
00294             if (Udot != 0)
00295                 delete Udot;
00296             if (Udotdot != 0)
00297                 delete Udotdot;
00298             
00299             Ut = 0; Utdot = 0; Utdotdot = 0;
00300             U = 0; Udot = 0; Udotdot = 0;
00301 
00302             return -1;
00303         }
00304     }        
00305     
00306     // now go through and populate U, Udot and Udotdot by iterating through
00307     // the DOF_Groups and getting the last committed velocity and accel
00308     DOF_GrpIter &theDOFs = myModel->getDOFs();
00309     DOF_Group *dofPtr;
00310     
00311     while ((dofPtr = theDOFs()) != 0)  {
00312         const ID &id = dofPtr->getID();
00313         int idSize = id.Size();
00314         
00315         int i;
00316         const Vector &disp = dofPtr->getCommittedDisp();        
00317         for (i=0; i < idSize; i++)  {
00318             int loc = id(i);
00319             if (loc >= 0)  {
00320                 (*U)(loc) = disp(i);            
00321             }
00322         }
00323         
00324         const Vector &vel = dofPtr->getCommittedVel();
00325         for (i=0; i < idSize; i++)  {
00326             int loc = id(i);
00327             if (loc >= 0)  {
00328                 (*Udot)(loc) = vel(i);
00329             }
00330         }
00331         
00332         const Vector &accel = dofPtr->getCommittedAccel();      
00333         for (i=0; i < idSize; i++)  {
00334             int loc = id(i);
00335             if (loc >= 0)  {
00336                 (*Udotdot)(loc) = accel(i);
00337             }
00338         }        
00339     }    
00340     
00341     return 0;
00342 }
00343 
00344 
00345 int Collocation::update(const Vector &deltaU)
00346 {
00347     AnalysisModel *theModel = this->getAnalysisModelPtr();
00348     if (theModel == 0)  {
00349         opserr << "WARNING Collocation::update() - no AnalysisModel set\n";
00350         return -1;
00351     }   
00352     
00353     // check domainChanged() has been called, i.e. Ut will not be zero
00354     if (Ut == 0)  {
00355         opserr << "WARNING Collocation::update() - domainChange() failed or not called\n";
00356         return -2;
00357     }   
00358     
00359     // check deltaU is of correct size
00360     if (deltaU.Size() != U->Size())  {
00361         opserr << "WARNING Collocation::update() - Vectors of incompatible size ";
00362         opserr << " expecting " << U->Size() << " obtained " << deltaU.Size() << endln;
00363         return -3;
00364     }
00365     
00366     // determine the response at t+theta*deltaT
00367     (*U) += deltaU;
00368 
00369     Udot->addVector(1.0, deltaU, c2);
00370     
00371     Udotdot->addVector(1.0, deltaU, c3);
00372     
00373     // update the response at the DOFs
00374     theModel->setResponse(*U,*Udot,*Udotdot);        
00375     if (theModel->updateDomain() < 0)  {
00376         opserr << "Collocation::update() - failed to update the domain\n";
00377         return -4;
00378     }
00379     
00380     return 0;
00381 }    
00382 
00383 
00384 int Collocation::commit(void)
00385 {
00386     
00387     AnalysisModel *theModel = this->getAnalysisModelPtr();
00388     if (theModel == 0)  {
00389         opserr << "WARNING Collocation::commit() - no AnalysisModel set\n";
00390         return -1;
00391     }     
00392         
00393     // determine response quantities at t+deltaT
00394     Udotdot->addVector(1.0/theta, *Utdotdot, (theta-1.0)/theta);
00395     
00396     (*Udot) = *Utdot;
00397     double a1 = deltaT*(1.0 - gamma);
00398     double a2 = deltaT*gamma;
00399     Udot->addVector(1.0, *Utdotdot, a1);
00400     Udot->addVector(1.0, *Udotdot, a2);
00401     
00402     (*U) = *Ut;
00403     U->addVector(1.0, *Utdot, deltaT);
00404     double a3 = deltaT*deltaT*(0.5 - beta);
00405     double a4 = deltaT*deltaT*beta;
00406     U->addVector(1.0, *Utdotdot, a3);
00407     U->addVector(1.0, *Udotdot, a4);
00408 
00409     // update the response at the DOFs
00410     theModel->setResponse(*U,*Udot,*Udotdot);        
00411 //    if (theModel->updateDomain() < 0)  {
00412 //        opserr << "Collocation::commit() - failed to update the domain\n";
00413 //        return -4;
00414 //    }
00415     
00416     // set the time to be t+delta t
00417     double time = theModel->getCurrentDomainTime();
00418     time += (1.0-theta)*deltaT;
00419     theModel->setCurrentDomainTime(time);
00420     
00421     return theModel->commitDomain();
00422 }
00423 
00424 
00425 int Collocation::sendSelf(int cTag, Channel &theChannel)
00426 {
00427     static Vector data(5);
00428     data(0) = theta;
00429     data(1) = alphaM;
00430     data(2) = betaK;
00431     data(3) = betaKi;
00432     data(4) = betaKc;
00433     
00434     if (theChannel.sendVector(this->getDbTag(), cTag, data) < 0)  {
00435         opserr << "WARNING Collocation::sendSelf() - failed to send the data\n";
00436         return -1;
00437     }
00438 
00439     return 0;
00440 }
00441 
00442 
00443 int Collocation::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
00444 {
00445     Vector data(5);
00446     if (theChannel.recvVector(this->getDbTag(), cTag, data) < 0)  {
00447         opserr << "WARNING Collocation::recvSelf() - could not receive data\n";
00448         return -1;
00449     }
00450     
00451     theta  = data(0);
00452     alphaM = data(1);
00453     betaK  = data(2);
00454     betaKi = data(3);
00455     betaKc = data(4);
00456     
00457     return 0;
00458 }
00459 
00460 
00461 void Collocation::Print(OPS_Stream &s, int flag)
00462 {
00463     AnalysisModel *theModel = this->getAnalysisModelPtr();
00464     if (theModel != 0)  {
00465         double currentTime = theModel->getCurrentDomainTime();
00466         s << "\t Collocation - currentTime: " << currentTime << endln;
00467         s << "  theta: " << theta << endln;
00468         s << "  Rayleigh Damping - alphaM: " << alphaM;
00469         s << "  betaK: " << betaK << "  betaKi: " << betaKi << endln;       
00470     } else 
00471         s << "\t Collocation - no associated AnalysisModel\n";
00472 }

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