20 #include <H5Tpublic.h>
29 #ifdef MOAB_HAVE_TEMPESTREMAP
30 #include "STLStringHelper.h"
53 #ifdef MOAB_HAVE_TEMPESTREMAP
54 struct TempestMapAppData
57 std::map< std::string, moab::TempestOnlineMap* > weightMaps;
93 std::map< int, ParCommGraph* > pgraph;
96 #ifdef MOAB_HAVE_TEMPESTREMAP
97 TempestMapAppData tempestData;
98 std::map< std::string, std::string > metadataMap;
118 std::vector< ParallelComm* > pcomms;
151 for(
int i = 0; i < 4; i++ )
167 MPI_Initialized( &flagInit );
215 std::string name( app_name );
219 std::cout <<
" application " << name <<
" already registered \n";
220 return moab::MB_FAILURE;
227 MPI_Comm_rank( *comm, &rankHere );
229 if( !rankHere ) std::cout <<
" application " << name <<
" with ID = " << *pid <<
" is registered now \n";
232 std::cout <<
" convention for external application is to have its id positive \n";
233 return moab::MB_FAILURE;
238 std::cout <<
" external application with comp id " << *compid <<
" is already registered\n";
239 return moab::MB_FAILURE;
251 int index = pco->
get_id();
252 assert( index == *pid );
256 context.pcomms.push_back( pco );
271 app_data.
name = name;
273 #ifdef MOAB_HAVE_TEMPESTREMAP
274 app_data.tempestData.remapper = NULL;
301 assert( app_name !=
nullptr );
302 std::string name( app_name );
311 ccomm = MPI_Comm_f2c( (MPI_Fint)*comm );
338 rankHere = pco->
rank();
341 std::cout <<
" application with ID: " << *pid <<
" global id: " << data.
global_id <<
" name: " << data.
name
342 <<
" is de-registered now \n";
349 fileents.
insert( fileSet );
353 #ifdef MOAB_HAVE_TEMPESTREMAP
354 if( data.tempestData.remapper )
delete data.tempestData.remapper;
355 if( data.tempestData.weightMaps.size() ) data.tempestData.weightMaps.clear();
363 std::map< int, ParCommGraph* >& pargs = data.pgraph;
366 for( std::map< int, ParCommGraph* >::iterator mt = pargs.begin(); mt != pargs.end(); mt++ )
372 if( pco )
delete pco;
386 if( !adj_ents_left.
empty() )
390 vertices =
subtract( vertices, conn_verts );
395 std::map< std::string, int >::iterator mit;
399 int pidx = mit->second;
408 std::map< int, int >::iterator mit1;
412 int pidx = mit1->second;
441 std::string& separator,
442 std::vector< std::string >& list_tag_names )
446 while( ( pos = input_names.find( separator ) ) != std::string::npos )
448 token = input_names.substr( 0, pos );
449 if( !token.empty() ) list_tag_names.push_back( token );
451 input_names.erase( 0, pos + separator.length() );
453 if( !input_names.empty() )
456 list_tag_names.push_back( input_names );
462 int* num_global_vertices,
463 int* num_global_elements,
468 IMOAB_ASSERT( strlen( filename ),
"Invalid filename length." );
470 #ifdef MOAB_HAVE_HDF5
471 std::string filen( filename );
476 if( num_global_vertices ) *num_global_vertices = 0;
477 if( num_global_elements ) *num_global_elements = 0;
478 if( num_dimension ) *num_dimension = 0;
479 if( num_parts ) *num_parts = 0;
483 unsigned long max_id;
486 file =
mhdf_openFile( filen.c_str(), 0, &max_id, -1, &status );
490 fprintf( stderr,
"%s: %s\n", filename,
mhdf_message( &status ) );
491 return moab::MB_FAILURE;
499 fprintf( stderr,
"%s: %s\n", filename,
mhdf_message( &status ) );
500 return moab::MB_FAILURE;
504 if( num_global_vertices ) *num_global_vertices = (int)data->
nodes.
count;
513 edges += ent_d->
count;
518 faces += ent_d->
count;
523 faces += ent_d->
count;
528 faces += ent_d->
count;
533 regions += ent_d->
count;
538 regions += ent_d->
count;
543 regions += ent_d->
count;
548 regions += ent_d->
count;
553 regions += ent_d->
count;
558 regions += ent_d->
count;
563 regions += ent_d->
count;
567 if( num_parts ) *num_parts = data->
numEntSets[0];
572 if( num_dimension ) *num_dimension = 1;
573 if( num_global_elements ) *num_global_elements = edges;
578 if( num_dimension ) *num_dimension = 2;
579 if( num_global_elements ) *num_global_elements = faces;
584 if( num_dimension ) *num_dimension = 3;
585 if( num_global_elements ) *num_global_elements = regions;
593 std::cout << filename
594 <<
": Please reconfigure with HDF5. Cannot retrieve header information for file "
595 "formats other than a h5m file.\n";
596 if( num_global_vertices ) *num_global_vertices = 0;
597 if( num_global_elements ) *num_global_elements = 0;
598 if( num_dimension ) *num_dimension = 0;
599 if( num_parts ) *num_parts = 0;
608 int* num_ghost_layers )
611 IMOAB_ASSERT( strlen( filename ),
"Invalid filename length." );
614 std::ostringstream newopts;
615 if( read_options ) newopts << read_options;
623 std::string opts( ( read_options ? read_options :
"" ) );
624 std::string pcid(
"PARALLEL_COMM=" );
625 std::size_t found = opts.find( pcid );
627 if( found != std::string::npos )
629 std::cerr <<
" cannot specify PARALLEL_COMM option, it is implicit \n";
630 return moab::MB_FAILURE;
635 std::string filen( filename );
636 std::string::size_type idx = filen.rfind(
'.' );
638 if( idx != std::string::npos )
640 std::string extension = filen.substr( idx + 1 );
641 if( ( extension == std::string(
"h5m" ) ) || ( extension == std::string(
"nc" ) ) )
642 newopts <<
";;PARALLEL_COMM=" << *pid;
645 if( *num_ghost_layers >= 1 )
649 std::string pcid2(
"PARALLEL_GHOSTS=" );
650 std::size_t found2 = opts.find( pcid2 );
652 if( found2 != std::string::npos )
654 std::cout <<
" PARALLEL_GHOSTS option is already specified, ignore passed "
655 "number of layers \n";
662 newopts <<
";PARALLEL_GHOSTS=3.0." << *num_ghost_layers <<
".3";
668 IMOAB_ASSERT( *num_ghost_layers == 0,
"Cannot provide ghost layers in serial." );
676 std::ostringstream outfile;
678 int rank =
context.pcomms[*pid]->rank();
679 int nprocs =
context.pcomms[*pid]->size();
680 outfile <<
"TaskMesh_n" << nprocs <<
"." << rank <<
".h5m";
682 outfile <<
"TaskMesh_n1.0.h5m";
696 IMOAB_ASSERT( strlen( filename ),
"Invalid filename length." );
701 std::ostringstream newopts;
703 std::string write_opts( ( write_options ? write_options :
"" ) );
704 std::string pcid(
"PARALLEL_COMM=" );
705 std::size_t found = write_opts.find( pcid );
707 if( found != std::string::npos )
709 std::cerr <<
" cannot specify PARALLEL_COMM option, it is implicit \n";
710 return moab::MB_FAILURE;
714 std::string pw(
"PARALLEL=WRITE_PART" );
715 found = write_opts.find( pw );
717 if( found != std::string::npos )
719 newopts <<
"PARALLEL_COMM=" << *pid <<
";";
724 if( write_options ) newopts << write_options;
726 std::vector< Tag > copyTagList = data.
tagList;
728 std::string gid_name_tag(
"GLOBAL_ID" );
731 if( data.
tagMap.find( gid_name_tag ) == data.
tagMap.end() )
734 copyTagList.push_back( gid );
737 std::string pp_name_tag(
"PARALLEL_PARTITION" );
740 if( data.
tagMap.find( pp_name_tag ) == data.
tagMap.end() )
744 if( ptag ) copyTagList.push_back( ptag );
757 IMOAB_ASSERT( strlen( prefix ),
"Invalid prefix string length." );
759 std::ostringstream file_name;
760 int rank = 0,
size = 1;
762 rank =
context.pcomms[*pid]->rank();
765 file_name << prefix <<
"_" <<
size <<
"_" << rank <<
".h5m";
840 int local[2], global[2];
843 MPI_Allreduce( local, global, 2, MPI_INT, MPI_SUM, pco->
comm() );
873 int* num_visible_vertices,
874 int* num_visible_elements,
875 int* num_visible_blocks,
876 int* num_visible_surfaceBC,
877 int* num_visible_vertexBC )
885 if( num_visible_elements )
889 num_visible_elements[0] =
static_cast< int >( data.
owned_elems.
size() );
890 num_visible_elements[1] =
static_cast< int >( data.
ghost_elems.
size() );
892 if( num_visible_vertices )
894 num_visible_vertices[2] =
static_cast< int >( data.
all_verts.
size() );
897 num_visible_vertices[0] = num_visible_vertices[2] - num_visible_vertices[1];
900 if( num_visible_blocks )
906 num_visible_blocks[0] = num_visible_blocks[2];
907 num_visible_blocks[1] = 0;
910 if( num_visible_surfaceBC )
915 num_visible_surfaceBC[2] = 0;
920 for(
int i = 0; i < numNeuSets; i++ )
929 Range adjPrimaryEnts;
932 num_visible_surfaceBC[2] += (int)adjPrimaryEnts.
size();
936 num_visible_surfaceBC[0] = num_visible_surfaceBC[2];
937 num_visible_surfaceBC[1] = 0;
940 if( num_visible_vertexBC )
945 num_visible_vertexBC[2] = 0;
948 for(
int i = 0; i < numDiriSets; i++ )
954 num_visible_vertexBC[2] += (int)verts.
size();
957 num_visible_vertexBC[0] = num_visible_vertexBC[2];
958 num_visible_vertexBC[1] = 0;
971 IMOAB_ASSERT( *vertices_length ==
static_cast< int >( verts.
size() ),
"Invalid vertices length provided" );
979 assert( vertices_length && *vertices_length );
980 assert( visible_global_rank_ID );
992 if( i != *vertices_length )
994 return moab::MB_FAILURE;
1000 if( (
int)verts.
size() != *vertices_length )
1002 return moab::MB_FAILURE;
1007 visible_global_rank_ID[i] = 0;
1020 if( *coords_length != 3 * (
int)verts.
size() )
1022 return moab::MB_FAILURE;
1036 if( *block_length != (
int)matSets.
size() )
1038 return moab::MB_FAILURE;
1046 for(
unsigned i = 0; i < matSets.
size(); i++ )
1048 matIdx[global_block_IDs[i]] = i;
1056 int* vertices_per_element,
1057 int* num_elements_in_block )
1059 assert( global_block_ID );
1062 std::map< int, int >::iterator it = matMap.find( *global_block_ID );
1064 if( it == matMap.end() )
1066 return moab::MB_FAILURE;
1069 int blockIndex = matMap[*global_block_ID];
1076 return moab::MB_FAILURE;
1083 return moab::MB_FAILURE;
1090 *vertices_per_element = num_verts;
1091 *num_elements_in_block = (int)blo_elems.
size();
1097 int* num_visible_elements,
1103 #ifdef MOAB_HAVE_MPI
1113 #ifdef MOAB_HAVE_MPI
1138 return moab::MB_FAILURE;
1141 if( -1 >= *num_visible_elements )
1143 return moab::MB_FAILURE;
1146 block_IDs[index] = valMatTag;
1155 int* connectivity_length,
1156 int* element_connectivity )
1158 assert( global_block_ID );
1159 assert( connectivity_length );
1162 std::map< int, int >& matMap = data.
matIndex;
1163 std::map< int, int >::iterator it = matMap.find( *global_block_ID );
1165 if( it == matMap.end() )
1167 return moab::MB_FAILURE;
1170 int blockIndex = matMap[*global_block_ID];
1172 std::vector< EntityHandle > elems;
1178 return moab::MB_FAILURE;
1181 std::vector< EntityHandle > vconnect;
1184 if( *connectivity_length != (
int)vconnect.size() )
1186 return moab::MB_FAILURE;
1189 for(
int i = 0; i < *connectivity_length; i++ )
1195 return moab::MB_FAILURE;
1198 element_connectivity[i] = inx;
1206 int* connectivity_length,
1207 int* element_connectivity )
1209 assert( elem_index );
1210 assert( connectivity_length );
1213 assert( ( *elem_index >= 0 ) && ( *elem_index < (
int)data.
primary_elems.
size() ) );
1222 if( *connectivity_length < num_nodes )
1224 return moab::MB_FAILURE;
1227 for(
int i = 0; i < num_nodes; i++ )
1233 return moab::MB_FAILURE;
1236 element_connectivity[i] = index;
1239 *connectivity_length = num_nodes;
1246 int* num_elements_in_block,
1247 int* element_ownership )
1249 assert( global_block_ID );
1250 assert( num_elements_in_block );
1254 std::map< int, int >::iterator it = matMap.find( *global_block_ID );
1256 if( it == matMap.end() )
1258 return moab::MB_FAILURE;
1261 int blockIndex = matMap[*global_block_ID];
1269 return moab::MB_FAILURE;
1272 if( *num_elements_in_block != (
int)elems.
size() )
1274 return moab::MB_FAILURE;
1278 #ifdef MOAB_HAVE_MPI
1284 #ifdef MOAB_HAVE_MPI
1287 element_ownership[i] = 0;
1296 int* num_elements_in_block,
1300 assert( global_block_ID );
1301 assert( num_elements_in_block );
1304 std::map< int, int >& matMap = data.
matIndex;
1306 std::map< int, int >::iterator it = matMap.find( *global_block_ID );
1308 if( it == matMap.end() )
1310 return moab::MB_FAILURE;
1313 int blockIndex = matMap[*global_block_ID];
1320 return moab::MB_FAILURE;
1323 if( *num_elements_in_block != (
int)elems.
size() )
1325 return moab::MB_FAILURE;
1331 for(
int i = 0; i < *num_elements_in_block; i++ )
1335 if( -1 == local_element_ID[i] )
1337 return moab::MB_FAILURE;
1345 int* surface_BC_length,
1347 int* reference_surface_ID,
1348 int* boundary_condition_value )
1359 for(
int i = 0; i < numNeuSets; i++ )
1371 Range adjPrimaryEnts;
1389 if( -1 == local_element_ID[index] )
1391 return moab::MB_FAILURE;
1398 boundary_condition_value[index] = neuVal;
1404 if( index != *surface_BC_length )
1406 return moab::MB_FAILURE;
1413 int* vertex_BC_length,
1415 int* boundary_condition_value )
1424 for(
int i = 0; i < numDiriSets; i++ )
1443 if( -1 == local_vertex_ID[index] )
1445 return moab::MB_FAILURE;
1448 boundary_condition_value[index] = diriVal;
1453 if( *vertex_BC_length != index )
1455 return moab::MB_FAILURE;
1464 int* components_per_entity,
1468 if( *tag_type < 0 || *tag_type > 5 )
1470 return moab::MB_FAILURE;
1475 void* defaultVal = NULL;
1476 int* defInt =
new int[*components_per_entity];
1477 double* defDouble =
new double[*components_per_entity];
1480 for(
int i = 0; i < *components_per_entity; i++ )
1483 defDouble[i] = -1e+10;
1492 defaultVal = defInt;
1498 defaultVal = defDouble;
1504 defaultVal = defHandle;
1510 defaultVal = defInt;
1516 defaultVal = defDouble;
1522 defaultVal = defHandle;
1529 return moab::MB_FAILURE;
1536 std::string tag_name( tag_storage_name );
1540 std::vector< std::string > tagNames;
1541 std::string separator(
":" );
1546 int already_defined_tags = 0;
1548 for(
size_t i = 0; i < tagNames.size(); i++ )
1555 std::map< std::string, Tag >& mTags = data.
tagMap;
1556 std::map< std::string, Tag >::iterator mit = mTags.find( tagNames[i].c_str() );
1558 if( mit == mTags.end() )
1561 mTags[tagNames[i]] = tagHandle;
1563 *tag_index = (int)data.
tagList.size();
1564 data.
tagList.push_back( tagHandle );
1567 already_defined_tags++;
1571 data.
tagMap[tagNames[i]] = tagHandle;
1572 *tag_index = (int)data.
tagList.size();
1573 data.
tagList.push_back( tagHandle );
1577 rval = moab::MB_FAILURE;
1582 #ifdef MOAB_HAVE_MPI
1584 rankHere = pco->
rank();
1586 if( !rankHere && already_defined_tags )
1587 std::cout <<
" application with ID: " << *pid <<
" global id: " << data.
global_id <<
" name: " << data.
name
1588 <<
" has " << already_defined_tags <<
" already defined tags out of " << tagNames.size()
1598 int* num_tag_storage_length,
1600 int* tag_storage_data )
1602 std::string tag_name( tag_storage_name );
1607 if( data.
tagMap.find( tag_name ) == data.
tagMap.end() )
1609 return moab::MB_FAILURE;
1622 return moab::MB_FAILURE;
1628 if( *ent_type == 0 )
1637 int nents_to_be_set = *num_tag_storage_length / tagLength;
1639 if( nents_to_be_set > (
int)ents_to_set->
size() || nents_to_be_set < 1 )
1641 return moab::MB_FAILURE;
1655 int* num_tag_storage_length,
1657 int* tag_storage_data )
1660 std::string tag_name( tag_storage_name );
1664 if( data.
tagMap.find( tag_name ) == data.
tagMap.end() )
1666 return moab::MB_FAILURE;
1679 return moab::MB_FAILURE;
1685 if( *ent_type == 0 )
1694 int nents_to_get = *num_tag_storage_length / tagLength;
1696 if( nents_to_get > (
int)ents_to_get->
size() || nents_to_get < 1 )
1698 return moab::MB_FAILURE;
1712 int* num_tag_storage_length,
1714 double* tag_storage_data )
1717 std::string tag_names( tag_storage_names );
1719 std::vector< std::string > tagNames;
1720 std::vector< Tag > tagHandles;
1721 std::string separator(
":" );
1725 Range* ents_to_set = NULL;
1727 if( *ent_type == 0 )
1731 else if( *ent_type == 1 )
1736 int nents_to_be_set = (int)( *ents_to_set ).
size();
1739 for(
size_t i = 0; i < tagNames.size(); i++ )
1741 if( data.
tagMap.find( tagNames[i] ) == data.
tagMap.end() )
1743 return moab::MB_FAILURE;
1756 return moab::MB_FAILURE;
1760 if( position + tagLength * nents_to_be_set > *num_tag_storage_length )
1761 return moab::MB_FAILURE;
1765 position = position + tagLength * nents_to_be_set;
1772 int* num_tag_storage_length,
1774 double* tag_storage_data,
1778 std::string tag_names( tag_storage_names );
1780 std::vector< std::string > tagNames;
1781 std::vector< Tag > tagHandles;
1782 std::string separator(
":" );
1786 Range* ents_to_set = NULL;
1788 if( *ent_type == 0 )
1792 else if( *ent_type == 1 )
1797 int nents_to_be_set = (int)( *ents_to_set ).
size();
1800 std::vector< int > gids;
1801 gids.resize( nents_to_be_set );
1807 std::map< int, EntityHandle > eh_by_gid;
1811 eh_by_gid[gids[i]] = *it;
1814 std::vector< int > tagLengths( tagNames.size() );
1815 std::vector< Tag > tagList;
1816 size_t total_tag_len = 0;
1817 for(
size_t i = 0; i < tagNames.size(); i++ )
1819 if( data.
tagMap.find( tagNames[i] ) == data.
tagMap.end() )
1821 MB_SET_ERR( moab::MB_FAILURE,
"tag missing" );
1825 tagList.push_back( tag );
1830 total_tag_len += tagLength;
1831 tagLengths[i] = tagLength;
1837 MB_SET_ERR( moab::MB_FAILURE,
"tag not double type" );
1841 #ifdef MOAB_HAVE_MPI
1843 unsigned num_procs = pco->
size();
1844 if( num_procs > 1 ) serial =
false;
1849 assert( total_tag_len * nents_to_be_set - *num_tag_storage_length == 0 );
1851 for(
int i = 0; i < nents_to_be_set; i++ )
1853 int gid = globalIds[i];
1856 int indexInTagValues = 0;
1857 for(
size_t j = 0; j < tagList.size(); j++ )
1859 indexInTagValues += i * tagLengths[j];
1862 indexInTagValues += ( nents_to_be_set - i ) * tagLengths[j];
1866 #ifdef MOAB_HAVE_MPI
1872 int nbLocalVals = *num_tag_storage_length / ( (int)tagNames.size() );
1873 assert( nbLocalVals * tagNames.size() - *num_tag_storage_length == 0 );
1875 TLsend.
initialize( 2, 0, 0, total_tag_len, nbLocalVals );
1879 int indexInRealLocal = 0;
1880 for(
int i = 0; i < nbLocalVals; i++ )
1883 int marker = globalIds[i];
1884 int to_proc = marker % num_procs;
1885 int n = TLsend.
get_n();
1886 TLsend.
vi_wr[2 * n] = to_proc;
1887 TLsend.
vi_wr[2 * n + 1] = marker;
1888 int indexInTagValues = 0;
1890 for(
size_t j = 0; j < tagList.size(); j++ )
1892 indexInTagValues += i * tagLengths[j];
1893 for(
int k = 0; k < tagLengths[j]; k++ )
1895 TLsend.
vr_wr[indexInRealLocal++] = tag_storage_data[indexInTagValues + k];
1897 indexInTagValues += ( nbLocalVals - i ) * tagLengths[j];
1902 assert( nbLocalVals * total_tag_len - indexInRealLocal == 0 );
1907 TLreq.
initialize( 2, 0, 0, 0, nents_to_be_set );
1909 for(
int i = 0; i < nents_to_be_set; i++ )
1912 int marker = gids[i];
1913 int to_proc = marker % num_procs;
1914 int n = TLreq.
get_n();
1915 TLreq.
vi_wr[2 * n] = to_proc;
1916 TLreq.
vi_wr[2 * n + 1] = marker;
1927 sort_buffer.buffer_init( TLreq.
get_n() );
1928 TLreq.
sort( 1, &sort_buffer );
1929 sort_buffer.
reset();
1930 sort_buffer.buffer_init( TLsend.
get_n() );
1931 TLsend.
sort( 1, &sort_buffer );
1932 sort_buffer.
reset();
1940 TLBack.
initialize( 3, 0, 0, total_tag_len, 0 );
1943 int n1 = TLreq.
get_n();
1944 int n2 = TLsend.
get_n();
1946 int indexInTLreq = 0;
1947 int indexInTLsend = 0;
1948 if( n1 > 0 && n2 > 0 )
1951 while( indexInTLreq < n1 && indexInTLsend < n2 )
1953 int currentValue1 = TLreq.
vi_rd[2 * indexInTLreq + 1];
1954 int currentValue2 = TLsend.
vi_rd[2 * indexInTLsend + 1];
1955 if( currentValue1 < currentValue2 )
1963 if( currentValue1 > currentValue2 )
1971 while( indexInTLreq + size1 < n1 && currentValue1 == TLreq.
vi_rd[2 * ( indexInTLreq + size1 ) + 1] )
1973 while( indexInTLsend + size2 < n2 && currentValue2 == TLsend.
vi_rd[2 * ( indexInTLsend + size2 ) + 1] )
1976 for(
int i1 = 0; i1 < size1; i1++ )
1978 for(
int i2 = 0; i2 < size2; i2++ )
1981 int n = TLBack.
get_n();
1983 TLBack.
vi_wr[3 * n] = TLreq.
vi_rd[2 * ( indexInTLreq + i1 )];
1985 TLBack.
vi_wr[3 * n + 1] = currentValue1;
1986 TLBack.
vi_wr[3 * n + 2] = TLsend.
vi_rd[2 * ( indexInTLsend + i2 )];
1988 for(
size_t k = 0; k < total_tag_len; k++ )
1990 TLBack.
vr_rd[total_tag_len * n + k] =
1991 TLsend.
vr_rd[total_tag_len * indexInTLsend + k];
1995 indexInTLreq += size1;
1996 indexInTLsend += size2;
2002 n1 = TLBack.
get_n();
2003 double* ptrVal = &TLBack.
vr_rd[0];
2004 for(
int i = 0; i < n1; i++ )
2006 int gid = TLBack.
vi_rd[3 * i + 1];
2010 for(
size_t j = 0; j < tagList.size(); j++ )
2014 ptrVal += tagLengths[j];
2023 int* num_tag_storage_length,
2025 double* tag_storage_data )
2029 std::string tag_names( tag_storage_names );
2031 std::vector< std::string > tagNames;
2032 std::vector< Tag > tagHandles;
2033 std::string separator(
":" );
2039 Range* ents_to_get = NULL;
2041 if( *ent_type == 0 )
2045 else if( *ent_type == 1 )
2049 int nents_to_get = (int)ents_to_get->
size();
2051 for(
size_t i = 0; i < tagNames.size(); i++ )
2053 if( data.
tagMap.find( tagNames[i] ) == data.
tagMap.end() )
2055 return moab::MB_FAILURE;
2068 return moab::MB_FAILURE;
2071 if( position + nents_to_get * tagLength > *num_tag_storage_length )
2072 return moab::MB_FAILURE;
2075 position = position + nents_to_get * tagLength;
2083 #ifdef MOAB_HAVE_MPI
2086 std::vector< Tag > tags;
2088 for(
int i = 0; i < *num_tag; i++ )
2090 if( tag_indices[i] < 0 || tag_indices[i] >= (
int)data.
tagList.size() )
2092 return moab::MB_FAILURE;
2095 tags.push_back( data.
tagList[tag_indices[i]] );
2098 if( *ent_type == 0 )
2102 else if( *ent_type == 1 )
2108 return moab::MB_FAILURE;
2119 int k = *pid + *num_tag + *tag_indices + *ent_type;
2128 #ifdef MOAB_HAVE_MPI
2132 if( *tag_index < 0 || *tag_index >= (
int)data.
tagList.size() )
2134 return moab::MB_FAILURE;
2139 if( *ent_type == 0 )
2143 else if( *ent_type == 1 )
2149 return moab::MB_FAILURE;
2160 int k = *pid + *tag_index + *ent_type;
2168 int* num_adjacent_elements,
2180 if( *num_adjacent_elements < (
int)adjs.
size() )
2182 return moab::MB_FAILURE;
2185 *num_adjacent_elements = (int)adjs.
size();
2187 for(
int i = 0; i < *num_adjacent_elements; i++ )
2211 return moab::MB_FAILURE;
2214 int nverts = *coords_len / *
dim;
2228 int* num_nodes_per_element,
2238 EntityType mbtype = (EntityType)( *type );
2247 for(
int j = 0; j < *num_elem * ( *num_nodes_per_element ); j++ )
2249 array[j] = connectivity[j] + firstVertex - 1;
2252 Range new_elems( actual_start_handle, actual_start_handle + *num_elem - 1 );
2263 int set_no = *block_ID;
2264 const void* setno_ptr = &set_no;
2269 if( MB_FAILURE == rval || sets.
empty() )
2281 block_set = sets[0];
2301 if( NULL != num_global_verts )
2305 if( NULL != num_global_elems )
2313 #ifdef MOAB_HAVE_MPI
2339 return moab::MB_FAILURE;
2356 std::cout <<
" can't get par part tag.\n";
2357 return moab::MB_FAILURE;
2360 int rank = pco->
rank();
2372 if( *num_ghost_layers <= 0 )
2389 ErrCode iMOAB_SendMesh(
iMOAB_AppID pid, MPI_Comm* join, MPI_Group* receivingGroup,
int* rcompid,
int* method )
2391 assert( join !=
nullptr );
2392 assert( receivingGroup !=
nullptr );
2393 assert( rcompid !=
nullptr );
2400 MPI_Comm global = ( data.
is_fortran ? MPI_Comm_f2c( *
reinterpret_cast< MPI_Fint*
>( join ) ) : *join );
2401 MPI_Group recvGroup =
2402 ( data.
is_fortran ? MPI_Group_f2c( *
reinterpret_cast< MPI_Fint*
>( receivingGroup ) ) : *receivingGroup );
2403 MPI_Comm sender = pco->
comm();
2405 MPI_Group senderGroup;
2406 ierr = MPI_Comm_group( sender, &senderGroup );
2407 if( ierr != 0 )
return moab::MB_FAILURE;
2418 int sender_rank = -1;
2419 MPI_Comm_rank( sender, &sender_rank );
2424 std::vector< int > number_elems_per_part;
2428 if( owned.
size() == 0 )
2437 int local_owned_elem = (int)owned.
size();
2439 int rank = pco->
rank();
2440 number_elems_per_part.resize(
size );
2441 number_elems_per_part[rank] = local_owned_elem;
2442 #if( MPI_VERSION >= 2 )
2444 ierr = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &number_elems_per_part[0], 1, MPI_INT, sender );
2447 std::vector< int > all_tmp(
size );
2448 ierr = MPI_Allgather( &number_elems_per_part[rank], 1, MPI_INT, &all_tmp[0], 1, MPI_INT, sender );
2449 number_elems_per_part = all_tmp;
2455 return moab::MB_FAILURE;
2476 MPI_Group_free( &senderGroup );
2480 ErrCode iMOAB_ReceiveMesh(
iMOAB_AppID pid, MPI_Comm* join, MPI_Group* sendingGroup,
int* scompid )
2482 assert( join !=
nullptr );
2483 assert( sendingGroup !=
nullptr );
2484 assert( scompid !=
nullptr );
2489 MPI_Comm receive = pco->
comm();
2492 MPI_Comm global = ( data.
is_fortran ? MPI_Comm_f2c( *
reinterpret_cast< MPI_Fint*
>( join ) ) : *join );
2493 MPI_Group sendGroup =
2494 ( data.
is_fortran ? MPI_Group_f2c( *
reinterpret_cast< MPI_Fint*
>( sendingGroup ) ) : *sendingGroup );
2498 MPI_Group receiverGroup;
2499 int ierr = MPI_Comm_group( receive, &receiverGroup );
CHK_MPI_ERR( ierr );
2508 int receiver_rank = -1;
2509 MPI_Comm_rank( receive, &receiver_rank );
2513 std::vector< int > pack_array;
2517 int current_receiver = cgraph->
receiver( receiver_rank );
2519 std::vector< int > senders_local;
2522 while( n < pack_array.size() )
2524 if( current_receiver == pack_array[n] )
2526 for(
int j = 0; j < pack_array[n + 1]; j++ )
2528 senders_local.push_back( pack_array[n + 2 + j] );
2534 n = n + 2 + pack_array[n + 1];
2538 std::cout <<
" receiver " << current_receiver <<
" at rank " << receiver_rank <<
" will receive from "
2539 << senders_local.size() <<
" tasks: ";
2541 for(
int k = 0; k < (int)senders_local.size(); k++ )
2543 std::cout <<
" " << senders_local[k];
2549 if( senders_local.empty() )
2551 std::cout <<
" we do not have any senders for receiver rank " << receiver_rank <<
"\n";
2567 if( (
int)senders_local.size() >= 2 )
2578 std::cout <<
"current_receiver " << current_receiver <<
" local verts: " << local_verts.
size() <<
"\n";
2582 rval = mm.merge_using_integer_tag( local_verts, idtag );
MB_CHK_ERR( rval );
2588 std::cout <<
"after merging: new verts: " << new_verts.
size() <<
"\n";
2606 if( NULL != densePartTag &&
MB_SUCCESS == rval )
2610 std::vector< int > vals;
2611 int rank = pco->
rank();
2612 vals.resize( local_verts.
size(), rank );
2624 std::cout <<
" can't get par part tag.\n";
2625 return moab::MB_FAILURE;
2628 int rank = pco->
rank();
2640 MPI_Group_free( &receiverGroup );
2648 std::map< int, ParCommGraph* >::iterator mt = data.pgraph.find( *context_id );
2649 if( mt == data.pgraph.end() )
2651 return moab::MB_FAILURE;
2659 MPI_Comm global = ( data.
is_fortran ? MPI_Comm_f2c( *
reinterpret_cast< MPI_Fint*
>( join ) ) : *join );
2668 #ifdef MOAB_HAVE_TEMPESTREMAP
2669 if( data.tempestData.remapper != NULL )
2683 if( 0 != cover_set )
2689 std::string tag_name( tag_storage_name );
2694 std::vector< std::string > tagNames;
2695 std::vector< Tag > tagHandles;
2696 std::string separator(
":" );
2698 for(
size_t i = 0; i < tagNames.size(); i++ )
2702 if(
MB_SUCCESS != rval || NULL == tagHandle )
2704 std::cout <<
" can't get tag handle for tag named:" << tagNames[i].c_str() <<
" at index " << i <<
"\n";
MB_CHK_SET_ERR( rval,
"can't get tag handle" );
2706 tagHandles.push_back( tagHandle );
2721 std::map< int, ParCommGraph* >::iterator mt = data.pgraph.find( *context_id );
2722 if( mt == data.pgraph.end() )
2724 return moab::MB_FAILURE;
2728 MPI_Comm global = ( data.
is_fortran ? MPI_Comm_f2c( *
reinterpret_cast< MPI_Fint*
>( join ) ) : *join );
2738 if( 0 != cover_set )
2749 #ifdef MOAB_HAVE_TEMPESTREMAP
2750 if( data.tempestData.remapper != NULL )
2763 std::string tag_name( tag_storage_name );
2766 std::vector< std::string > tagNames;
2767 std::vector< Tag > tagHandles;
2768 std::string separator(
":" );
2770 for(
size_t i = 0; i < tagNames.size(); i++ )
2774 if(
MB_SUCCESS != rval || NULL == tagHandle )
2776 std::cout <<
" can't get tag handle for tag named:" << tagNames[i].c_str() <<
" at index " << i <<
"\n";
MB_CHK_SET_ERR( rval,
"can't get tag handle" );
2778 tagHandles.push_back( tagHandle );
2782 std::cout << pco->
rank() <<
". Looking to receive data for tags: " << tag_name
2783 <<
" and file set = " << ( data.
file_set ) <<
"\n";
2790 std::cout << pco->
rank() <<
". Looking to receive data for tags: " << tag_name <<
"\n";
2801 std::map< int, ParCommGraph* >::iterator mt = data.pgraph.find( *context_id );
2802 if( mt == data.pgraph.end() )
return moab::MB_FAILURE;
2804 mt->second->release_send_buffers();
2827 int localRank = 0, numProcs = 1;
2829 bool isFortran =
false;
2830 if( *pid1 >= 0 ) isFortran = isFortran ||
context.
appDatas[*pid1].is_fortran;
2831 if( *pid2 >= 0 ) isFortran = isFortran ||
context.
appDatas[*pid2].is_fortran;
2833 MPI_Comm global = ( isFortran ? MPI_Comm_f2c( *
reinterpret_cast< MPI_Fint*
>( join ) ) : *join );
2834 MPI_Group srcGroup = ( isFortran ? MPI_Group_f2c( *
reinterpret_cast< MPI_Fint*
>( group1 ) ) : *group1 );
2835 MPI_Group tgtGroup = ( isFortran ? MPI_Group_f2c( *
reinterpret_cast< MPI_Fint*
>( group2 ) ) : *group2 );
2837 MPI_Comm_rank( global, &localRank );
2838 MPI_Comm_size( global, &numProcs );
2843 bool already_exists =
false;
2847 std::map< int, ParCommGraph* >::iterator mt = data.pgraph.find( *comp2 );
2848 if( mt != data.pgraph.end() ) already_exists =
true;
2853 std::map< int, ParCommGraph* >::iterator mt = data.pgraph.find( *comp1 );
2854 if( mt != data.pgraph.end() ) already_exists =
true;
2857 if( already_exists )
2861 std::cout <<
" parcomgraph already existing between components " << *comp1 <<
" and " << *comp2
2862 <<
". Do not compute again\n";
2870 if( *pid1 >= 0 ) cgraph =
new ParCommGraph( global, srcGroup, tgtGroup, *comp1, *comp2 );
2872 if( *pid2 >= 0 ) cgraph_rev =
new ParCommGraph( global, tgtGroup, srcGroup, *comp2, *comp1 );
2889 int lenTagType1 = 1;
2890 if( 1 == *type1 || 1 == *type2 )
2897 std::vector< int > valuesComp1;
2904 #ifdef MOAB_HAVE_TEMPESTREMAP
2905 if( data1.tempestData.remapper != NULL )
2911 Range ents_of_interest;
2916 valuesComp1.resize( ents_of_interest.
size() * lenTagType1 );
2919 else if( *type1 == 2 )
2922 valuesComp1.resize( ents_of_interest.
size() );
2925 else if( *type1 == 3 )
2928 valuesComp1.resize( ents_of_interest.
size() );
2937 std::set< int > uniq( valuesComp1.begin(), valuesComp1.end() );
2938 TLcomp1.
resize( uniq.size() );
2939 for( std::set< int >::iterator sit = uniq.begin(); sit != uniq.end(); sit++ )
2943 int to_proc = marker % numProcs;
2944 int n = TLcomp1.
get_n();
2945 TLcomp1.
vi_wr[2 * n] = to_proc;
2946 TLcomp1.
vi_wr[2 * n + 1] = marker;
2952 pc.crystal_router()->gs_transfer( 1, TLcomp1,
2956 std::stringstream ff1;
2957 ff1 <<
"TLcomp1_" << localRank <<
".txt";
2961 sort_buffer.buffer_init( TLcomp1.
get_n() );
2962 TLcomp1.
sort( 1, &sort_buffer );
2963 sort_buffer.
reset();
2972 std::vector< int > valuesComp2;
2978 #ifdef MOAB_HAVE_TEMPESTREMAP
2979 if( data2.tempestData.remapper != NULL )
2986 Range ents_of_interest;
2991 valuesComp2.resize( ents_of_interest.
size() * lenTagType1 );
2994 else if( *type2 == 2 )
2997 valuesComp2.resize( ents_of_interest.
size() );
3000 else if( *type2 == 3 )
3003 valuesComp2.resize( ents_of_interest.
size() );
3011 std::set< int > uniq( valuesComp2.begin(), valuesComp2.end() );
3012 TLcomp2.
resize( uniq.size() );
3013 for( std::set< int >::iterator sit = uniq.begin(); sit != uniq.end(); sit++ )
3017 int to_proc = marker % numProcs;
3018 int n = TLcomp2.
get_n();
3019 TLcomp2.
vi_wr[2 * n] = to_proc;
3020 TLcomp2.
vi_wr[2 * n + 1] = marker;
3024 pc.crystal_router()->gs_transfer( 1, TLcomp2,
3028 std::stringstream ff2;
3029 ff2 <<
"TLcomp2_" << localRank <<
".txt";
3032 sort_buffer.buffer_reserve( TLcomp2.
get_n() );
3033 TLcomp2.
sort( 1, &sort_buffer );
3034 sort_buffer.
reset();
3056 int n1 = TLcomp1.
get_n();
3057 int n2 = TLcomp2.
get_n();
3059 int indexInTLComp1 = 0;
3060 int indexInTLComp2 = 0;
3061 if( n1 > 0 && n2 > 0 )
3064 while( indexInTLComp1 < n1 && indexInTLComp2 < n2 )
3066 int currentValue1 = TLcomp1.
vi_rd[2 * indexInTLComp1 + 1];
3067 int currentValue2 = TLcomp2.
vi_rd[2 * indexInTLComp2 + 1];
3068 if( currentValue1 < currentValue2 )
3076 if( currentValue1 > currentValue2 )
3084 while( indexInTLComp1 + size1 < n1 && currentValue1 == TLcomp1.
vi_rd[2 * ( indexInTLComp1 + size1 ) + 1] )
3086 while( indexInTLComp2 + size2 < n2 && currentValue2 == TLcomp2.
vi_rd[2 * ( indexInTLComp2 + size2 ) + 1] )
3089 for(
int i1 = 0; i1 < size1; i1++ )
3091 for(
int i2 = 0; i2 < size2; i2++ )
3094 int n = TLBackToComp1.
get_n();
3096 TLBackToComp1.
vi_wr[3 * n] =
3097 TLcomp1.
vi_rd[2 * ( indexInTLComp1 + i1 )];
3099 TLBackToComp1.
vi_wr[3 * n + 1] = currentValue1;
3100 TLBackToComp1.
vi_wr[3 * n + 2] = TLcomp2.
vi_rd[2 * ( indexInTLComp2 + i2 )];
3101 n = TLBackToComp2.
get_n();
3103 TLBackToComp2.
vi_wr[3 * n] =
3104 TLcomp2.
vi_rd[2 * ( indexInTLComp2 + i2 )];
3105 TLBackToComp2.
vi_wr[3 * n + 1] = currentValue1;
3106 TLBackToComp2.
vi_wr[3 * n + 2] = TLcomp1.
vi_rd[2 * ( indexInTLComp1 + i1 )];
3110 indexInTLComp1 += size1;
3111 indexInTLComp2 += size2;
3114 pc.crystal_router()->gs_transfer( 1, TLBackToComp1, 0 );
3115 pc.crystal_router()->gs_transfer( 1, TLBackToComp2, 0 );
3124 std::stringstream f1;
3125 f1 <<
"TLBack1_" << localRank <<
".txt";
3128 sort_buffer.buffer_reserve( TLBackToComp1.
get_n() );
3129 TLBackToComp1.
sort( 1, &sort_buffer );
3130 sort_buffer.
reset();
3145 std::stringstream f2;
3146 f2 <<
"TLBack2_" << localRank <<
".txt";
3149 sort_buffer.buffer_reserve( TLBackToComp2.
get_n() );
3150 TLBackToComp2.
sort( 2, &sort_buffer );
3151 sort_buffer.
reset();
3169 std::vector< Tag > tagsList;
3172 if( !tag || rval !=
MB_SUCCESS )
return moab::MB_FAILURE;
3173 tagsList.push_back( tag );
3175 if( tag && rval ==
MB_SUCCESS ) tagsList.push_back( tag );
3177 if( tag && rval ==
MB_SUCCESS ) tagsList.push_back( tag );
3178 double tol = 1.0e-9;
3214 std::cout <<
" can't get par part tag.\n";
3215 return moab::MB_FAILURE;
3218 int rank = pco->
rank();
3224 #ifdef MOAB_HAVE_TEMPESTREMAP
3231 ErrCode iMOAB_CoverageGraph( MPI_Comm* join,
3242 std::vector< int > srcSenders;
3243 std::vector< int > receivers;
3246 int default_context_id = -1;
3247 bool is_fortran_context =
false;
3253 default_context_id = *migr_id;
3255 is_fortran_context =
context.
appDatas[*pid_src].is_fortran || is_fortran_context;
3259 srcSenders = sendGraph->
senders();
3262 std::cout <<
"senders: " << srcSenders.size() <<
" first sender: " << srcSenders[0] << std::endl;
3267 if( *pid_migr >= 0 )
3269 is_fortran_context =
context.
appDatas[*pid_migr].is_fortran || is_fortran_context;
3271 default_context_id = *src_id;
3275 srcSenders = recvGraph->
senders();
3278 std::cout <<
"receivers: " << receivers.size() <<
" first receiver: " << receivers[0] << std::endl;
3282 if( *pid_intx >= 0 ) is_fortran_context =
context.
appDatas[*pid_intx].is_fortran || is_fortran_context;
3295 MPI_Comm global = ( is_fortran_context ? MPI_Comm_f2c( *
reinterpret_cast< MPI_Fint*
>( join ) ) : *join );
3296 int currentRankInJointComm = -1;
3297 ierr = MPI_Comm_rank( global, ¤tRankInJointComm );
CHK_MPI_ERR( ierr );
3301 if( find( receivers.begin(), receivers.end(), currentRankInJointComm ) !=
3305 if( *pid_intx >= (
int)
context.
appDatas.size() )
return moab::MB_FAILURE;
3308 Tag parentTag, orgSendProcTag;
3311 if( !parentTag )
return moab::MB_FAILURE;
3314 if( !orgSendProcTag )
return moab::MB_FAILURE;
3322 std::map< int, std::set< int > > idsFromProcs;
3326 int gidCell, origProc;
3332 if( origProc < 0 )
continue;
3334 std::set< int >& setInts = idsFromProcs[origProc];
3335 setInts.insert( gidCell );
3344 assert( *pid_intx >= 0 );
3356 int gidCell, origProc;
3364 std::set< int >& setInts = idsFromProcs[origProc];
3365 setInts.insert( gidCell );
3370 std::ofstream dbfile;
3371 std::stringstream outf;
3372 outf <<
"idsFromProc_0" << currentRankInJointComm <<
".txt";
3373 dbfile.open( outf.str().c_str() );
3374 dbfile <<
"Writing this to a file.\n";
3376 dbfile <<
" map size:" << idsFromProcs.size()
3380 for( std::map<
int, std::set< int > >::iterator mt = idsFromProcs.begin(); mt != idsFromProcs.end(); mt++ )
3382 std::set< int >& setIds = mt->second;
3383 dbfile <<
"from id: " << mt->first <<
" receive " << setIds.size() <<
" cells \n";
3385 for( std::set< int >::iterator st = setIds.begin(); st != setIds.end(); st++ )
3388 dbfile <<
" " << valueID;
3390 if( counter % 10 == 0 ) dbfile <<
"\n";
3396 if( NULL != recvGraph )
3403 assert( *pid_intx >= 0 );
3409 for( std::map<
int, std::set< int > >::iterator mit = idsFromProcs.begin(); mit != idsFromProcs.end(); mit++ )
3411 int procToSendTo = mit->first;
3412 std::set< int >& idSet = mit->second;
3413 for( std::set< int >::iterator sit = idSet.begin(); sit != idSet.end(); sit++ )
3415 int n = TLcovIDs.
get_n();
3417 TLcovIDs.
vi_wr[2 * n] = procToSendTo;
3418 TLcovIDs.
vi_wr[2 * n + 1] = *sit;
3424 pc.crystal_router()->gs_transfer( 1, TLcovIDs,
3429 if( NULL != sendGraph )
3443 assert( prefix && strlen( prefix ) );
3446 std::string prefix_str( prefix );
3448 if( NULL != cgraph )
3452 std::cout <<
" cannot find ParCommGraph on app with pid " << *pid <<
" name: " <<
context.
appDatas[*pid].name
3453 <<
" context: " << *context_id <<
"\n";
3462 #ifdef MOAB_HAVE_TEMPESTREMAP
3464 #ifdef MOAB_HAVE_NETCDF
3465 ErrCode iMOAB_LoadMappingWeightsFromFile(
3474 bool row_based_partition =
true;
3475 if( *col_or_row == 1 ) row_based_partition =
false;
3481 TempestMapAppData& tdata = data_intx.tempestData;
3484 if( tdata.remapper == NULL )
3487 #ifdef MOAB_HAVE_MPI
3493 tdata.remapper->meshValidate =
true;
3494 tdata.remapper->constructEdgeMap =
true;
3496 tdata.remapper->initialize(
false );
3502 tdata.weightMaps[std::string( solution_weights_identifier )] =
new moab::TempestOnlineMap( tdata.remapper );
3506 assert( weightMap != NULL );
3517 int lenTagType1 = 1;
3525 std::vector< int > dofValues;
3536 dofValues.resize( ents_of_interest.
size() * lenTagType1 );
3538 ndofPerEl = lenTagType1;
3540 else if( *type == 2 )
3543 dofValues.resize( ents_of_interest.
size() );
3546 else if( *type == 3 )
3549 dofValues.resize( ents_of_interest.
size() );
3557 std::vector< int > orderDofs( dofValues.begin(), dofValues.end() );
3559 std::sort( orderDofs.begin(), orderDofs.end() );
3560 orderDofs.erase( std::unique( orderDofs.begin(), orderDofs.end() ), orderDofs.end() );
3565 if( row_based_partition )
3567 tdata.pid_dest = pid_cpl;
3574 tdata.pid_src = pid_cpl;
3582 std::vector< int > tmp_owned_ids;
3589 ErrCode iMOAB_WriteMappingWeightsToFile(
3594 assert( solution_weights_identifier && strlen( solution_weights_identifier ) );
3595 assert( remap_weights_filename && strlen( remap_weights_filename ) );
3601 TempestMapAppData& tdata = data_intx.tempestData;
3604 assert( tdata.remapper != NULL );
3608 assert( weightMap != NULL );
3610 std::string filename = std::string( remap_weights_filename );
3612 std::string metadataStr = data_intx.metadataMap[std::string( solution_weights_identifier )];
3613 std::map<std::string, std::string> attrMap;
3614 attrMap[
"title"] =
"MOAB-TempestRemap Online Regridding Weight Generator";
3615 attrMap[
"normalization"] =
"ovarea";
3616 attrMap[
"map_aPb"] = filename;
3618 const std::string delim =
";";
3619 size_t pos = 0, index = 0;
3620 std::vector<std::string> stringAttr(3);
3622 while (( pos = metadataStr.find (delim)) != std::string::npos)
3624 std::string token1 = metadataStr.substr(0, pos);
3625 if ( token1.size() > 0 || index == 0 )
3626 stringAttr[index++] = token1;
3627 metadataStr.erase(0, pos + delim.length());
3629 stringAttr[index] = metadataStr;
3631 attrMap[
"remap_options"] = stringAttr[0];
3632 attrMap[
"methodorder_b"] = stringAttr[1];
3633 attrMap[
"methodorder_a"] = stringAttr[2];
3634 attrMap[
"concave_a"] =
"false";
3635 attrMap[
"concave_b"] =
"false";
3636 attrMap[
"bubble"] =
"true";
3645 #ifdef MOAB_HAVE_MPI
3649 MPI_Comm* jointcomm,
3657 assert( jointcomm );
3660 bool is_fortran =
false;
3661 if( *pid1 >= 0 ) is_fortran =
context.
appDatas[*pid1].is_fortran || is_fortran;
3662 if( *pid2 >= 0 ) is_fortran =
context.
appDatas[*pid2].is_fortran || is_fortran;
3663 if( *pid3 >= 0 ) is_fortran =
context.
appDatas[*pid3].is_fortran || is_fortran;
3665 MPI_Comm joint_communicator =
3666 ( is_fortran ? MPI_Comm_f2c( *
reinterpret_cast< MPI_Fint*
>( jointcomm ) ) : *jointcomm );
3667 MPI_Group group_first = ( is_fortran ? MPI_Group_f2c( *
reinterpret_cast< MPI_Fint*
>( groupA ) ) : *groupA );
3668 MPI_Group group_second = ( is_fortran ? MPI_Group_f2c( *
reinterpret_cast< MPI_Fint*
>( groupB ) ) : *groupB );
3671 int localRank = 0, numProcs = 1;
3674 MPI_Comm_rank( joint_communicator, &localRank );
3675 MPI_Comm_size( joint_communicator, &numProcs );
3681 if( *pid1 >= 0 ) cgraph =
new ParCommGraph( joint_communicator, group_first, group_second, *comp1, *comp2 );
3684 if( *pid3 >= 0 ) cgraph_rev =
new ParCommGraph( joint_communicator, group_second, group_first, *comp2, *comp1 );
3705 int lenTagType1 = 1;
3713 std::vector< int > valuesComp1;
3716 Range ents_of_interest;
3727 valuesComp1.resize( ents_of_interest.
size() * lenTagType1 );
3730 else if( *type == 2 )
3733 valuesComp1.resize( ents_of_interest.
size() );
3736 else if( *type == 3 )
3739 valuesComp1.resize( ents_of_interest.
size() );
3748 std::set< int > uniq( valuesComp1.begin(), valuesComp1.end() );
3749 TLcomp1.
resize( uniq.size() );
3750 for( std::set< int >::iterator sit = uniq.begin(); sit != uniq.end(); sit++ )
3754 int to_proc = marker % numProcs;
3755 int n = TLcomp1.
get_n();
3756 TLcomp1.
vi_wr[2 * n] = to_proc;
3757 TLcomp1.
vi_wr[2 * n + 1] = marker;
3763 pc.crystal_router()->gs_transfer( 1, TLcomp1,
3767 std::stringstream ff1;
3768 ff1 <<
"TLcomp1_" << localRank <<
".txt";
3772 sort_buffer.buffer_init( TLcomp1.
get_n() );
3773 TLcomp1.
sort( 1, &sort_buffer );
3774 sort_buffer.
reset();
3786 std::vector< int > valuesComp2;
3790 TempestMapAppData& tdata = data2.tempestData;
3792 assert( tdata.weightMaps.size() == 1 );
3794 weightMap = tdata.weightMaps.begin()->second;
3797 if( *direction == 1 )
3803 else if( *direction == 2 )
3811 std::set< int > uniq( valuesComp2.begin(), valuesComp2.end() );
3812 TLcomp2.
resize( uniq.size() );
3813 for( std::set< int >::iterator sit = uniq.begin(); sit != uniq.end(); sit++ )
3817 int to_proc = marker % numProcs;
3818 int n = TLcomp2.
get_n();
3819 TLcomp2.
vi_wr[2 * n] = to_proc;
3820 TLcomp2.
vi_wr[2 * n + 1] = marker;
3824 pc.crystal_router()->gs_transfer( 1, TLcomp2,
3829 std::stringstream ff2;
3830 ff2 <<
"TLcomp2_" << localRank <<
".txt";
3833 sort_buffer.buffer_reserve( TLcomp2.
get_n() );
3834 TLcomp2.
sort( 1, &sort_buffer );
3835 sort_buffer.
reset();
3857 int n1 = TLcomp1.
get_n();
3858 int n2 = TLcomp2.
get_n();
3860 int indexInTLComp1 = 0;
3861 int indexInTLComp2 = 0;
3862 if( n1 > 0 && n2 > 0 )
3865 while( indexInTLComp1 < n1 && indexInTLComp2 < n2 )
3867 int currentValue1 = TLcomp1.
vi_rd[2 * indexInTLComp1 + 1];
3868 int currentValue2 = TLcomp2.
vi_rd[2 * indexInTLComp2 + 1];
3869 if( currentValue1 < currentValue2 )
3877 if( currentValue1 > currentValue2 )
3885 while( indexInTLComp1 + size1 < n1 && currentValue1 == TLcomp1.
vi_rd[2 * ( indexInTLComp1 + size1 ) + 1] )
3887 while( indexInTLComp2 + size2 < n2 && currentValue2 == TLcomp2.
vi_rd[2 * ( indexInTLComp2 + size2 ) + 1] )
3890 for(
int i1 = 0; i1 < size1; i1++ )
3892 for(
int i2 = 0; i2 < size2; i2++ )
3895 int n = TLBackToComp1.
get_n();
3897 TLBackToComp1.
vi_wr[3 * n] =
3898 TLcomp1.
vi_rd[2 * ( indexInTLComp1 + i1 )];
3900 TLBackToComp1.
vi_wr[3 * n + 1] = currentValue1;
3901 TLBackToComp1.
vi_wr[3 * n + 2] = TLcomp2.
vi_rd[2 * ( indexInTLComp2 + i2 )];
3902 n = TLBackToComp2.
get_n();
3904 TLBackToComp2.
vi_wr[3 * n] =
3905 TLcomp2.
vi_rd[2 * ( indexInTLComp2 + i2 )];
3906 TLBackToComp2.
vi_wr[3 * n + 1] = currentValue1;
3907 TLBackToComp2.
vi_wr[3 * n + 2] = TLcomp1.
vi_rd[2 * ( indexInTLComp1 + i1 )];
3911 indexInTLComp1 += size1;
3912 indexInTLComp2 += size2;
3915 pc.crystal_router()->gs_transfer( 1, TLBackToComp1, 0 );
3916 pc.crystal_router()->gs_transfer( 1, TLBackToComp2, 0 );
3928 std::stringstream f1;
3929 f1 <<
"TLBack1_" << localRank <<
".txt";
3934 rval = cgraph->
set_split_ranges( *comp1, TLBackToComp1, valuesComp1, lenTagType1, ents_of_interest, *type );
MB_CHK_ERR( rval );
3947 int size_tuple = 2 + ( ( *type != 1 ) ? 0 : lenTagType1 ) + 1 +
3953 pc.crystal_router()->gs_transfer( 1, TLv, 0 );
3954 if( *type != 2 ) pc.crystal_router()->gs_transfer( 1, TLc, 0 );
3960 Range primary_ents3;
3961 std::vector< int > values_entities;
3966 if( 1 == *type ) ndofPerEl = (int)( sqrt( lenTagType1 ) );
3968 assert( *pid2 >= 0 );
3970 TempestMapAppData& tdata = dataIntx.tempestData;
3973 if( 1 == *direction )
3975 tdata.pid_src = pid3;
3983 tdata.pid_dest = pid3;
4019 bool validate =
true;
4021 double radius_source = 1.0;
4022 double radius_target = 1.0;
4023 const double epsrel = ReferenceTolerance;
4024 constexpr
double boxeps = 1.e-1;
4025 constexpr
bool gnomonic =
false;
4031 #ifdef MOAB_HAVE_MPI
4036 TempestMapAppData& tdata = data_intx.tempestData;
4039 bool is_parallel =
false, is_root =
true;
4041 #ifdef MOAB_HAVE_MPI
4044 rank = pco_intx->
rank();
4045 is_parallel = ( pco_intx->
size() > 1 );
4046 is_root = ( rank == 0 );
4052 outputFormatter.set_prefix(
"[iMOAB_ComputeMeshIntersectionOnSphere]: " );
4058 ComputeSphereRadius( pid_src, &radius_source );
4059 ComputeSphereRadius( pid_tgt, &radius_target );
4062 outputFormatter.printf( 0,
"Radius of spheres: source = %12.14f, and target = %12.14f\n", radius_source,
4068 double defaultradius = 1.0;
4069 if( fabs( radius_source - radius_target ) > 1e-10 )
4076 #ifdef MOAB_HAVE_TEMPESTREMAP
4083 bool use_kdtree_search =
false;
4084 double srctgt_areas[2], srctgt_areas_glb[2];
4091 srctgt_areas[0] = areaAdaptor.area_on_sphere(
context.
MBI, data_src.
file_set, defaultradius );
4094 outputFormatter.printf( 0,
"The red set contains %d vertices and %d elements \n", rintxverts.
size(),
4095 rintxelems.
size() );
4103 srctgt_areas[1] = areaAdaptor.area_on_sphere(
context.
MBI, data_tgt.
file_set, defaultradius );
4106 outputFormatter.printf( 0,
"The blue set contains %d vertices and %d elements \n", bintxverts.
size(),
4107 bintxelems.
size() );
4109 #ifdef MOAB_HAVE_MPI
4110 MPI_Allreduce( &srctgt_areas[0], &srctgt_areas_glb[0], 2, MPI_DOUBLE, MPI_SUM, pco_intx->
comm() );
4112 srctgt_areas_glb[0] = srctgt_areas[0];
4113 srctgt_areas_glb[1] = srctgt_areas[1];
4115 use_kdtree_search = ( srctgt_areas_glb[0] < srctgt_areas_glb[1] );
4120 tdata.pid_src = pid_src;
4121 tdata.pid_dest = pid_tgt;
4124 #ifdef MOAB_HAVE_MPI
4129 tdata.remapper->meshValidate =
true;
4130 tdata.remapper->constructEdgeMap =
true;
4133 tdata.remapper->initialize(
false );
4142 rval = tdata.remapper->ConstructCoveringSet( epsrel, 1.0, 1.0, boxeps,
false, gnomonic );
MB_CHK_ERR( rval );
4146 rval = tdata.remapper->ComputeOverlapMesh( use_kdtree_search,
false );
MB_CHK_ERR( rval );
4153 local_area = areaAdaptor.area_on_sphere(
context.
MBI, data_intx.
file_set, radius_source );
4155 global_areas[0] = srctgt_areas_glb[0];
4156 global_areas[1] = srctgt_areas_glb[1];
4158 #ifdef MOAB_HAVE_MPI
4161 MPI_Reduce( &local_area, &global_areas[2], 1, MPI_DOUBLE, MPI_SUM, 0, pco_intx->
comm() );
4165 global_areas[2] = local_area;
4168 global_areas[2] = local_area;
4172 outputFormatter.printf( 0,
4173 "initial area: source mesh = %12.14f, target mesh = %12.14f, "
4174 "overlap mesh = %12.14f\n",
4175 global_areas[0], global_areas[1], global_areas[2] );
4176 outputFormatter.printf( 0,
" relative error w.r.t source = %12.14e, and target = %12.14e\n",
4177 fabs( global_areas[0] - global_areas[2] ) / global_areas[0],
4178 fabs( global_areas[1] - global_areas[2] ) / global_areas[1] );
4190 double radius_source = 1.0;
4191 double radius_target = 1.0;
4192 const double epsrel = ReferenceTolerance;
4193 const double boxeps = 1.e-8;
4199 #ifdef MOAB_HAVE_MPI
4204 TempestMapAppData& tdata = data_intx.tempestData;
4207 #ifdef MOAB_HAVE_MPI
4218 ComputeSphereRadius( pid_src, &radius_source );
4219 ComputeSphereRadius( pid_tgt, &radius_target );
4230 std::cout <<
"The red set contains " << rintxverts.
size() <<
" vertices and " << rintxelems.
size()
4240 std::cout <<
"The blue set contains " << bintxverts.
size() <<
" vertices and " << bintxelems.
size()
4248 tdata.pid_src = pid_src;
4249 tdata.pid_dest = pid_tgt;
4254 #ifdef MOAB_HAVE_MPI
4259 tdata.remapper->meshValidate =
true;
4260 tdata.remapper->constructEdgeMap =
true;
4263 tdata.remapper->initialize(
false );
4270 if( fabs( radius_source - radius_target ) > 1e-10 )
4280 rval = tdata.remapper->ConstructCoveringSet( epsrel, 1.0, 1.0, boxeps,
false );
MB_CHK_ERR( rval );
4282 #ifdef MOAB_HAVE_MPI
4291 rval = tdata.remapper->ComputeGlobalLocalMaps();
MB_CHK_ERR( rval );
4296 ErrCode iMOAB_ComputeScalarProjectionWeights(
4300 int* disc_order_source,
4302 int* disc_order_target,
4305 int* fMonotoneTypeID,
4307 int* fInverseDistanceMap,
4308 int* fNoConservation,
4315 assert( disc_order_source && disc_order_target && *disc_order_source > 0 && *disc_order_target > 0 );
4316 assert( solution_weights_identifier && strlen( solution_weights_identifier ) );
4317 assert( disc_method_source && strlen( disc_method_source ) );
4318 assert( disc_method_target && strlen( disc_method_target ) );
4319 assert( source_solution_tag_dof_name && strlen( source_solution_tag_dof_name ) );
4320 assert( target_solution_tag_dof_name && strlen( target_solution_tag_dof_name ) );
4324 TempestMapAppData& tdata = data_intx.tempestData;
4328 tdata.weightMaps[std::string( solution_weights_identifier )] =
new moab::TempestOnlineMap( tdata.remapper );
4332 assert( weightMap != NULL );
4334 GenerateOfflineMapAlgorithmOptions mapOptions;
4335 mapOptions.nPin = *disc_order_source;
4336 mapOptions.nPout = *disc_order_target;
4337 mapOptions.fSourceConcave =
false;
4338 mapOptions.fTargetConcave =
false;
4340 mapOptions.strMethod =
"";
4341 if( fv_method ) mapOptions.strMethod += std::string( fv_method ) +
";";
4342 if( fMonotoneTypeID )
4344 switch( *fMonotoneTypeID )
4347 mapOptions.fMonotone =
false;
4350 mapOptions.strMethod +=
"mono3;";
4353 mapOptions.strMethod +=
"mono2;";
4356 mapOptions.fMonotone =
true;
4360 mapOptions.fMonotone =
false;
4362 mapOptions.fNoBubble = ( fNoBubble ? *fNoBubble : false );
4363 mapOptions.fNoConservation = ( fNoConservation ? *fNoConservation > 0 : false );
4364 mapOptions.fNoCorrectAreas =
false;
4365 mapOptions.fNoCheck = !( fValidate ? *fValidate : true );
4366 if( fVolumetric && *fVolumetric ) mapOptions.strMethod +=
"volumetric;";
4367 if( fInverseDistanceMap && *fInverseDistanceMap ) mapOptions.strMethod +=
"invdist;";
4369 std::string metadataStr = mapOptions.strMethod +
";" +
4370 std::string( disc_method_source ) +
":" + std::to_string(*disc_order_source) +
":" + std::string( source_solution_tag_dof_name ) +
";" +
4371 std::string( disc_method_target ) +
":" + std::to_string(*disc_order_target) +
":" + std::string( target_solution_tag_dof_name );
4373 data_intx.metadataMap[std::string( solution_weights_identifier )] = metadataStr;
4379 std::string( disc_method_source ),
4380 std::string( disc_method_target ),
4382 std::string( source_solution_tag_dof_name ),
4383 std::string( target_solution_tag_dof_name )
4389 ErrCode iMOAB_ApplyScalarProjectionWeights(
4395 assert( solution_weights_identifier && strlen( solution_weights_identifier ) );
4396 assert( source_solution_tag_name && strlen( source_solution_tag_name ) );
4397 assert( target_solution_tag_name && strlen( target_solution_tag_name ) );
4403 TempestMapAppData& tdata = data_intx.tempestData;
4407 if( !tdata.weightMaps.count( std::string( solution_weights_identifier ) ) )
4412 std::vector< std::string > srcNames;
4413 std::vector< std::string > tgtNames;
4414 std::vector< Tag > srcTagHandles;
4415 std::vector< Tag > tgtTagHandles;
4416 std::string separator(
":" );
4417 std::string src_name( source_solution_tag_name );
4418 std::string tgt_name( target_solution_tag_name );
4421 if( srcNames.size() != tgtNames.size() )
4423 std::cout <<
" error in parsing source and target tag names. \n";
4424 return moab::MB_FAILURE;
4427 for(
size_t i = 0; i < srcNames.size(); i++ )
4431 if(
MB_SUCCESS != rval || NULL == tagHandle )
4433 return moab::MB_FAILURE;
4435 srcTagHandles.push_back( tagHandle );
4437 if(
MB_SUCCESS != rval || NULL == tagHandle )
4439 return moab::MB_FAILURE;
4441 tgtTagHandles.push_back( tagHandle );
4444 std::vector< double > solSTagVals;
4445 std::vector< double > solTTagVals;
4455 solSTagVals.resize( covSrcEnts.
size(), -1.0 );
4469 solTTagVals.resize( tgtEnts.
size(), -1.0 );
4495 for(
size_t i = 0; i < srcTagHandles.size(); i++ )
4498 Tag ssolnTag = srcTagHandles[i];
4499 Tag tsolnTag = tgtTagHandles[i];
4517 Tag ssolnTag = srcTagHandles[ivar];
4518 std::stringstream sstr;
4519 sstr <<
"covsrcTagData_" << *pid_intersection <<
"_" << ivar <<
"_" << pco_intx->
rank() <<
".txt";
4520 std::ofstream output_file( sstr.str().c_str() );
4521 for(
unsigned i = 0; i < sents.
size(); ++i )
4524 std::vector< double > locsolSTagVals( 16 );
4527 for(
unsigned j = 0; j < 16; ++j )
4528 output_file << locsolSTagVals[j] <<
" ";
4530 output_file.flush();
4531 output_file.close();
4534 std::stringstream sstr;
4535 sstr <<
"outputSrcDest_" << *pid_intersection <<
"_" << ivar <<
"_" << pco_intx->
rank() <<
".h5m";
4541 std::stringstream sstr;
4542 sstr <<
"outputCovSrcDest_" << *pid_intersection <<
"_" << ivar <<
"_" << pco_intx->
rank() <<
".h5m";
4549 std::stringstream sstr;
4550 sstr <<
"covMesh_" << *pid_intersection <<
"_" << pco_intx->
rank() <<
".vtk";
4556 std::stringstream sstr;
4557 sstr <<
"tgtMesh_" << *pid_intersection <<
"_" << pco_intx->
rank() <<
".vtk";
4563 std::stringstream sstr;
4564 sstr <<
"colvector_" << *pid_intersection <<
"_" << ivar <<
"_" << pco_intx->
rank() <<
".txt";
4565 std::ofstream output_file( sstr.str().c_str() );
4566 for(
unsigned i = 0; i < solSTagVals.size(); ++i )
4567 output_file << i <<
" " << weightMap->col_dofmap[i] <<
" " << weightMap->
col_gdofmap[i] <<
" "
4568 << solSTagVals[i] <<
"\n";
4569 output_file.flush();
4570 output_file.close();