Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
h5mtoscrip.cpp
Go to the documentation of this file.
1 // Usage:
2 // tools/h5mtoscrip -w map_atm_to_ocn.h5m -s map_atm_to_ocn.nc --coords
3 //
4 #include <iostream>
5 #include <exception>
6 #include <cmath>
7 #include <cassert>
8 #include <vector>
9 #include <string>
10 #include <fstream>
11 #include <iomanip>
12 
13 #include "moab/ProgOptions.hpp"
14 #include "moab/Core.hpp"
15 
16 #ifdef MOAB_HAVE_MPI
17 #include "moab_mpi.h"
18 #endif
19 
20 #ifndef MOAB_HAVE_TEMPESTREMAP
21 #error Tool requires compilation with TempestRemap dependency
22 #endif
23 
24 // TempestRemap includes
25 #include "OfflineMap.h"
26 #include "netcdfcpp.h"
27 #include "NetCDFUtilities.h"
28 #include "DataArray2D.h"
29 
30 using namespace moab;
31 
32 template < typename T >
33 ErrorCode get_vartag_data( moab::Interface* mbCore,
34  Tag tag,
35  moab::Range& sets,
36  int& out_data_size,
37  std::vector< T >& data )
38 {
39  int* tag_sizes = new int[sets.size()];
40  const void** tag_data = (const void**)new void*[sets.size()];
41 
42  MB_CHK_SET_ERR( mbCore->tag_get_by_ptr( tag, sets, tag_data, tag_sizes ), "Getting matrix rows failed" );
43 
44  out_data_size = 0;
45  for( unsigned is = 0; is < sets.size(); ++is )
46  out_data_size += tag_sizes[is];
47 
48  data.resize( out_data_size );
49  int ioffset = 0;
50  for( unsigned index = 0; index < sets.size(); index++ )
51  {
52  T* m_vals = (T*)tag_data[index];
53  for( int k = 0; k < tag_sizes[index]; k++ )
54  {
55  data[ioffset++] = m_vals[k];
56  }
57  }
58 
59  return moab::MB_SUCCESS;
60 }
61 
62 void ReadFileMetaData( std::string& metaFilename, std::map< std::string, std::string >& metadataVals )
63 {
64  std::ifstream metafile;
65  std::string line;
66 
67  metafile.open( metaFilename.c_str() );
68  metadataVals["Title"] = "MOAB-TempestRemap (MBTR) Offline Regridding Weight Converter (h5mtoscrip)";
69  std::string key, value;
70  while( std::getline( metafile, line ) )
71  {
72  size_t lastindex = line.find_last_of( "=" );
73  key = line.substr( 0, lastindex - 1 );
74  value = line.substr( lastindex + 2, line.length() );
75 
76  metadataVals[std::string( key )] = std::string( value );
77  }
78  metafile.close();
79 }
80 
81 int main( int argc, char* argv[] )
82 {
83  int dimension = 2;
84  NcError error2( NcError::verbose_nonfatal );
85  std::stringstream sstr;
86  ProgOptions opts;
87  std::string h5mfilename, scripfile;
88  bool noMap = false;
89  bool writeXYCoords = false;
90 
91 #ifdef MOAB_HAVE_MPI
92  MPI_Init( &argc, &argv );
93 #endif
94 
95  opts.addOpt< std::string >( "weights,w", "h5m remapping weights filename", &h5mfilename );
96  opts.addOpt< std::string >( "scrip,s", "Output SCRIP map filename", &scripfile );
97  opts.addOpt< int >( "dim,d", "Dimension of entities to use for partitioning", &dimension );
98  opts.addOpt< void >( "mesh,m", "Only convert the mesh and exclude the remap weight details", &noMap );
99  opts.addOpt< void >( "coords,c", "Write the center and vertex coordinates in lat/lon format", &writeXYCoords );
100 
101  opts.parseCommandLine( argc, argv );
102 
103  if( h5mfilename.empty() || scripfile.empty() )
104  {
105  opts.printHelp();
106  exit( 1 );
107  }
108 
109  moab::Interface* mbCore = new( std::nothrow ) moab::Core;
110 
111  if( NULL == mbCore )
112  {
113  return 1;
114  }
115 
116  // Set the read options for parallel file loading
117  const std::string partition_set_name = "PARALLEL_PARTITION";
118  const std::string global_id_name = "GLOBAL_ID";
119 
120  // Load file
121  MB_CHK_ERR( mbCore->load_mesh( h5mfilename.c_str() ) );
122 
123  try
124  {
125  // Temporarily change rval reporting
126  NcError error_temp( NcError::verbose_fatal );
127 
128  // Open an output file
129  NcFile ncMap( scripfile.c_str(), NcFile::Replace, NULL, 0, NcFile::Offset64Bits );
130  if( !ncMap.is_valid() )
131  {
132  _EXCEPTION1( "Unable to open output map file \"%s\"", scripfile.c_str() );
133  }
134 
135  {
136  // NetCDF-SCRIP Global Attributes
137  std::map< std::string, std::string > mapAttributes;
138  size_t lastindex = h5mfilename.find_last_of( "." );
139  std::stringstream sstr;
140  sstr << h5mfilename.substr( 0, lastindex ) << ".meta";
141  std::string metaFilename = sstr.str();
142  ReadFileMetaData( metaFilename, mapAttributes );
143  mapAttributes["Command"] =
144  "Converted with MOAB:h5mtoscrip with --w=" + h5mfilename + " and --s=" + scripfile;
145 
146  // Add global attributes
147  std::map< std::string, std::string >::const_iterator iterAttributes = mapAttributes.begin();
148  for( ; iterAttributes != mapAttributes.end(); iterAttributes++ )
149  {
150 
151  std::cout << iterAttributes->first << " -- " << iterAttributes->second << std::endl;
152  ncMap.add_att( iterAttributes->first.c_str(), iterAttributes->second.c_str() );
153  }
154  std::cout << "\n";
155  }
156 
157  Tag globalIDTag, materialSetTag;
158  globalIDTag = mbCore->globalId_tag();
159  // materialSetTag = mbCore->material_tag();
160  MB_CHK_ERR( mbCore->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, materialSetTag, MB_TAG_SPARSE ) );
161 
162  // Get sets entities, by type
163  moab::Range meshsets;
164  MB_CHK_ERR( mbCore->get_entities_by_type_and_tag( 0, MBENTITYSET, &globalIDTag, NULL, 1, meshsets,
165  moab::Interface::UNION, true ) );
166 
167  moab::EntityHandle rootset = 0;
168  ///////////////////////////////////////////////////////////////////////////
169  // The metadata in H5M file contains the following data:
170  //
171  // 1. n_a: Total source entities: (number of elements in source mesh)
172  // 2. n_b: Total target entities: (number of elements in target mesh)
173  // 3. nv_a: Max edge size of elements in source mesh
174  // 4. nv_b: Max edge size of elements in target mesh
175  // 5. maxrows: Number of rows in remap weight matrix
176  // 6. maxcols: Number of cols in remap weight matrix
177  // 7. nnz: Number of total nnz in sparse remap weight matrix
178  // 8. np_a: The order of the field description on the source mesh: >= 1
179  // 9. np_b: The order of the field description on the target mesh: >= 1
180  // 10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2
181  // = dGLL]
182  // 11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2
183  // = dGLL]
184  // 12. conserved: Flag to specify whether the remap operator has conservation constraints:
185  // [0, 1]
186  // 13. monotonicity: Flags to specify whether the remap operator has monotonicity
187  // constraints: [0, 1, 2]
188  //
189  ///////////////////////////////////////////////////////////////////////////
190  Tag smatMetadataTag;
191  int smat_metadata_glb[13];
192  MB_CHK_ERR( mbCore->tag_get_handle( "SMAT_DATA", 13, MB_TYPE_INTEGER, smatMetadataTag, MB_TAG_SPARSE ) );
193  MB_CHK_ERR( mbCore->tag_get_data( smatMetadataTag, &rootset, 1, smat_metadata_glb ) );
194  // std::cout << "Number of mesh sets is " << meshsets.size() << std::endl;
195 
196 #define DTYPE( a ) \
197  { \
198  ( ( ( a ) == 0 ) ? "FV" : ( ( ( a ) == 1 ) ? "cGLL" : "dGLL" ) ) \
199  }
200  // Map dimensions
201  int nA = smat_metadata_glb[0];
202  int nB = smat_metadata_glb[1];
203  int nVA = smat_metadata_glb[2];
204  int nVB = smat_metadata_glb[3];
205  int nDofB = smat_metadata_glb[4];
206  int nDofA = smat_metadata_glb[5];
207  int NNZ = smat_metadata_glb[6];
208  int nOrdA = smat_metadata_glb[7];
209  int nOrdB = smat_metadata_glb[8];
210  int nBasA = smat_metadata_glb[9];
211  std::string methodA = DTYPE( nBasA );
212  int nBasB = smat_metadata_glb[10];
213  std::string methodB = DTYPE( nBasB );
214  int bConserved = smat_metadata_glb[11];
215  int bMonotonicity = smat_metadata_glb[12];
216 
217  EntityHandle source_mesh = 0, target_mesh = 0, overlap_mesh = 0;
218  for( unsigned im = 0; im < meshsets.size(); ++im )
219  {
220  moab::Range elems;
221  MB_CHK_ERR( mbCore->get_entities_by_dimension( meshsets[im], 2, elems ) );
222  if( elems.size() - nA == 0 && source_mesh == 0 )
223  source_mesh = meshsets[im];
224  else if( elems.size() - nB == 0 && target_mesh == 0 )
225  target_mesh = meshsets[im];
226  else if( overlap_mesh == 0 )
227  overlap_mesh = meshsets[im];
228  else
229  continue;
230  }
231 
232  Tag srcIDTag, srcAreaTag, tgtIDTag, tgtAreaTag;
233  MB_CHK_ERR( mbCore->tag_get_handle( "SourceGIDS", srcIDTag ) );
234  MB_CHK_ERR( mbCore->tag_get_handle( "SourceAreas", srcAreaTag ) );
235  MB_CHK_ERR( mbCore->tag_get_handle( "TargetGIDS", tgtIDTag ) );
236  MB_CHK_ERR( mbCore->tag_get_handle( "TargetAreas", tgtAreaTag ) );
237  Tag smatRowdataTag, smatColdataTag, smatValsdataTag;
238  MB_CHK_ERR( mbCore->tag_get_handle( "SMAT_ROWS", smatRowdataTag ) );
239  MB_CHK_ERR( mbCore->tag_get_handle( "SMAT_COLS", smatColdataTag ) );
240  MB_CHK_ERR( mbCore->tag_get_handle( "SMAT_VALS", smatValsdataTag ) );
241  Tag srcCenterLon, srcCenterLat, tgtCenterLon, tgtCenterLat;
242  MB_CHK_ERR( mbCore->tag_get_handle( "SourceCoordCenterLon", srcCenterLon ) );
243  MB_CHK_ERR( mbCore->tag_get_handle( "SourceCoordCenterLat", srcCenterLat ) );
244  MB_CHK_ERR( mbCore->tag_get_handle( "TargetCoordCenterLon", tgtCenterLon ) );
245  MB_CHK_ERR( mbCore->tag_get_handle( "TargetCoordCenterLat", tgtCenterLat ) );
246  Tag srcVertexLon, srcVertexLat, tgtVertexLon, tgtVertexLat;
247  MB_CHK_ERR( mbCore->tag_get_handle( "SourceCoordVertexLon", srcVertexLon ) );
248  MB_CHK_ERR( mbCore->tag_get_handle( "SourceCoordVertexLat", srcVertexLat ) );
249  MB_CHK_ERR( mbCore->tag_get_handle( "TargetCoordVertexLon", tgtVertexLon ) );
250  MB_CHK_ERR( mbCore->tag_get_handle( "TargetCoordVertexLat", tgtVertexLat ) );
251 
252  // Get sets entities, by type
253  moab::Range sets;
254  // MB_CHK_ERR( mbCore->get_entities_by_type(0, MBENTITYSET, sets) );
255  MB_CHK_ERR( mbCore->get_entities_by_type_and_tag( 0, MBENTITYSET, &smatRowdataTag, NULL, 1, sets,
256  moab::Interface::UNION, true ) );
257 
258  std::vector< int > src_gids, tgt_gids;
259  std::vector< double > src_areas, tgt_areas;
260  int srcID_size, tgtID_size, srcArea_size, tgtArea_size;
261  MB_CHK_SET_ERR( get_vartag_data( mbCore, srcIDTag, sets, srcID_size, src_gids ),
262  "Getting source mesh IDs failed" );
263  MB_CHK_SET_ERR( get_vartag_data( mbCore, tgtIDTag, sets, tgtID_size, tgt_gids ),
264  "Getting target mesh IDs failed" );
265  MB_CHK_SET_ERR( get_vartag_data( mbCore, srcAreaTag, sets, srcArea_size, src_areas ),
266  "Getting source mesh areas failed" );
267  MB_CHK_SET_ERR( get_vartag_data( mbCore, tgtAreaTag, sets, tgtArea_size, tgt_areas ),
268  "Getting target mesh areas failed" );
269 
270  assert( srcArea_size == srcID_size );
271  assert( tgtArea_size == tgtID_size );
272 
273  std::vector< double > src_glob_areas( nDofA, 0.0 ), tgt_glob_areas( nDofB, 0.0 );
274  for( int i = 0; i < srcArea_size; ++i )
275  {
276  // printf("%d/%d: %d = Found ID %d and area %5.6e\n", i, srcArea_size, nDofA,
277  // src_gids[i], src_areas[i]);
278  assert( i < srcID_size );
279  assert( src_gids[i] < nDofA );
280  if( src_areas[i] > src_glob_areas[src_gids[i]] ) src_glob_areas[src_gids[i]] = src_areas[i];
281  }
282  for( int i = 0; i < tgtArea_size; ++i )
283  {
284  // printf("%d/%d: %d = Found ID %d and area %5.6e\n", i, tgtArea_size, nDofB,
285  // tgt_gids[i], tgt_areas[i]);
286  assert( i < tgtID_size );
287  assert( tgt_gids[i] < nDofB );
288  if( tgt_areas[i] > tgt_glob_areas[tgt_gids[i]] ) tgt_glob_areas[tgt_gids[i]] = tgt_areas[i];
289  }
290 
291  // Write output dimensions entries
292  int nSrcGridDims = 1;
293  int nDstGridDims = 1;
294 
295  NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
296  NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );
297 
298  NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
299  NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );
300 
301  if( nA == nDofA )
302  {
303  varSrcGridDims->put( &nA, 1 );
304  varSrcGridDims->add_att( "name0", "num_elem" );
305  }
306  else
307  {
308  varSrcGridDims->put( &nDofA, 1 );
309  varSrcGridDims->add_att( "name1", "num_dof" );
310  }
311 
312  if( nB == nDofB )
313  {
314  varDstGridDims->put( &nB, 1 );
315  varDstGridDims->add_att( "name0", "num_elem" );
316  }
317  else
318  {
319  varDstGridDims->put( &nDofB, 1 );
320  varDstGridDims->add_att( "name1", "num_dof" );
321  }
322 
323  // Source and Target mesh resolutions
324  NcDim* dimNA = ncMap.add_dim( "n_a", nDofA );
325  NcDim* dimNB = ncMap.add_dim( "n_b", nDofB );
326 
327  // Source and Target verticecs per elements
328  const int nva = ( nA == nDofA ? nVA : 1 );
329  const int nvb = ( nB == nDofB ? nVB : 1 );
330  NcDim* dimNVA = ncMap.add_dim( "nv_a", nva );
331  NcDim* dimNVB = ncMap.add_dim( "nv_b", nvb );
332 
333  // Source and Target verticecs per elements
334  // NcDim * dimNEA = ncMap.add_dim("ne_a", nA);
335  // NcDim * dimNEB = ncMap.add_dim("ne_b", nB);
336 
337  if( writeXYCoords )
338  {
339  // Write coordinates
340  NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA /*dimNA*/ );
341  NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB /*dimNB*/ );
342 
343  NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA /*dimNA*/ );
344  NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB /*dimNB*/ );
345 
346  NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA /*dimNA*/, dimNVA );
347  NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB /*dimNB*/, dimNVB );
348 
349  NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA /*dimNA*/, dimNVA );
350  NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB /*dimNB*/, dimNVB );
351 
352  varYCA->add_att( "units", "degrees" );
353  varYCB->add_att( "units", "degrees" );
354 
355  varXCA->add_att( "units", "degrees" );
356  varXCB->add_att( "units", "degrees" );
357 
358  varYVA->add_att( "units", "degrees" );
359  varYVB->add_att( "units", "degrees" );
360 
361  varXVA->add_att( "units", "degrees" );
362  varXVB->add_att( "units", "degrees" );
363 
364  std::vector< double > src_centerlat, src_centerlon;
365  int srccenter_size;
366  MB_CHK_SET_ERR( get_vartag_data( mbCore, srcCenterLat, sets, srccenter_size, src_centerlat ),
367  "Getting source mesh areas failed" );
368  MB_CHK_SET_ERR( get_vartag_data( mbCore, srcCenterLon, sets, srccenter_size, src_centerlon ),
369  "Getting target mesh areas failed" );
370  std::vector< double > src_glob_centerlat( nDofA, 0.0 ), src_glob_centerlon( nDofA, 0.0 );
371 
372  for( int i = 0; i < srccenter_size; ++i )
373  {
374  assert( i < srcID_size );
375  assert( src_gids[i] < nDofA );
376 
377  src_glob_centerlat[src_gids[i]] = src_centerlat[i];
378  src_glob_centerlon[src_gids[i]] = src_centerlon[i];
379  }
380 
381  std::vector< double > tgt_centerlat, tgt_centerlon;
382  int tgtcenter_size;
383  MB_CHK_SET_ERR( get_vartag_data( mbCore, tgtCenterLat, sets, tgtcenter_size, tgt_centerlat ),
384  "Getting source mesh areas failed" );
385  MB_CHK_SET_ERR( get_vartag_data( mbCore, tgtCenterLon, sets, tgtcenter_size, tgt_centerlon ),
386  "Getting target mesh areas failed" );
387  std::vector< double > tgt_glob_centerlat( nDofB, 0.0 ), tgt_glob_centerlon( nDofB, 0.0 );
388  for( int i = 0; i < tgtcenter_size; ++i )
389  {
390  assert( i < tgtID_size );
391  assert( tgt_gids[i] < nDofB );
392 
393  tgt_glob_centerlat[tgt_gids[i]] = tgt_centerlat[i];
394  tgt_glob_centerlon[tgt_gids[i]] = tgt_centerlon[i];
395  }
396 
397  varYCA->put( &( src_glob_centerlat[0] ), nDofA );
398  varYCB->put( &( tgt_glob_centerlat[0] ), nDofB );
399  varXCA->put( &( src_glob_centerlon[0] ), nDofA );
400  varXCB->put( &( tgt_glob_centerlon[0] ), nDofB );
401 
402  src_centerlat.clear();
403  src_centerlon.clear();
404  tgt_centerlat.clear();
405  tgt_centerlon.clear();
406 
407  DataArray2D< double > src_glob_vertexlat( nDofA, nva ), src_glob_vertexlon( nDofA, nva );
408  if( nva > 1 )
409  {
410  std::vector< double > src_vertexlat, src_vertexlon;
411  int srcvertex_size;
412  MB_CHK_SET_ERR( get_vartag_data( mbCore, srcVertexLat, sets, srcvertex_size, src_vertexlat ),
413  "Getting source mesh areas failed" );
414  MB_CHK_SET_ERR( get_vartag_data( mbCore, srcVertexLon, sets, srcvertex_size, src_vertexlon ),
415  "Getting target mesh areas failed" );
416  int offset = 0;
417  for( unsigned vIndex = 0; vIndex < src_gids.size(); ++vIndex )
418  {
419  for( int vNV = 0; vNV < nva; ++vNV )
420  {
421  assert( offset < srcvertex_size );
422  src_glob_vertexlat[src_gids[vIndex]][vNV] = src_vertexlat[offset];
423  src_glob_vertexlon[src_gids[vIndex]][vNV] = src_vertexlon[offset];
424  offset++;
425  }
426  }
427  }
428 
429  DataArray2D< double > tgt_glob_vertexlat( nDofB, nvb ), tgt_glob_vertexlon( nDofB, nvb );
430  if( nvb > 1 )
431  {
432  std::vector< double > tgt_vertexlat, tgt_vertexlon;
433  int tgtvertex_size;
434  MB_CHK_SET_ERR( get_vartag_data( mbCore, tgtVertexLat, sets, tgtvertex_size, tgt_vertexlat ),
435  "Getting source mesh areas failed" );
436  MB_CHK_SET_ERR( get_vartag_data( mbCore, tgtVertexLon, sets, tgtvertex_size, tgt_vertexlon ),
437  "Getting target mesh areas failed" );
438  int offset = 0;
439  for( unsigned vIndex = 0; vIndex < tgt_gids.size(); ++vIndex )
440  {
441  for( int vNV = 0; vNV < nvb; ++vNV )
442  {
443  assert( offset < tgtvertex_size );
444  tgt_glob_vertexlat[tgt_gids[vIndex]][vNV] = tgt_vertexlat[offset];
445  tgt_glob_vertexlon[tgt_gids[vIndex]][vNV] = tgt_vertexlon[offset];
446  offset++;
447  }
448  }
449  }
450 
451  varYVA->put( &( src_glob_vertexlat[0][0] ), nDofA, nva );
452  varYVB->put( &( tgt_glob_vertexlat[0][0] ), nDofB, nvb );
453 
454  varXVA->put( &( src_glob_vertexlon[0][0] ), nDofA, nva );
455  varXVB->put( &( tgt_glob_vertexlon[0][0] ), nDofB, nvb );
456  }
457 
458  // Write areas
459  NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
460  varAreaA->put( &( src_glob_areas[0] ), nDofA );
461  // varAreaA->add_att("units", "steradians");
462 
463  NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
464  varAreaB->put( &( tgt_glob_areas[0] ), nDofB );
465  // varAreaB->add_att("units", "steradians");
466 
467  std::vector< int > mat_rows, mat_cols;
468  std::vector< double > mat_vals;
469  int row_sizes, col_sizes, val_sizes;
470  MB_CHK_SET_ERR( get_vartag_data( mbCore, smatRowdataTag, sets, row_sizes, mat_rows ),
471  "Getting matrix row data failed" );
472  assert( row_sizes == NNZ );
473  MB_CHK_SET_ERR( get_vartag_data( mbCore, smatColdataTag, sets, col_sizes, mat_cols ),
474  "Getting matrix col data failed" );
475  assert( col_sizes == NNZ );
476  MB_CHK_SET_ERR( get_vartag_data( mbCore, smatValsdataTag, sets, val_sizes, mat_vals ),
477  "Getting matrix values failed" );
478  assert( val_sizes == NNZ );
479 
480  // Let us form the matrix in-memory and consolidate shared DoF rows from shared-process
481  // contributions
482  SparseMatrix< double > mapMatrix;
483 
484  for( int innz = 0; innz < NNZ; ++innz )
485  {
486 #ifdef VERBOSE
487  if( fabs( mapMatrix( mat_rows[innz], mat_cols[innz] ) ) > 1e-12 )
488  {
489  printf( "Adding to existing loc: (%d, %d) = %12.8f\n", mat_rows[innz], mat_cols[innz],
490  mapMatrix( mat_rows[innz], mat_cols[innz] ) );
491  }
492 #endif
493  mapMatrix( mat_rows[innz], mat_cols[innz] ) += mat_vals[innz];
494  }
495 
496  // Write SparseMatrix entries
497  DataArray1D< int > vecRow;
498  DataArray1D< int > vecCol;
499  DataArray1D< double > vecS;
500 
501  mapMatrix.GetEntries( vecRow, vecCol, vecS );
502 
503  int nS = vecS.GetRows();
504 
505  // Print more information about what we are converting:
506  // Source elements/vertices/type (Discretization ?)
507  // Target elements/vertices/type (Discretization ?)
508  // Overlap elements/types
509  // Rmeapping weights matrix: rows/cols/NNZ
510  // Output the number of sets
511  printf( "Primary sets: %15zu\n", sets.size() );
512  printf( "Original NNZ: %18d\n", NNZ );
513  printf( "Consolidated Total NNZ: %8d\n", nS );
514  printf( "Conservative weights ? %6d\n", ( bConserved > 0 ) );
515  printf( "Monotone weights ? %10d\n", ( bMonotonicity > 0 ) );
516 
517  printf( "\n--------------------------------------------------------------\n" );
518  printf( "%20s %21s %15s\n", "Description", "Source", "Target" );
519  printf( "--------------------------------------------------------------\n" );
520 
521  printf( "%25s %15d %15d\n", "Number of elements:", nA, nB );
522  printf( "%25s %15d %15d\n", "Number of DoFs:", nDofA, nDofB );
523  printf( "%25s %15d %15d\n", "Maximum vertex/element:", nVA, nVB );
524  printf( "%25s %15s %15s\n", "Discretization type:", methodA.c_str(), methodB.c_str() );
525  printf( "%25s %15d %15d\n", "Discretization order:", nOrdA, nOrdB );
526 
527  // Calculate and write fractional coverage arrays
528  {
529  DataArray1D< double > dFracA( nDofA );
530  DataArray1D< double > dFracB( nDofB );
531 
532  for( int i = 0; i < nS; i++ )
533  {
534  // std::cout << i << " - mat_vals = " << mat_vals[i] << " dFracA = " << mat_vals[i]
535  // / src_glob_areas[vecCol[i]] * tgt_glob_areas[vecRow[i]] << std::endl;
536  dFracA[vecCol[i]] += vecS[i] / src_glob_areas[vecCol[i]] * tgt_glob_areas[vecRow[i]];
537  dFracB[vecRow[i]] += vecS[i];
538  }
539 
540  NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
541  varFracA->put( &( dFracA[0] ), nDofA );
542  varFracA->add_att( "name", "fraction of target coverage of source dof" );
543  varFracA->add_att( "units", "unitless" );
544 
545  NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
546  varFracB->put( &( dFracB[0] ), nDofB );
547  varFracB->add_att( "name", "fraction of source coverage of target dof" );
548  varFracB->add_att( "units", "unitless" );
549  }
550 
551  // Write out data
552  NcDim* dimNS = ncMap.add_dim( "n_s", nS );
553 
554  NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
555  varRow->add_att( "name", "sparse matrix target dof index" );
556  varRow->add_att( "first_index", "1" );
557 
558  NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
559  varCol->add_att( "name", "sparse matrix source dof index" );
560  varCol->add_att( "first_index", "1" );
561 
562  NcVar* varS = ncMap.add_var( "S", ncDouble, dimNS );
563  varS->add_att( "name", "sparse matrix coefficient" );
564 
565  // Increment vecRow and vecCol: make it 1-based
566  for( int i = 0; i < nS; i++ )
567  {
568  vecRow[i]++;
569  vecCol[i]++;
570  }
571 
572  varRow->set_cur( (long)0 );
573  varRow->put( &( vecRow[0] ), nS );
574 
575  varCol->set_cur( (long)0 );
576  varCol->put( &( vecCol[0] ), nS );
577 
578  varS->set_cur( (long)0 );
579  varS->put( &( vecS[0] ), nS );
580 
581  ncMap.close();
582 
583  // MB_CHK_ERR( mbCore->write_file(scripfile.c_str()) );
584  }
585  catch( std::exception& e )
586  {
587  std::cout << " exception caught during tree initialization " << e.what() << std::endl;
588  }
589  delete mbCore;
590 
591 #ifdef MOAB_HAVE_MPI
592  MPI_Finalize();
593 #endif
594 
595  exit( 0 );
596 }