fathom/moab-docs-develop/visuMapVtk_8cpp_source.html

00001 /*
00002  * visuMapVtk.cpp
00003  * this tool will take a source file, target file (h5m) and a map file in nc format,
00004  *
00005  * example of usage:
00006  * ./mbvisumap -s source.h5m -t target.h5m -m map.nc -b startSourceID -e endSourceID  -c startTargetID -f endTargetID
00007  * will associate row i in map with a partial mesh
00008  *  will associate row j in map with a partial mesh
00009  *
00010  * can be built only if netcdf and hdf5 and eigen3 are available
00011  *
00012  *
00013  */
00014 #include "moab/MOABConfig.h"
00015
00016 #ifndef MOAB_HAVE_EIGEN3
00017 #error compareMaps tool requires eigen3 configuration
00018 #endif
00019
00020 #include "moab/ProgOptions.hpp"
00021 #include "moab/Core.hpp"
00022 #include "moab/Range.hpp"
00023
00024 #include "netcdf.h"
00025 #include
00026 #include
00027 #include
00028 #include
00029
00030 #define ERR_NC( e )                                \
00031     {                                              \
00032         printf( "Error: %s\n", nc_strerror( e ) ); \
00033         exit( 2 );                                 \
00034     }
00035
00036 // copy from ReadNCDF.cpp some useful macros for reading from a netcdf file (exodus?)
00037 // ncFile is an integer initialized when opening the nc file in read mode
00038
00039 int ncFile1;
00040
00041 #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
00042
00043 #define GET_DIM1( ncdim, name, val )                            \
00044     {                                                           \
00045         int gdfail = nc_inq_dimid( ncFile1, name, &( ncdim ) ); \
00046         if( NC_NOERR == gdfail )                                \
00047         {                                                       \
00048             size_t tmp_val;                                     \
00049             gdfail = nc_inq_dimlen( ncFile1, ncdim, &tmp_val ); \
00050             if( NC_NOERR != gdfail )                            \
00051             {                                                   \
00052                 ERR_NC( gdfail )                                \
00053             }                                                   \
00054             else                                                \
00055                 ( val ) = tmp_val;                              \
00056         }                                                       \
00057         else                                                    \
00058             ( val ) = 0;                                        \
00059     }
00060
00061 #define GET_DIMB1( ncdim, name, varname, id, val ) \
00062     INS_ID( name, varname, id );                   \
00063     GET_DIM1( ncdim, name, val );
00064
00065 #define GET_VAR1( name, id, dims )                                     \
00066     {                                                                  \
00067         ( id )     = -1;                                               \
00068         int gvfail = nc_inq_varid( ncFile1, name, &( id ) );           \
00069         if( NC_NOERR == gvfail )                                       \
00070         {                                                              \
00071             int ndims;                                                 \
00072             gvfail = nc_inq_varndims( ncFile1, id, &ndims );           \
00073             if( NC_NOERR == gvfail )                                   \
00074             {                                                          \
00075                 ( dims ).resize( ndims );                              \
00076                 gvfail = nc_inq_vardimid( ncFile1, id, &( dims )[0] ); \
00077                 if( NC_NOERR != gvfail )                               \
00078                 {                                                      \
00079                     ERR_NC( gvfail )                                   \
00080                 }                                                      \
00081             }                                                          \
00082         }                                                              \
00083     }
00084
00085 #define GET_1D_INT_VAR1( name, id, vals )                                               \
00086     {                                                                                   \
00087         GET_VAR1( name, id, vals );                                                     \
00088         if( -1 != ( id ) )                                                              \
00089         {                                                                               \
00090             size_t ntmp;                                                                \
00091             int ivfail = nc_inq_dimlen( ncFile1, ( vals )[0], &ntmp );                  \
00092             if( NC_NOERR != ivfail )                                                    \
00093             {                                                                           \
00094                 ERR_NC( ivfail )                                                        \
00095             }                                                                           \
00096             ( vals ).resize( ntmp );                                                    \
00097             size_t ntmp1 = 0;                                                           \
00098             ivfail       = nc_get_vara_int( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
00099             if( NC_NOERR != ivfail )                                                    \
00100             {                                                                           \
00101                 ERR_NC( ivfail )                                                        \
00102             }                                                                           \
00103         }                                                                               \
00104     }
00105
00106 #define GET_1D_DBL_VAR1( name, id, vals )                                                  \
00107     {                                                                                      \
00108         std::vector< int > dum_dims;                                                       \
00109         GET_VAR1( name, id, dum_dims );                                                    \
00110         if( -1 != ( id ) )                                                                 \
00111         {                                                                                  \
00112             size_t ntmp;                                                                   \
00113             int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp );                     \
00114             if( NC_NOERR != dvfail )                                                       \
00115             {                                                                              \
00116                 ERR_NC( dvfail )                                                           \
00117             }                                                                              \
00118             ( vals ).resize( ntmp );                                                       \
00119             size_t ntmp1 = 0;                                                              \
00120             dvfail       = nc_get_vara_double( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
00121             if( NC_NOERR != dvfail )                                                       \
00122             {                                                                              \
00123                 ERR_NC( dvfail )                                                           \
00124             }                                                                              \
00125         }                                                                                  \
00126     }
00127
00128 using namespace moab;
00129 using namespace std;
00130
00131 int main( int argc, char* argv[] )
00132 {
00133
00134     ProgOptions opts;
00135     int dimSource = 2;  // for FV meshes is 2; for SE meshes, use fine mesh, dim will be 0
00136     int dimTarget = 2;  //
00137     int otype     = 0;
00138
00139     std::string inputfile1, inputSource, inputTarget;
00140     opts.addOpt< std::string >( "map,m", "input map ", &inputfile1 );
00141     opts.addOpt< std::string >( "source,s", "source mesh", &inputSource );
00142     opts.addOpt< std::string >( "target,t", "target mesh", &inputTarget );
00143     opts.addOpt< int >( "dimSource,d", "dimension of source  ", &dimSource );
00144     opts.addOpt< int >( "dimTarget,g", "dimension of target  ", &dimTarget );
00145     opts.addOpt< int >( "typeOutput,o", " output type vtk(0), h5m(1)", &otype );
00146
00147     int startSourceID = -1, endSourceID = -1, startTargetID = -1, endTargetID = -1;
00148     opts.addOpt< int >( "startSourceID,b", "start source id ", &startSourceID );
00149     opts.addOpt< int >( "endSourceID,e", "end source id ", &endSourceID );
00150     opts.addOpt< int >( "startTargetID,c", "start target id ", &startTargetID );
00151     opts.addOpt< int >( "endTargetID,f", "end target id ", &endTargetID );
00152     //  -b startSourceID -e endSourceID  -c startTargetID -f endTargetID
00153
00154     opts.parseCommandLine( argc, argv );
00155
00156     std::string extension = ".vtk";
00157     if( 1 == otype ) extension = ".h5m";
00158     // Open netcdf/exodus file
00159     int fail = nc_open( inputfile1.c_str(), 0, &ncFile1 );
00160     if( NC_NOWRITE != fail )
00161     {
00162         ERR_NC( fail )
00163     }
00164
00165     std::cout << " opened " << inputfile1 << " for map 1 \n";
00166
00167     int temp_dim;
00168     int na1, nb1, ns1;
00169     GET_DIM1( temp_dim, "n_a", na1 );
00170     GET_DIM1( temp_dim, "n_b", nb1 );
00171     GET_DIM1( temp_dim, "n_s", ns1 );
00172     std::cout << " n_a, n_b, n_s : " << na1 << ", " << nb1 << ", " << ns1 << " for map 1 \n";
00173     std::vector< int > col1( ns1 ), row1( ns1 );
00174     std::vector< double > val1( ns1 );
00175     int idrow1, idcol1, ids1;
00176     GET_1D_INT_VAR1( "row", idrow1, row1 );
00177     GET_1D_INT_VAR1( "col", idcol1, col1 );
00178     GET_1D_DBL_VAR1( "S", ids1, val1 );
00179     // first matrix
00180     typedef Eigen::Triplet< double > Triplet;
00181     std::vector< Triplet > tripletList;
00182     tripletList.reserve( ns1 );
00183     for( int iv = 0; iv < ns1; iv++ )
00184     {
00185         // all row and column indices are 1-based in the map file.
00186         // inside Eigen3, we will use 0-based; then we will have to add back 1 when
00187         //     we dump out the files
00188         tripletList.push_back( Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
00189     }
00190     Eigen::SparseMatrix< double > weight1( nb1, na1 );
00191
00192     weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
00193     weight1.makeCompressed();
00194     // we read the matrix; now read moab source and target
00195     Core core;
00196     Interface* mb = &core;
00197     ErrorCode rval;
00198     Tag gtag = mb->globalId_tag();
00199     EntityHandle sourceSet, targetSet;
00200     // those are the maps from global ids to the moab entity handles corresponding to those global ids
00201     //   which are corresponding to the global DOFs
00202     map< int, EntityHandle > sourceHandles;
00203     map< int, EntityHandle > targetHandles;
00204     rval = mb->create_meshset( MESHSET_SET, sourceSet );MB_CHK_SET_ERR( rval, "can't create source mesh set" );
00205     rval = mb->create_meshset( MESHSET_SET, targetSet );MB_CHK_SET_ERR( rval, "can't create target mesh set" );
00206     const char* readopts = "";
00207     rval                 = mb->load_file( inputSource.c_str(), &sourceSet, readopts );MB_CHK_SET_ERR( rval, "Failed to read" );
00208     rval = mb->load_file( inputTarget.c_str(), &targetSet, readopts );MB_CHK_SET_ERR( rval, "Failed to read" );
00209     Range sRange;
00210     rval = mb->get_entities_by_dimension( sourceSet, dimSource, sRange );MB_CHK_SET_ERR( rval, "Failed to get sRange" );
00211     vector< int > sids;
00212     sids.resize( sRange.size() );
00213     rval = mb->tag_get_data( gtag, sRange, &sids[0] );MB_CHK_SET_ERR( rval, "Failed to get ids for srange" );
00214     // all global ids are 1 based in the source file, and they correspond to the dofs in the map file
00215     for( size_t i = 0; i < sids.size(); i++ )
00216     {
00217         EntityHandle eh    = sRange[i];
00218         int gid            = sids[i];
00219         sourceHandles[gid] = eh;
00220     }
00221     Range tRange;
00222     rval = mb->get_entities_by_dimension( targetSet, dimTarget, tRange );MB_CHK_SET_ERR( rval, "Failed to get tRange" );
00223     vector< int > tids;
00224     tids.resize( tRange.size() );
00225     rval = mb->tag_get_data( gtag, tRange, &tids[0] );MB_CHK_SET_ERR( rval, "Failed to get ids for trange" );
00226     // all global ids are 1 based in the target file, and they correspond to the dofs in the map file
00227     for( size_t i = 0; i < tids.size(); i++ )
00228     {
00229         EntityHandle eh    = tRange[i];
00230         int gid            = tids[i];
00231         targetHandles[gid] = eh;
00232     }
00233     // a dense tag for weights
00234     Tag wtag;
00235     double defVal = 0;
00236
00237     std::string name_map = inputfile1;
00238     // strip last 3 chars (.nc extension)
00239     name_map.erase( name_map.begin() + name_map.length() - 3, name_map.end() );
00240     // if path , remove from name
00241     size_t pos = name_map.rfind( '/', name_map.length() );
00242     if( pos != std::string::npos ) name_map = name_map.erase( 0, pos + 1 );
00243
00244     rval = mb->tag_get_handle( "weight", 1, MB_TYPE_DOUBLE, wtag, MB_TAG_CREAT | MB_TAG_DENSE, &defVal );MB_CHK_SET_ERR( rval, "Failed to create weight" );
00245     EntityHandle partialSet;
00246     rval = mb->create_meshset( MESHSET_SET, partialSet );MB_CHK_SET_ERR( rval, "can't create partial set" );
00247     // how to get a complete row in sparse matrix? Or a complete column ?
00248     for( int col = startSourceID - 1; col <= endSourceID - 1; col++ )
00249     {
00250         Range targetEnts;  // will find entries for column col-1 in sparse matrix weight1, and its entries
00251         // will assign a dense tag with values, and write out the file
00252         if( col < 0 ) continue;
00253         Eigen::SparseVector< double > colVect = weight1.col( col );
00254         // the row indices correspond to target cells
00255         for( Eigen::SparseVector< double >::InnerIterator it( colVect ); it; ++it )
00256         {
00257             double weight = it.value();  // == vec[ it.index() ]
00258                                          // we add back the 1 that we subtract
00259             int globalIdRow = it.index() + 1;
00260             EntityHandle th = targetHandles[globalIdRow];
00261             targetEnts.insert( th );
00262             rval = mb->tag_set_data( wtag, &th, 1, &weight );MB_CHK_SET_ERR( rval, "Failed to set weight tag on target" );
00263         }
00264
00265         if( dimTarget == 0 )
00266         {
00267             Range adjCells;
00268             rval = mb->get_adjacencies( targetEnts, 2, false, adjCells, Interface::UNION );MB_CHK_SET_ERR( rval, " can't get adj cells " );
00269             targetEnts.merge( adjCells );
00270         }
00271
00272         rval = mb->add_entities( partialSet, targetEnts );MB_CHK_SET_ERR( rval, "Failed to add target entities to partial set" );
00273         // write now the set in a numbered file
00274         std::stringstream fff;
00275         fff << name_map << "_column" << col + 1 << extension;
00276         rval = mb->write_mesh( fff.str().c_str(), &partialSet, 1 );MB_CHK_ERR( rval );
00277         // remove from partial set the entities it has
00278         rval = mb->clear_meshset( &partialSet, 1 );MB_CHK_SET_ERR( rval, "Failed to empty partial set" );
00279     }
00280
00281     // how to get a complete row in sparse matrix?
00282     for( int row = startTargetID - 1; row <= endTargetID - 1; row++ )
00283     {
00284         Range sourceEnts;  // will find entries for row in sparse matrix weight1, and its entries
00285         // will assign a dense tag with values, and write out the file
00286         if( row < 0 ) continue;
00287         Eigen::SparseVector< double > rowVect = weight1.row( row );
00288         // the row indices correspond to target cells
00289         for( Eigen::SparseVector< double >::InnerIterator it( rowVect ); it; ++it )
00290         {
00291             double weight   = it.value();  // == vec[ it.index() ]
00292             int globalIdCol = it.index() + 1;
00293             EntityHandle sh = sourceHandles[globalIdCol];
00294             sourceEnts.insert( sh );
00295             rval = mb->tag_set_data( wtag, &sh, 1, &weight );MB_CHK_SET_ERR( rval, "Failed to set weight tag on source" );
00296         }
00297         if( dimSource == 0 )
00298         {
00299             Range adjCells;
00300             rval = mb->get_adjacencies( sourceEnts, 2, false, adjCells, Interface::UNION );MB_CHK_SET_ERR( rval, " can't get adj cells " );
00301             sourceEnts.merge( adjCells );
00302         }
00303         rval = mb->add_entities( partialSet, sourceEnts );MB_CHK_SET_ERR( rval, "Failed to add source entities" );
00304         // write now the set in a numbered file
00305         std::stringstream fff;
00306         fff << name_map << "_row" << row + 1 << extension;
00307         rval = mb->write_mesh( fff.str().c_str(), &partialSet, 1 );MB_CHK_ERR( rval );
00308         rval = mb->clear_meshset( &partialSet, 1 );MB_CHK_SET_ERR( rval, "Failed to empty partial set" );
00309     }
00310     return 0;
00311 }