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 <ArcLength1.h>
00043 #include <AnalysisModel.h>
00044 #include <LinearSOE.h>
00045 #include <Vector.h>
00046 #include <Channel.h>
00047 #include <math.h>
00048
00049 ArcLength1::ArcLength1(double arcLength, double alpha)
00050 :StaticIntegrator(INTEGRATOR_TAGS_ArcLength1),
00051 arcLength2(arcLength*arcLength), alpha2(alpha*alpha),
00052 deltaUhat(0), deltaUbar(0), deltaU(0), deltaUstep(0),
00053 phat(0), deltaLambdaStep(0.0), currentLambda(0.0),
00054 signLastDeltaLambdaStep(1)
00055 {
00056
00057 }
00058
00059 ArcLength1::~ArcLength1()
00060 {
00061
00062 if (deltaUhat != 0)
00063 delete deltaUhat;
00064 if (deltaU != 0)
00065 delete deltaU;
00066 if (deltaUstep != 0)
00067 delete deltaUstep;
00068 if (deltaUbar != 0)
00069 delete deltaUbar;
00070 if (phat != 0)
00071 delete phat;
00072 }
00073
00074 int
00075 ArcLength1::newStep(void)
00076 {
00077
00078 AnalysisModel *theModel = this->getAnalysisModelPtr();
00079 LinearSOE *theLinSOE = this->getLinearSOEPtr();
00080 if (theModel == 0 || theLinSOE == 0) {
00081 opserr << "WARNING ArcLength1::newStep() ";
00082 opserr << "No AnalysisModel or LinearSOE has been set\n";
00083 return -1;
00084 }
00085
00086
00087 currentLambda = theModel->getCurrentDomainTime();
00088
00089 if (deltaLambdaStep < 0)
00090 signLastDeltaLambdaStep = -1;
00091 else
00092 signLastDeltaLambdaStep = +1;
00093
00094
00095 this->formTangent();
00096 theLinSOE->setB(*phat);
00097 theLinSOE->solve();
00098 (*deltaUhat) = theLinSOE->getX();
00099 Vector &dUhat = *deltaUhat;
00100
00101
00102 double dLambda = sqrt(arcLength2/((dUhat^dUhat)+alpha2));
00103 dLambda *= signLastDeltaLambdaStep;
00104
00105 deltaLambdaStep = dLambda;
00106 currentLambda += dLambda;
00107
00108
00109 (*deltaU) = dUhat;
00110 (*deltaU) *= dLambda;
00111 (*deltaUstep) = (*deltaU);
00112
00113
00114 theModel->incrDisp(*deltaU);
00115 theModel->applyLoadDomain(currentLambda);
00116 theModel->updateDomain();
00117
00118 return 0;
00119 }
00120
00121 int
00122 ArcLength1::update(const Vector &dU)
00123 {
00124 AnalysisModel *theModel = this->getAnalysisModelPtr();
00125 LinearSOE *theLinSOE = this->getLinearSOEPtr();
00126 if (theModel == 0 || theLinSOE == 0) {
00127 opserr << "WARNING ArcLength1::update() ";
00128 opserr << "No AnalysisModel or LinearSOE has been set\n";
00129 return -1;
00130 }
00131
00132 (*deltaUbar) = dU;
00133
00134
00135 theLinSOE->setB(*phat);
00136 theLinSOE->solve();
00137 (*deltaUhat) = theLinSOE->getX();
00138
00139
00140 double a = (*deltaUstep)^(*deltaUbar);
00141 double b = (*deltaUstep)^(*deltaUhat) + alpha2*deltaLambdaStep;
00142 if (b == 0) {
00143 opserr << "ArcLength1::update() - zero denominator,";
00144 opserr << " alpha was set to 0.0 and zero reference load\n";
00145 return -1;
00146 }
00147 double dLambda = -a/b;
00148
00149
00150 (*deltaU) = (*deltaUbar);
00151 deltaU->addVector(1.0, *deltaUhat,dLambda);
00152
00153
00154 (*deltaUstep) += *deltaU;
00155 deltaLambdaStep += dLambda;
00156 currentLambda += dLambda;
00157
00158
00159 theModel->incrDisp(*deltaU);
00160 theModel->applyLoadDomain(currentLambda);
00161 theModel->updateDomain();
00162
00163
00164 theLinSOE->setX(*deltaU);
00165
00166 return 0;
00167 }
00168
00169
00170
00171 int
00172 ArcLength1::domainChanged(void)
00173 {
00174
00175 AnalysisModel *theModel = this->getAnalysisModelPtr();
00176 LinearSOE *theLinSOE = this->getLinearSOEPtr();
00177 if (theModel == 0 || theLinSOE == 0) {
00178 opserr << "WARNING ArcLength1::update() ";
00179 opserr << "No AnalysisModel or LinearSOE has been set\n";
00180 return -1;
00181 }
00182 int size = theModel->getNumEqn();
00183
00184 if (deltaUhat == 0 || deltaUhat->Size() != size) {
00185 if (deltaUhat != 0)
00186 delete deltaUhat;
00187 deltaUhat = new Vector(size);
00188 if (deltaUhat == 0 || deltaUhat->Size() != size) {
00189 opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00190 opserr << " deltaUhat Vector of size " << size << endln;
00191 exit(-1);
00192 }
00193 }
00194
00195 if (deltaUbar == 0 || deltaUbar->Size() != size) {
00196 if (deltaUbar != 0)
00197 delete deltaUbar;
00198 deltaUbar = new Vector(size);
00199 if (deltaUbar == 0 || deltaUbar->Size() != size) {
00200 opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00201 opserr << " deltaUbar Vector of size " << size << endln;
00202 exit(-1);
00203 }
00204 }
00205
00206
00207 if (deltaU == 0 || deltaU->Size() != size) {
00208 if (deltaU != 0)
00209 delete deltaU;
00210 deltaU = new Vector(size);
00211 if (deltaU == 0 || deltaU->Size() != size) {
00212 opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00213 opserr << " deltaU Vector of size " << size << endln;
00214 exit(-1);
00215 }
00216 }
00217
00218 if (deltaUstep == 0 || deltaUstep->Size() != size) {
00219 if (deltaUstep != 0)
00220 delete deltaUstep;
00221 deltaUstep = new Vector(size);
00222 if (deltaUstep == 0 || deltaUstep->Size() != size) {
00223 opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00224 opserr << " deltaUstep Vector of size " << size << endln;
00225 exit(-1);
00226 }
00227 }
00228
00229 if (phat == 0 || phat->Size() != size) {
00230 if (phat != 0)
00231 delete phat;
00232 phat = new Vector(size);
00233 if (phat == 0 || phat->Size() != size) {
00234 opserr << "FATAL ArcLength1::domainChanged() - ran out of memory for";
00235 opserr << " phat Vector of size " << size << endln;
00236 exit(-1);
00237 }
00238 }
00239
00240
00241
00242
00243 currentLambda = theModel->getCurrentDomainTime();
00244 currentLambda += 1.0;
00245 theModel->applyLoadDomain(currentLambda);
00246 this->formUnbalance();
00247 (*phat) = theLinSOE->getB();
00248 currentLambda -= 1.0;
00249 theModel->setCurrentDomainTime(currentLambda);
00250
00251 return 0;
00252 }
00253
00254 int
00255 ArcLength1::sendSelf(int cTag,
00256 Channel &theChannel)
00257 {
00258 Vector data(5);
00259 data(0) = arcLength2;
00260 data(1) = alpha2;
00261 data(2) = deltaLambdaStep;
00262 data(3) = currentLambda;
00263 data(4) = signLastDeltaLambdaStep;
00264
00265 if (theChannel.sendVector(this->getDbTag(), cTag, data) < 0) {
00266 opserr << "ArcLength1::sendSelf() - failed to send the data\n";
00267 return -1;
00268 }
00269 return 0;
00270 }
00271
00272
00273 int
00274 ArcLength1::recvSelf(int cTag,
00275 Channel &theChannel, FEM_ObjectBroker &theBroker)
00276 {
00277 Vector data(5);
00278 if (theChannel.recvVector(this->getDbTag(), cTag, data) < 0) {
00279 opserr << "ArcLength1::sendSelf() - failed to send the data\n";
00280 return -1;
00281 }
00282
00283
00284 arcLength2 = data(0);
00285 alpha2 = data(1);
00286 deltaLambdaStep = data(2);
00287 currentLambda = data(3);
00288 signLastDeltaLambdaStep = data(4);
00289 return 0;
00290 }
00291
00292 void
00293 ArcLength1::Print(OPS_Stream &s, int flag)
00294 {
00295 AnalysisModel *theModel = this->getAnalysisModelPtr();
00296 if (theModel != 0) {
00297 double cLambda = theModel->getCurrentDomainTime();
00298 s << "\t ArcLength1 - currentLambda: " << cLambda;
00299 s << " ArcLength1: " << sqrt(arcLength2) << " alpha: ";
00300 s << sqrt(alpha2) << endln;
00301 } else
00302 s << "\t ArcLength1 - no associated AnalysisModel\n";
00303 }
00304
00305
00306
00307
00308
00309
00310
00311