 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
#include <iostream>#include <exception>#include <cmath>#include <cassert>#include <vector>#include <string>#include <fstream>#include <iomanip>#include "moab/ProgOptions.hpp"#include "moab/Core.hpp"#include "OfflineMap.h"#include "netcdfcpp.h"#include "NetCDFUtilities.h"#include "DataArray2D.h"
Include dependency graph for h5mtoscrip.cpp:Go to the source code of this file.
Defines | |
| #define | DTYPE(a) |
Functions | |
| template<typename T > | |
| ErrorCode | get_vartag_data (moab::Interface *mbCore, Tag tag, moab::Range &sets, int &out_data_size, std::vector< T > &data) |
| void | ReadFileMetaData (std::string &metaFilename, std::map< std::string, std::string > &metadataVals) |
| int | main (int argc, char *argv[]) |
| #define DTYPE | ( | a | ) |
{ \
( ( ( a ) == 0 ) ? "FV" : ( ( ( a ) == 1 ) ? "cGLL" : "dGLL" ) ) \
}
Referenced by main().
| ErrorCode get_vartag_data | ( | moab::Interface * | mbCore, |
| Tag | tag, | ||
| moab::Range & | sets, | ||
| int & | out_data_size, | ||
| std::vector< T > & | data | ||
| ) |
Definition at line 33 of file h5mtoscrip.cpp.
References ErrorCode, MB_CHK_SET_ERR, MB_SUCCESS, moab::Range::size(), T, and moab::Interface::tag_get_by_ptr().
Referenced by main().
{
int* tag_sizes = new int[sets.size()];
const void** tag_data = (const void**)new void*[sets.size()];
ErrorCode rval = mbCore->tag_get_by_ptr( tag, sets, tag_data, tag_sizes );MB_CHK_SET_ERR( rval, "Getting matrix rows failed" );
out_data_size = 0;
for( unsigned is = 0; is < sets.size(); ++is )
out_data_size += tag_sizes[is];
data.resize( out_data_size );
int ioffset = 0;
for( unsigned index = 0; index < sets.size(); index++ )
{
T* m_vals = (T*)tag_data[index];
for( int k = 0; k < tag_sizes[index]; k++ )
{
data[ioffset++] = m_vals[k];
}
}
return moab::MB_SUCCESS;
}
| int main | ( | int | argc, |
| char * | argv[] | ||
| ) |
Definition at line 81 of file h5mtoscrip.cpp.
References ProgOptions::addOpt(), DTYPE, ErrorCode, moab::Interface::get_entities_by_dimension(), moab::Interface::get_entities_by_type_and_tag(), get_vartag_data(), moab::Interface::globalId_tag(), moab::Interface::load_mesh(), MB_CHK_ERR, MB_CHK_SET_ERR, MB_TAG_SPARSE, MB_TYPE_INTEGER, MBENTITYSET, ProgOptions::parseCommandLine(), ProgOptions::printHelp(), ReadFileMetaData(), moab::Range::size(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), and moab::Interface::UNION.
{
moab::ErrorCode rval;
int dimension = 2;
NcError error2( NcError::verbose_nonfatal );
std::stringstream sstr;
ProgOptions opts;
std::string h5mfilename, scripfile;
bool noMap = false;
bool writeXYCoords = false;
#ifdef MOAB_HAVE_MPI
MPI_Init( &argc, &argv );
#endif
opts.addOpt< std::string >( "weights,w", "h5m remapping weights filename", &h5mfilename );
opts.addOpt< std::string >( "scrip,s", "Output SCRIP map filename", &scripfile );
opts.addOpt< int >( "dim,d", "Dimension of entities to use for partitioning", &dimension );
opts.addOpt< void >( "mesh,m", "Only convert the mesh and exclude the remap weight details", &noMap );
opts.addOpt< void >( "coords,c", "Write the center and vertex coordinates in lat/lon format", &writeXYCoords );
opts.parseCommandLine( argc, argv );
if( h5mfilename.empty() || scripfile.empty() )
{
opts.printHelp();
exit( 1 );
}
moab::Interface* mbCore = new( std::nothrow ) moab::Core;
if( NULL == mbCore )
{
return 1;
}
// Set the read options for parallel file loading
const std::string partition_set_name = "PARALLEL_PARTITION";
const std::string global_id_name = "GLOBAL_ID";
// Load file
rval = mbCore->load_mesh( h5mfilename.c_str() );MB_CHK_ERR( rval );
try
{
// Temporarily change rval reporting
NcError error_temp( NcError::verbose_fatal );
// Open an output file
NcFile ncMap( scripfile.c_str(), NcFile::Replace, NULL, 0, NcFile::Offset64Bits );
if( !ncMap.is_valid() )
{
_EXCEPTION1( "Unable to open output map file \"%s\"", scripfile.c_str() );
}
{
// NetCDF-SCRIP Global Attributes
std::map< std::string, std::string > mapAttributes;
size_t lastindex = h5mfilename.find_last_of( "." );
std::stringstream sstr;
sstr << h5mfilename.substr( 0, lastindex ) << ".meta";
std::string metaFilename = sstr.str();
ReadFileMetaData( metaFilename, mapAttributes );
mapAttributes["Command"] =
"Converted with MOAB:h5mtoscrip with --w=" + h5mfilename + " and --s=" + scripfile;
// Add global attributes
std::map< std::string, std::string >::const_iterator iterAttributes = mapAttributes.begin();
for( ; iterAttributes != mapAttributes.end(); iterAttributes++ )
{
std::cout << iterAttributes->first << " -- " << iterAttributes->second << std::endl;
ncMap.add_att( iterAttributes->first.c_str(), iterAttributes->second.c_str() );
}
std::cout << "\n";
}
Tag globalIDTag, materialSetTag;
globalIDTag = mbCore->globalId_tag();
// materialSetTag = mbCore->material_tag();
rval = mbCore->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, materialSetTag, MB_TAG_SPARSE );MB_CHK_ERR( rval );
// Get sets entities, by type
moab::Range meshsets;
rval = mbCore->get_entities_by_type_and_tag( 0, MBENTITYSET, &globalIDTag, NULL, 1, meshsets,
moab::Interface::UNION, true );MB_CHK_ERR( rval );
moab::EntityHandle rootset = 0;
///////////////////////////////////////////////////////////////////////////
// The metadata in H5M file contains the following data:
//
// 1. n_a: Total source entities: (number of elements in source mesh)
// 2. n_b: Total target entities: (number of elements in target mesh)
// 3. nv_a: Max edge size of elements in source mesh
// 4. nv_b: Max edge size of elements in target mesh
// 5. maxrows: Number of rows in remap weight matrix
// 6. maxcols: Number of cols in remap weight matrix
// 7. nnz: Number of total nnz in sparse remap weight matrix
// 8. np_a: The order of the field description on the source mesh: >= 1
// 9. np_b: The order of the field description on the target mesh: >= 1
// 10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2
// = dGLL]
// 11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2
// = dGLL]
// 12. conserved: Flag to specify whether the remap operator has conservation constraints:
// [0, 1]
// 13. monotonicity: Flags to specify whether the remap operator has monotonicity
// constraints: [0, 1, 2]
//
///////////////////////////////////////////////////////////////////////////
Tag smatMetadataTag;
int smat_metadata_glb[13];
rval = mbCore->tag_get_handle( "SMAT_DATA", 13, MB_TYPE_INTEGER, smatMetadataTag, MB_TAG_SPARSE );MB_CHK_ERR( rval );
rval = mbCore->tag_get_data( smatMetadataTag, &rootset, 1, smat_metadata_glb );MB_CHK_ERR( rval );
// std::cout << "Number of mesh sets is " << meshsets.size() << std::endl;
#define DTYPE( a ) \
{ \
( ( ( a ) == 0 ) ? "FV" : ( ( ( a ) == 1 ) ? "cGLL" : "dGLL" ) ) \
}
// Map dimensions
int nA = smat_metadata_glb[0];
int nB = smat_metadata_glb[1];
int nVA = smat_metadata_glb[2];
int nVB = smat_metadata_glb[3];
int nDofB = smat_metadata_glb[4];
int nDofA = smat_metadata_glb[5];
int NNZ = smat_metadata_glb[6];
int nOrdA = smat_metadata_glb[7];
int nOrdB = smat_metadata_glb[8];
int nBasA = smat_metadata_glb[9];
std::string methodA = DTYPE( nBasA );
int nBasB = smat_metadata_glb[10];
std::string methodB = DTYPE( nBasB );
int bConserved = smat_metadata_glb[11];
int bMonotonicity = smat_metadata_glb[12];
EntityHandle source_mesh = 0, target_mesh = 0, overlap_mesh = 0;
for( unsigned im = 0; im < meshsets.size(); ++im )
{
moab::Range elems;
rval = mbCore->get_entities_by_dimension( meshsets[im], 2, elems );MB_CHK_ERR( rval );
if( elems.size() - nA == 0 && source_mesh == 0 )
source_mesh = meshsets[im];
else if( elems.size() - nB == 0 && target_mesh == 0 )
target_mesh = meshsets[im];
else if( overlap_mesh == 0 )
overlap_mesh = meshsets[im];
else
continue;
}
Tag srcIDTag, srcAreaTag, tgtIDTag, tgtAreaTag;
rval = mbCore->tag_get_handle( "SourceGIDS", srcIDTag );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "SourceAreas", srcAreaTag );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "TargetGIDS", tgtIDTag );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "TargetAreas", tgtAreaTag );MB_CHK_ERR( rval );
Tag smatRowdataTag, smatColdataTag, smatValsdataTag;
rval = mbCore->tag_get_handle( "SMAT_ROWS", smatRowdataTag );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "SMAT_COLS", smatColdataTag );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "SMAT_VALS", smatValsdataTag );MB_CHK_ERR( rval );
Tag srcCenterLon, srcCenterLat, tgtCenterLon, tgtCenterLat;
rval = mbCore->tag_get_handle( "SourceCoordCenterLon", srcCenterLon );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "SourceCoordCenterLat", srcCenterLat );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "TargetCoordCenterLon", tgtCenterLon );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "TargetCoordCenterLat", tgtCenterLat );MB_CHK_ERR( rval );
Tag srcVertexLon, srcVertexLat, tgtVertexLon, tgtVertexLat;
rval = mbCore->tag_get_handle( "SourceCoordVertexLon", srcVertexLon );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "SourceCoordVertexLat", srcVertexLat );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "TargetCoordVertexLon", tgtVertexLon );MB_CHK_ERR( rval );
rval = mbCore->tag_get_handle( "TargetCoordVertexLat", tgtVertexLat );MB_CHK_ERR( rval );
// Get sets entities, by type
moab::Range sets;
// rval = mbCore->get_entities_by_type(0, MBENTITYSET, sets);MB_CHK_ERR(rval);
rval = mbCore->get_entities_by_type_and_tag( 0, MBENTITYSET, &smatRowdataTag, NULL, 1, sets,
moab::Interface::UNION, true );MB_CHK_ERR( rval );
std::vector< int > src_gids, tgt_gids;
std::vector< double > src_areas, tgt_areas;
int srcID_size, tgtID_size, srcArea_size, tgtArea_size;
rval = get_vartag_data( mbCore, srcIDTag, sets, srcID_size, src_gids );MB_CHK_SET_ERR( rval, "Getting source mesh IDs failed" );
rval = get_vartag_data( mbCore, tgtIDTag, sets, tgtID_size, tgt_gids );MB_CHK_SET_ERR( rval, "Getting target mesh IDs failed" );
rval = get_vartag_data( mbCore, srcAreaTag, sets, srcArea_size, src_areas );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
rval = get_vartag_data( mbCore, tgtAreaTag, sets, tgtArea_size, tgt_areas );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
assert( srcArea_size == srcID_size );
assert( tgtArea_size == tgtID_size );
std::vector< double > src_glob_areas( nDofA, 0.0 ), tgt_glob_areas( nDofB, 0.0 );
for( int i = 0; i < srcArea_size; ++i )
{
// printf("%d/%d: %d = Found ID %d and area %5.6e\n", i, srcArea_size, nDofA,
// src_gids[i], src_areas[i]);
assert( i < srcID_size );
assert( src_gids[i] < nDofA );
if( src_areas[i] > src_glob_areas[src_gids[i]] ) src_glob_areas[src_gids[i]] = src_areas[i];
}
for( int i = 0; i < tgtArea_size; ++i )
{
// printf("%d/%d: %d = Found ID %d and area %5.6e\n", i, tgtArea_size, nDofB,
// tgt_gids[i], tgt_areas[i]);
assert( i < tgtID_size );
assert( tgt_gids[i] < nDofB );
if( tgt_areas[i] > tgt_glob_areas[tgt_gids[i]] ) tgt_glob_areas[tgt_gids[i]] = tgt_areas[i];
}
// Write output dimensions entries
int nSrcGridDims = 1;
int nDstGridDims = 1;
NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );
NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );
if( nA == nDofA )
{
varSrcGridDims->put( &nA, 1 );
varSrcGridDims->add_att( "name0", "num_elem" );
}
else
{
varSrcGridDims->put( &nDofA, 1 );
varSrcGridDims->add_att( "name1", "num_dof" );
}
if( nB == nDofB )
{
varDstGridDims->put( &nB, 1 );
varDstGridDims->add_att( "name0", "num_elem" );
}
else
{
varDstGridDims->put( &nDofB, 1 );
varDstGridDims->add_att( "name1", "num_dof" );
}
// Source and Target mesh resolutions
NcDim* dimNA = ncMap.add_dim( "n_a", nDofA );
NcDim* dimNB = ncMap.add_dim( "n_b", nDofB );
// Source and Target verticecs per elements
const int nva = ( nA == nDofA ? nVA : 1 );
const int nvb = ( nB == nDofB ? nVB : 1 );
NcDim* dimNVA = ncMap.add_dim( "nv_a", nva );
NcDim* dimNVB = ncMap.add_dim( "nv_b", nvb );
// Source and Target verticecs per elements
// NcDim * dimNEA = ncMap.add_dim("ne_a", nA);
// NcDim * dimNEB = ncMap.add_dim("ne_b", nB);
if( writeXYCoords )
{
// Write coordinates
NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA /*dimNA*/ );
NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB /*dimNB*/ );
NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA /*dimNA*/ );
NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB /*dimNB*/ );
NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA /*dimNA*/, dimNVA );
NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB /*dimNB*/, dimNVB );
NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA /*dimNA*/, dimNVA );
NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB /*dimNB*/, dimNVB );
varYCA->add_att( "units", "degrees" );
varYCB->add_att( "units", "degrees" );
varXCA->add_att( "units", "degrees" );
varXCB->add_att( "units", "degrees" );
varYVA->add_att( "units", "degrees" );
varYVB->add_att( "units", "degrees" );
varXVA->add_att( "units", "degrees" );
varXVB->add_att( "units", "degrees" );
std::vector< double > src_centerlat, src_centerlon;
int srccenter_size;
rval = get_vartag_data( mbCore, srcCenterLat, sets, srccenter_size, src_centerlat );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
rval = get_vartag_data( mbCore, srcCenterLon, sets, srccenter_size, src_centerlon );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
std::vector< double > src_glob_centerlat( nDofA, 0.0 ), src_glob_centerlon( nDofA, 0.0 );
for( int i = 0; i < srccenter_size; ++i )
{
assert( i < srcID_size );
assert( src_gids[i] < nDofA );
src_glob_centerlat[src_gids[i]] = src_centerlat[i];
src_glob_centerlon[src_gids[i]] = src_centerlon[i];
}
std::vector< double > tgt_centerlat, tgt_centerlon;
int tgtcenter_size;
rval = get_vartag_data( mbCore, tgtCenterLat, sets, tgtcenter_size, tgt_centerlat );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
rval = get_vartag_data( mbCore, tgtCenterLon, sets, tgtcenter_size, tgt_centerlon );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
std::vector< double > tgt_glob_centerlat( nDofB, 0.0 ), tgt_glob_centerlon( nDofB, 0.0 );
for( int i = 0; i < tgtcenter_size; ++i )
{
assert( i < tgtID_size );
assert( tgt_gids[i] < nDofB );
tgt_glob_centerlat[tgt_gids[i]] = tgt_centerlat[i];
tgt_glob_centerlon[tgt_gids[i]] = tgt_centerlon[i];
}
varYCA->put( &( src_glob_centerlat[0] ), nDofA );
varYCB->put( &( tgt_glob_centerlat[0] ), nDofB );
varXCA->put( &( src_glob_centerlon[0] ), nDofA );
varXCB->put( &( tgt_glob_centerlon[0] ), nDofB );
src_centerlat.clear();
src_centerlon.clear();
tgt_centerlat.clear();
tgt_centerlon.clear();
DataArray2D< double > src_glob_vertexlat( nDofA, nva ), src_glob_vertexlon( nDofA, nva );
if( nva > 1 )
{
std::vector< double > src_vertexlat, src_vertexlon;
int srcvertex_size;
rval = get_vartag_data( mbCore, srcVertexLat, sets, srcvertex_size, src_vertexlat );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
rval = get_vartag_data( mbCore, srcVertexLon, sets, srcvertex_size, src_vertexlon );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
int offset = 0;
for( unsigned vIndex = 0; vIndex < src_gids.size(); ++vIndex )
{
for( int vNV = 0; vNV < nva; ++vNV )
{
assert( offset < srcvertex_size );
src_glob_vertexlat[src_gids[vIndex]][vNV] = src_vertexlat[offset];
src_glob_vertexlon[src_gids[vIndex]][vNV] = src_vertexlon[offset];
offset++;
}
}
}
DataArray2D< double > tgt_glob_vertexlat( nDofB, nvb ), tgt_glob_vertexlon( nDofB, nvb );
if( nvb > 1 )
{
std::vector< double > tgt_vertexlat, tgt_vertexlon;
int tgtvertex_size;
rval = get_vartag_data( mbCore, tgtVertexLat, sets, tgtvertex_size, tgt_vertexlat );MB_CHK_SET_ERR( rval, "Getting source mesh areas failed" );
rval = get_vartag_data( mbCore, tgtVertexLon, sets, tgtvertex_size, tgt_vertexlon );MB_CHK_SET_ERR( rval, "Getting target mesh areas failed" );
int offset = 0;
for( unsigned vIndex = 0; vIndex < tgt_gids.size(); ++vIndex )
{
for( int vNV = 0; vNV < nvb; ++vNV )
{
assert( offset < tgtvertex_size );
tgt_glob_vertexlat[tgt_gids[vIndex]][vNV] = tgt_vertexlat[offset];
tgt_glob_vertexlon[tgt_gids[vIndex]][vNV] = tgt_vertexlon[offset];
offset++;
}
}
}
varYVA->put( &( src_glob_vertexlat[0][0] ), nDofA, nva );
varYVB->put( &( tgt_glob_vertexlat[0][0] ), nDofB, nvb );
varXVA->put( &( src_glob_vertexlon[0][0] ), nDofA, nva );
varXVB->put( &( tgt_glob_vertexlon[0][0] ), nDofB, nvb );
}
// Write areas
NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
varAreaA->put( &( src_glob_areas[0] ), nDofA );
// varAreaA->add_att("units", "steradians");
NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
varAreaB->put( &( tgt_glob_areas[0] ), nDofB );
// varAreaB->add_att("units", "steradians");
std::vector< int > mat_rows, mat_cols;
std::vector< double > mat_vals;
int row_sizes, col_sizes, val_sizes;
rval = get_vartag_data( mbCore, smatRowdataTag, sets, row_sizes, mat_rows );MB_CHK_SET_ERR( rval, "Getting matrix row data failed" );
assert( row_sizes == NNZ );
rval = get_vartag_data( mbCore, smatColdataTag, sets, col_sizes, mat_cols );MB_CHK_SET_ERR( rval, "Getting matrix col data failed" );
assert( col_sizes == NNZ );
rval = get_vartag_data( mbCore, smatValsdataTag, sets, val_sizes, mat_vals );MB_CHK_SET_ERR( rval, "Getting matrix values failed" );
assert( val_sizes == NNZ );
// Let us form the matrix in-memory and consolidate shared DoF rows from shared-process
// contributions
SparseMatrix< double > mapMatrix;
for( int innz = 0; innz < NNZ; ++innz )
{
#ifdef VERBOSE
if( fabs( mapMatrix( mat_rows[innz], mat_cols[innz] ) ) > 1e-12 )
{
printf( "Adding to existing loc: (%d, %d) = %12.8f\n", mat_rows[innz], mat_cols[innz],
mapMatrix( mat_rows[innz], mat_cols[innz] ) );
}
#endif
mapMatrix( mat_rows[innz], mat_cols[innz] ) += mat_vals[innz];
}
// Write SparseMatrix entries
DataArray1D< int > vecRow;
DataArray1D< int > vecCol;
DataArray1D< double > vecS;
mapMatrix.GetEntries( vecRow, vecCol, vecS );
int nS = vecS.GetRows();
// Print more information about what we are converting:
// Source elements/vertices/type (Discretization ?)
// Target elements/vertices/type (Discretization ?)
// Overlap elements/types
// Rmeapping weights matrix: rows/cols/NNZ
// Output the number of sets
printf( "Primary sets: %15zu\n", sets.size() );
printf( "Original NNZ: %18d\n", NNZ );
printf( "Consolidated Total NNZ: %8d\n", nS );
printf( "Conservative weights ? %6d\n", ( bConserved > 0 ) );
printf( "Monotone weights ? %10d\n", ( bMonotonicity > 0 ) );
printf( "\n--------------------------------------------------------------\n" );
printf( "%20s %21s %15s\n", "Description", "Source", "Target" );
printf( "--------------------------------------------------------------\n" );
printf( "%25s %15d %15d\n", "Number of elements:", nA, nB );
printf( "%25s %15d %15d\n", "Number of DoFs:", nDofA, nDofB );
printf( "%25s %15d %15d\n", "Maximum vertex/element:", nVA, nVB );
printf( "%25s %15s %15s\n", "Discretization type:", methodA.c_str(), methodB.c_str() );
printf( "%25s %15d %15d\n", "Discretization order:", nOrdA, nOrdB );
// Calculate and write fractional coverage arrays
{
DataArray1D< double > dFracA( nDofA );
DataArray1D< double > dFracB( nDofB );
for( int i = 0; i < nS; i++ )
{
// std::cout << i << " - mat_vals = " << mat_vals[i] << " dFracA = " << mat_vals[i]
// / src_glob_areas[vecCol[i]] * tgt_glob_areas[vecRow[i]] << std::endl;
dFracA[vecCol[i]] += vecS[i] / src_glob_areas[vecCol[i]] * tgt_glob_areas[vecRow[i]];
dFracB[vecRow[i]] += vecS[i];
}
NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
varFracA->put( &( dFracA[0] ), nDofA );
varFracA->add_att( "name", "fraction of target coverage of source dof" );
varFracA->add_att( "units", "unitless" );
NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
varFracB->put( &( dFracB[0] ), nDofB );
varFracB->add_att( "name", "fraction of source coverage of target dof" );
varFracB->add_att( "units", "unitless" );
}
// Write out data
NcDim* dimNS = ncMap.add_dim( "n_s", nS );
NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
varRow->add_att( "name", "sparse matrix target dof index" );
varRow->add_att( "first_index", "1" );
NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
varCol->add_att( "name", "sparse matrix source dof index" );
varCol->add_att( "first_index", "1" );
NcVar* varS = ncMap.add_var( "S", ncDouble, dimNS );
varS->add_att( "name", "sparse matrix coefficient" );
// Increment vecRow and vecCol: make it 1-based
for( int i = 0; i < nS; i++ )
{
vecRow[i]++;
vecCol[i]++;
}
varRow->set_cur( (long)0 );
varRow->put( &( vecRow[0] ), nS );
varCol->set_cur( (long)0 );
varCol->put( &( vecCol[0] ), nS );
varS->set_cur( (long)0 );
varS->put( &( vecS[0] ), nS );
ncMap.close();
// rval = mbCore->write_file(scripfile.c_str());MB_CHK_ERR(rval);
}
catch( std::exception& e )
{
std::cout << " exception caught during tree initialization " << e.what() << std::endl;
}
delete mbCore;
#ifdef MOAB_HAVE_MPI
MPI_Finalize();
#endif
exit( 0 );
}
| void ReadFileMetaData | ( | std::string & | metaFilename, |
| std::map< std::string, std::string > & | metadataVals | ||
| ) |
Definition at line 62 of file h5mtoscrip.cpp.
Referenced by main().
{
std::ifstream metafile;
std::string line;
metafile.open( metaFilename.c_str() );
metadataVals["Title"] = "MOAB-TempestRemap (MBTR) Offline Regridding Weight Converter (h5mtoscrip)";
std::string key, value;
while( std::getline( metafile, line ) )
{
size_t lastindex = line.find_last_of( "=" );
key = line.substr( 0, lastindex - 1 );
value = line.substr( lastindex + 2, line.length() );
metadataVals[std::string( key )] = std::string( value );
}
metafile.close();
}
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.