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