|
MOAB: Mesh Oriented datABase
(version 5.2.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* 00002 * addncdata.cpp 00003 * this tool will take an existing h5m file and add data from an nc type file 00004 * will support mainly showing the data associated with unstructured meshes (Homme, MPAS) with Visit 00005 * 00006 * example of usage: 00007 * ./mbaddnc -i wholeFineATM.h5m -n surfdata_ne11np4_simyr1850_c160614.nc -o 00008 * whole_LONGXY_surfdata.h5m -v LONGXY 00009 * 00010 * Basically, will output a new h5m file (whole_LONGXY_surfdata.h5m), which has an extra tag, 00011 * corresponding to the variable LONGXY from the file surfdata_ne11np4_simyr1850_c160614.nc; 00012 * matching is based on the global ids between what we think is the order on the original file 00013 * (wholeFineATM.h5m) and the order of surfdata_ne11np4_simyr1850_c160614.nc 00014 * 00015 * file wholeFineATM.h5m is obtained from a coupled run in e3sm, with the ne 11, np 4, 00016 * 00017 * add an option to output the nc data to the original coarse ATM file, the one that also has the 00018 * GLOBAL_DOFS tag with the global DOFs of the GLL points 00019 */ 00020 00021 #include "moab/ProgOptions.hpp" 00022 #include "moab/Core.hpp" 00023 00024 #include "netcdf.h" 00025 #include <math.h> 00026 #include <sstream> 00027 00028 using namespace moab; 00029 00030 // copy from ReadNCDF.cpp some useful macros for reading from a netcdf file (exodus?) 00031 // ncFile is an integer initialized when opening the nc file in read mode 00032 00033 int ncFile; 00034 00035 #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id ) 00036 00037 #define GET_DIM( ncdim, name, val ) \ 00038 { \ 00039 int gdfail = nc_inq_dimid( ncFile, name, &ncdim ); \ 00040 if( NC_NOERR == gdfail ) \ 00041 { \ 00042 size_t tmp_val; \ 00043 gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val ); \ 00044 if( NC_NOERR != gdfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); } \ 00045 else \ 00046 val = tmp_val; \ 00047 } \ 00048 else \ 00049 val = 0; \ 00050 } 00051 00052 #define GET_DIMB( ncdim, name, varname, id, val ) \ 00053 INS_ID( name, varname, id ); \ 00054 GET_DIM( ncdim, name, val ); 00055 00056 #define GET_VAR( name, id, dims ) \ 00057 { \ 00058 id = -1; \ 00059 int gvfail = nc_inq_varid( ncFile, name, &id ); \ 00060 if( NC_NOERR == gvfail ) \ 00061 { \ 00062 int ndims; \ 00063 gvfail = nc_inq_varndims( ncFile, id, &ndims ); \ 00064 if( NC_NOERR == gvfail ) \ 00065 { \ 00066 dims.resize( ndims ); \ 00067 gvfail = nc_inq_vardimid( ncFile, id, &dims[0] ); \ 00068 if( NC_NOERR != gvfail ) \ 00069 { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get variable dimension IDs" ); } \ 00070 } \ 00071 } \ 00072 } 00073 00074 #define GET_1D_INT_VAR( name, id, vals ) \ 00075 { \ 00076 GET_VAR( name, id, vals ); \ 00077 if( -1 != id ) \ 00078 { \ 00079 size_t ntmp; \ 00080 int ivfail = nc_inq_dimlen( ncFile, vals[0], &ntmp ); \ 00081 if( NC_NOERR != ivfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); } \ 00082 vals.resize( ntmp ); \ 00083 size_t ntmp1 = 0; \ 00084 ivfail = nc_get_vara_int( ncFile, id, &ntmp1, &ntmp, &vals[0] ); \ 00085 if( NC_NOERR != ivfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << name ); } \ 00086 } \ 00087 } 00088 00089 #define GET_1D_DBL_VAR( name, id, vals ) \ 00090 { \ 00091 std::vector< int > dum_dims; \ 00092 GET_VAR( name, id, dum_dims ); \ 00093 if( -1 != id ) \ 00094 { \ 00095 size_t ntmp; \ 00096 int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp ); \ 00097 if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); } \ 00098 vals.resize( ntmp ); \ 00099 size_t ntmp1 = 0; \ 00100 dvfail = nc_get_vara_double( ncFile, id, &ntmp1, &ntmp, &vals[0] ); \ 00101 if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << name ); } \ 00102 } \ 00103 } 00104 00105 #define GET_1D_FLT_VAR( name, id, vals ) \ 00106 { \ 00107 std::vector< int > dum_dims; \ 00108 GET_VAR( name, id, dum_dims ); \ 00109 if( -1 != id ) \ 00110 { \ 00111 size_t ntmp; \ 00112 int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp ); \ 00113 if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); } \ 00114 vals.resize( ntmp ); \ 00115 size_t ntmp1 = 0; \ 00116 dvfail = nc_get_vara_float( ncFile, id, &ntmp1, &ntmp, &vals[0] ); \ 00117 if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << name ); } \ 00118 } \ 00119 } 00120 00121 int main( int argc, char* argv[] ) 00122 { 00123 00124 ProgOptions opts; 00125 00126 std::string inputfile, outfile( "out.h5m" ), netcdfFile, variable_name, sefile_name; 00127 00128 opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile ); 00129 opts.addOpt< std::string >( "netcdfFile,n", "netcdf file aligned with the mesh input file", &netcdfFile ); 00130 opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile ); 00131 00132 opts.addOpt< std::string >( "var,v", "variable to extract and add to output file", &variable_name ); 00133 00134 opts.addOpt< std::string >( "sefile,s", "spectral elements file (coarse SE mesh)", &sefile_name ); 00135 opts.parseCommandLine( argc, argv ); 00136 00137 ErrorCode rval; 00138 Core* mb = new Core(); 00139 00140 rval = mb->load_file( inputfile.c_str() );MB_CHK_SET_ERR( rval, "can't load input file" ); 00141 00142 std::cout << " opened " << inputfile << " with initial h5m data.\n"; 00143 // open the netcdf file, and see if it has that variable we are looking for 00144 00145 Range nodes; 00146 rval = mb->get_entities_by_dimension( 0, 0, nodes );MB_CHK_SET_ERR( rval, "can't get nodes" ); 00147 00148 Range edges; 00149 rval = mb->get_entities_by_dimension( 0, 1, edges );MB_CHK_SET_ERR( rval, "can't get edges" ); 00150 00151 Range cells; 00152 rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells" ); 00153 00154 std::cout << " it has " << nodes.size() << " vertices " << edges.size() << " edges " << cells.size() << " cells\n"; 00155 00156 // construct maps between global id and handles 00157 std::map< int, EntityHandle > vGidHandle; 00158 std::map< int, EntityHandle > eGidHandle; 00159 std::map< int, EntityHandle > cGidHandle; 00160 std::vector< int > gids; 00161 Tag gid; 00162 rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get global id tag" ); 00163 gids.resize( nodes.size() ); 00164 rval = mb->tag_get_data( gid, nodes, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on vertices" ); 00165 int i = 0; 00166 for( Range::iterator vit = nodes.begin(); vit != nodes.end(); vit++ ) 00167 { 00168 vGidHandle[gids[i++]] = *vit; 00169 } 00170 00171 gids.resize( edges.size() ); 00172 rval = mb->tag_get_data( gid, edges, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on edges" ); 00173 i = 0; 00174 for( Range::iterator vit = edges.begin(); vit != edges.end(); vit++ ) 00175 { 00176 eGidHandle[gids[i++]] = *vit; 00177 } 00178 00179 gids.resize( cells.size() ); 00180 rval = mb->tag_get_data( gid, cells, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on cells" ); 00181 i = 0; 00182 for( Range::iterator vit = cells.begin(); vit != cells.end(); vit++ ) 00183 { 00184 cGidHandle[gids[i++]] = *vit; 00185 } 00186 00187 // Open netcdf/exodus file 00188 int fail = nc_open( netcdfFile.c_str(), 0, &ncFile ); 00189 if( NC_NOWRITE != fail ) 00190 { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf II file " << netcdfFile ); } 00191 00192 std::cout << " opened " << netcdfFile << " with new data \n"; 00193 std::vector< int > dims; 00194 int nc_var; 00195 00196 size_t recs; 00197 char recname[NC_MAX_NAME + 1]; 00198 00199 std::cout << " looking for variable " << variable_name << "\n"; 00200 GET_VAR( variable_name.c_str(), nc_var, dims ); 00201 std::cout << " it has " << dims.size() << " dimensions\n"; 00202 00203 int dimIndex = -1; // index of the dimension of interest 00204 bool vertex_data = false; 00205 bool cell_data = false; 00206 for( size_t j = 0; j < dims.size(); j++ ) 00207 { 00208 fail = nc_inq_dim( ncFile, dims[j], recname, &recs ); 00209 std::string name_dim( recname ); 00210 std::cout << " dimension index " << j << " in file: " << dims[j] << " name: " << name_dim << " recs:" << recs 00211 << "\n"; 00212 if( recs == nodes.size() ) 00213 { 00214 dimIndex = j; 00215 vertex_data = true; 00216 } 00217 if( recs == cells.size() ) 00218 { 00219 dimIndex = j; 00220 cell_data = true; 00221 } 00222 } 00223 int otherDim = 1 - dimIndex; // used only if 2 dimensions ; could be 0 or 1; 00224 size_t size_tag = 1; // good for one dimension 00225 std::vector< int > evals; // size of size_tag 00226 std::vector< double > dvals; // size of size_tag 00227 nc_type dataType; 00228 // read the variable, and set it to the tag 00229 fail = nc_inq_vartype( ncFile, nc_var, &dataType ); 00230 DataType mbtype = MB_TYPE_DOUBLE; 00231 bool float_var = false; 00232 if( NC_INT == dataType ) 00233 mbtype = MB_TYPE_INTEGER; 00234 else if( NC_DOUBLE != dataType && NC_FLOAT != dataType ) 00235 MB_CHK_SET_ERR( MB_FAILURE, "unknown type" ); 00236 00237 if( NC_FLOAT == dataType ) float_var = true; 00238 00239 int time_id = -1; 00240 fail = nc_inq_varid( ncFile, "time", &time_id ); 00241 std::vector< float > times; 00242 if( NC_NOERR == fail ) 00243 { 00244 int ii; 00245 GET_1D_FLT_VAR( "time", ii, times ); 00246 } 00247 Tag newTag; 00248 00249 if( ( dims.size() >= 1 && dims.size() <= 2 ) && ( vertex_data || cell_data ) ) 00250 { 00251 00252 if( dims.size() == 2 ) { fail = nc_inq_dim( ncFile, dims[otherDim], recname, &size_tag ); } 00253 00254 int def_val = 0; 00255 rval = mb->tag_get_handle( variable_name.c_str(), (int)size_tag, mbtype, newTag, MB_TAG_CREAT | MB_TAG_DENSE, 00256 &def_val );MB_CHK_SET_ERR( rval, "can't define new tag" ); 00257 00258 if( NC_INT == dataType ) 00259 { 00260 std::vector< int > vals; 00261 if( 1 == dims.size() ) 00262 { 00263 GET_1D_INT_VAR( variable_name.c_str(), dims[0], vals ); 00264 if( vertex_data ) 00265 { 00266 for( size_t k = 0; k < vals.size(); k++ ) 00267 { 00268 EntityHandle vh = vGidHandle[k + 1]; 00269 rval = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00270 } 00271 } 00272 else // cell_data 00273 { 00274 for( size_t k = 0; k < vals.size(); k++ ) 00275 { 00276 EntityHandle ch = cGidHandle[k + 1]; // cell handle 00277 rval = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on cell" ); 00278 } 00279 } 00280 } 00281 else // dims.size() == 2 00282 { 00283 // Single var for all coords 00284 size_t start[2] = { 0, 0 }, count[2] = { 1, 1 }; 00285 count[dimIndex] = recs; 00286 count[otherDim] = size_tag; 00287 vals.resize( recs * size_tag ); 00288 fail = nc_get_vara_int( ncFile, nc_var, start, count, &vals[0] ); 00289 evals.resize( size_tag ); 00290 00291 if( vertex_data ) 00292 { 00293 for( size_t k = 0; k < recs; k++ ) 00294 { 00295 EntityHandle vh = vGidHandle[k + 1]; 00296 size_t start_in_vals = k * size_tag, stride = 1; 00297 for( size_t j = 0; j < size_tag; j++ ) 00298 evals[j] = vals[start_in_vals + j * stride]; 00299 rval = mb->tag_set_data( newTag, &vh, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00300 } 00301 } 00302 else // cell_data 00303 { 00304 for( size_t k = 0; k < recs; k++ ) 00305 { 00306 EntityHandle ch = cGidHandle[k + 1]; // cell handle 00307 size_t start_in_vals = k * size_tag, stride = 1; 00308 for( size_t j = 0; j < size_tag; j++ ) 00309 evals[j] = vals[start_in_vals + j * stride]; 00310 rval = mb->tag_set_data( newTag, &ch, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" ); 00311 } 00312 } 00313 } 00314 } 00315 else 00316 { 00317 std::vector< double > vals; 00318 if( 1 == dims.size() ) 00319 { 00320 GET_1D_DBL_VAR( variable_name.c_str(), dims[0], vals ); 00321 if( vertex_data ) 00322 { 00323 for( size_t k = 0; k < vals.size(); k++ ) 00324 { 00325 EntityHandle vh = vGidHandle[k + 1]; 00326 rval = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00327 } 00328 } 00329 else // cell_data 00330 { 00331 for( size_t k = 0; k < vals.size(); k++ ) 00332 { 00333 EntityHandle ch = cGidHandle[k + 1]; // cell handle 00334 rval = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00335 } 00336 } 00337 } 00338 else // dims.size() == 2 00339 { 00340 // Single var for all coords 00341 size_t start[2] = { 0, 0 }, count[2] = { 1, 1 }; 00342 count[dimIndex] = recs; 00343 count[otherDim] = size_tag; 00344 vals.resize( recs * size_tag ); 00345 fail = nc_get_vara_double( ncFile, nc_var, start, count, &vals[0] ); 00346 dvals.resize( size_tag ); 00347 00348 if( vertex_data ) 00349 { 00350 for( size_t k = 0; k < recs; k++ ) 00351 { 00352 EntityHandle vh = vGidHandle[k + 1]; 00353 size_t start_in_vals = k * size_tag, stride = 1; 00354 for( size_t j = 0; j < size_tag; j++ ) 00355 dvals[j] = vals[start_in_vals + j * stride]; 00356 rval = mb->tag_set_data( newTag, &vh, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00357 } 00358 } 00359 else // cell_data 00360 { 00361 for( size_t k = 0; k < recs; k++ ) 00362 { 00363 EntityHandle ch = cGidHandle[k + 1]; // cell handle 00364 size_t start_in_vals = k * size_tag, stride = 1; 00365 for( size_t j = 0; j < size_tag; j++ ) 00366 dvals[j] = vals[start_in_vals + j * stride]; 00367 rval = mb->tag_set_data( newTag, &ch, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" ); 00368 } 00369 } 00370 } 00371 } 00372 } 00373 else if( ( dims.size() == 3 ) && vertex_data && dimIndex == 2 && 00374 mbtype == MB_TYPE_DOUBLE ) // the last one is the vertex 00375 { 00376 // the case when the last dim is ncol (for homme type mesh)) 00377 size_t dim0, dim1; // dim 2 is ncol.. 00378 char recname0[NC_MAX_NAME + 1], recname1[NC_MAX_NAME + 1]; 00379 fail = nc_inq_dim( ncFile, dims[0], recname0, &dim0 ); 00380 fail = nc_inq_dim( ncFile, dims[1], recname1, &dim1 ); 00381 std::string name0( recname0 ); 00382 std::string name1( recname1 ); 00383 std::string timestr( "time" ); 00384 size_t start[3] = { 0, 0, 0 }; 00385 size_t count[3] = { 1, 0, 0 }; 00386 count[1] = dim1; 00387 count[2] = nodes.size(); 00388 // create a few tags, with name inserted 00389 if( name0.compare( "time" ) == 0 ) 00390 { 00391 00392 std::vector< double > dvalues( dim1 * count[2] ); 00393 std::vector< float > fvalues( dim1 * count[2] ); 00394 std::vector< double > transp( dim1 * count[2] ); 00395 // get the variable time 00396 for( size_t k = 0; k < dim0; k++ ) 00397 { 00398 // create a tag for each time, and 00399 std::stringstream tag_name; 00400 tag_name << variable_name << "_t" << times[k]; 00401 std::vector< double > defvals( dim1, 0. ); 00402 rval = mb->tag_get_handle( tag_name.str().c_str(), (int)dim1, mbtype, newTag, 00403 MB_TAG_CREAT | MB_TAG_DENSE, &defvals[0] );MB_CHK_SET_ERR( rval, "can't define new tag" ); 00404 start[0] = k; 00405 00406 if( float_var ) 00407 { 00408 fail = nc_get_vara_float( ncFile, nc_var, start, count, &fvalues[0] ); 00409 // now arrange them in a tag, transpose data 00410 for( size_t ii = 0; ii < dim1; ii++ ) 00411 for( size_t j = 0; j < count[2]; j++ ) 00412 { 00413 transp[j * dim1 + ii] = fvalues[ii * count[2] + j]; 00414 } 00415 } 00416 else // double 00417 { 00418 fail = nc_get_vara_double( ncFile, nc_var, start, count, &dvalues[0] ); 00419 // now arrange them in a tag, transpose data 00420 for( size_t ii = 0; ii < dim1; ii++ ) 00421 for( size_t j = 0; j < count[2]; j++ ) 00422 { 00423 transp[j * dim1 + ii] = dvalues[ii * count[2] + j]; 00424 } 00425 } 00426 for( size_t ii = 0; ii < nodes.size(); ii++ ) 00427 { 00428 EntityHandle vh = vGidHandle[ii + 1]; 00429 rval = mb->tag_set_data( newTag, &vh, 1, &transp[ii * dim1] );MB_CHK_SET_ERR( rval, "can't set tag on nodes" ); 00430 } 00431 } 00432 } 00433 } 00434 rval = mb->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" ); 00435 std::cout << " wrote file " << outfile << "\n"; 00436 00437 // now, if s option, load the coarse mesh and put data on each element, according to a matrix 00438 if( !sefile_name.empty() ) 00439 { 00440 // load the file, check for GLOBAL_DOFS tag, and create a new tag with the data associated 00441 Core* mb2 = new Core(); 00442 rval = mb2->load_file( sefile_name.c_str() );MB_CHK_SET_ERR( rval, "can't load spectral element file" ); 00443 std::cout << " loaded spectral file " << sefile_name << "\n"; 00444 // look for GLOBAL_DOFS tag 00445 Tag gdofeTag; 00446 rval = mb2->tag_get_handle( "GLOBAL_DOFS", gdofeTag );MB_CHK_SET_ERR( rval, "file does not have GLOBAL_DOFS file" ); 00447 int sizeTag; 00448 rval = mb2->tag_get_length( gdofeTag, sizeTag );MB_CHK_SET_ERR( rval, "can't get size of tag" ); 00449 int np = (int)sqrt( 1.0 * sizeTag ); 00450 std::cout << " size of tag: " << sizeTag << " np = " << np << "\n"; 00451 std::vector< int > gdofs; 00452 Range cells2; 00453 rval = mb2->get_entities_by_dimension( 0, 2, cells2 );MB_CHK_SET_ERR( rval, "can't get cells on spectral mesh" ); 00454 gdofs.resize( cells2.size() * sizeTag ); 00455 rval = mb2->tag_get_data( gdofeTag, cells2, &gdofs[0] );MB_CHK_SET_ERR( rval, "can't get global dofs tag" ); 00456 // create a new tag for element data arranged as DFIELD 00457 00458 std::vector< double > dfield; 00459 dfield.resize( sizeTag, 0.0 ); 00460 Tag newTag2; 00461 rval = mb2->tag_get_handle( variable_name.c_str(), (int)sizeTag, mbtype, newTag2, MB_TAG_CREAT | MB_TAG_DENSE, 00462 &dfield[0] );MB_CHK_SET_ERR( rval, "can't define new tag" ); 00463 00464 int i1 = 0; // index in the gdofs array, per element 00465 00466 // get the tag values from the other moab core, for newTag 00467 int dataTagLen; 00468 rval = mb->tag_get_length( newTag, dataTagLen );MB_CHK_SET_ERR( rval, "can't get size of newTag" ); 00469 // 00470 std::vector< double > oldData; 00471 oldData.resize( dataTagLen * nodes.size() ); // 00472 00473 // get the "old" values 00474 rval = mb->tag_get_data( newTag, nodes, &oldData[0] );MB_CHK_SET_ERR( rval, "can't get old values" ); 00475 for( Range::iterator it = cells2.begin(); it != cells2.end(); it++ ) 00476 { 00477 EntityHandle cel = *it; 00478 // gdofs per element are gdofs[i:i+np*np]; 00479 for( int k = 0; k < sizeTag; k++ ) 00480 { 00481 int gdof = gdofs[i1 + k]; 00482 EntityHandle node = vGidHandle[gdof]; 00483 int index = nodes.index( node ); 00484 // get last layer of data 00485 double val = oldData[index * dataTagLen + dataTagLen - 1]; // 00486 dfield[k] = val; 00487 } 00488 i1 = i1 + sizeTag; 00489 rval = mb2->tag_set_data( newTag2, &cel, 1, &dfield[0] );MB_CHK_SET_ERR( rval, "can't set new tag" ); 00490 } 00491 00492 // write the appended file with the new field: 00493 rval = mb2->write_file( "atm2.h5m" );MB_CHK_SET_ERR( rval, "can't write new spectral file" ); 00494 std::cout << " wrote file atm2.h5m \n"; 00495 } 00496 return 0; 00497 }
1.7.6.1 
