172 MPI_Init( &argc, &argv );
173 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
176 #ifdef MOAB_HAVE_TEMPESTREMAP
177 bool tempestin =
false, tempestout =
false;
185 #if ( defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_TEMPESTREMAP ) )
189 bool append_rank =
false;
190 bool percent_rank_subst =
false;
191 bool file_written =
false;
193 std::list< std::string >::iterator j;
194 bool dims[4] = {
false,
false,
false,
false };
195 const char* format = NULL;
196 std::list< std::string > in;
199 std::set< int > geom[4], mesh[4];
200 std::vector< EntityHandle > set_list;
201 std::vector< std::string > write_opts, read_opts;
202 std::string metis_partition_file;
203 #ifdef MOAB_HAVE_TEMPESTREMAP
204 std::string globalid_tag_name;
205 int spectral_order = 1;
206 bool unitscaling =
false;
211 const char*
const geom_names[] = {
"VERTEX",
"CURVE",
"SURFACE",
"VOLUME" };
215 bool print_times =
false;
216 bool generate[] = {
false,
false,
false };
218 bool parallel =
false, resolve_shared =
false, exchange_ghosts =
false;
219 for( i = 1; i < argc; i++ )
223 if( do_flag && argv[i][0] ==
'-' )
225 if( !argv[i][1] || ( argv[i][1] !=
'M' && argv[i][2] ) )
usage_error( argv[0] );
253 percent_rank_subst =
true;
257 if( argv[i][2] ==
'1' || argv[i][2] ==
'2' ) resolve_shared =
true;
258 if( argv[i][2] ==
'2' ) exchange_ghosts =
true;
261 #ifdef MOAB_HAVE_TEMPESTREMAP
275 dims[argv[i][1] -
'0'] =
true;
280 if( i == argc || argv[i][0] ==
'-' )
282 std::cerr <<
"Expected argument following " << argv[i - 1] << std::endl;
285 if( argv[i - 1][1] ==
'I' )
287 dim = atoi( argv[i] );
288 if( dim < 1 || dim > 2 )
290 std::cerr <<
"Invalid dimension value following -I" << std::endl;
293 generate[dim] =
true;
297 switch( argv[i - 1][1] )
300 read_opts.push_back( std::string(
"SAT_FILE=" ) + argv[i] );
308 write_opts.push_back( argv[i] );
312 read_opts.push_back( argv[i] );
315 #ifdef MOAB_HAVE_TEMPESTREMAP
317 globalid_tag_name = std::string( argv[i] );
321 spectral_order = atoi( argv[i] );
350 metis_partition_file = argv[i];
354 std::cerr <<
"Invalid option: " << argv[i] << std::endl;
359 std::cerr <<
"Invalid flag or flag value: " << argv[i - 1] <<
" " << argv[i] << std::endl;
367 in.push_back( argv[i] );
372 std::cerr <<
"No output file name specified." << std::endl;
381 std::ostringstream mod;
382 mod << out <<
"." << proc_id;
386 if( percent_rank_subst )
388 for( j = in.begin(); j != in.end(); ++j )
394 std::string read_options, write_options;
397 read_opts.push_back(
"PARALLEL=READ_PART" );
398 read_opts.push_back(
"PARTITION=PARALLEL_PARTITION" );
399 if( !append_rank && !percent_rank_subst ) write_opts.push_back(
"PARALLEL=WRITE_PART" );
401 if( resolve_shared ) read_opts.push_back(
"PARALLEL_RESOLVE_SHARED_ENTS" );
402 if( exchange_ghosts ) read_opts.push_back(
"PARALLEL_GHOSTS=3.0.1" );
412 if( !metis_partition_file.empty() )
414 if( ( in.size() != 1 ) || ( proc_id != 0 ) )
416 std::cerr <<
" mpas partition allows only one input file, in serial conversion\n";
427 #ifdef MOAB_HAVE_TEMPESTREMAP
428 if( tempestin && in.size() > 1 )
430 std::cerr <<
" we can read only one tempest files at a time\n";
437 if( tempestin or tempestout )
446 bool use_overlap_context =
false;
447 Tag srcParentTag, tgtParentTag;
463 std::cerr <<
"Failed to create file set." << std::endl;
470 for( j = in.begin(); j != in.end(); ++j )
472 std::string inFileName = *j;
476 #ifdef MOAB_HAVE_TEMPESTREMAP
485 #ifdef MOAB_HAVE_NETCDF
493 tempestMesh->RemoveZeroEdges();
494 tempestMesh->RemoveCoincidentNodes();
498 MB_CHK_SET_ERR( result,
"can't convert the TempestRemap mesh to MOAB" );
507 MB_CHK_SET_ERR( result,
"can't get entities from the TempestRemap source mesh" );
509 MB_CHK_SET_ERR( result,
"can't add the TempestRemap mesh to the file set" );
519 NcFile ncInput( inFileName.c_str(), NcFile::ReadOnly );
521 NcError error_temp( NcError::silent_nonfatal );
523 NcAtt* attRectilinear = ncInput.get_att(
"rectilinear" );
524 NcVar* varGridDims = ncInput.get_var(
"grid_dims" );
527 std::vector< int > vecDimSizes( 3, 0 );
536 if( attRectilinear !=
nullptr )
539 NcAtt* attRectilinearDim0Size = ncInput.get_att(
"rectilinear_dim0_size" );
540 NcAtt* attRectilinearDim1Size = ncInput.get_att(
"rectilinear_dim1_size" );
542 if( attRectilinearDim0Size ==
nullptr )
544 _EXCEPTIONT(
"Missing attribute \"rectilinear_dim0_size\"" );
546 if( attRectilinearDim1Size ==
nullptr )
548 _EXCEPTIONT(
"Missing attribute \"rectilinear_dim1_size\"" );
551 int nDim0Size = attRectilinearDim0Size->as_int( 0 );
552 int nDim1Size = attRectilinearDim1Size->as_int( 0 );
555 NcAtt* attRectilinearDim0Name = ncInput.get_att(
"rectilinear_dim0_name" );
556 NcAtt* attRectilinearDim1Name = ncInput.get_att(
"rectilinear_dim1_name" );
558 if( attRectilinearDim0Name ==
nullptr )
560 _EXCEPTIONT(
"Missing attribute \"rectilinear_dim0_name\"" );
562 if( attRectilinearDim1Name ==
nullptr )
564 _EXCEPTIONT(
"Missing attribute \"rectilinear_dim1_name\"" );
567 std::string strDim0Name = attRectilinearDim0Name->as_string( 0 );
568 std::string strDim1Name = attRectilinearDim1Name->as_string( 0 );
570 std::map< std::string, int > vecDimNameSizes;
572 vecDimNameSizes[strDim0Name] = nDim0Size;
573 vecDimNameSizes[strDim1Name] = nDim1Size;
575 vecDimSizes[1] = vecDimNameSizes[
"lat"];
576 vecDimSizes[2] = vecDimNameSizes[
"lon"];
578 else if( varGridDims !=
nullptr )
581 NcDim* dimGridRank = varGridDims->get_dim( 0 );
582 if( dimGridRank == NULL )
584 _EXCEPTIONT(
"Variable \"grid_dims\" has no dimensions" );
587 int gridrank = dimGridRank->size();
590 varGridDims->get( &( gridsizes[0] ), dimGridRank->size() );
601 vecDimSizes[1] = elems.
size();
607 vecDimSizes[1] = gridsizes[0];
608 vecDimSizes[2] = gridsizes[1];
620 vecDimSizes[1] = elems.
size();
625 result =
gMB->
tag_set_data( rectilinearTag, &mSet, 1, vecDimSizes.data() );
628 switch( vecDimSizes[0] )
631 printf(
"Cubed-Sphere mesh: %d (elements), %d (vertices)\n", vecDimSizes[1], vecDimSizes[2] );
634 printf(
"Rectilinear RLL mesh: (lon) %d X (lat) %d\n", vecDimSizes[2], vecDimSizes[1] );
637 printf(
"Icosahedral (triangular) mesh: %d (elements), %d (vertices)\n", vecDimSizes[1],
642 printf(
"Polygonal mesh: %d (elements), %d (vertices)\n", vecDimSizes[1], vecDimSizes[2] );
648 const size_t nOverlapFaces = tempestMesh->faces.size();
649 if( tempestMesh->vecSourceFaceIx.size() == nOverlapFaces &&
650 tempestMesh->vecSourceFaceIx.size() == nOverlapFaces )
653 use_overlap_context =
true;
666 std::vector< int > gids( faces.
size() ), srcpar( faces.
size() ), tgtpar( faces.
size() );
670 for(
unsigned ii = 0; ii < faces.
size(); ++ii )
672 srcpar[ii] = tempestMesh->vecSourceFaceIx[gids[ii] - 1];
673 tgtpar[ii] = tempestMesh->vecTargetFaceIx[gids[ii] - 1];
687 "Reading TempestRemap/SCRIP meshes (-T tempest input) requires NetCDF, which is "
688 "disabled in this MOAB build" );
691 else if( tempestout )
697 std::vector< int > metadata( 2 );
713 use_overlap_context =
true;
736 assert( metadata.size() );
737 std::cout <<
"Converting a RLL mesh with rectilinear dimension: " << metadata[0] <<
" X "
738 << metadata[1] << std::endl;
748 result =
gMB->
load_file( j->c_str(), &file_set, read_options.c_str() );
750 result =
gMB->
load_file( j->c_str(), &file_set, read_options.c_str() );
754 std::cerr <<
"Failed to load \"" << *j <<
"\"." << std::endl;
755 std::cerr <<
"Error code: " <<
gMB->
get_error_string( result ) <<
" (" << result <<
")" << std::endl;
758 std::cerr <<
"Error message: " << message << std::endl;
764 if( !proc_id ) std::cerr <<
"Read \"" << *j <<
"\"" << std::endl;
765 if( print_times && !proc_id )
write_times( std::cout );
769 bool have_geom =
false;
770 for( dim = 0; dim <= 3; ++dim )
772 if( !geom[dim].empty() ) have_geom =
true;
777 bool have_sets = have_geom;
785 std::cerr <<
"No ID tag defined." << std::endl;
791 std::cerr <<
"No geometry tag defined." << std::endl;
800 Tag tags[] = { id_tag, dim_tag };
801 const void* vals[] = { &id_val, &dim };
802 for( dim = 0; dim <= 3; ++dim )
804 int init_count = set_list.size();
805 for( std::set< int >::iterator iter = geom[dim].begin(); iter != geom[dim].end(); ++iter )
813 std::cerr << geom_names[dim] <<
" " << id_val <<
" not found.\n";
815 std::copy( range.
begin(), range.
end(), std::back_inserter( set_list ) );
819 std::cout <<
"Found " << ( set_list.size() - init_count ) <<
' ' << geom_names[dim] <<
" sets"
825 for( i = 0; i < 4; ++i )
829 if( mesh[i].empty() )
continue;
837 std::cerr <<
"Tag not found: " << mesh_tag_names[i] << std::endl;
842 int init_count = set_list.size();
843 for( std::set< int >::iterator iter = mesh[i].begin(); iter != mesh[i].end(); ++iter )
846 const void* vals[] = { &*iter };
851 std::cerr << mesh_tag_names[i] <<
" " << *iter <<
" not found.\n";
853 std::copy( range.
begin(), range.
end(), std::back_inserter( set_list ) );
857 std::cout <<
"Found " << ( set_list.size() - init_count ) <<
' ' << mesh_tag_names[i] <<
" sets"
863 for( dim = 1; dim < 4; ++dim )
864 if( dims[dim] ) bydim =
true;
869 if( generate[1] && !dims[1] )
871 std::cerr <<
"Warning: Request to generate 1D internal entities but not export them." << std::endl;
874 if( generate[2] && !dims[2] )
876 std::cerr <<
"Warning: Request to generate 2D internal entities but not export them." << std::endl;
882 if( generate[1] || generate[2] )
889 num_sets = set_list.size();
892 for( i = 0; i < num_sets; ++i )
894 Range dim3, dim2, adj;
914 Range dead_entities, tmp_range;
915 for( dim = 1; dim <= 3; ++dim )
917 if( dims[dim] )
continue;
919 dead_entities.
merge( tmp_range );
925 while( !empty_sets.
empty() )
928 dead_entities.
merge( empty_sets );
931 empty_sets =
subtract( tmp_range, dead_entities );
938 if( have_sets && set_list.empty() )
940 std::cerr <<
"Nothing to write." << std::endl;
948 if( !metis_partition_file.empty() )
953 std::cerr <<
"Failed to process partition file \"" << metis_partition_file <<
"\"." << std::endl;
963 std::cout <<
"Found " << set_list.size() <<
" specified sets to write (total)." << std::endl;
965 std::cout <<
"No sets specifed. Writing entire mesh." << std::endl;
970 #ifdef MOAB_HAVE_TEMPESTREMAP
980 int ntot_elements = 0, nelements = faces.
size();
982 int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, pcomm->
comm() );
983 if( ierr != 0 )
MB_CHK_SET_ERR( MB_FAILURE,
"MPI_Allreduce failed to get total source elements" );
985 ntot_elements = nelements;
989 std::vector< int > gids( faces.
size() );
993 if( faces.
size() > 1 && gids[0] == gids[1] && !use_overlap_context )
1009 if( spectral_order > 1 && globalid_tag_name.size() > 1 )
1013 result = remapper->
GenerateMeshMetadata( *tempestMesh, ntot_elements, faces, NULL, globalid_tag_name,
1021 const int nP = spectral_order;
1022 const int dofsPerElem = nP * nP;
1023 const std::string dTagName =
"D" + globalid_tag_name;
1033 std::vector< int > elemGids( faces.
size() );
1037 std::vector< int > dofIDs( dofsPerElem );
1038 for(
size_t ie = 0; ie < faces.
size(); ++ie )
1041 const int elemIdx = elemGids[ie] - 1;
1042 for(
int j = 0; j < nP; ++j )
1043 for(
int i = 0; i < nP; ++i )
1044 dofIDs[j * nP + i] = elemIdx * dofsPerElem + j * nP + i + 1;
1052 std::cout <<
"Generated discontinuous DoF tag \"" << dTagName <<
"\" with " << dofsPerElem
1053 <<
" DoFs/element (" << ntot_elements * dofsPerElem <<
" total)\n";
1061 if( use_overlap_context &&
false )
1063 const int nOverlapFaces = faces.
size();
1065 tempestMesh->vecSourceFaceIx.resize( nOverlapFaces );
1066 tempestMesh->vecTargetFaceIx.resize( nOverlapFaces );
1067 result =
gMB->
tag_get_data( srcParentTag, faces, &tempestMesh->vecSourceFaceIx[0] );
1069 result =
gMB->
tag_get_data( tgtParentTag, faces, &tempestMesh->vecTargetFaceIx[0] );
1073 #ifdef MOAB_HAVE_NETCDF
1074 tempestMesh->Write( out, NcFile::Netcdf4 );
1075 file_written =
true;
1078 "Writing TempestRemap meshes (-T tempest output) requires NetCDF, which is "
1079 "disabled in this MOAB build" );
1089 result =
gMB->
write_file( out.c_str(), format, write_options.c_str(), &set_list[0], set_list.size() );
1095 result =
gMB->
write_file( out.c_str(), format, write_options.c_str(), &file_set, 1 );
1098 std::cerr <<
"Failed to write \"" << out <<
"\"." << std::endl;
1099 std::cerr <<
"Error code: " <<
gMB->
get_error_string( result ) <<
" (" << result <<
")" << std::endl;
1100 std::string message;
1102 std::cerr <<
"Error message: " << message << std::endl;
1103 #ifdef MOAB_HAVE_MPI
1110 if( !proc_id ) std::cerr <<
"Wrote \"" << out <<
"\"" << std::endl;
1111 if( print_times && !proc_id )
write_times( std::cout );
1113 #ifdef MOAB_HAVE_MPI