C++ RCP::NeighborhoodList方法代码示例

Teuchos::RCP<PeridigmNS::NeighborhoodData> PeridigmNS::Block::createNeighborhoodDataFromGlobalNeighborhoodData(Teuchos::RCP<const Epetra_BlockMap> globalOverlapScalarPointMap,
                                                                                                               Teuchos::RCP<const PeridigmNS::NeighborhoodData> globalNeighborhoodData)
{
  int numOwnedPoints = ownedScalarPointMap->NumMyElements();
  int* ownedPointGlobalIDs = ownedScalarPointMap->MyGlobalElements();

  vector<int> ownedIDs(numOwnedPoints);
  vector<int> neighborhoodList;
  vector<int> neighborhoodPtr(numOwnedPoints);

  int* const globalNeighborhoodList = globalNeighborhoodData->NeighborhoodList();
  int* const globalNeighborhoodPtr = globalNeighborhoodData->NeighborhoodPtr();

  // Create the neighborhoodList and neighborhoodPtr for this block.
  // All the IDs in the neighborhoodList and neighborhoodPtr are local IDs into 
  // the block-specific overlap map.

  for(int i=0 ; i<numOwnedPoints ; ++i){
    neighborhoodPtr[i] = (int)(neighborhoodList.size());
    int globalID = ownedPointGlobalIDs[i];
    ownedIDs[i] = overlapScalarPointMap->LID(globalID);
    int globalNeighborhoodListIndex = globalNeighborhoodPtr[globalOverlapScalarPointMap->LID(globalID)];
    int numNeighbors = globalNeighborhoodList[globalNeighborhoodListIndex++];
    neighborhoodList.push_back(numNeighbors);
    for(int j=0 ; j<numNeighbors ; ++j){
      int globalNeighborID = globalOverlapScalarPointMap->GID(globalNeighborhoodList[globalNeighborhoodListIndex++]);
      neighborhoodList.push_back( overlapScalarPointMap->LID(globalNeighborID) );
    }
  }

  // create the NeighborhoodData for this block

  Teuchos::RCP<PeridigmNS::NeighborhoodData> blockNeighborhoodData = Teuchos::rcp(new PeridigmNS::NeighborhoodData);
  blockNeighborhoodData->SetNumOwned(ownedIDs.size());
  if(ownedIDs.size() > 0){
    memcpy(blockNeighborhoodData->OwnedIDs(), 
           &ownedIDs.at(0),
           ownedIDs.size()*sizeof(int));
  }
  if(neighborhoodPtr.size() > 0){
    memcpy(blockNeighborhoodData->NeighborhoodPtr(), 
           &neighborhoodPtr.at(0),
           neighborhoodPtr.size()*sizeof(int));
  }
  blockNeighborhoodData->SetNeighborhoodListSize(neighborhoodList.size());
  if(neighborhoodList.size() > 0){
    memcpy(blockNeighborhoodData->NeighborhoodList(),
           &neighborhoodList.at(0),
           neighborhoodList.size()*sizeof(int));
  }

  return blockNeighborhoodData;
}

void 
PeridigmNS::ModelEvaluator::evalJacobian(Teuchos::RCP<Workset> workset) const
{
  const double dt = workset->timeStep;
  std::vector<PeridigmNS::Block>::iterator blockIt;
  PeridigmNS::Material::JacobianType jacobianType = *(workset->jacobianType);
  PeridigmNS::SerialMatrix& jacobian = *(workset->jacobian);

  // ---- Compute the Tangent Stiffness Matrix ----

  for(blockIt = workset->blocks->begin() ; blockIt != workset->blocks->end() ; blockIt++){

    Teuchos::RCP<PeridigmNS::NeighborhoodData> neighborhoodData = blockIt->getNeighborhoodData();
    const int numOwnedPoints = neighborhoodData->NumOwnedPoints();
    const int* ownedIDs = neighborhoodData->OwnedIDs();
    const int* neighborhoodList = neighborhoodData->NeighborhoodList();
    Teuchos::RCP<PeridigmNS::DataManager> dataManager = blockIt->getDataManager();
    Teuchos::RCP<const PeridigmNS::Material> materialModel = blockIt->getMaterialModel();

    materialModel->computeJacobian(dt, 
                                   numOwnedPoints,
                                   ownedIDs,
                                   neighborhoodList,
                                   *dataManager,
                                   jacobian,
                                   jacobianType);
  }
}

void PeridigmNS::Compute_Number_Of_Neighbors::initialize( Teuchos::RCP< std::vector<PeridigmNS::Block> > blocks ) {

  std::vector<PeridigmNS::Block>::iterator blockIt;
  for(blockIt = blocks->begin() ; blockIt != blocks->end() ; blockIt++){
    Teuchos::RCP<PeridigmNS::NeighborhoodData> neighborhoodData = blockIt->getNeighborhoodData();
    const int numOwnedPoints = neighborhoodData->NumOwnedPoints();
    const int* neighborhoodList = neighborhoodData->NeighborhoodList();
    double *numberOfNeighbors;
    blockIt->getData(m_numberOfNeighborsFieldId, PeridigmField::STEP_NONE)->ExtractView(&numberOfNeighbors);

    int neighborhoodListIndex = 0;
    for(int iID=0 ; iID<numOwnedPoints ; ++iID){
      int numNeighbors = neighborhoodList[neighborhoodListIndex++];
      numberOfNeighbors[iID] = numNeighbors;
      neighborhoodListIndex += numNeighbors;
    }
  }
}

void EvaluateForce<EvalT, Traits>::evaluateFields(typename Traits::EvalData cellData)
{
  const double dt = *cellData.timeStep;

  std::vector<PeridigmNS::Block>::iterator blockIt;
  for(blockIt = cellData.blocks->begin() ; blockIt != cellData.blocks->end() ; blockIt++){

    Teuchos::RCP<PeridigmNS::NeighborhoodData> neighborhoodData = blockIt->getNeighborhoodData();
    const int numOwnedPoints = neighborhoodData->NumOwnedPoints();
    const int* ownedIDs = neighborhoodData->OwnedIDs();
    const int* neighborhoodList = neighborhoodData->NeighborhoodList();
    Teuchos::RCP<PeridigmNS::DataManager> dataManager = blockIt->getDataManager();
    Teuchos::RCP<const PeridigmNS::Material> materialModel = blockIt->getMaterialModel();

    materialModel->computeForce(dt, 
                                numOwnedPoints,
                                ownedIDs,
                                neighborhoodList,
                                *dataManager);
  }
}

Teuchos::RCP<PeridigmNS::NeighborhoodData>
PeridigmNS::TextFileDiscretization::filterBonds(Teuchos::RCP<PeridigmNS::NeighborhoodData> unfilteredNeighborhoodData)
{
  // Set up a block bonding matrix, which defines whether or not bonds should be formed across blocks
  int numBlocks = getNumBlocks();
  std::vector< std::vector<bool> > blockBondingMatrix(numBlocks);
  for(int i=0 ; i<numBlocks ; ++i){
    blockBondingMatrix[i].resize(numBlocks, true);
  }

  if(bondFilterCommand == "None"){
    // All blocks are bonded, the blockBondingMatrix is unchanged
    return unfilteredNeighborhoodData;
  }
  else if(bondFilterCommand == "All"){
    // No blocks are bonded, the blockBondingMatrix is the identity matrix
    for(int i=0 ; i<numBlocks ; ++i){
      for(int j=0 ; j<numBlocks ; ++j){
        if(i != j)
          blockBondingMatrix[i][j] = false;
      }
    }
  }
  else{
    string msg = "**** Error, unrecognized value for \"Omit Bonds Between Blocks\":  ";
    msg += bondFilterCommand + "\n";
    msg += "**** Valid options are:  All, None\n";
    TEUCHOS_TEST_FOR_EXCEPT_MSG(true, msg);
  }

  // Create an overlap vector containing the block IDs of each cell
  Teuchos::RCP<const Epetra_BlockMap> ownedMap = getGlobalOwnedMap(1);
  Teuchos::RCP<const Epetra_BlockMap> overlapMap = getGlobalOverlapMap(1);
  Epetra_Vector blockIDs(*overlapMap);
  Epetra_Import importer(*overlapMap, *ownedMap);
  Teuchos::RCP<Epetra_Vector> ownedBlockIDs = getBlockID();
  blockIDs.Import(*ownedBlockIDs, importer, Insert);

  // Apply the block bonding matrix and create a new NeighborhoodData
  Teuchos::RCP<PeridigmNS::NeighborhoodData> neighborhoodData = Teuchos::rcp(new PeridigmNS::NeighborhoodData);
  neighborhoodData->SetNumOwned(unfilteredNeighborhoodData->NumOwnedPoints());
  memcpy(neighborhoodData->OwnedIDs(), unfilteredNeighborhoodData->OwnedIDs(), neighborhoodData->NumOwnedPoints()*sizeof(int));
  vector<int> neighborhoodListVec;
  neighborhoodListVec.reserve(unfilteredNeighborhoodData->NeighborhoodListSize());
  int* const neighborhoodPtr = neighborhoodData->NeighborhoodPtr();

  int numOwnedPoints = neighborhoodData->NumOwnedPoints();
  int* const unfilteredNeighborhoodList = unfilteredNeighborhoodData->NeighborhoodList();
  int unfilteredNeighborhoodListIndex(0);
  for(int iID=0 ; iID<numOwnedPoints ; ++iID){
    int blockID = static_cast<int>(blockIDs[iID]);
	int numUnfilteredNeighbors = unfilteredNeighborhoodList[unfilteredNeighborhoodListIndex++];
    unsigned int numNeighborsIndex = neighborhoodListVec.size();
    neighborhoodListVec.push_back(-1); // placeholder for number of neighbors
    int numNeighbors = 0;
	for(int iNID=0 ; iNID<numUnfilteredNeighbors ; ++iNID){
      int unfilteredNeighborID = unfilteredNeighborhoodList[unfilteredNeighborhoodListIndex++];
      int unfilteredNeighborBlockID = static_cast<int>(blockIDs[unfilteredNeighborID]);
      if(blockBondingMatrix[blockID-1][unfilteredNeighborBlockID-1] == true){
        neighborhoodListVec.push_back(unfilteredNeighborID);
        numNeighbors += 1;
      }
    }
    neighborhoodListVec[numNeighborsIndex] = numNeighbors;
    neighborhoodPtr[iID] = numNeighborsIndex;
  }

  neighborhoodData->SetNeighborhoodListSize(neighborhoodListVec.size());
  memcpy(neighborhoodData->NeighborhoodList(), &neighborhoodListVec[0], neighborhoodListVec.size()*sizeof(int));

  return neighborhoodData;
}

//.........这里部分代码省略.........
    TEST_FLOATING_EQUALITY((*initialX)[4],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[5],  0.25, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[6],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[7],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[8],  0.75, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[9],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[10], 0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[11], 0.75, 1.0e-16);
  }
  if(rank == 1){
    TEST_FLOATING_EQUALITY((*initialX)[0],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[1],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[2],  0.75, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[3],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[4],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[5],  0.75, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[6],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[7],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[8],  0.25, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[9],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[10], 0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[11], 0.25, 1.0e-16);
  }

  // check cell volumes
  Teuchos::RCP<Epetra_Vector> volume = discretization->getCellVolume();
  TEST_ASSERT(volume->MyLength() == 4);
  TEST_ASSERT(volume->GlobalLength() == 8);
  for(int i=0 ; i<volume->MyLength() ; ++i)
    TEST_FLOATING_EQUALITY((*volume)[i], 0.125, 1.0e-16);

  // check the neighbor lists
  Teuchos::RCP<PeridigmNS::NeighborhoodData> neighborhoodData = discretization->getNeighborhoodData();
  TEST_ASSERT(neighborhoodData->NumOwnedPoints() == 4);
  int* ownedIds = neighborhoodData->OwnedIDs();
  TEST_ASSERT(ownedIds[0] == 0);
  TEST_ASSERT(ownedIds[1] == 1);
  TEST_ASSERT(ownedIds[2] == 2);
  TEST_ASSERT(ownedIds[3] == 3);
  TEST_ASSERT(neighborhoodData->NeighborhoodListSize() == 16);
  int* neighborhood = neighborhoodData->NeighborhoodList();
  int* neighborhoodPtr = neighborhoodData->NeighborhoodPtr();
  // remember, these are local IDs on each processor, 
  // which includes both owned and ghost nodes (confusing!)
  if(rank == 0){
    TEST_ASSERT(neighborhoodPtr[0] == 0);
    TEST_ASSERT(neighborhood[0]    == 3);
    TEST_ASSERT(neighborhood[1]    == 4);
    TEST_ASSERT(neighborhood[2]    == 1);
    TEST_ASSERT(neighborhood[3]    == 2);

    TEST_ASSERT(neighborhoodPtr[1] == 4);
    TEST_ASSERT(neighborhood[4]    == 3);
    TEST_ASSERT(neighborhood[5]    == 0);
    TEST_ASSERT(neighborhood[6]    == 5);
    TEST_ASSERT(neighborhood[7]    == 3);

    TEST_ASSERT(neighborhoodPtr[2] == 8);
    TEST_ASSERT(neighborhood[8]    == 3);
    TEST_ASSERT(neighborhood[9]    == 0);
    TEST_ASSERT(neighborhood[10]   == 6);
    TEST_ASSERT(neighborhood[11]   == 3);

    TEST_ASSERT(neighborhoodPtr[3] == 12);
    TEST_ASSERT(neighborhood[12]   == 3);
    TEST_ASSERT(neighborhood[13]   == 1);
    TEST_ASSERT(neighborhood[14]   == 2);
    TEST_ASSERT(neighborhood[15]   == 7);
  }
  if(rank == 1){
    TEST_ASSERT(neighborhoodPtr[0] == 0);
    TEST_ASSERT(neighborhood[0]    == 3);
    TEST_ASSERT(neighborhood[1]    == 2);
    TEST_ASSERT(neighborhood[2]    == 6);
    TEST_ASSERT(neighborhood[3]    == 1);

    TEST_ASSERT(neighborhoodPtr[1] == 4);
    TEST_ASSERT(neighborhood[4]    == 3);
    TEST_ASSERT(neighborhood[5]    == 3);
    TEST_ASSERT(neighborhood[6]    == 0);
    TEST_ASSERT(neighborhood[7]    == 7);

    TEST_ASSERT(neighborhoodPtr[2] == 8);
    TEST_ASSERT(neighborhood[8]    == 3);
    TEST_ASSERT(neighborhood[9]    == 4);
    TEST_ASSERT(neighborhood[10]   == 3);
    TEST_ASSERT(neighborhood[11]   == 0);

    TEST_ASSERT(neighborhoodPtr[3] == 12);
    TEST_ASSERT(neighborhood[12]   == 3);
    TEST_ASSERT(neighborhood[13]   == 2);
    TEST_ASSERT(neighborhood[14]   == 5);
    TEST_ASSERT(neighborhood[15]   == 1);
  }
}

double PeridigmNS::ComputeCriticalTimeStep(const Epetra_Comm& comm, PeridigmNS::Block& block){

  Teuchos::RCP<PeridigmNS::NeighborhoodData> neighborhoodData = block.getNeighborhoodData();
  const int numOwnedPoints = neighborhoodData->NumOwnedPoints();
  const int* ownedIDs = neighborhoodData->OwnedIDs();
  const int* neighborhoodList = neighborhoodData->NeighborhoodList();
  Teuchos::RCP<const PeridigmNS::Material> materialModel = block.getMaterialModel();

  double density = materialModel()->Density();
  double bulkModulus = materialModel()->BulkModulus();

  double horizon(0.0);
  string blockName = block.getName();
  PeridigmNS::HorizonManager& horizonManager = PeridigmNS::HorizonManager::self();
  bool blockHasConstantHorizon = horizonManager.blockHasConstantHorizon(blockName);
  if(blockHasConstantHorizon)
    horizon = horizonManager.getBlockConstantHorizonValue(blockName);

  double *cellVolume, *x;
  PeridigmNS::FieldManager& fieldManager = PeridigmNS::FieldManager::self();
  block.getData(fieldManager.getFieldId("Volume"), PeridigmField::STEP_NONE)->ExtractView(&cellVolume);
  block.getData(fieldManager.getFieldId("Model_Coordinates"), PeridigmField::STEP_NONE)->ExtractView(&x);

  const double pi = boost::math::constants::pi<double>();
  double springConstant(0.0);
  if(blockHasConstantHorizon)
    springConstant = 18.0*bulkModulus/(pi*horizon*horizon*horizon*horizon);

  double minCriticalTimeStep = 1.0e50;

  int neighborhoodListIndex = 0;
  for(int iID=0 ; iID<numOwnedPoints ; ++iID){

    double timestepDenominator = 0.0;
    int nodeID = ownedIDs[iID];
    double X[3] = { x[nodeID*3], x[nodeID*3+1], x[nodeID*3+2] };
	int numNeighbors = neighborhoodList[neighborhoodListIndex++];

    if(!blockHasConstantHorizon){
      double delta = horizonManager.evaluateHorizon(blockName, X[0], X[1], X[2]);
      springConstant = 18.0*bulkModulus/(pi*delta*delta*delta*delta);
    }

    for(int iNID=0 ; iNID<numNeighbors ; ++iNID){
      int neighborID = neighborhoodList[neighborhoodListIndex++];
      double neighborVolume = cellVolume[neighborID];
      double initialDistance = sqrt( (X[0] - x[neighborID*3  ])*(X[0] - x[neighborID*3  ]) +
                                     (X[1] - x[neighborID*3+1])*(X[1] - x[neighborID*3+1]) +
                                     (X[2] - x[neighborID*3+2])*(X[2] - x[neighborID*3+2]) );

      // Issue a warning if the bond length is very very small (as in zero)
      static bool warningGiven = false;
      if(!warningGiven && initialDistance < 1.0e-50){
        cout << "\nWarning:  Possible zero length bond detected (length = " << initialDistance << ")." << endl;
        cout << "            Bonds of length zero are not valid, the input mesh may contain coincident nodes.\n" << endl;
        warningGiven = true;
      }

      timestepDenominator += neighborVolume*springConstant/initialDistance;
    }

    double criticalTimeStep = 1.0e50;
    if(numNeighbors > 0)
      criticalTimeStep = sqrt(2.0*density/timestepDenominator);
    if(criticalTimeStep < minCriticalTimeStep)
      minCriticalTimeStep = criticalTimeStep;
  }

  // Find the minimum time step for this block across all processors
  double globalMinCriticalTimeStep;
  comm.MinAll(&minCriticalTimeStep, &globalMinCriticalTimeStep, 1);

  return globalMinCriticalTimeStep;
}