C++ MeshPtr::activeElements方法代码示例

 double getMinCellMeasure(){
   double minMeasure = 1e7;
   vector<ElementPtr> elems = _mesh->activeElements();
   for (vector<ElementPtr>::iterator elemIt = elems.begin();elemIt!=elems.end();elemIt++){
     minMeasure = min(minMeasure, _mesh->getCellMeasure((*elemIt)->cellID()));
   }
   return minMeasure;
 }

  vector<int> getMinCellSizeCellIDs(){
    double minMeasure = getMinCellMeasure();
    vector<int> minMeasureCellIDs;
    vector<ElementPtr> elems = _mesh->activeElements();
    for (vector<ElementPtr>::iterator elemIt = elems.begin();elemIt!=elems.end();elemIt++){
      if (minMeasure <= _mesh->getCellMeasure((*elemIt)->cellID())){
	minMeasureCellIDs.push_back((*elemIt)->cellID());
      }
    }
    return minMeasureCellIDs;
  }

// builds a [0,2]x[0,1] L shaped domain with 2 main blocks
MeshPtr MeshUtilities::buildLongRampMesh(double rampHeight, Teuchos::RCP< BilinearForm > bilinearForm, int H1Order, int pTest){

  MeshPtr mesh;
  // L-shaped domain for double ramp problem
  vector<double> A(2), B(2), C(2), D(2), E(2), F(2), G(2), H(2);
  A[0] = -1.0; A[1] = 0.0;
  B[0] = 0.0; B[1] = 0.0;
  C[0] = 1.0; C[1] = 0.0;
  D[0] = 2.0; D[1] = 0.0 + rampHeight;
  E[0] = 2.0; E[1] = 1.0;
  F[0] = 1.0; F[1] = 1.0;
  G[0] = 0.0; G[1] = 1.0;
  H[0] = -1.0; H[1] = 1.0;
  vector<vector<double> > vertices;
  vertices.push_back(A); int A_index = 0;
  vertices.push_back(B); int B_index = 1;
  vertices.push_back(C); int C_index = 2;
  vertices.push_back(D); int D_index = 3;
  vertices.push_back(E); int E_index = 4;
  vertices.push_back(F); int F_index = 5;
  vertices.push_back(G); int G_index = 6;
  vertices.push_back(H); int H_index = 7;
  vector< vector<unsigned> > elementVertices;
  vector<unsigned> el1, el2, el3;
  // left patch:
  el1.push_back(A_index); el1.push_back(B_index); el1.push_back(G_index); el1.push_back(H_index);
  // center patch:
  el2.push_back(B_index); el2.push_back(C_index); el2.push_back(F_index); el2.push_back(G_index);
  // right:
  el3.push_back(C_index); el3.push_back(D_index); el3.push_back(E_index); el3.push_back(F_index);

  elementVertices.push_back(el1);
  elementVertices.push_back(el2);
  elementVertices.push_back(el3);
  int pToAdd = pTest-H1Order;
  mesh = Teuchos::rcp( new Mesh(vertices, elementVertices, bilinearForm, H1Order, pToAdd) );  
  vector<ElementPtr> elems = mesh->activeElements();
  vector<GlobalIndexType> cellsToRefine;
  for (vector<ElementPtr>::iterator elemIt = elems.begin();elemIt!=elems.end();elemIt++){
    cellsToRefine.push_back((*elemIt)->cellID());
  }
  mesh->hRefine(cellsToRefine, RefinementPattern::regularRefinementPatternQuad()); // refine all cells (we know they're quads)
  return mesh;
}

