C++ AnalysisModel::getDOFs方法代码示例

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

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

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 CollocationHSIncrReduct::domainChanged()
{
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theLinSOE = this->getLinearSOE();
    const Vector &x = theLinSOE->getX();
    int size = x.Size();
    
    // 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 (scaledDeltaU != 0)
            delete scaledDeltaU;
        
        // 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);
        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 ||
            scaledDeltaU == 0 || scaledDeltaU->Size() != size)  {
            
            opserr << "CollocationHSIncrReduct::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 (scaledDeltaU != 0)
                delete scaledDeltaU;
            
            Ut = 0; Utdot = 0; Utdotdot = 0;
            U = 0; Udot = 0; Udotdot = 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 = theModel->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);
            }
        }
        
        const Vector &accel = dofPtr->getCommittedAccel();
        for (i=0; i < idSize; i++)  {
            int loc = id(i);
            if (loc >= 0)  {
                (*Udotdot)(loc) = accel(i);
//.........这里部分代码省略.........

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

//.........这里部分代码省略.........
        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
PlainNumberer::numberDOF(int lastDOF)
{
    int eqnNumber = 0; // start equation number = 0

    // get a pointer to the model & check its not null
    AnalysisModel *theModel = this->getAnalysisModelPtr();
    Domain *theDomain = 0;
    if (theModel != 0) theDomain = theModel->getDomainPtr();

    if (theModel == 0 || theDomain == 0) {
	opserr << "WARNING PlainNumberer::numberDOF(int) -";
	opserr << " - no AnalysisModel - has setLinks() been invoked?\n";
	return -1;
    }
    
    if (lastDOF != -1) {
	opserr << "WARNING PlainNumberer::numberDOF(int lastDOF):";
	opserr << " does not use the lastDOF as requested\n";
    }
    
    // iterate throgh  the DOFs first time setting -2 values
    DOF_GrpIter &theDOFs = theModel->getDOFs();
    DOF_Group *dofPtr;
    
    while ((dofPtr = theDOFs()) != 0) {
	const ID &theID = dofPtr->getID();
	for (int i=0; i<theID.Size(); i++)
	    if (theID(i) == -2) 
	      dofPtr->setID(i,eqnNumber++);
    }

    // iterate throgh  the DOFs second time setting -3 values
    DOF_GrpIter &moreDOFs = theModel->getDOFs();
    
    while ((dofPtr = moreDOFs()) != 0) {
	const ID &theID = dofPtr->getID();
	for (int i=0; i<theID.Size(); i++)
	    if (theID(i) == -3) dofPtr->setID(i,eqnNumber++);
    }

    // iterate through the DOFs one last time setting any -4 values
    DOF_GrpIter &tDOFs = theModel->getDOFs();
    while ((dofPtr = tDOFs()) != 0) {
    	const ID &theID = dofPtr->getID();
    	int have4s = 0;
	for (int i=0; i<theID.Size(); i++)
	    if (theID(i) == -4) have4s = 1;

	if (have4s == 1) {
		int nodeID = dofPtr->getNodeTag();
		// loop through the MP_Constraints to see if any of the
		// DOFs are constrained, note constraint matrix must be diagonal
		// with 1's on the diagonal
		MP_ConstraintIter &theMPs = theDomain->getMPs();
		MP_Constraint *mpPtr;
		while ((mpPtr = theMPs()) != 0 ) {
			// note keep looping over all in case multiple constraints
			// are used to constrain a node -- can't assume intelli user
	    		if (mpPtr->getNodeConstrained() == nodeID) {
	    			int nodeRetained = mpPtr->getNodeRetained();
	    			Node *nodeRetainedPtr = theDomain->getNode(nodeRetained);
	    			DOF_Group *retainedDOF = nodeRetainedPtr->getDOF_GroupPtr();
	    			const ID&retainedDOFIDs = retainedDOF->getID();
	    			const ID&constrainedDOFs = mpPtr->getConstrainedDOFs();
	    			const ID&retainedDOFs = mpPtr->getRetainedDOFs();
	    			for (int i=0; i<constrainedDOFs.Size(); i++) {
	    				int dofC = constrainedDOFs(i);
	    				int dofR = retainedDOFs(i);
	    				int dofID = retainedDOFIDs(dofR);
	    				dofPtr->setID(dofC, dofID);
	    			}
	    		}
		}		
	}	
    }

    eqnNumber--;
    int numEqn = eqnNumber - START_EQN_NUMBER +1;
	
    // iterate through the FE_Element getting them to set their IDs
    FE_EleIter &theEle = theModel->getFEs();
    FE_Element *elePtr;
    while ((elePtr = theEle()) != 0)
	elePtr->setID();

    // set the numOfEquation in the Model
    theModel->setNumEqn(numEqn);

    return numEqn;
}

void
ArpackSolver::myMv(int n, double *v, double *result)
{
  Vector x(v, n);
  Vector y(result,n);
    
  bool mDiagonal = theArpackSOE->mDiagonal;

  if (mDiagonal == true) {

    int Msize = theArpackSOE->Msize;
    double *M = theArpackSOE->M;

    /* for output
    DataFileStream dataStream("M.txt");
    dataStream.open();
    for (int i=0; i<n; i++)
      dataStream << M[i] << endln;
    dataStream.close();
    */

    if (n <= Msize) {
      for (int i=0; i<n; i++)
	result[i] = M[i]*v[i];
    } else {
      opserr << "ArpackSolver::myMv() n > Msize!\n";
      return;
    }

  } else {

    y.Zero();

    AnalysisModel *theAnalysisModel = theArpackSOE->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();
    while ((dofPtr = theDofs()) != 0) {
      const Vector &a = dofPtr->getM_Force(x,1.0);      
      y.Assemble(a, dofPtr->getID(), 1.0);
    }
  }

  // if paallel we have to merge the results
  int processID = theArpackSOE->processID;
  if (processID != -1) {
    Channel **theChannels = theArpackSOE->theChannels;
    int numChannels = theArpackSOE->numChannels;
    if (processID != 0) {
      theChannels[0]->sendVector(0, 0, y);
      theChannels[0]->recvVector(0, 0, y);
    } else {
      Vector other(workArea, n);
      // recv contribution from remote & add
      for (int i=0; i<numChannels; i++) {
	theChannels[i]->recvVector(0,0,other);
	y += other;
      }
      // send result back
      for (int i=0; i<numChannels; i++) {
	theChannels[i]->sendVector(0,0,y);
      }
    }
  }
}

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

int CentralDifference::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 element & 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)  {
        
        if (Utm1 != 0)
            delete Utm1;
        if (Ut != 0)
            delete Ut;
        if (Utdot != 0)
            delete Utdot;
        if (Utdotdot != 0)
            delete Utdotdot;
        if (Udot != 0)
            delete Udot;
        if (Udotdot != 0)
            delete Udotdot;
        
        // create the new
        Utm1 = new Vector(size);
        Ut = new Vector(size);
        Utdot = new Vector(size);
        Utdotdot = new Vector(size);
        Udot = new Vector(size);
        Udotdot = new Vector(size);
        
        // check we obtained the new
        if (Utm1 == 0 || Utm1->Size() != size ||
            Ut == 0 || Ut->Size() != size ||
            Utdot == 0 || Utdot->Size() != size ||
            Utdotdot == 0 || Utdotdot->Size() != size ||
            Udot == 0 || Udot->Size() != size ||
            Udotdot == 0 || Udotdot->Size() != size)  {
            
            opserr << "CentralDifference::domainChanged - ran out of memory\n";
            
            // delete the old
            if (Utm1 != 0)
                delete Utm1;
            if (Ut != 0)
                delete Ut;
            if (Utdot != 0)
                delete Utdot;
            if (Utdotdot != 0)
                delete Utdotdot;
            if (Udot != 0)
                delete Udot;
            if (Udotdot != 0)
                delete Udotdot;
            
            Utm1 = 0;
            Ut = 0; Utdot = 0; Utdotdot = 0;
            Udot = 0; Udotdot = 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);		
            }
        }
        
        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);
            }
        }
