C++ AnalysisModel类代码示例

int
CentralDifferenceAlternative::newStep(double _deltaT)
{
  updateCount = 0;
  
  deltaT = _deltaT;

  if (deltaT <= 0.0) {
    opserr << "CentralDifference::newStep() - error in variable\n";
    opserr << "dT = " << deltaT << endln;
    return -2;	
  }

  AnalysisModel *theModel = this->getAnalysisModel();
  double time = theModel->getCurrentDomainTime();
  theModel->applyLoadDomain(time);

  return 0;
}

int AlphaOS::commit(void)
{
    AnalysisModel *theModel = this->getAnalysisModel();
    if (theModel == 0)  {
        opserr << "WARNING AlphaOS::commit() - no AnalysisModel set\n";
        return -1;
    }
    
    // set the time to be t+deltaT
    double time = theModel->getCurrentDomainTime();
    time += (1.0-alpha)*deltaT;
    theModel->setCurrentDomainTime(time);
    
    // update the displacements in the elements
    if (updElemDisp == true)
        theModel->updateDomain();
    
    return theModel->commitDomain();
}

int 
TransientIntegrator::formTangent(int statFlag)
{
    int result = 0;
    statusFlag = statFlag;

    LinearSOE *theLinSOE = this->getLinearSOE();
    AnalysisModel *theModel = this->getAnalysisModel();
    if (theLinSOE == 0 || theModel == 0) {
	opserr << "WARNING TransientIntegrator::formTangent() ";
	opserr << "no LinearSOE or AnalysisModel has been set\n";
	return -1;
    }
    
    // the loops to form and add the tangents are broken into two for 
    // efficiency when performing parallel computations
    
    theLinSOE->zeroA();

    // loop through the DOF_Groups and add the unbalance
    DOF_GrpIter &theDOFs = theModel->getDOFs();
    DOF_Group *dofPtr;
    
    while ((dofPtr = theDOFs()) != 0) {
	if (theLinSOE->addA(dofPtr->getTangent(this),dofPtr->getID()) <0) {
	    opserr << "TransientIntegrator::formTangent() - failed to addA:dof\n";
	    result = -1;
	}
    }    

    // loop through the FE_Elements getting them to add the tangent    
    FE_EleIter &theEles2 = theModel->getFEs();    
    FE_Element *elePtr;    
    while((elePtr = theEles2()) != 0)     {
	if (theLinSOE->addA(elePtr->getTangent(this),elePtr->getID()) < 0) {
	    opserr << "TransientIntegrator::formTangent() - failed to addA:ele\n";
	    result = -2;
	}
    }

    return result;
}

int TRBDF2::update(const Vector &deltaU)
{
    AnalysisModel *theModel = this->getAnalysisModel();
    if (theModel == 0)  {
        opserr << "WARNING TRBDF2::update() - no AnalysisModel set\n";
        return -1;
    }	
    
    // check domainChanged() has been called, i.e. Ut will not be zero
    if (Ut == 0)  {
        opserr << "WARNING TRBDF2::update() - domainChange() failed or not called\n";
        return -2;
    }	
    
    // check deltaU is of correct size
    if (deltaU.Size() != U->Size())  {
        opserr << "WARNING TRBDF2::update() - Vectors of incompatible size ";
        opserr << " expecting " << U->Size() << " obtained " << deltaU.Size() << endln;
        return -3;
    }
    
    //  determine the response at t+deltaT
    if (step == 0)  {
        (*U) += deltaU;
        Udot->addVector(1.0, deltaU, c2);
        Udotdot->addVector(1.0, deltaU, c3);
    } else  {
        (*U) += deltaU;
        Udot->addVector(1.0, deltaU, c2);
        Udotdot->addVector(1.0, deltaU, c3);
    }

    // update the response at the DOFs
    theModel->setResponse(*U,*Udot,*Udotdot);
    if (theModel->updateDomain() < 0)  {
      opserr << "TRBDF2::update() - failed to update the domain\n";
      return -4;
    }
    
    return 0;
}    

