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