C++ Domain::getNode方法代码示例

PenaltyMP_FE::PenaltyMP_FE(int tag, Domain &theDomain, 
			   MP_Constraint &TheMP, double Alpha)
:FE_Element(tag, 2,(TheMP.getConstrainedDOFs()).Size()+
 (TheMP.getRetainedDOFs()).Size()),
 theMP(&TheMP), theConstrainedNode(0) , theRetainedNode(0),
 tang(0), resid(0), C(0), alpha(Alpha)
{
    
    int size;
    const ID &id1 = theMP->getConstrainedDOFs();
    size = id1.Size();
    const ID &id2 = theMP->getRetainedDOFs();    
    size += id2.Size();

    tang = new Matrix(size,size);
    resid = new Vector(size);
    C = new Matrix(id1.Size(),size);

    if (tang == 0 || resid == 0 || C == 0 ||
	tang->noCols() != size || C->noCols() != size || 
	resid->Size() != size) {
	opserr << "FATAL PenaltyMP_FE::PenaltyMP_FE() - out of memory\n";
	exit(-1);
    }
	    
    theRetainedNode = theDomain.getNode(theMP->getNodeRetained());    
    theConstrainedNode = theDomain.getNode(theMP->getNodeConstrained());

    if (theRetainedNode == 0 || theConstrainedNode == 0) {
	opserr << "FATAL PenaltyMP_FE::PenaltyMP_FE() - Constrained or Retained";
	opserr << " Node does not exist in Domain\n";
	opserr << theMP->getNodeRetained() << " " << theMP->getNodeConstrained() << endln;
	exit(-1);
    }	


    // set up the dof groups tags
    DOF_Group *dofGrpPtr = 0;
    dofGrpPtr = theRetainedNode->getDOF_GroupPtr();
    if (dofGrpPtr != 0) 
	myDOF_Groups(0) = dofGrpPtr->getTag();	    
    else 
	opserr << "WARNING PenaltyMP_FE::PenaltyMP_FE() - node no Group yet?\n"; 
    dofGrpPtr = theConstrainedNode->getDOF_GroupPtr();
    if (dofGrpPtr != 0) 
	myDOF_Groups(1) = dofGrpPtr->getTag();	        
    else
	opserr << "WARNING PenaltyMP_FE::PenaltyMP_FE() - node no Group yet?\n"; 
    
    
    if (theMP->isTimeVarying() == false) {
	this->determineTangent();
	// we can free up the space taken by C as it is no longer needed
	if (C != 0)
	    delete C;
	C = 0;
    }
}

int
MeshRegion::setRayleighDampingFactors(double alpham, double betak, double betak0, double betakc)
{
  alphaM = alpham;
  betaK  = betak;
  betaK0 = betak0;
  betaKc = betakc;

  // now set the damping factors at the nodes & elements
  Domain *theDomain = this->getDomain();
  if (theDomain == 0) {
    opserr << "MeshRegion::setRayleighDampingFactors() - no domain yet set\n";
    return -1;
  }
  if (theElements != 0) {
    for (int i=0; i<theElements->Size(); i++) {
      int eleTag = (*theElements)(i);
      Element *theEle = theDomain->getElement(eleTag);
      if (theEle != 0)
	theEle->setRayleighDampingFactors(alphaM, betaK, betaK0, betaKc);
    }
  }

  if (theNodes != 0) {
    for (int i=0; i<theNodes->Size(); i++) {
      int nodTag = (*theNodes)(i);
      Node *theNode = theDomain->getNode(nodTag);
      if (theNode != 0)
	theNode->setRayleighDampingFactor(alphaM);
    }
  }

  return 0;
}

int ExpControlSP::applyConstraint(double loadFactor)
{
    // on first call
    if (theNode == 0 || theNodeResponse == 0)  {
        Domain *theDomain = this->getDomain();
        
        theNode = theDomain->getNode(nodeTag);
        if (theNode == 0)
            return -1;
        
        theNodeResponse = new Vector(theNode->getNumberDOF());
        if (theNodeResponse == 0)
            return -2;
    }
    
    // set the responses at the node
    // disp response is the responsibility of constraint handler
    //*theNodeResponse = theNode->getTrialDisp();
    //(*theNodeResponse)(dofNumber) = dispFact * (*ctrlDisp);
    //theNode->setTrialDisp(*theNodeResponse);
    
    if (ctrlVel != 0 && velFact != 0.0)  {
        *theNodeResponse = theNode->getTrialVel();
        (*theNodeResponse)(dofNumber) = velFact * (*ctrlVel);
        theNode->setTrialVel(*theNodeResponse);
    }
    
    if (ctrlAccel != 0 && accelFact != 0.0)  {
        *theNodeResponse = theNode->getTrialAccel();
        (*theNodeResponse)(dofNumber) = accelFact * (*ctrlAccel);
        theNode->setTrialAccel(*theNodeResponse);
    }
    
    return 0;
}

Node*
Pressure_Constraint::getPressureNode()
{
    if(pval != 0) return 0;

    Domain* theDomain = this->getDomain();
    if(theDomain == 0) {
        opserr<<"WARNING: domain has not been set";
        opserr<<" -- Pressure_Constraint::getPressureNode\n";
        return 0;
    }
    return theDomain->getNode(pTag);
}

LagrangeSP_FE::LagrangeSP_FE(int tag, Domain &theDomain, SP_Constraint &TheSP,
			     DOF_Group &theGroup, double Alpha)
:FE_Element(tag, 2,2),
 alpha(Alpha), tang(0), resid(0), theSP(&TheSP), theDofGroup(&theGroup)
{
    // create a Matrix and a Vector for the tangent and residual
    tang = new Matrix(2,2);
    resid = new Vector(2);
    if ((tang == 0) || (tang->noCols() == 0) || (resid == 0) ||
	(resid->Size() == 0)) {
	opserr << "WARNING LagrangeSP_FE::LagrangeSP_FE()";
	opserr << "- ran out of memory\n";
	exit(-1);
    }

    // zero the Matrix and Vector
    resid->Zero();
    tang->Zero();

    theNode = theDomain.getNode(theSP->getNodeTag());    
    if (theNode == 0) {
	opserr << "WARNING LagrangeSP_FE::LagrangeSP_FE()";
	opserr << "- no asscoiated Node\n";
	exit(-1);
    }

    // set the tangent
    (*tang)(0,1) = alpha;
    (*tang)(1,0) = alpha;
    
    // set the myDOF_Groups tags indicating the attached id's of the
    // DOF_Group objects
    DOF_Group *theNodesDOFs = theNode->getDOF_GroupPtr();
    if (theNodesDOFs == 0) {
	opserr << "WARNING LagrangeSP_FE::LagrangeSP_FE()";
	opserr << " - no DOF_Group with Constrained Node\n";
	exit(-1);	
    }    

    myDOF_Groups(0) = theNodesDOFs->getTag();
    myDOF_Groups(1) = theDofGroup->getTag();
}

double
Pressure_Constraint::getPressure(int last)
{
    if (pval != 0) {
	return pval[0];
    }

    Domain* theDomain = this->getDomain();
    if(theDomain == 0) {
        opserr<<"WARNING: domain has not been set";
        opserr<<" -- Pressure_Constraint::getPressureNode\n";
        return 0;
    }
    Node* pNode = theDomain->getNode(pTag);
    if(pNode == 0) return 0.0;
    const Vector& vel = pNode->getVel();
    if(last == 1) {
        if(vel.Size()==0) return 0.0;
        return vel(0);
    }
    return 0.0;
}

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

