mesh_communication_global_indices.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/libmesh_common.h"
26 #include "libmesh/mesh_base.h"
27 #include "libmesh/mesh_tools.h"
29 #include "libmesh/parallel.h"
31 #include "libmesh/parallel_sort.h"
32 #include "libmesh/elem.h"
33 #include "libmesh/elem_range.h"
34 #include "libmesh/node_range.h"
35 #ifdef LIBMESH_HAVE_LIBHILBERT
36 # include "hilbert.h"
37 #endif
38 
39 #ifdef LIBMESH_HAVE_LIBHILBERT
40 namespace { // anonymous namespace for helper functions
41 
42  using namespace libMesh;
43 
44  // Utility function to map (x,y,z) in [bbox.min, bbox.max]^3 into
45  // [0,max_inttype]^3 for computing Hilbert keys
46  void get_hilbert_coords (const Point &p,
47  const MeshTools::BoundingBox &bbox,
48  CFixBitVec icoords[3])
49  {
50  static const Hilbert::inttype max_inttype = static_cast<Hilbert::inttype>(-1);
51 
52  const long double // put (x,y,z) in [0,1]^3 (don't divide by 0)
53  x = ((bbox.first(0) == bbox.second(0)) ? 0. :
54  (p(0)-bbox.first(0))/(bbox.second(0)-bbox.first(0))),
55 
56 #if LIBMESH_DIM > 1
57  y = ((bbox.first(1) == bbox.second(1)) ? 0. :
58  (p(1)-bbox.first(1))/(bbox.second(1)-bbox.first(1))),
59 #else
60  y = 0.,
61 #endif
62 
63 #if LIBMESH_DIM > 2
64  z = ((bbox.first(2) == bbox.second(2)) ? 0. :
65  (p(2)-bbox.first(2))/(bbox.second(2)-bbox.first(2)));
66 #else
67  z = 0.;
68 #endif
69 
70  // (icoords) in [0,max_inttype]^3
71  icoords[0] = static_cast<Hilbert::inttype>(x*max_inttype);
72  icoords[1] = static_cast<Hilbert::inttype>(y*max_inttype);
73  icoords[2] = static_cast<Hilbert::inttype>(z*max_inttype);
74  }
75 
76 
77 
78  // Compute the hilbert index
79  template <typename T>
80  Hilbert::HilbertIndices
81  get_hilbert_index (const T *p,
82  const MeshTools::BoundingBox &bbox)
83  {
84  static const unsigned int sizeof_inttype = sizeof(Hilbert::inttype);
85 
86  Hilbert::HilbertIndices index;
87  CFixBitVec icoords[3];
88  Hilbert::BitVecType bv;
89  get_hilbert_coords (*p, bbox, icoords);
90  Hilbert::coordsToIndex (icoords, 8*sizeof_inttype, 3, bv);
91  index = bv;
92 
93  return index;
94  }
95 
96  template <>
97  Hilbert::HilbertIndices
98  get_hilbert_index (const Elem *e,
99  const MeshTools::BoundingBox &bbox)
100  {
101  static const unsigned int sizeof_inttype = sizeof(Hilbert::inttype);
102 
103  Hilbert::HilbertIndices index;
104  CFixBitVec icoords[3];
105  Hilbert::BitVecType bv;
106  get_hilbert_coords (e->centroid(), bbox, icoords);
107  Hilbert::coordsToIndex (icoords, 8*sizeof_inttype, 3, bv);
108  index = bv;
109 
110  return index;
111  }
112 
113 
114 
115  // Compute the hilbert index
116  Hilbert::HilbertIndices
117  get_hilbert_index (const Point &p,
118  const MeshTools::BoundingBox &bbox)
119  {
120  static const unsigned int sizeof_inttype = sizeof(Hilbert::inttype);
121 
122  Hilbert::HilbertIndices index;
123  CFixBitVec icoords[3];
124  Hilbert::BitVecType bv;
125  get_hilbert_coords (p, bbox, icoords);
126  Hilbert::coordsToIndex (icoords, 8*sizeof_inttype, 3, bv);
127  index = bv;
128 
129  return index;
130  }
131 
132  // Helper class for threaded Hilbert key computation
133  class ComputeHilbertKeys
134  {
135  public:
136  ComputeHilbertKeys (const MeshTools::BoundingBox &bbox,
137  std::vector<Hilbert::HilbertIndices> &keys) :
138  _bbox(bbox),
139  _keys(keys)
140  {}
141 
142  // computes the hilbert index for a node
143  void operator() (const ConstNodeRange &range) const
144  {
145  dof_id_type pos = range.first_idx();
146  for (ConstNodeRange::const_iterator it = range.begin(); it!=range.end(); ++it)
147  {
148  const Node* node = (*it);
149  libmesh_assert(node);
150  libmesh_assert_less (pos, _keys.size());
151  _keys[pos++] = get_hilbert_index (*node, _bbox);
152  }
153  }
154 
155  // computes the hilbert index for an element
156  void operator() (const ConstElemRange &range) const
157  {
158  dof_id_type pos = range.first_idx();
159  for (ConstElemRange::const_iterator it = range.begin(); it!=range.end(); ++it)
160  {
161  const Elem* elem = (*it);
162  libmesh_assert(elem);
163  libmesh_assert_less (pos, _keys.size());
164  _keys[pos++] = get_hilbert_index (elem->centroid(), _bbox);
165  }
166  }
167 
168  private:
169  const MeshTools::BoundingBox &_bbox;
170  std::vector<Hilbert::HilbertIndices> &_keys;
171  };
172 }
173 #endif
174 
175 
176 namespace libMesh
177 {
178 
179 
180 // ------------------------------------------------------------
181 // MeshCommunication class members
182 #if defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
184 {
185  START_LOG ("assign_global_indices()", "MeshCommunication");
186 
187  // This method determines partition-agnostic global indices
188  // for nodes and elements.
189 
190  // Algorithm:
191  // (1) compute the Hilbert key for each local node/element
192  // (2) perform a parallel sort of the Hilbert key
193  // (3) get the min/max value on each processor
194  // (4) determine the position in the global ranking for
195  // each local object
196 
197  const Parallel::Communicator &communicator (mesh.comm());
198 
199  // Global bounding box
202 
203  //-------------------------------------------------------------
204  // (1) compute Hilbert keys
205  std::vector<Hilbert::HilbertIndices>
206  node_keys, elem_keys;
207 
208  {
209  // Nodes first
210  {
212  mesh.local_nodes_end());
213  node_keys.resize (nr.size());
214  Threads::parallel_for (nr, ComputeHilbertKeys (bbox, node_keys));
215 
216 #if 0
217  // It's O(N^2) to check that these keys don't duplicate before the
218  // sort...
220  for (std::size_t i = 0; i != node_keys.size(); ++i, ++nodei)
221  {
223  for (std::size_t j = 0; j != i; ++j, ++nodej)
224  {
225  if (node_keys[i] == node_keys[j])
226  {
227  libMesh::err << "Error: nodes with duplicate Hilbert keys!" <<
228  std::endl;
229  CFixBitVec icoords[3], jcoords[3];
230  get_hilbert_coords(**nodej, bbox, jcoords);
231  libMesh::err <<
232  "node " << (*nodej)->id() << ", " <<
233  *(Point*)(*nodej) << " has HilbertIndices " <<
234  node_keys[j] << std::endl;
235  get_hilbert_coords(**nodei, bbox, icoords);
236  libMesh::err <<
237  "node " << (*nodei)->id() << ", " <<
238  *(Point*)(*nodei) << " has HilbertIndices " <<
239  node_keys[i] << std::endl;
240  libmesh_error();
241  }
242  }
243  }
244 #endif
245  }
246 
247  // Elements next
248  {
250  mesh.local_elements_end());
251  elem_keys.resize (er.size());
252  Threads::parallel_for (er, ComputeHilbertKeys (bbox, elem_keys));
253 
254 #if 0
255  // For elements, the keys can be (and in the case of TRI, are
256  // expected to be) duplicates, but only if the elements are at
257  // different levels
259  for (std::size_t i = 0; i != elem_keys.size(); ++i, ++elemi)
260  {
262  for (std::size_t j = 0; j != i; ++j, ++elemj)
263  {
264  if ((elem_keys[i] == elem_keys[j]) &&
265  ((*elemi)->level() == (*elemj)->level()))
266  {
267  libMesh::err <<
268  "Error: level " << (*elemi)->level() <<
269  " elements with duplicate Hilbert keys!" <<
270  std::endl;
271  libMesh::err <<
272  "level " << (*elemj)->level() << " elem\n" <<
273  (**elemj) << " centroid " <<
274  (*elemj)->centroid() << " has HilbertIndices " <<
275  elem_keys[j] << " or " <<
276  get_hilbert_index((*elemj)->centroid(), bbox) <<
277  std::endl;
278  libMesh::err <<
279  "level " << (*elemi)->level() << " elem\n" <<
280  (**elemi) << " centroid " <<
281  (*elemi)->centroid() << " has HilbertIndices " <<
282  elem_keys[i] << " or " <<
283  get_hilbert_index((*elemi)->centroid(), bbox) <<
284  std::endl;
285  libmesh_error();
286  }
287  }
288  }
289 #endif
290  }
291  } // done computing Hilbert keys
292 
293 
294 
295  //-------------------------------------------------------------
296  // (2) parallel sort the Hilbert keys
298  node_keys);
299  node_sorter.sort(); /* done with node_keys */ //node_keys.clear();
300 
301  const std::vector<Hilbert::HilbertIndices> &my_node_bin =
302  node_sorter.bin();
303 
305  elem_keys);
306  elem_sorter.sort(); /* done with elem_keys */ //elem_keys.clear();
307 
308  const std::vector<Hilbert::HilbertIndices> &my_elem_bin =
309  elem_sorter.bin();
310 
311 
312 
313  //-------------------------------------------------------------
314  // (3) get the max value on each processor
315  std::vector<Hilbert::HilbertIndices>
316  node_upper_bounds(communicator.size()),
317  elem_upper_bounds(communicator.size());
318 
319  { // limit scope of temporaries
320  std::vector<Hilbert::HilbertIndices> recvbuf(2*communicator.size());
321  std::vector<unsigned short int> /* do not use a vector of bools here since it is not always so! */
322  empty_nodes (communicator.size()),
323  empty_elem (communicator.size());
324  std::vector<Hilbert::HilbertIndices> my_max(2);
325 
326  communicator.allgather (static_cast<unsigned short int>(my_node_bin.empty()), empty_nodes);
327  communicator.allgather (static_cast<unsigned short int>(my_elem_bin.empty()), empty_elem);
328 
329  if (!my_node_bin.empty()) my_max[0] = my_node_bin.back();
330  if (!my_elem_bin.empty()) my_max[1] = my_elem_bin.back();
331 
332  communicator.allgather (my_max, /* identical_buffer_sizes = */ true);
333 
334  // Be cereful here. The *_upper_bounds will be used to find the processor
335  // a given object belongs to. So, if a processor contains no objects (possible!)
336  // then copy the bound from the lower processor id.
337  for (processor_id_type p=0; p<communicator.size(); p++)
338  {
339  node_upper_bounds[p] = my_max[2*p+0];
340  elem_upper_bounds[p] = my_max[2*p+1];
341 
342  if (p > 0) // default hilbert index value is the OK upper bound for processor 0.
343  {
344  if (empty_nodes[p]) node_upper_bounds[p] = node_upper_bounds[p-1];
345  if (empty_elem[p]) elem_upper_bounds[p] = elem_upper_bounds[p-1];
346  }
347  }
348  }
349 
350 
351 
352  //-------------------------------------------------------------
353  // (4) determine the position in the global ranking for
354  // each local object
355  {
356  //----------------------------------------------
357  // Nodes first -- all nodes, not just local ones
358  {
359  // Request sets to send to each processor
360  std::vector<std::vector<Hilbert::HilbertIndices> >
361  requested_ids (communicator.size());
362  // Results to gather from each processor
363  std::vector<std::vector<dof_id_type> >
364  filled_request (communicator.size());
365 
366  {
369 
370  // build up list of requests
371  for (; it != end; ++it)
372  {
373  const Node* node = (*it);
374  libmesh_assert(node);
375  const Hilbert::HilbertIndices hi =
376  get_hilbert_index (*node, bbox);
377  const processor_id_type pid =
378  libmesh_cast_int<processor_id_type>
379  (std::distance (node_upper_bounds.begin(),
380  std::lower_bound(node_upper_bounds.begin(),
381  node_upper_bounds.end(),
382  hi)));
383 
384  libmesh_assert_less (pid, communicator.size());
385 
386  requested_ids[pid].push_back(hi);
387  }
388  }
389 
390  // The number of objects in my_node_bin on each processor
391  std::vector<dof_id_type> node_bin_sizes(communicator.size());
392  communicator.allgather (static_cast<dof_id_type>(my_node_bin.size()), node_bin_sizes);
393 
394  // The offset of my first global index
395  dof_id_type my_offset = 0;
396  for (processor_id_type pid=0; pid<communicator.rank(); pid++)
397  my_offset += node_bin_sizes[pid];
398 
399  // start with pid=0, so that we will trade with ourself
400  for (processor_id_type pid=0; pid<communicator.size(); pid++)
401  {
402  // Trade my requests with processor procup and procdown
403  const processor_id_type procup = (communicator.rank() + pid) %
404  communicator.size();
405  const processor_id_type procdown = (communicator.size() +
406  communicator.rank() - pid) %
407  communicator.size();
408 
409  std::vector<Hilbert::HilbertIndices> request_to_fill;
410  communicator.send_receive(procup, requested_ids[procup],
411  procdown, request_to_fill);
412 
413  // Fill the requests
414  std::vector<dof_id_type> global_ids; global_ids.reserve(request_to_fill.size());
415  for (std::size_t idx=0; idx<request_to_fill.size(); idx++)
416  {
417  const Hilbert::HilbertIndices &hi = request_to_fill[idx];
418  libmesh_assert_less_equal (hi, node_upper_bounds[communicator.rank()]);
419 
420  // find the requested index in my node bin
421  std::vector<Hilbert::HilbertIndices>::const_iterator pos =
422  std::lower_bound (my_node_bin.begin(), my_node_bin.end(), hi);
423  libmesh_assert (pos != my_node_bin.end());
424  libmesh_assert_equal_to (*pos, hi);
425 
426  // Finally, assign the global index based off the position of the index
427  // in my array, properly offset.
428  global_ids.push_back (std::distance(my_node_bin.begin(), pos) + my_offset);
429  }
430 
431  // and trade back
432  communicator.send_receive (procdown, global_ids,
433  procup, filled_request[procup]);
434  }
435 
436  // We now have all the filled requests, so we can loop through our
437  // nodes once and assign the global index to each one.
438  {
439  std::vector<std::vector<dof_id_type>::const_iterator>
440  next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
441  for (processor_id_type pid=0; pid<communicator.size(); pid++)
442  next_obj_on_proc.push_back(filled_request[pid].begin());
443 
444  {
446  const MeshBase::node_iterator end = mesh.nodes_end();
447 
448  for (; it != end; ++it)
449  {
450  Node* node = (*it);
451  libmesh_assert(node);
452  const Hilbert::HilbertIndices hi =
453  get_hilbert_index (*node, bbox);
454  const processor_id_type pid =
455  libmesh_cast_int<processor_id_type>
456  (std::distance (node_upper_bounds.begin(),
457  std::lower_bound(node_upper_bounds.begin(),
458  node_upper_bounds.end(),
459  hi)));
460 
461  libmesh_assert_less (pid, communicator.size());
462  libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());
463 
464  const dof_id_type global_index = *next_obj_on_proc[pid];
465  libmesh_assert_less (global_index, mesh.n_nodes());
466  node->set_id() = global_index;
467 
468  ++next_obj_on_proc[pid];
469  }
470  }
471  }
472  }
473 
474  //---------------------------------------------------
475  // elements next -- all elements, not just local ones
476  {
477  // Request sets to send to each processor
478  std::vector<std::vector<Hilbert::HilbertIndices> >
479  requested_ids (communicator.size());
480  // Results to gather from each processor
481  std::vector<std::vector<dof_id_type> >
482  filled_request (communicator.size());
483 
484  {
487 
488  for (; it != end; ++it)
489  {
490  const Elem* elem = (*it);
491  libmesh_assert(elem);
492  const Hilbert::HilbertIndices hi =
493  get_hilbert_index (elem->centroid(), bbox);
494  const processor_id_type pid =
495  libmesh_cast_int<processor_id_type>
496  (std::distance (elem_upper_bounds.begin(),
497  std::lower_bound(elem_upper_bounds.begin(),
498  elem_upper_bounds.end(),
499  hi)));
500 
501  libmesh_assert_less (pid, communicator.size());
502 
503  requested_ids[pid].push_back(hi);
504  }
505  }
506 
507  // The number of objects in my_elem_bin on each processor
508  std::vector<dof_id_type> elem_bin_sizes(communicator.size());
509  communicator.allgather (static_cast<dof_id_type>(my_elem_bin.size()), elem_bin_sizes);
510 
511  // The offset of my first global index
512  dof_id_type my_offset = 0;
513  for (processor_id_type pid=0; pid<communicator.rank(); pid++)
514  my_offset += elem_bin_sizes[pid];
515 
516  // start with pid=0, so that we will trade with ourself
517  for (processor_id_type pid=0; pid<communicator.size(); pid++)
518  {
519  // Trade my requests with processor procup and procdown
520  const processor_id_type procup = (communicator.rank() + pid) %
521  communicator.size();
522  const processor_id_type procdown = (communicator.size() +
523  communicator.rank() - pid) %
524  communicator.size();
525 
526  std::vector<Hilbert::HilbertIndices> request_to_fill;
527  communicator.send_receive(procup, requested_ids[procup],
528  procdown, request_to_fill);
529 
530  // Fill the requests
531  std::vector<dof_id_type> global_ids; global_ids.reserve(request_to_fill.size());
532  for (std::size_t idx=0; idx<request_to_fill.size(); idx++)
533  {
534  const Hilbert::HilbertIndices &hi = request_to_fill[idx];
535  libmesh_assert_less_equal (hi, elem_upper_bounds[communicator.rank()]);
536 
537  // find the requested index in my elem bin
538  std::vector<Hilbert::HilbertIndices>::const_iterator pos =
539  std::lower_bound (my_elem_bin.begin(), my_elem_bin.end(), hi);
540  libmesh_assert (pos != my_elem_bin.end());
541  libmesh_assert_equal_to (*pos, hi);
542 
543  // Finally, assign the global index based off the position of the index
544  // in my array, properly offset.
545  global_ids.push_back (std::distance(my_elem_bin.begin(), pos) + my_offset);
546  }
547 
548  // and trade back
549  communicator.send_receive (procdown, global_ids,
550  procup, filled_request[procup]);
551  }
552 
553  // We now have all the filled requests, so we can loop through our
554  // elements once and assign the global index to each one.
555  {
556  std::vector<std::vector<dof_id_type>::const_iterator>
557  next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
558  for (processor_id_type pid=0; pid<communicator.size(); pid++)
559  next_obj_on_proc.push_back(filled_request[pid].begin());
560 
561  {
564 
565  for (; it != end; ++it)
566  {
567  Elem* elem = (*it);
568  libmesh_assert(elem);
569  const Hilbert::HilbertIndices hi =
570  get_hilbert_index (elem->centroid(), bbox);
571  const processor_id_type pid =
572  libmesh_cast_int<processor_id_type>
573  (std::distance (elem_upper_bounds.begin(),
574  std::lower_bound(elem_upper_bounds.begin(),
575  elem_upper_bounds.end(),
576  hi)));
577 
578  libmesh_assert_less (pid, communicator.size());
579  libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());
580 
581  const dof_id_type global_index = *next_obj_on_proc[pid];
582  libmesh_assert_less (global_index, mesh.n_elem());
583  elem->set_id() = global_index;
584 
585  ++next_obj_on_proc[pid];
586  }
587  }
588  }
589  }
590  }
591 
592  STOP_LOG ("assign_global_indices()", "MeshCommunication");
593 }
594 #else // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
596 {
597 }
598 #endif // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
599 
600 
601 
602 #if defined(LIBMESH_HAVE_LIBHILBERT) && defined(LIBMESH_HAVE_MPI)
603 template <typename ForwardIterator>
605  const MeshTools::BoundingBox &bbox,
606  const ForwardIterator &begin,
607  const ForwardIterator &end,
608  std::vector<dof_id_type> &index_map) const
609 {
610  START_LOG ("find_global_indices()", "MeshCommunication");
611 
612  // This method determines partition-agnostic global indices
613  // for nodes and elements.
614 
615  // Algorithm:
616  // (1) compute the Hilbert key for each local node/element
617  // (2) perform a parallel sort of the Hilbert key
618  // (3) get the min/max value on each processor
619  // (4) determine the position in the global ranking for
620  // each local object
621  index_map.clear();
622  index_map.reserve(std::distance (begin, end));
623 
624  //-------------------------------------------------------------
625  // (1) compute Hilbert keys
626  // These aren't trivial to compute, and we will need them again.
627  // But the binsort will sort the input vector, trashing the order
628  // that we'd like to rely on. So, two vectors...
629  std::vector<Hilbert::HilbertIndices>
630  sorted_hilbert_keys,
631  hilbert_keys;
632  sorted_hilbert_keys.reserve(index_map.capacity());
633  hilbert_keys.reserve(index_map.capacity());
634  {
635  START_LOG("compute_hilbert_indices()", "MeshCommunication");
636  for (ForwardIterator it=begin; it!=end; ++it)
637  {
638  const Hilbert::HilbertIndices hi(get_hilbert_index (*it, bbox));
639  hilbert_keys.push_back(hi);
640 
641  if ((*it)->processor_id() == communicator.rank())
642  sorted_hilbert_keys.push_back(hi);
643 
644  // someone needs to take care of unpartitioned objects!
645  if ((communicator.rank() == 0) &&
646  ((*it)->processor_id() == DofObject::invalid_processor_id))
647  sorted_hilbert_keys.push_back(hi);
648  }
649  STOP_LOG("compute_hilbert_indices()", "MeshCommunication");
650  }
651 
652  //-------------------------------------------------------------
653  // (2) parallel sort the Hilbert keys
654  START_LOG ("parallel_sort()", "MeshCommunication");
655  Parallel::Sort<Hilbert::HilbertIndices> sorter (communicator,
656  sorted_hilbert_keys);
657  sorter.sort();
658  STOP_LOG ("parallel_sort()", "MeshCommunication");
659  const std::vector<Hilbert::HilbertIndices> &my_bin = sorter.bin();
660 
661  // The number of objects in my_bin on each processor
662  std::vector<unsigned int> bin_sizes(communicator.size());
663  communicator.allgather (static_cast<unsigned int>(my_bin.size()), bin_sizes);
664 
665  // The offset of my first global index
666  unsigned int my_offset = 0;
667  for (unsigned int pid=0; pid<communicator.rank(); pid++)
668  my_offset += bin_sizes[pid];
669 
670  //-------------------------------------------------------------
671  // (3) get the max value on each processor
672  std::vector<Hilbert::HilbertIndices>
673  upper_bounds(1);
674 
675  if (!my_bin.empty())
676  upper_bounds[0] = my_bin.back();
677 
678  communicator.allgather (upper_bounds, /* identical_buffer_sizes = */ true);
679 
680  // Be cereful here. The *_upper_bounds will be used to find the processor
681  // a given object belongs to. So, if a processor contains no objects (possible!)
682  // then copy the bound from the lower processor id.
683  for (unsigned int p=1; p<communicator.size(); p++)
684  if (!bin_sizes[p]) upper_bounds[p] = upper_bounds[p-1];
685 
686 
687  //-------------------------------------------------------------
688  // (4) determine the position in the global ranking for
689  // each local object
690  {
691  //----------------------------------------------
692  // all objects, not just local ones
693 
694  // Request sets to send to each processor
695  std::vector<std::vector<Hilbert::HilbertIndices> >
696  requested_ids (communicator.size());
697  // Results to gather from each processor
698  std::vector<std::vector<unsigned int> >
699  filled_request (communicator.size());
700 
701  // build up list of requests
702  std::vector<Hilbert::HilbertIndices>::const_iterator hi =
703  hilbert_keys.begin();
704 
705  for (ForwardIterator it = begin; it != end; ++it)
706  {
707  libmesh_assert (hi != hilbert_keys.end());
708  const processor_id_type pid =
709  libmesh_cast_int<processor_id_type>
710  (std::distance (upper_bounds.begin(),
711  std::lower_bound(upper_bounds.begin(),
712  upper_bounds.end(),
713  *hi)));
714 
715  libmesh_assert_less (pid, communicator.size());
716 
717  requested_ids[pid].push_back(*hi);
718 
719  hi++;
720  // go ahead and put pid in index_map, that way we
721  // don't have to repeat the std::lower_bound()
722  index_map.push_back(pid);
723  }
724 
725  // start with pid=0, so that we will trade with ourself
726  std::vector<Hilbert::HilbertIndices> request_to_fill;
727  std::vector<unsigned int> global_ids;
728  for (unsigned int pid=0; pid<communicator.size(); pid++)
729  {
730  // Trade my requests with processor procup and procdown
731  const unsigned int procup = (communicator.rank() + pid) %
732  communicator.size();
733  const unsigned int procdown = (communicator.size() +
734  communicator.rank() - pid) %
735  communicator.size();
736 
737  communicator.send_receive(procup, requested_ids[procup],
738  procdown, request_to_fill);
739 
740  // Fill the requests
741  global_ids.clear(); global_ids.reserve(request_to_fill.size());
742  for (unsigned int idx=0; idx<request_to_fill.size(); idx++)
743  {
744  const Hilbert::HilbertIndices &hilbert_indices = request_to_fill[idx];
745  libmesh_assert_less_equal (hilbert_indices, upper_bounds[communicator.rank()]);
746 
747  // find the requested index in my node bin
748  std::vector<Hilbert::HilbertIndices>::const_iterator pos =
749  std::lower_bound (my_bin.begin(), my_bin.end(), hilbert_indices);
750  libmesh_assert (pos != my_bin.end());
751  libmesh_assert_equal_to (*pos, hilbert_indices);
752 
753  // Finally, assign the global index based off the position of the index
754  // in my array, properly offset.
755  global_ids.push_back (std::distance(my_bin.begin(), pos) + my_offset);
756  }
757 
758  // and trade back
759  communicator.send_receive (procdown, global_ids,
760  procup, filled_request[procup]);
761  }
762 
763  // We now have all the filled requests, so we can loop through our
764  // nodes once and assign the global index to each one.
765  {
766  std::vector<std::vector<unsigned int>::const_iterator>
767  next_obj_on_proc; next_obj_on_proc.reserve(communicator.size());
768  for (unsigned int pid=0; pid<communicator.size(); pid++)
769  next_obj_on_proc.push_back(filled_request[pid].begin());
770 
771  unsigned int cnt=0;
772  for (ForwardIterator it = begin; it != end; ++it, cnt++)
773  {
774  const unsigned int pid = index_map[cnt];
775 
776  libmesh_assert_less (pid, communicator.size());
777  libmesh_assert (next_obj_on_proc[pid] != filled_request[pid].end());
778 
779  const unsigned int global_index = *next_obj_on_proc[pid];
780  index_map[cnt] = global_index;
781 
782  ++next_obj_on_proc[pid];
783  }
784  }
785  }
786 
787  STOP_LOG ("find_global_indices()", "MeshCommunication");
788 }
789 #else // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
790 template <typename ForwardIterator>
792  const MeshTools::BoundingBox &,
793  const ForwardIterator &begin,
794  const ForwardIterator &end,
795  std::vector<dof_id_type> &index_map) const
796 {
797  index_map.clear();
798  index_map.reserve(std::distance (begin, end));
799  unsigned int index = 0;
800  for (ForwardIterator it=begin; it!=end; ++it)
801  index_map.push_back(index++);
802 }
803 #endif // LIBMESH_HAVE_LIBHILBERT, LIBMESH_HAVE_MPI
804 
805 
806 
807 //------------------------------------------------------------------
808 template void MeshCommunication::find_global_indices<MeshBase::const_node_iterator> (const Parallel::Communicator &,
809  const MeshTools::BoundingBox &,
812  std::vector<dof_id_type> &) const;
813 
814 template void MeshCommunication::find_global_indices<MeshBase::const_element_iterator> (const Parallel::Communicator &,
815  const MeshTools::BoundingBox &,
818  std::vector<dof_id_type> &) const;
819 template void MeshCommunication::find_global_indices<MeshBase::node_iterator> (const Parallel::Communicator &,
820  const MeshTools::BoundingBox &,
821  const MeshBase::node_iterator &,
822  const MeshBase::node_iterator &,
823  std::vector<dof_id_type> &) const;
824 
825 template void MeshCommunication::find_global_indices<MeshBase::element_iterator> (const Parallel::Communicator &,
826  const MeshTools::BoundingBox &,
829  std::vector<dof_id_type> &) const;
830 
831 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo