MOAB: Mesh Oriented datABase  (version 5.4.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,
00352                            ParallelComm* pc,
00353                            EntityHandle fset,
00354                            int* level_degrees,
00355                            const int num_levels,
00356                            bool output )
00357 {
00358     ErrorCode error;
00359 
00360     // Get the range of entities in the initial mesh
00361     Range init_ents[4];
00362 
00363     int dim[3] = { 1, 2, 3 };
00364 
00365     if( output )
00366     {
00367         //  int inents = init_ents.size();
00368         std::stringstream file;
00369 #ifdef MOAB_HAVE_MPI
00370         file << "MESH_LEVEL_0.h5m";
00371 #else
00372         file << "MESH_LEVEL_0.vtk";
00373 #endif
00374         std::string str         = file.str();
00375         const char* output_file = str.c_str();
00376 #ifdef MOAB_HAVE_MPI
00377         error = mb->write_file( output_file, 0, ";;PARALLEL=WRITE_PART" );CHECK_ERR( error );
00378 #else
00379         error = mb->write_file( output_file, 0, NULL );CHECK_ERR( error );
00380 #endif
00381     }
00382 
00383     // Create an hm object and generate the hierarchy
00384     std::cout << "Creating a hm object" << std::endl;
00385 
00386 #ifdef MOAB_HAVE_MPI
00387     Range averts, aedges, afaces, acells;
00388     error = mb->get_entities_by_dimension( fset, 0, averts );MB_CHK_ERR( error );
00389     error = mb->get_entities_by_dimension( fset, 1, aedges );MB_CHK_ERR( error );
00390     error = mb->get_entities_by_dimension( fset, 2, afaces );MB_CHK_ERR( error );
00391     error = mb->get_entities_by_dimension( fset, 3, acells );MB_CHK_ERR( error );
00392 
00393     /* filter based on parallel status */
00394     if( pc )
00395     {
00396         error = pc->filter_pstatus( averts, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[0] );MB_CHK_ERR( error );
00397         error = pc->filter_pstatus( aedges, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[1] );MB_CHK_ERR( error );
00398         error = pc->filter_pstatus( afaces, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[2] );MB_CHK_ERR( error );
00399         error = pc->filter_pstatus( acells, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[3] );MB_CHK_ERR( error );
00400     }
00401     else
00402     {
00403         init_ents[0] = averts;
00404         init_ents[1] = aedges;
00405         init_ents[2] = afaces;
00406         init_ents[3] = acells;
00407     }
00408 #else
00409     error = mb->get_entities_by_dimension( fset, 0, init_ents[0] );CHECK_ERR( error );
00410     error = mb->get_entities_by_dimension( fset, 1, init_ents[1] );CHECK_ERR( error );
00411     error = mb->get_entities_by_dimension( fset, 2, init_ents[2] );CHECK_ERR( error );
00412     error = mb->get_entities_by_dimension( fset, 3, init_ents[3] );CHECK_ERR( error );
00413 #endif
00414 
00415     NestedRefine uref( dynamic_cast< Core* >( mb ), pc, fset );
00416     std::vector< EntityHandle > set;
00417 
00418     std::cout << "Starting hierarchy generation" << std::endl;
00419     bool opt = true;
00420     // bool opt = false;
00421     error = uref.generate_mesh_hierarchy( num_levels, level_degrees, set, opt );CHECK_ERR( error );
00422     std::cout << "Finished hierarchy generation in " << uref.timeall.tm_total << "  secs" << std::endl;
00423 #ifdef MOAB_HAVE_MPI
00424     if( pc )
00425     {
00426         std::cout << "Time taken for refinement " << uref.timeall.tm_refine << "  secs" << std::endl;
00427         std::cout << "Time taken for resolving shared interface " << uref.timeall.tm_resolve << "  secs" << std::endl;
00428     }
00429 #endif
00430 
00431     // error = uref.exchange_ghosts(set, 1); CHECK_ERR(error);
00432 
00433     std::cout << std::endl;
00434     std::cout << "Mesh size for level 0  :: inverts = " << init_ents[0].size() << ", inedges = " << init_ents[1].size()
00435               << ", infaces = " << init_ents[2].size() << ", incells = " << init_ents[3].size() << std::endl;
00436 
00437     Range prev_ents[4];
00438     for( int i = 0; i < 4; i++ )
00439         prev_ents[i] = init_ents[i];
00440 
00441     // Loop over each mesh level and check its topological properties
00442     for( int l = 0; l < num_levels; l++ )
00443     {
00444         Range all_ents;
00445         error = mb->get_entities_by_handle( set[l + 1], all_ents );CHECK_ERR( error );
00446 
00447         Range ents[4];
00448         for( int k = 0; k < 4; k++ )
00449             ents[k] = all_ents.subset_by_dimension( k );
00450 
00451         if( ents[0].empty() || all_ents.empty() ) std::cout << "Something is not right" << std::endl;
00452 
00453         std::cout << std::endl;
00454         std::cout << "Mesh size for level " << l + 1 << "  :: nverts = " << ents[0].size()
00455                   << ", nedges = " << ents[1].size() << ", nfaces = " << ents[2].size()
00456                   << ", ncells = " << ents[3].size() << std::endl;
00457 
00458         // Check if the number of new entities created are correct.
00459 
00460         for( int type = 1; type < 3; type++ )
00461         {
00462             int factor = 1;
00463             if( !ents[type + 1].empty() )
00464             {
00465 
00466                 for( int p = 0; p <= l; p++ )
00467                 {
00468                     for( int d = 0; d < dim[type]; d++ )
00469                         factor *= level_degrees[p];
00470                 }
00471                 int expected_nents = factor * init_ents[type + 1].size();
00472                 CHECK_EQUAL( expected_nents, (int)ents[type + 1].size() );
00473             }
00474         }
00475 
00476         // Check adjacencies
00477         error = test_adjacencies( mb, &uref, all_ents );CHECK_ERR( error );
00478 
00479         // Check interlevel child-parent query between previous and current level
00480         for( int type = 1; type < 3; type++ )
00481         {
00482             if( !prev_ents[type + 1].empty() )
00483             {
00484                 for( Range::iterator e = prev_ents[type + 1].begin(); e != prev_ents[type + 1].end(); e++ )
00485                 {
00486                     std::vector< EntityHandle > children;
00487                     error = uref.parent_to_child( *e, l, l + 1, children );CHECK_ERR( error );
00488                     for( int i = 0; i < (int)children.size(); i++ )
00489                     {
00490                         EntityHandle parent;
00491                         error = uref.child_to_parent( children[i], l + 1, l, &parent );CHECK_ERR( error );
00492                         assert( parent == *e );
00493                     }
00494                 }
00495             }
00496         }
00497 
00498         for( int i = 0; i < 4; i++ )
00499             prev_ents[i] = ents[i];
00500 
00501         // Print out the boundary vertices
00502         /*  EntityHandle bnd_set;
00503           error = mb->create_meshset(MESHSET_SET, bnd_set); MB_CHK_ERR(error);
00504           std::vector<EntityHandle> vbnd;
00505 
00506           for (int k=0; k<3; k++){
00507               std::cout<<"Finding boundary of dimension "<<k<<" with size
00508           "<<ents[k].size()<<std::endl; if (ents[k].size() != 0){ for (Range::iterator v =
00509           ents[k].begin(); v != ents[k].end(); v++)
00510                     {
00511                       EntityHandle ent = *v;
00512                       bool bnd = uref.is_entity_on_boundary(ent);
00513                       if (bnd)
00514                         vbnd.push_back(*v);
00515                     }
00516                 }
00517             }
00518 
00519           std::cout<<"vbnd.size = "<<vbnd.size()<<std::endl;
00520           error = mb->add_entities(bnd_set, &vbnd[0], (int)vbnd.size()); MB_CHK_ERR(error);
00521           if (output)
00522             {
00523               std::stringstream file;
00524               file <<  "VBND_LEVEL_" <<l+1<<".vtk";
00525               std::string str = file.str();
00526               const char* output_file = str.c_str();
00527               char * write_opts = NULL;
00528               error = mb->write_file(output_file, 0, write_opts, &bnd_set, 1); CHECK_ERR(error);
00529             }*/
00530 
00531         // Print out the mesh
00532         if( output )
00533         {
00534             std::stringstream file;
00535 
00536 #ifdef MOAB_HAVE_MPI
00537             file << "MESH_LEVEL_" << l + 1 << ".h5m";
00538 #else
00539             file << "MESH_LEVEL_" << l + 1 << ".vtk";
00540 #endif
00541             std::string str         = file.str();
00542             const char* output_file = str.c_str();
00543 #ifdef MOAB_HAVE_MPI
00544             error = mb->write_file( output_file, 0, ";;PARALLEL=WRITE_PART", &set[l + 1], 1 );CHECK_ERR( error );
00545 #else
00546             error = mb->write_file( output_file, 0, NULL, &set[l + 1], 1 );CHECK_ERR( error );
00547 #endif
00548         }
00549     }
00550 
00551     // Check interlevel child-parent query between initial and most refined mesh
00552     for( int type = 1; type < 3; type++ )
00553     {
00554         if( !init_ents[type + 1].empty() )
00555         {
00556             for( Range::iterator e = init_ents[type + 1].begin(); e != init_ents[type + 1].end(); e++ )
00557             {
00558                 std::vector< EntityHandle > children;
00559                 error = uref.parent_to_child( *e, 0, num_levels, children );CHECK_ERR( error );
00560                 for( int i = 0; i < (int)children.size(); i++ )
00561                 {
00562                     EntityHandle parent;
00563                     error = uref.child_to_parent( children[i], num_levels, 0, &parent );CHECK_ERR( error );
00564                     assert( parent == *e );
00565                 }
00566             }
00567         }
00568     }
00569 
00570     // Print out the whole hierarchy into a single file
00571     /* if (output)
00572        {
00573          std::stringstream file;
00574          file <<  "MESH_HIERARCHY.vtk";
00575          std::string str = file.str();
00576          const char* output_file = str.c_str();
00577          error = mb->write_file(output_file); CHECK_ERR(error);
00578        }*/
00579 
00580     return MB_SUCCESS;
00581 }
00582 
00583 ErrorCode create_single_entity( Interface* mbImpl, EntityType type )
00584 {
00585     ErrorCode error;
00586     if( type == MBEDGE )
00587     {
00588         const double coords[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 };
00589         const size_t num_vtx  = sizeof( coords ) / sizeof( double ) / 3;
00590 
00591         const int conn[]       = { 0, 1 };
00592         const size_t num_elems = sizeof( conn ) / sizeof( conn[0] ) / 2;
00593 
00594         std::cout << "Specify verts and ents" << std::endl;
00595 
00596         EntityHandle verts[num_vtx], edges[num_elems];
00597         for( size_t i = 0; i < num_vtx; ++i )
00598         {
00599             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00600         }
00601 
00602         std::cout << "Created vertices" << std::endl;
00603 
00604         for( size_t i = 0; i < num_elems; ++i )
00605         {
00606             EntityHandle c[2];
00607             c[0] = verts[conn[0]];
00608             c[1] = verts[conn[1]];
00609 
00610             error = mbImpl->create_element( MBEDGE, c, 2, edges[i] );CHECK_ERR( error );
00611         }
00612 
00613         std::cout << "Created ents" << std::endl;
00614     }
00615     else if( type == MBTRI )
00616     {
00617         const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0 };
00618         const size_t num_vtx  = sizeof( coords ) / sizeof( double ) / 3;
00619 
00620         const int conn[]       = { 0, 1, 2 };
00621         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3;
00622 
00623         EntityHandle verts[num_vtx], faces[num_elems];
00624         for( size_t i = 0; i < num_vtx; ++i )
00625         {
00626             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );
00627             if( error != MB_SUCCESS ) return error;
00628         }
00629 
00630         for( size_t i = 0; i < num_elems; ++i )
00631         {
00632             EntityHandle c[3];
00633             for( int j = 0; j < 3; j++ )
00634                 c[j] = verts[conn[3 * i + j]];
00635 
00636             error = mbImpl->create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error );
00637         }
00638     }
00639     else if( type == MBQUAD )
00640     {
00641         const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0 };
00642         const size_t num_vtx  = sizeof( coords ) / sizeof( double ) / 3;
00643 
00644         const int conn[]       = { 0, 1, 2, 3 };
00645         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3;
00646 
00647         EntityHandle verts[num_vtx], faces[num_elems];
00648         for( size_t i = 0; i < num_vtx; ++i )
00649         {
00650             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00651         }
00652 
00653         for( size_t i = 0; i < num_elems; ++i )
00654         {
00655             EntityHandle c[4];
00656             for( int j = 0; j < 4; j++ )
00657                 c[j] = verts[conn[j]];
00658 
00659             error = mbImpl->create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error );
00660         }
00661     }
00662     else if( type == MBTET )
00663     {
00664         const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 };
00665 
00666         const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
00667 
00668         const int conn[]       = { 0, 1, 2, 3 };
00669         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
00670 
00671         EntityHandle verts[num_vtx], cells[num_elems];
00672         for( size_t i = 0; i < num_vtx; ++i )
00673         {
00674             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00675         }
00676 
00677         for( size_t i = 0; i < num_elems; ++i )
00678         {
00679             EntityHandle c[4];
00680             for( int j = 0; j < 4; j++ )
00681                 c[j] = verts[conn[j]];
00682 
00683             error = mbImpl->create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error );
00684         }
00685     }
00686     else if( type == MBHEX )
00687     {
00688         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 };
00689         const size_t num_vtx  = sizeof( coords ) / sizeof( double ) / 3;
00690 
00691         const int conn[]       = { 0, 1, 2, 3, 4, 5, 6, 7 };
00692         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8;
00693 
00694         EntityHandle verts[num_vtx], cells[num_elems];
00695         for( size_t i = 0; i < num_vtx; ++i )
00696         {
00697             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00698         }
00699 
00700         for( size_t i = 0; i < num_elems; ++i )
00701         {
00702             EntityHandle c[8];
00703             for( int j = 0; j < 8; j++ )
00704                 c[j] = verts[conn[j]];
00705 
00706             error = mbImpl->create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error );
00707         }
00708     }
00709     return MB_SUCCESS;
00710 }
00711 
00712 ErrorCode create_mesh( Interface* mbImpl, EntityType type )
00713 {
00714     ErrorCode error;
00715     if( type == MBEDGE )
00716     {
00717         const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0 };
00718         const size_t num_vtx  = sizeof( coords ) / sizeof( double ) / 3;
00719 
00720         const int conn[]       = { 1, 0, 0, 3, 2, 0, 0, 4 };
00721         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 2;
00722 
00723         EntityHandle verts[num_vtx], edges[num_elems];
00724         for( size_t i = 0; i < num_vtx; ++i )
00725         {
00726             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00727         }
00728 
00729         for( size_t i = 0; i < num_elems; ++i )
00730         {
00731             EntityHandle c[2];
00732             c[0] = verts[conn[2 * i]];
00733             c[1] = verts[conn[2 * i + 1]];
00734 
00735             error = mbImpl->create_element( MBEDGE, c, 2, edges[i] );CHECK_ERR( error );
00736         }
00737     }
00738     else if( type == MBTRI )
00739     {
00740         const double coords[] = { 0, 0, 0, 1, -1, 0, 1, 1, 0, -1, 1, 0, -1, -1, 0, 0, 0, 1 };
00741 
00742         const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
00743 
00744         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 };
00745 
00746         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3;
00747 
00748         EntityHandle verts[num_vtx], faces[num_elems];
00749         for( size_t i = 0; i < num_vtx; ++i )
00750         {
00751             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00752         }
00753 
00754         for( size_t i = 0; i < num_elems; ++i )
00755         {
00756             EntityHandle c[3];
00757             for( int j = 0; j < 3; j++ )
00758                 c[j] = verts[conn[3 * i + j]];
00759 
00760             error = mbImpl->create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error );
00761         }
00762     }
00763     else if( type == MBQUAD )
00764     {
00765         const double coords[] = { 0,  0,  0, 1, 0,  0, 1, 1,  0, 0, 1, 0, -1, 1, 0, -1, 0,  0,
00766                                   -1, -1, 0, 0, -1, 0, 1, -1, 0, 0, 1, 1, 0,  0, 1, 0,  -1, 1 };
00767 
00768         const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
00769 
00770         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 };
00771 
00772         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
00773 
00774         EntityHandle verts[num_vtx], faces[num_elems];
00775         for( size_t i = 0; i < num_vtx; ++i )
00776         {
00777             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00778         }
00779 
00780         for( size_t i = 0; i < num_elems; ++i )
00781         {
00782             EntityHandle c[4];
00783             for( int j = 0; j < 4; j++ )
00784                 c[j] = verts[conn[4 * i + j]];
00785 
00786             error = mbImpl->create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error );
00787         }
00788     }
00789     else if( type == MBTET )
00790     {
00791         const double coords[] = { 0, -1, 0, 0, 2, 0, 1, 0, 0, -0.5, 0, 0, 0, 0, 1, 0, 0, -2 };
00792 
00793         const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
00794 
00795         const int conn[] = { 0, 2, 1, 4, 3, 0, 1, 4, 5, 2, 1, 0, 5, 0, 1, 3 };
00796 
00797         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
00798 
00799         EntityHandle verts[num_vtx], cells[num_elems];
00800         for( size_t i = 0; i < num_vtx; ++i )
00801         {
00802             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00803         }
00804 
00805         for( size_t i = 0; i < num_elems; ++i )
00806         {
00807             EntityHandle c[4];
00808             for( int j = 0; j < 4; j++ )
00809                 c[j] = verts[conn[4 * i + j]];
00810 
00811             error = mbImpl->create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error );
00812         }
00813     }
00814     else if( type == MBHEX )
00815     {
00816         const double coords[] = { 0, -1, 0,  1, -1, 0,  1, 1, 0,  0, 1, 0,  -1, 1, 0,  -1, -1, 0,
00817                                   0, -1, 1,  1, -1, 1,  1, 1, 1,  0, 1, 1,  -1, 1, 1,  -1, -1, 1,
00818                                   0, -1, -1, 1, -1, -1, 1, 1, -1, 0, 1, -1, -1, 1, -1, -1, -1, -1 };
00819         const size_t num_vtx  = sizeof( coords ) / sizeof( double ) / 3;
00820 
00821         const int conn[]       = { 0,  1,  2,  3,  6, 7, 8, 9, 5,  0,  3,  4,  11, 6, 9, 10,
00822                              12, 13, 14, 15, 0, 1, 2, 3, 17, 12, 15, 16, 5,  0, 3, 4 };
00823         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8;
00824 
00825         EntityHandle verts[num_vtx], cells[num_elems];
00826         for( size_t i = 0; i < num_vtx; ++i )
00827         {
00828             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00829         }
00830 
00831         for( size_t i = 0; i < num_elems; ++i )
00832         {
00833             EntityHandle c[8];
00834             for( int j = 0; j < 8; j++ )
00835                 c[j] = verts[conn[8 * i + j]];
00836 
00837             error = mbImpl->create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error );
00838         }
00839     }
00840     return MB_SUCCESS;
00841 }
00842 
00843 ErrorCode create_simple_mesh( Interface* mbImpl, EntityType type )
00844 {
00845     ErrorCode error;
00846     if( type == MBEDGE )
00847     {
00848         const double coords[] = { 0, 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 0, 4, 0, 0,
00849                                   4, 1, 0, 3, 1, 0, 2, 1, 0, 1, 1, 0, 0, 1, 0 };
00850         const size_t num_vtx  = sizeof( coords ) / sizeof( double ) / 3;
00851 
00852         const int conn[]       = { 1, 0, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 0 };
00853         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 2;
00854 
00855         EntityHandle verts[num_vtx], edges[num_elems];
00856         for( size_t i = 0; i < num_vtx; ++i )
00857         {
00858             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00859         }
00860 
00861         for( size_t i = 0; i < num_elems; ++i )
00862         {
00863             EntityHandle c[2];
00864             c[0] = verts[conn[2 * i]];
00865             c[1] = verts[conn[2 * i + 1]];
00866 
00867             error = mbImpl->create_element( MBEDGE, c, 2, edges[i] );CHECK_ERR( error );
00868         }
00869     }
00870     else if( type == MBTRI )
00871     {
00872         const double coords[] = { 0, 0,    0, 1, 0,    0,  2, 0,   0,  2.5, 1,   0,  1.5, 1,   0,  0.5, 1,
00873                                   0, -0.5, 1, 0, -0.5, -1, 0, 0.5, -1, 0,   1.5, -1, 0,   2.5, -1, 0 };
00874 
00875         const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
00876 
00877         const int conn[] = { 0, 5, 6, 0, 1, 5, 1, 4, 5, 1, 2, 4, 2, 3,  4,
00878                              7, 8, 0, 8, 1, 0, 8, 9, 1, 9, 2, 1, 9, 10, 2 };
00879 
00880         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3;
00881 
00882         EntityHandle verts[num_vtx], faces[num_elems];
00883         for( size_t i = 0; i < num_vtx; ++i )
00884         {
00885             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00886         }
00887 
00888         for( size_t i = 0; i < num_elems; ++i )
00889         {
00890             EntityHandle c[3];
00891             for( int j = 0; j < 3; j++ )
00892                 c[j] = verts[conn[3 * i + j]];
00893 
00894             error = mbImpl->create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error );
00895         }
00896     }
00897     else if( type == MBQUAD )
00898     {
00899         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,
00900                                   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 };
00901 
00902         const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
00903 
00904         const int conn[] = { 0, 1, 7,  6,  1, 2, 8,  7,  2, 3, 9,  8,  3, 4,  10, 9,  4,  5,  11, 10,
00905                              6, 7, 13, 12, 7, 8, 14, 13, 8, 9, 15, 14, 9, 10, 16, 15, 10, 11, 17, 16 };
00906 
00907         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
00908 
00909         EntityHandle verts[num_vtx], faces[num_elems];
00910         for( size_t i = 0; i < num_vtx; ++i )
00911         {
00912             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00913         }
00914 
00915         for( size_t i = 0; i < num_elems; ++i )
00916         {
00917             EntityHandle c[4];
00918             for( int j = 0; j < 4; j++ )
00919                 c[j] = verts[conn[4 * i + j]];
00920 
00921             error = mbImpl->create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error );
00922         }
00923     }
00924     else if( type == MBTET )
00925     {
00926         const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1 };
00927 
00928         const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
00929 
00930         const int conn[] = { 0, 1, 2, 5, 3, 0, 2, 5, 4, 1, 0, 5, 4, 0, 3, 5 };
00931 
00932         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
00933 
00934         EntityHandle verts[num_vtx], cells[num_elems];
00935         for( size_t i = 0; i < num_vtx; ++i )
00936         {
00937             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00938         }
00939 
00940         for( size_t i = 0; i < num_elems; ++i )
00941         {
00942             EntityHandle c[4];
00943             for( int j = 0; j < 4; j++ )
00944                 c[j] = verts[conn[4 * i + j]];
00945 
00946             error = mbImpl->create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error );
00947         }
00948     }
00949     else if( type == MBHEX )
00950     {
00951         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,
00952                                   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 };
00953         const size_t num_vtx  = sizeof( coords ) / sizeof( double ) / 3;
00954 
00955         const int conn[]       = { 0, 1, 4, 3, 9,  10, 13, 12, 1, 2, 5, 4, 10, 11, 14, 13,
00956                              3, 4, 7, 6, 12, 13, 16, 15, 4, 5, 8, 7, 13, 14, 17, 16 };
00957         const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8;
00958 
00959         EntityHandle verts[num_vtx], cells[num_elems];
00960         for( size_t i = 0; i < num_vtx; ++i )
00961         {
00962             error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
00963         }
00964 
00965         for( size_t i = 0; i < num_elems; ++i )
00966         {
00967             EntityHandle c[8];
00968             for( int j = 0; j < 8; j++ )
00969                 c[j] = verts[conn[8 * i + j]];
00970 
00971             error = mbImpl->create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error );
00972         }
00973     }
00974     return MB_SUCCESS;
00975 }
00976 
00977 ErrorCode test_entities( int mesh_type, EntityType type, int* level_degrees, int num_levels, bool output )
00978 {
00979     ErrorCode error;
00980     Core mb;
00981     Interface* mbimpl = &mb;
00982 
00983     // Create entities
00984     if( mesh_type == 1 )
00985     {
00986         error = create_single_entity( mbimpl, type );
00987         if( error != MB_SUCCESS ) return error;
00988         std::cout << "Entity created successfully" << std::endl;
00989     }
00990     else if( mesh_type == 2 )
00991     {
00992         error = create_mesh( mbimpl, type );
00993         if( error != MB_SUCCESS ) return error;
00994         std::cout << "Small mesh created successfully" << std::endl;
00995     }
00996     else if( mesh_type == 3 )
00997     {
00998         error = create_simple_mesh( mbimpl, type );
00999         if( error != MB_SUCCESS ) return error;
01000         std::cout << "Small simple mesh created successfully" << std::endl;
01001     }
01002 
01003     // Generate hierarchy
01004     error = refine_entities( mbimpl, 0, 0, level_degrees, num_levels, output );
01005     if( error != MB_SUCCESS ) return error;
01006 
01007     return MB_SUCCESS;
01008 }
01009 ErrorCode test_1D()
01010 {
01011     ErrorCode error;
01012 
01013     std::cout << "Testing EDGE" << std::endl;
01014     EntityType type = MBEDGE;
01015 
01016     std::cout << "Testing single entity" << std::endl;
01017     int deg[3] = { 2, 3, 5 };
01018     int len    = sizeof( deg ) / sizeof( int );
01019     error      = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
01020 
01021     std::cout << std::endl;
01022     std::cout << "Testing a small mesh" << std::endl;
01023     error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
01024 
01025     std::cout << std::endl;
01026     std::cout << "Testing a small simple mesh" << std::endl;
01027     int degree[4] = { 5, 5, 2, 2 };
01028     len           = sizeof( degree ) / sizeof( int );
01029     error         = test_entities( 3, type, degree, len, false );CHECK_ERR( error );
01030 
01031     return MB_SUCCESS;
01032 }
01033 
01034 ErrorCode test_2D()
01035 {
01036     ErrorCode error;
01037 
01038     std::cout << "Testing TRI" << std::endl;
01039     EntityType type = MBTRI;
01040 
01041     std::cout << "Testing single entity" << std::endl;
01042     int deg[3] = { 2, 3, 5 };
01043     int len    = sizeof( deg ) / sizeof( int );
01044     error      = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
01045 
01046     std::cout << std::endl;
01047     std::cout << "Testing a small mesh" << std::endl;
01048     error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
01049 
01050     std::cout << std::endl;
01051     std::cout << "Testing a small simple mesh" << std::endl;
01052     int degree[2] = { 5, 2 };
01053     int length    = sizeof( degree ) / sizeof( int );
01054     error         = test_entities( 3, type, degree, length, false );CHECK_ERR( error );
01055 
01056     std::cout << std::endl;
01057     std::cout << "Testing QUAD" << std::endl;
01058     type = MBQUAD;
01059 
01060     std::cout << "Testing single entity" << std::endl;
01061     error = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
01062 
01063     std::cout << std::endl;
01064     std::cout << "Testing a small mesh" << std::endl;
01065     error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
01066 
01067     std::cout << std::endl;
01068     std::cout << "Testing a small simple mesh" << std::endl;
01069     error = test_entities( 3, type, degree, length, false );CHECK_ERR( error );
01070 
01071     return MB_SUCCESS;
01072 }
01073 
01074 ErrorCode test_3D()
01075 {
01076     ErrorCode error;
01077 
01078     std::cout << "Testing TET" << std::endl;
01079     EntityType type = MBTET;
01080     int deg[2]      = { 2, 3 };
01081     int len         = sizeof( deg ) / sizeof( int );
01082 
01083     std::cout << "Testing single entity" << std::endl;
01084     error = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
01085 
01086     std::cout << std::endl;
01087     std::cout << "Testing a small mesh" << std::endl;
01088     error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
01089 
01090     std::cout << std::endl;
01091     std::cout << "Testing a small simple mesh" << std::endl;
01092     int degree[4] = { 2, 2, 2, 2 };
01093     int length    = sizeof( degree ) / sizeof( int );
01094     error         = test_entities( 3, type, degree, length, false );CHECK_ERR( error );
01095 
01096     std::cout << std::endl;
01097     std::cout << "Testing HEX" << std::endl;
01098     type = MBHEX;
01099 
01100     std::cout << "Testing single entity" << std::endl;
01101     error = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
01102 
01103     std::cout << std::endl;
01104     std::cout << "Testing a small mesh" << std::endl;
01105     error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
01106 
01107     std::cout << std::endl;
01108     std::cout << "Testing a small simple mesh" << std::endl;
01109     error = test_entities( 3, type, degree, length, false );CHECK_ERR( error );
01110 
01111     return MB_SUCCESS;
01112 }
01113 
01114 ErrorCode test_mesh( const char* filename, int* level_degrees, int num_levels )
01115 {
01116     Core moab;
01117     Interface* mbImpl = &moab;
01118     ParallelComm* pc  = NULL;
01119     EntityHandle fileset;
01120 
01121     ErrorCode error;
01122     error = mbImpl->create_meshset( moab::MESHSET_SET, fileset );CHECK_ERR( error );
01123 
01124 #ifdef MOAB_HAVE_MPI
01125     MPI_Comm comm = MPI_COMM_WORLD;
01126     EntityHandle partnset;
01127     error = mbImpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
01128     pc = moab::ParallelComm::get_pcomm( mbImpl, partnset, &comm );
01129 
01130     int procs = 1;
01131     MPI_Comm_size( MPI_COMM_WORLD, &procs );
01132 
01133     if( procs > 1 )
01134     {
01135         read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS";
01136 
01137         error = mbImpl->load_file( filename, &fileset, read_options.c_str() );CHECK_ERR( error );
01138 
01139         // DBG
01140         std::set< unsigned int > shprocs;
01141         error = pc->get_comm_procs( shprocs );CHECK_ERR( error );
01142         std::cout << "#sprocs = " << shprocs.size() << std::endl;
01143         // DBG
01144     }
01145     else if( procs == 1 )
01146     {
01147 #endif
01148         error = mbImpl->load_file( filename, &fileset );CHECK_ERR( error );
01149 #ifdef MOAB_HAVE_MPI
01150     }
01151 #endif
01152 
01153     // Generate hierarchy
01154     error = refine_entities( &moab, pc, fileset, level_degrees, num_levels, false );CHECK_ERR( error );
01155 
01156     return MB_SUCCESS;
01157 }
01158 
01159 int main( int argc, char* argv[] )
01160 {
01161 #ifdef MOAB_HAVE_MPI
01162     MPI_Init( &argc, &argv );
01163 
01164     int nprocs, rank;
01165     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
01166     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
01167 #endif
01168 
01169     ErrorCode result;
01170     if( argc == 1 )
01171     {
01172         result = test_1D();
01173         handle_error_code( result, number_tests_failed, number_tests_successful );
01174         std::cout << "\n";
01175 
01176         result = test_2D();
01177         handle_error_code( result, number_tests_failed, number_tests_successful );
01178         std::cout << "\n";
01179 
01180         result = test_3D();
01181         handle_error_code( result, number_tests_failed, number_tests_successful );
01182         std::cout << "\n";
01183     }
01184     else if( argc == 2 )
01185     {
01186         const char* filename = argv[1];
01187         int deg[1]           = { 2 };
01188         int len              = sizeof( deg ) / sizeof( int );
01189         result               = test_mesh( filename, deg, len );
01190         handle_error_code( result, number_tests_failed, number_tests_successful );
01191         std::cout << "\n";
01192     }
01193 
01194 #ifdef MOAB_HAVE_MPI
01195     MPI_Finalize();
01196 #endif
01197 
01198     return number_tests_failed;
01199 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines