MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /*This unit test is for the uniform refinement capability based on AHF datastructures*/ 00002 #include <iostream> 00003 #include <string> 00004 #include <sstream> 00005 #include <ctime> 00006 #include <vector> 00007 #include <algorithm> 00008 #include "moab/Core.hpp" 00009 #include "moab/Range.hpp" 00010 #include "moab/MeshTopoUtil.hpp" 00011 #include "moab/ReadUtilIface.hpp" 00012 #include "moab/NestedRefine.hpp" 00013 #include "TestUtil.hpp" 00014 00015 #ifdef MOAB_HAVE_MPI 00016 #include "moab/ParallelComm.hpp" 00017 #include "MBParallelConventions.h" 00018 #include "ReadParallel.hpp" 00019 #include "moab/FileOptions.hpp" 00020 #include "MBTagConventions.hpp" 00021 #include "moab_mpi.h" 00022 #endif 00023 00024 using namespace moab; 00025 00026 #ifdef MOAB_HAVE_MPI 00027 std::string read_options; 00028 #endif 00029 00030 int number_tests_successful = 0; 00031 int number_tests_failed = 0; 00032 00033 void handle_error_code( ErrorCode rv, int& number_failed, int& number_successful ) 00034 { 00035 if( rv == MB_SUCCESS ) 00036 { 00037 #ifdef MOAB_HAVE_MPI 00038 int rank = 0; 00039 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 00040 if( rank == 0 ) std::cout << "Success"; 00041 #else 00042 std::cout << "Success"; 00043 #endif 00044 number_successful++; 00045 } 00046 else 00047 { 00048 std::cout << "Failure"; 00049 number_failed++; 00050 } 00051 } 00052 00053 ErrorCode test_adjacencies( Interface* mbImpl, NestedRefine* nr, Range all_ents ) 00054 { 00055 MeshTopoUtil mtu( mbImpl ); 00056 ErrorCode error; 00057 Range verts, edges, faces, cells; 00058 verts = all_ents.subset_by_dimension( 0 ); 00059 edges = all_ents.subset_by_dimension( 1 ); 00060 faces = all_ents.subset_by_dimension( 2 ); 00061 cells = all_ents.subset_by_dimension( 3 ); 00062 00063 std::vector< EntityHandle > adjents; 00064 Range mbents, ahfents; 00065 00066 // Update the moab native data structures 00067 ReadUtilIface* read_face; 00068 error = mbImpl->query_interface( read_face );CHECK_ERR( error ); 00069 std::vector< EntityHandle > ents, conn; 00070 00071 if( !edges.empty() ) 00072 { 00073 conn.clear(); 00074 ents.clear(); 00075 00076 for( Range::iterator it = edges.begin(); it != edges.end(); it++ ) 00077 ents.push_back( *it ); 00078 // std::copy(edges.begin(), edges.end(), ents.begin()); 00079 error = mbImpl->get_connectivity( &ents[0], (int)ents.size(), conn );CHECK_ERR( error ); 00080 error = read_face->update_adjacencies( *ents.begin(), (int)ents.size(), 2, &conn[0] );CHECK_ERR( error ); 00081 } 00082 00083 if( !faces.empty() ) 00084 { 00085 conn.clear(); 00086 ents.clear(); 00087 00088 for( Range::iterator it = faces.begin(); it != faces.end(); it++ ) 00089 ents.push_back( *it ); 00090 00091 // std::copy(faces.begin(), faces.end(), ents.begin()); 00092 error = mbImpl->get_connectivity( &ents[0], 1, conn );CHECK_ERR( error ); 00093 int nvF = conn.size(); 00094 conn.clear(); 00095 error = mbImpl->get_connectivity( &ents[0], (int)ents.size(), conn );CHECK_ERR( error ); 00096 error = read_face->update_adjacencies( *ents.begin(), (int)ents.size(), nvF, &conn[0] );CHECK_ERR( error ); 00097 } 00098 00099 if( !cells.empty() ) 00100 { 00101 conn.clear(); 00102 ents.clear(); 00103 00104 for( Range::iterator it = cells.begin(); it != cells.end(); it++ ) 00105 ents.push_back( *it ); 00106 // std::copy(cells.begin(), cells.end(), ents.begin()); 00107 error = mbImpl->get_connectivity( &ents[0], 1, conn );CHECK_ERR( error ); 00108 int nvF = conn.size(); 00109 conn.clear(); 00110 error = mbImpl->get_connectivity( &ents[0], (int)ents.size(), conn );CHECK_ERR( error ); 00111 error = read_face->update_adjacencies( *ents.begin(), (int)ents.size(), nvF, &conn[0] );CHECK_ERR( error ); 00112 } 00113 00114 if( !edges.empty() ) 00115 { 00116 // 1D Queries // 00117 // IQ1: For every vertex, obtain incident edges 00118 for( Range::iterator i = verts.begin(); i != verts.end(); ++i ) 00119 { 00120 adjents.clear(); 00121 mbents.clear(); 00122 ahfents.clear(); 00123 error = nr->get_adjacencies( *i, 1, adjents );CHECK_ERR( error ); 00124 error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error ); 00125 std::sort( adjents.begin(), adjents.end() ); 00126 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00127 00128 if( ahfents.size() != mbents.size() ) 00129 { 00130 std::cout << "ahf results = " << std::endl; 00131 ahfents.print(); 00132 std::cout << "native results = " << std::endl; 00133 mbents.print(); 00134 } 00135 00136 CHECK_EQUAL( adjents.size(), mbents.size() ); 00137 mbents = subtract( mbents, ahfents ); 00138 if( ahfents.size() != mbents.size() ) 00139 { 00140 // std::cout<<"ahf results = "<<std::endl; 00141 // ahfents.print(); 00142 // std::cout<<"native results = "<<std::endl; 00143 // mbents.print(); 00144 } 00145 // CHECK_EQUAL(adjents.size(),mbents.size()); 00146 CHECK( !mbents.size() ); 00147 } 00148 00149 // NQ1: For every edge, obtain neighbor edges 00150 for( Range::iterator i = edges.begin(); i != edges.end(); ++i ) 00151 { 00152 adjents.clear(); 00153 mbents.clear(); 00154 ahfents.clear(); 00155 error = nr->get_adjacencies( *i, 1, adjents );CHECK_ERR( error ); 00156 error = mtu.get_bridge_adjacencies( *i, 0, 1, mbents );CHECK_ERR( error ); 00157 CHECK_EQUAL( adjents.size(), mbents.size() ); 00158 std::sort( adjents.begin(), adjents.end() ); 00159 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00160 mbents = subtract( mbents, ahfents ); 00161 CHECK( !mbents.size() ); 00162 } 00163 } 00164 00165 if( !faces.empty() ) 00166 { 00167 // IQ21: For every vertex, obtain incident faces 00168 for( Range::iterator i = verts.begin(); i != verts.end(); ++i ) 00169 { 00170 adjents.clear(); 00171 mbents.clear(); 00172 ahfents.clear(); 00173 error = nr->get_adjacencies( *i, 2, adjents );CHECK_ERR( error ); 00174 error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error ); 00175 CHECK_EQUAL( adjents.size(), mbents.size() ); 00176 std::sort( adjents.begin(), adjents.end() ); 00177 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00178 mbents = subtract( mbents, ahfents ); 00179 CHECK( !mbents.size() ); 00180 } 00181 00182 // IQ22: For every edge, obtain incident faces 00183 if( !edges.empty() ) 00184 { 00185 for( Range::iterator i = edges.begin(); i != edges.end(); ++i ) 00186 { 00187 adjents.clear(); 00188 mbents.clear(); 00189 ahfents.clear(); 00190 error = nr->get_adjacencies( *i, 2, adjents );CHECK_ERR( error ); 00191 error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error ); 00192 CHECK_EQUAL( adjents.size(), mbents.size() ); 00193 std::sort( adjents.begin(), adjents.end() ); 00194 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00195 mbents = subtract( mbents, ahfents ); 00196 CHECK( !mbents.size() ); 00197 } 00198 } 00199 00200 // NQ2: For every face, obtain neighbor faces 00201 for( Range::iterator i = faces.begin(); i != faces.end(); ++i ) 00202 { 00203 adjents.clear(); 00204 mbents.clear(); 00205 ahfents.clear(); 00206 error = nr->get_adjacencies( *i, 2, adjents );CHECK_ERR( error ); 00207 error = mtu.get_bridge_adjacencies( *i, 1, 2, mbents );CHECK_ERR( error ); 00208 CHECK_EQUAL( adjents.size(), mbents.size() ); 00209 std::sort( adjents.begin(), adjents.end() ); 00210 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00211 mbents = subtract( mbents, ahfents ); 00212 CHECK( !mbents.size() ); 00213 } 00214 00215 if( !edges.empty() ) 00216 { 00217 for( Range::iterator i = faces.begin(); i != faces.end(); ++i ) 00218 { 00219 adjents.clear(); 00220 mbents.clear(); 00221 ahfents.clear(); 00222 error = nr->get_adjacencies( *i, 1, adjents );CHECK_ERR( error ); 00223 error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error ); 00224 CHECK_EQUAL( adjents.size(), mbents.size() ); 00225 std::sort( adjents.begin(), adjents.end() ); 00226 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00227 mbents = subtract( mbents, ahfents ); 00228 CHECK( !mbents.size() ); 00229 } 00230 } 00231 } 00232 00233 if( !cells.empty() ) 00234 { 00235 // IQ 31: For every vertex, obtain incident cells 00236 for( Range::iterator i = verts.begin(); i != verts.end(); ++i ) 00237 { 00238 adjents.clear(); 00239 mbents.clear(); 00240 ahfents.clear(); 00241 error = nr->get_adjacencies( *i, 3, adjents );CHECK_ERR( error ); 00242 error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error ); 00243 00244 if( adjents.size() != mbents.size() ) 00245 { 00246 // std::cout<<"ahf results = "<<std::endl; 00247 // ahfents.print(); 00248 // std::cout<<"native results = "<<std::endl; 00249 // mbents.print(); 00250 } 00251 00252 CHECK_EQUAL( adjents.size(), mbents.size() ); 00253 std::sort( adjents.begin(), adjents.end() ); 00254 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00255 mbents = subtract( mbents, ahfents ); 00256 CHECK( !mbents.size() ); 00257 } 00258 00259 if( !edges.empty() ) 00260 { 00261 for( Range::iterator i = edges.begin(); i != edges.end(); ++i ) 00262 { 00263 adjents.clear(); 00264 mbents.clear(); 00265 ahfents.clear(); 00266 error = nr->get_adjacencies( *i, 3, adjents );CHECK_ERR( error ); 00267 error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error ); 00268 CHECK_EQUAL( adjents.size(), mbents.size() ); 00269 std::sort( adjents.begin(), adjents.end() ); 00270 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00271 mbents = subtract( mbents, ahfents ); 00272 CHECK( !mbents.size() ); 00273 } 00274 } 00275 00276 if( !faces.empty() ) 00277 { 00278 for( Range::iterator i = faces.begin(); i != faces.end(); ++i ) 00279 { 00280 adjents.clear(); 00281 mbents.clear(); 00282 ahfents.clear(); 00283 error = nr->get_adjacencies( *i, 3, adjents );CHECK_ERR( error ); 00284 error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error ); 00285 CHECK_EQUAL( adjents.size(), mbents.size() ); 00286 std::sort( adjents.begin(), adjents.end() ); 00287 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00288 mbents = subtract( mbents, ahfents ); 00289 CHECK( !mbents.size() ); 00290 } 00291 } 00292 00293 // NQ3: For every cell, obtain neighbor cells 00294 for( Range::iterator i = cells.begin(); i != cells.end(); ++i ) 00295 { 00296 adjents.clear(); 00297 mbents.clear(); 00298 ahfents.clear(); 00299 error = nr->get_adjacencies( *i, 3, adjents );CHECK_ERR( error ); 00300 error = mtu.get_bridge_adjacencies( *i, 2, 3, mbents );CHECK_ERR( error ); 00301 CHECK_EQUAL( adjents.size(), mbents.size() ); 00302 std::sort( adjents.begin(), adjents.end() ); 00303 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00304 mbents = subtract( mbents, ahfents ); 00305 CHECK( !mbents.size() ); 00306 } 00307 00308 if( !edges.empty() ) 00309 { 00310 for( Range::iterator i = cells.begin(); i != cells.end(); ++i ) 00311 { 00312 adjents.clear(); 00313 mbents.clear(); 00314 ahfents.clear(); 00315 error = nr->get_adjacencies( *i, 1, adjents );CHECK_ERR( error ); 00316 error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error ); 00317 CHECK_EQUAL( adjents.size(), mbents.size() ); 00318 std::sort( adjents.begin(), adjents.end() ); 00319 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00320 mbents = subtract( mbents, ahfents ); 00321 CHECK( !mbents.size() ); 00322 } 00323 } 00324 00325 if( !faces.empty() ) 00326 { 00327 for( Range::iterator i = cells.begin(); i != cells.end(); ++i ) 00328 { 00329 adjents.clear(); 00330 mbents.clear(); 00331 ahfents.clear(); 00332 error = nr->get_adjacencies( *i, 2, adjents );CHECK_ERR( error ); 00333 error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error ); 00334 CHECK_EQUAL( adjents.size(), mbents.size() ); 00335 std::sort( adjents.begin(), adjents.end() ); 00336 std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) ); 00337 mbents = subtract( mbents, ahfents ); 00338 CHECK( !mbents.size() ); 00339 } 00340 } 00341 } 00342 00343 return MB_SUCCESS; 00344 } 00345 00346 ErrorCode refine_entities( Interface* mb, ParallelComm* pc, EntityHandle fset, int* level_degrees, const int num_levels, 00347 bool output ) 00348 { 00349 ErrorCode error; 00350 00351 // Get the range of entities in the initial mesh 00352 Range init_ents[4]; 00353 00354 int dim[3] = { 1, 2, 3 }; 00355 00356 if( output ) 00357 { 00358 // int inents = init_ents.size(); 00359 std::stringstream file; 00360 #ifdef MOAB_HAVE_MPI 00361 file << "MESH_LEVEL_0.h5m"; 00362 #else 00363 file << "MESH_LEVEL_0.vtk"; 00364 #endif 00365 std::string str = file.str(); 00366 const char* output_file = str.c_str(); 00367 #ifdef MOAB_HAVE_MPI 00368 error = mb->write_file( output_file, 0, ";;PARALLEL=WRITE_PART" );CHECK_ERR( error ); 00369 #else 00370 error = mb->write_file( output_file, 0, NULL );CHECK_ERR( error ); 00371 #endif 00372 } 00373 00374 // Create an hm object and generate the hierarchy 00375 std::cout << "Creating a hm object" << std::endl; 00376 00377 #ifdef MOAB_HAVE_MPI 00378 Range averts, aedges, afaces, acells; 00379 error = mb->get_entities_by_dimension( fset, 0, averts );MB_CHK_ERR( error ); 00380 error = mb->get_entities_by_dimension( fset, 1, aedges );MB_CHK_ERR( error ); 00381 error = mb->get_entities_by_dimension( fset, 2, afaces );MB_CHK_ERR( error ); 00382 error = mb->get_entities_by_dimension( fset, 3, acells );MB_CHK_ERR( error ); 00383 00384 /* filter based on parallel status */ 00385 if( pc ) 00386 { 00387 error = pc->filter_pstatus( averts, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[0] );MB_CHK_ERR( error ); 00388 error = pc->filter_pstatus( aedges, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[1] );MB_CHK_ERR( error ); 00389 error = pc->filter_pstatus( afaces, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[2] );MB_CHK_ERR( error ); 00390 error = pc->filter_pstatus( acells, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[3] );MB_CHK_ERR( error ); 00391 } 00392 else 00393 { 00394 init_ents[0] = averts; 00395 init_ents[1] = aedges; 00396 init_ents[2] = afaces; 00397 init_ents[3] = acells; 00398 } 00399 #else 00400 error = mb->get_entities_by_dimension( fset, 0, init_ents[0] );CHECK_ERR( error ); 00401 error = mb->get_entities_by_dimension( fset, 1, init_ents[1] );CHECK_ERR( error ); 00402 error = mb->get_entities_by_dimension( fset, 2, init_ents[2] );CHECK_ERR( error ); 00403 error = mb->get_entities_by_dimension( fset, 3, init_ents[3] );CHECK_ERR( error ); 00404 #endif 00405 00406 NestedRefine uref( dynamic_cast< Core* >( mb ), pc, fset ); 00407 std::vector< EntityHandle > set; 00408 00409 std::cout << "Starting hierarchy generation" << std::endl; 00410 bool opt = true; 00411 // bool opt = false; 00412 error = uref.generate_mesh_hierarchy( num_levels, level_degrees, set, opt );CHECK_ERR( error ); 00413 std::cout << "Finished hierarchy generation in " << uref.timeall.tm_total << " secs" << std::endl; 00414 #ifdef MOAB_HAVE_MPI 00415 if( pc ) 00416 { 00417 std::cout << "Time taken for refinement " << uref.timeall.tm_refine << " secs" << std::endl; 00418 std::cout << "Time taken for resolving shared interface " << uref.timeall.tm_resolve << " secs" << std::endl; 00419 } 00420 #endif 00421 00422 // error = uref.exchange_ghosts(set, 1); CHECK_ERR(error); 00423 00424 std::cout << std::endl; 00425 std::cout << "Mesh size for level 0 :: inverts = " << init_ents[0].size() << ", inedges = " << init_ents[1].size() 00426 << ", infaces = " << init_ents[2].size() << ", incells = " << init_ents[3].size() << std::endl; 00427 00428 Range prev_ents[4]; 00429 for( int i = 0; i < 4; i++ ) 00430 prev_ents[i] = init_ents[i]; 00431 00432 // Loop over each mesh level and check its topological properties 00433 for( int l = 0; l < num_levels; l++ ) 00434 { 00435 Range all_ents; 00436 error = mb->get_entities_by_handle( set[l + 1], all_ents );CHECK_ERR( error ); 00437 00438 Range ents[4]; 00439 for( int k = 0; k < 4; k++ ) 00440 ents[k] = all_ents.subset_by_dimension( k ); 00441 00442 if( ents[0].empty() || all_ents.empty() ) std::cout << "Something is not right" << std::endl; 00443 00444 std::cout << std::endl; 00445 std::cout << "Mesh size for level " << l + 1 << " :: nverts = " << ents[0].size() 00446 << ", nedges = " << ents[1].size() << ", nfaces = " << ents[2].size() 00447 << ", ncells = " << ents[3].size() << std::endl; 00448 00449 // Check if the number of new entities created are correct. 00450 00451 for( int type = 1; type < 3; type++ ) 00452 { 00453 int factor = 1; 00454 if( !ents[type + 1].empty() ) 00455 { 00456 00457 for( int p = 0; p <= l; p++ ) 00458 { 00459 for( int d = 0; d < dim[type]; d++ ) 00460 factor *= level_degrees[p]; 00461 } 00462 int expected_nents = factor * init_ents[type + 1].size(); 00463 CHECK_EQUAL( expected_nents, (int)ents[type + 1].size() ); 00464 } 00465 } 00466 00467 // Check adjacencies 00468 error = test_adjacencies( mb, &uref, all_ents );CHECK_ERR( error ); 00469 00470 // Check interlevel child-parent query between previous and current level 00471 for( int type = 1; type < 3; type++ ) 00472 { 00473 if( !prev_ents[type + 1].empty() ) 00474 { 00475 for( Range::iterator e = prev_ents[type + 1].begin(); e != prev_ents[type + 1].end(); e++ ) 00476 { 00477 std::vector< EntityHandle > children; 00478 error = uref.parent_to_child( *e, l, l + 1, children );CHECK_ERR( error ); 00479 for( int i = 0; i < (int)children.size(); i++ ) 00480 { 00481 EntityHandle parent; 00482 error = uref.child_to_parent( children[i], l + 1, l, &parent );CHECK_ERR( error ); 00483 assert( parent == *e ); 00484 } 00485 } 00486 } 00487 } 00488 00489 for( int i = 0; i < 4; i++ ) 00490 prev_ents[i] = ents[i]; 00491 00492 // Print out the boundary vertices 00493 /* EntityHandle bnd_set; 00494 error = mb->create_meshset(MESHSET_SET, bnd_set); MB_CHK_ERR(error); 00495 std::vector<EntityHandle> vbnd; 00496 00497 for (int k=0; k<3; k++){ 00498 std::cout<<"Finding boundary of dimension "<<k<<" with size 00499 "<<ents[k].size()<<std::endl; if (ents[k].size() != 0){ for (Range::iterator v = 00500 ents[k].begin(); v != ents[k].end(); v++) 00501 { 00502 EntityHandle ent = *v; 00503 bool bnd = uref.is_entity_on_boundary(ent); 00504 if (bnd) 00505 vbnd.push_back(*v); 00506 } 00507 } 00508 } 00509 00510 std::cout<<"vbnd.size = "<<vbnd.size()<<std::endl; 00511 error = mb->add_entities(bnd_set, &vbnd[0], (int)vbnd.size()); MB_CHK_ERR(error); 00512 if (output) 00513 { 00514 std::stringstream file; 00515 file << "VBND_LEVEL_" <<l+1<<".vtk"; 00516 std::string str = file.str(); 00517 const char* output_file = str.c_str(); 00518 char * write_opts = NULL; 00519 error = mb->write_file(output_file, 0, write_opts, &bnd_set, 1); CHECK_ERR(error); 00520 }*/ 00521 00522 // Print out the mesh 00523 if( output ) 00524 { 00525 std::stringstream file; 00526 00527 #ifdef MOAB_HAVE_MPI 00528 file << "MESH_LEVEL_" << l + 1 << ".h5m"; 00529 #else 00530 file << "MESH_LEVEL_" << l + 1 << ".vtk"; 00531 #endif 00532 std::string str = file.str(); 00533 const char* output_file = str.c_str(); 00534 #ifdef MOAB_HAVE_MPI 00535 error = mb->write_file( output_file, 0, ";;PARALLEL=WRITE_PART", &set[l + 1], 1 );CHECK_ERR( error ); 00536 #else 00537 error = mb->write_file( output_file, 0, NULL, &set[l + 1], 1 );CHECK_ERR( error ); 00538 #endif 00539 } 00540 } 00541 00542 // Check interlevel child-parent query between initial and most refined mesh 00543 for( int type = 1; type < 3; type++ ) 00544 { 00545 if( !init_ents[type + 1].empty() ) 00546 { 00547 for( Range::iterator e = init_ents[type + 1].begin(); e != init_ents[type + 1].end(); e++ ) 00548 { 00549 std::vector< EntityHandle > children; 00550 error = uref.parent_to_child( *e, 0, num_levels, children );CHECK_ERR( error ); 00551 for( int i = 0; i < (int)children.size(); i++ ) 00552 { 00553 EntityHandle parent; 00554 error = uref.child_to_parent( children[i], num_levels, 0, &parent );CHECK_ERR( error ); 00555 assert( parent == *e ); 00556 } 00557 } 00558 } 00559 } 00560 00561 // Print out the whole hierarchy into a single file 00562 /* if (output) 00563 { 00564 std::stringstream file; 00565 file << "MESH_HIERARCHY.vtk"; 00566 std::string str = file.str(); 00567 const char* output_file = str.c_str(); 00568 error = mb->write_file(output_file); CHECK_ERR(error); 00569 }*/ 00570 00571 return MB_SUCCESS; 00572 } 00573 00574 ErrorCode create_single_entity( Interface* mbImpl, EntityType type ) 00575 { 00576 ErrorCode error; 00577 if( type == MBEDGE ) 00578 { 00579 const double coords[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 }; 00580 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00581 00582 const int conn[] = { 0, 1 }; 00583 const size_t num_elems = sizeof( conn ) / sizeof( conn[0] ) / 2; 00584 00585 std::cout << "Specify verts and ents" << std::endl; 00586 00587 EntityHandle verts[num_vtx], edges[num_elems]; 00588 for( size_t i = 0; i < num_vtx; ++i ) 00589 { 00590 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00591 } 00592 00593 std::cout << "Created vertices" << std::endl; 00594 00595 for( size_t i = 0; i < num_elems; ++i ) 00596 { 00597 EntityHandle c[2]; 00598 c[0] = verts[conn[0]]; 00599 c[1] = verts[conn[1]]; 00600 00601 error = mbImpl->create_element( MBEDGE, c, 2, edges[i] );CHECK_ERR( error ); 00602 } 00603 00604 std::cout << "Created ents" << std::endl; 00605 } 00606 else if( type == MBTRI ) 00607 { 00608 const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0 }; 00609 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00610 00611 const int conn[] = { 0, 1, 2 }; 00612 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3; 00613 00614 EntityHandle verts[num_vtx], faces[num_elems]; 00615 for( size_t i = 0; i < num_vtx; ++i ) 00616 { 00617 error = mbImpl->create_vertex( coords + 3 * i, verts[i] ); 00618 if( error != MB_SUCCESS ) return error; 00619 } 00620 00621 for( size_t i = 0; i < num_elems; ++i ) 00622 { 00623 EntityHandle c[3]; 00624 for( int j = 0; j < 3; j++ ) 00625 c[j] = verts[conn[3 * i + j]]; 00626 00627 error = mbImpl->create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error ); 00628 } 00629 } 00630 else if( type == MBQUAD ) 00631 { 00632 const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0 }; 00633 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00634 00635 const int conn[] = { 0, 1, 2, 3 }; 00636 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3; 00637 00638 EntityHandle verts[num_vtx], faces[num_elems]; 00639 for( size_t i = 0; i < num_vtx; ++i ) 00640 { 00641 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00642 } 00643 00644 for( size_t i = 0; i < num_elems; ++i ) 00645 { 00646 EntityHandle c[4]; 00647 for( int j = 0; j < 4; j++ ) 00648 c[j] = verts[conn[j]]; 00649 00650 error = mbImpl->create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error ); 00651 } 00652 } 00653 else if( type == MBTET ) 00654 { 00655 const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 }; 00656 00657 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00658 00659 const int conn[] = { 0, 1, 2, 3 }; 00660 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4; 00661 00662 EntityHandle verts[num_vtx], cells[num_elems]; 00663 for( size_t i = 0; i < num_vtx; ++i ) 00664 { 00665 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00666 } 00667 00668 for( size_t i = 0; i < num_elems; ++i ) 00669 { 00670 EntityHandle c[4]; 00671 for( int j = 0; j < 4; j++ ) 00672 c[j] = verts[conn[j]]; 00673 00674 error = mbImpl->create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error ); 00675 } 00676 } 00677 else if( type == MBHEX ) 00678 { 00679 const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1 }; 00680 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00681 00682 const int conn[] = { 0, 1, 2, 3, 4, 5, 6, 7 }; 00683 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8; 00684 00685 EntityHandle verts[num_vtx], cells[num_elems]; 00686 for( size_t i = 0; i < num_vtx; ++i ) 00687 { 00688 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00689 } 00690 00691 for( size_t i = 0; i < num_elems; ++i ) 00692 { 00693 EntityHandle c[8]; 00694 for( int j = 0; j < 8; j++ ) 00695 c[j] = verts[conn[j]]; 00696 00697 error = mbImpl->create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error ); 00698 } 00699 } 00700 return MB_SUCCESS; 00701 } 00702 00703 ErrorCode create_mesh( Interface* mbImpl, EntityType type ) 00704 { 00705 ErrorCode error; 00706 if( type == MBEDGE ) 00707 { 00708 const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0 }; 00709 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00710 00711 const int conn[] = { 1, 0, 0, 3, 2, 0, 0, 4 }; 00712 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 2; 00713 00714 EntityHandle verts[num_vtx], edges[num_elems]; 00715 for( size_t i = 0; i < num_vtx; ++i ) 00716 { 00717 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00718 } 00719 00720 for( size_t i = 0; i < num_elems; ++i ) 00721 { 00722 EntityHandle c[2]; 00723 c[0] = verts[conn[2 * i]]; 00724 c[1] = verts[conn[2 * i + 1]]; 00725 00726 error = mbImpl->create_element( MBEDGE, c, 2, edges[i] );CHECK_ERR( error ); 00727 } 00728 } 00729 else if( type == MBTRI ) 00730 { 00731 const double coords[] = { 0, 0, 0, 1, -1, 0, 1, 1, 0, -1, 1, 0, -1, -1, 0, 0, 0, 1 }; 00732 00733 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00734 00735 const int conn[] = { 0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 1, 0, 5, 3, 0, 2, 5, 0, 4, 5, 0, 5, 1 }; 00736 00737 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3; 00738 00739 EntityHandle verts[num_vtx], faces[num_elems]; 00740 for( size_t i = 0; i < num_vtx; ++i ) 00741 { 00742 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00743 } 00744 00745 for( size_t i = 0; i < num_elems; ++i ) 00746 { 00747 EntityHandle c[3]; 00748 for( int j = 0; j < 3; j++ ) 00749 c[j] = verts[conn[3 * i + j]]; 00750 00751 error = mbImpl->create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error ); 00752 } 00753 } 00754 else if( type == MBQUAD ) 00755 { 00756 const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, -1, 1, 0, -1, 0, 0, 00757 -1, -1, 0, 0, -1, 0, 1, -1, 0, 0, 1, 1, 0, 0, 1, 0, -1, 1 }; 00758 00759 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00760 00761 const int conn[] = { 0, 1, 2, 3, 0, 3, 4, 5, 7, 8, 1, 0, 6, 7, 0, 5, 0, 3, 9, 10, 0, 10, 11, 7 }; 00762 00763 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4; 00764 00765 EntityHandle verts[num_vtx], faces[num_elems]; 00766 for( size_t i = 0; i < num_vtx; ++i ) 00767 { 00768 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00769 } 00770 00771 for( size_t i = 0; i < num_elems; ++i ) 00772 { 00773 EntityHandle c[4]; 00774 for( int j = 0; j < 4; j++ ) 00775 c[j] = verts[conn[4 * i + j]]; 00776 00777 error = mbImpl->create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error ); 00778 } 00779 } 00780 else if( type == MBTET ) 00781 { 00782 const double coords[] = { 0, -1, 0, 0, 2, 0, 1, 0, 0, -0.5, 0, 0, 0, 0, 1, 0, 0, -2 }; 00783 00784 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00785 00786 const int conn[] = { 0, 2, 1, 4, 3, 0, 1, 4, 5, 2, 1, 0, 5, 0, 1, 3 }; 00787 00788 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4; 00789 00790 EntityHandle verts[num_vtx], cells[num_elems]; 00791 for( size_t i = 0; i < num_vtx; ++i ) 00792 { 00793 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00794 } 00795 00796 for( size_t i = 0; i < num_elems; ++i ) 00797 { 00798 EntityHandle c[4]; 00799 for( int j = 0; j < 4; j++ ) 00800 c[j] = verts[conn[4 * i + j]]; 00801 00802 error = mbImpl->create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error ); 00803 } 00804 } 00805 else if( type == MBHEX ) 00806 { 00807 const double coords[] = { 0, -1, 0, 1, -1, 0, 1, 1, 0, 0, 1, 0, -1, 1, 0, -1, -1, 0, 00808 0, -1, 1, 1, -1, 1, 1, 1, 1, 0, 1, 1, -1, 1, 1, -1, -1, 1, 00809 0, -1, -1, 1, -1, -1, 1, 1, -1, 0, 1, -1, -1, 1, -1, -1, -1, -1 }; 00810 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00811 00812 const int conn[] = { 0, 1, 2, 3, 6, 7, 8, 9, 5, 0, 3, 4, 11, 6, 9, 10, 00813 12, 13, 14, 15, 0, 1, 2, 3, 17, 12, 15, 16, 5, 0, 3, 4 }; 00814 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8; 00815 00816 EntityHandle verts[num_vtx], cells[num_elems]; 00817 for( size_t i = 0; i < num_vtx; ++i ) 00818 { 00819 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00820 } 00821 00822 for( size_t i = 0; i < num_elems; ++i ) 00823 { 00824 EntityHandle c[8]; 00825 for( int j = 0; j < 8; j++ ) 00826 c[j] = verts[conn[8 * i + j]]; 00827 00828 error = mbImpl->create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error ); 00829 } 00830 } 00831 return MB_SUCCESS; 00832 } 00833 00834 ErrorCode create_simple_mesh( Interface* mbImpl, EntityType type ) 00835 { 00836 ErrorCode error; 00837 if( type == MBEDGE ) 00838 { 00839 const double coords[] = { 0, 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 0, 4, 0, 0, 00840 4, 1, 0, 3, 1, 0, 2, 1, 0, 1, 1, 0, 0, 1, 0 }; 00841 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00842 00843 const int conn[] = { 1, 0, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 0 }; 00844 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 2; 00845 00846 EntityHandle verts[num_vtx], edges[num_elems]; 00847 for( size_t i = 0; i < num_vtx; ++i ) 00848 { 00849 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00850 } 00851 00852 for( size_t i = 0; i < num_elems; ++i ) 00853 { 00854 EntityHandle c[2]; 00855 c[0] = verts[conn[2 * i]]; 00856 c[1] = verts[conn[2 * i + 1]]; 00857 00858 error = mbImpl->create_element( MBEDGE, c, 2, edges[i] );CHECK_ERR( error ); 00859 } 00860 } 00861 else if( type == MBTRI ) 00862 { 00863 const double coords[] = { 0, 0, 0, 1, 0, 0, 2, 0, 0, 2.5, 1, 0, 1.5, 1, 0, 0.5, 1, 00864 0, -0.5, 1, 0, -0.5, -1, 0, 0.5, -1, 0, 1.5, -1, 0, 2.5, -1, 0 }; 00865 00866 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00867 00868 const int conn[] = { 00869 0, 5, 6, 0, 1, 5, 1, 4, 5, 1, 2, 4, 2, 3, 4, 7, 8, 0, 8, 1, 0, 8, 9, 1, 9, 2, 1, 9, 10, 2 00870 }; 00871 00872 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3; 00873 00874 EntityHandle verts[num_vtx], faces[num_elems]; 00875 for( size_t i = 0; i < num_vtx; ++i ) 00876 { 00877 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00878 } 00879 00880 for( size_t i = 0; i < num_elems; ++i ) 00881 { 00882 EntityHandle c[3]; 00883 for( int j = 0; j < 3; j++ ) 00884 c[j] = verts[conn[3 * i + j]]; 00885 00886 error = mbImpl->create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error ); 00887 } 00888 } 00889 else if( type == MBQUAD ) 00890 { 00891 const double coords[] = { 0, 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 0, 4, 0, 0, 5, 0, 0, 0, 1, 0, 1, 1, 0, 2, 1, 0, 00892 3, 1, 0, 4, 1, 0, 5, 1, 0, 0, 2, 0, 1, 2, 0, 2, 2, 0, 3, 2, 0, 4, 2, 0, 5, 2, 0 }; 00893 00894 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00895 00896 const int conn[] = { 0, 1, 7, 6, 1, 2, 8, 7, 2, 3, 9, 8, 3, 4, 10, 9, 4, 5, 11, 10, 00897 6, 7, 13, 12, 7, 8, 14, 13, 8, 9, 15, 14, 9, 10, 16, 15, 10, 11, 17, 16 }; 00898 00899 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4; 00900 00901 EntityHandle verts[num_vtx], faces[num_elems]; 00902 for( size_t i = 0; i < num_vtx; ++i ) 00903 { 00904 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00905 } 00906 00907 for( size_t i = 0; i < num_elems; ++i ) 00908 { 00909 EntityHandle c[4]; 00910 for( int j = 0; j < 4; j++ ) 00911 c[j] = verts[conn[4 * i + j]]; 00912 00913 error = mbImpl->create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error ); 00914 } 00915 } 00916 else if( type == MBTET ) 00917 { 00918 const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1 }; 00919 00920 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00921 00922 const int conn[] = { 0, 1, 2, 5, 3, 0, 2, 5, 4, 1, 0, 5, 4, 0, 3, 5 }; 00923 00924 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4; 00925 00926 EntityHandle verts[num_vtx], cells[num_elems]; 00927 for( size_t i = 0; i < num_vtx; ++i ) 00928 { 00929 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00930 } 00931 00932 for( size_t i = 0; i < num_elems; ++i ) 00933 { 00934 EntityHandle c[4]; 00935 for( int j = 0; j < 4; j++ ) 00936 c[j] = verts[conn[4 * i + j]]; 00937 00938 error = mbImpl->create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error ); 00939 } 00940 } 00941 else if( type == MBHEX ) 00942 { 00943 const double coords[] = { 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 1, 0, 1, 1, 0, 2, 1, 0, 0, 2, 0, 1, 2, 0, 2, 2, 0, 00944 0, 0, 1, 1, 0, 1, 2, 0, 1, 0, 1, 1, 1, 1, 1, 2, 1, 1, 0, 2, 1, 1, 2, 1, 2, 2, 1 }; 00945 const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3; 00946 00947 const int conn[] = { 0, 1, 4, 3, 9, 10, 13, 12, 1, 2, 5, 4, 10, 11, 14, 13, 00948 3, 4, 7, 6, 12, 13, 16, 15, 4, 5, 8, 7, 13, 14, 17, 16 }; 00949 const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8; 00950 00951 EntityHandle verts[num_vtx], cells[num_elems]; 00952 for( size_t i = 0; i < num_vtx; ++i ) 00953 { 00954 error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error ); 00955 } 00956 00957 for( size_t i = 0; i < num_elems; ++i ) 00958 { 00959 EntityHandle c[8]; 00960 for( int j = 0; j < 8; j++ ) 00961 c[j] = verts[conn[8 * i + j]]; 00962 00963 error = mbImpl->create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error ); 00964 } 00965 } 00966 return MB_SUCCESS; 00967 } 00968 00969 ErrorCode test_entities( int mesh_type, EntityType type, int* level_degrees, int num_levels, bool output ) 00970 { 00971 ErrorCode error; 00972 Core mb; 00973 Interface* mbimpl = &mb; 00974 00975 // Create entities 00976 if( mesh_type == 1 ) 00977 { 00978 error = create_single_entity( mbimpl, type ); 00979 if( error != MB_SUCCESS ) return error; 00980 std::cout << "Entity created successfully" << std::endl; 00981 } 00982 else if( mesh_type == 2 ) 00983 { 00984 error = create_mesh( mbimpl, type ); 00985 if( error != MB_SUCCESS ) return error; 00986 std::cout << "Small mesh created successfully" << std::endl; 00987 } 00988 else if( mesh_type == 3 ) 00989 { 00990 error = create_simple_mesh( mbimpl, type ); 00991 if( error != MB_SUCCESS ) return error; 00992 std::cout << "Small simple mesh created successfully" << std::endl; 00993 } 00994 00995 // Generate hierarchy 00996 error = refine_entities( mbimpl, 0, 0, level_degrees, num_levels, output ); 00997 if( error != MB_SUCCESS ) return error; 00998 00999 return MB_SUCCESS; 01000 } 01001 ErrorCode test_1D() 01002 { 01003 ErrorCode error; 01004 01005 std::cout << "Testing EDGE" << std::endl; 01006 EntityType type = MBEDGE; 01007 01008 std::cout << "Testing single entity" << std::endl; 01009 int deg[3] = { 2, 3, 5 }; 01010 int len = sizeof( deg ) / sizeof( int ); 01011 error = test_entities( 1, type, deg, len, false );CHECK_ERR( error ); 01012 01013 std::cout << std::endl; 01014 std::cout << "Testing a small mesh" << std::endl; 01015 error = test_entities( 2, type, deg, len, false );CHECK_ERR( error ); 01016 01017 std::cout << std::endl; 01018 std::cout << "Testing a small simple mesh" << std::endl; 01019 int degree[4] = { 5, 5, 2, 2 }; 01020 len = sizeof( degree ) / sizeof( int ); 01021 error = test_entities( 3, type, degree, len, false );CHECK_ERR( error ); 01022 01023 return MB_SUCCESS; 01024 } 01025 01026 ErrorCode test_2D() 01027 { 01028 ErrorCode error; 01029 01030 std::cout << "Testing TRI" << std::endl; 01031 EntityType type = MBTRI; 01032 01033 std::cout << "Testing single entity" << std::endl; 01034 int deg[3] = { 2, 3, 5 }; 01035 int len = sizeof( deg ) / sizeof( int ); 01036 error = test_entities( 1, type, deg, len, false );CHECK_ERR( error ); 01037 01038 std::cout << std::endl; 01039 std::cout << "Testing a small mesh" << std::endl; 01040 error = test_entities( 2, type, deg, len, false );CHECK_ERR( error ); 01041 01042 std::cout << std::endl; 01043 std::cout << "Testing a small simple mesh" << std::endl; 01044 int degree[2] = { 5, 2 }; 01045 int length = sizeof( degree ) / sizeof( int ); 01046 error = test_entities( 3, type, degree, length, false );CHECK_ERR( error ); 01047 01048 std::cout << std::endl; 01049 std::cout << "Testing QUAD" << std::endl; 01050 type = MBQUAD; 01051 01052 std::cout << "Testing single entity" << std::endl; 01053 error = test_entities( 1, type, deg, len, false );CHECK_ERR( error ); 01054 01055 std::cout << std::endl; 01056 std::cout << "Testing a small mesh" << std::endl; 01057 error = test_entities( 2, type, deg, len, false );CHECK_ERR( error ); 01058 01059 std::cout << std::endl; 01060 std::cout << "Testing a small simple mesh" << std::endl; 01061 error = test_entities( 3, type, degree, length, false );CHECK_ERR( error ); 01062 01063 return MB_SUCCESS; 01064 } 01065 01066 ErrorCode test_3D() 01067 { 01068 ErrorCode error; 01069 01070 std::cout << "Testing TET" << std::endl; 01071 EntityType type = MBTET; 01072 int deg[2] = { 2, 3 }; 01073 int len = sizeof( deg ) / sizeof( int ); 01074 01075 std::cout << "Testing single entity" << std::endl; 01076 error = test_entities( 1, type, deg, len, false );CHECK_ERR( error ); 01077 01078 std::cout << std::endl; 01079 std::cout << "Testing a small mesh" << std::endl; 01080 error = test_entities( 2, type, deg, len, false );CHECK_ERR( error ); 01081 01082 std::cout << std::endl; 01083 std::cout << "Testing a small simple mesh" << std::endl; 01084 int degree[4] = { 2, 2, 2, 2 }; 01085 int length = sizeof( degree ) / sizeof( int ); 01086 error = test_entities( 3, type, degree, length, false );CHECK_ERR( error ); 01087 01088 std::cout << std::endl; 01089 std::cout << "Testing HEX" << std::endl; 01090 type = MBHEX; 01091 01092 std::cout << "Testing single entity" << std::endl; 01093 error = test_entities( 1, type, deg, len, false );CHECK_ERR( error ); 01094 01095 std::cout << std::endl; 01096 std::cout << "Testing a small mesh" << std::endl; 01097 error = test_entities( 2, type, deg, len, false );CHECK_ERR( error ); 01098 01099 std::cout << std::endl; 01100 std::cout << "Testing a small simple mesh" << std::endl; 01101 error = test_entities( 3, type, degree, length, false );CHECK_ERR( error ); 01102 01103 return MB_SUCCESS; 01104 } 01105 01106 ErrorCode test_mesh( const char* filename, int* level_degrees, int num_levels ) 01107 { 01108 Core moab; 01109 Interface* mbImpl = &moab; 01110 ParallelComm* pc = NULL; 01111 EntityHandle fileset; 01112 01113 ErrorCode error; 01114 error = mbImpl->create_meshset( moab::MESHSET_SET, fileset );CHECK_ERR( error ); 01115 01116 #ifdef MOAB_HAVE_MPI 01117 MPI_Comm comm = MPI_COMM_WORLD; 01118 EntityHandle partnset; 01119 error = mbImpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error ); 01120 pc = moab::ParallelComm::get_pcomm( mbImpl, partnset, &comm ); 01121 01122 int procs = 1; 01123 MPI_Comm_size( MPI_COMM_WORLD, &procs ); 01124 01125 if( procs > 1 ) 01126 { 01127 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS"; 01128 01129 error = mbImpl->load_file( filename, &fileset, read_options.c_str() );CHECK_ERR( error ); 01130 01131 // DBG 01132 std::set< unsigned int > shprocs; 01133 error = pc->get_comm_procs( shprocs );CHECK_ERR( error ); 01134 std::cout << "#sprocs = " << shprocs.size() << std::endl; 01135 // DBG 01136 } 01137 else if( procs == 1 ) 01138 { 01139 #endif 01140 error = mbImpl->load_file( filename, &fileset );CHECK_ERR( error ); 01141 #ifdef MOAB_HAVE_MPI 01142 } 01143 #endif 01144 01145 // Generate hierarchy 01146 error = refine_entities( &moab, pc, fileset, level_degrees, num_levels, false );CHECK_ERR( error ); 01147 01148 return MB_SUCCESS; 01149 } 01150 01151 int main( int argc, char* argv[] ) 01152 { 01153 #ifdef MOAB_HAVE_MPI 01154 MPI_Init( &argc, &argv ); 01155 01156 int nprocs, rank; 01157 MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); 01158 MPI_Comm_rank( MPI_COMM_WORLD, &rank ); 01159 #endif 01160 01161 ErrorCode result; 01162 if( argc == 1 ) 01163 { 01164 result = test_1D(); 01165 handle_error_code( result, number_tests_failed, number_tests_successful ); 01166 std::cout << "\n"; 01167 01168 result = test_2D(); 01169 handle_error_code( result, number_tests_failed, number_tests_successful ); 01170 std::cout << "\n"; 01171 01172 result = test_3D(); 01173 handle_error_code( result, number_tests_failed, number_tests_successful ); 01174 std::cout << "\n"; 01175 } 01176 else if( argc == 2 ) 01177 { 01178 const char* filename = argv[1]; 01179 int deg[1] = { 2 }; 01180 int len = sizeof( deg ) / sizeof( int ); 01181 result = test_mesh( filename, deg, len ); 01182 handle_error_code( result, number_tests_failed, number_tests_successful ); 01183 std::cout << "\n"; 01184 } 01185 01186 #ifdef MOAB_HAVE_MPI 01187 MPI_Finalize(); 01188 #endif 01189 01190 return number_tests_failed; 01191 }