Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
compareFiles.cpp File Reference
#include "moab/ProgOptions.hpp"
#include "moab/Core.hpp"
#include <cmath>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <fstream>
+ Include dependency graph for compareFiles.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 32 of file compareFiles.cpp.

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 }

References ProgOptions::addOpt(), moab::Range::begin(), moab::Core::delete_entities(), moab::Range::end(), ErrorCode, moab::Core::get_entities_by_dimension(), moab::Range::index(), moab::Core::load_file(), mb, MB_CHK_SET_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, ProgOptions::parseCommandLine(), moab::Range::size(), moab::Core::tag_get_data(), moab::Core::tag_get_handle(), moab::Core::tag_get_length(), moab::Core::tag_set_data(), and moab::Core::write_file().