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