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