 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 }
1.9.1.