PottsMeshGenerator.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
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 
00029 #include "PottsMeshGenerator.hpp"
00030 
00031 template<unsigned DIM>
00032 PottsMeshGenerator<DIM>::PottsMeshGenerator(unsigned numNodesAcross, unsigned numElementsAcross, unsigned elementWidth,
00033                                             unsigned numNodesUp, unsigned numElementsUp, unsigned elementHeight,
00034                                             unsigned numNodesDeep, unsigned numElementsDeep, unsigned elementDepth,
00035                                             bool startAtBottomLeft, bool isPeriodicInX, bool isPeriodicInY ,bool isPeriodicInZ)
00036 {
00037     assert(numElementsAcross > 0);
00038     assert(numElementsUp > 0);
00039     assert(elementWidth > 0);
00040     assert(elementHeight > 0);
00041     assert(numElementsDeep > 0);
00042     assert(numNodesDeep > 0);
00043     assert(elementDepth > 0);
00044 
00045     assert(numElementsAcross*elementWidth <= numNodesAcross);
00046     assert(numElementsUp*elementHeight <= numNodesUp);
00047     assert(numElementsDeep*elementDepth <= numNodesDeep);
00048 
00049     std::vector<Node<DIM>*> nodes;
00050     std::vector<PottsElement<DIM>*>  elements;
00051     std::vector<std::set<unsigned> > moore_neighbours;
00052     std::vector<std::set<unsigned> > von_neumann_neighbours;
00053 
00054     unsigned num_nodes = numNodesAcross*numNodesUp*numNodesDeep;
00055 
00056     unsigned node_index = 0;
00057     unsigned node_indices[elementWidth*elementHeight*elementDepth];
00058     unsigned element_index;
00059 
00060     unsigned index_offset = 0;
00061 
00062     if (!startAtBottomLeft) // Elements in centre of mesh
00063     {
00064         // Calculate the width of the medium on the edge and offset the node index so that the elements are in the centre of the mesh.
00065         unsigned across_gap = (numNodesAcross -  numElementsAcross*elementWidth)/2;
00066         unsigned up_gap = (numNodesUp -  numElementsUp*elementHeight)/2;
00067         unsigned deep_gap = (numNodesDeep -  numElementsDeep*elementDepth)/2;
00068 
00069         index_offset = deep_gap*numNodesAcross*numNodesUp + up_gap*numNodesAcross + across_gap;
00070     }
00071 
00072     /*
00073      * Create the nodes, row by row, from the bottom up
00074      * On the first and last row we have numNodesAcross nodes, all of which are boundary
00075      * nodes. On each interior row we have numNodesAcross nodes, the first and last nodes
00076      * are boundary nodes.
00077      */
00078     for (unsigned k=0; k<numNodesDeep; k++)
00079     {
00080         for (unsigned j=0; j<numNodesUp; j++)
00081         {
00082             for (unsigned i=0; i<numNodesAcross; i++)
00083             {
00084                 bool is_boundary_node=false;
00085                 if (DIM==2)
00086                 {
00087                     is_boundary_node = (j==0 || j==numNodesUp-1 || (i==0 && !isPeriodicInX) || (i==numNodesAcross-1 && !isPeriodicInX) ) ? true : false;
00088                 }
00089                 if (DIM==3)
00090                 {
00091                     is_boundary_node = (j==0 || j==numNodesUp-1 || (i==0 && !isPeriodicInX) || (i==numNodesAcross-1 && !isPeriodicInX) || k==0 || k==numNodesDeep-1) ? true : false;
00092                 }
00093                 Node<DIM>* p_node = new Node<DIM>(node_index, is_boundary_node, i, j, k);
00094                 nodes.push_back(p_node);
00095                 node_index++;
00096             }
00097         }
00098     }
00099     assert(nodes.size()==num_nodes);
00100 
00101     /*
00102      * Create the elements. The array node_indices contains the
00103      * global node indices, in increasing order.
00104      */
00105     for (unsigned n=0; n<numElementsDeep; n++)
00106     {
00107         for (unsigned j=0; j<numElementsUp; j++)
00108         {
00109             for (unsigned i=0; i<numElementsAcross; i++)
00110             {
00111                 for (unsigned m=0; m<elementDepth; m++)
00112                 {
00113                     for (unsigned l=0; l<elementHeight; l++)
00114                     {
00115                         for (unsigned k=0; k<elementWidth; k++)
00116                         {
00117                             node_indices[m*elementHeight*elementWidth + l*elementWidth + k] = n*elementDepth*numNodesUp*numNodesAcross +
00118                                                                                               j*elementHeight*numNodesAcross +
00119                                                                                               i*elementWidth +
00120                                                                                               m*numNodesAcross*numNodesUp +
00121                                                                                               l*numNodesAcross +
00122                                                                                               k + index_offset;
00123                         }
00124                     }
00125                 }
00126                 std::vector<Node<DIM>*> element_nodes;
00127                 for (unsigned k=0; k<elementDepth*elementHeight*elementWidth; k++)
00128                 {
00129                    element_nodes.push_back(nodes[node_indices[k]]);
00130                 }
00131 
00132                 element_index = n*numElementsAcross*numElementsUp + j*numElementsAcross + i;
00133                 PottsElement<DIM>* p_element = new PottsElement<DIM>(element_index, element_nodes);
00134                 elements.push_back(p_element);
00135             }
00136         }
00137     }
00138 
00139     /*
00140      * Create the neighborhoods of each node
00141      */
00142 
00143     moore_neighbours.resize(num_nodes);
00144     von_neumann_neighbours.resize(num_nodes);
00145 
00146     for (unsigned node_index=0; node_index<num_nodes; node_index++)
00147     {
00148         // Clear the set of neighbouring node indices
00149 
00150         moore_neighbours[node_index].clear();
00151 
00152         switch (DIM)
00153         {
00154             case 2:
00155             {
00156                 assert(DIM == 2);
00157                 /*
00158                  * This stores the available neighbours using the following numbering:
00159                  *
00160                  *  1-----0-----7
00161                  *  |     |     |
00162                  *  |     |     |
00163                  *  2-----x-----6
00164                  *  |     |     |
00165                  *  |     |     |
00166                  *  3-----4-----5
00167                  */
00168 
00169                 // Create a vector of possible neighbouring node indices
00170                 std::vector<unsigned> moore_neighbour_indices_vector(8, node_index);
00171                 moore_neighbour_indices_vector[0] += numNodesAcross;
00172                 moore_neighbour_indices_vector[1] += numNodesAcross - 1;
00173                 moore_neighbour_indices_vector[2] -= 1;
00174                 moore_neighbour_indices_vector[3] -= numNodesAcross + 1;
00175                 moore_neighbour_indices_vector[4] -= numNodesAcross;
00176                 moore_neighbour_indices_vector[5] -= numNodesAcross - 1;
00177                 moore_neighbour_indices_vector[6] += 1;
00178                 moore_neighbour_indices_vector[7] += numNodesAcross + 1;
00179 
00180                 // Work out whether this node lies on any edge of the mesh
00181                 bool on_south_edge = (node_index < numNodesAcross);
00182                 bool on_north_edge = (node_index > numNodesAcross*(numNodesUp - 1) - 1);
00183                 bool on_west_edge = (node_index%numNodesAcross == 0);
00184                 bool on_east_edge = (node_index%numNodesAcross == numNodesAcross - 1);
00185 
00186                 if(isPeriodicInX)
00187                 {
00188                     if(on_west_edge)
00189                     {
00190                         moore_neighbour_indices_vector[1] = node_index + 2*numNodesAcross - 1;
00191                         moore_neighbour_indices_vector[2] = node_index + numNodesAcross - 1;
00192                         moore_neighbour_indices_vector[3] = node_index - 1;
00193                     }
00194 
00195                     if(on_east_edge)
00196                     {
00197                         moore_neighbour_indices_vector[5] = node_index - 2*numNodesAcross + 1;
00198                         moore_neighbour_indices_vector[6] = node_index - numNodesAcross + 1;
00199                         moore_neighbour_indices_vector[7] = node_index + 1;
00200                     }
00201                 }
00202 
00203                 if(isPeriodicInY)
00204                 {
00205                     if(on_north_edge)
00206                     {
00207                         moore_neighbour_indices_vector[0] = node_index - numNodesAcross*(numNodesUp-1);
00208                         moore_neighbour_indices_vector[1] = moore_neighbour_indices_vector[0] - 1;
00209                         moore_neighbour_indices_vector[7] = moore_neighbour_indices_vector[0] + 1;
00210 
00211                         if(on_west_edge)
00212                         {
00213                             moore_neighbour_indices_vector[1] = moore_neighbour_indices_vector[1] + numNodesAcross;
00214                         }
00215                         if(on_east_edge)
00216                         {
00217                             moore_neighbour_indices_vector[7] = moore_neighbour_indices_vector[7] - numNodesAcross;
00218                         }
00219                     }
00220 
00221                     if(on_south_edge)
00222                     {
00223                         moore_neighbour_indices_vector[4] = node_index + numNodesAcross*(numNodesUp-1);
00224                         moore_neighbour_indices_vector[3] = moore_neighbour_indices_vector[4] - 1;
00225                         moore_neighbour_indices_vector[5] = moore_neighbour_indices_vector[4] + 1;
00226 
00227                         if(on_west_edge)
00228                         {
00229                             moore_neighbour_indices_vector[3] = moore_neighbour_indices_vector[3] + numNodesAcross;
00230                         }
00231                         if(on_east_edge)
00232                         {
00233                             moore_neighbour_indices_vector[5] = moore_neighbour_indices_vector[5] - numNodesAcross;
00234                         }
00235                     }
00236                 }
00237 
00238 
00239                 if(isPeriodicInX)
00240                 {
00241                     on_east_edge = false;
00242                     on_west_edge = false;
00243                 }
00244                 if(isPeriodicInY)
00245                 {
00246                     on_south_edge = false;
00247                     on_north_edge = false;
00248                 }
00249 
00250                 // Create a vector of booleans for which neighbours are available
00251                 // Use the order N, NW, W, SW, S, SE, E, NE
00252                 std::vector<bool> available_neighbours = std::vector<bool>(8, true);
00253                 available_neighbours[0] = !on_north_edge;
00254                 available_neighbours[1] = !(on_north_edge || on_west_edge);
00255                 available_neighbours[2] = !on_west_edge;
00256                 available_neighbours[3] = !(on_south_edge || on_west_edge);
00257                 available_neighbours[4] = !on_south_edge;
00258                 available_neighbours[5] = !(on_south_edge || on_east_edge);
00259                 available_neighbours[6] = !on_east_edge;
00260                 available_neighbours[7] = !(on_north_edge || on_east_edge);
00261 
00262                 // Using neighbour_indices_vector and available_neighbours, store the indices of all available neighbours to the set all_neighbours
00263                 for (unsigned i=0; i<8; i++)
00264                 {
00265                     if (available_neighbours[i])
00266                     {
00267                         assert(moore_neighbour_indices_vector[i] < nodes.size());
00268                         moore_neighbours[node_index].insert(moore_neighbour_indices_vector[i]);
00269                     }
00270                 }
00271                 break;
00272             }
00273             case 3:
00274             {
00275                 assert(DIM ==3 );
00276                 /*
00277                  * This stores the available neighbours using the following numbering:
00278                  *                      FRONT           BACK
00279                  *  1-----0-----7   10-----9---16   19----18---25
00280                  *  |     |     |   |     |     |   |     |     |
00281                  *  |     |     |   |     |     |   |     |     |
00282                  *  2-----x-----6   11----8-----15  20----17---24
00283                  *  |     |     |   |     |     |   |     |     |
00284                  *  |     |     |   |     |     |   |     |     |
00285                  *  3-----4-----5   12----13----14  21---22----23
00286                  */
00287 
00288                 // Create a vector of possible neighbouring node indices
00289                 std::vector<unsigned> moore_neighbour_indices_vector(26, node_index);
00290                 moore_neighbour_indices_vector[0] += numNodesAcross;
00291                 moore_neighbour_indices_vector[1] += numNodesAcross - 1;
00292                 moore_neighbour_indices_vector[2] -= 1;
00293                 moore_neighbour_indices_vector[3] -= numNodesAcross + 1;
00294                 moore_neighbour_indices_vector[4] -= numNodesAcross;
00295                 moore_neighbour_indices_vector[5] -= numNodesAcross - 1;
00296                 moore_neighbour_indices_vector[6] += 1;
00297                 moore_neighbour_indices_vector[7] += numNodesAcross + 1;
00298                 moore_neighbour_indices_vector[8] -= numNodesAcross*numNodesUp;
00299                 for( unsigned i=9; i<17; i++)
00300                 {
00301                     moore_neighbour_indices_vector[i] = moore_neighbour_indices_vector[i-9]-numNodesAcross*numNodesUp;
00302                 }
00303                 moore_neighbour_indices_vector[17] += numNodesAcross*numNodesUp;
00304                 for (unsigned i=18; i<26; i++)
00305                 {
00306                     moore_neighbour_indices_vector[i]=moore_neighbour_indices_vector[i-18]+numNodesAcross*numNodesUp;
00307                 }
00308 
00309                 // Work out whether this node lies on any edge of the mesh
00310                 bool on_south_edge = (node_index%(numNodesAcross*numNodesUp)<numNodesAcross);
00311                 bool on_north_edge = (node_index%(numNodesAcross*numNodesUp)>(numNodesAcross*numNodesUp-numNodesAcross-1));
00312                 bool on_west_edge = (node_index%numNodesAcross == 0);
00313                 bool on_east_edge = (node_index%numNodesAcross == numNodesAcross - 1);
00314                 bool on_front_edge = (node_index < numNodesAcross*numNodesUp);
00315                 bool on_back_edge = (node_index > numNodesAcross*numNodesUp*(numNodesDeep-1)-1);
00316 
00317                 if(isPeriodicInX)
00318                 {
00319                     if(on_west_edge)
00320                     {
00321                         moore_neighbour_indices_vector[1] = node_index + 2*numNodesAcross - 1;
00322                         moore_neighbour_indices_vector[2] = node_index + numNodesAcross - 1;
00323                         moore_neighbour_indices_vector[3] = node_index - 1;
00324 
00325                         moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[1] - numNodesAcross*numNodesUp;
00326                         moore_neighbour_indices_vector[11] = moore_neighbour_indices_vector[2] - numNodesAcross*numNodesUp;
00327                         moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[3] - numNodesAcross*numNodesUp;
00328 
00329                         moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[1] + numNodesAcross*numNodesUp;
00330                         moore_neighbour_indices_vector[20] = moore_neighbour_indices_vector[2] + numNodesAcross*numNodesUp;
00331                         moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[3] + numNodesAcross*numNodesUp;
00332                     }
00333 
00334                     if(on_east_edge)
00335                     {
00336                         moore_neighbour_indices_vector[5] = node_index - 2*numNodesAcross + 1;
00337                         moore_neighbour_indices_vector[6] = node_index - numNodesAcross + 1;
00338                         moore_neighbour_indices_vector[7] = node_index + 1;
00339 
00340                         moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[5] - numNodesAcross*numNodesUp;
00341                         moore_neighbour_indices_vector[15] = moore_neighbour_indices_vector[6] - numNodesAcross*numNodesUp;
00342                         moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[7] - numNodesAcross*numNodesUp;
00343 
00344                         moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[5] + numNodesAcross*numNodesUp;
00345                         moore_neighbour_indices_vector[24] = moore_neighbour_indices_vector[6] + numNodesAcross*numNodesUp;
00346                         moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[7] + numNodesAcross*numNodesUp;
00347                     }
00348                 }
00349 
00350                 if(isPeriodicInY)
00351                 {
00352                     if(on_north_edge)
00353                     {
00354                         moore_neighbour_indices_vector[0] = node_index - numNodesAcross*(numNodesUp-1);
00355                         moore_neighbour_indices_vector[1] = moore_neighbour_indices_vector[0] - 1;
00356                         moore_neighbour_indices_vector[7] = moore_neighbour_indices_vector[0] + 1;
00357 
00358                         moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[1] - numNodesAcross*numNodesUp;
00359                         moore_neighbour_indices_vector[9] = moore_neighbour_indices_vector[0] - numNodesAcross*numNodesUp;
00360                         moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[7] - numNodesAcross*numNodesUp;
00361 
00362                         moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[1] + numNodesAcross*numNodesUp;
00363                         moore_neighbour_indices_vector[18] = moore_neighbour_indices_vector[0] + numNodesAcross*numNodesUp;
00364                         moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[7] + numNodesAcross*numNodesUp;
00365 
00366                         if(on_west_edge)
00367                         {
00368                             moore_neighbour_indices_vector[1] = moore_neighbour_indices_vector[1] + numNodesAcross;
00369                             moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[10] + numNodesAcross;
00370                             moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[19] + numNodesAcross;
00371                         }
00372                         if(on_east_edge)
00373                         {
00374                             moore_neighbour_indices_vector[7] = moore_neighbour_indices_vector[7] - numNodesAcross;
00375                             moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[16] - numNodesAcross;
00376                             moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[25] - numNodesAcross;
00377                         }
00378                     }
00379 
00380                     if(on_south_edge)
00381                     {
00382                         moore_neighbour_indices_vector[4] = node_index + numNodesAcross*(numNodesUp-1);
00383                         moore_neighbour_indices_vector[3] = moore_neighbour_indices_vector[4] - 1;
00384                         moore_neighbour_indices_vector[5] = moore_neighbour_indices_vector[4] + 1;
00385 
00386                         moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[3] - numNodesAcross*numNodesUp;
00387                         moore_neighbour_indices_vector[13] = moore_neighbour_indices_vector[4] - numNodesAcross*numNodesUp;
00388                         moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[5] - numNodesAcross*numNodesUp;
00389 
00390                         moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[3] + numNodesAcross*numNodesUp;
00391                         moore_neighbour_indices_vector[22] = moore_neighbour_indices_vector[4] + numNodesAcross*numNodesUp;
00392                         moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[5] + numNodesAcross*numNodesUp;
00393 
00394                         if(on_west_edge)
00395                         {
00396                             moore_neighbour_indices_vector[3] = moore_neighbour_indices_vector[3] + numNodesAcross;
00397                             moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[12] + numNodesAcross;
00398                             moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[21] + numNodesAcross;
00399                         }
00400                         if(on_east_edge)
00401                         {
00402                             moore_neighbour_indices_vector[5] = moore_neighbour_indices_vector[5] - numNodesAcross;
00403                             moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[14] - numNodesAcross;
00404                             moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[23] - numNodesAcross;
00405                         }
00406                     }
00407                 }
00408 
00409                 if(isPeriodicInZ)
00410                 {
00411                     if(on_back_edge)
00412                     {
00413                         moore_neighbour_indices_vector[17] = node_index - numNodesAcross*numNodesUp*(numNodesDeep-1);
00414                         moore_neighbour_indices_vector[20] = moore_neighbour_indices_vector[17] - 1;
00415                         moore_neighbour_indices_vector[24] = moore_neighbour_indices_vector[17] + 1;
00416 
00417                         moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[20] - numNodesAcross;
00418                         moore_neighbour_indices_vector[22] = moore_neighbour_indices_vector[17] - numNodesAcross;
00419                         moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[24] - numNodesAcross;
00420 
00421                         moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[20] + numNodesAcross;
00422                         moore_neighbour_indices_vector[18] = moore_neighbour_indices_vector[17] + numNodesAcross;
00423                         moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[24] + numNodesAcross;
00424 
00425                         if(on_west_edge)
00426                         {
00427                             moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[19] + numNodesAcross;
00428                             moore_neighbour_indices_vector[20] = moore_neighbour_indices_vector[20] + numNodesAcross;
00429                             moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[21] + numNodesAcross;
00430                         }
00431                         if(on_east_edge)
00432                         {
00433                             moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[23] - numNodesAcross;
00434                             moore_neighbour_indices_vector[24] = moore_neighbour_indices_vector[24] - numNodesAcross;
00435                             moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[25] - numNodesAcross;
00436                         }
00437 
00438                         if(on_south_edge)
00439                         {
00440                             moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[21] + numNodesAcross*numNodesUp;
00441                             moore_neighbour_indices_vector[22] = moore_neighbour_indices_vector[22] + numNodesAcross*numNodesUp;
00442                             moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[23] + numNodesAcross*numNodesUp;
00443                         }
00444 
00445                         if(on_north_edge)
00446                         {
00447                             moore_neighbour_indices_vector[18] = moore_neighbour_indices_vector[18] - numNodesAcross*numNodesUp;
00448                             moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[19] - numNodesAcross*numNodesUp;
00449                             moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[25] - numNodesAcross*numNodesUp;
00450                         }
00451                     }
00452 
00453                     if(on_front_edge)
00454                     {
00455                         moore_neighbour_indices_vector[8] = node_index + numNodesAcross*numNodesUp*(numNodesDeep-1);
00456                         moore_neighbour_indices_vector[11] = moore_neighbour_indices_vector[8] - 1;
00457                         moore_neighbour_indices_vector[15] = moore_neighbour_indices_vector[8] + 1;
00458 
00459                         moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[11] - numNodesAcross;
00460                         moore_neighbour_indices_vector[13] = moore_neighbour_indices_vector[8] - numNodesAcross;
00461                         moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[15] - numNodesAcross;
00462 
00463                         moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[11] + numNodesAcross;
00464                         moore_neighbour_indices_vector[9] = moore_neighbour_indices_vector[8] + numNodesAcross;
00465                         moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[15] + numNodesAcross;
00466 
00467                         if(on_west_edge)
00468                         {
00469                             moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[10] + numNodesAcross;
00470                             moore_neighbour_indices_vector[11] = moore_neighbour_indices_vector[11] + numNodesAcross;
00471                             moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[12] + numNodesAcross;
00472                         }
00473                         if(on_east_edge)
00474                         {
00475                             moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[14] - numNodesAcross;
00476                             moore_neighbour_indices_vector[15] = moore_neighbour_indices_vector[15] - numNodesAcross;
00477                             moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[16] - numNodesAcross;
00478                         }
00479 
00480                         if(on_south_edge)
00481                         {
00482                             moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[12] + numNodesAcross*numNodesUp;
00483                             moore_neighbour_indices_vector[13] = moore_neighbour_indices_vector[13] + numNodesAcross*numNodesUp;
00484                             moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[14] + numNodesAcross*numNodesUp;
00485                         }
00486 
00487                         if(on_north_edge)
00488                         {
00489                             moore_neighbour_indices_vector[9] = moore_neighbour_indices_vector[9] - numNodesAcross*numNodesUp;
00490                             moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[10] - numNodesAcross*numNodesUp;
00491                             moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[16] - numNodesAcross*numNodesUp;
00492                         }
00493                     }
00494 
00495 
00496                 }
00497 
00498                 if(isPeriodicInX)
00499                 {
00500                     on_east_edge = false;
00501                     on_west_edge = false;
00502                 }
00503                 if(isPeriodicInY)
00504                 {
00505                     on_south_edge = false;
00506                     on_north_edge = false;
00507                 }
00508                 if(isPeriodicInZ)
00509                 {
00510                     on_front_edge = false;
00511                     on_back_edge = false;
00512                 }
00513 
00514 
00515                 // Create a vector of booleans for which neighbours are available
00516                 // Use the order N, NW, W, SW, S, SE, E, NE
00517                 std::vector<bool> available_neighbours = std::vector<bool>(26, true);
00518                 available_neighbours[0] = !on_north_edge;
00519                 available_neighbours[1] = !(on_north_edge || on_west_edge);
00520                 available_neighbours[2] = !on_west_edge;
00521                 available_neighbours[3] = !(on_south_edge || on_west_edge);
00522                 available_neighbours[4] = !on_south_edge;
00523                 available_neighbours[5] = !(on_south_edge || on_east_edge);
00524                 available_neighbours[6] = !on_east_edge;
00525                 available_neighbours[7] = !(on_north_edge || on_east_edge);
00526                 available_neighbours[8] = !(on_front_edge);
00527                 for (unsigned i=9; i<17; i++)
00528                 {
00529                     available_neighbours[i] = (available_neighbours[i-9] && !(on_front_edge));
00530                 }
00531                 available_neighbours[17] = !(on_back_edge);
00532                 for (unsigned i=18; i<26; i++)
00533                 {
00534                     available_neighbours[i] = (available_neighbours[i-18] && !(on_back_edge));
00535                 }
00536 
00537                 // Using neighbour_indices_vector and available_neighbours, store the indices of all available neighbours to the set all_neighbours
00538                 for (unsigned i=0; i<26; i++)
00539                 {
00540                     if (available_neighbours[i] && moore_neighbour_indices_vector[i] < numNodesAcross*numNodesUp*numNodesDeep)
00541                     {
00542                         assert(moore_neighbour_indices_vector[i] < nodes.size());
00543                         moore_neighbours[node_index].insert(moore_neighbour_indices_vector[i]);
00544                     }
00545                 }
00546                 break;
00547             }
00548             default:
00549                 NEVER_REACHED;
00550         }
00551 
00552         // Clear the set of neighbouring node indices
00553         von_neumann_neighbours[node_index].clear();
00554 
00555         switch (DIM)
00556         {
00557             case 2:
00558             {
00559                 assert(DIM == 2);
00560                 /*
00561                  * This stores the available neighbours using the following numbering:
00562                  *
00563                  *        0
00564                  *        |
00565                  *        |
00566                  *  1-----x-----3
00567                  *        |
00568                  *        |
00569                  *        2
00570                  */
00571                 // Create a vector of possible neighbouring node indices
00572                 std::vector<unsigned> von_neumann_neighbour_indices_vector(4, node_index);
00573                 von_neumann_neighbour_indices_vector[0] += numNodesAcross;
00574                 von_neumann_neighbour_indices_vector[1] -= 1;
00575                 von_neumann_neighbour_indices_vector[2] -= numNodesAcross;
00576                 von_neumann_neighbour_indices_vector[3] += 1;
00577 
00578                 // Work out whether this node lies on any edge of the mesh
00579                 bool on_south_edge = (node_index < numNodesAcross);
00580                 bool on_north_edge = (node_index > numNodesAcross*(numNodesUp - 1) - 1);
00581                 bool on_west_edge = (node_index%numNodesAcross == 0);
00582                 bool on_east_edge = (node_index%numNodesAcross == numNodesAcross - 1);
00583 
00584 
00585 
00586                 if(isPeriodicInX)
00587                 {
00588                     if(on_west_edge)
00589                     {
00590                         von_neumann_neighbour_indices_vector[1] = node_index + numNodesAcross - 1;
00591                         on_west_edge = false;
00592                     }
00593 
00594                     if(on_east_edge)
00595                     {
00596                         von_neumann_neighbour_indices_vector[3] = node_index - numNodesAcross + 1;
00597                         on_east_edge = false;
00598                     }
00599                 }
00600 
00601                 if(isPeriodicInY)
00602                 {
00603                     if(on_north_edge)
00604                     {
00605                         von_neumann_neighbour_indices_vector[0] = node_index - numNodesAcross*(numNodesUp-1);
00606                         on_north_edge = false;
00607                     }
00608 
00609                     if(on_south_edge)
00610                     {
00611                         von_neumann_neighbour_indices_vector[2] = node_index + numNodesAcross*(numNodesUp-1);
00612                         on_south_edge = false;
00613                     }
00614                 }
00615 
00616                 // Create a vector of booleans for which neighbours are available
00617                 // Use the order N, W, S, E
00618                 std::vector<bool> available_neighbours = std::vector<bool>(2*DIM, true);
00619                 available_neighbours[0] = !on_north_edge;
00620                 available_neighbours[1] = !on_west_edge;
00621                 available_neighbours[2] = !on_south_edge;
00622                 available_neighbours[3] = !on_east_edge;
00623 
00624                 // Using von_neumann_neighbour_indices_vector and available_neighbours, store the indices of all available neighbours to the set all_neighbours
00625                 for (unsigned i=0; i<4; i++)
00626                 {
00627                     if (available_neighbours[i])
00628                     {
00629                         assert(von_neumann_neighbour_indices_vector[i] < nodes.size());
00630                         von_neumann_neighbours[node_index].insert(von_neumann_neighbour_indices_vector[i]);
00631                     }
00632                 }
00633                 break;
00634             }
00635             case 3:
00636             {
00637                 assert(DIM == 3);
00638 
00639                 /*
00640                  * This stores the available neighbours using the following numbering:
00641                  *
00642                  *        0  5
00643                  *        | /
00644                  *        |/
00645                  *  1-----x-----3
00646                  *      / |
00647                  *     /  |
00648                  *    4   2
00649                  */
00650                 // Create a vector of possible neighbouring node indices
00651                 std::vector<unsigned> von_neumann_neighbour_indices_vector(6, node_index);
00652                 von_neumann_neighbour_indices_vector[0] += numNodesAcross;
00653                 von_neumann_neighbour_indices_vector[1] -= 1;
00654                 von_neumann_neighbour_indices_vector[2] -= numNodesAcross;
00655                 von_neumann_neighbour_indices_vector[3] += 1;
00656                 von_neumann_neighbour_indices_vector[4] -= numNodesAcross*numNodesUp;
00657                 von_neumann_neighbour_indices_vector[5] += numNodesAcross*numNodesUp;
00658 
00659                 // Work out whether this node lies on any edge of the mesh
00660                 bool on_south_edge = (node_index%(numNodesAcross*numNodesUp)<numNodesAcross);
00661                 bool on_north_edge = (node_index%(numNodesAcross*numNodesUp)>(numNodesAcross*numNodesUp-numNodesAcross-1));
00662                 bool on_west_edge = (node_index%numNodesAcross== 0);
00663                 bool on_east_edge = (node_index%numNodesAcross == numNodesAcross - 1);
00664                 bool on_front_edge = (node_index < numNodesAcross*numNodesUp);
00665                 bool on_back_edge = (node_index > numNodesAcross*numNodesUp*(numNodesDeep-1)-1);
00666 
00667 
00668                 if(isPeriodicInX)
00669                 {
00670                     if(on_west_edge)
00671                     {
00672                         von_neumann_neighbour_indices_vector[1] = node_index + numNodesAcross - 1;
00673                         on_west_edge = false;
00674                     }
00675 
00676                     if(on_east_edge)
00677                     {
00678                         von_neumann_neighbour_indices_vector[3] = node_index - numNodesAcross + 1;
00679                         on_east_edge = false;
00680                     }
00681                 }
00682 
00683                 if(isPeriodicInY)
00684                 {
00685                     if(on_north_edge)
00686                     {
00687                         von_neumann_neighbour_indices_vector[0] = node_index - numNodesAcross*(numNodesUp-1);
00688                         on_north_edge = false;
00689                     }
00690 
00691                     if(on_south_edge)
00692                     {
00693                         von_neumann_neighbour_indices_vector[2] = node_index + numNodesAcross*(numNodesUp-1);
00694                         on_south_edge = false;
00695                     }
00696                 }
00697 
00698                 if(isPeriodicInZ)
00699                 {
00700                     if(on_front_edge)
00701                     {
00702                         von_neumann_neighbour_indices_vector[4] = node_index + numNodesAcross*numNodesUp*(numNodesDeep-1);
00703                         on_front_edge = false;
00704                     }
00705 
00706                     if(on_back_edge)
00707                     {
00708                         von_neumann_neighbour_indices_vector[5] = node_index - numNodesAcross*numNodesUp*(numNodesDeep-1);
00709                         on_back_edge = false;
00710                     }
00711                 }
00712 
00713 
00714 
00715 
00716                 // Create a vector of booleans for which neighbours are available
00717                 // Use the order N, W, S, E, F, B
00718                 std::vector<bool> available_neighbours = std::vector<bool>(2*DIM, true);
00719                 available_neighbours[0] = !on_north_edge;
00720                 available_neighbours[1] = !on_west_edge;
00721                 available_neighbours[2] = !on_south_edge;
00722                 available_neighbours[3] = !on_east_edge;
00723                 available_neighbours[4] = !on_front_edge;
00724                 available_neighbours[5] = !on_back_edge;
00725 
00726                 // Using von_neumann_neighbour_indices_vector and available_neighbours, store the indices of all available neighbours to the set all_neighbours
00727                 for (unsigned i=0; i<6; i++)
00728                 {
00729                     if (available_neighbours[i] && von_neumann_neighbour_indices_vector[i]<numNodesAcross*numNodesUp*numNodesDeep)
00730                     {
00731                         assert(von_neumann_neighbour_indices_vector[i] < nodes.size());
00732                         von_neumann_neighbours[node_index].insert(von_neumann_neighbour_indices_vector[i]);
00733                     }
00734                 }
00735                 break;
00736             }
00737             default:
00738                 NEVER_REACHED;
00739         }
00740     }
00741 
00742 
00743     mpMesh = new PottsMesh<DIM>(nodes, elements, von_neumann_neighbours, moore_neighbours);
00744 }
00745 
00746 template<unsigned DIM>
00747 PottsMeshGenerator<DIM>::~PottsMeshGenerator()
00748 {
00749     delete mpMesh;
00750 }
00751 
00752 template<unsigned DIM>
00753 PottsMesh<DIM>* PottsMeshGenerator<DIM>::GetMesh()
00754 {
00755     return mpMesh;
00756 }
00757 
00759 // Explicit instantiation
00761 
00762 template class PottsMeshGenerator<1>;
00763 template class PottsMeshGenerator<2>;
00764 template class PottsMeshGenerator<3>;
Generated on Thu Dec 22 13:00:04 2011 for Chaste by  doxygen 1.6.3