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

TEUCHOS_UNIT_TEST(PdQuickGridDiscretization_MPI_np2, SimpleTensorProductMeshTest) {

  Teuchos::RCP<Epetra_Comm> comm;
  comm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));

  int numProcs = comm->NumProc();
  int rank     = comm->MyPID();

  TEST_COMPARE(numProcs, ==, 2);

  if(numProcs != 2){
     std::cerr << "Unit test runtime ERROR: utPeridigm_PdQuickGridDiscretization_MPI_np2 only makes sense on 2 processors" << std::endl;
     return;
  }

  RCP<ParameterList> discParams = rcp(new ParameterList);

  // create a 2x2x2 discretization
  // specify a spherical neighbor search with the horizon a tad longer than the mesh spacing
  discParams->set("Type", "PdQuickGrid");
  discParams->set("NeighborhoodType", "Spherical");
  ParameterList& quickGridParams = discParams->sublist("TensorProduct3DMeshGenerator");
  quickGridParams.set("Type", "PdQuickGrid");
  quickGridParams.set("X Origin", 0.0);
  quickGridParams.set("Y Origin", 0.0);
  quickGridParams.set("Z Origin", 0.0);
  quickGridParams.set("X Length", 1.0);
  quickGridParams.set("Y Length", 1.0);
  quickGridParams.set("Z Length", 1.0);
  quickGridParams.set("Number Points X", 2);
  quickGridParams.set("Number Points Y", 2);
  quickGridParams.set("Number Points Z", 2);

  // initialize the horizon manager and set the horizon to 0.501
  ParameterList blockParameterList;
  ParameterList& blockParams = blockParameterList.sublist("My Block");
  blockParams.set("Block Names", "block_1");
  blockParams.set("Horizon", 0.501);
  PeridigmNS::HorizonManager::self().loadHorizonInformationFromBlockParameters(blockParameterList);

  // create the discretization
  RCP<PdQuickGridDiscretization> discretization =
    rcp(new PdQuickGridDiscretization(comm, discParams));

  // sanity check, calling with a dimension other than 1 or 3 should throw an exception
  TEST_THROW(discretization->getGlobalOwnedMap(0), Teuchos::Exceptions::InvalidParameter);
  TEST_THROW(discretization->getGlobalOwnedMap(2), Teuchos::Exceptions::InvalidParameter);
  TEST_THROW(discretization->getGlobalOwnedMap(4), Teuchos::Exceptions::InvalidParameter);

  // basic checks on the 1d map
  Teuchos::RCP<const Epetra_BlockMap> map = discretization->getGlobalOwnedMap(1);
  TEST_ASSERT(map->NumGlobalElements() == 8);
  TEST_ASSERT(map->NumMyElements() == 4);
  TEST_ASSERT(map->ElementSize() == 1);
  TEST_ASSERT(map->IndexBase() == 0);
  TEST_ASSERT(map->UniqueGIDs() == true);
  int* myGlobalElements = map->MyGlobalElements();
  if(rank == 0){
    TEST_ASSERT(myGlobalElements[0] == 0);
    TEST_ASSERT(myGlobalElements[1] == 2);
    TEST_ASSERT(myGlobalElements[2] == 4);
    TEST_ASSERT(myGlobalElements[3] == 6);
  }
  if(rank == 1){
    TEST_ASSERT(myGlobalElements[0] == 5);
    TEST_ASSERT(myGlobalElements[1] == 7);
    TEST_ASSERT(myGlobalElements[2] == 1);
    TEST_ASSERT(myGlobalElements[3] == 3);
  }

  // check the 1d overlap map
  // for this simple discretization, everything should be ghosted on both processors
  Teuchos::RCP<const Epetra_BlockMap> overlapMap = discretization->getGlobalOverlapMap(1);
  TEST_ASSERT(overlapMap->NumGlobalElements() == 16);
  TEST_ASSERT(overlapMap->NumMyElements() == 8);
  TEST_ASSERT(overlapMap->ElementSize() == 1);
  TEST_ASSERT(overlapMap->IndexBase() == 0);
  TEST_ASSERT(overlapMap->UniqueGIDs() == false);
  myGlobalElements = overlapMap->MyGlobalElements();
  if(rank == 0){
    TEST_ASSERT(myGlobalElements[0] == 0);
    TEST_ASSERT(myGlobalElements[1] == 2);
    TEST_ASSERT(myGlobalElements[2] == 4);
    TEST_ASSERT(myGlobalElements[3] == 6);
    TEST_ASSERT(myGlobalElements[4] == 1);
    TEST_ASSERT(myGlobalElements[5] == 3);
    TEST_ASSERT(myGlobalElements[6] == 5);
    TEST_ASSERT(myGlobalElements[7] == 7);
  }
  if(rank == 1){
    TEST_ASSERT(myGlobalElements[0] == 5);
    TEST_ASSERT(myGlobalElements[1] == 7);
    TEST_ASSERT(myGlobalElements[2] == 1);
    TEST_ASSERT(myGlobalElements[3] == 3);
    TEST_ASSERT(myGlobalElements[4] == 0);
    TEST_ASSERT(myGlobalElements[5] == 2);
    TEST_ASSERT(myGlobalElements[6] == 4);
    TEST_ASSERT(myGlobalElements[7] == 6);
  }

//.........这里部分代码省略.........

void PeridigmNS::Block::createMapsFromGlobalMaps(Teuchos::RCP<const Epetra_BlockMap> globalOwnedScalarPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOverlapScalarPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOwnedVectorPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOverlapVectorPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOwnedScalarBondMap,
                                                 Teuchos::RCP<const Epetra_Vector>   globalBlockIds,
                                                 Teuchos::RCP<const PeridigmNS::NeighborhoodData> globalNeighborhoodData,
                                                 Teuchos::RCP<const PeridigmNS::NeighborhoodData> globalContactNeighborhoodData)
{
  double* globalBlockIdsPtr;
  globalBlockIds->ExtractView(&globalBlockIdsPtr);

  // Create a list of all the on-processor elements that are part of this block

  vector<int> IDs;
  IDs.reserve(globalOverlapScalarPointMap->NumMyElements()); // upper bound
  vector<int> bondIDs;
  bondIDs.reserve(globalOverlapScalarPointMap->NumMyElements());
  vector<int> bondElementSize;
  bondElementSize.reserve(globalOwnedScalarPointMap->NumMyElements());

  for(int iLID=0 ; iLID<globalOwnedScalarPointMap->NumMyElements() ; ++iLID){
    if(globalBlockIdsPtr[iLID] == blockID) {
      int globalID = globalOwnedScalarPointMap->GID(iLID);
      IDs.push_back(globalID);
    }
  }

  // Record the size of these elements in the bond map
  // Note that if an element has no bonds, it has no entry in the bondMap
  // So, the bond map and the scalar map can have a different number of entries (different local IDs)

  for(int iLID=0 ; iLID<globalOwnedScalarBondMap->NumMyElements() ; ++iLID){
    int globalID = globalOwnedScalarBondMap->GID(iLID);
    int localID = globalOwnedScalarPointMap->LID(globalID);
    if(globalBlockIdsPtr[localID] == blockID){
      bondIDs.push_back(globalID);
      bondElementSize.push_back(globalOwnedScalarBondMap->ElementSize(iLID));
    }
  }

  // Create the owned scalar point map, the owned vector point map, and the owned scalar bond map

  int numGlobalElements = -1;
  int numMyElements = IDs.size();
  int* myGlobalElements = 0;
  if(numMyElements > 0)
    myGlobalElements = &IDs.at(0);
  int elementSize = 1;
  int indexBase = 0;
  ownedScalarPointMap =
    Teuchos::rcp(new Epetra_BlockMap(numGlobalElements, numMyElements, myGlobalElements, elementSize, indexBase, globalOwnedScalarPointMap->Comm()));

  elementSize = 3;
  ownedVectorPointMap =
    Teuchos::rcp(new Epetra_BlockMap(numGlobalElements, numMyElements, myGlobalElements, elementSize, indexBase, globalOwnedScalarPointMap->Comm()));

  numMyElements = bondElementSize.size();
  myGlobalElements = 0;
  int* elementSizeList = 0;
  if(numMyElements > 0){
    myGlobalElements = &bondIDs.at(0);
    elementSizeList = &bondElementSize.at(0);
  }
  ownedScalarBondMap =
    Teuchos::rcp(new Epetra_BlockMap(numGlobalElements, numMyElements, myGlobalElements, elementSizeList, indexBase, globalOwnedScalarPointMap->Comm()));

  // Create a list of nodes that need to be ghosted (both across material boundaries and across processor boundaries)
  set<int> ghosts;

  // Check the neighborhood list for things that need to be ghosted
  int* const globalNeighborhoodList = globalNeighborhoodData->NeighborhoodList();
  int globalNeighborhoodListIndex = 0;
  for(int iLID=0 ; iLID<globalNeighborhoodData->NumOwnedPoints() ; ++iLID){
    int numNeighbors = globalNeighborhoodList[globalNeighborhoodListIndex++];
    if(globalBlockIdsPtr[iLID] == blockID) {
      for(int i=0 ; i<numNeighbors ; ++i){
        int neighborGlobalID = globalOverlapScalarPointMap->GID( globalNeighborhoodList[globalNeighborhoodListIndex + i] );
        ghosts.insert(neighborGlobalID);
      }
    }
    globalNeighborhoodListIndex += numNeighbors;
  }

  // Check the contact neighborhood list for things that need to be ghosted
  if(!globalContactNeighborhoodData.is_null()){
    int* const globalContactNeighborhoodList = globalContactNeighborhoodData->NeighborhoodList();
    int globalContactNeighborhoodListIndex = 0;
    for(int iLID=0 ; iLID<globalContactNeighborhoodData->NumOwnedPoints() ; ++iLID){
      int numNeighbors = globalContactNeighborhoodList[globalContactNeighborhoodListIndex++];
      if(globalBlockIdsPtr[iLID] == blockID) {
        for(int i=0 ; i<numNeighbors ; ++i){
          int neighborGlobalID = globalOverlapScalarPointMap->GID( globalContactNeighborhoodList[globalContactNeighborhoodListIndex + i] );
          ghosts.insert(neighborGlobalID);
        }
      }
      globalContactNeighborhoodListIndex += numNeighbors;
    }
  }

//.........这里部分代码省略.........