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
00037
00038
00039
00040
00041
00042 #include <DisplacementControl.h>
00043 #include <AnalysisModel.h>
00044 #include <LinearSOE.h>
00045 #include <Vector.h>
00046 #include <Channel.h>
00047 #include <math.h>
00048 #include <Domain.h>
00049 #include <Node.h>
00050 #include <DOF_Group.h>
00051 #include <ID.h>
00052
00053 DisplacementControl::DisplacementControl(int node, int dof,
00054 double increment,
00055 Domain *domain,
00056 int numIncr,
00057 double min, double max)
00058 :StaticIntegrator(INTEGRATOR_TAGS_DisplacementControl),
00059 theNode(node), theDof(dof), theIncrement(increment), theDomain(domain),
00060 theDofID(0),
00061 deltaUhat(0), deltaUbar(0), deltaU(0), deltaUstep(0),
00062 phat(0), deltaLambdaStep(0.0), currentLambda(0.0),
00063 specNumIncrStep(numIncr), numIncrLastStep(numIncr),
00064 minIncrement(min), maxIncrement(max)
00065 {
00066
00067 if (numIncr == 0) {
00068 opserr << "WARNING DisplacementControl::DisplacementControl() -";
00069 opserr << " numIncr set to 0, 1 assumed\n";
00070 specNumIncrStep = 1.0;
00071 numIncrLastStep = 1.0;
00072 }
00073 }
00074
00075 DisplacementControl::~DisplacementControl()
00076 {
00077
00078 if (deltaUhat != 0)
00079 delete deltaUhat;
00080 if (deltaU != 0)
00081 delete deltaU;
00082 if (deltaUstep != 0)
00083 delete deltaUstep;
00084 if (deltaUbar != 0)
00085 delete deltaUbar;
00086 if (phat != 0)
00087 delete phat;
00088 }
00089
00090 int
00091 DisplacementControl::newStep(void)
00092 {
00093
00094 AnalysisModel *theModel = this->getAnalysisModelPtr();
00095 LinearSOE *theLinSOE = this->getLinearSOEPtr();
00096 if (theModel == 0 || theLinSOE == 0) {
00097 opserr << "WARNING DisplacementControl::newStep() ";
00098 opserr << "No AnalysisModel or LinearSOE has been set\n";
00099 return -1;
00100 }
00101
00102
00103
00104 double factor = specNumIncrStep/numIncrLastStep;
00105 theIncrement *=factor;
00106
00107 if (theIncrement < minIncrement)
00108 theIncrement = minIncrement;
00109 else if (theIncrement > maxIncrement)
00110 theIncrement = maxIncrement;
00111
00112
00113
00114 currentLambda = theModel->getCurrentDomainTime();
00115
00116
00117 this->formTangent();
00118 theLinSOE->setB(*phat);
00119
00120 if (theLinSOE->solve() < 0) {
00121 opserr << "DisplacementControl::newStep(void) - failed in solver\n";
00122 return -1;
00123 }
00124
00125 (*deltaUhat) = theLinSOE->getX();
00126 Vector &dUhat = *deltaUhat;
00127
00128 double dUahat = dUhat(theDofID);
00129 if (dUahat == 0.0) {
00130 opserr << "WARNING DisplacementControl::newStep() ";
00131 opserr << "dUahat is zero -- zero reference displacement at control node DOF\n";
00132 return -1;
00133 }
00134
00135
00136 double dLambda = theIncrement/dUahat;
00137
00138 deltaLambdaStep = dLambda;
00139 currentLambda += dLambda;
00140
00141
00142
00143 (*deltaU) = dUhat;
00144 (*deltaU) *= dLambda;
00145 (*deltaUstep) = (*deltaU);
00146
00147
00148 theModel->incrDisp(*deltaU);
00149 theModel->applyLoadDomain(currentLambda);
00150 if (theModel->updateDomain() < 0) {
00151 opserr << "DisplacementControl::newStep - model failed to update for new dU\n";
00152 return -1;
00153 }
00154
00155 numIncrLastStep = 0;
00156
00157 return 0;
00158 }
00159
00160 int
00161 DisplacementControl::update(const Vector &dU)
00162 {
00163 AnalysisModel *theModel = this->getAnalysisModelPtr();
00164 LinearSOE *theLinSOE = this->getLinearSOEPtr();
00165 if (theModel == 0 || theLinSOE == 0) {
00166 opserr << "WARNING DisplacementControl::update() ";
00167 opserr << "No AnalysisModel or LinearSOE has been set\n";
00168 return -1;
00169 }
00170
00171 (*deltaUbar) = dU;
00172 double dUabar = (*deltaUbar)(theDofID);
00173
00174
00175 theLinSOE->setB(*phat);
00176 theLinSOE->solve();
00177 (*deltaUhat) = theLinSOE->getX();
00178
00179 double dUahat = (*deltaUhat)(theDofID);
00180 if (dUahat == 0.0) {
00181 opserr << "WARNING DisplacementControl::update() ";
00182 opserr << "dUahat is zero -- zero reference displacement at control node DOF\n";
00183 return -1;
00184 }
00185
00186
00187 double dLambda = -dUabar/dUahat;
00188
00189
00190 (*deltaU) = (*deltaUbar);
00191 deltaU->addVector(1.0, *deltaUhat,dLambda);
00192
00193
00194 (*deltaUstep) += *deltaU;
00195 deltaLambdaStep += dLambda;
00196 currentLambda += dLambda;
00197
00198
00199 theModel->incrDisp(*deltaU);
00200 theModel->applyLoadDomain(currentLambda);
00201 if (theModel->updateDomain() < 0) {
00202 opserr << "DisplacementControl::update - model failed to update for new dU\n";
00203 return -1;
00204 }
00205
00206
00207
00208 theLinSOE->setX(*deltaU);
00209
00210 numIncrLastStep++;
00211
00212 return 0;
00213 }
00214
00215
00216
00217 int
00218 DisplacementControl::domainChanged(void)
00219 {
00220
00221 AnalysisModel *theModel = this->getAnalysisModelPtr();
00222 LinearSOE *theLinSOE = this->getLinearSOEPtr();
00223 if (theModel == 0 || theLinSOE == 0) {
00224 opserr << "WARNING DisplacementControl::update() ";
00225 opserr << "No AnalysisModel or LinearSOE has been set\n";
00226 return -1;
00227 }
00228 int size = theModel->getNumEqn();
00229
00230 if (deltaUhat == 0 || deltaUhat->Size() != size) {
00231 if (deltaUhat != 0)
00232 delete deltaUhat;
00233 deltaUhat = new Vector(size);
00234 if (deltaUhat == 0 || deltaUhat->Size() != size) {
00235 opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
00236 opserr << " deltaUhat Vector of size " << size << endln;
00237 exit(-1);
00238 }
00239 }
00240
00241 if (deltaUbar == 0 || deltaUbar->Size() != size) {
00242 if (deltaUbar != 0)
00243 delete deltaUbar;
00244 deltaUbar = new Vector(size);
00245 if (deltaUbar == 0 || deltaUbar->Size() != size) {
00246 opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
00247 opserr << " deltaUbar Vector of size " << size << endln;
00248 exit(-1);
00249 }
00250 }
00251
00252 if (deltaU == 0 || deltaU->Size() != size) {
00253 if (deltaU != 0)
00254 delete deltaU;
00255 deltaU = new Vector(size);
00256 if (deltaU == 0 || deltaU->Size() != size) {
00257 opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
00258 opserr << " deltaU Vector of size " << size << endln;
00259 exit(-1);
00260 }
00261 }
00262
00263 if (deltaUstep == 0 || deltaUstep->Size() != size) {
00264 if (deltaUstep != 0)
00265 delete deltaUstep;
00266 deltaUstep = new Vector(size);
00267 if (deltaUstep == 0 || deltaUstep->Size() != size) {
00268 opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
00269 opserr << " deltaUstep Vector of size " << size << endln;
00270 exit(-1);
00271 }
00272 }
00273
00274 if (phat == 0 || phat->Size() != size) {
00275 if (phat != 0)
00276 delete phat;
00277 phat = new Vector(size);
00278 if (phat == 0 || phat->Size() != size) {
00279 opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
00280 opserr << " phat Vector of size " << size << endln;
00281 exit(-1);
00282 }
00283 }
00284
00285
00286
00287
00288 currentLambda = theModel->getCurrentDomainTime();
00289 currentLambda += 1.0;
00290 theModel->applyLoadDomain(currentLambda);
00291 this->formUnbalance();
00292 (*phat) = theLinSOE->getB();
00293 currentLambda -= 1.0;
00294 theModel->setCurrentDomainTime(currentLambda);
00295
00296
00297
00298 int haveLoad = 0;
00299 for (int i=0; i<size; i++)
00300 if ( (*phat)(i) != 0.0 ) {
00301 haveLoad = 1;
00302 i = size;
00303 }
00304
00305 if (haveLoad == 0) {
00306 opserr << "WARNING DisplacementControl::domainChanged() - zero reference load";
00307 return -1;
00308 }
00309
00310
00311
00312
00313 Node *theNodePtr = theDomain->getNode(theNode);
00314 DOF_Group *theGroup = theNodePtr->getDOF_GroupPtr();
00315 const ID &theID = theGroup->getID();
00316 theDofID = theID(theDof);
00317
00318 return 0;
00319 }
00320
00321 int
00322 DisplacementControl::sendSelf(int cTag,
00323 Channel &theChannel)
00324 {
00325
00326 return 0;
00327 }
00328
00329
00330 int
00331 DisplacementControl::recvSelf(int cTag,
00332 Channel &theChannel, FEM_ObjectBroker &theBroker)
00333 {
00334
00335 return 0;
00336 }
00337
00338 void
00339 DisplacementControl::Print(OPS_Stream &s, int flag)
00340 {
00341
00342 }