Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
compareFiles.cpp
Go to the documentation of this file.
00001 /*
00002  * compareFiles.cpp
00003  * this tool will take two existing h5m files, for the same mesh;
00004  *  they will both have the same GLOBAL_IDs for the elements, but the entity handles can be
00005  *  very different (depending on how the mesh was partitioned, and saved in parallel)
00006  *
00007  *  will compare then the difference between tags, and store the result on one of the files (saved
00008  * again)
00009  *
00010  *
00011  * example of usage:
00012  * ./mbcmpfiles -i file1.h5m -j file2.h5m -n <tag_name>  -o out.file
00013  *
00014  *
00015  * Basically, will output a new h5m file (out.file), which has an extra tag, corresponding to the
00016  * difference between the 2 values
00017  *
00018  */
00019 
00020 #include "moab/ProgOptions.hpp"
00021 #include "moab/Core.hpp"
00022 
00023 #include <cmath>
00024 #include <sstream>
00025 #include <iostream>
00026 #include <iomanip>
00027 #include <fstream>
00028 
00029 using namespace moab;
00030 using namespace std;
00031 
00032 int main( int argc, char* argv[] )
00033 {
00034 
00035     ProgOptions opts;
00036 
00037     std::string inputfile1( "fTargetIntx.h5m" ), inputfile2( "ocn_proj.h5m" ), outfile( "out.h5m" );
00038 
00039     std::string tag_name( "a2oTAG_proj" );
00040 
00041     opts.addOpt< std::string >( "input1,i", "input mesh filename 1", &inputfile1 );
00042     opts.addOpt< std::string >( "input2,j", "input mesh filename 2", &inputfile2 );
00043     opts.addOpt< std::string >( "tagname,n", "tag to compare", &tag_name );
00044     opts.addOpt< std::string >( "outfile,o", "output file", &outfile );
00045 
00046     opts.parseCommandLine( argc, argv );
00047 
00048     ErrorCode rval;
00049     Core* mb = new Core();
00050 
00051     rval = mb->load_file( inputfile1.c_str() );MB_CHK_SET_ERR( rval, "can't load input file 1" );
00052 
00053     Core* mb2 = new Core();
00054     rval      = mb2->load_file( inputfile2.c_str() );MB_CHK_SET_ERR( rval, "can't load input file 2" );
00055 
00056     std::cout << " opened " << inputfile1 << " and " << inputfile2 << " with initial h5m data.\n";
00057     // open the netcdf file, and see if it has that variable we are looking for
00058 
00059     Range nodes;
00060     rval = mb->get_entities_by_dimension( 0, 0, nodes );MB_CHK_SET_ERR( rval, "can't get nodes" );
00061 
00062     Range edges;
00063     rval = mb->get_entities_by_dimension( 0, 1, edges );MB_CHK_SET_ERR( rval, "can't get edges" );
00064 
00065     Range cells;
00066     rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells" );
00067 
00068     std::cout << inputfile1 << " has " << nodes.size() << " vertices " << edges.size() << " edges " << cells.size()
00069               << " cells\n";
00070 
00071     // construct maps between global id and handles
00072     std::map< int, EntityHandle > vGidHandle;
00073     std::map< int, EntityHandle > eGidHandle;
00074     std::map< int, EntityHandle > cGidHandle;
00075     std::vector< int > gids;
00076     Tag gid;
00077     rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get global id tag" );
00078     gids.resize( nodes.size() );
00079     rval = mb->tag_get_data( gid, nodes, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on vertices" );
00080     int i = 0;
00081     for( Range::iterator vit = nodes.begin(); vit != nodes.end(); vit++ )
00082     {
00083         vGidHandle[gids[i++]] = *vit;
00084     }
00085 
00086     gids.resize( edges.size() );
00087     rval = mb->tag_get_data( gid, edges, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on edges" );
00088     i = 0;
00089     for( Range::iterator vit = edges.begin(); vit != edges.end(); vit++ )
00090     {
00091         eGidHandle[gids[i++]] = *vit;
00092     }
00093 
00094     gids.resize( cells.size() );
00095     rval = mb->tag_get_data( gid, cells, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on cells" );
00096     i = 0;
00097     for( Range::iterator vit = cells.begin(); vit != cells.end(); vit++ )
00098     {
00099         cGidHandle[gids[i++]] = *vit;
00100     }
00101 
00102     Range nodes2;
00103     rval = mb2->get_entities_by_dimension( 0, 0, nodes2 );MB_CHK_SET_ERR( rval, "can't get nodes2" );
00104 
00105     Range edges2;
00106     rval = mb2->get_entities_by_dimension( 0, 1, edges2 );MB_CHK_SET_ERR( rval, "can't get edges2" );
00107 
00108     Range cells2;
00109     rval = mb2->get_entities_by_dimension( 0, 2, cells2 );MB_CHK_SET_ERR( rval, "can't get cells2" );
00110 
00111     std::cout << inputfile2 << " has " << nodes2.size() << " vertices " << edges2.size() << " edges " << cells2.size()
00112               << " cells\n";
00113 
00114     // construct maps between global id and handles
00115     std::map< int, EntityHandle > vGidHandle2;
00116     std::map< int, EntityHandle > eGidHandle2;
00117     std::map< int, EntityHandle > cGidHandle2;
00118 
00119     Tag gid2;
00120     rval = mb2->tag_get_handle( "GLOBAL_ID", gid2 );MB_CHK_SET_ERR( rval, "can't get global id tag2" );
00121     gids.resize( nodes2.size() );
00122     rval = mb2->tag_get_data( gid2, nodes2, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on vertices2" );
00123 
00124     i = 0;
00125     for( Range::iterator vit = nodes2.begin(); vit != nodes2.end(); vit++ )
00126     {
00127         vGidHandle2[gids[i++]] = *vit;
00128     }
00129 
00130     gids.resize( edges2.size() );
00131     rval = mb2->tag_get_data( gid2, edges2, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on edges2" );
00132     i = 0;
00133     for( Range::iterator vit = edges2.begin(); vit != edges2.end(); vit++ )
00134     {
00135         eGidHandle2[gids[i++]] = *vit;
00136     }
00137 
00138     gids.resize( cells2.size() );
00139     rval = mb2->tag_get_data( gid2, cells2, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on cells2" );
00140     i = 0;
00141     for( Range::iterator vit = cells2.begin(); vit != cells2.end(); vit++ )
00142     {
00143         cGidHandle2[gids[i++]] = *vit;
00144     }
00145 
00146     Tag tag;
00147     rval = mb->tag_get_handle( tag_name.c_str(), tag );MB_CHK_SET_ERR( rval, "can't get tag on file 1" );
00148 
00149     int len_tag = 0;
00150     rval        = mb->tag_get_length( tag, len_tag );MB_CHK_SET_ERR( rval, "can't get tag length on tag" );
00151     std::cout << "length tag : " << len_tag << "\n";
00152 
00153     if( cells.size() != cells2.size() )
00154     {
00155         std::cout << " meshes are different between 2 files, cells.size do not agree \n";
00156         exit( 1 );
00157     }
00158     std::vector< double > vals;
00159     vals.resize( len_tag * cells.size() );
00160     rval = mb->tag_get_data( tag, cells, &vals[0] );MB_CHK_SET_ERR( rval, "can't get tag data" );
00161 
00162     Tag tag2;
00163     rval = mb2->tag_get_handle( tag_name.c_str(), tag2 );MB_CHK_SET_ERR( rval, "can't get tag on file 2" );
00164     std::vector< double > vals2;
00165     vals2.resize( len_tag * cells2.size() );
00166     rval = mb2->tag_get_data( tag2, cells2, &vals2[0] );MB_CHK_SET_ERR( rval, "can't get tag data on file 2" );
00167 
00168     rval = mb->delete_entities( edges );MB_CHK_SET_ERR( rval, "can't delete edges from file 1" );
00169 
00170     std::string new_tag_name = tag_name + "_2";
00171     Tag newTag, newTagDiff;
00172     double def_val = -1000;
00173     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" );
00174 
00175     std::string tag_name_diff = tag_name + "_diff";
00176     rval = mb->tag_get_handle( tag_name_diff.c_str(), 1, MB_TYPE_DOUBLE, newTagDiff, MB_TAG_CREAT | MB_TAG_DENSE,
00177                                &def_val );MB_CHK_SET_ERR( rval, "can't define new tag diff" );
00178     i = 0;
00179     for( Range::iterator c2it = cells2.begin(); c2it != cells2.end(); c2it++ )
00180     {
00181         double val2 = vals2[i];
00182         int id2     = gids[i];
00183         i++;
00184         EntityHandle c1 = cGidHandle[id2];
00185 
00186         rval = mb->tag_set_data( newTag, &c1, 1, &val2 );MB_CHK_SET_ERR( rval, "can't set new tag" );
00187         int indx    = cells.index( c1 );
00188         double diff = vals[indx] - val2;
00189         rval        = mb->tag_set_data( newTagDiff, &c1, 1, &diff );MB_CHK_SET_ERR( rval, "can't set new tag" );
00190     }
00191 
00192     rval = mb->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
00193     std::cout << " wrote file " << outfile << "\n";
00194     return 0;
00195 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines