MOAB: Mesh Oriented datABase  (version 5.2.1)
HalfFacetRep.cpp
Go to the documentation of this file.
00001 /**
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 #include "moab/HalfFacetRep.hpp"
00017 #include "Internals.hpp"
00018 #include <iostream>
00019 #include <assert.h>
00020 #include <vector>
00021 #include <map>
00022 #include "MBTagConventions.hpp"
00023 #include "moab/ScdInterface.hpp"
00024 #ifdef MOAB_HAVE_MPI
00025 #include "moab/ParallelComm.hpp"
00026 #endif
00027 
00028 namespace moab
00029 {
00030 
00031 inline EntityHandle CREATE_HALFFACET( const unsigned lid, const EntityID id )
00032 {
00033     assert( id <= MB_END_ID && lid < 6 );
00034     return ( ( (HFacet)lid ) << MB_ID_WIDTH ) | id;
00035 }
00036 inline EntityID FID_FROM_HALFFACET( HFacet handle )
00037 {
00038     return ( handle & MB_ID_MASK );
00039 }
00040 inline int LID_FROM_HALFFACET( HFacet handle )
00041 {
00042     return static_cast< int >( handle >> MB_ID_WIDTH );
00043 }
00044 
00045 HalfFacetRep::HalfFacetRep( Core* impl, ParallelComm* comm, moab::EntityHandle rset, bool filter_ghosts )
00046     : thismeshtype( CURVE ), mb( impl ), pcomm( comm ), _rset( rset ), _filterghost( filter_ghosts )
00047 {
00048     assert( NULL != impl );
00049     mInitAHFmaps = false;
00050     chk_mixed    = false;
00051     is_mixed     = false;
00052 }
00053 
00054 HalfFacetRep::~HalfFacetRep() {}
00055 
00056 MESHTYPE HalfFacetRep::get_mesh_type( int nverts, int nedges, int nfaces, int ncells )
00057 {
00058     MESHTYPE mesh_type = CURVE;
00059 
00060     if( nverts && nedges && ( !nfaces ) && ( !ncells ) )
00061         mesh_type = CURVE;
00062     else if( nverts && !nedges && nfaces && !ncells )
00063         mesh_type = SURFACE;
00064     else if( nverts && nedges && nfaces && !ncells )
00065         mesh_type = SURFACE_MIXED;
00066     else if( nverts && !nedges && !nfaces && ncells )
00067         mesh_type = VOLUME;
00068     else if( nverts && nedges && !nfaces && ncells )
00069         mesh_type = VOLUME_MIXED_1;
00070     else if( nverts && !nedges && nfaces && ncells )
00071         mesh_type = VOLUME_MIXED_2;
00072     else if( nverts && nedges && nfaces && ncells )
00073         mesh_type = VOLUME_MIXED;
00074 
00075     return mesh_type;
00076 }
00077 
00078 const HalfFacetRep::adj_matrix HalfFacetRep::adjMatrix[7] = {
00079     // Stores the adjacency matrix for each mesh type.
00080     // CURVE
00081     { { { 0, 1, 0, 0 }, { 1, 1, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
00082 
00083     // SURFACE
00084     { { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 1, 0 }, { 0, 0, 0, 0 } } },
00085 
00086     // SURFACE_MIXED
00087     { { { 0, 1, 1, 0 }, { 1, 1, 1, 0 }, { 1, 1, 1, 0 }, { 0, 0, 0, 0 } } },
00088 
00089     // VOLUME
00090     { { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 1 } } },
00091 
00092     // VOLUME_MIXED_1
00093     { { { 0, 1, 0, 1 }, { 1, 1, 0, 1 }, { 0, 0, 0, 0 }, { 1, 1, 0, 1 } } },
00094 
00095     // VOLUME_MIXED_2
00096     { { { 0, 0, 1, 1 }, { 0, 0, 0, 0 }, { 1, 0, 1, 1 }, { 1, 0, 1, 1 } } },
00097 
00098     // VOLUME_MIXED
00099     { { { 0, 1, 1, 1 }, { 1, 1, 1, 1 }, { 1, 1, 1, 1 }, { 1, 1, 1, 1 } } }
00100 };
00101 
00102 int HalfFacetRep::get_index_for_meshtype( MESHTYPE mesh_type )
00103 {
00104     int index = 0;
00105     switch( mesh_type )
00106     {
00107         case CURVE:
00108             index = 0;
00109             break;
00110         case SURFACE:
00111             index = 1;
00112             break;
00113         case SURFACE_MIXED:
00114             index = 2;
00115             break;
00116         case VOLUME:
00117             index = 3;
00118             break;
00119         case VOLUME_MIXED_1:
00120             index = 4;
00121             break;
00122         case VOLUME_MIXED_2:
00123             index = 5;
00124             break;
00125         case VOLUME_MIXED:
00126             index = 6;
00127             break;
00128     }
00129     return index;
00130 }
00131 
00132 bool HalfFacetRep::check_mixed_entity_type()
00133 {
00134     if( !chk_mixed )
00135     {
00136         chk_mixed = true;
00137 
00138         ErrorCode error;
00139         Range felems, celems;
00140 
00141         error = mb->get_entities_by_dimension( this->_rset, 2, felems );MB_CHK_ERR( error );
00142 
00143         if( felems.size() )
00144         {
00145             Range tris, quad, poly;
00146             tris = felems.subset_by_type( MBTRI );
00147             quad = felems.subset_by_type( MBQUAD );
00148             poly = felems.subset_by_type( MBPOLYGON );
00149             if( ( tris.size() && quad.size() ) || ( tris.size() && poly.size() ) || ( quad.size() && poly.size() ) )
00150                 is_mixed = true;
00151             if( poly.size() ) is_mixed = true;
00152 
00153             if( is_mixed ) return is_mixed;
00154         }
00155 
00156         error = mb->get_entities_by_dimension( this->_rset, 3, celems );MB_CHK_ERR( error );
00157         if( celems.size() )
00158         {
00159             Range tet, pyr, prism, hex, polyhed;
00160             tet     = celems.subset_by_type( MBTET );
00161             pyr     = celems.subset_by_type( MBPYRAMID );
00162             prism   = celems.subset_by_type( MBPRISM );
00163             hex     = celems.subset_by_type( MBHEX );
00164             polyhed = celems.subset_by_type( MBPOLYHEDRON );
00165             if( ( tet.size() && pyr.size() ) || ( tet.size() && prism.size() ) || ( tet.size() && hex.size() ) ||
00166                 ( tet.size() && polyhed.size() ) || ( pyr.size() && prism.size() ) || ( pyr.size() && hex.size() ) ||
00167                 ( pyr.size() && polyhed.size() ) || ( prism.size() && hex.size() ) ||
00168                 ( prism.size() && polyhed.size() ) || ( hex.size() && polyhed.size() ) )
00169                 is_mixed = true;
00170 
00171             if( polyhed.size() ) is_mixed = true;
00172         }
00173 
00174         ScdInterface* scdi = NULL;
00175         error              = mb->query_interface( scdi );MB_CHK_ERR( error );
00176         if( scdi )
00177         {
00178             Range boxes;
00179             error = scdi->find_boxes( boxes );MB_CHK_ERR( error );
00180 
00181             if( !boxes.empty() ) is_mixed = true;
00182         }
00183     }
00184     return is_mixed;
00185 }
00186 
00187 /*******************************************************
00188  * initialize                                          *
00189  ******************************************************/
00190 
00191 ErrorCode HalfFacetRep::initialize()
00192 {
00193     ErrorCode error;
00194 
00195     if( !mInitAHFmaps )
00196     {
00197         mInitAHFmaps = true;
00198 #ifdef MOAB_HAVE_MPI
00199         if( pcomm && _filterghost )
00200         {
00201             moab::Range _averts, _aedgs, _afacs, _acels;
00202             error = mb->get_entities_by_dimension( this->_rset, 0, _averts, true );MB_CHK_ERR( error );
00203             error = mb->get_entities_by_dimension( this->_rset, 1, _aedgs, true );MB_CHK_ERR( error );
00204             error = mb->get_entities_by_dimension( this->_rset, 2, _afacs, true );MB_CHK_ERR( error );
00205             error = mb->get_entities_by_dimension( this->_rset, 3, _acels, true );MB_CHK_ERR( error );
00206 
00207             // filter based on parallel status
00208             error = pcomm->filter_pstatus( _averts, PSTATUS_GHOST, PSTATUS_NOT, -1, &_verts );MB_CHK_ERR( error );
00209             error = pcomm->filter_pstatus( _aedgs, PSTATUS_GHOST, PSTATUS_NOT, -1, &_edges );MB_CHK_ERR( error );
00210             error = pcomm->filter_pstatus( _afacs, PSTATUS_GHOST, PSTATUS_NOT, -1, &_faces );MB_CHK_ERR( error );
00211             error = pcomm->filter_pstatus( _acels, PSTATUS_GHOST, PSTATUS_NOT, -1, &_cells );MB_CHK_ERR( error );
00212         }
00213         else
00214         {
00215             error = mb->get_entities_by_dimension( this->_rset, 0, _verts, true );MB_CHK_ERR( error );
00216             error = mb->get_entities_by_dimension( this->_rset, 1, _edges, true );MB_CHK_ERR( error );
00217             error = mb->get_entities_by_dimension( this->_rset, 2, _faces, true );MB_CHK_ERR( error );
00218             error = mb->get_entities_by_dimension( this->_rset, 3, _cells, true );MB_CHK_ERR( error );
00219         }
00220 #else
00221         error = mb->get_entities_by_dimension( this->_rset, 0, _verts, true );MB_CHK_ERR( error );
00222         error = mb->get_entities_by_dimension( this->_rset, 1, _edges, true );MB_CHK_ERR( error );
00223         error = mb->get_entities_by_dimension( this->_rset, 2, _faces, true );MB_CHK_ERR( error );
00224         error = mb->get_entities_by_dimension( this->_rset, 3, _cells, true );MB_CHK_ERR( error );
00225 
00226 #endif
00227 
00228         int nverts = _verts.size();
00229         int nedges = _edges.size();
00230         int nfaces = _faces.size();
00231         int ncells = _cells.size();
00232 
00233         MESHTYPE mesh_type = get_mesh_type( nverts, nedges, nfaces, ncells );
00234         thismeshtype       = mesh_type;
00235 
00236         // Initialize mesh type specific maps
00237         if( thismeshtype == CURVE )
00238         {
00239             error = init_curve();MB_CHK_ERR( error );
00240         }
00241         else if( thismeshtype == SURFACE )
00242         {
00243             error = init_surface();MB_CHK_ERR( error );
00244         }
00245         else if( thismeshtype == SURFACE_MIXED )
00246         {
00247             error = init_curve();MB_CHK_ERR( error );
00248             error = init_surface();MB_CHK_ERR( error );
00249         }
00250         else if( thismeshtype == VOLUME )
00251         {
00252             error = init_volume();MB_CHK_ERR( error );
00253         }
00254         else if( thismeshtype == VOLUME_MIXED_1 )
00255         {
00256             error = init_curve();MB_CHK_ERR( error );
00257             error = init_volume();MB_CHK_ERR( error );
00258         }
00259         else if( thismeshtype == VOLUME_MIXED_2 )
00260         {
00261             error = init_surface();MB_CHK_ERR( error );
00262             error = init_volume();MB_CHK_ERR( error );
00263         }
00264         else if( thismeshtype == VOLUME_MIXED )
00265         {
00266             error = init_curve();MB_CHK_ERR( error );
00267             error = init_surface();MB_CHK_ERR( error );
00268             error = init_volume();MB_CHK_ERR( error );
00269         }
00270     }
00271     return MB_SUCCESS;
00272 }
00273 
00274 ErrorCode HalfFacetRep::deinitialize()
00275 {
00276     return MB_SUCCESS;
00277 }
00278 
00279 ErrorCode HalfFacetRep::init_curve()
00280 {
00281     ErrorCode error;
00282 
00283     int nv = ID_FROM_HANDLE( *( _verts.end() - 1 ) );
00284     int ne = ID_FROM_HANDLE( *( _edges.end() - 1 ) );
00285 
00286     v2hv.resize( nv, 0 );
00287     sibhvs.resize( ne * 2, 0 );
00288 
00289     error = determine_sibling_halfverts( _verts, _edges );MB_CHK_ERR( error );
00290     error = determine_incident_halfverts( _edges );MB_CHK_ERR( error );
00291 
00292     return MB_SUCCESS;
00293 }
00294 
00295 ErrorCode HalfFacetRep::init_surface()
00296 {
00297     ErrorCode error;
00298     EntityType ftype = mb->type_from_handle( *_faces.begin() );
00299     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
00300 
00301     int nv = ID_FROM_HANDLE( *( _verts.end() - 1 ) );
00302     int nf = ID_FROM_HANDLE( *( _faces.end() - 1 ) );
00303 
00304     v2he.resize( nv, 0 );
00305     sibhes.resize( nf * nepf, 0 );
00306 
00307     // Construct ahf maps
00308     error = determine_sibling_halfedges( _faces );MB_CHK_ERR( error );
00309     error = determine_incident_halfedges( _faces );MB_CHK_ERR( error );
00310 
00311     // Initialize queues for storing face and local id's during local search
00312     for( int i = 0; i < MAXSIZE; i++ )
00313     {
00314         queue_fid[i]  = 0;
00315         queue_lid[i]  = 0;
00316         trackfaces[i] = 0;
00317     }
00318 
00319     return MB_SUCCESS;
00320 }
00321 
00322 ErrorCode HalfFacetRep::init_volume()
00323 {
00324     ErrorCode error;
00325 
00326     // Initialize std::map between cell-types and their index in lConnMap3D
00327     cell_index[MBTET]     = 0;
00328     cell_index[MBPYRAMID] = 1;
00329     cell_index[MBPRISM]   = 2;
00330     cell_index[MBHEX]     = 3;
00331 
00332     int index = get_index_in_lmap( *_cells.begin() );
00333     int nfpc  = lConnMap3D[index].num_faces_in_cell;
00334     int nv    = ID_FROM_HANDLE( *( _verts.end() - 1 ) );
00335     int nc    = ID_FROM_HANDLE( *( _cells.end() - 1 ) );
00336     ;
00337 
00338     v2hf.resize( nv, 0 );
00339     sibhfs.resize( nc * nfpc, 0 );
00340 
00341     // Construct the maps
00342     error = determine_sibling_halffaces( _cells );MB_CHK_ERR( error );
00343     error = determine_incident_halffaces( _cells );MB_CHK_ERR( error );
00344 
00345     // Initialize queues for storing face and local id's during local search
00346     for( int i = 0; i < MAXSIZE; i++ )
00347     {
00348         Stkcells[i]   = 0;
00349         cellq[i]      = 0;
00350         trackcells[i] = 0;
00351     }
00352 
00353     return MB_SUCCESS;
00354 }
00355 
00356 //////////////////////////////////////////////////
00357 ErrorCode HalfFacetRep::print_tags( int dim )
00358 {
00359     if( dim == 1 )
00360     {
00361         EntityHandle start_edge = *_edges.begin();
00362         std::cout << "start_edge = " << start_edge << std::endl;
00363         std::cout << "<SIBHVS_EID,SIBHVS_LVID>" << std::endl;
00364 
00365         for( Range::iterator i = _edges.begin(); i != _edges.end(); ++i )
00366         {
00367             EntityHandle eid[2];
00368             int lvid[2];
00369             int eidx   = ID_FROM_HANDLE( *i ) - 1;
00370             HFacet hf1 = sibhvs[2 * eidx];
00371             HFacet hf2 = sibhvs[2 * eidx + 1];
00372             eid[0]     = fid_from_halfacet( hf1, MBEDGE );
00373             eid[1]     = fid_from_halfacet( hf2, MBEDGE );
00374             lvid[0]    = lid_from_halffacet( hf1 );
00375             lvid[1]    = lid_from_halffacet( hf2 );
00376             std::cout << "Entity = " << *i << " :: <" << eid[0] << "," << lvid[0] << ">"
00377                       << "      "
00378                       << "<" << eid[1] << "," << lvid[1] << ">" << std::endl;
00379         }
00380 
00381         std::cout << "<V2HV_EID, V2HV_LVID>" << std::endl;
00382 
00383         for( Range::iterator i = _verts.begin(); i != _verts.end(); ++i )
00384         {
00385             int vidx         = ID_FROM_HANDLE( *i ) - 1;
00386             HFacet hf        = v2hv[vidx];
00387             EntityHandle eid = fid_from_halfacet( hf, MBEDGE );
00388             int lvid         = lid_from_halffacet( hf );
00389             std::cout << "Vertex = " << *i << " :: <" << eid << "," << lvid << ">" << std::endl;
00390         }
00391     }
00392     else if( dim == 2 )
00393     {
00394         EntityType ftype        = mb->type_from_handle( *_faces.begin() );
00395         int nepf                = lConnMap2D[ftype - 2].num_verts_in_face;
00396         EntityHandle start_face = *_faces.begin();
00397         std::cout << "start_face = " << start_face << std::endl;
00398         std::cout << "<SIBHES_FID,SIBHES_LEID>" << std::endl;
00399 
00400         for( Range::iterator i = _faces.begin(); i != _faces.end(); ++i )
00401         {
00402             int fidx = ID_FROM_HANDLE( *i ) - 1;
00403             std::cout << "Entity = " << *i;
00404             for( int j = 0; j < nepf; j++ )
00405             {
00406                 HFacet hf        = sibhes[nepf * fidx + j];
00407                 EntityHandle sib = fid_from_halfacet( hf, ftype );
00408                 int lid          = lid_from_halffacet( hf );
00409                 std::cout << " :: <" << sib << "," << lid << ">"
00410                           << "       ";
00411             }
00412             std::cout << std::endl;
00413         }
00414 
00415         std::cout << "<V2HE_FID, V2HE_LEID>" << std::endl;
00416 
00417         for( Range::iterator i = _verts.begin(); i != _verts.end(); ++i )
00418         {
00419             int vidx         = ID_FROM_HANDLE( *i ) - 1;
00420             HFacet hf        = v2he[vidx];
00421             EntityHandle fid = fid_from_halfacet( hf, ftype );
00422             int lid          = lid_from_halffacet( hf );
00423             std::cout << "Vertex = " << *i << " :: <" << fid << "," << lid << ">" << std::endl;
00424         }
00425     }
00426     else if( dim == 3 )
00427     {
00428         EntityType ctype = mb->type_from_handle( *_cells.begin() );
00429 
00430         int index               = get_index_in_lmap( *_cells.begin() );
00431         int nfpc                = lConnMap3D[index].num_faces_in_cell;
00432         EntityHandle start_cell = *_cells.begin();
00433         std::cout << "start_cell = " << start_cell << std::endl;
00434         std::cout << "<SIBHES_CID,SIBHES_LFID>" << std::endl;
00435 
00436         for( Range::iterator i = _cells.begin(); i != _cells.end(); ++i )
00437         {
00438             int cidx = ID_FROM_HANDLE( *i ) - 1;
00439             std::cout << "Entity = " << *i;
00440             for( int j = 0; j < nfpc; j++ )
00441             {
00442                 HFacet hf        = sibhfs[nfpc * cidx + j];
00443                 EntityHandle sib = fid_from_halfacet( hf, ctype );
00444                 int lid          = lid_from_halffacet( hf );
00445                 std::cout << " :: <" << sib << "," << lid << ">"
00446                           << "       ";
00447             }
00448             std::cout << std::endl;
00449         }
00450 
00451         std::cout << "<V2HF_CID, V2HF_LFID>" << std::endl;
00452         EntityHandle cid;
00453         int lid;
00454 
00455         for( Range::iterator i = _verts.begin(); i != _verts.end(); ++i )
00456         {
00457             int vidx  = ID_FROM_HANDLE( *i ) - 1;
00458             HFacet hf = v2hf[vidx];
00459 
00460             if( hf == 0 && ( v2hfs.find( *i ) != v2hfs.end() ) )
00461             {
00462                 std::pair< std::multimap< EntityHandle, HFacet >::iterator,
00463                            std::multimap< EntityHandle, HFacet >::iterator >
00464                     it_hfs;
00465                 it_hfs = v2hfs.equal_range( *i );
00466 
00467                 for( std::multimap< EntityHandle, HFacet >::iterator it = it_hfs.first; it != it_hfs.second; ++it )
00468                 {
00469                     cid = fid_from_halfacet( it->second, ctype );
00470                     lid = lid_from_halffacet( hf );
00471 
00472                     std::cout << "Vertex = " << *i << " :: <" << cid << "," << lid << ">" << std::endl;
00473                 }
00474             }
00475             else
00476             {
00477                 cid = fid_from_halfacet( hf, ctype );
00478                 lid = lid_from_halffacet( hf );
00479                 std::cout << "Vertex = " << *i << " :: <" << cid << "," << lid << ">" << std::endl;
00480             }
00481         }
00482     }
00483     return MB_SUCCESS;
00484 }
00485 
00486 /**********************************************************
00487  *      User interface for adjacency functions            *
00488  ********************************************************/
00489 
00490 ErrorCode HalfFacetRep::get_adjacencies( const EntityHandle source_entity, const unsigned int target_dimension,
00491                                          std::vector< EntityHandle >& target_entities )
00492 {
00493 
00494     ErrorCode error;
00495 
00496     unsigned int source_dimension = mb->dimension_from_handle( source_entity );
00497     assert( ( source_dimension <= target_dimension ) || ( source_dimension > target_dimension ) );
00498 
00499     if( mInitAHFmaps == false )
00500     {
00501         error = initialize();MB_CHK_ERR( error );
00502     }
00503 
00504     int mindex       = get_index_for_meshtype( thismeshtype );
00505     int adj_possible = adjMatrix[mindex].val[source_dimension][target_dimension];
00506 
00507     if( adj_possible )
00508     {
00509         if( source_dimension < target_dimension )
00510         {
00511             error = get_up_adjacencies( source_entity, target_dimension, target_entities );MB_CHK_ERR( error );
00512         }
00513         else if( source_dimension == target_dimension )
00514         {
00515             error = get_neighbor_adjacencies( source_entity, target_entities );MB_CHK_ERR( error );
00516         }
00517         else
00518         {
00519             error = get_down_adjacencies( source_entity, target_dimension, target_entities );MB_CHK_ERR( error );
00520         }
00521     }
00522     else
00523         return MB_SUCCESS;
00524 
00525     return MB_SUCCESS;
00526 }
00527 
00528 ErrorCode HalfFacetRep::get_up_adjacencies( EntityHandle ent, int out_dim, std::vector< EntityHandle >& adjents,
00529                                             std::vector< int >* lids )
00530 {
00531     ErrorCode error;
00532     int in_dim = mb->dimension_from_handle( ent );
00533     assert( ( in_dim >= 0 && in_dim <= 2 ) && ( out_dim > in_dim ) );
00534 
00535     if( in_dim == 0 )
00536     {
00537         if( out_dim == 1 )
00538         {
00539             error = get_up_adjacencies_1d( ent, adjents, lids );MB_CHK_ERR( error );
00540         }
00541         else if( out_dim == 2 )
00542         {
00543             error = get_up_adjacencies_vert_2d( ent, adjents );MB_CHK_ERR( error );
00544         }
00545         else if( out_dim == 3 )
00546         {
00547             error = get_up_adjacencies_vert_3d( ent, adjents );MB_CHK_ERR( error );
00548         }
00549     }
00550 
00551     else if( ( in_dim == 1 ) && ( out_dim == 2 ) )
00552     {
00553         error = get_up_adjacencies_2d( ent, adjents, lids );MB_CHK_ERR( error );
00554     }
00555     else if( ( in_dim == 1 ) && ( out_dim == 3 ) )
00556     {
00557         error = get_up_adjacencies_edg_3d( ent, adjents, lids );MB_CHK_ERR( error );
00558     }
00559     else if( ( in_dim == 2 ) && ( out_dim == 3 ) )
00560     {
00561         error = get_up_adjacencies_face_3d( ent, adjents, lids );MB_CHK_ERR( error );
00562     }
00563     return MB_SUCCESS;
00564 }
00565 
00566 ErrorCode HalfFacetRep::get_neighbor_adjacencies( EntityHandle ent, std::vector< EntityHandle >& adjents )
00567 {
00568     ErrorCode error;
00569     int in_dim = mb->dimension_from_handle( ent );
00570     assert( in_dim >= 1 && in_dim <= 3 );
00571 
00572     if( in_dim == 1 )
00573     {
00574         error = get_neighbor_adjacencies_1d( ent, adjents );MB_CHK_ERR( error );
00575     }
00576 
00577     else if( in_dim == 2 )
00578     {
00579         error = get_neighbor_adjacencies_2d( ent, adjents );MB_CHK_ERR( error );
00580     }
00581     else if( in_dim == 3 )
00582     {
00583         error = get_neighbor_adjacencies_3d( ent, adjents );MB_CHK_ERR( error );
00584     }
00585     return MB_SUCCESS;
00586 }
00587 
00588 ErrorCode HalfFacetRep::get_down_adjacencies( EntityHandle ent, int out_dim, std::vector< EntityHandle >& adjents )
00589 {
00590     ErrorCode error;
00591     int in_dim = mb->dimension_from_handle( ent );
00592     assert( ( in_dim >= 2 && in_dim <= 3 ) && ( out_dim < in_dim ) );
00593 
00594     if( ( in_dim == 2 ) && ( out_dim == 1 ) )
00595     {
00596         error = get_down_adjacencies_2d( ent, adjents );MB_CHK_ERR( error );
00597     }
00598     else if( ( in_dim == 3 ) && ( out_dim == 1 ) )
00599     {
00600         error = get_down_adjacencies_edg_3d( ent, adjents );MB_CHK_ERR( error );
00601     }
00602     else if( ( in_dim == 3 ) && ( out_dim == 2 ) )
00603     {
00604         error = get_down_adjacencies_face_3d( ent, adjents );MB_CHK_ERR( error );
00605     }
00606     return MB_SUCCESS;
00607 }
00608 
00609 ErrorCode HalfFacetRep::count_subentities( Range& edges, Range& faces, Range& cells, int* nedges, int* nfaces )
00610 {
00611     ErrorCode error;
00612     if( edges.size() && !faces.size() && !cells.size() )
00613     {
00614         nedges[0] = edges.size();
00615         nfaces[0] = 0;
00616     }
00617     else if( faces.size() && !cells.size() )
00618     {
00619         nedges[0] = find_total_edges_2d( faces );
00620         nfaces[0] = 0;
00621     }
00622     else if( cells.size() )
00623     {
00624         error = find_total_edges_faces_3d( cells, nedges, nfaces );MB_CHK_ERR( error );
00625     }
00626     return MB_SUCCESS;
00627 }
00628 
00629 /********************************************************
00630  * 1D: sibhvs, v2hv, incident and neighborhood queries   *
00631  *********************************************************/
00632 ErrorCode HalfFacetRep::determine_sibling_halfverts( Range& verts, Range& edges )
00633 {
00634     ErrorCode error;
00635 
00636     // Step 1: Create an index list storing the starting position for each vertex
00637     int nv = verts.size();
00638     std::vector< int > is_index( nv + 1 );
00639     for( int i = 0; i < nv + 1; i++ )
00640         is_index[i] = 0;
00641 
00642     for( Range::iterator eid = edges.begin(); eid != edges.end(); ++eid )
00643     {
00644         const EntityHandle* conn;
00645         int num_conn = 0;
00646         error        = mb->get_connectivity( *eid, conn, num_conn, true );MB_CHK_ERR( error );
00647 
00648         int index = verts.index( conn[0] );
00649         is_index[index + 1] += 1;
00650         index = verts.index( conn[1] );
00651         is_index[index + 1] += 1;
00652     }
00653     is_index[0] = 0;
00654 
00655     for( int i = 0; i < nv; i++ )
00656         is_index[i + 1] = is_index[i] + is_index[i + 1];
00657 
00658     // Step 2: Define two arrays v2hv_eid, v2hv_lvid storing every half-facet on a vertex
00659     std::vector< EntityHandle > v2hv_map_eid( 2 * edges.size() );
00660     std::vector< int > v2hv_map_lvid( 2 * edges.size() );
00661 
00662     for( Range::iterator eid = edges.begin(); eid != edges.end(); ++eid )
00663     {
00664         const EntityHandle* conn;
00665         int num_conn = 0;
00666         error        = mb->get_connectivity( *eid, conn, num_conn, true );MB_CHK_ERR( error );
00667 
00668         for( int j = 0; j < 2; j++ )
00669         {
00670             int v                      = verts.index( conn[j] );
00671             v2hv_map_eid[is_index[v]]  = *eid;
00672             v2hv_map_lvid[is_index[v]] = j;
00673             is_index[v] += 1;
00674         }
00675     }
00676 
00677     for( int i = nv - 2; i >= 0; i-- )
00678         is_index[i + 1] = is_index[i];
00679     is_index[0] = 0;
00680 
00681     // Step 3: Fill up sibling half-verts map
00682     for( Range::iterator vid = verts.begin(); vid != verts.end(); ++vid )
00683     {
00684         int v    = verts.index( *vid );
00685         int last = is_index[v + 1] - 1;
00686         if( last > is_index[v] )
00687         {
00688             EntityHandle prev_eid = v2hv_map_eid[last];
00689             int prev_lvid         = v2hv_map_lvid[last];
00690 
00691             for( int i = is_index[v]; i <= last; i++ )
00692             {
00693                 EntityHandle cur_eid = v2hv_map_eid[i];
00694                 int cur_lvid         = v2hv_map_lvid[i];
00695 
00696                 int pidx                     = ID_FROM_HANDLE( prev_eid ) - 1;
00697                 sibhvs[2 * pidx + prev_lvid] = create_halffacet( cur_eid, cur_lvid );
00698 
00699                 prev_eid  = cur_eid;
00700                 prev_lvid = cur_lvid;
00701             }
00702         }
00703     }
00704 
00705     return MB_SUCCESS;
00706 }
00707 
00708 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
00709 ErrorCode HalfFacetRep::determine_incident_halfverts( Range& edges )
00710 {
00711     ErrorCode error;
00712 
00713     for( Range::iterator e_it = edges.begin(); e_it != edges.end(); ++e_it )
00714     {
00715         EntityHandle cur_eid = *e_it;
00716         const EntityHandle* conn;
00717         int num_conn = 0;
00718         error        = mb->get_connectivity( *e_it, conn, num_conn, true );MB_CHK_ERR( error );
00719 
00720         for( int i = 0; i < 2; ++i )
00721         {
00722             EntityHandle v = conn[i];
00723             int vidx       = ID_FROM_HANDLE( v ) - 1;
00724 
00725             HFacet hf        = v2hv[vidx];
00726             EntityHandle eid = fid_from_halfacet( hf, MBEDGE );
00727             if( eid == 0 ) { v2hv[vidx] = create_halffacet( cur_eid, i ); }
00728         }
00729     }
00730 
00731     return MB_SUCCESS;
00732 }
00733 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
00734 ErrorCode HalfFacetRep::get_up_adjacencies_1d( EntityHandle vid, std::vector< EntityHandle >& adjents,
00735                                                std::vector< int >* lvids )
00736 {
00737     adjents.clear();
00738     adjents.reserve( 20 );
00739 
00740     if( lvids != NULL ) lvids->reserve( 20 );
00741 
00742     int vidx  = ID_FROM_HANDLE( vid ) - 1;
00743     HFacet hf = v2hv[vidx];
00744 
00745     EntityHandle start_eid = fid_from_halfacet( hf, MBEDGE );
00746     int start_lid          = lid_from_halffacet( hf );
00747 
00748     EntityHandle eid = start_eid;
00749     int lid          = start_lid;
00750 
00751     if( eid != 0 )
00752     {
00753         adjents.push_back( eid );
00754         if( lvids != NULL ) lvids->push_back( lid );
00755 
00756         while( eid != 0 )
00757         {
00758             int eidx   = ID_FROM_HANDLE( eid ) - 1;
00759             HFacet shf = sibhvs[2 * eidx + lid];
00760             eid        = fid_from_halfacet( shf, MBEDGE );
00761             lid        = lid_from_halffacet( shf );
00762 
00763             if( ( !eid ) || ( eid == start_eid ) ) break;
00764 
00765             adjents.push_back( eid );
00766             if( lvids != NULL ) lvids->push_back( lid );
00767         }
00768     }
00769 
00770     return MB_SUCCESS;
00771 }
00772 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
00773 ErrorCode HalfFacetRep::get_neighbor_adjacencies_1d( EntityHandle eid, std::vector< EntityHandle >& adjents )
00774 {
00775     adjents.clear();
00776     adjents.reserve( 20 );
00777 
00778     EntityHandle sibhv_eid;
00779     int sibhv_lid;
00780     int eidx = ID_FROM_HANDLE( eid ) - 1;
00781 
00782     for( int lid = 0; lid < 2; ++lid )
00783     {
00784         HFacet shf = sibhvs[2 * eidx + lid];
00785         sibhv_eid  = fid_from_halfacet( shf, MBEDGE );
00786         sibhv_lid  = lid_from_halffacet( shf );
00787 
00788         if( sibhv_eid != 0 )
00789         {
00790             adjents.push_back( sibhv_eid );
00791 
00792             eidx                = ID_FROM_HANDLE( sibhv_eid ) - 1;
00793             HFacet nhf          = sibhvs[2 * eidx + sibhv_lid];
00794             EntityHandle hv_eid = fid_from_halfacet( nhf, MBEDGE );
00795             int hv_lid          = lid_from_halffacet( nhf );
00796 
00797             while( hv_eid != 0 )
00798             {
00799                 if( hv_eid != eid ) adjents.push_back( hv_eid );
00800 
00801                 eidx              = ID_FROM_HANDLE( hv_eid ) - 1;
00802                 HFacet hf         = sibhvs[2 * eidx + hv_lid];
00803                 EntityHandle edge = fid_from_halfacet( hf, MBEDGE );
00804                 if( edge == sibhv_eid ) break;
00805 
00806                 hv_eid = edge;
00807                 hv_lid = lid_from_halffacet( hf );
00808             }
00809         }
00810     }
00811 
00812     return MB_SUCCESS;
00813 }
00814 
00815 /*******************************************************
00816  * 2D: sibhes, v2he, incident and neighborhood queries  *
00817  ********************************************************/
00818 const HalfFacetRep::LocalMaps2D HalfFacetRep::lConnMap2D[2] = {
00819     // Triangle
00820     { 3, { 1, 2, 0 }, { 2, 0, 1 } },
00821     // Quad
00822     { 4, { 1, 2, 3, 0 }, { 3, 0, 1, 2 } }
00823 };
00824 
00825 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
00826 ErrorCode HalfFacetRep::determine_sibling_halfedges( Range& faces )
00827 {
00828     ErrorCode error;
00829     EntityHandle start_face = *faces.begin();
00830     EntityType ftype        = mb->type_from_handle( start_face );
00831     int nfaces              = faces.size();
00832     int nepf                = lConnMap2D[ftype - 2].num_verts_in_face;
00833 
00834     // Step 1: Create an index list storing the starting position for each vertex
00835     int nv = _verts.size();
00836     std::vector< int > is_index( nv + 1 );
00837     for( int i = 0; i < nv + 1; i++ )
00838         is_index[i] = 0;
00839 
00840     int index;
00841 
00842     for( Range::iterator fid = faces.begin(); fid != faces.end(); ++fid )
00843     {
00844         const EntityHandle* conn;
00845         error = mb->get_connectivity( *fid, conn, nepf, true );MB_CHK_ERR( error );
00846 
00847         for( int i = 0; i < nepf; i++ )
00848         {
00849             index = _verts.index( conn[i] );
00850             is_index[index + 1] += 1;
00851         }
00852     }
00853     is_index[0] = 0;
00854 
00855     for( int i = 0; i < nv; i++ )
00856         is_index[i + 1] = is_index[i] + is_index[i + 1];
00857 
00858     // Step 2: Define two arrays v2hv_eid, v2hv_lvid storing every half-facet on a vertex
00859     std::vector< EntityHandle > v2nv( nepf * nfaces );
00860     std::vector< EntityHandle > v2he_map_fid( nepf * nfaces );
00861     std::vector< int > v2he_map_leid( nepf * nfaces );
00862 
00863     for( Range::iterator fid = faces.begin(); fid != faces.end(); ++fid )
00864     {
00865         const EntityHandle* conn;
00866         error = mb->get_connectivity( *fid, conn, nepf, true );MB_CHK_ERR( error );
00867 
00868         for( int j = 0; j < nepf; j++ )
00869         {
00870             int v                      = _verts.index( conn[j] );
00871             int nidx                   = lConnMap2D[ftype - 2].next[j];
00872             v2nv[is_index[v]]          = conn[nidx];
00873             v2he_map_fid[is_index[v]]  = *fid;
00874             v2he_map_leid[is_index[v]] = j;
00875             is_index[v] += 1;
00876         }
00877     }
00878 
00879     for( int i = nv - 2; i >= 0; i-- )
00880         is_index[i + 1] = is_index[i];
00881     is_index[0] = 0;
00882 
00883     // Step 3: Fill up sibling half-verts map
00884     for( Range::iterator fid = faces.begin(); fid != faces.end(); ++fid )
00885     {
00886         const EntityHandle* conn;
00887         error = mb->get_connectivity( *fid, conn, nepf, true );MB_CHK_ERR( error );
00888 
00889         int fidx = ID_FROM_HANDLE( *fid ) - 1;
00890         for( int k = 0; k < nepf; k++ )
00891         {
00892             HFacet hf           = sibhes[nepf * fidx + k];
00893             EntityHandle sibfid = fid_from_halfacet( hf, ftype );
00894 
00895             if( sibfid != 0 ) continue;
00896 
00897             int nidx = lConnMap2D[ftype - 2].next[k];
00898             int v    = _verts.index( conn[k] );
00899             int vn   = _verts.index( conn[nidx] );
00900 
00901             EntityHandle first_fid = *fid;
00902             int first_leid         = k;
00903 
00904             EntityHandle prev_fid = *fid;
00905             int prev_leid         = k;
00906 
00907             for( index = is_index[vn]; index <= is_index[vn + 1] - 1; index++ )
00908             {
00909                 if( v2nv[index] == conn[k] )
00910                 {
00911                     EntityHandle cur_fid = v2he_map_fid[index];
00912                     int cur_leid         = v2he_map_leid[index];
00913 
00914                     int pidx                        = ID_FROM_HANDLE( prev_fid ) - 1;
00915                     sibhes[nepf * pidx + prev_leid] = create_halffacet( cur_fid, cur_leid );
00916 
00917                     prev_fid  = cur_fid;
00918                     prev_leid = cur_leid;
00919                 }
00920             }
00921 
00922             for( index = is_index[v]; index <= is_index[v + 1] - 1; index++ )
00923             {
00924                 if( ( v2nv[index] == conn[nidx] ) && ( v2he_map_fid[index] != *fid ) )
00925                 {
00926 
00927                     EntityHandle cur_fid = v2he_map_fid[index];
00928                     int cur_leid         = v2he_map_leid[index];
00929 
00930                     int pidx                        = ID_FROM_HANDLE( prev_fid ) - 1;
00931                     sibhes[nepf * pidx + prev_leid] = create_halffacet( cur_fid, cur_leid );
00932 
00933                     prev_fid  = cur_fid;
00934                     prev_leid = cur_leid;
00935                 }
00936             }
00937 
00938             if( prev_fid != first_fid )
00939             {
00940 
00941                 int pidx                        = ID_FROM_HANDLE( prev_fid ) - 1;
00942                 sibhes[nepf * pidx + prev_leid] = create_halffacet( first_fid, first_leid );
00943             }
00944         }
00945     }
00946 
00947     return MB_SUCCESS;
00948 }
00949 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
00950 ErrorCode HalfFacetRep::determine_incident_halfedges( Range& faces )
00951 {
00952     ErrorCode error;
00953     EntityType ftype = mb->type_from_handle( *faces.begin() );
00954     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
00955 
00956     std::vector< char > markEdges( nepf * faces.size(), 0 );
00957 
00958     for( Range::iterator it = faces.begin(); it != faces.end(); ++it )
00959     {
00960         EntityHandle fid = *it;
00961         const EntityHandle* conn;
00962         error = mb->get_connectivity( fid, conn, nepf, true );MB_CHK_ERR( error );
00963 
00964         for( int i = 0; i < nepf; ++i )
00965         {
00966             EntityHandle v = conn[i];
00967             int vidx       = ID_FROM_HANDLE( v ) - 1;
00968             HFacet hf      = v2he[vidx];
00969 
00970             if( hf == 0 && ( v2hes.empty() || ( v2hes.find( v ) == v2hes.end() ) ) )
00971             {
00972                 // This is the first time a half-facet is assigned to a vertex.
00973                 HFacet nwhf = 0;
00974                 error       = mark_halfedges( v, fid, i, faces, markEdges, nwhf );MB_CHK_ERR( error );
00975 
00976                 if( nwhf == 0 ) nwhf = create_halffacet( fid, i );
00977 
00978                 v2he[vidx] = nwhf;
00979             }
00980             else if( hf != 0 && !markEdges[nepf * faces.index( fid ) + i] )
00981             {
00982                 // This is the first time a non-manifold vertex is encountered. Copy the existing he
00983                 // in v2he[v] to the multimap.
00984                 v2hes.insert( std::pair< EntityHandle, HFacet >( v, hf ) );
00985                 HFacet nwhf = 0;
00986                 error       = mark_halfedges( v, fid, i, faces, markEdges, nwhf );MB_CHK_ERR( error );
00987 
00988                 if( nwhf == 0 ) nwhf = create_halffacet( fid, i );
00989 
00990                 v2hes.insert( std::pair< EntityHandle, HFacet >( v, nwhf ) );
00991                 v2he[vidx] = 0;
00992             }
00993             else if( hf == 0 && ( !v2hes.empty() ) && ( v2hes.find( v ) != v2hes.end() ) &&
00994                      !markEdges[nepf * faces.index( fid ) + i] )
00995             {
00996                 // This is check if reached if the vertex is non-manifold and has encountered a
00997                 // half-facet to a new component.
00998                 HFacet nwhf = 0;
00999                 error       = mark_halfedges( v, fid, i, faces, markEdges, nwhf );MB_CHK_ERR( error );
01000 
01001                 if( nwhf == 0 ) nwhf = create_halffacet( fid, i );
01002 
01003                 v2hes.insert( std::pair< EntityHandle, HFacet >( v, nwhf ) );
01004             }
01005         }
01006     }
01007 
01008     // error = print_tags(2);
01009 
01010     return MB_SUCCESS;
01011 }
01012 
01013 ///////////////////////////////////////////////////////////////////
01014 ErrorCode HalfFacetRep::mark_halfedges( EntityHandle vid, EntityHandle he_fid, int he_lid, Range& faces,
01015                                         std::vector< char >& markHEdgs, HFacet& bnd_hf )
01016 {
01017     ErrorCode error;
01018     EntityType ftype = mb->type_from_handle( he_fid );
01019     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
01020 
01021     int qsize = 0, count = -1;
01022     int num_qvals = 0;
01023 
01024     error = gather_halfedges( vid, he_fid, he_lid, &qsize, &count );MB_CHK_ERR( error );
01025 
01026     while( num_qvals < qsize )
01027     {
01028         EntityHandle curfid = queue_fid[num_qvals];
01029         int curlid          = queue_lid[num_qvals];
01030         num_qvals += 1;
01031 
01032         int fidx = ID_FROM_HANDLE( curfid ) - 1;
01033 
01034         const EntityHandle* conn;
01035         error = mb->get_connectivity( curfid, conn, nepf, true );MB_CHK_ERR( error );
01036 
01037         if( !markHEdgs[nepf * faces.index( curfid ) + curlid] && ( conn[curlid] == vid ) )
01038         {
01039             markHEdgs[nepf * faces.index( curfid ) + curlid] = 1;
01040             HFacet hf                                        = sibhes[nepf * fidx + curlid];
01041             EntityHandle sibfid                              = fid_from_halfacet( hf, ftype );
01042             if( sibfid == 0 ) bnd_hf = create_halffacet( curfid, curlid );
01043         }
01044 
01045         EntityHandle he2_fid = 0;
01046         int he2_lid          = 0;
01047         error                = another_halfedge( vid, curfid, curlid, &he2_fid, &he2_lid );MB_CHK_ERR( error );
01048 
01049         if( !markHEdgs[nepf * faces.index( curfid ) + he2_lid] && ( conn[he2_lid] == vid ) )
01050         {
01051             markHEdgs[nepf * faces.index( curfid ) + he2_lid] = 1;
01052             HFacet hf                                         = sibhes[nepf * fidx + he2_lid];
01053             EntityHandle sibfid                               = fid_from_halfacet( hf, ftype );
01054             if( sibfid == 0 ) bnd_hf = create_halffacet( he2_fid, he2_lid );
01055         }
01056 
01057         bool val = find_match_in_array( he2_fid, trackfaces, count );
01058 
01059         if( val ) continue;
01060 
01061         count += 1;
01062         trackfaces[count] = he2_fid;
01063 
01064         error = get_up_adjacencies_2d( he2_fid, he2_lid, &qsize, &count );MB_CHK_ERR( error );
01065     }
01066 
01067     // Change the visited faces to false, also empty the queue
01068     for( int i = 0; i <= qsize; i++ )
01069     {
01070         queue_fid[i] = 0;
01071         queue_lid[i] = 0;
01072     }
01073 
01074     for( int i = 0; i <= count; i++ )
01075         trackfaces[i] = 0;
01076     return MB_SUCCESS;
01077 }
01078 
01079 ///////////////////////////////////////////////////////////////////
01080 ErrorCode HalfFacetRep::get_up_adjacencies_vert_2d( EntityHandle vid, std::vector< EntityHandle >& adjents )
01081 {
01082     ErrorCode error;
01083     EntityType ftype = mb->type_from_handle( *_faces.begin() );
01084 
01085     int vidx  = ID_FROM_HANDLE( vid ) - 1;
01086     HFacet hf = v2he[vidx];
01087 
01088     std::vector< EntityHandle > start_fids;
01089     std::vector< int > start_lids;
01090 
01091     if( hf == 0 && ( v2hes.find( vid ) != v2hes.end() ) )
01092     {
01093         std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
01094             it_hes;
01095         it_hes = v2hes.equal_range( vid );
01096 
01097         for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
01098         {
01099             start_fids.push_back( fid_from_halfacet( it->second, ftype ) );
01100             start_lids.push_back( lid_from_halffacet( it->second ) );
01101         }
01102     }
01103     else if( hf != 0 )
01104     {
01105         start_fids.push_back( fid_from_halfacet( hf, ftype ) );
01106         start_lids.push_back( lid_from_halffacet( hf ) );
01107     }
01108 
01109     if( start_fids.empty() ) return MB_SUCCESS;
01110 
01111     int qsize = 0, count = -1;
01112     int num_qvals = 0;
01113 
01114     adjents.reserve( (int)start_fids.size() );
01115 
01116     for( int i = 0; i < (int)start_fids.size(); i++ )
01117     {
01118         adjents.push_back( start_fids[i] );
01119         error = gather_halfedges( vid, start_fids[i], start_lids[i], &qsize, &count );MB_CHK_ERR( error );
01120     }
01121 
01122     while( num_qvals < qsize )
01123     {
01124         EntityHandle curfid = queue_fid[num_qvals];
01125         int curlid          = queue_lid[num_qvals];
01126         num_qvals += 1;
01127 
01128         EntityHandle he2_fid = 0;
01129         int he2_lid          = 0;
01130         error                = another_halfedge( vid, curfid, curlid, &he2_fid, &he2_lid );MB_CHK_ERR( error );
01131 
01132         bool val = find_match_in_array( he2_fid, trackfaces, count );
01133 
01134         if( val ) continue;
01135 
01136         count += 1;
01137         trackfaces[count] = he2_fid;
01138 
01139         error = get_up_adjacencies_2d( he2_fid, he2_lid, &qsize, &count );MB_CHK_ERR( error );
01140 
01141         adjents.push_back( he2_fid );
01142     }
01143 
01144     // Change the visited faces to false, also empty the queue
01145     for( int i = 0; i <= qsize; i++ )
01146     {
01147         queue_fid[i] = 0;
01148         queue_lid[i] = 0;
01149     }
01150 
01151     for( int i = 0; i <= count; i++ )
01152         trackfaces[i] = 0;
01153 
01154     return MB_SUCCESS;
01155 }
01156 
01157 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
01158 ErrorCode HalfFacetRep::get_up_adjacencies_2d( EntityHandle eid, std::vector< EntityHandle >& adjents,
01159                                                std::vector< int >* leids )
01160 {
01161 
01162     // Given an explicit edge eid, find the incident faces.
01163     ErrorCode error;
01164     EntityHandle he_fid = 0;
01165     int he_lid          = 0;
01166 
01167     // Step 1: Given an explicit edge, find a corresponding half-edge from the surface mesh.
01168     bool found = find_matching_halfedge( eid, &he_fid, &he_lid );
01169 
01170     // Step 2: If there is a corresponding half-edge, collect all sibling half-edges and store the
01171     // incident faces.
01172     if( found )
01173     {
01174         error = get_up_adjacencies_2d( he_fid, he_lid, true, adjents, leids );MB_CHK_ERR( error );
01175     }
01176 
01177     return MB_SUCCESS;
01178 }
01179 
01180 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
01181 ErrorCode HalfFacetRep::get_up_adjacencies_2d( EntityHandle fid, int leid, bool add_inent,
01182                                                std::vector< EntityHandle >& adj_ents, std::vector< int >* adj_leids,
01183                                                std::vector< int >* adj_orients )
01184 {
01185     // Given an implicit half-edge <fid, leid>, find the incident half-edges.
01186     ErrorCode error;
01187 
01188     EntityType ftype = mb->type_from_handle( fid );
01189     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
01190 
01191     if( !fid ) return MB_FAILURE;
01192     adj_ents.reserve( 20 );
01193 
01194     bool local_id = false;
01195     bool orient   = false;
01196     if( adj_leids != NULL )
01197     {
01198         local_id = true;
01199         adj_leids->reserve( 20 );
01200     }
01201     if( adj_orients != NULL )
01202     {
01203         orient = true;
01204         adj_orients->reserve( 20 );
01205     }
01206 
01207     if( add_inent )
01208     {
01209         adj_ents.push_back( fid );
01210         if( local_id ) adj_leids->push_back( leid );
01211     }
01212 
01213     EntityHandle fedge[2] = { 0, 0 };
01214 
01215     if( orient )
01216     {
01217         // get connectivity and match their directions
01218         const EntityHandle* fid_conn;
01219         error = mb->get_connectivity( fid, fid_conn, nepf, true );MB_CHK_ERR( error );
01220 
01221         int nidx = lConnMap2D[ftype - 2].next[leid];
01222         fedge[0] = fid_conn[leid];
01223         fedge[1] = fid_conn[nidx];
01224     }
01225 
01226     int fidx            = ID_FROM_HANDLE( fid ) - 1;
01227     HFacet hf           = sibhes[nepf * fidx + leid];
01228     EntityHandle curfid = fid_from_halfacet( hf, ftype );
01229     int curlid          = lid_from_halffacet( hf );
01230 
01231     while( ( curfid != fid ) && ( curfid != 0 ) )
01232     {  // Should not go into the loop when no sibling exists
01233         adj_ents.push_back( curfid );
01234 
01235         if( local_id ) adj_leids->push_back( curlid );
01236 
01237         if( orient )
01238         {
01239             // get connectivity and match their directions
01240             const EntityHandle* conn;
01241             error = mb->get_connectivity( curfid, conn, nepf, true );MB_CHK_ERR( error );
01242 
01243             int nidx = lConnMap2D[ftype - 2].next[curlid];
01244 
01245             if( ( fedge[0] == conn[curlid] ) && ( fedge[1] == conn[nidx] ) )
01246                 adj_orients->push_back( 1 );
01247             else if( ( fedge[1] == conn[curlid] ) && ( fedge[0] == conn[nidx] ) )
01248                 adj_orients->push_back( 0 );
01249         }
01250 
01251         int cidx = ID_FROM_HANDLE( curfid ) - 1;
01252         hf       = sibhes[nepf * cidx + curlid];
01253         curfid   = fid_from_halfacet( hf, ftype );
01254         curlid   = lid_from_halffacet( hf );
01255     }
01256 
01257     return MB_SUCCESS;
01258 }
01259 
01260 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
01261 ErrorCode HalfFacetRep::get_up_adjacencies_2d( EntityHandle fid, int lid, int* qsize, int* count )
01262 {
01263     EntityType ftype = mb->type_from_handle( fid );
01264     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
01265 
01266     int fidx            = ID_FROM_HANDLE( fid ) - 1;
01267     HFacet hf           = sibhes[nepf * fidx + lid];
01268     EntityHandle curfid = fid_from_halfacet( hf, ftype );
01269     int curlid          = lid_from_halffacet( hf );
01270 
01271     if( curfid == 0 )
01272     {
01273         int index      = 0;
01274         bool found_ent = find_match_in_array( fid, queue_fid, qsize[0] - 1, true, &index );
01275 
01276         if( ( !found_ent ) || ( ( found_ent ) && ( queue_lid[index] != lid ) ) )
01277         {
01278             queue_fid[qsize[0]] = fid;
01279             queue_lid[qsize[0]] = lid;
01280             qsize[0] += 1;
01281         }
01282     }
01283 
01284     while( ( curfid != fid ) && ( curfid != 0 ) )
01285     {
01286         bool val = find_match_in_array( curfid, trackfaces, count[0] );
01287 
01288         if( !val )
01289         {
01290             queue_fid[qsize[0]] = curfid;
01291             queue_lid[qsize[0]] = curlid;
01292             qsize[0] += 1;
01293         }
01294 
01295         int cidx = ID_FROM_HANDLE( curfid ) - 1;
01296         hf       = sibhes[nepf * cidx + curlid];
01297         curfid   = fid_from_halfacet( hf, ftype );
01298         curlid   = lid_from_halffacet( hf );
01299     }
01300 
01301     return MB_SUCCESS;
01302 }
01303 
01304 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
01305 bool HalfFacetRep::find_matching_halfedge( EntityHandle eid, EntityHandle* hefid, int* helid )
01306 {
01307     ErrorCode error;
01308     EntityType ftype = mb->type_from_handle( *_faces.begin() );
01309 
01310     const EntityHandle* conn;
01311     int num_conn = 0;
01312     error        = mb->get_connectivity( eid, conn, num_conn, true );MB_CHK_ERR( error );
01313 
01314     EntityHandle vid = conn[0];
01315     int vidx         = ID_FROM_HANDLE( conn[0] ) - 1;
01316     HFacet hf        = v2he[vidx];
01317 
01318     if( hf == 0 )
01319     {
01320         vidx = ID_FROM_HANDLE( conn[1] ) - 1;
01321         hf   = v2he[vidx];
01322 
01323         if( hf == 0 )  // The edge is either a dangling edge or attached to two non-manifold
01324                        // vertices
01325             return MB_SUCCESS;
01326 
01327         vid = conn[1];
01328     }
01329 
01330     EntityHandle fid = fid_from_halfacet( hf, ftype );
01331     int lid          = lid_from_halffacet( hf );
01332 
01333     bool found = false;
01334     int qsize = 0, count = -1;
01335 
01336     error = gather_halfedges( vid, fid, lid, &qsize, &count );MB_CHK_ERR( error );
01337 
01338     found = collect_and_compare( vid, conn, &qsize, &count, hefid, helid );MB_CHK_ERR( error );
01339 
01340     // Change the visited faces to false
01341     for( int i = 0; i < qsize; i++ )
01342     {
01343         queue_fid[i] = 0;
01344         queue_lid[i] = 0;
01345     }
01346 
01347     for( int i = 0; i <= count; i++ )
01348         trackfaces[i] = 0;
01349 
01350     return found;
01351 }
01352 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
01353 ErrorCode HalfFacetRep::gather_halfedges( EntityHandle vid, EntityHandle he_fid, int he_lid, int* qsize, int* count )
01354 {
01355     ErrorCode error;
01356     EntityHandle he2_fid = 0;
01357     int he2_lid          = 0;
01358 
01359     error = another_halfedge( vid, he_fid, he_lid, &he2_fid, &he2_lid );MB_CHK_ERR( error );
01360 
01361     queue_fid[*qsize] = he_fid;
01362     queue_lid[*qsize] = he_lid;
01363     *qsize += 1;
01364 
01365     queue_fid[*qsize] = he2_fid;
01366     queue_lid[*qsize] = he2_lid;
01367     *qsize += 1;
01368 
01369     *count += 1;
01370     trackfaces[*count] = he_fid;
01371 
01372     error = get_up_adjacencies_2d( he_fid, he_lid, qsize, count );MB_CHK_ERR( error );
01373     error = get_up_adjacencies_2d( he2_fid, he2_lid, qsize, count );MB_CHK_ERR( error );
01374 
01375     return MB_SUCCESS;
01376 }
01377 
01378 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
01379 ErrorCode HalfFacetRep::another_halfedge( EntityHandle vid, EntityHandle he_fid, int he_lid, EntityHandle* he2_fid,
01380                                           int* he2_lid )
01381 {
01382     ErrorCode error;
01383     EntityType ftype = mb->type_from_handle( he_fid );
01384     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
01385 
01386     const EntityHandle* conn;
01387     error = mb->get_connectivity( he_fid, conn, nepf, true );MB_CHK_ERR( error );
01388 
01389     *he2_fid = he_fid;
01390     if( conn[he_lid] == vid )
01391         *he2_lid = lConnMap2D[ftype - 2].prev[he_lid];
01392     else
01393         *he2_lid = lConnMap2D[ftype - 2].next[he_lid];
01394 
01395     return MB_SUCCESS;
01396 }
01397 
01398 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
01399 bool HalfFacetRep::collect_and_compare( const EntityHandle vid, const EntityHandle* edg_vert, int* qsize, int* count,
01400                                         EntityHandle* he_fid, int* he_lid )
01401 {
01402     ErrorCode error;
01403     EntityType ftype = mb->type_from_handle( *_faces.begin() );
01404     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
01405 
01406     bool found    = false;
01407     int num_qvals = 0, counter = 0;
01408 
01409     while( num_qvals < *qsize && counter < MAXSIZE )
01410     {
01411         EntityHandle curfid = queue_fid[num_qvals];
01412         int curlid          = queue_lid[num_qvals];
01413         num_qvals += 1;
01414 
01415         const EntityHandle* conn;
01416         error = mb->get_connectivity( curfid, conn, nepf, true );MB_CHK_ERR( error );
01417 
01418         int id = lConnMap2D[ftype - 2].next[curlid];
01419         if( ( ( conn[curlid] == edg_vert[0] ) && ( conn[id] == edg_vert[1] ) ) ||
01420             ( ( conn[curlid] == edg_vert[1] ) && ( conn[id] == edg_vert[0] ) ) )
01421         {
01422             *he_fid = curfid;
01423             *he_lid = curlid;
01424             found   = true;
01425             break;
01426         }
01427 
01428         bool val = find_match_in_array( curfid, trackfaces, count[0] );
01429 
01430         if( val ) continue;
01431 
01432         count[0] += 1;
01433         trackfaces[*count] = curfid;
01434 
01435         EntityHandle he2_fid;
01436         int he2_lid;
01437         error = another_halfedge( vid, curfid, curlid, &he2_fid, &he2_lid );MB_CHK_ERR( error );
01438         error = get_up_adjacencies_2d( he2_fid, he2_lid, qsize, count );MB_CHK_ERR( error );
01439 
01440         counter += 1;
01441     }
01442     return found;
01443 }
01444 
01445 ///////////////////////////////////////////////////////////////////////////////////////////////////
01446 ErrorCode HalfFacetRep::get_neighbor_adjacencies_2d( EntityHandle fid, std::vector< EntityHandle >& adjents )
01447 {
01448     ErrorCode error;
01449 
01450     if( fid != 0 )
01451     {
01452         EntityType ftype = mb->type_from_handle( fid );
01453         int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
01454 
01455         for( int lid = 0; lid < nepf; ++lid )
01456         {
01457             error = get_up_adjacencies_2d( fid, lid, false, adjents );MB_CHK_ERR( error );
01458         }
01459     }
01460 
01461     return MB_SUCCESS;
01462 }
01463 
01464 /////////////////////////////////////////////////////////////////////////////////////////////////
01465 ErrorCode HalfFacetRep::get_down_adjacencies_2d( EntityHandle fid, std::vector< EntityHandle >& adjents )
01466 {
01467     // Returns explicit edges, if any, of the face
01468     ErrorCode error;
01469     adjents.reserve( 10 );
01470     EntityType ftype = mb->type_from_handle( fid );
01471     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
01472 
01473     const EntityHandle* conn;
01474     error = mb->get_connectivity( fid, conn, nepf, true );MB_CHK_ERR( error );
01475 
01476     std::vector< EntityHandle > temp;
01477 
01478     // Loop over 2 vertices
01479     for( int i = 0; i < 2; i++ )
01480     {
01481         // Choose the two adjacent vertices for a triangle, and two opposite vertices for a quad
01482         int l;
01483         if( ftype == MBTRI )
01484             l = i;
01485         else
01486             l = 2 * i;
01487 
01488         // Get the current, next and prev vertices
01489         int nidx           = lConnMap2D[ftype - 2].next[l];
01490         int pidx           = lConnMap2D[ftype - 2].prev[l];
01491         EntityHandle v     = conn[l];
01492         EntityHandle vnext = conn[nidx];
01493         EntityHandle vprev = conn[pidx];
01494 
01495         // Get incident edges on v
01496         error = get_up_adjacencies_1d( v, temp );MB_CHK_ERR( error );
01497 
01498         // Loop over the incident edges and check if its end vertices match those in the face
01499         for( int k = 0; k < (int)temp.size(); k++ )
01500         {
01501             const EntityHandle* econn;
01502             int num_conn = 0;
01503             error        = mb->get_connectivity( temp[k], econn, num_conn, true );MB_CHK_ERR( error );
01504 
01505             if( ( econn[0] == v && econn[1] == vnext ) || ( econn[0] == v && econn[1] == vprev ) ||
01506                 ( econn[0] == vnext && econn[1] == v ) || ( econn[0] == vprev && econn[1] == v ) )
01507             {
01508                 bool found = false;
01509                 for( int j = 0; j < (int)adjents.size(); j++ )
01510                 {
01511                     if( adjents[j] == temp[k] )
01512                     {
01513                         found = true;
01514                         break;
01515                     }
01516                 }
01517                 if( !found ) adjents.push_back( temp[k] );
01518             }
01519         }
01520     }
01521 
01522     return MB_SUCCESS;
01523 }
01524 
01525 ////////////////////////////////////////////////////////////////////////////////////////////////
01526 int HalfFacetRep::find_total_edges_2d( Range& faces )
01527 {
01528     ErrorCode error;
01529     EntityType ftype = mb->type_from_handle( *faces.begin() );
01530     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
01531     int nfaces       = faces.size();
01532 
01533     int total_edges = nepf * nfaces;
01534 
01535     std::vector< int > trackF( total_edges, 0 );
01536     std::vector< EntityHandle > adj_fids;
01537     std::vector< int > adj_lids;
01538 
01539     for( Range::iterator f = faces.begin(); f != faces.end(); ++f )
01540     {
01541         for( int l = 0; l < nepf; l++ )
01542         {
01543 
01544             adj_fids.clear();
01545             adj_lids.clear();
01546 
01547             int id = nepf * ( faces.index( *f ) ) + l;
01548             if( !trackF[id] )
01549             {
01550                 error = get_up_adjacencies_2d( *f, l, false, adj_fids, &adj_lids );MB_CHK_ERR( error );
01551 
01552                 total_edges -= adj_fids.size();
01553 
01554                 for( int i = 0; i < (int)adj_fids.size(); i++ )
01555                     trackF[nepf * ( faces.index( adj_fids[i] ) ) + adj_lids[i]] = 1;
01556             };
01557         };
01558     };
01559 
01560     return total_edges;
01561 }
01562 
01563 ErrorCode HalfFacetRep::get_face_edges( EntityHandle fid, std::vector< EntityHandle >& edges )
01564 {
01565     ErrorCode error;
01566     edges.clear();
01567 
01568     EntityType ftype = mb->type_from_handle( fid );
01569     int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
01570 
01571     std::vector< EntityHandle > conn;
01572     error = mb->get_connectivity( &fid, 1, conn );MB_CHK_ERR( error );
01573 
01574     for( int i = 0; i < nepf; i++ )
01575     {
01576         EntityHandle v0 = conn[i];
01577         EntityHandle v1 = conn[lConnMap2D[ftype - 2].next[i]];
01578 
01579         std::vector< EntityHandle > e0, e1, ecom;
01580         error = get_up_adjacencies_1d( v0, e0 );MB_CHK_ERR( error );
01581         error = get_up_adjacencies_1d( v1, e1 );MB_CHK_ERR( error );
01582 
01583         std::sort( e0.begin(), e0.end() );
01584         std::sort( e1.begin(), e1.end() );
01585         std::set_intersection( e0.begin(), e0.end(), e1.begin(), e1.end(), std::back_inserter( ecom ) );
01586         assert( ecom.size() == 1 || ecom.size() == 0 );
01587         if( ecom.size() == 0 )
01588             edges.push_back( 0 );
01589         else
01590             edges.push_back( ecom[0] );
01591     }
01592 
01593     return MB_SUCCESS;
01594 }
01595 
01596 /*******************************************************
01597  * 3D: sibhfs, v2hf, incident and neighborhood queries  *
01598  ********************************************************/
01599 
01600 int HalfFacetRep::get_index_in_lmap( EntityHandle cid )
01601 {
01602     EntityType type = mb->type_from_handle( cid );
01603     int index       = cell_index.find( type )->second;
01604     return index;
01605 }
01606 
01607 const HalfFacetRep::LocalMaps3D HalfFacetRep::lConnMap3D[4] = {
01608     // Tet
01609     { 4,
01610       6,
01611       4,
01612       { 3, 3, 3, 3 },
01613       { { 0, 1, 3 }, { 1, 2, 3 }, { 2, 0, 3 }, { 0, 2, 1 } },
01614       { 3, 3, 3, 3 },
01615       { { 0, 2, 3 }, { 0, 1, 3 }, { 1, 2, 3 }, { 0, 1, 2 } },
01616       { { 0, 1 }, { 1, 2 }, { 2, 0 }, { 0, 3 }, { 1, 3 }, { 2, 3 } },
01617       { { 3, 0 }, { 3, 1 }, { 3, 2 }, { 0, 2 }, { 0, 1 }, { 1, 2 } },
01618       { { 0, 4, 3 }, { 1, 5, 4 }, { 2, 3, 5 }, { 2, 1, 0 } },
01619       { { -1, 0, 2, 3 }, { 0, -1, 1, 4 }, { 2, 1, -1, 5 }, { 3, 4, 5, -1 } },
01620       { 3, 0, 1, 2 },
01621       { 0, 1 },
01622       { { 3, 1, 2, 3 }, { 2, 2, 3 }, { 1, 3 } } },
01623 
01624     // Pyramid: Note: In MOAB pyramid follows the CGNS convention. Look up src/MBCNArrays.hpp
01625     { 5,
01626       8,
01627       5,
01628       { 4, 3, 3, 3, 3 },
01629       { { 0, 3, 2, 1 }, { 0, 1, 4 }, { 1, 2, 4 }, { 2, 3, 4 }, { 3, 0, 4 } },
01630       { 3, 3, 3, 3, 4 },
01631       { { 0, 1, 4 }, { 0, 1, 2 }, { 0, 2, 3 }, { 0, 3, 4 }, { 1, 2, 3, 4 } },
01632       { { 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 0 }, { 0, 4 }, { 1, 4 }, { 2, 4 }, { 3, 4 } },
01633       { { 0, 1 }, { 0, 2 }, { 0, 3 }, { 0, 4 }, { 1, 4 }, { 1, 2 }, { 2, 3 }, { 3, 4 } },
01634       { { 3, 2, 1, 0 }, { 0, 5, 4 }, { 1, 6, 5 }, { 2, 7, 6 }, { 3, 4, 7 } },
01635       { { -1, 0, -1, 3, 4 }, { 0, -1, 1, -1, 5 }, { -1, 1, -1, 2, 6 }, { 3, -1, 2, -1, 7 }, { 4, 5, 6, 7, -1 } },
01636       { 3, 4, 2, 0 },
01637       { 0, 4 },
01638       { { 4, 0, 1, 2, 3 }, { 2, 1, 3 }, { 2, 1, 3 } } },
01639 
01640     // Prism
01641     { 6,
01642       9,
01643       5,
01644       { 4, 4, 4, 3, 3 },
01645       { { 0, 1, 4, 3 }, { 1, 2, 5, 4 }, { 0, 3, 5, 2 }, { 0, 2, 1 }, { 3, 4, 5 } },
01646       { 3, 3, 3, 3, 3, 3 },
01647       { { 0, 2, 3 }, { 0, 1, 3 }, { 1, 2, 3 }, { 0, 2, 4 }, { 0, 1, 4 }, { 1, 4, 2 } },
01648       { { 0, 1 }, { 1, 2 }, { 2, 0 }, { 0, 3 }, { 1, 4 }, { 2, 5 }, { 3, 4 }, { 4, 5 }, { 5, 3 } },
01649       { { 0, 3 }, { 1, 3 }, { 2, 3 }, { 0, 2 }, { 0, 1 }, { 1, 2 }, { 0, 4 }, { 1, 4 }, { 2, 4 } },
01650       { { 0, 4, 6, 3 }, { 1, 5, 7, 4 }, { 2, 3, 8, 5 }, { 2, 1, 0 }, { 6, 7, 8 } },
01651       { { -1, 0, 2, 3, -1, -1 },
01652         { 0, -1, 1, -1, 4, -1 },
01653         { 2, 1, -1, -1, -1, 5 },
01654         { 3, -1, -1, -1, 6, 8 },
01655         { -1, 4, -1, 6, -1, 7 },
01656         { -1, -1, 5, 8, 7, -1 } },
01657       { 4, 0, 5, 4, 1 },
01658       { 0, 5 },
01659       { { 3, 1, 2, 3 }, { 3, 2, 4, 3 }, { 2, 3, 1 }, { 1, 2 } } },
01660 
01661     // Hex
01662     { 8,
01663       12,
01664       6,
01665       { 4, 4, 4, 4, 4, 4 },
01666       { { 0, 1, 5, 4 }, { 1, 2, 6, 5 }, { 2, 3, 7, 6 }, { 3, 0, 4, 7 }, { 0, 3, 2, 1 }, { 4, 5, 6, 7 } },
01667       { 3, 3, 3, 3, 3, 3, 3, 3 },
01668       { { 0, 3, 4 }, { 0, 1, 4 }, { 1, 2, 4 }, { 2, 3, 4 }, { 0, 3, 5 }, { 0, 1, 5 }, { 1, 2, 5 }, { 2, 3, 5 } },
01669       { { 0, 1 },
01670         { 1, 2 },
01671         { 2, 3 },
01672         { 3, 0 },
01673         { 0, 4 },
01674         { 1, 5 },
01675         { 2, 6 },
01676         { 3, 7 },
01677         { 4, 5 },
01678         { 5, 6 },
01679         { 6, 7 },
01680         { 7, 4 } },
01681       { { 0, 4 },
01682         { 1, 4 },
01683         { 2, 4 },
01684         { 3, 4 },
01685         { 0, 3 },
01686         { 0, 1 },
01687         { 1, 2 },
01688         { 2, 3 },
01689         { 0, 5 },
01690         { 1, 5 },
01691         { 2, 5 },
01692         { 3, 5 } },
01693       { { 0, 5, 8, 4 }, { 1, 6, 9, 5 }, { 2, 7, 10, 6 }, { 3, 4, 11, 7 }, { 3, 2, 1, 0 }, { 8, 9, 10, 11 } },
01694       { { -1, 0, -1, 3, 4, -1, -1, -1 },
01695         { 0, -1, 1, -1, -1, 5, -1, -1 },
01696         { -1, 1, -1, 2, -1, -1, 6, -1 },
01697         { 3, -1, 2, -1, -1, -1, -1, 7 },
01698         { 4, -1, -1, -1, -1, 8, -1, 11 },
01699         { -1, 5, -1, -1, 8, -1, 9, -1 },
01700         { -1, -1, 6, -1, -1, 9, -1, 10 },
01701         { -1, -1, -1, 7, 11, -1, 10, -1 } },
01702       { 4, 0, 2, 5, 7 },
01703       { 0, 6 },
01704       { { 3, 1, 3, 4 }, { 3, 1, 3, 6 }, { 3, 1, 4, 6 }, { 3, 3, 6, 4 } } }
01705 };
01706 
01707 //////////////////////////////////////////////////////////////////////////////////
01708 ErrorCode HalfFacetRep::determine_sibling_halffaces( Range& cells )
01709 {
01710     ErrorCode error;
01711     EntityHandle start_cell = *cells.begin();
01712     EntityType ctype        = mb->type_from_handle( start_cell );
01713     int index               = get_index_in_lmap( start_cell );
01714     int nvpc                = lConnMap3D[index].num_verts_in_cell;
01715     int nfpc                = lConnMap3D[index].num_faces_in_cell;
01716 
01717     // Step 1: Create an index list storing the starting position for each vertex
01718     int nv = _verts.size();
01719     std::vector< int > is_index( nv + 1 );
01720     for( int i = 0; i < nv + 1; i++ )
01721         is_index[i] = 0;
01722 
01723     int vindex;
01724 
01725     for( Range::iterator cid = cells.begin(); cid != cells.end(); ++cid )
01726     {
01727         const EntityHandle* conn;
01728         error = mb->get_connectivity( *cid, conn, nvpc, true );MB_CHK_ERR( error );
01729 
01730         for( int i = 0; i < nfpc; ++i )
01731         {
01732             int nvF        = lConnMap3D[index].hf2v_num[i];
01733             EntityHandle v = 0;
01734             for( int k = 0; k < nvF; k++ )
01735             {
01736                 int id = lConnMap3D[index].hf2v[i][k];
01737                 if( v <= conn[id] ) v = conn[id];
01738             }
01739             vindex = _verts.index( v );
01740             is_index[vindex + 1] += 1;
01741         }
01742     }
01743     is_index[0] = 0;
01744 
01745     for( int i = 0; i < nv; i++ )
01746         is_index[i + 1] = is_index[i] + is_index[i + 1];
01747 
01748     // Step 2: Define four arrays v2hv_eid, v2hv_lvid storing every half-facet on a vertex
01749     std::vector< EntityHandle > v2oe_v1( is_index[nv] );
01750     std::vector< EntityHandle > v2oe_v2( is_index[nv] );
01751     std::vector< EntityHandle > v2hf_map_cid( is_index[nv] );
01752     std::vector< int > v2hf_map_lfid( is_index[nv] );
01753 
01754     for( Range::iterator cid = cells.begin(); cid != cells.end(); ++cid )
01755     {
01756         const EntityHandle* conn;
01757         error = mb->get_connectivity( *cid, conn, nvpc, true );MB_CHK_ERR( error );
01758 
01759         for( int i = 0; i < nfpc; i++ )
01760         {
01761             int nvF = lConnMap3D[index].hf2v_num[i];
01762             std::vector< EntityHandle > vs( nvF );
01763             EntityHandle vmax = 0;
01764             int lv            = -1;
01765             for( int k = 0; k < nvF; k++ )
01766             {
01767                 int id = lConnMap3D[index].hf2v[i][k];
01768                 vs[k]  = conn[id];
01769                 if( vmax <= conn[id] )
01770                 {
01771                     vmax = conn[id];
01772                     lv   = k;
01773                 }
01774             }
01775 
01776             int nidx = lConnMap2D[nvF - 3].next[lv];
01777             int pidx = lConnMap2D[nvF - 3].prev[lv];
01778 
01779             int v                      = _verts.index( vmax );
01780             v2oe_v1[is_index[v]]       = vs[nidx];
01781             v2oe_v2[is_index[v]]       = vs[pidx];
01782             v2hf_map_cid[is_index[v]]  = *cid;
01783             v2hf_map_lfid[is_index[v]] = i;
01784             is_index[v] += 1;
01785         }
01786     }
01787 
01788     for( int i = nv - 2; i >= 0; i-- )
01789         is_index[i + 1] = is_index[i];
01790     is_index[0] = 0;
01791 
01792     // Step 3: Fill up sibling half-verts map
01793     for( Range::iterator cid = cells.begin(); cid != cells.end(); ++cid )
01794     {
01795         const EntityHandle* conn;
01796         error = mb->get_connectivity( *cid, conn, nvpc, true );MB_CHK_ERR( error );
01797 
01798         int cidx = ID_FROM_HANDLE( *cid ) - 1;
01799         for( int i = 0; i < nfpc; i++ )
01800         {
01801             HFacet hf           = sibhfs[nfpc * cidx + i];
01802             EntityHandle sibcid = fid_from_halfacet( hf, ctype );
01803 
01804             if( sibcid != 0 ) continue;
01805 
01806             int nvF = lConnMap3D[index].hf2v_num[i];
01807             std::vector< EntityHandle > vs( nvF );
01808             EntityHandle vmax = 0;
01809             int lv            = -1;
01810             for( int k = 0; k < nvF; k++ )
01811             {
01812                 int id = lConnMap3D[index].hf2v[i][k];
01813                 vs[k]  = conn[id];
01814                 if( vmax <= conn[id] )
01815                 {
01816                     vmax = conn[id];
01817                     lv   = k;
01818                 }
01819             }
01820             if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
01821             int nidx = lConnMap2D[nvF - 3].next[lv];
01822             int pidx = lConnMap2D[nvF - 3].prev[lv];
01823 
01824             int v           = _verts.index( vmax );
01825             EntityHandle v1 = vs[pidx];
01826             EntityHandle v2 = vs[nidx];
01827 
01828             for( int ind = is_index[v]; ind <= is_index[v + 1] - 1; ind++ )
01829             {
01830                 if( ( v2oe_v1[ind] == v1 ) && ( v2oe_v2[ind] == v2 ) )
01831                 {
01832                     // Map to opposite hf
01833                     EntityHandle cur_cid = v2hf_map_cid[ind];
01834                     int cur_lfid         = v2hf_map_lfid[ind];
01835 
01836                     sibhfs[nfpc * cidx + i] = create_halffacet( cur_cid, cur_lfid );
01837 
01838                     int scidx                       = ID_FROM_HANDLE( cur_cid ) - 1;
01839                     sibhfs[nfpc * scidx + cur_lfid] = create_halffacet( *cid, i );
01840                 }
01841             }
01842         }
01843     }
01844 
01845     return MB_SUCCESS;
01846 }
01847 
01848 ErrorCode HalfFacetRep::determine_incident_halffaces( Range& cells )
01849 {
01850     ErrorCode error;
01851     int index = get_index_in_lmap( *cells.begin() );
01852     int nvpc  = lConnMap3D[index].num_verts_in_cell;
01853 
01854     std::multimap< EntityHandle, EntityHandle > comps;
01855     HFacet nwhf;
01856 
01857     for( Range::iterator cid = cells.begin(); cid != cells.end(); ++cid )
01858     {
01859         EntityHandle cell = *cid;
01860         const EntityHandle* conn;
01861         error = mb->get_connectivity( *cid, conn, nvpc, true );MB_CHK_ERR( error );
01862 
01863         for( int i = 0; i < nvpc; ++i )
01864         {
01865             EntityHandle v = conn[i];
01866             int vidx       = ID_FROM_HANDLE( v ) - 1;
01867             HFacet hf      = v2hf[vidx];
01868 
01869             bool found = find_cell_in_component( v, cell, comps );
01870 
01871             if( hf == 0 && !found && ( v2hfs.empty() || ( v2hfs.find( v ) == v2hfs.end() ) ) )
01872             {
01873                 nwhf  = 0;
01874                 error = add_cells_of_single_component( v, cell, lConnMap3D[index].v2hf[i][0], comps, nwhf );MB_CHK_ERR( error );
01875 
01876                 v2hf[vidx] = nwhf;
01877             }
01878             else if( hf != 0 && !found )
01879             {
01880                 nwhf  = 0;
01881                 error = add_cells_of_single_component( v, cell, lConnMap3D[index].v2hf[i][0], comps, nwhf );MB_CHK_ERR( error );
01882 
01883                 v2hfs.insert( std::pair< EntityHandle, HFacet >( v, hf ) );
01884                 v2hfs.insert( std::pair< EntityHandle, HFacet >( v, nwhf ) );
01885                 v2hf[vidx] = 0;
01886             }
01887             else if( hf == 0 && !found && ( !v2hfs.empty() ) && ( v2hfs.find( v ) != v2hfs.end() ) )
01888             {
01889                 nwhf  = 0;
01890                 error = add_cells_of_single_component( v, cell, lConnMap3D[index].v2hf[i][0], comps, nwhf );MB_CHK_ERR( error );
01891                 v2hfs.insert( std::pair< EntityHandle, HFacet >( v, nwhf ) );
01892             }
01893         }
01894     }
01895 
01896     // error = print_tags(3);
01897 
01898     return MB_SUCCESS;
01899 }
01900 
01901 ErrorCode HalfFacetRep::determine_border_vertices( Range& cells, Tag isborder )
01902 {
01903     ErrorCode error;
01904     EntityHandle start_cell = *cells.begin();
01905     EntityType ctype        = mb->type_from_handle( *cells.begin() );
01906     int index               = get_index_in_lmap( start_cell );
01907     int nvpc                = lConnMap3D[index].num_verts_in_cell;
01908     int nfpc                = lConnMap3D[index].num_faces_in_cell;
01909 
01910     int val = 1;
01911 
01912     for( Range::iterator t = cells.begin(); t != cells.end(); ++t )
01913     {
01914 
01915         const EntityHandle* conn;
01916         error = mb->get_connectivity( *t, conn, nvpc, true );MB_CHK_ERR( error );
01917 
01918         int cidx = ID_FROM_HANDLE( *t ) - 1;
01919         for( int i = 0; i < nfpc; ++i )
01920         {
01921             HFacet hf            = sibhfs[nfpc * cidx + i];
01922             EntityHandle sib_cid = fid_from_halfacet( hf, ctype );
01923 
01924             if( sib_cid == 0 )
01925             {
01926                 int nvF = lConnMap3D[index].hf2v_num[i];
01927 
01928                 for( int j = 0; j < nvF; ++j )
01929                 {
01930                     int ind = lConnMap3D[index].hf2v[i][j];
01931                     error   = mb->tag_set_data( isborder, &conn[ind], 1, &val );MB_CHK_ERR( error );
01932                 }
01933             }
01934         }
01935     }
01936 
01937     return MB_SUCCESS;
01938 }
01939 
01940 /////////////////////////////////////////////////////////////////////////
01941 ErrorCode HalfFacetRep::add_cells_of_single_component( EntityHandle vid, EntityHandle curcid, int curlid,
01942                                                        std::multimap< EntityHandle, EntityHandle >& comps, HFacet& hf )
01943 {
01944     ErrorCode error;
01945     EntityType ctype = mb->type_from_handle( curcid );
01946     int index        = get_index_in_lmap( curcid );
01947     int nvpc         = lConnMap3D[index].num_verts_in_cell;
01948     int nfpc         = lConnMap3D[index].num_faces_in_cell;
01949 
01950     int Stksize = 0, count = -1;
01951     Stkcells[0] = curcid;
01952 
01953     hf = create_halffacet( curcid, curlid );
01954 
01955     EntityHandle cur_cid;
01956     while( Stksize >= 0 )
01957     {
01958         cur_cid = Stkcells[Stksize];
01959         Stksize -= 1;
01960 
01961         bool found = find_match_in_array( cur_cid, trackcells, count );
01962         if( !found )
01963         {
01964             count += 1;
01965             trackcells[count] = cur_cid;
01966 
01967             // Add the current cell
01968             comps.insert( std::pair< EntityHandle, EntityHandle >( vid, cur_cid ) );
01969         }
01970 
01971         // Connectivity of the cell
01972         const EntityHandle* conn;
01973         error = mb->get_connectivity( cur_cid, conn, nvpc, true );MB_CHK_ERR( error );
01974 
01975         // Local id of vid in the cell and the half-faces incident on it
01976         int lv = -1;
01977         for( int i = 0; i < nvpc; ++i )
01978         {
01979             if( conn[i] == vid )
01980             {
01981                 lv = i;
01982                 break;
01983             }
01984         }
01985         if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
01986         int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
01987         int cidx      = ID_FROM_HANDLE( cur_cid ) - 1;
01988 
01989         // Add new cells into the stack
01990         EntityHandle ngb;
01991         HFacet hf_ngb;
01992         for( int i = 0; i < nhf_thisv; ++i )
01993         {
01994             int ind = lConnMap3D[index].v2hf[lv][i];
01995             hf_ngb  = sibhfs[nfpc * cidx + ind];
01996             ngb     = fid_from_halfacet( hf_ngb, ctype );
01997 
01998             if( ngb )
01999             {
02000                 bool found_ent = find_match_in_array( ngb, trackcells, count );
02001 
02002                 if( !found_ent )
02003                 {
02004                     Stksize += 1;
02005                     Stkcells[Stksize] = ngb;
02006                 }
02007             }
02008             else
02009                 hf = create_halffacet( cur_cid, ind );
02010         }
02011     }
02012 
02013     // Change the visited faces to false
02014     for( int i = 0; i < Stksize; i++ )
02015         Stkcells[i] = 0;
02016 
02017     for( int i = 0; i <= count; i++ )
02018         trackcells[i] = 0;
02019 
02020     return MB_SUCCESS;
02021 }
02022 
02023 bool HalfFacetRep::find_cell_in_component( EntityHandle vid, EntityHandle cell,
02024                                            std::multimap< EntityHandle, EntityHandle >& comps )
02025 {
02026     bool found = false;
02027 
02028     if( comps.empty() ) return found;
02029 
02030     std::pair< std::multimap< EntityHandle, EntityHandle >::iterator,
02031                std::multimap< EntityHandle, EntityHandle >::iterator >
02032         rit;
02033 
02034     rit = comps.equal_range( vid );
02035 
02036     for( std::multimap< EntityHandle, EntityHandle >::iterator it = rit.first; it != rit.second; ++it )
02037     {
02038         if( it->second == cell )
02039         {
02040             found = true;
02041             break;
02042         }
02043     }
02044 
02045     return found;
02046 }
02047 
02048 ////////////////////////////////////////////////////////////////////////
02049 ErrorCode HalfFacetRep::get_up_adjacencies_vert_3d( EntityHandle vid, std::vector< EntityHandle >& adjents )
02050 {
02051     ErrorCode error;
02052     adjents.reserve( 20 );
02053     EntityType ctype = mb->type_from_handle( *_cells.begin() );
02054 
02055     // Obtain a half-face/s incident on v
02056     int vidx  = ID_FROM_HANDLE( vid ) - 1;
02057     HFacet hf = v2hf[vidx];
02058 
02059     std::vector< EntityHandle > start_cells;
02060     if( hf == 0 && ( v2hfs.find( vid ) != v2hfs.end() ) )  // Vertex is non-manifold
02061     {
02062         std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
02063             it_hes;
02064         it_hes = v2hfs.equal_range( vid );
02065 
02066         for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
02067         {
02068             start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
02069         }
02070     }
02071     else if( hf != 0 )
02072         start_cells.push_back( fid_from_halfacet( hf, ctype ) );
02073 
02074     if( start_cells.empty() ) return MB_SUCCESS;
02075 
02076     int index = get_index_in_lmap( start_cells[0] );
02077     int nvpc  = lConnMap3D[index].num_verts_in_cell;
02078     int nfpc  = lConnMap3D[index].num_faces_in_cell;
02079 
02080     for( int i = 0; i < (int)start_cells.size(); i++ )
02081         cellq[i] = start_cells[i];
02082 
02083     int qsize = start_cells.size();
02084     EntityHandle cur_cid;
02085     int num_qvals = 0;
02086 
02087     while( num_qvals < qsize )
02088     {
02089         cur_cid = cellq[num_qvals];
02090         num_qvals += 1;
02091 
02092         // Add the current cell to output adj vector
02093         adjents.push_back( cur_cid );
02094 
02095         // Connectivity of the cell
02096         const EntityHandle* conn;
02097         error = mb->get_connectivity( cur_cid, conn, nvpc, true );MB_CHK_ERR( error );
02098 
02099         // Local id of vid in the cell and the half-faces incident on it
02100         int lv = -1;
02101         for( int i = 0; i < nvpc; ++i )
02102         {
02103             if( conn[i] == vid )
02104             {
02105                 lv = i;
02106                 break;
02107             }
02108         };
02109         if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
02110         // Number of local half-faces incident on the current vertex
02111         int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
02112         int cidx      = ID_FROM_HANDLE( cur_cid ) - 1;
02113 
02114         // Add new cells into the stack
02115         EntityHandle ngb;
02116         for( int i = 0; i < nhf_thisv; ++i )
02117         {
02118             int ind = lConnMap3D[index].v2hf[lv][i];
02119             hf      = sibhfs[nfpc * cidx + ind];
02120             ngb     = fid_from_halfacet( hf, ctype );
02121 
02122             if( ngb )
02123             {
02124                 bool found_ent = find_match_in_array( ngb, cellq, qsize - 1 );
02125 
02126                 if( !found_ent )
02127                 {
02128                     cellq[qsize] = ngb;
02129                     qsize += 1;
02130                 }
02131             }
02132         }
02133     }
02134 
02135     // Change the visited faces to false
02136     for( int i = 0; i < qsize; i++ )
02137         cellq[i] = 0;
02138 
02139     return MB_SUCCESS;
02140 }
02141 
02142 ErrorCode HalfFacetRep::get_up_adjacencies_edg_3d( EntityHandle eid, std::vector< EntityHandle >& adjents,
02143                                                    std::vector< int >* leids )
02144 {
02145     ErrorCode error;
02146     EntityType ctype        = mb->type_from_handle( *_cells.begin() );
02147     EntityHandle start_cell = *_cells.begin();
02148     int index               = get_index_in_lmap( start_cell );
02149     int nvpc                = lConnMap3D[index].num_verts_in_cell;
02150     int nfpc                = lConnMap3D[index].num_faces_in_cell;
02151 
02152     adjents.reserve( 20 );
02153     if( leids != NULL ) leids->reserve( 20 );
02154 
02155     // Find the edge vertices
02156     const EntityHandle* econn;
02157     int num_conn = 0;
02158     error        = mb->get_connectivity( eid, econn, num_conn, true );MB_CHK_ERR( error );
02159 
02160     EntityHandle v_start = econn[0], v_end = econn[1];
02161     int v1idx = ID_FROM_HANDLE( v_start ) - 1;
02162     int v2idx = ID_FROM_HANDLE( v_end ) - 1;
02163 
02164     // Find an half-facets incident to each end vertex of the edge
02165     std::vector< EntityHandle > start_cells;
02166     HFacet hf1 = v2hf[v1idx];
02167     HFacet hf2 = v2hf[v2idx];
02168 
02169     if( hf1 == 0 && !v2hfs.empty() )
02170     {
02171         std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
02172             it_hes;
02173         it_hes = v2hfs.equal_range( v_start );
02174 
02175         for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
02176         {
02177             start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
02178         }
02179     }
02180     else if( hf1 != 0 )
02181     {
02182         start_cells.push_back( fid_from_halfacet( hf1, ctype ) );
02183     }
02184 
02185     if( hf2 == 0 && !v2hfs.empty() )
02186     {
02187         std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
02188             it_hes;
02189         it_hes = v2hfs.equal_range( v_end );
02190 
02191         for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
02192         {
02193             start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
02194         }
02195     }
02196     else if( hf2 != 0 )
02197     {
02198         start_cells.push_back( fid_from_halfacet( hf2, ctype ) );
02199     }
02200 
02201     if( start_cells.empty() ) return MB_SUCCESS;
02202 
02203     std::sort( start_cells.begin(), start_cells.end() );
02204     std::vector< EntityHandle >::iterator last = std::unique( start_cells.begin(), start_cells.end() );
02205     start_cells.erase( last, start_cells.end() );
02206 
02207     for( int i = 0; i < (int)start_cells.size(); i++ )
02208         cellq[i] = start_cells[i];
02209 
02210     int qsize     = start_cells.size();
02211     int num_qvals = 0;
02212 
02213     while( num_qvals < qsize )
02214     {
02215         EntityHandle cell_id = cellq[num_qvals];
02216         num_qvals += 1;
02217 
02218         const EntityHandle* conn;
02219         error = mb->get_connectivity( cell_id, conn, nvpc, true );MB_CHK_ERR( error );
02220 
02221         int lv0 = -1, lv1 = -1, lv = -1;
02222 
02223         // locate v_origin in poped out tet, check if v_end is in
02224         for( int i = 0; i < nvpc; i++ )
02225         {
02226             if( v_start == conn[i] )
02227             {
02228                 lv0 = i;
02229                 lv  = lv0;
02230             }
02231             else if( v_end == conn[i] )
02232             {
02233                 lv1 = i;
02234                 lv  = lv1;
02235             }
02236         }
02237 
02238         if( ( lv0 >= 0 ) && ( lv1 >= 0 ) )
02239         {
02240             adjents.push_back( cell_id );
02241             if( leids != NULL ) leids->push_back( lConnMap3D[index].lookup_leids[lv0][lv1] );
02242         }
02243 
02244         // push back new found unchecked incident tets of v_start
02245         int cidx = ID_FROM_HANDLE( cell_id ) - 1;
02246         if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
02247         int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
02248 
02249         for( int i = 0; i < nhf_thisv; i++ )
02250         {
02251             int ind          = lConnMap3D[index].v2hf[lv][i];
02252             HFacet hf        = sibhfs[nfpc * cidx + ind];
02253             EntityHandle ngb = fid_from_halfacet( hf, ctype );
02254 
02255             if( ngb )
02256             {
02257                 bool found_ent = find_match_in_array( ngb, &cellq[0], qsize - 1 );
02258 
02259                 if( !found_ent )
02260                 {
02261                     cellq[qsize] = ngb;
02262                     qsize += 1;
02263                 }
02264             }
02265         }
02266     }
02267 
02268     for( int i = 0; i < qsize; i++ )
02269         cellq[i] = 0;
02270 
02271     return MB_SUCCESS;
02272 }
02273 
02274 ErrorCode HalfFacetRep::get_up_adjacencies_edg_3d( EntityHandle cid, int leid, std::vector< EntityHandle >& adjents,
02275                                                    std::vector< int >* leids, std::vector< int >* adj_orients )
02276 {
02277     ErrorCode error;
02278     EntityType ctype = mb->type_from_handle( cid );
02279     int index        = get_index_in_lmap( cid );
02280     int nvpc         = lConnMap3D[index].num_verts_in_cell;
02281     int nfpc         = lConnMap3D[index].num_faces_in_cell;
02282     adjents.clear();
02283     adjents.reserve( 20 );
02284 
02285     if( leids != NULL )
02286     {
02287         leids->clear();
02288         leids->reserve( 20 );
02289     }
02290     if( adj_orients != NULL )
02291     {
02292         adj_orients->clear();
02293         adj_orients->reserve( 20 );
02294     }
02295 
02296     const EntityHandle* econn;
02297     error = mb->get_connectivity( cid, econn, nvpc, true );MB_CHK_ERR( error );
02298 
02299     // Get the end vertices of the edge <cid,leid>
02300     int id               = lConnMap3D[index].e2v[leid][0];
02301     EntityHandle v_start = econn[id];
02302     id                   = lConnMap3D[index].e2v[leid][1];
02303     EntityHandle v_end   = econn[id];
02304 
02305     int v1idx = ID_FROM_HANDLE( v_start ) - 1;
02306     int v2idx = ID_FROM_HANDLE( v_end ) - 1;
02307 
02308     // Find an half-facets incident to each end vertex of the edge
02309     std::vector< EntityHandle > start_cells;
02310     HFacet hf1 = v2hf[v1idx];
02311     HFacet hf2 = v2hf[v2idx];
02312 
02313     if( hf1 == 0 && !v2hfs.empty() )
02314     {
02315         std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
02316             it_hes;
02317         it_hes = v2hfs.equal_range( v_start );
02318 
02319         for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
02320         {
02321             start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
02322         }
02323     }
02324     else if( hf1 != 0 )
02325     {
02326         start_cells.push_back( fid_from_halfacet( hf1, ctype ) );
02327     }
02328 
02329     if( hf2 == 0 && !v2hfs.empty() )
02330     {
02331         std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
02332             it_hes;
02333         it_hes = v2hfs.equal_range( v_end );
02334 
02335         for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
02336         {
02337             start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
02338         }
02339     }
02340     else if( hf2 != 0 )
02341     {
02342         start_cells.push_back( fid_from_halfacet( hf2, ctype ) );
02343     }
02344 
02345     if( start_cells.empty() ) return MB_SUCCESS;
02346 
02347     std::sort( start_cells.begin(), start_cells.end() );
02348     std::vector< EntityHandle >::iterator last = std::unique( start_cells.begin(), start_cells.end() );
02349     start_cells.erase( last, start_cells.end() );
02350 
02351     for( int i = 0; i < (int)start_cells.size(); i++ )
02352         cellq[i] = start_cells[i];
02353 
02354     int qsize     = start_cells.size();
02355     int num_qvals = 0;
02356 
02357     while( num_qvals < qsize )
02358     {
02359         EntityHandle cell_id = cellq[num_qvals];
02360         num_qvals += 1;
02361 
02362         const EntityHandle* conn;
02363         error = mb->get_connectivity( cell_id, conn, nvpc, true );MB_CHK_ERR( error );
02364 
02365         int lv0 = -1, lv1 = -1, lv = -1;
02366 
02367         // locate v_origin in poped out tet, check if v_end is in
02368         for( int i = 0; i < nvpc; i++ )
02369         {
02370             if( v_start == conn[i] )
02371             {
02372                 lv0 = i;
02373                 lv  = lv0;
02374             }
02375             else if( v_end == conn[i] )
02376             {
02377                 lv1 = i;
02378                 lv  = lv1;
02379             }
02380         }
02381 
02382         if( ( lv0 >= 0 ) && ( lv1 >= 0 ) )
02383         {
02384             adjents.push_back( cell_id );
02385             if( leids != NULL ) leids->push_back( lConnMap3D[index].lookup_leids[lv0][lv1] );
02386 
02387             if( adj_orients != NULL )
02388             {
02389                 int cur_leid = lConnMap3D[index].lookup_leids[lv0][lv1];
02390                 int id1      = lConnMap3D[index].e2v[cur_leid][0];
02391                 int id2      = lConnMap3D[index].e2v[cur_leid][1];
02392                 if( ( v_start == conn[id1] ) && ( v_end == conn[id2] ) )
02393                     adj_orients->push_back( 1 );
02394                 else if( ( v_start == conn[id2] ) && ( v_end == conn[id1] ) )
02395                     adj_orients->push_back( 0 );
02396             }
02397         }
02398         if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
02399         // push back new found unchecked incident tets of v_start
02400         int cidx      = ID_FROM_HANDLE( cell_id ) - 1;
02401         int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
02402 
02403         for( int i = 0; i < nhf_thisv; i++ )
02404         {
02405             int ind          = lConnMap3D[index].v2hf[lv][i];
02406             HFacet hf        = sibhfs[nfpc * cidx + ind];
02407             EntityHandle ngb = fid_from_halfacet( hf, ctype );
02408 
02409             if( ngb )
02410             {
02411                 bool found_ent = find_match_in_array( ngb, &cellq[0], qsize - 1 );
02412 
02413                 if( !found_ent )
02414                 {
02415                     cellq[qsize] = ngb;
02416                     qsize += 1;
02417                 }
02418             }
02419         }
02420     }
02421 
02422     for( int i = 0; i < qsize; i++ )
02423         cellq[i] = 0;
02424 
02425     return MB_SUCCESS;
02426 }
02427 
02428 ErrorCode HalfFacetRep::get_up_adjacencies_edg_3d_comp( EntityHandle cid, int leid,
02429                                                         std::vector< EntityHandle >& adjents, std::vector< int >* leids,
02430                                                         std::vector< int >* adj_orients )
02431 {
02432     ErrorCode error;
02433     EntityType ctype = mb->type_from_handle( cid );
02434     int index        = get_index_in_lmap( cid );
02435     int nvpc         = lConnMap3D[index].num_verts_in_cell;
02436     int nfpc         = lConnMap3D[index].num_faces_in_cell;
02437     adjents.clear();
02438     adjents.reserve( 20 );
02439 
02440     if( leids != NULL )
02441     {
02442         leids->clear();
02443         leids->reserve( 20 );
02444     }
02445     if( adj_orients != NULL )
02446     {
02447         adj_orients->clear();
02448         adj_orients->reserve( 20 );
02449     }
02450 
02451     const EntityHandle* econn;
02452     error = mb->get_connectivity( cid, econn, nvpc );MB_CHK_ERR( error );
02453 
02454     // Get the end vertices of the edge <cid,leid>
02455     int id               = lConnMap3D[index].e2v[leid][0];
02456     EntityHandle v_start = econn[id];
02457     id                   = lConnMap3D[index].e2v[leid][1];
02458     EntityHandle v_end   = econn[id];
02459 
02460     int v1idx = ID_FROM_HANDLE( v_start ) - 1;
02461     int v2idx = ID_FROM_HANDLE( v_end ) - 1;
02462 
02463     // Find an half-facets incident to each end vertex of the edge
02464     std::vector< EntityHandle > start_cells;
02465     HFacet hf1 = v2hf[v1idx];
02466     HFacet hf2 = v2hf[v2idx];
02467 
02468     if( ( hf1 == 0 ) && ( v2hfs.find( v_start ) != v2hfs.end() ) && ( hf2 == 0 ) &&
02469         ( v2hfs.find( v_end ) != v2hfs.end() ) )
02470     {
02471         std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
02472             it_hes;
02473         it_hes = v2hfs.equal_range( v_start );
02474 
02475         for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
02476         {
02477             start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
02478         }
02479     }
02480     else
02481         return MB_SUCCESS;
02482 
02483     if( start_cells.empty() ) return MB_SUCCESS;
02484 
02485     // std::sort(start_cells.begin(), start_cells.end());
02486     //  std::vector<EntityHandle>::iterator last = std::unique(start_cells.begin(),
02487     //  start_cells.end()); start_cells.erase(last, start_cells.end());
02488 
02489     for( int c = 0; c < (int)start_cells.size(); c++ )
02490     {
02491         cellq[0] = start_cells[c];
02492 
02493         int qsize     = 1;
02494         int num_qvals = 0;
02495 
02496         while( num_qvals < qsize )
02497         {
02498             EntityHandle cell_id = cellq[num_qvals];
02499             num_qvals += 1;
02500 
02501             const EntityHandle* conn;
02502             error = mb->get_connectivity( cell_id, conn, nvpc );MB_CHK_ERR( error );
02503 
02504             int lv0 = -1, lv1 = -1, lv = -1;
02505 
02506             // locate v_origin in poped out tet, check if v_end is in
02507             for( int i = 0; i < nvpc; i++ )
02508             {
02509                 if( v_start == conn[i] )
02510                 {
02511                     lv0 = i;
02512                     lv  = lv0;
02513                 }
02514                 else if( v_end == conn[i] )
02515                 {
02516                     lv1 = i;
02517                     lv  = lv1;
02518                 }
02519             }
02520 
02521             if( ( lv0 >= 0 ) && ( lv1 >= 0 ) )
02522             {
02523                 adjents.push_back( cell_id );
02524                 if( leids != NULL ) leids->push_back( lConnMap3D[index].lookup_leids[lv0][lv1] );
02525 
02526                 if( adj_orients != NULL )
02527                 {
02528                     int cur_leid = lConnMap3D[index].lookup_leids[lv0][lv1];
02529                     int id1      = lConnMap3D[index].e2v[cur_leid][0];
02530                     int id2      = lConnMap3D[index].e2v[cur_leid][1];
02531                     if( ( v_start == conn[id1] ) && ( v_end == conn[id2] ) )
02532                         adj_orients->push_back( 1 );
02533                     else if( ( v_start == conn[id2] ) && ( v_end == conn[id1] ) )
02534                         adj_orients->push_back( 0 );
02535                 }
02536 
02537                 for( int i = 0; i < qsize; i++ )
02538                     cellq[i] = 0;
02539 
02540                 break;
02541             }
02542 
02543             // push back new found unchecked incident tets of v_start
02544             int cidx      = ID_FROM_HANDLE( cell_id ) - 1;
02545             int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
02546 
02547             for( int i = 0; i < nhf_thisv; i++ )
02548             {
02549                 int ind          = lConnMap3D[index].v2hf[lv][i];
02550                 HFacet hf        = sibhfs[nfpc * cidx + ind];
02551                 EntityHandle ngb = fid_from_halfacet( hf, ctype );
02552 
02553                 if( ngb )
02554                 {
02555                     bool found_ent = find_match_in_array( ngb, &cellq[0], qsize - 1 );
02556 
02557                     if( !found_ent )
02558                     {
02559                         cellq[qsize] = ngb;
02560                         qsize += 1;
02561                     }
02562                 }
02563             }
02564         }
02565 
02566         for( int i = 0; i < qsize; i++ )
02567             cellq[i] = 0;
02568     }
02569 
02570     return MB_SUCCESS;
02571 }
02572 
02573 ErrorCode HalfFacetRep::get_up_adjacencies_face_3d( EntityHandle fid, std::vector< EntityHandle >& adjents,
02574                                                     std::vector< int >* lfids )
02575 {
02576     ErrorCode error;
02577 
02578     EntityHandle cid = 0;
02579     int lid          = 0;
02580     bool found       = find_matching_halfface( fid, &cid, &lid );
02581 
02582     if( found )
02583     {
02584         error = get_up_adjacencies_face_3d( cid, lid, adjents, lfids );MB_CHK_ERR( error );
02585     }
02586 
02587     return MB_SUCCESS;
02588 }
02589 
02590 ErrorCode HalfFacetRep::get_up_adjacencies_face_3d( EntityHandle cid, int lfid, std::vector< EntityHandle >& adjents,
02591                                                     std::vector< int >* lfids )
02592 {
02593 
02594     EntityHandle start_cell = *_cells.begin();
02595     EntityType ctype        = mb->type_from_handle( start_cell );
02596     int index               = get_index_in_lmap( start_cell );
02597     int nfpc                = lConnMap3D[index].num_faces_in_cell;
02598 
02599     adjents.reserve( 4 );
02600     adjents.push_back( cid );
02601 
02602     if( lfids != NULL )
02603     {
02604         lfids->reserve( 4 );
02605         lfids->push_back( lfid );
02606     }
02607 
02608     int cidx            = ID_FROM_HANDLE( cid ) - 1;
02609     HFacet hf           = sibhfs[nfpc * cidx + lfid];
02610     EntityHandle sibcid = fid_from_halfacet( hf, ctype );
02611     int siblid          = lid_from_halffacet( hf );
02612 
02613     if( sibcid != 0 )
02614     {
02615         adjents.push_back( sibcid );
02616         if( lfids != NULL ) lfids->push_back( siblid );
02617     }
02618 
02619     return MB_SUCCESS;
02620 }
02621 
02622 bool HalfFacetRep::find_matching_implicit_edge_in_cell( EntityHandle eid, std::vector< EntityHandle >& cid,
02623                                                         std::vector< int >& leid )
02624 {
02625     ErrorCode error;
02626     EntityType ctype        = mb->type_from_handle( *_cells.begin() );
02627     EntityHandle start_cell = *_cells.begin();
02628     int index               = get_index_in_lmap( start_cell );
02629     int nvpc                = lConnMap3D[index].num_verts_in_cell;
02630     int nfpc                = lConnMap3D[index].num_faces_in_cell;
02631 
02632     // Find the edge vertices
02633     const EntityHandle* econn;
02634     int num_conn = 0;
02635     error        = mb->get_connectivity( eid, econn, num_conn, true );MB_CHK_ERR( error );
02636 
02637     EntityHandle v_start = econn[0], v_end = econn[1];
02638     int v1idx = ID_FROM_HANDLE( v_start ) - 1;
02639     int v2idx = ID_FROM_HANDLE( v_end ) - 1;
02640 
02641     // Find an incident cell to v_start
02642     std::vector< EntityHandle > start_cells;
02643     HFacet hf1 = v2hf[v1idx];
02644     HFacet hf2 = v2hf[v2idx];
02645 
02646     int ncomps1 = 0, ncomps2 = 0, ncomp;
02647     if( hf1 == 0 && !v2hfs.empty() )
02648     {
02649         std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
02650             it_hes;
02651         it_hes = v2hfs.equal_range( v_start );
02652 
02653         for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
02654         {
02655             start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
02656             ncomps1 += 1;
02657         }
02658     }
02659     else if( hf1 != 0 )
02660     {
02661         start_cells.push_back( fid_from_halfacet( hf1, ctype ) );
02662         ncomps1 += 1;
02663     }
02664 
02665     if( hf2 == 0 && !v2hfs.empty() )
02666     {
02667         std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
02668             it_hes;
02669         it_hes = v2hfs.equal_range( v_end );
02670 
02671         for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
02672         {
02673             start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
02674             ncomps2 += 1;
02675         }
02676     }
02677     else if( hf2 != 0 )
02678     {
02679         start_cells.push_back( fid_from_halfacet( hf2, ctype ) );
02680         ncomps2 += 1;
02681     }
02682 
02683     ncomp = std::min( ncomps1, ncomps2 );
02684 
02685     bool found = false;
02686     if( start_cells.empty() ) return found;
02687 
02688     for( int i = 0; i < (int)start_cells.size(); i++ )
02689         cellq[i] = start_cells[i];
02690 
02691     int qsize     = start_cells.size();
02692     int num_qvals = 0;
02693 
02694     while( num_qvals < qsize )
02695     {
02696         EntityHandle cell_id = cellq[num_qvals];
02697         num_qvals += 1;
02698 
02699         const EntityHandle* conn;
02700         error = mb->get_connectivity( cell_id, conn, nvpc, true );MB_CHK_ERR( error );
02701 
02702         int lv0 = -1, lv1 = -1, lv = -1;
02703 
02704         // locate v_origin in poped out tet, check if v_end is in
02705         for( int i = 0; i < nvpc; i++ )
02706         {
02707             if( v_start == conn[i] )
02708             {
02709                 lv0 = i;
02710                 lv  = lv0;
02711             }
02712             else if( v_end == conn[i] )
02713             {
02714                 lv1 = i;
02715                 lv  = lv1;
02716             }
02717         }
02718 
02719         if( ( lv0 >= 0 ) && ( lv1 >= 0 ) )
02720         {
02721             found = true;
02722             cid.push_back( cell_id );
02723             leid.push_back( lConnMap3D[index].lookup_leids[lv0][lv1] );
02724 
02725             if( (int)cid.size() == ncomp ) break;
02726         }
02727 
02728         // push back new found unchecked incident tets of v_start
02729         int cidx      = ID_FROM_HANDLE( cell_id ) - 1;
02730         int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
02731         if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
02732         for( int i = 0; i < nhf_thisv; i++ )
02733         {
02734             int ind          = lConnMap3D[index].v2hf[lv][i];
02735             HFacet hf        = sibhfs[nfpc * cidx + ind];
02736             EntityHandle ngb = fid_from_halfacet( hf, ctype );
02737 
02738             if( ngb )
02739             {
02740                 bool found_ent = find_match_in_array( ngb, &cellq[0], qsize - 1 );
02741 
02742                 if( !found_ent )
02743                 {
02744                     cellq[qsize] = ngb;
02745                     qsize += 1;
02746                 }
02747             }
02748         }
02749     }
02750 
02751     for( int i = 0; i < qsize; i++ )
02752         cellq[i] = 0;
02753 
02754     return found;
02755 }
02756 
02757 bool HalfFacetRep::find_matching_halfface( EntityHandle fid, EntityHandle* cid, int* lid )
02758 {
02759     ErrorCode error;
02760     EntityHandle start_cell = *_cells.begin();
02761     EntityType ctype        = mb->type_from_handle( start_cell );
02762     int index               = get_index_in_lmap( start_cell );
02763     int nvpc                = lConnMap3D[index].num_verts_in_cell;
02764     int nfpc                = lConnMap3D[index].num_faces_in_cell;
02765     EntityType ftype        = mb->type_from_handle( fid );
02766     int nvF                 = lConnMap2D[ftype - 2].num_verts_in_face;
02767 
02768     const EntityHandle* fid_verts;
02769     error = mb->get_connectivity( fid, fid_verts, nvF, true );MB_CHK_ERR( error );
02770 
02771     std::vector< EntityHandle > start_cells;
02772     int vidx, locfv0 = -1;
02773     HFacet hf = 0;
02774 
02775     for( int i = 0; i < nvF; i++ )
02776     {
02777         vidx = ID_FROM_HANDLE( fid_verts[i] ) - 1;
02778         hf   = v2hf[vidx];
02779         if( hf != 0 )
02780         {
02781             start_cells.push_back( fid_from_halfacet( hf, ctype ) );
02782             locfv0 = i;
02783             break;
02784         }
02785         else if( hf == 0 && v2hfs.find( fid_verts[i] ) != v2hfs.end() )
02786         {
02787             std::pair< std::multimap< EntityHandle, HFacet >::iterator,
02788                        std::multimap< EntityHandle, HFacet >::iterator >
02789                 it_hfs;
02790             it_hfs = v2hfs.equal_range( fid_verts[i] );
02791 
02792             for( std::multimap< EntityHandle, HFacet >::iterator it = it_hfs.first; it != it_hfs.second; ++it )
02793             {
02794                 start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
02795             }
02796             locfv0 = i;
02797             break;
02798         }
02799     }
02800 
02801     if( start_cells.empty() ) return false;
02802 
02803     EntityHandle cur_cid;
02804     bool found = false;
02805 
02806     int Stksize = 0, count = -1;
02807 
02808     for( int i = 0; i < (int)start_cells.size(); i++ )
02809         Stkcells[i] = start_cells[i];
02810 
02811     Stksize = start_cells.size() - 1;
02812 
02813     while( Stksize >= 0 )
02814     {
02815         cur_cid = Stkcells[Stksize];
02816         Stksize -= 1;
02817         count += 1;
02818         trackcells[count] = cur_cid;
02819 
02820         const EntityHandle* conn;
02821         error = mb->get_connectivity( cur_cid, conn, nvpc, true );MB_CHK_ERR( error );
02822 
02823         int lv[4] = { -1, -1, -1, -1 };
02824         int cnt   = 0;
02825         for( int i = 0; i < nvpc; i++ )
02826         {
02827             for( int j = 0; j < nvF; j++ )
02828                 if( conn[i] == fid_verts[j] )
02829                 {
02830                     lv[j] = i;
02831                     cnt += 1;
02832                 }
02833         }
02834         if( cnt == nvF )  // All face verts are part of the cell
02835         {
02836             found         = true;
02837             int nhf_thisv = lConnMap3D[index].v2hf_num[lv[locfv0]];
02838             int lfid      = -1;
02839             for( int i = 0; i < nhf_thisv; ++i )
02840             {
02841                 lfid     = lConnMap3D[index].v2hf[lv[locfv0]][i];
02842                 int lcnt = 0;
02843                 for( int j = 0; j < nvF; ++j )
02844                 {
02845                     for( int k = 0; k < nvF; ++k )
02846                     {
02847                         if( lv[k] == lConnMap3D[index].hf2v[lfid][j] ) lcnt += 1;
02848                     }
02849                 }
02850                 if( lcnt == nvF ) break;
02851             }
02852             cid[0] = cur_cid;
02853             lid[0] = lfid;
02854 
02855             break;
02856         }
02857         else
02858         {
02859             // Add other cells that are incident on fid_verts[0]
02860             if( locfv0 < 0 || lv[locfv0] < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
02861             int nhf_thisv = lConnMap3D[index].v2hf_num[lv[locfv0]];
02862             int cidx      = ID_FROM_HANDLE( cur_cid ) - 1;
02863 
02864             // Add new cells into the stack
02865             EntityHandle ngb;
02866             for( int i = 0; i < nhf_thisv; ++i )
02867             {
02868                 int ind = lConnMap3D[index].v2hf[lv[locfv0]][i];
02869                 hf      = sibhfs[nfpc * cidx + ind];
02870                 ngb     = fid_from_halfacet( hf, ctype );
02871 
02872                 if( ngb )
02873                 {
02874 
02875                     bool found_ent = find_match_in_array( ngb, trackcells, count );
02876 
02877                     if( !found_ent )
02878                     {
02879                         Stksize += 1;
02880                         Stkcells[Stksize] = ngb;
02881                     }
02882                 }
02883             }
02884         }
02885     }
02886 
02887     // Change the visited faces to false
02888     for( int i = 0; i < Stksize; i++ )
02889         Stkcells[i] = 0;
02890 
02891     for( int i = 0; i <= count; i++ )
02892         trackcells[i] = 0;
02893 
02894     return found;
02895 }
02896 
02897 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
02898 ErrorCode HalfFacetRep::get_neighbor_adjacencies_3d( EntityHandle cid, std::vector< EntityHandle >& adjents )
02899 {
02900     adjents.reserve( 20 );
02901     EntityType ctype = mb->type_from_handle( cid );
02902     int index        = get_index_in_lmap( cid );
02903     int nfpc         = lConnMap3D[index].num_faces_in_cell;
02904     int cidx         = ID_FROM_HANDLE( cid ) - 1;
02905 
02906     if( cid != 0 )
02907     {
02908         for( int lfid = 0; lfid < nfpc; ++lfid )
02909         {
02910             HFacet hf           = sibhfs[nfpc * cidx + lfid];
02911             EntityHandle sibcid = fid_from_halfacet( hf, ctype );
02912             if( sibcid != 0 ) adjents.push_back( sibcid );
02913         }
02914     }
02915 
02916     return MB_SUCCESS;
02917 }
02918 
02919 /////////////////////////////////////////////////////////////////////////////////////////////////
02920 ErrorCode HalfFacetRep::get_down_adjacencies_edg_3d( EntityHandle cid, std::vector< EntityHandle >& adjents )
02921 {
02922     // TODO: Try intersection without using std templates
02923     // Returns explicit edges, if any, of the face
02924     ErrorCode error;
02925     adjents.reserve( 20 );
02926     int index = get_index_in_lmap( cid );
02927     int nvpc  = lConnMap3D[index].num_verts_in_cell;
02928 
02929     const EntityHandle* conn;
02930     error = mb->get_connectivity( cid, conn, nvpc, true );MB_CHK_ERR( error );
02931 
02932     // Gather all the incident edges on each vertex of the face
02933     int ns = lConnMap3D[index].search_everts[0];
02934     std::vector< EntityHandle > temp;
02935     for( int i = 0; i < ns; i++ )
02936     {
02937         temp.clear();
02938         int lv0 = lConnMap3D[index].search_everts[i + 1];
02939         error   = get_up_adjacencies_1d( conn[lv0], temp );MB_CHK_ERR( error );
02940 
02941         int nle   = lConnMap3D[index].v2le[i][0];
02942         int count = 0;
02943         for( int j = 0; j < (int)temp.size(); j++ )
02944         {
02945             const EntityHandle* econn;
02946             int nvpe = 0;
02947             error    = mb->get_connectivity( temp[j], econn, nvpe, true );MB_CHK_ERR( error );
02948 
02949             for( int k = 0; k < nle; k++ )
02950             {
02951                 int lv1 = lConnMap3D[index].v2le[i][k + 1];
02952                 if( ( ( econn[0] == conn[lv0] ) && ( econn[1] == conn[lv1] ) ) ||
02953                     ( ( econn[1] == conn[lv0] ) && ( econn[0] == conn[lv1] ) ) )
02954                 {
02955                     adjents.push_back( temp[j] );
02956                     count += 1;
02957                 }
02958             }
02959             if( count == nle ) break;
02960         }
02961     }
02962     return MB_SUCCESS;
02963 }
02964 
02965 ErrorCode HalfFacetRep::get_down_adjacencies_face_3d( EntityHandle cid, std::vector< EntityHandle >& adjents )
02966 {
02967     // Returns explicit face, if any of the cell
02968     ErrorCode error;
02969     adjents.reserve( 10 );
02970     int index = get_index_in_lmap( cid );
02971     int nvpc  = lConnMap3D[index].num_verts_in_cell;
02972     int nfpc  = lConnMap3D[index].num_faces_in_cell;
02973 
02974     // Get the connectivity of the input cell
02975     const EntityHandle* conn;
02976     error = mb->get_connectivity( cid, conn, nvpc, true );MB_CHK_ERR( error );
02977 
02978     // Collect all the half-faces of the cell
02979     EntityHandle half_faces[6][4];
02980     for( int i = 0; i < nfpc; i++ )
02981     {
02982         int nvf = lConnMap3D[index].hf2v_num[i];
02983         for( int j = 0; j < nvf; j++ )
02984         {
02985             int ind          = lConnMap3D[index].hf2v[i][j];
02986             half_faces[i][j] = conn[ind];
02987         }
02988     }
02989 
02990     // Add two vertices for which the upward adjacencies are computed
02991     int search_verts[2];
02992     search_verts[0] = lConnMap3D[index].search_fverts[0];
02993     search_verts[1] = lConnMap3D[index].search_fverts[1];
02994 
02995     std::vector< EntityHandle > temp;
02996     temp.reserve( 20 );
02997     for( int i = 0; i < 2; i++ )
02998     {
02999         // Get the incident faces on the local vertex
03000         int lv = search_verts[i];
03001         temp.clear();
03002         error = get_up_adjacencies_vert_2d( conn[lv], temp );MB_CHK_ERR( error );
03003 
03004         if( temp.size() == 0 ) continue;
03005 
03006         // Get the half-faces incident on the local vertex and match it with the obtained faces
03007         int nhfthisv = lConnMap3D[index].v2hf_num[lv];
03008         for( int k = 0; k < (int)temp.size(); k++ )
03009         {
03010             const EntityHandle* fid_verts;
03011             int fsize = 0;
03012             error     = mb->get_connectivity( temp[k], fid_verts, fsize, true );MB_CHK_ERR( error );
03013 
03014             for( int j = 0; j < nhfthisv; j++ )
03015             {
03016                 // Collect all the vertices of this half-face
03017                 int idx = lConnMap3D[index].v2hf[lv][j];
03018                 int nvF = lConnMap3D[index].hf2v_num[idx];
03019 
03020                 if( fsize != nvF ) continue;
03021 
03022                 int direct, offset;
03023                 bool they_match = CN::ConnectivityMatch( &half_faces[idx][0], &fid_verts[0], nvF, direct, offset );
03024 
03025                 if( they_match )
03026                 {
03027                     bool found = false;
03028                     for( int p = 0; p < (int)adjents.size(); p++ )
03029                     {
03030                         if( adjents[p] == temp[k] )
03031                         {
03032                             found = true;
03033                             break;
03034                         }
03035                     }
03036                     if( !found ) adjents.push_back( temp[k] );
03037                 }
03038             }
03039         }
03040     }
03041 
03042     return MB_SUCCESS;
03043 }
03044 ////////////////////////////////////////////////////////////////////////////////
03045 ErrorCode HalfFacetRep::find_total_edges_faces_3d( Range cells, int* nedges, int* nfaces )
03046 {
03047     ErrorCode error;
03048     int index       = get_index_in_lmap( *cells.begin() );
03049     int nepc        = lConnMap3D[index].num_edges_in_cell;
03050     int nfpc        = lConnMap3D[index].num_faces_in_cell;
03051     int ncells      = cells.size();
03052     int total_edges = nepc * ncells;
03053     int total_faces = nfpc * ncells;
03054 
03055     std::vector< int > trackE( total_edges, 0 );
03056     std::vector< int > trackF( total_faces, 0 );
03057 
03058     std::vector< EntityHandle > inc_cids, sib_cids;
03059     std::vector< int > inc_leids, sib_lfids;
03060 
03061     for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
03062     {
03063         // Count edges
03064         for( int i = 0; i < nepc; i++ )
03065         {
03066             inc_cids.clear();
03067             inc_leids.clear();
03068 
03069             int id = nepc * ( cells.index( *it ) ) + i;
03070             if( !trackE[id] )
03071             {
03072                 error = get_up_adjacencies_edg_3d( *it, i, inc_cids, &inc_leids );MB_CHK_ERR( error );
03073 
03074                 total_edges -= inc_cids.size() - 1;
03075                 for( int j = 0; j < (int)inc_cids.size(); j++ )
03076                     trackE[nepc * ( cells.index( inc_cids[j] ) ) + inc_leids[j]] = 1;
03077             }
03078         }
03079 
03080         // Count faces
03081         for( int i = 0; i < nfpc; i++ )
03082         {
03083             sib_cids.clear();
03084             sib_lfids.clear();
03085 
03086             int id = nfpc * ( cells.index( *it ) ) + i;
03087             if( !trackF[id] )
03088             {
03089                 error = get_up_adjacencies_face_3d( *it, i, sib_cids, &sib_lfids );MB_CHK_ERR( error );
03090 
03091                 if( sib_cids.size() == 1 ) continue;
03092 
03093                 total_faces -= sib_cids.size() - 1;
03094                 trackF[nfpc * ( cells.index( sib_cids[1] ) ) + sib_lfids[1]] = 1;
03095             }
03096         }
03097     }
03098 
03099     nedges[0] = total_edges;
03100     nfaces[0] = total_faces;
03101 
03102     return MB_SUCCESS;
03103 }
03104 
03105 bool HalfFacetRep::find_match_in_array( EntityHandle ent, EntityHandle* ent_list, int count, bool get_index,
03106                                         int* index )
03107 {
03108     bool found = false;
03109     for( int i = 0; i <= count; i++ )
03110     {
03111         if( ent == ent_list[i] )
03112         {
03113             found = true;
03114             if( get_index ) *index = i;
03115             break;
03116         }
03117     }
03118 
03119     return found;
03120 }
03121 
03122 ErrorCode HalfFacetRep::get_half_facet_in_comp( EntityHandle cid, int leid, std::vector< EntityHandle >& ents,
03123                                                 std::vector< int >& lids, std::vector< int >& lfids )
03124 {
03125     ErrorCode error;
03126     ents.clear();
03127     lids.clear();
03128     EntityType ctype = mb->type_from_handle( cid );
03129     int index        = get_index_in_lmap( *_cells.begin() );
03130     int nfpc         = lConnMap3D[index].num_faces_in_cell;
03131     int nvpc         = lConnMap3D[index].num_verts_in_cell;
03132 
03133     // Get all incident cells
03134     std::vector< EntityHandle > adjents;
03135     std::vector< int > adjlids;
03136     error = get_up_adjacencies_edg_3d( cid, leid, adjents, &adjlids );MB_CHK_ERR( error );
03137 
03138     // Get the end vertices of the edge <cid,leid>
03139     const EntityHandle* econn;
03140     error = mb->get_connectivity( cid, econn, nvpc, true );MB_CHK_ERR( error );
03141     int id              = lConnMap3D[index].e2v[leid][0];
03142     EntityHandle vstart = econn[id];
03143     id                  = lConnMap3D[index].e2v[leid][1];
03144     EntityHandle vend   = econn[id];
03145 
03146     std::vector< int > mark( adjents.size(), 0 );
03147     int count = 0;
03148 
03149     for( int k = 0; k < (int)adjents.size(); k++ )
03150     {
03151 
03152         if( mark[k] != 0 ) continue;
03153 
03154         count += 1;
03155         mark[k] = count;
03156 
03157         // Loop over each half-face incident on this local edge
03158         for( int i = 0; i < 2; i++ )
03159         {
03160             EntityHandle cur_cell = adjents[k];
03161             int cur_leid          = adjlids[k];
03162 
03163             int lface = i;
03164 
03165             while( true )
03166             {
03167                 int cidx = ID_FROM_HANDLE( cur_cell ) - 1;
03168                 int lfid = lConnMap3D[index].e2hf[cur_leid][lface];
03169 
03170                 HFacet hf = sibhfs[nfpc * cidx + lfid];
03171                 cur_cell  = fid_from_halfacet( hf, ctype );
03172                 lfid      = lid_from_halffacet( hf );
03173 
03174                 // Check if loop reached starting cell or a boundary
03175                 if( ( cur_cell == adjents[k] ) || ( cur_cell == 0 ) ) break;
03176 
03177                 const EntityHandle* sib_conn;
03178                 error = mb->get_connectivity( cur_cell, sib_conn, nvpc, true );MB_CHK_ERR( error );
03179 
03180                 // Find the local edge id wrt to sibhf
03181                 int nv_curF = lConnMap3D[index].hf2v_num[lfid];
03182                 int lv0 = -1, lv1 = -1, idx = -1;
03183                 for( int j = 0; j < nv_curF; j++ )
03184                 {
03185                     idx = lConnMap3D[index].hf2v[lfid][j];
03186                     if( vstart == sib_conn[idx] ) lv0 = idx;
03187                     if( vend == sib_conn[idx] ) lv1 = idx;
03188                 }
03189 
03190                 assert( ( lv0 >= 0 ) && ( lv1 >= 0 ) );
03191                 cur_leid = lConnMap3D[index].lookup_leids[lv0][lv1];
03192 
03193                 int chk_lfid = lConnMap3D[index].e2hf[cur_leid][0];
03194 
03195                 if( lfid == chk_lfid )
03196                     lface = 1;
03197                 else
03198                     lface = 0;
03199 
03200                 int ind   = std::find( adjents.begin(), adjents.end(), cur_cell ) - adjents.begin();
03201                 mark[ind] = count;
03202             }
03203 
03204             // Loop back
03205             if( cur_cell != 0 ) break;
03206         }
03207     }
03208 
03209     // Loop over again to find cells on the boundary
03210     for( int c = 0; c < count; c++ )
03211     {
03212         for( int i = 0; i < (int)adjents.size(); i++ )
03213         {
03214             if( mark[i] == c + 1 )
03215             {
03216                 int cidx = ID_FROM_HANDLE( adjents[i] ) - 1;
03217                 for( int j = 0; j < nfpc; j++ )
03218                 {
03219                     HFacet hf         = sibhfs[nfpc * cidx + j];
03220                     EntityHandle cell = fid_from_halfacet( hf, ctype );
03221                     if( cell == 0 )
03222                     {
03223                         ents.push_back( adjents[i] );
03224                         lids.push_back( adjlids[i] );
03225                         lfids.push_back( j );
03226                         break;
03227                     }
03228                 }
03229             }
03230         }
03231     }
03232 
03233     return MB_SUCCESS;
03234 }
03235 
03236 ///////////////////////////////////////////////////////////////////////////////////////////
03237 ErrorCode HalfFacetRep::resize_hf_maps( EntityHandle start_vert, int nverts, EntityHandle start_edge, int nedges,
03238                                         EntityHandle start_face, int nfaces, EntityHandle start_cell, int ncells )
03239 {
03240     int nwsz = 0, insz = 0;
03241     if( nedges )
03242     {
03243         if( ID_FROM_HANDLE( ( *( _edges.end() - 1 ) + 1 ) ) != ID_FROM_HANDLE( start_edge ) )
03244             nwsz = ( ID_FROM_HANDLE( start_edge ) - ID_FROM_HANDLE( *_edges.end() ) + nedges ) * 2;
03245         else
03246             nwsz = nedges * 2;
03247         insz = sibhvs.size();
03248         sibhvs.resize( insz + nwsz, 0 );
03249 
03250         if( v2hv.empty() )
03251         {
03252             if( ( !v2he.empty() ) )
03253                 insz = v2he.size();
03254             else if( ( !v2hf.empty() ) )
03255                 insz = v2hf.size();
03256             else
03257                 MB_SET_ERR( MB_FAILURE, "Trying to resize ahf maps for a mesh with no edges, faces and cells" );
03258         }
03259         else
03260             insz = v2hv.size();
03261 
03262         if( ID_FROM_HANDLE( *( _verts.end() - 1 ) + 1 ) != ID_FROM_HANDLE( start_vert ) )
03263             nwsz = ID_FROM_HANDLE( start_vert ) - ID_FROM_HANDLE( *_verts.end() ) + nverts;
03264         else
03265             nwsz = nverts;
03266         v2hv.resize( insz + nwsz, 0 );
03267     }
03268 
03269     if( nfaces )
03270     {
03271         EntityType ftype = mb->type_from_handle( *_faces.begin() );
03272         int nepf         = lConnMap2D[ftype - 2].num_verts_in_face;
03273 
03274         if( ID_FROM_HANDLE( ( *( _faces.end() - 1 ) + 1 ) ) != ID_FROM_HANDLE( start_face ) )
03275             nwsz = ( ID_FROM_HANDLE( start_face ) - ID_FROM_HANDLE( *_faces.end() ) + nfaces ) * nepf;
03276         else
03277             nwsz = nfaces * nepf;
03278         insz = sibhes.size();
03279         sibhes.resize( insz + nwsz, 0 );
03280 
03281         if( ID_FROM_HANDLE( *( _verts.end() - 1 ) + 1 ) != ID_FROM_HANDLE( start_vert ) )
03282             nwsz = ID_FROM_HANDLE( start_vert ) - ID_FROM_HANDLE( *_verts.end() ) + nverts;
03283         else
03284             nwsz = nverts;
03285         insz = v2he.size();
03286         v2he.resize( insz + nwsz, 0 );
03287     }
03288 
03289     if( ncells )
03290     {
03291         int index = get_index_in_lmap( *_cells.begin() );
03292         int nfpc  = lConnMap3D[index].num_faces_in_cell;
03293 
03294         if( ID_FROM_HANDLE( ( *( _cells.end() - 1 ) + 1 ) ) != ID_FROM_HANDLE( start_cell ) )
03295             nwsz = ( ID_FROM_HANDLE( start_cell ) - ID_FROM_HANDLE( *_cells.end() ) + ncells ) * nfpc;
03296         else
03297             nwsz = ncells * nfpc;
03298         insz = sibhfs.size();
03299         sibhfs.resize( insz + nwsz, 0 );
03300 
03301         if( ID_FROM_HANDLE( *( _verts.end() - 1 ) + 1 ) != ID_FROM_HANDLE( start_vert ) )
03302             nwsz = ID_FROM_HANDLE( start_vert ) - ID_FROM_HANDLE( *_verts.end() ) + nverts;
03303         else
03304             nwsz = nverts;
03305         insz = v2hf.size();
03306         v2hf.resize( insz + nwsz, 0 );
03307     }
03308 
03309     return MB_SUCCESS;
03310 }
03311 
03312 bool HalfFacetRep::check_nonmanifold_vertices( EntityType type, EntityHandle vid )
03313 {
03314     bool status = false;
03315     if( type == MBTRI || type == MBQUAD )
03316     {
03317         HFacet hf = v2he[ID_FROM_HANDLE( vid ) - 1];
03318         if( hf == 0 && ( v2hes.find( vid ) != v2hes.end() ) ) status = true;
03319     }
03320     else if( type == MBTET || type == MBHEX )
03321     {
03322         HFacet hf = v2hf[ID_FROM_HANDLE( vid ) - 1];
03323         if( hf == 0 && ( v2hfs.find( vid ) != v2hfs.end() ) ) status = true;
03324     }
03325     else
03326         MB_SET_ERR( MB_FAILURE, "Requesting non-manifold vertex checks for either (1) 1D mesh or "
03327                                 "(2) not-implemented entity types" );
03328 
03329     return status;
03330 }
03331 
03332 ErrorCode HalfFacetRep::get_sibling_map( EntityType type, EntityHandle ent, EntityHandle* sib_entids, int* sib_lids,
03333                                          int num_halffacets )
03334 {
03335 
03336     if( type == MBEDGE )
03337     {
03338         if( num_halffacets != 2 ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halfvertices." );
03339 
03340         int eidx = ID_FROM_HANDLE( ent ) - 1;
03341         for( int i = 0; i < 2; i++ )
03342         {
03343             HFacet hf     = sibhvs[2 * eidx + i];
03344             sib_entids[i] = fid_from_halfacet( hf, MBEDGE );
03345             sib_lids[i]   = lid_from_halffacet( hf );
03346         }
03347     }
03348     else if( type == MBTRI || type == MBQUAD )
03349     {
03350         int nepf = lConnMap2D[type - 2].num_verts_in_face;
03351 
03352         if( num_halffacets != nepf ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halfedges." );
03353 
03354         int fidx = ID_FROM_HANDLE( ent ) - 1;
03355         for( int i = 0; i < nepf; i++ )
03356         {
03357             HFacet hf     = sibhes[nepf * fidx + i];
03358             sib_entids[i] = fid_from_halfacet( hf, type );
03359             sib_lids[i]   = lid_from_halffacet( hf );
03360         }
03361     }
03362     else
03363     {
03364         int idx  = get_index_in_lmap( *_cells.begin() );
03365         int nfpc = lConnMap3D[idx].num_faces_in_cell;
03366 
03367         if( num_halffacets != nfpc ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halffaces." );
03368 
03369         int cidx = ID_FROM_HANDLE( ent ) - 1;
03370         for( int i = 0; i < nfpc; i++ )
03371         {
03372             HFacet hf     = sibhfs[nfpc * cidx + i];
03373             sib_entids[i] = fid_from_halfacet( hf, type );
03374             sib_lids[i]   = lid_from_halffacet( hf );
03375         }
03376     }
03377     return MB_SUCCESS;
03378 }
03379 
03380 ErrorCode HalfFacetRep::get_sibling_map( EntityType type, EntityHandle ent, int lid, EntityHandle& sib_entid,
03381                                          int& sib_lid )
03382 {
03383 
03384     if( type == MBEDGE )
03385     {
03386         int eidx  = ID_FROM_HANDLE( ent ) - 1;
03387         HFacet hf = sibhvs[2 * eidx + lid];
03388         sib_entid = fid_from_halfacet( hf, MBEDGE );
03389         sib_lid   = lid_from_halffacet( hf );
03390     }
03391     else if( type == MBTRI || type == MBQUAD )
03392     {
03393         int nepf  = lConnMap2D[type - 2].num_verts_in_face;
03394         int fidx  = ID_FROM_HANDLE( ent ) - 1;
03395         HFacet hf = sibhes[nepf * fidx + lid];
03396         sib_entid = fid_from_halfacet( hf, type );
03397         sib_lid   = lid_from_halffacet( hf );
03398     }
03399     else
03400     {
03401         int idx   = get_index_in_lmap( *_cells.begin() );
03402         int nfpc  = lConnMap3D[idx].num_faces_in_cell;
03403         int cidx  = ID_FROM_HANDLE( ent ) - 1;
03404         HFacet hf = sibhfs[nfpc * cidx + lid];
03405         sib_entid = fid_from_halfacet( hf, type );
03406         sib_lid   = lid_from_halffacet( hf );
03407     }
03408     return MB_SUCCESS;
03409 }
03410 
03411 ErrorCode HalfFacetRep::set_sibling_map( EntityType type, EntityHandle ent, EntityHandle* set_entids, int* set_lids,
03412                                          int num_halffacets )
03413 {
03414 
03415     if( type == MBEDGE )
03416     {
03417         if( num_halffacets != 2 ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halfvertices" );
03418 
03419         int eidx = ID_FROM_HANDLE( ent ) - 1;
03420         for( int i = 0; i < 2; i++ )
03421         {
03422             sibhvs[2 * eidx + i] = create_halffacet( set_entids[i], set_lids[i] );
03423         }
03424     }
03425     else if( type == MBTRI || type == MBQUAD )
03426     {
03427         int nepf = lConnMap2D[type - 2].num_verts_in_face;
03428         if( num_halffacets != nepf ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halfedges." );
03429 
03430         int fidx = ID_FROM_HANDLE( ent ) - 1;
03431         for( int i = 0; i < nepf; i++ )
03432         {
03433             sibhes[nepf * fidx + i] = create_halffacet( set_entids[i], set_lids[i] );
03434         }
03435     }
03436     else
03437     {
03438         int idx  = get_index_in_lmap( *_cells.begin() );
03439         int nfpc = lConnMap3D[idx].num_faces_in_cell;
03440         if( num_halffacets != nfpc ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halffaces." );
03441 
03442         int cidx = ID_FROM_HANDLE( ent ) - 1;
03443         for( int i = 0; i < nfpc; i++ )
03444         {
03445             sibhfs[nfpc * cidx + i] = create_halffacet( set_entids[i], set_lids[i] );
03446         }
03447     }
03448 
03449     return MB_SUCCESS;
03450 }
03451 
03452 ErrorCode HalfFacetRep::set_sibling_map( EntityType type, EntityHandle ent, int lid, EntityHandle& set_entid,
03453                                          int& set_lid )
03454 {
03455 
03456     if( type == MBEDGE )
03457     {
03458         int eidx               = ID_FROM_HANDLE( ent ) - 1;
03459         sibhvs[2 * eidx + lid] = create_halffacet( set_entid, set_lid );
03460     }
03461     else if( type == MBTRI || type == MBQUAD )
03462     {
03463         int nepf                  = lConnMap2D[type - 2].num_verts_in_face;
03464         int fidx                  = ID_FROM_HANDLE( ent ) - 1;
03465         sibhes[nepf * fidx + lid] = create_halffacet( set_entid, set_lid );
03466     }
03467     else
03468     {
03469         int idx                   = get_index_in_lmap( *_cells.begin() );
03470         int nfpc                  = lConnMap3D[idx].num_faces_in_cell;
03471         int cidx                  = ID_FROM_HANDLE( ent ) - 1;
03472         sibhfs[nfpc * cidx + lid] = create_halffacet( set_entid, set_lid );
03473     }
03474 
03475     return MB_SUCCESS;
03476 }
03477 
03478 ErrorCode HalfFacetRep::get_incident_map( EntityType type, EntityHandle vid, std::vector< EntityHandle >& inci_entid,
03479                                           std::vector< int >& inci_lid )
03480 {
03481     inci_entid.clear();
03482     inci_lid.clear();
03483 
03484     if( type == MBEDGE )
03485     {
03486         HFacet hf = v2hv[ID_FROM_HANDLE( vid ) - 1];
03487         inci_entid.push_back( fid_from_halfacet( hf, type ) );
03488         inci_lid.push_back( lid_from_halffacet( hf ) );
03489     }
03490     else if( type == MBTRI || type == MBQUAD )
03491     {
03492         HFacet hf = v2he[ID_FROM_HANDLE( vid ) - 1];
03493         if( hf == 0 && ( v2hes.find( vid ) != v2hes.end() ) )
03494         {
03495             std::pair< std::multimap< EntityHandle, HFacet >::iterator,
03496                        std::multimap< EntityHandle, HFacet >::iterator >
03497                 it_hes;
03498             it_hes = v2hes.equal_range( vid );
03499 
03500             for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
03501             {
03502                 inci_entid.push_back( fid_from_halfacet( it->second, type ) );
03503                 inci_lid.push_back( lid_from_halffacet( it->second ) );
03504             }
03505         }
03506         else if( hf != 0 )
03507         {
03508             inci_entid.push_back( fid_from_halfacet( hf, type ) );
03509             inci_lid.push_back( lid_from_halffacet( hf ) );
03510         }
03511         else if( hf == 0 && ( v2hes.find( vid ) == v2hes.end() ) )
03512         {
03513             inci_entid.push_back( fid_from_halfacet( hf, type ) );
03514             inci_lid.push_back( lid_from_halffacet( hf ) );
03515         }
03516     }
03517     else
03518     {
03519         HFacet hf = v2hf[ID_FROM_HANDLE( vid ) - 1];
03520         if( hf == 0 && ( v2hfs.find( vid ) != v2hfs.end() ) )
03521         {
03522             std::pair< std::multimap< EntityHandle, HFacet >::iterator,
03523                        std::multimap< EntityHandle, HFacet >::iterator >
03524                 it_hes;
03525             it_hes = v2hfs.equal_range( vid );
03526 
03527             for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
03528             {
03529                 inci_entid.push_back( fid_from_halfacet( it->second, type ) );
03530                 inci_lid.push_back( lid_from_halffacet( it->second ) );
03531             }
03532         }
03533         else if( hf != 0 )
03534         {
03535             inci_entid.push_back( fid_from_halfacet( hf, type ) );
03536             inci_lid.push_back( lid_from_halffacet( hf ) );
03537         }
03538         else if( hf == 0 && ( v2hfs.find( vid ) == v2hfs.end() ) )
03539         {
03540             inci_entid.push_back( fid_from_halfacet( hf, type ) );
03541             inci_lid.push_back( lid_from_halffacet( hf ) );
03542         }
03543     }
03544 
03545     return MB_SUCCESS;
03546 }
03547 
03548 ErrorCode HalfFacetRep::set_incident_map( EntityType type, EntityHandle vid, std::vector< EntityHandle >& set_entid,
03549                                           std::vector< int >& set_lid )
03550 {
03551     if( type == MBEDGE ) { v2hv[ID_FROM_HANDLE( vid ) - 1] = create_halffacet( set_entid[0], set_lid[0] ); }
03552     else if( type == MBTRI || type == MBQUAD )
03553     {
03554         if( set_entid.size() == 1 )
03555             v2he[ID_FROM_HANDLE( vid ) - 1] = create_halffacet( set_entid[0], set_lid[0] );
03556         else
03557         {
03558             HFacet hf = 0;
03559             for( int i = 0; i < (int)set_entid.size(); i++ )
03560             {
03561                 hf = create_halffacet( set_entid[i], set_lid[i] );
03562                 v2hes.insert( std::pair< EntityHandle, HFacet >( vid, hf ) );
03563             }
03564         }
03565     }
03566     else
03567     {
03568         if( set_entid.size() == 1 )
03569             v2hf[ID_FROM_HANDLE( vid ) - 1] = create_halffacet( set_entid[0], set_lid[0] );
03570         else
03571         {
03572             HFacet hf = v2hf[ID_FROM_HANDLE( vid ) - 1];
03573             if( hf != 0 ) { v2hf[ID_FROM_HANDLE( vid ) - 1] = 0; }
03574             for( int i = 0; i < (int)set_entid.size(); i++ )
03575             {
03576                 hf = create_halffacet( set_entid[i], set_lid[i] );
03577                 v2hfs.insert( std::pair< EntityHandle, HFacet >( vid, hf ) );
03578             }
03579         }
03580     }
03581 
03582     return MB_SUCCESS;
03583 }
03584 
03585 ErrorCode HalfFacetRep::get_entity_ranges( Range& verts, Range& edges, Range& faces, Range& cells )
03586 {
03587     verts = _verts;
03588     edges = _edges;
03589     faces = _faces;
03590     cells = _cells;
03591     return MB_SUCCESS;
03592 }
03593 
03594 ErrorCode HalfFacetRep::update_entity_ranges( EntityHandle fileset )
03595 {
03596     ErrorCode error;
03597 
03598     error = mb->get_entities_by_dimension( fileset, 0, _verts, true );MB_CHK_ERR( error );
03599 
03600     error = mb->get_entities_by_dimension( fileset, 1, _edges, true );MB_CHK_ERR( error );
03601 
03602     error = mb->get_entities_by_dimension( fileset, 2, _faces, true );MB_CHK_ERR( error );
03603 
03604     error = mb->get_entities_by_dimension( fileset, 3, _cells, true );MB_CHK_ERR( error );
03605 
03606     return MB_SUCCESS;
03607 }
03608 
03609 void HalfFacetRep::get_memory_use( unsigned long long& entity_total, unsigned long long& memory_total )
03610 {
03611     entity_total = memory_total = 0;
03612     // 1D
03613     if( !v2hv.empty() ) entity_total += v2hv.capacity() * sizeof( HFacet ) + sizeof( v2hv );
03614     if( !sibhvs.empty() ) entity_total += sibhvs.capacity() * sizeof( HFacet ) + sizeof( sibhvs );
03615 
03616     // 2D
03617     if( !v2he.empty() ) entity_total += v2he.capacity() * sizeof( HFacet ) + sizeof( v2he );
03618     if( !sibhes.empty() ) entity_total += sibhes.capacity() * sizeof( HFacet ) + sizeof( sibhes );
03619 
03620     // 3D
03621     if( !v2hf.empty() ) entity_total += v2hf.capacity() * sizeof( HFacet ) + sizeof( v2hf );
03622     if( !sibhfs.empty() ) entity_total += sibhfs.capacity() * sizeof( HFacet ) + sizeof( sibhfs );
03623 
03624     memory_total = entity_total;
03625 }
03626 
03627 HFacet HalfFacetRep::create_halffacet( EntityHandle handle, int lid )
03628 {
03629     EntityID fid = ID_FROM_HANDLE( handle );
03630     return CREATE_HALFFACET( lid, fid );
03631 }
03632 
03633 EntityHandle HalfFacetRep::fid_from_halfacet( const HFacet facet, EntityType type )
03634 {
03635     EntityID id         = FID_FROM_HALFFACET( facet );
03636     EntityHandle handle = 0;
03637     if( id == 0 ) return handle;
03638 
03639     ErrorCode error = mb->handle_from_id( type, id, handle );MB_CHK_ERR( error );
03640     return handle;
03641 }
03642 
03643 int HalfFacetRep::lid_from_halffacet( const HFacet facet )
03644 {
03645     if( facet == 0 )
03646         return 0;
03647     else
03648         return LID_FROM_HALFFACET( facet );
03649 }
03650 
03651 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines