1 #ifndef DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH
2 #define DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH
11 #include <dune/common/parallel/mpihelper.hh>
12 #include <dune/common/exceptions.hh>
14 #include <dune/geometry/referenceelements.hh>
25 #ifdef PARMETIS_MAJOR_VERSION
34 template<
class Gr
idView>
35 struct ParMetisGridPartitioner {
38 #if PARMETIS_MAJOR_VERSION > 3
39 typedef idx_t idx_type;
40 typedef ::real_t real_type;
43 typedef float real_type;
44 #endif // PARMETIS_MAJOR_VERSION > 3
46 typedef typename GridView::template Codim<0>::Iterator ElementIterator;
47 typedef typename GridView::template Codim<0>::template Partition<Interior_Partition>::Iterator InteriorElementIterator;
65 static std::vector<unsigned> partition(
const GridView& gv,
const Dune::MPIHelper& mpihelper) {
66 const unsigned numElements = gv.size(0);
68 std::vector<unsigned> part(numElements);
74 idx_type ncommonnodes = 2;
75 idx_type options[4] = {0, 0, 0, 0};
77 idx_type nparts = mpihelper.size();
78 std::vector<real_type> tpwgts(ncon*nparts, 1./nparts);
79 std::vector<real_type> ubvec(ncon, 1.05);
82 std::vector<idx_type> elmdist(nparts+1);
84 std::fill(elmdist.begin()+1, elmdist.end(), gv.size(0));
88 std::vector<idx_type> eptr, eind;
90 eptr.push_back(numVertices);
92 for (InteriorElementIterator eIt = gv.template begin<0, Interior_Partition>(); eIt != gv.template end<0, Interior_Partition>(); ++eIt) {
93 const size_t curNumVertices = ReferenceElements<double, dimension>::general(eIt->type()).size(dimension);
95 numVertices += curNumVertices;
96 eptr.push_back(numVertices);
98 for (
size_t k = 0; k < curNumVertices; ++k)
99 eind.push_back(gv.indexSet().subIndex(*eIt, k, dimension));
103 if (0 == mpihelper.rank()) {
104 MPI_Comm comm = Dune::MPIHelper::getLocalCommunicator();
106 #if PARMETIS_MAJOR_VERSION >= 4
109 ParMETIS_V3_PartMeshKway(elmdist.data(), eptr.data(), eind.data(), NULL, &wgtflag, &numflag,
110 &ncon, &ncommonnodes, &nparts, tpwgts.data(), ubvec.data(),
111 options, &edgecut, reinterpret_cast<idx_type*>(part.data()), &comm);
113 #if PARMETIS_MAJOR_VERSION >= 4
115 DUNE_THROW(Dune::Exception,
"ParMETIS returned an error code.");
133 static std::vector<unsigned> repartition(
const GridView& gv,
const Dune::MPIHelper& mpihelper, real_type& itr = 1000) {
136 GlobalIndexSet<GridView> globalIndex(gv,0);
138 int numElements = std::distance(gv.template begin<0, Interior_Partition>(),
139 gv.template end<0, Interior_Partition>());
141 std::vector<unsigned> interiorPart(numElements);
144 idx_type wgtflag = 0;
145 idx_type numflag = 0;
147 idx_type options[4] = {0, 0, 0, 0};
149 idx_type nparts = mpihelper.size();
150 std::vector<real_type> tpwgts(ncon*nparts, 1./nparts);
151 std::vector<real_type> ubvec(ncon, 1.05);
153 MPI_Comm comm = Dune::MPIHelper::getCommunicator();
156 std::vector<int> offset(gv.comm().size());
157 std::fill(offset.begin(), offset.end(), 0);
159 gv.comm().template allgather<int>(&numElements, 1, offset.data());
162 std::vector<idx_type> vtxdist(gv.comm().size()+1);
165 for (
unsigned int i=1; i<vtxdist.size(); ++i)
166 vtxdist[i] = vtxdist[i-1] + offset[i-1];
169 std::vector<idx_type> xadj, adjncy;
172 for (InteriorElementIterator eIt = gv.template begin<0, Interior_Partition>(); eIt != gv.template end<0, Interior_Partition>(); ++eIt)
174 size_t numNeighbors = 0;
176 for (IntersectionIterator iIt = gv.template ibegin(*eIt); iIt != gv.template iend(*eIt); ++iIt) {
177 if (iIt->neighbor()) {
178 adjncy.push_back(globalIndex.index(iIt->outside()));
184 xadj.push_back(xadj.back() + numNeighbors);
187 #if PARMETIS_MAJOR_VERSION >= 4
190 ParMETIS_V3_AdaptiveRepart(vtxdist.data(), xadj.data(), adjncy.data(), NULL, NULL, NULL,
191 &wgtflag, &numflag, &ncon, &nparts, tpwgts.data(), ubvec.data(),
192 &itr, options, &edgecut, reinterpret_cast<idx_type*>(interiorPart.data()), &comm);
194 #if PARMETIS_MAJOR_VERSION >= 4
196 DUNE_THROW(Dune::Exception,
"ParMETIS returned error code " << OK);
205 typedef MultipleCodimMultipleGeomTypeMapper<GridView, MCMGElementLayout> ElementMapper;
206 ElementMapper elementMapper(gv);
208 std::vector<unsigned int> part(gv.size(0));
209 std::fill(part.begin(), part.end(), 0);
211 for (InteriorElementIterator eIt = gv.template begin<0, Interior_Partition>();
212 eIt != gv.template end<0, Interior_Partition>();
215 part[elementMapper.index(*eIt)] = interiorPart[c++];
223 #else // PARMETIS_MAJOR_VERSION
224 #warning "You seem to be using the ParMETIS emulation layer of scotch, which does not work with this file."
227 #else // HAVE_PARMETIS
228 #warning "PARMETIS was not found, please check your configuration"
231 #endif // DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH