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