Chaste  Release::2018.1
DistributedBoxCollection.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 #include "DistributedBoxCollection.hpp"
36 #include "Exception.hpp"
37 #include "MathsCustomFunctions.hpp"
38 #include "Warnings.hpp"
39 
40 // Static member for "fudge factor" is instantiated here
41 template<unsigned DIM>
42 const double DistributedBoxCollection<DIM>::msFudge = 5e-14;
43 
44 template<unsigned DIM>
45 DistributedBoxCollection<DIM>::DistributedBoxCollection(double boxWidth, c_vector<double, 2*DIM> domainSize, bool isPeriodicInX, int localRows)
46  : mBoxWidth(boxWidth),
47  mIsPeriodicInX(isPeriodicInX),
48  mAreLocalBoxesSet(false),
49  mCalculateNodeNeighbours(true)
50 {
51  // Periodicity only works in 2d
52  if (isPeriodicInX)
53  {
54  assert(DIM==2); // LCOV_EXCL_LINE
55  }
56 
57  // If the domain size is not 'divisible' (i.e. fmod(width, box_size) > 0.0) we swell the domain to enforce this.
58  for (unsigned i=0; i<DIM; i++)
59  {
60  double r = fmod((domainSize[2*i+1]-domainSize[2*i]), boxWidth);
61  if (r > 0.0)
62  {
63  domainSize[2*i+1] += boxWidth - r;
64  }
65  }
66 
67  mDomainSize = domainSize;
68 
69  // Calculate the number of boxes in each direction.
70  mNumBoxesEachDirection = scalar_vector<unsigned>(DIM, 0u);
71 
72  for (unsigned i=0; i<DIM; i++)
73  {
74  double counter = mDomainSize(2*i);
75  while (counter + msFudge < mDomainSize(2*i+1))
76  {
78  counter += mBoxWidth;
79  }
80  }
81 
82  // Make sure there are enough boxes for the number of processes.
84  {
85  WARNING("There are more processes than convenient for the domain/mesh/box size. The domain size has been swollen.")
88  }
89 
90  // Make a distributed vector factory to split the rows of boxes between processes.
92 
93  // Calculate how many boxes in a row / face. A useful piece of data in the class.
94  mNumBoxes = 1u;
95  for (unsigned dim=0; dim<DIM; dim++)
96  {
98  }
99 
101 
102  unsigned num_local_boxes = mNumBoxesInAFace * GetNumLocalRows();
103 
106 
107  // Create the correct number of boxes and set up halos
108  mBoxes.resize(num_local_boxes);
109  SetupHaloBoxes();
110 }
111 
112 template<unsigned DIM>
114 {
115  delete mpDistributedBoxStackFactory;
116 }
117 
118 template<unsigned DIM>
120 {
121  for (unsigned i=0; i<mBoxes.size(); i++)
122  {
123  mBoxes[i].ClearNodes();
124  }
125  for (unsigned i=0; i<mHaloBoxes.size(); i++)
126  {
127  mHaloBoxes[i].ClearNodes();
128  }
129 }
130 
131 template<unsigned DIM>
133 {
134  // Get top-most and bottom-most value of Distributed Box Stack.
135  unsigned hi = mpDistributedBoxStackFactory->GetHigh();
136  unsigned lo = mpDistributedBoxStackFactory->GetLow();
137 
138  // If I am not the top-most process, add halo structures above.
139  if (!PetscTools::AmTopMost())
140  {
141  for (unsigned i=0; i < mNumBoxesInAFace; i++)
142  {
143  Box<DIM> new_box;
144  mHaloBoxes.push_back(new_box);
145 
146  unsigned global_index = hi * mNumBoxesInAFace + i;
147  mHaloBoxesMapping[global_index] = mHaloBoxes.size()-1;
148  mHalosRight.push_back(global_index - mNumBoxesInAFace);
149  }
150  }
151 
152  // If I am not the bottom-most process, add halo structures below.
153  if (!PetscTools::AmMaster())
154  {
155  for (unsigned i=0; i< mNumBoxesInAFace; i++)
156  {
157  Box<DIM> new_box;
158  mHaloBoxes.push_back(new_box);
159 
160  unsigned global_index = (lo - 1) * mNumBoxesInAFace + i;
161  mHaloBoxesMapping[global_index] = mHaloBoxes.size() - 1;
162 
163  mHalosLeft.push_back(global_index + mNumBoxesInAFace);
164  }
165  }
166 }
167 
168 template<unsigned DIM>
170 {
171  mHaloNodesLeft.clear();
172  for (unsigned i=0; i<mHalosLeft.size(); i++)
173  {
174  for (typename std::set<Node<DIM>* >::iterator iter=this->rGetBox(mHalosLeft[i]).rGetNodesContained().begin();
175  iter!=this->rGetBox(mHalosLeft[i]).rGetNodesContained().end();
176  iter++)
177  {
178  mHaloNodesLeft.push_back((*iter)->GetIndex());
179  }
180  }
181 
182  // Send right
183  mHaloNodesRight.clear();
184  for (unsigned i=0; i<mHalosRight.size(); i++)
185  {
186  for (typename std::set<Node<DIM>* >::iterator iter=this->rGetBox(mHalosRight[i]).rGetNodesContained().begin();
187  iter!=this->rGetBox(mHalosRight[i]).rGetNodesContained().end();
188  iter++)
189  {
190  mHaloNodesRight.push_back((*iter)->GetIndex());
191  }
192  }
193 }
194 
195 template<unsigned DIM>
197 {
198  return mpDistributedBoxStackFactory->GetHigh() - mpDistributedBoxStackFactory->GetLow();
199 }
200 
201 template<unsigned DIM>
203 {
204  return (!(globalIndex<mMinBoxIndex) && !(mMaxBoxIndex<globalIndex));
205 }
206 
207 template<unsigned DIM>
209 {
210  bool is_halo_right = ((globalIndex > mMaxBoxIndex) && !(globalIndex > mMaxBoxIndex + mNumBoxesInAFace));
211  bool is_halo_left = ((globalIndex < mMinBoxIndex) && !(globalIndex < mMinBoxIndex - mNumBoxesInAFace));
212 
213  return (PetscTools::IsParallel() && (is_halo_right || is_halo_left));
214 }
215 
216 template<unsigned DIM>
218 {
219  bool is_on_boundary = !(globalIndex < mMaxBoxIndex - mNumBoxesInAFace) || (globalIndex < mMinBoxIndex + mNumBoxesInAFace);
220 
221  return (PetscTools::IsSequential() || !(is_on_boundary));
222 }
223 
224 template<unsigned DIM>
225 unsigned DistributedBoxCollection<DIM>::CalculateGlobalIndex(c_vector<unsigned, DIM> gridIndices)
226 {
228  unsigned global_index;
229 
230  switch (DIM)
231  {
232  case 1:
233  {
234  global_index = gridIndices(0);
235  break;
236  }
237  case 2:
238  {
239  global_index = gridIndices(0) +
240  gridIndices(1) * mNumBoxesEachDirection(0);
241  break;
242  }
243  case 3:
244  {
245  global_index = gridIndices(0) +
246  gridIndices(1) * mNumBoxesEachDirection(0) +
247  gridIndices(2) * mNumBoxesEachDirection(0) * mNumBoxesEachDirection(1);
248  break;
249  }
250  default:
251  {
253  }
254  }
255  return global_index;
256 }
257 
258 template<unsigned DIM>
260 {
261  // Get the location of the node
262  c_vector<double, DIM> location = pNode->rGetLocation();
263  return CalculateContainingBox(location);
264 }
265 
266 template<unsigned DIM>
267 unsigned DistributedBoxCollection<DIM>::CalculateContainingBox(c_vector<double, DIM>& rLocation)
268 {
269  // The node must lie inside the boundary of the box collection
270  for (unsigned i=0; i<DIM; i++)
271  {
272  if ((rLocation[i] < mDomainSize(2*i)) || !(rLocation[i] < mDomainSize(2*i+1)))
273  {
274  EXCEPTION("The point provided is outside all of the boxes");
275  }
276  }
277 
278  // Compute the containing box index in each dimension
279  c_vector<unsigned, DIM> containing_box_indices = scalar_vector<unsigned>(DIM, 0u);
280  for (unsigned i=0; i<DIM; i++)
281  {
282  double box_counter = mDomainSize(2*i);
283  while (!((box_counter + mBoxWidth) > rLocation[i] + msFudge))
284  {
285  containing_box_indices[i]++;
286  box_counter += mBoxWidth;
287  }
288  }
289 
290  // Use these to compute the index of the containing box
291  unsigned containing_box_index = 0;
292  for (unsigned i=0; i<DIM; i++)
293  {
294  unsigned temp = 1;
295  for (unsigned j=0; j<i; j++)
296  {
297  temp *= mNumBoxesEachDirection(j);
298  }
299  containing_box_index += temp*containing_box_indices[i];
300  }
301 
302  // This index must be less than the total number of boxes
303  assert(containing_box_index < mNumBoxes);
304 
305  return containing_box_index;
306 }
307 
308 template<unsigned DIM>
309 c_vector<unsigned, DIM> DistributedBoxCollection<DIM>::CalculateGridIndices(unsigned globalIndex)
310 {
311  c_vector<unsigned, DIM> grid_indices;
312 
313  switch (DIM)
314  {
315  case 1:
316  {
317  grid_indices(0) = globalIndex;
318  break;
319  }
320  case 2:
321  {
322  unsigned num_x = mNumBoxesEachDirection(0);
323  grid_indices(0) = globalIndex % num_x;
324  grid_indices(1) = (globalIndex - grid_indices(0)) / num_x;
325  break;
326  }
327  case 3:
328  {
329  unsigned num_x = mNumBoxesEachDirection(0);
330  unsigned num_xy = mNumBoxesEachDirection(0)*mNumBoxesEachDirection(1);
331  grid_indices(0) = globalIndex % num_x;
332  grid_indices(1) = (globalIndex % num_xy - grid_indices(0)) / num_x;
333  grid_indices(2) = globalIndex / num_xy;
334  break;
335  }
336  default:
337  {
339  }
340  }
341 
342  return grid_indices;
343 }
344 
345 template<unsigned DIM>
347 {
348  // Check first for local ownership
349  if (!(boxIndex<mMinBoxIndex) && !(mMaxBoxIndex<boxIndex))
350  {
351  return mBoxes[boxIndex-mMinBoxIndex];
352  }
353 
354  // If normal execution reaches this point then the box does not belong to the process so we will check for a halo box
355  return rGetHaloBox(boxIndex);
356 }
357 
358 template<unsigned DIM>
360 {
361  assert(IsHaloBox(boxIndex));
362 
363  unsigned local_index = mHaloBoxesMapping.find(boxIndex)->second;
364 
365  return mHaloBoxes[local_index];
366 }
367 
368 template<unsigned DIM>
370 {
371  return mNumBoxes;
372 }
373 
374 template<unsigned DIM>
376 {
377  return mBoxes.size();
378 }
379 
380 template<unsigned DIM>
381 c_vector<double, 2*DIM> DistributedBoxCollection<DIM>::rGetDomainSize() const
382 {
383  return mDomainSize;
384 }
385 
386 template<unsigned DIM>
388 {
389  return mAreLocalBoxesSet;
390 }
391 
392 template<unsigned DIM>
394 {
395  return mBoxWidth;
396 }
397 
398 template<unsigned DIM>
400 {
401  return mIsPeriodicInX;
402 }
403 
404 template<unsigned DIM>
406 {
407  return mpDistributedBoxStackFactory->GetHigh() - mpDistributedBoxStackFactory->GetLow();
408 }
409 
410 template<unsigned DIM>
411 int DistributedBoxCollection<DIM>::LoadBalance(std::vector<int> localDistribution)
412 {
413  MPI_Status status;
414 
415  int proc_right = (PetscTools::AmTopMost()) ? MPI_PROC_NULL : (int)PetscTools::GetMyRank() + 1;
416  int proc_left = (PetscTools::AmMaster()) ? MPI_PROC_NULL : (int)PetscTools::GetMyRank() - 1;
417 
418  // A variable that will return the new number of rows.
419  int new_rows = localDistribution.size();
420 
424  unsigned rows_on_left_process = 0;
425  std::vector<int> node_distr_on_left_process;
426 
427  unsigned num_local_rows = localDistribution.size();
428 
429  MPI_Send(&num_local_rows, 1, MPI_UNSIGNED, proc_right, 123, PETSC_COMM_WORLD);
430  MPI_Recv(&rows_on_left_process, 1, MPI_UNSIGNED, proc_left, 123, PETSC_COMM_WORLD, &status);
431 
432  node_distr_on_left_process.resize(rows_on_left_process > 0 ? rows_on_left_process : 1);
433 
434  MPI_Send(&localDistribution[0], num_local_rows, MPI_INT, proc_right, 123, PETSC_COMM_WORLD);
435  MPI_Recv(&node_distr_on_left_process[0], rows_on_left_process, MPI_INT, proc_left, 123, PETSC_COMM_WORLD, &status);
436 
440  int local_load = 0;
441  for (unsigned i=0; i<localDistribution.size(); i++)
442  {
443  local_load += localDistribution[i];
444  }
445  int load_on_left_proc = 0;
446  for (unsigned i=0; i<node_distr_on_left_process.size(); i++)
447  {
448  load_on_left_proc += node_distr_on_left_process[i];
449  }
450 
451  if (!PetscTools::AmMaster())
452  {
453  // Calculate (Difference in load with a shift) - (Difference in current loads) for a move left and right of the boundary
454  // This code uses integer arithmetic in order to avoid the rounding errors associated with doubles
455  int local_to_left_sq = (local_load - load_on_left_proc) * (local_load - load_on_left_proc);
456  int delta_left = ( (local_load + node_distr_on_left_process[node_distr_on_left_process.size() - 1]) - (load_on_left_proc - node_distr_on_left_process[node_distr_on_left_process.size() - 1]) );
457  delta_left = delta_left*delta_left - local_to_left_sq;
458 
459  int delta_right = ( (local_load - localDistribution[0]) - (load_on_left_proc + localDistribution[0]));
460  delta_right = delta_right*delta_right - local_to_left_sq;
461 
462  // If a delta is negative we should accept that change. If both are negative choose the largest change.
463  int local_change = 0;
464  bool move_left = (!(delta_left > 0) && (node_distr_on_left_process.size() > 1));
465  if (move_left)
466  {
467  local_change = 1;
468  }
469 
470  bool move_right = !(delta_right > 0) && (localDistribution.size() > 2);
471  if (move_right)
472  {
473  local_change = -1;
474  }
475 
476  if (move_left && move_right)
477  {
478  local_change = (fabs((double)delta_right) > fabs((double)delta_left)) ? -1 : 1;
479  }
480 
481  // Update the number of local rows.
482  new_rows += local_change;
483 
484  // Send the result of the calculation back to the left processes.
485  MPI_Send(&local_change, 1, MPI_INT, proc_left, 123, PETSC_COMM_WORLD);
486  }
487 
488  // Receive changes from right hand process.
489  int remote_change = 0;
490  MPI_Recv(&remote_change, 1, MPI_INT, proc_right, 123, PETSC_COMM_WORLD, &status);
491 
492  // Update based on change or right/top boundary
493  new_rows -= remote_change;
494 
495  return new_rows;
496 }
497 
498 template<unsigned DIM>
500 {
501  if (mAreLocalBoxesSet)
502  {
503  EXCEPTION("Local Boxes Are Already Set");
504  }
505  else
506  {
507  switch (DIM)
508  {
509  case 1:
510  {
511  // We only need to look for neighbours in the current and successive boxes plus some others for halos
512  mLocalBoxes.clear();
513 
514  // Iterate over the global box indices
515  for (unsigned global_index = mMinBoxIndex; global_index<mMaxBoxIndex+1; global_index++)
516  {
517  std::set<unsigned> local_boxes;
518 
519  // Insert the current box
520  local_boxes.insert(global_index);
521 
522  // Set some bools to find out where we are
523  bool right = (global_index==mNumBoxesEachDirection(0)-1);
524  bool left = (global_index == 0);
525  bool proc_left = (global_index == mpDistributedBoxStackFactory->GetLow());
526 
527  // If we're not at the right-most box, then insert the box to the right
528  if (!right)
529  {
530  local_boxes.insert(global_index+1);
531  }
532  // If we're on a left process boundary and not on process 0, insert the (halo) box to the left
533  if (proc_left && !left)
534  {
535  local_boxes.insert(global_index-1);
536  }
537 
538  mLocalBoxes.push_back(local_boxes);
539  }
540  break;
541  }
542  case 2:
543  {
544  // We only need to look for neighbours in the current box and half the neighbouring boxes plus some others for halos
545  mLocalBoxes.clear();
546 
547  for (unsigned global_index = mMinBoxIndex; global_index<mMaxBoxIndex+1; global_index++)
548  {
549  std::set<unsigned> local_boxes;
550 
551  // Set up bools to find out where we are
552  bool left = (global_index%mNumBoxesEachDirection(0) == 0);
553  bool right = (global_index%mNumBoxesEachDirection(0) == mNumBoxesEachDirection(0)-1);
554  bool top = !(global_index < mNumBoxesEachDirection(0)*mNumBoxesEachDirection(1) - mNumBoxesEachDirection(0));
555  bool bottom = (global_index < mNumBoxesEachDirection(0));
556  bool bottom_proc = (CalculateGridIndices(global_index)[1] == mpDistributedBoxStackFactory->GetLow());
557 
558  // Insert the current box
559  local_boxes.insert(global_index);
560 
561  // If we're on the bottom of the process boundary, but not the bottom of the domain add boxes below
562  if (!bottom && bottom_proc)
563  {
564  local_boxes.insert(global_index - mNumBoxesEachDirection(0));
565  if (!left)
566  {
567  local_boxes.insert(global_index - mNumBoxesEachDirection(0) - 1);
568  }
569  if (!right)
570  {
571  local_boxes.insert(global_index - mNumBoxesEachDirection(0) + 1);
572  }
573  }
574 
575  // If we're not at the top of the domain insert boxes above
576  if (!top)
577  {
578  local_boxes.insert(global_index + mNumBoxesEachDirection(0));
579 
580  if (!right)
581  {
582  local_boxes.insert(global_index + mNumBoxesEachDirection(0) + 1);
583  }
584  if (!left)
585  {
586  local_boxes.insert(global_index + mNumBoxesEachDirection(0) - 1);
587  }
588  else if ((global_index % mNumBoxesEachDirection(0) == 0) && (mIsPeriodicInX)) // If we're on the left edge but its periodic include the box on the far right and up one.
589  {
590  local_boxes.insert(global_index + 2 * mNumBoxesEachDirection(0) - 1);
591  }
592  }
593 
594  // If we're not on the far right hand side inseryt box to the right
595  if (!right)
596  {
597  local_boxes.insert(global_index + 1);
598  }
599  // If we're on the right edge but it's periodic include the box on the far left of the domain
600  else if ((global_index % mNumBoxesEachDirection(0) == mNumBoxesEachDirection(0)-1) && (mIsPeriodicInX))
601  {
602  local_boxes.insert(global_index - mNumBoxesEachDirection(0) + 1);
603  // If we're also not on the top-most row, then insert the box above- on the far left of the domain
604  if (global_index < mBoxes.size() - mNumBoxesEachDirection(0))
605  {
606  local_boxes.insert(global_index + 1);
607  }
608  }
609 
610  mLocalBoxes.push_back(local_boxes);
611  }
612  break;
613  }
614  case 3:
615  {
616  // We only need to look for neighbours in the current box and half the neighbouring boxes plus some others for halos
617  mLocalBoxes.clear();
618  unsigned num_boxes_xy = mNumBoxesEachDirection(0)*mNumBoxesEachDirection(1);
619 
620 
621  for (unsigned global_index = mMinBoxIndex; global_index<mMaxBoxIndex+1; global_index++)
622  {
623  std::set<unsigned> local_boxes;
624 
625  // Set some bools to find out where we are
626  bool top = !(global_index % num_boxes_xy < num_boxes_xy - mNumBoxesEachDirection(0));
627  bool bottom = (global_index % num_boxes_xy < mNumBoxesEachDirection(0));
628  bool left = (global_index % mNumBoxesEachDirection(0) == 0);
629  bool right = (global_index % mNumBoxesEachDirection(0) == mNumBoxesEachDirection(0) - 1);
630  bool front = (global_index < num_boxes_xy);
631  bool back = !(global_index < num_boxes_xy*mNumBoxesEachDirection(2) - num_boxes_xy);
632  bool proc_front = (CalculateGridIndices(global_index)[2] == mpDistributedBoxStackFactory->GetLow());
633  bool proc_back = (CalculateGridIndices(global_index)[2] == mpDistributedBoxStackFactory->GetHigh()-1);
634 
635  // Insert the current box
636  local_boxes.insert(global_index);
637 
638  // If we're not on the front face, add appropriate boxes on the closer face
639  if (!front)
640  {
641  // If we're not on the top of the domain
642  if (!top)
643  {
644  local_boxes.insert( global_index - num_boxes_xy + mNumBoxesEachDirection(0) );
645  if (!left)
646  {
647  local_boxes.insert( global_index - num_boxes_xy + mNumBoxesEachDirection(0) - 1);
648  }
649  if (!right)
650  {
651  local_boxes.insert( global_index - num_boxes_xy + mNumBoxesEachDirection(0) + 1);
652  }
653  }
654  if (!right)
655  {
656  local_boxes.insert( global_index - num_boxes_xy + 1);
657  }
658 
659  // If we are on the front of the process we have to add extra boxes as they are halos.
660  if (proc_front)
661  {
662  local_boxes.insert( global_index - num_boxes_xy );
663 
664  if (!left)
665  {
666  local_boxes.insert( global_index - num_boxes_xy - 1);
667  }
668  if (!bottom)
669  {
670  local_boxes.insert( global_index - num_boxes_xy - mNumBoxesEachDirection(0));
671 
672  if (!left)
673  {
674  local_boxes.insert( global_index - num_boxes_xy - mNumBoxesEachDirection(0) - 1);
675  }
676  if (!right)
677  {
678  local_boxes.insert( global_index - num_boxes_xy - mNumBoxesEachDirection(0) + 1);
679  }
680  }
681 
682  }
683  }
684  if (!right)
685  {
686  local_boxes.insert( global_index + 1);
687  }
688  // If we're not on the very top add boxes above
689  if (!top)
690  {
691  local_boxes.insert( global_index + mNumBoxesEachDirection(0));
692 
693  if (!right)
694  {
695  local_boxes.insert( global_index + mNumBoxesEachDirection(0) + 1);
696  }
697  if (!left)
698  {
699  local_boxes.insert( global_index + mNumBoxesEachDirection(0) - 1);
700  }
701  }
702 
703  // If we're not on the back add boxes behind
704  if (!back)
705  {
706  local_boxes.insert(global_index + num_boxes_xy);
707 
708  if (!right)
709  {
710  local_boxes.insert(global_index + num_boxes_xy + 1);
711  }
712  if (!top)
713  {
714  local_boxes.insert(global_index + num_boxes_xy + mNumBoxesEachDirection(0));
715  if (!right)
716  {
717  local_boxes.insert(global_index + num_boxes_xy + mNumBoxesEachDirection(0) + 1);
718  }
719  if (!left)
720  {
721  local_boxes.insert(global_index + num_boxes_xy + mNumBoxesEachDirection(0) - 1);
722  }
723  }
724  // If we are on the back proc we should make sure we get everything in the face further back
725  if (proc_back)
726  {
727  if (!left)
728  {
729  local_boxes.insert(global_index + num_boxes_xy - 1);
730  }
731  if (!bottom)
732  {
733  local_boxes.insert(global_index + num_boxes_xy - mNumBoxesEachDirection(0));
734 
735  if (!left)
736  {
737  local_boxes.insert(global_index + num_boxes_xy - mNumBoxesEachDirection(0) - 1);
738  }
739  if (!right)
740  {
741  local_boxes.insert(global_index + num_boxes_xy - mNumBoxesEachDirection(0) + 1);
742  }
743  }
744  }
745  }
746  mLocalBoxes.push_back(local_boxes);
747  }
748 
749  break;
750  }
751  default:
753  }
754  mAreLocalBoxesSet=true;
755  }
756 }
757 
758 template<unsigned DIM>
760 {
761  mAreLocalBoxesSet = true;
762  switch (DIM)
763  {
764  case 1:
765  {
766  for (unsigned i=mMinBoxIndex; i<mMaxBoxIndex+1; i++)
767  {
768  std::set<unsigned> local_boxes;
769 
770  local_boxes.insert(i);
771 
772  // add the two neighbours
773  if (i!=0)
774  {
775  local_boxes.insert(i-1);
776  }
777  if (i+1 != mNumBoxesEachDirection(0))
778  {
779  local_boxes.insert(i+1);
780  }
781 
782  mLocalBoxes.push_back(local_boxes);
783  }
784  break;
785  }
786  case 2:
787  {
788  mLocalBoxes.clear();
789 
790  unsigned M = mNumBoxesEachDirection(0);
791  unsigned N = mNumBoxesEachDirection(1);
792 
793  std::vector<bool> is_xmin(N*M); // far left
794  std::vector<bool> is_xmax(N*M); // far right
795  std::vector<bool> is_ymin(N*M); // bottom
796  std::vector<bool> is_ymax(N*M); // top
797 
798  for (unsigned i=0; i<M*N; i++)
799  {
800  is_xmin[i] = (i%M==0);
801  is_xmax[i] = ((i+1)%M==0);
802  is_ymin[i] = (i%(M*N)<M);
803  is_ymax[i] = (i%(M*N)>=(N-1)*M);
804  }
805 
806  for (unsigned i=mMinBoxIndex; i<mMaxBoxIndex+1; i++)
807  {
808  std::set<unsigned> local_boxes;
809 
810  local_boxes.insert(i);
811 
812  // add the box to the left
813  if (!is_xmin[i])
814  {
815  local_boxes.insert(i-1);
816  }
817  else // Add Periodic Box if needed
818  {
819  if (mIsPeriodicInX)
820  {
821  local_boxes.insert(i+M-1);
822  }
823  }
824 
825  // add the box to the right
826  if (!is_xmax[i])
827  {
828  local_boxes.insert(i+1);
829  }
830  else // Add Periodic Box if needed
831  {
832  if (mIsPeriodicInX)
833  {
834  local_boxes.insert(i-M+1);
835  }
836  }
837 
838  // add the one below
839  if (!is_ymin[i])
840  {
841  local_boxes.insert(i-M);
842  }
843 
844  // add the one above
845  if (!is_ymax[i])
846  {
847  local_boxes.insert(i+M);
848  }
849 
850  // add the four corner boxes
851 
852  if ((!is_xmin[i]) && (!is_ymin[i]))
853  {
854  local_boxes.insert(i-1-M);
855  }
856  if ((!is_xmin[i]) && (!is_ymax[i]))
857  {
858  local_boxes.insert(i-1+M);
859  }
860  if ((!is_xmax[i]) && (!is_ymin[i]))
861  {
862  local_boxes.insert(i+1-M);
863  }
864  if ((!is_xmax[i]) && (!is_ymax[i]))
865  {
866  local_boxes.insert(i+1+M);
867  }
868 
869  // Add periodic corner boxes if needed
870  if (mIsPeriodicInX)
871  {
872  if ((is_xmin[i]) && (!is_ymin[i]))
873  {
874  local_boxes.insert(i-1);
875  }
876  if ((is_xmin[i]) && (!is_ymax[i]))
877  {
878  local_boxes.insert(i-1+2*M);
879  }
880  if ((is_xmax[i]) && (!is_ymin[i]))
881  {
882  local_boxes.insert(i+1-2*M);
883  }
884  if ((is_xmax[i]) && (!is_ymax[i]))
885  {
886  local_boxes.insert(i+1);
887  }
888  }
889 
890  mLocalBoxes.push_back(local_boxes);
891  }
892  break;
893  }
894  case 3:
895  {
896  mLocalBoxes.clear();
897 
898  unsigned M = mNumBoxesEachDirection(0);
899  unsigned N = mNumBoxesEachDirection(1);
900  unsigned P = mNumBoxesEachDirection(2);
901 
902  std::vector<bool> is_xmin(N*M*P); // far left
903  std::vector<bool> is_xmax(N*M*P); // far right
904  std::vector<bool> is_ymin(N*M*P); // nearest
905  std::vector<bool> is_ymax(N*M*P); // furthest
906  std::vector<bool> is_zmin(N*M*P); // bottom layer
907  std::vector<bool> is_zmax(N*M*P); // top layer
908 
909  for (unsigned i=0; i<M*N*P; i++)
910  {
911  is_xmin[i] = (i%M==0);
912  is_xmax[i] = ((i+1)%M==0);
913  is_ymin[i] = (i%(M*N)<M);
914  is_ymax[i] = (i%(M*N)>=(N-1)*M);
915  is_zmin[i] = (i<M*N);
916  is_zmax[i] = (i>=M*N*(P-1));
917  }
918 
919  for (unsigned i=mMinBoxIndex; i<mMaxBoxIndex+1; i++)
920  {
921  std::set<unsigned> local_boxes;
922 
923  // add itself as a local box
924  local_boxes.insert(i);
925 
926  // now add all 26 other neighbours.....
927 
928  // add the box left
929  if (!is_xmin[i])
930  {
931  local_boxes.insert(i-1);
932 
933  // plus some others towards the left
934  if (!is_ymin[i])
935  {
936  local_boxes.insert(i-1-M);
937  }
938 
939  if (!is_ymax[i])
940  {
941  local_boxes.insert(i-1+M);
942  }
943 
944  if (!is_zmin[i])
945  {
946  local_boxes.insert(i-1-M*N);
947  }
948 
949  if (!is_zmax[i])
950  {
951  local_boxes.insert(i-1+M*N);
952  }
953  }
954 
955  // add the box to the right
956  if (!is_xmax[i])
957  {
958  local_boxes.insert(i+1);
959 
960  // plus some others towards the right
961  if (!is_ymin[i])
962  {
963  local_boxes.insert(i+1-M);
964  }
965 
966  if (!is_ymax[i])
967  {
968  local_boxes.insert(i+1+M);
969  }
970 
971  if (!is_zmin[i])
972  {
973  local_boxes.insert(i+1-M*N);
974  }
975 
976  if (!is_zmax[i])
977  {
978  local_boxes.insert(i+1+M*N);
979  }
980  }
981 
982  // add the boxes next along the y axis
983  if (!is_ymin[i])
984  {
985  local_boxes.insert(i-M);
986 
987  // and more in this plane
988  if (!is_zmin[i])
989  {
990  local_boxes.insert(i-M-M*N);
991  }
992 
993  if (!is_zmax[i])
994  {
995  local_boxes.insert(i-M+M*N);
996  }
997  }
998 
999  // add the boxes next along the y axis
1000  if (!is_ymax[i])
1001  {
1002  local_boxes.insert(i+M);
1003 
1004  // and more in this plane
1005  if (!is_zmin[i])
1006  {
1007  local_boxes.insert(i+M-M*N);
1008  }
1009 
1010  if (!is_zmax[i])
1011  {
1012  local_boxes.insert(i+M+M*N);
1013  }
1014  }
1015 
1016  // add the box directly above
1017  if (!is_zmin[i])
1018  {
1019  local_boxes.insert(i-N*M);
1020  }
1021 
1022  // add the box directly below
1023  if (!is_zmax[i])
1024  {
1025  local_boxes.insert(i+N*M);
1026  }
1027 
1028  // finally, the 8 corners are left
1029 
1030  if ((!is_xmin[i]) && (!is_ymin[i]) && (!is_zmin[i]))
1031  {
1032  local_boxes.insert(i-1-M-M*N);
1033  }
1034 
1035  if ((!is_xmin[i]) && (!is_ymin[i]) && (!is_zmax[i]))
1036  {
1037  local_boxes.insert(i-1-M+M*N);
1038  }
1039 
1040  if ((!is_xmin[i]) && (!is_ymax[i]) && (!is_zmin[i]))
1041  {
1042  local_boxes.insert(i-1+M-M*N);
1043  }
1044 
1045  if ((!is_xmin[i]) && (!is_ymax[i]) && (!is_zmax[i]))
1046  {
1047  local_boxes.insert(i-1+M+M*N);
1048  }
1049 
1050  if ((!is_xmax[i]) && (!is_ymin[i]) && (!is_zmin[i]))
1051  {
1052  local_boxes.insert(i+1-M-M*N);
1053  }
1054 
1055  if ((!is_xmax[i]) && (!is_ymin[i]) && (!is_zmax[i]))
1056  {
1057  local_boxes.insert(i+1-M+M*N);
1058  }
1059 
1060  if ((!is_xmax[i]) && (!is_ymax[i]) && (!is_zmin[i]))
1061  {
1062  local_boxes.insert(i+1+M-M*N);
1063  }
1064 
1065  if ((!is_xmax[i]) && (!is_ymax[i]) && (!is_zmax[i]))
1066  {
1067  local_boxes.insert(i+1+M+M*N);
1068  }
1069 
1070  mLocalBoxes.push_back(local_boxes);
1071  }
1072  break;
1073  }
1074  default:
1075  NEVER_REACHED;
1076  }
1077 }
1078 
1079 template<unsigned DIM>
1080 std::set<unsigned>& DistributedBoxCollection<DIM>::rGetLocalBoxes(unsigned boxIndex)
1081 {
1082  // Make sure the box is locally owned
1083  assert(!(boxIndex < mMinBoxIndex) && !(mMaxBoxIndex<boxIndex));
1084  return mLocalBoxes[boxIndex-mMinBoxIndex];
1085 }
1086 
1087 template<unsigned DIM>
1089 {
1090  unsigned index = CalculateContainingBox(pNode);
1091 
1092  return IsBoxOwned(index);
1093 }
1094 
1095 template<unsigned DIM>
1096 bool DistributedBoxCollection<DIM>::IsOwned(c_vector<double, DIM>& location)
1097 {
1098  unsigned index = CalculateContainingBox(location);
1099 
1100  return IsBoxOwned(index);
1101 }
1102 
1103 template<unsigned DIM>
1105 {
1106  unsigned box_index = CalculateContainingBox(pNode);
1107  unsigned containing_process = PetscTools::GetMyRank();
1108 
1109  if (box_index > mMaxBoxIndex)
1110  {
1111  containing_process++;
1112  }
1113  else if (box_index < mMinBoxIndex)
1114  {
1115  containing_process--;
1116  }
1117 
1118  return containing_process;
1119 }
1120 
1121 template<unsigned DIM>
1123 {
1124  return mHaloNodesRight;
1125 }
1126 
1127 template<unsigned DIM>
1129 {
1130  return mHaloNodesLeft;
1131 }
1132 
1133 template<unsigned DIM>
1135 {
1136  mCalculateNodeNeighbours = calculateNodeNeighbours;
1137 }
1138 
1139 template<unsigned DIM>
1140 void DistributedBoxCollection<DIM>::CalculateNodePairs(std::vector<Node<DIM>*>& rNodes, std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
1141 {
1142  rNodePairs.clear();
1143 
1144  // Create an empty neighbours set for each node
1145  for (unsigned i=0; i<rNodes.size(); i++)
1146  {
1147  // Get the box containing this node as only nodes on this process have NodeAttributes
1148  // and therefore Neighbours setup.
1149  unsigned box_index = CalculateContainingBox(rNodes[i]);
1150 
1151  if (IsBoxOwned(box_index))
1152  {
1153  rNodes[i]->ClearNeighbours();
1154  }
1155  }
1156 
1157  for (unsigned box_index=mMinBoxIndex; box_index<=mMaxBoxIndex; box_index++)
1158  {
1159  AddPairsFromBox(box_index, rNodePairs);
1160  }
1161 
1162  if (mCalculateNodeNeighbours)
1163  {
1164  for (unsigned i = 0; i < rNodes.size(); i++)
1165  {
1166  // Get the box containing this node as only nodes on this process have NodeAttributes
1167  // and therefore Neighbours setup.
1168  unsigned box_index = CalculateContainingBox(rNodes[i]);
1169 
1170  if (IsBoxOwned(box_index))
1171  {
1172  rNodes[i]->RemoveDuplicateNeighbours();
1173  }
1174  }
1175  }
1176 }
1177 
1178 template<unsigned DIM>
1179 void DistributedBoxCollection<DIM>::CalculateInteriorNodePairs(std::vector<Node<DIM>*>& rNodes, std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
1180 {
1181  rNodePairs.clear();
1182 
1183  // Create an empty neighbours set for each node
1184  for (unsigned i=0; i<rNodes.size(); i++)
1185  {
1186  // Get the box containing this node as only nodes on this process have NodeAttributes
1187  // and therefore Neighbours setup.
1188  unsigned box_index = CalculateContainingBox(rNodes[i]);
1189 
1190  if (IsBoxOwned(box_index))
1191  {
1192  rNodes[i]->ClearNeighbours();
1193  rNodes[i]->SetNeighboursSetUp(false);
1194  }
1195  }
1196 
1197  for (unsigned box_index=mMinBoxIndex; box_index<=mMaxBoxIndex; box_index++)
1198  {
1199  if (IsInteriorBox(box_index))
1200  {
1201  AddPairsFromBox(box_index, rNodePairs);
1202  }
1203  }
1204 
1205  if (mCalculateNodeNeighbours)
1206  {
1207  for (unsigned i = 0; i < rNodes.size(); i++)
1208  {
1209  // Get the box containing this node as only nodes on this process have NodeAttributes
1210  // and therefore Neighbours setup.
1211  unsigned box_index = CalculateContainingBox(rNodes[i]);
1212 
1213  if (IsBoxOwned(box_index))
1214  {
1215  rNodes[i]->RemoveDuplicateNeighbours();
1216  rNodes[i]->SetNeighboursSetUp(true);
1217  }
1218  }
1219  }
1220 }
1221 
1222 template<unsigned DIM>
1223 void DistributedBoxCollection<DIM>::CalculateBoundaryNodePairs(std::vector<Node<DIM>*>& rNodes, std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
1224 {
1225  for (unsigned box_index=mMinBoxIndex; box_index<=mMaxBoxIndex; box_index++)
1226  {
1227  if (!IsInteriorBox(box_index))
1228  {
1229  AddPairsFromBox(box_index, rNodePairs);
1230  }
1231  }
1232 
1233  if (mCalculateNodeNeighbours)
1234  {
1235  for (unsigned i = 0; i < rNodes.size(); i++)
1236  {
1237  // Get the box containing this node as only nodes on this process have NodeAttributes
1238  // and therefore Neighbours setup.
1239  unsigned box_index = CalculateContainingBox(rNodes[i]);
1240 
1241  if (IsBoxOwned(box_index))
1242  {
1243  rNodes[i]->RemoveDuplicateNeighbours();
1244  rNodes[i]->SetNeighboursSetUp(true);
1245  }
1246  }
1247  }
1248 }
1249 
1250 template<unsigned DIM>
1252  std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
1253 {
1254  // Get the box
1255  Box<DIM>& r_box = rGetBox(boxIndex);
1256 
1257  // Get the set of nodes in this box
1258  const std::set< Node<DIM>* >& r_contained_nodes = r_box.rGetNodesContained();
1259 
1260  // Get the local boxes to this box
1261  const std::set<unsigned>& local_boxes_indices = rGetLocalBoxes(boxIndex);
1262 
1263  // Loop over all the local boxes
1264  for (std::set<unsigned>::iterator box_iter = local_boxes_indices.begin();
1265  box_iter != local_boxes_indices.end();
1266  box_iter++)
1267  {
1268  Box<DIM>* p_neighbour_box;
1269 
1270  // Establish whether box is locally owned or halo.
1271  if (IsBoxOwned(*box_iter))
1272  {
1273  p_neighbour_box = &mBoxes[*box_iter - mMinBoxIndex];
1274  }
1275  else // Assume it is a halo.
1276  {
1277  p_neighbour_box = &mHaloBoxes[mHaloBoxesMapping[*box_iter]];
1278  }
1279  assert(p_neighbour_box);
1280 
1281  // Get the set of nodes contained in this box
1282  std::set< Node<DIM>* >& r_contained_neighbour_nodes = p_neighbour_box->rGetNodesContained();
1283 
1284  // Loop over these nodes
1285  for (typename std::set<Node<DIM>*>::iterator neighbour_node_iter = r_contained_neighbour_nodes.begin();
1286  neighbour_node_iter != r_contained_neighbour_nodes.end();
1287  ++neighbour_node_iter)
1288  {
1289  // Get the index of the other node
1290  unsigned other_node_index = (*neighbour_node_iter)->GetIndex();
1291 
1292  // Loop over nodes in this box
1293  for (typename std::set<Node<DIM>*>::iterator node_iter = r_contained_nodes.begin();
1294  node_iter != r_contained_nodes.end();
1295  ++node_iter)
1296  {
1297  unsigned node_index = (*node_iter)->GetIndex();
1298 
1299  // If we're in the same box, then take care not to store the node pair twice
1300  if (*box_iter != boxIndex || other_node_index > node_index)
1301  {
1302  rNodePairs.push_back(std::pair<Node<DIM>*, Node<DIM>*>((*node_iter), (*neighbour_node_iter)));
1303  if (mCalculateNodeNeighbours)
1304  {
1305  (*node_iter)->AddNeighbour(other_node_index);
1306  (*neighbour_node_iter)->AddNeighbour(node_index);
1307  }
1308  }
1309 
1310  }
1311  }
1312  }
1313 }
1314 
1315 template<unsigned DIM>
1317 {
1318  std::vector<int> cell_numbers(mpDistributedBoxStackFactory->GetHigh() - mpDistributedBoxStackFactory->GetLow(), 0);
1319 
1320  for (unsigned global_index=mMinBoxIndex; global_index<=mMaxBoxIndex; global_index++)
1321  {
1322  c_vector<unsigned, DIM> coords = CalculateGridIndices(global_index);
1323  unsigned location_in_vector = coords[DIM-1] - mpDistributedBoxStackFactory->GetLow();
1324  unsigned local_index = global_index - mMinBoxIndex;
1325  cell_numbers[location_in_vector] += mBoxes[local_index].rGetNodesContained().size();
1326  }
1327 
1328  return cell_numbers;
1329 }
1330 
1332 
1333 template class DistributedBoxCollection<1>;
1334 template class DistributedBoxCollection<2>;
1335 template class DistributedBoxCollection<3>;
std::vector< Box< DIM > > mBoxes
bool IsBoxOwned(unsigned globalIndex)
std::vector< unsigned > & rGetHaloNodesRight()
unsigned CalculateContainingBox(Node< DIM > *pNode)
void CalculateInteriorNodePairs(std::vector< Node< DIM > * > &rNodes, std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs)
Definition: Node.hpp:58
bool IsInteriorBox(unsigned globalIndex)
c_vector< double, 2 *DIM > rGetDomainSize() const
c_vector< double, 2 *DIM > mDomainSize
#define EXCEPTION(message)
Definition: Exception.hpp:143
void SetCalculateNodeNeighbours(bool calculateNodeNeighbours)
unsigned CalculateGlobalIndex(c_vector< unsigned, DIM > gridIndices)
void CalculateBoundaryNodePairs(std::vector< Node< DIM > * > &rNodes, std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs)
static bool AmMaster()
Definition: PetscTools.cpp:120
std::set< Node< DIM > * > & rGetNodesContained()
Definition: Box.cpp:56
std::vector< int > CalculateNumberOfNodesInEachStrip()
#define NEVER_REACHED
Definition: Exception.hpp:206
unsigned GetProcessOwningNode(Node< DIM > *pNode)
std::set< unsigned > & rGetLocalBoxes(unsigned boxIndex)
static bool IsSequential()
Definition: PetscTools.cpp:91
c_vector< unsigned, DIM > mNumBoxesEachDirection
std::vector< unsigned > & rGetHaloNodesLeft()
bool IsOwned(Node< DIM > *pNode)
DistributedBoxCollection(double boxWidth, c_vector< double, 2 *DIM > domainSize, bool isPeriodicInX=false, int localRows=PETSC_DECIDE)
void CalculateNodePairs(std::vector< Node< DIM > * > &rNodes, std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs)
int LoadBalance(std::vector< int > localDistribution)
Definition: Box.hpp:50
static bool IsParallel()
Definition: PetscTools.cpp:97
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition: Node.cpp:139
c_vector< unsigned, DIM > CalculateGridIndices(unsigned globalIndex)
bool IsHaloBox(unsigned globalIndex)
Box< DIM > & rGetHaloBox(unsigned boxIndex)
static bool AmTopMost()
Definition: PetscTools.cpp:126
static unsigned GetNumProcs()
Definition: PetscTools.cpp:108
static unsigned GetMyRank()
Definition: PetscTools.cpp:114
Box< DIM > & rGetBox(unsigned boxIndex)
DistributedVectorFactory * mpDistributedBoxStackFactory
void AddPairsFromBox(unsigned boxIndex, std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs)