Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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,
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines