MOAB: Mesh Oriented datABase  (version 5.2.1)
read_ucd_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 "TagInfo.hpp"
00005 #include "MBTagConventions.hpp"
00006 
00007 using namespace moab;
00008 
00009 std::string example    = TestDir + "/io/homme3x3458.t.3.nc";
00010 std::string conn_fname = TestDir + "/io/HommeMapping.nc";
00011 
00012 #ifdef MOAB_HAVE_MPI
00013 #include "moab_mpi.h"
00014 #include "moab/ParallelComm.hpp"
00015 #endif
00016 
00017 void test_read_all();
00018 void test_read_onevar();
00019 void test_read_onetimestep();
00020 void test_read_nomesh();
00021 void test_read_novars();
00022 void test_read_coord_vars();  // Test reading coordinate variables
00023 void test_gather_onevar();    // Test gather set with one variable
00024 void test_read_conn();        // Test reading connectivity file only
00025 
00026 void get_options( std::string& opts );
00027 
00028 const int levels = 3;
00029 
00030 int main( int argc, char* argv[] )
00031 {
00032     int result = 0;
00033 
00034 #ifdef MOAB_HAVE_MPI
00035     int fail = MPI_Init( &argc, &argv );
00036     if( fail ) return 1;
00037 #else
00038     argv[0]   = argv[argc - argc];  // To remove the warnings in serial mode about unused variables
00039 #endif
00040 
00041     result += RUN_TEST( test_read_all );
00042     result += RUN_TEST( test_read_onevar );
00043     result += RUN_TEST( test_read_onetimestep );
00044     result += RUN_TEST( test_read_nomesh );
00045     result += RUN_TEST( test_read_novars );
00046     result += RUN_TEST( test_read_coord_vars );
00047     result += RUN_TEST( test_gather_onevar );
00048     result += RUN_TEST( test_read_conn );
00049 
00050 #ifdef MOAB_HAVE_MPI
00051     fail = MPI_Finalize();
00052     if( fail ) return 1;
00053 #endif
00054 
00055     return result;
00056 }
00057 
00058 void test_read_all()
00059 {
00060     Core moab;
00061     Interface& mb = moab;
00062 
00063     std::string opts;
00064     get_options( opts );
00065 
00066     ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
00067 
00068     // Check for proper tags
00069     Tag Ttag0, Ttag1;
00070     rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );CHECK_ERR( rval );
00071 
00072     rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
00073 }
00074 
00075 void test_read_onevar()
00076 {
00077     Core moab;
00078     Interface& mb = moab;
00079 
00080     std::string opts;
00081     get_options( opts );
00082 
00083     // Read mesh and read vertex variable T at all timesteps
00084     opts += std::string( ";VARIABLE=T" );
00085     ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
00086 
00087 #ifdef MOAB_HAVE_MPI
00088     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00089     int procs           = pcomm->proc_config().proc_size();
00090 #else
00091     int procs = 1;
00092 #endif
00093 
00094     // Make check runs this test on one processor
00095     if( 1 == procs )
00096     {
00097         // Check for proper tags
00098         Tag Ttag0, Ttag1;
00099         rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );CHECK_ERR( rval );
00100         rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
00101 
00102         // Get vertices
00103         Range verts;
00104         rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
00105         CHECK_EQUAL( (size_t)3458, verts.size() );
00106 
00107         // Get all values of tag T0
00108         int count;
00109         void* Tbuf;
00110         rval = mb.tag_iterate( Ttag0, verts.begin(), verts.end(), count, Tbuf );CHECK_ERR( rval );
00111         CHECK_EQUAL( (size_t)count, verts.size() );
00112 
00113         const double eps = 0.0001;
00114         double* data     = (double*)Tbuf;
00115 
00116         // Check first level values on 4 strategically selected vertices
00117         CHECK_REAL_EQUAL( 233.1136, data[0 * levels], eps );     // First vert
00118         CHECK_REAL_EQUAL( 236.1505, data[1728 * levels], eps );  // Median vert
00119         CHECK_REAL_EQUAL( 235.7722, data[1729 * levels], eps );  // Median vert
00120         CHECK_REAL_EQUAL( 234.0416, data[3457 * levels], eps );  // Last vert
00121     }
00122 }
00123 
00124 void test_read_onetimestep()
00125 {
00126     Core moab;
00127     Interface& mb = moab;
00128 
00129     std::string opts;
00130     get_options( opts );
00131 
00132     opts += std::string( ";TIMESTEP=1" );
00133     ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
00134 
00135     // Check for proper tags
00136     Tag Ttag0, Ttag1;
00137     rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );
00138     CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
00139 
00140     rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
00141 }
00142 
00143 void test_read_nomesh()
00144 {
00145     Core moab;
00146     Interface& mb = moab;
00147 
00148     // Need a set for nomesh to work right
00149     EntityHandle file_set;
00150     ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
00151 
00152     std::string orig, opts;
00153     get_options( orig );
00154 
00155     opts = orig + std::string( ";TIMESTEP=0" );
00156     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00157 
00158     // Check for proper tag
00159     Tag Ttag0, Ttag1;
00160     rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );CHECK_ERR( rval );
00161 
00162     rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );
00163     CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
00164 
00165     // Now read 2nd timestep with nomesh option
00166     opts = orig + std::string( ";TIMESTEP=1;NOMESH" );
00167     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00168 
00169     // Check for proper tag
00170     rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
00171 }
00172 
00173 void test_read_novars()
00174 {
00175     Core moab;
00176     Interface& mb = moab;
00177 
00178     // Need a set for nomesh to work right
00179     EntityHandle set;
00180     ErrorCode rval = mb.create_meshset( MESHSET_SET, set );CHECK_ERR( rval );
00181 
00182     std::string orig, opts;
00183     get_options( orig );
00184 
00185     opts = orig + std::string( ";NOMESH;VARIABLE=" );
00186     rval = mb.load_file( example.c_str(), &set, opts.c_str() );CHECK_ERR( rval );
00187 
00188     opts = orig + std::string( ";VARIABLE=;TIMESTEP=0" );
00189     rval = mb.load_file( example.c_str(), &set, opts.c_str() );CHECK_ERR( rval );
00190 
00191     // Check for proper tag
00192     Tag Ttag0, Ttag1;
00193     rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );
00194     CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
00195 
00196     opts = orig + std::string( ";VARIABLE=T;TIMESTEP=0;NOMESH" );
00197     rval = mb.load_file( example.c_str(), &set, opts.c_str() );CHECK_ERR( rval );
00198 
00199     rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );CHECK_ERR( rval );
00200 
00201     rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );
00202     CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
00203 
00204     // Now read 2nd timestep with nomesh option
00205     opts = orig + std::string( ";TIMESTEP=1;NOMESH" );
00206     rval = mb.load_file( example.c_str(), &set, opts.c_str() );CHECK_ERR( rval );
00207 
00208     // Check for proper tag
00209     rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
00210 }
00211 
00212 void test_read_coord_vars()
00213 {
00214     Core moab;
00215     Interface& mb = moab;
00216 
00217     EntityHandle file_set;
00218     ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
00219 
00220     std::string orig, opts;
00221     get_options( orig );
00222 
00223     opts = orig + std::string( ";NOMESH;VARIABLE=" );
00224     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00225 
00226     std::string tag_name;
00227     int var_len;
00228     Tag var_tag;
00229     const void* var_data;
00230 
00231     // Check tag for regular coordinate variable lev
00232     tag_name = "lev";
00233     var_len  = 0;
00234     rval     = mb.tag_get_handle( tag_name.c_str(), var_len, MB_TYPE_OPAQUE, var_tag, MB_TAG_SPARSE | MB_TAG_VARLEN );CHECK_ERR( rval );
00235     CHECK_EQUAL( true, var_tag->variable_length() );
00236     CHECK_EQUAL( MB_TYPE_DOUBLE, var_tag->get_data_type() );
00237 
00238     // Check lev tag size and values on file_set
00239     rval = mb.tag_get_by_ptr( var_tag, &file_set, 1, &var_data, &var_len );CHECK_ERR( rval );
00240     CHECK_EQUAL( levels, var_len );
00241     double* lev_val  = (double*)var_data;
00242     const double eps = 1e-10;
00243     CHECK_REAL_EQUAL( 3.54463800000002, lev_val[0], eps );
00244     CHECK_REAL_EQUAL( 13.9672100000001, lev_val[levels - 1], eps );
00245 
00246     // Check tag for dummy coordinate variable ncol
00247     tag_name = "ncol";
00248     var_len  = 0;
00249     rval     = mb.tag_get_handle( tag_name.c_str(), var_len, MB_TYPE_OPAQUE, var_tag, MB_TAG_SPARSE | MB_TAG_VARLEN );CHECK_ERR( rval );
00250     CHECK_EQUAL( true, var_tag->variable_length() );
00251     CHECK_EQUAL( MB_TYPE_INTEGER, var_tag->get_data_type() );
00252 
00253     // Check ncol tag size and values on file_set
00254     rval = mb.tag_get_by_ptr( var_tag, &file_set, 1, &var_data, &var_len );CHECK_ERR( rval );
00255     CHECK_EQUAL( 1, var_len );
00256     int* ncol_val = (int*)var_data;
00257     CHECK_EQUAL( 3458, ncol_val[0] );
00258 
00259     // Create another file set
00260     EntityHandle file_set2;
00261     rval = mb.create_meshset( MESHSET_SET, file_set2 );CHECK_ERR( rval );
00262 
00263     // Read file again with file_set2
00264     rval = mb.load_file( example.c_str(), &file_set2, opts.c_str() );CHECK_ERR( rval );
00265 
00266     // Check tag for regular coordinate lev
00267     tag_name = "lev";
00268     var_len  = 0;
00269     rval     = mb.tag_get_handle( tag_name.c_str(), var_len, MB_TYPE_OPAQUE, var_tag, MB_TAG_SPARSE | MB_TAG_VARLEN );CHECK_ERR( rval );
00270     CHECK_EQUAL( true, var_tag->variable_length() );
00271     CHECK_EQUAL( MB_TYPE_DOUBLE, var_tag->get_data_type() );
00272 
00273     // Check lev tag size and values on file_set2
00274     rval = mb.tag_get_by_ptr( var_tag, &file_set2, 1, &var_data, &var_len );CHECK_ERR( rval );
00275     CHECK_EQUAL( levels, var_len );
00276     lev_val = (double*)var_data;
00277     CHECK_REAL_EQUAL( 3.54463800000002, lev_val[0], eps );
00278     CHECK_REAL_EQUAL( 13.9672100000001, lev_val[levels - 1], eps );
00279 
00280     // Check tag for dummy coordinate variable ncol
00281     tag_name = "ncol";
00282     var_len  = 0;
00283     rval     = mb.tag_get_handle( tag_name.c_str(), var_len, MB_TYPE_OPAQUE, var_tag, MB_TAG_SPARSE | MB_TAG_VARLEN );CHECK_ERR( rval );
00284     CHECK_EQUAL( true, var_tag->variable_length() );
00285     CHECK_EQUAL( MB_TYPE_INTEGER, var_tag->get_data_type() );
00286 
00287     // Check ncol tag size and values on file_set2
00288     rval = mb.tag_get_by_ptr( var_tag, &file_set2, 1, &var_data, &var_len );CHECK_ERR( rval );
00289     CHECK_EQUAL( 1, var_len );
00290     ncol_val = (int*)var_data;
00291     CHECK_EQUAL( 3458, ncol_val[0] );
00292 }
00293 
00294 void test_gather_onevar()
00295 {
00296     Core moab;
00297     Interface& mb = moab;
00298 
00299     EntityHandle file_set;
00300     ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
00301 
00302     std::string opts;
00303     get_options( opts );
00304 
00305     // Read vertex variable T and create gather set on processor 0
00306     opts += ";VARIABLE=T;GATHER_SET=0";
00307 #ifdef MOAB_HAVE_MPI
00308     opts += ";PARALLEL_RESOLVE_SHARED_ENTS";
00309 #endif
00310     rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
00311 
00312 #ifdef MOAB_HAVE_MPI
00313     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00314     int rank            = pcomm->proc_config().proc_rank();
00315 
00316     Range verts, verts_owned;
00317     rval = mb.get_entities_by_type( file_set, MBVERTEX, verts );CHECK_ERR( rval );
00318 
00319     // Get local owned vertices
00320     rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &verts_owned );CHECK_ERR( rval );
00321 
00322     EntityHandle gather_set = 0;
00323     if( 0 == rank )
00324     {
00325         // Get gather set
00326         ReadUtilIface* readUtilIface;
00327         mb.query_interface( readUtilIface );
00328         rval = readUtilIface->get_gather_set( gather_set );CHECK_ERR( rval );
00329         assert( gather_set != 0 );
00330     }
00331 
00332     Tag Ttag0, gid_tag;
00333     rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0, MB_TAG_DENSE );CHECK_ERR( rval );
00334 
00335     gid_tag = mb.globalId_tag();
00336 
00337     pcomm->gather_data( verts_owned, Ttag0, gid_tag, gather_set, 0 );
00338 
00339     if( 0 == rank )
00340     {
00341         // Get gather set vertices
00342         Range gather_set_verts;
00343         rval = mb.get_entities_by_type( gather_set, MBVERTEX, gather_set_verts );CHECK_ERR( rval );
00344         CHECK_EQUAL( (size_t)3458, gather_set_verts.size() );
00345 
00346         // Get T0 tag values on 4 strategically selected gather set vertices
00347         double T0_val[4 * levels];
00348         EntityHandle vert_ents[] = { gather_set_verts[0], gather_set_verts[1728], gather_set_verts[1729],
00349                                      gather_set_verts[3457] };
00350         rval                     = mb.tag_get_data( Ttag0, vert_ents, 4, T0_val );CHECK_ERR( rval );
00351 
00352         const double eps = 0.001;
00353 
00354         // Check first level values
00355         CHECK_REAL_EQUAL( 233.1136, T0_val[0 * levels], eps );  // First vert
00356         CHECK_REAL_EQUAL( 236.1505, T0_val[1 * levels], eps );  // Median vert
00357         CHECK_REAL_EQUAL( 235.7722, T0_val[2 * levels], eps );  // Median vert
00358         CHECK_REAL_EQUAL( 234.0416, T0_val[3 * levels], eps );  // Last vert
00359     }
00360 #endif
00361 }
00362 
00363 void test_read_conn()
00364 {
00365     Core moab;
00366     Interface& mb = moab;
00367 
00368     std::string opts;
00369     get_options( opts );
00370 
00371     ErrorCode rval = mb.load_file( conn_fname.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
00372 
00373 #ifdef MOAB_HAVE_MPI
00374     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00375     int procs           = pcomm->proc_config().proc_size();
00376 #else
00377     int procs = 1;
00378 #endif
00379 
00380     // Make check runs this test on one processor
00381     if( 1 == procs )
00382     {
00383         // Get vertices
00384         Range verts;
00385         rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
00386         CHECK_EQUAL( (size_t)3458, verts.size() );
00387 
00388         // Get cells
00389         Range cells;
00390         rval = mb.get_entities_by_type( 0, MBQUAD, cells );CHECK_ERR( rval );
00391         CHECK_EQUAL( (size_t)3456, cells.size() );
00392     }
00393 }
00394 
00395 void get_options( std::string& opts )
00396 {
00397 #ifdef MOAB_HAVE_MPI
00398     // Use parallel options
00399     opts = std::string( ";;PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL" );
00400 #else
00401     opts      = std::string( ";;" );
00402 #endif
00403 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines