addConnec.cpp

Go to the documentation of this file.

    1 /*
    2  *
    3  *
    4 /** @example addConnec.cpp  Add connectivity to point cloud data, for better view in visit
    5  * this tool will take an existing h5m fine point cloud data (phys grid or land pc)
    6  *   and add 2d cells from a fine atm mesh
    7  *
    8  * example of usage:
    9  * ./addConnec -i wholeFineATM.h5m -s wholeLND_proj01.h5m -o wholeLndM.h5m
   10  *
   11  * Basically, will output a new h5m file (wholeLndM.h5m), which has also cell connectivity copied from fine atm
   12  *  matching is based on the global ids
   13  * between what we think is the order on the original file (wholeFineATM.h5m)
   14  *
   15  *  file  wholeFineATM.h5m is obtained from a coupled run in e3sm, with the ne 11, np 4,
   16  *
   17  */
   18 #include "moab/ProgOptions.hpp"
   19 #include "moab/Core.hpp"
   20  
   21 #include <cmath>
   22 #include <sstream>
   23  
   24 using namespace moab;
   25  
   26 int main( int argc, char* argv[] )
   27 {
   28  
   29     ProgOptions opts;
   30  
   31     std::string inputfile, outfile( "out.h5m" ), sourcefile;
   32  
   33     opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile );
   34     opts.addOpt< std::string >( "source,s", "h5m file aligned with the mesh input file", &sourcefile );
   35     opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
   36  
   37     opts.parseCommandLine( argc, argv );
   38  
   39     std::cout << "input mesh file: " << inputfile << "\n";
   40     std::cout << "source file: " << sourcefile << "\n";
   41     std::cout << "output file: " << outfile << "\n";
   42  
   43     if( inputfile.empty() )
   44     {
   45         opts.printHelp();
   46         return 0;
   47     }
   48     ErrorCode rval;
   49     Core* mb = new Core();
   50  
   51     rval = mb->load_file( inputfile.c_str() );MB_CHK_SET_ERR( rval, "can't load input file" );
   52  
   53     Core* mb2 = new Core();
   54     rval      = mb2->load_file( sourcefile.c_str() );MB_CHK_SET_ERR( rval, "can't load source file" );
   55  
   56     // get vertices on ini mesh; get global id on ini mesh
   57     // get global id on source mesh
   58     Tag gid;
   59     Tag gid2;
   60     rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get GLOBAL_ID tag on ini mesh " );
   61     rval = mb2->tag_get_handle( "GLOBAL_ID", gid2 );MB_CHK_SET_ERR( rval, "can't get GLOBAL_ID tag on source mesh " );
   62  
   63     // get vertices on ini mesh; build
   64     Range iniVerts;
   65     rval = mb->get_entities_by_dimension( 0, 0, iniVerts );MB_CHK_SET_ERR( rval, "can't get verts on initial mesh " );
   66  
   67     std::vector< int > gids;
   68     gids.resize( iniVerts.size() );
   69     rval = mb->tag_get_data( gid, iniVerts, &( gids[0] ) );MB_CHK_SET_ERR( rval, "can't get gid on initial verts " );
   70     // build now the map
   71     std::map< int, EntityHandle > fromGidToEh;
   72     int i = 0;
   73     for( Range::iterator vit = iniVerts.begin(); vit != iniVerts.end(); ++vit, i++ )
   74     {
   75         fromGidToEh[gids[i]] = *vit;
   76     }
   77     // now get the source verts, and tags, and set it on new mesh
   78  
   79     Range sourceVerts;
   80     rval = mb2->get_entities_by_dimension( 0, 0, sourceVerts );MB_CHK_SET_ERR( rval, "can't get verts on source mesh " );
   81     std::vector< int > gids2;
   82     gids2.resize( sourceVerts.size() );
   83     rval = mb2->tag_get_data( gid2, sourceVerts, &( gids2[0] ) );MB_CHK_SET_ERR( rval, "can't get gid2 on cloud mesh " );
   84     std::map< int, EntityHandle > fromGid2ToEh;
   85     i = 0;
   86     for( Range::iterator vit = sourceVerts.begin(); vit != sourceVerts.end(); ++vit, i++ )
   87     {
   88         fromGid2ToEh[gids2[i]] = *vit;
   89     }
   90  
   91     Range usedVerts;
   92     for( size_t i = 0; i < gids2.size(); i++ )
   93     {
   94         int globalId          = gids2[i];
   95         EntityHandle usedVert = fromGidToEh[globalId];
   96         usedVerts.insert( usedVert );
   97     }
   98     Range cells;
   99     rval = mb->get_adjacencies( usedVerts, 2, false, cells, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adj cells " );
  100  
  101     // now create a new cell in mb2 for each one in mb
  102     for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
  103     {
  104         EntityHandle cell        = *cit;
  105         int nnodes               = 0;
  106         const EntityHandle* conn = NULL;
  107         rval                     = mb->get_connectivity( cell, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity " );
  108  
  109         std::vector< EntityHandle > nconn( nnodes );
  110         bool goodCell = true;
  111         for( int i = 0; i < nnodes; i++ )
  112         {
  113             EntityHandle v = conn[i];
  114             int index      = iniVerts.index( v );
  115             if( index < 0 ) goodCell = false;
  116             int id = gids[index];
  117             if( fromGid2ToEh.find( id ) == fromGid2ToEh.end() )
  118                 goodCell = false;
  119             else
  120                 nconn[i] = fromGid2ToEh[id];
  121         }
  122         if( !goodCell ) continue;
  123         EntityType type = MBTRI;
  124         if( nnodes == 3 )
  125             type = MBTRI;
  126         else if( nnodes == 4 )
  127             type = MBQUAD;
  128         else if( nnodes > 4 )
  129             type = MBPOLYGON;
  130         EntityHandle nel;
  131         rval = mb2->create_element( type, &nconn[0], nnodes, nel );MB_CHK_SET_ERR( rval, "can't create new cell " );
  132     }
  133     // save file
  134     rval = mb2->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
  135  
  136     delete mb;
  137     delete mb2;
  138  
  139     return 0;
  140 }