BoxCollection.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 #include "BoxCollection.hpp"
00029 #include "PetscTools.hpp"
00030 
00032 // Box methods
00034 
00035 
00036 template<unsigned DIM>
00037 Box<DIM>::Box(c_vector<double, 2*DIM>& rMinAndMaxValues)
00038 {
00039     mMinAndMaxValues = rMinAndMaxValues;
00040 }
00041 
00042 template<unsigned DIM>
00043 c_vector<double, 2*DIM>& Box<DIM>::rGetMinAndMaxValues()
00044 {
00045     return mMinAndMaxValues;
00046 }
00047 
00048 template<unsigned DIM>
00049 void Box<DIM>::AddNode(Node<DIM>* pNode)
00050 {
00051     mNodesContained.insert(pNode);
00052 }
00053 
00054 template<unsigned DIM>
00055 void Box<DIM>::RemoveNode(Node<DIM>* pNode)
00056 {
00057     mNodesContained.erase(pNode);
00058 }
00059 
00060 template<unsigned DIM>
00061 std::set< Node<DIM>* >& Box<DIM>::rGetNodesContained()
00062 {
00063     return mNodesContained;
00064 }
00065 
00066 template<unsigned DIM>
00067 void Box<DIM>::AddElement(Element<DIM,DIM>* pElement)
00068 {
00069     mElementsContained.insert(pElement);
00070 }
00071 
00072 template<unsigned DIM>
00073 std::set< Element<DIM,DIM>* >& Box<DIM>::rGetElementsContained()
00074 {
00075     return mElementsContained;
00076 }
00077 
00078 
00080 // BoxCollection methods
00082 
00083 
00084 template<unsigned DIM>
00085 BoxCollection<DIM>::BoxCollection(double boxWidth, c_vector<double, 2*DIM> domainSize)
00086     : mDomainSize(domainSize),
00087       mBoxWidth(boxWidth)
00088 {
00089     /*
00090      * Start by calculating the number of boxes in each direction and total number of boxes.
00091      * Also create a helper vector of coefficients, whose first entry is 1 and whose i-th
00092      * entry (for i>1) is the i-th partial product of the vector mNumBoxesEachDirection. This
00093      * vector of coefficients will be used in the next code block to compute how many boxes
00094      * along in each dimension each box, given its index.
00095      */
00096     unsigned num_boxes = 1;
00097     std::vector<unsigned> coefficients;
00098     coefficients.push_back(1);
00099 
00100     for (unsigned i=0; i<DIM; i++)
00101     {
00102         mNumBoxesEachDirection(i) = (unsigned)((domainSize(2*i+1) - domainSize(2*i))/boxWidth) + 1;
00103         num_boxes *= mNumBoxesEachDirection(i);
00104         coefficients.push_back(coefficients[i]*mNumBoxesEachDirection(i));
00105     }
00106 
00107     for (unsigned box_index=0; box_index<num_boxes; box_index++)
00108     {
00109         /*
00110          * The code block below computes how many boxes along in each dimension the
00111          * current box is and stores this information in the second, ..., (DIM+1)th
00112          * entries of the vector current_box_indices. The first entry of
00113          * current_box_indices is zero to ensure that the for loop works correctly.
00114          *
00115          * Our convention is that in 3D the index of each box, box_index, is related
00116          * to its indices (i,j,k) by
00117          *
00118          * box_index = i + mNumBoxesEachDirection(0)*j
00119          *               + mNumBoxesEachDirection(0)*mNumBoxesEachDirection(1)*k,
00120          *
00121          * while in 2D, box_index is related to (i,j) by
00122          *
00123          * box_index = i + mNumBoxesEachDirection(0)*j
00124          *
00125          * and in 1D we simply have box_index = i.
00126          */
00127         c_vector<unsigned, DIM+1> current_box_indices;
00128         current_box_indices[0] = 0;
00129 
00130         for (unsigned i=0; i<DIM; i++)
00131         {
00132             unsigned temp = 0;
00133             for (unsigned j=1; j<i; j++)
00134             {
00135                 temp += coefficients[j]*current_box_indices[j-1];
00136             }
00137             current_box_indices[i+1] = (box_index%coefficients[i+1] - temp)/coefficients[i];
00138         }
00139 
00140         /*
00141          * We now use the information stores in current_box_indices to construct the
00142          * Box, which we add to mBoxes.
00143          */
00144         c_vector<double, 2*DIM> box_coords;
00145         for (unsigned l=0; l<DIM; l++)
00146         {
00147             box_coords(2*l) = domainSize(2*l) + current_box_indices(l+1)*boxWidth;
00148             box_coords(2*l+1) = domainSize(2*l) + (current_box_indices(l+1)+1)*boxWidth;
00149         }
00150 
00151         Box<DIM> new_box(box_coords);
00152         mBoxes.push_back(new_box);
00153     }
00154 
00155     // Check that we have the correct number of boxes
00156     assert(num_boxes == mBoxes.size());
00157 }
00158 
00159 template<unsigned DIM>
00160 unsigned BoxCollection<DIM>::CalculateContainingBox(Node<DIM>* pNode)
00161 {
00162     // Get the location of the node
00163     c_vector<double, DIM> location = pNode->rGetLocation();
00164     return CalculateContainingBox(location);
00165 }
00166 
00167 
00168 template<unsigned DIM>
00169 unsigned BoxCollection<DIM>::CalculateContainingBox(c_vector<double, DIM>& rLocation)
00170 {
00171     // The node must lie inside the boundary of the box collection
00172     for (unsigned i=0; i<DIM; i++)
00173     {
00174         if( (rLocation[i] < mDomainSize(2*i)) || (rLocation[i] > mDomainSize(2*i+1)) )
00175         {
00176             EXCEPTION("The point provided in outside all of the boxes");
00177         }
00178     }
00179 
00180     // Compute the containing box index in each dimension
00181     c_vector<unsigned, DIM> containing_box_indices;
00182     for (unsigned i=0; i<DIM; i++)
00183     {
00184         containing_box_indices[i] = (unsigned) floor((rLocation[i] - mDomainSize(2*i))/mBoxWidth);
00185     }
00186 
00187     // Use these to compute the index of the containing box
00188     unsigned containing_box_index = 0;
00189     for (unsigned i=0; i<DIM; i++)
00190     {
00191         unsigned temp = 1;
00192         for (unsigned j=0; j<i; j++)
00193         {
00194             temp *= mNumBoxesEachDirection(j);
00195         }
00196         containing_box_index += temp*containing_box_indices[i];
00197     }
00198 
00199     // This index must be less than the number of boxes
00200     assert(containing_box_index < mBoxes.size());
00201 
00202     return containing_box_index;
00203 }
00204 
00205 template<unsigned DIM>
00206 Box<DIM>& BoxCollection<DIM>::rGetBox(unsigned boxIndex)
00207 {
00208     assert(boxIndex < mBoxes.size());
00209     return mBoxes[boxIndex];
00210 }
00211 
00212 template<unsigned DIM>
00213 unsigned BoxCollection<DIM>::GetNumBoxes()
00214 {
00215     return mBoxes.size();
00216 }
00217 
00218 template<unsigned DIM>
00219 void BoxCollection<DIM>::SetupLocalBoxesHalfOnly()
00220 {
00221     switch (DIM)
00222     {
00223         case 1:
00224         {
00225             // We only need to look for neighbours in the current and successive boxes
00226             mLocalBoxes.clear();
00227             for (unsigned box_index=0; box_index<mBoxes.size(); box_index++)
00228             {
00229                 std::set<unsigned> local_boxes;
00230 
00231                 // Insert the current box
00232                 local_boxes.insert(box_index);
00233 
00234                 // If we're not at the right-most box, then insert the box to the right
00235                 if (box_index < mNumBoxesEachDirection(0)-1)
00236                 {
00237                     local_boxes.insert(box_index+1);
00238                 }
00239                 mLocalBoxes.push_back(local_boxes);
00240             }
00241             break;
00242         }
00243         case 2:
00244         {
00245             // We only need to look for neighbours in the current box and half the neighbouring boxes
00246             mLocalBoxes.clear();
00247             for (unsigned box_index=0; box_index<mBoxes.size(); box_index++)
00248             {
00249                 std::set<unsigned> local_boxes;
00250 
00251                 // Insert the current box
00252                 local_boxes.insert(box_index);
00253 
00254                 // If we're not on the top-most row, then insert the box above
00255                 if (box_index < mBoxes.size() - mNumBoxesEachDirection(0))
00256                 {
00257                     local_boxes.insert(box_index + mNumBoxesEachDirection(0));
00258 
00259                     // If we're also not on the left-most column, then insert the box above-left
00260                     if (box_index % mNumBoxesEachDirection(0) != 0)
00261                     {
00262                         local_boxes.insert(box_index + mNumBoxesEachDirection(0) - 1);
00263                     }
00264                 }
00265                 // If we're not on the right-most column, then insert the box to the right
00266                 if (box_index % mNumBoxesEachDirection(0) != mNumBoxesEachDirection(0)-1)
00267                 {
00268                     local_boxes.insert(box_index + 1);
00269 
00270                     // If we're also not on the top-most row, then insert the box above-right
00271                     if (box_index < mBoxes.size() - mNumBoxesEachDirection(0))
00272                     {
00273                         local_boxes.insert(box_index + mNumBoxesEachDirection(0) + 1);
00274                     }
00275                 }
00276 
00277                 mLocalBoxes.push_back(local_boxes);
00278             }
00279             break;
00280         }
00281         case 3:
00282         {
00283             // We only need to look for neighbours in the current box and half the neighbouring boxes
00284             mLocalBoxes.clear();
00285             unsigned num_boxes_xy = mNumBoxesEachDirection(0)*mNumBoxesEachDirection(1);
00286 
00287             for (unsigned box_index=0; box_index<mBoxes.size(); box_index++)
00288             {
00289                 std::set<unsigned> local_boxes;
00290 
00291                 // Insert the current box
00292                 local_boxes.insert(box_index);
00293 
00294                 // If we're not on the far face (y max), then insert the far box
00295                 if (box_index % num_boxes_xy < num_boxes_xy - mNumBoxesEachDirection(0))
00296                 {
00297                     local_boxes.insert(box_index + mNumBoxesEachDirection(0));
00298 
00299                     // If we're also not on the left face (x min), then insert the box to the left
00300                     if (box_index % mNumBoxesEachDirection(0) != 0)
00301                     {
00302                         local_boxes.insert(box_index + mNumBoxesEachDirection(0) - 1);
00303                     }
00304                 }
00305                 // If we're not on the right face (x max), then insert the box to the right
00306                 if (box_index % mNumBoxesEachDirection(0) != mNumBoxesEachDirection(0)-1)
00307                 {
00308                     local_boxes.insert(box_index + 1);
00309 
00310                     // If we're also not on the far face (y max) row, then insert the box to the far-right
00311                     if (box_index % num_boxes_xy < num_boxes_xy - mNumBoxesEachDirection(0))
00312                     {
00313                         local_boxes.insert(box_index + mNumBoxesEachDirection(0) + 1);
00314                     }
00315                 }
00316                 // If we're not on the top face (z max), then insert the box above
00317                 if (box_index < mBoxes.size() - num_boxes_xy)
00318                 {
00319                     local_boxes.insert(box_index + num_boxes_xy);
00320 
00321                     // If we're also not on the far face (y max), then insert the above-far box
00322                     if (box_index % num_boxes_xy < num_boxes_xy - mNumBoxesEachDirection(0))
00323                     {
00324                         local_boxes.insert(box_index + num_boxes_xy + mNumBoxesEachDirection(0));
00325 
00326                         // If we're also not on the left face (x min), then insert the box to the above-left
00327                         if (box_index % mNumBoxesEachDirection(0) != 0)
00328                         {
00329                             local_boxes.insert(box_index + num_boxes_xy + mNumBoxesEachDirection(0) - 1);
00330                         }
00331                     }
00332                     // If we're also not on the right face (x max), then insert the box to the above-right
00333                     if (box_index % mNumBoxesEachDirection(0) != mNumBoxesEachDirection(0)-1)
00334                     {
00335                         local_boxes.insert(box_index + num_boxes_xy + 1);
00336 
00337                         // If we're also not on the far face (y max) row, then insert the box to the above-far-right
00338                         if (box_index % num_boxes_xy < num_boxes_xy - mNumBoxesEachDirection(0))
00339                         {
00340                             local_boxes.insert(box_index + num_boxes_xy + mNumBoxesEachDirection(0) + 1);
00341                         }
00342                     }
00343                 }
00344                 // If we're not on the bottom face (z min), then insert the box above
00345                 if (box_index >= num_boxes_xy)
00346                 {
00347                     local_boxes.insert(box_index - num_boxes_xy);
00348 
00349                     // If we're also not on the far face (y max), then insert the below-far box
00350                     if (box_index % num_boxes_xy < num_boxes_xy - mNumBoxesEachDirection(0))
00351                     {
00352                         local_boxes.insert(box_index - num_boxes_xy + mNumBoxesEachDirection(0));
00353 
00354                         // If we're also not on the left face (x min), then insert the box to the below-left
00355                         if (box_index % mNumBoxesEachDirection(0) != 0)
00356                         {
00357                             local_boxes.insert(box_index - num_boxes_xy + mNumBoxesEachDirection(0) - 1);
00358                         }
00359                     }
00360                     // If we're also not on the right face (x max), then insert the box to the below-right
00361                     if (box_index % mNumBoxesEachDirection(0) != mNumBoxesEachDirection(0)-1)
00362                     {
00363                         local_boxes.insert(box_index - num_boxes_xy + 1);
00364 
00365                         // If we're also not on the far face (y max) row, then insert the box to the below-far-right
00366                         if (box_index % num_boxes_xy < num_boxes_xy - mNumBoxesEachDirection(0))
00367                         {
00368                             local_boxes.insert(box_index - num_boxes_xy + mNumBoxesEachDirection(0) + 1);
00369                         }
00370                     }
00371                 }
00372                 mLocalBoxes.push_back(local_boxes);
00373             }
00374             break;
00375         }
00376         default:
00377             NEVER_REACHED;
00378     }
00379 }
00380 
00381 
00382 
00383 template<unsigned DIM>
00384 void BoxCollection<DIM>::SetupAllLocalBoxes()
00385 {
00386     switch (DIM)
00387     {
00388         case 1:
00389         {
00390             for (unsigned i=0; i<mBoxes.size(); i++)
00391             {
00392                 std::set<unsigned> local_boxes;
00393 
00394                 local_boxes.insert(i);
00395              
00396                 // add the two neighbours   
00397                 if(i!=0)
00398                 {
00399                     local_boxes.insert(i-1);
00400                 }
00401                 if(i+1 != mNumBoxesEachDirection(0))
00402                 {
00403                     local_boxes.insert(i+1);
00404                 }
00405                 
00406                 mLocalBoxes.push_back(local_boxes);
00407             }
00408             break;
00409         }
00410         case 2:
00411         {
00412             mLocalBoxes.clear();
00413             
00414             unsigned M = mNumBoxesEachDirection(0);
00415             unsigned N = mNumBoxesEachDirection(1);
00416 
00417             std::vector<bool> is_xmin(N*M); // far left
00418             std::vector<bool> is_xmax(N*M); // far right
00419             std::vector<bool> is_ymin(N*M); // bottom
00420             std::vector<bool> is_ymax(N*M); // top
00421             
00422             for(unsigned i=0; i<M*N; i++)
00423             {
00424                 is_xmin[i] = (i%M==0);
00425                 is_xmax[i] = ((i+1)%M==0);
00426                 is_ymin[i] = (i%(M*N)<M);
00427                 is_ymax[i] = (i%(M*N)>=(N-1)*M);
00428             }
00429 
00430             for (unsigned i=0; i<mBoxes.size(); i++)
00431             {
00432                 std::set<unsigned> local_boxes;
00433 
00434                 local_boxes.insert(i);
00435                 
00436                 // add the box to the left
00437                 if(!is_xmin[i])
00438                 {
00439                     local_boxes.insert(i-1);
00440                 }
00441         
00442                 // add the box to the right
00443                 if(!is_xmax[i])
00444                 {
00445                     local_boxes.insert(i+1);
00446                 }
00447                 
00448                 // add the one below
00449                 if(!is_ymin[i])
00450                 {
00451                     local_boxes.insert(i-M);
00452                 }
00453                 
00454                 // add the one above
00455                 if(!is_ymax[i])
00456                 {
00457                     local_boxes.insert(i+M);
00458                 }
00459                 
00460                 // add the four corner boxes
00461         
00462                 if( (!is_xmin[i]) && (!is_ymin[i]) )
00463                 {
00464                     local_boxes.insert(i-1-M);
00465                 }
00466         
00467                 if( (!is_xmin[i]) && (!is_ymax[i]) )
00468                 {
00469                     local_boxes.insert(i-1+M);
00470                 }
00471         
00472                 if( (!is_xmax[i]) && (!is_ymin[i]) )
00473                 {
00474                     local_boxes.insert(i+1-M);
00475                 }
00476         
00477                 if( (!is_xmax[i]) && (!is_ymax[i]) )
00478                 {
00479                     local_boxes.insert(i+1+M);
00480                 }
00481 
00482                 mLocalBoxes.push_back(local_boxes);
00483             }
00484             break;
00485         }
00486         case 3:
00487         {
00488             mLocalBoxes.clear();
00489             
00490             unsigned M = mNumBoxesEachDirection(0);
00491             unsigned N = mNumBoxesEachDirection(1);
00492             unsigned P = mNumBoxesEachDirection(2);
00493 
00494             std::vector<bool> is_xmin(N*M*P); // far left
00495             std::vector<bool> is_xmax(N*M*P); // far right
00496             std::vector<bool> is_ymin(N*M*P); // nearest
00497             std::vector<bool> is_ymax(N*M*P); // furthest
00498             std::vector<bool> is_zmin(N*M*P); // bottom layer
00499             std::vector<bool> is_zmax(N*M*P); // top layer
00500     
00501             for(unsigned i=0; i<M*N*P; i++)
00502             {
00503                 is_xmin[i] = (i%M==0);
00504                 is_xmax[i] = ((i+1)%M==0);
00505                 is_ymin[i] = (i%(M*N)<M);
00506                 is_ymax[i] = (i%(M*N)>=(N-1)*M);
00507                 is_zmin[i] = (i<M*N);
00508                 is_zmax[i] = (i>=M*N*(P-1));
00509             }
00510 
00511             for (unsigned i=0; i<mBoxes.size(); i++)
00512             {
00513                 std::set<unsigned> local_boxes;
00514 
00515                 // add itself as a local box
00516                 local_boxes.insert(i);
00517                 
00518                 // now add all 26 other neighbours.....
00519                 
00520                 // add the box left
00521                 if(!is_xmin[i])
00522                 {
00523                     local_boxes.insert(i-1);
00524                     
00525                     // plus some others towards the left
00526                     if(!is_ymin[i])
00527                     {
00528                         local_boxes.insert(i-1-M);
00529                     }
00530         
00531                     if(!is_ymax[i])
00532                     {
00533                         local_boxes.insert(i-1+M);
00534                     }
00535         
00536                     if(!is_zmin[i])
00537                     {
00538                         local_boxes.insert(i-1-M*N);
00539                     }
00540         
00541                     if(!is_zmax[i])
00542                     {
00543                         local_boxes.insert(i-1+M*N);
00544                     }
00545                 }
00546         
00547                 // add the box to the right
00548                 if(!is_xmax[i])
00549                 {
00550                     local_boxes.insert(i+1);
00551 
00552                     // plus some others towards the right
00553                     if(!is_ymin[i])
00554                     {
00555                         local_boxes.insert(i+1-M);
00556                     }
00557         
00558                     if(!is_ymax[i])
00559                     {
00560                         local_boxes.insert(i+1+M);
00561                     }
00562         
00563                     if(!is_zmin[i])
00564                     {
00565                         local_boxes.insert(i+1-M*N);
00566                     }
00567         
00568                     if(!is_zmax[i])
00569                     {
00570                         local_boxes.insert(i+1+M*N);
00571                     }
00572                 }
00573                 
00574                 // add the boxes next along the y axis
00575                 if(!is_ymin[i])
00576                 {
00577                     local_boxes.insert(i-M);
00578                     
00579                     // and more in this plane
00580                     if(!is_zmin[i])
00581                     {
00582                         local_boxes.insert(i-M-M*N);
00583                     }
00584         
00585                     if(!is_zmax[i])
00586                     {
00587                         local_boxes.insert(i-M+M*N);
00588                     }
00589                 }
00590                 
00591                 // add the boxes next along the y axis
00592                 if(!is_ymax[i])
00593                 {
00594                     local_boxes.insert(i+M);
00595         
00596                     // and more in this plane
00597                     if(!is_zmin[i])
00598                     {
00599                         local_boxes.insert(i+M-M*N);
00600                     }
00601         
00602                     if(!is_zmax[i])
00603                     {
00604                         local_boxes.insert(i+M+M*N);
00605                     }
00606                 }
00607                 
00608                 // add the box directly above
00609                 if(!is_zmin[i])
00610                 {
00611                     local_boxes.insert(i-N*M);
00612                 }            
00613                 
00614                 // add the box directly below
00615                 if(!is_zmax[i])
00616                 {
00617                     local_boxes.insert(i+N*M);
00618                 }
00619         
00620                 // finally, the 8 corners are left
00621         
00622                 if( (!is_xmin[i]) && (!is_ymin[i]) && (!is_zmin[i]) )
00623                 {
00624                     local_boxes.insert(i-1-M-M*N);
00625                 }
00626         
00627                 if( (!is_xmin[i]) && (!is_ymin[i]) && (!is_zmax[i]) )
00628                 {
00629                     local_boxes.insert(i-1-M+M*N);
00630                 }
00631         
00632                 if( (!is_xmin[i]) && (!is_ymax[i]) && (!is_zmin[i]) )
00633                 {
00634                     local_boxes.insert(i-1+M-M*N);
00635                 }
00636         
00637                 if( (!is_xmin[i]) && (!is_ymax[i]) && (!is_zmax[i]) )
00638                 {
00639                     local_boxes.insert(i-1+M+M*N);
00640                 }
00641         
00642                 if( (!is_xmax[i]) && (!is_ymin[i]) && (!is_zmin[i]) )
00643                 {
00644                     local_boxes.insert(i+1-M-M*N);
00645                 }
00646         
00647                 if( (!is_xmax[i]) && (!is_ymin[i]) && (!is_zmax[i]) )
00648                 {
00649                     local_boxes.insert(i+1-M+M*N);
00650                 }
00651         
00652                 if( (!is_xmax[i]) && (!is_ymax[i]) && (!is_zmin[i]) )
00653                 {
00654                     local_boxes.insert(i+1+M-M*N);
00655                 }
00656         
00657                 if( (!is_xmax[i]) && (!is_ymax[i]) && (!is_zmax[i]) )
00658                 {
00659                     local_boxes.insert(i+1+M+M*N);
00660                 }
00661 
00662                 mLocalBoxes.push_back(local_boxes);
00663             }
00664             break;
00665         }
00666         default:
00667             NEVER_REACHED;
00668     }
00669 }
00670 
00671 
00672 
00673 
00674 
00675 template<unsigned DIM>
00676 std::set<unsigned> BoxCollection<DIM>::GetLocalBoxes(unsigned boxIndex)
00677 {
00678     assert(boxIndex < mLocalBoxes.size());
00679     return mLocalBoxes[boxIndex];
00680 }
00681 
00682 template<unsigned DIM>
00683 void BoxCollection<DIM>::CalculateNodePairs(std::vector<Node<DIM>*>& rNodes, std::set<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs)
00684 {
00685     rNodePairs.clear();
00686     for (unsigned node_index=0; node_index<rNodes.size(); node_index++)
00687     {
00688         // Get the box containing this node
00689         unsigned box_index = CalculateContainingBox(rNodes[node_index]);
00690 
00691         // Get the local boxes to this node
00692         std::set<unsigned> local_boxes_indices = GetLocalBoxes(box_index);
00693 
00694         // Loop over all the local boxes
00695         for (std::set<unsigned>::iterator box_iter = local_boxes_indices.begin();
00696              box_iter != local_boxes_indices.end();
00697              box_iter++)
00698         {
00699             // Get the set of nodes contained in this box
00700             std::set< Node<DIM>* >& r_contained_nodes = mBoxes[*box_iter].rGetNodesContained();
00701 
00702             // Loop over these nodes
00703             for (typename std::set<Node<DIM>*>::iterator node_iter = r_contained_nodes.begin();
00704                  node_iter != r_contained_nodes.end();
00705                  ++node_iter)
00706             {
00707                 // Get the index of the other node
00708                 unsigned other_node_index = (*node_iter)->GetIndex();
00709 
00710                 // If we're in the same box, then take care not to store the node pair twice
00711                 if (*box_iter == box_index)
00712                 {
00713                     if (other_node_index > node_index)
00714                     {
00715                         rNodePairs.insert(std::pair<Node<DIM>*, Node<DIM>*>(rNodes[node_index], rNodes[other_node_index]));
00716                     }
00717                 }
00718                 else
00719                 {
00720                     rNodePairs.insert(std::pair<Node<DIM>*, Node<DIM>*>(rNodes[node_index], rNodes[other_node_index]));
00721                 }
00722             }
00723         }
00724     }
00725 }
00726 
00728 // Explicit instantiation
00730 
00731 template class Box<1>;
00732 template class Box<2>;
00733 template class Box<3>;
00734 template class BoxCollection<1>;
00735 template class BoxCollection<2>;
00736 template class BoxCollection<3>;

Generated on Mon Nov 1 12:35:23 2010 for Chaste by  doxygen 1.5.5