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