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