|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
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 <cmath> 00026 #include <sstream> 00027 #include <map> 00028 #include <Eigen/Sparse> 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 }
1.7.6.1 
