MOAB: Mesh Oriented datABase  (version 5.4.1)
addConnec.cpp
Go to the documentation of this file.
00001 /*
00002  *
00003  *
00004 /** @example addConnec.cpp  Add connectivity to point cloud data, for better view in visit
00005  * this tool will take an existing h5m fine point cloud data (phys grid or land pc)
00006  *   and add 2d cells from a fine atm mesh
00007  *
00008  * example of usage:
00009  * ./addConnec -i wholeFineATM.h5m -s wholeLND_proj01.h5m -o wholeLndM.h5m
00010  *
00011  * Basically, will output a new h5m file (wholeLndM.h5m), which has also cell connectivity copied from fine atm
00012  *  matching is based on the global ids
00013  * between what we think is the order on the original file (wholeFineATM.h5m)
00014  *
00015  *  file  wholeFineATM.h5m is obtained from a coupled run in e3sm, with the ne 11, np 4,
00016  *
00017  */
00018 #include "moab/ProgOptions.hpp"
00019 #include "moab/Core.hpp"
00020 
00021 #include <cmath>
00022 #include <sstream>
00023 
00024 using namespace moab;
00025 
00026 int main( int argc, char* argv[] )
00027 {
00028 
00029     ProgOptions opts;
00030 
00031     std::string inputfile, outfile( "out.h5m" ), sourcefile;
00032 
00033     opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile );
00034     opts.addOpt< std::string >( "source,s", "h5m file aligned with the mesh input file", &sourcefile );
00035     opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
00036 
00037     opts.parseCommandLine( argc, argv );
00038 
00039     std::cout << "input mesh file: " << inputfile << "\n";
00040     std::cout << "source file: " << sourcefile << "\n";
00041     std::cout << "output file: " << outfile << "\n";
00042 
00043     if( inputfile.empty() )
00044     {
00045         opts.printHelp();
00046         return 0;
00047     }
00048     ErrorCode rval;
00049     Core* mb = new Core();
00050 
00051     rval = mb->load_file( inputfile.c_str() );MB_CHK_SET_ERR( rval, "can't load input file" );
00052 
00053     Core* mb2 = new Core();
00054     rval      = mb2->load_file( sourcefile.c_str() );MB_CHK_SET_ERR( rval, "can't load source file" );
00055 
00056     // get vertices on ini mesh; get global id on ini mesh
00057     // get global id on source mesh
00058     Tag gid;
00059     Tag gid2;
00060     rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get GLOBAL_ID tag on ini mesh " );
00061     rval = mb2->tag_get_handle( "GLOBAL_ID", gid2 );MB_CHK_SET_ERR( rval, "can't get GLOBAL_ID tag on source mesh " );
00062 
00063     // get vertices on ini mesh; build
00064     Range iniVerts;
00065     rval = mb->get_entities_by_dimension( 0, 0, iniVerts );MB_CHK_SET_ERR( rval, "can't get verts on initial mesh " );
00066 
00067     std::vector< int > gids;
00068     gids.resize( iniVerts.size() );
00069     rval = mb->tag_get_data( gid, iniVerts, &( gids[0] ) );MB_CHK_SET_ERR( rval, "can't get gid on initial verts " );
00070     // build now the map
00071     std::map< int, EntityHandle > fromGidToEh;
00072     int i = 0;
00073     for( Range::iterator vit = iniVerts.begin(); vit != iniVerts.end(); ++vit, i++ )
00074     {
00075         fromGidToEh[gids[i]] = *vit;
00076     }
00077     // now get the source verts, and tags, and set it on new mesh
00078 
00079     Range sourceVerts;
00080     rval = mb2->get_entities_by_dimension( 0, 0, sourceVerts );MB_CHK_SET_ERR( rval, "can't get verts on source mesh " );
00081     std::vector< int > gids2;
00082     gids2.resize( sourceVerts.size() );
00083     rval = mb2->tag_get_data( gid2, sourceVerts, &( gids2[0] ) );MB_CHK_SET_ERR( rval, "can't get gid2 on cloud mesh " );
00084     std::map< int, EntityHandle > fromGid2ToEh;
00085     i = 0;
00086     for( Range::iterator vit = sourceVerts.begin(); vit != sourceVerts.end(); ++vit, i++ )
00087     {
00088         fromGid2ToEh[gids2[i]] = *vit;
00089     }
00090 
00091     Range usedVerts;
00092     for( size_t i = 0; i < gids2.size(); i++ )
00093     {
00094         int globalId          = gids2[i];
00095         EntityHandle usedVert = fromGidToEh[globalId];
00096         usedVerts.insert( usedVert );
00097     }
00098     Range cells;
00099     rval = mb->get_adjacencies( usedVerts, 2, false, cells, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adj cells " );
00100 
00101     // now create a new cell in mb2 for each one in mb
00102     for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
00103     {
00104         EntityHandle cell        = *cit;
00105         int nnodes               = 0;
00106         const EntityHandle* conn = NULL;
00107         rval                     = mb->get_connectivity( cell, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity " );
00108 
00109         std::vector< EntityHandle > nconn( nnodes );
00110         bool goodCell = true;
00111         for( int i = 0; i < nnodes; i++ )
00112         {
00113             EntityHandle v = conn[i];
00114             int index      = iniVerts.index( v );
00115             if( index < 0 ) goodCell = false;
00116             int id = gids[index];
00117             if( fromGid2ToEh.find( id ) == fromGid2ToEh.end() )
00118                 goodCell = false;
00119             else
00120                 nconn[i] = fromGid2ToEh[id];
00121         }
00122         if( !goodCell ) continue;
00123         EntityType type = MBTRI;
00124         if( nnodes == 3 )
00125             type = MBTRI;
00126         else if( nnodes == 4 )
00127             type = MBQUAD;
00128         else if( nnodes > 4 )
00129             type = MBPOLYGON;
00130         EntityHandle nel;
00131         rval = mb2->create_element( type, &nconn[0], nnodes, nel );MB_CHK_SET_ERR( rval, "can't create new cell " );
00132     }
00133     // save file
00134     rval = mb2->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
00135 
00136     delete mb;
00137     delete mb2;
00138 
00139     return 0;
00140 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines