C++ mesh::BulkData类代码示例

    void setup_elem_entities_and_connectivity_tables(const stk::mesh::BulkData& bulk, const stk::mesh::Selector& selector)
    {
        const stk::mesh::BucketVector& elementBuckets = bulk.get_buckets(stk::topology::ELEM_RANK, selector);
        unsigned numElementBuckets = elementBuckets.size();

        connBucketOffsets = DeviceViewIntType("D_connBucketOffsets", numElementBuckets);
        hostConnBucketOffsets = Kokkos::create_mirror_view(connBucketOffsets);

        elemBucketOffsets = DeviceViewIntType("D_elemBucketOffsets", numElementBuckets);
        hostElemBucketOffsets = Kokkos::create_mirror_view(elemBucketOffsets);

        elemsPerBucket = DeviceViewIntType("D_elemsPerBucket", numElementBuckets);
        hostElemsPerBucket = Kokkos::create_mirror_view(elemsPerBucket);

        nodesPerElement = DeviceViewIntType("D_nodesPerElement", numElementBuckets);
        hostNodesPerElement = Kokkos::create_mirror_view(nodesPerElement);

        unsigned numConnectivities = 0;
	unsigned numElements = 0;
	
        for (unsigned elemBucketIndex = 0; elemBucketIndex < numElementBuckets; ++elemBucketIndex)
	{
           const stk::mesh::Bucket& bucket = *elementBuckets[elemBucketIndex];

           unsigned numNodesPerElem = bucket.topology().num_nodes();
	   unsigned numElemsInBucket = bucket.size();

           hostConnBucketOffsets(elemBucketIndex) = numConnectivities;
           hostElemBucketOffsets(elemBucketIndex) = numElements;
	   hostNodesPerElement(elemBucketIndex) = numNodesPerElem;
	   hostElemsPerBucket(elemBucketIndex) = numElemsInBucket;
	   
    	   numConnectivities += numNodesPerElem*numElemsInBucket;
	   numElements += numElemsInBucket;
        }

	Kokkos::deep_copy(elemsPerBucket, hostElemsPerBucket);
	Kokkos::deep_copy(nodesPerElement, hostNodesPerElement);
	Kokkos::deep_copy(connBucketOffsets, hostConnBucketOffsets);
	Kokkos::deep_copy(elemBucketOffsets, hostElemBucketOffsets);

        elementNodeConnectivity = DeviceViewFlatConnectivityType("DElementNodeConnectivity", numConnectivities);
        hostElementNodeConnectivity =  Kokkos::create_mirror_view(elementNodeConnectivity);

        elemEntities = DeviceViewEntitiesType("DElemEntities", numElements);
        hostElemEntities = Kokkos::create_mirror_view(elemEntities);

        for (unsigned elemBucketIndex = 0; elemBucketIndex < numElementBuckets; ++elemBucketIndex)
        {
	    unsigned connBucketOffset = hostConnBucketOffsets(elemBucketIndex);
	    unsigned elemBucketOffset = hostElemBucketOffsets(elemBucketIndex);

            const stk::mesh::Bucket& bucket = *elementBuckets[elemBucketIndex];
            unsigned nodesPerElem = bucket.topology().num_nodes();

            for(unsigned elemIndex = 0; elemIndex < bucket.size(); ++elemIndex)
            {
	        unsigned connElemOffset = elemIndex*nodesPerElem + connBucketOffset;

                stk::mesh::Entity element = bucket[elemIndex];
                hostElemEntities(elemBucketOffset + elemIndex) = element;

                const stk::mesh::Entity * elemNodes = bulk.begin_nodes(element);
                for(unsigned iNode = 0; iNode < nodesPerElem; ++iNode)
                {
		    unsigned nodeOffset = connElemOffset + iNode;
                    hostElementNodeConnectivity(nodeOffset) = elemNodes[iNode];
                }
            }
        }

        Kokkos::deep_copy(elementNodeConnectivity, hostElementNodeConnectivity);
        Kokkos::deep_copy(elemEntities, hostElemEntities);
    }

stk::mesh::Selector get_owned_or_shared_selector(const stk::mesh::BulkData & bulkData)
{
    return bulkData.mesh_meta_data().locally_owned_part() | bulkData.mesh_meta_data().globally_shared_part();
}

void pack_selected_value_for_par_info(stk::CommSparse &comm, int procRank, const stk::mesh::BulkData& bulkData, stk::mesh::Entity local_element, stk::mesh::Selector sel)
{
    if(sel(bulkData.bucket(local_element)))
        comm.send_buffer(procRank).pack<int64_t>(bulkData.identifier(local_element));
}

void verify_unbuildable_element(stk::mesh::BulkData &bulk,
                                const stk::topology topo,
                                const stk::mesh::EntityIdVector & elem_node_ids,
                                const stk::mesh::EntityIdVector & side_ids,
                                const std::vector < std::vector < unsigned > > &gold_side_node_ids,
                                bool *sides_connectibility_check,
                                const stk::mesh::EntityIdVector & edge_ids,
                                const std::vector < std::vector < unsigned > > &gold_edge_node_ids,
                                bool *edges_connectibility_check)
{
  stk::mesh::EntityId element_id[1] = {1};
  stk::mesh::MetaData &meta = bulk.mesh_meta_data();
  stk::mesh::Part &elem_part = meta.declare_part_with_topology("elem_part", topo);

  meta.commit();
  bulk.modification_begin();

  stk::mesh::Entity elem = stk::mesh::declare_element(bulk, elem_part, element_id[0], elem_node_ids);

  stk::mesh::EntityVector side_nodes;
  uint num_sides = topo.num_sides();
  stk::topology::rank_t sub_topo_rank = topo.side_rank();

  for(uint i = 0; i < num_sides; ++i)
  {
    stk::topology sub_topo = topo.side_topology(i);
    side_nodes.clear();

    stk::mesh::Entity side = bulk.declare_entity(sub_topo_rank, side_ids[i], meta.get_topology_root_part(sub_topo));

    for (uint j = 0; j < sub_topo.num_nodes(); ++j)
    {
      stk::mesh::Entity side_node = bulk.get_entity(stk::topology::NODE_RANK, gold_side_node_ids[i][j]);
      side_nodes.push_back(side_node);
      bulk.declare_relation(side, side_node, j);
    }

    std::pair<stk::mesh::ConnectivityOrdinal, stk::mesh::Permutation> ordinalAndPermutation
      = stk::mesh::get_ordinal_and_permutation(bulk, elem, sub_topo_rank, side_nodes);

    if (sides_connectibility_check[i])
    {
      EXPECT_NE(ordinalAndPermutation.first, stk::mesh::ConnectivityOrdinal::INVALID_CONNECTIVITY_ORDINAL);
      EXPECT_NE(ordinalAndPermutation.second, stk::mesh::Permutation::INVALID_PERMUTATION);
    }
    else
    {
      EXPECT_EQ(ordinalAndPermutation.first, stk::mesh::ConnectivityOrdinal::INVALID_CONNECTIVITY_ORDINAL);
      EXPECT_EQ(ordinalAndPermutation.second, stk::mesh::Permutation::INVALID_PERMUTATION);
    }
  }

  if (edge_ids.empty()) {
    bulk.modification_end();
    return;
  }

  stk::mesh::EntityVector edge_nodes;
  uint num_edges = topo.num_edges();

  for(uint i = 0; i < num_edges; ++i)
  {
    edge_nodes.clear();

    stk::mesh::Entity edge = bulk.declare_entity(stk::topology::EDGE_RANK, edge_ids[i],
                                                 meta.get_topology_root_part(topo.edge_topology()));

    for (uint j = 0; j < topo.edge_topology().num_nodes(); ++j)
    {
      stk::mesh::Entity edge_node = bulk.get_entity(stk::topology::NODE_RANK, gold_edge_node_ids[i][j]);
      edge_nodes.push_back(edge_node);
      bulk.declare_relation(edge, edge_node, j);
    }

    std::pair<stk::mesh::ConnectivityOrdinal, stk::mesh::Permutation> ordinalAndPermutation
      = stk::mesh::get_ordinal_and_permutation(bulk, elem, stk::topology::EDGE_RANK, edge_nodes);

    if (edges_connectibility_check[i])
    {
      EXPECT_NE(ordinalAndPermutation.first, stk::mesh::ConnectivityOrdinal::INVALID_CONNECTIVITY_ORDINAL);
      EXPECT_NE(ordinalAndPermutation.second, stk::mesh::Permutation::INVALID_PERMUTATION);
    }
    else
    {
      EXPECT_EQ(ordinalAndPermutation.first, stk::mesh::ConnectivityOrdinal::INVALID_CONNECTIVITY_ORDINAL);
      EXPECT_EQ(ordinalAndPermutation.second, stk::mesh::Permutation::INVALID_PERMUTATION);
    }
  }

  bulk.modification_end();
}

 ElemElemGraphTester(stk::mesh::BulkData& bulkData)
   : ElemElemGraph(bulkData, bulkData.mesh_meta_data().universal_part()) {};

void Gear::mesh( stk::mesh::BulkData & M )
{
  stk::mesh::EntityRank element_rank = stk::topology::ELEMENT_RANK;
  stk::mesh::EntityRank side_rank = m_mesh_meta_data.side_rank();

  M.modification_begin();

  m_mesh = & M ;

  const unsigned p_size = M.parallel_size();
  const unsigned p_rank = M.parallel_rank();

  std::vector<size_t> counts ;
  stk::mesh::comm_mesh_counts(M, counts);

  // max_id is no longer available from comm_mesh_stats.
  // If we assume uniform numbering from 1.., then max_id
  // should be equal to counts...
  const stk::mesh::EntityId node_id_base = counts[ stk::topology::NODE_RANK ] + 1 ;
  const stk::mesh::EntityId elem_id_base = counts[ element_rank ] + 1 ;

  const unsigned long elem_id_gear_max =
    m_angle_num * ( m_rad_num - 1 ) * ( m_z_num - 1 );

  std::vector<stk::mesh::Part*> elem_parts ;
  std::vector<stk::mesh::Part*> face_parts ;
  std::vector<stk::mesh::Part*> node_parts ;

  {
    stk::mesh::Part * const p_gear = & m_gear ;
    stk::mesh::Part * const p_surf = & m_surf ;

    elem_parts.push_back( p_gear );
    face_parts.push_back( p_surf );
  }

  for ( unsigned ia = 0 ; ia < m_angle_num ; ++ia ) {
    for ( unsigned ir = 0 ; ir < m_rad_num - 1 ; ++ir ) {
      for ( unsigned iz = 0 ; iz < m_z_num - 1 ; ++iz ) {

        stk::mesh::EntityId elem_id_gear = identifier( m_z_num-1 , m_rad_num-1 , iz , ir , ia );

        if ( ( ( elem_id_gear * p_size ) / elem_id_gear_max ) == p_rank ) {

          stk::mesh::EntityId elem_id = elem_id_base + elem_id_gear ;

          // Create the node and set the model_coordinates

          const size_t ia_1 = ( ia + 1 ) % m_angle_num ;
          const size_t ir_1 = ir + 1 ;
          const size_t iz_1 = iz + 1 ;

          stk::mesh::Entity node[8] ;

          node[0] = create_node( node_parts, node_id_base, iz  , ir  , ia_1 );
          node[1] = create_node( node_parts, node_id_base, iz_1, ir  , ia_1 );
          node[2] = create_node( node_parts, node_id_base, iz_1, ir  , ia   );
          node[3] = create_node( node_parts, node_id_base, iz  , ir  , ia   );
          node[4] = create_node( node_parts, node_id_base, iz  , ir_1, ia_1 );
          node[5] = create_node( node_parts, node_id_base, iz_1, ir_1, ia_1 );
          node[6] = create_node( node_parts, node_id_base, iz_1, ir_1, ia   );
          node[7] = create_node( node_parts, node_id_base, iz  , ir_1, ia   );
#if 0 /* VERIFY_CENTROID */

          // Centroid of the element for verification

          const double TWO_PI = 2.0 * acos( -1.0 );
          const double angle = m_ang_inc * (0.5 + ia);
          const double z = m_center[2] + m_z_min + m_z_inc * (0.5 + iz);

          double c[3] = { 0 , 0 , 0 };

          for ( size_t j = 0 ; j < 8 ; ++j ) {
            double * const coord_data = field_data( m_model_coord , *node[j] );
            c[0] += coord_data[0] ;
            c[1] += coord_data[1] ;
            c[2] += coord_data[2] ;
          }
          c[0] /= 8 ; c[1] /= 8 ; c[2] /= 8 ;
          c[0] -= m_center[0] ;
          c[1] -= m_center[1] ;

          double val_a = atan2( c[1] , c[0] );
          if ( val_a < 0 ) { val_a += TWO_PI ; }
          const double err_a = angle - val_a ;
          const double err_z = z - c[2] ;

          const double eps = 100 * std::numeric_limits<double>::epsilon();

          if ( err_z < - eps || eps < err_z ||
               err_a < - eps || eps < err_a ) {
            std::string msg ;
            msg.append("problem setup element centroid error" );
            throw std::logic_error( msg );
          }
#endif

          stk::mesh::Entity elem =
            M.declare_entity( element_rank, elem_id, elem_parts );

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

    /**
     * Algorithm:
     *   1. for each element, loop over its edges, for each edge's node, loop over elements attached to node
     *   1a.   for each neighboring element, check if current edge is invalid
     */
    bool TopologyVerifier::isTopologyBad(stk::mesh::BulkData& bulk) //, stk::mesh::Part& mesh_part )
    {
      const stk::mesh::fem::FEMMetaData& meta = stk::mesh::fem::FEMMetaData::get(bulk);

      stk::mesh::Field<double, stk::mesh::Cartesian> *coord_field =
        meta.get_field<stk::mesh::Field<double, stk::mesh::Cartesian> >("coordinates");

      //mesh::Selector select_owned( meta_data.locally_owned_part() );

      const std::vector<mesh::Bucket*> & buckets = bulk.buckets(meta.element_rank() );

      for ( std::vector<mesh::Bucket *>::const_iterator ik = buckets.begin() ; ik != buckets.end() ; ++ik )
        {
          // if ( select_owned( **ik ) ) {...

          const mesh::Bucket & bucket = **ik ;

          // Number of elems in this bucket of elems and elem field data
          const unsigned number_elems = bucket.size();

          double * elem_node_data = field_data( *coord_field , bucket.begin() );
          //double * elem_centroid_data = field_data( elem_centroid_field , bucket.begin() );

          // FIXME
          if (0) { elem_node_data[0]++;}

#if 1
          const CellTopologyData * const bucket_cell_topo = stk::percept::PerceptMesh::get_cell_topology(bucket);
          int bucket_shardsId = ShardsInterfaceTable::s_singleton.lookupShardsId(bucket_cell_topo->name);
#endif

          //if (0) { std::cout << bucket_cell_topo->name; }
          if (0) { std::cout << "bucket_shardsId= " << bucket_shardsId << " name= " << bucket_cell_topo->name <<  std::endl; }

          if (0) { std::cout << "number_elems= " << number_elems << std::endl;}

          for ( unsigned i = 0 ; i < number_elems ; ++i)
            {
              mesh::Entity & elem = bucket[i] ;
              bool isDuplicateNode = isTopologyBad(elem);
              if (isDuplicateNode)
                {
                  std::cout << "duplicate node found: elem = " << elem << std::endl;
                  return true;
                }
              if (0) std::cout << "elemOfBucket= " << elem << std::endl;
              const mesh::PairIterRelation elem_nodes = elem.relations( stk::mesh::fem::FEMMetaData::NODE_RANK );

              //const CellTopologyData * const cell_topo = stk::percept::PerceptMesh::get_cell_topology(elem);
              const CellTopologyData * const cell_topo = stk::percept::PerceptMesh::get_cell_topology(elem);
              int shardsId = ShardsInterfaceTable::s_singleton.lookupShardsId(cell_topo->name);
              if (0) { std::cout << "shardsId= " << shardsId << " name= " << cell_topo->name <<  std::endl; }

              for (unsigned iedgeOrd = 0; iedgeOrd < cell_topo->edge_count; iedgeOrd++)
                {
                  //const CellTopologyData_Subcell& edge =

                  unsigned in0 = cell_topo->edge[iedgeOrd].node[0];
                  unsigned in1 = cell_topo->edge[iedgeOrd].node[1];

                  MyEdge<unsigned> potential_bad_edge(elem_nodes[in0].entity()->identifier(), elem_nodes[in1].entity()->identifier());
                  //if (potential_bad_edge.getId0() == 3 && potential_bad_edge.getId1() == 6)
                  if (0) std::cout << "potential_bad_edge: " << potential_bad_edge.getId0() << " " << potential_bad_edge.getId1() << std::endl;
                }
              for (unsigned iedgeOrd = 0; iedgeOrd < cell_topo->edge_count; iedgeOrd++)
                {
                  //const CellTopologyData_Subcell& edge =

                  unsigned in0 = cell_topo->edge[iedgeOrd].node[0];
                  unsigned in1 = cell_topo->edge[iedgeOrd].node[1];

                  MyEdge<unsigned> potential_bad_edge(elem_nodes[in0].entity()->identifier(), elem_nodes[in1].entity()->identifier());
                  //if (potential_bad_edge.getId0() == 3 && potential_bad_edge.getId1() == 6)
                  if (0) std::cout << "potential_bad_edge: " << potential_bad_edge.getId0() << " " << potential_bad_edge.getId1() << std::endl;

                  if (0 && MyEdge<unsigned>(3,6) == potential_bad_edge)
                    {
                      std::cout << "bad edge" << std::endl;
                    }

                  for (unsigned inodeOnPotBadEdge = 0; inodeOnPotBadEdge < 2; inodeOnPotBadEdge++)
                    {
                      unsigned inodeOnPotBadEdgeInElem = cell_topo->edge[iedgeOrd].node[inodeOnPotBadEdge];

                      const mesh::PairIterRelation node_elems = elem_nodes[inodeOnPotBadEdgeInElem].entity()->relations( meta.element_rank() );
                      unsigned num_elems_on_node = node_elems.size();

                      for (unsigned iele = 0; iele < num_elems_on_node; iele++)
                        {
                          mesh::Entity & elemOnNode = *node_elems[iele].entity();
                          const mesh::PairIterRelation elemOnNode_nodes = elemOnNode.relations( stk::mesh::fem::FEMMetaData::NODE_RANK );

                          const CellTopologyData * const local_cell_topo = stk::percept::PerceptMesh::get_cell_topology(elemOnNode);
                          int local_shardsId = ShardsInterfaceTable::s_singleton.lookupShardsId(local_cell_topo->name);
                          //if (1) { std::cout << "shardsId= " << shardsId << " name= " << cell_topo->name <<  std::endl; }
//.........这里部分代码省略.........

void fill_pre_req_data(
  ElemDataRequests& dataNeeded,
  const stk::mesh::BulkData& bulkData,
  stk::topology topo,
  stk::mesh::Entity elem,
  const stk::mesh::FieldBase* coordField,
  ScratchViews& prereqData)
{
  int nodesPerElem = topo.num_nodes();

  MasterElement *meSCS = dataNeeded.get_cvfem_surface_me();
  MasterElement *meSCV = dataNeeded.get_cvfem_volume_me();
  prereqData.elemNodes = bulkData.begin_nodes(elem);

  const FieldSet& neededFields = dataNeeded.get_fields();
  for(const FieldInfo& fieldInfo : neededFields) {
    stk::mesh::EntityRank fieldEntityRank = fieldInfo.field->entity_rank();
    unsigned scalarsDim1 = fieldInfo.scalarsDim1;
    bool isTensorField = fieldInfo.scalarsDim2 > 1;

    if (fieldEntityRank==stk::topology::ELEM_RANK) {
      if (isTensorField) {
        SharedMemView<double**>& shmemView = prereqData.get_scratch_view_2D(*fieldInfo.field);
        gather_elem_tensor_field(*fieldInfo.field, elem, scalarsDim1, fieldInfo.scalarsDim2, shmemView);
      }
      else {
        SharedMemView<double*>& shmemView = prereqData.get_scratch_view_1D(*fieldInfo.field);
        unsigned len = shmemView.dimension(0);
        double* fieldDataPtr = static_cast<double*>(stk::mesh::field_data(*fieldInfo.field, elem));
        for(unsigned i=0; i<len; ++i) {
          shmemView(i) = fieldDataPtr[i];
        }
      }
    }
    else if (fieldEntityRank == stk::topology::NODE_RANK) {
      if (isTensorField) {
        SharedMemView<double***>& shmemView3D = prereqData.get_scratch_view_3D(*fieldInfo.field);
        gather_elem_node_tensor_field(*fieldInfo.field, nodesPerElem, scalarsDim1, fieldInfo.scalarsDim2, bulkData.begin_nodes(elem), shmemView3D);
      }
      else {
        if (scalarsDim1 == 1) {
          SharedMemView<double*>& shmemView1D = prereqData.get_scratch_view_1D(*fieldInfo.field);
          gather_elem_node_field(*fieldInfo.field, nodesPerElem, prereqData.elemNodes, shmemView1D);
        }
        else {
          SharedMemView<double**>& shmemView2D = prereqData.get_scratch_view_2D(*fieldInfo.field);
          if (scalarsDim1 == 3) {
            gather_elem_node_field_3D(*fieldInfo.field, nodesPerElem, prereqData.elemNodes, shmemView2D);
          }
          else {
            gather_elem_node_field(*fieldInfo.field, nodesPerElem, scalarsDim1, prereqData.elemNodes, shmemView2D);
          }
        }
      }
    }
    else {
      ThrowRequireMsg(false, "Only node and element fields supported currently.");
    }
  } 
      
  SharedMemView<double**>* coordsView = nullptr;
  if (coordField != nullptr) {
    coordsView = &prereqData.get_scratch_view_2D(*coordField);
  }

  const std::set<ELEM_DATA_NEEDED>& dataEnums = dataNeeded.get_data_enums();
  double error = 0;
  for(ELEM_DATA_NEEDED data : dataEnums) {
    switch(data)
    {
      case SCS_AREAV:
         ThrowRequireMsg(meSCS != nullptr, "ERROR, meSCS needs to be non-null if SCS_AREAV is requested.");
         ThrowRequireMsg(coordsView != nullptr, "ERROR, coords null but SCS_AREAV requested.");
         meSCS->determinant(1, &((*coordsView)(0,0)), &prereqData.scs_areav(0,0), &error);
         break;
      case SCS_GRAD_OP:
         ThrowRequireMsg(meSCS != nullptr, "ERROR, meSCS needs to be non-null if SCS_GRAD_OP is requested.");
         ThrowRequireMsg(coordsView != nullptr, "ERROR, coords null but SCS_GRAD_OP requested.");
         meSCS->grad_op(1, &((*coordsView)(0,0)), &prereqData.dndx(0,0,0), &prereqData.deriv(0), &prereqData.det_j(0), &error);
         break;
      case SCS_SHIFTED_GRAD_OP:
        ThrowRequireMsg(meSCS != nullptr, "ERROR, meSCS needs to be non-null if SCS_GRAD_OP is requested.");
        ThrowRequireMsg(coordsView != nullptr, "ERROR, coords null but SCS_GRAD_OP requested.");
        meSCS->shifted_grad_op(1, &((*coordsView)(0,0)), &prereqData.dndx_shifted(0,0,0), &prereqData.deriv(0), &prereqData.det_j(0), &error);
        break;
      case SCS_GIJ:
         ThrowRequireMsg(meSCS != nullptr, "ERROR, meSCS needs to be non-null if SCS_GIJ is requested.");
         ThrowRequireMsg(coordsView != nullptr, "ERROR, coords null but SCS_GIJ requested.");
         meSCS->gij(&((*coordsView)(0,0)), &prereqData.gijUpper(0,0,0), &prereqData.gijLower(0,0,0), &prereqData.deriv(0));
         break;
      case SCV_VOLUME:
         ThrowRequireMsg(meSCV != nullptr, "ERROR, meSCV needs to be non-null if SCV_VOLUME is requested.");
         ThrowRequireMsg(coordsView != nullptr, "ERROR, coords null but SCV_VOLUME requested.");
         meSCV->determinant(1, &((*coordsView)(0,0)), &prereqData.scv_volume(0), &error);
         break;
      default: break;
    }
  }
}

inline
void setup2Block2HexMesh(stk::mesh::BulkData& bulk)
{
//
//   proc 0      proc 1
//             |
//    block_1  |  block_2
//             |
//      8----7 |    7----12
//     /    /| |   /    / |
//    5----6 3 |  6----11 10
//    | 1  |/  |  | 2  | /
//    1----2   |  2----9
//             |
//             |
//             |
//
//shared nodes 2, 3, 6, 7
//

  if (bulk.parallel_size() > 2) {
    return;
  }

  stk::mesh::MetaData& meta = bulk.mesh_meta_data();

  stk::topology hex = stk::topology::HEX_8;
  stk::mesh::Part& block_1 = meta.declare_part_with_topology("block_1", hex);
  stk::mesh::Part& block_2 = meta.declare_part_with_topology("block_2", hex);
  meta.commit();

  bulk.modification_begin();

  stk::mesh::EntityIdVector elem1_nodes {1, 2, 3, 4, 5, 6, 7, 8};
  stk::mesh::EntityIdVector elem2_nodes {2, 9, 10, 3, 6, 11, 12, 7};

  stk::mesh::EntityId elemId = 1;
  if (bulk.parallel_rank() == 0) {
    stk::mesh::declare_element(bulk, block_1, elemId, elem1_nodes);
  }
  if (bulk.parallel_rank() == 1 || bulk.parallel_size() == 1) {
    elemId = 2;
    stk::mesh::declare_element(bulk, block_2, elemId, elem2_nodes);
  }
  if(bulk.parallel_rank() == 0 && bulk.parallel_size() == 2)
  {
    bulk.add_node_sharing(bulk.get_entity(stk::topology::NODE_RANK , 2), 1);
    bulk.add_node_sharing(bulk.get_entity(stk::topology::NODE_RANK , 3), 1);
    bulk.add_node_sharing(bulk.get_entity(stk::topology::NODE_RANK , 6), 1);
    bulk.add_node_sharing(bulk.get_entity(stk::topology::NODE_RANK , 7), 1);
  }
  if(bulk.parallel_rank() == 1 && bulk.parallel_size() == 2)
  {
    bulk.add_node_sharing(bulk.get_entity(stk::topology::NODE_RANK , 2), 0);
    bulk.add_node_sharing(bulk.get_entity(stk::topology::NODE_RANK , 3), 0);
    bulk.add_node_sharing(bulk.get_entity(stk::topology::NODE_RANK , 6), 0);
    bulk.add_node_sharing(bulk.get_entity(stk::topology::NODE_RANK , 7), 0);
  }

  bulk.modification_end();
}

void assemble_elem_matrices_and_vectors(stk::mesh::BulkData& mesh,
                                        stk::mesh::FieldBase& field,
                                        stk::linsys::DofMapper& dof_mapper,
                                        fei::Matrix& matrix,
                                        fei::Vector& rhs)
{
  stk::mesh::fem::FEMMetaData &fem = stk::mesh::fem::FEMMetaData::get(mesh);
  const stk::mesh::EntityRank element_rank = fem.element_rank();

  const std::vector<stk::mesh::Bucket*>& mesh_buckets = mesh.buckets(element_rank);

  std::vector<stk::mesh::Bucket*> part_buckets;
  stk::mesh::Selector select_owned(stk::mesh::MetaData::get(mesh).locally_owned_part());
  stk::mesh::get_buckets(select_owned, mesh_buckets, part_buckets);

  int field_id = dof_mapper.get_field_id(field);

  stk::mesh::Entity& first_entity = *(part_buckets[0]->begin());
  stk::mesh::PairIterRelation rel = first_entity.relations(stk::mesh::fem::FEMMetaData::NODE_RANK);
  int num_nodes_per_elem = rel.second - rel.first;

  fei::SharedPtr<fei::MatrixGraph> matgraph = matrix.getMatrixGraph();
  int pattern_id = matgraph->definePattern(num_nodes_per_elem, stk::mesh::fem::FEMMetaData::NODE_RANK, field_id);

  std::vector<int> node_ids(num_nodes_per_elem);

  const int field_size = dof_mapper.get_fei_VectorSpace()->getFieldSize(field_id);
  const int matsize = num_nodes_per_elem*field_size*num_nodes_per_elem*field_size;
  const int vecsize = num_nodes_per_elem*field_size;

  std::vector<double> elem_matrix_1d(matsize, 0);
  std::vector<double*> elem_matrix_2d(vecsize);

  std::vector<double> elem_vector(vecsize, 0);

  for(size_t i=0; i<elem_matrix_2d.size(); ++i) {
    elem_matrix_2d[i] = &elem_matrix_1d[i*vecsize];
  }

  //fill our dummy elem-matrix:
  //This dummy matrix will be the same for every element. A real application
  //would form a different elem-matrix for each element.
  for(size_t i=0; i<elem_matrix_2d.size(); ++i) {
    double* row = elem_matrix_2d[i];
    if (i>=1) row[i-1] = -1;
    row[i] = 2;
    if (i<elem_matrix_2d.size()-1) row[i+1] = -1;

    elem_vector[i] = 1;
  }

  std::vector<int> eqn_indices(vecsize);

  for(size_t i=0; i<part_buckets.size(); ++i) {
    stk::mesh::Bucket::iterator
      b_iter = part_buckets[i]->begin(),
             b_end  = part_buckets[i]->end();
    for(; b_iter != b_end; ++b_iter) {
      stk::mesh::Entity& elem = *b_iter;
      rel = elem.relations(stk::mesh::fem::FEMMetaData::NODE_RANK);
      for(int j=0; rel.first != rel.second; ++rel.first, ++j) {
        node_ids[j] = rel.first->entity()->identifier();
      }

      matgraph->getPatternIndices(pattern_id, &node_ids[0], eqn_indices);

      matrix.sumIn(vecsize, &eqn_indices[0], vecsize, &eqn_indices[0],
                    &elem_matrix_2d[0]);
      rhs.sumIn(vecsize, &eqn_indices[0], &elem_vector[0]);
    }
  }
}

void use_case_7_generate_mesh(
  const std::string& mesh_options ,
  stk::mesh::BulkData & mesh ,
  const VectorFieldType & node_coord ,
  stk::mesh::Part & hex_block ,
  stk::mesh::Part & quad_shell_block )
{
  mesh.modification_begin();

  const unsigned parallel_size = mesh.parallel_size();
  const unsigned parallel_rank = mesh.parallel_rank();

  double t = 0 ;
  size_t num_hex = 0 ;
  size_t num_shell = 0 ;
  size_t num_nodes = 0 ;
  size_t num_block = 0 ;
  int error_flag = 0 ;

  try {

    Iogn::GeneratedMesh gmesh( mesh_options, parallel_size, parallel_rank );

    num_nodes = gmesh.node_count_proc();
    num_block = gmesh.block_count();

    t = stk::wall_time();

    std::vector<int> node_map( num_nodes , 0 );

    gmesh.node_map( node_map );

    {

      for ( size_t i = 1 ; i <= num_block ; ++i ) {
        const size_t                        num_elem = gmesh.element_count_proc(i);
        const std::pair<std::string,int> top_info = gmesh.topology_type(i);

	std::vector<int> elem_map( num_elem , 0 );
        std::vector<int> elem_conn( num_elem * top_info.second );

	gmesh.element_map( i, elem_map );
        gmesh.connectivity( i , elem_conn );

        if ( top_info.second == 8 ) {

          for ( size_t j = 0 ; j < num_elem ; ++j ) {

            const int * const local_node_id = & elem_conn[ j * 8 ] ;

            const stk::mesh::EntityId node_id[8] = {
              local_node_id[0] ,
              local_node_id[1] ,
              local_node_id[2] ,
              local_node_id[3] ,
              local_node_id[4] ,
              local_node_id[5] ,
              local_node_id[6] ,
              local_node_id[7]
            };

            const stk::mesh::EntityId elem_id = elem_map[ j ];

            stk::mesh::fem::declare_element( mesh , hex_block , elem_id , node_id );

            ++num_hex ;
          }
        }
        else if ( top_info.second == 4 ) {

          for ( size_t j = 0 ; j < num_elem ; ++j ) {

            const int * const local_node_id = & elem_conn[ j * 4 ] ;

            const stk::mesh::EntityId node_id[4] = {
              local_node_id[0] ,
              local_node_id[1] ,
              local_node_id[2] ,
              local_node_id[3]
            };

            const stk::mesh::EntityId elem_id = elem_map[ j ];

            stk::mesh::fem::declare_element( mesh , quad_shell_block , elem_id , node_id );

            ++num_shell ;
          }
        }
      }
    }

    std::vector<double> node_coordinates( 3 * node_map.size() );

    gmesh.coordinates( node_coordinates );

    if ( 3 * node_map.size() != node_coordinates.size() ) {
      std::ostringstream msg ;
      msg << "  P" << mesh.parallel_rank()
          << ": ERROR, node_map.size() = "
          << node_map.size()
//.........这里部分代码省略.........

void make_small_hybrid_mesh(stk::mesh::MetaData &meta, stk::mesh::BulkData &mesh,
                            bool user_attempt_no_induce = false, bool user_parts_force_no_induce = true)
{
    stk::ParallelMachine pm = MPI_COMM_WORLD;
    int p_size = stk::parallel_machine_size(pm);

    if(p_size > 2)
    {
        return;
    }

    const unsigned p_rank = mesh.parallel_rank();

    stk::mesh::Part * hexPart = &meta.get_topology_root_part(stk::topology::HEX_8);
    stk::mesh::Part * pyrPart = &meta.get_topology_root_part(stk::topology::PYRAMID_5);
    stk::mesh::Part * tetPart = &meta.get_topology_root_part(stk::topology::TET_4);

    if (user_attempt_no_induce)
    {
        hexPart = &meta.declare_part_with_topology("my_hex_part",stk::topology::HEX_8, user_parts_force_no_induce);
        pyrPart = &meta.declare_part_with_topology("my_pyr_part",stk::topology::PYRAMID_5, user_parts_force_no_induce);
        tetPart = &meta.declare_part_with_topology("my_tet_part",stk::topology::TET_4, user_parts_force_no_induce);

        EXPECT_EQ(user_parts_force_no_induce, hexPart->force_no_induce());
        EXPECT_EQ(user_parts_force_no_induce, pyrPart->force_no_induce());
        EXPECT_EQ(user_parts_force_no_induce, tetPart->force_no_induce());
    }

    meta.commit();

    const size_t numHex = 1;
    stk::mesh::EntityIdVector hexNodeIDs[] {
        { 1, 2, 3, 4, 5, 6, 7, 8 }
    };
    stk::mesh::EntityId hexElemIDs[] = { 1 };

    const size_t numPyr = 1;
    stk::mesh::EntityIdVector pyrNodeIDs[] {
        { 5, 6, 7, 8, 9 }
    };
    stk::mesh::EntityId pyrElemIDs[] = { 2 };

    const size_t numTet = 4;
    stk::mesh::EntityIdVector tetNodeIDs[] {
        { 7, 8, 9, 12 },
        { 6, 9, 10, 7 },
        { 7, 9, 10, 12 },
        { 7, 12, 10, 11 }
    };
    stk::mesh::EntityId tetElemIDs[] = { 3, 4, 5, 6 };

    // list of triplets: (owner-proc, shared-nodeID, sharing-proc)
    std::vector< std::vector<unsigned> > shared_nodeIDs_and_procs
    {
        { 0, 5, 1 },  // proc 0
        { 0, 6, 1 },
        { 0, 7, 1 },
        { 0, 8, 1 },
        { 1, 5, 0 },  // proc 1
        { 1, 6, 0 },
        { 1, 7, 0 },
        { 1, 8, 0 }
    };

    mesh.modification_begin();

    if (0 == p_rank) {
        for (size_t i = 0; i < numHex; ++i) {
          stk::mesh::declare_element(mesh, *hexPart, hexElemIDs[i], hexNodeIDs[i]);
        }
    }
    if ( (1 == p_rank) || (1 == p_size) )  { // setup the pyramids/tets for either np 2 or serial
        for (size_t i = 0; i < numPyr; ++i) {
          stk::mesh::declare_element(mesh, *pyrPart, pyrElemIDs[i], pyrNodeIDs[i]);
        }
        for (size_t i = 0; i < numTet; ++i) {
          stk::mesh::declare_element(mesh, *tetPart, tetElemIDs[i], tetNodeIDs[i]);
        }
    }

    if (p_size > 1)
    {
        for (size_t nodeIdx = 0, end = shared_nodeIDs_and_procs.size(); nodeIdx < end; ++nodeIdx) {
            if (p_rank == shared_nodeIDs_and_procs[nodeIdx][0]) {
                stk::mesh::EntityId nodeID = shared_nodeIDs_and_procs[nodeIdx][1];
                int sharingProc = shared_nodeIDs_and_procs[nodeIdx][2];
                stk::mesh::Entity node = mesh.get_entity(stk::topology::NODE_RANK, nodeID);
                mesh.add_node_sharing(node, sharingProc);
            }
        }
    }

    mesh.modification_end();
}

Intrepid::FieldContainer<double> STKMeshHelpers::extractEntityNodeCoordinates( 
    const Teuchos::Array<stk::mesh::Entity>& stk_entities, 
    const stk::mesh::BulkData& bulk_data,
    const int space_dim )
{
    // Cast the field.
    const stk::mesh::FieldBase* coord_field_base= 
	bulk_data.mesh_meta_data().coordinate_field();
    const stk::mesh::Field<double,FieldType>* coord_field =
	dynamic_cast<const stk::mesh::Field<double,FieldType>* >(
	    coord_field_base);

    // Allocate the coordinate array.
    int num_cells = stk_entities.size();
    int num_nodes = 0;
    stk::mesh::EntityRank stk_rank = stk::topology::INVALID_RANK;
    if ( num_cells > 0 )
    {
	stk_rank = bulk_data.entity_rank(stk_entities[0]);
	if ( stk::topology::NODE_RANK == stk_rank )
	{
	    num_nodes = 1;
	}
	else
	{
	    const stk::mesh::Entity* begin = 
		bulk_data.begin_nodes( stk_entities[0] );
	    const stk::mesh::Entity* end = 
		bulk_data.end_nodes( stk_entities[0] );
	    num_nodes = std::distance( begin, end );
	}
    }
    Intrepid::FieldContainer<double> coords( num_cells, num_nodes, space_dim );

    // Extract the coordinates.
    double* node_coords = 0;
    for ( int c = 0; c < num_cells; ++c )
    {
	if ( stk::topology::NODE_RANK == stk_rank )
	{
	    node_coords = stk::mesh::field_data( *coord_field, stk_entities[c] );
	    for ( int d = 0; d < space_dim; ++d )
	    {
		coords(c,0,d) = node_coords[d];
	    }
	}
	else
	{
	    const stk::mesh::Entity* begin = bulk_data.begin_nodes( stk_entities[c] );
	    DTK_REMEMBER(
		const stk::mesh::Entity* end = bulk_data.end_nodes( stk_entities[c] ) 
		);
	    DTK_CHECK( std::distance(begin,end) == num_nodes );
	    for ( int n = 0; n < num_nodes; ++n )
	    {
		node_coords = stk::mesh::field_data( *coord_field, begin[n] );
		for ( int d = 0; d < space_dim; ++d )
		{
		    coords(c,n,d) = node_coords[d];
		}
	    }
	}
    }

    return coords;
}

PromotedElementIO::PromotedElementIO(
  const ElementDescription& elem,
  const stk::mesh::MetaData& metaData,
  stk::mesh::BulkData& bulkData,
  const stk::mesh::PartVector& baseParts,
  const std::string& fileName,
  const VectorFieldType& coordField
) : elem_(elem),
    metaData_(metaData),
    bulkData_(bulkData),
    fileName_(fileName),
    coordinates_(coordField),
    nDim_(metaData.spatial_dimension())
{
  Ioss::Init::Initializer init_db;

  Ioss::PropertyManager properties;

  Ioss::Property intSizeAPI("INTEGER_SIZE_API", 8);
  properties.add(intSizeAPI);

  Ioss::Property intSizeDB("INTEGER_SIZE_DB", 8);
  properties.add(intSizeDB);

  databaseIO = Ioss::IOFactory::create(
        "exodus",
        fileName_,
        Ioss::WRITE_RESULTS,
        bulkData_.parallel(),
        properties
  );
  ThrowRequire(databaseIO != nullptr && databaseIO->ok(true));

  output_ = make_unique<Ioss::Region>(databaseIO, "HighOrderOutput"); //sink for databaseIO
  ThrowRequire(output_ != nullptr);

  const stk::mesh::BucketVector& elem_buckets = bulkData_.get_buckets(
    stk::topology::ELEM_RANK, stk::mesh::selectUnion(baseParts));

  size_t numSubElems = num_sub_elements(nDim_, elem_buckets, elem_.polyOrder);
  std::vector<stk::mesh::EntityId> subElemIds;

  bulkData.generate_new_ids(stk::topology::ELEM_RANK,  numSubElems,  subElemIds);
  ThrowRequire(subElemIds.size() == numSubElems);

  superElemParts_ = super_elem_part_vector(baseParts);
  ThrowRequireMsg(part_vector_is_valid_and_nonempty(superElemParts_),
    "Not all element parts have a super-element mirror");

  output_->begin_mode(Ioss::STATE_DEFINE_MODEL);
  write_node_block_definitions(superElemParts_);
  write_elem_block_definitions(superElemParts_);
  write_sideset_definitions(baseParts);
  output_->end_mode(Ioss::STATE_DEFINE_MODEL);

  output_->begin_mode(Ioss::STATE_MODEL);
  write_coordinate_list(superElemParts_);
  write_element_connectivity(superElemParts_, subElemIds);
  write_sideset_connectivity(baseParts);
  output_->end_mode(Ioss::STATE_MODEL);
}

void fixup_ghosted_to_shared_nodes(stk::mesh::BulkData & bulk)
{
    stk::mesh::EntityVector ghosted_nodes_that_are_now_shared;
    find_ghosted_nodes_that_need_to_be_shared(bulk, ghosted_nodes_that_are_now_shared);

    stk::CommSparse comm(bulk.parallel());

    for (int phase=0;phase<2;++phase)
    {
        for (size_t i = 0; i < ghosted_nodes_that_are_now_shared.size(); ++i)
        {
            stk::mesh::Entity node = ghosted_nodes_that_are_now_shared[i];
            int proc = bulk.parallel_owner_rank(node);
            comm.send_buffer(proc).pack<stk::mesh::EntityKey>(bulk.entity_key(node));
        }
        if (phase == 0 )
        {
            comm.allocate_buffers();
        }
        else
        {
            comm.communicate();
        }
    }

    stk::mesh::EntityVector sharedNodes;
    for (int process=0;process<bulk.parallel_size();++process)
    {
        while(comm.recv_buffer(process).remaining())
        {
            stk::mesh::EntityKey key;
            comm.recv_buffer(process).unpack<stk::mesh::EntityKey>(key);

            stk::mesh::Entity entity = bulk.get_entity(key);
            if ( bulk.state(entity) != stk::mesh::Deleted && bulk.is_valid(entity) )
            {
                bulk.add_node_sharing(entity, process);
                sharedNodes.push_back(entity);
            }
        }
    }
/////////////////////////

    stk::CommSparse commSecondStage(bulk.parallel());
    for (int phase=0;phase<2;++phase)
    {
        for (size_t i=0;i<sharedNodes.size();++i)
        {
            std::vector<int> procs;
            stk::mesh::EntityKey key = bulk.entity_key(sharedNodes[i]);
            bulk.comm_shared_procs(key, procs);
            for (size_t j=0;j<procs.size();++j)
            {
                if ( procs[j] != bulk.parallel_rank() )
                {
                    commSecondStage.send_buffer(procs[j]).pack<int>(bulk.parallel_rank()).pack<stk::mesh::EntityKey>(key);
                    for (size_t k=0;k<procs.size();++k)
                    {
                        commSecondStage.send_buffer(procs[j]).pack<int>(procs[k]).pack<stk::mesh::EntityKey>(key);
                    }
                }
            }
        }
        if (phase == 0 )
        {
            commSecondStage.allocate_buffers();
        }
        else
        {
            commSecondStage.communicate();
        }
    }

    for (int proc_that_sent_message=0;proc_that_sent_message<bulk.parallel_size();++proc_that_sent_message)
    {
        if ( proc_that_sent_message == bulk.parallel_rank() ) continue;
        while(commSecondStage.recv_buffer(proc_that_sent_message).remaining())
        {
            stk::mesh::EntityKey key;
            int sharingProc;
            commSecondStage.recv_buffer(proc_that_sent_message).unpack<int>(sharingProc).unpack<stk::mesh::EntityKey>(key);
            if ( sharingProc != bulk.parallel_rank() )
            {
                stk::mesh::Entity entity = bulk.get_entity(key);
                if ( bulk.state(entity) != stk::mesh::Deleted && bulk.is_valid(entity) && !bulk.in_shared(key, sharingProc) )
                {
                    bulk.add_node_sharing(entity, sharingProc);
                }
            }
        }
    }
}