ExtractLand.cpp

Go to the documentation of this file.

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