C++ DOF_Group::getID方法代码示例

int 
IncrementalIntegrator::doMv(const Vector &v, Vector &res) {

  int n = v.Size();
  if (isDiagonal == true) {
    for (int i=0; i<n; i++)
      res[i] = diagMass[i]*v[i];
    return 0;
  }

  res.Zero();

  // loop over the FE_Elements
  FE_Element *elePtr;
  FE_EleIter &theEles = theAnalysisModel->getFEs();    
  while((elePtr = theEles()) != 0) {
    const Vector &b = elePtr->getM_Force(v, 1.0);
    res.Assemble(b, elePtr->getID(), 1.0);
  }

  // loop over the DOF_Groups
  DOF_Group *dofPtr;
  DOF_GrpIter &theDofs = theAnalysisModel->getDOFs();
  while ((dofPtr = theDofs()) != 0) {
    const Vector &a = dofPtr->getM_Force(v, 1.0);      
    res.Assemble(a, dofPtr->getID(), 1.0);
  }
  return 0;
}

void
BandArpackSolver::myMv(int n, double *v, double *result)
{
    Vector x(v, n);
    Vector y(result,n);

    y.Zero();
    AnalysisModel *theAnalysisModel = theSOE->theModel;

    // loop over the FE_Elements
    FE_Element *elePtr;
    FE_EleIter &theEles = theAnalysisModel->getFEs();    
    while((elePtr = theEles()) != 0) {
      const Vector &b = elePtr->getM_Force(x, 1.0);
      y.Assemble(b, elePtr->getID(), 1.0);
    }

    // loop over the DOF_Groups
    DOF_Group *dofPtr;
    DOF_GrpIter &theDofs = theAnalysisModel->getDOFs();
    Integrator *theIntegrator = 0;
    while ((dofPtr = theDofs()) != 0) {
      const Vector &a = dofPtr->getM_Force(x,1.0);      
      y.Assemble(a,dofPtr->getID(),1.0);
    }
}

// void setID(int index, int value);
//	Method to set the correSPonding index of the ID to value.
int
LagrangeSP_FE::setID(void)
{
    int result = 0;

    // first determine the IDs in myID for those DOFs marked
    // as constrained DOFs, this is obtained from the DOF_Group
    // associated with the constrained node
    DOF_Group *theNodesDOFs = theNode->getDOF_GroupPtr();
    if (theNodesDOFs == 0) {
	opserr << "WARNING LagrangeSP_FE::setID(void)";
	opserr << " - no DOF_Group with Constrained Node\n";
	return -1;
    }    

    int restrainedDOF = theSP->getDOF_Number();
    const ID &theNodesID = theNodesDOFs->getID();
    
    if (restrainedDOF < 0 || restrainedDOF >= theNodesID.Size()) {
	opserr << "WARNING LagrangeSP_FE::setID(void)";
	opserr << " - restrained DOF invalid\n";
	return -2;
    }    	
    
    myID(0) = theNodesID(restrainedDOF);
    myID(1) = (theDofGroup->getID())(0);
    
    return result;
}

int
RitzIntegrator::formM()
{
	if (theAnalysisModel == 0 || theSOE == 0) {
		opserr << "WARNING RitzIntegrator::formM -";
		opserr << " no AnalysisModel or EigenSOE has been set\n";
		return -1;
	}

	// the loops to form and add the tangents are broken into two for 
	// efficiency when performing parallel computations

	// loop through the FE_Elements getting them to form the tangent
	// FE_EleIter &theEles1 = theAnalysisModel->getFEs();
	FE_Element *elePtr;

	flagK = 1;
	theSOE->zeroM();

	// while((elePtr = theEles1()) != 0) 
	//     elePtr->formTangent(this);

	// loop through the FE_Elements getting them to add the tangent    
	int result = 0;
	FE_EleIter &theEles2 = theAnalysisModel->getFEs();    
	while((elePtr = theEles2()) != 0) {     
		if (theSOE->addM(elePtr->getTangent(this), elePtr->getID()) < 0) {
			opserr << "WARNING RitzIntegrator::formM -";
			opserr << " failed in addM for ID " << elePtr->getID();	    
			result = -2;
		}
	}

	DOF_Group *dofPtr;
	DOF_GrpIter &theDofs = theAnalysisModel->getDOFs();    
	while((dofPtr = theDofs()) != 0) {
		//   	dofPtr->formTangent(this);
		if (theSOE->addM(dofPtr->getTangent(this),dofPtr->getID()) < 0) {
			opserr << "WARNING RitzIntegrator::formM -";
			opserr << " failed in addM for ID " << dofPtr->getID();	    
			result = -3;
		}
	}

	return result;    
}

