Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
DistributedBoxCollection.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35#include "DistributedBoxCollection.hpp"
36#include "Exception.hpp"
38#include "Warnings.hpp"
39
40// Static member for "fudge factor" is instantiated here
41template<unsigned DIM>
43
44template<unsigned DIM>
45DistributedBoxCollection<DIM>::DistributedBoxCollection(double boxWidth, c_vector<double, 2*DIM> domainSize, bool isPeriodicInX, bool isPeriodicInY, bool isPeriodicInZ, int localRows)
46 : mBoxWidth(boxWidth),
47 mIsPeriodicInX(isPeriodicInX),
48 mIsPeriodicInY(isPeriodicInY),
49 mIsPeriodicInZ(isPeriodicInZ),
50 mAreLocalBoxesSet(false),
51 mCalculateNodeNeighbours(true)
52{
53 // If the domain size is not 'divisible' (i.e. fmod(width, box_size) > 0.0) we swell the domain to enforce this.
54 for (unsigned i=0; i<DIM; i++)
55 {
56 double r = fmod((domainSize[2*i+1]-domainSize[2*i]), boxWidth);
57 if (r > 0.0)
58 {
59 domainSize[2*i+1] += boxWidth - r;
60 }
61 }
62
63 mDomainSize = domainSize;
64
65 // Set the periodicity across procs flag
66 mIsPeriodicAcrossProcs = (DIM==1 && mIsPeriodicInX) || (DIM==2 && mIsPeriodicInY) || (DIM==3 && mIsPeriodicInZ);
67
68 // Calculate the number of boxes in each direction.
69 mNumBoxesEachDirection = scalar_vector<unsigned>(DIM, 0u);
70
71 for (unsigned i=0; i<DIM; i++)
72 {
73 double counter = mDomainSize(2*i);
74 while (counter + msFudge < mDomainSize(2*i+1))
75 {
77 counter += mBoxWidth;
78 }
79 }
80
81 // Make sure there are enough boxes for the number of processes.
83 {
84 WARNING("There are more processes than convenient for the domain/mesh/box size. The domain size has been swollen.")
87 }
88
89 // Make a distributed vector factory to split the rows of boxes between processes.
91
92 // Calculate how many boxes in a row / face. A useful piece of data in the class.
93 mNumBoxes = 1u;
94 for (unsigned dim=0; dim<DIM; dim++)
95 {
97 }
98
100
101 unsigned num_local_boxes = mNumBoxesInAFace * GetNumLocalRows();
102
105
106 // Create the correct number of boxes and set up halos
107 mBoxes.resize(num_local_boxes);
109}
110
111template<unsigned DIM>
113{
114 delete mpDistributedBoxStackFactory;
115}
116
117template<unsigned DIM>
119{
120 for (unsigned i=0; i<mBoxes.size(); i++)
121 {
122 mBoxes[i].ClearNodes();
123 }
124 for (unsigned i=0; i<mHaloBoxes.size(); i++)
125 {
126 mHaloBoxes[i].ClearNodes();
127 }
128}
129
130template<unsigned DIM>
132{
133 // We don't need to do this if not parallel
135 {
136 return;
137 }
138
139 // Get top-most and bottom-most value of Distributed Box Stack.
140 unsigned hi = mpDistributedBoxStackFactory->GetHigh();
141 unsigned lo = mpDistributedBoxStackFactory->GetLow();
142
143 // If I am not the top-most process, add halo structures above.
145 {
146 for (unsigned i=0; i < mNumBoxesInAFace; i++)
147 {
148 Box<DIM> new_box;
149 mHaloBoxes.push_back(new_box);
150
151 unsigned global_index = hi * mNumBoxesInAFace + i;
152 mHaloBoxesMapping[global_index] = mHaloBoxes.size()-1;
153 mHalosRight.push_back(global_index - mNumBoxesInAFace);
154 }
155 }
156 // Otherwise if I am the top most and periodic in y (2d) or z (3d) add halo boxes for
157 // the base process
158 else if ( mIsPeriodicAcrossProcs )
159 {
160 for (unsigned i=0; i < mNumBoxesInAFace; i++)
161 {
162 Box<DIM> new_box;
163 mHaloBoxes.push_back(new_box);
164
165 mHaloBoxesMapping[i] = mHaloBoxes.size()-1;
167 mHalosRight.push_back( (hi-1)*mNumBoxesInAFace + i );
168 }
169 }
170
171 // If I am not the bottom-most process, add halo structures below.
173 {
174 for (unsigned i=0; i< mNumBoxesInAFace; i++)
175 {
176 Box<DIM> new_box;
177 mHaloBoxes.push_back(new_box);
178
179 unsigned global_index = (lo - 1) * mNumBoxesInAFace + i;
180 mHaloBoxesMapping[global_index] = mHaloBoxes.size() - 1;
181
182 mHalosLeft.push_back(global_index + mNumBoxesInAFace);
183 }
184 }
185 // Otherwise if I am the bottom most and periodic in y (2d) or z (3d) add halo boxes for
186 // the top process
187 else if ( mIsPeriodicAcrossProcs )
188 {
189 for (unsigned i=0; i < mNumBoxesInAFace; i++)
190 {
191 Box<DIM> new_box;
192 mHaloBoxes.push_back(new_box);
194 unsigned global_index = (mNumBoxesEachDirection(DIM-1) - 1) * mNumBoxesInAFace + i;
195 mHaloBoxesMapping[global_index] = mHaloBoxes.size()-1;
196
197 mHalosLeft.push_back( i );
198 }
199
200 }
201}
203template<unsigned DIM>
205{
206 mHaloNodesLeft.clear();
207 for (unsigned i=0; i<mHalosLeft.size(); i++)
209 for (typename std::set<Node<DIM>* >::iterator iter=this->rGetBox(mHalosLeft[i]).rGetNodesContained().begin();
210 iter!=this->rGetBox(mHalosLeft[i]).rGetNodesContained().end();
211 iter++)
212 {
213 mHaloNodesLeft.push_back((*iter)->GetIndex());
214 }
215 }
217 // Send right
218 mHaloNodesRight.clear();
219 for (unsigned i=0; i<mHalosRight.size(); i++)
220 {
221 for (typename std::set<Node<DIM>* >::iterator iter=this->rGetBox(mHalosRight[i]).rGetNodesContained().begin();
222 iter!=this->rGetBox(mHalosRight[i]).rGetNodesContained().end();
223 iter++)
224 {
225 mHaloNodesRight.push_back((*iter)->GetIndex());
226 }
227 }
229
230template<unsigned DIM>
232{
233 return mpDistributedBoxStackFactory->GetHigh() - mpDistributedBoxStackFactory->GetLow();
234}
235
236template<unsigned DIM>
238{
239 return (!(globalIndex<mMinBoxIndex) && !(mMaxBoxIndex<globalIndex));
240}
241
242template<unsigned DIM>
245 bool is_halo_right = ((globalIndex > mMaxBoxIndex) && !(globalIndex > mMaxBoxIndex + mNumBoxesInAFace));
246 bool is_halo_left = ((globalIndex < mMinBoxIndex) && !(globalIndex < mMinBoxIndex - mNumBoxesInAFace));
247
248 // Also need to check for periodic boxes
249 if ( mIsPeriodicAcrossProcs )
250 {
252 {
253 is_halo_right = (globalIndex < mNumBoxesInAFace);
254 }
255 if ( PetscTools::AmMaster() )
257 is_halo_left = (!(globalIndex < (mNumBoxesEachDirection[DIM-1]-1)*mNumBoxesInAFace) && (globalIndex < mNumBoxesEachDirection[DIM-1]*mNumBoxesInAFace ));
258 }
259 }
260
261 return (PetscTools::IsParallel() && (is_halo_right || is_halo_left));
262}
263
264template<unsigned DIM>
267 bool is_on_boundary = !(globalIndex < mMaxBoxIndex - mNumBoxesInAFace) || (globalIndex < mMinBoxIndex + mNumBoxesInAFace);
268
269 return (PetscTools::IsSequential() || !(is_on_boundary));
270}
272template<unsigned DIM>
273unsigned DistributedBoxCollection<DIM>::CalculateGlobalIndex(c_vector<unsigned, DIM> gridIndices)
274{
276 unsigned global_index;
277
278 switch (DIM)
279 {
280 case 1:
282 global_index = gridIndices(0);
283 break;
284 }
285 case 2:
287 global_index = gridIndices(0) +
288 gridIndices(1) * mNumBoxesEachDirection(0);
289 break;
290 }
291 case 3:
292 {
293 global_index = gridIndices(0) +
294 gridIndices(1) * mNumBoxesEachDirection(0) +
295 gridIndices(2) * mNumBoxesEachDirection(0) * mNumBoxesEachDirection(1);
296 break;
297 }
298 default:
299 {
301 }
303 return global_index;
304}
305
306template<unsigned DIM>
308{
309 // Get the location of the node
310 c_vector<double, DIM> location = pNode->rGetLocation();
311 return CalculateContainingBox(location);
312}
313
314template<unsigned DIM>
315unsigned DistributedBoxCollection<DIM>::CalculateContainingBox(c_vector<double, DIM>& rLocation)
317 // The node must lie inside the boundary of the box collection
318 for (unsigned i=0; i<DIM; i++)
319 {
320 if ((rLocation[i] < mDomainSize(2*i)) || !(rLocation[i] < mDomainSize(2*i+1)))
321 {
322 EXCEPTION("The point provided is outside all of the boxes");
324 }
325
326 // Compute the containing box index in each dimension
327 c_vector<unsigned, DIM> containing_box_indices = scalar_vector<unsigned>(DIM, 0u);
328 for (unsigned i=0; i<DIM; i++)
329 {
330 double box_counter = mDomainSize(2*i);
331 while (!((box_counter + mBoxWidth) > rLocation[i] + msFudge))
332 {
333 containing_box_indices[i]++;
334 box_counter += mBoxWidth;
335 }
337
338 // Use these to compute the index of the containing box
339 unsigned containing_box_index = 0;
340 for (unsigned i=0; i<DIM; i++)
341 {
342 unsigned temp = 1;
343 for (unsigned j=0; j<i; j++)
344 {
345 temp *= mNumBoxesEachDirection(j);
346 }
347 containing_box_index += temp*containing_box_indices[i];
349
350 // This index must be less than the total number of boxes
351 assert(containing_box_index < mNumBoxes);
352
353 return containing_box_index;
354}
355
356template<unsigned DIM>
357c_vector<unsigned, DIM> DistributedBoxCollection<DIM>::CalculateGridIndices(unsigned globalIndex)
358{
359 c_vector<unsigned, DIM> grid_indices;
360
361 switch (DIM)
363 case 1:
364 {
365 grid_indices(0) = globalIndex;
366 break;
368 case 2:
369 {
370 unsigned num_x = mNumBoxesEachDirection(0);
371 grid_indices(0) = globalIndex % num_x;
372 grid_indices(1) = (globalIndex - grid_indices(0)) / num_x;
373 break;
375 case 3:
376 {
377 unsigned num_x = mNumBoxesEachDirection(0);
378 unsigned num_xy = mNumBoxesEachDirection(0)*mNumBoxesEachDirection(1);
379 grid_indices(0) = globalIndex % num_x;
380 grid_indices(1) = (globalIndex % num_xy - grid_indices(0)) / num_x;
381 grid_indices(2) = globalIndex / num_xy;
382 break;
383 }
384 default:
387 }
388 }
389
390 return grid_indices;
391}
392
393template<unsigned DIM>
395{
396 // Check first for local ownership
397 if (!(boxIndex<mMinBoxIndex) && !(mMaxBoxIndex<boxIndex))
398 {
399 return mBoxes[boxIndex-mMinBoxIndex];
400 }
402 // If normal execution reaches this point then the box does not belong to the process so we will check for a halo box
403 return rGetHaloBox(boxIndex);
404}
405
406template<unsigned DIM>
408{
409 assert(IsHaloBox(boxIndex));
410
411 unsigned local_index = mHaloBoxesMapping.find(boxIndex)->second;
412
413 return mHaloBoxes[local_index];
414}
415
416template<unsigned DIM>
418{
419 return mNumBoxes;
420}
421
422template<unsigned DIM>
424{
425 return mBoxes.size();
426}
427
428template<unsigned DIM>
430{
431 return mDomainSize;
432}
433
434template<unsigned DIM>
436{
437 return mAreLocalBoxesSet;
438}
439
440template<unsigned DIM>
442{
443 return mBoxWidth;
444}
445
446template<unsigned DIM>
448{
449 return mIsPeriodicInX;
450}
451
452template<unsigned DIM>
454{
455 return mIsPeriodicInY;
456}
457
458template<unsigned DIM>
460{
461 return mIsPeriodicInZ;
462}
463
464template<unsigned DIM>
466{
467 return mIsPeriodicAcrossProcs;
468}
469
470template<unsigned DIM>
472{
473 c_vector<bool, DIM> periodic_dims = zero_vector<double>(DIM);
474
475 if constexpr (DIM == 1)
476 {
477 periodic_dims[0] = mIsPeriodicInX;
478 }
479 else if constexpr (DIM == 2)
480 {
481 periodic_dims[0] = mIsPeriodicInX;
482 periodic_dims[1] = mIsPeriodicInY;
483 }
484 else if constexpr (DIM == 3)
485 {
486 periodic_dims[0] = mIsPeriodicInX;
487 periodic_dims[1] = mIsPeriodicInY;
488 periodic_dims[2] = mIsPeriodicInZ;
489 }
490
491 return periodic_dims;
492}
493
494template<unsigned DIM>
496{
497 return mpDistributedBoxStackFactory->GetHigh() - mpDistributedBoxStackFactory->GetLow();
498}
499
500template<unsigned DIM>
501int DistributedBoxCollection<DIM>::LoadBalance(std::vector<int> localDistribution)
502{
503 MPI_Status status;
504
505 int proc_right = (PetscTools::AmTopMost()) ? MPI_PROC_NULL : (int)PetscTools::GetMyRank() + 1;
506 int proc_left = (PetscTools::AmMaster()) ? MPI_PROC_NULL : (int)PetscTools::GetMyRank() - 1;
507
508 // A variable that will return the new number of rows.
509 int new_rows = localDistribution.size();
510
514 unsigned rows_on_left_process = 0;
515 std::vector<int> node_distr_on_left_process;
516
517 unsigned num_local_rows = localDistribution.size();
518
519 MPI_Send(&num_local_rows, 1, MPI_UNSIGNED, proc_right, 123, PETSC_COMM_WORLD);
520 MPI_Recv(&rows_on_left_process, 1, MPI_UNSIGNED, proc_left, 123, PETSC_COMM_WORLD, &status);
521
522 node_distr_on_left_process.resize(rows_on_left_process > 0 ? rows_on_left_process : 1);
523
524 MPI_Send(&localDistribution[0], num_local_rows, MPI_INT, proc_right, 123, PETSC_COMM_WORLD);
525 MPI_Recv(&node_distr_on_left_process[0], rows_on_left_process, MPI_INT, proc_left, 123, PETSC_COMM_WORLD, &status);
526
530 int local_load = 0;
531 for (unsigned i=0; i<localDistribution.size(); i++)
532 {
533 local_load += localDistribution[i];
534 }
535 int load_on_left_proc = 0;
536 for (unsigned i=0; i<node_distr_on_left_process.size(); i++)
537 {
538 load_on_left_proc += node_distr_on_left_process[i];
539 }
540
542 {
543 // Calculate (Difference in load with a shift) - (Difference in current loads) for a move left and right of the boundary
544 // This code uses integer arithmetic in order to avoid the rounding errors associated with doubles
545 int local_to_left_sq = (local_load - load_on_left_proc) * (local_load - load_on_left_proc);
546 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]) );
547 delta_left = delta_left*delta_left - local_to_left_sq;
548
549 int delta_right = ( (local_load - localDistribution[0]) - (load_on_left_proc + localDistribution[0]));
550 delta_right = delta_right*delta_right - local_to_left_sq;
551
552 // If a delta is negative we should accept that change. If both are negative choose the largest change.
553 int local_change = 0;
554 bool move_left = (!(delta_left > 0) && (node_distr_on_left_process.size() > 1));
555 if (move_left)
556 {
557 local_change = 1;
558 }
559
560 bool move_right = !(delta_right > 0) && (localDistribution.size() > 2);
561 if (move_right)
562 {
563 local_change = -1;
564 }
565
566 if (move_left && move_right)
567 {
568 local_change = (fabs((double)delta_right) > fabs((double)delta_left)) ? -1 : 1;
569 }
570
571 // Update the number of local rows.
572 new_rows += local_change;
573
574 // Send the result of the calculation back to the left processes.
575 MPI_Send(&local_change, 1, MPI_INT, proc_left, 123, PETSC_COMM_WORLD);
576 }
577
578 // Receive changes from right hand process.
579 int remote_change = 0;
580 MPI_Recv(&remote_change, 1, MPI_INT, proc_right, 123, PETSC_COMM_WORLD, &status);
581
582 // Update based on change or right/top boundary
583 new_rows -= remote_change;
584
585 return new_rows;
586}
587
588template<unsigned DIM>
590{
591 if (mAreLocalBoxesSet)
592 {
593 EXCEPTION("Local Boxes Are Already Set");
594 }
595 else
596 {
597 switch (DIM)
598 {
599 case 1:
600 {
601 // We only need to look for neighbours in the current and successive boxes plus some others for halos
602 mLocalBoxes.clear();
603
604 // Iterate over the global box indices
605 for (unsigned global_index = mMinBoxIndex; global_index<mMaxBoxIndex+1; global_index++)
606 {
607 std::set<unsigned> local_boxes;
608
609 // Insert the current box
610 local_boxes.insert(global_index);
611
612 // Set some bools to find out where we are
613 bool right = (global_index==mNumBoxesEachDirection(0)-1);
614 bool left = (global_index == 0);
615 bool proc_left = (global_index == mpDistributedBoxStackFactory->GetLow());
616
617 // If we're not at the right-most box, then insert the box to the right
618 if (!right)
619 {
620 local_boxes.insert(global_index+1);
621 }
622 // If we're on a left process boundary and not on process 0, insert the (halo) box to the left
623 if (proc_left && !left)
624 {
625 local_boxes.insert(global_index-1);
626 }
627
628 mLocalBoxes.push_back(local_boxes);
629 }
630 break;
631 }
632 case 2:
633 {
634 // Clear the local boxes
635 mLocalBoxes.clear();
636
637 // Now we need to work out the min and max y indices
638 unsigned j_start = CalculateGridIndices(mMinBoxIndex)(1);
639 unsigned j_end = CalculateGridIndices(mMaxBoxIndex)(1);
640 // Determine the number of boxes in each direction
641 unsigned nI = mNumBoxesEachDirection(0);
642 unsigned nJ = mNumBoxesEachDirection(1);
643 // Locally store the bottom processor row
644 unsigned bottom_proc = mpDistributedBoxStackFactory->GetLow();
645 unsigned top_proc = mpDistributedBoxStackFactory->GetHigh();
646
647 // Outer loop: over y
648 for ( unsigned j = j_start; j < (j_end+1); j++ )
649 {
650 // Inner loop: over x
651 for ( unsigned i = 0; i < nI; i++ )
652 {
653 std::set<unsigned> local_boxes;
654 /* We want to add (for non-boundary)
655 (i-1,j+1) (i,j+1) (i+1,j+1)
656 (i,j) (i+1,j) */
657
658 int j_mod = j; // This is used to account for y periodic boundaries
659 // If we are on the bottom of the processor, we may need to add the row below
660 int dj = -1 * (int)(j == bottom_proc && (j > 0 || (mIsPeriodicInY && top_proc < nJ)) );
661 int j_mod_2 = 0;
662 // The min ensures we don't go above the top boundary
663 for (; dj < std::min((int)nJ-j_mod,(int)2); dj++ )
664 {
665 // We need to change to the top row if we are in the condition where dj == -1 and j == 0 which is only
666 // when we are periodic in y and the top row is on a different processor to the bottom row
667 if ( mIsPeriodicInY && j == 0 && dj < 1 && top_proc < nJ )
668 {
669 j_mod_2 = ( dj < 0 ) ? nJ : 0;
670 }
671
672 // The -1*dj ensures we get the upper left, the max ensures we don't hit the left boundary,
673 // the min ensures we don't hit the right boundary
674 int boxi = std::max((int) i-1*std::abs(dj),(int) 0);
675 for ( ; boxi < std::min((int)i+2,(int)nI); boxi++ )
676 {
677 local_boxes.insert( (j_mod+dj+j_mod_2)*nI + boxi );
678 }
679 // Add in the x periodicity at the right boundary, we then want: (0,j) and (0,j+1)
680 if ( i==(nI-1) && mIsPeriodicInX )
681 {
682 local_boxes.insert( (j_mod+dj+j_mod_2)*nI );
683 }
684 // If the y boundary is periodic, the new j level is the bottom row; we use jmod to adjust
685 // for this to adjust for dj = 1
686 if ( mIsPeriodicInY && j == (nJ-1) && dj == 0 )
687 {
688 j_mod = -1;
689 }
690 }
691 // Now add the left-upper box if x periodic and on the boundary
692 if ( mIsPeriodicInX && i == 0 )
693 {
694 if ( j < (nJ-1) )
695 {
696 local_boxes.insert( (j+2)*nI - 1 );
697 if ( j==0 && mIsPeriodicInY && top_proc < nJ )
698 {
699 local_boxes.insert( nI*nJ - 1 );
700 }
701 }
702 else if ( mIsPeriodicInY )
703 {
704 local_boxes.insert( nI-1 );
705 }
706 // If wee are on the bottom of the process need to add
707 // the box on the right in the row below
708 if ( j == bottom_proc && j > 0 )
709 {
710 local_boxes.insert( j*nI - 1 );
711 }
712 }
713
714 // Add to the local boxes
715 mLocalBoxes.push_back(local_boxes);
716 }
717 }
718
719 break;
720 }
721 case 3:
722 {
723 // Clear the local boxes
724 mLocalBoxes.clear();
725
726 // Now we need to work out the min and max z indices
727 unsigned k_start = CalculateGridIndices(mMinBoxIndex)(2);
728 unsigned k_end = CalculateGridIndices(mMaxBoxIndex)(2);
729
730 // Determine the number of boxes in each direction
731 unsigned nI = mNumBoxesEachDirection(0);
732 unsigned nJ = mNumBoxesEachDirection(1);
733 unsigned nK = mNumBoxesEachDirection(2);
734
735 // Work out the bottom/top processor
736 unsigned bottom_proc = mpDistributedBoxStackFactory->GetLow();
737 unsigned top_proc = mpDistributedBoxStackFactory->GetHigh();
738
739 // Outer loop: over z
740 for ( unsigned k = k_start; k <= k_end; k++ )
741 {
742 // Middle loop: over y
743 for ( unsigned j = 0; j < nJ; j++ )
744 {
745 // Inner loop: over x
746 for ( unsigned i = 0; i < nI; i++ )
747 {
748 std::set<unsigned> local_boxes;
749
750 // Same z level
751 unsigned z_offset = k*nI*nJ;
752 // (See case dim=2 for commented code of X-Y implementation)
753 int j_mod = (int)j;
754 for ( int dj = 0; dj < std::min((int)nJ-j_mod,2); dj++ )
755 {
756 for ( int boxi = std::max((int)i-1*dj,0);
757 boxi < std::min((int)i+2,(int)nI); boxi++ )
758 {
759 local_boxes.insert( z_offset + (j_mod+dj)*nI + boxi );
760 }
761 if ( i==(nI-1) && mIsPeriodicInX )
762 {
763 local_boxes.insert( z_offset + (j_mod+dj)*nI );
764 }
765 if ( mIsPeriodicInY && j == (nJ-1) )
766 {
767 j_mod = -1;
768 }
769 }
770 // Now add the left-upper box if x periodic and on the boundary
771 if ( mIsPeriodicInX && i == 0 )
772 {
773 if ( j < (nJ-1) )
774 {
775 local_boxes.insert( z_offset + (j+2)*nI - 1 );
776 }
777 else if ( mIsPeriodicInY )
778 {
779 local_boxes.insert( z_offset + nI-1 );
780 }
781 }
782
783 // Add all the surrounding boxes from the z level above if required
784 std::vector<unsigned> k_offset;
785 if ( k < (nK-1) )
786 {
787 k_offset.push_back(k+1);
788 }
789 if ( k == bottom_proc && k > 0 )
790 {
791 // Need the level below if we are on the lowest processor
792 k_offset.push_back(k-1);
793 }
794 if ( mIsPeriodicInZ && k == (nK-1) )
795 {
796 k_offset.push_back(0);
797 }
798 else if ( mIsPeriodicInZ && k==0 && (top_proc < nK) )
799 {
800 k_offset.push_back(nK-1);
801 }
802
803 for ( std::vector<unsigned>::iterator k_offset_it = k_offset.begin(); k_offset_it != k_offset.end(); ++k_offset_it )
804 {
805 z_offset = (*k_offset_it)*nI*nJ;
806 // Periodicity adjustments
807 int pX = (int) mIsPeriodicInX;
808 int pY = (int) mIsPeriodicInY;
809 for ( int boxi = std::max((int)i-1,-1*pX); boxi < std::min((int)i+2,(int)nI+pX); boxi++ )
810 {
811 for ( int boxj = std::max((int)j-1,-1*pY); boxj < std::min((int)j+2,(int)nJ+pY); boxj++ )
812 {
813 int box_to_add = z_offset + (boxj*(int)nI) + boxi;
814 // Check for x periodicity
815 // i==(nI-1) only when we are periodic and on the right,
816 // i==-1 only when we are periodic and on the left
817 box_to_add += ( (unsigned)(boxi==-1) - (unsigned)(boxi==(int)nI) )*nI;
818 // Check for y periodicity
819 // j==nJ only when we are periodic and on the right,
820 // j==-1 only when we are periodic and on the left
821 box_to_add += ( (unsigned)(boxj==-1) - (unsigned)(boxj==(int)nJ) )*nI*nJ;
822 local_boxes.insert( box_to_add );
823 }
824 }
825 }
826
827 // Add to the local boxes
828 mLocalBoxes.push_back(local_boxes);
829 }
830 }
831 }
832 break;
833 }
834 default:
836 }
837 mAreLocalBoxesSet=true;
838 }
839}
840
841template<unsigned DIM>
843{
844 mAreLocalBoxesSet = true;
845 switch (DIM)
846 {
847 case 1:
848 {
849 for (unsigned i=mMinBoxIndex; i<mMaxBoxIndex+1; i++)
850 {
851 std::set<unsigned> local_boxes;
852
853 local_boxes.insert(i);
854
855 // add the two neighbours
856 if (i!=0)
857 {
858 local_boxes.insert(i-1);
859 }
860 if (i+1 != mNumBoxesEachDirection(0))
861 {
862 local_boxes.insert(i+1);
863 }
864
865 mLocalBoxes.push_back(local_boxes);
866 }
867 break;
868 }
869 case 2:
870 {
871 mLocalBoxes.clear();
872
873 unsigned M = mNumBoxesEachDirection(0);
874 unsigned N = mNumBoxesEachDirection(1);
875
876 std::vector<bool> is_xmin(N*M); // far left
877 std::vector<bool> is_xmax(N*M); // far right
878 std::vector<bool> is_ymin(N*M); // bottom
879 std::vector<bool> is_ymax(N*M); // top
880
881 for (unsigned i=0; i<M*N; i++)
882 {
883 is_xmin[i] = (i%M==0);
884 is_xmax[i] = ((i+1)%M==0);
885 is_ymin[i] = (i%(M*N)<M);
886 is_ymax[i] = (i%(M*N)>=(N-1)*M);
887 }
888
889 for (unsigned i=mMinBoxIndex; i<mMaxBoxIndex+1; i++)
890 {
891 std::set<unsigned> local_boxes;
892
893 local_boxes.insert(i);
894
895 // add the box to the left
896 if (!is_xmin[i])
897 {
898 local_boxes.insert(i-1);
899 }
900 else // Add Periodic Box if needed
901 {
902 if (mIsPeriodicInX)
903 {
904 local_boxes.insert(i+M-1);
905 }
906 }
907
908 // add the box to the right
909 if (!is_xmax[i])
910 {
911 local_boxes.insert(i+1);
912 }
913 else // Add Periodic Box if needed
914 {
915 if (mIsPeriodicInX)
916 {
917 local_boxes.insert(i-M+1);
918 }
919 }
920
921 // add the one below
922 if (!is_ymin[i])
923 {
924 local_boxes.insert(i-M);
925 }
926 else // Add periodic box if needed
927 {
928 if (mIsPeriodicInY)
929 {
930 local_boxes.insert(i+(N-1)*M);
931 }
932 }
933
934 // add the one above
935 if (!is_ymax[i])
936 {
937 local_boxes.insert(i+M);
938 }
939 else // Add periodic box if needed
940 {
941 if (mIsPeriodicInY)
942 {
943 local_boxes.insert(i-(N-1)*M);
944 }
945 }
946
947 // add the four corner boxes
948
949 if ( (!is_xmin[i]) && (!is_ymin[i]) )
950 {
951 local_boxes.insert(i-1-M);
952 }
953 if ( (!is_xmin[i]) && (!is_ymax[i]) )
954 {
955 local_boxes.insert(i-1+M);
956 }
957 if ( (!is_xmax[i]) && (!is_ymin[i]) )
958 {
959 local_boxes.insert(i+1-M);
960 }
961 if ( (!is_xmax[i]) && (!is_ymax[i]) )
962 {
963 local_boxes.insert(i+1+M);
964 }
965
966 // Add periodic corner boxes if needed
967 if (mIsPeriodicInX)
968 {
969 if ( (is_xmin[i]) && (!is_ymin[i]) )
970 {
971 local_boxes.insert(i-1);
972 }
973 if ( (is_xmin[i]) && (!is_ymax[i]) )
974 {
975 local_boxes.insert(i-1+2*M);
976 }
977 if ( (is_xmax[i]) && (!is_ymin[i]) )
978 {
979 local_boxes.insert(i+1-2*M);
980 }
981 if ( (is_xmax[i]) && (!is_ymax[i]) )
982 {
983 local_boxes.insert(i+1);
984 }
985 }
986 if(mIsPeriodicInY)
987 {
988 if( (is_ymin[i]) && !(is_xmin[i]) )
989 {
990 local_boxes.insert(i+(N-1)*M-1);
991 }
992 if( (is_ymin[i]) && !(is_xmax[i]) )
993 {
994 local_boxes.insert(i+(N-1)*M+1);
995 }
996 if( (is_ymax[i]) && !(is_xmin[i]) )
997 {
998 local_boxes.insert(i-(N-1)*M-1);
999 }
1000 if( (is_ymax[i]) && !(is_xmax[i]) )
1001 {
1002 local_boxes.insert(i-(N-1)*M+1);
1003 }
1004 }
1005 if(mIsPeriodicInX && mIsPeriodicInY)
1006 {
1007 if( i==0 ) // Lower left corner
1008 {
1009 local_boxes.insert(M*N-1); // Add upper right corner
1010 }
1011 else if( i==(M-1) ) // Lower right corner
1012 {
1013 local_boxes.insert(M*(N-1)); // Add upper left corner
1014 }
1015 else if( i==(M*(N-1)) ) // Upper left corner
1016 {
1017 local_boxes.insert(M-1); // Add lower right corner
1018 }
1019 else if( i==(M*N-1) ) // Upper right corner
1020 {
1021 local_boxes.insert(0); // Lower left corner
1022 }
1023 }
1024
1025 mLocalBoxes.push_back(local_boxes);
1026 }
1027 break;
1028 }
1029 case 3:
1030 {
1031 mLocalBoxes.clear();
1032
1033 unsigned M = mNumBoxesEachDirection(0);
1034 unsigned N = mNumBoxesEachDirection(1);
1035 unsigned P = mNumBoxesEachDirection(2);
1036
1037 std::vector<bool> is_xmin(N*M*P); // far left
1038 std::vector<bool> is_xmax(N*M*P); // far right
1039 std::vector<bool> is_ymin(N*M*P); // nearest
1040 std::vector<bool> is_ymax(N*M*P); // furthest
1041 std::vector<bool> is_zmin(N*M*P); // bottom layer
1042 std::vector<bool> is_zmax(N*M*P); // top layer
1043
1044 for (unsigned i=0; i<M*N*P; i++)
1045 {
1046 is_xmin[i] = (i%M==0);
1047 is_xmax[i] = ((i+1)%M==0);
1048 is_ymin[i] = (i%(M*N)<M);
1049 is_ymax[i] = (i%(M*N)>=(N-1)*M);
1050 is_zmin[i] = (i<M*N);
1051 is_zmax[i] = (i>=M*N*(P-1));
1052 }
1053
1054 for (unsigned i=mMinBoxIndex; i<mMaxBoxIndex+1; i++)
1055 {
1056 std::set<unsigned> local_boxes;
1057
1058 // add itself as a local box
1059 local_boxes.insert(i);
1060
1061 // now add all 26 other neighbours.....
1062
1063 // add the box left
1064 if (!is_xmin[i])
1065 {
1066 local_boxes.insert(i-1);
1067
1068 // plus some others towards the left
1069 if (!is_ymin[i])
1070 {
1071 local_boxes.insert(i-1-M);
1072 }
1073
1074 if (!is_ymax[i])
1075 {
1076 local_boxes.insert(i-1+M);
1077 }
1078
1079 if (!is_zmin[i])
1080 {
1081 local_boxes.insert(i-1-M*N);
1082 }
1083
1084 if (!is_zmax[i])
1085 {
1086 local_boxes.insert(i-1+M*N);
1087 }
1088 }
1089
1090 // add the box to the right
1091 if (!is_xmax[i])
1092 {
1093 local_boxes.insert(i+1);
1094
1095 // plus some others towards the right
1096 if (!is_ymin[i])
1097 {
1098 local_boxes.insert(i+1-M);
1099 }
1100
1101 if (!is_ymax[i])
1102 {
1103 local_boxes.insert(i+1+M);
1104 }
1105
1106 if (!is_zmin[i])
1107 {
1108 local_boxes.insert(i+1-M*N);
1109 }
1110
1111 if (!is_zmax[i])
1112 {
1113 local_boxes.insert(i+1+M*N);
1114 }
1115 }
1116
1117 // add the boxes next along the y axis
1118 if (!is_ymin[i])
1119 {
1120 local_boxes.insert(i-M);
1121
1122 // and more in this plane
1123 if (!is_zmin[i])
1124 {
1125 local_boxes.insert(i-M-M*N);
1126 }
1127
1128 if (!is_zmax[i])
1129 {
1130 local_boxes.insert(i-M+M*N);
1131 }
1132 }
1133
1134 // add the boxes next along the y axis
1135 if (!is_ymax[i])
1136 {
1137 local_boxes.insert(i+M);
1138
1139 // and more in this plane
1140 if (!is_zmin[i])
1141 {
1142 local_boxes.insert(i+M-M*N);
1143 }
1144
1145 if (!is_zmax[i])
1146 {
1147 local_boxes.insert(i+M+M*N);
1148 }
1149 }
1150
1151 // add the box directly above
1152 if (!is_zmin[i])
1153 {
1154 local_boxes.insert(i-N*M);
1155 }
1156
1157 // add the box directly below
1158 if (!is_zmax[i])
1159 {
1160 local_boxes.insert(i+N*M);
1161 }
1162
1163 // finally, the 8 corners are left
1164
1165 if ( (!is_xmin[i]) && (!is_ymin[i]) && (!is_zmin[i]) )
1166 {
1167 local_boxes.insert(i-1-M-M*N);
1168 }
1169
1170 if ( (!is_xmin[i]) && (!is_ymin[i]) && (!is_zmax[i]) )
1171 {
1172 local_boxes.insert(i-1-M+M*N);
1173 }
1174
1175 if ( (!is_xmin[i]) && (!is_ymax[i]) && (!is_zmin[i]) )
1176 {
1177 local_boxes.insert(i-1+M-M*N);
1178 }
1179
1180 if ( (!is_xmin[i]) && (!is_ymax[i]) && (!is_zmax[i]) )
1181 {
1182 local_boxes.insert(i-1+M+M*N);
1183 }
1184
1185 if ( (!is_xmax[i]) && (!is_ymin[i]) && (!is_zmin[i]) )
1186 {
1187 local_boxes.insert(i+1-M-M*N);
1188 }
1189
1190 if ( (!is_xmax[i]) && (!is_ymin[i]) && (!is_zmax[i]) )
1191 {
1192 local_boxes.insert(i+1-M+M*N);
1193 }
1194
1195 if ( (!is_xmax[i]) && (!is_ymax[i]) && (!is_zmin[i]) )
1196 {
1197 local_boxes.insert(i+1+M-M*N);
1198 }
1199
1200 if ( (!is_xmax[i]) && (!is_ymax[i]) && (!is_zmax[i]) )
1201 {
1202 local_boxes.insert(i+1+M+M*N);
1203 }
1204
1205 // Now add the periodic boxes if any periodicity
1206 if (mIsPeriodicInX && ( is_xmin[i] || is_xmax[i]) )
1207 {
1208 // We are repeating the same steps on each z level, so easiest is
1209 // to make a vector of the z levels we want
1210 std::vector< int > z_i_offsets(1,0); // Add the current z box level
1211 if ( !is_zmin[i] )
1212 {
1213 z_i_offsets.push_back(-M*N); // Add the z box level below
1214 }
1215 if ( !is_zmax[i] )
1216 {
1217 z_i_offsets.push_back(M*N); // Add the z box level above
1218 }
1219
1220 // If we are on the left, add the nine on the right
1221 if ( is_xmin[i] )
1222 {
1223 // Loop over the z levels
1224 for ( std::vector<int>::iterator it = z_i_offsets.begin(); it != z_i_offsets.end(); it++ )
1225 {
1226 local_boxes.insert( i + (*it) + (M-1) ); // The right-most box on the same row
1227 // We also need to check for y boundaries
1228 if ( !is_ymin[i] )
1229 {
1230 local_boxes.insert( i + (*it) - 1 ); // The right-most box one row below
1231 }
1232 if ( !is_ymax[i] )
1233 {
1234 local_boxes.insert( i + (*it) + (2*M-1) ); // The right-most box one row above
1235 }
1236 }
1237 }
1238
1239 // If we are on the right, add the nine on the left
1240 else if ( is_xmax[i] )
1241 {
1242 // Loop over the z levels
1243 for ( std::vector<int>::iterator it = z_i_offsets.begin(); it != z_i_offsets.end(); it++ )
1244 {
1245 local_boxes.insert( i + (*it) - (M-1) ); // The left-most box on the same row
1246 // We also need to check for y boundaries
1247 if ( !is_ymin[i] )
1248 {
1249 local_boxes.insert( i + (*it) - (2*M-1) ); // The left-most box one row below
1250 }
1251 if ( !is_ymax[i] )
1252 {
1253 local_boxes.insert( i + (*it) + 1 ); // The left-most box one row below
1254 }
1255 }
1256 }
1257 }
1258
1259 if ( mIsPeriodicInY && (is_ymax[i] || is_ymin[i]) )
1260 {
1261 // We consider the current and upper z level and create a vector of the
1262 // opposite box indices that need to be added
1263 std::vector<unsigned> opp_box_i(0);
1264 if ( is_ymin[i] )
1265 {
1266 opp_box_i.push_back(i + (N-1)*M); // Current z level
1267 if ( !is_zmin[i] )
1268 {
1269 opp_box_i.push_back( i - M ); // z level below
1270 }
1271 if ( !is_zmax[i] )
1272 {
1273 opp_box_i.push_back(i + 2*M*N - M); // z level above
1274 }
1275 }
1276 else if ( is_ymax[i] )
1277 {
1278 opp_box_i.push_back( i - (N-1)*M ); // Current z level
1279 if ( !is_zmin[i] )
1280 {
1281 opp_box_i.push_back( i - 2*M*N + M ); // z level below
1282 }
1283 if ( !is_zmax[i] )
1284 {
1285 opp_box_i.push_back( i + M ); // z level above
1286 }
1287 }
1288
1289 // Now we add the different boxes, checking for left and right
1290 for ( std::vector<unsigned>::iterator it_opp_box = opp_box_i.begin(); it_opp_box != opp_box_i.end(); it_opp_box++ )
1291 {
1292 local_boxes.insert( *it_opp_box );
1293 if ( !is_xmin[i] )
1294 {
1295 local_boxes.insert( *it_opp_box - 1 );
1296 }
1297 if ( !is_xmax[i] )
1298 {
1299 local_boxes.insert( *it_opp_box + 1 );
1300 }
1301 }
1302 }
1303
1304 if ( mIsPeriodicInX && mIsPeriodicInY )
1305 {
1306 // Need to add the corners
1307 if ( is_xmin[i] && is_ymin[i] )
1308 {
1309 // Current z level
1310 local_boxes.insert(i+M*N-1);
1311 if ( !is_zmax[i] )
1312 {
1313 // Upper z level
1314 local_boxes.insert(i+2*M*N-1);
1315 }
1316 if ( !is_zmin[i] )
1317 {
1318 // Lower z level
1319 local_boxes.insert(i-1);
1320 }
1321 }
1322 if ( is_xmax[i] && is_ymin[i] )
1323 {
1324 local_boxes.insert(i + (N-2)*M + 1);
1325 if ( !is_zmax[i] )
1326 {
1327 local_boxes.insert(i + M*N + (N-2)*M + 1);
1328 }
1329 if ( !is_zmin[i] )
1330 {
1331 // Lower z level
1332 local_boxes.insert(i-2*M+1);
1333 }
1334 }
1335 if ( is_xmin[i] && is_ymax[i] )
1336 {
1337 local_boxes.insert(i + (N-2)*M - 1);
1338 if (!is_zmax[i])
1339 {
1340 // Upper z level
1341 local_boxes.insert(i - 2*M - 1);
1342 }
1343 if (!is_zmin[i])
1344 {
1345 // Lower z level
1346 local_boxes.insert(i-2*(N-1)*M-1);
1347 }
1348
1349 }
1350 if ( is_xmax[i] && is_ymax[i] )
1351 {
1352 if (!is_zmax[i])
1353 {
1354 // Upper z level
1355 local_boxes.insert(i + 1);
1356 }
1357 if (!is_zmin[i])
1358 {
1359 // Lower z level
1360 local_boxes.insert(i - 2*M*N + 1);
1361 }
1362
1363 }
1364 }
1365
1366 if (mIsPeriodicInZ && (is_zmin[i] || is_zmax[i]))
1367 {
1368 if ( is_zmin[i] )
1369 {
1370 // We need to add the top level
1371 unsigned above_box = i+(P-1)*M*N;
1372 local_boxes.insert(above_box);
1373 if (!is_xmin[i])
1374 {
1375 // Also add the boxes to the left and at the top
1376 local_boxes.insert(above_box-1);
1377 if (!is_ymax[i])
1378 {
1379 local_boxes.insert(above_box+M-1);
1380 }
1381
1382 if (!is_ymin[i])
1383 {
1384 local_boxes.insert(above_box-M-1);
1385 }
1386 }
1387 else if ( mIsPeriodicInX )
1388 {
1389 // Add x periodic box in top layer
1390 local_boxes.insert(above_box+M-1);
1391 if ( !is_ymin[i] )
1392 {
1393 local_boxes.insert(above_box+2*M-1);
1394 }
1395 else if ( mIsPeriodicInY )
1396 {
1397 // Add the xy periodic box in top layer if at y min or max
1398 local_boxes.insert(above_box+M*N-1);
1399 }
1400 if ( !is_ymax[i] )
1401 {
1402 local_boxes.insert(above_box+2*M-1);
1403 }
1404 else if ( mIsPeriodicInY )
1405 {
1406 // Add the xy periodic box in top layer if at y min or max
1407 local_boxes.insert(above_box-M*(N-2)-1);
1408 }
1409 }
1410
1411 if (!is_xmax[i])
1412 {
1413 // Also add the boxes to the left and at the top
1414 local_boxes.insert(above_box+1);
1415 if (!is_ymax[i])
1416 {
1417 local_boxes.insert(above_box+M+1);
1418 }
1419 if (!is_ymin[i])
1420 {
1421 local_boxes.insert(above_box-M+1);
1422 }
1423 }
1424 else if ( mIsPeriodicInX )
1425 {
1426 // Add x periodic box in top layer
1427 local_boxes.insert(above_box - M+1);
1428 if ( mIsPeriodicInY )
1429 {
1430 if ( is_ymin[i] )
1431 {
1432 // Add the xy periodic box in top layer if at y min or max
1433 local_boxes.insert(above_box+M*(N-2)+1);
1434 }
1435 else if ( is_ymax[i] )
1436 {
1437 // Add the xy periodic box in top layer if at y min or max
1438 local_boxes.insert(above_box-M*N+1);
1439 }
1440 }
1441 }
1442
1443
1444 if (!is_ymax[i])
1445 {
1446 local_boxes.insert(above_box+M);
1447 }
1448 else if ( mIsPeriodicInY )
1449 {
1450 local_boxes.insert(above_box-M*(N-1));
1451 if ( !is_xmin[i] )
1452 {
1453 local_boxes.insert(above_box-M*(N-1)-1);
1454 }
1455 if ( !is_xmax[i] )
1456 {
1457 local_boxes.insert(above_box-M*(N-1)+1);
1458 }
1459 }
1460 if (!is_ymin[i])
1461 {
1462 local_boxes.insert(above_box-M);
1463 }
1464 else if ( mIsPeriodicInY )
1465 {
1466 local_boxes.insert(above_box+M*(N-1));
1467 if ( !is_xmin[i] )
1468 {
1469 local_boxes.insert(above_box+M*(N-1)-1);
1470 }
1471 if ( !is_xmax[i] )
1472 {
1473 local_boxes.insert(above_box+M*(N-1)+1);
1474 }
1475 }
1476 }
1477 else if ( is_zmax[i] )
1478 {
1479 // We need to add the bottom level
1480 unsigned below_box = i-(P-1)*M*N;
1481 local_boxes.insert(below_box);
1482 if (!is_xmin[i])
1483 {
1484 // Also add the boxes to the left and at the top
1485 local_boxes.insert(below_box-1);
1486 if (!is_ymax[i])
1487 {
1488 local_boxes.insert(below_box+M-1);
1489 }
1490
1491 if (!is_ymin[i])
1492 {
1493 local_boxes.insert(below_box-M-1);
1494 }
1495 }
1496 else if ( mIsPeriodicInX )
1497 {
1498 // Add x periodic box in top layer
1499 local_boxes.insert(below_box+M-1);
1500 if ( mIsPeriodicInY )
1501 {
1502 if ( is_ymin[i] )
1503 {
1504 // Add the xy periodic box in top layer if at y min or max
1505 local_boxes.insert(below_box+M*N-1);
1506 }
1507 else if ( is_ymax[i] )
1508 {
1509 // Add the xy periodic box in top layer if at y min or max
1510 local_boxes.insert(below_box-M*(N-2)-1);
1511 }
1512 }
1513 }
1514
1515 if (!is_xmax[i])
1516 {
1517 // Also add the boxes to the left and at the top
1518 local_boxes.insert(below_box+1);
1519 if (!is_ymax[i])
1520 {
1521 local_boxes.insert(below_box+M+1);
1522 }
1523 if (!is_ymin[i])
1524 {
1525 local_boxes.insert(below_box-M+1);
1526 }
1527 }
1528 else if ( mIsPeriodicInX )
1529 {
1530 // Add x periodic box in top layer
1531 local_boxes.insert(below_box - M+1);
1532 if ( mIsPeriodicInY )
1533 {
1534 if ( is_ymin[i] )
1535 {
1536 // Add the xy periodic box in top layer if at y min or max
1537 local_boxes.insert(below_box+M*(N-2)+1);
1538 }
1539 else if ( is_ymax[i] )
1540 {
1541 // Add the xy periodic box in top layer if at y min or max
1542 local_boxes.insert(below_box-M*N+1);
1543 }
1544 }
1545 }
1546
1547
1548 if (!is_ymax[i])
1549 {
1550 local_boxes.insert(below_box+M);
1551 }
1552 else if ( mIsPeriodicInY )
1553 {
1554 local_boxes.insert(below_box-M*(N-1));
1555 if ( !is_xmin[i] )
1556 {
1557 local_boxes.insert(below_box-M*(N-1)+1);
1558 }
1559 if (!is_xmax[i])
1560 {
1561 local_boxes.insert(below_box-M*(N-1)-1);
1562 }
1563 }
1564 if (!is_ymin[i])
1565 {
1566 local_boxes.insert(below_box-M);
1567 }
1568 else if ( mIsPeriodicInY )
1569 {
1570 local_boxes.insert(below_box+M*(N-1));
1571 if ( !is_xmin[i] )
1572 {
1573 local_boxes.insert(below_box+M*(N-1)+1);
1574 }
1575 if (!is_xmax[i])
1576 {
1577 local_boxes.insert(below_box+M*(N-1)-1);
1578 }
1579 }
1580 }
1581 }
1582
1583 mLocalBoxes.push_back(local_boxes);
1584 }
1585 break;
1586 }
1587 default:
1589 }
1590}
1591
1592template<unsigned DIM>
1593std::set<unsigned>& DistributedBoxCollection<DIM>::rGetLocalBoxes(unsigned boxIndex)
1594{
1595 // Make sure the box is locally owned
1596 assert(!(boxIndex < mMinBoxIndex) && !(mMaxBoxIndex<boxIndex));
1597 return mLocalBoxes[boxIndex-mMinBoxIndex];
1598}
1599
1600template<unsigned DIM>
1602{
1603 unsigned index = CalculateContainingBox(pNode);
1604
1605 return IsBoxOwned(index);
1606}
1607
1608template<unsigned DIM>
1609bool DistributedBoxCollection<DIM>::IsOwned(c_vector<double, DIM>& location)
1610{
1611 unsigned index = CalculateContainingBox(location);
1612
1613 return IsBoxOwned(index);
1614}
1615
1616template<unsigned DIM>
1618{
1619 unsigned box_index = CalculateContainingBox(pNode);
1620 unsigned containing_process = PetscTools::GetMyRank();
1621
1622 if (box_index > mMaxBoxIndex)
1623 {
1624 containing_process++;
1625 }
1626 else if (box_index < mMinBoxIndex)
1627 {
1628 containing_process--;
1629 }
1630
1631 // Need a special case for periodicity
1632 if ( mIsPeriodicAcrossProcs )
1633 {
1634 if (PetscTools::AmMaster() && box_index > ((mNumBoxesEachDirection[DIM-1]-1)*mNumBoxesInAFace-1))
1635 {
1636 // It needs to move to the top process
1637 containing_process = PetscTools::GetNumProcs()-1;
1638 }
1639 else if (PetscTools::AmTopMost() && box_index < mNumBoxesInAFace)
1640 {
1641 // It needs to move to the bottom process
1642 containing_process = 0;
1643 }
1644 }
1645
1646 return containing_process;
1647}
1648
1649template<unsigned DIM>
1651{
1652 return mHaloNodesRight;
1653}
1654
1655template<unsigned DIM>
1657{
1658 return mHaloNodesLeft;
1659}
1660
1661template<unsigned DIM>
1663{
1664 mCalculateNodeNeighbours = calculateNodeNeighbours;
1665}
1666
1667template<unsigned DIM>
1668void DistributedBoxCollection<DIM>::CalculateNodePairs(std::vector<Node<DIM>*>& rNodes, std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
1669{
1670 rNodePairs.clear();
1671
1672 // Create an empty neighbours set for each node
1673 for (unsigned i=0; i<rNodes.size(); i++)
1674 {
1675 // Get the box containing this node as only nodes on this process have NodeAttributes
1676 // and therefore Neighbours setup.
1677 unsigned box_index = CalculateContainingBox(rNodes[i]);
1678
1679 if (IsBoxOwned(box_index))
1680 {
1681 rNodes[i]->ClearNeighbours();
1682 }
1683 }
1684
1685 for (unsigned box_index=mMinBoxIndex; box_index<=mMaxBoxIndex; box_index++)
1686 {
1687 AddPairsFromBox(box_index, rNodePairs);
1688 }
1689
1690 if (mCalculateNodeNeighbours)
1691 {
1692 for (unsigned i = 0; i < rNodes.size(); i++)
1693 {
1694 // Get the box containing this node as only nodes on this process have NodeAttributes
1695 // and therefore Neighbours setup.
1696 unsigned box_index = CalculateContainingBox(rNodes[i]);
1697
1698 if (IsBoxOwned(box_index))
1699 {
1700 rNodes[i]->RemoveDuplicateNeighbours();
1701 }
1702 }
1703 }
1704}
1705
1706template<unsigned DIM>
1707void DistributedBoxCollection<DIM>::CalculateInteriorNodePairs(std::vector<Node<DIM>*>& rNodes, std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
1708{
1709 rNodePairs.clear();
1710
1711 // Create an empty neighbours set for each node
1712 for (unsigned i=0; i<rNodes.size(); i++)
1713 {
1714 // Get the box containing this node as only nodes on this process have NodeAttributes
1715 // and therefore Neighbours setup.
1716 unsigned box_index = CalculateContainingBox(rNodes[i]);
1717
1718 if (IsBoxOwned(box_index))
1719 {
1720 rNodes[i]->ClearNeighbours();
1721 rNodes[i]->SetNeighboursSetUp(false);
1722 }
1723 }
1724
1725 for (unsigned box_index=mMinBoxIndex; box_index<=mMaxBoxIndex; box_index++)
1726 {
1727 if (IsInteriorBox(box_index))
1728 {
1729 AddPairsFromBox(box_index, rNodePairs);
1730 }
1731 }
1732
1733 if (mCalculateNodeNeighbours)
1734 {
1735 for (unsigned i = 0; i < rNodes.size(); i++)
1736 {
1737 // Get the box containing this node as only nodes on this process have NodeAttributes
1738 // and therefore Neighbours setup.
1739 unsigned box_index = CalculateContainingBox(rNodes[i]);
1740
1741 if (IsBoxOwned(box_index))
1742 {
1743 rNodes[i]->RemoveDuplicateNeighbours();
1744 rNodes[i]->SetNeighboursSetUp(true);
1745 }
1746 }
1747 }
1748}
1749
1750template<unsigned DIM>
1751void DistributedBoxCollection<DIM>::CalculateBoundaryNodePairs(std::vector<Node<DIM>*>& rNodes, std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
1752{
1753 for (unsigned box_index=mMinBoxIndex; box_index<=mMaxBoxIndex; box_index++)
1754 {
1755 if (!IsInteriorBox(box_index))
1756 {
1757 AddPairsFromBox(box_index, rNodePairs);
1758 }
1759 }
1760
1761 if (mCalculateNodeNeighbours)
1762 {
1763 for (unsigned i = 0; i < rNodes.size(); i++)
1764 {
1765 // Get the box containing this node as only nodes on this process have NodeAttributes
1766 // and therefore Neighbours setup.
1767 unsigned box_index = CalculateContainingBox(rNodes[i]);
1768
1769 if (IsBoxOwned(box_index))
1770 {
1771 rNodes[i]->RemoveDuplicateNeighbours();
1772 rNodes[i]->SetNeighboursSetUp(true);
1773 }
1774 }
1775 }
1776}
1777
1778template<unsigned DIM>
1780 std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
1781{
1782 // Get the box
1783 Box<DIM>& r_box = rGetBox(boxIndex);
1784
1785 // Get the set of nodes in this box
1786 const std::set< Node<DIM>* >& r_contained_nodes = r_box.rGetNodesContained();
1787
1788 // Get the local boxes to this box
1789 const std::set<unsigned>& local_boxes_indices = rGetLocalBoxes(boxIndex);
1790
1791 // Loop over all the local boxes
1792 for (std::set<unsigned>::iterator box_iter = local_boxes_indices.begin();
1793 box_iter != local_boxes_indices.end();
1794 box_iter++)
1795 {
1796 Box<DIM>* p_neighbour_box;
1797
1798 // Establish whether box is locally owned or halo.
1799 if (IsBoxOwned(*box_iter))
1800 {
1801 p_neighbour_box = &mBoxes[*box_iter - mMinBoxIndex];
1802 }
1803 else // Assume it is a halo.
1804 {
1805 p_neighbour_box = &mHaloBoxes[mHaloBoxesMapping[*box_iter]];
1806 }
1807 assert(p_neighbour_box);
1808
1809 // Get the set of nodes contained in this box
1810 std::set< Node<DIM>* >& r_contained_neighbour_nodes = p_neighbour_box->rGetNodesContained();
1811
1812 // Loop over these nodes
1813 for (typename std::set<Node<DIM>*>::iterator neighbour_node_iter = r_contained_neighbour_nodes.begin();
1814 neighbour_node_iter != r_contained_neighbour_nodes.end();
1815 ++neighbour_node_iter)
1816 {
1817 // Get the index of the other node
1818 unsigned other_node_index = (*neighbour_node_iter)->GetIndex();
1819
1820 // Loop over nodes in this box
1821 for (typename std::set<Node<DIM>*>::iterator node_iter = r_contained_nodes.begin();
1822 node_iter != r_contained_nodes.end();
1823 ++node_iter)
1824 {
1825 unsigned node_index = (*node_iter)->GetIndex();
1826
1827 // If we're in the same box, then take care not to store the node pair twice
1828 if (*box_iter != boxIndex || other_node_index > node_index)
1829 {
1830 rNodePairs.push_back(std::pair<Node<DIM>*, Node<DIM>*>((*node_iter), (*neighbour_node_iter)));
1831 if (mCalculateNodeNeighbours)
1832 {
1833 (*node_iter)->AddNeighbour(other_node_index);
1834 (*neighbour_node_iter)->AddNeighbour(node_index);
1835 }
1836 }
1837
1838 }
1839 }
1840 }
1841}
1842
1843template<unsigned DIM>
1845{
1846 std::vector<int> cell_numbers(mpDistributedBoxStackFactory->GetHigh() - mpDistributedBoxStackFactory->GetLow(), 0);
1847
1848 for (unsigned global_index=mMinBoxIndex; global_index<=mMaxBoxIndex; global_index++)
1849 {
1850 c_vector<unsigned, DIM> coords = CalculateGridIndices(global_index);
1851 unsigned location_in_vector = coords[DIM-1] - mpDistributedBoxStackFactory->GetLow();
1852 unsigned local_index = global_index - mMinBoxIndex;
1853 cell_numbers[location_in_vector] += mBoxes[local_index].rGetNodesContained().size();
1854 }
1855
1856 return cell_numbers;
1857}
1858
1860
1861template class DistributedBoxCollection<1>;
1862template class DistributedBoxCollection<2>;
1863template class DistributedBoxCollection<3>;
#define EXCEPTION(message)
#define NEVER_REACHED
Definition Box.hpp:51
std::set< Node< DIM > * > & rGetNodesContained()
Definition Box.cpp:56
void AddPairsFromBox(unsigned boxIndex, std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs)
std::vector< Box< DIM > > mBoxes
std::set< unsigned > & rGetLocalBoxes(unsigned boxIndex)
DistributedBoxCollection(double boxWidth, c_vector< double, 2 *DIM > domainSize, bool isPeriodicInX=false, bool isPeriodicInY=false, bool isPeriodicInZ=false, int localRows=PETSC_DECIDE)
Box< DIM > & rGetHaloBox(unsigned boxIndex)
c_vector< double, 2 *DIM > rGetDomainSize() const
int LoadBalance(std::vector< int > localDistribution)
void CalculateNodePairs(std::vector< Node< DIM > * > &rNodes, std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs)
void CalculateInteriorNodePairs(std::vector< Node< DIM > * > &rNodes, std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs)
unsigned GetProcessOwningNode(Node< DIM > *pNode)
void SetCalculateNodeNeighbours(bool calculateNodeNeighbours)
bool IsInteriorBox(unsigned globalIndex)
unsigned CalculateContainingBox(Node< DIM > *pNode)
std::vector< unsigned > & rGetHaloNodesLeft()
std::vector< int > CalculateNumberOfNodesInEachStrip()
bool IsBoxOwned(unsigned globalIndex)
std::vector< unsigned > & rGetHaloNodesRight()
bool IsHaloBox(unsigned globalIndex)
DistributedVectorFactory * mpDistributedBoxStackFactory
c_vector< bool, DIM > GetIsPeriodicAllDims() const
c_vector< double, 2 *DIM > mDomainSize
void CalculateBoundaryNodePairs(std::vector< Node< DIM > * > &rNodes, std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs)
c_vector< unsigned, DIM > mNumBoxesEachDirection
unsigned CalculateGlobalIndex(c_vector< unsigned, DIM > gridIndices)
c_vector< unsigned, DIM > CalculateGridIndices(unsigned globalIndex)
Box< DIM > & rGetBox(unsigned boxIndex)
Definition Node.hpp:59
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition Node.cpp:139
static bool AmMaster()
static bool AmTopMost()
static bool IsParallel()
static bool IsSequential()
static unsigned GetMyRank()
static unsigned GetNumProcs()