MOAB: Mesh Oriented datABase  (version 5.4.1)
read_mpas_nc.cpp
Go to the documentation of this file.
00001 #include "TestUtil.hpp"
00002 #include "moab/Core.hpp"
00003 #include "moab/ReadUtilIface.hpp"
00004 #include "MBTagConventions.hpp"
00005 
00006 using namespace moab;
00007 
00008 std::string example = TestDir + "unittest/io/mpasx1.642.t.2.nc";
00009 
00010 #ifdef MOAB_HAVE_MPI
00011 #include "moab_mpi.h"
00012 #include "moab/ParallelComm.hpp"
00013 #endif
00014 
00015 void test_read_all();
00016 void test_read_onevar();
00017 void test_read_onetimestep();
00018 void test_read_nomesh();
00019 void test_read_novars();
00020 void test_read_no_mixed_elements();  // Test read option NO_MIXED_ELEMENTS
00021 void test_read_no_edges();           // Test read option NO_EDGES
00022 void test_gather_onevar();           // Test gather set with one variable
00023 
00024 void get_options( std::string& opts );
00025 
00026 const double eps = 1e-20;
00027 
00028 int main( int argc, char* argv[] )
00029 {
00030     int result = 0;
00031 
00032 #ifdef MOAB_HAVE_MPI
00033     int fail = MPI_Init( &argc, &argv );
00034     if( fail ) return 1;
00035 #else
00036     argv[0]   = argv[argc - argc];  // To remove the warnings in serial mode about unused variables
00037 #endif
00038 
00039     result += RUN_TEST( test_read_all );
00040     result += RUN_TEST( test_read_onevar );
00041     result += RUN_TEST( test_read_onetimestep );
00042     result += RUN_TEST( test_read_nomesh );
00043     result += RUN_TEST( test_read_novars );
00044     result += RUN_TEST( test_read_no_mixed_elements );
00045     result += RUN_TEST( test_read_no_edges );
00046     result += RUN_TEST( test_gather_onevar );
00047 
00048 #ifdef MOAB_HAVE_MPI
00049     fail = MPI_Finalize();
00050     if( fail ) return 1;
00051 #endif
00052 
00053     return result;
00054 }
00055 
00056 void test_read_all()
00057 {
00058     Core moab;
00059     Interface& mb = moab;
00060 
00061     std::string opts;
00062     get_options( opts );
00063 
00064     // Read mesh and read all variables at all timesteps
00065     ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
00066 
00067 #ifdef MOAB_HAVE_MPI
00068     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00069     int procs           = pcomm->proc_config().proc_size();
00070 #else
00071     int procs = 1;
00072 #endif
00073 
00074     // Make check runs this test on one processor
00075     if( 1 == procs )
00076     {
00077         // For each tag, check two values
00078         double val[2];
00079 
00080         // Check tags for vertex variable vorticity
00081         Tag vorticity_tag0, vorticity_tag1;
00082         rval = mb.tag_get_handle( "vorticity0", 1, MB_TYPE_DOUBLE, vorticity_tag0 );CHECK_ERR( rval );
00083         rval = mb.tag_get_handle( "vorticity1", 1, MB_TYPE_DOUBLE, vorticity_tag1 );CHECK_ERR( rval );
00084 
00085         // Get vertices (1280 edges)
00086         Range verts;
00087         rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
00088         CHECK_EQUAL( (size_t)1280, verts.size() );
00089         CHECK_EQUAL( (size_t)1, verts.psize() );
00090 
00091         // Check vorticity tag values on first two vertices
00092         EntityHandle vert_ents[] = { verts[0], verts[1] };
00093         rval                     = mb.tag_get_data( vorticity_tag0, vert_ents, 2, val );CHECK_ERR( rval );
00094         CHECK_REAL_EQUAL( 1.1, val[0], eps );
00095         CHECK_REAL_EQUAL( 1.2, val[1], eps );
00096         rval = mb.tag_get_data( vorticity_tag1, vert_ents, 2, val );CHECK_ERR( rval );
00097         CHECK_REAL_EQUAL( 2.1, val[0], eps );
00098         CHECK_REAL_EQUAL( 2.2, val[1], eps );
00099 
00100         // Check tags for edge variable u
00101         Tag u_tag0, u_tag1;
00102         rval = mb.tag_get_handle( "u0", 1, MB_TYPE_DOUBLE, u_tag0 );CHECK_ERR( rval );
00103         rval = mb.tag_get_handle( "u1", 1, MB_TYPE_DOUBLE, u_tag1 );CHECK_ERR( rval );
00104 
00105         // Get edges (1920 edges)
00106         Range edges;
00107         rval = mb.get_entities_by_type( 0, MBEDGE, edges );CHECK_ERR( rval );
00108         CHECK_EQUAL( (size_t)1920, edges.size() );
00109         CHECK_EQUAL( (size_t)1, edges.psize() );
00110 
00111         // Check u tag values on two specified edges
00112         EntityHandle edge_ents[] = { edges[5], edges[6] };
00113         rval                     = mb.tag_get_data( u_tag0, edge_ents, 2, val );CHECK_ERR( rval );
00114         CHECK_REAL_EQUAL( 1.113138721544778, val[0], eps );
00115         CHECK_REAL_EQUAL( -1.113138721930009, val[1], eps );
00116         rval = mb.tag_get_data( u_tag1, edge_ents, 2, val );CHECK_ERR( rval );
00117         CHECK_REAL_EQUAL( 2.113138721544778, val[0], eps );
00118         CHECK_REAL_EQUAL( -2.113138721930009, val[1], eps );
00119 
00120         // Check tags for cell variable ke
00121         Tag ke_tag0, ke_tag1;
00122         rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );CHECK_ERR( rval );
00123         rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );CHECK_ERR( rval );
00124 
00125         // Get cells (12 pentagons and 630 hexagons)
00126         Range cells;
00127         rval = mb.get_entities_by_type( 0, MBPOLYGON, cells );CHECK_ERR( rval );
00128         CHECK_EQUAL( (size_t)642, cells.size() );
00129         CHECK_EQUAL( (size_t)1, cells.psize() );
00130 
00131         // Check ke tag values on first pentagon and first hexagon
00132         EntityHandle cell_ents[] = { cells[0], cells[12] };
00133         rval                     = mb.tag_get_data( ke_tag0, cell_ents, 2, val );CHECK_ERR( rval );
00134         CHECK_REAL_EQUAL( 15.001, val[0], eps );
00135         CHECK_REAL_EQUAL( 16.013, val[1], eps );
00136         rval = mb.tag_get_data( ke_tag1, cell_ents, 2, val );CHECK_ERR( rval );
00137         CHECK_REAL_EQUAL( 25.001, val[0], eps );
00138         CHECK_REAL_EQUAL( 26.013, val[1], eps );
00139     }
00140 }
00141 
00142 void test_read_onevar()
00143 {
00144     Core moab;
00145     Interface& mb = moab;
00146 
00147     std::string opts;
00148     get_options( opts );
00149 
00150     // Read mesh and read cell variable ke at all timesteps
00151     opts += ";VARIABLE=ke";
00152     ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
00153 
00154 #ifdef MOAB_HAVE_MPI
00155     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00156     int procs           = pcomm->proc_config().proc_size();
00157 #else
00158     int procs = 1;
00159 #endif
00160 
00161     // Make check runs this test on one processor
00162     if( 1 == procs )
00163     {
00164         // Check ke tags
00165         Tag ke_tag0, ke_tag1;
00166         rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );CHECK_ERR( rval );
00167         rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );CHECK_ERR( rval );
00168 
00169         // Get cells (12 pentagons and 630 hexagons)
00170         Range cells;
00171         rval = mb.get_entities_by_type( 0, MBPOLYGON, cells );CHECK_ERR( rval );
00172         CHECK_EQUAL( (size_t)642, cells.size() );
00173         CHECK_EQUAL( (size_t)1, cells.psize() );
00174 
00175         // Check ke tag values on 4 cells: first pentagon, last pentagon,
00176         // first hexagon, and last hexagon
00177         EntityHandle cell_ents[] = { cells[0], cells[11], cells[12], cells[641] };
00178         double ke0_val[4];
00179         rval = mb.tag_get_data( ke_tag0, cell_ents, 4, ke0_val );CHECK_ERR( rval );
00180         CHECK_REAL_EQUAL( 15.001, ke0_val[0], eps );
00181         CHECK_REAL_EQUAL( 15.012, ke0_val[1], eps );
00182         CHECK_REAL_EQUAL( 16.013, ke0_val[2], eps );
00183         CHECK_REAL_EQUAL( 16.642, ke0_val[3], eps );
00184         double ke1_val[4];
00185         rval = mb.tag_get_data( ke_tag1, cell_ents, 4, ke1_val );CHECK_ERR( rval );
00186         CHECK_REAL_EQUAL( 25.001, ke1_val[0], eps );
00187         CHECK_REAL_EQUAL( 25.012, ke1_val[1], eps );
00188         CHECK_REAL_EQUAL( 26.013, ke1_val[2], eps );
00189         CHECK_REAL_EQUAL( 26.642, ke1_val[3], eps );
00190     }
00191 }
00192 
00193 void test_read_onetimestep()
00194 {
00195     Core moab;
00196     Interface& mb = moab;
00197 
00198     std::string opts;
00199     get_options( opts );
00200 
00201     // Read mesh and read all variables at 2nd timestep
00202     opts += ";TIMESTEP=1";
00203     ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
00204 
00205     // Check vorticity tags
00206     Tag vorticity_tag0, vorticity_tag1;
00207     rval = mb.tag_get_handle( "vorticity0", 1, MB_TYPE_DOUBLE, vorticity_tag0 );
00208     // Tag vorticity0 should not exist
00209     CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
00210     rval = mb.tag_get_handle( "vorticity1", 1, MB_TYPE_DOUBLE, vorticity_tag1 );CHECK_ERR( rval );
00211 }
00212 
00213 void test_read_nomesh()
00214 {
00215     Core moab;
00216     Interface& mb = moab;
00217 
00218     // Need a file set for nomesh to work right
00219     EntityHandle file_set;
00220     ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
00221 
00222     std::string orig, opts;
00223     get_options( orig );
00224 
00225     // Read mesh and read all variables at 1st timestep
00226     opts = orig + ";TIMESTEP=0";
00227     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00228 
00229     // Check u tags
00230     Tag u_tag0, u_tag1;
00231     rval = mb.tag_get_handle( "u0", 1, MB_TYPE_DOUBLE, u_tag0 );CHECK_ERR( rval );
00232     // Tag u1 should not exist
00233     rval = mb.tag_get_handle( "u1", 1, MB_TYPE_DOUBLE, u_tag1 );
00234     CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
00235 
00236     // Read all variables at 2nd timestep 0, no need to read mesh
00237     opts = orig + ";TIMESTEP=1;NOMESH";
00238     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00239 
00240     // Check tag u1 again
00241     rval = mb.tag_get_handle( "u1", 1, MB_TYPE_DOUBLE, u_tag1 );
00242     // Tag u1 should exist at this time
00243     CHECK_ERR( rval );
00244 }
00245 
00246 void test_read_novars()
00247 {
00248     Core moab;
00249     Interface& mb = moab;
00250 
00251     // Need a file set for nomesh to work right
00252     EntityHandle file_set;
00253     ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
00254 
00255     std::string orig, opts;
00256     get_options( orig );CHECK_ERR( rval );
00257 
00258     // Read header info only, no mesh, no variables
00259     opts = orig + ";NOMESH;VARIABLE=";
00260     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00261 
00262     // Read mesh, but still no variables
00263     opts = orig + ";VARIABLE=;TIMESTEP=0";
00264     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00265 
00266     // Check ke tags
00267     Tag ke_tag0, ke_tag1;
00268     rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );
00269     // Tag ke0 should not exist
00270     CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
00271     rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );
00272     // Tag ke1 should not exist
00273     CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
00274 
00275     // Read ke at 1st timestep, no need to read mesh
00276     opts = orig + ";VARIABLE=ke;TIMESTEP=0;NOMESH";
00277     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00278 
00279     // Check ke tags again
00280     rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );
00281     // Tag ke0 should exist at this time
00282     CHECK_ERR( rval );
00283     // Tag ke1 should still not exist
00284     rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );
00285     CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
00286 
00287     // Read ke at 2nd timestep, no need to read mesh
00288     opts = orig + ";VARIABLE=ke;TIMESTEP=1;NOMESH";
00289     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00290 
00291     // Check tag ke1 again
00292     rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );
00293     // Tag ke1 should exist at this time
00294     CHECK_ERR( rval );
00295 
00296 #ifdef MOAB_HAVE_MPI
00297     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00298     int procs           = pcomm->proc_config().proc_size();
00299 #else
00300     int procs = 1;
00301 #endif
00302 
00303     // Make check runs this test on one processor
00304     if( 1 == procs )
00305     {
00306         // Get cells (12 pentagons and 630 hexagons)
00307         Range cells;
00308         rval = mb.get_entities_by_type( file_set, MBPOLYGON, cells );CHECK_ERR( rval );
00309         CHECK_EQUAL( (size_t)642, cells.size() );
00310         CHECK_EQUAL( (size_t)1, cells.psize() );
00311 
00312         // Check ke tag values on 4 cells: first pentagon, last pentagon,
00313         // first hexagon, and last hexagon
00314         EntityHandle cell_ents[] = { cells[0], cells[11], cells[12], cells[641] };
00315         double ke0_val[4];
00316         rval = mb.tag_get_data( ke_tag0, cell_ents, 4, ke0_val );CHECK_ERR( rval );
00317         CHECK_REAL_EQUAL( 15.001, ke0_val[0], eps );
00318         CHECK_REAL_EQUAL( 15.012, ke0_val[1], eps );
00319         CHECK_REAL_EQUAL( 16.013, ke0_val[2], eps );
00320         CHECK_REAL_EQUAL( 16.642, ke0_val[3], eps );
00321         double ke1_val[4];
00322         rval = mb.tag_get_data( ke_tag1, cell_ents, 4, ke1_val );CHECK_ERR( rval );
00323         CHECK_REAL_EQUAL( 25.001, ke1_val[0], eps );
00324         CHECK_REAL_EQUAL( 25.012, ke1_val[1], eps );
00325         CHECK_REAL_EQUAL( 26.013, ke1_val[2], eps );
00326         CHECK_REAL_EQUAL( 26.642, ke1_val[3], eps );
00327     }
00328 }
00329 
00330 void test_read_no_mixed_elements()
00331 {
00332     Core moab;
00333     Interface& mb = moab;
00334 
00335     std::string opts;
00336     get_options( opts );
00337 
00338     // Read mesh with no mixed elements and read all variables at all timesteps
00339     opts += ";NO_MIXED_ELEMENTS";
00340     ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
00341 
00342 #ifdef MOAB_HAVE_MPI
00343     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00344     int procs           = pcomm->proc_config().proc_size();
00345 #else
00346     int procs = 1;
00347 #endif
00348 
00349     // Make check runs this test on one processor
00350     if( 1 == procs )
00351     {
00352         // Check ke tags
00353         Tag ke_tag0, ke_tag1;
00354         rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0 );CHECK_ERR( rval );
00355         rval = mb.tag_get_handle( "ke1", 1, MB_TYPE_DOUBLE, ke_tag1 );CHECK_ERR( rval );
00356 
00357         // Get cells (12 pentagons and 630 hexagons)
00358         Range cells;
00359         rval = mb.get_entities_by_type( 0, MBPOLYGON, cells );CHECK_ERR( rval );
00360         CHECK_EQUAL( (size_t)642, cells.size() );
00361         // Only one group of cells (pentagons are padded to hexagons,
00362         // e.g. connectivity [1 2 3 4 5] => [1 2 3 4 5 5])
00363         CHECK_EQUAL( (size_t)1, cells.psize() );
00364 
00365         // Check ke tag values on 4 cells: first pentagon, last pentagon,
00366         // first hexagon, and last hexagon
00367         EntityHandle cell_ents[] = { cells[0], cells[11], cells[12], cells[641] };
00368         double ke0_val[4];
00369         rval = mb.tag_get_data( ke_tag0, cell_ents, 4, ke0_val );CHECK_ERR( rval );
00370         CHECK_REAL_EQUAL( 15.001, ke0_val[0], eps );
00371         CHECK_REAL_EQUAL( 15.012, ke0_val[1], eps );
00372         CHECK_REAL_EQUAL( 16.013, ke0_val[2], eps );
00373         CHECK_REAL_EQUAL( 16.642, ke0_val[3], eps );
00374         double ke1_val[4];
00375         rval = mb.tag_get_data( ke_tag1, cell_ents, 4, ke1_val );CHECK_ERR( rval );
00376         CHECK_REAL_EQUAL( 25.001, ke1_val[0], eps );
00377         CHECK_REAL_EQUAL( 25.012, ke1_val[1], eps );
00378         CHECK_REAL_EQUAL( 26.013, ke1_val[2], eps );
00379         CHECK_REAL_EQUAL( 26.642, ke1_val[3], eps );
00380     }
00381 }
00382 
00383 void test_read_no_edges()
00384 {
00385     Core moab;
00386     Interface& mb = moab;
00387 
00388     std::string opts;
00389     get_options( opts );
00390 
00391     opts += ";NO_EDGES;VARIABLE=";
00392     ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
00393 
00394     // Get edges
00395     Range edges;
00396     rval = mb.get_entities_by_type( 0, MBEDGE, edges );CHECK_ERR( rval );
00397     CHECK_EQUAL( (size_t)0, edges.size() );
00398 }
00399 
00400 void test_gather_onevar()
00401 {
00402     Core moab;
00403     Interface& mb = moab;
00404 
00405     EntityHandle file_set;
00406     ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
00407 
00408     std::string opts;
00409     get_options( opts );
00410 
00411     // Read cell variable ke and create gather set on processor 0
00412     opts += ";VARIABLE=ke;GATHER_SET=0";
00413     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00414 
00415 #ifdef MOAB_HAVE_MPI
00416     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00417     int rank            = pcomm->proc_config().proc_rank();
00418 
00419     Range cells, cells_owned;
00420     rval = mb.get_entities_by_type( file_set, MBPOLYGON, cells );CHECK_ERR( rval );
00421 
00422     // Get local owned cells
00423     rval = pcomm->filter_pstatus( cells, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &cells_owned );CHECK_ERR( rval );
00424 
00425     EntityHandle gather_set = 0;
00426     if( 0 == rank )
00427     {
00428         // Get gather set
00429         ReadUtilIface* readUtilIface;
00430         mb.query_interface( readUtilIface );
00431         rval = readUtilIface->get_gather_set( gather_set );CHECK_ERR( rval );
00432         assert( gather_set != 0 );
00433     }
00434 
00435     Tag ke_tag0, gid_tag;
00436     rval = mb.tag_get_handle( "ke0", 1, MB_TYPE_DOUBLE, ke_tag0, MB_TAG_DENSE );CHECK_ERR( rval );
00437 
00438     gid_tag = mb.globalId_tag();
00439 
00440     pcomm->gather_data( cells_owned, ke_tag0, gid_tag, gather_set, 0 );
00441 
00442     if( 0 == rank )
00443     {
00444         // Get gather set cells
00445         Range gather_set_cells;
00446         rval = mb.get_entities_by_type( gather_set, MBPOLYGON, gather_set_cells );CHECK_ERR( rval );
00447         CHECK_EQUAL( (size_t)642, gather_set_cells.size() );
00448         CHECK_EQUAL( (size_t)1, gather_set_cells.psize() );
00449 
00450         // Check ke0 tag values on 4 gather set cells: first pentagon, last pentagon,
00451         // first hexagon, and last hexagon
00452         double ke0_val[4];
00453         EntityHandle cell_ents[] = { gather_set_cells[0], gather_set_cells[11], gather_set_cells[12],
00454                                      gather_set_cells[641] };
00455         rval                     = mb.tag_get_data( ke_tag0, cell_ents, 4, ke0_val );CHECK_ERR( rval );
00456         CHECK_REAL_EQUAL( 15.001, ke0_val[0], eps );
00457         CHECK_REAL_EQUAL( 15.012, ke0_val[1], eps );
00458         CHECK_REAL_EQUAL( 16.013, ke0_val[2], eps );
00459         CHECK_REAL_EQUAL( 16.642, ke0_val[3], eps );
00460     }
00461 #endif
00462 }
00463 
00464 void get_options( std::string& opts )
00465 {
00466 #ifdef MOAB_HAVE_MPI
00467     // Use parallel options
00468     opts = std::string( ";;PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL" );
00469 #else
00470     opts      = std::string( ";;" );
00471 #endif
00472 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines