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 #include <CentralDifferenceAlternative.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
00045 CentralDifferenceAlternative::CentralDifferenceAlternative()
00046 :TransientIntegrator(INTEGRATOR_TAGS_CentralDifferenceAlternative),
00047 updateCount(0), Ut(0), Utp1(0), Udot(0), deltaT(0)
00048 {
00049
00050 }
00051
00052 CentralDifferenceAlternative::~CentralDifferenceAlternative()
00053 {
00054
00055 if (Ut != 0)
00056 delete Ut;
00057 if (Utp1 != 0)
00058 delete Utp1;
00059 if (Udot != 0)
00060 delete Udot;
00061 }
00062
00063 int
00064 CentralDifferenceAlternative::newStep(double _deltaT)
00065 {
00066 updateCount = 0;
00067
00068 deltaT = _deltaT;
00069
00070 if (deltaT <= 0.0) {
00071 opserr << "CentralDifference::newStep() - error in variable\n";
00072 opserr << "dT = " << deltaT << endln;
00073 return -2;
00074 }
00075
00076 AnalysisModel *theModel = this->getAnalysisModelPtr();
00077 double time = theModel->getCurrentDomainTime();
00078 theModel->applyLoadDomain(time);
00079
00080 return 0;
00081 }
00082
00083 int
00084 CentralDifferenceAlternative::formEleTangent(FE_Element *theEle)
00085 {
00086 theEle->zeroTangent();
00087 theEle->addMtoTang();
00088
00089 return 0;
00090 }
00091
00092 int
00093 CentralDifferenceAlternative::formNodTangent(DOF_Group *theDof)
00094 {
00095 theDof->zeroTangent();
00096 theDof->addMtoTang();
00097 return(0);
00098 }
00099
00100 int
00101 CentralDifferenceAlternative::domainChanged()
00102 {
00103 AnalysisModel *myModel = this->getAnalysisModelPtr();
00104 LinearSOE *theLinSOE = this->getLinearSOEPtr();
00105 const Vector &x = theLinSOE->getX();
00106 int size = x.Size();
00107
00108
00109 if (Ut == 0 || Ut->Size() != size) {
00110
00111
00112 if (Ut != 0)
00113 delete Ut;
00114 if (Utp1 != 0)
00115 delete Utp1;
00116 if (Udot != 0)
00117 delete Udot;
00118
00119
00120 Ut = new Vector(size);
00121 Utp1 = new Vector(size);
00122 Udot = new Vector(size);
00123
00124
00125 if (Ut == 0 || Ut->Size() != size ||
00126 Utp1 == 0 || Utp1->Size() != size ||
00127 Udot == 0 || Udot->Size() != size) {
00128
00129 opserr << "CentralDifferenceAlternative::domainChanged - ran out of memory\n";
00130
00131
00132 if (Ut != 0)
00133 delete Ut;
00134 if (Utp1 != 0)
00135 delete Utp1;
00136 if (Udot != 0)
00137 delete Udot;
00138
00139 Ut = 0; Utp1 = 0; Udot = 0;
00140 return -1;
00141 }
00142 }
00143
00144
00145
00146
00147 DOF_GrpIter &theDOFs = myModel->getDOFs();
00148 DOF_Group *dofPtr;
00149
00150 while ((dofPtr = theDOFs()) != 0) {
00151 const ID &id = dofPtr->getID();
00152 int idSize = id.Size();
00153 int i;
00154 const Vector &disp = dofPtr->getCommittedDisp();
00155 for (i=0; i < idSize; i++) {
00156 int loc = id(i);
00157 if (loc >= 0) {
00158 (*Ut)(loc) = disp(i);
00159 }
00160 }
00161
00162 const Vector &vel = dofPtr->getCommittedVel();
00163 for (i=0; i < idSize; i++) {
00164 int loc = id(i);
00165 if (loc >= 0) {
00166 (*Udot)(loc) = vel(i);
00167 }
00168 }
00169 }
00170
00171 return 0;
00172 }
00173
00174
00175 int
00176 CentralDifferenceAlternative::update(const Vector &X)
00177 {
00178 updateCount++;
00179 if (updateCount > 1) {
00180 opserr << "ERROR CentralDifferenceAlternative::update() - called more than once -";
00181 opserr << " Central Difference integraion schemes require a LINEAR solution algorithm\n";
00182 return -1;
00183 }
00184
00185 AnalysisModel *theModel = this->getAnalysisModelPtr();
00186
00187 if (theModel == 0) {
00188 opserr << "ERROR CentralDifferenceAlternative::update() - no AnalysisModel set\n";
00189 return -2;
00190 }
00191
00192
00193 if (Ut == 0) {
00194 opserr << "WARNING CentralDifferenceAlternative::update() - domainChange() failed or not called\n";
00195 return -2;
00196 }
00197
00198
00199 if (X.Size() != Ut->Size()) {
00200 opserr << "WARNING CentralDifferenceAlternative::update() - Vectors of incompatable size ";
00201 opserr << " expecting " << Ut->Size() << " obtained " << X.Size() << endln;
00202 return -3;
00203 }
00204
00205
00206
00207 Utp1->addVector(0.0, X, deltaT * deltaT);
00208 (*Utp1) += *Ut;
00209 Utp1->addVector(1.0, *Udot, deltaT);
00210
00211
00212 (*Udot) = *Utp1;
00213 (*Udot) -= *Ut;
00214 (*Udot) *= (1.0/deltaT);
00215
00216
00217 theModel->setDisp(*Utp1);
00218 theModel->setVel(*Udot);
00219 theModel->updateDomain();
00220
00221 return 0;
00222 }
00223
00224 int
00225 CentralDifferenceAlternative::commit(void)
00226 {
00227 AnalysisModel *theModel = this->getAnalysisModelPtr();
00228 if (theModel == 0) {
00229 opserr << "WARNING CentralDifferenceAlternative::commit() - no AnalysisModel set\n";
00230 return -1;
00231 }
00232
00233 *Ut = *Utp1;
00234
00235
00236 double time = theModel->getCurrentDomainTime() + deltaT;
00237 theModel->setCurrentDomainTime(time);
00238
00239 return theModel->commitDomain();
00240
00241 return 0;
00242 }
00243
00244 int
00245 CentralDifferenceAlternative::sendSelf(int cTag, Channel &theChannel)
00246 {
00247 return 0;
00248 }
00249
00250 int
00251 CentralDifferenceAlternative::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
00252 {
00253 return 0;
00254 }
00255
00256 void
00257 CentralDifferenceAlternative::Print(OPS_Stream &s, int flag)
00258 {
00259 AnalysisModel *theModel = this->getAnalysisModelPtr();
00260 if (theModel != 0) {
00261 double currentTime = theModel->getCurrentDomainTime();
00262 s << "\t CentralDifferenceAlternative - currentTime: " << currentTime;
00263 } else
00264 s << "\t CentralDifferenceAlternative - no associated AnalysisModel\n";
00265 }
00266