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