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