DisplacementControl.cpp

Go to the documentation of this file.
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.10 $
00022 // $Date: 2005/10/19 21:53:57 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/integrator/DisplacementControl.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/analysis/integrator/DisplacementControl.C
00027 // 
00028 // Written: fmk 
00029 // Created: 07/98
00030 // Revision: A
00031 //
00032 // Description: This file contains the class definition for DisplacementControl.
00033 // DisplacementControl is an algorithmic class for perfroming a static analysis
00034 // using the arc length scheme, that is within a load step the follwing
00035 // constraint is enforced: dU^TdU + alpha^2*dLambda^2 = DisplacementControl^2
00036 // where dU is change in nodal displacements for step, dLambda is
00037 // change in applied load and DisplacementControl is a control parameter.
00038 //
00039 // What: "@(#) DisplacementControl.C, revA"
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   // to avoid divide-by-zero error on first update() ensure numIncr != 0
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     // delete any vector object created
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     // get pointers to AnalysisModel and LinearSOE
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     // determine increment for this iteration
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     // get the current load factor
00114     currentLambda = theModel->getCurrentDomainTime();
00115 
00116     // determine dUhat
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     // determine delta lambda(1) == dlambda    
00136     double dLambda = theIncrement/dUahat;
00137 
00138     deltaLambdaStep = dLambda;
00139     currentLambda += dLambda;
00140  //   opserr << "DisplacementControl: " << dUahat  << " " << theDofID << endln;
00141  //   opserr << "DisplacementControl::newStep() : " << deltaLambdaStep << endln;
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     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; // have to do this as the SOE is gonna change
00172     double dUabar = (*deltaUbar)(theDofID);
00173     
00174     // determine dUhat    
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     // determine delta lambda(1) == dlambda    
00187     double dLambda = -dUabar/dUahat;
00188     
00189     // determine delta U(i)
00190     (*deltaU) = (*deltaUbar);    
00191     deltaU->addVector(1.0, *deltaUhat,dLambda);
00192     
00193     // update dU and dlambda
00194     (*deltaUstep) += *deltaU;
00195     deltaLambdaStep += dLambda;
00196     currentLambda += dLambda;
00197 
00198     // update the model
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     // set the X soln in linearSOE to be deltaU for convergence Test
00208     theLinSOE->setX(*deltaU);
00209 
00210     numIncrLastStep++;
00211 
00212     return 0;
00213 }
00214 
00215 
00216 
00217 int 
00218 DisplacementControl::domainChanged(void)
00219 {
00220     // we first create the Vectors needed
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(); // ask model in case N+1 space
00229 
00230     if (deltaUhat == 0 || deltaUhat->Size() != size) { // create new Vector
00231         if (deltaUhat != 0)
00232             delete deltaUhat;   // delete the old
00233         deltaUhat = new Vector(size);
00234         if (deltaUhat == 0 || deltaUhat->Size() != size) { // check got it
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) { // create new Vector
00242         if (deltaUbar != 0)
00243             delete deltaUbar;   // delete the old
00244         deltaUbar = new Vector(size);
00245         if (deltaUbar == 0 || deltaUbar->Size() != size) { // check got it
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) { // create new Vector
00253         if (deltaU != 0)
00254             delete deltaU;   // delete the old
00255         deltaU = new Vector(size);
00256         if (deltaU == 0 || deltaU->Size() != size) { // check got it
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     // now we have to determine phat
00286     // do this by incrementing lambda by 1, applying load
00287     // and getting phat from unbalance.
00288     currentLambda = theModel->getCurrentDomainTime();
00289     currentLambda += 1.0;
00290     theModel->applyLoadDomain(currentLambda);    
00291     this->formUnbalance(); // NOTE: this assumes unbalance at last was 0
00292     (*phat) = theLinSOE->getB();
00293     currentLambda -= 1.0;
00294     theModel->setCurrentDomainTime(currentLambda);    
00295 
00296 
00297     // check there is a reference load
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     // lastly we determine the id of the nodal dof
00311     // EXTRA CODE TO DO SOME ERROR CHECKING REQUIRED
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     // TO FINISH
00326   return 0;
00327 }
00328 
00329 
00330 int
00331 DisplacementControl::recvSelf(int cTag,
00332                     Channel &theChannel, FEM_ObjectBroker &theBroker)
00333 {
00334     // TO FINISH
00335   return 0;
00336 }
00337 
00338 void
00339 DisplacementControl::Print(OPS_Stream &s, int flag)
00340 {
00341     // TO FINISH    
00342 }

Generated on Mon Oct 23 15:04:58 2006 for OpenSees by doxygen 1.5.0