//.........这里部分代码省略.........

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

//.........这里部分代码省略.........
        Ualpha = new Vector(size);
        Ualphadot = new Vector(size);
        Ualphadotdot = new Vector(size);
        Utdothat = new Vector(size);
        
        // check we obtained the new
        if (alpha1 == 0 || alpha1->noRows() != size || alpha1->noCols() != size ||
            alpha3 == 0 || alpha3->noRows() != size || alpha3->noCols() != size ||
            Mhat == 0 || Mhat->noRows() != size || Mhat->noCols() != size ||
            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 ||
            Utdothat == 0 || Utdothat->Size() != size)  {
            
            opserr << "WARNING KRAlphaExplicit::domainChanged() - ";
            opserr << "ran out of memory\n";
            
            // delete the old
            if (alpha1 != 0)
                delete alpha1;
            if (alpha3 != 0)
                delete alpha3;
            if (Mhat != 0)
                delete Mhat;
            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 (Utdothat != 0)
                delete Utdothat;
            
            alpha1 = 0; alpha3 = 0; Mhat = 0;
            Ut = 0; Utdot = 0; Utdotdot = 0;
            U = 0; Udot = 0; Udotdot = 0;
            Ualpha = 0; Ualphadot = 0; Ualphadotdot = 0;
            Utdothat = 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();
        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);
            }
        }
    }
    
    // recalculate integration parameter matrices b/c domain changed
    initAlphaMatrices = 1;
    
    return 0;
}