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

//----------------------------------------------------------------------------
//
// Random fracture criterion function.
//
bool
RandomCriterion::computeFractureCriterion(stk::mesh::Entity& entity, double p) {

  // Fracture only defined on the boundary of the elements
  stk::mesh::EntityRank rank = entity.entity_rank();
  assert(rank == num_dim_ - 1);

  stk::mesh::PairIterRelation neighbor_elems =
    entity.relations(element_rank_);

  // Need an element on each side
  if(neighbor_elems.size() != 2)
    return false;

  bool is_open = false;

  // All we need to do is generate a number between 0 and 1
  double random = 0.5 + 0.5 * Teuchos::ScalarTraits<double>::random();

  if(random < p) {
    is_open = true;
  }

  return is_open;
}

bool
AAdapt::STKUnifUnrefineField::operator()(const stk::mesh::Entity element,
    stk::mesh::FieldBase* field,  const stk::mesh::BulkData& bulkData) {

  const stk::mesh::PairIterRelation elem_nodes = element.relations(stk::topology::NODE_RANK);
  unsigned num_node = elem_nodes.size();
  double* f_data = stk::percept::PerceptMesh::field_data_entity(field, element);
  Albany::AbstractSTKFieldContainer::VectorFieldType* coordField = m_eMesh.get_coordinates_field();

  bool found = true;

  for(unsigned inode = 0; inode < num_node; inode++) {
    stk::mesh::Entity node = * elem_nodes[ inode ].entity();
    double* coord_data = stk::percept::PerceptMesh::field_data(coordField, node);

    if(coord_data[0] < 0.0 || coord_data[1] < 0.0) { // || coord_data[2] > 1.1)
      found = false;
      break;
    }
  }

  if(found)
    f_data[0] = -1.0;

  else
    f_data[0] = 0.0;

  return false;  // don't terminate the loop
}

void SideConnector::connect_side_to_all_elements(stk::mesh::Entity sideEntity, stk::mesh::Entity elemEntity, int elemSide)
{
    connect_side_to_elem(sideEntity, elemEntity, elemSide);
    stk::mesh::impl::LocalId elemLocalId = m_entity_to_local_id[elemEntity.local_offset()];
    connect_side_to_coincident_elements(sideEntity, elemLocalId, elemSide);
    connect_side_to_adjacent_elements(sideEntity, elemLocalId, elemSide);
}

int IEdgeAdapter::markUnrefine(const stk::mesh::Entity& element)
{
    const CellTopologyData * const cell_topo_data = stk::percept::PerceptMesh::get_cell_topology(element);

    CellTopology cell_topo(cell_topo_data);
    const mesh::PairIterRelation elem_nodes = element.relations(stk::mesh::fem::FEMMetaData::NODE_RANK);

    VectorFieldType* coordField = m_eMesh.get_coordinates_field();

    unsigned numSubDimNeededEntities = 0;
    numSubDimNeededEntities = cell_topo_data->edge_count;

    bool unrefAllEdges = true;
    for (unsigned iSubDimOrd = 0; iSubDimOrd < numSubDimNeededEntities; iSubDimOrd++)
    {
        stk::mesh::Entity & node0 = *elem_nodes[cell_topo_data->edge[iSubDimOrd].node[0]].entity();
        stk::mesh::Entity & node1 = *elem_nodes[cell_topo_data->edge[iSubDimOrd].node[1]].entity();
        double * const coord0 = stk::mesh::field_data( *coordField , node0 );
        double * const coord1 = stk::mesh::field_data( *coordField , node1 );

        int markInfo = mark(element, iSubDimOrd, node0, node1, coord0, coord1, 0);
        bool do_unref = markInfo & DO_UNREFINE;
        if (!do_unref)
        {
            unrefAllEdges = false;
            break;
        }
    }
    if (unrefAllEdges)
        return -1;
    else
        return 0;
}

bool gather_field_data(unsigned expected_num_rel, const field_type & field,
  const stk::mesh::Entity & entity,
  typename stk::mesh::FieldTraits<field_type>::data_type * dst,
  stk::mesh::EntityRank entity_rank, const int dim)
{
  typedef typename stk::mesh::FieldTraits<field_type>::data_type T;

  stk::mesh::PairIterRelation rel = entity.relations(entity_rank);

  bool result = expected_num_rel == (unsigned) rel.size();

  if (result)
  {
    // Iterate over field data for each related entity and copy data
    // into src for one entity at a time
    T * const dst_end = dst + dim * expected_num_rel;
    for (; dst < dst_end; ++rel, dst += dim)
    {
      const T* src = field_data(field, *rel->entity());
      if (!src)
      {
        break;
      }
      // FIXME:
      if (dim == 1) stk::Copy<1>(dst, src);
      if (dim == 2) stk::Copy<2>(dst, src);
      if (dim == 3) stk::Copy<3>(dst, src);
    }
    result = dst == dst_end;
  }
  return result;
}

 KOKKOS_INLINE_FUNCTION unsigned get_bucket_id(stk::mesh::Entity entity) const
 {
     const stk::mesh::FastMeshIndex& meshIndex = constMeshIndices(entity.local_offset());
     return meshIndex.bucket_id;
 }

 STK_FUNCTION unsigned get_index(stk::mesh::Entity entity) const
 {
     const stk::mesh::FastMeshIndex& meshIndex = constMeshIndices(entity.local_offset());
     return meshIndex.bucket_id * bucketCapacity + meshIndex.bucket_ord;
 }

void TestLocalRefinerTri2::
refineMethodApply(NodeRegistry::ElementFunctionPrototype function, const stk::mesh::Entity& element, vector<NeededEntityType>& needed_entity_ranks)
{
    const CellTopologyData * const cell_topo_data = stk::percept::PerceptMesh::get_cell_topology(element);

    CellTopology cell_topo(cell_topo_data);
    const mesh::PairIterRelation elem_nodes = element.relations(stk::mesh::fem::FEMMetaData::NODE_RANK);

    VectorFieldType* coordField = m_eMesh.get_coordinates_field();

    for (unsigned ineed_ent=0; ineed_ent < needed_entity_ranks.size(); ineed_ent++)
    {
        unsigned numSubDimNeededEntities = 0;
        stk::mesh::EntityRank needed_entity_rank = needed_entity_ranks[ineed_ent].first;

        if (needed_entity_rank == m_eMesh.edge_rank())
        {
            numSubDimNeededEntities = cell_topo_data->edge_count;
        }
        else if (needed_entity_rank == m_eMesh.face_rank())
        {
            numSubDimNeededEntities = cell_topo_data->side_count;
        }
        else if (needed_entity_rank == m_eMesh.element_rank())
        {
            numSubDimNeededEntities = 1;
        }

        // see how many edges are already marked
        int num_marked=0;
        if (needed_entity_rank == m_eMesh.edge_rank())
        {
            for (unsigned iSubDimOrd = 0; iSubDimOrd < numSubDimNeededEntities; iSubDimOrd++)
            {
                bool is_empty = m_nodeRegistry->is_empty( element, needed_entity_rank, iSubDimOrd);
                if (!is_empty) ++num_marked;
            }
        }

        for (unsigned iSubDimOrd = 0; iSubDimOrd < numSubDimNeededEntities; iSubDimOrd++)
        {
            /// note: at this level of granularity we can do single edge refinement, hanging nodes, etc.
            //SubDimCell_SDSEntityType subDimEntity;
            //getSubDimEntity(subDimEntity, element, needed_entity_rank, iSubDimOrd);
            //bool is_empty = m_nodeRegistry->is_empty( element, needed_entity_rank, iSubDimOrd);
            //if(1||!is_empty)

            if (needed_entity_rank == m_eMesh.edge_rank())
            {
                stk::mesh::Entity & node0 = *elem_nodes[cell_topo_data->edge[iSubDimOrd].node[0]].entity();
                stk::mesh::Entity & node1 = *elem_nodes[cell_topo_data->edge[iSubDimOrd].node[1]].entity();
                double * const coord0 = stk::mesh::field_data( *coordField , node0 );
                double * const coord1 = stk::mesh::field_data( *coordField , node1 );

                // only refine diagonals of the split quads
                if (m_diagonals)
                {
                    if ( std::abs(coord0[0]-coord1[0]) > 1.e-3 && std::abs(coord0[1]-coord1[1]) > 1.e-3 )
                    {
                        (m_nodeRegistry ->* function)(element, needed_entity_ranks[ineed_ent], iSubDimOrd, true);
                    }
                }
                else
                {
                    if ( std::abs(coord0[0]-coord1[0]) < 1.e-3 && std::abs(coord0[1]-coord1[1]) > 1.e-3 )
                    {
                        (m_nodeRegistry ->* function)(element, needed_entity_ranks[ineed_ent], iSubDimOrd, true);
                    }
                }
            }

        } // iSubDimOrd
    } // ineed_ent
}

inline int Albany::STKDiscretization::gid(const stk::mesh::Entity& node) const
{ return node.identifier()-1; }

unsigned
Albany::STKDiscretization::determine_local_side_id( const stk::mesh::Entity & elem , stk::mesh::Entity & side ) {

  using namespace stk;

  const CellTopologyData * const elem_top = mesh::fem::get_cell_topology( elem ).getCellTopologyData();

  const mesh::PairIterRelation elem_nodes = elem.relations( mesh::fem::FEMMetaData::NODE_RANK );
  const mesh::PairIterRelation side_nodes = side.relations( mesh::fem::FEMMetaData::NODE_RANK );

  int side_id = -1 ;

  if(elem_nodes.size() == 0 || side_nodes.size() == 0){ // Node relations are not present, look at elem->face

    int elem_rank = elem.entity_rank();
    const mesh::PairIterRelation elem_sides = elem.relations( elem_rank - 1);

    for ( unsigned i = 0 ; i < elem_sides.size() ; ++i ) {

      const stk::mesh::Entity & elem_side = *elem_sides[i].entity();

      if(elem_side.identifier() == side.identifier()){ // Found the local side in the element

         side_id = static_cast<int>(i);

         return side_id;

      }

    }

    if ( side_id < 0 ) {
      std::ostringstream msg ;
      msg << "determine_local_side_id( " ;
      msg << elem_top->name ;
      msg << " , Element[ " ;
      msg << elem.identifier();
      msg << " ]{" ;
      for ( unsigned i = 0 ; i < elem_sides.size() ; ++i ) {
        msg << " " << elem_sides[i].entity()->identifier();
      }
      msg << " } , Side[ " ;
      msg << side.identifier();
      msg << " ] ) FAILED" ;
      throw std::runtime_error( msg.str() );
    }

  }
  else { // Conventional elem->node - side->node connectivity present

    for ( unsigned i = 0 ; side_id == -1 && i < elem_top->side_count ; ++i ) {
      const CellTopologyData & side_top = * elem_top->side[i].topology ;
      const unsigned     * side_map =   elem_top->side[i].node ;

      if ( side_nodes.size() == side_top.node_count ) {

        side_id = i ;

        for ( unsigned j = 0 ;
              side_id == static_cast<int>(i) && j < side_top.node_count ; ++j ) {

          mesh::Entity * const elem_node = elem_nodes[ side_map[j] ].entity();

          bool found = false ;

          for ( unsigned k = 0 ; ! found && k < side_top.node_count ; ++k ) {
            found = elem_node == side_nodes[k].entity();
          }

          if ( ! found ) { side_id = -1 ; }
        }
      }
    }

    if ( side_id < 0 ) {
      std::ostringstream msg ;
      msg << "determine_local_side_id( " ;
      msg << elem_top->name ;
      msg << " , Element[ " ;
      msg << elem.identifier();
      msg << " ]{" ;
      for ( unsigned i = 0 ; i < elem_nodes.size() ; ++i ) {
        msg << " " << elem_nodes[i].entity()->identifier();
      }
      msg << " } , Side[ " ;
      msg << side.identifier();
      msg << " ]{" ;
      for ( unsigned i = 0 ; i < side_nodes.size() ; ++i ) {
        msg << " " << side_nodes[i].entity()->identifier();
      }
      msg << " } ) FAILED" ;
      throw std::runtime_error( msg.str() );
    }
  }

  return static_cast<unsigned>(side_id) ;
}

    void TestLocalRefinerTet_N_2_1::
    refineMethodApply(NodeRegistry::ElementFunctionPrototype function, const stk::mesh::Entity& element, vector<NeededEntityType>& needed_entity_ranks)
    {
      //static int n_seq = 400;

      const CellTopologyData * const cell_topo_data = stk::percept::PerceptMesh::get_cell_topology(element);
                
      CellTopology cell_topo(cell_topo_data);
      const mesh::PairIterRelation elem_nodes = element.relations(stk::mesh::fem::FEMMetaData::NODE_RANK);

      //VectorFieldType* coordField = m_eMesh.get_coordinates_field();

      for (unsigned ineed_ent=0; ineed_ent < needed_entity_ranks.size(); ineed_ent++)
        {
          unsigned numSubDimNeededEntities = 0;
          stk::mesh::EntityRank needed_entity_rank = needed_entity_ranks[ineed_ent].first;

          if (needed_entity_rank == m_eMesh.edge_rank())
            {
              numSubDimNeededEntities = cell_topo_data->edge_count;
            }
          else if (needed_entity_rank == m_eMesh.face_rank())
            {
              numSubDimNeededEntities = cell_topo_data->side_count;
            }
          else if (needed_entity_rank == m_eMesh.element_rank())
            {
              numSubDimNeededEntities = 1;
            }

          // see how many edges are already marked
          int num_marked=0;
          if (needed_entity_rank == m_eMesh.edge_rank())
            {
              for (unsigned iSubDimOrd = 0; iSubDimOrd < numSubDimNeededEntities; iSubDimOrd++)
                {
                  bool is_empty = m_nodeRegistry->is_empty( element, needed_entity_rank, iSubDimOrd);
                  if (!is_empty) ++num_marked;
                }
            }

          for (unsigned iSubDimOrd = 0; iSubDimOrd < numSubDimNeededEntities; iSubDimOrd++)
            {
              /// note: at this level of granularity we can do single edge refinement, hanging nodes, etc.
              //SubDimCell_SDSEntityType subDimEntity;
              //getSubDimEntity(subDimEntity, element, needed_entity_rank, iSubDimOrd);
              //bool is_empty = m_nodeRegistry->is_empty( element, needed_entity_rank, iSubDimOrd);
              //if(1||!is_empty)

              if (needed_entity_rank == m_eMesh.edge_rank())
                {
#if 0
                  stk::mesh::Entity & node0 = *elem_nodes[cell_topo_data->edge[iSubDimOrd].node[0]].entity();
                  stk::mesh::Entity & node1 = *elem_nodes[cell_topo_data->edge[iSubDimOrd].node[1]].entity();
                  double * const coord0 = stk::mesh::field_data( *coordField , node0 );
                  double * const coord1 = stk::mesh::field_data( *coordField , node1 );
                  
                  // vertical line position
                  const double vx = 0.21;

                  // horizontal line position
                  const double vy = 1.21;

                  // choose to refine or not 
                  if (
                      ( std::fabs(coord0[0]-coord1[0]) > 1.e-3 &&
                        ( (coord0[0] < vx && vx < coord1[0]) || (coord1[0] < vx && vx < coord0[0]) )
                        )
                      ||
                      ( std::fabs(coord0[1]-coord1[1]) > 1.e-3 &&
                        ( (coord0[1] < vy && vy < coord1[1]) || (coord1[1] < vy && vy < coord0[1]) )
                        )
                      )
                    {
                      (m_nodeRegistry ->* function)(element, needed_entity_ranks[ineed_ent], iSubDimOrd, true);
                    }

#endif
                  // mark first m_edge_mark_bitcode edges 
                  std::cout << "tmp TestLocalRefinerTet_N_2_1 element.identifier() = " << element.identifier() 
                            << " m_edge_mark_bitcode= " << m_edge_mark_bitcode << std::endl;

                  if ( ((1 << iSubDimOrd) & m_edge_mark_bitcode) && element.identifier() == 1)
                    (m_nodeRegistry ->* function)(element, needed_entity_ranks[ineed_ent], iSubDimOrd, true);

                }

            } // iSubDimOrd
        } // ineed_ent
    }

      void 
      createNewElements(percept::PerceptMesh& eMesh, NodeRegistry& nodeRegistry, 
                        stk::mesh::Entity& element,  NewSubEntityNodesType& new_sub_entity_nodes, vector<stk::mesh::Entity *>::iterator& element_pool,
                        stk::mesh::FieldBase *proc_rank_field=0)
      {
        const CellTopologyData * const cell_topo_data = stk::percept::PerceptMesh::get_cell_topology(element);
        typedef boost::tuple<stk::mesh::EntityId, stk::mesh::EntityId, stk::mesh::EntityId, stk::mesh::EntityId> quad_tuple_type;
        static vector<quad_tuple_type> elems(4);

        CellTopology cell_topo(cell_topo_data);
        const stk::mesh::PairIterRelation elem_nodes = element.relations(stk::mesh::fem::FEMMetaData::NODE_RANK);

        //stk::mesh::Part & active = mesh->ActivePart();
        //stk::mesh::Part & quad4  = mesh->QuadPart();

        std::vector<stk::mesh::Part*> add_parts;
        std::vector<stk::mesh::Part*> remove_parts;

        //add_parts.push_back( &active );
        //FIXME 
        //add_parts.push_back( const_cast<mesh::Part*>( eMesh.getPart(m_toTopoPartName) ));
        add_parts = m_toParts;
        
        double tmp_x[3];
        for (int iedge = 0; iedge < 4; iedge++)
          {
            double * mp = midPoint(EDGE_COORD(iedge,0), EDGE_COORD(iedge,1), eMesh.get_spatial_dim(), tmp_x);

            if (!EDGE_N(iedge))
              {
                std::cout << "P[" << eMesh.get_rank() << " nid ## = 0 << " << std::endl;
              }
            eMesh.createOrGetNode(EDGE_N(iedge), mp);

          }

        nodeRegistry.makeCentroidCoords(*const_cast<stk::mesh::Entity *>(&element), m_eMesh.element_rank(), 0u);


// new_sub_entity_nodes[i][j]
#define CENTROID_N NN(m_primaryEntityRank,0)  

        elems[0] = quad_tuple_type(VERT_N(0), EDGE_N(0), CENTROID_N, EDGE_N(3));
        elems[1] = quad_tuple_type(VERT_N(1), EDGE_N(1), CENTROID_N, EDGE_N(0));
        elems[2] = quad_tuple_type(VERT_N(2), EDGE_N(2), CENTROID_N, EDGE_N(1));
        elems[3] = quad_tuple_type(VERT_N(3), EDGE_N(3), CENTROID_N, EDGE_N(2));

#undef CENTROID_N

        // write a diagram of the refinement pattern as a vtk file, or a latex/tikz/pgf file
#define WRITE_DIAGRAM 0
#if WRITE_DIAGRAM

#endif
        
        for (unsigned ielem=0; ielem < elems.size(); ielem++)
          {
            stk::mesh::Entity& newElement = *(*element_pool);

            if (proc_rank_field)
              {
                double *fdata = stk::mesh::field_data( *static_cast<const ScalarFieldType *>(proc_rank_field) , newElement );
                //fdata[0] = double(m_eMesh.get_rank());
                fdata[0] = double(newElement.owner_rank());
              }

            eMesh.get_bulk_data()->change_entity_parts( newElement, add_parts, remove_parts );

            {
              if (!elems[ielem].get<0>())
                {
                  std::cout << "P[" << eMesh.get_rank() << " nid = 0 << " << std::endl;
                  exit(1);
                }

            }
            eMesh.get_bulk_data()->declare_relation(newElement, eMesh.createOrGetNode(elems[ielem].get<0>()), 0);
            eMesh.get_bulk_data()->declare_relation(newElement, eMesh.createOrGetNode(elems[ielem].get<1>()), 1);
            eMesh.get_bulk_data()->declare_relation(newElement, eMesh.createOrGetNode(elems[ielem].get<2>()), 2);
            eMesh.get_bulk_data()->declare_relation(newElement, eMesh.createOrGetNode(elems[ielem].get<3>()), 3);

            set_parent_child_relations(eMesh, element, newElement, ielem);


            element_pool++;

          }

      }

//----------------------------------------------------------------------------
//
// Stress fracture criterion function.
//
bool
StressFracture::computeFractureCriterion(stk::mesh::Entity entity, double p) {
  // Fracture only defined on the boundary of the elements
  stk::mesh::EntityRank rank = entity.entity_rank();
  assert(rank == num_dim_ - 1);

  stk::mesh::PairIterRelation neighbor_elems =
    entity.relations(element_rank_);

  // Need an element on each side of the edge
  if(neighbor_elems.size() != 2)
    return false;

  // Note that these are element GIDs

  stk::mesh::EntityId elem_0_Id =
    neighbor_elems[0].entity()->identifier();
  stk::mesh::EntityId elem_1_Id =
    neighbor_elems[1].entity()->identifier();

  Albany::WsLIDList& elemGIDws = stk_.getElemGIDws();

  // Have two elements, one on each size of the edge (or
  // face). Check and see if the stresses are such that we want to
  // split the mesh here.
  //
  // Initial hack - GAH: if the average stress between two elements
  // is above the input value "Fracture Stress", split them at the
  // edge

  bool is_open = false;

  // Check criterion
  // If average between cells is above crit, split
  //  if (0.5 * (avg_stresses[elemGIDws[elem_0_Id].ws][elemGIDws[elem_0_Id].LID] +
  //    avg_stresses[elemGIDws[elem_1_Id].ws][elemGIDws[elem_1_Id].LID]) >= crit_stress){
  // If stress difference across face it above crit, split
  //  if (fabs(avg_stresses[elemGIDws[elem_0_Id].ws][elemGIDws[elem_0_Id].LID] -
  //    avg_stresses[elemGIDws[elem_1_Id].ws][elemGIDws[elem_1_Id].LID]) >= crit_stress){
  // Just split the doggone elements already!!!
  if(p == 5) {
    if((elem_0_Id - 1 == 35 && elem_1_Id - 1 == 140) ||
        (elem_1_Id - 1 == 35 && elem_0_Id - 1 == 140)) {

      is_open = true;

      std::cout << "Splitting elements " << elem_0_Id - 1 << " and " << elem_1_Id - 1 << std::endl;
      //std::cout << avg_stresses[elemGIDws[elem_0_Id].ws][elemGIDws[elem_0_Id].LID] << " " <<
      //    avg_stresses[elemGIDws[elem_1_Id].ws][elemGIDws[elem_1_Id].LID] << std::endl;

    }
  }

  else if(p == 10) {
    if((elem_0_Id - 1 == 42 && elem_1_Id - 1 == 147) ||
        (elem_1_Id - 1 == 42 && elem_0_Id - 1 == 147)) {

      is_open = true;

      std::cout << "Splitting elements " << elem_0_Id - 1 << " and " << elem_1_Id - 1 << std::endl;
    }
  }

  else if(p == 15) {
    if((elem_0_Id - 1 == 49 && elem_1_Id - 1 == 154) ||
        (elem_1_Id - 1 == 49 && elem_0_Id - 1 == 154)) {

      is_open = true;

      std::cout << "Splitting elements " << elem_0_Id - 1 << " and " << elem_1_Id - 1 << std::endl;
    }
  }

  return is_open;
}

      void 
      createNewElements(percept::PerceptMesh& eMesh, NodeRegistry& nodeRegistry, 
                        stk::mesh::Entity& element,  NewSubEntityNodesType& new_sub_entity_nodes, vector<stk::mesh::Entity *>::iterator& element_pool,
                        stk::mesh::FieldBase *proc_rank_field=0)
      {
        const CellTopologyData * const cell_topo_data = stk::percept::PerceptMesh::get_cell_topology(element);
        typedef boost::tuple<stk::mesh::EntityId, stk::mesh::EntityId> line_tuple_type;
        static vector<line_tuple_type> elems(2);

        CellTopology cell_topo(cell_topo_data);
        const stk::mesh::PairIterRelation elem_nodes = element.relations(stk::mesh::fem::FEMMetaData::NODE_RANK);

        std::vector<stk::mesh::Part*> add_parts;
        std::vector<stk::mesh::Part*> remove_parts;

        add_parts = m_toParts;
        
        unsigned num_nodes_on_edge = new_sub_entity_nodes[m_eMesh.edge_rank()][0].size();
        if (!num_nodes_on_edge)
          return;

        double coord_x[3];
        for (int iedge = 0; iedge < 1; iedge++)
          {
            //double * mp = midPoint(EDGE_COORD(iedge,0), EDGE_COORD(iedge,1), eMesh.get_spatial_dim(), coord_x);
            //double * mp = midPoint(FACE_COORD(iedge,0), FACE_COORD(iedge,1), eMesh.get_spatial_dim(), coord_x);
            double * mp = midPoint(VERT_COORD(0), VERT_COORD(1), eMesh.get_spatial_dim(), coord_x);

            if (!EDGE_N(iedge))
              {
                std::cout << "P[" << eMesh.get_rank() << " nid ## = 0  " << std::endl;
              }

            eMesh.createOrGetNode(EDGE_N(iedge), mp);
          }

        // FIXME
        nodeRegistry.makeCentroidCoords(*const_cast<stk::mesh::Entity *>(&element), m_primaryEntityRank, 0u);
        nodeRegistry.addToExistingParts(*const_cast<stk::mesh::Entity *>(&element), m_primaryEntityRank, 0u);

        nodeRegistry.interpolateFields(*const_cast<stk::mesh::Entity *>(&element), m_primaryEntityRank, 0u);

        Elem::CellTopology elem_celltopo = Elem::getCellTopology< FromTopology >();
        const Elem::RefinementTopology* ref_topo_p = Elem::getRefinementTopology(elem_celltopo);
        const Elem::RefinementTopology& ref_topo = *ref_topo_p;

#ifndef NDEBUG
        unsigned num_child = ref_topo.num_child();
        VERIFY_OP(num_child, == , 2, "createNewElements num_child problem");
        bool homogeneous_child = ref_topo.homogeneous_child();
        VERIFY_OP(homogeneous_child, ==, true, "createNewElements homogeneous_child");
#endif

        // new_sub_entity_nodes[i][j]
        //const UInt * const * child_nodes() const {
        //const UInt * child_node_0 = ref_topo.child_node(0);

        typedef Elem::StdMeshObjTopologies::RefTopoX RefTopoX;
        RefTopoX& l2 = Elem::StdMeshObjTopologies::RefinementTopologyExtra< FromTopology > ::refinement_topology;

#define CENTROID_N NN(m_primaryEntityRank,0)  

        for (unsigned iChild = 0; iChild < 2; iChild++)
          {
            unsigned EN[2];
            for (unsigned jNode = 0; jNode < 2; jNode++)
              {
                unsigned childNodeIdx = ref_topo.child_node(iChild)[jNode];
#ifndef NDEBUG
                unsigned childNodeIdxCheck = l2[childNodeIdx].ordinal_of_node;
                VERIFY_OP(childNodeIdx, ==, childNodeIdxCheck, "childNodeIdxCheck");
#endif
                unsigned inode=0;

                if (l2[childNodeIdx].rank_of_subcell == 0)
                  inode = VERT_N(l2[childNodeIdx].ordinal_of_subcell);
                else if (l2[childNodeIdx].rank_of_subcell == 1)
                  inode = EDGE_N(l2[childNodeIdx].ordinal_of_subcell);

                //                 else if (l2[childNodeIdx].rank_of_subcell == 2)
                //                   inode = CENTROID_N;

                EN[jNode] = inode;
              }
            elems[iChild] = line_tuple_type(EN[0], EN[1]);
          }

#undef CENTROID_N

        for (unsigned ielem=0; ielem < elems.size(); ielem++)
          {
            stk::mesh::Entity& newElement = *(*element_pool);

#if 0
            if (proc_rank_field && proc_rank_field->rank() == m_eMesh.edge_rank()) //&& m_eMesh.get_spatial_dim()==1)
              {
                double *fdata = stk::mesh::field_data( *static_cast<const ScalarFieldType *>(proc_rank_field) , newElement );
                //fdata[0] = double(m_eMesh.get_rank());
                fdata[0] = double(newElement.owner_rank());
              }
//.........这里部分代码省略.........

void IEdgeAdapter::
refineMethodApply(NodeRegistry::ElementFunctionPrototype function, const stk::mesh::Entity& element,
                  vector<NeededEntityType>& needed_entity_ranks)
{
    const CellTopologyData * const cell_topo_data = stk::percept::PerceptMesh::get_cell_topology(element);

    CellTopology cell_topo(cell_topo_data);
    const mesh::PairIterRelation elem_nodes = element.relations(stk::mesh::fem::FEMMetaData::NODE_RANK);

    VectorFieldType* coordField = m_eMesh.get_coordinates_field();

    for (unsigned ineed_ent=0; ineed_ent < needed_entity_ranks.size(); ineed_ent++)
    {
        unsigned numSubDimNeededEntities = 0;
        stk::mesh::EntityRank needed_entity_rank = needed_entity_ranks[ineed_ent].first;

        if (needed_entity_rank == m_eMesh.edge_rank())
        {
            numSubDimNeededEntities = cell_topo_data->edge_count;
        }
        else if (needed_entity_rank == m_eMesh.face_rank())
        {
            numSubDimNeededEntities = cell_topo_data->side_count;
            throw std::runtime_error("IEdgeAdapter::apply can't use IEdgeAdapter for RefinerPatterns that require face nodes");
        }
        else if (needed_entity_rank == m_eMesh.element_rank())
        {
            numSubDimNeededEntities = 1;
            throw std::runtime_error("IEdgeAdapter::apply can't use IEdgeAdapter for RefinerPatterns that require volume nodes");
        }

        // see how many edges are already marked
        int num_marked=0;
        std::vector<int> edge_marks(numSubDimNeededEntities,0);
        if (needed_entity_rank == m_eMesh.edge_rank())
        {
            for (unsigned iSubDimOrd = 0; iSubDimOrd < numSubDimNeededEntities; iSubDimOrd++)
            {
                bool is_empty = m_nodeRegistry->is_empty( element, needed_entity_rank, iSubDimOrd);
                if (!is_empty)
                {
                    edge_marks[iSubDimOrd] = 1;
                    ++num_marked;
                }
            }
        }

        for (unsigned iSubDimOrd = 0; iSubDimOrd < numSubDimNeededEntities; iSubDimOrd++)
        {
            if (needed_entity_rank == m_eMesh.edge_rank())
            {
                stk::mesh::Entity & node0 = *elem_nodes[cell_topo_data->edge[iSubDimOrd].node[0]].entity();
                stk::mesh::Entity & node1 = *elem_nodes[cell_topo_data->edge[iSubDimOrd].node[1]].entity();
                double * const coord0 = stk::mesh::field_data( *coordField , node0 );
                double * const coord1 = stk::mesh::field_data( *coordField , node1 );


                int markInfo = mark(element, iSubDimOrd, node0, node1, coord0, coord1, &edge_marks);

                bool needNodes = (DO_REFINE & markInfo);
                {
                    (m_nodeRegistry ->* function)(element, needed_entity_ranks[ineed_ent], iSubDimOrd, needNodes);
                }
            }

        } // iSubDimOrd
    } // ineed_ent
}