int PFEMIntegrator::newStep(double deltaT)
{

    if (deltaT <= 0.0)  {
        opserr << "PFEMIntegrator::newStep() - error in variable\n";
        opserr << "dT = " << deltaT << endln;
        return -2;	
    }

    // get a pointer to the AnalysisModel and Domain
    AnalysisModel *theModel = this->getAnalysisModel();
    if(theModel == 0) {
        opserr << "Analysis model has not been linked - PFEMIntegrator::newStep()\n";
        return -1;
    }
    Domain* theDomain = theModel->getDomainPtr();
    if(theDomain == 0) {
        opserr<<"WARNING: no domain is set for the model";
        opserr<<" -- PFEMIntegrator::newStep()\n";
        return -1;
    }
    
    // set the constants
    c1 = deltaT;
    c2 = 1.0;
    c3 = 1.0/deltaT;

    c4 = deltaT*deltaT;
    c5 = deltaT;
    c6 = 1.0;

    // check if domainchange() is called
    if (U == 0)  {
        opserr << "PFEMIntegrator::newStep() - domainChange() failed or hasn't been called\n";
        return -3;	
    }
    
    // set response at t to be that at t+deltaT of previous step
    (*Ut) = *U;        
    (*Utdot) = *Udot;  
    (*Utdotdot) = *Udotdot;
    
    // determinte new disps and accels
    U->addVector(1.0, *Utdot, deltaT);
    Udotdot->Zero();

    // set states
    theModel->setDisp(*U);
    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 << "PFEMIntegrator::newStep() - failed to update the domain\n";
        return -4;
    }
    
    return 0;
}

int CollocationHSIncrReduct::update(const Vector &deltaU)
{
    AnalysisModel *theModel = this->getAnalysisModel();
    if (theModel == 0)  {
        opserr << "WARNING CollocationHSIncrReduct::update() - no AnalysisModel set\n";
        return -1;
    }
    
    // check domainChanged() has been called, i.e. Ut will not be zero
    if (Ut == 0)  {
        opserr << "WARNING CollocationHSIncrReduct::update() - domainChange() failed or not called\n";
        return -2;
    }
    
    // check deltaU is of correct size
    if (deltaU.Size() != U->Size())  {
        opserr << "WARNING CollocationHSIncrReduct::update() - Vectors of incompatible size ";
        opserr << " expecting " << U->Size() << " obtained " << deltaU.Size() << endln;
        return -3;
    }
    
    // get scaled increment
    (*scaledDeltaU) = reduct*deltaU;
    
    // determine the response at t+theta*deltaT
    U->addVector(1.0, *scaledDeltaU, c1);
    
    Udot->addVector(1.0, *scaledDeltaU, c2);
    
    Udotdot->addVector(1.0, *scaledDeltaU, c3);
    
    // update the response at the DOFs
    theModel->setResponse(*U, *Udot, *Udotdot);
    if (theModel->updateDomain() < 0)  {
        opserr << "CollocationHSIncrReduct::update() - failed to update the domain\n";
        return -4;
    }
    
    return 0;
}

int NewmarkHSFixedNumIter::commit(void)
{
    AnalysisModel *theModel = this->getAnalysisModel();
    if (theModel == 0)  {
        opserr << "WARNING NewmarkHSFixedNumIter::commit() - no AnalysisModel set\n";
        return -1;
    }

    LinearSOE *theSOE = this->getLinearSOE();
    if (theSOE == 0)  {
        opserr << "WARNING NewmarkHSFixedNumIter::commit() - no LinearSOE set\n";
        return -2;
    }

    if (this->formTangent(statusFlag) < 0)  {
        opserr << "WARNING NewmarkHSFixedNumIter::commit() - "
               << "the Integrator failed in formTangent()\n";
        return -3;
    }

    if (theSOE->solve() < 0)  {
        opserr << "WARNING NewmarkHSFixedNumIter::commit() - "
               << "the LinearSysOfEqn failed in solve()\n";
        return -4;
    }
    const Vector &deltaU = theSOE->getX();

    //  determine the response at t+deltaT
    U->addVector(1.0, deltaU, c1);

    Udot->addVector(1.0, deltaU, c2);

    Udotdot->addVector(1.0, deltaU, c3);

    // update the response at the DOFs
    theModel->setResponse(*U,*Udot,*Udotdot);

    return theModel->commitDomain();
}

int XC::Collocation::commit(void)
{
    
    AnalysisModel *theModel = this->getAnalysisModelPtr();
    if (theModel == 0)  {
        std::cerr << "WARNING XC::Collocation::commit() - no XC::AnalysisModel set\n";
        return -1;
    }	  
        
    // determine response quantities at t+deltaT
    U.getDotDot().addVector(1.0/theta, Ut.getDotDot(), (theta-1.0)/theta);
    
    (U.getDot()) = Ut.getDot();
    double a1 = deltaT*(1.0 - gamma);
    double a2 = deltaT*gamma;
    U.getDot().addVector(1.0, Ut.getDotDot(), a1);
    U.getDot().addVector(1.0, U.getDotDot(), a2);
    
    U.get()= Ut.get();
    U.get().addVector(1.0, Ut.getDot(), deltaT);
    double a3 = deltaT*deltaT*(0.5 - beta);
    double a4 = deltaT*deltaT*beta;
    U.get().addVector(1.0, Ut.getDotDot(), a3);
    U.get().addVector(1.0, U.getDotDot(), a4);

    // update the response at the DOFs
    theModel->setResponse(U.get(),U.getDot(),U.getDotDot());        
//    if (theModel->updateDomain() < 0)  {
//        std::cerr << "XC::Collocation::commit() - failed to update the domain\n";
//        return -4;
//    }
    
    // set the time to be t+delta t
    double time= getCurrentModelTime();
    time += (1.0-theta)*deltaT;
    setCurrentModelTime(time);
    return commitModel();
  }

int XC::HHTHybridSimulation::commit(void)
  {
    AnalysisModel *theModel = this->getAnalysisModelPtr();
    if(theModel == 0)
      {
        std::cerr << "WARNING XC::HHTHybridSimulation::commit() - no XC::AnalysisModel set\n";
        return -1;
      }
    
    // update the response at the DOFs
    theModel->setResponse(U.get(),U.getDot(),U.getDotDot());
//    if (theModel->updateDomain() < 0)  {
//        std::cerr << "XC::HHTHybridSimulation::commit() - failed to update the domain\n";
//        return -4;
//    }
    
    // set the time to be t+deltaT
    double time = getCurrentModelTime();
    time+= (1.0-alphaF)*deltaT;
    setCurrentModelTime(time);

    return commitModel();
  }

int XC::Collocation::update(const Vector &deltaU)
  {
    AnalysisModel *theModel = this->getAnalysisModelPtr();
    if (theModel == 0)  {
        std::cerr << "WARNING XC::Collocation::update() - no XC::AnalysisModel set\n";
        return -1;
    }	
    
    // check domainChanged() has been called, i.e. Ut will not be zero
    if (Ut.get().Size() == 0)  {
        std::cerr << "WARNING XC::Collocation::update() - domainChange() failed or not called\n";
        return -2;
    }	
    
    // check deltaU is of correct size
    if (deltaU.Size() != U.get().Size())  {
        std::cerr << "WARNING XC::Collocation::update() - Vectors of incompatible size ";
        std::cerr << " expecting " << U.get().Size() << " obtained " << deltaU.Size() << std::endl;
        return -3;
    }
    
    // determine the response at t+theta*deltaT
    (U.get()) += deltaU;

    U.getDot().addVector(1.0, deltaU, c2);
    
    U.getDotDot().addVector(1.0, deltaU, c3);
    
    // update the response at the DOFs
    theModel->setResponse(U.get(),U.getDot(),U.getDotDot());        
    if(updateModel() < 0)
      {
        std::cerr << "XC::Collocation::update() - failed to update the domain\n";
        return -4;
      }
    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;
}

//! @brief Constructor.
//!
//! The constructor is responsible for constructing the graph given {\em
//! theModel}. It creates the vertices of the graph, one for every
//! equation  (each DOF that has not been constrained out by the
//! constraint handler)  in the model and adds all edges based on the
//! FE\_Element connectivity. For this reason the model must be fully
//! populated with the DOF\_Group and FE\_Element objects before
// the constructor is called.
XC::DOF_GroupGraph::DOF_GroupGraph(const AnalysisModel &theModel)
  :ModelGraph(theModel.getNumDOF_Groups()+START_VERTEX_NUM,theModel)
  {
    assert(myModel);
    const int numVertex= myModel->getNumDOF_Groups();

    if(numVertex>0)
      {
	
        const DOF_Group *dofGroupPtr= nullptr;

        // now create the vertices with a reference equal to the DOF_Group number.
        // and a tag which ranges from 0 through numVertex-1
        DOF_GrpConstIter &dofIter2 = myModel->getConstDOFs();
        while((dofGroupPtr = dofIter2()) != 0)
          {
	    const int DOF_GroupTag = dofGroupPtr->getTag();
            const int DOF_GroupNodeTag = dofGroupPtr->getNodeTag();
	    const int numDOF = dofGroupPtr->getNumFreeDOF();
            Vertex vrt(DOF_GroupTag, DOF_GroupNodeTag, 0, numDOF);
            this->addVertex(vrt);
          }


        // now add the edges, by looping over the Elements, getting their
        // IDs and adding edges between DOFs for equation numbers >= START_EQN_NUM
    
        const FE_Element *elePtr= nullptr;
        FE_EleConstIter &eleIter = myModel->getConstFEs();
        while((elePtr = eleIter()) != 0)
          {
	    const ID &id= elePtr->getDOFtags();
	    const int size = id.Size();
            for(int i=0; i<size; i++)
              {
	        int dof1 = id(i);
	        for(int j=0; j<size; j++)
                  if(i != j)
                    {
		      const int dof2 = id(j);
		      this->addEdge(dof1,dof2);
		    }
	      }
          }
      }
  }

int GeneralizedAlpha::commit(void)
{
    AnalysisModel *theModel = this->getAnalysisModel();
    if (theModel == 0)  {
        opserr << "WARNING GeneralizedAlpha::commit() - no AnalysisModel set\n";
        return -1;
    }	  
    
    // update the response at the DOFs
    theModel->setResponse(*U,*Udot,*Udotdot);
    if (theModel->updateDomain() < 0)  {
        opserr << "GeneralizedAlpha::commit() - failed to update the domain\n";
        return -4;
    }
    
    // set the time to be t+deltaT
    double time = theModel->getCurrentDomainTime();
    time += (1.0-alphaF)*deltaT;
    theModel->setCurrentDomainTime(time);

    return theModel->commitDomain();
}

//! @brief Handle the constraints.
//! 
//! Determines the number of FE\_Elements and DOF\_Groups needed from the
//! Domain (a one to one mapping between Elements and FE\_Elements and
//! Nodes and DOF\_Groups) Creates two arrays of pointers to store the
//! FE\_elements and DOF\_Groups, returning a warning message and a \f$-2\f$
//! or \f$-3\f$ if not enough memory is available for these arrays. Then the
//! object will iterate through the Nodes of the Domain, creating a
//! DOF\_Group for each node and setting the initial id for each dof to
//! \f$-2\f$ if no SFreedom\_Constraint exists for the dof, or \f$-1\f$ if a
//! SFreedom\_Constraint exists or \f$-3\f$ if the node identifier is in {\em
//! nodesToBeNumberedLast}. The object then iterates through the Elements
//! of the Domain creating a FE\_Element for each Element, if the Element
//! is a Subdomain setFE\_ElementPtr() is invoked on the Subdomain
//! with the new FE\_Element as the argument. If not enough memory is
//! available for any DOF\_Group or FE\_element a warning message is
//! printed and a \f$-4\f$ or \f$-5\f$ is returned. If any MFreedom\_Constraint
//! objects exist in the Domain a warning message is printed and \f$-6\f$ is
//! returned. If all is successful, the method returns the number of
//! degrees-of-freedom associated with the DOF\_Groups in {\em
//! nodesToBeNumberedLast}.
int XC::PlainHandler::handle(const ID *nodesLast)
  {
    // first check links exist to a Domain and an AnalysisModel object
    Domain *theDomain = this->getDomainPtr();
    AnalysisModel *theModel = this->getAnalysisModelPtr();
    Integrator *theIntegrator = this->getIntegratorPtr();    
    
    if((!theDomain) || (!theModel) || (!theIntegrator))
      {
        std::cerr << getClassName() << "::" << __FUNCTION__
                  << "; domain, model or integrator was not set.\n";
        return -1;
      }

    // initialse the DOF_Groups and add them to the AnalysisModel.
    //    : must of course set the initial IDs
    NodeIter &theNod= theDomain->getNodes();
    Node *nodPtr= nullptr;
    SFreedom_Constraint *spPtr= nullptr;
    DOF_Group *dofPtr= nullptr;

    int numDOF = 0;
    int count3 = 0;
    int countDOF =0;
    while((nodPtr = theNod()) != nullptr)
      {
        dofPtr= theModel->createDOF_Group(numDOF++, nodPtr);
        // initially set all the ID value to -2
        countDOF+= dofPtr->inicID(-2);

        // loop through the SFreedom_Constraints to see if any of the
        // DOFs are constrained, if so set initial XC::ID value to -1
        int nodeID = nodPtr->getTag();
        SFreedom_ConstraintIter &theSPs = theDomain->getConstraints().getDomainAndLoadPatternSPs();
        while((spPtr = theSPs()) != 0)
            if(spPtr->getNodeTag() == nodeID)
              {
                if(spPtr->isHomogeneous() == false)
                  std::cerr << getClassName() << "::" << __FUNCTION__
                            << ";  non-homogeneos constraint"
                            << " for node " << spPtr->getNodeTag()
                            << " homo assumed\n";
                const ID &id = dofPtr->getID();
                int dof = spPtr->getDOF_Number();                
                if(id(dof) == -2)
                  {
                        dofPtr->setID(spPtr->getDOF_Number(),-1);
                        countDOF--;        
                  }
                else
                  std::cerr << getClassName() << "::" << __FUNCTION__
                            << "; multiple single pointconstraints at DOF "
                            << dof << " for node " << spPtr->getNodeTag()
                            << std::endl;
              }

        // loop through the MFreedom_Constraints to see if any of the
        // DOFs are constrained, note constraint matrix must be diagonal
        // with 1's on the diagonal
        MFreedom_ConstraintIter &theMPs = theDomain->getConstraints().getMPs();
        MFreedom_Constraint *mpPtr;
        while((mpPtr = theMPs()) != 0)
          {
            if(mpPtr->getNodeConstrained() == nodeID)
              {
                if(mpPtr->isTimeVarying() == true)
                  std::cerr << getClassName() << "::" << __FUNCTION__
                            << ";  time-varying constraint"
                            << " for node " << nodeID
                            << " non-varying assumed\n";
                const Matrix &C = mpPtr->getConstraint();
                int numRows = C.noRows();
                int numCols = C.noCols();
                if(numRows != numCols)
                  std::cerr << getClassName() << "::" << __FUNCTION__
                            << " constraint matrix not diagonal,"
                            << " ignoring constraint for node "
                            << nodeID << std::endl;
                else
//.........这里部分代码省略.........

int 
HHT1::domainChanged()
{
  AnalysisModel *myModel = this->getAnalysisModel();
  LinearSOE *theLinSOE = this->getLinearSOE();
  const Vector &x = theLinSOE->getX();
  int size = x.Size();

  // if damping factors exist set them in the ele & node of the domain
  if (alphaM != 0.0 || betaK != 0.0 || betaKi != 0.0 || betaKc != 0.0)
    myModel->setRayleighDampingFactors(alphaM, betaK, betaKi, betaKc);

  // create the new Vector objects
  if (Ut == 0 || Ut->Size() != size) {

    // delete the old
    if (Ut != 0)
      delete Ut;
    if (Utdot != 0)
      delete Utdot;
    if (Utdotdot != 0)
      delete Utdotdot;
    if (U != 0)
      delete U;
    if (Udot != 0)
      delete Udot;
    if (Udotdot != 0)
      delete Udotdot;
    if (Ualpha != 0)
      delete Ualpha;
    if (Udotalpha != 0)
      delete Udotalpha;
    
    // create the new
    Ut = new Vector(size);
    Utdot = new Vector(size);
    Utdotdot = new Vector(size);
    U = new Vector(size);
    Udot = new Vector(size);
    Udotdot = new Vector(size);
    Ualpha = new Vector(size);
    Udotalpha = new Vector(size);

    // check we obtained the new
    if (Ut == 0 || Ut->Size() != size ||
	Utdot == 0 || Utdot->Size() != size ||
	Utdotdot == 0 || Utdotdot->Size() != size ||
	U == 0 || U->Size() != size ||
	Udot == 0 || Udot->Size() != size ||
	Udotdot == 0 || Udotdot->Size() != size ||
	Ualpha == 0 || Ualpha->Size() != size ||
	Udotalpha == 0 || Udotalpha->Size() != size) {
  
      opserr << "HHT1::domainChanged - ran out of memory\n";

      // delete the old
      if (Ut != 0)
	delete Ut;
      if (Utdot != 0)
	delete Utdot;
      if (Utdotdot != 0)
	delete Utdotdot;
      if (U != 0)
	delete U;
      if (Udot != 0)
	delete Udot;
      if (Udotdot != 0)
	delete Udotdot;
    if (Ualpha != 0)
      delete Ualpha;
    if (Udotalpha != 0)
      delete Udotalpha;

      Ut = 0; Utdot = 0; Utdotdot = 0;
      U = 0; Udot = 0; Udotdot = 0; Udotalpha=0; Ualpha =0;
      return -1;
    }
  }        
    
  // now go through and populate U, Udot and Udotdot by iterating through
  // the DOF_Groups and getting the last committed velocity and accel

  DOF_GrpIter &theDOFs = myModel->getDOFs();
  DOF_Group *dofPtr;

  while ((dofPtr = theDOFs()) != 0) {
    const ID &id = dofPtr->getID();
    int idSize = id.Size();


	int i;
    const Vector &disp = dofPtr->getCommittedDisp();	
    for (i=0; i < idSize; i++) {
      int loc = id(i);
      if (loc >= 0) {
 	(*U)(loc) = disp(i);		
      }
    }

    const Vector &vel = dofPtr->getCommittedVel();
//.........这里部分代码省略.........