Chaste  Release::3.4
PottsMeshGenerator.cpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "PottsMeshGenerator.hpp"
37 
38 #include <boost/scoped_array.hpp>
39 
40 template<unsigned DIM>
41 PottsMeshGenerator<DIM>::PottsMeshGenerator(unsigned numNodesAcross, unsigned numElementsAcross, unsigned elementWidth,
42  unsigned numNodesUp, unsigned numElementsUp, unsigned elementHeight,
43  unsigned numNodesDeep, unsigned numElementsDeep, unsigned elementDepth,
44  bool startAtBottomLeft, bool isPeriodicInX, bool isPeriodicInY ,bool isPeriodicInZ)
45 {
46  assert(numNodesAcross > 0);
47  assert(numNodesUp > 0);
48  assert(numNodesDeep > 0);
49 
50  assert(numElementsAcross*elementWidth <= numNodesAcross);
51  assert(numElementsUp*elementHeight <= numNodesUp);
52  assert(numElementsDeep*elementDepth <= numNodesDeep);
53 
54  std::vector<Node<DIM>*> nodes;
55  std::vector<PottsElement<DIM>*> elements;
56  std::vector<std::set<unsigned> > moore_neighbours;
57  std::vector<std::set<unsigned> > von_neumann_neighbours;
58 
59  unsigned num_nodes = numNodesAcross*numNodesUp*numNodesDeep;
60 
61  unsigned next_node_index = 0;
62  boost::scoped_array<unsigned> node_indices(new unsigned[elementWidth*elementHeight*elementDepth]);
63  unsigned element_index;
64 
65  unsigned index_offset = 0;
66 
67  if (!startAtBottomLeft) // Elements in centre of mesh
68  {
69  // 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.
70  unsigned across_gap = (numNodesAcross - numElementsAcross*elementWidth)/2;
71  unsigned up_gap = (numNodesUp - numElementsUp*elementHeight)/2;
72  unsigned deep_gap = (numNodesDeep - numElementsDeep*elementDepth)/2;
73 
74  index_offset = deep_gap*numNodesAcross*numNodesUp + up_gap*numNodesAcross + across_gap;
75  }
76 
77  /*
78  * Create the nodes, row by row, from the bottom up
79  * On the first and last row we have numNodesAcross nodes, all of which are boundary
80  * nodes. On each interior row we have numNodesAcross nodes, the first and last nodes
81  * are boundary nodes.
82  */
83  for (unsigned k=0; k<numNodesDeep; k++)
84  {
85  for (unsigned j=0; j<numNodesUp; j++)
86  {
87  for (unsigned i=0; i<numNodesAcross; i++)
88  {
89  bool is_boundary_node=false;
90  if (DIM==2)
91  {
92  is_boundary_node = (j==0 || j==numNodesUp-1 || (i==0 && !isPeriodicInX) || (i==numNodesAcross-1 && !isPeriodicInX) ) ? true : false;
93  }
94  if (DIM==3)
95  {
96  is_boundary_node = (j==0 || j==numNodesUp-1 || (i==0 && !isPeriodicInX) || (i==numNodesAcross-1 && !isPeriodicInX) || k==0 || k==numNodesDeep-1) ? true : false;
97  }
98  Node<DIM>* p_node = new Node<DIM>(next_node_index, is_boundary_node, i, j, k);
99  nodes.push_back(p_node);
100  next_node_index++;
101  }
102  }
103  }
104  assert(nodes.size()==num_nodes);
105 
106  /*
107  * Create the elements. The array node_indices contains the
108  * global node indices, in increasing order.
109  */
110  for (unsigned n=0; n<numElementsDeep; n++)
111  {
112  for (unsigned j=0; j<numElementsUp; j++)
113  {
114  for (unsigned i=0; i<numElementsAcross; i++)
115  {
116  for (unsigned m=0; m<elementDepth; m++)
117  {
118  for (unsigned l=0; l<elementHeight; l++)
119  {
120  for (unsigned k=0; k<elementWidth; k++)
121  {
122  node_indices[m*elementHeight*elementWidth + l*elementWidth + k] = n*elementDepth*numNodesUp*numNodesAcross +
123  j*elementHeight*numNodesAcross +
124  i*elementWidth +
125  m*numNodesAcross*numNodesUp +
126  l*numNodesAcross +
127  k + index_offset;
128  }
129  }
130  }
131  std::vector<Node<DIM>*> element_nodes;
132  for (unsigned k=0; k<elementDepth*elementHeight*elementWidth; k++)
133  {
134  element_nodes.push_back(nodes[node_indices[k]]);
135  }
136 
137  element_index = n*numElementsAcross*numElementsUp + j*numElementsAcross + i;
138  PottsElement<DIM>* p_element = new PottsElement<DIM>(element_index, element_nodes);
139  elements.push_back(p_element);
140  }
141  }
142  }
143 
144  /*
145  * Create the neighborhoods of each node
146  */
147 
148  moore_neighbours.resize(num_nodes);
149  von_neumann_neighbours.resize(num_nodes);
150 
151  for (unsigned node_index=0; node_index<num_nodes; node_index++)
152  {
153  // Clear the set of neighbouring node indices
154 
155  moore_neighbours[node_index].clear();
156 
157  switch (DIM)
158  {
159  case 2:
160  {
161  assert(DIM == 2);
162  /*
163  * This stores the available neighbours using the following numbering:
164  *
165  * 1-----0-----7
166  * | | |
167  * | | |
168  * 2-----x-----6
169  * | | |
170  * | | |
171  * 3-----4-----5
172  */
173 
174  // Create a vector of possible neighbouring node indices
175  std::vector<unsigned> moore_neighbour_indices_vector(8, node_index);
176  moore_neighbour_indices_vector[0] += numNodesAcross;
177  moore_neighbour_indices_vector[1] += numNodesAcross - 1;
178  moore_neighbour_indices_vector[2] -= 1;
179  moore_neighbour_indices_vector[3] -= numNodesAcross + 1;
180  moore_neighbour_indices_vector[4] -= numNodesAcross;
181  moore_neighbour_indices_vector[5] -= numNodesAcross - 1;
182  moore_neighbour_indices_vector[6] += 1;
183  moore_neighbour_indices_vector[7] += numNodesAcross + 1;
184 
185  // Work out whether this node lies on any edge of the mesh
186  bool on_south_edge = (node_index < numNodesAcross);
187  bool on_north_edge = ((int) node_index > (int)numNodesAcross*((int)numNodesUp - 1) - 1);
188  bool on_west_edge = (node_index%numNodesAcross == 0);
189  bool on_east_edge = (node_index%numNodesAcross == numNodesAcross - 1);
190 
191  if (isPeriodicInX)
192  {
193  if (on_west_edge)
194  {
195  moore_neighbour_indices_vector[1] = node_index + 2*numNodesAcross - 1;
196  moore_neighbour_indices_vector[2] = node_index + numNodesAcross - 1;
197  moore_neighbour_indices_vector[3] = node_index - 1;
198  }
199 
200  if (on_east_edge)
201  {
202  moore_neighbour_indices_vector[5] = node_index - 2*numNodesAcross + 1;
203  moore_neighbour_indices_vector[6] = node_index - numNodesAcross + 1;
204  moore_neighbour_indices_vector[7] = node_index + 1;
205  }
206  }
207 
208  if (isPeriodicInY)
209  {
210  if (on_north_edge)
211  {
212  moore_neighbour_indices_vector[0] = node_index - numNodesAcross*(numNodesUp-1);
213  moore_neighbour_indices_vector[1] = moore_neighbour_indices_vector[0] - 1;
214  moore_neighbour_indices_vector[7] = moore_neighbour_indices_vector[0] + 1;
215 
216  if (on_west_edge)
217  {
218  moore_neighbour_indices_vector[1] = moore_neighbour_indices_vector[1] + numNodesAcross;
219  }
220  if (on_east_edge)
221  {
222  moore_neighbour_indices_vector[7] = moore_neighbour_indices_vector[7] - numNodesAcross;
223  }
224  }
225 
226  if (on_south_edge)
227  {
228  moore_neighbour_indices_vector[4] = node_index + numNodesAcross*(numNodesUp-1);
229  moore_neighbour_indices_vector[3] = moore_neighbour_indices_vector[4] - 1;
230  moore_neighbour_indices_vector[5] = moore_neighbour_indices_vector[4] + 1;
231 
232  if (on_west_edge)
233  {
234  moore_neighbour_indices_vector[3] = moore_neighbour_indices_vector[3] + numNodesAcross;
235  }
236  if (on_east_edge)
237  {
238  moore_neighbour_indices_vector[5] = moore_neighbour_indices_vector[5] - numNodesAcross;
239  }
240  }
241  }
242 
243  if (isPeriodicInX)
244  {
245  on_east_edge = false;
246  on_west_edge = false;
247  }
248  if (isPeriodicInY)
249  {
250  on_south_edge = false;
251  on_north_edge = false;
252  }
253 
254  // Create a vector of booleans for which neighbours are available
255  // Use the order N, NW, W, SW, S, SE, E, NE
256  std::vector<bool> available_neighbours = std::vector<bool>(8, true);
257  available_neighbours[0] = !on_north_edge;
258  available_neighbours[1] = !(on_north_edge || on_west_edge);
259  available_neighbours[2] = !on_west_edge;
260  available_neighbours[3] = !(on_south_edge || on_west_edge);
261  available_neighbours[4] = !on_south_edge;
262  available_neighbours[5] = !(on_south_edge || on_east_edge);
263  available_neighbours[6] = !on_east_edge;
264  available_neighbours[7] = !(on_north_edge || on_east_edge);
265 
266  // Using neighbour_indices_vector and available_neighbours, store the indices of all available neighbours to the set all_neighbours
267  for (unsigned i=0; i<8; i++)
268  {
269  if (available_neighbours[i])
270  {
271  assert(moore_neighbour_indices_vector[i] < nodes.size());
272  moore_neighbours[node_index].insert(moore_neighbour_indices_vector[i]);
273  }
274  }
275  break;
276  }
277  case 3:
278  {
279  assert(DIM ==3);
280  /*
281  * This stores the available neighbours using the following numbering:
282  * FRONT BACK
283  * 1-----0-----7 10-----9---16 19----18---25
284  * | | | | | | | | |
285  * | | | | | | | | |
286  * 2-----x-----6 11----8-----15 20----17---24
287  * | | | | | | | | |
288  * | | | | | | | | |
289  * 3-----4-----5 12----13----14 21---22----23
290  */
291 
292  // Create a vector of possible neighbouring node indices
293  std::vector<unsigned> moore_neighbour_indices_vector(26, node_index);
294  moore_neighbour_indices_vector[0] += numNodesAcross;
295  moore_neighbour_indices_vector[1] += numNodesAcross - 1;
296  moore_neighbour_indices_vector[2] -= 1;
297  moore_neighbour_indices_vector[3] -= numNodesAcross + 1;
298  moore_neighbour_indices_vector[4] -= numNodesAcross;
299  moore_neighbour_indices_vector[5] -= numNodesAcross - 1;
300  moore_neighbour_indices_vector[6] += 1;
301  moore_neighbour_indices_vector[7] += numNodesAcross + 1;
302  moore_neighbour_indices_vector[8] -= numNodesAcross*numNodesUp;
303  for (unsigned i=9; i<17; i++)
304  {
305  moore_neighbour_indices_vector[i] = moore_neighbour_indices_vector[i-9]-numNodesAcross*numNodesUp;
306  }
307  moore_neighbour_indices_vector[17] += numNodesAcross*numNodesUp;
308  for (unsigned i=18; i<26; i++)
309  {
310  moore_neighbour_indices_vector[i]=moore_neighbour_indices_vector[i-18]+numNodesAcross*numNodesUp;
311  }
312 
313  // Work out whether this node lies on any edge of the mesh
314  bool on_south_edge = (node_index%(numNodesAcross*numNodesUp)<numNodesAcross);
315  bool on_north_edge = (node_index%(numNodesAcross*numNodesUp)>(numNodesAcross*numNodesUp-numNodesAcross-1));
316  bool on_west_edge = (node_index%numNodesAcross == 0);
317  bool on_east_edge = (node_index%numNodesAcross == numNodesAcross - 1);
318  bool on_front_edge = (node_index < numNodesAcross*numNodesUp);
319  bool on_back_edge = (node_index > numNodesAcross*numNodesUp*(numNodesDeep-1)-1);
320 
321  if (isPeriodicInX)
322  {
323  if (on_west_edge)
324  {
325  moore_neighbour_indices_vector[1] = node_index + 2*numNodesAcross - 1;
326  moore_neighbour_indices_vector[2] = node_index + numNodesAcross - 1;
327  moore_neighbour_indices_vector[3] = node_index - 1;
328 
329  moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[1] - numNodesAcross*numNodesUp;
330  moore_neighbour_indices_vector[11] = moore_neighbour_indices_vector[2] - numNodesAcross*numNodesUp;
331  moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[3] - numNodesAcross*numNodesUp;
332 
333  moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[1] + numNodesAcross*numNodesUp;
334  moore_neighbour_indices_vector[20] = moore_neighbour_indices_vector[2] + numNodesAcross*numNodesUp;
335  moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[3] + numNodesAcross*numNodesUp;
336  }
337 
338  if (on_east_edge)
339  {
340  moore_neighbour_indices_vector[5] = node_index - 2*numNodesAcross + 1;
341  moore_neighbour_indices_vector[6] = node_index - numNodesAcross + 1;
342  moore_neighbour_indices_vector[7] = node_index + 1;
343 
344  moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[5] - numNodesAcross*numNodesUp;
345  moore_neighbour_indices_vector[15] = moore_neighbour_indices_vector[6] - numNodesAcross*numNodesUp;
346  moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[7] - numNodesAcross*numNodesUp;
347 
348  moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[5] + numNodesAcross*numNodesUp;
349  moore_neighbour_indices_vector[24] = moore_neighbour_indices_vector[6] + numNodesAcross*numNodesUp;
350  moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[7] + numNodesAcross*numNodesUp;
351  }
352  }
353 
354  if (isPeriodicInY)
355  {
356  if (on_north_edge)
357  {
358  moore_neighbour_indices_vector[0] = node_index - numNodesAcross*(numNodesUp-1);
359  moore_neighbour_indices_vector[1] = moore_neighbour_indices_vector[0] - 1;
360  moore_neighbour_indices_vector[7] = moore_neighbour_indices_vector[0] + 1;
361 
362  moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[1] - numNodesAcross*numNodesUp;
363  moore_neighbour_indices_vector[9] = moore_neighbour_indices_vector[0] - numNodesAcross*numNodesUp;
364  moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[7] - numNodesAcross*numNodesUp;
365 
366  moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[1] + numNodesAcross*numNodesUp;
367  moore_neighbour_indices_vector[18] = moore_neighbour_indices_vector[0] + numNodesAcross*numNodesUp;
368  moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[7] + numNodesAcross*numNodesUp;
369 
370  if (on_west_edge)
371  {
372  moore_neighbour_indices_vector[1] = moore_neighbour_indices_vector[1] + numNodesAcross;
373  moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[10] + numNodesAcross;
374  moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[19] + numNodesAcross;
375  }
376  if (on_east_edge)
377  {
378  moore_neighbour_indices_vector[7] = moore_neighbour_indices_vector[7] - numNodesAcross;
379  moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[16] - numNodesAcross;
380  moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[25] - numNodesAcross;
381  }
382  }
383 
384  if (on_south_edge)
385  {
386  moore_neighbour_indices_vector[4] = node_index + numNodesAcross*(numNodesUp-1);
387  moore_neighbour_indices_vector[3] = moore_neighbour_indices_vector[4] - 1;
388  moore_neighbour_indices_vector[5] = moore_neighbour_indices_vector[4] + 1;
389 
390  moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[3] - numNodesAcross*numNodesUp;
391  moore_neighbour_indices_vector[13] = moore_neighbour_indices_vector[4] - numNodesAcross*numNodesUp;
392  moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[5] - numNodesAcross*numNodesUp;
393 
394  moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[3] + numNodesAcross*numNodesUp;
395  moore_neighbour_indices_vector[22] = moore_neighbour_indices_vector[4] + numNodesAcross*numNodesUp;
396  moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[5] + numNodesAcross*numNodesUp;
397 
398  if (on_west_edge)
399  {
400  moore_neighbour_indices_vector[3] = moore_neighbour_indices_vector[3] + numNodesAcross;
401  moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[12] + numNodesAcross;
402  moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[21] + numNodesAcross;
403  }
404  if (on_east_edge)
405  {
406  moore_neighbour_indices_vector[5] = moore_neighbour_indices_vector[5] - numNodesAcross;
407  moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[14] - numNodesAcross;
408  moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[23] - numNodesAcross;
409  }
410  }
411  }
412 
413  if (isPeriodicInZ)
414  {
415  if (on_back_edge)
416  {
417  moore_neighbour_indices_vector[17] = node_index - numNodesAcross*numNodesUp*(numNodesDeep-1);
418  moore_neighbour_indices_vector[20] = moore_neighbour_indices_vector[17] - 1;
419  moore_neighbour_indices_vector[24] = moore_neighbour_indices_vector[17] + 1;
420 
421  moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[20] - numNodesAcross;
422  moore_neighbour_indices_vector[22] = moore_neighbour_indices_vector[17] - numNodesAcross;
423  moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[24] - numNodesAcross;
424 
425  moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[20] + numNodesAcross;
426  moore_neighbour_indices_vector[18] = moore_neighbour_indices_vector[17] + numNodesAcross;
427  moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[24] + numNodesAcross;
428 
429  if (on_west_edge)
430  {
431  moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[19] + numNodesAcross;
432  moore_neighbour_indices_vector[20] = moore_neighbour_indices_vector[20] + numNodesAcross;
433  moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[21] + numNodesAcross;
434  }
435  if (on_east_edge)
436  {
437  moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[23] - numNodesAcross;
438  moore_neighbour_indices_vector[24] = moore_neighbour_indices_vector[24] - numNodesAcross;
439  moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[25] - numNodesAcross;
440  }
441 
442  if (on_south_edge)
443  {
444  moore_neighbour_indices_vector[21] = moore_neighbour_indices_vector[21] + numNodesAcross*numNodesUp;
445  moore_neighbour_indices_vector[22] = moore_neighbour_indices_vector[22] + numNodesAcross*numNodesUp;
446  moore_neighbour_indices_vector[23] = moore_neighbour_indices_vector[23] + numNodesAcross*numNodesUp;
447  }
448 
449  if (on_north_edge)
450  {
451  moore_neighbour_indices_vector[18] = moore_neighbour_indices_vector[18] - numNodesAcross*numNodesUp;
452  moore_neighbour_indices_vector[19] = moore_neighbour_indices_vector[19] - numNodesAcross*numNodesUp;
453  moore_neighbour_indices_vector[25] = moore_neighbour_indices_vector[25] - numNodesAcross*numNodesUp;
454  }
455  }
456 
457  if (on_front_edge)
458  {
459  moore_neighbour_indices_vector[8] = node_index + numNodesAcross*numNodesUp*(numNodesDeep-1);
460  moore_neighbour_indices_vector[11] = moore_neighbour_indices_vector[8] - 1;
461  moore_neighbour_indices_vector[15] = moore_neighbour_indices_vector[8] + 1;
462 
463  moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[11] - numNodesAcross;
464  moore_neighbour_indices_vector[13] = moore_neighbour_indices_vector[8] - numNodesAcross;
465  moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[15] - numNodesAcross;
466 
467  moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[11] + numNodesAcross;
468  moore_neighbour_indices_vector[9] = moore_neighbour_indices_vector[8] + numNodesAcross;
469  moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[15] + numNodesAcross;
470 
471  if (on_west_edge)
472  {
473  moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[10] + numNodesAcross;
474  moore_neighbour_indices_vector[11] = moore_neighbour_indices_vector[11] + numNodesAcross;
475  moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[12] + numNodesAcross;
476  }
477  if (on_east_edge)
478  {
479  moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[14] - numNodesAcross;
480  moore_neighbour_indices_vector[15] = moore_neighbour_indices_vector[15] - numNodesAcross;
481  moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[16] - numNodesAcross;
482  }
483 
484  if (on_south_edge)
485  {
486  moore_neighbour_indices_vector[12] = moore_neighbour_indices_vector[12] + numNodesAcross*numNodesUp;
487  moore_neighbour_indices_vector[13] = moore_neighbour_indices_vector[13] + numNodesAcross*numNodesUp;
488  moore_neighbour_indices_vector[14] = moore_neighbour_indices_vector[14] + numNodesAcross*numNodesUp;
489  }
490 
491  if (on_north_edge)
492  {
493  moore_neighbour_indices_vector[9] = moore_neighbour_indices_vector[9] - numNodesAcross*numNodesUp;
494  moore_neighbour_indices_vector[10] = moore_neighbour_indices_vector[10] - numNodesAcross*numNodesUp;
495  moore_neighbour_indices_vector[16] = moore_neighbour_indices_vector[16] - numNodesAcross*numNodesUp;
496  }
497  }
498  }
499 
500  if (isPeriodicInX)
501  {
502  on_east_edge = false;
503  on_west_edge = false;
504  }
505  if (isPeriodicInY)
506  {
507  on_south_edge = false;
508  on_north_edge = false;
509  }
510  if (isPeriodicInZ)
511  {
512  on_front_edge = false;
513  on_back_edge = false;
514  }
515 
516  // Create a vector of booleans for which neighbours are available
517  // Use the order N, NW, W, SW, S, SE, E, NE
518  std::vector<bool> available_neighbours = std::vector<bool>(26, true);
519  available_neighbours[0] = !on_north_edge;
520  available_neighbours[1] = !(on_north_edge || on_west_edge);
521  available_neighbours[2] = !on_west_edge;
522  available_neighbours[3] = !(on_south_edge || on_west_edge);
523  available_neighbours[4] = !on_south_edge;
524  available_neighbours[5] = !(on_south_edge || on_east_edge);
525  available_neighbours[6] = !on_east_edge;
526  available_neighbours[7] = !(on_north_edge || on_east_edge);
527  available_neighbours[8] = !(on_front_edge);
528  for (unsigned i=9; i<17; i++)
529  {
530  available_neighbours[i] = (available_neighbours[i-9] && !(on_front_edge));
531  }
532  available_neighbours[17] = !(on_back_edge);
533  for (unsigned i=18; i<26; i++)
534  {
535  available_neighbours[i] = (available_neighbours[i-18] && !(on_back_edge));
536  }
537 
538  // Using neighbour_indices_vector and available_neighbours, store the indices of all available neighbours to the set all_neighbours
539  for (unsigned i=0; i<26; i++)
540  {
541  if (available_neighbours[i] && moore_neighbour_indices_vector[i] < numNodesAcross*numNodesUp*numNodesDeep)
542  {
543  assert(moore_neighbour_indices_vector[i] < nodes.size());
544  moore_neighbours[node_index].insert(moore_neighbour_indices_vector[i]);
545  }
546  }
547  break;
548  }
549  default:
550  NEVER_REACHED;
551  }
552 
553  // Clear the set of neighbouring node indices
554  von_neumann_neighbours[node_index].clear();
555 
556  switch (DIM)
557  {
558  case 2:
559  {
560  assert(DIM == 2);
561  /*
562  * This stores the available neighbours using the following numbering:
563  *
564  * 0
565  * |
566  * |
567  * 1-----x-----3
568  * |
569  * |
570  * 2
571  */
572  // Create a vector of possible neighbouring node indices
573  std::vector<unsigned> von_neumann_neighbour_indices_vector(4, node_index);
574  von_neumann_neighbour_indices_vector[0] += numNodesAcross;
575  von_neumann_neighbour_indices_vector[1] -= 1;
576  von_neumann_neighbour_indices_vector[2] -= numNodesAcross;
577  von_neumann_neighbour_indices_vector[3] += 1;
578 
579  // Work out whether this node lies on any edge of the mesh
580  bool on_south_edge = (node_index < numNodesAcross);
581  bool on_north_edge = ((int)node_index > (int)numNodesAcross*((int)numNodesUp - 1) - 1);
582  bool on_west_edge = (node_index%numNodesAcross == 0);
583  bool on_east_edge = (node_index%numNodesAcross == numNodesAcross - 1);
584 
585  if (isPeriodicInX)
586  {
587  if (on_west_edge)
588  {
589  von_neumann_neighbour_indices_vector[1] = node_index + numNodesAcross - 1;
590  on_west_edge = false;
591  }
592 
593  if (on_east_edge)
594  {
595  von_neumann_neighbour_indices_vector[3] = node_index - numNodesAcross + 1;
596  on_east_edge = false;
597  }
598  }
599 
600  if (isPeriodicInY)
601  {
602  if (on_north_edge)
603  {
604  von_neumann_neighbour_indices_vector[0] = node_index - numNodesAcross*(numNodesUp-1);
605  on_north_edge = false;
606  }
607 
608  if (on_south_edge)
609  {
610  von_neumann_neighbour_indices_vector[2] = node_index + numNodesAcross*(numNodesUp-1);
611  on_south_edge = false;
612  }
613  }
614 
615  // Create a vector of booleans for which neighbours are available
616  // Use the order N, W, S, E
617  std::vector<bool> available_neighbours = std::vector<bool>(2*DIM, true);
618  available_neighbours[0] = !on_north_edge;
619  available_neighbours[1] = !on_west_edge;
620  available_neighbours[2] = !on_south_edge;
621  available_neighbours[3] = !on_east_edge;
622 
623  // Using von_neumann_neighbour_indices_vector and available_neighbours, store the indices of all available neighbours to the set all_neighbours
624  for (unsigned i=0; i<4; i++)
625  {
626  if (available_neighbours[i])
627  {
628  assert(von_neumann_neighbour_indices_vector[i] < nodes.size());
629  von_neumann_neighbours[node_index].insert(von_neumann_neighbour_indices_vector[i]);
630  }
631  }
632  break;
633  }
634  case 3:
635  {
636  assert(DIM == 3);
637 
638  /*
639  * This stores the available neighbours using the following numbering:
640  *
641  * 0 5
642  * | /
643  * |/
644  * 1-----x-----3
645  * / |
646  * / |
647  * 4 2
648  */
649  // Create a vector of possible neighbouring node indices
650  std::vector<unsigned> von_neumann_neighbour_indices_vector(6, node_index);
651  von_neumann_neighbour_indices_vector[0] += numNodesAcross;
652  von_neumann_neighbour_indices_vector[1] -= 1;
653  von_neumann_neighbour_indices_vector[2] -= numNodesAcross;
654  von_neumann_neighbour_indices_vector[3] += 1;
655  von_neumann_neighbour_indices_vector[4] -= numNodesAcross*numNodesUp;
656  von_neumann_neighbour_indices_vector[5] += numNodesAcross*numNodesUp;
657 
658  // Work out whether this node lies on any edge of the mesh
659  bool on_south_edge = (node_index%(numNodesAcross*numNodesUp)<numNodesAcross);
660  bool on_north_edge = (node_index%(numNodesAcross*numNodesUp)>(numNodesAcross*numNodesUp-numNodesAcross-1));
661  bool on_west_edge = (node_index%numNodesAcross== 0);
662  bool on_east_edge = (node_index%numNodesAcross == numNodesAcross - 1);
663  bool on_front_edge = (node_index < numNodesAcross*numNodesUp);
664  bool on_back_edge = (node_index > numNodesAcross*numNodesUp*(numNodesDeep-1)-1);
665 
666  if (isPeriodicInX)
667  {
668  if (on_west_edge)
669  {
670  von_neumann_neighbour_indices_vector[1] = node_index + numNodesAcross - 1;
671  on_west_edge = false;
672  }
673 
674  if (on_east_edge)
675  {
676  von_neumann_neighbour_indices_vector[3] = node_index - numNodesAcross + 1;
677  on_east_edge = false;
678  }
679  }
680 
681  if (isPeriodicInY)
682  {
683  if (on_north_edge)
684  {
685  von_neumann_neighbour_indices_vector[0] = node_index - numNodesAcross*(numNodesUp-1);
686  on_north_edge = false;
687  }
688 
689  if (on_south_edge)
690  {
691  von_neumann_neighbour_indices_vector[2] = node_index + numNodesAcross*(numNodesUp-1);
692  on_south_edge = false;
693  }
694  }
695 
696  if (isPeriodicInZ)
697  {
698  if (on_front_edge)
699  {
700  von_neumann_neighbour_indices_vector[4] = node_index + numNodesAcross*numNodesUp*(numNodesDeep-1);
701  on_front_edge = false;
702  }
703 
704  if (on_back_edge)
705  {
706  von_neumann_neighbour_indices_vector[5] = node_index - numNodesAcross*numNodesUp*(numNodesDeep-1);
707  on_back_edge = false;
708  }
709  }
710 
711  // Create a vector of booleans for which neighbours are available
712  // Use the order N, W, S, E, F, B
713  std::vector<bool> available_neighbours = std::vector<bool>(2*DIM, true);
714  available_neighbours[0] = !on_north_edge;
715  available_neighbours[1] = !on_west_edge;
716  available_neighbours[2] = !on_south_edge;
717  available_neighbours[3] = !on_east_edge;
718  available_neighbours[4] = !on_front_edge;
719  available_neighbours[5] = !on_back_edge;
720 
721  // Using von_neumann_neighbour_indices_vector and available_neighbours, store the indices of all available neighbours to the set all_neighbours
722  for (unsigned i=0; i<6; i++)
723  {
724  if (available_neighbours[i] && von_neumann_neighbour_indices_vector[i]<numNodesAcross*numNodesUp*numNodesDeep)
725  {
726  assert(von_neumann_neighbour_indices_vector[i] < nodes.size());
727  von_neumann_neighbours[node_index].insert(von_neumann_neighbour_indices_vector[i]);
728  }
729  }
730  break;
731  }
732  default:
733  NEVER_REACHED;
734  }
735  }
736 
737  mpMesh = new PottsMesh<DIM>(nodes, elements, von_neumann_neighbours, moore_neighbours);
738 }
739 
740 template<unsigned DIM>
741 PottsMeshGenerator<DIM>::~PottsMeshGenerator()
742 {
743  delete mpMesh;
744 }
745 
746 template<unsigned DIM>
747 PottsMesh<DIM>* PottsMeshGenerator<DIM>::GetMesh()
748 {
749  return mpMesh;
750 }
751 
752 // Explicit instantiation
753 template class PottsMeshGenerator<1>;
754 template class PottsMeshGenerator<2>;
755 template class PottsMeshGenerator<3>;
Node
Definition: Node.hpp:58
NEVER_REACHED
#define NEVER_REACHED
Definition: Exception.hpp:206
PottsMeshGenerator::GetMesh
virtual PottsMesh< DIM > * GetMesh()
Definition: PottsMeshGenerator.cpp:747