void ErrorPercentageRefinementStrategy::refine(bool printToConsole) {
  MeshPtr mesh = this->mesh();
  
  map<GlobalIndexType, double> energyError;
  if (_rieszRep.get() != NULL) {
    _rieszRep->computeRieszRep();
    energyError = _rieszRep->getNormsSquared();
    // take square roots:
    for (map<GlobalIndexType, double>::iterator energyEntryIt = energyError.begin();
         energyEntryIt != energyError.end(); energyEntryIt++) {
      energyEntryIt->second = sqrt( energyEntryIt->second );
    }
  } else {
    energyError = _solution->globalEnergyError();
  }
  vector< Teuchos::RCP< Element > > activeElements = mesh->activeElements();
  
  double maxError = 0.0;
  double totalEnergyErrorSquared = 0.0;
  
  map<GlobalIndexType, double> cellMeasures;
  set<GlobalIndexType> cellIDs = mesh->getActiveCellIDs();
  for (set<GlobalIndexType>::iterator cellIt=cellIDs.begin(); cellIt != cellIDs.end(); cellIt++) {
    int cellID = *cellIt;
    cellMeasures[cellID] = mesh->getCellMeasure(cellID);
  }
  
  map<double, vector<GlobalIndexType> > errorReverseLookup; // an easy way to sort by error amount
  
  for (vector< Teuchos::RCP< Element > >::iterator activeElemIt = activeElements.begin();
       activeElemIt != activeElements.end(); activeElemIt++) {
    Teuchos::RCP< Element > current_element = *(activeElemIt);
    int cellID = current_element->cellID();
    double cellEnergyError = energyError.find(cellID)->second;
    
    errorReverseLookup[cellEnergyError].push_back(cellID);
    
    double h = sqrt(cellMeasures[cellID]);
    if (h > _min_h) {
      maxError = max(cellEnergyError,maxError);
    }
    totalEnergyErrorSquared += cellEnergyError * cellEnergyError;
  }
  double totalEnergyError = sqrt(totalEnergyErrorSquared);
  if ( printToConsole && _reportPerCellErrors ) {
    cout << "per-cell Energy Error Squared for cells with > 0.1% of squared energy error\n";
    for (vector< Teuchos::RCP< Element > >::iterator activeElemIt = activeElements.begin();
         activeElemIt != activeElements.end(); activeElemIt++) {
      Teuchos::RCP< Element > current_element = *(activeElemIt);
      int cellID = current_element->cellID();
      double cellEnergyError = energyError.find(cellID)->second;
      double percent = (cellEnergyError*cellEnergyError) / totalEnergyErrorSquared * 100;
      if (percent > 0.1) {
        cout << cellID << ": " << cellEnergyError*cellEnergyError << " ( " << percent << " %)\n";
      }
    }
  }

  // record results prior to refinement
  RefinementResults results;
  setResults(results, mesh->numElements(), mesh->numGlobalDofs(), totalEnergyError);
  _results.push_back(results);
  
  vector<GlobalIndexType> cellsToRefine;
  vector<GlobalIndexType> cellsToPRefine;
  
  double errorSquaredThatWeCanIgnore = (1- _percentageThreshold) * totalEnergyErrorSquared;
  
  double errorSquaredEncounteredThusFar = 0;
  for (map<double, vector<GlobalIndexType> >::iterator errorEntryIt=errorReverseLookup.begin(); errorEntryIt != errorReverseLookup.end(); errorEntryIt++) {
    double error = errorEntryIt->first;
    vector<GlobalIndexType> cells = errorEntryIt->second;
    for (vector<GlobalIndexType>::iterator cellIt = cells.begin(); cellIt != cells.end(); cellIt++) {
      errorSquaredEncounteredThusFar += error * error;
      if (errorSquaredEncounteredThusFar > errorSquaredThatWeCanIgnore) {
        GlobalIndexType cellID = *cellIt;
        double h = sqrt(cellMeasures[cellID]);
        int p = mesh->cellPolyOrder(cellID);
        
        //      cout << "refining cellID " << cellID << endl;
        if (!_preferPRefinements) {
          if (h > _min_h) {
            cellsToRefine.push_back(cellID);
          } else {
            cellsToPRefine.push_back(cellID);
          }
        } else {
          if (p < _max_p) {
            cellsToPRefine.push_back(cellID);
          } else {
            cellsToRefine.push_back(cellID);
          }
        }
      }
    }
  }
  
  if (printToConsole) {
    if (cellsToRefine.size() > 0) Camellia::print("cells for h-refinement", cellsToRefine);
    if (cellsToPRefine.size() > 0) Camellia::print("cells for p-refinement", cellsToPRefine);
//.........这里部分代码省略.........