C++ STK_Interface类代码示例

void STK_ExodusReaderFactory::registerSidesets(STK_Interface & mesh,stk::io::MeshData & meshData) const
{
   using Teuchos::RCP;

   RCP<stk::mesh::fem::FEMMetaData> metaData = mesh.getMetaData();
   const stk::mesh::PartVector & parts = metaData->get_parts();

   stk::mesh::PartVector::const_iterator partItr;
   for(partItr=parts.begin();partItr!=parts.end();++partItr) {
      const stk::mesh::Part * part = *partItr;
      const stk::mesh::PartVector & subsets = part->subsets();
      // const CellTopologyData * ct = stk::mesh::fem::get_cell_topology(*part).getCellTopologyData();
      const CellTopologyData * ct = metaData->get_cell_topology(*part).getCellTopologyData();

      // if a side part ==> this is a sideset: now storage is recursive
      // on part contains all sub parts with consistent topology
      if(part->primary_entity_rank()==mesh.getSideRank() && ct==0 && subsets.size()>0) {
         TEUCHOS_TEST_FOR_EXCEPTION(subsets.size()!=1,std::runtime_error,
                            "STK_ExodusReaderFactory::registerSidesets error - part \"" << part->name() << 
                            "\" has more than one subset"); 

         // grab cell topology and name of subset part
         const stk::mesh::Part * ss_part = subsets[0];
         // const CellTopologyData * ss_ct = stk::mesh::fem::get_cell_topology(*ss_part).getCellTopologyData();
         const CellTopologyData * ss_ct = metaData->get_cell_topology(*ss_part).getCellTopologyData();
 
         // only add subset parts that have no topology
         if(ss_ct!=0) 
            mesh.addSideset(part->name(),ss_ct);
      }
   }
}

void buildLocalIds(const STK_Interface & mesh,
                   std::map<std::string,Teuchos::RCP<std::vector<std::size_t> > > & localIds)
{
   // defines ordering of blocks
   std::vector<std::string> blockIds;
   mesh.getElementBlockNames(blockIds);

   std::vector<std::string>::const_iterator idItr;
   for(idItr=blockIds.begin();idItr!=blockIds.end();++idItr) {
      std::string blockId = *idItr;

      localIds[blockId] = Teuchos::rcp(new std::vector<std::size_t>);
      std::vector<std::size_t> & localBlockIds = *localIds[blockId];

      // grab elements on this block
      std::vector<stk::mesh::Entity*> blockElmts;
      mesh.getMyElements(blockId,blockElmts);

      std::vector<stk::mesh::Entity*>::const_iterator itr;
      for(itr=blockElmts.begin();itr!=blockElmts.end();++itr)
         localBlockIds.push_back(mesh.elementLocalId(*itr));

      std::sort(localBlockIds.begin(),localBlockIds.end());
   }
}

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

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

    // add node and element information
    buildElements(mesh);

    // build edges and faces; fyi: addSides(mesh) builds only edges
    mesh.buildSubcells();
    mesh.buildLocalElementIDs();



    // now that edges are built, side and node sets can be added
    addSideSets(mesh);

    // set solution fields
    fillSolutionFieldData(mesh);

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

void SculptMeshFactory::buildElements(stk_classic::ParallelMachine parallelMach,STK_Interface & mesh) const
{
   struct MeshStorageStruct *mss = get_sculpt_mesh();
   int num_blocks  = mss->num_elem_blk;
  

   int *block_id = new int[num_blocks];
   //char ** element_types = new std::string[num_blocks];
   int *elements = new int[num_blocks];
   int *nodes_per_element = new int[num_blocks];
   int *element_attributes = new int[num_blocks];
   int **elmt_node_linkage = new int*[num_blocks];

   for(int b = 0; b < num_blocks; b++){
      block_id[b] = mss->block_id[b];
    //  element_types[b] = mss->element_types[b];
      elements[b] = mss->elements[b];
      nodes_per_element[b] = mss->nodes_per_element[b];
      element_attributes[b] = mss->element_attributes[b];
   }


   int elm_start = 1;
   mesh.beginModification();
   // build each block
   for(int ib=0;ib<num_blocks;ib++) {
     buildBlock(parallelMach,mesh, ib, block_id, elm_start, elements, nodes_per_element, element_attributes, elmt_node_linkage );
     elm_start += elements[ib];
   }
   mesh.endModification();
}

