C++ FE_Element类代码示例

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

int 
ShadowSubdomain::computeNodalResponse(void)    
{

  DomainDecompositionAnalysis *theDDA = this->getDDAnalysis();
  if (theDDA != 0 && theDDA->doesIndependentAnalysis() != true) {
    FE_Element *theFePtr = this->getFE_ElementPtr();

    if (theFePtr != 0) {

      const Vector &lastChange = theFePtr->getLastResponse();
      msgData(0) =  ShadowActorSubdomain_computeNodalResponse;
      msgData(1) = lastChange.Size();
      if (numDOF != msgData(1)) {
	opserr << "ShadowSubdomain::update(void)";
	opserr << " - numDOF " << numDOF << " and size of Vector ";
	opserr << msgData(1) << "do not agree?\n";
	numDOF = msgData(1);
      }
      this->sendID(msgData);
      Vector theChange(lastChange);
      this->sendVector(theChange);
    }
  }
  
  return 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;
}

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
ConstraintHandler::doneNumberingDOF(void)
{
  // iterate through the FE_Element getting them to set their IDs
  FE_EleIter &theEle = theAnalysisModelPtr->getFEs();
  FE_Element *elePtr;
  while ((elePtr = theEle()) != 0)
    elePtr->setID();
  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
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 
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 
TransformationConstraintHandler::doneNumberingDOF(void)
{
    // iterate through the DOF_Groups telling them that their ID has now been set
    AnalysisModel *theModel1=this->getAnalysisModelPtr();
    DOF_GrpIter &theDOFS = theModel1->getDOFs();
    DOF_Group *dofPtr;
    while ((dofPtr = theDOFS()) != 0) {
       dofPtr->doneID();
    }


    // iterate through the FE_Element getting them to set their IDs
    AnalysisModel *theModel=this->getAnalysisModelPtr();
    FE_EleIter &theEle = theModel->getFEs();
    FE_Element *elePtr;
    while ((elePtr = theEle()) != 0) {
      elePtr->setID();
    }

    return 0;
}

int 
TransformationConstraintHandler::enforceSPs(void)
{
    for (int i=1; i<=numConstrainedNodes; i++) {
	// upward cast - safe as i put it in this location
	TransformationDOF_Group *theDof  =
	    (TransformationDOF_Group *)theDOFs[numDOF-i];
	theDof->enforceSPs(1);
    }
    for (int k=1; k<=numConstrainedNodes; k++) {
	// upward cast - safe as i put it in this location
	TransformationDOF_Group *theDof  =
	    (TransformationDOF_Group *)theDOFs[numDOF-k];
	theDof->enforceSPs(0);
    }

    for (int j=0; j<numFE; j++) {
      FE_Element *theEle = theFEs[j];
      theEle->updateElement();
    }

    return 0;
}

Graph &
AnalysisModel::getDOFGroupGraph(void)
{
  if (myGroupGraph == 0) {
    int numVertex = this->getNumDOF_Groups();

    if (numVertex == 0) {
	opserr << "WARNING AnalysisMode::getGroupGraph";
	opserr << "  - 0 vertices, has the Domain been populated?\n";
	exit(-1);
    }	

    //    myGroupGraph = new Graph(numVertex);
    MapOfTaggedObjects *graphStorage = new MapOfTaggedObjects();
    myGroupGraph = new Graph(*graphStorage);

    if (numVertex == 0) {
	opserr << "WARNING AnalysisMode::getGroupGraph";
	opserr << "  - out of memory\n";
	exit(-1);
    }	
	
    DOF_Group *dofPtr;

    // now create the vertices with a reference equal to the DOF_Group number.
    // and a tag which ranges from 0 through numVertex-1

    DOF_GrpIter &dofIter2 = this->getDOFs();
    int count = START_VERTEX_NUM;
    while ((dofPtr = dofIter2()) != 0) {
	int DOF_GroupTag = dofPtr->getTag();
	int DOF_GroupNodeTag = dofPtr->getNodeTag();
	int numDOF = dofPtr->getNumFreeDOF();
	Vertex *vertexPtr = new Vertex(DOF_GroupTag, DOF_GroupNodeTag, 0, numDOF);

	if (vertexPtr == 0) {
	    opserr << "WARNING DOF_GroupGraph::DOF_GroupGraph";
	    opserr << " - Not Enough Memory to create ";
	    opserr << count << "th Vertex\n";
	    return *myGroupGraph;
	}
	
	myGroupGraph->addVertex(vertexPtr);
    }

    // now add the edges, by looping over the Elements, getting their
    // IDs and adding edges between DOFs for equation numbers >= START_EQN_NUM
    
    FE_Element *elePtr;
    FE_EleIter &eleIter = this->getFEs();

    while((elePtr = eleIter()) != 0) {
	const ID &id = elePtr->getDOFtags();
	int size = id.Size();
	for (int i=0; i<size; i++) {
	    int dof1 = id(i);
	    for (int j=0; j<size; j++) 
		if (i != j) {
		    int dof2 = id(j);
		    myGroupGraph->addEdge(dof1,dof2);
		}
	}
    }
  }

  return *myGroupGraph;
}

Graph &
AnalysisModel::getDOFGraph(void)
{
  if (myDOFGraph == 0) {
    int numVertex = this->getNumDOF_Groups();

    //    myDOFGraph = new Graph(numVertex);
    MapOfTaggedObjects *graphStorage = new MapOfTaggedObjects();
    myDOFGraph = new Graph(*graphStorage);

    //
    // create a vertex for each dof
    //
    
    DOF_Group *dofPtr =0;
    DOF_GrpIter &theDOFs = this->getDOFs();
    while ((dofPtr = theDOFs()) != 0) {
      const ID &id = dofPtr->getID();
      int size = id.Size();
      for (int i=0; i<size; i++) {
	int dofTag = id(i);
	if (dofTag >= START_EQN_NUM) {
	  Vertex *vertexPtr = myDOFGraph->getVertexPtr(dofTag);
	  if (vertexPtr == 0) {
	    Vertex *vertexPtr = new Vertex(dofTag, dofTag);      
	    if (vertexPtr == 0) {
	      opserr << "WARNING AnalysisModel::getDOFGraph";
	      opserr << " - Not Enough Memory to create " << i+1 << "th Vertex\n";
	      return *myDOFGraph;
	    }
	    if (myDOFGraph->addVertex(vertexPtr, false) == false) {
	      opserr << "WARNING AnalysisModel::getDOFGraph - error adding vertex\n";
	      return *myDOFGraph;
	    }
	  }
	}
      }
    }
    
    // now add the edges, by looping over the FE_elements, getting their
    // IDs and adding edges between DOFs for equation numbers >= START_EQN_NUM
    
    FE_Element *elePtr =0;
    FE_EleIter &eleIter = this->getFEs();
    int cnt = 0;
    
    while((elePtr = eleIter()) != 0) {
      const ID &id = elePtr->getID();
      cnt++;
      int size = id.Size();
      for (int i=0; i<size; i++) {
	int eqn1 = id(i);
	
	// if eqnNum of DOF is a valid eqn number add an edge
	// to all other DOFs with valid eqn numbers.
	
	if (eqn1 >=START_EQN_NUM) {
	  for (int j=i+1; j<size; j++) {
	    int eqn2 = id(j);
	    if (eqn2 >=START_EQN_NUM)
	      myDOFGraph->addEdge(eqn1-START_EQN_NUM+START_VERTEX_NUM,
				  eqn2-START_EQN_NUM+START_VERTEX_NUM);
	  }
	}
      }
    }
  }    

  return *myDOFGraph;
}

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

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