ParCommGraph.cpp

Go to the documentation of this file.

/*
 * ParCommGraph.cpp
 *
 */
 
#include "moab/ParCommGraph.hpp"
// we need to recompute adjacencies for merging to work
#include "moab/Core.hpp"
#include "AEntityFactory.hpp"
 
#ifdef MOAB_HAVE_ZOLTAN
#include "moab/ZoltanPartitioner.hpp"
#endif
 
// #define VERBOSE
// #define GRAPH_INFO
 
namespace moab
{
ParCommGraph::ParCommGraph( MPI_Comm joincomm, MPI_Group group1, MPI_Group group2, int coid1, int coid2 )
 : comm( joincomm ), compid1( coid1 ), compid2( coid2 )
{
 // find out the tasks from each group, in the joint communicator
 find_group_ranks( group1, comm, senderTasks );
 find_group_ranks( group2, comm, receiverTasks );
 
 rootSender = rootReceiver = false;
 rankInGroup1 = rankInGroup2 = rankInJoin = -1; // not initialized, or not part of the group
 
 int mpierr = MPI_Group_rank( group1, &rankInGroup1 );
 if( MPI_SUCCESS != mpierr || rankInGroup1 == MPI_UNDEFINED ) rankInGroup1 = -1;
 
 mpierr = MPI_Group_rank( group2, &rankInGroup2 );
 if( MPI_SUCCESS != mpierr || rankInGroup2 == MPI_UNDEFINED ) rankInGroup2 = -1;
 
 mpierr = MPI_Comm_rank( comm, &rankInJoin );
 if( MPI_SUCCESS != mpierr ) // it should be a fatal error
 rankInJoin = -1;
 
 mpierr = MPI_Comm_size( comm, &joinSize );
 if( MPI_SUCCESS != mpierr ) // it should be a fatal error
 joinSize = -1;
 
 if( 0 == rankInGroup1 ) rootSender = true;
 if( 0 == rankInGroup2 ) rootReceiver = true;
 graph_type = INITIAL_MIGRATE; // 0
 comm_graph = NULL;
 context_id = -1;
 cover_set = 0; // refers to nothing yet
}
 
// copy constructor will copy only few basic things; split ranges will not be copied
ParCommGraph::ParCommGraph( const ParCommGraph& src )
{
 comm = src.comm;
 senderTasks = src.senderTasks; // these are the sender tasks in joint comm
 receiverTasks = src.receiverTasks; // these are all the receiver tasks in joint comm
 rootSender = src.rootSender;
 rootReceiver = src.rootReceiver;
 rankInGroup1 = src.rankInGroup1;
 rankInGroup2 = src.rankInGroup2; // group 1 is sender, 2 is receiver
 rankInJoin = src.rankInJoin;
 joinSize = src.joinSize;
 compid1 = src.compid1;
 compid2 = src.compid2;
 comm_graph = NULL;
 graph_type = src.graph_type;
 context_id = src.context_id;
 cover_set = src.cover_set;
 return;
}
 
ParCommGraph::~ParCommGraph()
{
 // TODO Auto-generated destructor stub
}
 
// utility to find out the ranks of the processes of a group, with respect to a joint comm,
// which spans for sure the group
// it is used locally (in the constructor), but it can be used as a utility
void ParCommGraph::find_group_ranks( MPI_Group group, MPI_Comm joincomm, std::vector< int >& ranks )
{
 MPI_Group global_grp;
 MPI_Comm_group( joincomm, &global_grp );
 
 int grp_size;
 
 MPI_Group_size( group, &grp_size );
 std::vector< int > rks( grp_size );
 ranks.resize( grp_size );
 
 for( int i = 0; i < grp_size; i++ )
 rks[i] = i;
 
 MPI_Group_translate_ranks( group, grp_size, &rks[0], global_grp, &ranks[0] );
 MPI_Group_free( &global_grp );
 return;
}
 
ErrorCode ParCommGraph::compute_trivial_partition( std::vector< int >& numElemsPerTaskInGroup1 )
{
 
 recv_graph.clear();
 recv_sizes.clear();
 sender_graph.clear();
 sender_sizes.clear();
 
 if( numElemsPerTaskInGroup1.size() != senderTasks.size() )
 return MB_FAILURE; // each sender has a number of elements that it owns
 
 // first find out total number of elements to be sent from all senders
 int total_elems = 0;
 std::vector< int > accum;
 accum.push_back( 0 );
 
 int num_senders = (int)senderTasks.size();
 
 for( size_t k = 0; k < numElemsPerTaskInGroup1.size(); k++ )
 {
 total_elems += numElemsPerTaskInGroup1[k];
 accum.push_back( total_elems );
 }
 
 int num_recv = ( (int)receiverTasks.size() );
 // in trivial partition, every receiver should get about total_elems/num_receivers elements
 int num_per_receiver = (int)( total_elems / num_recv );
 int leftover = total_elems - num_per_receiver * num_recv;
 
 // so receiver k will receive [starts[k], starts[k+1] ) interval
 std::vector< int > starts;
 starts.resize( num_recv + 1 );
 starts[0] = 0;
 for( int k = 0; k < num_recv; k++ )
 {
 starts[k + 1] = starts[k] + num_per_receiver;
 if( k < leftover ) starts[k + 1]++;
 }
 
 // each sender will send to a number of receivers, based on how the
 // arrays starts[0:num_recv] and accum[0:sendr] overlap
 int lastUsedReceiverRank = 0; // first receiver was not treated yet
 for( int j = 0; j < num_senders; j++ )
 {
 // we could start the receiver loop with the latest receiver that received from previous
 // sender
 for( int k = lastUsedReceiverRank; k < num_recv; k++ )
 {
 // if overlap:
 if( starts[k] < accum[j + 1] && starts[k + 1] > accum[j] )
 {
 recv_graph[receiverTasks[k]].push_back( senderTasks[j] );
 sender_graph[senderTasks[j]].push_back( receiverTasks[k] );
 
 // we still need to decide what is the overlap
 int sizeOverlap = 1; // at least 1, for sure
 // 1
 if( starts[k] >= accum[j] ) // one end is starts[k]
 {
 if( starts[k + 1] >= accum[j + 1] ) // the other end is accum[j+1]
 sizeOverlap = accum[j + 1] - starts[k];
 else //
 sizeOverlap = starts[k + 1] - starts[k];
 }
 else // one end is accum[j]
 {
 if( starts[k + 1] >= accum[j + 1] ) // the other end is accum[j+1]
 sizeOverlap = accum[j + 1] - accum[j];
 else
 sizeOverlap = starts[k + 1] - accum[j];
 }
 recv_sizes[receiverTasks[k]].push_back( sizeOverlap ); // basically, task k will receive from
 // sender j, sizeOverlap elems
 sender_sizes[senderTasks[j]].push_back( sizeOverlap );
 if( starts[k] > accum[j + 1] )
 {
 lastUsedReceiverRank = k - 1; // so next k loop will start a little higher, we
 // probably finished with first few receivers (up
 // to receiver lastUsedReceiverRank)
 break; // break the k loop, we distributed all elements from sender j to some
 // receivers
 }
 }
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode ParCommGraph::pack_receivers_graph( std::vector< int >& packed_recv_array )
{
 // it will basically look at local data, to pack communication graph, each receiver task will
 // have to post receives for each sender task that will send data to it; the array will be
 // communicated to root receiver, and eventually distributed to receiver tasks
 
 /*
 * packed_array will have receiver, number of senders, then senders, etc
 */
 if( recv_graph.size() < receiverTasks.size() )
 {
 // big problem, we have empty partitions in receive
 std::cout << " WARNING: empty partitions, some receiver tasks will receive nothing.\n";
 }
 for( std::map< int, std::vector< int > >::iterator it = recv_graph.begin(); it != recv_graph.end(); it++ )
 {
 int recv = it->first;
 std::vector< int >& senders = it->second;
 packed_recv_array.push_back( recv );
 packed_recv_array.push_back( (int)senders.size() );
 
 for( int k = 0; k < (int)senders.size(); k++ )
 packed_recv_array.push_back( senders[k] );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode ParCommGraph::split_owned_range( int sender_rank, Range& owned )
{
 int senderTask = senderTasks[sender_rank];
 std::vector< int >& distribution = sender_sizes[senderTask];
 std::vector< int >& receivers = sender_graph[senderTask];
 if( distribution.size() != receivers.size() ) //
 return MB_FAILURE;
 
 Range current = owned; // get the full range first, then we will subtract stuff, for
 // the following ranges
 
 Range rleftover = current;
 for( size_t k = 0; k < receivers.size(); k++ )
 {
 Range newr;
 newr.insert( current.begin(), current.begin() + distribution[k] );
 split_ranges[receivers[k]] = newr;
 
 rleftover = subtract( current, newr );
 current = rleftover;
 }
 
 return MB_SUCCESS;
}
 
// use for this the corresponding tasks and sizes
ErrorCode ParCommGraph::split_owned_range( Range& owned )
{
 if( corr_tasks.size() != corr_sizes.size() ) //
 return MB_FAILURE;
 
 Range current = owned; // get the full range first, then we will subtract stuff, for
 // the following ranges
 
 Range rleftover = current;
 for( size_t k = 0; k < corr_tasks.size(); k++ )
 {
 Range newr;
 newr.insert( current.begin(), current.begin() + corr_sizes[k] );
 split_ranges[corr_tasks[k]] = newr;
 
 rleftover = subtract( current, newr );
 current = rleftover;
 }
 
 return MB_SUCCESS;
}
 
ErrorCode ParCommGraph::send_graph( MPI_Comm jcomm )
{
 if( is_root_sender() )
 {
 int ierr;
 // will need to build a communication graph, because each sender knows now to which receiver
 // to send data the receivers need to post receives for each sender that will send data to
 // them will need to gather on rank 0 on the sender comm, global ranks of sender with
 // receivers to send build communication matrix, each receiver will receive from what sender
 
 std::vector< int > packed_recv_array;
 ErrorCode rval = pack_receivers_graph( packed_recv_array );
 if( MB_SUCCESS != rval ) return rval;
 
 int size_pack_array = (int)packed_recv_array.size();
 comm_graph = new int[size_pack_array + 1];
 comm_graph[0] = size_pack_array;
 for( int k = 0; k < size_pack_array; k++ )
 comm_graph[k + 1] = packed_recv_array[k];
 // will add 2 requests
 /// use tag 10 to send size and tag 20 to send the packed array
 sendReqs.resize( 1 );
 // do not send the size in advance, because we use probe now
 /*ierr = MPI_Isend(&comm_graph[0], 1, MPI_INT, receiver(0), 10, jcomm, &sendReqs[0]); // we
 have to use global communicator if (ierr!=0) return MB_FAILURE;*/
 ierr = MPI_Isend( &comm_graph[1], size_pack_array, MPI_INT, receiver( 0 ), 20, jcomm,
 &sendReqs[0] ); // we have to use global communicator
 if( ierr != 0 ) return MB_FAILURE;
 }
 return MB_SUCCESS;
}
 
// pco has MOAB too get_moab()
// do we need to store "method" as a member variable ?
ErrorCode ParCommGraph::send_mesh_parts( MPI_Comm jcomm, ParallelComm* pco, Range& owned )
{
 
 ErrorCode rval;
 if( split_ranges.empty() ) // in trivial partition
 {
 rval = split_owned_range( rankInGroup1, owned );
 if( rval != MB_SUCCESS ) return rval;
 // we know this on the sender side:
 corr_tasks = sender_graph[senderTasks[rankInGroup1]]; // copy
 corr_sizes = sender_sizes[senderTasks[rankInGroup1]]; // another copy
 }
 
 int indexReq = 0;
 int ierr; // MPI error
 if( is_root_sender() ) indexReq = 1; // for sendReqs
 sendReqs.resize( indexReq + split_ranges.size() );
 for( std::map< int, Range >::iterator it = split_ranges.begin(); it != split_ranges.end(); it++ )
 {
 int receiver_proc = it->first;
 Range ents = it->second;
 
 // add necessary vertices too
 Range verts;
 rval = pco->get_moab()->get_adjacencies( ents, 0, false, verts, Interface::UNION );
 if( rval != MB_SUCCESS )
 {
 std::cout << " can't get adjacencies. for entities to send\n";
 return rval;
 }
 ents.merge( verts );
 ParallelComm::Buffer* buffer = new ParallelComm::Buffer( ParallelComm::INITIAL_BUFF_SIZE );
 buffer->reset_ptr( sizeof( int ) );
 rval = pco->pack_buffer( ents, false, true, false, -1, buffer );
 if( rval != MB_SUCCESS )
 {
 std::cout << " can't pack buffer for entities to send\n";
 return rval;
 }
 int size_pack = buffer->get_current_size();
 
 // TODO there could be an issue with endian things; check !!!!!
 // we are sending the size of the buffer first as an int!!!
 /// not anymore !
 /* ierr = MPI_Isend(buffer->mem_ptr, 1, MPI_INT, receiver_proc, 1, jcomm,
 &sendReqs[indexReq]); // we have to use global communicator if (ierr!=0) return MB_FAILURE;
 indexReq++;*/
 
 ierr = MPI_Isend( buffer->mem_ptr, size_pack, MPI_CHAR, receiver_proc, 2, jcomm,
 &sendReqs[indexReq] ); // we have to use global communicator
 if( ierr != 0 ) return MB_FAILURE;
 indexReq++;
 localSendBuffs.push_back( buffer );
 }
 return MB_SUCCESS;
}
 
// this is called on receiver side
ErrorCode ParCommGraph::receive_comm_graph( MPI_Comm jcomm, ParallelComm* pco, std::vector< int >& pack_array )
{
 // first, receive from sender_rank 0, the communication graph (matrix), so each receiver
 // knows what data to expect
 MPI_Comm receive = pco->comm();
 int size_pack_array, ierr;
 MPI_Status status;
 if( rootReceiver )
 {
 /*
 * MPI_Probe(
 int source,
 int tag,
 MPI_Comm comm,
 MPI_Status* status)
 *
 */
 ierr = MPI_Probe( sender( 0 ), 20, jcomm, &status );
 if( 0 != ierr )
 {
 std::cout << " MPI_Probe failure: " << ierr << "\n";
 return MB_FAILURE;
 }
 // get the count of data received from the MPI_Status structure
 ierr = MPI_Get_count( &status, MPI_INT, &size_pack_array );
 if( 0 != ierr )
 {
 std::cout << " MPI_Get_count failure: " << ierr << "\n";
 return MB_FAILURE;
 }
#ifdef VERBOSE
 std::cout << " receive comm graph size: " << size_pack_array << "\n";
#endif
 pack_array.resize( size_pack_array );
 ierr = MPI_Recv( &pack_array[0], size_pack_array, MPI_INT, sender( 0 ), 20, jcomm, &status );
 if( 0 != ierr ) return MB_FAILURE;
#ifdef VERBOSE
 std::cout << " receive comm graph ";
 for( int k = 0; k < (int)pack_array.size(); k++ )
 std::cout << " " << pack_array[k];
 std::cout << "\n";
#endif
 }
 
 // now broadcast this whole array to all receivers, so they know what to expect
 ierr = MPI_Bcast( &size_pack_array, 1, MPI_INT, 0, receive );
 if( 0 != ierr ) return MB_FAILURE;
 pack_array.resize( size_pack_array );
 ierr = MPI_Bcast( &pack_array[0], size_pack_array, MPI_INT, 0, receive );
 if( 0 != ierr ) return MB_FAILURE;
 return MB_SUCCESS;
}
 
ErrorCode ParCommGraph::receive_mesh( MPI_Comm jcomm,
 ParallelComm* pco,
 EntityHandle local_set,
 std::vector< int >& senders_local )
{
 ErrorCode rval;
 int ierr;
 MPI_Status status;
 // we also need to fill corresponding mesh info on the other side
 corr_tasks = senders_local;
 Range newEnts;
 
 Tag orgSendProcTag; // this will be a tag set on the received mesh, with info about from what
 // task / PE the
 // primary element came from, in the joint communicator ; this will be forwarded by coverage
 // mesh
 int defaultInt = -1; // no processor, so it was not migrated from somewhere else
 rval = pco->get_moab()->tag_get_handle( "orig_sending_processor", 1, MB_TYPE_INTEGER, orgSendProcTag,
 MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create original sending processor tag" );
 if( !senders_local.empty() )
 {
 for( size_t k = 0; k < senders_local.size(); k++ )
 {
 int sender1 = senders_local[k];
 // first receive the size of the buffer using probe
 /*
 * MPI_Probe(
 int source,
 int tag,
 MPI_Comm comm,
 MPI_Status* status)
 *
 */
 ierr = MPI_Probe( sender1, 2, jcomm, &status );
 if( 0 != ierr )
 {
 std::cout << " MPI_Probe failure in ParCommGraph::receive_mesh " << ierr << "\n";
 return MB_FAILURE;
 }
 // get the count of data received from the MPI_Status structure
 int size_pack;
 ierr = MPI_Get_count( &status, MPI_CHAR, &size_pack );
 if( 0 != ierr )
 {
 std::cout << " MPI_Get_count failure in ParCommGraph::receive_mesh " << ierr << "\n";
 return MB_FAILURE;
 }
 
 /* ierr = MPI_Recv (&size_pack, 1, MPI_INT, sender1, 1, jcomm, &status);
 if (0!=ierr) return MB_FAILURE;*/
 // now resize the buffer, then receive it
 ParallelComm::Buffer* buffer = new ParallelComm::Buffer( size_pack );
 // buffer->reserve(size_pack);
 
 ierr = MPI_Recv( buffer->mem_ptr, size_pack, MPI_CHAR, sender1, 2, jcomm, &status );
 if( 0 != ierr )
 {
 std::cout << " MPI_Recv failure in ParCommGraph::receive_mesh " << ierr << "\n";
 return MB_FAILURE;
 }
 // now unpack the buffer we just received
 Range entities;
 std::vector< std::vector< EntityHandle > > L1hloc, L1hrem;
 std::vector< std::vector< int > > L1p;
 std::vector< EntityHandle > L2hloc, L2hrem;
 std::vector< unsigned int > L2p;
 
 buffer->reset_ptr( sizeof( int ) );
 std::vector< EntityHandle > entities_vec( entities.size() );
 std::copy( entities.begin(), entities.end(), entities_vec.begin() );
 rval = pco->unpack_buffer( buffer->buff_ptr, false, -1, -1, L1hloc, L1hrem, L1p, L2hloc, L2hrem, L2p,
 entities_vec );
 delete buffer;
 if( MB_SUCCESS != rval ) return rval;
 
 std::copy( entities_vec.begin(), entities_vec.end(), range_inserter( entities ) );
 // we have to add them to the local set
 rval = pco->get_moab()->add_entities( local_set, entities );
 if( MB_SUCCESS != rval ) return rval;
 // corr_sizes is the size of primary entities received
 Range verts = entities.subset_by_dimension( 0 );
 Range local_primary_ents = subtract( entities, verts );
 if( local_primary_ents.empty() )
 {
 // it is possible that all ents sent were vertices (point cloud)
 // then consider primary entities the vertices
 local_primary_ents = verts;
 }
 else
 {
 // set a tag with the original sender for the primary entity
 // will be used later for coverage mesh
 std::vector< int > orig_senders( local_primary_ents.size(), sender1 );
 rval = pco->get_moab()->tag_set_data( orgSendProcTag, local_primary_ents, &orig_senders[0] );
 }
 corr_sizes.push_back( (int)local_primary_ents.size() );
 
 newEnts.merge( entities );
 // make these in split ranges
 split_ranges[sender1] = local_primary_ents;
 
#ifdef VERBOSE
 std::ostringstream partial_outFile;
 
 partial_outFile << "part_send_" << sender1 << "."
 << "recv" << rankInJoin << ".vtk";
 
 // the mesh contains ghosts too, but they are not part of mat/neumann set
 // write in serial the file, to see what tags are missing
 std::cout << " writing from receiver " << rankInJoin << " from sender " << sender1
 << " entities: " << entities.size() << std::endl;
 rval = pco->get_moab()->write_file( partial_outFile.str().c_str(), 0, 0, &local_set,
 1 ); // everything on local set received
 if( MB_SUCCESS != rval ) return rval;
#endif
 }
 }
 // make sure adjacencies are updated on the new elements
 
 if( newEnts.empty() )
 {
 std::cout << " WARNING: this task did not receive any entities \n";
 }
 // in order for the merging to work, we need to be sure that the adjacencies are updated
 // (created)
 Range local_verts = newEnts.subset_by_type( MBVERTEX );
 newEnts = subtract( newEnts, local_verts );
 Core* mb = (Core*)pco->get_moab();
 AEntityFactory* adj_fact = mb->a_entity_factory();
 if( !adj_fact->vert_elem_adjacencies() )
 adj_fact->create_vert_elem_adjacencies();
 else
 {
 for( Range::iterator it = newEnts.begin(); it != newEnts.end(); ++it )
 {
 EntityHandle eh = *it;
 const EntityHandle* conn = NULL;
 int num_nodes = 0;
 rval = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
 adj_fact->notify_create_entity( eh, conn, num_nodes );
 }
 }
 
 return MB_SUCCESS;
}
 
// VSM: Why is the communicator never used. Remove the argument ?
ErrorCode ParCommGraph::release_send_buffers()
{
 int ierr, nsize = (int)sendReqs.size();
 std::vector< MPI_Status > mult_status;
 mult_status.resize( sendReqs.size() );
 ierr = MPI_Waitall( nsize, &sendReqs[0], &mult_status[0] );
 
 if( ierr != 0 ) return MB_FAILURE;
 // now we can free all buffers
 delete[] comm_graph;
 comm_graph = NULL;
 std::vector< ParallelComm::Buffer* >::iterator vit;
 for( vit = localSendBuffs.begin(); vit != localSendBuffs.end(); ++vit )
 delete( *vit );
 localSendBuffs.clear();
 return MB_SUCCESS;
}
 
// again, will use the send buffers, for nonblocking sends;
// should be the receives non-blocking too?
ErrorCode ParCommGraph::send_tag_values( MPI_Comm jcomm,
 ParallelComm* pco,
 Range& owned,
 std::vector< Tag >& tag_handles )
{
 // basically, owned.size() needs to be equal to sum(corr_sizes)
 // get info about the tag size, type, etc
 int ierr;
 Core* mb = (Core*)pco->get_moab();
 // get info about the tag
 //! Get the size of the specified tag in bytes
 int total_bytes_per_entity = 0; // we need to know, to allocate buffers
 ErrorCode rval;
 std::vector< int > vect_bytes_per_tag;
#ifdef VERBOSE
 std::vector< int > tag_sizes;
#endif
 for( size_t i = 0; i < tag_handles.size(); i++ )
 {
 int bytes_per_tag;
 rval = mb->tag_get_bytes( tag_handles[i], bytes_per_tag );MB_CHK_ERR( rval );
 int tag_size1; // length
 rval = mb->tag_get_length( tag_handles[i], tag_size1 );MB_CHK_ERR( rval );
 if( graph_type == DOF_BASED )
 bytes_per_tag = bytes_per_tag / tag_size1; // we know we have one double per tag , per ID sent;
 // could be 8 for double, 4 for int, etc
 total_bytes_per_entity += bytes_per_tag;
 vect_bytes_per_tag.push_back( bytes_per_tag );
#ifdef VERBOSE
 int tag_size;
 rval = mb->tag_get_length( tag_handles[i], tag_size );MB_CHK_ERR( rval );
 tag_sizes.push_back( tag_size );
#endif
 }
 
 int indexReq = 0;
 if( graph_type == INITIAL_MIGRATE ) // original send
 {
 // use the buffers data structure to allocate memory for sending the tags
 sendReqs.resize( split_ranges.size() );
 
 for( std::map< int, Range >::iterator it = split_ranges.begin(); it != split_ranges.end(); it++ )
 {
 int receiver_proc = it->first;
 Range ents = it->second; // primary entities, with the tag data
 int size_buffer = 4 + total_bytes_per_entity *
 (int)ents.size(); // hopefully, below 2B; if more, we have a big problem ...
 ParallelComm::Buffer* buffer = new ParallelComm::Buffer( size_buffer );
 
 buffer->reset_ptr( sizeof( int ) );
 for( size_t i = 0; i < tag_handles.size(); i++ )
 {
 // copy tag data to buffer->buff_ptr, and send the buffer (we could have used
 // regular char arrays)
 rval = mb->tag_get_data( tag_handles[i], ents, (void*)( buffer->buff_ptr ) );MB_CHK_ERR( rval );
 // advance the butter
 buffer->buff_ptr += vect_bytes_per_tag[i] * ents.size();
 }
 *( (int*)buffer->mem_ptr ) = size_buffer;
 // int size_pack = buffer->get_current_size(); // debug check
 ierr = MPI_Isend( buffer->mem_ptr, size_buffer, MPI_CHAR, receiver_proc, 222, jcomm,
 &sendReqs[indexReq] ); // we have to use global communicator
 if( ierr != 0 ) return MB_FAILURE;
 indexReq++;
 localSendBuffs.push_back( buffer ); // we will release them after nonblocking sends are completed
 }
 }
 else if( graph_type == COVERAGE )
 {
 // we know that we will need to send some tag data in a specific order
 // first, get the ids of the local elements, from owned Range; arrange the buffer in order
 // of increasing global id
 Tag gidTag = mb->globalId_tag();
 std::vector< int > gids;
 gids.resize( owned.size() );
 rval = mb->tag_get_data( gidTag, owned, &gids[0] );MB_CHK_ERR( rval );
 std::map< int, EntityHandle > gidToHandle;
 size_t i = 0;
 for( Range::iterator it = owned.begin(); it != owned.end(); it++ )
 {
 EntityHandle eh = *it;
 gidToHandle[gids[i++]] = eh;
 }
 // now, pack the data and send it
 sendReqs.resize( involved_IDs_map.size() );
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin();
 mit != involved_IDs_map.end(); mit++ )
 {
 int receiver_proc = mit->first;
 std::vector< int >& eids = mit->second;
 int size_buffer = 4 + total_bytes_per_entity *
 (int)eids.size(); // hopefully, below 2B; if more, we have a big problem ...
 ParallelComm::Buffer* buffer = new ParallelComm::Buffer( size_buffer );
 buffer->reset_ptr( sizeof( int ) );
#ifdef VERBOSE
 std::ofstream dbfile;
 std::stringstream outf;
 outf << "from_" << rankInJoin << "_send_to_" << receiver_proc << ".txt";
 dbfile.open( outf.str().c_str() );
 dbfile << "from " << rankInJoin << " send to " << receiver_proc << "\n";
#endif
 // copy tag data to buffer->buff_ptr, and send the buffer (we could have used regular
 // char arrays) pack data by tag, to be consistent with above, even though we loop
 // through the entities for each tag
 
 for( std::vector< int >::iterator it = eids.begin(); it != eids.end(); it++ )
 {
 int eID = *it;
 EntityHandle eh = gidToHandle[eID];
 for( i = 0; i < tag_handles.size(); i++ )
 {
 rval = mb->tag_get_data( tag_handles[i], &eh, 1, (void*)( buffer->buff_ptr ) );
 if( rval != MB_SUCCESS )
 {
 delete buffer; // free parallel comm buffer first
 
 MB_SET_ERR( rval, "Tag get data failed" );
 }
#ifdef VERBOSE
 dbfile << "global ID " << eID << " local handle " << mb->id_from_handle( eh ) << " vals: ";
 double* vals = (double*)( buffer->buff_ptr );
 for( int kk = 0; kk < tag_sizes[i]; kk++ )
 {
 dbfile << " " << *vals;
 vals++;
 }
 dbfile << "\n";
#endif
 buffer->buff_ptr += vect_bytes_per_tag[i];
 }
 }
 
#ifdef VERBOSE
 dbfile.close();
#endif
 *( (int*)buffer->mem_ptr ) = size_buffer;
 // int size_pack = buffer->get_current_size(); // debug check
 ierr = MPI_Isend( buffer->mem_ptr, size_buffer, MPI_CHAR, receiver_proc, 222, jcomm,
 &sendReqs[indexReq] ); // we have to use global communicator
 if( ierr != 0 ) return MB_FAILURE;
 indexReq++;
 localSendBuffs.push_back( buffer ); // we will release them after nonblocking sends are completed
 }
 }
 else if( graph_type == DOF_BASED )
 {
 // need to fill up the buffer, in the order desired, send it
 // get all the tags, for all owned entities, and pack the buffers accordingly
 // we do not want to get the tags by entity, it may be too expensive
 std::vector< std::vector< double > > valuesTags;
 valuesTags.resize( tag_handles.size() );
 for( size_t i = 0; i < tag_handles.size(); i++ )
 {
 
 int bytes_per_tag;
 rval = mb->tag_get_bytes( tag_handles[i], bytes_per_tag );MB_CHK_ERR( rval );
 valuesTags[i].resize( owned.size() * bytes_per_tag / sizeof( double ) );
 // fill the whole array, we will pick up from here
 rval = mb->tag_get_data( tag_handles[i], owned, (void*)( &( valuesTags[i][0] ) ) );MB_CHK_ERR( rval );
 }
 // now, pack the data and send it
 sendReqs.resize( involved_IDs_map.size() );
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin();
 mit != involved_IDs_map.end(); mit++ )
 {
 int receiver_proc = mit->first;
 std::vector< int >& eids = mit->second;
 std::vector< int >& index_in_values = map_index[receiver_proc];
 std::vector< int >& index_ptr = map_ptr[receiver_proc]; // this is eids.size()+1;
 int size_buffer = 4 + total_bytes_per_entity *
 (int)eids.size(); // hopefully, below 2B; if more, we have a big problem ...
 ParallelComm::Buffer* buffer = new ParallelComm::Buffer( size_buffer );
 buffer->reset_ptr( sizeof( int ) );
#ifdef VERBOSE
 std::ofstream dbfile;
 std::stringstream outf;
 outf << "from_" << rankInJoin << "_send_to_" << receiver_proc << ".txt";
 dbfile.open( outf.str().c_str() );
 dbfile << "from " << rankInJoin << " send to " << receiver_proc << "\n";
#endif
 // copy tag data to buffer->buff_ptr, and send the buffer
 // pack data by tag, to be consistent with above
 int j = 0;
 for( std::vector< int >::iterator it = eids.begin(); it != eids.end(); it++, j++ )
 {
 int index_in_v = index_in_values[index_ptr[j]];
 for( size_t i = 0; i < tag_handles.size(); i++ )
 {
 // right now, move just doubles; but it could be any type of tag
 *( (double*)( buffer->buff_ptr ) ) = valuesTags[i][index_in_v];
 buffer->buff_ptr += 8; // we know we are working with doubles only !!!
 }
 };
 *( (int*)buffer->mem_ptr ) = size_buffer;
 // int size_pack = buffer->get_current_size(); // debug check
 ierr = MPI_Isend( buffer->mem_ptr, size_buffer, MPI_CHAR, receiver_proc, 222, jcomm,
 &sendReqs[indexReq] ); // we have to use global communicator
 if( ierr != 0 ) return MB_FAILURE;
 indexReq++;
 localSendBuffs.push_back( buffer ); // we will release them after nonblocking sends are completed
 }
 }
 return MB_SUCCESS;
}
 
ErrorCode ParCommGraph::receive_tag_values( MPI_Comm jcomm,
 ParallelComm* pco,
 Range& owned,
 std::vector< Tag >& tag_handles )
{
 // opposite to sending, we will use blocking receives
 int ierr;
 MPI_Status status;
 // basically, owned.size() needs to be equal to sum(corr_sizes)
 // get info about the tag size, type, etc
 Core* mb = (Core*)pco->get_moab();
 // get info about the tag
 //! Get the size of the specified tag in bytes
 ErrorCode rval;
 int total_bytes_per_entity = 0;
 std::vector< int > vect_bytes_per_tag;
#ifdef VERBOSE
 std::vector< int > tag_sizes;
#endif
 for( size_t i = 0; i < tag_handles.size(); i++ )
 {
 int bytes_per_tag;
 rval = mb->tag_get_bytes( tag_handles[i], bytes_per_tag );MB_CHK_ERR( rval );
 total_bytes_per_entity += bytes_per_tag;
 vect_bytes_per_tag.push_back( bytes_per_tag );
#ifdef VERBOSE
 int tag_size;
 rval = mb->tag_get_length( tag_handles[i], tag_size );MB_CHK_ERR( rval );
 tag_sizes.push_back( tag_size );
#endif
 }
 
 if( graph_type == INITIAL_MIGRATE )
 {
 // std::map<int, Range> split_ranges;
 // rval = split_owned_range ( owned);MB_CHK_ERR ( rval );
 
 // use the buffers data structure to allocate memory for receiving the tags
 for( std::map< int, Range >::iterator it = split_ranges.begin(); it != split_ranges.end(); it++ )
 {
 int sender_proc = it->first;
 Range ents = it->second; // primary entities, with the tag data, we will receive
 int size_buffer = 4 + total_bytes_per_entity *
 (int)ents.size(); // hopefully, below 2B; if more, we have a big problem ...
 ParallelComm::Buffer* buffer = new ParallelComm::Buffer( size_buffer );
 
 buffer->reset_ptr( sizeof( int ) );
 
 *( (int*)buffer->mem_ptr ) = size_buffer;
 // int size_pack = buffer->get_current_size(); // debug check
 
 ierr = MPI_Recv( buffer->mem_ptr, size_buffer, MPI_CHAR, sender_proc, 222, jcomm, &status );
 if( ierr != 0 ) return MB_FAILURE;
 // now set the tag
 // copy to tag
 
 for( size_t i = 0; i < tag_handles.size(); i++ )
 {
 rval = mb->tag_set_data( tag_handles[i], ents, (void*)( buffer->buff_ptr ) );
 buffer->buff_ptr += vect_bytes_per_tag[i] * ents.size();
 }
 delete buffer; // no need for it afterwards
 MB_CHK_ERR( rval );
 }
 }
 else if( graph_type == COVERAGE ) // receive buffer, then extract tag data, in a loop
 {
 // we know that we will need to receive some tag data in a specific order (by ids stored)
 // first, get the ids of the local elements, from owned Range; unpack the buffer in order
 Tag gidTag = mb->globalId_tag();
 std::vector< int > gids;
 gids.resize( owned.size() );
 rval = mb->tag_get_data( gidTag, owned, &gids[0] );MB_CHK_ERR( rval );
 std::map< int, EntityHandle > gidToHandle;
 size_t i = 0;
 for( Range::iterator it = owned.begin(); it != owned.end(); it++ )
 {
 EntityHandle eh = *it;
 gidToHandle[gids[i++]] = eh;
 }
 //
 // now, unpack the data and set it to the tag
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin();
 mit != involved_IDs_map.end(); mit++ )
 {
 int sender_proc = mit->first;
 std::vector< int >& eids = mit->second;
 int size_buffer = 4 + total_bytes_per_entity *
 (int)eids.size(); // hopefully, below 2B; if more, we have a big problem ...
 ParallelComm::Buffer* buffer = new ParallelComm::Buffer( size_buffer );
 buffer->reset_ptr( sizeof( int ) );
 *( (int*)buffer->mem_ptr ) = size_buffer; // this is really not necessary, it should receive this too
 
 // receive the buffer
 ierr = MPI_Recv( buffer->mem_ptr, size_buffer, MPI_CHAR, sender_proc, 222, jcomm, &status );
 if( ierr != 0 ) return MB_FAILURE;
// start copy
#ifdef VERBOSE
 std::ofstream dbfile;
 std::stringstream outf;
 outf << "recvFrom_" << sender_proc << "_on_proc_" << rankInJoin << ".txt";
 dbfile.open( outf.str().c_str() );
 dbfile << "recvFrom_" << sender_proc << " on proc " << rankInJoin << "\n";
#endif
 
 // copy tag data from buffer->buff_ptr
 // data is arranged by tag , and repeat the loop for each entity ()
 // maybe it should be arranged by entity now, not by tag (so one loop for entities,
 // outside)
 
 for( std::vector< int >::iterator it = eids.begin(); it != eids.end(); it++ )
 {
 int eID = *it;
 std::map< int, EntityHandle >::iterator mit2 = gidToHandle.find( eID );
 if( mit2 == gidToHandle.end() )
 {
 std::cout << " on rank: " << rankInJoin << " cannot find entity handle with global ID " << eID
 << "\n";
 return MB_FAILURE;
 }
 EntityHandle eh = mit2->second;
 for( i = 0; i < tag_handles.size(); i++ )
 {
 rval = mb->tag_set_data( tag_handles[i], &eh, 1, (void*)( buffer->buff_ptr ) );MB_CHK_ERR( rval );
#ifdef VERBOSE
 dbfile << "global ID " << eID << " local handle " << mb->id_from_handle( eh ) << " vals: ";
 double* vals = (double*)( buffer->buff_ptr );
 for( int kk = 0; kk < tag_sizes[i]; kk++ )
 {
 dbfile << " " << *vals;
 vals++;
 }
 dbfile << "\n";
#endif
 buffer->buff_ptr += vect_bytes_per_tag[i];
 }
 }
 
 // delete receive buffer
 delete buffer;
#ifdef VERBOSE
 dbfile.close();
#endif
 }
 }
 else if( graph_type == DOF_BASED )
 {
 // need to fill up the values for each tag, in the order desired, from the buffer received
 //
 // get all the tags, for all owned entities, and pack the buffers accordingly
 // we do not want to get the tags by entity, it may be too expensive
 std::vector< std::vector< double > > valuesTags;
 valuesTags.resize( tag_handles.size() );
 for( size_t i = 0; i < tag_handles.size(); i++ )
 {
 int bytes_per_tag;
 rval = mb->tag_get_bytes( tag_handles[i], bytes_per_tag );MB_CHK_ERR( rval );
 valuesTags[i].resize( owned.size() * bytes_per_tag / sizeof( double ) );
 // fill the whole array, we will pick up from here
 // we will fill this array, using data from received buffer
 // rval = mb->tag_get_data(owned, (void*)( &(valuesTags[i][0])) );MB_CHK_ERR ( rval );
 }
 // now, unpack the data and set the tags
 sendReqs.resize( involved_IDs_map.size() );
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin();
 mit != involved_IDs_map.end(); mit++ )
 {
 int sender_proc = mit->first;
 std::vector< int >& eids = mit->second;
 std::vector< int >& index_in_values = map_index[sender_proc];
 std::vector< int >& index_ptr = map_ptr[sender_proc]; // this is eids.size()+1;
 int size_buffer = 4 + total_bytes_per_entity *
 (int)eids.size(); // hopefully, below 2B; if more, we have a big problem ...
 ParallelComm::Buffer* buffer = new ParallelComm::Buffer( size_buffer );
 buffer->reset_ptr( sizeof( int ) );
 
 // receive the buffer
 ierr = MPI_Recv( buffer->mem_ptr, size_buffer, MPI_CHAR, sender_proc, 222, jcomm, &status );
 if( ierr != 0 ) return MB_FAILURE;
 // use the values in buffer to populate valuesTag arrays, fill it up!
 int j = 0;
 for( std::vector< int >::iterator it = eids.begin(); it != eids.end(); it++, j++ )
 {
 for( size_t i = 0; i < tag_handles.size(); i++ )
 {
 // right now, move just doubles; but it could be any type of tag
 double val = *( (double*)( buffer->buff_ptr ) );
 buffer->buff_ptr += 8; // we know we are working with doubles only !!!
 for( int k = index_ptr[j]; k < index_ptr[j + 1]; k++ )
 valuesTags[i][index_in_values[k]] = val;
 }
 }
 // we are done with the buffer in which we received tags, release / delete it
 delete buffer;
 }
 // now we populated the values for all tags; set now the tags!
 for( size_t i = 0; i < tag_handles.size(); i++ )
 {
 // we will fill this array, using data from received buffer
 rval = mb->tag_set_data( tag_handles[i], owned, (void*)( &( valuesTags[i][0] ) ) );MB_CHK_ERR( rval );
 }
 }
 return MB_SUCCESS;
}
/*
 * for example
 */
ErrorCode ParCommGraph::settle_send_graph( TupleList& TLcovIDs )
{
 // fill involved_IDs_map with data
 // will have "receiving proc" and global id of element
 int n = TLcovIDs.get_n();
 graph_type = COVERAGE; // do not rely only on involved_IDs_map.size(); this can be 0 in some cases
 for( int i = 0; i < n; i++ )
 {
 int to_proc = TLcovIDs.vi_wr[2 * i];
 int globalIdElem = TLcovIDs.vi_wr[2 * i + 1];
 involved_IDs_map[to_proc].push_back( globalIdElem );
 }
#ifdef VERBOSE
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin(); mit != involved_IDs_map.end();
 mit++ )
 {
 std::cout << " towards task " << mit->first << " send: " << mit->second.size() << " cells " << std::endl;
 for( size_t i = 0; i < mit->second.size(); i++ )
 {
 std::cout << " " << mit->second[i];
 }
 std::cout << std::endl;
 }
#endif
 return MB_SUCCESS;
}
 
// this will set involved_IDs_map will store all ids to be received from one sender task
void ParCommGraph::SetReceivingAfterCoverage(
 std::map< int, std::set< int > >& idsFromProcs ) // will make sense only on receivers, right now after cov
{
 for( std::map< int, std::set< int > >::iterator mt = idsFromProcs.begin(); mt != idsFromProcs.end(); mt++ )
 {
 int fromProc = mt->first;
 std::set< int >& setIds = mt->second;
 involved_IDs_map[fromProc].resize( setIds.size() );
 std::vector< int >& listIDs = involved_IDs_map[fromProc];
 size_t indx = 0;
 for( std::set< int >::iterator st = setIds.begin(); st != setIds.end(); st++ )
 {
 int valueID = *st;
 listIDs[indx++] = valueID;
 }
 }
 graph_type = COVERAGE;
 return;
}
//#define VERBOSE
void ParCommGraph::settle_comm_by_ids( int comp, TupleList& TLBackToComp, std::vector< int >& valuesComp )
{
 // settle comm graph on comp
 if( rootSender || rootReceiver ) std::cout << " settle comm graph by id on component " << comp << "\n";
 int n = TLBackToComp.get_n();
 // third_method = true; // do not rely only on involved_IDs_map.size(); this can be 0 in some
 // cases
 std::map< int, std::set< int > > uniqueIDs;
 
 for( int i = 0; i < n; i++ )
 {
 int to_proc = TLBackToComp.vi_wr[3 * i + 2];
 int globalId = TLBackToComp.vi_wr[3 * i + 1];
 uniqueIDs[to_proc].insert( globalId );
 }
 
 // Vector to store element
 // with respective present index
 std::vector< std::pair< int, int > > vp;
 vp.reserve( valuesComp.size() );
 
 // Inserting element in pair vector
 // to keep track of previous indexes in valuesComp
 for( int i = 0; i < (int)valuesComp.size(); ++i )
 {
 vp.push_back( std::make_pair( valuesComp[i], i ) );
 }
 // Sorting pair vector
 sort( vp.begin(), vp.end() );
 
 // vp[i].first, second
 
 // count now how many times some value appears in ordered (so in valuesComp)
 for( std::map< int, std::set< int > >::iterator it = uniqueIDs.begin(); it != uniqueIDs.end(); it++ )
 {
 int procId = it->first;
 std::set< int >& nums = it->second;
 std::vector< int >& indx = map_ptr[procId];
 std::vector< int >& indices = map_index[procId];
 indx.resize( nums.size() + 1 );
 int indexInVp = 0;
 int indexVal = 0;
 indx[0] = 0; // start from 0
 for( std::set< int >::iterator sst = nums.begin(); sst != nums.end(); sst++, indexVal++ )
 {
 int val = *sst;
 involved_IDs_map[procId].push_back( val );
 indx[indexVal + 1] = indx[indexVal];
 while( ( indexInVp < (int)valuesComp.size() ) && ( vp[indexInVp].first <= val ) ) // should be equal !
 {
 if( vp[indexInVp].first == val )
 {
 indx[indexVal + 1]++;
 indices.push_back( vp[indexInVp].second );
 }
 indexInVp++;
 }
 }
 }
#ifdef VERBOSE
 std::stringstream f1;
 std::ofstream dbfile;
 f1 << "Involve_" << comp << "_" << rankInJoin << ".txt";
 dbfile.open( f1.str().c_str() );
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin(); mit != involved_IDs_map.end();
 mit++ )
 {
 int corrTask = mit->first;
 std::vector< int >& corrIds = mit->second;
 std::vector< int >& indx = map_ptr[corrTask];
 std::vector< int >& indices = map_index[corrTask];
 
 dbfile << " towards proc " << corrTask << " \n";
 for( int i = 0; i < (int)corrIds.size(); i++ )
 {
 dbfile << corrIds[i] << " [" << indx[i] << "," << indx[i + 1] << ") : ";
 for( int j = indx[i]; j < indx[i + 1]; j++ )
 dbfile << indices[j] << " ";
 dbfile << "\n";
 }
 dbfile << " \n";
 }
 dbfile.close();
#endif
 
 graph_type = DOF_BASED;
 // now we need to fill back and forth information, needed to fill the arrays
 // for example, from spectral to involved_IDs_map, in case we want to send data from
 // spectral to phys
}
//#undef VERBOSE
// new partition calculation
ErrorCode ParCommGraph::compute_partition( ParallelComm* pco, Range& owned, int met )
{
 // we are on a task on sender, and need to compute a new partition;
 // primary cells need to be distributed to nb receivers tasks
 // first, we will use graph partitioner, with zoltan;
 // in the graph that we need to build, the first layer of ghosts is needed;
 // can we avoid that ? For example, we can find out from each boundary edge/face what is the
 // other cell (on the other side), then form the global graph, and call zoltan in parallel met 1
 // would be a geometric partitioner, and met 2 would be a graph partitioner for method 1 we do
 // not need any ghost exchange
 
 // find first edges that are shared
 if( owned.empty() )
 return MB_SUCCESS; // nothing to do? empty partition is not allowed, maybe we should return
 // error?
 Core* mb = (Core*)pco->get_moab();
 
 double t1, t2, t3;
 t1 = MPI_Wtime();
 int primaryDim = mb->dimension_from_handle( *owned.rbegin() );
 int interfaceDim = primaryDim - 1; // should be 1 or 2
 Range sharedEdges;
 ErrorCode rval;
 
 std::vector< int > shprocs( MAX_SHARING_PROCS );
 std::vector< EntityHandle > shhandles( MAX_SHARING_PROCS );
 
 Tag gidTag = mb->globalId_tag();
 int np;
 unsigned char pstatus;
 
 std::multimap< int, int > extraGraphEdges;
 // std::map<int, int> adjCellsId;
 std::map< int, int > extraCellsProc;
 // if method is 2, no need to do the exchange for adjacent cells across partition boundary
 // these maps above will be empty for method 2 (geometry)
 if( 1 == met )
 {
 rval = pco->get_shared_entities( /*int other_proc*/ -1, sharedEdges, interfaceDim,
 /*const bool iface*/ true );MB_CHK_ERR( rval );
 
#ifdef VERBOSE
 std::cout << " on sender task " << pco->rank() << " number of shared interface cells " << sharedEdges.size()
 << "\n";
#endif
 // find to what processors we need to send the ghost info about the edge
 // first determine the local graph; what elements are adjacent to each cell in owned range
 // cells that are sharing a partition interface edge, are identified first, and form a map
 TupleList TLe; // tuple list for cells
 TLe.initialize( 2, 0, 1, 0, sharedEdges.size() ); // send to, id of adj cell, remote edge
 TLe.enableWriteAccess();
 
 std::map< EntityHandle, int > edgeToCell; // from local boundary edge to adjacent cell id
 // will be changed after
 for( Range::iterator eit = sharedEdges.begin(); eit != sharedEdges.end(); eit++ )
 {
 EntityHandle edge = *eit;
 // get the adjacent cell
 Range adjEnts;
 rval = mb->get_adjacencies( &edge, 1, primaryDim, false, adjEnts );MB_CHK_ERR( rval );
 if( adjEnts.size() > 0 )
 {
 EntityHandle adjCell = adjEnts[0];
 int gid;
 rval = mb->tag_get_data( gidTag, &adjCell, 1, &gid );MB_CHK_ERR( rval );
 rval = pco->get_sharing_data( edge, &shprocs[0], &shhandles[0], pstatus, np );MB_CHK_ERR( rval );
 int n = TLe.get_n();
 TLe.vi_wr[2 * n] = shprocs[0];
 TLe.vi_wr[2 * n + 1] = gid;
 TLe.vul_wr[n] = shhandles[0]; // the remote edge corresponding to shared edge
 edgeToCell[edge] = gid; // store the map between edge and local id of adj cell
 TLe.inc_n();
 }
 }
 
#ifdef VERBOSE
 std::stringstream ff2;
 ff2 << "TLe_" << pco->rank() << ".txt";
 TLe.print_to_file( ff2.str().c_str() );
#endif
 // send the data to the other processors:
 ( pco->proc_config().crystal_router() )->gs_transfer( 1, TLe, 0 );
 // on receiver side, each local edge will have the remote cell adjacent to it!
 
 int ne = TLe.get_n();
 for( int i = 0; i < ne; i++ )
 {
 int sharedProc = TLe.vi_rd[2 * i]; // this info is coming from here, originally
 int remoteCellID = TLe.vi_rd[2 * i + 1]; // this is the id of the remote cell, on sharedProc
 EntityHandle localCell = TLe.vul_rd[i]; // this is now local edge/face on this proc
 int localCellId = edgeToCell[localCell]; // this is the local cell adjacent to edge/face
 // now, we will need to add to the graph the pair <localCellId, remoteCellID>
 std::pair< int, int > extraAdj = std::make_pair( localCellId, remoteCellID );
 extraGraphEdges.insert( extraAdj );
 // adjCellsId [edgeToCell[localCell]] = remoteCellID;
 extraCellsProc[remoteCellID] = sharedProc;
#ifdef VERBOSE
 std::cout << "local ID " << edgeToCell[localCell] << " remote cell ID: " << remoteCellID << "\n";
#endif
 }
 }
 t2 = MPI_Wtime();
 if( rootSender ) std::cout << " time preparing the input for Zoltan:" << t2 - t1 << " seconds. \n";
 // so adj cells ids; need to call zoltan for parallel partition
#ifdef MOAB_HAVE_ZOLTAN
 ZoltanPartitioner* mbZTool = new ZoltanPartitioner( mb, pco );
 if( 1 <= met ) // partition in zoltan, either graph or geometric partitioner
 {
 std::map< int, Range > distribution; // how to distribute owned elements by processors in receiving groups
 // in how many tasks do we want to be distributed?
 int numNewPartitions = (int)receiverTasks.size();
 Range primaryCells = owned.subset_by_dimension( primaryDim );
 rval = mbZTool->partition_owned_cells( primaryCells, extraGraphEdges, extraCellsProc, numNewPartitions,
 distribution, met );MB_CHK_ERR( rval );
 for( std::map< int, Range >::iterator mit = distribution.begin(); mit != distribution.end(); mit++ )
 {
 int part_index = mit->first;
 assert( part_index < numNewPartitions );
 split_ranges[receiverTasks[part_index]] = mit->second;
 }
 }
 // delete now the partitioner
 delete mbZTool;
#endif
 t3 = MPI_Wtime();
 if( rootSender ) std::cout << " time spent by Zoltan " << t3 - t2 << " seconds. \n";
 return MB_SUCCESS;
}
 
// after map read, we need to know what entities we need to send to receiver
ErrorCode ParCommGraph::set_split_ranges( int comp,
 TupleList& TLBackToComp1,
 std::vector< int >& valuesComp1,
 int lenTag,
 Range& ents_of_interest,
 int /*type*/ )
{
 // settle split_ranges // same role as partitioning
 if( rootSender ) std::cout << " find split_ranges on component " << comp << " according to read map \n";
 // Vector to store element
 // with respective present index
 int n = TLBackToComp1.get_n();
 // third_method = true; // do not rely only on involved_IDs_map.size(); this can be 0 in some
 // cases
 std::map< int, std::set< int > > uniqueIDs;
 
 for( int i = 0; i < n; i++ )
 {
 int to_proc = TLBackToComp1.vi_wr[3 * i + 2];
 int globalId = TLBackToComp1.vi_wr[3 * i + 1];
 uniqueIDs[to_proc].insert( globalId );
 }
 
 for( int i = 0; i < (int)ents_of_interest.size(); i++ )
 {
 EntityHandle ent = ents_of_interest[i];
 for( int j = 0; j < lenTag; j++ )
 {
 int marker = valuesComp1[i * lenTag + j];
 for( auto mit = uniqueIDs.begin(); mit != uniqueIDs.end(); mit++ )
 {
 int proc = mit->first;
 std::set< int >& setIds = mit->second;
 if( setIds.find( marker ) != setIds.end() )
 {
 split_ranges[proc].insert( ent );
 }
 }
 }
 }
 
 return MB_SUCCESS;
}
 
// new methods to migrate mesh after reading map
ErrorCode ParCommGraph::form_tuples_to_migrate_mesh( Interface* mb,
 TupleList& TLv,
 TupleList& TLc,
 int type,
 int lenTagType1 )
{
 // we will always need GlobalID tag
 Tag gidtag = mb->globalId_tag();
 // for Type1, we need GLOBAL_DOFS tag;
 Tag gds;
 ErrorCode rval;
 if( type == 1 )
 {
 rval = mb->tag_get_handle( "GLOBAL_DOFS", gds );
 }
 // find vertices to be sent, for each of the split_ranges procs
 std::map< int, Range > verts_to_proc;
 int numv = 0, numc = 0;
 for( auto it = split_ranges.begin(); it != split_ranges.end(); it++ )
 {
 int to_proc = it->first;
 Range verts;
 if( type != 2 )
 {
 rval = mb->get_connectivity( it->second, verts );MB_CHK_ERR( rval );
 numc += (int)it->second.size();
 }
 else
 verts = it->second;
 verts_to_proc[to_proc] = verts;
 numv += (int)verts.size();
 }
 // first vertices:
 TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, 3 real coordinates
 TLv.enableWriteAccess();
 // use the global id of vertices for connectivity
 for( auto it = verts_to_proc.begin(); it != verts_to_proc.end(); it++ )
 {
 int to_proc = it->first;
 Range& verts = it->second;
 for( Range::iterator vit = verts.begin(); vit != verts.end(); ++vit )
 {
 EntityHandle v = *vit;
 int n = TLv.get_n(); // current size of tuple list
 TLv.vi_wr[2 * n] = to_proc; // send to processor
 
 rval = mb->tag_get_data( gidtag, &v, 1, &( TLv.vi_wr[2 * n + 1] ) );MB_CHK_ERR( rval );
 rval = mb->get_coords( &v, 1, &( TLv.vr_wr[3 * n] ) );MB_CHK_ERR( rval );
 TLv.inc_n(); // increment tuple list size
 }
 }
 if( type == 2 ) return MB_SUCCESS; // no need to fill TLc
 // to proc, ID cell, gdstag, nbv, id conn,
 int size_tuple = 2 + ( ( type != 1 ) ? 0 : lenTagType1 ) + 1 + 10; // 10 is the max number of vertices in cell
 
 std::vector< int > gdvals;
 
 TLc.initialize( size_tuple, 0, 0, 0, numc ); // to proc, GLOBAL ID, 3 real coordinates
 TLc.enableWriteAccess();
 for( auto it = split_ranges.begin(); it != split_ranges.end(); it++ )
 {
 int to_proc = it->first;
 Range& cells = it->second;
 for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
 {
 EntityHandle cell = *cit;
 int n = TLc.get_n(); // current size of tuple list
 TLc.vi_wr[size_tuple * n] = to_proc;
 int current_index = 2;
 rval = mb->tag_get_data( gidtag, &cell, 1, &( TLc.vi_wr[size_tuple * n + 1] ) );MB_CHK_ERR( rval );
 if( 1 == type )
 {
 rval = mb->tag_get_data( gds, &cell, 1, &( TLc.vi_wr[size_tuple * n + current_index] ) );MB_CHK_ERR( rval );
 current_index += lenTagType1;
 }
 // now get connectivity
 const EntityHandle* conn = NULL;
 int nnodes = 0;
 rval = mb->get_connectivity( cell, conn, nnodes );MB_CHK_ERR( rval );
 // fill nnodes:
 TLc.vi_wr[size_tuple * n + current_index] = nnodes;
 rval = mb->tag_get_data( gidtag, conn, nnodes, &( TLc.vi_wr[size_tuple * n + current_index + 1] ) );MB_CHK_ERR( rval );
 TLc.inc_n(); // increment tuple list size
 }
 }
 return MB_SUCCESS;
}
ErrorCode ParCommGraph::form_mesh_from_tuples( Interface* mb,
 TupleList& TLv,
 TupleList& TLc,
 int type,
 int lenTagType1,
 EntityHandle fset,
 Range& primary_ents,
 std::vector< int >& values_entities )
{
 // might need to fill also the split_range things
 // we will always need GlobalID tag
 Tag gidtag = mb->globalId_tag();
 // for Type1, we need GLOBAL_DOFS tag;
 Tag gds;
 ErrorCode rval;
 std::vector< int > def_val( lenTagType1, 0 );
 if( type == 1 )
 {
 // we may need to create this tag
 rval = mb->tag_get_handle( "GLOBAL_DOFS", lenTagType1, MB_TYPE_INTEGER, gds, MB_TAG_CREAT | MB_TAG_DENSE,
 &def_val[0] );MB_CHK_ERR( rval );
 }
 
 std::map< int, EntityHandle > vertexMap; //
 Range verts;
 // always form vertices and add them to the fset;
 int n = TLv.get_n();
 EntityHandle vertex;
 for( int i = 0; i < n; i++ )
 {
 int gid = TLv.vi_rd[2 * i + 1];
 if( vertexMap.find( gid ) == vertexMap.end() )
 {
 // need to form this vertex
 rval = mb->create_vertex( &( TLv.vr_rd[3 * i] ), vertex );MB_CHK_ERR( rval );
 vertexMap[gid] = vertex;
 verts.insert( vertex );
 rval = mb->tag_set_data( gidtag, &vertex, 1, &gid );MB_CHK_ERR( rval );
 // if point cloud,
 }
 if( 2 == type ) // point cloud, add it to the split_ranges ?
 {
 split_ranges[TLv.vi_rd[2 * i]].insert( vertexMap[gid] );
 }
 // todo : when to add the values_entities ?
 }
 rval = mb->add_entities( fset, verts );MB_CHK_ERR( rval );
 if( 2 == type )
 {
 primary_ents = verts;
 values_entities.resize( verts.size() ); // just get the ids of vertices
 rval = mb->tag_get_data( gidtag, verts, &values_entities[0] );MB_CHK_ERR( rval );
 return MB_SUCCESS;
 }
 
 n = TLc.get_n();
 int size_tuple = 2 + ( ( type != 1 ) ? 0 : lenTagType1 ) + 1 + 10; // 10 is the max number of vertices in cell
 
 EntityHandle new_element;
 Range cells;
 std::map< int, EntityHandle > cellMap; // do no tcreate one if it already exists, maybe from other processes
 for( int i = 0; i < n; i++ )
 {
 int from_proc = TLc.vi_rd[size_tuple * i];
 int globalIdEl = TLc.vi_rd[size_tuple * i + 1];
 if( cellMap.find( globalIdEl ) == cellMap.end() ) // need to create the cell
 {
 int current_index = 2;
 if( 1 == type ) current_index += lenTagType1;
 int nnodes = TLc.vi_rd[size_tuple * i + current_index];
 std::vector< EntityHandle > conn;
 conn.resize( nnodes );
 for( int j = 0; j < nnodes; j++ )
 {
 conn[j] = vertexMap[TLc.vi_rd[size_tuple * i + current_index + j + 1]];
 }
 //
 EntityType entType = MBQUAD;
 if( nnodes > 4 ) entType = MBPOLYGON;
 if( nnodes < 4 ) entType = MBTRI;
 rval = mb->create_element( entType, &conn[0], nnodes, new_element );MB_CHK_SET_ERR( rval, "can't create new element " );
 cells.insert( new_element );
 cellMap[globalIdEl] = new_element;
 rval = mb->tag_set_data( gidtag, &new_element, 1, &globalIdEl );MB_CHK_SET_ERR( rval, "can't set global id tag on cell " );
 if( 1 == type )
 {
 // set the gds tag
 rval = mb->tag_set_data( gds, &new_element, 1, &( TLc.vi_rd[size_tuple * i + 2] ) );MB_CHK_SET_ERR( rval, "can't set gds tag on cell " );
 }
 }
 split_ranges[from_proc].insert( cellMap[globalIdEl] );
 }
 rval = mb->add_entities( fset, cells );MB_CHK_ERR( rval );
 primary_ents = cells;
 if( 1 == type )
 {
 values_entities.resize( lenTagType1 * primary_ents.size() );
 rval = mb->tag_get_data( gds, primary_ents, &values_entities[0] );MB_CHK_ERR( rval );
 }
 else // type == 3
 {
 values_entities.resize( primary_ents.size() ); // just get the ids !
 rval = mb->tag_get_data( gidtag, primary_ents, &values_entities[0] );MB_CHK_ERR( rval );
 }
 return MB_SUCCESS;
}
 
// at this moment, each sender task has split_ranges formed;
// we need to aggregate that info and send it to receiver
ErrorCode ParCommGraph::send_graph_partition( ParallelComm* pco, MPI_Comm jcomm )
{
 // first, accumulate the info to root of sender; use gatherv
 // first, accumulate number of receivers from each sender, to the root receiver
 int numberReceivers =
 (int)split_ranges.size(); // these are ranges of elements to be sent to each receiver, from this task
 int nSenders = (int)senderTasks.size(); //
 // this sender will have to send to this many receivers
 std::vector< int > displs( 1 ); // displacements for gatherv
 std::vector< int > counts( 1 );
 if( is_root_sender() )
 {
 displs.resize( nSenders + 1 );
 counts.resize( nSenders );
 }
 
 int ierr = MPI_Gather( &numberReceivers, 1, MPI_INT, &counts[0], 1, MPI_INT, 0, pco->comm() );
 if( ierr != MPI_SUCCESS ) return MB_FAILURE;
 // compute now displacements
 if( is_root_sender() )
 {
 displs[0] = 0;
 for( int k = 0; k < nSenders; k++ )
 {
 displs[k + 1] = displs[k] + counts[k];
 }
 }
 std::vector< int > buffer;
 if( is_root_sender() ) buffer.resize( displs[nSenders] ); // the last one will have the total count now
 
 std::vector< int > recvs;
 for( std::map< int, Range >::iterator mit = split_ranges.begin(); mit != split_ranges.end(); mit++ )
 {
 recvs.push_back( mit->first );
 }
 ierr =
 MPI_Gatherv( &recvs[0], numberReceivers, MPI_INT, &buffer[0], &counts[0], &displs[0], MPI_INT, 0, pco->comm() );
 if( ierr != MPI_SUCCESS ) return MB_FAILURE;
 
 // now form recv_graph map; points from the
 // now form the graph to be sent to the other side; only on root
 if( is_root_sender() )
 {
#ifdef GRAPH_INFO
 std::ofstream dbfileSender;
 std::stringstream outf;
 outf << "S_" << compid1 << "_R_" << compid2 << "_SenderGraph.txt";
 dbfileSender.open( outf.str().c_str() );
 dbfileSender << " number senders: " << nSenders << "\n";
 dbfileSender << " senderRank \treceivers \n";
 for( int k = 0; k < nSenders; k++ )
 {
 int indexInBuff = displs[k];
 int senderTask = senderTasks[k];
 dbfileSender << senderTask << "\t\t";
 for( int j = 0; j < counts[k]; j++ )
 {
 int recvTask = buffer[indexInBuff + j];
 dbfileSender << recvTask << " ";
 }
 dbfileSender << "\n";
 }
 dbfileSender.close();
#endif
 for( int k = 0; k < nSenders; k++ )
 {
 int indexInBuff = displs[k];
 int senderTask = senderTasks[k];
 for( int j = 0; j < counts[k]; j++ )
 {
 int recvTask = buffer[indexInBuff + j];
 recv_graph[recvTask].push_back( senderTask ); // this will be packed and sent to root receiver, with
 // nonblocking send
 }
 }
#ifdef GRAPH_INFO
 std::ofstream dbfile;
 std::stringstream outf2;
 outf2 << "S_" << compid1 << "_R_" << compid2 << "_RecvGraph.txt";
 dbfile.open( outf2.str().c_str() );
 dbfile << " number receivers: " << recv_graph.size() << "\n";
 dbfile << " receiverRank \tsenders \n";
 for( std::map< int, std::vector< int > >::iterator mit = recv_graph.begin(); mit != recv_graph.end(); mit++ )
 {
 int recvTask = mit->first;
 std::vector< int >& senders = mit->second;
 dbfile << recvTask << "\t\t";
 for( std::vector< int >::iterator vit = senders.begin(); vit != senders.end(); vit++ )
 dbfile << *vit << " ";
 dbfile << "\n";
 }
 dbfile.close();
#endif
 // this is the same as trivial partition
 ErrorCode rval = send_graph( jcomm );MB_CHK_ERR( rval );
 }
 
 return MB_SUCCESS;
}
// method to expose local graph info: sender id, receiver id, sizes of elements to send, after or
// before intersection
ErrorCode ParCommGraph::dump_comm_information( std::string prefix, int is_send )
{
 //
 if( -1 != rankInGroup1 && 1 == is_send ) // it is a sender task
 {
 std::ofstream dbfile;
 std::stringstream outf;
 outf << prefix << "_sender_" << rankInGroup1 << "_joint" << rankInJoin << "_type_" << (int)graph_type << ".txt";
 dbfile.open( outf.str().c_str() );
 
 if( graph_type == COVERAGE )
 {
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin();
 mit != involved_IDs_map.end(); mit++ )
 {
 int receiver_proc = mit->first;
 std::vector< int >& eids = mit->second;
 dbfile << "receiver: " << receiver_proc << " size:" << eids.size() << "\n";
 }
 }
 else if( graph_type == INITIAL_MIGRATE ) // just after migration
 {
 for( std::map< int, Range >::iterator mit = split_ranges.begin(); mit != split_ranges.end(); mit++ )
 {
 int receiver_proc = mit->first;
 Range& eids = mit->second;
 dbfile << "receiver: " << receiver_proc << " size:" << eids.size() << "\n";
 }
 }
 else if( graph_type == DOF_BASED ) // just after migration, or from computeGraph
 {
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin();
 mit != involved_IDs_map.end(); mit++ )
 {
 int receiver_proc = mit->first;
 dbfile << "receiver: " << receiver_proc << " size:" << mit->second.size() << "\n";
 }
 }
 dbfile.close();
 }
 if( -1 != rankInGroup2 && 0 == is_send ) // it is a receiver task
 {
 std::ofstream dbfile;
 std::stringstream outf;
 outf << prefix << "_receiver_" << rankInGroup2 << "_joint" << rankInJoin << "_type_" << (int)graph_type
 << ".txt";
 dbfile.open( outf.str().c_str() );
 
 if( graph_type == COVERAGE )
 {
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin();
 mit != involved_IDs_map.end(); mit++ )
 {
 int sender_proc = mit->first;
 std::vector< int >& eids = mit->second;
 dbfile << "sender: " << sender_proc << " size:" << eids.size() << "\n";
 }
 }
 else if( graph_type == INITIAL_MIGRATE ) // just after migration
 {
 for( std::map< int, Range >::iterator mit = split_ranges.begin(); mit != split_ranges.end(); mit++ )
 {
 int sender_proc = mit->first;
 Range& eids = mit->second;
 dbfile << "sender: " << sender_proc << " size:" << eids.size() << "\n";
 }
 }
 else if( graph_type == DOF_BASED ) // just after migration
 {
 for( std::map< int, std::vector< int > >::iterator mit = involved_IDs_map.begin();
 mit != involved_IDs_map.end(); mit++ )
 {
 int sender_proc = mit->first;
 dbfile << "receiver: " << sender_proc << " size:" << mit->second.size() << "\n";
 }
 }
 dbfile.close();
 }
 return MB_SUCCESS;
}
} // namespace moab