00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
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
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
00117 AnalysisModel *theModel = this->getAnalysisModelPtr();
00118
00119
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
00129 (*Ut) = *U;
00130 (*Utdot) = *Udot;
00131 (*Utdotdot) = *Udotdot;
00132
00133
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
00142 theModel->setDisp(*U);
00143 theModel->setVel(*Udot);
00144
00145
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
00154 Udotdot->Zero();
00155
00156
00157 theModel->setAccel(*Udotdot);
00158
00159 return 0;
00160 }
00161
00162
00163 int NewmarkExplicit::revertToLastStep()
00164 {
00165
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
00206 if (alphaM != 0.0 || betaK != 0.0 || betaKi != 0.0 || betaKc != 0.0)
00207 myModel->setRayleighDampingFactors(alphaM, betaK, betaKi, betaKc);
00208
00209
00210 if (U == 0 || U->Size() != size) {
00211
00212
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
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
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
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
00266
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
00320 if (Ut == 0) {
00321 opserr << "WARNING NewmarkExplicit::update() - domainChange() failed or not called\n";
00322 return -2;
00323 }
00324
00325
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
00333 Udot->addVector(1.0, aiPlusOne, c2);
00334
00335 (*Udotdot) = aiPlusOne;
00336
00337
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