C++ STK_Interface::getBulkData方法代码示例

void SculptMeshFactory::addNodeSets(STK_Interface & mesh) const
{
    mesh.beginModification();

    struct MeshStorageStruct *mss = get_sculpt_mesh();
    int num_node_sets  = mss->num_node_sets;
 

    if (num_node_sets) {
      int *node_set_id = new int[num_node_sets];
      int *num_nodes_in_node_set = new int[num_node_sets];
      int *num_df_in_node_set = new int[num_node_sets];
      int **node_set_nodes = new int*[num_node_sets];
      double **node_set_df = new double*[num_node_sets];

      for(int ict = 0; ict < num_node_sets; ict ++){
        node_set_id[ict] = mss->node_set_id[ict];
        num_nodes_in_node_set[ict] = mss->num_nodes_in_node_set[ict];
        num_df_in_node_set[ict] = mss->num_df_in_node_set[ict];
      }

      for(int i = 0; i < num_node_sets; i++) {
        node_set_nodes[i] = new int[num_nodes_in_node_set[i]];
        node_set_df[i] = NULL;
        if(num_nodes_in_node_set[i]) {
          for(int nct = 0; nct < num_nodes_in_node_set[i];nct ++){
            node_set_nodes[i][nct] = mss->node_set_nodes[i][nct];
          }
        }
      }
    

      for(int i = 0; i < num_node_sets; i++) {
    
        std::stringstream nodesetName;
        nodesetName << "Nodeset-" << mss->node_set_id[i];
        stk_classic::mesh::Part * nodeset = mesh.getNodeset(nodesetName.str());

        for( int j = 0; j < num_nodes_in_node_set[i]; j++ )
        {
           int node_id = node_set_nodes[i][j];
           Teuchos::RCP<stk_classic::mesh::BulkData> bulkData = mesh.getBulkData();
           if(machRank_==0)
           {  
              stk_classic::mesh::Entity * node = bulkData->get_entity(mesh.getNodeRank(),node_id);
              mesh.addEntityToNodeset(*node, nodeset);
           }
        }
      }

    }
    mesh.endModification();
}

void STK_ExodusReaderFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
{
   PANZER_FUNC_TIME_MONITOR("panzer::STK_ExodusReaderFactory::completeMeshConstruction()");

   using Teuchos::RCP;
   using Teuchos::rcp;

   if(not mesh.isInitialized())
      mesh.initialize(parallelMach);

   // grab mesh data pointer to build the bulk data
   stk::mesh::MetaData & metaData = stk::mesh::fem::FEMMetaData::get_meta_data(*mesh.getMetaData());
   stk::io::MeshData * meshData = 
         const_cast<stk::io::MeshData *>(metaData.get_attribute<stk::io::MeshData>());
         // if const_cast is wrong ... why does it feel so right?
         // I believe this is safe since we are basically hiding this object under the covers
         // until the mesh construction can be completed...below I cleanup the object myself.
   TEUCHOS_ASSERT(metaData.remove_attribute(meshData)); 
      // remove the MeshData attribute

   RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();

   // build mesh bulk data
   mesh.beginModification();
   stk::io::populate_bulk_data(*bulkData, *meshData);
   mesh.endModification();

   // put in a negative index and (like python) the restart will be from the back
   // (-1 is the last time step)
   int restartIndex = restartIndex_;
   if(restartIndex<0) {
     std::pair<int,double> lastTimeStep = meshData->m_input_region->get_max_time();
     restartIndex = 1+restartIndex+lastTimeStep.first;
   }

   // populate mesh fields with specific index
   stk::io::process_input_request(*meshData,*bulkData,restartIndex);

   mesh.buildSubcells();
   mesh.buildLocalElementIDs();

   if(restartIndex>0) // process_input_request is a no-op if restartIndex<=0 ... thus there would be no inital time
      mesh.setInitialStateTime(meshData->m_input_region->get_state_time(restartIndex));
   else
      mesh.setInitialStateTime(0.0); // no initial time to speak, might as well use 0.0

   // clean up mesh data object
   delete meshData;

   // calls Stk_MeshFactory::rebalance
   this->rebalance(mesh);
}

void STK_ExodusReaderFactory::completeMeshConstruction(STK_Interface & mesh,stk_classic::ParallelMachine parallelMach) const
{
   PANZER_FUNC_TIME_MONITOR("panzer::STK_ExodusReaderFactory::completeMeshConstruction()");

   using Teuchos::RCP;
   using Teuchos::rcp;

   if(not mesh.isInitialized())
      mesh.initialize(parallelMach);

   // grab mesh data pointer to build the bulk data
   stk_classic::mesh::MetaData & metaData = stk_classic::mesh::fem::FEMMetaData::get_meta_data(*mesh.getMetaData());
   stk_classic::io::MeshData * meshData = 
         const_cast<stk_classic::io::MeshData *>(metaData.get_attribute<stk_classic::io::MeshData>());
         // if const_cast is wrong ... why does it feel so right?
         // I believe this is safe since we are basically hiding this object under the covers
         // until the mesh construction can be completed...below I cleanup the object myself.
   TEUCHOS_ASSERT(metaData.remove_attribute(meshData)); 
      // remove the MeshData attribute

   RCP<stk_classic::mesh::BulkData> bulkData = mesh.getBulkData();

   // build mesh bulk data
   mesh.beginModification();
   stk_classic::io::populate_bulk_data(*bulkData, *meshData);

   // The following section of code is applicable if mesh scaling is
   // turned on from the input file.
   if (userMeshScaling_)
   {
     stk_classic::mesh::Field<double,stk_classic::mesh::Cartesian>* coord_field =
       metaData.get_field<stk_classic::mesh::Field<double, stk_classic::mesh::Cartesian> >("coordinates");

     std::vector<stk_classic::mesh::Bucket*> const all_node_buckets =
       bulkData->buckets(stk_classic::mesh::fem::FEMMetaData::NODE_RANK);

     stk_classic::mesh::Selector select_all_local = metaData.locally_owned_part() | metaData.globally_shared_part();
     std::vector<stk_classic::mesh::Bucket*> my_node_buckets;
     stk_classic::mesh::get_buckets(select_all_local, all_node_buckets, my_node_buckets);

     int mesh_dim = mesh.getDimension();

     // Scale the mesh
     for (size_t i=0; i < my_node_buckets.size(); ++i)
     {
       stk_classic::mesh::Bucket& b = *(my_node_buckets[i]);
       stk_classic::mesh::BucketArray<stk_classic::mesh::Field<double,stk_classic::mesh::Cartesian> > 
         coordinate_data(*coord_field, b);

       for (size_t j=0; j < b.size(); ++j) {

         int index = j;

         double inv_msf = 1.0/meshScaleFactor_;
         for (int k=0; k < mesh_dim; ++k)
           coordinate_data(k, index) = coordinate_data(k, index) * inv_msf;
       }
     }
   }

   mesh.endModification();

   // put in a negative index and (like python) the restart will be from the back
   // (-1 is the last time step)
   int restartIndex = restartIndex_;
   if(restartIndex<0) {
     std::pair<int,double> lastTimeStep = meshData->m_input_region->get_max_time();
     restartIndex = 1+restartIndex+lastTimeStep.first;
   }

   // populate mesh fields with specific index
   stk_classic::io::process_input_request(*meshData,*bulkData,restartIndex);

   mesh.buildSubcells();
   mesh.buildLocalElementIDs();

   if(restartIndex>0) // process_input_request is a no-op if restartIndex<=0 ... thus there would be no inital time
      mesh.setInitialStateTime(meshData->m_input_region->get_state_time(restartIndex));
   else
      mesh.setInitialStateTime(0.0); // no initial time to speak, might as well use 0.0

   // clean up mesh data object
   delete meshData;

   // calls Stk_MeshFactory::rebalance
   this->rebalance(mesh);
}

std::pair<Teuchos::RCP<std::vector<std::size_t> >,
          Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > >
getSideIdsAndCoords(const STK_Interface & mesh,
              const std::string & sideName, const std::string type_)
{
   Epetra_MpiComm Comm(mesh.getBulkData()->parallel());

   unsigned physicalDim = mesh.getDimension();
 
   // grab local IDs and coordinates on this side
   // and build local epetra vector
   //////////////////////////////////////////////////////////////////

   std::pair<Teuchos::RCP<std::vector<std::size_t> >,
             Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > > sidePair =
          getLocalSideIdsAndCoords(mesh,sideName,type_);

   std::vector<std::size_t> & local_side_ids = *sidePair.first;
   std::vector<Teuchos::Tuple<double,3> > & local_side_coords = *sidePair.second;
   int nodeCount = local_side_ids.size();

   // build local Epetra objects
   Epetra_Map idMap(-1,nodeCount,0,Comm);
   Teuchos::RCP<Epetra_IntVector> localIdVec = Teuchos::rcp(new Epetra_IntVector(idMap));
   Teuchos::RCP<Epetra_MultiVector> localCoordVec = Teuchos::rcp(new Epetra_MultiVector(idMap,physicalDim));

   // copy local Ids into Epetra vector
   for(std::size_t n=0;n<local_side_ids.size();n++) {
      std::size_t nodeId = local_side_ids[n];
      Teuchos::Tuple<double,3> & coords = local_side_coords[n];

      (*localIdVec)[n] = nodeId;
      for(unsigned d=0;d<physicalDim;d++)
         (*(*localCoordVec)(d))[n] = coords[d];
   }

   // fully distribute epetra vector across all processors 
   // (these are "distributed" or "dist" objects)
   //////////////////////////////////////////////////////////////

   int dist_nodeCount = idMap.NumGlobalElements();

   // build global epetra objects
   Epetra_LocalMap distMap(dist_nodeCount,0,Comm);
   Teuchos::RCP<Epetra_IntVector> distIdVec = Teuchos::rcp(new Epetra_IntVector(distMap));
   Teuchos::RCP<Epetra_MultiVector> distCoordVec = Teuchos::rcp(new Epetra_MultiVector(distMap,physicalDim));

   // export to the localVec object from the "vector" object
   Epetra_Import importer(distMap,idMap);
   TEUCHOS_ASSERT(distIdVec->Import(*localIdVec,importer,Insert)==0);
   TEUCHOS_ASSERT(distCoordVec->Import(*localCoordVec,importer,Insert)==0);

   // convert back to generic stl vector objects
   ///////////////////////////////////////////////////////////

   Teuchos::RCP<std::vector<std::size_t> > dist_side_ids
      = Teuchos::rcp(new std::vector<std::size_t>(dist_nodeCount));
   Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > dist_side_coords
      = Teuchos::rcp(new std::vector<Teuchos::Tuple<double,3> >(dist_nodeCount));

   // copy local Ids into Epetra vector
   for(std::size_t n=0;n<dist_side_ids->size();n++) {
      (*dist_side_ids)[n] = (*distIdVec)[n];

      Teuchos::Tuple<double,3> & coords = (*dist_side_coords)[n];
      for(unsigned d=0;d<physicalDim;d++)
         coords[d] = (*(*distCoordVec)(d))[n];
   }

   return std::make_pair(dist_side_ids,dist_side_coords);
}

std::pair<Teuchos::RCP<std::vector<std::size_t> >,
          Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > >
getLocalSideIdsAndCoords(const STK_Interface & mesh,
                         const std::string & sideName, const std::string type_)
{
   unsigned physicalDim = mesh.getDimension();
   
   Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
   Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();

   // grab nodes owned by requested side
   /////////////////////////////////////////////
   std::stringstream ss;
   ss << "Can't find part=\"" << sideName << "\"" << std::endl;
   stk::mesh::Part * side = metaData->get_part(sideName,ss.str().c_str());
   stk::mesh::Selector mySides = (*side) & metaData->locally_owned_part();

   stk::mesh::EntityRank rank;
   const STK_Interface::VectorFieldType * field = 0;
   unsigned int offset = 0;
   if(type_ == "coord"){
     rank = mesh.getNodeRank();
     field = & mesh.getCoordinatesField();
   } else if(type_ == "edge"){
     rank = mesh.getEdgeRank();
     field = & mesh.getEdgesField();
     offset = mesh.getMaxEntityId(mesh.getNodeRank());
   } else if(type_ == "face"){
     rank = mesh.getFaceRank();
     field = & mesh.getFacesField();
     offset = mesh.getMaxEntityId(mesh.getNodeRank())+mesh.getMaxEntityId(mesh.getEdgeRank());
   } else {
     ss << "Can't do BCs of type " << type_  << std::endl;
     TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error, ss.str())
   }

   std::vector<stk::mesh::Bucket*> const& nodeBuckets =
     bulkData->get_buckets(rank, mySides);

   // build id vector
   ////////////////////////////////////////////
   std::size_t nodeCount = 0;
   for(std::size_t b=0;b<nodeBuckets.size();b++)
      nodeCount += nodeBuckets[b]->size();

   Teuchos::RCP<std::vector<std::size_t> > sideIds
      = Teuchos::rcp(new std::vector<std::size_t>(nodeCount));
   Teuchos::RCP<std::vector<Teuchos::Tuple<double,3> > > sideCoords
      = Teuchos::rcp(new std::vector<Teuchos::Tuple<double,3> >(nodeCount));

   // loop over node buckets
   for(std::size_t b=0,index=0;b<nodeBuckets.size();b++) {
      stk::mesh::Bucket & bucket = *nodeBuckets[b];
      double const* array = stk::mesh::field_data(*field, bucket);

      for(std::size_t n=0;n<bucket.size();n++,index++) {
         (*sideIds)[index] = bulkData->identifier(bucket[n]) + offset;
         Teuchos::Tuple<double,3> & coord = (*sideCoords)[index];

         // copy coordinates into multi vector
         for(std::size_t d=0;d<physicalDim;d++)
            coord[d] = array[physicalDim*n + d];
      }
   }

   return std::make_pair(sideIds,sideCoords);
}