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 #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
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
00098 AnalysisModel *theModel = this->getAnalysisModelPtr();
00099
00100
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
00110 double time = theModel->getCurrentDomainTime();
00111 theModel->applyLoadDomain(time);
00112
00113
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
00120 theModel->setVel(*Utdot);
00121 theModel->setAccel(*Utdotdot);
00122
00123
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
00161 if (alphaM != 0.0 || betaK != 0.0 || betaKi != 0.0 || betaKc != 0.0)
00162 myModel->setRayleighDampingFactors(alphaM, betaK, betaKi, betaKc);
00163
00164
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
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
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
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
00221
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
00278 if (Ut == 0) {
00279 opserr << "WARNING CentralDifference::update() - domainChange() failed or not called\n";
00280 return -2;
00281 }
00282
00283
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
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
00301 (*Utm1) = *Ut;
00302 (*Ut) = U;
00303
00304
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
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 }