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