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

int
DisplacementControl::update(const Vector &dU)
{

  if (theDofID == -1) {
    opserr << "DisplacementControl::newStep() - domainChanged has not been called\n";
    return -1;
  } 
    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;
    }

    (*deltaUbar) = dU; // have to do this as the SOE is gonna change
    double dUabar = (*deltaUbar)(theDofID);
    
    // determine dUhat    
    theLinSOE->setB(*phat);
    theLinSOE->solve();
    (*deltaUhat) = theLinSOE->getX();    

    double dUahat = (*deltaUhat)(theDofID);
    if (dUahat == 0.0) {
	opserr << "WARNING DisplacementControl::update() ";
	opserr << "dUahat is zero -- zero reference displacement at control node DOF\n";
	return -1;
    }
    
    // determine delta lambda(1) == dlambda    
    double dLambda = -dUabar/dUahat;
    
    // determine delta U(i)
    (*deltaU) = (*deltaUbar);    
    deltaU->addVector(1.0, *deltaUhat,dLambda);
    
    // update dU and dlambda
    (*deltaUstep) += *deltaU;
    deltaLambdaStep += dLambda;
    currentLambda += dLambda;

    // update the model
    theModel->incrDisp(*deltaU);    
    theModel->applyLoadDomain(currentLambda);    
    if (theModel->updateDomain() < 0) {
      opserr << "DisplacementControl::update - model failed to update for new dU\n";
      return -1;
    }
	
    
    // set the X soln in linearSOE to be deltaU for convergence Test
    theLinSOE->setX(*deltaU);

    numIncrLastStep++;

    return 0;
}

int
ArcLength1::update(const Vector &dU)
{
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();    
    if (theModel == 0 || theLinSOE == 0) {
	opserr << "WARNING ArcLength1::update() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }

    (*deltaUbar) = dU; // have to do this as the SOE is gonna change

    // determine dUhat    
    theLinSOE->setB(*phat);
    theLinSOE->solve();
    (*deltaUhat) = theLinSOE->getX();    

    // determine delta lambda(i)
    double a = (*deltaUstep)^(*deltaUbar);
    double b = (*deltaUstep)^(*deltaUhat) + alpha2*deltaLambdaStep;
    if (b == 0) {
      opserr << "ArcLength1::update() - zero denominator,";
      opserr << " alpha was set to 0.0 and zero reference load\n";
      return -1;
    }
    double dLambda = -a/b;

    // determine delta U(i)
    (*deltaU) = (*deltaUbar);    
    deltaU->addVector(1.0, *deltaUhat,dLambda);
    
    // update dU and dlambda
    (*deltaUstep) += *deltaU;
    deltaLambdaStep += dLambda;
    currentLambda += dLambda;

    // update the model
    theModel->incrDisp(*deltaU);    
    theModel->applyLoadDomain(currentLambda);    
    theModel->updateDomain();
    
    // set the X soln in linearSOE to be deltaU for convergence Test
    theLinSOE->setX(*deltaU);

    return 0;
}

int
ArcLengthw::update(const Vector &dU)
{
    ofstream factor;
    factor.open("FS.dat",ios::app);

    factor<<"insideupdate"<<endln;
    //factor>>dU;
    factor<<"insideupdate1"<<endln;

    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();    
    if (theModel == 0 || theLinSOE == 0) {
	opserr << "WARNING ArcLengthw::update() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }

    (*deltaUbar) = dU; // have to do this as the SOE is gonna change

    // determine dUhat    
    theLinSOE->setB(*phat);
    theLinSOE->solve();
    (*deltaUhat) = theLinSOE->getX();    

    double dLambda = -((*phat)^(*deltaUbar))/((*phat)^(*deltaUhat));

    (*deltaU) = (*deltaUbar);    
    deltaU->addVector(1.0, *deltaUhat,dLambda);
    
    // update dU and dlambda
    (*deltaUstep) += *deltaU;
    deltaLambdaStep += dLambda;
    currentLambda += dLambda;

    // update the model
    theModel->incrDisp(*deltaU);    
    theModel->applyLoadDomain(currentLambda);    
    theModel->updateDomain();
    
    // set the X soln in linearSOE to be deltaU for convergence Test
    theLinSOE->setX(*deltaU);

    return 0;
}

int
LoadControl::update(const Vector &deltaU)
{
    AnalysisModel *myModel = this->getAnalysisModel();
    LinearSOE *theSOE = this->getLinearSOE();
    if (myModel == 0 || theSOE == 0) {
	opserr << "WARNING LoadControl::update() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }

    myModel->incrDisp(deltaU);    
    if (myModel->updateDomain() < 0) {
      opserr << "LoadControl::update - model failed to update for new dU\n";
      return -1;
    }

    // Set deltaU for the convergence test
    theSOE->setX(deltaU);

    numIncrLastStep++;

    return 0;
}

int
ArcLength::update(const Vector &dU)
{
    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;
    }

    (*deltaUbar) = dU; // have to do this as the SOE is gonna change

    // determine dUhat    
    theLinSOE->setB(*phat);
    theLinSOE->solve();

    (*deltaUhat) = theLinSOE->getX();    

    // determine the coeeficients of our quadratic equation
    double a = alpha2 + ((*deltaUhat)^(*deltaUhat));
    double b = alpha2*deltaLambdaStep 
      + ((*deltaUhat)^(*deltaUbar))
      + ((*deltaUstep)^(*deltaUhat));
    b *= 2.0;
    double c = 2*((*deltaUstep)^(*deltaUbar)) + ((*deltaUbar)^(*deltaUbar));
    // check for a solution to quadratic
    double b24ac = b*b - 4.0*a*c;
    if (b24ac < 0) {
      opserr << "ArcLength::update() - imaginary roots due to multiple instability";
      opserr << " directions - initial load increment was too large\n";
      opserr << "a: " << a << " b: " << b << " c: " << c << " b24ac: " << b24ac << endln;
      return -1;
    }			       
    double a2 = 2.0*a;
    if (a2 == 0.0) {
      opserr << "ArcLength::update() - zero denominator";
      opserr << " alpha was set to 0.0 and zero reference load\n";
      return -2;
    }			       

    // determine the roots of the quadratic
    double sqrtb24ac = sqrt(b24ac);
    double dlambda1 = (-b + sqrtb24ac)/a2;
    double dlambda2 = (-b - sqrtb24ac)/a2;

    double val = (*deltaUhat)^(*deltaUstep);
    double theta1 = ((*deltaUstep)^(*deltaUstep)) + ((*deltaUbar)^(*deltaUstep));
    //    double theta2 = theta1 + dlambda2*val;
    theta1 += dlambda1*val;

    // choose dLambda based on angle between incremental displacement before
    // and after this step -- want positive
    double dLambda;
    if (theta1 > 0)
      dLambda = dlambda1;
    else
      dLambda = dlambda2;


    // determine delta U(i)
    (*deltaU) = (*deltaUbar);    
    deltaU->addVector(1.0, *deltaUhat,dLambda);
    
    // update dU and dlambda
    (*deltaUstep) += *deltaU;
    deltaLambdaStep += dLambda;
    currentLambda += dLambda;

    // update the model
    theModel->incrDisp(*deltaU);    
    theModel->applyLoadDomain(currentLambda);    


    theModel->updateDomain();
    
    // set the X soln in linearSOE to be deltaU for convergence Test
    theLinSOE->setX(*deltaU);

    return 0;
}

//.........这里部分代码省略.........
  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;
  if (eta != 0.0)
    *x *= eta;
  theSOE.setX(*x);

  return 0;
}

//.........这里部分代码省略.........
{
  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;

  theSOE.setX(*x);

  return 0;
}

int
DisplacementPath::update(const Vector &dU)
{
    // opserr << " update is invoked " << endln;
	if (theDofID == -1) {
		opserr << "DisplacementControl::newStep() - domainChanged has not been called\n";
		return -1;
	}
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();    
    if (theModel == 0 || theLinSOE == 0) {
	opserr << "WARNING DisplacementPath::update() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }

    (*deltaUbar) = dU; // have to do this as the SOE is gonna change
    double dUabar = (*deltaUbar)(theDofID);
    
    // determine dUhat    
    theLinSOE->setB(*phat);
    theLinSOE->solve();
    (*deltaUhat) = theLinSOE->getX();    

	// add by zhong for check purpose
    //int size = deltaUhat->Size();
	//opserr << "\n size of deltaUhat = " << size << endln;
    //for (int i=0; i<size; i++) {
    //   opserr << " dektaUhat(i) = " << (*deltaUhat)(i) << endln;
    //}
	// finish here

    double dUahat = (*deltaUhat)(theDofID);
	//opserr << " theDofID = " << theDofID << endln;
	
	//opserr << "update( ) " << endln;
	//opserr << "dUahat = " << dUahat << endln;

    if (dUahat == 0.0) {
	opserr << "WARNING DisplacementPath::update() ";
	opserr << "dUahat is zero -- zero reference displacement at control node DOF\n";
	return -1;
    }
    
    // determine delta lambda(1) == dlambda    
    double dLambda = -dUabar/dUahat;

    // add by zhong
	//opserr << "\n dUahat = " << dUahat << endln;
	//opserr << " dUabar = " << dUabar << endln;
	//opserr << " dLambda = " << dLambda << endln;
	// finish
    
    // determine delta U(i)
    (*deltaU) = (*deltaUbar);    
    deltaU->addVector(1.0, *deltaUhat,dLambda);
    
    // update dU and dlambda
    (*deltaUstep) += *deltaU;
    deltaLambdaStep += dLambda;
    currentLambda += dLambda;

    // update the model
    theModel->incrDisp(*deltaU);    
    theModel->applyLoadDomain(currentLambda);    
    if (theModel->updateDomain() < 0) {
      opserr << "DisplacementPath::update - model failed to update for new dU\n";
      return -1;
    }
	
    
    // set the X soln in linearSOE to be deltaU for convergence Test
    theLinSOE->setX(*deltaU);


    return 0;
}