MOAB: Mesh Oriented datABase  (version 5.2.1)
adj_moab_test.cpp
Go to the documentation of this file.
00001 /*This function tests the AHF datastructures on CST meshes*/
00002 #include <iostream>
00003 #include <vector>
00004 #include <algorithm>
00005 #include "moab/Core.hpp"
00006 #include "moab/Range.hpp"
00007 #include "moab/MeshTopoUtil.hpp"
00008 #include "moab/HalfFacetRep.hpp"
00009 #include "TestUtil.hpp"
00010 
00011 #ifdef MOAB_HAVE_MPI
00012 #include "moab/ParallelComm.hpp"
00013 #include "MBParallelConventions.h"
00014 #include "moab/FileOptions.hpp"
00015 #include "MBTagConventions.hpp"
00016 #include "moab_mpi.h"
00017 #endif
00018 
00019 using namespace moab;
00020 
00021 #ifdef MOAB_HAVE_MPI
00022 std::string read_options;
00023 #endif
00024 
00025 int number_tests_successful = 0;
00026 int number_tests_failed     = 0;
00027 
00028 void handle_error_code( ErrorCode rv, int& number_failed, int& number_successful )
00029 {
00030     if( rv == MB_SUCCESS )
00031     {
00032 #ifdef MOAB_HAVE_MPI
00033         int rank = 0;
00034         MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00035         if( rank == 0 ) std::cout << "Success";
00036 #else
00037         std::cout << "Success";
00038 #endif
00039         number_successful++;
00040     }
00041     else
00042     {
00043         std::cout << "Failure";
00044         number_failed++;
00045     }
00046 }
00047 
00048 ErrorCode ahf_test( const char* filename )
00049 {
00050 
00051     Core moab;
00052     Interface* mbImpl = &moab;
00053     ParallelComm* pc  = NULL;
00054     MeshTopoUtil mtu( mbImpl );
00055     ErrorCode error;
00056     EntityHandle fileset;
00057     error = mbImpl->create_meshset( moab::MESHSET_SET, fileset );CHECK_ERR( error );
00058 
00059 #ifdef MOAB_HAVE_MPI
00060     int procs = 1;
00061     MPI_Comm_size( MPI_COMM_WORLD, &procs );
00062 
00063     if( procs > 1 )
00064     {
00065         read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
00066 
00067         error = mbImpl->load_file( filename, &fileset, read_options.c_str() );CHECK_ERR( error );
00068     }
00069     else if( procs == 1 )
00070     {
00071 #endif
00072         error = mbImpl->load_file( filename, &fileset );CHECK_ERR( error );
00073 #ifdef MOAB_HAVE_MPI
00074     }
00075 #endif
00076 
00077     /*Create ranges for handles of explicit elements of the mixed mesh*/
00078     Range verts, edges, faces, cells;
00079     error = mbImpl->get_entities_by_dimension( fileset, 0, verts );CHECK_ERR( error );
00080     error = mbImpl->get_entities_by_dimension( fileset, 1, edges );CHECK_ERR( error );
00081     error = mbImpl->get_entities_by_dimension( fileset, 2, faces );CHECK_ERR( error );
00082     error = mbImpl->get_entities_by_dimension( fileset, 3, cells );CHECK_ERR( error );
00083 
00084     // Create an ahf instance
00085 #ifdef MOAB_HAVE_MPI
00086     pc = ParallelComm::get_pcomm( mbImpl, 0 );
00087     if( !pc ) pc = new moab::ParallelComm( &moab, MPI_COMM_WORLD );
00088 #endif
00089 
00090     HalfFacetRep ahf( &moab, pc, fileset );
00091 
00092     // Call the initialize function which creates the maps for each dimension
00093     error = ahf.initialize();CHECK_ERR( error );
00094 
00095     std::cout << "Finished AHF initialization" << std::endl;
00096 
00097     // Perform queries
00098     std::vector< EntityHandle > adjents;
00099     Range mbents, ahfents;
00100 
00101     // 1D Queries //
00102     // IQ1: For every vertex, obtain incident edges
00103     if( edges.size() )
00104     {
00105         Range everts;
00106         error = mbImpl->get_connectivity( edges, everts, true );CHECK_ERR( error );
00107         for( Range::iterator i = everts.begin(); i != everts.end(); ++i )
00108         {
00109             adjents.clear();
00110             error = ahf.get_up_adjacencies( *i, 1, adjents );CHECK_ERR( error );
00111             mbents.clear();
00112             error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error );
00113 
00114             CHECK_EQUAL( mbents.size(), adjents.size() );
00115             std::sort( adjents.begin(), adjents.end() );
00116             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00117             mbents = subtract( mbents, ahfents );
00118             CHECK( !mbents.size() );
00119         }
00120     }
00121 
00122     // NQ1:  For every edge, obtain neighbor edges
00123     if( edges.size() )
00124     {
00125         for( Range::iterator i = edges.begin(); i != edges.end(); ++i )
00126         {
00127             adjents.clear();
00128             error = ahf.get_neighbor_adjacencies( *i, adjents );CHECK_ERR( error );
00129             mbents.clear();
00130             error = mtu.get_bridge_adjacencies( *i, 0, 1, mbents );CHECK_ERR( error );
00131 
00132             CHECK_EQUAL( mbents.size(), adjents.size() );
00133             std::sort( adjents.begin(), adjents.end() );
00134             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00135             mbents = subtract( mbents, ahfents );
00136             CHECK( !mbents.size() );
00137         }
00138     }
00139 
00140     std::cout << "Finished 1D queries" << std::endl;
00141 
00142     // 2D Queries
00143     // IQ21: For every vertex, obtain incident faces
00144     if( faces.size() )
00145     {
00146         Range fverts;
00147         error = mbImpl->get_connectivity( faces, fverts, true );CHECK_ERR( error );
00148         for( Range::iterator i = fverts.begin(); i != fverts.end(); ++i )
00149         {
00150             adjents.clear();
00151             error = ahf.get_up_adjacencies( *i, 2, adjents );CHECK_ERR( error );
00152             mbents.clear();
00153             error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error );
00154 
00155             CHECK_EQUAL( mbents.size(), adjents.size() );
00156             std::sort( adjents.begin(), adjents.end() );
00157             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00158             mbents = subtract( mbents, ahfents );
00159             CHECK( !mbents.size() );
00160         }
00161     }
00162 
00163     // IQ22: For every edge, obtain incident faces
00164     if( edges.size() && faces.size() )
00165     {
00166         for( Range::iterator i = edges.begin(); i != edges.end(); ++i )
00167         {
00168             adjents.clear();
00169             error = ahf.get_up_adjacencies( *i, 2, adjents );CHECK_ERR( error );
00170             mbents.clear();
00171             error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error );
00172 
00173             CHECK_EQUAL( mbents.size(), adjents.size() );
00174             std::sort( adjents.begin(), adjents.end() );
00175             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00176             mbents = subtract( mbents, ahfents );
00177             CHECK( !mbents.size() );
00178         }
00179     }
00180 
00181     // NQ2: For every face, obtain neighbor faces
00182     if( faces.size() )
00183     {
00184         for( Range::iterator i = faces.begin(); i != faces.end(); ++i )
00185         {
00186             adjents.clear();
00187             error = ahf.get_neighbor_adjacencies( *i, adjents );CHECK_ERR( error );
00188             mbents.clear();
00189             error = mtu.get_bridge_adjacencies( *i, 1, 2, mbents );CHECK_ERR( error );
00190 
00191             CHECK_EQUAL( mbents.size(), adjents.size() );
00192             std::sort( adjents.begin(), adjents.end() );
00193             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00194             mbents = subtract( mbents, ahfents );
00195             CHECK( !mbents.size() );
00196         }
00197     }
00198 
00199     // DQ 21: For every face, obtain its edges
00200     if( edges.size() && faces.size() )
00201     {
00202         for( Range::iterator i = faces.begin(); i != faces.end(); ++i )
00203         {
00204             adjents.clear();
00205             error = ahf.get_down_adjacencies( *i, 1, adjents );CHECK_ERR( error );
00206             mbents.clear();
00207             error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error );
00208 
00209             CHECK_EQUAL( mbents.size(), adjents.size() );
00210             std::sort( adjents.begin(), adjents.end() );
00211             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00212             mbents = subtract( mbents, ahfents );
00213             CHECK( !mbents.size() );
00214         }
00215     }
00216 
00217     std::cout << "Finished 2D queries: " << std::endl;
00218 
00219     // 3D Queries
00220     // IQ 31: For every vertex, obtain incident cells
00221     if( cells.size() )
00222     {
00223         Range cverts;
00224         error = mbImpl->get_connectivity( cells, cverts, true );CHECK_ERR( error );
00225         for( Range::iterator i = cverts.begin(); i != cverts.end(); ++i )
00226         {
00227             adjents.clear();
00228             error = ahf.get_up_adjacencies( *i, 3, adjents );CHECK_ERR( error );
00229             mbents.clear();
00230             error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error );
00231 
00232             CHECK_EQUAL( mbents.size(), adjents.size() );
00233             std::sort( adjents.begin(), adjents.end() );
00234             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00235             mbents = subtract( mbents, ahfents );
00236             CHECK( !mbents.size() );
00237         }
00238     }
00239 
00240     // IQ 32: For every edge, obtain incident cells
00241     if( edges.size() && cells.size() )
00242     {
00243         for( Range::iterator i = edges.begin(); i != edges.end(); ++i )
00244         {
00245             adjents.clear();
00246             error = ahf.get_up_adjacencies( *i, 3, adjents );CHECK_ERR( error );
00247             mbents.clear();
00248             error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error );
00249 
00250             if( adjents.size() != mbents.size() )
00251             {
00252                 //   std::cout<<"ahf results = "<<std::endl;
00253                 //  ahfents.print();
00254                 //   std::cout<<"native results = "<<std::endl;
00255                 //   mbents.print();
00256             }
00257 
00258             CHECK_EQUAL( mbents.size(), adjents.size() );
00259             std::sort( adjents.begin(), adjents.end() );
00260             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00261             mbents = subtract( mbents, ahfents );
00262             CHECK( !mbents.size() );
00263         }
00264     }
00265 
00266     // IQ33: For every face, obtain incident cells
00267     if( faces.size() && cells.size() )
00268     {
00269         for( Range::iterator i = faces.begin(); i != faces.end(); ++i )
00270         {
00271             adjents.clear();
00272             error = ahf.get_up_adjacencies( *i, 3, adjents );CHECK_ERR( error );
00273             mbents.clear();
00274             error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error );
00275 
00276             if( adjents.size() != mbents.size() )
00277             {
00278                 //   std::cout<<"ahf results = "<<std::endl;
00279                 //  ahfents.print();
00280                 //   std::cout<<"native results = "<<std::endl;
00281                 //   mbents.print();
00282             }
00283 
00284             CHECK_EQUAL( mbents.size(), adjents.size() );
00285             std::sort( adjents.begin(), adjents.end() );
00286             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00287             mbents = subtract( mbents, ahfents );
00288             CHECK( !mbents.size() );
00289         }
00290     }
00291 
00292     // NQ3: For every cell, obtain neighbor cells
00293     if( cells.size() )
00294     {
00295         for( Range::iterator i = cells.begin(); i != cells.end(); ++i )
00296         {
00297             adjents.clear();
00298             error = ahf.get_neighbor_adjacencies( *i, adjents );CHECK_ERR( error );
00299             mbents.clear();
00300             error = mtu.get_bridge_adjacencies( *i, 2, 3, mbents );CHECK_ERR( error );
00301 
00302             CHECK_EQUAL( mbents.size(), adjents.size() );
00303             std::sort( adjents.begin(), adjents.end() );
00304             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00305             mbents = subtract( mbents, ahfents );
00306             CHECK( !mbents.size() );
00307         }
00308     }
00309 
00310     // DQ 31: For every cell, obtain its edges
00311     if( edges.size() && cells.size() )
00312     {
00313         for( Range::iterator i = cells.begin(); i != cells.end(); ++i )
00314         {
00315             adjents.clear();
00316             error = ahf.get_down_adjacencies( *i, 1, adjents );CHECK_ERR( error );
00317             mbents.clear();
00318             error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error );
00319 
00320             CHECK_EQUAL( mbents.size(), adjents.size() );
00321             std::sort( adjents.begin(), adjents.end() );
00322             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00323             mbents = subtract( mbents, ahfents );
00324             CHECK( !mbents.size() );
00325         }
00326     }
00327 
00328     // DQ 32: For every cell, obtain its faces
00329     if( faces.size() && cells.size() )
00330     {
00331         for( Range::iterator i = cells.begin(); i != cells.end(); ++i )
00332         {
00333             adjents.clear();
00334             error = ahf.get_down_adjacencies( *i, 2, adjents );CHECK_ERR( error );
00335             mbents.clear();
00336             error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error );
00337 
00338             CHECK_EQUAL( mbents.size(), adjents.size() );
00339             std::sort( adjents.begin(), adjents.end() );
00340             std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
00341             mbents = subtract( mbents, ahfents );
00342             CHECK( !mbents.size() );
00343         }
00344     }
00345 
00346     std::cout << "Finished 3D queries" << std::endl;
00347 
00348     error = ahf.deinitialize();CHECK_ERR( error );
00349 
00350     return MB_SUCCESS;
00351 }
00352 
00353 int main( int argc, char* argv[] )
00354 {
00355 
00356 #ifdef MOAB_HAVE_MPI
00357     MPI_Init( &argc, &argv );
00358 
00359     int nprocs, rank;
00360     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00361     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00362 #endif
00363 
00364     std::string filename;
00365 #ifdef MOAB_HAVE_HDF5
00366 #ifdef MOAB_HAVE_AHF
00367     filename = TestDir + "/spectral.h5m";
00368 #else
00369     filename = TestDir + "/32hex_ef.h5m";
00370 #endif
00371 #else
00372     filename = TestDir + "/hexes_mixed.vtk";
00373 #endif
00374 
00375     if( argc == 1 )
00376     {
00377 #ifdef MOAB_HAVE_MPI
00378         if( rank == 0 ) std::cout << "Using default input file:" << filename << std::endl;
00379 #else
00380         std::cout << "Using default input file:" << filename << std::endl;
00381 #endif
00382     }
00383 
00384     else if( argc == 2 )
00385         filename = std::string( argv[1] );
00386     else
00387     {
00388         std::cerr << "Usage: " << argv[0] << " [filename]" << std::endl;
00389         return 1;
00390     }
00391 
00392     ErrorCode result;
00393 
00394 #ifdef MOAB_HAVE_MPI
00395     if( rank == 0 ) std::cout << " para_ahf_test: ";
00396 #else
00397     std::cout << "ahf_test:";
00398 #endif
00399 
00400     result = ahf_test( filename.c_str() );
00401     handle_error_code( result, number_tests_failed, number_tests_successful );
00402     std::cout << "\n";
00403 
00404 #ifdef MOAB_HAVE_MPI
00405     MPI_Finalize();
00406 #endif
00407 
00408     return number_tests_failed;
00409 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines