Method update

SRC/analysis/integrator/ArcLength.cpp:225–311

}

int
ArcLength::update(const Vector &dU)
{
  /// opserr<<" update: start"<<endln;

    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
//opserr<<"deltaUbar= "<<*deltaUbar<<endln;

//opserr<<"Update:   phat is = "<<*phat<<"////////////////////////////"<<endln;
    // determine dUhat    
    theLinSOE->setB(*phat);
    theLinSOE->solve();

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

    // determine the coeeficients of our quadratic equation
           a = alpha2 + ((*deltaUhat)^(*deltaUhat));
           b = alpha2*deltaLambdaStep 
      + ((*deltaUhat)^(*deltaUbar))
      + ((*deltaUstep)^(*deltaUhat));
      b *= 2.0;
      c = 2*((*deltaUstep)^(*deltaUbar)) + ((*deltaUbar)^(*deltaUbar));
    // check for a solution to quadratic
         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;
//opserr<<"squareRoot of b24ac= "<<sqrtb24ac<<endln;
    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;