 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
00001 // Usage:
00002 // tools/h5mtoscrip -w map_atm_to_ocn.h5m -s map_atm_to_ocn.nc --coords
00003 //
00004 #include
00005 #include
00006 #include
00007 #include
00008 #include
00009 #include
00010 #include
00011 #include
00012
00013 #include "moab/ProgOptions.hpp"
00014 #include "moab/Core.hpp"
00015
00016 #ifdef MOAB_HAVE_MPI
00017 #include "moab_mpi.h"
00018 #endif
00019
00020 #ifndef MOAB_HAVE_TEMPESTREMAP
00021 #error Tool requires compilation with TempestRemap dependency
00022 #endif
00023
00024 // TempestRemap includes
00025 #include "OfflineMap.h"
00026 #include "netcdfcpp.h"
00027 #include "NetCDFUtilities.h"
00028 #include "DataArray2D.h"
00029
00030 using namespace moab;
00031
00032 template < typename T >
00033 ErrorCode get_vartag_data( moab::Interface* mbCore,
00034 Tag tag,
00035 moab::Range& sets,
00036 int& out_data_size,
00037 std::vector< T >& data )
00038 {
00039 int* tag_sizes = new int[sets.size()];
00040 const void** tag_data = (const void**)new void*[sets.size()];
00041
00042 ErrorCode rval = mbCore->tag_get_by_ptr( tag, sets, tag_data, tag_sizes );MB_CHK_SET_ERR( rval, "Getting matrix rows failed" );
00043
00044 out_data_size = 0;
00045 for( unsigned is = 0; is < sets.size(); ++is )
00046 out_data_size += tag_sizes[is];
00047
00048 data.resize( out_data_size );
00049 int ioffset = 0;
00050 for( unsigned index = 0; index < sets.size(); index++ )
00051 {
00052 T* m_vals = (T*)tag_data[index];
00053 for( int k = 0; k < tag_sizes[index]; k++ )
00054 {
00055 data[ioffset++] = m_vals[k];
00056 }
00057 }
00058
00059 return moab::MB_SUCCESS;
00060 }
00061
00062 void ReadFileMetaData( std::string& metaFilename, std::map< std::string, std::string >& metadataVals )
00063 {
00064 std::ifstream metafile;
00065 std::string line;
00066
00067 metafile.open( metaFilename.c_str() );
00068 metadataVals["Title"] = "MOAB-TempestRemap (MBTR) Offline Regridding Weight Converter (h5mtoscrip)";
00069 std::string key, value;
00070 while( std::getline( metafile, line ) )
00071 {
00072 size_t lastindex = line.find_last_of( "=" );
00073 key = line.substr( 0, lastindex - 1 );
00074 value = line.substr( lastindex + 2, line.length() );
00075
00076 metadataVals[std::string( key )] = std::string( value );
00077 }
00078 metafile.close();
00079 }
00080
00081 int main( int argc, char* argv[] )
00082 {
00083 moab::ErrorCode rval;
00084 int dimension = 2;
00085 NcError error2( NcError::verbose_nonfatal );
00086 std::stringstream sstr;
00087 ProgOptions opts;
00088 std::string h5mfilename, scripfile;
00089 bool noMap = false;
00090 bool writeXYCoords = false;
00091
00092 #ifdef MOAB_HAVE_MPI
00093 MPI_Init( &argc, &argv );
00094 #endif
00095
00096 opts.addOpt< std::string >( "weights,w", "h5m remapping weights filename", &h5mfilename );
00097 opts.addOpt< std::string >( "scrip,s", "Output SCRIP map filename", &scripfile );
00098 opts.addOpt< int >( "dim,d", "Dimension of entities to use for partitioning", &dimension );
00099 opts.addOpt< void >( "mesh,m", "Only convert the mesh and exclude the remap weight details", &noMap );
00100 opts.addOpt< void >( "coords,c", "Write the center and vertex coordinates in lat/lon format", &writeXYCoords );
00101
00102 opts.parseCommandLine( argc, argv );
00103
00104 if( h5mfilename.empty() || scripfile.empty() )
00105 {
00106 opts.printHelp();
00107 exit( 1 );
00108 }
00109
00110 moab::Interface* mbCore = new( std::nothrow ) moab::Core;
00111
00112 if( NULL == mbCore )
00113 {
00114 return 1;
00115 }
00116
00117 // Set the read options for parallel file loading
00118 const std::string partition_set_name = "PARALLEL_PARTITION";
00119 const std::string global_id_name = "GLOBAL_ID";
00120
00121 // Load file
00122 rval = mbCore->load_mesh( h5mfilename.c_str() );MB_CHK_ERR( rval );
00123
00124 try
00125 {
00126
00127 // Temporarily change rval reporting
00128 NcError error_temp( NcError::verbose_fatal );
00129
00130 // Open an output file
00131 NcFile ncMap( scripfile.c_str(), NcFile::Replace, NULL, 0, NcFile::Offset64Bits );
00132 if( !ncMap.is_valid() )
00133 {
00134 _EXCEPTION1( "Unable to open output map file \"%s\"", scripfile.c_str() );
00135 }
00136
00137 {
00138 // NetCDF-SCRIP Global Attributes
00139 std::map< std::string, std::string > mapAttributes;
00140 size_t lastindex = h5mfilename.find_last_of( "." );
00141 std::stringstream sstr;
00142 sstr << h5mfilename.substr( 0, lastindex ) << ".meta";
00143 std::string metaFilename = sstr.str();
00144 ReadFileMetaData( metaFilename, mapAttributes );
00145 mapAttributes["Command"] =
00146 "Converted with MOAB:h5mtoscrip with --w=" + h5mfilename + " and --s=" + scripfile;
00147
00148 // Add global attributes
00149 std::map< std::string, std::string >::const_iterator iterAttributes = mapAttributes.begin();
00150 for( ; iterAttributes != mapAttributes.end(); iterAttributes++ )
00151 {
00152
00153 std::cout << iterAttributes->first << " -- " << iterAttributes->second << std::endl;
00154 ncMap.add_att( iterAttributes->first.c_str(), iterAttributes->second.c_str() );
00155 }
00156 std::cout << "\n";
00157 }
00158
00159 Tag globalIDTag, materialSetTag;
00160 globalIDTag = mbCore->globalId_tag();
00161 // materialSetTag = mbCore->material_tag();
00162 rval = mbCore->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, materialSetTag, MB_TAG_SPARSE );MB_CHK_ERR( rval );
00163
00164 // Get sets entities, by type
00165 moab::Range meshsets;
00166 rval = mbCore->get_entities_by_type_and_tag( 0, MBENTITYSET, &globalIDTag, NULL, 1, meshsets,
00167 moab::Interface::UNION, true );MB_CHK_ERR( rval );
00168
00169 moab::EntityHandle rootset = 0;
00170 ///////////////////////////////////////////////////////////////////////////
00171 // The metadata in H5M file contains the following data:
00172 //
00173 // 1. n_a: Total source entities: (number of elements in source mesh)
00174 // 2. n_b: Total target entities: (number of elements in target mesh)
00175 // 3. nv_a: Max edge size of elements in source mesh
00176 // 4. nv_b: Max edge size of elements in target mesh
00177 // 5. maxrows: Number of rows in remap weight matrix
00178 // 6. maxcols: Number of cols in remap weight matrix
00179 // 7. nnz: Number of total nnz in sparse remap weight matrix
00180 // 8. np_a: The order of the field description on the source mesh: >= 1
00181 // 9. np_b: The order of the field description on the target mesh: >= 1
00182 // 10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2
00183 // = dGLL]
00184 // 11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2
00185 // = dGLL]
00186 // 12. conserved: Flag to specify whether the remap operator has conservation constraints:
00187 // [0, 1]
00188 // 13. monotonicity: Flags to specify whether the remap operator has monotonicity
00189 // constraints: [0, 1, 2]
00190 //
00191 ///////////////////////////////////////////////////////////////////////////
00192 Tag smatMetadataTag;
00193 int smat_metadata_glb[13];
00194 rval = mbCore->tag_get_handle( "SMAT_DATA", 13, MB_TYPE_INTEGER, smatMetadataTag, MB_TAG_SPARSE );MB_CHK_ERR( rval );
00195 rval = mbCore->tag_get_data( smatMetadataTag, &rootset, 1, smat_metadata_glb );MB_CHK_ERR( rval );
00196 // std::cout << "Number of mesh sets is " << meshsets.size() << std::endl;
00197
00198 #define DTYPE( a ) \
00199 { \
00200 ( ( ( a ) == 0 ) ? "FV" : ( ( ( a ) == 1 ) ? "cGLL" : "dGLL" ) ) \
00201 }
00202 // Map dimensions
00203 int nA = smat_metadata_glb[0];
00204 int nB = smat_metadata_glb[1];
00205 int nVA = smat_metadata_glb[2];
00206 int nVB = smat_metadata_glb[3];
00207 int nDofB = smat_metadata_glb[4];
00208 int nDofA = smat_metadata_glb[5];
00209 int NNZ = smat_metadata_glb[6];
00210 int nOrdA = smat_metadata_glb[7];
00211 int nOrdB = smat_metadata_glb[8];
00212 int nBasA = smat_metadata_glb[9];
00213 std::string methodA = DTYPE( nBasA );
00214 int nBasB = smat_metadata_glb[10];
00215 std::string methodB = DTYPE( nBasB );
00216 int bConserved = smat_metadata_glb[11];
00217 int bMonotonicity = smat_metadata_glb[12];
00218
00219 EntityHandle source_mesh = 0, target_mesh = 0, overlap_mesh = 0;
00220 for( unsigned im = 0; im < meshsets.size(); ++im )
00221 {
00222 moab::Range elems;
00223 rval = mbCore->get_entities_by_dimension( meshsets[im], 2, elems );MB_CHK_ERR( rval );
00224 if( elems.size() - nA == 0 && source_mesh == 0 )
00225 source_mesh = meshsets[im];
00226 else if( elems.size() - nB == 0 && target_mesh == 0 )
00227 target_mesh = meshsets[im];
00228 else if( overlap_mesh == 0 )
00229 overlap_mesh = meshsets[im];
00230 else
00231 continue;
00232 }
00233
00234 Tag srcIDTag, srcAreaTag, tgtIDTag, tgtAreaTag;
00235 rval = mbCore->tag_get_handle( "SourceGIDS", srcIDTag );MB_CHK_ERR( rval );
00236 rval = mbCore->tag_get_handle( "SourceAreas", srcAreaTag );MB_CHK_ERR( rval );
00237 rval = mbCore->tag_get_handle( "TargetGIDS", tgtIDTag );MB_CHK_ERR( rval );
00238 rval = mbCore->tag_get_handle( "TargetAreas", tgtAreaTag );MB_CHK_ERR( rval );
00239 Tag smatRowdataTag, smatColdataTag, smatValsdataTag;
00240 rval = mbCore->tag_get_handle( "SMAT_ROWS", smatRowdataTag );MB_CHK_ERR( rval );
00241 rval = mbCore->tag_get_handle( "SMAT_COLS", smatColdataTag );MB_CHK_ERR( rval );
00242 rval = mbCore->tag_get_handle( "SMAT_VALS", smatValsdataTag );MB_CHK_ERR( rval );
00243 Tag srcCenterLon, srcCenterLat, tgtCenterLon, tgtCenterLat;
00244 rval = mbCore->tag_get_handle( "SourceCoordCenterLon", srcCenterLon );MB_CHK_ERR( rval );
00245 rval = mbCore->tag_get_handle( "SourceCoordCenterLat", srcCenterLat );MB_CHK_ERR( rval );
00246 rval = mbCore->tag_get_handle( "TargetCoordCenterLon", tgtCenterLon );MB_CHK_ERR( rval );
00247 rval = mbCore->tag_get_handle( "TargetCoordCenterLat", tgtCenterLat );MB_CHK_ERR( rval );
00248 Tag srcVertexLon, srcVertexLat, tgtVertexLon, tgtVertexLat;
00249 rval = mbCore->tag_get_handle( "SourceCoordVertexLon", srcVertexLon );MB_CHK_ERR( rval );
00250 rval = mbCore->tag_get_handle( "SourceCoordVertexLat", srcVertexLat );MB_CHK_ERR( rval );
00251 rval = mbCore->tag_get_handle( "TargetCoordVertexLon", tgtVertexLon );MB_CHK_ERR( rval );
00252 rval = mbCore->tag_get_handle( "TargetCoordVertexLat", tgtVertexLat );MB_CHK_ERR( rval );
00253
00254 // Get sets entities, by type
00255 moab::Range sets;
00256 // rval = mbCore->get_entities_by_type(0, MBENTITYSET, sets);MB_CHK_ERR(rval);
00257 rval = mbCore->get_entities_by_type_and_tag( 0, MBENTITYSET, &smatRowdataTag, NULL, 1, sets,
00258 moab::Interface::UNION, true );MB_CHK_ERR( rval );
00259
00260 std::vector< int > src_gids, tgt_gids;
00261 std::vector< double > src_areas, tgt_areas;
00262 int srcID_size, tgtID_size, srcArea_size, tgtArea_size;
00263 rval = get_vartag_data( mbCore, srcIDTag, sets, srcID_size, src_gids );MB_CHK_SET_ERR( rval, "Getting source mesh IDs failed" );
00264 rval = get_vartag_data( mbCore, tgtIDTag, sets, tgtID_size, tgt_gids );MB_CHK_SET_ERR( rval, "Getting target mesh IDs failed" );
00265 rval = get_vartag_data( mbCore, srcAreaTag, sets, srcArea_size, src_areas );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
00266 rval = get_vartag_data( mbCore, tgtAreaTag, sets, tgtArea_size, tgt_areas );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
00267
00268 assert( srcArea_size == srcID_size );
00269 assert( tgtArea_size == tgtID_size );
00270
00271 std::vector< double > src_glob_areas( nDofA, 0.0 ), tgt_glob_areas( nDofB, 0.0 );
00272 for( int i = 0; i < srcArea_size; ++i )
00273 {
00274 // printf("%d/%d: %d = Found ID %d and area %5.6e\n", i, srcArea_size, nDofA,
00275 // src_gids[i], src_areas[i]);
00276 assert( i < srcID_size );
00277 assert( src_gids[i] < nDofA );
00278 if( src_areas[i] > src_glob_areas[src_gids[i]] ) src_glob_areas[src_gids[i]] = src_areas[i];
00279 }
00280 for( int i = 0; i < tgtArea_size; ++i )
00281 {
00282 // printf("%d/%d: %d = Found ID %d and area %5.6e\n", i, tgtArea_size, nDofB,
00283 // tgt_gids[i], tgt_areas[i]);
00284 assert( i < tgtID_size );
00285 assert( tgt_gids[i] < nDofB );
00286 if( tgt_areas[i] > tgt_glob_areas[tgt_gids[i]] ) tgt_glob_areas[tgt_gids[i]] = tgt_areas[i];
00287 }
00288
00289 // Write output dimensions entries
00290 int nSrcGridDims = 1;
00291 int nDstGridDims = 1;
00292
00293 NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
00294 NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );
00295
00296 NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
00297 NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );
00298
00299 if( nA == nDofA )
00300 {
00301 varSrcGridDims->put( &nA, 1 );
00302 varSrcGridDims->add_att( "name0", "num_elem" );
00303 }
00304 else
00305 {
00306 varSrcGridDims->put( &nDofA, 1 );
00307 varSrcGridDims->add_att( "name1", "num_dof" );
00308 }
00309
00310 if( nB == nDofB )
00311 {
00312 varDstGridDims->put( &nB, 1 );
00313 varDstGridDims->add_att( "name0", "num_elem" );
00314 }
00315 else
00316 {
00317 varDstGridDims->put( &nDofB, 1 );
00318 varDstGridDims->add_att( "name1", "num_dof" );
00319 }
00320
00321 // Source and Target mesh resolutions
00322 NcDim* dimNA = ncMap.add_dim( "n_a", nDofA );
00323 NcDim* dimNB = ncMap.add_dim( "n_b", nDofB );
00324
00325 // Source and Target verticecs per elements
00326 const int nva = ( nA == nDofA ? nVA : 1 );
00327 const int nvb = ( nB == nDofB ? nVB : 1 );
00328 NcDim* dimNVA = ncMap.add_dim( "nv_a", nva );
00329 NcDim* dimNVB = ncMap.add_dim( "nv_b", nvb );
00330
00331 // Source and Target verticecs per elements
00332 // NcDim * dimNEA = ncMap.add_dim("ne_a", nA);
00333 // NcDim * dimNEB = ncMap.add_dim("ne_b", nB);
00334
00335 if( writeXYCoords )
00336 {
00337 // Write coordinates
00338 NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA /*dimNA*/ );
00339 NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB /*dimNB*/ );
00340
00341 NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA /*dimNA*/ );
00342 NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB /*dimNB*/ );
00343
00344 NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA /*dimNA*/, dimNVA );
00345 NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB /*dimNB*/, dimNVB );
00346
00347 NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA /*dimNA*/, dimNVA );
00348 NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB /*dimNB*/, dimNVB );
00349
00350 varYCA->add_att( "units", "degrees" );
00351 varYCB->add_att( "units", "degrees" );
00352
00353 varXCA->add_att( "units", "degrees" );
00354 varXCB->add_att( "units", "degrees" );
00355
00356 varYVA->add_att( "units", "degrees" );
00357 varYVB->add_att( "units", "degrees" );
00358
00359 varXVA->add_att( "units", "degrees" );
00360 varXVB->add_att( "units", "degrees" );
00361
00362 std::vector< double > src_centerlat, src_centerlon;
00363 int srccenter_size;
00364 rval = get_vartag_data( mbCore, srcCenterLat, sets, srccenter_size, src_centerlat );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
00365 rval = get_vartag_data( mbCore, srcCenterLon, sets, srccenter_size, src_centerlon );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
00366 std::vector< double > src_glob_centerlat( nDofA, 0.0 ), src_glob_centerlon( nDofA, 0.0 );
00367
00368 for( int i = 0; i < srccenter_size; ++i )
00369 {
00370 assert( i < srcID_size );
00371 assert( src_gids[i] < nDofA );
00372
00373 src_glob_centerlat[src_gids[i]] = src_centerlat[i];
00374 src_glob_centerlon[src_gids[i]] = src_centerlon[i];
00375 }
00376
00377 std::vector< double > tgt_centerlat, tgt_centerlon;
00378 int tgtcenter_size;
00379 rval = get_vartag_data( mbCore, tgtCenterLat, sets, tgtcenter_size, tgt_centerlat );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
00380 rval = get_vartag_data( mbCore, tgtCenterLon, sets, tgtcenter_size, tgt_centerlon );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
00381 std::vector< double > tgt_glob_centerlat( nDofB, 0.0 ), tgt_glob_centerlon( nDofB, 0.0 );
00382 for( int i = 0; i < tgtcenter_size; ++i )
00383 {
00384 assert( i < tgtID_size );
00385 assert( tgt_gids[i] < nDofB );
00386
00387 tgt_glob_centerlat[tgt_gids[i]] = tgt_centerlat[i];
00388 tgt_glob_centerlon[tgt_gids[i]] = tgt_centerlon[i];
00389 }
00390
00391 varYCA->put( &( src_glob_centerlat[0] ), nDofA );
00392 varYCB->put( &( tgt_glob_centerlat[0] ), nDofB );
00393 varXCA->put( &( src_glob_centerlon[0] ), nDofA );
00394 varXCB->put( &( tgt_glob_centerlon[0] ), nDofB );
00395
00396 src_centerlat.clear();
00397 src_centerlon.clear();
00398 tgt_centerlat.clear();
00399 tgt_centerlon.clear();
00400
00401 DataArray2D< double > src_glob_vertexlat( nDofA, nva ), src_glob_vertexlon( nDofA, nva );
00402 if( nva > 1 )
00403 {
00404 std::vector< double > src_vertexlat, src_vertexlon;
00405 int srcvertex_size;
00406 rval = get_vartag_data( mbCore, srcVertexLat, sets, srcvertex_size, src_vertexlat );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
00407 rval = get_vartag_data( mbCore, srcVertexLon, sets, srcvertex_size, src_vertexlon );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
00408 int offset = 0;
00409 for( unsigned vIndex = 0; vIndex < src_gids.size(); ++vIndex )
00410 {
00411 for( int vNV = 0; vNV < nva; ++vNV )
00412 {
00413 assert( offset < srcvertex_size );
00414 src_glob_vertexlat[src_gids[vIndex]][vNV] = src_vertexlat[offset];
00415 src_glob_vertexlon[src_gids[vIndex]][vNV] = src_vertexlon[offset];
00416 offset++;
00417 }
00418 }
00419 }
00420
00421 DataArray2D< double > tgt_glob_vertexlat( nDofB, nvb ), tgt_glob_vertexlon( nDofB, nvb );
00422 if( nvb > 1 )
00423 {
00424 std::vector< double > tgt_vertexlat, tgt_vertexlon;
00425 int tgtvertex_size;
00426 rval = get_vartag_data( mbCore, tgtVertexLat, sets, tgtvertex_size, tgt_vertexlat );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
00427 rval = get_vartag_data( mbCore, tgtVertexLon, sets, tgtvertex_size, tgt_vertexlon );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
00428 int offset = 0;
00429 for( unsigned vIndex = 0; vIndex < tgt_gids.size(); ++vIndex )
00430 {
00431 for( int vNV = 0; vNV < nvb; ++vNV )
00432 {
00433 assert( offset < tgtvertex_size );
00434 tgt_glob_vertexlat[tgt_gids[vIndex]][vNV] = tgt_vertexlat[offset];
00435 tgt_glob_vertexlon[tgt_gids[vIndex]][vNV] = tgt_vertexlon[offset];
00436 offset++;
00437 }
00438 }
00439 }
00440
00441 varYVA->put( &( src_glob_vertexlat[0][0] ), nDofA, nva );
00442 varYVB->put( &( tgt_glob_vertexlat[0][0] ), nDofB, nvb );
00443
00444 varXVA->put( &( src_glob_vertexlon[0][0] ), nDofA, nva );
00445 varXVB->put( &( tgt_glob_vertexlon[0][0] ), nDofB, nvb );
00446 }
00447
00448 // Write areas
00449 NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
00450 varAreaA->put( &( src_glob_areas[0] ), nDofA );
00451 // varAreaA->add_att("units", "steradians");
00452
00453 NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
00454 varAreaB->put( &( tgt_glob_areas[0] ), nDofB );
00455 // varAreaB->add_att("units", "steradians");
00456
00457 std::vector< int > mat_rows, mat_cols;
00458 std::vector< double > mat_vals;
00459 int row_sizes, col_sizes, val_sizes;
00460 rval = get_vartag_data( mbCore, smatRowdataTag, sets, row_sizes, mat_rows );MB_CHK_SET_ERR( rval, "Getting matrix row data failed" );
00461 assert( row_sizes == NNZ );
00462 rval = get_vartag_data( mbCore, smatColdataTag, sets, col_sizes, mat_cols );MB_CHK_SET_ERR( rval, "Getting matrix col data failed" );
00463 assert( col_sizes == NNZ );
00464 rval = get_vartag_data( mbCore, smatValsdataTag, sets, val_sizes, mat_vals );MB_CHK_SET_ERR( rval, "Getting matrix values failed" );
00465 assert( val_sizes == NNZ );
00466
00467 // Let us form the matrix in-memory and consolidate shared DoF rows from shared-process
00468 // contributions
00469 SparseMatrix< double > mapMatrix;
00470
00471 for( int innz = 0; innz < NNZ; ++innz )
00472 {
00473 #ifdef VERBOSE
00474 if( fabs( mapMatrix( mat_rows[innz], mat_cols[innz] ) ) > 1e-12 )
00475 {
00476 printf( "Adding to existing loc: (%d, %d) = %12.8f\n", mat_rows[innz], mat_cols[innz],
00477 mapMatrix( mat_rows[innz], mat_cols[innz] ) );
00478 }
00479 #endif
00480 mapMatrix( mat_rows[innz], mat_cols[innz] ) += mat_vals[innz];
00481 }
00482
00483 // Write SparseMatrix entries
00484 DataArray1D< int > vecRow;
00485 DataArray1D< int > vecCol;
00486 DataArray1D< double > vecS;
00487
00488 mapMatrix.GetEntries( vecRow, vecCol, vecS );
00489
00490 int nS = vecS.GetRows();
00491
00492 // Print more information about what we are converting:
00493 // Source elements/vertices/type (Discretization ?)
00494 // Target elements/vertices/type (Discretization ?)
00495 // Overlap elements/types
00496 // Rmeapping weights matrix: rows/cols/NNZ
00497 // Output the number of sets
00498 printf( "Primary sets: %15zu\n", sets.size() );
00499 printf( "Original NNZ: %18d\n", NNZ );
00500 printf( "Consolidated Total NNZ: %8d\n", nS );
00501 printf( "Conservative weights ? %6d\n", ( bConserved > 0 ) );
00502 printf( "Monotone weights ? %10d\n", ( bMonotonicity > 0 ) );
00503
00504 printf( "\n--------------------------------------------------------------\n" );
00505 printf( "%20s %21s %15s\n", "Description", "Source", "Target" );
00506 printf( "--------------------------------------------------------------\n" );
00507
00508 printf( "%25s %15d %15d\n", "Number of elements:", nA, nB );
00509 printf( "%25s %15d %15d\n", "Number of DoFs:", nDofA, nDofB );
00510 printf( "%25s %15d %15d\n", "Maximum vertex/element:", nVA, nVB );
00511 printf( "%25s %15s %15s\n", "Discretization type:", methodA.c_str(), methodB.c_str() );
00512 printf( "%25s %15d %15d\n", "Discretization order:", nOrdA, nOrdB );
00513
00514 // Calculate and write fractional coverage arrays
00515 {
00516 DataArray1D< double > dFracA( nDofA );
00517 DataArray1D< double > dFracB( nDofB );
00518
00519 for( int i = 0; i < nS; i++ )
00520 {
00521 // std::cout << i << " - mat_vals = " << mat_vals[i] << " dFracA = " << mat_vals[i]
00522 // / src_glob_areas[vecCol[i]] * tgt_glob_areas[vecRow[i]] << std::endl;
00523 dFracA[vecCol[i]] += vecS[i] / src_glob_areas[vecCol[i]] * tgt_glob_areas[vecRow[i]];
00524 dFracB[vecRow[i]] += vecS[i];
00525 }
00526
00527 NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
00528 varFracA->put( &( dFracA[0] ), nDofA );
00529 varFracA->add_att( "name", "fraction of target coverage of source dof" );
00530 varFracA->add_att( "units", "unitless" );
00531
00532 NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
00533 varFracB->put( &( dFracB[0] ), nDofB );
00534 varFracB->add_att( "name", "fraction of source coverage of target dof" );
00535 varFracB->add_att( "units", "unitless" );
00536 }
00537
00538 // Write out data
00539 NcDim* dimNS = ncMap.add_dim( "n_s", nS );
00540
00541 NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
00542 varRow->add_att( "name", "sparse matrix target dof index" );
00543 varRow->add_att( "first_index", "1" );
00544
00545 NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
00546 varCol->add_att( "name", "sparse matrix source dof index" );
00547 varCol->add_att( "first_index", "1" );
00548
00549 NcVar* varS = ncMap.add_var( "S", ncDouble, dimNS );
00550 varS->add_att( "name", "sparse matrix coefficient" );
00551
00552 // Increment vecRow and vecCol: make it 1-based
00553 for( int i = 0; i < nS; i++ )
00554 {
00555 vecRow[i]++;
00556 vecCol[i]++;
00557 }
00558
00559 varRow->set_cur( (long)0 );
00560 varRow->put( &( vecRow[0] ), nS );
00561
00562 varCol->set_cur( (long)0 );
00563 varCol->put( &( vecCol[0] ), nS );
00564
00565 varS->set_cur( (long)0 );
00566 varS->put( &( vecS[0] ), nS );
00567
00568 ncMap.close();
00569
00570 // rval = mbCore->write_file(scripfile.c_str());MB_CHK_ERR(rval);
00571 }
00572 catch( std::exception& e )
00573 {
00574 std::cout << " exception caught during tree initialization " << e.what() << std::endl;
00575 }
00576 delete mbCore;
00577
00578 #ifdef MOAB_HAVE_MPI
00579 MPI_Finalize();
00580 #endif
00581
00582 exit( 0 );
00583 }
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.