Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
copyPartition.cpp
Go to the documentation of this file.
00001 /** @example copyPartition.cpp  copy partition info on a mesh file from a point cloud
00002  * this tool will take an existing h5m  phys grid partition file (point cloud) and copy the
00003  * partition information on a pg2 mesh file, for better viewing with VisIt
00004  *
00005  * example of usage:
00006  * ./copyPartition -p AtmPhys.h5m -g wholeATM_PG2.h5m -o atm_PG2_part.h5m
00007  *
00008  * the *PG2" style atm mesh is available in E3SM only for pg2 runs, something like
00009  *  --res ne30pg2_r05_oECv3_ICG --compset A_WCYCL1850S_CMIP6
00010  *  or
00011  *  --res ne4pg2_ne4pg2 --compset FC5AV1C-L
00012  */
00013 #include "moab/ProgOptions.hpp"
00014 #include "moab/Core.hpp"
00015 
00016 #include <cmath>
00017 #include <sstream>
00018 
00019 using namespace moab;
00020 
00021 int main( int argc, char* argv[] )
00022 {
00023 
00024     ProgOptions opts;
00025 
00026     std::string physfile, pg2file, outfile;
00027 
00028     opts.addOpt< std::string >( "physgridfile,p", "phys grid filename", &physfile );
00029     opts.addOpt< std::string >( "pg2file,g", "pg2 mesh file", &pg2file );
00030     opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
00031 
00032     opts.parseCommandLine( argc, argv );
00033 
00034     std::cout << "phys grid cloud file: " << physfile << "\n";
00035     std::cout << "pg2 mesh file: " << pg2file << "\n";
00036     std::cout << "output file: " << outfile << "\n";
00037 
00038     if( physfile.empty() )
00039     {
00040         opts.printHelp();
00041         return 0;
00042     }
00043     ErrorCode rval;
00044     Core* mb = new Core();
00045 
00046     rval = mb->load_file( physfile.c_str() );MB_CHK_SET_ERR( rval, "can't load phys grid file" );
00047 
00048     Core* mb2 = new Core();
00049     rval      = mb2->load_file( pg2file.c_str() );MB_CHK_SET_ERR( rval, "can't load pg2 mesh file" );
00050 
00051     Tag globalIDTag1 = mb->globalId_tag();
00052     Tag parti;
00053     rval = mb->tag_get_handle( "partition", parti );MB_CHK_SET_ERR( rval, "can't get partition tag phys grid mesh " );
00054 
00055     Tag globalIDTag2 = mb2->globalId_tag();
00056 
00057     Range verts1;
00058     rval = mb->get_entities_by_dimension( 0, 0, verts1 );MB_CHK_SET_ERR( rval, "can't get vertices " );
00059 
00060     std::vector< int > partValues;
00061     partValues.resize( verts1.size() );
00062     rval = mb->tag_get_data( parti, verts1, &partValues[0] );MB_CHK_SET_ERR( rval, "can't get parts values on vertices " );
00063 
00064     Range cells;
00065     rval = mb2->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "can't get 2d cells " );
00066     std::vector< int > globalIdsCells;
00067     globalIdsCells.resize( cells.size() );
00068     rval = mb2->tag_get_data( globalIDTag2, cells, &globalIdsCells[0] );MB_CHK_SET_ERR( rval, "can't get global ids cells " );
00069 
00070     std::vector< int > globalIdsVerts;
00071     globalIdsVerts.resize( verts1.size() );
00072     rval = mb->tag_get_data( globalIDTag1, verts1, &globalIdsVerts[0] );MB_CHK_SET_ERR( rval, "can't get global ids cells " );
00073 
00074     Tag partTag;
00075     rval = mb2->tag_get_handle( "PARALLEL_PARTITION", partTag );MB_CHK_SET_ERR( rval, "can't partition tag " );
00076 
00077     Range sets;
00078     rval = mb2->get_entities_by_type_and_tag( 0, MBENTITYSET, &partTag, NULL, 1, sets );MB_CHK_ERR( rval );
00079 
00080     std::vector< int > setValues;
00081     setValues.resize( sets.size() );
00082     rval = mb2->tag_get_data( partTag, sets, &setValues[0] );MB_CHK_ERR( rval );
00083 
00084     std::map< int, EntityHandle > valToSet;
00085     int i = 0;
00086     for( Range::iterator st = sets.begin(); st != sets.end(); ++st, i++ )
00087     {
00088         valToSet[setValues[i]] = *st;
00089     }
00090     // now, every cell will be put into one set, by looking at the global id of cell
00091 
00092     std::map< int, EntityHandle > gidToCell;
00093     i = 0;
00094     for( Range::iterator it = cells.begin(); it != cells.end(); ++it, i++ )
00095     {
00096         gidToCell[globalIdsCells[i]] = *it;
00097     }
00098     // empty all sets
00099     rval = mb2->clear_meshset( sets );MB_CHK_ERR( rval );
00100 
00101     // look now at parti values for vertices, and their global ids
00102     for( i = 0; i < (int)verts1.size(); i++ )
00103     {
00104         int part          = partValues[i];
00105         int gid           = globalIdsVerts[i];
00106         EntityHandle set1 = valToSet[part];
00107         EntityHandle cell = gidToCell[gid];
00108         rval              = mb2->add_entities( set1, &cell, 1 );MB_CHK_ERR( rval );
00109     }
00110 
00111     rval = mb2->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
00112 
00113     delete mb;
00114     delete mb2;
00115 
00116     return 0;
00117 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines