MOAB: Mesh Oriented datABase  (version 5.2.1)
NestedRefine.cpp
Go to the documentation of this file.
00001 #include "moab/NestedRefine.hpp"
00002 #include "moab/NestedRefineTemplates.hpp"
00003 #include "moab/HalfFacetRep.hpp"
00004 #include "moab/CpuTimer.hpp"
00005 #include "moab/ReadUtilIface.hpp"
00006 #include "Internals.hpp"
00007 #include "MBTagConventions.hpp"
00008 
00009 #ifdef MOAB_HAVE_MPI
00010 #include "moab/ParallelComm.hpp"
00011 #include "moab/ParallelMergeMesh.hpp"
00012 #include "moab/Skinner.hpp"
00013 #endif
00014 
00015 #include <iostream>
00016 #include <assert.h>
00017 #include <vector>
00018 #include <limits>
00019 #include <cmath>
00020 
00021 namespace moab
00022 {
00023 
00024 NestedRefine::NestedRefine( Core* impl, ParallelComm* comm, EntityHandle rset )
00025     : mbImpl( impl ), pcomm( comm ), _rset( rset )
00026 {
00027     ErrorCode error;
00028     assert( NULL != impl );
00029 
00030 #ifdef MOAB_HAVE_MPI
00031     // Get the Parallel Comm instance to prepare all new sets to work in parallel
00032     // in case the user did not provide any arguments
00033     if( !comm ) pcomm = moab::ParallelComm::get_pcomm( mbImpl, 0 );
00034 #endif
00035     error = initialize();
00036     if( error != MB_SUCCESS )
00037     {
00038         std::cout << "Error initializing NestedRefine\n" << std::endl;
00039         exit( 1 );
00040     }
00041 }
00042 
00043 NestedRefine::~NestedRefine()
00044 {
00045 #ifdef MOAB_HAVE_AHF
00046     ahf = NULL;
00047 #else
00048     delete ahf;
00049 #endif
00050     delete tm;
00051 }
00052 
00053 ErrorCode NestedRefine::initialize()
00054 {
00055     ErrorCode error;
00056 
00057     tm = new CpuTimer();
00058     if( !tm ) return MB_MEMORY_ALLOCATION_FAILED;
00059 
00060 #ifdef MOAB_HAVE_AHF
00061     ahf = mbImpl->a_half_facet_rep();
00062 #else
00063     ahf = new HalfFacetRep( mbImpl, pcomm, _rset, true );
00064     if( !ahf ) return MB_MEMORY_ALLOCATION_FAILED;
00065 #endif
00066 
00067     // Check for mixed entity type
00068     bool chk_mixed = ahf->check_mixed_entity_type();
00069     if( chk_mixed ) MB_SET_ERR( MB_NOT_IMPLEMENTED, "Encountered a mesh with mixed entity types" );
00070 
00071     error = ahf->initialize();MB_CHK_ERR( error );
00072     error = ahf->get_entity_ranges( _inverts, _inedges, _infaces, _incells );MB_CHK_ERR( error );
00073 
00074     // Check for supported entity type
00075     if( !_incells.empty() )
00076     {
00077         EntityType type = mbImpl->type_from_handle( _incells[0] );
00078         if( type != MBTET && type != MBHEX )
00079             MB_SET_ERR( MB_FAILURE, "Not supported 3D entity types: MBPRISM, MBPYRAMID, MBKNIFE, MBPOLYHEDRON" );
00080 
00081         meshdim    = 3;
00082         elementype = type;
00083     }
00084     else if( !_infaces.empty() )
00085     {
00086         EntityType type = mbImpl->type_from_handle( _infaces[0] );
00087         if( type == MBPOLYGON ) MB_SET_ERR( MB_FAILURE, "Not supported 2D entity type: POLYGON" );
00088 
00089         meshdim    = 2;
00090         elementype = type;
00091     }
00092     else if( !_inedges.empty() )
00093     {
00094         meshdim    = 1;
00095         elementype = MBEDGE;
00096     }
00097     else
00098         MB_SET_ERR( MB_NOT_IMPLEMENTED, "Encountered a mixed-dimensional or invalid mesh" );
00099 
00100     // Initialize std::map to get indices of degrees.
00101     deg_index[2] = 0;
00102     deg_index[3] = 1;
00103     deg_index[5] = 2;
00104 
00105     // Set ghost flag to false
00106     hasghost = false;
00107     return MB_SUCCESS;
00108 }
00109 
00110 /************************************************************
00111  *     Interface Functions                                  *
00112  ************************************************************/
00113 
00114 ErrorCode NestedRefine::generate_mesh_hierarchy( int num_level, int* level_degrees,
00115                                                  std::vector< EntityHandle >& level_sets, bool optimize )
00116 {
00117     assert( num_level > 0 );
00118     nlevels = num_level;
00119 
00120     ErrorCode error;
00121     std::vector< moab::EntityHandle > hmsets( num_level );
00122 
00123     if( meshdim <= 2 )
00124     {
00125         for( int i = 0; i < num_level; i++ )
00126         {
00127             assert( ( level_degrees[i] == 2 ) || ( level_degrees[i] == 3 ) || ( level_degrees[i] == 5 ) );
00128             level_dsequence[i] = level_degrees[i];
00129         }
00130     }
00131     else
00132     {
00133         for( int i = 0; i < num_level; i++ )
00134         {
00135             assert( ( level_degrees[i] == 2 ) || ( level_degrees[i] == 3 ) );
00136             level_dsequence[i] = level_degrees[i];
00137         }
00138     }
00139 
00140     error = generate_hm( level_degrees, num_level, &hmsets[0], optimize );MB_CHK_ERR( error );
00141 
00142     // copy the entity handles
00143     level_sets.resize( num_level + 1 );
00144     level_sets[0] = _rset;
00145     for( int i = 0; i < num_level; i++ )
00146         level_sets[i + 1] = hmsets[i];
00147 
00148     return MB_SUCCESS;
00149 }
00150 
00151 ErrorCode NestedRefine::get_connectivity( EntityHandle ent, int level, std::vector< EntityHandle >& conn )
00152 {
00153     ErrorCode error;
00154     EntityType type = mbImpl->type_from_handle( ent );
00155     EntityHandle start_ent;
00156     if( !conn.empty() ) conn.clear();
00157     if( level > 0 )
00158     {
00159         if( type == MBEDGE )
00160         {
00161             conn.reserve( 2 );
00162             start_ent       = level_mesh[level - 1].start_edge;
00163             EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
00164             conn.push_back( level_mesh[level - 1].edge_conn[2 * offset] );
00165             conn.push_back( level_mesh[level - 1].edge_conn[2 * offset + 1] );
00166         }
00167         else if( type == MBTRI || type == MBQUAD )
00168         {
00169             int num_corners = ahf->lConnMap2D[type - 2].num_verts_in_face;
00170             conn.reserve( num_corners );
00171             start_ent       = level_mesh[level - 1].start_face;
00172             EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
00173 
00174             for( int i = 0; i < num_corners; i++ )
00175                 conn.push_back( level_mesh[level - 1].face_conn[num_corners * offset + i] );
00176         }
00177         else if( type == MBTET || type == MBHEX )
00178         {
00179             int index       = ahf->get_index_in_lmap( *_incells.begin() );
00180             int num_corners = ahf->lConnMap3D[index].num_verts_in_cell;
00181             conn.reserve( num_corners );
00182             start_ent       = level_mesh[level - 1].start_cell;
00183             EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
00184             for( int i = 0; i < num_corners; i++ )
00185                 conn.push_back( level_mesh[level - 1].cell_conn[num_corners * offset + i] );
00186         }
00187         else
00188             MB_SET_ERR( MB_FAILURE, "Requesting connectivity for an unsupported entity type" );
00189     }
00190     else
00191     {
00192         error = mbImpl->get_connectivity( &ent, 1, conn );MB_CHK_ERR( error );
00193     }
00194 
00195     return MB_SUCCESS;
00196 }
00197 
00198 ErrorCode NestedRefine::get_coordinates( EntityHandle* verts, int num_verts, int level, double* coords )
00199 {
00200     if( level > 0 )
00201     {
00202         EntityID vstart = ID_FROM_HANDLE( level_mesh[level - 1].start_vertex );
00203         for( int i = 0; i < num_verts; i++ )
00204         {
00205             const EntityHandle& vid = verts[i];
00206             EntityID offset         = ID_FROM_HANDLE( vid ) - vstart;
00207             coords[3 * i]           = level_mesh[level - 1].coordinates[0][offset];
00208             coords[3 * i + 1]       = level_mesh[level - 1].coordinates[1][offset];
00209             coords[3 * i + 2]       = level_mesh[level - 1].coordinates[2][offset];
00210         }
00211     }
00212     else
00213     {
00214         ErrorCode error;
00215         error = mbImpl->get_coords( verts, num_verts, coords );MB_CHK_ERR( error );
00216     }
00217 
00218     return MB_SUCCESS;
00219 }
00220 
00221 ErrorCode NestedRefine::get_adjacencies( const EntityHandle source_entity, const unsigned int target_dimension,
00222                                          std::vector< EntityHandle >& target_entities )
00223 
00224 {
00225     ErrorCode error;
00226     error = ahf->get_adjacencies( source_entity, target_dimension, target_entities );MB_CHK_ERR( error );
00227 
00228     return MB_SUCCESS;
00229 }
00230 
00231 ErrorCode NestedRefine::child_to_parent( EntityHandle child, int child_level, int parent_level, EntityHandle* parent )
00232 {
00233     assert( ( child_level > 0 ) && ( child_level > parent_level ) );
00234     EntityType type = mbImpl->type_from_handle( child );
00235     assert( type != MBVERTEX );
00236 
00237     int child_index;
00238     if( type == MBEDGE )
00239         child_index = child - level_mesh[child_level - 1].start_edge;
00240     else if( type == MBTRI || type == MBQUAD )
00241         child_index = child - level_mesh[child_level - 1].start_face;
00242     else if( type == MBTET || type == MBHEX )
00243         child_index = child - level_mesh[child_level - 1].start_cell;
00244     else
00245         MB_SET_ERR( MB_FAILURE, "Requesting parent for unsupported entity type" );
00246 
00247     int parent_index;
00248     int l = child_level - parent_level;
00249     for( int i = 0; i < l; i++ )
00250     {
00251         int d       = get_index_from_degree( level_dsequence[child_level - i - 1] );
00252         int nch     = refTemplates[type - 1][d].total_new_ents;
00253         child_index = child_index / nch;
00254     }
00255     parent_index = child_index;
00256 
00257     if( type == MBEDGE )
00258     {
00259         if( parent_level > 0 )
00260             *parent = level_mesh[parent_level - 1].start_edge + parent_index;
00261         else
00262             *parent = _inedges[parent_index];
00263     }
00264     else if( type == MBTRI || type == MBQUAD )
00265     {
00266         if( parent_level > 0 )
00267             *parent = level_mesh[parent_level - 1].start_face + parent_index;
00268         else
00269             *parent = _infaces[parent_index];
00270     }
00271     else if( type == MBTET || type == MBHEX )
00272     {
00273         if( parent_level > 0 )
00274             *parent = level_mesh[parent_level - 1].start_cell + parent_index;
00275         else
00276             *parent = _incells[parent_index];
00277     }
00278 
00279     return MB_SUCCESS;
00280 }
00281 
00282 ErrorCode NestedRefine::parent_to_child( EntityHandle parent, int parent_level, int child_level,
00283                                          std::vector< EntityHandle >& children )
00284 {
00285     assert( ( child_level > 0 ) && ( child_level > parent_level ) );
00286     EntityType type = mbImpl->type_from_handle( parent );
00287     assert( type != MBVERTEX );
00288 
00289     int parent_index;
00290     if( type == MBEDGE )
00291     {
00292         if( parent_level > 0 )
00293             parent_index = parent - level_mesh[parent_level - 1].start_edge;
00294         else
00295             parent_index = _inedges.index( parent );
00296     }
00297     else if( type == MBTRI || type == MBQUAD )
00298     {
00299         if( parent_level > 0 )
00300             parent_index = parent - level_mesh[parent_level - 1].start_face;
00301         else
00302             parent_index = _infaces.index( parent );
00303     }
00304     else if( type == MBTET || type == MBHEX )
00305     {
00306         if( parent_level > 0 )
00307             parent_index = parent - level_mesh[parent_level - 1].start_cell;
00308         else
00309             parent_index = _incells.index( parent );
00310     }
00311     else
00312         MB_SET_ERR( MB_FAILURE, "Requesting children for unsupported entity type" );
00313 
00314     int start, end;
00315     start = end = parent_index;
00316     for( int i = parent_level; i < child_level; i++ )
00317     {
00318         int d   = get_index_from_degree( level_dsequence[i] );
00319         int nch = refTemplates[type - 1][d].total_new_ents;
00320         start   = start * nch;
00321         end     = end * nch + nch - 1;
00322     }
00323 
00324     int num_child = end - start;
00325     children.reserve( num_child );
00326 
00327     for( int i = start; i <= end; i++ )
00328     {
00329         EntityHandle child;
00330         if( type == MBEDGE )
00331             child = level_mesh[child_level - 1].start_edge + i;
00332         else if( type == MBTRI || type == MBQUAD )
00333             child = level_mesh[child_level - 1].start_face + i;
00334         else if( type == MBTET || type == MBHEX )
00335             child = level_mesh[child_level - 1].start_cell + i;
00336 
00337         children.push_back( child );
00338     }
00339 
00340     return MB_SUCCESS;
00341 }
00342 
00343 ErrorCode NestedRefine::vertex_to_entities_up( EntityHandle vertex, int vert_level, int parent_level,
00344                                                std::vector< EntityHandle >& incident_entities )
00345 {
00346     assert( vert_level > parent_level );
00347     ErrorCode error;
00348 
00349     // Step 1: Get the incident entities at the current level
00350     std::vector< EntityHandle > inents;
00351     if( meshdim == 1 )
00352     {
00353         error = ahf->get_up_adjacencies_1d( vertex, inents );MB_CHK_ERR( error );
00354     }
00355     else if( meshdim == 2 )
00356     {
00357         error = ahf->get_up_adjacencies_vert_2d( vertex, inents );MB_CHK_ERR( error );
00358     }
00359     else if( meshdim == 3 )
00360     {
00361         error = ahf->get_up_adjacencies_vert_3d( vertex, inents );MB_CHK_ERR( error );
00362     }
00363 
00364     // Step 2: Loop over all the incident entities at the current level and gather their parents
00365     for( int i = 0; i < (int)inents.size(); i++ )
00366     {
00367         EntityHandle ent = inents[i];
00368         EntityHandle parent;
00369         error = child_to_parent( ent, vert_level, parent_level, &parent );MB_CHK_ERR( error );
00370         incident_entities.push_back( parent );
00371     }
00372 
00373     // Step 3: Sort and remove duplicates
00374     std::sort( incident_entities.begin(), incident_entities.end() );
00375     incident_entities.erase( std::unique( incident_entities.begin(), incident_entities.end() ),
00376                              incident_entities.end() );
00377 
00378     return MB_SUCCESS;
00379 }
00380 
00381 ErrorCode NestedRefine::vertex_to_entities_down( EntityHandle vertex, int vert_level, int child_level,
00382                                                  std::vector< EntityHandle >& incident_entities )
00383 {
00384     assert( vert_level < child_level );
00385     ErrorCode error;
00386 
00387     // Step 1: Get the incident entities at the current level
00388     std::vector< EntityHandle > inents;
00389     if( meshdim == 1 )
00390     {
00391         error = ahf->get_up_adjacencies_1d( vertex, inents );MB_CHK_ERR( error );
00392     }
00393     else if( meshdim == 2 )
00394     {
00395         error = ahf->get_up_adjacencies_vert_2d( vertex, inents );MB_CHK_ERR( error );
00396     }
00397     else if( meshdim == 3 )
00398     {
00399         error = ahf->get_up_adjacencies_vert_3d( vertex, inents );MB_CHK_ERR( error );
00400     }
00401 
00402     // Step 2: Loop over all the incident entities at the current level and gather their parents
00403     std::vector< EntityHandle > childs;
00404     for( int i = 0; i < (int)inents.size(); i++ )
00405     {
00406         childs.clear();
00407         EntityHandle ent = inents[i];
00408         error            = parent_to_child( ent, vert_level, child_level, childs );MB_CHK_ERR( error );
00409         for( int j = 0; j < (int)childs.size(); j++ )
00410             incident_entities.push_back( childs[j] );
00411     }
00412 
00413     return MB_SUCCESS;
00414 }
00415 
00416 ErrorCode NestedRefine::get_vertex_duplicates( EntityHandle vertex, int level, EntityHandle& dupvertex )
00417 {
00418     if( ( vertex - *_inverts.begin() ) > _inverts.size() )
00419         MB_SET_ERR( MB_FAILURE, "Requesting duplicates for non-coarse vertices" );
00420 
00421     dupvertex = level_mesh[level - 1].start_vertex + ( vertex - *_inverts.begin() );
00422 
00423     return MB_SUCCESS;
00424 }
00425 
00426 bool NestedRefine::is_entity_on_boundary( const EntityHandle& entity )
00427 {
00428     bool is_border  = false;
00429     EntityType type = mbImpl->type_from_handle( entity );
00430 
00431     if( type == MBVERTEX )
00432         is_border = is_vertex_on_boundary( entity );
00433     else if( type == MBEDGE )
00434         is_border = is_edge_on_boundary( entity );
00435     else if( type == MBTRI || type == MBQUAD )
00436         is_border = is_face_on_boundary( entity );
00437     else if( type == MBTET || type == MBHEX )
00438         is_border = is_cell_on_boundary( entity );
00439     else
00440         MB_SET_ERR( MB_FAILURE, "Requesting boundary information for unsupported entity type" );
00441 
00442     return is_border;
00443 }
00444 
00445 ErrorCode NestedRefine::exchange_ghosts( std::vector< EntityHandle >& lsets, int num_glayers )
00446 {
00447     ErrorCode error;
00448 
00449     if( hasghost ) return MB_SUCCESS;
00450 
00451     hasghost = true;
00452 #ifdef MOAB_HAVE_MPI
00453     error = pcomm->exchange_ghost_cells( meshdim, 0, num_glayers, 0, true, false );MB_CHK_ERR( error );
00454     {
00455         Range empty_range;
00456         error = pcomm->exchange_tags( GLOBAL_ID_TAG_NAME, empty_range );MB_CHK_ERR( error );
00457         // error = pcomm->assign_global_ids(lsets[i], 0, 1, false, true, false);MB_CHK_ERR(error);
00458     }
00459 #else
00460     MB_SET_ERR( MB_FAILURE, "Requesting ghost layers for a serial mesh" );
00461 #endif
00462 
00463     Range* lverts = new Range[lsets.size()];
00464     Range* lents  = new Range[lsets.size()];
00465     for( size_t i = 0; i < lsets.size(); i++ )
00466     {
00467         error = mbImpl->get_entities_by_dimension( lsets[i], meshdim, lents[i] );MB_CHK_ERR( error );
00468         error = mbImpl->get_connectivity( lents[i], lverts[i] );MB_CHK_ERR( error );
00469 
00470         for( int gl = 0; gl < num_glayers; gl++ )
00471         {
00472             error = mbImpl->get_adjacencies( lverts[i], meshdim, false, lents[i], Interface::UNION );MB_CHK_ERR( error );
00473             error = mbImpl->get_connectivity( lents[i], lverts[i] );MB_CHK_ERR( error );
00474         }
00475     }
00476     for( size_t i = 0; i < lsets.size(); i++ )
00477     {
00478         error = mbImpl->add_entities( lsets[i], lverts[i] );MB_CHK_ERR( error );
00479         error = mbImpl->add_entities( lsets[i], lents[i] );MB_CHK_ERR( error );
00480     }
00481 
00482     delete[] lverts;
00483     delete[] lents;
00484     return MB_SUCCESS;
00485 }
00486 
00487 ErrorCode NestedRefine::update_special_tags( int level, EntityHandle& lset )
00488 {
00489     assert( level > 0 && level < nlevels + 1 );
00490 
00491     ErrorCode error;
00492     std::vector< Tag > mtags( 3 );
00493 
00494     error = mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mtags[0] );MB_CHK_ERR( error );
00495     error = mbImpl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mtags[1] );MB_CHK_ERR( error );
00496     error = mbImpl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mtags[2] );MB_CHK_ERR( error );
00497 
00498     for( int i = 0; i < 3; i++ )
00499     {
00500         // Gather sets of a particular tag
00501         Range sets;
00502         error = mbImpl->get_entities_by_type_and_tag( _rset, MBENTITYSET, &mtags[i], NULL, 1, sets );MB_CHK_ERR( error );
00503 
00504         // Loop over all sets, gather entities in each set and add their children at all levels to
00505         // the set
00506         Range set_ents;
00507         Range::iterator set_it;
00508         std::vector< EntityHandle > childs;
00509 
00510         for( set_it = sets.begin(); set_it != sets.end(); ++set_it )
00511         {
00512             // Get the entities in the set, recursively
00513             set_ents.clear();
00514             childs.clear();
00515             error = mbImpl->get_entities_by_handle( *set_it, set_ents, true );MB_CHK_ERR( error );
00516 
00517             // Gather child entities at the input level
00518             for( Range::iterator sit = set_ents.begin(); sit != set_ents.end(); sit++ )
00519             {
00520                 EntityType type = mbImpl->type_from_handle( *sit );
00521                 if( type == MBVERTEX )
00522                 {
00523                     Range conn;
00524                     std::vector< EntityHandle > cents;
00525                     error = vertex_to_entities_down( *sit, 0, level, cents );MB_CHK_ERR( error );
00526                     error = mbImpl->get_connectivity( &cents[0], (int)cents.size(), conn, true );MB_CHK_ERR( error );
00527                     childs.insert( childs.end(), cents.begin(), cents.end() );
00528                 }
00529                 else
00530                 {
00531                     error = parent_to_child( *sit, 0, level, childs );MB_CHK_ERR( error );
00532                 }
00533 
00534                 std::sort( childs.begin(), childs.end() );
00535                 childs.erase( std::unique( childs.begin(), childs.end() ), childs.end() );
00536 
00537                 // Add child entities to tagged sets
00538                 error = mbImpl->add_entities( *set_it, &childs[0], childs.size() );MB_CHK_ERR( error );
00539             }
00540 
00541             // Remove the coarse entities
00542             error = mbImpl->remove_entities( *set_it, set_ents );MB_CHK_ERR( error );
00543 
00544             // Add
00545             error = mbImpl->add_entities( lset, &( *set_it ), 1 );MB_CHK_ERR( error );
00546         }
00547     }
00548     return MB_SUCCESS;
00549 }
00550 
00551 /***********************************************
00552  *  Basic functionalities: generate HM         *
00553  ***********************************************/
00554 
00555 ErrorCode NestedRefine::estimate_hm_storage( EntityHandle set, int level_degree, int cur_level, int hmest[4] )
00556 {
00557     ErrorCode error;
00558 
00559     // Obtain the size of input mesh.
00560     int nverts_prev, nedges_prev, nfaces_prev, ncells_prev;
00561     if( cur_level )
00562     {
00563         nverts_prev = level_mesh[cur_level - 1].num_verts;
00564         nedges_prev = level_mesh[cur_level - 1].num_edges;
00565         nfaces_prev = level_mesh[cur_level - 1].num_faces;
00566         ncells_prev = level_mesh[cur_level - 1].num_cells;
00567     }
00568     else
00569     {
00570         nverts_prev = _inverts.size();
00571         nedges_prev = _inedges.size();
00572         nfaces_prev = _infaces.size();
00573         ncells_prev = _incells.size();
00574     }
00575 
00576     // Estimate mesh size of current level mesh.
00577     int nedges = 0, nfaces = 0;
00578     error = count_subentities( set, cur_level - 1, &nedges, &nfaces );MB_CHK_ERR( error );
00579 
00580     int d      = get_index_from_degree( level_degree );
00581     int nverts = refTemplates[MBEDGE - 1][d].nv_edge * nedges;
00582     hmest[0]   = nverts_prev + nverts;
00583     hmest[1]   = nedges_prev * refTemplates[MBEDGE - 1][d].total_new_ents;
00584     hmest[2]   = 0;
00585     hmest[3]   = 0;
00586 
00587     int findex, cindex;
00588     if( nfaces_prev != 0 )
00589     {
00590         EntityHandle start_face;
00591         if( cur_level )
00592             start_face = level_mesh[cur_level - 1].start_face;
00593         else
00594             start_face = *_infaces.begin();
00595         findex   = mbImpl->type_from_handle( start_face ) - 1;
00596         hmest[2] = nfaces_prev * refTemplates[findex][d].total_new_ents;
00597 
00598         if( meshdim == 2 ) hmest[0] += refTemplates[findex][d].nv_face * nfaces_prev;
00599 
00600         if( meshdim == 3 ) hmest[1] += nfaces_prev * intFacEdg[findex - 1][d].nie;
00601     }
00602 
00603     if( ncells_prev != 0 )
00604     {
00605         cindex   = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
00606         hmest[3] = ncells_prev * refTemplates[cindex][d].total_new_ents;
00607 
00608         hmest[0] += refTemplates[cindex][d].nv_face * nfaces;
00609         hmest[0] += refTemplates[cindex][d].nv_cell * ncells_prev;
00610     }
00611 
00612     return MB_SUCCESS;
00613 }
00614 
00615 ErrorCode NestedRefine::create_hm_storage_single_level( EntityHandle* set, int cur_level, int estL[4] )
00616 {
00617     // Obtain chunks of memory for the current level. Add them to a particular meshset.
00618     EntityHandle set_handle;
00619     ErrorCode error = mbImpl->create_meshset( MESHSET_SET, set_handle );MB_CHK_SET_ERR( error, "Cannot create mesh for the current level" );
00620     *set = set_handle;
00621 
00622     ReadUtilIface* read_iface;
00623     error = mbImpl->query_interface( read_iface );MB_CHK_ERR( error );
00624 
00625     // Vertices
00626     error = read_iface->get_node_coords( 3, estL[0], 0, level_mesh[cur_level].start_vertex,
00627                                          level_mesh[cur_level].coordinates );MB_CHK_ERR( error );
00628     level_mesh[cur_level].num_verts = estL[0];
00629 
00630     Range newverts( level_mesh[cur_level].start_vertex, level_mesh[cur_level].start_vertex + estL[0] - 1 );
00631     error = mbImpl->add_entities( *set, newverts );MB_CHK_ERR( error );
00632     level_mesh[cur_level].verts = newverts;
00633 
00634     Tag gidtag;
00635     error = mbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, gidtag );MB_CHK_ERR( error );
00636     error = read_iface->assign_ids( gidtag, newverts, level_mesh[cur_level].start_vertex );MB_CHK_ERR( error );
00637 
00638     // Edges
00639     if( estL[1] )
00640     {
00641         error = read_iface->get_element_connect( estL[1], 2, MBEDGE, 0, level_mesh[cur_level].start_edge,
00642                                                  level_mesh[cur_level].edge_conn );MB_CHK_ERR( error );
00643         level_mesh[cur_level].num_edges = estL[1];
00644 
00645         Range newedges( level_mesh[cur_level].start_edge, level_mesh[cur_level].start_edge + estL[1] - 1 );
00646         error = mbImpl->add_entities( *set, newedges );MB_CHK_ERR( error );
00647         level_mesh[cur_level].edges = newedges;
00648     }
00649     else
00650         level_mesh[cur_level].num_edges = 0;
00651 
00652     // Faces
00653     if( estL[2] )
00654     {
00655         EntityType type = mbImpl->type_from_handle( *( _infaces.begin() ) );
00656         int nvpf        = ahf->lConnMap2D[type - 2].num_verts_in_face;
00657         error           = read_iface->get_element_connect( estL[2], nvpf, type, 0, level_mesh[cur_level].start_face,
00658                                                  level_mesh[cur_level].face_conn );MB_CHK_ERR( error );
00659         level_mesh[cur_level].num_faces = estL[2];
00660 
00661         Range newfaces( level_mesh[cur_level].start_face, level_mesh[cur_level].start_face + estL[2] - 1 );
00662         error = mbImpl->add_entities( *set, newfaces );MB_CHK_ERR( error );
00663         level_mesh[cur_level].faces = newfaces;
00664     }
00665     else
00666         level_mesh[cur_level].num_faces = 0;
00667 
00668     // Cells
00669     if( estL[3] )
00670     {
00671         EntityType type = mbImpl->type_from_handle( *( _incells.begin() ) );
00672         int index       = ahf->get_index_in_lmap( *_incells.begin() );
00673         int nvpc        = ahf->lConnMap3D[index].num_verts_in_cell;
00674         error           = read_iface->get_element_connect( estL[3], nvpc, type, 0, level_mesh[cur_level].start_cell,
00675                                                  level_mesh[cur_level].cell_conn );MB_CHK_ERR( error );
00676         level_mesh[cur_level].num_cells = estL[3];
00677 
00678         Range newcells( level_mesh[cur_level].start_cell, level_mesh[cur_level].start_cell + estL[3] - 1 );
00679         error = mbImpl->add_entities( *set, newcells );MB_CHK_ERR( error );
00680         level_mesh[cur_level].cells = newcells;
00681     }
00682     else
00683         level_mesh[cur_level].num_cells = 0;
00684 
00685     // Resize the ahf maps
00686     error = ahf->resize_hf_maps( level_mesh[cur_level].start_vertex, level_mesh[cur_level].num_verts,
00687                                  level_mesh[cur_level].start_edge, level_mesh[cur_level].num_edges,
00688                                  level_mesh[cur_level].start_face, level_mesh[cur_level].num_faces,
00689                                  level_mesh[cur_level].start_cell, level_mesh[cur_level].num_cells );MB_CHK_ERR( error );
00690 
00691     error = ahf->update_entity_ranges( *set );MB_CHK_ERR( error );
00692 
00693     // If the mesh type changes, then update the member variable in ahf to use the applicable
00694     // adjacency matrix
00695     MESHTYPE nwmesh = ahf->get_mesh_type( level_mesh[cur_level].num_verts, level_mesh[cur_level].num_edges,
00696                                           level_mesh[cur_level].num_faces, level_mesh[cur_level].num_cells );MB_CHK_ERR( error );
00697     if( ahf->thismeshtype != nwmesh ) ahf->thismeshtype = nwmesh;
00698 
00699     return MB_SUCCESS;
00700 }
00701 
00702 /**********************************
00703  *   Hierarchical Mesh Generation  *
00704  * *********************************/
00705 
00706 ErrorCode NestedRefine::generate_hm( int* level_degrees, int num_level, EntityHandle* hm_set, bool optimize )
00707 {
00708     ErrorCode error;
00709 
00710     Tag gidtag;
00711     error = mbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, gidtag );MB_CHK_ERR( error );
00712 
00713     nlevels = num_level;
00714 
00715     timeall.tm_total   = 0;
00716     timeall.tm_refine  = 0;
00717     timeall.tm_resolve = 0;
00718 
00719     for( int l = 0; l < num_level; l++ )
00720     {
00721         double tstart;
00722         tstart = tm->time_elapsed();
00723 
00724         // Estimate storage
00725         int hmest[4] = { 0, 0, 0, 0 };
00726         EntityHandle set;
00727         if( l )
00728             set = hm_set[l - 1];
00729         else
00730             set = _rset;
00731         error = estimate_hm_storage( set, level_degrees[l], l, hmest );MB_CHK_ERR( error );
00732 
00733         // Create arrays for storing the current level
00734         error = create_hm_storage_single_level( &hm_set[l], l, hmest );MB_CHK_ERR( error );
00735 
00736         // Copy the old vertices along with their coordinates
00737         error = copy_vertices_from_prev_level( l );MB_CHK_ERR( error );
00738 
00739         // Create the new entities and new vertices
00740         error = construct_hm_entities( l, level_degrees[l] );MB_CHK_ERR( error );
00741 
00742         timeall.tm_refine += tm->time_elapsed() - tstart;
00743 
00744         // Go into parallel communication
00745         if( !optimize )
00746         {
00747 #ifdef MOAB_HAVE_MPI
00748             if( pcomm && ( pcomm->size() > 1 ) )
00749             {
00750                 double tpstart = tm->time_elapsed();
00751                 error          = resolve_shared_ents_parmerge( l, hm_set[l] );MB_CHK_ERR( error );
00752                 timeall.tm_resolve += tm->time_elapsed() - tpstart;
00753             }
00754 #endif
00755         }
00756     }
00757 
00758     if( optimize )
00759     {
00760 #ifdef MOAB_HAVE_MPI
00761         if( pcomm && ( pcomm->size() > 1 ) )
00762         {
00763             double tpstart = tm->time_elapsed();
00764             error          = resolve_shared_ents_opt( hm_set, nlevels );MB_CHK_ERR( error );
00765             timeall.tm_resolve = tm->time_elapsed() - tpstart;
00766         }
00767 #endif
00768     }
00769     timeall.tm_total = timeall.tm_refine + timeall.tm_resolve;
00770 
00771     return MB_SUCCESS;
00772 }
00773 
00774 ErrorCode NestedRefine::construct_hm_entities( int cur_level, int deg )
00775 {
00776     ErrorCode error;
00777 
00778     // Generate mesh for current level by refining previous level.
00779     if( ahf->thismeshtype == CURVE )
00780     {
00781         error = construct_hm_1D( cur_level, deg );MB_CHK_ERR( error );
00782     }
00783     else if( ahf->thismeshtype == SURFACE || ahf->thismeshtype == SURFACE_MIXED )
00784     {
00785         error = construct_hm_2D( cur_level, deg );MB_CHK_ERR( error );
00786     }
00787     else
00788     {
00789         error = construct_hm_3D( cur_level, deg );MB_CHK_ERR( error );
00790     }
00791 
00792     return MB_SUCCESS;
00793 }
00794 
00795 ErrorCode NestedRefine::construct_hm_1D( int cur_level, int deg )
00796 {
00797     ErrorCode error;
00798     int nverts_prev, nents_prev;
00799     if( cur_level )
00800     {
00801         nverts_prev = level_mesh[cur_level - 1].num_verts;
00802         nents_prev  = level_mesh[cur_level - 1].num_edges;
00803     }
00804     else
00805     {
00806         nverts_prev = _inverts.size();
00807         nents_prev  = _inedges.size();
00808     }
00809 
00810     int d      = get_index_from_degree( deg );
00811     int vtotal = 2 + refTemplates[0][d].total_new_verts;
00812     std::vector< EntityHandle > vbuffer( vtotal );
00813 
00814     std::vector< EntityHandle > conn;
00815     int count_nents = 0;
00816     int count_verts = nverts_prev;
00817 
00818     // Step 1: Create the subentities via refinement of the previous mesh
00819     for( int eid = 0; eid < nents_prev; eid++ )
00820     {
00821         conn.clear();
00822 
00823         // EntityHandle of the working edge
00824         EntityHandle edge;
00825         if( cur_level )
00826             edge = level_mesh[cur_level - 1].start_edge + eid;
00827         else
00828             edge = _inedges[eid];  // Makes the assumption initial mesh is contiguous in memory
00829 
00830         error = get_connectivity( edge, cur_level, conn );MB_CHK_ERR( error );
00831 
00832         // Add the vertex handles to vbuffer for the current level for the working edge
00833 
00834         // Since the old vertices are copied first, their local indices do not change as new levels
00835         // are added. Clearly the local indices of the new vertices introduced in the current level
00836         // is still the same when the old vertices are copied. Thus, there is no need to explicitly
00837         // store another map between the old and duplicates in the subsequent levels. The second
00838         // part in the following sum is the local index in the previous level.
00839 
00840         // Add the corners to the vbuffer first.
00841 
00842         for( int i = 0; i < (int)conn.size(); i++ )
00843         {
00844             if( cur_level )
00845                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
00846             else
00847                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
00848         }
00849 
00850         // Adding rest of the entityhandles to working buffer for vertices.
00851         int num_new_verts = refTemplates[0][d].total_new_verts;
00852 
00853         for( int i = 0; i < num_new_verts; i++ )
00854         {
00855             vbuffer[i + 2] = level_mesh[cur_level].start_vertex + count_verts;
00856             count_verts += 1;
00857         }
00858 
00859         // Use the template to obtain the subentities
00860         int id1, id2;
00861         int etotal = refTemplates[0][d].total_new_ents;
00862         std::vector< EntityHandle > ent_buffer( etotal );
00863 
00864         for( int i = 0; i < etotal; i++ )
00865         {
00866             id1                                                      = refTemplates[0][d].ents_conn[i][0];
00867             id2                                                      = refTemplates[0][d].ents_conn[i][1];
00868             level_mesh[cur_level].edge_conn[2 * ( count_nents )]     = vbuffer[id1];
00869             level_mesh[cur_level].edge_conn[2 * ( count_nents ) + 1] = vbuffer[id2];
00870             ent_buffer[i]                                            = level_mesh[cur_level].start_edge + count_nents;
00871             count_nents += 1;
00872         };
00873 
00874         error = update_local_ahf( deg, MBEDGE, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
00875 
00876         // Compute the coordinates of the new vertices: Linear interpolation
00877         int idx;
00878         double xi;
00879         for( int i = 0; i < num_new_verts; i++ )
00880         {
00881             xi  = refTemplates[0][d].vert_nat_coord[i][0];
00882             idx = vbuffer[i + 2] - level_mesh[cur_level].start_vertex;  // index of new vertex in current level
00883             if( cur_level )
00884             {
00885                 id1 = conn[0] - level_mesh[cur_level - 1].start_vertex;  // index of old end vertices in current level
00886                 id2 = conn[1] - level_mesh[cur_level - 1].start_vertex;
00887             }
00888             else
00889             {
00890                 id1 = _inverts.index( conn[0] );
00891                 id2 = _inverts.index( conn[1] );
00892             }
00893 
00894             level_mesh[cur_level].coordinates[0][idx] =
00895                 ( 1 - xi ) * level_mesh[cur_level].coordinates[0][id1] + xi * level_mesh[cur_level].coordinates[0][id2];
00896             level_mesh[cur_level].coordinates[1][idx] =
00897                 ( 1 - xi ) * level_mesh[cur_level].coordinates[1][id1] + xi * level_mesh[cur_level].coordinates[1][id2];
00898             level_mesh[cur_level].coordinates[2][idx] =
00899                 ( 1 - xi ) * level_mesh[cur_level].coordinates[2][id1] + xi * level_mesh[cur_level].coordinates[2][id2];
00900         }
00901     }
00902 
00903     error = update_global_ahf( MBEDGE, cur_level, deg );MB_CHK_ERR( error );
00904 
00905     return MB_SUCCESS;
00906 }
00907 
00908 ErrorCode NestedRefine::construct_hm_1D( int cur_level, int deg, EntityType type,
00909                                          std::vector< EntityHandle >& trackverts )
00910 {
00911     ErrorCode error;
00912 
00913     int nedges_prev;
00914     if( cur_level )
00915         nedges_prev = level_mesh[cur_level - 1].num_edges;
00916     else
00917         nedges_prev = _inedges.size();
00918 
00919     int d      = get_index_from_degree( deg );
00920     int nve    = refTemplates[0][d].nv_edge;
00921     int vtotal = 2 + refTemplates[0][d].total_new_verts;
00922     int etotal = refTemplates[0][d].total_new_ents;
00923     int ne = 0, dim = 0, index = 0;
00924     if( type == MBTRI || type == MBQUAD )
00925     {
00926         index = type - 2;
00927         ne    = ahf->lConnMap2D[index].num_verts_in_face;
00928         dim   = 2;
00929     }
00930     else if( type == MBTET || type == MBHEX )
00931     {
00932         index = ahf->get_index_in_lmap( *( _incells.begin() ) );
00933         ne    = ahf->lConnMap3D[index].num_edges_in_cell;
00934         dim   = 3;
00935     }
00936 
00937     std::vector< EntityHandle > vbuffer( vtotal );
00938     std::vector< EntityHandle > ent_buffer( etotal );
00939 
00940     std::vector< EntityHandle > adjents, econn, fconn;
00941     std::vector< int > leids;
00942     int count_nents = 0;
00943 
00944     // Loop over all the edges and gather the vertices to be used for refinement
00945     for( int eid = 0; eid < nedges_prev; eid++ )
00946     {
00947         adjents.clear();
00948         leids.clear();
00949         econn.clear();
00950         fconn.clear();
00951         for( int i = 0; i < vtotal; i++ )
00952             vbuffer[i] = 0;
00953         for( int i = 0; i < etotal; i++ )
00954             ent_buffer[i] = 0;
00955 
00956         EntityHandle edge;
00957         if( cur_level )
00958             edge = level_mesh[cur_level - 1].start_edge + eid;
00959         else
00960             edge = _inedges[eid];
00961 
00962         error = get_connectivity( edge, cur_level, econn );MB_CHK_ERR( error );
00963 
00964         for( int i = 0; i < (int)econn.size(); i++ )
00965         {
00966             if( cur_level )
00967                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( econn[i] - level_mesh[cur_level - 1].start_vertex );
00968             else
00969                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( econn[i] - *_inverts.begin() );
00970         }
00971 
00972         int fid = -1, lid = -1, idx1 = -1, idx2 = -1;
00973 
00974         if( dim == 2 )
00975         {
00976             error = ahf->get_up_adjacencies_2d( edge, adjents, &leids );MB_CHK_ERR( error );
00977             if( cur_level )
00978                 fid = adjents[0] - level_mesh[cur_level - 1].start_face;
00979             else
00980                 fid = _infaces.index( adjents[0] );
00981 
00982             lid  = leids[0];
00983             idx1 = lid;
00984             idx2 = ahf->lConnMap2D[index].next[lid];
00985         }
00986         else if( dim == 3 )
00987         {
00988             error = ahf->get_up_adjacencies_edg_3d( edge, adjents, &leids );MB_CHK_ERR( error );
00989             if( cur_level )
00990                 fid = adjents[0] - level_mesh[cur_level - 1].start_cell;
00991             else
00992                 fid = _incells.index( adjents[0] );
00993 
00994             lid  = leids[0];
00995             idx1 = ahf->lConnMap3D[index].e2v[lid][0];
00996             idx2 = ahf->lConnMap3D[index].e2v[lid][1];
00997         }
00998 
00999         error = get_connectivity( adjents[0], cur_level, fconn );MB_CHK_ERR( error );
01000 
01001         bool orient = false;
01002         if( ( fconn[idx1] == econn[0] ) && ( fconn[idx2] == econn[1] ) ) orient = true;
01003 
01004         if( orient )
01005         {
01006             for( int j = 0; j < nve; j++ )
01007                 vbuffer[j + 2] = trackverts[fid * ne * nve + nve * lid + j];
01008         }
01009         else
01010         {
01011             for( int j = 0; j < nve; j++ )
01012                 vbuffer[( nve - j - 1 ) + 2] = trackverts[fid * ne * nve + nve * lid + j];
01013         }
01014 
01015         // Use the template to obtain the subentities
01016         int id1, id2;
01017 
01018         for( int i = 0; i < etotal; i++ )
01019         {
01020             id1                                                      = refTemplates[0][d].ents_conn[i][0];
01021             id2                                                      = refTemplates[0][d].ents_conn[i][1];
01022             level_mesh[cur_level].edge_conn[2 * ( count_nents )]     = vbuffer[id1];
01023             level_mesh[cur_level].edge_conn[2 * ( count_nents ) + 1] = vbuffer[id2];
01024             ent_buffer[i]                                            = level_mesh[cur_level].start_edge + count_nents;
01025 
01026             count_nents += 1;
01027         };
01028 
01029         error = update_local_ahf( deg, MBEDGE, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
01030     }
01031 
01032     error = update_global_ahf_1D_sub( cur_level, deg );MB_CHK_ERR( error );
01033 
01034     return MB_SUCCESS;
01035 }
01036 
01037 ErrorCode NestedRefine::construct_hm_2D( int cur_level, int deg )
01038 {
01039     ErrorCode error;
01040     int nverts_prev, nents_prev;
01041     if( cur_level )
01042     {
01043         nverts_prev = level_mesh[cur_level - 1].num_verts;
01044         nents_prev  = level_mesh[cur_level - 1].num_faces;
01045     }
01046     else
01047     {
01048         nverts_prev = _inverts.size();
01049         nents_prev  = _infaces.size();
01050     }
01051 
01052     // Create some book-keeping arrays over the old mesh to avoid introducing duplicate vertices and
01053     // calculating vertices more than once.
01054     EntityType ftype = mbImpl->type_from_handle( *_infaces.begin() );
01055     int nepf         = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
01056     int findex       = ftype - 1;
01057 
01058     int d      = get_index_from_degree( deg );
01059     int tnv    = refTemplates[findex][d].total_new_verts;
01060     int vtotal = nepf + tnv;
01061     std::vector< EntityHandle > vbuffer( vtotal );
01062     int etotal = refTemplates[findex][d].total_new_ents;
01063     std::vector< EntityHandle > ent_buffer( etotal );
01064 
01065     int nve = refTemplates[findex][d].nv_edge;
01066     std::vector< EntityHandle > trackvertsF( nents_prev * nepf * nve, 0 );
01067     int cur_nverts = level_mesh[cur_level].num_verts;
01068     std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
01069 
01070     int count_nverts = nverts_prev;
01071     int count_nents  = 0;
01072     std::vector< EntityHandle > conn, cur_conn;
01073 
01074     // Step 1: Create the subentities via refinement of the previous mesh
01075     for( int fid = 0; fid < nents_prev; fid++ )
01076     {
01077         conn.clear();
01078         cur_conn.clear();
01079         for( int i = 0; i < vtotal; i++ )
01080             vbuffer[i] = 0;
01081         for( int i = 0; i < etotal; i++ )
01082             ent_buffer[i] = 0;
01083 
01084         // EntityHandle of the working face
01085         EntityHandle face;
01086         if( cur_level )
01087             face = level_mesh[cur_level - 1].start_face + fid;
01088         else
01089             face = _infaces[fid];
01090 
01091         error = get_connectivity( face, cur_level, conn );MB_CHK_ERR( error );
01092 
01093         // Step 1: Add vertices from the current level for the working face that will be used for
01094         // subdivision.
01095         // Add the corners to vbuffer
01096         for( int i = 0; i < (int)conn.size(); i++ )
01097         {
01098             if( cur_level )
01099                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
01100             else
01101                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
01102 
01103             cur_conn.push_back( vbuffer[i] );
01104         }
01105 
01106         // Gather vertices already added to tracking array due to refinement of the sibling faces
01107 
01108         for( int i = 0; i < nepf; i++ )
01109         {
01110             for( int j = 0; j < nve; j++ )
01111             {
01112                 int id      = refTemplates[findex][d].vert_on_edges[i][j];
01113                 vbuffer[id] = trackvertsF[fid * nve * nepf + nve * i + j];
01114             }
01115         }
01116 
01117         // Add the remaining vertex handles to vbuffer for the current level for the working face
01118         for( int i = 0; i < tnv; i++ )
01119         {
01120             if( !vbuffer[i + nepf] )
01121             {
01122                 vbuffer[i + nepf] = level_mesh[cur_level].start_vertex + count_nverts;
01123                 count_nverts += 1;
01124             }
01125         }
01126 
01127         // Step 2: Create the subentities using the template and the vbuffer
01128         int idx;
01129         for( int i = 0; i < etotal; i++ )
01130         {
01131             for( int k = 0; k < nepf; k++ )
01132             {
01133                 idx                                                     = refTemplates[findex][d].ents_conn[i][k];
01134                 level_mesh[cur_level].face_conn[nepf * count_nents + k] = vbuffer[idx];
01135             }
01136             ent_buffer[i] = level_mesh[cur_level].start_face + count_nents;
01137             count_nents += 1;
01138         }
01139 
01140         // Step 3: Update the local AHF maps
01141         error = update_local_ahf( deg, ftype, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
01142 
01143         // Step 4: Add the new vertices to the tracking array
01144         int id;
01145 
01146         for( int i = 0; i < nepf; i++ )
01147         {
01148             // Add the vertices to trackvertsF for fid
01149             for( int j = 0; j < nve; j++ )
01150             {
01151                 id                                          = refTemplates[findex][d].vert_on_edges[i][j];
01152                 trackvertsF[fid * nepf * nve + nve * i + j] = vbuffer[id];
01153             }
01154 
01155             std::vector< EntityHandle > sibfids;
01156             std::vector< int > sibleids;
01157             std::vector< int > siborient;
01158 
01159             // Add the vertices to trackvertsF for siblings of fid, if any.
01160             error = ahf->get_up_adjacencies_2d( face, i, false, sibfids, &sibleids, &siborient );MB_CHK_ERR( error );
01161 
01162             if( !sibfids.size() ) continue;
01163 
01164             for( int s = 0; s < (int)sibfids.size(); s++ )
01165             {
01166                 int sibid;
01167                 if( cur_level )
01168                     sibid = sibfids[s] - level_mesh[cur_level - 1].start_face;
01169                 else
01170                     sibid = sibfids[s] - *_infaces.begin();
01171 
01172                 if( siborient[s] )  // Same half-edge direction as the current half-edge
01173                 {
01174                     for( int j = 0; j < nve; j++ )
01175                     {
01176                         id = refTemplates[findex][d].vert_on_edges[i][j];
01177                         trackvertsF[sibid * nepf * nve + nve * sibleids[s] + j] = vbuffer[id];
01178                     }
01179                 }
01180                 else
01181                 {
01182                     for( int j = 0; j < nve; j++ )
01183                     {
01184                         id = refTemplates[findex][d].vert_on_edges[i][nve - j - 1];
01185                         trackvertsF[sibid * nepf * nve + nve * sibleids[s] + j] = vbuffer[id];
01186                     }
01187                 }
01188             }
01189         }
01190 
01191         // Step 5: Compute the coordinates of the new vertices, avoids computing more than once via
01192         // the flag_verts array.
01193         std::vector< double > corner_coords( nepf * 3 );
01194         error = get_coordinates( &cur_conn[0], nepf, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
01195 
01196         error = compute_coordinates( cur_level, deg, ftype, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
01197                                      nverts_prev );MB_CHK_ERR( error );
01198     }
01199 
01200     // Step 6: Update the global maps
01201     error = update_global_ahf( ftype, cur_level, deg );MB_CHK_ERR( error );
01202 
01203     // Step 7: If edges exists, refine them.
01204     if( !_inedges.empty() )
01205     {
01206         error = construct_hm_1D( cur_level, deg, ftype, trackvertsF );MB_CHK_ERR( error );
01207     }
01208 
01209     return MB_SUCCESS;
01210 }
01211 
01212 ErrorCode NestedRefine::construct_hm_2D( int cur_level, int deg, EntityType type,
01213                                          std::vector< EntityHandle >& trackvertsE,
01214                                          std::vector< EntityHandle >& trackvertsF )
01215 {
01216     ErrorCode error;
01217 
01218     EntityType ftype = MBTRI;
01219     if( type == MBHEX ) ftype = MBQUAD;
01220 
01221     int d      = get_index_from_degree( deg );
01222     int findex = ftype - 1;
01223     int cidx   = ahf->get_index_in_lmap( *( _incells.begin() ) );
01224 
01225     int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
01226     int nepc = ahf->lConnMap3D[cidx].num_edges_in_cell;
01227     int nfpc = ahf->lConnMap3D[cidx].num_faces_in_cell;
01228 
01229     int tnv    = refTemplates[findex][d].total_new_verts;
01230     int nve    = refTemplates[findex][d].nv_edge;
01231     int nvf    = refTemplates[findex][d].nv_face;
01232     int vtotal = nepf + tnv;
01233     int etotal = refTemplates[findex][d].total_new_ents;
01234 
01235     std::vector< EntityHandle > vbuffer( vtotal );
01236     std::vector< EntityHandle > ent_buffer( etotal );
01237 
01238     std::vector< EntityHandle > adjents, fconn, cconn;
01239     std::vector< int > leids;
01240     int count_nents = 0;
01241 
01242     int nents_prev, ecount;
01243     if( cur_level )
01244     {
01245         nents_prev = level_mesh[cur_level - 1].num_faces;
01246         ecount     = level_mesh[cur_level - 1].num_edges * refTemplates[MBEDGE - 1][d].total_new_ents;
01247         ;
01248     }
01249     else
01250     {
01251         nents_prev = _infaces.size();
01252         ecount     = _inedges.size() * refTemplates[MBEDGE - 1][d].total_new_ents;
01253         ;
01254     }
01255 
01256     // Step 1: Create the subentities via refinement of the previous mesh
01257     for( int it = 0; it < nents_prev; it++ )
01258     {
01259         fconn.clear();
01260         cconn.clear();
01261         adjents.clear();
01262         leids.clear();
01263         for( int i = 0; i < vtotal; i++ )
01264             vbuffer[i] = 0;
01265         for( int i = 0; i < etotal; i++ )
01266             ent_buffer[i] = 0;
01267 
01268         // EntityHandle of the working face
01269         EntityHandle face;
01270         if( cur_level )
01271             face = level_mesh[cur_level - 1].start_face + it;
01272         else
01273             face = _infaces[it];
01274 
01275         error = get_connectivity( face, cur_level, fconn );MB_CHK_ERR( error );
01276 
01277         // Add the new handles for old connectivity in the buffer
01278         for( int i = 0; i < (int)fconn.size(); i++ )
01279         {
01280             if( cur_level )
01281                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( fconn[i] - level_mesh[cur_level - 1].start_vertex );
01282             else
01283                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( fconn[i] - *_inverts.begin() );
01284         }
01285 
01286         // Add handles for vertices on edges and faces from the already refined cell
01287         int fid, lid;
01288         error = ahf->get_up_adjacencies_face_3d( face, adjents, &leids );MB_CHK_ERR( error );
01289 
01290         if( cur_level )
01291             fid = adjents[0] - level_mesh[cur_level - 1].start_cell;
01292         else
01293             fid = _incells.index( adjents[0] );
01294 
01295         lid = leids[0];
01296 
01297         error = get_connectivity( adjents[0], cur_level, cconn );MB_CHK_ERR( error );
01298 
01299         // Find the orientation w.r.t the half-face and then add vertices properly.
01300         std::vector< EntityHandle > fac_conn( nepf );
01301         std::vector< EntityHandle > lfac_conn( nepf );
01302         for( int j = 0; j < nepf; j++ )
01303         {
01304             fac_conn[j]  = fconn[j];
01305             int id       = ahf->lConnMap3D[cidx].hf2v[lid][j];
01306             lfac_conn[j] = cconn[id];
01307         }
01308 
01309         std::vector< int > le_idx, indices;
01310 
01311         error = reorder_indices( deg, &fac_conn[0], &lfac_conn[0], nepf, le_idx, indices );MB_CHK_ERR( error );
01312 
01313         // Add the existing vertices on edges of the already refined cell to the vbuffer
01314         for( int j = 0; j < nepf; j++ )
01315         {
01316             int id  = le_idx[j];                              // Corresponding local edge
01317             int idx = ahf->lConnMap3D[cidx].f2leid[lid][id];  // Local edge in the cell
01318 
01319             // Get the orientation of the local edge of the face wrt the corresponding local edge in
01320             // the cell
01321             bool eorient = false;
01322             int fnext    = ahf->lConnMap2D[ftype - 2].next[j];
01323             int idx1     = ahf->lConnMap3D[cidx].e2v[idx][0];
01324             int idx2     = ahf->lConnMap3D[cidx].e2v[idx][1];
01325             if( ( fconn[j] == cconn[idx1] ) && ( fconn[fnext] == cconn[idx2] ) ) eorient = true;
01326 
01327             if( eorient )
01328             {
01329                 for( int k = 0; k < nve; k++ )
01330                 {
01331                     int ind      = refTemplates[findex][d].vert_on_edges[j][k];
01332                     vbuffer[ind] = trackvertsE[fid * nepc * nve + nve * idx + k];
01333                 }
01334             }
01335             else
01336             {
01337                 for( int k = 0; k < nve; k++ )
01338                 {
01339                     int ind      = refTemplates[findex][d].vert_on_edges[j][nve - k - 1];
01340                     vbuffer[ind] = trackvertsE[fid * nepc * nve + nve * idx + k];
01341                 }
01342             }
01343         }
01344 
01345         // Add the existing vertices on the face of the refine cell to vbuffer
01346         if( nvf )
01347         {
01348             for( int k = 0; k < nvf; k++ )
01349             {
01350                 int ind      = refTemplates[findex][d].vert_on_faces[0][k];
01351                 vbuffer[ind] = trackvertsF[fid * nfpc * nvf + nvf * lid + indices[k] - 1];
01352             }
01353         }
01354 
01355         // Create the subentities using the template and the vbuffer
01356         for( int i = 0; i < etotal; i++ )
01357         {
01358             for( int k = 0; k < nepf; k++ )
01359             {
01360                 int idx                                                 = refTemplates[findex][d].ents_conn[i][k];
01361                 level_mesh[cur_level].face_conn[nepf * count_nents + k] = vbuffer[idx];
01362             }
01363             ent_buffer[i] = level_mesh[cur_level].start_face + count_nents;
01364             count_nents += 1;
01365         }
01366 
01367         error = update_local_ahf( deg, ftype, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
01368 
01369         // Create the interior edges
01370         int id1, id2;
01371 
01372         int ne = intFacEdg[ftype - 2][d].nie;
01373         for( int i = 0; i < ne; i++ )
01374         {
01375             id1                                                 = intFacEdg[ftype - 2][d].ieconn[i][0];
01376             id2                                                 = intFacEdg[ftype - 2][d].ieconn[i][1];
01377             level_mesh[cur_level].edge_conn[2 * ( ecount )]     = vbuffer[id1];
01378             level_mesh[cur_level].edge_conn[2 * ( ecount ) + 1] = vbuffer[id2];
01379             ecount += 1;
01380         }
01381     }
01382 
01383     // Step 6: Update the global maps
01384     error = update_global_ahf_2D_sub( cur_level, deg );MB_CHK_ERR( error );
01385 
01386     // Step 7: Update the hf-maps for the edges
01387     error = update_ahf_1D( cur_level );MB_CHK_ERR( error );
01388 
01389     return MB_SUCCESS;
01390 }
01391 
01392 ErrorCode NestedRefine::construct_hm_3D( int cur_level, int deg )
01393 {
01394     ErrorCode error;
01395     EntityType type = mbImpl->type_from_handle( *( _incells.begin() ) );
01396     if( type == MBTET )
01397     {
01398         error = subdivide_tets( cur_level, deg );MB_CHK_ERR( error );
01399     }
01400     else
01401     {
01402         error = subdivide_cells( type, cur_level, deg );MB_CHK_ERR( error );
01403     }
01404 
01405     return MB_SUCCESS;
01406 }
01407 
01408 ErrorCode NestedRefine::subdivide_cells( EntityType type, int cur_level, int deg )
01409 {
01410     ErrorCode error;
01411     int nverts_prev, nents_prev;
01412     if( cur_level )
01413     {
01414         nverts_prev = level_mesh[cur_level - 1].num_verts;
01415         nents_prev  = level_mesh[cur_level - 1].num_cells;
01416     }
01417     else
01418     {
01419         nverts_prev = _inverts.size();
01420         nents_prev  = _incells.size();
01421     }
01422 
01423     // Create some book-keeping arrays over the parent mesh to avoid introducing duplicate vertices
01424     int cindex  = type - 1;
01425     int d       = get_index_from_degree( deg );
01426     int ne      = refTemplates[cindex][d].nv_edge;
01427     int nvf     = refTemplates[cindex][d].nv_face;
01428     int nvtotal = refTemplates[cindex][d].total_new_verts;
01429 
01430     int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
01431     int nvpc  = ahf->lConnMap3D[index].num_verts_in_cell;
01432     int nepc  = ahf->lConnMap3D[index].num_edges_in_cell;
01433     int nfpc  = ahf->lConnMap3D[index].num_faces_in_cell;
01434 
01435     int vtotal = nvpc + nvtotal;
01436     std::vector< EntityHandle > vbuffer( vtotal );
01437 
01438     std::vector< EntityHandle > trackvertsC_edg( nepc * ne * nents_prev, 0 );
01439     std::vector< EntityHandle > trackvertsC_face( nfpc * nvf * nents_prev, 0 );
01440 
01441     int cur_nverts = level_mesh[cur_level].num_verts;
01442     std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
01443 
01444     int count_nverts = nverts_prev;
01445     int count_ents   = 0;
01446     std::vector< EntityHandle > conn, cur_conn;
01447 
01448     // Step 1: Create the subentities via refinement of the previous mesh
01449     for( int cid = 0; cid < nents_prev; cid++ )
01450     {
01451         conn.clear();
01452         cur_conn.clear();
01453         for( int i = 0; i < vtotal; i++ )
01454             vbuffer[i] = 0;
01455 
01456         // EntityHandle of the working cell
01457         EntityHandle cell;
01458         if( cur_level )
01459             cell = level_mesh[cur_level - 1].start_cell + cid;
01460         else
01461             cell = _incells[cid];
01462 
01463         error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
01464 
01465         // Step 1: Add vertices from the current level for the working face that will be used for
01466         // subdivision.
01467         // Add the corners to vbuffer
01468         for( int i = 0; i < (int)conn.size(); i++ )
01469         {
01470             if( cur_level )
01471                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
01472             else
01473                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
01474 
01475             cur_conn.push_back( vbuffer[i] );
01476         }
01477 
01478         // Gather vertices already added to tracking array due to refinement of the sibling cells
01479         for( int i = 0; i < nepc; i++ )
01480         {
01481             for( int j = 0; j < ne; j++ )
01482             {
01483                 int idx      = refTemplates[cindex][d].vert_on_edges[i][j];
01484                 vbuffer[idx] = trackvertsC_edg[cid * nepc * ne + ne * i + j];
01485             }
01486         }
01487 
01488         // Add remaining new vertex handles
01489         for( int i = 0; i < nfpc; i++ )
01490         {
01491             for( int j = 0; j < nvf; j++ )
01492             {
01493                 int idx      = refTemplates[cindex][d].vert_on_faces[i][j];
01494                 vbuffer[idx] = trackvertsC_face[cid * nfpc * nvf + nvf * i + j];
01495             }
01496         }
01497 
01498         // Add the remaining vertex handles to vbuffer for the current level for the working cell
01499         for( int i = 0; i < nvtotal; i++ )
01500         {
01501             if( !vbuffer[i + nvpc] )
01502             {
01503                 vbuffer[i + nvpc] = level_mesh[cur_level].start_vertex + count_nverts;
01504                 count_nverts += 1;
01505             }
01506         }
01507 
01508         // Step 2: Use the template to obtain the subentities. The coordinates and local ahf maps
01509         // are also constructed. Connectivity of the children
01510         int etotal = refTemplates[type - 1][d].total_new_ents;
01511         std::vector< EntityHandle > ent_buffer( etotal );
01512 
01513         for( int i = 0; i < etotal; i++ )
01514         {
01515             for( int k = 0; k < nvpc; k++ )
01516             {
01517                 int idx                                                = refTemplates[type - 1][d].ents_conn[i][k];
01518                 level_mesh[cur_level].cell_conn[nvpc * count_ents + k] = vbuffer[idx];
01519             }
01520             ent_buffer[i] = level_mesh[cur_level].start_cell + count_ents;
01521             count_ents += 1;
01522         }
01523 
01524         // Step 3: Update local ahf maps
01525         error = update_local_ahf( deg, type, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
01526 
01527         // Step 4: Update tracking information
01528         error = update_tracking_verts( cell, cur_level, deg, trackvertsC_edg, trackvertsC_face, &vbuffer[0] );MB_CHK_ERR( error );
01529 
01530         // Step 5: Coordinates of the new vertices
01531         std::vector< double > corner_coords( nvpc * 3 );
01532         error = get_coordinates( &cur_conn[0], nvpc, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
01533 
01534         error = compute_coordinates( cur_level, deg, type, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
01535                                      nverts_prev );MB_CHK_ERR( error );
01536     }
01537 
01538     // error = ahf->print_tags(3);
01539 
01540     // Step 6: Update the global maps
01541     error = update_global_ahf( type, cur_level, deg );MB_CHK_ERR( error );
01542 
01543     // Step 7: If edges exists, refine them as well.
01544     if( level_mesh[cur_level].num_edges != 0 )
01545     {
01546         error = construct_hm_1D( cur_level, deg, type, trackvertsC_edg );MB_CHK_ERR( error );
01547     }
01548 
01549     // Step 8: If faces exists, refine them as well.
01550     if( !_infaces.empty() )
01551     {
01552         error = construct_hm_2D( cur_level, deg, type, trackvertsC_edg, trackvertsC_face );MB_CHK_ERR( error );
01553     }
01554 
01555     // error = ahf->print_tags(3);
01556 
01557     return MB_SUCCESS;
01558 }
01559 
01560 ErrorCode NestedRefine::subdivide_tets( int cur_level, int deg )
01561 {
01562     ErrorCode error;
01563     int nverts_prev, nents_prev;
01564     if( cur_level )
01565     {
01566         nverts_prev = level_mesh[cur_level - 1].num_verts;
01567         nents_prev  = level_mesh[cur_level - 1].num_cells;
01568     }
01569     else
01570     {
01571         nverts_prev = _inverts.size();
01572         nents_prev  = _incells.size();
01573     }
01574 
01575     EntityType type = MBTET;
01576     int cindex      = type - 1;
01577     int d           = get_index_from_degree( deg );
01578     int ne          = refTemplates[cindex][d].nv_edge;
01579     int nvf         = refTemplates[cindex][d].nv_face;
01580     int nvtotal     = refTemplates[cindex][d].total_new_verts;
01581 
01582     int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
01583     int nvpc  = ahf->lConnMap3D[index].num_verts_in_cell;
01584     int nepc  = ahf->lConnMap3D[index].num_edges_in_cell;
01585     int nfpc  = ahf->lConnMap3D[index].num_faces_in_cell;
01586 
01587     // Create vertex buffer
01588     int vtotal = nvpc + nvtotal;
01589     std::vector< EntityHandle > vbuffer( vtotal );
01590 
01591     // Create book-keeping arrays over the parent mesh to avoid introducing duplicate vertices
01592     std::vector< EntityHandle > trackvertsC_edg( nepc * ne * nents_prev, 0 );
01593     std::vector< EntityHandle > trackvertsC_face( nfpc * nvf * nents_prev, 0 );
01594 
01595     int cur_nverts = level_mesh[cur_level].num_verts;
01596     std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
01597     std::vector< int > cell_patterns( nents_prev, 0 );
01598 
01599     int count_nverts = nverts_prev;
01600     int count_ents   = 0;
01601     std::vector< EntityHandle > conn, cur_conn;
01602 
01603     // Step 1: Create the subentities via refinement of the previous mesh
01604     for( int cid = 0; cid < nents_prev; cid++ )
01605     {
01606         conn.clear();
01607         cur_conn.clear();
01608         for( int i = 0; i < vtotal; i++ )
01609             vbuffer[i] = 0;
01610 
01611         // EntityHandle of the working cell
01612         EntityHandle cell;
01613         if( cur_level )
01614             cell = level_mesh[cur_level - 1].start_cell + cid;
01615         else
01616             cell = _incells[cid];
01617 
01618         error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
01619 
01620         // Step 1: Add vertices from the current level for the working face that will be used for
01621         // subdivision.
01622         // Add the corners to vbuffer
01623         for( int i = 0; i < (int)conn.size(); i++ )
01624         {
01625             if( cur_level )
01626                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
01627             else
01628                 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
01629 
01630             cur_conn.push_back( vbuffer[i] );
01631         }
01632 
01633         // Gather vertices already added to tracking array due to refinement of the sibling cells
01634         for( int i = 0; i < nepc; i++ )
01635         {
01636             for( int j = 0; j < ne; j++ )
01637             {
01638                 int idx      = refTemplates[cindex][d].vert_on_edges[i][j];
01639                 vbuffer[idx] = trackvertsC_edg[cid * nepc * ne + ne * i + j];
01640             }
01641         }
01642 
01643         // Add remaining new vertex handles
01644         for( int i = 0; i < nfpc; i++ )
01645         {
01646             for( int j = 0; j < nvf; j++ )
01647             {
01648                 int idx      = refTemplates[cindex][d].vert_on_faces[i][j];
01649                 vbuffer[idx] = trackvertsC_face[cid * nfpc * nvf + nvf * i + j];
01650             }
01651         }
01652 
01653         // Add the remaining vertex handles to vbuffer for the current level for the working cell
01654         for( int i = 0; i < nvtotal; i++ )
01655         {
01656             if( !vbuffer[i + nvpc] )
01657             {
01658                 vbuffer[i + nvpc] = level_mesh[cur_level].start_vertex + count_nverts;
01659                 count_nverts += 1;
01660             }
01661         }
01662 
01663         // Step 2: Coordinates of the new vertices
01664         std::vector< double > corner_coords( nvpc * 3 );
01665         error = get_coordinates( &cur_conn[0], nvpc, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
01666 
01667         error = compute_coordinates( cur_level, deg, type, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
01668                                      nverts_prev );MB_CHK_ERR( error );
01669 
01670         // Step 3: Choose the tet refine pattern to be used for this tet
01671         int diag           = find_shortest_diagonal_octahedron( cur_level, deg, &vbuffer[0] );
01672         int pat_id         = diag + 2;
01673         cell_patterns[cid] = pat_id;
01674 
01675         // Step 4: Use the template to obtain the subentities. The coordinates and local ahf maps
01676         // are also constructed. Connectivity of the children
01677         int etotal = refTemplates[pat_id][d].total_new_ents;
01678         std::vector< EntityHandle > ent_buffer( etotal );
01679 
01680         for( int i = 0; i < etotal; i++ )
01681         {
01682             for( int k = 0; k < nvpc; k++ )
01683             {
01684                 int idx                                                = refTemplates[pat_id][d].ents_conn[i][k];
01685                 level_mesh[cur_level].cell_conn[nvpc * count_ents + k] = vbuffer[idx];
01686             }
01687             ent_buffer[i] = level_mesh[cur_level].start_cell + count_ents;
01688             count_ents += 1;
01689         }
01690 
01691         // Step 5: Update local ahf maps
01692         error = update_local_ahf( deg, MBTET, pat_id, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
01693 
01694         // Step 6: Update tracking information
01695         error = update_tracking_verts( cell, cur_level, deg, trackvertsC_edg, trackvertsC_face, &vbuffer[0] );MB_CHK_ERR( error );
01696     }
01697 
01698     // Step 7: Update the global maps
01699     //  error = update_global_ahf(cur_level, deg, cell_patterns); MB_CHK_ERR(error);
01700     error = update_global_ahf( type, cur_level, deg, &cell_patterns );MB_CHK_ERR( error );
01701 
01702     // Step 8: If edges exists, refine them as well.
01703     if( level_mesh[cur_level].num_edges != 0 )
01704     {
01705         error = construct_hm_1D( cur_level, deg, type, trackvertsC_edg );MB_CHK_ERR( error );
01706     }
01707 
01708     // Step 9: If faces exists, refine them as well.
01709     if( !_infaces.empty() )
01710     {
01711         error = construct_hm_2D( cur_level, deg, type, trackvertsC_edg, trackvertsC_face );MB_CHK_ERR( error );
01712     }
01713 
01714     return MB_SUCCESS;
01715 }
01716 
01717 ErrorCode NestedRefine::compute_coordinates( int cur_level, int deg, EntityType type, EntityHandle* vbuffer, int vtotal,
01718                                              double* corner_coords, std::vector< int >& vflag, int nverts_prev )
01719 {
01720     EntityHandle vstart = level_mesh[cur_level].start_vertex;
01721     int d               = get_index_from_degree( deg );
01722 
01723     if( type == MBTRI )
01724     {
01725         double xi, eta, N[3];
01726         int findex = mbImpl->type_from_handle( *( _infaces.begin() ) ) - 1;
01727 
01728         for( int i = 3; i < vtotal; i++ )
01729         {
01730             if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
01731 
01732             xi   = refTemplates[findex][d].vert_nat_coord[i - 3][0];
01733             eta  = refTemplates[findex][d].vert_nat_coord[i - 3][1];
01734             N[0] = 1 - xi - eta;
01735             N[1] = xi;
01736             N[2] = eta;
01737 
01738             double x = 0, y = 0, z = 0;
01739             for( int j = 0; j < 3; j++ )
01740             {
01741                 x += N[j] * corner_coords[3 * j];
01742                 y += N[j] * corner_coords[3 * j + 1];
01743                 z += N[j] * corner_coords[3 * j + 2];
01744             }
01745 
01746             level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
01747             level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
01748             level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
01749             vflag[vbuffer[i] - vstart - nverts_prev]                  = 1;
01750         }
01751     }
01752     else if( type == MBQUAD )
01753     {
01754         double xi, eta, N[4];
01755         int findex = mbImpl->type_from_handle( *( _infaces.begin() ) ) - 1;
01756 
01757         for( int i = 4; i < vtotal; i++ )
01758         {
01759             if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
01760 
01761             xi   = refTemplates[findex][d].vert_nat_coord[i - 4][0];
01762             eta  = refTemplates[findex][d].vert_nat_coord[i - 4][1];
01763             N[0] = ( 1 - xi ) * ( 1 - eta ) / 4;
01764             N[1] = ( 1 + xi ) * ( 1 - eta ) / 4;
01765             N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
01766 
01767             double x = 0, y = 0, z = 0;
01768             for( int j = 0; j < 4; j++ )
01769             {
01770                 x += N[j] * corner_coords[3 * j];
01771                 y += N[j] * corner_coords[3 * j + 1];
01772                 z += N[j] * corner_coords[3 * j + 2];
01773             }
01774 
01775             level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
01776             level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
01777             level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
01778             vflag[vbuffer[i] - vstart - nverts_prev]                  = 1;
01779         }
01780     }
01781     else if( type == MBTET )
01782     {
01783         double xi, eta, mu, N[4];
01784         int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
01785 
01786         for( int i = 4; i < vtotal; i++ )
01787         {
01788             if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
01789 
01790             xi  = refTemplates[cindex][d].vert_nat_coord[i - 4][0];
01791             eta = refTemplates[cindex][d].vert_nat_coord[i - 4][1];
01792             mu  = refTemplates[cindex][d].vert_nat_coord[i - 4][2];
01793 
01794             N[0] = 1 - xi - eta - mu;
01795             N[1] = xi;
01796             N[2] = eta, N[3] = mu;
01797 
01798             double x = 0, y = 0, z = 0;
01799             for( int j = 0; j < 4; j++ )
01800             {
01801                 x += N[j] * corner_coords[3 * j];
01802                 y += N[j] * corner_coords[3 * j + 1];
01803                 z += N[j] * corner_coords[3 * j + 2];
01804             }
01805 
01806             level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
01807             level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
01808             level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
01809             vflag[vbuffer[i] - vstart - nverts_prev]                  = 1;
01810         }
01811     }
01812     else if( type == MBPRISM )
01813     {
01814         double xi, eta, mu, N[6];
01815         int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
01816 
01817         for( int i = 6; i < vtotal; i++ )
01818         {
01819             if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
01820 
01821             xi  = refTemplates[cindex][d].vert_nat_coord[i - 6][0];
01822             eta = refTemplates[cindex][d].vert_nat_coord[i - 6][1];
01823             mu  = refTemplates[cindex][d].vert_nat_coord[i - 6][2];
01824 
01825             N[0] = ( 1 - xi - eta ) * ( 1 - mu ), N[1] = xi * ( 1 - mu ), N[2] = eta * ( 1 - mu ),
01826             N[3] = ( 1 - xi - eta ) * ( 1 + mu ), N[4] = xi * ( 1 + mu ), N[5] = eta * ( 1 + mu );
01827 
01828             double x = 0, y = 0, z = 0;
01829             for( int j = 0; j < 6; j++ )
01830             {
01831                 x += N[j] * corner_coords[3 * j];
01832                 y += N[j] * corner_coords[3 * j + 1];
01833                 z += N[j] * corner_coords[3 * j + 2];
01834             }
01835 
01836             level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
01837             level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
01838             level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
01839             vflag[vbuffer[i] - vstart - nverts_prev]                  = 1;
01840         }
01841     }
01842     else if( type == MBHEX )
01843     {
01844         double xi, eta, mu, N[8];
01845         double d1, d2, d3, s1, s2, s3;
01846         int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
01847 
01848         for( int i = 8; i < vtotal; i++ )
01849         {
01850 
01851             if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
01852 
01853             xi  = refTemplates[cindex][d].vert_nat_coord[i - 8][0];
01854             eta = refTemplates[cindex][d].vert_nat_coord[i - 8][1];
01855             mu  = refTemplates[cindex][d].vert_nat_coord[i - 8][2];
01856 
01857             d1   = 1 - xi;
01858             d2   = 1 - eta;
01859             d3   = 1 - mu;
01860             s1   = 1 + xi;
01861             s2   = 1 + eta;
01862             s3   = 1 + mu;
01863             N[0] = ( d1 * d2 * d3 ) / 8;
01864             N[1] = ( s1 * d2 * d3 ) / 8;
01865             N[2] = ( s1 * s2 * d3 ) / 8;
01866             N[3] = ( d1 * s2 * d3 ) / 8;
01867             N[4] = ( d1 * d2 * s3 ) / 8;
01868             N[5] = ( s1 * d2 * s3 ) / 8;
01869             N[6] = ( s1 * s2 * s3 ) / 8;
01870             N[7] = ( d1 * s2 * s3 ) / 8;
01871 
01872             double x = 0, y = 0, z = 0;
01873             for( int j = 0; j < 8; j++ )
01874             {
01875                 x += N[j] * corner_coords[3 * j];
01876                 y += N[j] * corner_coords[3 * j + 1];
01877                 z += N[j] * corner_coords[3 * j + 2];
01878             }
01879 
01880             level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
01881             level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
01882             level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
01883 
01884             vflag[vbuffer[i] - vstart - nverts_prev] = 1;
01885         }
01886     }
01887     return MB_SUCCESS;
01888 }
01889 /**********************************
01890  *      Parallel Communication       *
01891  * ********************************/
01892 
01893 #ifdef MOAB_HAVE_MPI
01894 ErrorCode NestedRefine::resolve_shared_ents_parmerge( int level, EntityHandle levelset )
01895 {
01896     // TEMP: Add the adjacencies for MOAB-native DS
01897     // NOTE (VSM): This is expensive since it creates a doubly
01898     // redundant copy of the adjacency data in both MOAB-native
01899     // and AHF. Need to fix this with AHF optimized branch.
01900     ErrorCode error;
01901     ReadUtilIface* read_iface;
01902     error = mbImpl->query_interface( read_iface );MB_CHK_ERR( error );
01903     if( level_mesh[level].num_edges != 0 )
01904     {
01905         error = read_iface->update_adjacencies( level_mesh[level].start_edge, level_mesh[level].num_edges, 2,
01906                                                 level_mesh[level].edge_conn );MB_CHK_ERR( error );
01907     }
01908     if( level_mesh[level].num_faces != 0 )
01909     {
01910         EntityType type = mbImpl->type_from_handle( *( _infaces.begin() ) );
01911         int nvpf        = ahf->lConnMap2D[type - 2].num_verts_in_face;
01912         error = read_iface->update_adjacencies( level_mesh[level].start_face, level_mesh[level].num_faces, nvpf,
01913                                                 level_mesh[level].face_conn );MB_CHK_ERR( error );
01914     }
01915     if( level_mesh[level].num_cells != 0 )
01916     {
01917         int index = ahf->get_index_in_lmap( *_incells.begin() );
01918         int nvpc  = ahf->lConnMap3D[index].num_verts_in_cell;
01919         error     = read_iface->update_adjacencies( level_mesh[level].start_cell, level_mesh[level].num_cells, nvpc,
01920                                                 level_mesh[level].cell_conn );MB_CHK_ERR( error );
01921     }
01922 
01923     if( pcomm->size() > 1 )
01924     {
01925 
01926         // get all entities on the rootset
01927         moab::Range vtxs, edgs, facs, elms;
01928         error = mbImpl->get_entities_by_dimension( levelset, 0, vtxs, false );MB_CHK_ERR( error );
01929         error = mbImpl->get_entities_by_dimension( levelset, 1, edgs, false );MB_CHK_ERR( error );
01930         error = mbImpl->get_entities_by_dimension( levelset, 2, facs, false );MB_CHK_ERR( error );
01931         error = mbImpl->get_entities_by_dimension( levelset, 3, elms, false );MB_CHK_ERR( error );
01932 
01933         // set the parallel partition tag data
01934         moab::Tag part_tag;
01935         int partid = pcomm->rank(), dum_id = -1;
01936         error = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, part_tag,
01937                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id );MB_CHK_ERR( error );
01938         error = mbImpl->tag_set_data( part_tag, &levelset, 1, &partid );MB_CHK_ERR( error );
01939         //
01940         // Now that we have the local piece of the mesh refined consistently,
01941         // call parallel merge instead of resolved_shared to stitch together the meshes
01942         // and to handle parallel communication of remote proc/entity-handle pairs for
01943         // shared + non-owned interfaces entities.
01944         //
01945         // TODO: This needs to be replaced by the following scheme in the future as
01946         // an optimization step.
01947         //   > Assign global IDs consistently for shared entities so that parallel
01948         //   > resolve shared ents can happen out of the box. This is the fastest option.
01949 
01950         ParallelMergeMesh pm( pcomm, 1e-08 );
01951         error = pm.merge( levelset, true );MB_CHK_ERR( error );
01952 
01953         //
01954         // Parallel Communication complete - all entities resolved
01955         //
01956     }
01957     return MB_SUCCESS;
01958 }
01959 
01960 ErrorCode NestedRefine::resolve_shared_ents_opt( EntityHandle* hm_set, int num_levels )
01961 {
01962     assert( pcomm->size() > 1 );
01963 
01964     ErrorCode error;
01965 
01966     // Step 1A: Pre-processing:  setting the parallel partition tag data
01967     for( int i = 0; i < num_levels; i++ )
01968     {
01969         Tag part_tag;
01970         int partid = pcomm->rank(), dum_id = -1;
01971         error = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, part_tag,
01972                                         moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id );MB_CHK_ERR( error );
01973         error = mbImpl->tag_set_data( part_tag, &hm_set[i], 1, &partid );MB_CHK_ERR( error );
01974     }
01975 
01976     // Step 1B: Pre-processing: gather all shared entities and list entities shared with each
01977     // sharing processor
01978 
01979     // All shared processors
01980     std::set< unsigned int > shprocs;
01981     error = pcomm->get_comm_procs( shprocs );MB_CHK_ERR( error );
01982 
01983     std::vector< int > sharedprocs;
01984     for( std::set< unsigned int >::iterator it = shprocs.begin(); it != shprocs.end(); it++ )
01985         sharedprocs.push_back( *it );
01986     int nprocs = sharedprocs.size();
01987 
01988     // Create buffer variables storing the entities to be sent and received
01989     std::vector< std::vector< int > > nsharedEntsperproc( nprocs );
01990     std::vector< std::vector< EntityHandle > > localBuffs( nprocs );
01991     std::vector< std::vector< EntityHandle > > remlocalBuffs( nprocs );
01992     std::vector< std::vector< EntityHandle > > remoteBuffs( nprocs );
01993 
01994     int i;
01995     Range sharedentities;
01996 
01997     for( i = 0; i < nprocs; i++ )
01998     {
01999         // List of shared entities at the coarsest level
02000         sharedentities.clear();
02001         error = pcomm->get_shared_entities( sharedprocs[i], sharedentities, -1, true );MB_CHK_ERR( error );
02002 
02003         // Get the list shared edges and vertices that are not part of the shared edges
02004         Range allEnts;
02005         error = collect_shared_entities_by_dimension( sharedentities, allEnts );MB_CHK_ERR( error );
02006 
02007         Range V0, E0, F0;
02008         V0 = allEnts.subset_by_dimension( 0 );
02009         E0 = allEnts.subset_by_dimension( 1 );
02010         F0 = allEnts.subset_by_dimension( 2 );
02011 
02012         // Step 2A: Prepare msg to be sent:
02013         //
02014         // FList = <F, FC, FE> where F, FC and FE are vectors containing the faces, their
02015         // connectivities and edges. F = <F_0, F_1, F_2, ..., F_L> where
02016         //       F_0 is the list of shared faces at the coarsest level,
02017         //       F_i are the list of children faces at subsequent levels.
02018         // FC vector contains the connectivities of F_i's and FE vector contains the bounding edges
02019         // of F_i's.
02020         //
02021         // EList = <E, EC> where E and EC are vectors containing the edges and their connectivities.
02022         // E = <E_0, E_1, ..., E_L> where
02023         //     E_0 is the list of shared edges at the coarsest level and
02024         //     E_i are the list of children edges at subsequent levels.
02025         // The EC vector contains the connectivities of E_i's.
02026         //
02027         // VList = <V> = <V_0, V_1, ...., V_L> where V_0 are shared vertices at the coarsest level
02028         //        and V_i are the duplicates in the subsequent levels.
02029 
02030         std::vector< EntityHandle > locFList, remFList;
02031         std::vector< EntityHandle > locEList, remEList;
02032         std::vector< EntityHandle > locVList, remVList;
02033 
02034         // collect faces
02035         if( !F0.empty() )
02036         {
02037             error = collect_FList( sharedprocs[i], F0, locFList, remFList );MB_CHK_ERR( error );
02038         }
02039 
02040         // collect edges
02041         if( !E0.empty() )
02042         {
02043             error = collect_EList( sharedprocs[i], E0, locEList, remEList );MB_CHK_ERR( error );
02044         }
02045 
02046         // collect vertices
02047         if( !V0.empty() )
02048         {
02049             error = collect_VList( sharedprocs[i], V0, locVList, remVList );MB_CHK_ERR( error );
02050         }
02051 
02052         // Step 2B: Add data to NR local buffer to be sent
02053         std::vector< int > msgsz;
02054         msgsz.push_back( F0.size() );
02055         msgsz.push_back( locFList.size() );
02056         msgsz.push_back( E0.size() );
02057         msgsz.push_back( locEList.size() );
02058         msgsz.push_back( V0.size() );
02059         msgsz.push_back( locVList.size() );
02060         nsharedEntsperproc[i].insert( nsharedEntsperproc[i].end(), msgsz.begin(), msgsz.end() );
02061 
02062         if( !F0.empty() )
02063         {
02064             localBuffs[i].insert( localBuffs[i].end(), locFList.begin(), locFList.end() );
02065             remlocalBuffs[i].insert( remlocalBuffs[i].end(), remFList.begin(), remFList.end() );
02066         }
02067         if( !E0.empty() )
02068         {
02069             localBuffs[i].insert( localBuffs[i].end(), locEList.begin(), locEList.end() );
02070             remlocalBuffs[i].insert( remlocalBuffs[i].end(), remEList.begin(), remEList.end() );
02071         }
02072         if( !V0.empty() )
02073         {
02074             localBuffs[i].insert( localBuffs[i].end(), locVList.begin(), locVList.end() );
02075             remlocalBuffs[i].insert( remlocalBuffs[i].end(), remVList.begin(), remVList.end() );
02076         }
02077     }
02078 
02079     // Step 3: Send and receive the remote collection of child ents
02080     error = pcomm->send_recv_entities( sharedprocs, nsharedEntsperproc, remlocalBuffs, remoteBuffs );MB_CHK_ERR( error );
02081 
02082     // Step 5: Resolve shared child entities and update parallel tags
02083     std::multimap< EntityHandle, int > rprocs;
02084     std::multimap< EntityHandle, EntityHandle > rhandles;
02085 
02086     error = decipher_remote_handles( sharedprocs, nsharedEntsperproc, localBuffs, remoteBuffs, rprocs, rhandles );MB_CHK_ERR( error );
02087 
02088     // Step 6: Update pcomm tags
02089     error = update_parallel_tags( rprocs, rhandles );MB_CHK_ERR( error );
02090 
02091     return MB_SUCCESS;
02092 }
02093 
02094 ErrorCode NestedRefine::collect_shared_entities_by_dimension( Range sharedEnts, Range& allEnts )
02095 {
02096     ErrorCode error;
02097 
02098     Range F0, E0, V0, E0all, V0all;
02099     std::vector< EntityHandle > ents;
02100 
02101     F0    = sharedEnts.subset_by_dimension( 2 );
02102     E0all = sharedEnts.subset_by_dimension( 1 );
02103     V0all = sharedEnts.subset_by_dimension( 0 );
02104 
02105     if( !F0.empty() )
02106     {
02107         Range edges, verts;
02108 
02109         for( Range::iterator it = F0.begin(); it != F0.end(); it++ )
02110         {
02111             ents.clear();
02112             error = ahf->get_adjacencies( *it, 1, ents );MB_CHK_ERR( error );
02113             std::copy( ents.begin(), ents.end(), range_inserter( edges ) );
02114             ents.clear();
02115             error = mbImpl->get_connectivity( &( *it ), 0, ents );MB_CHK_ERR( error );
02116             std::copy( ents.begin(), ents.end(), range_inserter( verts ) );
02117         }
02118 
02119         E0 = subtract( E0all, edges );
02120         if( !E0.empty() )
02121         {
02122             for( Range::iterator it = E0.begin(); it != E0.end(); it++ )
02123             {
02124                 ents.clear();
02125                 error = mbImpl->get_connectivity( &( *it ), 0, ents );MB_CHK_ERR( error );
02126                 std::copy( ents.begin(), ents.end(), range_inserter( verts ) );
02127             }
02128         }
02129         V0 = subtract( V0all, verts );
02130     }
02131     else if( !E0all.empty() )
02132     {
02133         Range verts;
02134         for( Range::iterator it = E0all.begin(); it != E0all.end(); it++ )
02135         {
02136             ents.clear();
02137             error = mbImpl->get_connectivity( &( *it ), 1, ents );MB_CHK_ERR( error );
02138             std::copy( ents.begin(), ents.end(), range_inserter( verts ) );
02139         }
02140         E0 = E0all;
02141         V0 = subtract( V0all, verts );
02142     }
02143     else if( !V0all.empty() )
02144     {
02145         V0 = V0all;
02146     }
02147     else
02148         MB_SET_ERR( MB_FAILURE, "Trying to pack unsupported sub-entities for shared interface entities" );
02149 
02150     if( !F0.empty() ) std::copy( F0.begin(), F0.end(), range_inserter( allEnts ) );
02151     if( !E0.empty() ) std::copy( E0.begin(), E0.end(), range_inserter( allEnts ) );
02152     if( !V0.empty() ) std::copy( V0.begin(), V0.end(), range_inserter( allEnts ) );
02153 
02154     return MB_SUCCESS;
02155 }
02156 
02157 ErrorCode NestedRefine::collect_FList( int to_proc, Range faces, std::vector< EntityHandle >& FList,
02158                                        std::vector< EntityHandle >& RList )
02159 {
02160     ErrorCode error;
02161 
02162     FList.clear();
02163     std::vector< EntityHandle > F, FC, FE, lF, lFC, lFE, rF, rFC, rFE;
02164     std::vector< EntityHandle > childEnts, conn, fedges;
02165 
02166     for( Range::iterator it = faces.begin(); it != faces.end(); it++ )
02167     {
02168         EntityHandle face = *it;
02169         conn.clear();
02170         fedges.clear();
02171         error = mbImpl->get_connectivity( &face, 1, conn );MB_CHK_ERR( error );
02172         error = ahf->get_face_edges( *it, fedges );MB_CHK_ERR( error );
02173 
02174         // add local handles
02175         lF.push_back( *it );
02176         lFC.insert( lFC.end(), conn.begin(), conn.end() );
02177         lFE.insert( lFE.end(), fedges.begin(), fedges.end() );
02178 
02179         // replace local handles with remote handles on to_proc
02180         EntityHandle rval = 0;
02181         error             = pcomm->get_remote_handles( &face, &rval, 1, to_proc );MB_CHK_ERR( error );
02182         rF.push_back( rval );
02183 
02184         for( int i = 0; i < (int)conn.size(); i++ )
02185         {
02186             error = pcomm->get_remote_handles( &conn[i], &rval, 1, to_proc );MB_CHK_ERR( error );
02187             rFC.push_back( rval );
02188 
02189             error = pcomm->get_remote_handles( &fedges[i], &rval, 1, to_proc );MB_CHK_ERR( error );
02190             rFE.push_back( rval );
02191         }
02192     }
02193 
02194     for( int l = 0; l < nlevels; l++ )
02195     {
02196         for( Range::iterator it = faces.begin(); it != faces.end(); it++ )
02197         {
02198             childEnts.clear();
02199             error = parent_to_child( *it, 0, l + 1, childEnts );MB_CHK_ERR( error );
02200 
02201             for( int i = 0; i < (int)childEnts.size(); i++ )
02202             {
02203                 conn.clear();
02204                 fedges.clear();
02205                 error = mbImpl->get_connectivity( &childEnts[i], 1, conn );MB_CHK_ERR( error );
02206                 error = ahf->get_face_edges( childEnts[i], fedges );MB_CHK_ERR( error );
02207 
02208                 F.push_back( childEnts[i] );
02209                 FC.insert( FC.end(), conn.begin(), conn.end() );
02210                 FE.insert( FE.end(), fedges.begin(), fedges.end() );
02211             }
02212         }
02213     }
02214 
02215     FList.insert( FList.end(), lF.begin(), lF.end() );
02216     FList.insert( FList.end(), F.begin(), F.end() );
02217     FList.insert( FList.end(), lFC.begin(), lFC.end() );
02218     FList.insert( FList.end(), FC.begin(), FC.end() );
02219     FList.insert( FList.end(), lFE.begin(), lFE.end() );
02220     FList.insert( FList.end(), FE.begin(), FE.end() );
02221 
02222     RList.insert( RList.end(), rF.begin(), rF.end() );
02223     RList.insert( RList.end(), F.begin(), F.end() );
02224     RList.insert( RList.end(), rFC.begin(), rFC.end() );
02225     RList.insert( RList.end(), FC.begin(), FC.end() );
02226     RList.insert( RList.end(), rFE.begin(), rFE.end() );
02227     RList.insert( RList.end(), FE.begin(), FE.end() );
02228 
02229     return MB_SUCCESS;
02230 }
02231 
02232 ErrorCode NestedRefine::collect_EList( int to_proc, Range edges, std::vector< EntityHandle >& EList,
02233                                        std::vector< EntityHandle >& RList )
02234 {
02235     ErrorCode error;
02236     EList.clear();
02237     std::vector< EntityHandle > E, EC, lE, lEC, rE, rEC;
02238     std::vector< EntityHandle > childEnts, conn;
02239 
02240     // Add the edges and their connectivities at the coarsest level first.
02241     for( Range::iterator it = edges.begin(); it != edges.end(); it++ )
02242     {
02243         EntityHandle edg = *it;
02244         conn.clear();
02245         error = mbImpl->get_connectivity( &edg, 1, conn );MB_CHK_ERR( error );
02246 
02247         // add local handles
02248         lE.push_back( edg );
02249         lEC.insert( lEC.end(), conn.begin(), conn.end() );
02250 
02251         // replace local handles with remote handle on to_proc
02252         EntityHandle rval = 0;
02253         error             = pcomm->get_remote_handles( &edg, &rval, 1, to_proc );MB_CHK_ERR( error );
02254         rE.push_back( rval );
02255         error = pcomm->get_remote_handles( &conn[0], &rval, 1, to_proc );MB_CHK_ERR( error );
02256         rEC.push_back( rval );
02257         error = pcomm->get_remote_handles( &conn[1], &rval, 1, to_proc );MB_CHK_ERR( error );
02258         rEC.push_back( rval );
02259     }
02260 
02261     // Add the edges and their connectivities at subsequent levels.
02262     for( int l = 0; l < nlevels; l++ )
02263     {
02264         for( Range::iterator it = edges.begin(); it != edges.end(); it++ )
02265         {
02266             childEnts.clear();
02267             error = parent_to_child( *it, 0, l + 1, childEnts );MB_CHK_ERR( error );
02268 
02269             for( int i = 0; i < (int)childEnts.size(); i++ )
02270             {
02271                 conn.clear();
02272                 error = mbImpl->get_connectivity( &childEnts[i], 1, conn );MB_CHK_ERR( error );
02273                 E.push_back( childEnts[i] );
02274                 EC.insert( EC.end(), conn.begin(), conn.end() );
02275             }
02276         }
02277     }
02278 
02279     EList.insert( EList.end(), lE.begin(), lE.end() );
02280     EList.insert( EList.end(), E.begin(), E.end() );
02281     EList.insert( EList.end(), lEC.begin(), lEC.end() );
02282     EList.insert( EList.end(), EC.begin(), EC.end() );
02283 
02284     RList.insert( RList.end(), rE.begin(), rE.end() );
02285     RList.insert( RList.end(), E.begin(), E.end() );
02286     RList.insert( RList.end(), rEC.begin(), rEC.end() );
02287     RList.insert( RList.end(), EC.begin(), EC.end() );
02288 
02289     return MB_SUCCESS;
02290 }
02291 
02292 ErrorCode NestedRefine::collect_VList( int to_proc, Range verts, std::vector< EntityHandle >& VList,
02293                                        std::vector< EntityHandle >& RList )
02294 {
02295     ErrorCode error;
02296     std::vector< EntityHandle > V, lV, rV;
02297     VList.clear();
02298 
02299     // Add the vertices at the coarsest level first.
02300     for( Range::iterator it = verts.begin(); it != verts.end(); it++ )
02301     {
02302         EntityHandle v = *it;
02303 
02304         lV.push_back( v );
02305         EntityHandle rval = 0;
02306         error             = pcomm->get_remote_handles( &v, &rval, 1, to_proc );MB_CHK_ERR( error );
02307 
02308         rV.push_back( rval );
02309     }
02310 
02311     // Add the vertices at the subsequent levels .
02312     for( int l = 0; l < nlevels; l++ )
02313     {
02314         for( Range::iterator it = verts.begin(); it != verts.end(); it++ )
02315         {
02316             EntityHandle dupvert = 0;
02317             error                = get_vertex_duplicates( *it, l + 1, dupvert );MB_CHK_ERR( error );
02318             V.push_back( dupvert );
02319         }
02320     }
02321 
02322     // local vertex handles at the coarsest level
02323     VList.insert( VList.end(), lV.begin(), lV.end() );
02324     VList.insert( VList.end(), V.begin(), V.end() );
02325 
02326     // remote vertex handles at the coarsest level
02327     RList.insert( RList.end(), rV.begin(), rV.end() );
02328     RList.insert( RList.end(), V.begin(), V.end() );
02329 
02330     return MB_SUCCESS;
02331 }
02332 
02333 ErrorCode NestedRefine::decipher_remote_handles( std::vector< int >& sharedprocs,
02334                                                  std::vector< std::vector< int > >& auxinfo,
02335                                                  std::vector< std::vector< EntityHandle > >& localbuffers,
02336                                                  std::vector< std::vector< EntityHandle > >& remotebuffers,
02337                                                  std::multimap< EntityHandle, int >& remProcs,
02338                                                  std::multimap< EntityHandle, EntityHandle >& remHandles )
02339 {
02340     ErrorCode error;
02341 
02342     int i;
02343     int nprocs = sharedprocs.size();
02344 
02345     for( i = 0; i < nprocs; i++ )
02346     {
02347         std::vector< int > msgsz;
02348         for( int j = 0; j < (int)auxinfo[i].size(); j++ )
02349         {
02350             msgsz.push_back( auxinfo[i][j] );
02351         }
02352 
02353         if( msgsz[0] != 0 )  // Faces
02354         {
02355             // Get the local and remote face handles from the buffers
02356             std::vector< EntityHandle > LFList, RFList;
02357             LFList.insert( LFList.end(), localbuffers[i].begin(), localbuffers[i].begin() + msgsz[1] );
02358             RFList.insert( RFList.end(), remotebuffers[i].begin(), remotebuffers[i].begin() + msgsz[1] );
02359 
02360             error = decipher_remote_handles_face( sharedprocs[i], msgsz[0], LFList, RFList, remProcs, remHandles );MB_CHK_ERR( error );
02361 
02362             if( msgsz[2] != 0 )  // Edges
02363             {
02364                 std::vector< EntityHandle > LEList, REList;
02365                 LEList.insert( LEList.end(), localbuffers[i].begin() + msgsz[1] + 1,
02366                                localbuffers[i].begin() + msgsz[1] + msgsz[3] );
02367                 REList.insert( REList.end(), remotebuffers[i].begin() + msgsz[1] + 1,
02368                                remotebuffers[i].begin() + msgsz[1] + msgsz[3] );
02369 
02370                 error = decipher_remote_handles_edge( sharedprocs[i], msgsz[2], LEList, REList, remProcs, remHandles );MB_CHK_ERR( error );
02371 
02372                 if( msgsz[4] != 0 )  // Vertices
02373                 {
02374                     // Get the local and remote face handles from the buffers
02375                     std::vector< EntityHandle > LVList, RVList;
02376                     LVList.insert( LVList.end(), localbuffers[i].begin() + msgsz[1] + msgsz[3] + 1,
02377                                    localbuffers[i].begin() + msgsz[1] + msgsz[3] + msgsz[5] );
02378                     RVList.insert( RVList.end(), remotebuffers[i].begin() + msgsz[1] + msgsz[3] + 1,
02379                                    remotebuffers[i].begin() + msgsz[1] + msgsz[3] + msgsz[5] );
02380 
02381                     error = decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs,
02382                                                             remHandles );MB_CHK_ERR( error );
02383                 }
02384             }
02385             else if( msgsz[4] != 0 )  // Vertices
02386             {
02387                 // Get the local and remote face handles from the buffers
02388                 std::vector< EntityHandle > LVList, RVList;
02389                 LVList.insert( LVList.end(), localbuffers[i].begin() + msgsz[1] + 1,
02390                                localbuffers[i].begin() + msgsz[1] + msgsz[5] );
02391                 RVList.insert( RVList.end(), remotebuffers[i].begin() + msgsz[1] + 1,
02392                                remotebuffers[i].begin() + msgsz[1] + msgsz[5] );
02393 
02394                 error =
02395                     decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs, remHandles );MB_CHK_ERR( error );
02396             }
02397         }
02398 
02399         else if( msgsz[2] != 0 )  // Edges
02400         {
02401             // Get the local and remote face handles from the buffers
02402             std::vector< EntityHandle > LEList, REList;
02403             LEList.insert( LEList.end(), localbuffers[i].begin(), localbuffers[i].begin() + msgsz[3] );
02404             REList.insert( REList.end(), remotebuffers[i].begin(), remotebuffers[i].begin() + msgsz[3] );
02405 
02406             error = decipher_remote_handles_edge( sharedprocs[i], msgsz[2], LEList, REList, remProcs, remHandles );MB_CHK_ERR( error );
02407 
02408             if( msgsz[4] != 0 )  // Vertices
02409             {
02410                 // Get the local and remote face handles from the buffers
02411                 std::vector< EntityHandle > LVList, RVList;
02412                 LVList.insert( LVList.end(), localbuffers[i].begin() + msgsz[3] + 1,
02413                                localbuffers[i].begin() + msgsz[3] + msgsz[5] );
02414                 RVList.insert( RVList.end(), remotebuffers[i].begin() + msgsz[3] + 1,
02415                                remotebuffers[i].begin() + msgsz[3] + msgsz[5] );
02416 
02417                 error =
02418                     decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs, remHandles );MB_CHK_ERR( error );
02419             }
02420         }
02421 
02422         else if( msgsz[4] != 0 )  // Vertices
02423         {
02424             // Get the local and remote face handles from the buffers
02425             std::vector< EntityHandle > LVList, RVList;
02426             LVList.insert( LVList.end(), localbuffers[i].begin(), localbuffers[i].end() );
02427             RVList.insert( RVList.end(), remotebuffers[i].begin(), remotebuffers[i].end() );
02428 
02429             error = decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs, remHandles );MB_CHK_ERR( error );
02430         }
02431         else
02432             MB_SET_ERR( MB_FAILURE, "Trying to decipher entities other than verts, edges, faces" );
02433     }
02434 
02435     return MB_SUCCESS;
02436 }
02437 
02438 ErrorCode NestedRefine::decipher_remote_handles_face( int shared_proc, int numfaces,
02439                                                       std::vector< EntityHandle >& localFaceList,
02440                                                       std::vector< EntityHandle >& remFaceList,
02441                                                       std::multimap< EntityHandle, int >& remProcs,
02442                                                       std::multimap< EntityHandle, EntityHandle >& remHandles )
02443 {
02444     ErrorCode error;
02445 
02446     for( int i = 0; i < numfaces; i++ )
02447     {
02448         // find local and remote handles of the coarsest face
02449         EntityHandle Lface = localFaceList[i];
02450         int Rface_idx =
02451             ( std::find( remFaceList.begin(), remFaceList.begin() + numfaces - 1, Lface ) ) - remFaceList.begin();
02452 
02453         // get connectivities of the local and remote coarsest faces
02454         std::vector< EntityHandle > Lface_conn, Rface_conn;
02455         error = get_data_from_buff( 2, 1, 0, i, numfaces, localFaceList, Lface_conn );MB_CHK_ERR( error );
02456         error = get_data_from_buff( 2, 1, 0, Rface_idx, numfaces, remFaceList, Rface_conn );MB_CHK_ERR( error );
02457 
02458         // find the combination difference between local and remote coarsest face
02459         std::vector< int > cmap;
02460         int comb = 0;
02461         int nvF  = (int)Lface_conn.size();
02462         error    = reorder_indices( &Lface_conn[0], &Rface_conn[0], nvF, &cmap[0], comb );MB_CHK_ERR( error );
02463 
02464         // go into loop over all levels
02465         std::vector< EntityHandle > lchildents, lparents;
02466         std::vector< EntityHandle > lcents, rcents;
02467         int lidx, ridx;
02468 
02469         for( int l = 0; l < nlevels; l++ )
02470         {
02471             lchildents.clear();
02472             lparents.clear();
02473             lcents.clear();
02474             rcents.clear();
02475 
02476             // obtain children at the current level
02477             error = get_data_from_buff( 2, 0, l + 1, i, numfaces, localFaceList, lchildents );MB_CHK_ERR( error );
02478 
02479             // obtain parents at the previous level
02480             if( l == 0 ) { lparents.push_back( Lface ); }
02481             else
02482             {
02483                 error = get_data_from_buff( 2, 0, l, i, numfaces, localFaceList, lparents );MB_CHK_ERR( error );
02484             }
02485 
02486             //#children at the previous level and the current level
02487             EntityType ftype = mbImpl->type_from_handle( Lface );
02488             int d            = get_index_from_degree( level_dsequence[l] );
02489             int nch          = refTemplates[ftype - 1][d].total_new_ents;
02490             std::vector< int > fmap;
02491             error = reorder_indices( level_dsequence[l], nvF, comb, &fmap[0] );MB_CHK_ERR( error );
02492 
02493             // loop over all the lparents
02494             for( int j = 0; j < (int)lparents.size(); j++ )
02495             {
02496                 // list local childrent at the current level
02497                 lidx  = std::find( localFaceList.begin(), localFaceList.end(), lparents[j] ) - localFaceList.begin();
02498                 error = get_data_from_buff( 2, 0, l + 1, lidx, numfaces, localFaceList, lcents );MB_CHK_ERR( error );
02499 
02500                 // find the corresponding remote of lparent and its children
02501                 EntityHandle rparent = 0;
02502                 error = check_for_parallelinfo( lparents[j], shared_proc, remHandles, remProcs, rparent );MB_CHK_ERR( error );
02503                 ridx  = std::find( remFaceList.begin(), remFaceList.end(), rparent ) - remFaceList.begin();
02504                 error = get_data_from_buff( 2, 0, l + 1, ridx, numfaces, remFaceList, rcents );MB_CHK_ERR( error );
02505 
02506                 // match up local face with remote handles according to cmap
02507                 std::vector< EntityHandle > lconn, rconn, ledg, redg;
02508                 for( int k = 0; k < nch; k++ )
02509                 {
02510                     lconn.clear();
02511                     rconn.clear();
02512                     ledg.clear();
02513                     redg.clear();
02514 
02515                     // matching the local children face handles with their remotes
02516                     bool found = check_for_parallelinfo( lcents[k], shared_proc, remProcs );
02517                     if( !found )
02518                     {
02519                         remProcs.insert( std::pair< EntityHandle, int >( lcents[k], shared_proc ) );
02520                         remHandles.insert( std::pair< EntityHandle, EntityHandle >( lcents[k], rcents[fmap[k]] ) );
02521                     }
02522 
02523                     // find indices of the matched child face
02524                     lidx = std::find( localFaceList.begin(), localFaceList.end(), lcents[k] ) - localFaceList.begin();
02525                     ridx = std::find( remFaceList.begin(), remFaceList.end(), rcents[fmap[k]] ) - remFaceList.begin();
02526 
02527                     // find bounding edges and connectivity of the matched child face
02528                     error = get_data_from_buff( 2, 2, l + 1, lidx, numfaces, localFaceList, ledg );MB_CHK_ERR( error );
02529                     error = get_data_from_buff( 2, 1, l + 1, lidx, numfaces, localFaceList, lconn );MB_CHK_ERR( error );
02530                     error = get_data_from_buff( 2, 2, l + 1, ridx, numfaces, remFaceList, redg );MB_CHK_ERR( error );
02531                     error = get_data_from_buff( 2, 1, l + 1, ridx, numfaces, remFaceList, rconn );MB_CHK_ERR( error );
02532 
02533                     // now match the handles of the bounding edges and the vertices using
02534                     // combination difference
02535                     for( int m = 0; m < (int)ledg.size(); m++ )
02536                     {
02537                         found = check_for_parallelinfo( ledg[m], shared_proc, remProcs );
02538 
02539                         if( !found )
02540                         {
02541                             remProcs.insert( std::pair< EntityHandle, int >( ledg[m], shared_proc ) );
02542                             remHandles.insert( std::pair< EntityHandle, EntityHandle >( ledg[m], redg[cmap[m]] ) );
02543                         }
02544 
02545                         found = check_for_parallelinfo( lconn[m], shared_proc, remProcs );
02546 
02547                         if( !found )
02548                         {
02549                             remProcs.insert( std::pair< EntityHandle, int >( lconn[m], shared_proc ) );
02550                             remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[m], rconn[cmap[m]] ) );
02551                         }
02552                     }
02553                 }
02554             }
02555         }
02556     }
02557 
02558     return MB_SUCCESS;
02559 }
02560 
02561 ErrorCode NestedRefine::decipher_remote_handles_edge( int shared_proc, int numedges,
02562                                                       std::vector< EntityHandle >& localEdgeList,
02563                                                       std::vector< EntityHandle >& remEdgeList,
02564                                                       std::multimap< EntityHandle, int >& remProcs,
02565                                                       std::multimap< EntityHandle, EntityHandle >& remHandles )
02566 {
02567     ErrorCode error;
02568 
02569     for( int i = 0; i < numedges; i++ )
02570     {
02571         EntityHandle Ledge = localEdgeList[i];
02572         int Redge_idx =
02573             ( std::find( remEdgeList.begin(), remEdgeList.begin() + numedges - 1, Ledge ) ) - remEdgeList.begin();
02574 
02575         std::vector< EntityHandle > Ledge_conn, Redge_conn;
02576         error = get_data_from_buff( 1, 1, 0, i, numedges, localEdgeList, Ledge_conn );MB_CHK_ERR( error );
02577         error = get_data_from_buff( 1, 1, 0, Redge_idx, numedges, remEdgeList, Redge_conn );MB_CHK_ERR( error );
02578 
02579         bool orient = true;
02580         if( ( Ledge_conn[0] == Redge_conn[1] ) && ( Ledge_conn[1] == Redge_conn[0] ) ) orient = false;
02581 
02582         if( orient ) assert( ( Ledge_conn[0] == Redge_conn[0] ) && ( Ledge_conn[1] == Redge_conn[1] ) );
02583 
02584         std::vector< EntityHandle > lchildEdgs, rchildEdgs, lconn, rconn;
02585         for( int l = 0; l < nlevels; l++ )
02586         {
02587             lchildEdgs.clear();
02588             rchildEdgs.clear();
02589             error = get_data_from_buff( 1, 0, l + 1, i, numedges, localEdgeList, lchildEdgs );MB_CHK_ERR( error );
02590             error = get_data_from_buff( 1, 0, l + 1, Redge_idx, numedges, remEdgeList, rchildEdgs );MB_CHK_ERR( error );
02591 
02592             int nchd = lchildEdgs.size();
02593             if( orient )
02594             {
02595                 for( int j = 0; j < nchd; j++ )
02596                 {
02597                     // match entityhandles of child edges
02598                     bool found = check_for_parallelinfo( lchildEdgs[j], shared_proc, remProcs );
02599                     if( !found )
02600                     {
02601                         remProcs.insert( std::pair< EntityHandle, int >( lchildEdgs[j], shared_proc ) );
02602                         remHandles.insert( std::pair< EntityHandle, EntityHandle >( lchildEdgs[j], rchildEdgs[j] ) );
02603                     }
02604 
02605                     // match entityhandles of child vertices
02606                     lconn.clear();
02607                     rconn.clear();
02608 
02609                     int lidx =
02610                         std::find( localEdgeList.begin(), localEdgeList.end(), lchildEdgs[j] ) - localEdgeList.begin();
02611                     int ridx = std::find( remEdgeList.begin(), remEdgeList.end(), rchildEdgs[j] ) - remEdgeList.begin();
02612 
02613                     error = get_data_from_buff( 1, 1, l + 1, lidx, numedges, localEdgeList, lconn );MB_CHK_ERR( error );
02614                     error = get_data_from_buff( 1, 1, l + 1, ridx, numedges, remEdgeList, rconn );MB_CHK_ERR( error );
02615 
02616                     found = check_for_parallelinfo( lconn[0], shared_proc, remProcs );
02617                     if( !found )
02618                     {
02619                         remProcs.insert( std::pair< EntityHandle, int >( lconn[0], shared_proc ) );
02620                         remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[0], rconn[0] ) );
02621                     }
02622                     found = check_for_parallelinfo( lconn[1], shared_proc, remProcs );
02623                     if( !found )
02624                     {
02625                         remProcs.insert( std::pair< EntityHandle, int >( lconn[1], shared_proc ) );
02626                         remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[1], rconn[0] ) );
02627                     }
02628                 }
02629             }
02630 
02631             else
02632             {
02633                 for( int j = 0; j < nchd; j++ )
02634                 {
02635                     // match entityhandles of child edges
02636                     bool found = check_for_parallelinfo( lchildEdgs[j], shared_proc, remProcs );
02637 
02638                     if( !found )
02639                     {
02640                         remProcs.insert( std::pair< EntityHandle, int >( lchildEdgs[j], shared_proc ) );
02641                         remHandles.insert(
02642                             std::pair< EntityHandle, EntityHandle >( lchildEdgs[j], rchildEdgs[nchd - j - 1] ) );
02643                     }
02644 
02645                     // match entityhandles of child vertices
02646                     lconn.clear();
02647                     rconn.clear();
02648 
02649                     int lidx =
02650                         std::find( localEdgeList.begin(), localEdgeList.end(), lchildEdgs[j] ) - localEdgeList.begin();
02651                     int ridx = std::find( remEdgeList.begin(), remEdgeList.end(), rchildEdgs[nchd - j - 1] ) -
02652                                remEdgeList.begin();
02653 
02654                     error = get_data_from_buff( 1, 1, l + 1, lidx, numedges, localEdgeList, lconn );MB_CHK_ERR( error );
02655                     error = get_data_from_buff( 1, 1, l + 1, ridx, numedges, remEdgeList, rconn );MB_CHK_ERR( error );
02656                     found = check_for_parallelinfo( lconn[0], shared_proc, remProcs );
02657 
02658                     if( !found )
02659                     {
02660                         remProcs.insert( std::pair< EntityHandle, int >( lconn[0], shared_proc ) );
02661                         remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[0], rconn[1] ) );
02662                     }
02663 
02664                     found = check_for_parallelinfo( lconn[1], shared_proc, remProcs );
02665                     if( !found )
02666                     {
02667                         remProcs.insert( std::pair< EntityHandle, int >( lconn[1], shared_proc ) );
02668                         remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[1], rconn[1] ) );
02669                     }
02670                 }
02671             }
02672         }
02673     }
02674 
02675     return MB_SUCCESS;
02676 }
02677 
02678 ErrorCode NestedRefine::decipher_remote_handles_vertex( int shared_proc, int numverts,
02679                                                         std::vector< EntityHandle >& localVertexList,
02680                                                         std::vector< EntityHandle >& remVertexList,
02681                                                         std::multimap< EntityHandle, int >& remProcs,
02682                                                         std::multimap< EntityHandle, EntityHandle >& remHandles )
02683 {
02684     // LVList = <V> where V = <LV0, LV1, ..,LVL>
02685     // RVList = <V> where V = <LV0', RV1, ..,RVL>, LV0' is the local handles of coarsest vertices
02686     // but in the order in which they appear on the remote/shared proc
02687 
02688     ErrorCode error;
02689 
02690     for( int i = 0; i < numverts; i++ )
02691     {
02692         EntityHandle v = localVertexList[i];
02693         int Rvert_idx =
02694             ( std::find( remVertexList.begin(), remVertexList.begin() + numverts - 1, v ) ) - remVertexList.begin();
02695 
02696         std::vector< EntityHandle > lverts, rverts;
02697         for( int l = 0; l < nlevels; l++ )
02698         {
02699             error = get_data_from_buff( 0, 0, l + 1, i, numverts, localVertexList, lverts );MB_CHK_ERR( error );
02700             error = get_data_from_buff( 0, 0, l + 1, Rvert_idx, numverts, remVertexList, rverts );MB_CHK_ERR( error );
02701 
02702             bool found = check_for_parallelinfo( lverts[0], shared_proc, remProcs );
02703             if( !found )
02704             {
02705                 remProcs.insert( std::pair< EntityHandle, int >( lverts[0], shared_proc ) );
02706                 remHandles.insert( std::pair< EntityHandle, EntityHandle >( lverts[0], rverts[0] ) );
02707             }
02708         }
02709     }
02710 
02711     return MB_SUCCESS;
02712 }
02713 
02714 ErrorCode NestedRefine::update_parallel_tags( std::multimap< EntityHandle, int >& remProcs,
02715                                               std::multimap< EntityHandle, EntityHandle >& remHandles )
02716 {
02717     ErrorCode error;
02718 
02719     std::vector< int > rprocs;
02720     std::vector< EntityHandle > rhandles;
02721 
02722     std::multimap< EntityHandle, int >::iterator it;
02723     std::pair< std::multimap< EntityHandle, int >::iterator, std::multimap< EntityHandle, int >::iterator > it_procs;
02724     std::pair< std::multimap< EntityHandle, EntityHandle >::iterator,
02725                std::multimap< EntityHandle, EntityHandle >::iterator >
02726         it_handles;
02727 
02728     it = remProcs.begin();
02729     while( it != remProcs.end() )
02730     {
02731         rprocs.clear();
02732         rhandles.clear();
02733 
02734         EntityHandle entity = it->first;
02735         it_procs            = remProcs.equal_range( entity );
02736         it_handles          = remHandles.equal_range( entity );
02737 
02738         for( std::multimap< EntityHandle, int >::iterator pit = it_procs.first; pit != it_procs.second; pit++ )
02739             rprocs.push_back( pit->second );
02740         for( std::multimap< EntityHandle, EntityHandle >::iterator pit = it_handles.first; pit != it_handles.second;
02741              pit++ )
02742             rhandles.push_back( pit->second );
02743 
02744         error = pcomm->update_remote_data( entity, rprocs, rhandles );MB_CHK_ERR( error );
02745 
02746         it = remProcs.upper_bound( it->first );
02747     }
02748 
02749     return MB_SUCCESS;
02750 }
02751 
02752 ErrorCode NestedRefine::get_data_from_buff( int listtype, int datatype, int level, int entity_index, int nentities,
02753                                             std::vector< EntityHandle >& buffer, std::vector< EntityHandle >& data )
02754 {
02755     /**
02756      * listtype = 2, 1, 0 for FList (faces), EList (edge) and VList (vertices), respectively
02757      *  datatype  =  0(entityhandles of entities at level 0 and their children upto nlevels, in the
02758      *case of vertices its the duplicate vertices at subsequent levels), =  1(connectivity for the
02759      *above), =  2(subentities, only case is edges from FList) entity_index = integer index of the
02760      *entity in the buffer nentities = number of entities at the coarsest level i.e., |F0| or |E0|
02761      *or |V0|
02762      *
02763      * Two types of queries:
02764      * 1) Given an entity, find its children at a particular level.
02765      * 2) Given any entity at any level , find its connectivity (or bounding edges)
02766      *
02767      **/
02768 
02769     data.clear();
02770 
02771     if( listtype == 2 )  // FList
02772     {
02773         // FList = <F, FC, FE> where F = <F0, F1, ..,FL>, FC = <FC0, FC1, .., FCL> and FE = <FE0,
02774         // FE1, .., FEL>
02775         if( datatype == 0 )  // requesting child entities
02776         {
02777             EntityType ftype = mbImpl->type_from_handle( buffer[0] );
02778 
02779             int start, end, toadd, prev;
02780             start = end = entity_index;
02781             toadd = prev = nentities;
02782             for( int i = 0; i < level; i++ )
02783             {
02784                 int d   = get_index_from_degree( level_dsequence[i] );
02785                 int nch = refTemplates[ftype - 1][d].total_new_ents;
02786                 start   = start * nch;
02787                 end     = end * nch + nch - 1;
02788                 if( i < level - 1 )
02789                 {
02790                     prev = prev * nch;
02791                     toadd += prev;
02792                 }
02793             }
02794 
02795             start += toadd;
02796             end += toadd;
02797 
02798             int num_child = end - start + 1;
02799             data.reserve( num_child );
02800 
02801             for( int i = start; i <= end; i++ )
02802             {
02803                 EntityHandle child = buffer[i];
02804                 data.push_back( child );
02805             }
02806         }
02807         else if( datatype == 1 )  // requesting connectivity of an entity
02808         {
02809             EntityType ftype = mbImpl->type_from_handle( buffer[0] );
02810             int nepf         = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
02811 
02812             int toadd = nentities, prev = nentities;
02813             for( int i = 0; i < nlevels; i++ )
02814             {
02815                 int d   = get_index_from_degree( level_dsequence[i] );
02816                 int nch = refTemplates[ftype - 1][d].total_new_ents;
02817                 prev    = prev * nch;
02818                 toadd += prev;
02819             }
02820 
02821             for( int i = 0; i < nepf; i++ )
02822                 data.push_back( buffer[toadd + nepf * entity_index + i] );
02823         }
02824         else if( datatype == 2 )  // requesting bounding edges of an entity
02825         {
02826             EntityType ftype = mbImpl->type_from_handle( buffer[0] );
02827             int nepf         = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
02828 
02829             int toadd = nentities, prev = nentities;
02830             for( int i = 0; i < nlevels; i++ )
02831             {
02832                 int d   = get_index_from_degree( level_dsequence[i] );
02833                 int nch = refTemplates[ftype - 1][d].total_new_ents;
02834                 prev    = prev * nch;
02835                 toadd += prev;
02836             }
02837             toadd += toadd * nepf;
02838 
02839             for( int i = 0; i < nepf; i++ )
02840                 data.push_back( buffer[toadd + nepf * entity_index + i] );
02841         }
02842         else
02843             MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer for faces" );
02844     }
02845     else if( listtype == 1 )  // EList
02846     {
02847         // EList = <E, EC> where E = <E0, E1, ..,EL> and EC = <EC0, EC1, .., ECL>
02848         if( datatype == 0 )  // requesting child entities
02849         {
02850             int start, end, toadd, prev;
02851             start = end = entity_index;
02852             toadd = prev = nentities;
02853             for( int i = 0; i < level; i++ )
02854             {
02855                 int d   = get_index_from_degree( level_dsequence[i] );
02856                 int nch = refTemplates[MBEDGE - 1][d].total_new_ents;
02857                 start   = start * nch;
02858                 end     = end * nch + nch - 1;
02859 
02860                 if( i < level - 1 )
02861                 {
02862                     prev = prev * nch;
02863                     toadd += prev;
02864                 }
02865             }
02866 
02867             start += toadd;
02868             end += toadd;
02869 
02870             int num_child = end - start + 1;
02871             data.reserve( num_child );
02872 
02873             for( int i = start; i <= end; i++ )
02874             {
02875                 EntityHandle child = buffer[i];
02876                 data.push_back( child );
02877             }
02878         }
02879         else if( datatype == 1 )  // requesting connectivities of child entities
02880         {
02881             int toadd = nentities, prev = nentities;
02882             for( int i = 0; i < nlevels; i++ )
02883             {
02884                 int d   = get_index_from_degree( level_dsequence[i] );
02885                 int nch = refTemplates[MBEDGE - 1][d].total_new_ents;
02886                 prev    = prev * nch;
02887                 toadd += prev;
02888             }
02889 
02890             data.push_back( buffer[toadd + 2 * entity_index] );
02891             data.push_back( buffer[toadd + 2 * entity_index + 1] );
02892         }
02893         else
02894             MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer for edges" );
02895     }
02896     else if( listtype == 0 )  // VList
02897     {
02898         // VList = <V> where V = <V0, V1, ..,VL>
02899         if( datatype == 0 )
02900         {
02901             int idx = level * nentities + entity_index;
02902             data.push_back( buffer[idx] );
02903         }
02904         else
02905             MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer for vertices" );
02906     }
02907     else
02908         MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer" );
02909 
02910     return MB_SUCCESS;
02911 }
02912 
02913 bool NestedRefine::check_for_parallelinfo( EntityHandle entity, int proc, std::multimap< EntityHandle, int >& remProcs )
02914 {
02915     bool found = false;
02916 
02917     std::pair< std::multimap< EntityHandle, int >::iterator, std::multimap< EntityHandle, int >::iterator > it_hes;
02918     it_hes = remProcs.equal_range( entity );
02919 
02920     for( std::multimap< EntityHandle, int >::iterator it = it_hes.first; it != it_hes.second; ++it )
02921     {
02922         if( it->second == proc )
02923         {
02924             found = true;
02925             break;
02926         }
02927     }
02928 
02929     return found;
02930 }
02931 
02932 ErrorCode NestedRefine::check_for_parallelinfo( EntityHandle entity, int proc,
02933                                                 std::multimap< EntityHandle, EntityHandle >& remHandles,
02934                                                 std::multimap< EntityHandle, int >& remProcs, EntityHandle& rhandle )
02935 {
02936     // shared procs for given entity
02937     std::pair< std::multimap< EntityHandle, int >::iterator, std::multimap< EntityHandle, int >::iterator > it_ps;
02938     it_ps = remProcs.equal_range( entity );
02939 
02940     // shared remote handles for given entity
02941     std::pair< std::multimap< EntityHandle, EntityHandle >::iterator,
02942                std::multimap< EntityHandle, EntityHandle >::iterator >
02943         it_hs;
02944     it_hs = remHandles.equal_range( entity );
02945 
02946     std::multimap< EntityHandle, int >::iterator itp;
02947     std::multimap< EntityHandle, EntityHandle >::iterator ith;
02948 
02949     for( itp = it_ps.first, ith = it_hs.first; itp != it_ps.second; ++itp, ++ith )
02950     {
02951         if( itp->second == proc )
02952         {
02953             rhandle = ith->second;
02954             break;
02955         }
02956     }
02957 
02958     return MB_SUCCESS;
02959 }
02960 
02961 #endif
02962 
02963 /**********************************
02964  *          Update AHF maps           *
02965  * ********************************/
02966 
02967 ErrorCode NestedRefine::update_local_ahf( int deg, EntityType type, int pat_id, EntityHandle* vbuffer,
02968                                           EntityHandle* ent_buffer, int etotal )
02969 {
02970     ErrorCode error;
02971     int nhf = 0, nv = 0, total_new_verts = 0;
02972     int d = get_index_from_degree( deg );
02973 
02974     // Get the number of half-facets
02975     if( type == MBEDGE )
02976     {
02977         nhf             = 2;
02978         nv              = 2;
02979         total_new_verts = refTemplates[0][d].total_new_verts;
02980     }
02981     else if( type == MBTRI || type == MBQUAD )
02982     {
02983         nhf             = ahf->lConnMap2D[type - 2].num_verts_in_face;
02984         nv              = nhf;
02985         total_new_verts = refTemplates[pat_id][d].total_new_verts;
02986     }
02987     else if( type == MBTET || type == MBHEX )
02988     {
02989         int index       = ahf->get_index_in_lmap( *_incells.begin() );
02990         nhf             = ahf->lConnMap3D[index].num_faces_in_cell;
02991         nv              = ahf->lConnMap3D[index].num_verts_in_cell;
02992         total_new_verts = refTemplates[pat_id][d].total_new_verts;
02993     }
02994 
02995     std::vector< EntityHandle > ent;
02996     std::vector< int > lid;
02997 
02998     // Update the vertex to half-facet map
02999     for( int i = 0; i < total_new_verts; i++ )
03000     {
03001         ent.clear();
03002         lid.clear();
03003         EntityHandle vid = vbuffer[i + nv];
03004         error            = ahf->get_incident_map( type, vid, ent, lid );MB_CHK_ERR( error );
03005 
03006         if( ent[0] ) continue;
03007 
03008         int id = refTemplates[pat_id][d].v2hf[i + nv][0] - 1;
03009         ent[0] = ent_buffer[id];
03010         lid[0] = refTemplates[pat_id][d].v2hf[i + nv][1];
03011 
03012         error = ahf->set_incident_map( type, vid, ent, lid );MB_CHK_ERR( error );
03013     }
03014 
03015     // Update the sibling half-facet map
03016     for( int i = 0; i < etotal; i++ )
03017     {
03018         std::vector< EntityHandle > sib_entids( nhf );
03019         std::vector< int > sib_lids( nhf );
03020 
03021         error = ahf->get_sibling_map( type, ent_buffer[i], &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
03022 
03023         for( int l = 0; l < nhf; l++ )
03024         {
03025             if( sib_entids[l] ) continue;
03026 
03027             // Fill out the sibling values
03028             int id = refTemplates[pat_id][d].ents_opphfs[i][2 * l];
03029             if( id )
03030             {
03031                 sib_entids[l] = ent_buffer[id - 1];
03032                 sib_lids[l]   = refTemplates[pat_id][d].ents_opphfs[i][2 * l + 1];
03033             }
03034             else
03035             {
03036                 sib_entids[l] = 0;
03037                 sib_lids[l]   = 0;
03038             }
03039         }
03040 
03041         error = ahf->set_sibling_map( type, ent_buffer[i], &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
03042 
03043         for( int l = 0; l < nhf; l++ )
03044         {
03045             if( sib_entids[l] )
03046             {
03047 
03048                 EntityHandle set_entid = ent_buffer[i];
03049                 int set_lid            = l;
03050 
03051                 error = ahf->set_sibling_map( type, sib_entids[l], sib_lids[l], set_entid, set_lid );MB_CHK_ERR( error );
03052             }
03053         }
03054     }
03055     return MB_SUCCESS;
03056 }
03057 
03058 ErrorCode NestedRefine::update_local_ahf( int deg, EntityType type, EntityHandle* vbuffer, EntityHandle* ent_buffer,
03059                                           int etotal )
03060 {
03061     ErrorCode error;
03062     assert( type != MBTET );
03063     error = update_local_ahf( deg, type, type - 1, vbuffer, ent_buffer, etotal );MB_CHK_ERR( error );
03064 
03065     return MB_SUCCESS;
03066 }
03067 
03068 ErrorCode NestedRefine::update_global_ahf( EntityType type, int cur_level, int deg, std::vector< int >* pattern_ids )
03069 {
03070 
03071     ErrorCode error;
03072 
03073     // Get the number of half-facets and number of children of each type
03074     if( type == MBEDGE )
03075     {
03076         assert( pattern_ids == NULL );
03077         error = update_global_ahf_1D( cur_level, deg );MB_CHK_ERR( error );
03078     }
03079     else if( type == MBTRI || type == MBQUAD )
03080     {
03081         assert( pattern_ids == NULL );
03082         error = update_global_ahf_2D( cur_level, deg );MB_CHK_ERR( error );
03083     }
03084     else if( type == MBHEX )
03085     {
03086         assert( pattern_ids == NULL );
03087         error = update_global_ahf_3D( cur_level, deg );MB_CHK_ERR( error );
03088     }
03089     else if( type == MBTET )
03090     {
03091         assert( pattern_ids != NULL );
03092         error = update_global_ahf_3D( cur_level, deg, pattern_ids );MB_CHK_ERR( error );
03093     }
03094     else
03095         MB_SET_ERR( MB_NOT_IMPLEMENTED, "Requesting AHF update for an unsupported mesh entity type" );
03096 
03097     return MB_SUCCESS;
03098 }
03099 
03100 /*ErrorCode NestedRefine::update_global_ahf(int cur_level, int deg, std::vector<int> &pattern_ids)
03101 {
03102   ErrorCode error;
03103   error = update_global_ahf_3D(cur_level, deg, pattern_ids); MB_CHK_ERR(error);
03104 
03105   return MB_SUCCESS;
03106 }*/
03107 
03108 ErrorCode NestedRefine::update_global_ahf_1D( int cur_level, int deg )
03109 {
03110     ErrorCode error;
03111     int d = get_index_from_degree( deg );
03112     int nhf, nchilds, nverts_prev, nents_prev;
03113     nhf     = 2;
03114     nchilds = refTemplates[0][d].total_new_ents;
03115     if( cur_level )
03116     {
03117         nverts_prev = level_mesh[cur_level - 1].num_verts;
03118         nents_prev  = level_mesh[cur_level - 1].num_edges;
03119     }
03120     else
03121     {
03122         nverts_prev = _inverts.size();
03123         nents_prev  = _inedges.size();
03124     }
03125 
03126     std::vector< EntityHandle > inci_ent, child_ents;
03127     std::vector< int > inci_lid, child_lids;
03128 
03129     // Update the vertex to half-facet maps for duplicate vertices
03130     for( int i = 0; i < nverts_prev; i++ )
03131     {
03132         inci_ent.clear();
03133         inci_lid.clear();
03134         child_ents.clear();
03135         child_lids.clear();
03136 
03137         // Vertex id in the previous mesh and the current one
03138         EntityHandle vid;
03139         if( cur_level )
03140             vid = level_mesh[cur_level - 1].start_vertex + i;
03141         else
03142             vid = _inverts[i];
03143         EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
03144 
03145         // Get the incident half-vert in the previous mesh
03146         error = ahf->get_incident_map( MBEDGE, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
03147 
03148         // Obtain the corresponding incident child in the current mesh
03149         int lvid = get_local_vid( vid, inci_ent[0], cur_level - 1 );
03150         if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
03151         int chid = refTemplates[0][d].v2hf[lvid][0] - 1;
03152 
03153         int pid;
03154         if( cur_level )
03155             pid = inci_ent[0] - level_mesh[cur_level - 1].start_edge;
03156         else
03157             pid = inci_ent[0] - *_inedges.begin();
03158 
03159         int ind = nchilds * pid;
03160 
03161         child_ents.push_back( level_mesh[cur_level].start_edge + ind + chid );
03162         child_lids.push_back( refTemplates[0][d].v2hf[lvid][1] );
03163 
03164         error = ahf->set_incident_map( MBEDGE, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
03165     }
03166 
03167     // Update the sibling half-facet maps across entities
03168     for( int i = 0; i < nents_prev; i++ )
03169     {
03170         EntityHandle ent;
03171         if( cur_level )
03172             ent = level_mesh[cur_level - 1].start_edge + i;
03173         else
03174             ent = _inedges[i];
03175 
03176         std::vector< EntityHandle > sib_entids( nhf );
03177         std::vector< int > sib_lids( nhf );
03178 
03179         error = ahf->get_sibling_map( MBEDGE, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
03180 
03181         int id, idx;
03182 
03183         for( int l = 0; l < nhf; l++ )
03184         {
03185             if( !sib_entids[l] ) continue;
03186 
03187             // Find the child incident on the half-facet
03188             id                     = refTemplates[0][d].ents_on_pent[l][1] - 1;
03189             idx                    = nchilds * i;
03190             EntityHandle child_ent = level_mesh[cur_level].start_edge + idx + id;
03191             int ch_lid             = l;
03192 
03193             // Find the sibling of the child
03194             std::vector< EntityHandle > sib_childs( nhf );
03195             std::vector< int > sib_chlids( nhf );
03196 
03197             error = ahf->get_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
03198 
03199             // If the sibling already exists, dont do anything
03200             if( sib_childs[ch_lid] ) continue;
03201 
03202             // Get the correponding child of the sibling of the current parent
03203             int psib;
03204             if( cur_level )
03205                 psib = sib_entids[l] - level_mesh[cur_level - 1].start_edge;
03206             else
03207                 psib = sib_entids[l] - *_inedges.begin();
03208 
03209             int plid = sib_lids[l];
03210 
03211             id  = refTemplates[0][d].ents_on_pent[plid][1] - 1;
03212             idx = nchilds * psib;
03213 
03214             EntityHandle psib_child = level_mesh[cur_level].start_edge + idx + id;
03215             int psib_chlid          = plid;
03216 
03217             // Set the siblings
03218             sib_childs[ch_lid] = psib_child;
03219             sib_chlids[ch_lid] = psib_chlid;
03220 
03221             error = ahf->set_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
03222         }
03223     }
03224 
03225     return MB_SUCCESS;
03226 }
03227 
03228 ErrorCode NestedRefine::update_global_ahf_1D_sub( int cur_level, int deg )
03229 {
03230     ErrorCode error;
03231     int d = get_index_from_degree( deg );
03232     int nhf, nchilds, nents_prev;
03233     nhf     = 2;
03234     nchilds = refTemplates[0][d].total_new_ents;
03235     if( cur_level ) { nents_prev = level_mesh[cur_level - 1].num_edges; }
03236     else
03237     {
03238         nents_prev = _inedges.size();
03239     }
03240 
03241     // Update the sibling half-facet maps across entities
03242 
03243     std::vector< EntityHandle > conn;
03244     for( int i = 0; i < nents_prev; i++ )
03245     {
03246         EntityHandle ent;
03247         if( cur_level )
03248             ent = level_mesh[cur_level - 1].start_edge + i;
03249         else
03250             ent = _inedges[i];
03251 
03252         // Set incident hv maps
03253         conn.clear();
03254         error = get_connectivity( ent, cur_level, conn );MB_CHK_ERR( error );
03255 
03256         std::vector< EntityHandle > inci_ent, child_ents;
03257         std::vector< int > inci_lid, child_lids;
03258         for( int j = 0; j < 2; j++ )
03259         {
03260             inci_ent.clear();
03261             inci_lid.clear();
03262             child_ents.clear();
03263             child_lids.clear();
03264 
03265             // Get the entityhandle of the vertex from previous level in the current level
03266             EntityHandle cur_vid;
03267             if( cur_level )
03268                 cur_vid = level_mesh[cur_level].start_vertex + ( conn[j] - level_mesh[cur_level - 1].start_vertex );
03269             else
03270                 cur_vid = level_mesh[cur_level].start_vertex + ( conn[j] - *_inverts.begin() );
03271 
03272             // Obtain the incident half-facet. If exists, then no need to assign another
03273             error = ahf->get_incident_map( MBEDGE, cur_vid, inci_ent, inci_lid );MB_CHK_ERR( error );
03274             if( inci_ent[0] != 0 ) continue;
03275 
03276             // Get the incident half-facet on the old vertex
03277             error = ahf->get_incident_map( MBEDGE, conn[j], inci_ent, inci_lid );MB_CHK_ERR( error );
03278 
03279             // Obtain the corresponding incident child in the current mesh
03280             int lvid = get_local_vid( conn[j], inci_ent[0], cur_level - 1 );
03281             if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
03282             int chid = refTemplates[0][d].v2hf[lvid][0] - 1;
03283 
03284             int pid;
03285             if( cur_level )
03286                 pid = inci_ent[0] - level_mesh[cur_level - 1].start_edge;
03287             else
03288                 pid = inci_ent[0] - *_inedges.begin();
03289 
03290             int ind = nchilds * pid;
03291 
03292             child_ents.push_back( level_mesh[cur_level].start_edge + ind + chid );
03293             child_lids.push_back( refTemplates[0][d].v2hf[lvid][1] );
03294 
03295             error = ahf->set_incident_map( MBEDGE, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
03296         }
03297 
03298         std::vector< EntityHandle > sib_entids( nhf );
03299         std::vector< int > sib_lids( nhf );
03300 
03301         error = ahf->get_sibling_map( MBEDGE, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
03302 
03303         int id, idx;
03304 
03305         for( int l = 0; l < nhf; l++ )
03306         {
03307             if( !sib_entids[l] ) continue;
03308 
03309             // Find the child incident on the half-facet
03310             id                     = refTemplates[0][d].ents_on_pent[l][1] - 1;
03311             idx                    = nchilds * i;
03312             EntityHandle child_ent = level_mesh[cur_level].start_edge + idx + id;
03313             int ch_lid             = l;
03314 
03315             // Find the sibling of the child
03316             std::vector< EntityHandle > sib_childs( nhf );
03317             std::vector< int > sib_chlids( nhf );
03318 
03319             error = ahf->get_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
03320 
03321             // If the sibling already exists, dont do anything
03322             if( sib_childs[ch_lid] ) continue;
03323 
03324             // Get the correponding child of the sibling of the current parent
03325             int psib;
03326             if( cur_level )
03327                 psib = sib_entids[l] - level_mesh[cur_level - 1].start_edge;
03328             else
03329                 psib = sib_entids[l] - *_inedges.begin();
03330 
03331             int plid = sib_lids[l];
03332 
03333             id  = refTemplates[0][d].ents_on_pent[plid][1] - 1;
03334             idx = nchilds * psib;
03335 
03336             EntityHandle psib_child = level_mesh[cur_level].start_edge + idx + id;
03337             int psib_chlid          = plid;
03338 
03339             // Set the siblings
03340             sib_childs[ch_lid] = psib_child;
03341             sib_chlids[ch_lid] = psib_chlid;
03342 
03343             error = ahf->set_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
03344         }
03345     }
03346 
03347     return MB_SUCCESS;
03348 }
03349 
03350 ErrorCode NestedRefine::update_ahf_1D( int cur_level )
03351 {
03352     ErrorCode error;
03353     error = ahf->determine_sibling_halfverts( level_mesh[cur_level].verts, level_mesh[cur_level].edges );MB_CHK_ERR( error );
03354 
03355     error = ahf->determine_incident_halfverts( level_mesh[cur_level].edges );MB_CHK_ERR( error );
03356 
03357     return MB_SUCCESS;
03358 }
03359 
03360 ErrorCode NestedRefine::update_global_ahf_2D( int cur_level, int deg )
03361 {
03362     ErrorCode error;
03363 
03364     EntityType type = mbImpl->type_from_handle( *_infaces.begin() );
03365     int nhf, nchilds, nverts_prev, nents_prev;
03366 
03367     nhf     = ahf->lConnMap2D[type - 2].num_verts_in_face;
03368     int d   = get_index_from_degree( deg );
03369     nchilds = refTemplates[type - 1][d].total_new_ents;
03370 
03371     if( cur_level )
03372     {
03373         nverts_prev = level_mesh[cur_level - 1].num_verts;
03374         nents_prev  = level_mesh[cur_level - 1].num_faces;
03375     }
03376     else
03377     {
03378         nverts_prev = _inverts.size();
03379         nents_prev  = _infaces.size();
03380     }
03381 
03382     std::vector< EntityHandle > inci_ent, child_ents;
03383     std::vector< int > inci_lid, child_lids;
03384 
03385     // Update the vertex to half-edge maps for old/duplicate vertices
03386     for( int i = 0; i < nverts_prev; i++ )
03387     {
03388         inci_ent.clear();
03389         inci_lid.clear();
03390         child_ents.clear();
03391         child_lids.clear();
03392 
03393         // Vertex id in the previous mesh
03394         EntityHandle vid;
03395         if( cur_level )
03396             vid = level_mesh[cur_level - 1].start_vertex + i;
03397         else
03398             vid = _inverts[i];
03399         EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
03400 
03401         // Get the incident half-vert in the previous mesh
03402         error = ahf->get_incident_map( type, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
03403 
03404         // Obtain the corresponding incident child in the current mesh
03405         for( int j = 0; j < (int)inci_ent.size(); j++ )
03406         {
03407             int lvid = get_local_vid( vid, inci_ent[j], cur_level - 1 );
03408             if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
03409             int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
03410 
03411             int pid;
03412             if( cur_level )
03413                 pid = inci_ent[j] - level_mesh[cur_level - 1].start_face;
03414             else
03415                 pid = inci_ent[j] - *_infaces.begin();
03416 
03417             int ind = nchilds * pid;
03418 
03419             child_ents.push_back( level_mesh[cur_level].start_face + ind + chid );
03420             child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
03421         }
03422         error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
03423     }
03424 
03425     EntityHandle fedge[2];
03426 
03427     // Update the sibling half-facet maps across entities
03428     for( int i = 0; i < nents_prev; i++ )
03429     {
03430         EntityHandle ent;
03431         if( cur_level )
03432             ent = level_mesh[cur_level - 1].start_face + i;
03433         else
03434             ent = _infaces[i];
03435 
03436         std::vector< EntityHandle > fid_conn;
03437         error = get_connectivity( ent, cur_level, fid_conn );
03438         if( MB_SUCCESS != error ) return error;
03439 
03440         std::vector< EntityHandle > sib_entids( nhf );
03441         std::vector< int > sib_lids( nhf );
03442 
03443         error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
03444 
03445         int id, idx;
03446 
03447         for( int l = 0; l < nhf; l++ )
03448         {
03449             if( !sib_entids[l] ) continue;
03450 
03451             int nidx = ahf->lConnMap2D[type - 2].next[l];
03452             fedge[0] = fid_conn[l];
03453             fedge[1] = fid_conn[nidx];
03454 
03455             EntityHandle sfid = sib_entids[l];
03456             int slid          = sib_lids[l];
03457 
03458             std::vector< EntityHandle > conn;
03459             error = get_connectivity( sfid, cur_level, conn );
03460             if( MB_SUCCESS != error ) return error;
03461 
03462             bool orient = true;
03463             nidx        = ahf->lConnMap2D[type - 2].next[slid];
03464             if( ( fedge[1] == conn[slid] ) && ( fedge[0] == conn[nidx] ) ) orient = false;
03465 
03466             if( orient ) assert( ( fedge[0] == conn[slid] ) && ( fedge[1] == conn[nidx] ) );
03467 
03468             // Find the childrens incident on the half-facet
03469             int nch = refTemplates[type - 1][d].ents_on_pent[l][0];
03470             idx     = nchilds * i;
03471 
03472             // Loop over all the incident childrens
03473             for( int k = 0; k < nch; k++ )
03474             {
03475                 id                     = refTemplates[type - 1][d].ents_on_pent[l][k + 1] - 1;
03476                 EntityHandle child_ent = level_mesh[cur_level].start_face + idx + id;
03477                 int child_lid          = l;
03478 
03479                 // Find the sibling of the child
03480                 EntityHandle child_sibent;
03481                 int child_siblid;
03482                 error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
03483 
03484                 if( child_sibent != 0 ) continue;
03485 
03486                 // Get the correponding child of the sibling of the current parent
03487                 int psib;
03488                 if( cur_level )
03489                     psib = sfid - level_mesh[cur_level - 1].start_face;
03490                 else
03491                     psib = sfid - *_infaces.begin();
03492 
03493                 int plid = slid;
03494 
03495                 if( orient )
03496                     id = refTemplates[type - 1][d].ents_on_pent[plid][k + 1] - 1;
03497                 else
03498                     id = refTemplates[type - 1][d].ents_on_pent[plid][nch - k] - 1;
03499 
03500                 int sidx = nchilds * psib;
03501 
03502                 EntityHandle psib_child = level_mesh[cur_level].start_face + sidx + id;
03503                 int psib_chlid          = plid;
03504 
03505                 // Set the siblings
03506                 error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
03507             }
03508         }
03509     }
03510 
03511     return MB_SUCCESS;
03512 }
03513 
03514 ErrorCode NestedRefine::update_global_ahf_2D_sub( int cur_level, int deg )
03515 {
03516     ErrorCode error;
03517     int d           = get_index_from_degree( deg );
03518     EntityType type = mbImpl->type_from_handle( *_infaces.begin() );
03519     int nhf, nchilds, nents_prev;
03520     nhf     = ahf->lConnMap2D[type - 2].num_verts_in_face;
03521     nchilds = refTemplates[type - 1][d].total_new_ents;
03522 
03523     if( cur_level )
03524         nents_prev = level_mesh[cur_level - 1].num_faces;
03525     else
03526         nents_prev = _infaces.size();
03527 
03528     EntityHandle fedge[2];
03529 
03530     // Update the sibling half-facet maps across entities
03531     for( int i = 0; i < nents_prev; i++ )
03532     {
03533         EntityHandle ent;
03534         if( cur_level )
03535             ent = level_mesh[cur_level - 1].start_face + i;
03536         else
03537             ent = _infaces[i];
03538 
03539         std::vector< EntityHandle > fid_conn;
03540         error = get_connectivity( ent, cur_level, fid_conn );
03541         if( MB_SUCCESS != error ) return error;
03542 
03543         std::vector< EntityHandle > inci_ent, child_ents;
03544         std::vector< int > inci_lid, child_lids;
03545 
03546         // Set incident half-edges
03547         for( int j = 0; j < nhf; j++ )
03548         {
03549             inci_ent.clear();
03550             inci_lid.clear();
03551             child_ents.clear();
03552             child_lids.clear();
03553             EntityHandle cur_vid;
03554             if( cur_level )
03555                 cur_vid = level_mesh[cur_level].start_vertex + ( fid_conn[j] - level_mesh[cur_level - 1].start_vertex );
03556             else
03557                 cur_vid = level_mesh[cur_level].start_vertex + ( fid_conn[j] - *_inverts.begin() );
03558 
03559             // Obtain the incident half-facet. If exists, then no need to assign another
03560             error = ahf->get_incident_map( type, cur_vid, inci_ent, inci_lid );MB_CHK_ERR( error );
03561             if( inci_ent[0] != 0 ) continue;
03562 
03563             // Get the incident half-facet on the old vertex
03564             error = ahf->get_incident_map( type, fid_conn[j], inci_ent, inci_lid );MB_CHK_ERR( error );
03565 
03566             // Obtain the corresponding incident child in the current mesh
03567             for( int k = 0; k < (int)inci_ent.size(); k++ )
03568             {
03569                 int lvid = get_local_vid( fid_conn[j], inci_ent[k], cur_level - 1 );
03570                 if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
03571                 int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
03572 
03573                 int pid;
03574                 if( cur_level )
03575                     pid = inci_ent[k] - level_mesh[cur_level - 1].start_face;
03576                 else
03577                     pid = inci_ent[k] - *_infaces.begin();
03578 
03579                 int ind = nchilds * pid;
03580 
03581                 child_ents.push_back( level_mesh[cur_level].start_face + ind + chid );
03582                 child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
03583             }
03584 
03585             error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
03586         }
03587 
03588         // Set sibling half-edges
03589         std::vector< EntityHandle > sib_entids( nhf );
03590         std::vector< int > sib_lids( nhf );
03591 
03592         error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
03593 
03594         int id, idx;
03595 
03596         for( int l = 0; l < nhf; l++ )
03597         {
03598             if( !sib_entids[l] ) continue;
03599 
03600             int nidx = ahf->lConnMap2D[type - 2].next[l];
03601             fedge[0] = fid_conn[l];
03602             fedge[1] = fid_conn[nidx];
03603 
03604             EntityHandle sfid = sib_entids[l];
03605             int slid          = sib_lids[l];
03606 
03607             std::vector< EntityHandle > conn;
03608             error = get_connectivity( sfid, cur_level, conn );MB_CHK_ERR( error );
03609 
03610             assert( (int)conn.size() > nidx && (int)conn.size() > slid );
03611 
03612             bool orient = true;
03613             nidx        = ahf->lConnMap2D[type - 2].next[slid];
03614             if( ( fedge[1] == conn[slid] ) && ( fedge[0] == conn[nidx] ) ) orient = false;
03615 
03616             if( orient ) assert( ( fedge[0] == conn[slid] ) && ( fedge[1] == conn[nidx] ) );
03617 
03618             // Find the childrens incident on the half-facet
03619             int nch = refTemplates[type - 1][d].ents_on_pent[l][0];
03620             idx     = nchilds * i;
03621 
03622             // Loop over all the incident childrens
03623             for( int k = 0; k < nch; k++ )
03624             {
03625                 id                     = refTemplates[type - 1][d].ents_on_pent[l][k + 1] - 1;
03626                 EntityHandle child_ent = level_mesh[cur_level].start_face + idx + id;
03627                 int child_lid          = l;
03628 
03629                 // Find the sibling of the child
03630                 EntityHandle child_sibent;
03631                 int child_siblid;
03632                 error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
03633 
03634                 if( child_sibent != 0 ) continue;
03635 
03636                 // Get the correponding child of the sibling of the current parent
03637                 int psib;
03638                 if( cur_level )
03639                     psib = sfid - level_mesh[cur_level - 1].start_face;
03640                 else
03641                     psib = sfid - *_infaces.begin();
03642 
03643                 int plid = slid;
03644 
03645                 if( orient )
03646                     id = refTemplates[type - 1][d].ents_on_pent[plid][k + 1] - 1;
03647                 else
03648                     id = refTemplates[type - 1][d].ents_on_pent[plid][nch - k] - 1;
03649 
03650                 int sidx = nchilds * psib;
03651 
03652                 EntityHandle psib_child = level_mesh[cur_level].start_face + sidx + id;
03653                 int psib_chlid          = plid;
03654 
03655                 // Set the siblings
03656                 error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
03657             }
03658         }
03659     }
03660 
03661     return MB_SUCCESS;
03662 }
03663 
03664 ErrorCode NestedRefine::update_global_ahf_3D( int cur_level, int deg, std::vector< int >* pattern_ids )
03665 {
03666     ErrorCode error;
03667     int nvpc, ne, nhf, nchilds, nverts_prev, nents_prev;
03668 
03669     EntityType type = mbImpl->type_from_handle( *_incells.begin() );
03670     int index       = ahf->get_index_in_lmap( *_incells.begin() );
03671     int d           = get_index_from_degree( deg );
03672 
03673     nhf     = ahf->lConnMap3D[index].num_faces_in_cell;
03674     ne      = ahf->lConnMap3D[index].num_edges_in_cell;
03675     nvpc    = ahf->lConnMap3D[index].num_verts_in_cell;
03676     nchilds = refTemplates[type - 1][d].total_new_ents;
03677 
03678     if( cur_level )
03679     {
03680         nverts_prev = level_mesh[cur_level - 1].num_verts;
03681         nents_prev  = level_mesh[cur_level - 1].num_cells;
03682     }
03683     else
03684     {
03685         nverts_prev = _inverts.size();
03686         nents_prev  = _incells.size();
03687     }
03688 
03689     std::vector< EntityHandle > inci_ent, child_ents;
03690     std::vector< int > inci_lid, child_lids;
03691 
03692     // Step 1: Update the V2HF maps for old/duplicate vertices
03693     for( int i = 0; i < nverts_prev; i++ )
03694     {
03695         inci_ent.clear();
03696         inci_lid.clear();
03697         child_ents.clear();
03698         child_lids.clear();
03699 
03700         // Vertex id in the previous mesh
03701         EntityHandle vid;
03702         if( cur_level )
03703             vid = level_mesh[cur_level - 1].start_vertex + i;
03704         else
03705             vid = _inverts[i];
03706         EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
03707 
03708         // Get the incident half-vert in the previous mesh
03709         error = ahf->get_incident_map( type, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
03710 
03711         // Obtain the corresponding incident child in the current mesh
03712         for( int j = 0; j < (int)inci_ent.size(); j++ )
03713         {
03714             int lvid = get_local_vid( vid, inci_ent[j], cur_level - 1 );
03715             if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
03716             int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
03717 
03718             int pid;
03719             if( cur_level )
03720                 pid = inci_ent[j] - level_mesh[cur_level - 1].start_cell;
03721             else
03722                 pid = inci_ent[j] - *_incells.begin();
03723 
03724             int ind = nchilds * pid;
03725 
03726             //  EntityHandle child_ent = level_mesh[cur_level].start_cell + ind+chid ;
03727             // int child_lid = refTemplates[type-1][d].v2hf[lvid][1];
03728             child_ents.push_back( level_mesh[cur_level].start_cell + ind + chid );
03729             child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
03730         }
03731 
03732         error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
03733     }
03734 
03735     //  error = ahf->determine_incident_halffaces( level_mesh[cur_level].cells);MB_CHK_ERR(error);
03736 
03737     // Step 2: Update SIBHFS maps
03738     for( int i = 0; i < nents_prev; i++ )
03739     {
03740         EntityHandle ent;
03741         if( cur_level )
03742             ent = level_mesh[cur_level - 1].start_cell + i;
03743         else
03744             ent = _incells[i];
03745 
03746         std::vector< EntityHandle > sib_entids( nhf );
03747         std::vector< int > sib_lids( nhf );
03748 
03749         error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
03750 
03751         int id, idx;
03752 
03753         for( int l = 0; l < nhf; l++ )
03754         {
03755 
03756             if( !sib_entids[l] ) continue;
03757 
03758             // Get the number of children incident on this half-face
03759             int pat_id;
03760             if( type == MBTET )
03761                 pat_id = ( *pattern_ids )[i];
03762             else
03763                 pat_id = type - 1;
03764             int nch = refTemplates[pat_id][d].ents_on_pent[l][0];
03765 
03766             // Get the order of children indices incident on this half-face
03767             std::vector< int > id_sib( nch );
03768             for( int k = 0; k < nch; k++ )
03769                 id_sib[k] = 0;
03770 
03771             error = reorder_indices( cur_level, deg, ent, l, sib_entids[l], sib_lids[l], 1, &id_sib[0] );MB_CHK_ERR( error );
03772 
03773             // Get the parent index of the sibling cell
03774             int psib;
03775             if( cur_level )
03776                 psib = sib_entids[l] - level_mesh[cur_level - 1].start_cell;
03777             else
03778                 psib = sib_entids[l] - *_incells.begin();
03779 
03780             int plid = sib_lids[l];
03781             int sidx = nchilds * psib;
03782             int sibpat_id;
03783             if( type == MBTET )
03784                 sibpat_id = ( *pattern_ids )[psib];
03785             else
03786                 sibpat_id = type - 1;
03787 
03788             // Loop over all the childs incident on the working half-face
03789             idx = nchilds * i;
03790 
03791             for( int k = 0; k < nch; k++ )
03792             {
03793                 id                     = refTemplates[pat_id][d].ents_on_pent[l][k + 1] - 1;
03794                 EntityHandle child_ent = level_mesh[cur_level].start_cell + idx + id;
03795                 int child_lid          = l;
03796 
03797                 // Find the sibling of the working child
03798                 EntityHandle child_sibent;
03799                 int child_siblid;
03800                 error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
03801 
03802                 if( child_sibent != 0 ) continue;
03803 
03804                 // Get the correponding child of the sibling of the current parent
03805                 // We have already computed the order the children on incident corresponding to the
03806                 // working half-face
03807                 id = refTemplates[sibpat_id][d].ents_on_pent[plid][id_sib[k]] - 1;
03808 
03809                 EntityHandle psib_child = level_mesh[cur_level].start_cell + sidx + id;
03810                 int psib_chlid          = plid;
03811 
03812                 // Set the siblings of children incident on current half-face
03813                 error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
03814 
03815                 // Set the sibling of the sibling of the children to the children
03816                 error = ahf->set_sibling_map( type, psib_child, psib_chlid, child_ent, child_lid );MB_CHK_ERR( error );
03817             }
03818         }
03819 
03820         // Loop over edges to check if there are any non-manifold edges. If there are then the v2hfs
03821         // map should be updated for the new vertices on it.
03822         const EntityHandle* conn;
03823         error = mbImpl->get_connectivity( ent, conn, nvpc );MB_CHK_ERR( error );
03824 
03825         int nv = refTemplates[type - 1][d].nv_edge;  //#verts on each edge
03826         for( int l = 0; l < ne; l++ )
03827         {
03828             id                   = ahf->lConnMap3D[index].e2v[l][0];
03829             EntityHandle v_start = conn[id];
03830             id                   = ahf->lConnMap3D[index].e2v[l][1];
03831             EntityHandle v_end   = conn[id];
03832 
03833             bool visited = false;
03834 
03835             std::vector< EntityHandle > inci_ent1, inci_ent2;
03836             std::vector< int > inci_lid1, inci_lid2;
03837             error = ahf->get_incident_map( type, v_start, inci_ent1, inci_lid1 );MB_CHK_ERR( error );
03838             error = ahf->get_incident_map( type, v_end, inci_ent2, inci_lid2 );MB_CHK_ERR( error );
03839 
03840             if( inci_ent1.size() > 1 && inci_ent2.size() > 1 )
03841             {
03842                 std::vector< EntityHandle > cell_comps;
03843                 std::vector< int > leid_comps;
03844 
03845                 error = ahf->get_up_adjacencies_edg_3d_comp( ent, l, cell_comps, &leid_comps );MB_CHK_ERR( error );
03846 
03847                 int ncomps = cell_comps.size();
03848                 std::vector< EntityHandle > edgverts;
03849                 std::vector< EntityHandle > compchildents( nv * ncomps );
03850                 std::vector< int > compchildlfids( nv * ncomps );
03851 
03852                 for( int s = 0; s < nv * ncomps; s++ )
03853                 {
03854                     compchildents[s]  = 0;
03855                     compchildlfids[s] = 0;
03856                 }
03857 
03858                 for( int j = 0; j < (int)cell_comps.size(); j++ )
03859                 {
03860                     int ind;
03861                     if( cur_level )
03862                         ind = level_mesh[cur_level - 1].cells.index( cell_comps[j] );
03863                     else
03864                         ind = _incells.index( cell_comps[j] );
03865 
03866                     for( int k = 0; k < nv; k++ )
03867                     {
03868                         int chid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k] - 1;
03869                         int lfid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k + 1];
03870                         int lvid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k + 2];
03871 
03872                         EntityHandle childcell = level_mesh[cur_level].start_cell + ind * nchilds + chid;
03873 
03874                         const EntityHandle* econn;
03875                         error = mbImpl->get_connectivity( childcell, econn, nvpc );MB_CHK_ERR( error );
03876 
03877                         EntityHandle vert = econn[lvid];
03878 
03879                         if( ahf->check_nonmanifold_vertices( type, vert ) )
03880                         {
03881                             visited = true;
03882                             break;
03883                         }
03884 
03885                         if( edgverts.empty() )
03886                         {
03887                             edgverts.push_back( vert );
03888                             compchildents[0]  = childcell;
03889                             compchildlfids[0] = lfid;
03890                         }
03891                         else
03892                         {
03893                             std::vector< EntityHandle >::iterator it;
03894                             it       = find( edgverts.begin(), edgverts.end(), vert );
03895                             int indx = it - edgverts.begin();
03896 
03897                             if( it == edgverts.end() )
03898                             {
03899                                 edgverts.push_back( vert );
03900                                 compchildents[k * ncomps]  = childcell;
03901                                 compchildlfids[k * ncomps] = lfid;
03902                             }
03903                             else
03904                             {
03905                                 compchildents[indx * ncomps + j]  = childcell;
03906                                 compchildlfids[indx * ncomps + j] = lfid;
03907                             }
03908                         }
03909                     }
03910                 }
03911 
03912                 if( visited ) { break; }
03913 
03914                 // Set the incident half-facet map
03915                 for( int k = 0; k < nv; k++ )
03916                 {
03917                     std::vector< EntityHandle > set_childents;
03918                     std::vector< int > set_childlfids;
03919                     for( int j = 0; j < ncomps; j++ )
03920                     {
03921                         set_childents.push_back( compchildents[k * ncomps + j] );
03922                         set_childlfids.push_back( compchildlfids[k * ncomps + j] );
03923                     }
03924 
03925                     error = ahf->set_incident_map( type, edgverts[k], set_childents, set_childlfids );MB_CHK_ERR( error );
03926                 }
03927             }
03928         }
03929     }
03930 
03931     return MB_SUCCESS;
03932 }
03933 
03934 ErrorCode NestedRefine::get_lid_inci_child( EntityType type, int deg, int lfid, int leid, std::vector< int >& child_ids,
03935                                             std::vector< int >& child_lvids )
03936 {
03937     int index = ahf->get_index_in_lmap( *_incells.begin() );
03938     int d     = get_index_from_degree( deg );
03939 
03940     // int lv0 = ahf->lConnMap3D[index].e2v[leid][0];
03941     //  int lv1 = ahf->lConnMap3D[index].e2v[leid][1];
03942     int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
03943 
03944     int nv  = refTemplates[type - 1][d].nv_edge;
03945     int nch = refTemplates[type - 1][d].ents_on_pent[lfid][0];
03946 
03947     for( int i = 0; i < nch; i++ )
03948     {
03949         int id = refTemplates[type - 1][d].ents_on_pent[lfid][i + 1] - 1;
03950         for( int j = 0; j < nvpc; j++ )
03951         {
03952             int lv = refTemplates[type - 1][d].ents_conn[id][j];
03953             for( int k = 0; k < nv; k++ )
03954             {
03955                 if( lv == refTemplates[type - 1][d].vert_on_edges[leid][k] )
03956                 {
03957                     child_ids.push_back( id );
03958                     child_lvids.push_back( j );
03959                 }
03960             }
03961         }
03962     }
03963 
03964     return MB_SUCCESS;
03965 }
03966 
03967 /* **********************************
03968  *  *          Boundary Functions      *
03969  ************************************/
03970 
03971 bool NestedRefine::is_vertex_on_boundary( const EntityHandle& vertex )
03972 {
03973     ErrorCode error;
03974     EntityHandle sibents[27];
03975     int siblids[27];
03976     std::vector< EntityHandle > ent;
03977     std::vector< int > lid;
03978 
03979     int nhf;
03980     if( elementype == MBEDGE )
03981         nhf = 2;
03982     else if( ( elementype == MBTRI ) || ( elementype == MBQUAD ) )
03983         nhf = ahf->lConnMap2D[elementype - 2].num_verts_in_face;
03984     else if( ( elementype == MBTET ) || ( elementype == MBHEX ) )
03985     {
03986         int idx = ahf->get_index_in_lmap( *_incells.begin() );
03987         nhf     = ahf->lConnMap3D[idx].num_faces_in_cell;
03988     }
03989     else
03990         MB_SET_ERR( MB_FAILURE, "Requesting vertex boundary information for an unsupported entity type" );
03991 
03992     error = ahf->get_incident_map( elementype, vertex, ent, lid );MB_CHK_ERR( error );
03993     error = ahf->get_sibling_map( elementype, ent[0], &sibents[0], &siblids[0], nhf );MB_CHK_ERR( error );
03994 
03995     return ( sibents[lid[0]] == 0 );
03996 }
03997 
03998 bool NestedRefine::is_edge_on_boundary( const EntityHandle& entity )
03999 {
04000     ErrorCode error;
04001     bool is_border = false;
04002     if( meshdim == 1 )  // The edge has a vertex on the boundary in the curve mesh
04003     {
04004         EntityHandle sibents[2];
04005         int siblids[2];
04006         error = ahf->get_sibling_map( MBEDGE, entity, &sibents[0], &siblids[0], 2 );MB_CHK_ERR( error );
04007         for( int i = 0; i < 2; i++ )
04008         {
04009             if( sibents[i] == 0 )
04010             {
04011                 is_border = true;
04012                 break;
04013             }
04014         }
04015     }
04016     else if( meshdim == 2 )  // The edge is on the boundary of the 2d mesh
04017     {
04018         std::vector< EntityHandle > adjents;
04019         error = ahf->get_up_adjacencies_2d( entity, adjents );MB_CHK_ERR( error );
04020         if( adjents.size() == 1 ) is_border = true;
04021     }
04022     else if( meshdim == 3 )  // The edge is part of a face on the boundary of the 3d mesh
04023     {
04024         std::vector< EntityHandle > adjents;
04025         std::vector< int > leids;
04026         error = ahf->get_up_adjacencies_edg_3d( entity, adjents, &leids );MB_CHK_ERR( error );
04027         assert( !adjents.empty() );
04028 
04029         int index = ahf->get_index_in_lmap( adjents[0] );
04030         int nhf   = ahf->lConnMap3D[index].num_faces_in_cell;
04031 
04032         for( int i = 0; i < (int)adjents.size(); i++ )
04033         {
04034             EntityHandle sibents[6];
04035             int siblids[6];
04036             error = ahf->get_sibling_map( elementype, adjents[0], &sibents[0], &siblids[0], nhf );MB_CHK_ERR( error );
04037             for( int k = 0; k < 2; k++ )
04038             {
04039                 int hf = ahf->lConnMap3D[index].e2hf[leids[0]][k];
04040                 if( sibents[hf] == 0 )
04041                 {
04042                     is_border = true;
04043                     break;
04044                 }
04045             }
04046         }
04047     }
04048     return is_border;
04049 }
04050 
04051 bool NestedRefine::is_face_on_boundary( const EntityHandle& entity )
04052 {
04053     ErrorCode error;
04054     bool is_border = false;
04055 
04056     if( meshdim == 1 )
04057         MB_SET_ERR( MB_FAILURE, "Requesting boundary information for a face entity type on a curve mesh" );
04058     else if( meshdim == 2 )  // The face has a local edge on the boundary of the 2d mesh
04059     {
04060         EntityHandle sibents[4];
04061         int siblids[4];
04062         int nepf = ahf->lConnMap2D[elementype - 2].num_verts_in_face;
04063         error    = ahf->get_sibling_map( elementype, entity, &sibents[0], &siblids[0], nepf );MB_CHK_ERR( error );
04064 
04065         for( int i = 0; i < nepf; i++ )
04066         {
04067             if( sibents[i] == 0 )
04068             {
04069                 is_border = true;
04070                 break;
04071             }
04072         }
04073     }
04074     else if( meshdim == 3 )  // The face lies on the boundary of the 3d mesh
04075     {
04076         std::vector< EntityHandle > adjents;
04077         error = ahf->get_up_adjacencies_face_3d( entity, adjents );MB_CHK_ERR( error );
04078         if( adjents.size() == 1 ) is_border = true;
04079     }
04080     return is_border;
04081 }
04082 
04083 bool NestedRefine::is_cell_on_boundary( const EntityHandle& entity )
04084 {
04085     if( meshdim != 3 )
04086         MB_SET_ERR( MB_FAILURE, "Requesting boundary information for a cell entity type on a curve or surface mesh" );
04087 
04088     bool is_border = false;
04089     int index      = ahf->get_index_in_lmap( *_incells.begin() );
04090     int nfpc       = ahf->lConnMap3D[index].num_faces_in_cell;
04091     EntityHandle sibents[6];
04092     int siblids[6];
04093 
04094     ErrorCode error = ahf->get_sibling_map( elementype, entity, &sibents[0], &siblids[0], nfpc );MB_CHK_ERR( error );
04095 
04096     for( int i = 0; i < nfpc; i++ )
04097     {
04098         if( sibents[i] == 0 )
04099         {
04100             is_border = true;
04101             break;
04102         }
04103     }
04104     return is_border;
04105 }
04106 
04107 /* **********************************
04108  *          Helper Functions              *
04109  ************************************/
04110 ErrorCode NestedRefine::copy_vertices_from_prev_level( int cur_level )
04111 {
04112     ErrorCode error;
04113 
04114     if( cur_level )
04115     {
04116         int nverts_prev = level_mesh[cur_level - 1].num_verts;
04117         for( int i = 0; i < nverts_prev; i++ )
04118         {
04119             level_mesh[cur_level].coordinates[0][i] = level_mesh[cur_level - 1].coordinates[0][i];
04120             level_mesh[cur_level].coordinates[1][i] = level_mesh[cur_level - 1].coordinates[1][i];
04121             level_mesh[cur_level].coordinates[2][i] = level_mesh[cur_level - 1].coordinates[2][i];
04122         }
04123     }
04124     else  // Copy the vertices from the input mesh
04125     {
04126         int nverts_in = _inverts.size();
04127         std::vector< double > vcoords( 3 * nverts_in );
04128         error = mbImpl->get_coords( _inverts, &vcoords[0] );MB_CHK_ERR( error );
04129 
04130         for( int i = 0; i < nverts_in; i++ )
04131         {
04132             level_mesh[cur_level].coordinates[0][i] = vcoords[3 * i];
04133             level_mesh[cur_level].coordinates[1][i] = vcoords[3 * i + 1];
04134             level_mesh[cur_level].coordinates[2][i] = vcoords[3 * i + 2];
04135         }
04136     }
04137     return MB_SUCCESS;
04138     // To add: Map from old vertices to new duplicates: NOT NEEDED
04139 }
04140 
04141 ErrorCode NestedRefine::update_tracking_verts( EntityHandle cid, int cur_level, int deg,
04142                                                std::vector< EntityHandle >& trackvertsC_edg,
04143                                                std::vector< EntityHandle >& trackvertsC_face, EntityHandle* vbuffer )
04144 {
04145     // The vertices in the vbuffer are added to appropriate edges and faces of cells that are
04146     // incident on the working cell.
04147     ErrorCode error;
04148 
04149     EntityHandle cstart_prev;
04150     if( cur_level )
04151         cstart_prev = level_mesh[cur_level - 1].start_cell;
04152     else
04153         cstart_prev = *_incells.begin();
04154 
04155     EntityType cell_type = mbImpl->type_from_handle( cstart_prev );
04156     int cindex           = cell_type - 1;
04157     int d                = get_index_from_degree( deg );
04158 
04159     int nve = refTemplates[cindex][d].nv_edge;
04160     int nvf = refTemplates[cindex][d].nv_face;
04161 
04162     int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
04163     int nepc  = ahf->lConnMap3D[index].num_edges_in_cell;
04164     int nfpc  = ahf->lConnMap3D[index].num_faces_in_cell;
04165 
04166     // Step 1: Add the vertices on an edge of the working cell to tracking array of incident cells.
04167     for( int i = 0; i < nepc; i++ )
04168     {
04169         // Add the vertices to edges of the current cell
04170         for( int j = 0; j < nve; j++ )
04171         {
04172             int id  = refTemplates[cindex][d].vert_on_edges[i][j];
04173             int idx = cid - cstart_prev;
04174             int aid = idx * nve * nepc + nve * i + j;
04175 
04176             if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
04177         }
04178 
04179         // Obtain all the incident cells
04180         std::vector< EntityHandle > inc_cids;
04181         std::vector< int > inc_leids, inc_orient;
04182 
04183         error = ahf->get_up_adjacencies_edg_3d( cid, i, inc_cids, &inc_leids, &inc_orient );MB_CHK_ERR( error );
04184 
04185         if( inc_cids.size() == 1 ) continue;
04186 
04187         // Add the vertices to the edges of the incident cells
04188         for( int k = 0; k < (int)inc_cids.size(); k++ )
04189         {
04190             if( inc_cids[k] == cid ) continue;
04191 
04192             int idx = inc_cids[k] - cstart_prev;
04193 
04194             if( inc_orient[k] )  // Same edge direction as the current edge
04195             {
04196                 for( int j = 0; j < nve; j++ )
04197                 {
04198                     int id  = refTemplates[cindex][d].vert_on_edges[i][j];
04199                     int aid = idx * nve * nepc + nve * inc_leids[k] + j;
04200 
04201                     if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
04202                 }
04203             }
04204             else
04205             {
04206                 for( int j = 0; j < nve; j++ )
04207                 {
04208                     int id  = refTemplates[cindex][d].vert_on_edges[i][nve - j - 1];
04209                     int aid = idx * nve * nepc + nve * inc_leids[k] + j;
04210 
04211                     if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
04212                 }
04213             }
04214         }
04215     }
04216 
04217     // Step 2: Add the vertices on a face of the working cell to tracking array of incident cells.
04218     if( nvf )
04219     {
04220 
04221         for( int i = 0; i < nfpc; i++ )
04222         {
04223             // Add vertices to the tracking array of vertices on faces for the current cell
04224             std::vector< EntityHandle > face_vbuf( nvf, 0 );
04225             for( int j = 0; j < nvf; j++ )
04226             {
04227                 int id  = refTemplates[cindex][d].vert_on_faces[i][j];
04228                 int idx = cid - cstart_prev;
04229                 int aid = idx * nvf * nfpc + nvf * i + j;
04230 
04231                 if( !trackvertsC_face[aid] ) trackvertsC_face[aid] = vbuffer[id];
04232 
04233                 face_vbuf[j] = vbuffer[id];
04234             }
04235 
04236             // Obtain all the incident cells
04237             std::vector< EntityHandle > sib_cids;
04238             std::vector< int > sib_lfids;
04239             error = ahf->get_up_adjacencies_face_3d( cid, i, sib_cids, &sib_lfids );MB_CHK_ERR( error );
04240 
04241             if( sib_cids.size() == 1 ) continue;
04242 
04243             // Reorder the vertex local ids incident on the half-face
04244             std::vector< int > id_sib( nvf );
04245             for( int k = 0; k < nvf; k++ )
04246                 id_sib[k] = 0;
04247 
04248             error = reorder_indices( cur_level, deg, sib_cids[1], sib_lfids[1], cid, i, 0, &id_sib[0] );MB_CHK_ERR( error );
04249 
04250             // Add vertices to the tracking array of vertices on faces for the sibling cell of the
04251             // current cell
04252             for( int j = 0; j < nvf; j++ )
04253             {
04254                 int idx = sib_cids[1] - cstart_prev;
04255                 int aid = idx * nvf * nfpc + nvf * sib_lfids[1] + j;
04256 
04257                 if( !trackvertsC_face[aid] ) trackvertsC_face[aid] = face_vbuf[id_sib[j] - 1];
04258             }
04259         }
04260     }
04261     return MB_SUCCESS;
04262 }
04263 
04264 ErrorCode NestedRefine::reorder_indices( int cur_level, int deg, EntityHandle cell, int lfid, EntityHandle sib_cell,
04265                                          int sib_lfid, int index, int* id_sib )
04266 {
04267     // Reorders the indices of either vertices or children cell local ids to match with order of the
04268     // given cell and a local face. index = 0 : vertices,
04269     //           = 1 : face
04270 
04271     assert( deg == 2 || deg == 3 );
04272 
04273     ErrorCode error;
04274     int idx = ahf->get_index_in_lmap( *_incells.begin() );
04275     int nvF = ahf->lConnMap3D[idx].hf2v_num[lfid];
04276     int nco = permutation[nvF - 3].num_comb;
04277 
04278     if( !index && ( ( nvF == 3 && deg == 3 ) || ( nvF == 4 && deg == 2 ) ) ) { id_sib[0] = 1; }
04279     else
04280     {
04281         // Get connectivity of the cell and its sibling cell
04282         std::vector< EntityHandle > conn, sib_conn;
04283         error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
04284 
04285         error = get_connectivity( sib_cell, cur_level, sib_conn );MB_CHK_ERR( error );
04286 
04287         // Get the connectivity of the local face in the cell and its sibling
04288         std::vector< EntityHandle > lface( nvF );
04289         std::vector< EntityHandle > lface_sib( nvF );
04290         for( int i = 0; i < nvF; i++ )
04291         {
04292             int id   = ahf->lConnMap3D[idx].hf2v[lfid][i];
04293             lface[i] = conn[id];
04294 
04295             id           = ahf->lConnMap3D[idx].hf2v[sib_lfid][i];
04296             lface_sib[i] = sib_conn[id];
04297         }
04298 
04299         // Find the combination
04300         int c = 0;
04301         for( int i = 0; i < nco; i++ )
04302         {
04303             int count = 0;
04304             for( int j = 0; j < nvF; j++ )
04305             {
04306                 int id = permutation[nvF - 3].comb[i][j];
04307                 if( lface[j] == lface_sib[id] ) count += 1;
04308             }
04309 
04310             if( count == nvF )
04311             {
04312                 c = i;
04313                 break;
04314             }
04315         }
04316 
04317         if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
04318 
04319         // Get the ordered indices
04320         if( ( ( !index ) && ( nvF == 4 ) && ( deg == 3 ) ) || ( deg == 2 ) )
04321         {
04322             for( int i = 0; i < 4; i++ )
04323                 id_sib[i] = permutation[nvF - 3].porder2[c][i];
04324         }
04325         else
04326         {
04327             for( int i = 0; i < 9; i++ )
04328                 id_sib[i] = permutation[nvF - 3].porder3[c][i];
04329         }
04330     }
04331 
04332     return MB_SUCCESS;
04333 }
04334 
04335 ErrorCode NestedRefine::reorder_indices( int deg, EntityHandle* face1_conn, EntityHandle* face2_conn, int nvF,
04336                                          std::vector< int >& lemap, std::vector< int >& vidx, int* leorient )
04337 {
04338     // Given the connectivities of two faces, get the permuted indices w.r.t first face.
04339     // Step 1: First find the orientation
04340     int nco = permutation[nvF - 3].num_comb;
04341     int c   = 0;
04342     for( int i = 0; i < nco; i++ )
04343     {
04344         int count = 0;
04345         for( int j = 0; j < nvF; j++ )
04346         {
04347             int id = permutation[nvF - 3].comb[i][j];
04348             if( face1_conn[j] == face2_conn[id] ) count += 1;
04349         }
04350 
04351         if( count == nvF )
04352         {
04353             c = i;
04354             break;
04355         }
04356     }
04357 
04358     if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
04359 
04360     // Add the corresponding local edges
04361     lemap.reserve( nvF );
04362     for( int i = 0; i < nvF; i++ )
04363     {
04364         lemap.push_back( permutation[nvF - 3].lemap[c][i] );
04365     }
04366     if( leorient ) leorient[0] = permutation[nvF - 3].orient[c];
04367 
04368     if( nvF == 3 && deg == 2 ) return MB_SUCCESS;
04369 
04370     if( ( nvF == 3 && deg == 3 ) || ( nvF == 4 && deg == 2 ) ) { vidx.push_back( 1 ); }
04371     else if( nvF == 4 && deg == 3 )
04372     {
04373         for( int i = 0; i < 4; i++ )
04374             vidx.push_back( permutation[nvF - 3].porder2[c][i] );
04375     }
04376 
04377     return MB_SUCCESS;
04378 }
04379 
04380 ErrorCode NestedRefine::reorder_indices( int deg, int nvF, int comb, int* childfid_map )
04381 {
04382     // Given connectivities of two faces and a degree, get the permuted indices of the children
04383     // faces w.r.t first face.
04384 
04385     assert( deg == 2 || deg == 3 );
04386 
04387     // Get the ordered indices
04388     if( deg == 2 )
04389     {
04390         for( int i = 0; i < 4; i++ )
04391             childfid_map[i] = permutation[nvF - 3].porder2[comb][i];
04392     }
04393     else
04394     {
04395         for( int i = 0; i < 9; i++ )
04396             childfid_map[i] = permutation[nvF - 3].porder3[comb][i];
04397     }
04398 
04399     return MB_SUCCESS;
04400 }
04401 
04402 ErrorCode NestedRefine::reorder_indices( EntityHandle* face1_conn, EntityHandle* face2_conn, int nvF, int* conn_map,
04403                                          int& comb, int* orient )
04404 {
04405     // Given connectivities of two faces and a degree, get the permuted indices of the children
04406     // faces w.r.t first face.
04407 
04408     // Step 1: First find the combination
04409     int nco = permutation[nvF - 3].num_comb;
04410     int c   = 0;
04411     for( int i = 0; i < nco; i++ )
04412     {
04413         int count = 0;
04414         for( int j = 0; j < nvF; j++ )
04415         {
04416             int id = permutation[nvF - 3].comb[i][j];
04417             if( face1_conn[j] == face2_conn[id] ) count += 1;
04418         }
04419 
04420         if( count == nvF )
04421         {
04422             c = i;
04423             break;
04424         }
04425     }
04426 
04427     if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
04428 
04429     comb = c;
04430 
04431     if( orient ) orient[0] = permutation[nvF - 3].orient[c];
04432 
04433     for( int j = 0; j < nvF; j++ )
04434     {
04435         conn_map[j] = permutation[nvF - 3].comb[c][j];
04436     }
04437 
04438     return MB_SUCCESS;
04439 }
04440 
04441 ErrorCode NestedRefine::count_subentities( EntityHandle set, int cur_level, int* nedges, int* nfaces )
04442 {
04443     ErrorCode error;
04444 
04445     if( cur_level >= 0 )
04446     {
04447         Range edges, faces, cells;
04448 
04449         error = mbImpl->get_entities_by_dimension( set, 1, edges );MB_CHK_ERR( error );
04450 
04451         error = mbImpl->get_entities_by_dimension( set, 2, faces );MB_CHK_ERR( error );
04452 
04453         error = mbImpl->get_entities_by_dimension( set, 3, cells );MB_CHK_ERR( error );
04454 
04455         error = ahf->count_subentities( edges, faces, cells, nedges, nfaces );MB_CHK_ERR( error );
04456     }
04457     else
04458     {
04459         error = ahf->count_subentities( _inedges, _infaces, _incells, nedges, nfaces );MB_CHK_ERR( error );
04460     }
04461 
04462     return MB_SUCCESS;
04463 }
04464 
04465 ErrorCode NestedRefine::get_octahedron_corner_coords( int cur_level, int deg, EntityHandle* vbuffer, double* ocoords )
04466 {
04467     int lid[6] = { 0, 0, 0, 0, 0, 0 };
04468 
04469     if( deg == 2 )
04470     {
04471         lid[0] = 5;
04472         lid[1] = 8;
04473         lid[2] = 9;
04474         lid[3] = 6;
04475         lid[4] = 4;
04476         lid[5] = 7;
04477     }
04478     else if( deg == 3 )
04479     {
04480         lid[0] = 19;
04481         lid[1] = 16;
04482         lid[2] = 18;
04483         lid[3] = 9;
04484         lid[4] = 4;
04485         lid[5] = 10;
04486     }
04487 
04488     EntityHandle vstart = level_mesh[cur_level].start_vertex;
04489 
04490     for( int i = 0; i < 6; i++ )
04491     {
04492         EntityHandle vid   = vbuffer[lid[i]];
04493         ocoords[3 * i]     = level_mesh[cur_level].coordinates[0][vid - vstart];
04494         ocoords[3 * i + 1] = level_mesh[cur_level].coordinates[1][vid - vstart];
04495         ocoords[3 * i + 2] = level_mesh[cur_level].coordinates[2][vid - vstart];
04496     }
04497 
04498     return MB_SUCCESS;
04499 }
04500 
04501 int NestedRefine::find_shortest_diagonal_octahedron( int cur_level, int deg, EntityHandle* vbuffer )
04502 {
04503     ErrorCode error;
04504     double coords[18];
04505     error = get_octahedron_corner_coords( cur_level, deg, vbuffer, coords );
04506     if( error != MB_SUCCESS ) MB_SET_ERR( MB_FAILURE, "Error in obtaining octahedron corner coordinates" );
04507 
04508     int diag_map[6] = { 1, 3, 2, 4, 5, 0 };
04509     double length   = std::numeric_limits< double >::max();
04510 
04511     int diag = 0;
04512     double x, y, z;
04513     x = y = z = 0;
04514 
04515     for( int d = 0; d < 3; d++ )
04516     {
04517         int id1     = diag_map[2 * d];
04518         int id2     = diag_map[2 * d + 1];
04519         x           = coords[3 * id1] - coords[3 * id2];
04520         y           = coords[3 * id1 + 1] - coords[3 * id2 + 1];
04521         z           = coords[3 * id1 + 2] - coords[3 * id2 + 2];
04522         double dlen = sqrt( x * x + y * y + z * z );
04523         if( dlen < length )
04524         {
04525             length = dlen;
04526             diag   = d + 1;
04527         }
04528     }
04529 
04530     return diag;
04531 }
04532 
04533 int NestedRefine::get_local_vid( EntityHandle vid, EntityHandle ent, int level )
04534 {
04535     ErrorCode error;
04536     // Given a vertex, find its local id in the given entity
04537     std::vector< EntityHandle > conn;
04538 
04539     error = get_connectivity( ent, level + 1, conn );
04540     if( error != MB_SUCCESS ) MB_SET_ERR( MB_FAILURE, "Error in getting connectivity of the requested entity" );
04541 
04542     int lid = -1;
04543     for( int i = 0; i < (int)conn.size(); i++ )
04544     {
04545         if( conn[i] == vid )
04546         {
04547             lid = i;
04548             break;
04549         }
04550     }
04551     if( lid < 0 ) MB_SET_ERR( MB_FAILURE, "Error in getting local vertex id in the given entity" );
04552     return lid;
04553 }
04554 
04555 int NestedRefine::get_index_from_degree( int degree )
04556 {
04557     int d = deg_index.find( degree )->second;
04558     return d;
04559 }
04560 
04561 /*
04562 ErrorCode NestedRefine::print_maps_1D(int level)
04563 {
04564   ErrorCode error;
04565   int nv, ne;
04566   nv = level_mesh[level].num_verts;
04567   ne = level_mesh[level].num_edges;
04568 
04569   EntityHandle start_edge = level_mesh[level].start_edge;
04570 
04571   //V2HV
04572     std::cout<<"<V2HV_EID, V2HV_LVID>"<<std::endl;
04573   for (int i=0; i<nv; i++)
04574     {
04575       EntityHandle eid=0;
04576       int lvid=0;
04577       EntityHandle vid = level_mesh[level].start_vertex+i;
04578       error = ahf->get_incident_map(MBEDGE, vid, eid, lvid); MB_CHK_ERR(error);
04579 
04580       std::cout<<"For vertex = "<<vid<<"::Incident halfvertex "<<eid<<"  "<<lvid<<std::endl;
04581     }
04582 
04583   //SIBHVS
04584   std::cout<<"start_edge = "<<start_edge<<std::endl;
04585   std::cout<<"<SIBHVS_EID,SIBHVS_LVID>"<<std::endl;
04586   for (int i=0; i<ne; i++)
04587     {
04588       EntityHandle ent = start_edge+i;
04589 
04590       EntityHandle eid[2];  int lvid[2];
04591       error = ahf->get_sibling_map(MBEDGE, ent, &eid[0], &lvid[0], 2); MB_CHK_ERR(error);
04592       std::cout<<"<"<<eid[0]<<","<<lvid[0]<<">"<<" "<<"<"<<eid[1]<<","<<lvid[1]<<">"<<std::endl;
04593     }
04594 
04595   return MB_SUCCESS;
04596 }
04597 
04598 ErrorCode NestedRefine::print_maps_2D(int level, EntityType type)
04599 {
04600   ErrorCode error;
04601   int nv, nf;
04602   nv = level_mesh[level].num_verts;
04603   nf = level_mesh[level].num_faces;
04604 
04605   EntityHandle start_face = level_mesh[level].start_face;
04606 
04607   //V2HV
04608     std::cout<<"<V2HE_FID, V2HE_LEID>"<<std::endl;
04609   for (int i=0; i<nv; i++)
04610     {
04611       EntityHandle fid=0;
04612       int leid=0;
04613       EntityHandle vid = level_mesh[level].start_vertex+i;
04614       error = ahf->get_incident_map(type, vid, fid, leid); MB_CHK_ERR(error);
04615 
04616       std::cout<<"For vertex = "<<vid<<"::Incident halfedge "<<fid<<"  "<<leid<<std::endl;
04617     }
04618 
04619   //SIBHES
04620   std::cout<<"start_face = "<<start_face<<std::endl;
04621   std::cout<<"<SIBHES_FID,SIBHES_LEID>"<<std::endl;
04622   EntityType ftype = mbImpl->type_from_handle(*_infaces.begin());
04623   int nepf = ahf->lConnMap2D[ftype-2].num_verts_in_face;
04624 
04625   EntityHandle *fid = new EntityHandle[nepf];
04626   int *leid = new int[nepf];
04627 
04628   for (int i=0; i<nf; i++)
04629     {
04630       for (int j=0; j<nepf; j++)
04631         {
04632           fid[j] = 0;
04633           leid[j] = 0;
04634         }
04635 
04636       EntityHandle ent = start_face+i;
04637       error = ahf->get_sibling_map(type, ent, fid, leid, nepf); MB_CHK_ERR(error);
04638 
04639       for (int j=0; j<nepf; j++){
04640           std::cout<<"<"<<fid[j]<<","<<leid[j]<<">"<<"      ";
04641       }
04642       std::cout<<std::endl;
04643     }
04644 
04645   delete [] fid;
04646   delete [] leid;
04647 
04648   return MB_SUCCESS;
04649 }
04650 
04651 ErrorCode NestedRefine::print_maps_3D(int level, EntityType type)
04652 {
04653   ErrorCode error;
04654   int nv, nc;
04655   nv = level_mesh[level].num_verts;
04656   nc = level_mesh[level].num_cells;
04657   EntityHandle start_cell = level_mesh[level].start_cell;
04658 
04659   //V2HF
04660     std::cout<<"<V2HF_CID, V2HF_LFID>"<<std::endl;
04661   for (int i=0; i<nv; i++)
04662     {
04663       EntityHandle cid=0;
04664       int lfid=0;
04665       EntityHandle vid = level_mesh[level].start_vertex+i;
04666       error = ahf->get_incident_map(type, vid, cid, lfid); MB_CHK_ERR(error);
04667 
04668       std::cout<<"For vertex = "<<vid<<"::Incident halfface "<<cid<<"  "<<lfid<<std::endl;
04669     }
04670 
04671   //SIBHFS
04672   std::cout<<"start_cell = "<<start_cell<<std::endl;
04673   std::cout<<"<SIBHFS_CID,SIBHFS_LFID>"<<std::endl;
04674   int index = ahf->get_index_in_lmap(start_cell);
04675   int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
04676 
04677   EntityHandle *cid = new EntityHandle[nfpc];
04678   int *lfid = new int[nfpc];
04679   for (int i=0; i<nc; i++)
04680     {
04681       for (int k=0; k<nfpc; k++)
04682         {
04683           cid[k] = 0;
04684           lfid[k] = 0;
04685         }
04686 
04687       EntityHandle ent = start_cell+i;
04688       error = ahf->get_sibling_map(type, ent, cid, lfid, nfpc); MB_CHK_ERR(error);
04689 
04690       for (int j=0; j<nfpc; j++){
04691           std::cout<<"<"<<cid[j]<<","<<lfid[j]<<">"<<"      ";
04692       }
04693       std::cout<<std::endl;
04694     }
04695 
04696   delete [] cid;
04697   delete [] lfid;
04698 
04699   return MB_SUCCESS;
04700 }
04701 
04702 */
04703 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines