dune-grid  3.0-git
parmetisgridpartitioner.hh
Go to the documentation of this file.
1 #ifndef DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH
2 #define DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH
3 
8 #include <algorithm>
9 #include <vector>
10 
11 #include <dune/common/parallel/mpihelper.hh>
12 #include <dune/common/exceptions.hh>
13 
14 #include <dune/geometry/referenceelements.hh>
15 
18 
19 #if HAVE_PARMETIS
20 
21 #include <parmetis.h>
22 
23 // only enable for ParMETIS because the implementation uses functions that
24 // are not emulated by scotch
25 #ifdef PARMETIS_MAJOR_VERSION
26 
27 namespace Dune
28 {
29 
34  template<class GridView>
35  struct ParMetisGridPartitioner {
36 
37  // define index type as provided by ParMETIS
38 #if PARMETIS_MAJOR_VERSION > 3
39  typedef idx_t idx_type;
40  typedef ::real_t real_type;
41 #else
42  typedef int idx_type;
43  typedef float real_type;
44 #endif // PARMETIS_MAJOR_VERSION > 3
45 
46  typedef typename GridView::template Codim<0>::Iterator ElementIterator;
47  typedef typename GridView::template Codim<0>::template Partition<Interior_Partition>::Iterator InteriorElementIterator;
48  typedef typename GridView::IntersectionIterator IntersectionIterator;
49 
50  enum {
51  dimension = GridView::dimension
52  };
53 
54 
65  static std::vector<unsigned> partition(const GridView& gv, const Dune::MPIHelper& mpihelper) {
66  const unsigned numElements = gv.size(0);
67 
68  std::vector<unsigned> part(numElements);
69 
70  // Setup parameters for ParMETIS
71  idx_type wgtflag = 0; // we don't use weights
72  idx_type numflag = 0; // we are using C-style arrays
73  idx_type ncon = 1; // number of balance constraints
74  idx_type ncommonnodes = 2; // number of nodes elements must have in common to be considered adjacent to each other
75  idx_type options[4] = {0, 0, 0, 0}; // use default values for random seed, output and coupling
76  idx_type edgecut; // will store number of edges cut by partition
77  idx_type nparts = mpihelper.size(); // number of parts equals number of processes
78  std::vector<real_type> tpwgts(ncon*nparts, 1./nparts); // load per subdomain and weight (same load on every process)
79  std::vector<real_type> ubvec(ncon, 1.05); // weight tolerance (same weight tolerance for every weight there is)
80 
81  // The difference elmdist[i+1] - elmdist[i] is the number of nodes that are on process i
82  std::vector<idx_type> elmdist(nparts+1);
83  elmdist[0] = 0;
84  std::fill(elmdist.begin()+1, elmdist.end(), gv.size(0)); // all elements are on process zero
85 
86  // Create and fill arrays "eptr", where eptr[i] is the number of vertices that belong to the i-th element, and
87  // "eind" contains the vertex-numbers of the i-the element in eind[eptr[i]] to eind[eptr[i+1]-1]
88  std::vector<idx_type> eptr, eind;
89  int numVertices = 0;
90  eptr.push_back(numVertices);
91 
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);
94 
95  numVertices += curNumVertices;
96  eptr.push_back(numVertices);
97 
98  for (size_t k = 0; k < curNumVertices; ++k)
99  eind.push_back(gv.indexSet().subIndex(*eIt, k, dimension));
100  }
101 
102  // Partition mesh using ParMETIS
103  if (0 == mpihelper.rank()) {
104  MPI_Comm comm = Dune::MPIHelper::getLocalCommunicator();
105 
106 #if PARMETIS_MAJOR_VERSION >= 4
107  const int OK =
108 #endif
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);
112 
113 #if PARMETIS_MAJOR_VERSION >= 4
114  if (OK != METIS_OK)
115  DUNE_THROW(Dune::Exception, "ParMETIS returned an error code.");
116 #endif
117  }
118 
119  return part;
120  }
121 
133  static std::vector<unsigned> repartition(const GridView& gv, const Dune::MPIHelper& mpihelper, real_type& itr = 1000) {
134 
135  // Create global index map
136  GlobalIndexSet<GridView> globalIndex(gv,0);
137 
138  int numElements = std::distance(gv.template begin<0, Interior_Partition>(),
139  gv.template end<0, Interior_Partition>());
140 
141  std::vector<unsigned> interiorPart(numElements);
142 
143  // Setup parameters for ParMETIS
144  idx_type wgtflag = 0; // we don't use weights
145  idx_type numflag = 0; // we are using C-style arrays
146  idx_type ncon = 1; // number of balance constraints
147  idx_type options[4] = {0, 0, 0, 0}; // use default values for random seed, output and coupling
148  idx_type edgecut; // will store number of edges cut by partition
149  idx_type nparts = mpihelper.size(); // number of parts equals number of processes
150  std::vector<real_type> tpwgts(ncon*nparts, 1./nparts); // load per subdomain and weight (same load on every process)
151  std::vector<real_type> ubvec(ncon, 1.05); // weight tolerance (same weight tolerance for every weight there is)
152 
153  MPI_Comm comm = Dune::MPIHelper::getCommunicator();
154 
155  // Make the number of interior elements of each processor available to all processors
156  std::vector<int> offset(gv.comm().size());
157  std::fill(offset.begin(), offset.end(), 0);
158 
159  gv.comm().template allgather<int>(&numElements, 1, offset.data());
160 
161  // The difference vtxdist[i+1] - vtxdist[i] is the number of elements that are on process i
162  std::vector<idx_type> vtxdist(gv.comm().size()+1);
163  vtxdist[0] = 0;
164 
165  for (unsigned int i=1; i<vtxdist.size(); ++i)
166  vtxdist[i] = vtxdist[i-1] + offset[i-1];
167 
168  // Set up element adjacency lists
169  std::vector<idx_type> xadj, adjncy;
170  xadj.push_back(0);
171 
172  for (InteriorElementIterator eIt = gv.template begin<0, Interior_Partition>(); eIt != gv.template end<0, Interior_Partition>(); ++eIt)
173  {
174  size_t numNeighbors = 0;
175 
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()));
179 
180  ++numNeighbors;
181  }
182  }
183 
184  xadj.push_back(xadj.back() + numNeighbors);
185  }
186 
187 #if PARMETIS_MAJOR_VERSION >= 4
188  const int OK =
189 #endif
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);
193 
194 #if PARMETIS_MAJOR_VERSION >= 4
195  if (OK != METIS_OK)
196  DUNE_THROW(Dune::Exception, "ParMETIS returned error code " << OK);
197 #endif
198 
199  // At this point, interiorPart contains a target rank for each interior element, and they are sorted
200  // by the order in which the grid view traverses them. Now we need to do two things:
201  // a) Add additional dummy entries for the ghost elements
202  // b) Use the element index for the actual ordering. Since there may be different types of elements,
203  // we cannot use the index set directly, but have to go through a Mapper.
204 
205  typedef MultipleCodimMultipleGeomTypeMapper<GridView, MCMGElementLayout> ElementMapper;
206  ElementMapper elementMapper(gv);
207 
208  std::vector<unsigned int> part(gv.size(0));
209  std::fill(part.begin(), part.end(), 0);
210  unsigned int c = 0;
211  for (InteriorElementIterator eIt = gv.template begin<0, Interior_Partition>();
212  eIt != gv.template end<0, Interior_Partition>();
213  ++eIt)
214  {
215  part[elementMapper.index(*eIt)] = interiorPart[c++];
216  }
217  return part;
218  }
219  };
220 
221 } // namespace Dune
222 
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."
225 #endif
226 
227 #else // HAVE_PARMETIS
228 #warning "PARMETIS was not found, please check your configuration"
229 #endif
230 
231 #endif // DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH
Dune
Include standard header files.
Definition: agrid.hh:59
Dune::GridView::IntersectionIterator
Traits ::IntersectionIterator IntersectionIterator
type of the intersection iterator
Definition: common/gridview.hh:87
mcmgmapper.hh
Mapper for multiple codim and multiple geometry types.
Dune::GridView::dimension
The dimension of the grid.
Definition: common/gridview.hh:128
globalindexset.hh
Provides a globally unique index for all entities of a distributed Dune grid.