Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
CompareH5mFiles.cpp
Go to the documentation of this file.
1 /*
2  * CompareH5mFiles.cpp
3  *
4  * Created on: Sep 26, 2023
5  */
6 
7 #include "moab/Core.hpp"
8 #include "moab/ProgOptions.hpp"
9 #include "moab/ReadUtilIface.hpp"
10 #include "math.h"
11 #include <map>
12 #include <iostream>
13 #include <cassert>
14 
15 using namespace moab;
16 using namespace std;
17 
18 int main( int argc, char** argv )
19 {
20 
21  std::string file1( "NE4pg2_imp/run/OcnCplAftMm25.h5m" );
22  std::string file2( "NE4pg2_imp2/run/OcnCplAftMm25.h5m" );
23 
24  ProgOptions opts;
25 
26  int dim = 2;
27  opts.addOpt< std::string >( "file1,f", "first file", &file1 );
28  opts.addOpt< std::string >( "file2,g", "second file", &file2 );
29 
30  opts.addOpt<void>( "global,G", "use global id for correspondence" );
31 
32 
33  opts.parseCommandLine( argc, argv );
34 
35  bool use_global_id = opts.numOptSet("global") > 0;
36 
37  ErrorCode rval;
38  Core* mb = new Core();
39 
40  rval = mb->load_file( file1.c_str() );MB_CHK_SET_ERR( rval, "can't load file1" );
41 
42  Core* mb2 = new Core();
43  rval = mb2->load_file( file2.c_str() );MB_CHK_SET_ERR( rval, "can't load file2" );
44 
45  std::vector< Tag > list1;
46  rval = mb->tag_get_tags( list1 );MB_CHK_SET_ERR( rval, "can't get tags 1" );
47 
48  Tag gidTag = mb->globalId_tag();
49  Tag gidTag2 = mb2->globalId_tag();
50 
51  Range cells1;
52  rval = mb->get_entities_by_dimension( 0, dim, cells1 );MB_CHK_SET_ERR( rval, "can't get cells 1" );
53  if( cells1.size() == 0 )
54  {
55  dim = 0;
56  rval = mb->get_entities_by_dimension( 0, dim, cells1 );MB_CHK_SET_ERR( rval, "can't get cells 1" );
57  }
58 
59 
60  Range cells2;
61  rval = mb2->get_entities_by_dimension( 0, dim, cells2 );MB_CHK_SET_ERR( rval, "can't get cells 2" );
62 
63  if( cells1.size() != cells2.size() ) MB_CHK_SET_ERR( MB_FAILURE, "different size models " );
64  vector< double > vals1, vals2;
65  vals1.resize( cells1.size() );
66  vals2.resize( cells2.size() );
67 
68  std::vector<int> gids, gids2;
69  std::map<int, int> gidMap2;
70  if (use_global_id)
71  {
72  gids.resize(cells1.size());
73  rval = mb->tag_get_data( gidTag, cells1, &gids[0] );MB_CHK_SET_ERR( rval, "can't get ids for cells1" );
74  gids2.resize(cells2.size());
75  rval = mb2->tag_get_data( gidTag2, cells2, &gids2[0] );MB_CHK_SET_ERR( rval, "can't get ids for cells2" );
76  for (int i=0; i<cells2.size(); i++)
77  {
78  gidMap2[ gids2[i] ] = i;
79  }
80  }
81 
82  int k = 0; // number of different fields
83  int k1 = 0; // number of exactly the same fields
84  std::cout << " compare files: " << file1 << " and " << file2 << " dimension entity: " << dim << "\n";
85  std::vector< std::string > same_fields;
86  std::vector<Tag> diffTags;
87  for( size_t i = 0; i < list1.size(); i++ )
88  {
89  Tag tag = list1[i];
90  std::string name;
91  rval = mb->tag_get_name( tag, name );MB_CHK_SET_ERR( rval, "can't get tag name" );
92  DataType type;
93  rval = mb->tag_get_data_type( tag, type );MB_CHK_SET_ERR( rval, "can't get tag data type" );
94  if( MB_TYPE_DOUBLE != type ) continue;
95  TagType tag_type;
96  rval = mb->tag_get_type( tag, tag_type );MB_CHK_SET_ERR( rval, "can't get tag type" );
97  if( MB_TAG_DENSE != tag_type ) continue;
98  int length = 0;
99  rval = mb->tag_get_length( tag, length );MB_CHK_SET_ERR( rval, "can't get tag length" );
100  if( 1 != length ) continue;
101  Tag tag2;
102  rval = mb2->tag_get_handle( name.c_str(), tag2 );MB_CHK_SET_ERR( rval, "can't get tag on second model" );
103  rval = mb->tag_get_data( tag, cells1, &vals1[0] );MB_CHK_SET_ERR( rval, "can't get values on tag on model 1" );
104  rval = mb2->tag_get_data( tag2, cells2, &vals2[0] );MB_CHK_SET_ERR( rval, "can't get values on tag on model 2" );
105  double minv1, maxv1, minv2, maxv2;
106  if( vals1.size() > 0 )
107  {
108  minv1 = maxv1 = vals1[0];
109  }
110  if( vals2.size() > 0 )
111  {
112  minv2 = maxv2 = vals2[0];
113  }
114  // compute the difference
115  double sum = 0, value1, value2;
116  for( int j = 0; j < vals1.size(); j++ )
117  {
118  value1 = vals1[j];
119  if (use_global_id)
120  {
121 
122  int index2 = gidMap2[ gids[j] ];
123  value2 = vals2[index2];
124 
125  }
126  else
127  {
128  value2 = vals2[j];
129 
130 
131  }
132  sum += fabs( value1 - value2 );
133  if( value1 < minv1 ) minv1 = value1;
134  if( value1 > maxv1 ) maxv1 = value1;
135  if( value2 < minv2 ) minv2 = value2;
136  if( value2 > maxv2 ) maxv2 = value2;
137  }
138 
139  if( sum > 0. )
140  {
141  std::cout << " tag: " << name << " \t difference : " << sum << " \t min/max (" << minv1 << "/" << maxv1
142  << ") \t (" << minv2 << "/" << maxv2 << ") \n";
143  k++;
144  if (use_global_id)
145  {
146  for (int j=0; j< vals1.size(); j++)
147  {
148  int index2 = gidMap2[ gids[j] ];
149  vals1[j] -= vals2[index2];
150  }
151  }
152  else
153  {
154  for (int j=0; j< vals1.size(); j++)
155  {
156  vals1[j] -= vals2[j];
157  }
158  }
159  std::string diffTagName = name+"_diff";
160  Tag newTag;
161  rval = mb->tag_get_handle( diffTagName.c_str(), 1, MB_TYPE_DOUBLE, newTag, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_ERR( rval );
162  rval = mb->tag_set_data(newTag, cells1, &vals1[0]);MB_CHK_ERR( rval );
163  diffTags.push_back(newTag);
164 
165  }
166  else
167  {
168  same_fields.push_back( name );
169  k1++;
170  }
171  }
172  if (k>0)
173  {
174  rval = mb->write_file("diff_tags.h5m", 0, 0, 0, 0, &diffTags[0], diffTags.size() ); MB_CHK_ERR( rval );
175  }
176  std::cout << " different fields:" << k << " \n exactly the same fields:" << k1 << "\n";
177  for( size_t i = 0; i < same_fields.size(); i++ )
178  {
179  std::cout << " " << same_fields[i];
180  }
181  std::cout << "\n";
182  // Cr
183  delete mb;
184  delete mb2;
185 }