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