|
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 bool use_time = false; 00268 int time_id = -1; 00269 std::vector< float > times; 00270 fail = nc_inq_varid( ncFile, "time", &time_id ); 00271 if( NC_NOERR == fail ) 00272 { 00273 use_time = true; 00274 int ii; 00275 GET_1D_FLT_VAR( "time", ii, times ); 00276 } 00277 Tag newTag; 00278 00279 if( ( dims.size() >= 1 && dims.size() <= 2 ) && ( vertex_data || cell_data ) ) 00280 { 00281 00282 if( dims.size() == 2 ) 00283 { 00284 fail = nc_inq_dim( ncFile, dims[otherDim], recname, &size_tag ); 00285 if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" ); 00286 } 00287 00288 int def_val = 0; 00289 rval = mb->tag_get_handle( variable_name.c_str(), (int)size_tag, mbtype, newTag, MB_TAG_CREAT | MB_TAG_DENSE, 00290 &def_val );MB_CHK_SET_ERR( rval, "can't define new tag" ); 00291 00292 if( NC_INT == dataType ) 00293 { 00294 std::vector< int > vals; 00295 if( 1 == dims.size() ) 00296 { 00297 GET_1D_INT_VAR( variable_name.c_str(), dims[0], vals ); 00298 if( vertex_data ) 00299 { 00300 for( size_t k = 0; k < vals.size(); k++ ) 00301 { 00302 EntityHandle vh = vGidHandle[k + 1]; 00303 rval = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00304 } 00305 } 00306 else // cell_data 00307 { 00308 for( size_t k = 0; k < vals.size(); k++ ) 00309 { 00310 EntityHandle ch = cGidHandle[k + 1]; // cell handle 00311 rval = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on cell" ); 00312 } 00313 } 00314 } 00315 else // dims.size() == 2 00316 { 00317 // Single var for all coords 00318 size_t start[2] = { 0, 0 }, count[2] = { 1, 1 }; 00319 count[dimIndex] = recs; 00320 count[otherDim] = size_tag; 00321 vals.resize( recs * size_tag ); 00322 fail = nc_get_vara_int( ncFile, nc_var, start, count, &vals[0] ); 00323 evals.resize( size_tag ); 00324 00325 if( vertex_data ) 00326 { 00327 for( size_t k = 0; k < recs; k++ ) 00328 { 00329 EntityHandle vh = vGidHandle[k + 1]; 00330 size_t start_in_vals = k * size_tag, stride = 1; 00331 for( size_t j = 0; j < size_tag; j++ ) 00332 evals[j] = vals[start_in_vals + j * stride]; 00333 rval = mb->tag_set_data( newTag, &vh, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00334 } 00335 } 00336 else // cell_data 00337 { 00338 for( size_t k = 0; k < recs; k++ ) 00339 { 00340 EntityHandle ch = cGidHandle[k + 1]; // cell handle 00341 size_t start_in_vals = k * size_tag, stride = 1; 00342 for( size_t j = 0; j < size_tag; j++ ) 00343 evals[j] = vals[start_in_vals + j * stride]; 00344 rval = mb->tag_set_data( newTag, &ch, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" ); 00345 } 00346 } 00347 } 00348 } 00349 else 00350 { 00351 std::vector< double > vals; 00352 if( 1 == dims.size() ) 00353 { 00354 GET_1D_DBL_VAR( variable_name.c_str(), dims[0], vals ); 00355 if( vertex_data ) 00356 { 00357 for( size_t k = 0; k < vals.size(); k++ ) 00358 { 00359 EntityHandle vh = vGidHandle[k + 1]; 00360 rval = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00361 } 00362 } 00363 else // cell_data 00364 { 00365 for( size_t k = 0; k < vals.size(); k++ ) 00366 { 00367 EntityHandle ch = cGidHandle[k + 1]; // cell handle 00368 rval = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00369 } 00370 } 00371 } 00372 else // dims.size() == 2 00373 { 00374 // Single var for all coords 00375 size_t start[2] = { 0, 0 }, count[2] = { 1, 1 }; 00376 count[dimIndex] = recs; 00377 count[otherDim] = size_tag; 00378 vals.resize( recs * size_tag ); 00379 fail = nc_get_vara_double( ncFile, nc_var, start, count, &vals[0] ); 00380 dvals.resize( size_tag ); 00381 00382 if( vertex_data ) 00383 { 00384 for( size_t k = 0; k < recs; k++ ) 00385 { 00386 EntityHandle vh = vGidHandle[k + 1]; 00387 size_t start_in_vals = k * size_tag, stride = 1; 00388 for( size_t j = 0; j < size_tag; j++ ) 00389 dvals[j] = vals[start_in_vals + j * stride]; 00390 rval = mb->tag_set_data( newTag, &vh, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" ); 00391 } 00392 } 00393 else // cell_data 00394 { 00395 for( size_t k = 0; k < recs; k++ ) 00396 { 00397 EntityHandle ch = cGidHandle[k + 1]; // cell handle 00398 size_t start_in_vals = k * size_tag, stride = 1; 00399 for( size_t j = 0; j < size_tag; j++ ) 00400 dvals[j] = vals[start_in_vals + j * stride]; 00401 rval = mb->tag_set_data( newTag, &ch, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" ); 00402 } 00403 } 00404 } 00405 } 00406 } 00407 else if( ( dims.size() == 3 ) && vertex_data && dimIndex == 2 && 00408 mbtype == MB_TYPE_DOUBLE && use_time) // the last one is the vertex 00409 { 00410 // the case when the last dim is ncol (for homme type mesh)) 00411 size_t dim0, dim1; // dim 2 is ncol.. 00412 char recname0[NC_MAX_NAME + 1], recname1[NC_MAX_NAME + 1]; 00413 fail = nc_inq_dim( ncFile, dims[0], recname0, &dim0 ); 00414 if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" ); 00415 fail = nc_inq_dim( ncFile, dims[1], recname1, &dim1 ); 00416 if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" ); 00417 std::string name0( recname0 ); 00418 std::string name1( recname1 ); 00419 std::string timestr( "time" ); 00420 size_t start[3] = { 0, 0, 0 }; 00421 size_t count[3] = { 1, 0, 0 }; 00422 count[1] = dim1; 00423 count[2] = nodes.size(); 00424 // create a few tags, with name inserted 00425 if( name0.compare( "time" ) == 0 ) 00426 { 00427 00428 std::vector< double > dvalues( dim1 * count[2] ); 00429 std::vector< float > fvalues( dim1 * count[2] ); 00430 std::vector< double > transp( dim1 * count[2] ); 00431 // get the variable time 00432 for( size_t k = 0; k < dim0; k++ ) 00433 { 00434 // create a tag for each time, and 00435 std::stringstream tag_name; 00436 tag_name << variable_name << "_t" << times[k]; 00437 std::vector< double > defvals( dim1, 0. ); 00438 rval = mb->tag_get_handle( tag_name.str().c_str(), (int)dim1, mbtype, newTag, 00439 MB_TAG_CREAT | MB_TAG_DENSE, &defvals[0] );MB_CHK_SET_ERR( rval, "can't define new tag" ); 00440 start[0] = k; 00441 00442 if( float_var ) 00443 { 00444 fail = nc_get_vara_float( ncFile, nc_var, start, count, &fvalues[0] ); 00445 // now arrange them in a tag, transpose data 00446 for( size_t ii = 0; ii < dim1; ii++ ) 00447 for( size_t j = 0; j < count[2]; j++ ) 00448 { 00449 transp[j * dim1 + ii] = fvalues[ii * count[2] + j]; 00450 } 00451 } 00452 else // double 00453 { 00454 fail = nc_get_vara_double( ncFile, nc_var, start, count, &dvalues[0] ); 00455 // now arrange them in a tag, transpose data 00456 for( size_t ii = 0; ii < dim1; ii++ ) 00457 for( size_t j = 0; j < count[2]; j++ ) 00458 { 00459 transp[j * dim1 + ii] = dvalues[ii * count[2] + j]; 00460 } 00461 } 00462 for( size_t ii = 0; ii < nodes.size(); ii++ ) 00463 { 00464 EntityHandle vh = vGidHandle[ii + 1]; 00465 rval = mb->tag_set_data( newTag, &vh, 1, &transp[ii * dim1] );MB_CHK_SET_ERR( rval, "can't set tag on nodes" ); 00466 } 00467 } 00468 } 00469 } 00470 rval = mb->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" ); 00471 std::cout << " wrote file " << outfile << "\n"; 00472 00473 // now, if s option, load the coarse mesh and put data on each element, according to a matrix 00474 if( !sefile_name.empty() ) 00475 { 00476 // load the file, check for GLOBAL_DOFS tag, and create a new tag with the data associated 00477 Core* mb2 = new Core(); 00478 rval = mb2->load_file( sefile_name.c_str() );MB_CHK_SET_ERR( rval, "can't load spectral element file" ); 00479 std::cout << " loaded spectral file " << sefile_name << "\n"; 00480 // look for GLOBAL_DOFS tag 00481 Tag gdofeTag; 00482 rval = mb2->tag_get_handle( "GLOBAL_DOFS", gdofeTag );MB_CHK_SET_ERR( rval, "file does not have GLOBAL_DOFS file" ); 00483 int sizeTag; 00484 rval = mb2->tag_get_length( gdofeTag, sizeTag );MB_CHK_SET_ERR( rval, "can't get size of tag" ); 00485 int np = (int)sqrt( 1.0 * sizeTag ); 00486 std::cout << " size of tag: " << sizeTag << " np = " << np << "\n"; 00487 std::vector< int > gdofs; 00488 Range cells2; 00489 rval = mb2->get_entities_by_dimension( 0, 2, cells2 );MB_CHK_SET_ERR( rval, "can't get cells on spectral mesh" ); 00490 gdofs.resize( cells2.size() * sizeTag ); 00491 rval = mb2->tag_get_data( gdofeTag, cells2, &gdofs[0] );MB_CHK_SET_ERR( rval, "can't get global dofs tag" ); 00492 // create a new tag for element data arranged as DFIELD 00493 00494 std::vector< double > dfield; 00495 dfield.resize( sizeTag, 0.0 ); 00496 Tag newTag2; 00497 rval = mb2->tag_get_handle( variable_name.c_str(), (int)sizeTag, mbtype, newTag2, MB_TAG_CREAT | MB_TAG_DENSE, 00498 &dfield[0] );MB_CHK_SET_ERR( rval, "can't define new tag" ); 00499 00500 int i1 = 0; // index in the gdofs array, per element 00501 00502 // get the tag values from the other moab core, for newTag 00503 int dataTagLen; 00504 rval = mb->tag_get_length( newTag, dataTagLen );MB_CHK_SET_ERR( rval, "can't get size of newTag" ); 00505 // 00506 std::vector< double > oldData; 00507 oldData.resize( dataTagLen * nodes.size() ); // 00508 00509 // get the "old" values 00510 rval = mb->tag_get_data( newTag, nodes, &oldData[0] );MB_CHK_SET_ERR( rval, "can't get old values" ); 00511 for( Range::iterator it = cells2.begin(); it != cells2.end(); it++ ) 00512 { 00513 EntityHandle cel = *it; 00514 // gdofs per element are gdofs[i:i+np*np]; 00515 for( int k = 0; k < sizeTag; k++ ) 00516 { 00517 int gdof = gdofs[i1 + k]; 00518 EntityHandle node = vGidHandle[gdof]; 00519 int index = nodes.index( node ); 00520 // get last layer of data 00521 double val = oldData[index * dataTagLen + dataTagLen - 1]; // 00522 dfield[k] = val; 00523 } 00524 i1 = i1 + sizeTag; 00525 rval = mb2->tag_set_data( newTag2, &cel, 1, &dfield[0] );MB_CHK_SET_ERR( rval, "can't set new tag" ); 00526 } 00527 00528 // write the appended file with the new field: 00529 rval = mb2->write_file( "atm2.h5m" );MB_CHK_SET_ERR( rval, "can't write new spectral file" ); 00530 std::cout << " wrote file atm2.h5m \n"; 00531 } 00532 return 0; 00533 } 00534
1.7.6.1 
