C++ LinearSOE::getB方法代码示例

int
SecantAccelerator3::accelerate(Vector &vStar, LinearSOE &theSOE, 
			       IncrementalIntegrator &theIntegrator)
{
  // Current right hand side
  const Vector &rNew  = theSOE.getB();

  // Store for next iteration
  *r_1 = vStar;

  // No acceleration on first iteration
  if (iteration == 0) {
    // do nothing
  }
  else {
    // Store gamma in rOld ... \gamma = R_i - R_{i-1}
    rOld->addVector(-1.0, rNew, 1.0);

    // Store alpha in r_1 ... \alpha = r_i - r_{i-1}
    r_1->addVector(-1.0, vStar, 1.0);
    
    double den = 1.0 / ((*vOld)^(*rOld));
    double C   = ((*vOld)^rNew) * den;
    double A   = 1.0-C;
    double B   = -C - ((*r_1)^rNew) * den + C*((*r_1)^(*rOld)) * den;
    double D   = -C - A*(vStar^(*rOld))*den;
    double DA  = D/A;

    //opserr << "D = " << D << endln;
    //opserr << "B = " << B << ", C = " << C << endln;
    //opserr << "B+C = " << B+C << endln << endln;

    // Check "cut-out" criteria
    if (cutOut && (A > R1 || A < 1.0/R1 || DA > R2 || DA < -0.5*R2)) {
      // do nothing
      //opserr << "SecantAccelerator3::accelerate() -- cut out, A = " << A
      //	   << ", D/A = " << DA << endln;
    }
    else {
      vStar.addVector(A, *vOld, B);
      vStar.addVector(1.0, *rOld, C);
    }
  }

  // Store old values for next iteration
  *rOld = rNew;
  *vOld = vStar;

  iteration++;

  return 0; 
}

int 
AlgorithmIncrements::record(int cTag, double timeStamp)
{
  LinearSOE *theSOE = theAlgo->getLinearSOEptr();
  if (theSOE == 0)
    return 0;


  const Vector &B = theSOE->getB();
  const Vector &X = theSOE->getX();

  if (fileName != 0) {
    if (cTag == 0) {

      if (theFile)
	theFile.close();

      numRecord = 0;

      theFile.open(fileName, ios::out);
      if (!theFile) {
	opserr << "WARNING - AlgorithmIncrements::record()";
	opserr << " - could not open file " << fileName << endln;
      } 
    }
    if (theFile) {
      numRecord ++;
      int i;
      for (i=0; X.Size(); i++) theFile << X(i);
      for (i=0; X.Size(); i++) theFile << B(i);
    }
  }

  if (displayRecord == true)
    return this->plotData(X, B);

  return 0;
}

int 
BFGS::solveCurrentStep(void)
{
 
    // set up some pointers and check they are valid
    // NOTE this could be taken away if we set Ptrs as protecetd in superclass

    AnalysisModel   *theAnaModel = this->getAnalysisModelPtr();

    IncrementalIntegrator *theIntegrator = this->getIncrementalIntegratorPtr();

    LinearSOE  *theSOE = this->getLinearSOEptr();

    if ((theAnaModel == 0) || (theIntegrator == 0) || (theSOE == 0)
	|| (theTest == 0)){
	opserr << "WARNING BFGS::solveCurrentStep() - setLinks() has";
	opserr << " not been called - or no ConvergenceTest has been set\n";
	return -5;
    }	

    // set itself as the ConvergenceTest objects EquiSolnAlgo
    theTest->setEquiSolnAlgo(*this);
    if (theTest->start() < 0) {
      opserr << "BFGS::solveCurrentStep() -";
      opserr << "the ConvergenceTest object failed in start()\n";
      return -3;
    }

    localTest->setEquiSolnAlgo(*this);

    if (rdotz == 0)
       rdotz = new double[numberLoops+3];

    if (sdotr == 0)
	sdotr = new double[numberLoops+3];


    int result = -1;
    int count = 0;
    do {

      // opserr << "      BFGS -- Forming New Tangent" << endln;

      //form the initial tangent
      if (theIntegrator->formTangent(tangent) < 0){
         opserr << "WARNING BFGS::solveCurrentStep() -";
         opserr << "the Integrator failed in formTangent()\n";
         return -1; 
      }

      //form the initial residual 
      if (theIntegrator->formUnbalance() < 0) {
        opserr << "WARNING BFGS::solveCurrentStep() -";
        opserr << "the Integrator failed in formUnbalance()\n";	
      }	    

      //solve
      if (theSOE->solve() < 0) {
	  opserr << "WARNING BFGS::solveCurrentStep() -";
	  opserr << "the LinearSysOfEqn failed in solve()\n";	
	  return -3;
	}	    

      //update
      if ( theIntegrator->update(theSOE->getX() ) < 0) {
	opserr << "WARNING BFGS::solveCurrentStep() -";
	opserr << "the Integrator failed in update()\n";	
	return -4;
      }	        


      //    int systemSize = ( theSOE->getB() ).Size();
      int systemSize = theSOE->getNumEqn( );

      //temporary vector
      if (temp == 0 )
	temp = new Vector(systemSize);

      //initial displacement increment
      if ( s[1] == 0 ) 
	s[1] = new Vector(systemSize);

      *s[1] = theSOE->getX( );

      if ( residOld == 0 ) 
	residOld = new Vector(systemSize);

      *residOld = theSOE->getB( ) ;
      *residOld *= (-1.0 );

      //form the residual again
      if (theIntegrator->formUnbalance() < 0) {
        opserr << "WARNING BFGS::solveCurrentStep() -";
        opserr << "the Integrator failed in formUnbalance()\n";	
      }	    

      if ( residNew == 0 ) 
	residNew = new Vector(systemSize);
 
      if ( du == 0 ) 
//.........这里部分代码省略.........

int 
ArcLength::domainChanged(void)
{
    // we first create the Vectors needed
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();    
    if (theModel == 0 || theLinSOE == 0) {
	opserr << "WARNING ArcLength::update() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }    
    int size = theModel->getNumEqn(); // ask model in case N+1 space

    if (deltaUhat == 0 || deltaUhat->Size() != size) { // create new Vector
	if (deltaUhat != 0)
	    delete deltaUhat;   // delete the old
	deltaUhat = new Vector(size);
	if (deltaUhat == 0 || deltaUhat->Size() != size) { // check got it
	    opserr << "FATAL ArcLength::domainChanged() - ran out of memory for";
	    opserr << " deltaUhat Vector of size " << size << endln;
	    exit(-1);
	}
    }

    if (deltaUbar == 0 || deltaUbar->Size() != size) { // create new Vector
	if (deltaUbar != 0)
	    delete deltaUbar;   // delete the old
	deltaUbar = new Vector(size);
	if (deltaUbar == 0 || deltaUbar->Size() != size) { // check got it
	    opserr << "FATAL ArcLength::domainChanged() - ran out of memory for";
	    opserr << " deltaUbar Vector of size " << size << endln;
	    exit(-1);
	}
    }

    
    if (deltaU == 0 || deltaU->Size() != size) { // create new Vector
	if (deltaU != 0)
	    delete deltaU;   // delete the old
	deltaU = new Vector(size);
	if (deltaU == 0 || deltaU->Size() != size) { // check got it
	    opserr << "FATAL ArcLength::domainChanged() - ran out of memory for";
	    opserr << " deltaU Vector of size " << size << endln;
	    exit(-1);
	}
    }

    if (deltaUstep == 0 || deltaUstep->Size() != size) { 
	if (deltaUstep != 0)
	    delete deltaUstep;  
	deltaUstep = new Vector(size);
	if (deltaUstep == 0 || deltaUstep->Size() != size) { 
	    opserr << "FATAL ArcLength::domainChanged() - ran out of memory for";
	    opserr << " deltaUstep Vector of size " << size << endln;
	    exit(-1);
	}
    }

    if (phat == 0 || phat->Size() != size) { 
	if (phat != 0)
	    delete phat;  
	phat = new Vector(size);
	if (phat == 0 || phat->Size() != size) { 
	    opserr << "FATAL ArcLength::domainChanged() - ran out of memory for";
	    opserr << " phat Vector of size " << size << endln;
	    exit(-1);
	}
    }    

    // now we have to determine phat
    // do this by incrementing lambda by 1, applying load
    // and getting phat from unbalance.
    currentLambda = theModel->getCurrentDomainTime();
    currentLambda += 1.0;
    theModel->applyLoadDomain(currentLambda);    
    this->formUnbalance(); // NOTE: this assumes unbalance at last was 0
    (*phat) = theLinSOE->getB();
    currentLambda -= 1.0;
    theModel->setCurrentDomainTime(currentLambda);    
    

    // check there is a reference load
    int haveLoad = 0;
    for (int i=0; i<size; i++)
      if ( (*phat)(i) != 0.0 ) {
	haveLoad = 1;
	i = size;
      }

    if (haveLoad == 0) {
      opserr << "WARNING ArcLength::domainChanged() - zero reference load";
      return -1;
    }

    return 0;
}

int 
SecantLineSearch::search(double s0, 
				 double s1, 
				 LinearSOE &theSOE, 
				 IncrementalIntegrator &theIntegrator)
{
  double r0 = 0.0;

  if ( s0 != 0.0 ) 
    r0 = fabs( s1 / s0 );
	
  if  (r0 <= tolerance )
    return 0; // Line Search Not Required Residual Decrease Less Than Tolerance

  if (s1 == s0)
    return 0;  // Secant will have a divide-by-zero if continue

  // set some variables
  double eta    = 1.0;
  double s      = s1;
  double etaJ   = 1.0;
  double etaJm1 = 0.0;
  double sJ     = s1;
  double sJm1   = s0;
  double r = r0;

  const Vector &dU = theSOE.getX();

  if (printFlag == 0) {
    opserr << "Secant Line Search - initial: "
	 << "      eta(0) : " << eta << " , Ratio |s/s0| = " << r0 << endln;
  }

  // perform the secant iterations:
  //
  //                eta(j+1) = eta(j) -  s(j) * (eta(j-1)-eta(j))
  //                                     ------------------------
  //                                           s(j-1) - s(j)

  int count = 0; //intial value of iteration counter 
  while ( r > tolerance  &&  count < maxIter ) {
    
    count++;

    eta = etaJ - sJ * (etaJm1-etaJ) / (sJm1 - sJ);

    //-- want to put limits on eta and stop solution blowing up
    if (eta > maxEta)  eta = maxEta;
    if (r   > r0    )  eta =  1.0;
    if (eta < minEta)  eta = minEta;
    
    //update the incremental difference in response and determine new unbalance
    if (eta == etaJ) 
      break; // no change in response

    *x = dU;
    *x *= eta-etaJ;
	    
    if (theIntegrator.update(*x) < 0) {
      opserr << "WARNING SecantLineSearch::search() -";
      opserr << "the Integrator failed in update()\n";	
      return -1;
    }
    
    if (theIntegrator.formUnbalance() < 0) {
      opserr << "WARNING SecantLineSearch::search() -";
      opserr << "the Integrator failed in formUnbalance()\n";	
      return -2;
    }	

    //new residual
    const Vector &ResidJ = theSOE.getB();
    
    //new value of s
    s = dU ^ ResidJ;
    
    //new value of r 
    r = fabs( s / s0 ); 

    if (printFlag == 0) {
      opserr << "Secant Line Search - iteration: " << count 
	   << " , eta(j) : " << eta << " , Ratio |sj/s0| = " << r << endln;
    }

    if (etaJ == eta)
      count = maxIter;

    // set variables for next iteration
    etaJm1 = etaJ;
    etaJ = eta;
    sJm1 = sJ;
    sJ = s;

    if (sJm1 == sJ)
      count = maxIter;
    
  } //end while

  // set X in the SOE for the revised dU, needed for convergence tests
  *x = dU;
//.........这里部分代码省略.........

int
KrylovAccelerator2::accelerate(Vector &vStar, LinearSOE &theSOE, 
			      IncrementalIntegrator &theIntegrator)
{
  const Vector &R = theSOE.getB();

  int k = dimension;

  // Store residual for differencing at next iteration
  *(Av[k]) = R;

  // If subspace is not empty
  if (dimension > 0) {

    // Compute Av_k = f(y_{k-1}) - f(y_k) = r_{k-1} - r_k
    Av[k-1]->addVector(1.0, R, -1.0);
    
    int i,j;
    
    // Put subspace vectors into AvData
    Matrix A(AvData, numEqns, k);
    for (i = 0; i < k; i++) {
      Vector &Ai = *(Av[i]);
      for (j = 0; j < numEqns; j++)
	A(j,i) = Ai(j);
    }

    for (i = 0; i < k; i++) {
      for (int j = i+1; j < k; j++) {
	double sum = 0.0;
	double sumi = 0.0;
	double sumj = 0.0;
	for (int ii = 0; ii < numEqns; ii++) {
	  sum += A(ii,i)*A(ii,j);
	  sumi += A(ii,i)*A(ii,i);
	  sumj += A(ii,j)*A(ii,j);
	}
	sumi = sqrt(sumi);
	sumj = sqrt(sumj);
	sum = sum/(sumi*sumj);
	//if (fabs(sum) > 0.99)
	  //opserr << sum << ' ' << i << ' ' << j << "   ";
      }
    }

    // Put residual vector into rData (need to save r for later!)
    Vector B(rData, numEqns);
    B = R;
    
    // No transpose
    char *trans = "N";
    
    // The number of right hand side vectors
    int nrhs = 1;
    
    // Leading dimension of the right hand side vector
    int ldb = (numEqns > k) ? numEqns : k;
    
    // Subroutine error flag
    int info = 0;
    
    // Call the LAPACK least squares subroutine
#ifdef _WIN32
    unsigned int sizeC = 1;
    DGELS(trans, &sizeC, &numEqns, &k, &nrhs, AvData, &numEqns,
	  rData, &ldb, work, &lwork, &info);
#else
    //SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
    //		      $                  INFO )

    dgels_(trans, &numEqns, &k, &nrhs, AvData, &numEqns,
	   rData, &ldb, work, &lwork, &info);
#endif
    
    // Check for error returned by subroutine
    if (info < 0) {
      opserr << "WARNING KrylovAccelerator2::accelerate() - \n";
      opserr << "error code " << info << " returned by LAPACK dgels\n";
      return info;
    }
    
    Vector Q(numEqns);
    Q = R;

    // Compute the correction vector
    double cj;
    for (j = 0; j < k; j++) {
      
      // Solution to least squares is written to rData
      cj = rData[j];
      
      // Compute w_{k+1} = c_1 v_1 + ... + c_k v_k
      vStar.addVector(1.0, *(v[j]), cj);

      // Compute least squares residual
      // q_{k+1} = r_k - (c_1 Av_1 + ... + c_k Av_k)
      Q.addVector(1.0, *(Av[j]), -cj);
    }

    theSOE.setB(Q);
//.........这里部分代码省略.........

int 
NewtonLineSearch::solveCurrentStep(void)
{
    // set up some pointers and check they are valid
    // NOTE this could be taken away if we set Ptrs as protecetd in superclass
    AnalysisModel   *theAnaModel = this->getAnalysisModelPtr();
    IncrementalIntegrator *theIntegrator = this->getIncrementalIntegratorPtr();
    LinearSOE  *theSOE = this->getLinearSOEptr();

    if ((theAnaModel == 0) || (theIntegrator == 0) || (theSOE == 0)
	|| (theTest == 0)){
	opserr << "WARNING NewtonLineSearch::solveCurrentStep() - setLinks() has";
	opserr << " not been called - or no ConvergenceTest has been set\n";
	return -5;
    }	

    theLineSearch->newStep(*theSOE);

    // set itself as the ConvergenceTest objects EquiSolnAlgo
    theTest->setEquiSolnAlgo(*this);
    if (theTest->start() < 0) {
      opserr << "NewtonLineSearch::solveCurrentStep() -";
      opserr << "the ConvergenceTest object failed in start()\n";
      return -3;
    }

    if (theIntegrator->formUnbalance() < 0) {
      opserr << "WARNING NewtonLineSearch::solveCurrentStep() -";
      opserr << "the Integrator failed in formUnbalance()\n";	
      return -2;
    }	    

    int result = -1;
    do {

	//residual at this iteration before next solve 
	const Vector &Resid0 = theSOE->getB() ;
	
	//form the tangent
        if (theIntegrator->formTangent() < 0){
	    opserr << "WARNING NewtonLineSearch::solveCurrentStep() -";
	    opserr << "the Integrator failed in formTangent()\n";
	    return -1;
	}		    
	
	//solve 
	if (theSOE->solve() < 0) {
	    opserr << "WARNING NewtonLineSearch::solveCurrentStep() -";
	    opserr << "the LinearSysOfEqn failed in solve()\n";	
	    return -3;
	}	    


	//line search direction 
	const Vector &dx0 = theSOE->getX() ;

	//intial value of s
	double s0 = - (dx0 ^ Resid0) ; 

	if (theIntegrator->update(theSOE->getX()) < 0) {
	    opserr << "WARNING NewtonLineSearch::solveCurrentStep() -";
	    opserr << "the Integrator failed in update()\n";	
	    return -4;
	}	        

	if (theIntegrator->formUnbalance() < 0) {
	    opserr << "WARNING NewtonLineSearch::solveCurrentStep() -";
	    opserr << "the Integrator failed in formUnbalance()\n";	
	    return -2;
	}	

	// do a line search only if convergence criteria not met
	theOtherTest->start();
	result = theOtherTest->test();

	if (result < 1) {
	  //new residual 
	  const Vector &Resid = theSOE->getB() ;
	  
	  //new value of s 
	  double s = - ( dx0 ^ Resid ) ;
	  
	  if (theLineSearch != 0)
	    theLineSearch->search(s0, s, *theSOE, *theIntegrator);
	}

	this->record(0);
	  
	result = theTest->test();

    } while (result == -1);

    if (result == -2) {
      opserr << "NewtonLineSearch::solveCurrentStep() -";
      opserr << "the ConvergenceTest object failed in test()\n";
      return -3;
    }

    // note - if postive result we are returning what the convergence test returned
    // which should be the number of iterations
//.........这里部分代码省略.........

int 
InitialInterpolatedLineSearch::search(double s0, 
				      double s1, 
				      LinearSOE &theSOE, 
				      IncrementalIntegrator &theIntegrator)
{
  double s = s1;

  //intialize r = ratio of residuals 
  double r0 = 0.0;

  if ( s0 != 0.0 ) 
    r0 = fabs( s / s0 );
	
  if  (r0 <= tolerance )
    return 0; // Line Search Not Required Residual Decrease Less Than Tolerance

  double eta = 1.0;     //initial value of line search parameter
  double etaPrev = 1.0;
  double r = r0;

  const Vector &dU = theSOE.getX();

  int count = 0; //intial value of iteration counter 

  if (printFlag == 0) {
    opserr << "InitialInterpolated Line Search - initial       "
	 << "    eta : " << eta 
	 << " , Ratio |s/s0| = " << r0 << endln;
  }    


  // Solution procedure follows the one in Crissfields book.
  // (M.A. Crissfield, Nonlinear Finite Element Analysis of Solid and Structures, Wiley. 97).
  // NOTE: it is not quite linear interpolation/false-position/regula-falsi as eta(0) = 0.0
  // does not change. uses eta(i) = eta(i-1)*s0
  //                                -----------
  //                                s0 - s(i-1)  to compute eta(i)


  //  opserr << "dU: " << dU;

  while ( r > tolerance  &&  count < maxIter ) {

    count++;
    
    eta *=  s0 / (s0 - s); 

    //-- want to put limits on eta(i)
    if (eta > maxEta)  eta = maxEta;
    if (r   > r0    )  eta =  1.0;
    if (eta < minEta)  eta = minEta;

    if (eta == etaPrev)
      break; // no change in response break

    //dx = ( eta * dx0 ); 
    *x = dU;
    *x *= eta-etaPrev;
	    
    if (theIntegrator.update(*x) < 0) {
      opserr << "WARNInG InitialInterpolatedLineSearch::search() -";
      opserr << "the Integrator failed in update()\n";	
      return -1;
    }
    
    if (theIntegrator.formUnbalance() < 0) {
      opserr << "WARNInG InitialInterpolatedLineSearch::search() -";
      opserr << "the Integrator failed in formUnbalance()\n";	
      return -2;
    }	

    //new residual
    const Vector &ResidI = theSOE.getB();
    
    //new value of s
    s = dU ^ ResidI;

    //new value of r 
    r = fabs( s / s0 ); 

    if (printFlag == 0) {
      opserr << "InitialInterpolated Line Search - iteration: " << count 
	   << " , eta(j) : " << eta
	   << " , Ratio |sj/s0| = " << r << endln;
    }    

    // reset the variables, also check not just hitting bounds over and over
    if (eta == etaPrev)
      count = maxIter;
    else
      etaPrev = eta;

  } //end while

  // set X in the SOE for the revised dU, needed for convergence tests
  *x = dU;

  if (eta != 0.0)
    *x *= eta;
//.........这里部分代码省略.........

int HHTGeneralized_TP::newStep(double _deltaT)
{
    if (beta == 0 || gamma == 0 )  {
        opserr << "HHTGeneralized_TP::newStep() - error in variable\n";
        opserr << "gamma = " << gamma << " beta = " << beta << endln;
        return -1;
    }
    
    deltaT = _deltaT;
    if (deltaT <= 0.0)  {
        opserr << "HHTGeneralized_TP::newStep() - error in variable\n";
        opserr << "dT = " << deltaT << endln;
        return -2;
    }
    
    // get a pointer to the LinearSOE and the AnalysisModel
    LinearSOE *theLinSOE = this->getLinearSOE();
    AnalysisModel *theModel = this->getAnalysisModel();
    if (theLinSOE == 0 || theModel == 0)  {
        opserr << "WARNING HHT_TP::newStep() - ";
        opserr << "no LinearSOE or AnalysisModel has been set\n";
        return -3;
    }
    
    // set the constants
    c1 = 1.0;
    c2 = gamma/(beta*deltaT);
    c3 = 1.0/(beta*deltaT*deltaT);
       
    if (U == 0)  {
        opserr << "HHTGeneralized_TP::newStep() - domainChange() failed or hasn't been called\n";
        return -4;
    }
    
    // set response at t to be that at t+deltaT of previous step
    (*Ut) = *U;
    (*Utdot) = *Udot;
    (*Utdotdot) = *Udotdot;
    
    // get unbalance at t and store it
    alphaM = (1.0 - alphaI);
    alphaD = alphaR = alphaP = (1.0 - alphaF);
    this->TransientIntegrator::formUnbalance();
    (*Put) = theLinSOE->getB();
    
    // determine new velocities and accelerations at t+deltaT
    double a1 = (1.0 - gamma/beta);
    double a2 = deltaT*(1.0 - 0.5*gamma/beta);
    Udot->addVector(a1, *Utdotdot, a2);
    
    double a3 = -1.0/(beta*deltaT);
    double a4 = 1.0 - 0.5/beta;  
    Udotdot->addVector(a4, *Utdot, a3);
    
    // set the trial response quantities
    theModel->setVel(*Udot);
    theModel->setAccel(*Udotdot);
    
    // increment the time to t+deltaT and apply the load
    double time = theModel->getCurrentDomainTime();
    time += deltaT;
    if (theModel->updateDomain(time, deltaT) < 0)  {
        opserr << "HHTGeneralized_TP::newStep() - failed to update the domain\n";
        return -5;
    }
    
    // modify constants for subsequent iterations
    alphaM = alphaI;
    alphaD = alphaR = alphaP = alphaF;
    
    return 0;
}

int 
DisplacementControl::domainChanged(void)
{
    // we first create the Vectors needed
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();    
    if (theModel == 0 || theLinSOE == 0) {
	opserr << "WARNING DisplacementControl::update() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }    
    int size = theModel->getNumEqn(); // ask model in case N+1 space

    if (deltaUhat == 0 || deltaUhat->Size() != size) { // create new Vector
	if (deltaUhat != 0)
	    delete deltaUhat;   // delete the old
	deltaUhat = new Vector(size);
	if (deltaUhat == 0 || deltaUhat->Size() != size) { // check got it
	    opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
	    opserr << " deltaUhat Vector of size " << size << endln;
	    exit(-1);
	}
    }

    if (deltaUbar == 0 || deltaUbar->Size() != size) { // create new Vector
	if (deltaUbar != 0)
	    delete deltaUbar;   // delete the old
	deltaUbar = new Vector(size);
	if (deltaUbar == 0 || deltaUbar->Size() != size) { // check got it
	    opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
	    opserr << " deltaUbar Vector of size " << size << endln;
	    exit(-1);
	}
    }
    
    if (deltaU == 0 || deltaU->Size() != size) { // create new Vector
	if (deltaU != 0)
	    delete deltaU;   // delete the old
	deltaU = new Vector(size);
	if (deltaU == 0 || deltaU->Size() != size) { // check got it
	    opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
	    opserr << " deltaU Vector of size " << size << endln;
	    exit(-1);
	}
    }

    if (deltaUstep == 0 || deltaUstep->Size() != size) { 
	if (deltaUstep != 0)
	    delete deltaUstep;  
	deltaUstep = new Vector(size);
	if (deltaUstep == 0 || deltaUstep->Size() != size) { 
	    opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
	    opserr << " deltaUstep Vector of size " << size << endln;
	    exit(-1);
	}
    }

    if (phat == 0 || phat->Size() != size) { 
	if (phat != 0)
	    delete phat;  
	phat = new Vector(size);
	if (phat == 0 || phat->Size() != size) { 
	    opserr << "FATAL DisplacementControl::domainChanged() - ran out of memory for";
	    opserr << " phat Vector of size " << size << endln;
	    exit(-1);
	}
    }    

    // now we have to determine phat
    // do this by incrementing lambda by 1, applying load
    // and getting phat from unbalance.
    currentLambda = theModel->getCurrentDomainTime();
    currentLambda += 1.0;
    theModel->applyLoadDomain(currentLambda);    
    this->formUnbalance(); // NOTE: this assumes unbalance at last was 0
    (*phat) = theLinSOE->getB();
    currentLambda -= 1.0;
    theModel->setCurrentDomainTime(currentLambda);    


    // check there is a reference load
    int haveLoad = 0;
    for (int i=0; i<size; i++)
      if ( (*phat)(i) != 0.0 ) {
	haveLoad = 1;
	i = size;
      }

    if (haveLoad == 0) {
      opserr << "WARNING DisplacementControl::domainChanged() - zero reference load";
      return -1;
    }

    // lastly we determine the id of the nodal dof
    // EXTRA CODE TO DO SOME ERROR CHECKING REQUIRED
    
    Node *theNodePtr = theDomain->getNode(theNode);
	if (theNodePtr == 0) {
		opserr << "DisplacementControl::domainChanged - no node\n";
		return -1;
//.........这里部分代码省略.........