int 
PFEMIntegrator::saveSensitivity(const Vector & dVNew,int gradNum,int numGrads)
{
    // Recover sensitivity results from previous step
    int vectorSize = U->Size();
    Vector dUn(vectorSize);
    dVn.resize(vectorSize); dVn.Zero();

    AnalysisModel *myModel = this->getAnalysisModel();
    DOF_GrpIter &theDOFs = myModel->getDOFs();
    DOF_Group *dofPtr;
    while ((dofPtr = theDOFs()) != 0) {
	  
        const ID &id = dofPtr->getID();
        int idSize = id.Size();
        const Vector &dispSens = dofPtr->getDispSensitivity(gradNumber);	
        for (int i=0; i < idSize; i++) {
	    int loc = id(i);
	    if (loc >= 0) {
                dUn(loc) = dispSens(i);		
	    }
        }
	  
        const Vector &velSens = dofPtr->getVelSensitivity(gradNumber);
        for (int i=0; i < idSize; i++) {
	    int loc = id(i);
	    if (loc >= 0) {
                dVn(loc) = velSens(i);
	    }
        }

    }


    // Compute new acceleration and velocity vectors:
    Vector dUNew(vectorSize);
    Vector dANew(vectorSize);

    // dudotdot = 1/dt*dv{n+1} - 1/dt*dvn
    dANew.addVector(0.0, dVNew, c3);
    dANew.addVector(1.0, dVn, -c3);

    // du       = dun + dt*dv{n+1}
    dUNew.addVector(0.0, dVNew, c1);
    dUNew.addVector(1.0, dUn, 1.0);

    // Now we can save vNew, vdotNew and vdotdotNew
    DOF_GrpIter &theDOFGrps = myModel->getDOFs();
    DOF_Group 	*dofPtr1;
    while ( (dofPtr1 = theDOFGrps() ) != 0)  {
        dofPtr1->saveSensitivity(dUNew,dVNew,dANew,gradNum,numGrads);
    }
	
    return 0;
}

int XC::Subdomain::buildMap(void) const
  {
    if(mapBuilt == false)
      {
        // determine the mapping between local dof and subdomain ana dof
        int numDOF = this->getNumDOF();
        if(map == nullptr)
          map = new ID(numDOF);
        if(map->Size() != numDOF)
          {
            delete map;
            map = new ID(numDOF);
          }

        //int numExt = theAnalysis->getNumExternalEqn();
        int numInt = theAnalysis->getNumInternalEqn();

        const ID &theExtNodes = this->getExternalNodes();
        int numExtNodes = theExtNodes.Size();
        int locInMap =0;
        for(int i=0; i<numExtNodes; i++)
          {
            Node *nodePtr= const_cast<Node *>(this->getNode(theExtNodes(i)));
            int numNodeDOF = nodePtr->getNumberDOF();
            DOF_Group *theDOF = nodePtr->getDOF_GroupPtr();
            const ID &theLocalID = theDOF->getID();
            for(int j=0; j<numNodeDOF; j++)
              {
                int locInSubdomainExt = theLocalID(j)-numInt;
                (*map)(locInMap)=locInSubdomainExt;
                locInMap++;
              }
          }
        mapBuilt = true;

        if(mappedVect == nullptr)
          mappedVect = new Vector(numDOF);
        if(mappedVect->Size() != numDOF)
          {
            delete mappedVect;
            mappedVect = new Vector(numDOF);
          }
        if(mappedMatrix == nullptr)
          mappedMatrix = new Matrix(numDOF,numDOF);
        if(mappedMatrix->noRows() != numDOF)
          {
            delete mappedMatrix;
            mappedMatrix = new Matrix(numDOF,numDOF);
          }
      }
    return 0;
  }

int 
IncrementalIntegrator::addModalDampingForce(void)
{
  int res = 0;
  
  if (modalDampingValues == 0)
    return 0;

  int numModes = modalDampingValues->Size();
  const Vector &eigenvalues = theAnalysisModel->getEigenvalues();
  
  if (eigenvalues.Size() < numModes) 
    numModes = eigenvalues.Size();

  Vector dampingForces(theSOE->getNumEqn());

  dampingForces.Zero();

  for (int i=0; i<numModes; i++) {

    DOF_GrpIter &theDOFs1 = theAnalysisModel->getDOFs();
    DOF_Group *dofPtr;
    double beta = 0.0;
    double eigenvalue = eigenvalues(i); // theEigenSOE->getEigenvalue(i+1);
    double wn = 0.;
    if (eigenvalue > 0)
      wn = sqrt(eigenvalue);

    while ((dofPtr = theDOFs1()) != 0) { 
      beta += dofPtr->getDampingBetaFactor(i, (*modalDampingValues)(i), wn);
    }
    
    DOF_GrpIter &theDOFs2 = theAnalysisModel->getDOFs();
    while ((dofPtr = theDOFs2()) != 0) { 
      if (theSOE->addB(dofPtr->getDampingBetaForce(i, beta),dofPtr->getID()) <0) {
	opserr << "WARNING IncrementalIntegrator::failed in dofPtr";
	res = -1;
      }    
    }
  }

  return res;
}

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
NewmarkSensitivityIntegrator::formEleResidual(FE_Element *theEle)
{

	if (sensitivityFlag == 0) {  // NO SENSITIVITY ANALYSIS

		this->Newmark::formEleResidual(theEle);

	}
	else {  // (ASSEMBLE ALL TERMS)

		theEle->zeroResidual();

		// Compute the time-stepping parameters on the form
		// udotdot = a1*ui+1 + a2*ui + a3*udoti + a4*udotdoti
		// udot    = a5*ui+1 + a6*ui + a7*udoti + a8*udotdoti
		// (see p. 166 of Chopra)

		// The constants are:
		// a1 = 1.0/(beta*dt*dt)
		// a2 = -1.0/(beta*dt*dt)
		// a3 = -1.0/beta*dt
		// a4 = 1.0 - 1.0/(2.0*beta)
		// a5 = gamma/(beta*dt)
		// a6 = -gamma/(beta*dt)
		// a7 = 1.0 - gamma/beta
		// a8 = 1.0 - gamma/(2.0*beta)

		// We can make use of the data members c2 and c3 of this class. 
		// As long as disp==true, they are defined as:
		// c2 = gamma/(beta*dt)
		// c3 = 1.0/(beta*dt*dt)

		// So, the constants can be computed as follows:
		if (displ==false) {
			opserr << "ERROR: Newmark::formEleResidual() -- the implemented"
				<< " scheme only works if the displ variable is set to true." << endln;
		}
		double a2 = -c3;
		double a3 = -c2/gamma;
		double a4 = 1.0 - 1.0/(2.0*beta);
		double a6 = -c2;
		double a7 = 1.0 - gamma/beta;
		double dt = gamma/(beta*c2);
		double a8 = dt*(1.0 - gamma/(2.0*beta));


		// Obtain sensitivity vectors from previous step
		int vectorSize = U->Size();
		Vector V(vectorSize);
		Vector Vdot(vectorSize);
		Vector Vdotdot(vectorSize);
		int i, loc;

		AnalysisModel *myModel = this->getAnalysisModel();
		DOF_GrpIter &theDOFs = myModel->getDOFs();
		DOF_Group *dofPtr;
		while ((dofPtr = theDOFs()) != 0) {

			const ID &id = dofPtr->getID();
			int idSize = id.Size();
			const Vector &dispSens = dofPtr->getDispSensitivity(gradNumber);	
			for (i=0; i < idSize; i++) {
				loc = id(i);
				if (loc >= 0) {
					V(loc) = dispSens(i);		
				}
			}

			const Vector &velSens = dofPtr->getVelSensitivity(gradNumber);
			for (i=0; i < idSize; i++) {
				loc = id(i);
				if (loc >= 0) {
					Vdot(loc) = velSens(i);
				}
			}

			const Vector &accelSens = dofPtr->getAccSensitivity(gradNumber);	
			for (i=0; i < idSize; i++) {
				loc = id(i);
				if (loc >= 0) {
					Vdotdot(loc) = accelSens(i);
				}
			}
		}


		// Pre-compute the vectors involving a2, a3, etc.
		//Vector tmp1 = V*a2 + Vdot*a3 + Vdotdot*a4;
		Vector tmp1(vectorSize);
		tmp1.addVector(0.0, V, a2);
		tmp1.addVector(1.0, Vdot, a3);
		tmp1.addVector(1.0, Vdotdot, a4);
		//Vector tmp2 = V*a6 + Vdot*a7 + Vdotdot*a8;
		Vector tmp2(vectorSize);
		tmp2.addVector(0.0, V, a6);
		tmp2.addVector(1.0, Vdot, a7);
		tmp2.addVector(1.0, Vdotdot, a8);

		if (massMatrixMultiplicator == 0)
//.........这里部分代码省略.........

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


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

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


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

int XC::ParallelNumberer::numberDOF(int lastDOF)
  {
    int result = 0;

    // get a pointer to the model & check its not null
    AnalysisModel *theModel = this->getAnalysisModelPtr();
    Domain *theDomain = 0;
    if(theModel) theDomain = theModel->getDomainPtr();
  
    if(theModel == 0 || theDomain == 0)
      {
        std::cerr << "WARNING XC::ParallelNumberer::numberDOF(int) -";
        std::cerr << " - no AnalysisModel.\n";
        return -1;
      }
  
    if(lastDOF != -1)
      {
        std::cerr << "WARNING XC::ParallelNumberer::numberDOF(int lastDOF):";
        std::cerr << " does not use the lastDOF as requested\n";
      }

    Graph &theGraph= theModel->getDOFGroupGraph();

    // if subdomain, collect graph, send it off, get 
    // ID back containing dof tags & start id numbers.
    if(processID != 0)
      {
        CommParameters cp(0,*theChannels[0]);
        const int numVertex = theGraph.getNumVertex();

        /*
        static XC::ID test(2); test(0) = processID; test(1) = 25;
        theChannel->recvID(0, 0, test);
        */

        cp.sendMovable(theGraph,DistributedObj::getDbTagData(),CommMetaData(1));

        // recv iD
        ID theID(2*numVertex);
        cp.receiveID(theID,DistributedObj::getDbTagData(),CommMetaData(2));

        // set vertex numbering based on ID received
        for(int i=0; i<numVertex; i ++)
          {
            const int vertexTag= theID(i);
            int startID= theID(i+numVertex);
            //Vertex *vertexPtr = theGraph.getVertexPtr(vertexTag);
            const int dofTag= vertexTag;
            DOF_Group *dofPtr= theModel->getDOF_GroupPtr(dofTag);
            if(!dofPtr)
              {
                std::cerr << "WARNING ParallelNumberer::numberDOF - ";
                std::cerr << "DOF_Group " << dofTag << "not in XC::AnalysisModel!\n";
                result= -4;
              }
            else
              {
                const ID &theDOFID= dofPtr->getID();
                //std::cerr << "P: " << processID << " dofTag: " << dofTag << " " << "start: " << startID << " " << theDOFID;
                const int idSize= theDOFID.Size();
                for(int j=0; j<idSize; j++)
                  if(theDOFID(j) == -2 || theDOFID(j) == -3) dofPtr->setID(j, startID++);
              }
            //const ID &theDOFID= dofPtr->getID();
          }
        cp.sendID(theID,DistributedObj::getDbTagData(),CommMetaData(2));
      } 
    else
      {
        // if XC::main domain, collect graphs from all subdomains,
        // merge into 1, number this one, send to subdomains the
        // id containing dof tags & start id's.

        // for P0 domain determine original vertex and ref tags
        const int numVertex= theGraph.getNumVertex(); 
        const int numVertexP0= numVertex;

        ID vertexTags(numVertex);
        ID vertexRefs(numVertex);
        Vertex *vertexPtr;
        int loc= 0;
        VertexIter &theVertices= theGraph.getVertices();
        while((vertexPtr= theVertices()) != 0)
          {
            vertexTags[loc]= vertexPtr->getTag();
            vertexRefs[loc]= vertexPtr->getRef();
            loc++;
          }
    
        const int numChannels= theChannels.size();
	std::vector<ID> theSubdomainIDs(numChannels);
        FEM_ObjectBroker theBroker;

        // for each subdomain we receive graph, create an XC::ID (to store
        // subdomain graph to merged graph vertex mapping and the final
        // subdoain graph vertex to startDOF mapping) and finally merge the
        // subdomain graph

        for(int j=0; j<numChannels; j++)
//.........这里部分代码省略.........

int 
DirectIntegrationAnalysis::eigen(int numMode, bool generalized)
{
    if (theAnalysisModel == 0 || theEigenSOE == 0) {
      opserr << "WARNING DirectIntegrationAnalysis::eigen() - no EigenSOE has been set\n";
      return -1;
    }

    int result = 0;
    Domain *the_Domain = this->getDomainPtr();

    result = theAnalysisModel->eigenAnalysis(numMode, generalized);

    int stamp = the_Domain->hasDomainChanged();

    if (stamp != domainStamp) {
      domainStamp = stamp;
      
      result = this->domainChanged();
      
      if (result < 0) {
	opserr << "DirectIntegrationAnalysis::eigen() - domainChanged failed";
	return -1;
      }	
    }

    //
    // zero A and M
    //
    theEigenSOE->zeroA();
    theEigenSOE->zeroM();

    //
    // form K
    //

    FE_EleIter &theEles = theAnalysisModel->getFEs();    
    FE_Element *elePtr;

    while((elePtr = theEles()) != 0) {
      elePtr->zeroTangent();
      elePtr->addKtToTang(1.0);
      if (theEigenSOE->addA(elePtr->getTangent(0), elePtr->getID()) < 0) {
	opserr << "WARNING DirectIntegrationAnalysis::eigen() -";
	opserr << " failed in addA for ID " << elePtr->getID();	    
	result = -2;
      }
    }

    //
    // if generalized is true, form M
    //

    if (generalized == true) {
      int result = 0;
      FE_EleIter &theEles2 = theAnalysisModel->getFEs();    
      while((elePtr = theEles2()) != 0) {     
	elePtr->zeroTangent();
	elePtr->addMtoTang(1.0);
	if (theEigenSOE->addM(elePtr->getTangent(0), elePtr->getID()) < 0) {
	  opserr << "WARNING DirectIntegrationAnalysis::eigen() -";
	  opserr << " failed in addA for ID " << elePtr->getID();	    
	  result = -2;
	}
      }
      
      DOF_Group *dofPtr;
      DOF_GrpIter &theDofs = theAnalysisModel->getDOFs();    
      while((dofPtr = theDofs()) != 0) {
	dofPtr->zeroTangent();
	dofPtr->addMtoTang(1.0);
	if (theEigenSOE->addM(dofPtr->getTangent(0),dofPtr->getID()) < 0) {
	  opserr << "WARNING DirectIntegrationAnalysis::eigen() -";
	  opserr << " failed in addM for ID " << dofPtr->getID();	    
	  result = -3;
	}
      }
    }
    
    // 
    // solve for the eigen values & vectors
    //

    if (theEigenSOE->solve(numMode, generalized) < 0) {
	opserr << "WARNING DirectIntegrationAnalysis::eigen() - EigenSOE failed in solve()\n";
	return -4;
    }

    //
    // now set the eigenvalues and eigenvectors in the model
    //

    theAnalysisModel->setNumEigenvectors(numMode);
    Vector theEigenvalues(numMode);
    for (int i = 1; i <= numMode; i++) {
      theEigenvalues[i-1] = theEigenSOE->getEigenvalue(i);
      theAnalysisModel->setEigenvector(i, theEigenSOE->getEigenvector(i));
    }    
    theAnalysisModel->setEigenvalues(theEigenvalues);
    
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        Ualphadot = new Vector(size);
        Ualphadotdot = new Vector(size);
        Utm1 = new Vector(size);
        Utm2 = new Vector(size);
        scaledDeltaU = 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 ||
            Ualphadot == 0 || Ualphadot->Size() != size ||
            Ualphadotdot == 0 || Ualphadotdot->Size() != size ||
            Utm1 == 0 || Utm1->Size() != size ||
            Utm2 == 0 || Utm2->Size() != size ||
            scaledDeltaU == 0 || scaledDeltaU->Size() != size)  {
            
            opserr << "HHTHSFixedNumIter::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 (Ualphadot != 0)
                delete Ualphadot;
            if (Ualphadotdot != 0)
                delete Ualphadotdot;
            if (Utm1 != 0)
                delete Utm1;
            if (Utm2 != 0)
                delete Utm2;
            if (scaledDeltaU != 0)
                delete scaledDeltaU;
            
            Ut = 0; Utdot = 0; Utdotdot = 0;
            U = 0; Udot = 0; Udotdot = 0;
            Ualpha = 0; Ualphadot = 0; Ualphadotdot = 0;
            Utm1 = 0; Utm2 = 0; scaledDeltaU = 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)  {
                (*Utm1)(loc) = disp(i);
                (*Ut)(loc) = disp(i);
                (*U)(loc) = disp(i);
            }
        }
        
        const Vector &vel = dofPtr->getCommittedVel();
        for (i=0; i < idSize; i++)  {
            int loc = id(i);
            if (loc >= 0)  {
                (*Udot)(loc) = vel(i);
            }
        }
        
        const Vector &accel = dofPtr->getCommittedAccel();	
        for (i=0; i < idSize; i++)  {
            int loc = id(i);
            if (loc >= 0)  {
                (*Udotdot)(loc) = accel(i);
            }
        }
    }
    
    if (polyOrder == 2)
        opserr << "\nWARNING: HHTHSFixedNumIter::domainChanged() - assuming Ut-1 = Ut\n";
    else if (polyOrder == 3)
        opserr << "\nWARNING: HHTHSFixedNumIter::domainChanged() - assuming Ut-2 = Ut-1 = Ut\n";
    
    return 0;
}

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

    DOF_Group *theGroup = theNodePtr->getDOF_GroupPtr();
	if (theGroup == 0) {
	   return 0;
	}
    const ID &theID = theGroup->getID();
    theDofID = theID(theDof);
    
    return 0;
}

