MOAB: Mesh Oriented datABase  (version 5.4.1)
ExtractLand.cpp
Go to the documentation of this file.
00001 /** @example ExtractLand.cpp
00002  * this tool will take an existing h5m  pg2 mesh (fv)  file and a land point cloud mesh
00003  *
00004  * will extract the corresponding land mesh
00005  *
00006  * example of usage:
00007  * ./ExtractLand -p wholeATM_PG2.h5m -l wholeLnd.h5m -o LandMesh.h5m
00008  *
00009  * the *PG2" style atm mesh is available in E3SM only for pg2 runs, something like
00010  *  --res ne30pg2_r05_oECv3_ICG --compset A_WCYCL1850S_CMIP6
00011  *  or
00012  *  --res ne4pg2_ne4pg2 --compset FC5AV1C-L
00013  */
00014 #include "moab/ProgOptions.hpp"
00015 #include "moab/Core.hpp"
00016 
00017 #include <cmath>
00018 #include <sstream>
00019 
00020 using namespace moab;
00021 
00022 int main( int argc, char* argv[] )
00023 {
00024 
00025     ProgOptions opts;
00026 
00027     std::string  pg2file, lndfile, outfile;
00028 
00029     opts.addOpt< std::string >( "land,l", "phys grid filename", &lndfile );
00030     opts.addOpt< std::string >( "pg2file,p", "pg2 mesh file", &pg2file );
00031     opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
00032 
00033     opts.parseCommandLine( argc, argv );
00034 
00035     std::cout << " land file " << lndfile << "\n";
00036     std::cout << "pg2 mesh file: " << pg2file << "\n";
00037     std::cout << "output file: " << outfile << "\n";
00038 
00039     if( lndfile.empty() || pg2file.empty() || outfile.empty() )
00040     {
00041         opts.printHelp();
00042         return 0;
00043     }
00044     ErrorCode rval;
00045     Core* mb = new Core();
00046 
00047     rval = mb->load_file( lndfile.c_str() );MB_CHK_SET_ERR( rval, "can't load land pc file" );
00048 
00049     Core* mb2 = new Core();
00050     rval      = mb2->load_file( pg2file.c_str() );MB_CHK_SET_ERR( rval, "can't load pg2 mesh file" );
00051 
00052     Tag globalIDTag1 = mb->globalId_tag();
00053 
00054     Tag globalIDTag2 = mb2->globalId_tag();
00055 
00056     Range verts1;
00057     rval = mb->get_entities_by_dimension( 0, 0, verts1 );MB_CHK_SET_ERR( rval, "can't get vertices " );
00058 
00059     Range cells;
00060     rval = mb2->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "can't get 2d cells " );
00061 
00062     std::vector< int > globalIdsCells;
00063     globalIdsCells.resize( cells.size() );
00064     rval = mb2->tag_get_data( globalIDTag2, cells, &globalIdsCells[0] );MB_CHK_SET_ERR( rval, "can't get global ids cells " );
00065 
00066     std::vector< int > globalIdsVerts;
00067     globalIdsVerts.resize( verts1.size() );
00068     rval = mb->tag_get_data( globalIDTag1, verts1, &globalIdsVerts[0] );MB_CHK_SET_ERR( rval, "can't get global ids cells " );
00069 
00070     // now, every cell will be put into one set, by looking at the global id of cell
00071 
00072     std::map< int, EntityHandle > gidToCell;
00073     int i = 0;
00074     for( Range::iterator it = cells.begin(); it != cells.end(); ++it, i++ )
00075     {
00076         gidToCell[globalIdsCells[i]] = *it;
00077     }
00078     // create a new file set with land cells
00079     EntityHandle fileSet;
00080     rval = mb2->create_meshset( MESHSET_SET, fileSet );MB_CHK_SET_ERR( rval, "Error creating file set" );
00081     // empty all sets
00082 
00083     Range landCells;
00084     // look now at gid values for vertices
00085     for( i = 0; i < (int)verts1.size(); i++ )
00086     {
00087         int gid           = globalIdsVerts[i];
00088         landCells.insert(gidToCell[gid]);
00089     }
00090 
00091     rval = mb2->add_entities(fileSet, landCells);;MB_CHK_SET_ERR( rval, "can't add land cells" );
00092 
00093     rval = mb2->write_file( outfile.c_str(), 0, 0, &fileSet, 1 );MB_CHK_SET_ERR( rval, "can't write file" );
00094 
00095     // write the original with mask 0/1, default -1
00096     Tag mask;
00097     double def_val=-1.;
00098     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" );
00099     for( Range::iterator it = cells.begin(); it != cells.end(); ++it, i++ )
00100     {
00101         EntityHandle cell = *it;
00102         // set to 0
00103         double val = 0.;
00104         rval = mb2->tag_set_data(mask, &cell, 1, &val); MB_CHK_SET_ERR( rval, "can't set mask tag" );
00105     }
00106 
00107     for( Range::iterator it = landCells.begin(); it != landCells.end(); ++it, i++ )
00108     {
00109         EntityHandle cell = *it;
00110         // set to 0
00111         double val = 1.;
00112         rval = mb2->tag_set_data(mask, &cell, 1, &val); MB_CHK_SET_ERR( rval, "can't set mask tag" );
00113     }
00114     rval = mb2->delete_entities(&fileSet, 1); MB_CHK_SET_ERR( rval, "can't delete set" );
00115     rval = mb2->write_file( "AtmWithLandMask.h5m" );MB_CHK_SET_ERR( rval, "can't write file" );
00116     delete mb;
00117     delete mb2;
00118 
00119     return 0;
00120 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines