parallel/MeshParalleliser.cpp

#include "DAPTA.hpp"

#include <vector>

namespace DAPTA { // Define namespace DAPTA

template <typename Traits>
MeshParalleliser<Traits>::MeshParalleliser(MPI_Comm _comm, Mesh<Traits> *_mesh)
   : communicator(_comm), mesh(_mesh) {
   assert(_comm != MPI_COMM_NULL);
    assert(_mesh != NULL);
   decomposer = new DomainDecomposerType();
   
   std::vector<VertexHandle> elemCentres = mesh->GetCentres();
   decomposer->loadVertices(elemCentres);
   int numProcs, myRank;
   MPI_Comm_size(communicator, &numProcs);
   MPI_Comm_rank(communicator, &myRank);
   
   std::vector<int> myElems = decomposer->decompose(numProcs, myRank);
   
   mesh->SetActive(myElems);
   mesh->newPrefix = (myRank+1);
   
   // TODO: Code to make mesh "delete" the elements it doesn't need anymore
}

template <typename Traits>
MeshParalleliser<Traits>::~MeshParalleliser() {
    delete coupler;
    delete exchanger;
   delete decomposer;
}

template <typename Traits>
void MeshParalleliser<Traits>::initialiseCommunication() {
   std::set<EdgeElementType> tmpEdgesSet;
   for(typename Mesh<Traits>::ElementMapIterator e_iter=mesh->ElementMapBegin(); e_iter!=mesh->ElementMapEnd(); ++e_iter) {
      std::vector<EdgeElementType> edges_vec=(*e_iter)->GetAllEdges();
      std::copy(edges_vec.begin(), edges_vec.end(), std::inserter(tmpEdgesSet, tmpEdgesSet.end()));
   }
   std::copy(tmpEdgesSet.begin(), tmpEdgesSet.end(), std::back_inserter(edgesConnected));
   
   std::set<VertexHandle> tmpVertexSet(mesh->VertexMapBegin(), mesh->VertexMapEnd());
   verticesConnected.reserve(tmpVertexSet.size());
   std::copy(tmpVertexSet.begin(), tmpVertexSet.end(), std::back_inserter(verticesConnected));
   
   std::vector<typename VertexType::CoordType> coordinates;
   for(typename Mesh<Traits>::VertexMapIterator v_iter=mesh->VertexMapBegin(); v_iter!=mesh->VertexMapEnd(); ++v_iter) {
      std::copy((*v_iter)->coords.begin(), (*v_iter)->coords.end(), std::back_inserter(coordinates));
   }
   
   coupler = new pico::DomainCoupler<VertexType>(communicator);
   coupler->registerSubDomain(pico::RED, &(coordinates[0]), coordinates.size()/VertexType::dimension, 0.0);
   coupler->prepareCommunicationsIntra(pico::RED);
   
   std::vector<int> globalIDs;
   globalIDs.reserve(verticesConnected.size());
   std::transform(verticesConnected.begin(), verticesConnected.end(), std::back_inserter(globalIDs), GetID<VertexType>());
   
   exchanger = new pico::FieldExchanger<VertexType>(coupler, globalIDs.begin(), globalIDs.size());
}

template <typename Traits>
void MeshParalleliser<Traits>::exchange() {
   /*int myRank;
   MPI_Comm_rank(communicator, &myRank);
   
   std::vector<int> test;
   for(int i=0; i<verticesConnected.size(); i++)
      test.push_back(verticesConnected.at(i)->globalID);
   exchanger->template exchange<1, int, std::plus>(&(test[0]), test.size());
   for(int i=0; i<verticesConnected.size(); i++) {
      std::cout << "client(" << myRank << "): " << verticesConnected.at(i)->globalID << ":" << test.at(i) << std::endl;
   }
   
   MPI_Barrier(MPI_COMM_WORLD);*/
   
   int pad_size=8;
   std::vector<int> link_states;
   link_states.reserve(edgesConnected.size() * pad_size);
   
   for(_EdgeIt edge_iter=edgesConnected.begin(); edge_iter!=edgesConnected.end(); ++edge_iter) {
      mesh->edges[*edge_iter]->GetLinkStateArray(std::back_inserter(link_states));
      
      for(int i=0; i<(pad_size-mesh->edges[*edge_iter]->size()); i++) {
         link_states.push_back(0);
      }
   }
   
   exchanger->template exchange<8, int, LinkStateJoin>(&(link_states[0]), link_states.size());
   
   int i=0;
   for(_EdgeIt edge_iter=edgesConnected.begin(); edge_iter!=edgesConnected.end(); ++edge_iter) {
      mesh->edges[*edge_iter]->SetLinkStateArray(link_states.begin()+(i*pad_size));
      i++;
   }
}

} // namespace DAPTA

---