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

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 AlphaOSGeneralized::formElementResidual(void)
{
    // calculate Residual Force
    AnalysisModel *theModel = this->getAnalysisModel();
    LinearSOE *theSOE = this->getLinearSOE();
    
    // loop through the FE_Elements and add the residual
    FE_Element *elePtr;
    FE_EleIter &theEles = theModel->getFEs();
    while((elePtr = theEles()) != 0)  {
        if (theSOE->addB(elePtr->getResidual(this), elePtr->getID()) < 0)  {
            opserr << "WARNING AlphaOSGeneralized::formElementResidual -";
            opserr << " failed in addB for ID " << elePtr->getID();
            return -1;
        }
        if (alphaF < 1.0)  {
            if (statusFlag == CURRENT_TANGENT)  {
                if (theSOE->addB(elePtr->getK_Force(*Ut-*Upt), elePtr->getID(), alphaF-1.0) < 0)  {
                    opserr << "WARNING AlphaOSGeneralized::formElementResidual -";
                    opserr << " failed in addB for ID " << elePtr->getID();
                    return -2;
                }
            } else if (statusFlag == INITIAL_TANGENT)  {
                if (theSOE->addB(elePtr->getKi_Force(*Ut-*Upt), elePtr->getID(), alphaF-1.0) < 0)  {
                    opserr << "WARNING AlphaOSGeneralized::formElementResidual -";
                    opserr << " failed in addB for ID " << elePtr->getID();
                    return -2;
                }
            }
        }
    }

    return 0;
}

int 
PFEMIntegrator::commitSensitivity(int gradNum, int numGrads)
{

    // Loop through the FE_Elements and set unconditional sensitivities
    AnalysisModel *theAnalysisModel = this->getAnalysisModel();
    FE_Element *elePtr;
    FE_EleIter &theEles = theAnalysisModel->getFEs();    
    while((elePtr = theEles()) != 0) {
        elePtr->commitSensitivity(gradNum, numGrads);
    }

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

//.........这里部分代码省略.........
        // we use graph numberer if one was provided in constructor,
        // otherwise we number based on subdomains (all in subdomain 1 numbered first, 
        // then  those in 2 not in 1 and so on till done.
        //    GraphNumberer *theNumberer= this->getGraphNumbererPtr();

        ID theOrderedRefs(theGraph.getNumVertex());

        if(theNumberer)
          {
            // use the supplied graph numberer to number the merged graph
            theOrderedRefs= theNumberer->number(theGraph, lastDOF);     
          }
        else
          {
            // assign numbers based on the subdomains

            int loc= 0;
            for(int l=0; l<numChannels; l++)
              {
                const ID &theSubdomain= theSubdomainIDs[l];
                int numVertexSubdomain= theSubdomain.Size()/2;

                for(int i=0; i<numVertexSubdomain; i++)
                  {
                    const int vertexTagMerged= theSubdomain(i+numVertexSubdomain);
                    //  int refTag= vertexRefs[vertexTags.getLocation(vertexTagMerged)];
                    if(theOrderedRefs.getLocation(vertexTagMerged) == -1)
                      theOrderedRefs[loc++]= vertexTagMerged;
                  }
              }

            // now order those not yet ordered in p0
            for(int j=0; j<numVertexP0; j++)
              {
                int refTagP0= vertexTags[j];
                if(theOrderedRefs.getLocation(refTagP0) == -1)
                  theOrderedRefs[loc++]= refTagP0;
              }        
          }
    int count= 0;
    for(int i=0; i<theOrderedRefs.Size(); i++)
      {
        int vertexTag= theOrderedRefs(i);
        //      int vertexTag= vertexTags[vertexRefs.getLocation(tag)];
        Vertex *vertexPtr= theGraph.getVertexPtr(vertexTag);
        int numDOF= vertexPtr->getColor();
        vertexPtr->setTmp(count);
        count += numDOF;
      }

    // number own dof's
    for(int i=0; i<numVertexP0; i++  ) {
      int vertexTag= vertexTags(i);
      Vertex *vertexPtr= theGraph.getVertexPtr(vertexTag);

      int startID= vertexPtr->getTmp();
      int dofTag= vertexTag;
      DOF_Group *dofPtr;        
      dofPtr= theModel->getDOF_GroupPtr(dofTag);
      if(dofPtr == 0) {
        std::cerr << "WARNING XC::ParallelNumberer::numberDOF - ";
        std::cerr << "DOF_Group (P0) " << dofTag << "not in XC::AnalysisModel!\n";
        result= -4;
      } else {
        const ID &theDOFID= dofPtr->getID();
        int idSize= theDOFID.Size();
        for(int j=0; j<idSize; j++)
          if(theDOFID(j) == -2 || theDOFID(j) == -3) dofPtr->setID(j, startID++);
      }
    }

    // now given the ordered refs we determine the mapping for each subdomain
    // and send the id with the information back to the subdomain, which it uses to order
    // it's own graph
    for(int k=0; k<numChannels; k++)
     {
        CommParameters cp(0,*theChannels[k]);
        ID &theSubdomain= theSubdomainIDs[k];
        int numVertexSubdomain= theSubdomain.Size()/2;

        for(int i=0; i<numVertexSubdomain; i++)
          {
            int vertexTagMerged= theSubdomain[numVertexSubdomain+i];
            Vertex *vertexPtr= theGraph.getVertexPtr(vertexTagMerged);
            int startDOF= vertexPtr->getTmp();
            theSubdomain[i+numVertexSubdomain]= startDOF;
          }
        cp.sendID(theSubdomain,DistributedObj::getDbTagData(),CommMetaData(4));
        cp.receiveID(theSubdomain,DistributedObj::getDbTagData(),CommMetaData(4));
      }      
  }

  // iterate through the XC::FE_Element getting them to set their IDs
  FE_EleIter &theEle= theModel->getFEs();
  FE_Element *elePtr;
  while ((elePtr= theEle()) != 0)
    elePtr->setID();
  
  return result;
}

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