void *
OPS_RotationShearCurve(void)
{
  if (shearCurveCount == 0) {
    opserr << "RotationShearCurve limit curve - Written by MRL UT Austin Copyright 2012 -  Use at your Own Peril \n";
    shearCurveCount++;
  }

  int argc = OPS_GetNumRemainingInputArgs();

  if (!(argc == 9 || argc == 23)) { 
    opserr << "WARNING RotationShearCurve -- insufficient arguments\n";
    opserr << "For direct input of shear curve parameters and degrading slope want:\n\n";
    opserr << "limitCurve RotationShearCurve crvTag? eleTag? \n";
    opserr << "ndI? ndJ? rotAxis? Vn? Vr? Kdeg? rotLim? \n" << endln;
    
    opserr << "OR for calibrated shear curve and degrading slope want:\n\n";
    opserr << "limitCurve RotationShearCurve crvTag? eleTag?\n";
    opserr << "ndI? ndJ? rotAxis? Vn? Vr? Kdeg? defType?\n";
    opserr << "b? d? h? L? st? As? Acc? ld? db? rhot? f'c?\n";
    opserr << "fy? fyt? delta?\n" << endln;
    
    return 0;
  }

  int iTagData[2];	  
  int iNodeData[3];	  
  double dKdegData[3];  
  double dRotLimData[1];
  int iTypeData[1];	  
  double dPropData[14]; 
  
  int numData;
  numData = 2;
  if (OPS_GetIntInput(&numData, iTagData) != 0) {
    opserr << "WARNING RotationShearCurve -- invalid crvTag? eleTag?\n" << endln;
    return 0;
  }
  int eleTag = iTagData[1];
  Domain *theDomain = 0;
  theDomain = OPS_GetDomain();
  if (theDomain == 0) {
    opserr << "WARNING RotationShearCurve -- Pointer to Domain was not returned\n" << endln;
    return 0;
  }
  Element *theElement = 0;
  theElement = theDomain->getElement(eleTag);
  if (theElement == 0) {
    opserr << "WARNING RotationShearCurve -- Element with tag " << iTagData[1] << " does not exist for shear curve tag " << iTagData[0] << endln << endln;
    return 0;
  }
  numData = 3;
  if (OPS_GetIntInput(&numData, iNodeData) != 0) {
    opserr << "WARNING RotationShearCurve -- invalid ndI? ndJ? rotAxis?\n" << endln;
    return 0;
  }
  Node *theNodeI = 0;
  theNodeI = theDomain->getNode(iNodeData[0]);
  if (theNodeI == 0) {
    opserr << "WARNING RotationShearCurve -- Node with tag " << iNodeData[0] << " does not exist for shear curve tag " << iTagData[0] << endln << endln;
    return 0;
  }	
  Node *theNodeJ = 0;
  theNodeJ = theDomain->getNode(iNodeData[1]);
  if (theNodeJ == 0) {
    opserr << "WARNING RotationShearCurve -- Node with tag " << iNodeData[1] << " does not exist for shear curve tag " << iTagData[0] << endln << endln;
    return 0;
  }	
  if (iNodeData[2] < 3 || iNodeData[2] > 6) {
    opserr << "WARNING RotationShearCurve -- rotAxis is invalid\n";
    opserr << "rotAxis = 3 -- Rotation about z-axis - 2D\n";
    opserr << "rotAxis = 4 -- Rotation about x-axis - 3D\n";
    opserr << "rotAxis = 5 -- Rotation about y-axis - 3D\n";
    opserr << "rotAxis = 6 -- Rotation about z-axis - 3D\n" << endln;
    return 0;
  }
  if (argc == 9) {
    numData = 3;
    if (OPS_GetDoubleInput(&numData, dKdegData) != 0) {
      opserr << "WARNING RotationShearCurve -- invalid Vn? Vr? Kdeg?\n" << endln;
      return 0;
    }
    if (dKdegData[0] != -1 && !(dKdegData[0] > 0)) {
      opserr << "WARNING RotationShearCurve --  Vn input is invalid\n";
      opserr << "Vn = -1 -- Shear critical limit is not used\n";
      opserr << "Vn > 0 -- Shear critical limit is the input value\n" << endln;
      return 0;
    }
    if (dKdegData[1] < -1) {
      opserr << "WARNING RotationShearCurve -- Vr input is invalid\n";
      opserr << "Vr = -1 -- Residual shear strength = 0.2*(maximum shear at failure)\n";
      opserr << "-1 < Vr < 0 -- Residual shear strength = Vr*(maximum shear at failure)\n";
      opserr << "Vr >= 0 -- Residual shear strength is the input value\n" << endln;
      return 0;
    }
    if (dKdegData[2] >= 0) {
      opserr << "WARNING RotationShearCurve -- Kdeg input is invalid\n";
      opserr << "The degrading slope must be less than zero\n" << endln;
      return 0;
    }
//.........这里部分代码省略.........

int
LagrangeConstraintHandler::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 == 0) || (theModel == 0) || (theIntegrator == 0)) {
	opserr << "WARNING LagrangeConstraintHandler::handle() - ";
	opserr << " setLinks() has not been called\n";
	return -1;
    }

    // get number ofelements and nodes in the domain 
    // and init the theFEs and theDOFs arrays

    int numConstraints = 0;
    SP_ConstraintIter &theSPss = theDomain->getDomainAndLoadPatternSPs();
    SP_Constraint *spPtr;
    while ((spPtr = theSPss()) != 0)
      numConstraints++;

    numConstraints += theDomain->getNumMPs();

    //create a DOF_Group for each Node and add it to the AnalysisModel.
    //    : must of course set the initial IDs
    NodeIter &theNod = theDomain->getNodes();
    Node *nodPtr;
    MP_Constraint *mpPtr;    
    DOF_Group *dofPtr;
    
    int numDofGrp = 0;
    int count3 = 0;
    int countDOF =0;
    while ((nodPtr = theNod()) != 0) {
	if ((dofPtr = new DOF_Group(numDofGrp++, nodPtr)) == 0) {
	    opserr << "WARNING LagrangeConstraintHandler::handle() ";
	    opserr << "- ran out of memory";
	    opserr << " creating DOF_Group " << numDofGrp++ << endln;	
	    return -4;    		
	}
	// initially set all the ID value to -2
	
	const ID &id = dofPtr->getID();
	for (int j=0; j < id.Size(); j++) {
	    dofPtr->setID(j,-2);
	    countDOF++;
	}

	nodPtr->setDOF_GroupPtr(dofPtr);
	theModel->addDOF_Group(dofPtr);
    }

    // create the FE_Elements for the Elements and add to the AnalysisModel
    ElementIter &theEle = theDomain->getElements();
    Element *elePtr;

    int numFeEle = 0;
    FE_Element *fePtr;
    while ((elePtr = theEle()) != 0) {

      // only create an FE_Element for a subdomain element if it does not
      // do independent analysis .. then subdomain part of this analysis so create
      // an FE_element & set subdomain to point to it.
      if (elePtr->isSubdomain() == true) {
	Subdomain *theSub = (Subdomain *)elePtr;
	if (theSub->doesIndependentAnalysis() == false) {
	  if ((fePtr = new FE_Element(numFeEle++, elePtr)) == 0) {
	    opserr << "WARNING PlainHandler::handle() - ran out of memory";
	    opserr << " creating FE_Element " << elePtr->getTag() << endln; 
	    return -5;
	  }		

	  theModel->addFE_Element(fePtr);
	  theSub->setFE_ElementPtr(fePtr);

	} //  if (theSub->doesIndependentAnalysis() == false) {

      } else {

	// just a regular element .. create an FE_Element for it & add to AnalysisModel
	if ((fePtr = new FE_Element(numFeEle++, elePtr)) == 0) {
	  opserr << "WARNING PlainHandler::handle() - ran out of memory";
	  opserr << " creating FE_Element " << elePtr->getTag() << endln; 
	  return -5;
	}
	
	theModel->addFE_Element(fePtr);
      }
    }

    // create the LagrangeSP_FE for the SP_Constraints and 
    // add to the AnalysisModel

    SP_ConstraintIter &theSPs = theDomain->getDomainAndLoadPatternSPs();
    while ((spPtr = theSPs()) != 0) {
	if ((dofPtr = new LagrangeDOF_Group(numDofGrp++, *spPtr)) == 0) {
	    opserr << "WARNING LagrangeConstraintHandler::handle()";
	    opserr << " - ran out of memory";
//.........这里部分代码省略.........

RigidBeam::RigidBeam(Domain &theDomain, int nR, int nC) {


    // get a pointer to the retained and constrained nodes - make sure they exist
    Node *nodeR = theDomain.getNode(nR);
    if (nodeR == 0) {
        opserr << "RigidBeam::RigidBeam - retained Node" <<  nR <<  "not in domain\n";
        return;
    }
    Node *nodeC = theDomain.getNode(nC);
    if (nodeR == 0) {
        opserr << "RigidBeam::RigidBeam - constrained Node" <<  nC <<  "not in domain\n";
        return;
    }

    // get the coordinates of the two nodes - check dimensions are the same FOR THE MOMENT
    const Vector &crdR = nodeR->getCrds();
    const Vector &crdC = nodeC->getCrds();
    int dimR = crdR.Size();
    int dimC = crdC.Size();
    if (dimR != dimC) {
        opserr << "RigidBeam::RigidBeam - mismatch in dimension "  <<
               "between constrained Node " <<  nC <<  " and Retained node" << nR << endln;
        return;
    }

    // check the number of dof at each node is the same
    int numDOF = nodeR->getNumberDOF();
    if (numDOF != nodeC->getNumberDOF()) {
        opserr << "RigidBeam::RigidBeam - mismatch in numDOF "  <<
               "between constrained Node " <<  nC <<  " and Retained node" << nR << endln;
        return;
    }

    // check the number of dof at the nodes >= dimension of problem
    if(numDOF < dimR) {
        opserr << "RigidBeam::RigidBeam - numDOF at nodes " <<
               nR << " and " <<  nC <<  "must be >= dimension of problem\n";
        return;
    }


    // create the ID to identify the constrained dof
    ID id(numDOF);

    // construct the tranformation matrix Ccr, where  Uc = Ccr Ur & set the diag, Ccr = I
    Matrix mat(numDOF,numDOF);
    mat.Zero();

    // set the values
    for (int i=0; i<numDOF; i++) {
        mat(i,i) = 1.0;
        id(i) = i;
    }

    // if there are rotational dof - we must modify Ccr DONE ASSUMING SMALL ROTATIONS
    if (dimR != numDOF) {
        if (dimR == 2 && numDOF == 3) {
            double deltaX = crdC(0) - crdR(0);
            double deltaY = crdC(1) - crdR(1);
            mat(0,2) = -deltaY;
            mat(1,2) = deltaX;
        } else if (dimR == 3 && numDOF == 6) {
            double deltaX = crdC(0) - crdR(0);
            double deltaY = crdC(1) - crdR(1);
            double deltaZ = crdC(2) - crdR(2);
            // rotation about z/3 axis
            mat(0,5) = -deltaY;
            mat(1,5) = deltaX;

            // rotation about y/2 axis
            mat(0,4) = deltaZ;
            mat(2,4) = -deltaX;

            // rotation about x/1 axis
            mat(1,3) = -deltaZ;
            mat(2,3) = deltaY;
        } else { // not valid
            opserr << "RigidBeam::RigidBeam -  for nodes " <<
                   nR << "and " << nC <<  "nodes do not have valid numDOF for their dimension\n";
            return;
        }

    }

    // create the MP_Constraint
    MP_Constraint *newC = new MP_Constraint(nR, nC, mat, id, id);

    if (newC == 0) {
        opserr << "RigidBeam::RigidBeam - for nodes " << nC << " and " << nR << ", out of memory\n";
    } else {
        // add the constraint to the domain
        if (theDomain.addMP_Constraint(newC) == false) {
            opserr << "RigidBeam::RigidBeam - for nodes " << nC << " and " << nR << ", could not add to domain\n";
            delete newC;
        }
    }
}

int
TransformationConstraintHandler::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 == 0) || (theModel == 0) || (theIntegrator == 0)) {
	opserr << "WARNING TransformationConstraintHandler::handle() - ";
	opserr << " setLinks() has not been called\n";
	return -1;
    }
    
    // get number ofelements and nodes in the domain 
    // and init the theFEs and theDOFs arrays
    int numMPConstraints = theDomain->getNumMPs();

    //    int numSPConstraints = theDomain->getNumSPs();    
    int numSPConstraints = 0;
    SP_ConstraintIter &theSP1s = theDomain->getDomainAndLoadPatternSPs();
    SP_Constraint *theSP1; 
    while ((theSP1 = theSP1s()) != 0) 
	numSPConstraints++;
    
    numDOF = 0;
    ID transformedNode(0, 64);

    int i;
    
    // create an ID of constrained node tags in MP_Constraints
    ID constrainedNodesMP(0, numMPConstraints);
    MP_Constraint **mps =0;
    if (numMPConstraints != 0) {
	mps = new MP_Constraint *[numMPConstraints];
	if (mps == 0) {
	    opserr << "WARNING TransformationConstraintHandler::handle() - ";
	    opserr << "ran out of memory for MP_Constraints"; 
	    opserr << " array of size " << numMPConstraints << endln;
	    return -3;	    
	}
	MP_ConstraintIter &theMPs = theDomain->getMPs();
	MP_Constraint *theMP; 
	int index = 0;
	while ((theMP = theMPs()) != 0) {
	  int nodeConstrained = theMP->getNodeConstrained();
	  if (transformedNode.getLocation(nodeConstrained) < 0)
	    transformedNode[numDOF++] = nodeConstrained;
	  constrainedNodesMP[index] = nodeConstrained;
	  mps[index] = theMP;
	  index++;
	}	
    }

    // create an ID of constrained node tags in SP_Constraints
    ID constrainedNodesSP(0, numSPConstraints);;
    SP_Constraint **sps =0;
    if (numSPConstraints != 0) {
	sps = new SP_Constraint *[numSPConstraints];
	if (sps == 0) {
	    opserr << "WARNING TransformationConstraintHandler::handle() - ";
	    opserr << "ran out of memory for SP_Constraints"; 
	    opserr << " array of size " << numSPConstraints << endln;
	    if (mps != 0) delete [] mps;
	    if (sps != 0) delete [] sps;
	    return -3;	    
	}
	SP_ConstraintIter &theSPs = theDomain->getDomainAndLoadPatternSPs();
	SP_Constraint *theSP; 
	int index = 0;
	while ((theSP = theSPs()) != 0) {
	  int constrainedNode = theSP->getNodeTag();
	  if (transformedNode.getLocation(constrainedNode) < 0)
	    transformedNode[numDOF++] = constrainedNode;	    
	  constrainedNodesSP[index] = constrainedNode;
	  sps[index] = theSP;
	  index++;
	}	
    }

    // create an array for the DOF_Groups and zero it
    if ((numDOF != 0) && ((theDOFs = new DOF_Group *[numDOF]) == 0)) {
	opserr << "WARNING TransformationConstraintHandler::handle() - ";
        opserr << "ran out of memory for DOF_Groups";
	opserr << " array of size " << numDOF << endln;
	return -3;    
    }    
    for (i=0; i<numDOF; i++) theDOFs[i] = 0;

    //create a DOF_Group for each Node and add it to the AnalysisModel.
    //    :must of course set the initial IDs
    NodeIter &theNod = theDomain->getNodes();
    Node *nodPtr;

    int numDofGrp = 0;
    int count3 = 0;
    int countDOF =0;
    
    numConstrainedNodes = 0;
    numDOF = 0;
//.........这里部分代码省略.........

int
ImposedMotionSP::applyConstraint(double time)
{
  // on first 
  if (theGroundMotion == 0 || theNode == 0 || theNodeResponse == 0) {
    Domain *theDomain = this->getDomain();

    theNode = theDomain->getNode(nodeTag);
    if (theNode == 0) {
      opserr << "ImposedMotionSP::applyConstraint() - node " << nodeTag << " does not exist\n";
      return -1;
    }

    int numNodeDOF = theNode->getNumberDOF();

    if (dofNumber < 0 || numNodeDOF <= dofNumber) {
      opserr << "ImposedMotionSP::applyConstraint() - dof number " << dofNumber++ << " at node " << nodeTag << " not valid\n";
      return -2;
    }

    theNodeResponse = new Vector(numNodeDOF);
    if (theNodeResponse == 0) {
      opserr << "ImposedMotionSP::applyConstraint() - out of memory\n";
      return -2;
    }
    
    LoadPattern *theLoadPattern = theDomain->getLoadPattern(patternTag);
    if (theLoadPattern == 0)
      return -3;

    theGroundMotion = theLoadPattern->getMotion(groundMotionTag);
    if (theGroundMotion == 0)
      return -4;
  }

  if (theNodeResponse == 0) 
    return -1;


  
  // now get the response from the ground motion
  theGroundMotionResponse = theGroundMotion->getDispVelAccel(time);
  
  //
  // now set the responses at the node
  //
  
  /* ***********************************************************
   * disp response the responsibility of constraint handler
   
   *theNodeResponse = theNode->getTrialDisp();
   (*theNodeResponse)(dofNumber) = theGroundMotionResponse(0);
   theNode->setTrialDisp(*theNodeResponse);
  *************************************************************/
  
  *theNodeResponse = theNode->getTrialVel();
  (*theNodeResponse)(dofNumber) = theGroundMotionResponse(1);
  theNode->setTrialVel(*theNodeResponse);    
  
  *theNodeResponse = theNode->getTrialAccel();
  (*theNodeResponse)(dofNumber) = theGroundMotionResponse(2);
  theNode->setTrialAccel(*theNodeResponse);        
  
  return 0;
}

