 |
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 /*
2 * compareFiles.cpp
3 * this tool will take two existing h5m files, for the same mesh;
4 * they will both have the same GLOBAL_IDs for the elements, but the entity handles can be
5 * very different (depending on how the mesh was partitioned, and saved in parallel)
6 *
7 * will compare then the difference between tags, and store the result on one of the files (saved
8 * again)
9 *
10 *
11 * example of usage:
12 * ./mbcmpfiles -i file1.h5m -j file2.h5m -n <tag_name> -o out.file
13 *
14 *
15 * Basically, will output a new h5m file (out.file), which has an extra tag, corresponding to the
16 * difference between the 2 values
17 *
18 */
19
20 #include "moab/ProgOptions.hpp"
21 #include "moab/Core.hpp"
22
23 #include <cmath>
24 #include <sstream>
25 #include <iostream>
26 #include <iomanip>
27 #include <fstream>
28
29 using namespace moab;
30 using namespace std;
31
32 int main( int argc, char* argv[] )
33 {
34
35 ProgOptions opts;
36
37 std::string inputfile1( "fTargetIntx.h5m" ), inputfile2( "ocn_proj.h5m" ), outfile( "out.h5m" );
38
39 std::string tag_name( "a2oTAG_proj" );
40
41 opts.addOpt< std::string >( "input1,i", "input mesh filename 1", &inputfile1 );
42 opts.addOpt< std::string >( "input2,j", "input mesh filename 2", &inputfile2 );
43 opts.addOpt< std::string >( "tagname,n", "tag to compare", &tag_name );
44 opts.addOpt< std::string >( "outfile,o", "output file", &outfile );
45
46 opts.parseCommandLine( argc, argv );
47
48 ErrorCode rval;
49 Core* mb = new Core();
50
51 rval = mb->load_file( inputfile1.c_str() );MB_CHK_SET_ERR( rval, "can't load input file 1" );
52
53 Core* mb2 = new Core();
54 rval = mb2->load_file( inputfile2.c_str() );MB_CHK_SET_ERR( rval, "can't load input file 2" );
55
56 std::cout << " opened " << inputfile1 << " and " << inputfile2 << " with initial h5m data.\n";
57 // open the netcdf file, and see if it has that variable we are looking for
58
59 Range nodes;
60 rval = mb->get_entities_by_dimension( 0, 0, nodes );MB_CHK_SET_ERR( rval, "can't get nodes" );
61
62 Range edges;
63 rval = mb->get_entities_by_dimension( 0, 1, edges );MB_CHK_SET_ERR( rval, "can't get edges" );
64
65 Range cells;
66 rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells" );
67
68 std::cout << inputfile1 << " has " << nodes.size() << " vertices " << edges.size() << " edges " << cells.size()
69 << " cells\n";
70
71 // construct maps between global id and handles
72 std::map< int, EntityHandle > vGidHandle;
73 std::map< int, EntityHandle > eGidHandle;
74 std::map< int, EntityHandle > cGidHandle;
75 std::vector< int > gids;
76 Tag gid;
77 rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get global id tag" );
78 gids.resize( nodes.size() );
79 rval = mb->tag_get_data( gid, nodes, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on vertices" );
80 int i = 0;
81 for( Range::iterator vit = nodes.begin(); vit != nodes.end(); vit++ )
82 {
83 vGidHandle[gids[i++]] = *vit;
84 }
85
86 gids.resize( edges.size() );
87 rval = mb->tag_get_data( gid, edges, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on edges" );
88 i = 0;
89 for( Range::iterator vit = edges.begin(); vit != edges.end(); vit++ )
90 {
91 eGidHandle[gids[i++]] = *vit;
92 }
93
94 gids.resize( cells.size() );
95 rval = mb->tag_get_data( gid, cells, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on cells" );
96 i = 0;
97 for( Range::iterator vit = cells.begin(); vit != cells.end(); vit++ )
98 {
99 cGidHandle[gids[i++]] = *vit;
100 }
101
102 Range nodes2;
103 rval = mb2->get_entities_by_dimension( 0, 0, nodes2 );MB_CHK_SET_ERR( rval, "can't get nodes2" );
104
105 Range edges2;
106 rval = mb2->get_entities_by_dimension( 0, 1, edges2 );MB_CHK_SET_ERR( rval, "can't get edges2" );
107
108 Range cells2;
109 rval = mb2->get_entities_by_dimension( 0, 2, cells2 );MB_CHK_SET_ERR( rval, "can't get cells2" );
110
111 std::cout << inputfile2 << " has " << nodes2.size() << " vertices " << edges2.size() << " edges " << cells2.size()
112 << " cells\n";
113
114 // construct maps between global id and handles
115 std::map< int, EntityHandle > vGidHandle2;
116 std::map< int, EntityHandle > eGidHandle2;
117 std::map< int, EntityHandle > cGidHandle2;
118
119 Tag gid2;
120 rval = mb2->tag_get_handle( "GLOBAL_ID", gid2 );MB_CHK_SET_ERR( rval, "can't get global id tag2" );
121 gids.resize( nodes2.size() );
122 rval = mb2->tag_get_data( gid2, nodes2, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on vertices2" );
123
124 i = 0;
125 for( Range::iterator vit = nodes2.begin(); vit != nodes2.end(); vit++ )
126 {
127 vGidHandle2[gids[i++]] = *vit;
128 }
129
130 gids.resize( edges2.size() );
131 rval = mb2->tag_get_data( gid2, edges2, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on edges2" );
132 i = 0;
133 for( Range::iterator vit = edges2.begin(); vit != edges2.end(); vit++ )
134 {
135 eGidHandle2[gids[i++]] = *vit;
136 }
137
138 gids.resize( cells2.size() );
139 rval = mb2->tag_get_data( gid2, cells2, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on cells2" );
140 i = 0;
141 for( Range::iterator vit = cells2.begin(); vit != cells2.end(); vit++ )
142 {
143 cGidHandle2[gids[i++]] = *vit;
144 }
145
146 Tag tag;
147 rval = mb->tag_get_handle( tag_name.c_str(), tag );MB_CHK_SET_ERR( rval, "can't get tag on file 1" );
148
149 int len_tag = 0;
150 rval = mb->tag_get_length( tag, len_tag );MB_CHK_SET_ERR( rval, "can't get tag length on tag" );
151 std::cout << "length tag : " << len_tag << "\n";
152
153 if( cells.size() != cells2.size() )
154 {
155 std::cout << " meshes are different between 2 files, cells.size do not agree \n";
156 exit( 1 );
157 }
158 std::vector< double > vals;
159 vals.resize( len_tag * cells.size() );
160 rval = mb->tag_get_data( tag, cells, &vals[0] );MB_CHK_SET_ERR( rval, "can't get tag data" );
161
162 Tag tag2;
163 rval = mb2->tag_get_handle( tag_name.c_str(), tag2 );MB_CHK_SET_ERR( rval, "can't get tag on file 2" );
164 std::vector< double > vals2;
165 vals2.resize( len_tag * cells2.size() );
166 rval = mb2->tag_get_data( tag2, cells2, &vals2[0] );MB_CHK_SET_ERR( rval, "can't get tag data on file 2" );
167
168 rval = mb->delete_entities( edges );MB_CHK_SET_ERR( rval, "can't delete edges from file 1" );
169
170 std::string new_tag_name = tag_name + "_2";
171 Tag newTag, newTagDiff;
172 double def_val = -1000;
173 rval = mb->tag_get_handle( new_tag_name.c_str(), 1, MB_TYPE_DOUBLE, newTag, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );MB_CHK_SET_ERR( rval, "can't define new tag" );
174
175 std::string tag_name_diff = tag_name + "_diff";
176 rval = mb->tag_get_handle( tag_name_diff.c_str(), 1, MB_TYPE_DOUBLE, newTagDiff, MB_TAG_CREAT | MB_TAG_DENSE,
177 &def_val );MB_CHK_SET_ERR( rval, "can't define new tag diff" );
178 i = 0;
179 for( Range::iterator c2it = cells2.begin(); c2it != cells2.end(); c2it++ )
180 {
181 double val2 = vals2[i];
182 int id2 = gids[i];
183 i++;
184 EntityHandle c1 = cGidHandle[id2];
185
186 rval = mb->tag_set_data( newTag, &c1, 1, &val2 );MB_CHK_SET_ERR( rval, "can't set new tag" );
187 int indx = cells.index( c1 );
188 double diff = vals[indx] - val2;
189 rval = mb->tag_set_data( newTagDiff, &c1, 1, &diff );MB_CHK_SET_ERR( rval, "can't set new tag" );
190 }
191
192 rval = mb->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
193 std::cout << " wrote file " << outfile << "\n";
194 return 0;
195 }
1.9.1.