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