My title

00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.3 $
00022 // $Date: 2001/05/30 01:42:35 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/integrator/MinUnbalDispNorm.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/analysis/integrator/MinUnbalDispNorm.C
00027 // 
00028 // Written: fmk 
00029 // Created: 07/99
00030 // Revision: A
00031 //
00032 // What: "@(#) MinUnbalDispNorm.C, revA"
00033 
00034 
00035 #include <MinUnbalDispNorm.h>
00036 #include <AnalysisModel.h>
00037 #include <LinearSOE.h>
00038 #include <Vector.h>
00039 #include <Channel.h>
00040 #include <iostream.h>
00041 #include <math.h>
00042 
00043 MinUnbalDispNorm::MinUnbalDispNorm(double lambda1, int specNumIter,
00044        double min, double max, int signFirstStep)
00045 :StaticIntegrator(INTEGRATOR_TAGS_MinUnbalDispNorm),
00046  dLambda1LastStep(lambda1), 
00047  specNumIncrStep(specNumIter), numIncrLastStep(specNumIter),
00048  deltaUhat(0), deltaUbar(0), deltaU(0), deltaUstep(0), 
00049  phat(0), deltaLambdaStep(0.0), currentLambda(0.0), 
00050  signLastDeltaLambdaStep(1),
00051  dLambda1min(min), dLambda1max(max), signLastDeterminant(1), signFirstStepMethod(signFirstStep)
00052 {
00053   // to avoid divide-by-zero error on first update() ensure numIncr != 0
00054   if (specNumIncrStep == 0) {
00055     cerr << "WARNING LoadControl::LoadControl() - numIncr set to 0, 1 assumed\n";
00056     specNumIncrStep = 1.0;
00057     numIncrLastStep = 1.0;
00058   }
00059 }
00060 
00061 MinUnbalDispNorm::~MinUnbalDispNorm()
00062 {
00063     // delete any vector object created
00064     if (deltaUhat != 0)
00065  delete deltaUhat;
00066     if (deltaU != 0)
00067  delete deltaU;
00068     if (deltaUstep != 0)
00069  delete deltaUstep;
00070     if (phat != 0)
00071  delete phat;
00072 }
00073 
00074 int
00075 MinUnbalDispNorm::newStep(void)
00076 {
00077     // get pointers to AnalysisModel and LinearSOE
00078     AnalysisModel *theModel = this->getAnalysisModelPtr();
00079     LinearSOE *theLinSOE = this->getLinearSOEPtr();    
00080     if (theModel == 0 || theLinSOE == 0) {
00081  cerr << "WARNING MinUnbalDispNorm::newStep() ";
00082  cerr << "No AnalysisModel or LinearSOE has been set\n";
00083  return -1;
00084     }
00085 
00086     // get the current load factor
00087     currentLambda = theModel->getCurrentDomainTime();
00088 
00089 
00090     // determine dUhat
00091     this->formTangent();
00092     theLinSOE->setB(*phat);
00093     theLinSOE->solve();
00094     (*deltaUhat) = theLinSOE->getX();
00095     Vector &dUhat = *deltaUhat;
00096 
00097     // determine delta lambda(1) == dlambda
00098     double factor = specNumIncrStep/numIncrLastStep;
00099     double dLambda = dLambda1LastStep*factor;
00100 
00101     // check aaint min and max values specified in constructor
00102     if (dLambda < dLambda1min)
00103       dLambda = dLambda1min;
00104     else if (dLambda > dLambda1max)
00105       dLambda = dLambda1max;
00106 
00107     dLambda1LastStep = dLambda;
00108 
00109 
00110     if (signFirstStepMethod == SIGN_LAST_STEP) {
00111       if (deltaLambdaStep < 0)
00112  signLastDeltaLambdaStep = -1;
00113       else
00114  signLastDeltaLambdaStep = +1;
00115       
00116       dLambda *= signLastDeltaLambdaStep; // base sign of load change
00117                                         // on what was happening last step
00118     } else {
00119     
00120       double det = theLinSOE->getDeterminant();
00121       int signDeterminant = 1;
00122       if (det < 0)
00123  signDeterminant = -1;
00124     
00125       dLambda *= signDeterminant * signLastDeterminant;
00126     
00127       signLastDeterminant = signDeterminant;
00128     }
00129 
00130     /*
00131     double work = (*phat)^(dUhat);
00132     int signCurrentWork = 1;
00133     if (work < 0) signCurrentWork = -1;
00134 
00135     if (signCurrentWork != signLastDeltaStep)
00136     */
00137 
00138     deltaLambdaStep = dLambda;
00139     currentLambda += dLambda;
00140     numIncrLastStep = 0;
00141 
00142     // determine delta U(1) == dU
00143     (*deltaU) = dUhat;
00144     (*deltaU) *= dLambda;
00145     (*deltaUstep) = (*deltaU);
00146 
00147     // update model with delta lambda and delta U
00148     theModel->incrDisp(*deltaU);    
00149     theModel->applyLoadDomain(currentLambda);    
00150     theModel->updateDomain();
00151 
00152     return 0;
00153 }
00154 
00155 int
00156 MinUnbalDispNorm::update(const Vector &dU)
00157 {
00158     AnalysisModel *theModel = this->getAnalysisModelPtr();
00159     LinearSOE *theLinSOE = this->getLinearSOEPtr();    
00160     if (theModel == 0 || theLinSOE == 0) {
00161  cerr << "WARNING MinUnbalDispNorm::update() ";
00162  cerr << "No AnalysisModel or LinearSOE has been set\n";
00163  return -1;
00164     }
00165 
00166     (*deltaUbar) = dU; // have to do this as the SOE is gonna change
00167 
00168     // determine dUhat    
00169     theLinSOE->setB(*phat);
00170     theLinSOE->solve();
00171     (*deltaUhat) = theLinSOE->getX();    
00172 
00173     // determine delta lambda(i)
00174     double a = (*deltaUhat)^(*deltaUbar);
00175     double b = (*deltaUhat)^(*deltaUhat);
00176     if (b == 0) {
00177       cerr << "MinUnbalDispNorm::update() - zero denominator\n";
00178       return -1;
00179     }
00180 
00181     double dLambda = -a/b;
00182     
00183     // determine delta U(i)
00184     (*deltaU) = (*deltaUbar);    
00185     deltaU->addVector(1.0, *deltaUhat,dLambda);
00186     
00187     // update dU and dlambda
00188     (*deltaUstep) += *deltaU;
00189     deltaLambdaStep += dLambda;
00190     currentLambda += dLambda;
00191 
00192     // update the model
00193     theModel->incrDisp(*deltaU);    
00194     theModel->applyLoadDomain(currentLambda);    
00195     theModel->updateDomain();
00196     
00197     // set the X soln in linearSOE to be deltaU for convergence Test
00198     for (int i=0; i<deltaU->Size(); i++)
00199  theLinSOE->setX(i, (*deltaU)(i));
00200 
00201     numIncrLastStep++;
00202     return 0;
00203 }
00204 
00205 
00206 
00207 int 
00208 MinUnbalDispNorm::domainChanged(void)
00209 {
00210     // we first create the Vectors needed
00211     AnalysisModel *theModel = this->getAnalysisModelPtr();
00212     LinearSOE *theLinSOE = this->getLinearSOEPtr();    
00213     if (theModel == 0 || theLinSOE == 0) {
00214  cerr << "WARNING MinUnbalDispNorm::update() ";
00215  cerr << "No AnalysisModel or LinearSOE has been set\n";
00216  return -1;
00217     }    
00218     int size = theModel->getNumEqn(); // ask model in case N+1 space
00219 
00220     if (deltaUhat == 0 || deltaUhat->Size() != size) { // create new Vector
00221  if (deltaUhat != 0)
00222      delete deltaUhat;   // delete the old
00223  deltaUhat = new Vector(size);
00224  if (deltaUhat == 0 || deltaUhat->Size() != size) { // check got it
00225      cerr << "FATAL MinUnbalDispNorm::domainChanged() - ran out of memory for";
00226      cerr << " deltaUhat Vector of size " << size << endl;
00227      exit(-1);
00228  }
00229     }
00230 
00231     if (deltaUbar == 0 || deltaUbar->Size() != size) { // create new Vector
00232  if (deltaUbar != 0)
00233      delete deltaUbar;   // delete the old
00234  deltaUbar = new Vector(size);
00235  if (deltaUbar == 0 || deltaUbar->Size() != size) { // check got it
00236      cerr << "FATAL MinUnbalDispNorm::domainChanged() - ran out of memory for";
00237      cerr << " deltaUbar Vector of size " << size << endl;
00238      exit(-1);
00239  }
00240     }
00241 
00242     
00243     if (deltaU == 0 || deltaU->Size() != size) { // create new Vector
00244  if (deltaU != 0)
00245      delete deltaU;   // delete the old
00246  deltaU = new Vector(size);
00247  if (deltaU == 0 || deltaU->Size() != size) { // check got it
00248      cerr << "FATAL MinUnbalDispNorm::domainChanged() - ran out of memory for";
00249      cerr << " deltaU Vector of size " << size << endl;
00250      exit(-1);
00251  }
00252     }
00253 
00254     if (deltaUstep == 0 || deltaUstep->Size() != size) { 
00255  if (deltaUstep != 0)
00256      delete deltaUstep;  
00257  deltaUstep = new Vector(size);
00258  if (deltaUstep == 0 || deltaUstep->Size() != size) { 
00259      cerr << "FATAL MinUnbalDispNorm::domainChanged() - ran out of memory for";
00260      cerr << " deltaUstep Vector of size " << size << endl;
00261      exit(-1);
00262  }
00263     }
00264 
00265     if (phat == 0 || phat->Size() != size) { 
00266  if (phat != 0)
00267      delete phat;  
00268  phat = new Vector(size);
00269  if (phat == 0 || phat->Size() != size) { 
00270      cerr << "FATAL MinUnbalDispNorm::domainChanged() - ran out of memory for";
00271      cerr << " phat Vector of size " << size << endl;
00272      exit(-1);
00273  }
00274     }    
00275 
00276     // now we have to determine phat
00277     // do this by incrementing lambda by 1, applying load
00278     // and getting phat from unbalance.
00279     currentLambda = theModel->getCurrentDomainTime();
00280     currentLambda += 1.0;
00281     theModel->applyLoadDomain(currentLambda);    
00282     this->formUnbalance(); // NOTE: this assumes unbalance at last was 0
00283     (*phat) = theLinSOE->getB();
00284     currentLambda -= 1.0;
00285     theModel->setCurrentDomainTime(currentLambda);    
00286 
00287 
00288     // check there is a reference load
00289     int haveLoad = 0;
00290     for (int i=0; i<size; i++)
00291       if ( (*phat)(i) != 0.0 ) {
00292  haveLoad = 1;
00293  i = size;
00294       }
00295 
00296     if (haveLoad == 0) {
00297       cerr << "WARNING ArcLength::domainChanged() - zero reference load";
00298       return -1;
00299     }
00300 
00301     
00302     return 0;
00303 }
00304 
00305 int
00306 MinUnbalDispNorm::sendSelf(int cTag,
00307       Channel &theChannel)
00308 {
00309   Vector data(8);
00310   data(0) = dLambda1LastStep;
00311   data(1) = specNumIncrStep;
00312   data(2) = numIncrLastStep;
00313   data(3) = deltaLambdaStep;
00314   data(4) = currentLambda;
00315   if (signLastDeltaLambdaStep == 1)
00316     data(5)  = 1.0;
00317   else
00318     data(5) = 0.0;
00319   data(6) = dLambda1min;
00320   data(7) = dLambda1max;
00321 
00322   if (theChannel.sendVector(this->getDbTag(), cTag, data) < 0) {
00323       cerr << "MinUnbalDispNorm::sendSelf() - failed to send the data\n";
00324       return -1;
00325   }
00326   return 0;
00327 }
00328 
00329 
00330 int
00331 MinUnbalDispNorm::recvSelf(int cTag,
00332       Channel &theChannel, FEM_ObjectBroker &theBroker)
00333 {
00334   Vector data(8);
00335   if (theChannel.recvVector(this->getDbTag(), cTag, data) < 0) {
00336       cerr << "MinUnbalDispNorm::sendSelf() - failed to send the data\n";
00337       return -1;
00338   }      
00339 
00340   // set the data
00341 
00342   dLambda1LastStep = data(0);
00343   specNumIncrStep = data(1);
00344   numIncrLastStep = data(2);
00345   deltaLambdaStep = data(3);
00346   currentLambda = data(4);
00347   if (data(5)== 1.0)
00348     signLastDeltaLambdaStep = 1;
00349   else
00350     signLastDeltaLambdaStep = -1;
00351   dLambda1min = data(6);
00352   dLambda1max = data(7);
00353 
00354   return 0;
00355 }
00356 
00357 void
00358 MinUnbalDispNorm::Print(ostream &s, int flag)
00359 {
00360     AnalysisModel *theModel = this->getAnalysisModelPtr();
00361     if (theModel != 0) {
00362  double cLambda = theModel->getCurrentDomainTime();
00363  s << "\t MinUnbalDispNorm - currentLambda: " << cLambda;
00364     } else 
00365  s << "\t MinUnbalDispNorm - no associated AnalysisModel\n";
00366 }