AddFieldtoPG2.cpp

Go to the documentation of this file.

    1 #include "moab/ProgOptions.hpp"
    2 #include "moab/Core.hpp"
    3 #include <iostream>
    4  
    5 using namespace moab;
    6 using namespace std;
    7  
    8 int main( int argc, char* argv[] )
    9 {
   10  
   11     ProgOptions opts;
   12  
   13     std::string inputfile, outfile( "out.h5m" ), physgridfile, variable_name;
   14  
   15     opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile );
   16     opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
   17     opts.addOpt< std::string >( "phys,p", "phys grid solution filename", &physgridfile );
   18     opts.addOpt< std::string >( "var,v", "variable to extract and add to output file", &variable_name );
   19  
   20     opts.parseCommandLine( argc, argv );
   21  
   22     if( inputfile.empty() )
   23     {
   24         opts.printHelp();
   25         return 0;
   26     }
   27     ErrorCode rval;
   28     Core* mb = new Core();
   29  
   30     EntityHandle fset1, fset2;
   31     rval = mb->create_meshset( MESHSET_SET, fset1 );MB_CHK_SET_ERR( rval, "can't create mesh set" );
   32     rval = mb->load_file( inputfile.c_str(), &fset1 );MB_CHK_SET_ERR( rval, "can't load input file" );
   33  
   34     cout << " opened " << inputfile << " with initial h5m data.\n";
   35  
   36     rval = mb->create_meshset( MESHSET_SET, fset2 );MB_CHK_SET_ERR( rval, "can't create mesh set" );
   37     rval = mb->load_file( physgridfile.c_str(), &fset2 );MB_CHK_SET_ERR( rval, "can't load phys grid file" );
   38  
   39     Tag tagv;
   40     rval = mb->tag_get_handle( variable_name.c_str(), tagv );MB_CHK_SET_ERR( rval, "can't get tag handle" );
   41  
   42     Tag gitag = mb->globalId_tag();
   43  
   44     Range verts;  // from phys grid
   45     rval = mb->get_entities_by_dimension( fset2, 0, verts );MB_CHK_SET_ERR( rval, "can't get vertices" );
   46     std::vector< int > gids;
   47     gids.resize( verts.size() );
   48     rval = mb->tag_get_data( gitag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "can't get gi tag values" );
   49     std::vector< double > valsTag;
   50     valsTag.resize( verts.size() );
   51     rval = mb->tag_get_data( tagv, verts, &valsTag[0] );MB_CHK_SET_ERR( rval, "can't get tag vals" );
   52     Range cells;
   53  
   54     rval = mb->get_entities_by_dimension( fset1, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells" );
   55  
   56     std::map< int, double > valsByID;
   57     for( int i = 0; i < (int)gids.size(); i++ )
   58         valsByID[gids[i]] = valsTag[i];
   59  
   60     // set now cells values
   61     std::vector< int > cellsIds;
   62     cellsIds.resize( cells.size() );
   63     rval = mb->tag_get_data( gitag, cells, &cellsIds[0] );MB_CHK_SET_ERR( rval, "can't get cells ids" );
   64     for( int i = 0; i < (int)cells.size(); i++ )
   65     {
   66         valsTag[i] = valsByID[cellsIds[i]];
   67     }
   68     rval = mb->tag_set_data( tagv, cells, &valsTag[0] );MB_CHK_SET_ERR( rval, "can't set  cells tags" );
   69  
   70     rval = mb->write_file( outfile.c_str(), 0, 0, &fset1, 1 );MB_CHK_SET_ERR( rval, "can't write file" );
   71  
   72     return 0;
   73 }