MOAB: Mesh Oriented datABase  (version 5.5.0)
read_ucd_nc.cpp
Go to the documentation of this file.
1 #include "TestUtil.hpp"
2 #include "moab/Core.hpp"
3 #include "moab/ReadUtilIface.hpp"
4 #include "TagInfo.hpp"
5 #include "MBTagConventions.hpp"
6 
7 using namespace moab;
8 
9 std::string example = TestDir + "unittest/io/homme3x3458.t.3.nc";
10 std::string conn_fname = TestDir + "unittest/io/HommeMapping.nc";
11 
12 #ifdef MOAB_HAVE_MPI
13 #include "moab_mpi.h"
14 #include "moab/ParallelComm.hpp"
15 #endif
16 
17 void test_read_all();
18 void test_read_onevar();
19 void test_read_onetimestep();
20 void test_read_nomesh();
21 void test_read_novars();
22 void test_read_coord_vars(); // Test reading coordinate variables
23 void test_gather_onevar(); // Test gather set with one variable
24 void test_read_conn(); // Test reading connectivity file only
25 
26 void get_options( std::string& opts );
27 
28 const int levels = 3;
29 
30 int main( int argc, char* argv[] )
31 {
32  int result = 0;
33 
34 #ifdef MOAB_HAVE_MPI
35  int fail = MPI_Init( &argc, &argv );
36  if( fail ) return 1;
37 #else
38  argv[0] = argv[argc - argc]; // To remove the warnings in serial mode about unused variables
39 #endif
40 
41  result += RUN_TEST( test_read_all );
42  result += RUN_TEST( test_read_onevar );
43  result += RUN_TEST( test_read_onetimestep );
44  result += RUN_TEST( test_read_nomesh );
45  result += RUN_TEST( test_read_novars );
46  result += RUN_TEST( test_read_coord_vars );
47  result += RUN_TEST( test_gather_onevar );
48  result += RUN_TEST( test_read_conn );
49 
50 #ifdef MOAB_HAVE_MPI
51  fail = MPI_Finalize();
52  if( fail ) return 1;
53 #endif
54 
55  return result;
56 }
57 
58 void test_read_all()
59 {
60  Core moab;
61  Interface& mb = moab;
62 
63  std::string opts;
64  get_options( opts );
65 
66  ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
67 
68  // Check for proper tags
69  Tag Ttag0, Ttag1;
70  rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );CHECK_ERR( rval );
71 
72  rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
73 }
74 
75 void test_read_onevar()
76 {
77  Core moab;
78  Interface& mb = moab;
79 
80  std::string opts;
81  get_options( opts );
82 
83  // Read mesh and read vertex variable T at all timesteps
84  opts += std::string( ";VARIABLE=T" );
85  ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
86 
87 #ifdef MOAB_HAVE_MPI
88  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
89  int procs = pcomm->proc_config().proc_size();
90 #else
91  int procs = 1;
92 #endif
93 
94  // Make check runs this test on one processor
95  if( 1 == procs )
96  {
97  // Check for proper tags
98  Tag Ttag0, Ttag1;
99  rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );CHECK_ERR( rval );
100  rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
101 
102  // Get vertices
103  Range verts;
104  rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
105  CHECK_EQUAL( (size_t)3458, verts.size() );
106 
107  // Get all values of tag T0
108  int count;
109  void* Tbuf;
110  rval = mb.tag_iterate( Ttag0, verts.begin(), verts.end(), count, Tbuf );CHECK_ERR( rval );
111  CHECK_EQUAL( (size_t)count, verts.size() );
112 
113  const double eps = 0.0001;
114  double* data = (double*)Tbuf;
115 
116  // Check first level values on 4 strategically selected vertices
117  CHECK_REAL_EQUAL( 233.1136, data[0 * levels], eps ); // First vert
118  CHECK_REAL_EQUAL( 236.1505, data[1728 * levels], eps ); // Median vert
119  CHECK_REAL_EQUAL( 235.7722, data[1729 * levels], eps ); // Median vert
120  CHECK_REAL_EQUAL( 234.0416, data[3457 * levels], eps ); // Last vert
121  }
122 }
123 
124 void test_read_onetimestep()
125 {
126  Core moab;
127  Interface& mb = moab;
128 
129  std::string opts;
130  get_options( opts );
131 
132  opts += std::string( ";TIMESTEP=1" );
133  ErrorCode rval = mb.load_file( example.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
134 
135  // Check for proper tags
136  Tag Ttag0, Ttag1;
137  rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );
138  CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
139 
140  rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
141 }
142 
143 void test_read_nomesh()
144 {
145  Core moab;
146  Interface& mb = moab;
147 
148  // Need a set for nomesh to work right
149  EntityHandle file_set;
150  ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
151 
152  std::string orig, opts;
153  get_options( orig );
154 
155  opts = orig + std::string( ";TIMESTEP=0" );
156  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
157 
158  // Check for proper tag
159  Tag Ttag0, Ttag1;
160  rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );CHECK_ERR( rval );
161 
162  rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );
163  CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
164 
165  // Now read 2nd timestep with nomesh option
166  opts = orig + std::string( ";TIMESTEP=1;NOMESH" );
167  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
168 
169  // Check for proper tag
170  rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
171 }
172 
173 void test_read_novars()
174 {
175  Core moab;
176  Interface& mb = moab;
177 
178  // Need a set for nomesh to work right
179  EntityHandle set;
180  ErrorCode rval = mb.create_meshset( MESHSET_SET, set );CHECK_ERR( rval );
181 
182  std::string orig, opts;
183  get_options( orig );
184 
185  opts = orig + std::string( ";NOMESH;VARIABLE=" );
186  rval = mb.load_file( example.c_str(), &set, opts.c_str() );CHECK_ERR( rval );
187 
188  opts = orig + std::string( ";VARIABLE=;TIMESTEP=0" );
189  rval = mb.load_file( example.c_str(), &set, opts.c_str() );CHECK_ERR( rval );
190 
191  // Check for proper tag
192  Tag Ttag0, Ttag1;
193  rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );
194  CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
195 
196  opts = orig + std::string( ";VARIABLE=T;TIMESTEP=0;NOMESH" );
197  rval = mb.load_file( example.c_str(), &set, opts.c_str() );CHECK_ERR( rval );
198 
199  rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0 );CHECK_ERR( rval );
200 
201  rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );
202  CHECK_EQUAL( rval, MB_TAG_NOT_FOUND );
203 
204  // Now read 2nd timestep with nomesh option
205  opts = orig + std::string( ";TIMESTEP=1;NOMESH" );
206  rval = mb.load_file( example.c_str(), &set, opts.c_str() );CHECK_ERR( rval );
207 
208  // Check for proper tag
209  rval = mb.tag_get_handle( "T1", levels, MB_TYPE_DOUBLE, Ttag1 );CHECK_ERR( rval );
210 }
211 
212 void test_read_coord_vars()
213 {
214  Core moab;
215  Interface& mb = moab;
216 
217  EntityHandle file_set;
218  ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
219 
220  std::string orig, opts;
221  get_options( orig );
222 
223  opts = orig + std::string( ";NOMESH;VARIABLE=" );
224  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
225 
226  std::string tag_name;
227  int var_len;
228  Tag var_tag;
229  const void* var_data;
230 
231  // Check tag for regular coordinate variable lev
232  tag_name = "lev";
233  var_len = 0;
234  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 );
235  CHECK_EQUAL( true, var_tag->variable_length() );
236  CHECK_EQUAL( MB_TYPE_DOUBLE, var_tag->get_data_type() );
237 
238  // Check lev tag size and values on file_set
239  rval = mb.tag_get_by_ptr( var_tag, &file_set, 1, &var_data, &var_len );CHECK_ERR( rval );
240  CHECK_EQUAL( levels, var_len );
241  double* lev_val = (double*)var_data;
242  const double eps = 1e-10;
243  CHECK_REAL_EQUAL( 3.54463800000002, lev_val[0], eps );
244  CHECK_REAL_EQUAL( 13.9672100000001, lev_val[levels - 1], eps );
245 
246  // Check tag for dummy coordinate variable ncol
247  tag_name = "ncol";
248  var_len = 0;
249  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 );
250  CHECK_EQUAL( true, var_tag->variable_length() );
251  CHECK_EQUAL( MB_TYPE_INTEGER, var_tag->get_data_type() );
252 
253  // Check ncol tag size and values on file_set
254  rval = mb.tag_get_by_ptr( var_tag, &file_set, 1, &var_data, &var_len );CHECK_ERR( rval );
255  CHECK_EQUAL( 1, var_len );
256  int* ncol_val = (int*)var_data;
257  CHECK_EQUAL( 3458, ncol_val[0] );
258 
259  // Create another file set
260  EntityHandle file_set2;
261  rval = mb.create_meshset( MESHSET_SET, file_set2 );CHECK_ERR( rval );
262 
263  // Read file again with file_set2
264  rval = mb.load_file( example.c_str(), &file_set2, opts.c_str() );CHECK_ERR( rval );
265 
266  // Check tag for regular coordinate lev
267  tag_name = "lev";
268  var_len = 0;
269  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 );
270  CHECK_EQUAL( true, var_tag->variable_length() );
271  CHECK_EQUAL( MB_TYPE_DOUBLE, var_tag->get_data_type() );
272 
273  // Check lev tag size and values on file_set2
274  rval = mb.tag_get_by_ptr( var_tag, &file_set2, 1, &var_data, &var_len );CHECK_ERR( rval );
275  CHECK_EQUAL( levels, var_len );
276  lev_val = (double*)var_data;
277  CHECK_REAL_EQUAL( 3.54463800000002, lev_val[0], eps );
278  CHECK_REAL_EQUAL( 13.9672100000001, lev_val[levels - 1], eps );
279 
280  // Check tag for dummy coordinate variable ncol
281  tag_name = "ncol";
282  var_len = 0;
283  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 );
284  CHECK_EQUAL( true, var_tag->variable_length() );
285  CHECK_EQUAL( MB_TYPE_INTEGER, var_tag->get_data_type() );
286 
287  // Check ncol tag size and values on file_set2
288  rval = mb.tag_get_by_ptr( var_tag, &file_set2, 1, &var_data, &var_len );CHECK_ERR( rval );
289  CHECK_EQUAL( 1, var_len );
290  ncol_val = (int*)var_data;
291  CHECK_EQUAL( 3458, ncol_val[0] );
292 }
293 
294 void test_gather_onevar()
295 {
296  Core moab;
297  Interface& mb = moab;
298 
299  EntityHandle file_set;
300  ErrorCode rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
301 
302  std::string opts;
303  get_options( opts );
304 
305  // Read vertex variable T and create gather set on processor 0
306  opts += ";VARIABLE=T;GATHER_SET=0";
307 #ifdef MOAB_HAVE_MPI
308  opts += ";PARALLEL_RESOLVE_SHARED_ENTS";
309 #endif
310  rval = mb.load_file( example.c_str(), &file_set, opts.c_str() );CHECK_ERR( rval );
311 
312 #ifdef MOAB_HAVE_MPI
313  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
314  int rank = pcomm->proc_config().proc_rank();
315 
316  Range verts, verts_owned;
317  rval = mb.get_entities_by_type( file_set, MBVERTEX, verts );CHECK_ERR( rval );
318 
319  // Get local owned vertices
320  rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &verts_owned );CHECK_ERR( rval );
321 
322  EntityHandle gather_set = 0;
323  if( 0 == rank )
324  {
325  // Get gather set
326  ReadUtilIface* readUtilIface;
327  mb.query_interface( readUtilIface );
328  rval = readUtilIface->get_gather_set( gather_set );CHECK_ERR( rval );
329  assert( gather_set != 0 );
330  }
331 
332  Tag Ttag0, gid_tag;
333  rval = mb.tag_get_handle( "T0", levels, MB_TYPE_DOUBLE, Ttag0, MB_TAG_DENSE );CHECK_ERR( rval );
334 
335  gid_tag = mb.globalId_tag();
336 
337  pcomm->gather_data( verts_owned, Ttag0, gid_tag, gather_set, 0 );
338 
339  if( 0 == rank )
340  {
341  // Get gather set vertices
342  Range gather_set_verts;
343  rval = mb.get_entities_by_type( gather_set, MBVERTEX, gather_set_verts );CHECK_ERR( rval );
344  CHECK_EQUAL( (size_t)3458, gather_set_verts.size() );
345 
346  // Get T0 tag values on 4 strategically selected gather set vertices
347  double T0_val[4 * levels];
348  EntityHandle vert_ents[] = { gather_set_verts[0], gather_set_verts[1728], gather_set_verts[1729],
349  gather_set_verts[3457] };
350  rval = mb.tag_get_data( Ttag0, vert_ents, 4, T0_val );CHECK_ERR( rval );
351 
352  const double eps = 0.001;
353 
354  // Check first level values
355  CHECK_REAL_EQUAL( 233.1136, T0_val[0 * levels], eps ); // First vert
356  CHECK_REAL_EQUAL( 236.1505, T0_val[1 * levels], eps ); // Median vert
357  CHECK_REAL_EQUAL( 235.7722, T0_val[2 * levels], eps ); // Median vert
358  CHECK_REAL_EQUAL( 234.0416, T0_val[3 * levels], eps ); // Last vert
359  }
360 #endif
361 }
362 
363 void test_read_conn()
364 {
365  Core moab;
366  Interface& mb = moab;
367 
368  std::string opts;
369  get_options( opts );
370 
371  ErrorCode rval = mb.load_file( conn_fname.c_str(), NULL, opts.c_str() );CHECK_ERR( rval );
372 
373 #ifdef MOAB_HAVE_MPI
374  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
375  int procs = pcomm->proc_config().proc_size();
376 #else
377  int procs = 1;
378 #endif
379 
380  // Make check runs this test on one processor
381  if( 1 == procs )
382  {
383  // Get vertices
384  Range verts;
385  rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
386  CHECK_EQUAL( (size_t)3458, verts.size() );
387 
388  // Get cells
389  Range cells;
390  rval = mb.get_entities_by_type( 0, MBQUAD, cells );CHECK_ERR( rval );
391  CHECK_EQUAL( (size_t)3456, cells.size() );
392  }
393 }
394 
395 void get_options( std::string& opts )
396 {
397 #ifdef MOAB_HAVE_MPI
398  // Use parallel options
399  opts = std::string( ";;PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL" );
400 #else
401  opts = std::string( ";;" );
402 #endif
403 }