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