copyPartition.cpp

Go to the documentation of this file.

    1 /** @example copyPartition.cpp  copy partition info on a mesh file from a point cloud
    2  * this tool will take an existing h5m  phys grid partition file (point cloud) and copy the
    3  * partition information on a pg2 mesh file, for better viewing with VisIt
    4  *
    5  * example of usage:
    6  * ./copyPartition -p AtmPhys.h5m -g wholeATM_PG2.h5m -o atm_PG2_part.h5m
    7  *
    8  * the *PG2" style atm mesh is available in E3SM only for pg2 runs, something like
    9  *  --res ne30pg2_r05_oECv3_ICG --compset A_WCYCL1850S_CMIP6
   10  *  or
   11  *  --res ne4pg2_ne4pg2 --compset FC5AV1C-L
   12  */
   13 #include "moab/ProgOptions.hpp"
   14 #include "moab/Core.hpp"
   15  
   16 #include <cmath>
   17 #include <sstream>
   18  
   19 using namespace moab;
   20  
   21 int main( int argc, char* argv[] )
   22 {
   23  
   24     ProgOptions opts;
   25  
   26     std::string physfile, pg2file, outfile;
   27  
   28     opts.addOpt< std::string >( "physgridfile,p", "phys grid filename", &physfile );
   29     opts.addOpt< std::string >( "pg2file,g", "pg2 mesh file", &pg2file );
   30     opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
   31  
   32     opts.parseCommandLine( argc, argv );
   33  
   34     std::cout << "phys grid cloud file: " << physfile << "\n";
   35     std::cout << "pg2 mesh file: " << pg2file << "\n";
   36     std::cout << "output file: " << outfile << "\n";
   37  
   38     if( physfile.empty() )
   39     {
   40         opts.printHelp();
   41         return 0;
   42     }
   43     ErrorCode rval;
   44     Core* mb = new Core();
   45  
   46     rval = mb->load_file( physfile.c_str() );MB_CHK_SET_ERR( rval, "can't load phys grid file" );
   47  
   48     Core* mb2 = new Core();
   49     rval      = mb2->load_file( pg2file.c_str() );MB_CHK_SET_ERR( rval, "can't load pg2 mesh file" );
   50  
   51     Tag globalIDTag1 = mb->globalId_tag();
   52     Tag parti;
   53     rval = mb->tag_get_handle( "partition", parti );MB_CHK_SET_ERR( rval, "can't get partition tag phys grid mesh " );
   54  
   55     Tag globalIDTag2 = mb2->globalId_tag();
   56  
   57     Range verts1;
   58     rval = mb->get_entities_by_dimension( 0, 0, verts1 );MB_CHK_SET_ERR( rval, "can't get vertices " );
   59  
   60     std::vector< int > partValues;
   61     partValues.resize( verts1.size() );
   62     rval = mb->tag_get_data( parti, verts1, &partValues[0] );MB_CHK_SET_ERR( rval, "can't get parts values on vertices " );
   63  
   64     Range cells;
   65     rval = mb2->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "can't get 2d cells " );
   66     std::vector< int > globalIdsCells;
   67     globalIdsCells.resize( cells.size() );
   68     rval = mb2->tag_get_data( globalIDTag2, cells, &globalIdsCells[0] );MB_CHK_SET_ERR( rval, "can't get global ids cells " );
   69  
   70     std::vector< int > globalIdsVerts;
   71     globalIdsVerts.resize( verts1.size() );
   72     rval = mb->tag_get_data( globalIDTag1, verts1, &globalIdsVerts[0] );MB_CHK_SET_ERR( rval, "can't get global ids cells " );
   73  
   74     Tag partTag;
   75     rval = mb2->tag_get_handle( "PARALLEL_PARTITION", partTag );MB_CHK_SET_ERR( rval, "can't partition tag " );
   76  
   77     Range sets;
   78     rval = mb2->get_entities_by_type_and_tag( 0, MBENTITYSET, &partTag, NULL, 1, sets );MB_CHK_ERR( rval );
   79  
   80     std::vector< int > setValues;
   81     setValues.resize( sets.size() );
   82     rval = mb2->tag_get_data( partTag, sets, &setValues[0] );MB_CHK_ERR( rval );
   83  
   84     std::map< int, EntityHandle > valToSet;
   85     int i = 0;
   86     for( Range::iterator st = sets.begin(); st != sets.end(); ++st, i++ )
   87     {
   88         valToSet[setValues[i]] = *st;
   89     }
   90     // now, every cell will be put into one set, by looking at the global id of cell
   91  
   92     std::map< int, EntityHandle > gidToCell;
   93     i = 0;
   94     for( Range::iterator it = cells.begin(); it != cells.end(); ++it, i++ )
   95     {
   96         gidToCell[globalIdsCells[i]] = *it;
   97     }
   98     // empty all sets
   99     rval = mb2->clear_meshset( sets );MB_CHK_ERR( rval );
  100  
  101     // look now at parti values for vertices, and their global ids
  102     for( i = 0; i < (int)verts1.size(); i++ )
  103     {
  104         int part          = partValues[i];
  105         int gid           = globalIdsVerts[i];
  106         EntityHandle set1 = valToSet[part];
  107         EntityHandle cell = gidToCell[gid];
  108         rval              = mb2->add_entities( set1, &cell, 1 );MB_CHK_ERR( rval );
  109     }
  110  
  111     rval = mb2->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
  112  
  113     delete mb;
  114     delete mb2;
  115  
  116     return 0;
  117 }