Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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().