164 MPI_Init( &argc, &argv );
165 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
168 #ifdef MOAB_HAVE_TEMPESTREMAP
169 bool tempestin =
false, tempestout =
false;
177 #if( defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_TEMPESTREMAP ) )
181 bool append_rank =
false;
182 bool percent_rank_subst =
false;
183 bool file_written =
false;
185 std::list< std::string >::iterator j;
186 bool dims[4] = {
false,
false,
false,
false };
187 const char* format = NULL;
188 std::list< std::string > in;
191 std::set< int > geom[4], mesh[4];
192 std::vector< EntityHandle > set_list;
193 std::vector< std::string > write_opts, read_opts;
194 std::string metis_partition_file;
195 #ifdef MOAB_HAVE_TEMPESTREMAP
196 std::string globalid_tag_name;
197 int spectral_order = 1;
198 bool unitscaling =
false;
203 const char*
const geom_names[] = {
"VERTEX",
"CURVE",
"SURFACE",
"VOLUME" };
207 bool print_times =
false;
208 bool generate[] = {
false,
false,
false };
210 bool parallel =
false, resolve_shared =
false, exchange_ghosts =
false;
211 for( i = 1; i < argc; i++ )
215 if( do_flag && argv[i][0] ==
'-' )
217 if( !argv[i][1] || ( argv[i][1] !=
'M' && argv[i][2] ) )
usage_error( argv[0] );
245 percent_rank_subst =
true;
249 if( argv[i][2] ==
'1' || argv[i][2] ==
'2' ) resolve_shared =
true;
250 if( argv[i][2] ==
'2' ) exchange_ghosts =
true;
253 #ifdef MOAB_HAVE_TEMPESTREMAP
267 dims[argv[i][1] -
'0'] =
true;
272 if( i == argc || argv[i][0] ==
'-' )
274 std::cerr <<
"Expected argument following " << argv[i - 1] << std::endl;
277 if( argv[i - 1][1] ==
'I' )
279 dim = atoi( argv[i] );
280 if( dim < 1 || dim > 2 )
282 std::cerr <<
"Invalid dimension value following -I" << std::endl;
285 generate[
dim] =
true;
289 switch( argv[i - 1][1] )
292 read_opts.push_back( std::string(
"SAT_FILE=" ) + argv[i] );
300 write_opts.push_back( argv[i] );
304 read_opts.push_back( argv[i] );
307 #ifdef MOAB_HAVE_TEMPESTREMAP
309 globalid_tag_name = std::string( argv[i] );
313 spectral_order = atoi( argv[i] );
342 metis_partition_file = argv[i];
346 std::cerr <<
"Invalid option: " << argv[i] << std::endl;
351 std::cerr <<
"Invalid flag or flag value: " << argv[i - 1] <<
" " << argv[i] << std::endl;
359 in.push_back( argv[i] );
364 std::cerr <<
"No output file name specified." << std::endl;
373 std::ostringstream mod;
374 mod << out <<
"." << proc_id;
378 if( percent_rank_subst )
380 for( j = in.begin(); j != in.end(); ++j )
386 std::string read_options, write_options;
389 read_opts.push_back(
"PARALLEL=READ_PART" );
390 read_opts.push_back(
"PARTITION=PARALLEL_PARTITION" );
391 if( !append_rank && !percent_rank_subst ) write_opts.push_back(
"PARALLEL=WRITE_PART" );
393 if( resolve_shared ) read_opts.push_back(
"PARALLEL_RESOLVE_SHARED_ENTS" );
394 if( exchange_ghosts ) read_opts.push_back(
"PARALLEL_GHOSTS=3.0.1" );
404 if( !metis_partition_file.empty() )
406 if( ( in.size() != 1 ) || ( proc_id != 0 ) )
408 std::cerr <<
" mpas partition allows only one input file, in serial conversion\n";
419 #ifdef MOAB_HAVE_TEMPESTREMAP
420 if( tempestin && in.size() > 1 )
422 std::cerr <<
" we can read only one tempest files at a time\n";
429 if( tempestin or tempestout )
438 bool use_overlap_context =
false;
439 Tag srcParentTag, tgtParentTag;
442 for( j = in.begin(); j != in.end(); ++j )
444 std::string inFileName = *j;
448 #ifdef MOAB_HAVE_TEMPESTREMAP
463 tempestMesh->RemoveZeroEdges();
464 tempestMesh->RemoveCoincidentNodes();
475 NcFile ncInput( inFileName.c_str(), NcFile::ReadOnly );
477 NcError error_temp( NcError::silent_nonfatal );
479 NcAtt* attRectilinear = ncInput.get_att(
"rectilinear" );
482 std::vector< int > vecDimSizes( 3, 0 );
490 if( attRectilinear !=
nullptr )
493 NcAtt* attRectilinearDim0Size = ncInput.get_att(
"rectilinear_dim0_size" );
494 NcAtt* attRectilinearDim1Size = ncInput.get_att(
"rectilinear_dim1_size" );
496 if( attRectilinearDim0Size ==
nullptr )
498 _EXCEPTIONT(
"Missing attribute \"rectilinear_dim0_size\"" );
500 if( attRectilinearDim1Size ==
nullptr )
502 _EXCEPTIONT(
"Missing attribute \"rectilinear_dim1_size\"" );
505 int nDim0Size = attRectilinearDim0Size->as_int( 0 );
506 int nDim1Size = attRectilinearDim1Size->as_int( 0 );
509 NcAtt* attRectilinearDim0Name = ncInput.get_att(
"rectilinear_dim0_name" );
510 NcAtt* attRectilinearDim1Name = ncInput.get_att(
"rectilinear_dim1_name" );
512 if( attRectilinearDim0Name ==
nullptr )
514 _EXCEPTIONT(
"Missing attribute \"rectilinear_dim0_name\"" );
516 if( attRectilinearDim1Name ==
nullptr )
518 _EXCEPTIONT(
"Missing attribute \"rectilinear_dim1_name\"" );
521 std::string strDim0Name = attRectilinearDim0Name->as_string( 0 );
522 std::string strDim1Name = attRectilinearDim1Name->as_string( 0 );
524 std::map< std::string, int > vecDimNameSizes;
526 vecDimNameSizes[strDim0Name] = nDim0Size;
527 vecDimNameSizes[strDim1Name] = nDim1Size;
529 vecDimSizes[1] = vecDimNameSizes[
"lat"];
530 vecDimSizes[2] = vecDimNameSizes[
"lon"];
532 printf(
"Rectilinear RLL mesh size: (lat) %d X (lon) %d\n", vecDimSizes[1], vecDimSizes[2] );
547 vecDimSizes[1] = elems.
size();
550 switch( vecDimSizes[0] )
553 printf(
"Cubed-Sphere mesh: %d (elems), %d (nodes)\n", vecDimSizes[1], vecDimSizes[2] );
556 printf(
"Icosahedral (triangular) mesh: %d (elems), %d (nodes)\n", vecDimSizes[1],
561 printf(
"Polygonal mesh: %d (elems), %d (nodes)\n", vecDimSizes[1], vecDimSizes[2] );
572 const size_t nOverlapFaces = tempestMesh->faces.size();
573 if( tempestMesh->vecSourceFaceIx.size() == nOverlapFaces &&
574 tempestMesh->vecSourceFaceIx.size() == nOverlapFaces )
577 use_overlap_context =
true;
588 std::vector< int > gids( faces.
size() ), srcpar( faces.
size() ), tgtpar( faces.
size() );
591 for(
unsigned ii = 0; ii < faces.
size(); ++ii )
593 srcpar[ii] = tempestMesh->vecSourceFaceIx[gids[ii] - 1];
594 tgtpar[ii] = tempestMesh->vecTargetFaceIx[gids[ii] - 1];
605 else if( tempestout )
611 std::vector< int > metadata( 2 );
625 use_overlap_context =
true;
634 result = reorder_tool.handle_order_from_int_tag( srcParentTag, -1, order );
MB_CHK_ERR( result );
635 result = reorder_tool.reorder_entities( order );
MB_CHK_ERR( result );
643 assert( metadata.size() );
644 std::cout <<
"Converting a RLL mesh with rectilinear dimension: " << metadata[0] <<
" X "
645 << metadata[1] << std::endl;
654 result =
gMB->
load_file( j->c_str(), 0, read_options.c_str() );
656 result =
gMB->
load_file( j->c_str(), 0, read_options.c_str() );
660 std::cerr <<
"Failed to load \"" << *j <<
"\"." << std::endl;
661 std::cerr <<
"Error code: " <<
gMB->
get_error_string( result ) <<
" (" << result <<
")" << std::endl;
664 std::cerr <<
"Error message: " << message << std::endl;
670 if( !proc_id ) std::cerr <<
"Read \"" << *j <<
"\"" << std::endl;
671 if( print_times && !proc_id )
write_times( std::cout );
675 bool have_geom =
false;
678 if( !geom[
dim].empty() ) have_geom =
true;
683 bool have_sets = have_geom;
691 std::cerr <<
"No ID tag defined." << std::endl;
697 std::cerr <<
"No geometry tag defined." << std::endl;
706 Tag tags[] = { id_tag, dim_tag };
707 const void* vals[] = { &id_val, &
dim };
710 int init_count = set_list.size();
711 for( std::set< int >::iterator iter = geom[
dim].begin(); iter != geom[
dim].end(); ++iter )
719 std::cerr << geom_names[
dim] <<
" " << id_val <<
" not found.\n";
721 std::copy( range.
begin(), range.
end(), std::back_inserter( set_list ) );
725 std::cout <<
"Found " << ( set_list.size() - init_count ) <<
' ' << geom_names[
dim] <<
" sets"
731 for( i = 0; i < 4; ++i )
735 if( mesh[i].empty() )
continue;
743 std::cerr <<
"Tag not found: " << mesh_tag_names[i] << std::endl;
748 int init_count = set_list.size();
749 for( std::set< int >::iterator iter = mesh[i].begin(); iter != mesh[i].end(); ++iter )
752 const void* vals[] = { &*iter };
757 std::cerr << mesh_tag_names[i] <<
" " << *iter <<
" not found.\n";
759 std::copy( range.
begin(), range.
end(), std::back_inserter( set_list ) );
763 std::cout <<
"Found " << ( set_list.size() - init_count ) <<
' ' << mesh_tag_names[i] <<
" sets"
770 if( dims[
dim] ) bydim =
true;
775 if( generate[1] && !dims[1] )
777 std::cerr <<
"Warning: Request to generate 1D internal entities but not export them." << std::endl;
780 if( generate[2] && !dims[2] )
782 std::cerr <<
"Warning: Request to generate 2D internal entities but not export them." << std::endl;
788 if( generate[1] || generate[2] )
795 num_sets = set_list.size();
798 for( i = 0; i < num_sets; ++i )
800 Range dim3, dim2, adj;
820 Range dead_entities, tmp_range;
823 if( dims[
dim] )
continue;
825 dead_entities.
merge( tmp_range );
831 while( !empty_sets.
empty() )
834 dead_entities.
merge( empty_sets );
837 empty_sets =
subtract( tmp_range, dead_entities );
844 if( have_sets && set_list.empty() )
846 std::cerr <<
"Nothing to write." << std::endl;
854 if( !metis_partition_file.empty() )
859 std::cerr <<
"Failed to process partition file \"" << metis_partition_file <<
"\"." << std::endl;
869 std::cout <<
"Found " << set_list.size() <<
" specified sets to write (total)." << std::endl;
871 std::cout <<
"No sets specifed. Writing entire mesh." << std::endl;
876 #ifdef MOAB_HAVE_TEMPESTREMAP
885 int ntot_elements = 0, nelements = faces.
size();
887 int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, pcomm->
comm() );
888 if( ierr != 0 )
MB_CHK_SET_ERR( MB_FAILURE,
"MPI_Allreduce failed to get total source elements" );
890 ntot_elements = nelements;
894 std::vector< int > gids( faces.
size() );
897 if( faces.
size() > 1 && gids[0] == gids[1] && !use_overlap_context )
910 if( spectral_order > 1 && globalid_tag_name.size() > 1 )
912 result = remapper->
GenerateMeshMetadata( *tempestMesh, ntot_elements, faces, NULL, globalid_tag_name,
920 if( use_overlap_context &&
false )
922 const int nOverlapFaces = faces.
size();
924 tempestMesh->vecSourceFaceIx.resize( nOverlapFaces );
925 tempestMesh->vecTargetFaceIx.resize( nOverlapFaces );
930 tempestMesh->Write( out, NcFile::Netcdf4 );
940 result =
gMB->
write_file( out.c_str(), format, write_options.c_str(), &set_list[0], set_list.size() );
942 result =
gMB->
write_file( out.c_str(), format, write_options.c_str() );
945 std::cerr <<
"Failed to write \"" << out <<
"\"." << std::endl;
946 std::cerr <<
"Error code: " <<
gMB->
get_error_string( result ) <<
" (" << result <<
")" << std::endl;
949 std::cerr <<
"Error message: " << message << std::endl;
957 if( !proc_id ) std::cerr <<
"Wrote \"" << out <<
"\"" << std::endl;
958 if( print_times && !proc_id )
write_times( std::cout );