void LineMeshFactory::buildBlock(stk::ParallelMachine parallelMach,int xBlock,STK_Interface & mesh) const
{
   // grab this processors rank and machine size
   std::pair<int,int> sizeAndStartX = determineXElemSizeAndStart(xBlock,machSize_,machRank_);

   int myXElems_start = sizeAndStartX.first;
   int myXElems_end  = myXElems_start+sizeAndStartX.second;
   int totalXElems = nXElems_*xBlocks_;

   double deltaX = (xf_-x0_)/double(totalXElems);

   // build the nodes
   std::vector<double> coord(1,0.0);
   for(int nx=myXElems_start;nx<myXElems_end+1;++nx) {
      coord[0] = double(nx)*deltaX+x0_;
      mesh.addNode(nx+1,coord);
   }

   std::stringstream blockName;
   blockName << "eblock-" << xBlock;
   stk::mesh::Part * block = mesh.getElementBlockPart(blockName.str());

   // build the elements
   for(int nx=myXElems_start;nx<myXElems_end;++nx) {
      stk::mesh::EntityId gid = nx+1;
      std::vector<stk::mesh::EntityId> nodes(2);
      nodes[0] = nx+1;
      nodes[1] = nodes[0]+1;

      RCP<ElementDescriptor> ed = rcp(new ElementDescriptor(gid,nodes));
      mesh.addElement(ed,block);
   }
}

void LineMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
{
   mesh.beginModification();
      // build each block
      for(int xBlock=0;xBlock<xBlocks_;xBlock++) {
         buildBlock(parallelMach,xBlock,mesh);
      }
   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 CubeTetMeshFactory::buildElements(stk_classic::ParallelMachine parallelMach,STK_Interface & mesh) const
{
   mesh.beginModification();
      // build each block
      for(int xBlock=0;xBlock<xBlocks_;xBlock++) {
         for(int yBlock=0;yBlock<yBlocks_;yBlock++) {
            for(int zBlock=0;zBlock<zBlocks_;zBlock++) {
               buildBlock(parallelMach,xBlock,yBlock,zBlock,mesh);
            }
         }
      }
   mesh.endModification();
}

void MultiBlockMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
{
   mesh.beginModification();

   if(machRank_==0) {
      buildBlock(parallelMach,0,0,mesh);
   }
   else if(machRank_==1) {
      buildBlock(parallelMach,0,1,mesh);
   }
   else TEUCHOS_ASSERT(false);

   mesh.endModification();
}

  void
  CustomMeshFactory::fillSolutionFieldData(STK_Interface &mesh) const
  {
    for (int blk=0;blk<NumBlocks_;++blk) {

      std::stringstream block_id;
      block_id << "eblock-" << blk;
      
      // elements in this processor for this block
      std::vector<stk_classic::mesh::Entity*> elements;    
      mesh.getMyElements(block_id.str(), elements);

      // size of elements in the current block
      std::size_t n_elements = elements.size();
      
      // build local element index
      std::vector<std::size_t> local_ids;
      for (std::vector<stk_classic::mesh::Entity*>::const_iterator
             itr=elements.begin();itr!=elements.end();++itr) 
        local_ids.push_back(mesh.elementLocalId(*itr));

      // re-index solution fields in the same order of local_ids
      std::vector<double> charge_density_by_local_ids, electric_potential_by_local_ids;
      for (std::vector<stk_classic::mesh::Entity*>::const_iterator
             itr=elements.begin();itr!=elements.end();++itr) {
        int q = (*itr)->identifier() - OffsetToGlobalElementIDs_;
        for (int k=0;k<8;++k) {
          int loc = q*8 + k;
          charge_density_by_local_ids.push_back(ChargeDensity_[loc]);
          electric_potential_by_local_ids.push_back(ElectricPotential_[loc]);
        }
      }

      // wrap the buffer with a proper container
      FieldContainer charge_density(n_elements, 8, &charge_density_by_local_ids[0]),
        electric_potential(n_elements, 8, &electric_potential_by_local_ids[0]);

      // write out to stk mesh
      mesh.setSolutionFieldData("CHARGE_DENSITY",
                                block_id.str(),
                                local_ids,
                                charge_density, 1.0);
      
      mesh.setSolutionFieldData("ELECTRIC_POTENTIAL",
                                block_id.str(),
                                local_ids,
                                electric_potential, 1.0);
    }
  }

void SculptMeshFactory::buildNodes( stk_classic::ParallelMachine paralleMach, STK_Interface &mesh ) const
{
   struct MeshStorageStruct *mss = get_sculpt_mesh();
   int num_nodes = mss->num_nodes;


  int dimensionality = 3;

  if (num_nodes){
    int global_node_numbers;
    for(int ict = 0; ict < num_nodes; ict ++){
      global_node_numbers = mss->global_node_numbers[ict];
      std::vector<double> coord(3, 0.0);
      coord[0] = mss->coord[0*num_nodes+ict];
      coord[1] = mss->coord[1*num_nodes+ict];
      coord[2] = mss->coord[2*num_nodes+ict];
      mesh.addNode(global_node_numbers, coord );

      //std::cout<<"Node "<<global_node_numbers<<": ( "<<coord[0]<<", "<<coord[1]<<", "<<coord[2]<<" )"<<std::endl;      

    }
  }
 

} 

void MultiBlockMeshFactory::buildBlock(stk::ParallelMachine parallelMach,int xBlock,int yBlock,STK_Interface & mesh) const
{
   int myXElems_start = (machRank_==0 ? 0 : 2);
   int myXElems_end  = (machRank_==0 ? 1 : 3);
   int myYElems_start = 0;
   int myYElems_end  = 1;
   int totalXElems = nXElems_*2;
   int totalYElems = nYElems_*1;

   double x0_ = 0.0;
   double xf_ = 1.0;
   double y0_ = 0.0;
   double yf_ = 1.0;
   double deltaX = (xf_-x0_)/double(totalXElems);
   double deltaY = (yf_-y0_)/double(totalYElems);
 
   std::vector<double> coord(2,0.0);

   // build the nodes
   for(int nx=myXElems_start;nx<myXElems_end+1;++nx) {
      coord[0] = double(nx)*deltaX+x0_;
      for(int ny=myYElems_start;ny<myYElems_end+1;++ny) {
         coord[1] = double(ny)*deltaY+y0_;

         mesh.addNode(ny*(totalXElems+1)+nx+1,coord);
      }
   }

   std::stringstream blockName;
   blockName << "eblock-" << xBlock << "_" << yBlock;
   stk::mesh::Part * block = mesh.getElementBlockPart(blockName.str());

   // build the elements
   for(int nx=myXElems_start;nx<myXElems_end;++nx) {
      for(int ny=myYElems_start;ny<myYElems_end;++ny) {
         stk::mesh::EntityId gid = totalXElems*ny+nx+1;
         std::vector<stk::mesh::EntityId> nodes(4);
         nodes[0] = nx+1+ny*(totalXElems+1);
         nodes[1] = nodes[0]+1;
         nodes[2] = nodes[1]+(totalXElems+1);
         nodes[3] = nodes[2]-1;

         RCP<ElementDescriptor> ed = rcp(new ElementDescriptor(gid,nodes));
         mesh.addElement(ed,block);
      }
   }
}

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

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

   buildElements(parallelMach,mesh);

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

   addSideSets(mesh);
   addNodeSets(mesh);

   this->rebalance(mesh);
}

  void
  CustomMeshFactory::buildMetaData(STK_Interface &mesh) const
  {
    typedef shards::Hexahedron<8> HexTopo;
    const CellTopologyData *ctd = shards::getCellTopologyData<HexTopo>();
    const CellTopologyData *side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);

    for (int blk=0;blk<NumBlocks_;++blk) {
      std::stringstream block_id;
      block_id << "eblock-" << blk;

      // add element blocks
      mesh.addElementBlock(block_id.str(),ctd);

      mesh.addSolutionField("CHARGE_DENSITY",block_id.str());
      mesh.addSolutionField("ELECTRIC_POTENTIAL",block_id.str());
    }

    mesh.addSideset("left",  side_ctd);
    mesh.addSideset("right", side_ctd);
    mesh.addSideset("top",   side_ctd);
    mesh.addSideset("bottom",side_ctd);
    mesh.addSideset("front", side_ctd);
    mesh.addSideset("back",  side_ctd);

    mesh.addSideset("wall",  side_ctd);
  }

void CubeTetMeshFactory::buildMetaData(stk_classic::ParallelMachine parallelMach, STK_Interface & mesh) const
{
   typedef shards::Tetrahedron<4> TetTopo;
   const CellTopologyData * ctd = shards::getCellTopologyData<TetTopo>();
   const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);

   // build meta data
   //mesh.setDimension(2);
   for(int bx=0;bx<xBlocks_;bx++) {
      for(int by=0;by<yBlocks_;by++) {
         for(int bz=0;bz<zBlocks_;bz++) {

            std::stringstream ebPostfix;
            ebPostfix << "-" << bx << "_" << by << "_" << bz;

            // add element blocks
            mesh.addElementBlock("eblock"+ebPostfix.str(),ctd);
         }
      }
   }

   // add sidesets 
   mesh.addSideset("left",side_ctd);
   mesh.addSideset("right",side_ctd);
   mesh.addSideset("top",side_ctd);
   mesh.addSideset("bottom",side_ctd);
   mesh.addSideset("front",side_ctd);
   mesh.addSideset("back",side_ctd);
}