RigidRod::RigidRod(Domain &theDomain, int nR, int nC) {

    
    // get a pointer to the retained node and constrained nodes - ensure these exist
    Node *nodeR = theDomain.getNode(nR);
    if (nodeR == 0) {
      opserr << "RigidRod::RigidRod - retained Node" <<  nR <<  "not in domain\n";
      return;
    }
    Node *nodeC = theDomain.getNode(nC);
    if (nodeR == 0) {
      opserr << "RigidRod::RigidRod - constrained Node" <<  nC <<  "not in domain\n";
      return;
    }

    // get the coordinates of the two nodes - check dimensions are the same
    const Vector &crdR = nodeR->getCrds();
    const Vector &crdC = nodeC->getCrds();
    int dimR = crdR.Size();
    int dimC = crdC.Size();
    if (dimR != dimC) {
      opserr << "RigidRod::RigidRod - mismatch in dimension " <<
	"between constrained Node " <<  nC <<  " and Retained node " << nR << endln;
      return;
    }
    
    // check the number of dof at each node is the same 
    int numDOF = nodeR->getNumberDOF();
    if (numDOF != nodeC->getNumberDOF()){ 
      opserr << "RigidRod::RigidRod - mismatch in numDOF " <<
	"between constrained Node " <<  nC <<  " and Retained node " << nR << endln;
      return;
    }

    // check the number of dof at the nodes >= dimension of problem
    if(numDOF < dimR){    
      opserr << "RigidRod::RigidRod - numDOF at nodes " << nR << " and " << nC <<
	"must be >= dimension of problem\n";

      return;
    }

    
    // create the ID to identify the constrained dof 
    ID id(dimR);

    // construct the tranformation matrix Ccr, where  Uc = Ccr Ur & set the diag
    Matrix mat(dimR,dimR);
    mat.Zero();

    // set the values
    for (int i=0; i<dimR; i++) {
      mat(i,i) = 1.0;
      id(i) = i;
    }

    // create the MP_Constraint
    MP_Constraint *newC = new MP_Constraint(nR, nC, mat, id, id);
					    
    if (newC == 0) {
      opserr << "RigidRod::RigidRod - for nodes " << nR << " and " << nC << " out of memory\n";
      exit(-1);
    } else {
      // add the constraint to the domain
      if (theDomain.addMP_Constraint(newC) == false) {
	opserr << "RigidRod::RigidRod - for nodes " << nC << " and " << nR << " could not add to domain\n",
	delete newC;
      }
    }
}

//  TransformationFE(Element *, Integrator *theIntegrator);
//	construictor that take the corresponding model element.
TransformationFE::TransformationFE(int tag, Element *ele)
:FE_Element(tag, ele), theDOFs(0), numSPs(0), theSPs(0), modID(0), 
  modTangent(0), modResidual(0), numGroups(0), numTransformedDOF(0)
{
  // set number of original dof at ele
    numOriginalDOF = ele->getNumDOF();

    // create the array of pointers to DOF_Groups
    const ID &nodes = ele->getExternalNodes();
    Domain *theDomain = ele->getDomain();
    int numNodes = nodes.Size();
    theDOFs = new DOF_Group *[numNodes];
    if (theDOFs == 0) {
	opserr << "FATAL TransformationFE::TransformationFE() - out of memory craeting ";
	opserr << "array of size : " << numNodes << " for storage of DOF_Group\n";
	exit(-1);
    }

    numGroups = numNodes;

    // now fill the array of DOF_Group pointers
    for (int i=0; i<numNodes; i++) {
	Node *theNode = theDomain->getNode(nodes(i));
	if (theNode == 0) {
	    opserr << "FATAL TransformationFE::TransformationFE() - no Node with tag: ";
	    opserr << nodes(i) << " in the domain\n";;
	    exit(-1);
	}
	DOF_Group *theDofGroup = theNode->getDOF_GroupPtr();
	if (theDofGroup == 0) {
	    opserr << "FATAL TransformationFE::TransformationFE() - no DOF_Group : ";
	    opserr << " associated with node: " << nodes(i) << " in the domain\n";;
	    exit(-1);
	}	
	theDOFs[i] = theDofGroup;
    }

    // see if theTransformation array is big enough
    // if not delete the old and create a new one
    if (numNodes > sizeTransformations) {
	if (theTransformations != 0) 
	    delete [] theTransformations;
	
	theTransformations = new Matrix *[numNodes];
	if (theTransformations == 0) {
	    opserr << "FATAL TransformationFE::TransformationFE() - out of memory ";
	    opserr << "for array of pointers for Transformation matrices of size ";
	    opserr << numNodes;
	    exit(-1);
	}		    
	sizeTransformations = numNodes;
    }	

    // if this is the first element of this type create the arrays for 
    // modified tangent and residual matrices
    if (numTransFE == 0) {

	modMatrices = new Matrix *[MAX_NUM_DOF+1];
	modVectors  = new Vector *[MAX_NUM_DOF+1];
	dataBuffer = new double[MAX_NUM_DOF*MAX_NUM_DOF];
	localKbuffer = new double[MAX_NUM_DOF*MAX_NUM_DOF];
	dofData      = new int[MAX_NUM_DOF];
	sizeBuffer = MAX_NUM_DOF*MAX_NUM_DOF;
	
	if (modMatrices == 0 || modVectors == 0 || dataBuffer == 0 ||
	    localKbuffer == 0 || dofData == 0) {
	    opserr << "TransformationFE::TransformationFE(Element *) ";
	    opserr << " ran out of memory";	    
	}
	for (int i=0; i<MAX_NUM_DOF; i++) {
	    modMatrices[i] = 0;
	    modVectors[i] = 0;
	}
    }

    // increment the number of transformations
    numTransFE++;
}

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