Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
mbaddchunk.cpp
Go to the documentation of this file.
1 /*
2  * addncdata.cpp
3  * this tool will take an existing h5m file and add data from some chunk description files
4  * generated from e3sm runs;
5  * will support mainly showing the chunks in ViSit
6  *
7  * example of usage:
8  * ./mbaddchunk -i penta3d.h5m -n chunks_on_proc.txt -o penta3d_ch.h5m
9  *
10  *
11  * Basically, will output a new h5m file (penta3d_ch.h5m), which has 2 extra tags, corresponding to
12  * the chunks number and processors it sits on
13  *
14  *
15  * file penta3d.h5m is obtained from the pentagons file, with a python script
16  *
17  */
18 
19 #include "moab/ProgOptions.hpp"
20 #include "moab/Core.hpp"
21 
22 #include <cmath>
23 #include <sstream>
24 #include <iostream>
25 #include <iomanip>
26 #include <fstream>
27 
28 using namespace moab;
29 using namespace std;
30 
31 int main( int argc, char* argv[] )
32 {
33 
34  ProgOptions opts;
35 
36  std::string inputfile( "penta3d.h5m" ), outfile( "penta3d_ch.h5m" ), chunkfile_name, gsmapfile;
37 
38  opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile );
39  opts.addOpt< std::string >( "chunkFile,n", "chunk file from cam run", &chunkfile_name );
40  opts.addOpt< std::string >( "gsMAPfile,g", "gsmap file", &gsmapfile );
41 
42  opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
43 
44  opts.parseCommandLine( argc, argv );
45 
46  ErrorCode rval;
47  Core* mb = new Core();
48 
49  rval = mb->load_file( inputfile.c_str() );MB_CHK_SET_ERR( rval, "can't load input file" );
50 
51  std::cout << " opened " << inputfile << " with initial h5m data.\n";
52  // open the netcdf file, and see if it has that variable we are looking for
53 
54  Range nodes;
55  rval = mb->get_entities_by_dimension( 0, 0, nodes );MB_CHK_SET_ERR( rval, "can't get nodes" );
56 
57  Range edges;
58  rval = mb->get_entities_by_dimension( 0, 1, edges );MB_CHK_SET_ERR( rval, "can't get edges" );
59 
60  Range cells;
61  rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells" );
62 
63  std::cout << " it has " << nodes.size() << " vertices " << edges.size() << " edges " << cells.size() << " cells\n";
64 
65  // construct maps between global id and handles
66  std::map< int, EntityHandle > vGidHandle;
67  std::map< int, EntityHandle > eGidHandle;
68  std::map< int, EntityHandle > cGidHandle;
69  std::vector< int > gids;
70  Tag gid;
71  rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get global id tag" );
72  gids.resize( nodes.size() );
73  rval = mb->tag_get_data( gid, nodes, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on vertices" );
74  int i = 0;
75  for( Range::iterator vit = nodes.begin(); vit != nodes.end(); vit++ )
76  {
77  vGidHandle[gids[i++]] = *vit;
78  }
79 
80  gids.resize( edges.size() );
81  rval = mb->tag_get_data( gid, edges, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on edges" );
82  i = 0;
83  for( Range::iterator vit = edges.begin(); vit != edges.end(); vit++ )
84  {
85  eGidHandle[gids[i++]] = *vit;
86  }
87 
88  gids.resize( cells.size() );
89  rval = mb->tag_get_data( gid, cells, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on cells" );
90  i = 0;
91  for( Range::iterator vit = cells.begin(); vit != cells.end(); vit++ )
92  {
93  cGidHandle[gids[i++]] = *vit;
94  }
95 
96  if( chunkfile_name.length() > 0 )
97  {
98 
99  // Open chunk file
100  ifstream inFile;
101 
102  inFile.open( chunkfile_name.c_str() );
103  if( !inFile )
104  {
105  cout << "Unable to open chunk file";
106  exit( 1 ); // terminate with error
107  }
108  Tag pTag, cTag;
109  int def_val = -1;
110  rval = mb->tag_get_handle( "ProcID", 1, MB_TYPE_INTEGER, pTag, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );MB_CHK_SET_ERR( rval, "can't define processor tag" );
111  rval = mb->tag_get_handle( "ChunkID", 1, MB_TYPE_INTEGER, cTag, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );MB_CHK_SET_ERR( rval, "can't define chunk tag" );
112 
113  int proc, lcid, ncols;
114  while( inFile >> proc )
115  {
116  inFile >> lcid >> ncols;
117  int Gid;
118  for( i = 0; i < ncols; i++ )
119  {
120  inFile >> Gid;
121  EntityHandle cell = cGidHandle[Gid];
122  rval = mb->tag_set_data( pTag, &cell, 1, &proc );MB_CHK_SET_ERR( rval, "can't set proc tag" );
123  rval = mb->tag_set_data( cTag, &cell, 1, &lcid );MB_CHK_SET_ERR( rval, "can't set chunk tag" );
124  }
125  }
126 
127  inFile.close();
128  }
129 
130  if( gsmapfile.length() > 0 )
131  {
132 
133  // Open chunk file
134  ifstream inFile;
135 
136  inFile.open( gsmapfile.c_str() );
137  if( !inFile )
138  {
139  cout << "Unable to open gsmap file";
140  exit( 1 ); // terminate with error
141  }
142  Tag pTag, cTag;
143  int def_val = -1;
144  std::string procTagName = gsmapfile + "_proc";
145  rval =
146  mb->tag_get_handle( procTagName.c_str(), 1, MB_TYPE_INTEGER, pTag, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );MB_CHK_SET_ERR( rval, "can't define processor tag" );
147  std::string segTagName = gsmapfile + "_seg";
148  rval =
149  mb->tag_get_handle( segTagName.c_str(), 1, MB_TYPE_INTEGER, cTag, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );MB_CHK_SET_ERR( rval, "can't define segment tag" );
150 
151  int compid, ngseg, gsize;
152  inFile >> compid >> ngseg >> gsize;
153  for( i = 1; i <= ngseg; i++ )
154  {
155  int start, len, pe;
156  inFile >> start >> len >> pe;
157  int Gid;
158  for( int j = 0; j < len; j++ )
159  {
160  Gid = start + j;
161  EntityHandle cell = cGidHandle[Gid];
162  rval = mb->tag_set_data( pTag, &cell, 1, &pe );MB_CHK_SET_ERR( rval, "can't set proc tag" );
163  rval = mb->tag_set_data( cTag, &cell, 1, &i );MB_CHK_SET_ERR( rval, "can't set segment tag" );
164  }
165  }
166 
167  inFile.close();
168  }
169 
170  rval = mb->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
171  std::cout << " wrote file " << outfile << "\n";
172  return 0;
173 }