ExtractLand.cpp File Reference

#include "moab/ProgOptions.hpp"
#include "moab/Core.hpp"
#include <cmath>
#include <sstream>

Include dependency graph for ExtractLand.cpp:

Go to the source code of this file.

Functions
int	main (int argc, char *argv[])

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Examples

ExtractLand.cpp.

Definition at line 22 of file ExtractLand.cpp.

   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 }

References moab::Core::add_entities(), ProgOptions::addOpt(), moab::Range::begin(), moab::Core::create_meshset(), moab::Core::delete_entities(), moab::Range::end(), ErrorCode, moab::Core::get_entities_by_dimension(), moab::Core::globalId_tag(), moab::Range::insert(), moab::Core::load_file(), mb, MB_CHK_SET_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MESHSET_SET, ProgOptions::parseCommandLine(), ProgOptions::printHelp(), moab::Range::size(), moab::Core::tag_get_data(), moab::Core::tag_get_handle(), moab::Core::tag_set_data(), and moab::Core::write_file().