int 
CentralDifferenceNoDamping::domainChanged()
{
  AnalysisModel *myModel = this->getAnalysisModel();
  LinearSOE *theLinSOE = this->getLinearSOE();
  const Vector &x = theLinSOE->getX();
  int size = x.Size();
  
  // create the new Vector objects
  if (U == 0 || U->Size() != size) {

    // delete the old
    if (U != 0)
      delete U;
    if (Udot != 0)
      delete Udot;
    if (Udotdot != 0)
      delete Udotdot;

    // create the new
    U = new Vector(size);
    Udot = new Vector(size);
    Udotdot = new Vector(size);

    // cheack we obtained the new
    if (U == 0 || U->Size() != size ||
	Udot == 0 || Udot->Size() != size ||
	Udotdot == 0 || Udotdot->Size() != size) {
      
      opserr << "CentralDifferenceNoDamping::domainChanged - ran out of memory\n";

      // delete the old
      if (U != 0)
	delete U;
      if (Udot != 0)
	delete U;
      if (Udotdot != 0)
	delete Udot;

      U = 0; Udot = 0; Udotdot = 0;
      return -1;
    }
  }        
    
  // now go through and populate U and Udot 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();
    for (i=0; i < idSize; i++)  {
      int loc = id(i);
      if (loc >= 0)  {
	(*Udot)(loc) = vel(i);
      }
    }
  }    
  
  return 0;
}