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

int HHTGeneralized_TP::formUnbalance()
{
    // get a pointer to the LinearSOE and the AnalysisModel
    LinearSOE *theLinSOE = this->getLinearSOE();
    AnalysisModel *theModel = this->getAnalysisModel();
    if (theLinSOE == 0 || theModel == 0)  {
        opserr << "WARNING HHTGeneralized_TP::formUnbalance() - ";
        opserr << "no LinearSOE or AnalysisModel has been set\n";
        return -1;
    }
    
    theLinSOE->setB(*Put);
    
    // do modal damping
    const Vector *modalValues = theModel->getModalDampingFactors();
    if (modalValues != 0)  {
        this->addModalDampingForce(modalValues);
    }
    
    if (this->formElementResidual() < 0)  {
        opserr << "WARNING HHTGeneralized_TP::formUnbalance ";
        opserr << " - this->formElementResidual failed\n";
        return -2;
    }
    
    if (this->formNodalUnbalance() < 0)  {
        opserr << "WARNING HHTGeneralized_TP::formUnbalance ";
        opserr << " - this->formNodalUnbalance failed\n";
        return -3;
    }
    
    return 0;
}

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
ArcLength::newStep(void)
{
    // get pointers to AnalysisModel and LinearSOE
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();    
    if (theModel == 0 || theLinSOE == 0) {
	opserr << "WARNING ArcLength::newStep() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }

    // get the current load factor
    currentLambda = theModel->getCurrentDomainTime();

    if (deltaLambdaStep < 0)
	signLastDeltaLambdaStep = -1;
    else
	signLastDeltaLambdaStep = +1;

    // determine dUhat
    this->formTangent();
    theLinSOE->setB(*phat);
    if (theLinSOE->solve() < 0) {
      opserr << "ArcLength::newStep(void) - failed in solver\n";
      return -1;
    }

    (*deltaUhat) = theLinSOE->getX();
    Vector &dUhat = *deltaUhat;
    
    // determine delta lambda(1) == dlambda
    double dLambda = sqrt(arcLength2/((dUhat^dUhat)+alpha2));
    dLambda *= signLastDeltaLambdaStep; // base sign of load change
                                        // on what was happening last step
    deltaLambdaStep = dLambda;
    currentLambda += dLambda;

    // determine delta U(1) == dU
    (*deltaU) = dUhat;
    (*deltaU) *= dLambda;
    (*deltaUstep) = (*deltaU);

    // update model with delta lambda and delta U
    theModel->incrDisp(*deltaU);    
    theModel->applyLoadDomain(currentLambda);    
    theModel->updateDomain();

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

    // get pointers to AnalysisModel and LinearSOE
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();    
    if (theModel == 0 || theLinSOE == 0) {
	opserr << "WARNING DisplacementControl::newStep() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }


    // determine increment for this iteration
    double factor = specNumIncrStep/numIncrLastStep;
    theIncrement *=factor;

    if (theIncrement < minIncrement)
      theIncrement = minIncrement;
    else if (theIncrement > maxIncrement)
      theIncrement = maxIncrement;


    // get the current load factor
    currentLambda = theModel->getCurrentDomainTime();

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

    if (theLinSOE->solve() < 0) {
      opserr << "DisplacementControl::newStep(void) - failed in solver\n";
      return -1;
    }

    (*deltaUhat) = theLinSOE->getX();
    Vector &dUhat = *deltaUhat;

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

    deltaLambdaStep = dLambda;
    currentLambda += dLambda;
 //   opserr << "DisplacementControl: " << dUahat  << " " << theDofID << endln;
 //   opserr << "DisplacementControl::newStep() : " << deltaLambdaStep << endln;
    // determine delta U(1) == dU
    (*deltaU) = dUhat;
    (*deltaU) *= dLambda;
    (*deltaUstep) = (*deltaU);

    // update model with delta lambda and delta U
    theModel->incrDisp(*deltaU);    
    theModel->applyLoadDomain(currentLambda);    
    if (theModel->updateDomain() < 0) {
      opserr << "DisplacementControl::newStep - model failed to update for new dU\n";
      return -1;
    }

    numIncrLastStep = 0;

    return 0;
}

int
DisplacementPath::newStep(void)
{
	if (theDofID == -1) {
		opserr << "DisplacementPath::newStep() - domainChanged has not been called\n";
		return -1;
	}

    // get pointers to AnalysisModel and LinearSOE
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();    
    if (theModel == 0 || theLinSOE == 0) {
	opserr << "WARNING DisplacementPath::newStep() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }

	// check theIncrementVector Vector
	if ( theIncrementVector == 0 ) {
		opserr << "DisplacementPath::newStep() - no theIncrementVector associated with object\n";
		return -2;
	}


    // determine increment for this iteration
    if (currentStep < theIncrementVector->Size()) {
		theCurrentIncrement = (*theIncrementVector)(currentStep);
	}
	else {
		theCurrentIncrement = 0.0;
		opserr << "DisplacementPath::newStep() - reach the end of specified load path\n";
		opserr << " - setting theCurrentIncrement = 0.0\n";
	}


    // get the current load factor
    currentLambda = theModel->getCurrentDomainTime();


    // determine dUhat and dUabar
    this->formTangent();
    this->formUnbalance();

	(*deltaUbar) = theLinSOE->getX();
	double dUabar = (*deltaUbar)(theDofID);


    theLinSOE->setB(*phat);

    if (theLinSOE->solve() < 0) {
      opserr << "DisplacementControl::newStep(void) - failed in solver\n";
      return -1;
    }
    
	
	(*deltaUhat) = theLinSOE->getX();
    Vector &dUhat = *deltaUhat;

    double dUahat = dUhat(theDofID);

	//opserr << " newStep( ) " << endln;
    //opserr << " theDofID = " << theDofID << endln;
	//opserr << "dUahat = " << dUahat << endln;

    if (dUahat == 0.0) {
	opserr << "WARNING DisplacementPath::newStep() ";
	opserr << "dUahat is zero -- zero reference displacement at control node DOF\n";
    
	opserr << "currentStep = " << currentStep << endln; // add by zhong
	opserr << " theCurrentIncrement = " << theCurrentIncrement << endln; // zhong

	return -1;
    }
    

    // determine delta lambda(1) == dlambda    
    double dLambda = (theCurrentIncrement-dUabar)/dUahat;

    deltaLambdaStep = dLambda;
    currentLambda += dLambda;
 //   opserr << "DisplacementPath: " << dUahat  << " " << theDofID << endln;
 //   opserr << "DisplacementPath::newStep() : " << deltaLambdaStep << endln;
    // determine delta U(1) == dU
    (*deltaU) = dUhat;
    (*deltaU) *= dLambda;
    (*deltaUstep) = (*deltaU);

    // update model with delta lambda and delta U
    theModel->incrDisp(*deltaU);    
    theModel->applyLoadDomain(currentLambda);    
    if (theModel->updateDomain() < 0) {
      opserr << "DisplacementPath::newStep - model failed to update for new dU\n";
      return -1;
    }

    currentStep++;

    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;
}

int
ArcLengthw::newStep(void)
{
    ofstream factor;
    factor.open("factor.dat",ios::app);
    // get pointers to AnalysisModel and LinearSOE
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();
	if (theModel == 0 || theLinSOE == 0) {
	opserr << "WARNING ArcLengthw::newStep() ";
	opserr << "No AnalysisModel or LinearSOE has been set\n";
	return -1;
    }

    // get the current load factor
    currentLambda = theModel->getCurrentDomainTime();
 
    factor<<"currentLambda"<<endln;
    factor<<currentLambda<<endln;

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

    factor<<"dUhat"<<endln;
    //factor>>dUhat;

    int size = dUhat.Size();
    int i = 0;
    double sum = 0.0;
    int Ji_1 = 0;
    double dLambda = 0.0;
    factor<<"dWibefore"<<endln;
    factor<<dWi<<endln;
    factor<<"*phat"<<endln;
    //factor>>*phat;
    factor<<"dUhat"<<endln;
    //factor>>dUhat;
    factor<<"iFactor"<<endln;
    factor<<iFactor<<endln;
    factor<<"Jd"<<endln;
    factor<<Jd<<endln;
    factor<<"iflag"<<endln;
    factor<<iflag<<endln;
    double dJd = Jd;
    double dJi_1 = 1.0;
    if( iflag == 0 ){
       dWi = ( (*phat) ^ dUhat ) * iFactor * iFactor;
       dLambda = iFactor; 
       iflag = 1;  
    }
    else if( iflag == 1 ){
	   Ji_1 = 10; //theAlgo->getNumIteration();
       dJi_1 = Ji_1;
       dWi = dWi * pow(( dJd / dJi_1 ),0.01);
       dLambda = dWi / ( (*phat)^(dUhat) );
    }
    if( Ji_1 >0){
    factor<<"Jd/Ji-1"<<endln;
    factor<<dJd/dJi_1<<endln;
    }
    factor<<"iflag"<<endln;
    factor<<iflag<<endln;

    factor<<"Ji_1"<<endln;
    factor<<Ji_1<<endln;

    factor<<"dWi"<<endln;
    factor<<dWi<<endln;
    
    deltaLambdaStep = dLambda;
    currentLambda += dLambda;

    (*deltaU) = dUhat;
    (*deltaU) *= dLambda;
    (*deltaUstep) = (*deltaU);

    // update model with delta lambda and delta U
    theModel->incrDisp(*deltaU);    
    theModel->applyLoadDomain(currentLambda);    
    theModel->updateDomain();

    return 0;
}

int 
NewmarkSensitivityIntegrator::formSensitivityRHS(int passedGradNumber)
{
	sensitivityFlag = 1;


	// Set a couple of data members
	gradNumber = passedGradNumber;

	// Get pointer to the SOE
	LinearSOE *theSOE = this->getLinearSOE();


	// Possibly set the independent part of the RHS
	if (assemblyFlag != 0) {
		theSOE->setB(independentRHS);
	}

	// Get the analysis model
	AnalysisModel *theModel = this->getAnalysisModel();



	// Randomness in external load (including randomness in time series)
	// Get domain
	Domain *theDomain = theModel->getDomainPtr();

	// Loop through nodes to zero the unbalaced load
	Node *nodePtr;
	NodeIter &theNodeIter = theDomain->getNodes();
	while ((nodePtr = theNodeIter()) != 0)
	nodePtr->zeroUnbalancedLoad();


	// Loop through load patterns to add external load sensitivity
	LoadPattern *loadPatternPtr;
	LoadPatternIter &thePatterns = theDomain->getLoadPatterns();
	double time;
	while((loadPatternPtr = thePatterns()) != 0) {
		time = theDomain->getCurrentTime();
		loadPatternPtr->applyLoadSensitivity(time);
	}


	// Randomness in element/material contributions
	// Loop through FE elements
	FE_Element *elePtr;
	FE_EleIter &theEles = theModel->getFEs();    
	while((elePtr = theEles()) != 0) {
		theSOE->addB(  elePtr->getResidual(this),  elePtr->getID()  );
	}


	// Loop through DOF groups (IT IS IMPORTANT THAT THIS IS DONE LAST!)
	DOF_Group *dofPtr;
	DOF_GrpIter &theDOFs = theModel->getDOFs();
	while((dofPtr = theDOFs()) != 0) {
		theSOE->addB(  dofPtr->getUnbalance(this),  dofPtr->getID()  );
	}


	// Reset the sensitivity flag
	sensitivityFlag = 0;

	return 0;
}