MOAB: Mesh Oriented datABase  (version 5.2.1)
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 #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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines