C++ Entity::relations方法代码示例

//----------------------------------------------------------------------------
//
// 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 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;
}

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
}

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;
}

    /// modeled after code from Mesquite::IdealWeightMeanRatio::evaluate()
    bool JacobianUtil::operator()(double& m,  PerceptMesh& eMesh, stk::mesh::Entity& element, stk::mesh::FieldBase *coord_field,
                                  const CellTopologyData * topology_data )
    {
      MsqError err;
      static MsqMatrix<3,3> J;

      static Vector3D n;			// Surface normal for 2D objects

      // Prism and Hex element descriptions
      static const int locs_prism[6][4] = {{0, 1, 2, 3}, {1, 2, 0, 4},
                                           {2, 0, 1, 5}, {3, 5, 4, 0},
                                           {4, 3, 5, 1}, {5, 4, 3, 2}};
      static const int locs_hex[8][4] = {{0, 1, 3, 4}, {1, 2, 0, 5},
                                         {2, 3, 1, 6}, {3, 0, 2, 7},
                                         {4, 7, 5, 0}, {5, 4, 6, 1},
                                         {6, 5, 7, 2}, {7, 6, 4, 3}};
      int i=0;

      static const Vector3D d_con(1.0, 1.0, 1.0);

      bool metric_valid = false;
      if (!topology_data) topology_data = stk::percept::PerceptMesh::get_cell_topology(element);

      stk::mesh::PairIterRelation v_i = element.relations(eMesh.node_rank());
      m_num_nodes = v_i.size();

      //#define VERTEX_2D(vi) vector_2D( eMesh.field_data(coord_field, *vi.entity() ) )
      //#define VERTEX_3D(vi) vector_3D( eMesh.field_data(coord_field, *vi.entity() ) )

#define VERTEX_2D(vi) vector_2D( stk::mesh::field_data( *static_cast<const VectorFieldType *>(coord_field) , *vi.entity() ) )
#define VERTEX_3D(vi) vector_3D( stk::mesh::field_data( *static_cast<const VectorFieldType *>(coord_field) , *vi.entity() ) )


      switch(topology_data->key) 
        {
        case shards::Triangle<3>::key:
          n[0] = 0; n[1] = 0; n[2] = 1;
          mCoords[0] = VERTEX_2D(v_i[0]);
          mCoords[1] = VERTEX_2D(v_i[1]);
          mCoords[2] = VERTEX_2D(v_i[2]);
          metric_valid = jacobian_matrix_2D(m, J, mCoords, n, d_con);
          for (i = 0; i < 4; i++) { m_detJ[i] = m; m_J[i] = J; }
          break;
    
        case shards::Quadrilateral<4>::key:
          n[0] = 0; n[1] = 0; n[2] = 1;
          for (i = 0; i < 4; ++i) {
            mCoords[0] = VERTEX_2D(v_i[locs_hex[i][0]]);
            mCoords[1] = VERTEX_2D(v_i[locs_hex[i][1]]);
            mCoords[2] = VERTEX_2D(v_i[locs_hex[i][2]]);
            metric_valid = jacobian_matrix_2D(m_detJ[i], m_J[i], mCoords, n, d_con);
          }
          m = average_metrics(m_detJ, 4, err); MSQ_ERRZERO(err);
          break;

        case shards::Tetrahedron<4>::key:
          mCoords[0] = VERTEX_3D(v_i[0]);
          mCoords[1] = VERTEX_3D(v_i[1]);
          mCoords[2] = VERTEX_3D(v_i[2]);
          mCoords[3] = VERTEX_3D(v_i[3]);
          metric_valid = jacobian_matrix_3D(m, J, mCoords, n, d_con);
          for (i = 0; i < 4; i++) { m_detJ[i] = m; m_J[i] = J; }
          break;

        case shards::Pyramid<5>::key:
          for (i = 0; i < 4; ++i) {
            mCoords[0] = VERTEX_3D(v_i[ i     ]);
            mCoords[1] = VERTEX_3D(v_i[(i+1)%4]);
            mCoords[2] = VERTEX_3D(v_i[(i+3)%4]);
            mCoords[3] = VERTEX_3D(v_i[ 4     ]);
            metric_valid = jacobian_matrix_3D(m_detJ[i], m_J[i], mCoords, n, d_con);
          }
          // FIXME
          m_J[4] = (m_J[0]+m_J[1]+m_J[2]+m_J[3]);
          m_J[4] *= 0.25;
          m_detJ[4] = det(m_J[4]);
          m = average_metrics(m_detJ, 5, err); MSQ_ERRZERO(err);
          break;

        case shards::Wedge<6>::key:
          for (i = 0; i < 6; ++i) {
            mCoords[0] = VERTEX_3D(v_i[locs_prism[i][0]]);
            mCoords[1] = VERTEX_3D(v_i[locs_prism[i][1]]);
            mCoords[2] = VERTEX_3D(v_i[locs_prism[i][2]]);
            mCoords[3] = VERTEX_3D(v_i[locs_prism[i][3]]);
            metric_valid = jacobian_matrix_3D(m_detJ[i], m_J[i], mCoords, n, d_con);
          }
          m = average_metrics(m_detJ, 6, err); MSQ_ERRZERO(err);
          break;

        case shards::Hexahedron<8>::key:
          for (i = 0; i < 8; ++i) {
            mCoords[0] = VERTEX_3D(v_i[locs_hex[i][0]]);
            mCoords[1] = VERTEX_3D(v_i[locs_hex[i][1]]);
            mCoords[2] = VERTEX_3D(v_i[locs_hex[i][2]]);
            mCoords[3] = VERTEX_3D(v_i[locs_hex[i][3]]);
            metric_valid = jacobian_matrix_3D(m_detJ[i], m_J[i], mCoords, n, d_con);
          }
          m = average_metrics(m_detJ, 8, err); MSQ_ERRZERO(err);
//.........这里部分代码省略.........

      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());
              }
//.........这里部分代码省略.........

//----------------------------------------------------------------------------
//
// 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 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
}

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
}

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++;

          }

      }