metis_partitioner.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ Includes -----------------------------------
21 
22 // Local Includes -----------------------------------
23 #include "libmesh/libmesh_config.h"
24 #include "libmesh/mesh_base.h"
27 #include "libmesh/elem.h"
29 #include "libmesh/error_vector.h"
30 #include "libmesh/vectormap.h"
32 
33 #ifdef LIBMESH_HAVE_METIS
34 // MIPSPro 7.4.2 gets confused about these nested namespaces
35 # ifdef __sgi
36 # include <cstdarg>
37 # endif
38  namespace Metis {
39  extern "C" {
40 # include "libmesh/ignore_warnings.h"
41 # include "metis.h"
43  }
44  }
45 #else
46 # include "libmesh/sfc_partitioner.h"
47 #endif
48 
49 
50 
51 
52 namespace libMesh
53 {
54 
55 
56 // ------------------------------------------------------------
57 // MetisPartitioner implementation
59  const unsigned int n_pieces)
60 {
61  libmesh_assert_greater (n_pieces, 0);
62  libmesh_assert (mesh.is_serial());
63 
64  // Check for an easy return
65  if (n_pieces == 1)
66  {
67  this->single_partition (mesh);
68  return;
69  }
70 
71 // What to do if the Metis library IS NOT present
72 #ifndef LIBMESH_HAVE_METIS
73 
74  libmesh_here();
75  libMesh::err << "ERROR: The library has been built without" << std::endl
76  << "Metis support. Using a space-filling curve" << std::endl
77  << "partitioner instead!" << std::endl;
78 
79  SFCPartitioner sfcp;
80 
81  sfcp.partition (mesh, n_pieces);
82 
83 // What to do if the Metis library IS present
84 #else
85 
86  START_LOG("partition()", "MetisPartitioner");
87 
88  const dof_id_type n_active_elem = mesh.n_active_elem();
89 
90  // build the graph
91  // std::vector<int> options(5);
92  std::vector<int> vwgt(n_active_elem);
93  std::vector<int> part(n_active_elem);
94 
95  int
96  n = static_cast<int>(n_active_elem), // number of "nodes" (elements)
97  // in the graph
98 // wgtflag = 2, // weights on vertices only,
99 // // none on edges
100 // numflag = 0, // C-style 0-based numbering
101  nparts = static_cast<int>(n_pieces), // number of subdomains to create
102  edgecut = 0; // the numbers of edges cut by the
103  // resulting partition
104 
105  // Set the options
106  // options[0] = 0; // use default options
107 
108  // Metis will only consider the active elements.
109  // We need to map the active element ids into a
110  // contiguous range. Further, we want the unique range indexing to be
111  // independednt of the element ordering, otherwise a circular dependency
112  // can result in which the partitioning depends on the ordering which
113  // depends on the partitioning...
114  vectormap<dof_id_type, dof_id_type> global_index_map;
115  global_index_map.reserve (n_active_elem);
116 
117  {
118  std::vector<dof_id_type> global_index;
119 
122 
125  it, end, global_index);
126 
127  libmesh_assert_equal_to (global_index.size(), n_active_elem);
128 
129  for (std::size_t cnt=0; it != end; ++it)
130  {
131  const Elem *elem = *it;
132 
133  global_index_map.insert (std::make_pair(elem->id(), global_index[cnt++]));
134  }
135  libmesh_assert_equal_to (global_index_map.size(), n_active_elem);
136  }
137 
138 
139  // Invoke METIS, but only on processor 0.
140  // Then broadcast the resulting decomposition
141  if (mesh.processor_id() == 0)
142  {
143  METIS_CSR_Graph csr_graph;
144 
145  csr_graph.offsets.resize(n_active_elem+1, 0);
146 
147  // Local scope for these
148  {
149  // build the graph in CSR format. Note that
150  // the edges in the graph will correspond to
151  // face neighbors
152  std::vector<const Elem*> neighbors_offspring;
153 
155  const MeshBase::element_iterator elem_end = mesh.active_elements_end();
156 
157  std::size_t graph_size=0;
158 
159  // (1) first pass - get the row sizes for each element by counting the number
160  // of face neighbors. Also populate the vwght array if necessary
161  for (; elem_it != elem_end; ++elem_it)
162  {
163  const Elem* elem = *elem_it;
164 
165  const dof_id_type elem_global_index =
166  global_index_map[elem->id()];
167 
168  libmesh_assert_less (elem_global_index, vwgt.size());
169 
170  // maybe there is a better weight?
171  // The weight is used to define what a balanced graph is
172  if(!_weights)
173  vwgt[elem_global_index] = elem->n_nodes();
174  else
175  vwgt[elem_global_index] = static_cast<int>((*_weights)[elem->id()]);
176 
177  unsigned int num_neighbors = 0;
178 
179  // Loop over the element's neighbors. An element
180  // adjacency corresponds to a face neighbor
181  for (unsigned int ms=0; ms<elem->n_neighbors(); ms++)
182  {
183  const Elem* neighbor = elem->neighbor(ms);
184 
185  if (neighbor != NULL)
186  {
187  // If the neighbor is active treat it
188  // as a connection
189  if (neighbor->active())
190  num_neighbors++;
191 
192 #ifdef LIBMESH_ENABLE_AMR
193 
194  // Otherwise we need to find all of the
195  // neighbor's children that are connected to
196  // us and add them
197  else
198  {
199  // The side of the neighbor to which
200  // we are connected
201  const unsigned int ns =
202  neighbor->which_neighbor_am_i (elem);
203  libmesh_assert_less (ns, neighbor->n_neighbors());
204 
205  // Get all the active children (& grandchildren, etc...)
206  // of the neighbor.
207  neighbor->active_family_tree (neighbors_offspring);
208 
209  // Get all the neighbor's children that
210  // live on that side and are thus connected
211  // to us
212  for (unsigned int nc=0; nc<neighbors_offspring.size(); nc++)
213  {
214  const Elem* child =
215  neighbors_offspring[nc];
216 
217  // This does not assume a level-1 mesh.
218  // Note that since children have sides numbered
219  // coincident with the parent then this is a sufficient test.
220  if (child->neighbor(ns) == elem)
221  {
222  libmesh_assert (child->active());
223  num_neighbors++;
224  }
225  }
226  }
227 
228 #endif /* ifdef LIBMESH_ENABLE_AMR */
229 
230  }
231  }
232 
233  csr_graph.prep_n_nonzeros(elem_global_index, num_neighbors);
234  graph_size += num_neighbors;
235  }
236 
237  csr_graph.prepare_for_use();
238 
239  // (2) second pass - fill the compressed adjacency array
240  elem_it = mesh.active_elements_begin();
241 
242  for (; elem_it != elem_end; ++elem_it)
243  {
244  const Elem* elem = *elem_it;
245 
246  const dof_id_type elem_global_index =
247  global_index_map[elem->id()];
248 
249  unsigned int connection=0;
250 
251  // Loop over the element's neighbors. An element
252  // adjacency corresponds to a face neighbor
253  for (unsigned int ms=0; ms<elem->n_neighbors(); ms++)
254  {
255  const Elem* neighbor = elem->neighbor(ms);
256 
257  if (neighbor != NULL)
258  {
259  // If the neighbor is active treat it
260  // as a connection
261  if (neighbor->active())
262  csr_graph(elem_global_index, connection++) = global_index_map[neighbor->id()];
263 
264 #ifdef LIBMESH_ENABLE_AMR
265 
266  // Otherwise we need to find all of the
267  // neighbor's children that are connected to
268  // us and add them
269  else
270  {
271  // The side of the neighbor to which
272  // we are connected
273  const unsigned int ns =
274  neighbor->which_neighbor_am_i (elem);
275  libmesh_assert_less (ns, neighbor->n_neighbors());
276 
277  // Get all the active children (& grandchildren, etc...)
278  // of the neighbor.
279  neighbor->active_family_tree (neighbors_offspring);
280 
281  // Get all the neighbor's children that
282  // live on that side and are thus connected
283  // to us
284  for (unsigned int nc=0; nc<neighbors_offspring.size(); nc++)
285  {
286  const Elem* child =
287  neighbors_offspring[nc];
288 
289  // This does not assume a level-1 mesh.
290  // Note that since children have sides numbered
291  // coincident with the parent then this is a sufficient test.
292  if (child->neighbor(ns) == elem)
293  {
294  libmesh_assert (child->active());
295 
296  csr_graph(elem_global_index, connection++) = global_index_map[child->id()];
297  }
298  }
299  }
300 
301 #endif /* ifdef LIBMESH_ENABLE_AMR */
302 
303  }
304  }
305  }
306 
307  // We create a non-empty vals for a disconnected graph, to
308  // work around a segfault from METIS.
309  libmesh_assert_equal_to (csr_graph.vals.size(),
310  std::max(graph_size,std::size_t(1)));
311  } // done building the graph
312 
313  int ncon = 1;
314 
315  // Select which type of partitioning to create
316 
317  // Use recursive if the number of partitions is less than or equal to 8
318  if (n_pieces <= 8)
319  Metis::METIS_PartGraphRecursive(&n, &ncon, &csr_graph.offsets[0], &csr_graph.vals[0], &vwgt[0], NULL,
320  NULL, &nparts, NULL, NULL, NULL,
321  &edgecut, &part[0]);
322 
323  // Otherwise use kway
324  else
325  Metis::METIS_PartGraphKway(&n, &ncon, &csr_graph.offsets[0], &csr_graph.vals[0], &vwgt[0], NULL,
326  NULL, &nparts, NULL, NULL, NULL,
327  &edgecut, &part[0]);
328 
329  } // end processor 0 part
330 
331  // Broadcase the resutling partition
332  mesh.comm().broadcast(part);
333 
334  // Assign the returned processor ids. The part array contains
335  // the processor id for each active element, but in terms of
336  // the contiguous indexing we defined above
337  {
340 
341  for (; it!=end; ++it)
342  {
343  Elem* elem = *it;
344 
345  libmesh_assert (global_index_map.count(elem->id()));
346 
347  const dof_id_type elem_global_index =
348  global_index_map[elem->id()];
349 
350  libmesh_assert_less (elem_global_index, part.size());
351  const processor_id_type elem_procid =
352  static_cast<processor_id_type>(part[elem_global_index]);
353 
354  elem->processor_id() = elem_procid;
355  }
356  }
357 
358  STOP_LOG("partition()", "MetisPartitioner");
359 #endif
360 }
361 
362 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo