Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
addPCdata.cpp
Go to the documentation of this file.
1 /** @example addPCdata.cpp Add point cloud data
2  * this tool will take an existing h5m fine atm mesh file and add data from an h5m type file with
3  * point cloud mesh will support mainly showing the data with Visit
4  *
5  * example of usage:
6  * ./mbaddpcdata -i wholeFineATM.h5m -s wholeLND_proj01.h5m -o atm_a2l.h5m -v a2lTbot_proj
7  *
8  * it should work also for pentagon file style data
9  * ./mbaddpcdata -i <MOABsource>/MeshFiles/unittest/penta3d.h5m -s wholeLND_proj01.h5m -o
10  * atm_a2l.h5m -v a2lTbot_proj -p 1
11  *
12  * Basically, will output a new h5m file (atm_a2l.h5m), which has an extra tag, corresponding to the
13  * variable a2lTbot_proj from the file wholeLND_proj01.h5m; matching is based on the global ids
14  * between what we think is the order on the original file (wholeFineATM.h5m) and the order of
15  * surfdata_ne11np4_simyr1850_c160614.nc
16  *
17  * file wholeFineATM.h5m is obtained from a coupled run in e3sm, with the ne 11, np 4,
18  *
19  */
20 #include "moab/ProgOptions.hpp"
21 #include "moab/Core.hpp"
22 
23 #include <cmath>
24 #include <sstream>
25 
26 using namespace moab;
27 
28 int main( int argc, char* argv[] )
29 {
30 
31  ProgOptions opts;
32 
33  std::string inputfile, outfile( "out.h5m" ), sourcefile, variable_name;
34  int dual_mesh = 0;
35  opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile );
36  opts.addOpt< std::string >( "source,s", "h5m file aligned with the mesh input file", &sourcefile );
37  opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
38 
39  opts.addOpt< std::string >( "var,v", "variable to extract and add to output file", &variable_name );
40  opts.addOpt< int >( "pentagon,p", "switch for dual mesh ", &dual_mesh );
41 
42  opts.parseCommandLine( argc, argv );
43 
44  std::cout << "input file: " << inputfile << "\n";
45  std::cout << "source file: " << sourcefile << "\n";
46  std::cout << "variable to copy: " << variable_name << "\n";
47  std::cout << "output file: " << outfile << "\n";
48 
49  if( inputfile.empty() )
50  {
51  opts.printHelp();
52  return 0;
53  }
54  ErrorCode rval;
55  Core* mb = new Core();
56 
57  rval = mb->load_file( inputfile.c_str() );MB_CHK_SET_ERR( rval, "can't load input file" );
58 
59  Core* mb2 = new Core();
60  rval = mb2->load_file( sourcefile.c_str() );MB_CHK_SET_ERR( rval, "can't load source file" );
61 
62  Tag sourceTag;
63  rval = mb2->tag_get_handle( variable_name.c_str(), sourceTag );MB_CHK_SET_ERR( rval, "can't get tag from file" );
64 
65  int sizeTag = 0;
66  rval = mb2->tag_get_length( sourceTag, sizeTag );MB_CHK_SET_ERR( rval, "can't get size of tag " );
67 
68  DataType type;
69  rval = mb2->tag_get_data_type( sourceTag, type );MB_CHK_SET_ERR( rval, "can't get type of tag " );
70 
71  int sizeInBytes = 0;
72  rval = mb2->tag_get_bytes( sourceTag, sizeInBytes );MB_CHK_SET_ERR( rval, "can't get size in bytes of tag " );
73 
74  Tag newTag;
75  /*
76  * ErrorCode tag_get_handle( const char* name,
77  int size,
78  DataType type,
79  Tag& tag_handle,
80  unsigned flags = 0,
81  const void* default_value = 0 )
82  */
83  void* defVal;
84  int defInt = -100000;
85  double defDouble = -1.e10;
86  if( type == MB_TYPE_INTEGER )
87  {
88  defVal = (void*)( &defInt );
89  }
90  else if( type == MB_TYPE_DOUBLE )
91  {
92  defVal = (void*)( &defDouble );
93  }
94 
95  rval = mb->tag_get_handle( variable_name.c_str(), sizeTag, type, newTag, MB_TAG_DENSE | MB_TAG_CREAT, defVal );MB_CHK_SET_ERR( rval, "can't create new tag " );
96 
97  // get vertices on ini mesh; get global id on ini mesh
98  // get global id on source mesh
99  Tag gid;
100  Tag gid2;
101  rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get GLOBAL_ID tag on ini mesh " );
102  rval = mb2->tag_get_handle( "GLOBAL_ID", gid2 );MB_CHK_SET_ERR( rval, "can't get GLOBAL_ID tag on source mesh " );
103 
104  // get vertices on ini mesh; build
105  Range iniVerts;
106  rval = mb->get_entities_by_dimension( 0, 0, iniVerts );MB_CHK_SET_ERR( rval, "can't get verts on initial mesh " );
107  if( 0 != dual_mesh )
108  {
109  // verts will be polygons
110  rval = mb->get_entities_by_dimension( 0, 2, iniVerts );MB_CHK_SET_ERR( rval, "can't get polygons on initial mesh " );
111  }
112  std::vector< int > gids;
113  gids.resize( iniVerts.size() );
114  rval = mb->tag_get_data( gid, iniVerts, &( gids[0] ) );MB_CHK_SET_ERR( rval, "can't get gid on initial verts " );
115  // build now the map
116  std::map< int, EntityHandle > fromGidToEh;
117  int i = 0;
118  for( Range::iterator vit = iniVerts.begin(); vit != iniVerts.end(); vit++, i++ )
119  {
120  fromGidToEh[gids[i]] = *vit;
121  }
122  // now get the source verts, and tags, and set it on new mesh
123 
124  char* valTag = new char[sizeInBytes];
125 
126  std::cout << " size of tag in bytes:" << sizeInBytes << "\n";
127  Range sourceVerts;
128  rval = mb2->get_entities_by_dimension( 0, 0, sourceVerts );MB_CHK_SET_ERR( rval, "can't get verts on source mesh " );
129  for( Range::iterator sit = sourceVerts.begin(); sit != sourceVerts.end(); sit++ )
130  {
131  int globalId = 0;
132  EntityHandle sourceHandle = *sit;
133  rval = mb2->tag_get_data( gid2, &sourceHandle, 1, &globalId );MB_CHK_SET_ERR( rval, "can't get id on source mesh " );
134  // iniVert could be a polygon, actually
135  EntityHandle iniVert = fromGidToEh[globalId];
136  // find the value of source tag
137  rval = mb2->tag_get_data( sourceTag, &sourceHandle, 1, (void*)valTag );MB_CHK_SET_ERR( rval, "can't get value on source tag " );
138 
139  rval = mb->tag_set_data( newTag, &iniVert, 1, (void*)valTag );MB_CHK_SET_ERR( rval, "can't set value on initial mesh, new tag " );
140  }
141 
142  // save file
143  rval = mb->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
144 
145  delete[] valTag;
146  delete mb;
147  delete mb2;
148 
149  return 0;
150 }