rectangular_quadmesh.template.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#ifndef OOMPH_RECTANGULAR_QUADMESH_TEMPLATE_CC
27#define OOMPH_RECTANGULAR_QUADMESH_TEMPLATE_CC
28
29
30// OOMPH-LIB Headers
31#include "rectangular_quadmesh.template.h"
32
33
34namespace oomph
35{
36 //===========================================================================
37 /// Generic mesh construction. This function contains the "guts" of the
38 /// mesh generation process, including all the tedious loops, counting
39 /// and spacing functions. The function should be called in all constuctors
40 /// of any derived classes.
41 //===========================================================================
42 template<class ELEMENT>
43 void RectangularQuadMesh<ELEMENT>::build_mesh(TimeStepper* time_stepper_pt)
44 {
45 // Mesh can only be built with 2D Qelements.
46 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
47
48 // Set the number of boundaries
49 set_nboundary(4);
50
51 // Allocate the store for the elements
52 Element_pt.resize(Nx * Ny);
53 // Allocate the memory for the first element
54 Element_pt[0] = new ELEMENT;
55
56 // Read out the number of linear points in the element
57 Np = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
58
59 // Can now allocate the store for the nodes
60 Node_pt.resize((1 + (Np - 1) * Nx) * (1 + (Np - 1) * Ny));
61
62 // Set up geometrical data
63 unsigned long node_count = 0;
64
65 // Now assign the topology
66 // Boundaries are numbered 0 1 2 3 from the bottom proceeding anticlockwise
67 // Pinned value are denoted by an integer value 1
68 // Thus if a node is on two boundaries, ORing the values of the
69 // boundary conditions will give the most restrictive case (pinning)
70
71 // FIRST ELEMENT (lower left corner)
72
73 // Set the corner node
74 // Allocate memory for the node
75 Node_pt[node_count] =
76 finite_element_pt(0)->construct_boundary_node(0, time_stepper_pt);
77
78 // Set the position of the node
79 Node_pt[node_count]->x(0) = x_spacing_function(0, 0, 0, 0);
80 Node_pt[node_count]->x(1) = y_spacing_function(0, 0, 0, 0);
81
82 // Push the node back onto boundaries
83 add_boundary_node(0, Node_pt[node_count]);
84 add_boundary_node(3, Node_pt[node_count]);
85
86 // Increment the node number
87 node_count++;
88
89 // Loop over the other nodes in the first row
90 for (unsigned l2 = 1; l2 < Np; l2++)
91 {
92 // Allocate memory for the nodes
93 Node_pt[node_count] =
94 finite_element_pt(0)->construct_boundary_node(l2, time_stepper_pt);
95
96 // Set the position of the node
97 Node_pt[node_count]->x(0) = x_spacing_function(0, l2, 0, 0);
98 Node_pt[node_count]->x(1) = y_spacing_function(0, l2, 0, 0);
99
100 // Push the node back onto boundaries
101 add_boundary_node(0, Node_pt[node_count]);
102
103 // If we only have one column then the RHS node is on the right boundary
104 if ((Nx == 1) && (l2 == (Np - 1)))
105 {
106 add_boundary_node(1, Node_pt[node_count]);
107 }
108
109 // Increment the node number
110 node_count++;
111 }
112
113 // Loop over the other node columns
114 for (unsigned l1 = 1; l1 < Np; l1++)
115 {
116 // Allocate memory for the nodes
117 Node_pt[node_count] =
118 finite_element_pt(0)->construct_boundary_node(l1 * Np, time_stepper_pt);
119
120 // Set the position of the node
121 Node_pt[node_count]->x(0) = x_spacing_function(0, 0, 0, l1);
122 Node_pt[node_count]->x(1) = y_spacing_function(0, 0, 0, l1);
123
124 // Push the node back onto boundaries
125 add_boundary_node(3, Node_pt[node_count]);
126
127 // If we only have one row, then the top node is on the top boundary
128 if ((Ny == 1) && (l1 == (Np - 1)))
129 {
130 add_boundary_node(2, Node_pt[node_count]);
131 }
132
133 // Increment the node number
134 node_count++;
135
136 // Loop over the other nodes in the row
137 for (unsigned l2 = 1; l2 < Np; l2++)
138 {
139 // Allocate the memory for the node
140 // If it lies on a boundary make a boundary node
141 if (((Nx == 1) && (l2 == (Np - 1))) || ((Ny == 1) && (l1 == (Np - 1))))
142 {
143 Node_pt[node_count] = finite_element_pt(0)->construct_boundary_node(
144 l1 * Np + l2, time_stepper_pt);
145 }
146 // Otherwise its a normal node
147 else
148 {
149 Node_pt[node_count] =
150 finite_element_pt(0)->construct_node(l1 * Np + l2, time_stepper_pt);
151 }
152
153 // Set the position of the node
154 Node_pt[node_count]->x(0) = x_spacing_function(0, l2, 0, l1);
155 Node_pt[node_count]->x(1) = y_spacing_function(0, l2, 0, l1);
156
157 // If we only have one column then the RHS node is on the right boundary
158 if ((Nx == 1) && l2 == (Np - 1))
159 {
160 add_boundary_node(1, Node_pt[node_count]);
161 }
162 // If we only have one row, then the top node is on the top boundary
163 if ((Ny == 1) && (l1 == (Np - 1)))
164 {
165 add_boundary_node(2, Node_pt[node_count]);
166 }
167
168 // Increment the node number
169 node_count++;
170 }
171 }
172 // END OF FIRST ELEMENT
173
174 // CENTRE OF FIRST ROW OF ELEMENTS
175 // Now loop over the first row of elements, apart from final element
176 for (unsigned j = 1; j < (Nx - 1); j++)
177 {
178 // Allocate memory for new element
179 Element_pt[j] = new ELEMENT;
180 // Do first row of nodes
181 // First column of nodes is same as neighbouring element
182 finite_element_pt(j)->node_pt(0) =
183 finite_element_pt(j - 1)->node_pt((Np - 1));
184 // New nodes for other columns
185 for (unsigned l2 = 1; l2 < Np; l2++)
186 {
187 // Allocate memory for the nodes
188 Node_pt[node_count] =
189 finite_element_pt(j)->construct_boundary_node(l2, time_stepper_pt);
190
191 // Set the position of the node
192 Node_pt[node_count]->x(0) = x_spacing_function(j, l2, 0, 0);
193 Node_pt[node_count]->x(1) = y_spacing_function(j, l2, 0, 0);
194
195 // Push the node back onto boundaries
196 add_boundary_node(0, Node_pt[node_count]);
197
198 // Increment the node number
199 node_count++;
200 }
201
202 // Do the rest of the nodes
203 for (unsigned l1 = 1; l1 < Np; l1++)
204 {
205 // First column of nodes is same as neighbouring element
206 finite_element_pt(j)->node_pt(l1 * Np) =
207 finite_element_pt(j - 1)->node_pt(l1 * Np + (Np - 1));
208 // New nodes for other columns
209 for (unsigned l2 = 1; l2 < Np; l2++)
210 {
211 // Allocate memory for the nodes
212 // If we only have one row, this node could be on the boundary
213 if ((Ny == 1) && (l1 == (Np - 1)))
214 {
215 Node_pt[node_count] = finite_element_pt(j)->construct_boundary_node(
216 l1 * Np + l2, time_stepper_pt);
217 }
218 // Otherwise create a normal node
219 else
220 {
221 Node_pt[node_count] = finite_element_pt(j)->construct_node(
222 l1 * Np + l2, time_stepper_pt);
223 }
224
225 // Set the position of the node
226 Node_pt[node_count]->x(0) = x_spacing_function(j, l2, 0, l1);
227 Node_pt[node_count]->x(1) = y_spacing_function(j, l2, 0, l1);
228
229 // If we only have one row, then the top node is on the top boundary
230 if ((Ny == 1) && (l1 == (Np - 1)))
231 {
232 add_boundary_node(2, Node_pt[node_count]);
233 }
234
235 // Increment the node number
236 node_count++;
237 }
238 }
239 }
240 // END OF CENTRE OF FIRST ROW OF ELEMENTS
241
242 // FINAL ELEMENT IN FIRST ROW (lower right corner)
243 // Only allocate if there is more than one element in the row
244 if (Nx > 1)
245 {
246 // Allocate memory for element
247 Element_pt[Nx - 1] = new ELEMENT;
248
249 // First column of nodes is same as neighbouring element
250 finite_element_pt(Nx - 1)->node_pt(0) =
251 finite_element_pt(Nx - 2)->node_pt(Np - 1);
252
253 // New middle nodes
254 for (unsigned l2 = 1; l2 < (Np - 1); l2++)
255 {
256 // Allocate memory for node
257 Node_pt[node_count] =
258 finite_element_pt(Nx - 1)->construct_boundary_node(l2,
259 time_stepper_pt);
260
261 // Set the position of the node
262 Node_pt[node_count]->x(0) = x_spacing_function(Nx - 1, l2, 0, 0);
263 Node_pt[node_count]->x(1) = y_spacing_function(Nx - 1, l2, 0, 0);
264
265 // Push the node back onto boundaries
266 add_boundary_node(0, Node_pt[node_count]);
267
268 // Increment the node number
269 node_count++;
270 }
271
272 // Allocate memory for a new boundary node
273 Node_pt[node_count] = finite_element_pt(Nx - 1)->construct_boundary_node(
274 Np - 1, time_stepper_pt);
275
276 // If required make it periodic from the node on the other side
277 if (Xperiodic == true)
278 {
279 Node_pt[node_count]->make_periodic(finite_element_pt(0)->node_pt(0));
280 }
281
282 // Set the position of the node
283 Node_pt[node_count]->x(0) = x_spacing_function(Nx - 1, Np - 1, 0, 0);
284 Node_pt[node_count]->x(1) = y_spacing_function(Nx - 1, Np - 1, 0, 0);
285
286 // Push the node back onto boundaries
287 add_boundary_node(0, Node_pt[node_count]);
288 add_boundary_node(1, Node_pt[node_count]);
289
290 // Increment the node number
291 node_count++;
292
293 // Do the rest of the nodes
294 for (unsigned l1 = 1; l1 < Np; l1++)
295 {
296 // First column of nodes is same as neighbouring element
297 finite_element_pt(Nx - 1)->node_pt(l1 * Np) =
298 finite_element_pt(Nx - 2)->node_pt(l1 * Np + (Np - 1));
299
300 // New node for middle column
301 for (unsigned l2 = 1; l2 < (Np - 1); l2++)
302 {
303 // Allocate memory for node
304 // If it lies on the boundary, create a boundary node
305 if ((Ny == 1) && (l1 == (Np - 1)))
306 {
307 Node_pt[node_count] =
308 finite_element_pt(Nx - 1)->construct_boundary_node(
309 l1 * Np + l2, time_stepper_pt);
310 }
311 // Otherwise create a normal node
312 else
313 {
314 Node_pt[node_count] = finite_element_pt(Nx - 1)->construct_node(
315 l1 * Np + l2, time_stepper_pt);
316 }
317
318 // Set the position of the node
319 Node_pt[node_count]->x(0) = x_spacing_function(Nx - 1, l2, 0, l1);
320 Node_pt[node_count]->x(1) = y_spacing_function(Nx - 1, l2, 0, l1);
321
322 // If we only have one row, then the top node is on the top boundary
323 if ((Ny == 1) && (l1 == (Np - 1)))
324 {
325 add_boundary_node(2, Node_pt[node_count]);
326 }
327
328 // Increment the node number
329 node_count++;
330 }
331
332 // Allocate memory for a new boundary node
333 Node_pt[node_count] =
334 finite_element_pt(Nx - 1)->construct_boundary_node(l1 * Np + (Np - 1),
335 time_stepper_pt);
336
337 // If required make it periodic from the corresponding node on the other
338 // side
339 if (Xperiodic == true)
340 {
341 Node_pt[node_count]->make_periodic(
342 finite_element_pt(0)->node_pt(l1 * Np));
343 }
344
345 // Set the position of the node
346 Node_pt[node_count]->x(0) = x_spacing_function(Nx - 1, Np - 1, 0, l1);
347 Node_pt[node_count]->x(1) = y_spacing_function(Nx - 1, Np - 1, 0, l1);
348
349 // Push the node back onto boundaries
350 add_boundary_node(1, Node_pt[node_count]);
351
352 // If we only have one row, then the top node is on the top boundary
353 if ((Ny == 1) && (l1 == (Np - 1)))
354 {
355 add_boundary_node(2, Node_pt[node_count]);
356 }
357
358 // Increment the node number
359 node_count++;
360 }
361 }
362 // END OF FIRST ROW OF ELEMENTS
363
364 // ALL CENTRAL ELEMENT ROWS
365 // Loop over remaining element rows
366 for (unsigned i = 1; i < (Ny - 1); i++)
367 {
368 // Set the first element in the row
369 // Allocate memory for element
370 Element_pt[Nx * i] = new ELEMENT;
371
372 // The first row of nodes is copied from the element below
373 for (unsigned l2 = 0; l2 < Np; l2++)
374 {
375 finite_element_pt(Nx * i)->node_pt(l2) =
376 finite_element_pt(Nx * (i - 1))->node_pt((Np - 1) * Np + l2);
377 }
378
379 // Other rows are new nodes
380 for (unsigned l1 = 1; l1 < Np; l1++)
381 {
382 // First column of nodes
383 // Allocate memory for node
384 Node_pt[node_count] =
385 finite_element_pt(Nx * i)->construct_boundary_node(l1 * Np,
386 time_stepper_pt);
387
388 // Set the position of the node
389 Node_pt[node_count]->x(0) = x_spacing_function(0, 0, i, l1);
390 Node_pt[node_count]->x(1) = y_spacing_function(0, 0, i, l1);
391
392 // Push the node back onto boundaries
393 add_boundary_node(3, Node_pt[node_count]);
394
395 // Increment the node number
396 node_count++;
397
398 // Now do the other columns
399 for (unsigned l2 = 1; l2 < Np; l2++)
400 {
401 // Allocate memory for node
402 // If we only have one column, the node could be on the boundary
403 if ((Nx == 1) && (l2 == (Np - 1)))
404 {
405 Node_pt[node_count] =
406 finite_element_pt(Nx * i)->construct_boundary_node(
407 l1 * Np + l2, time_stepper_pt);
408 }
409 else
410 {
411 Node_pt[node_count] = finite_element_pt(Nx * i)->construct_node(
412 l1 * Np + l2, time_stepper_pt);
413 }
414
415 // Set the position of the node
416 Node_pt[node_count]->x(0) = x_spacing_function(0, l2, i, l1);
417 Node_pt[node_count]->x(1) = y_spacing_function(0, l2, i, l1);
418
419 // If we only have one column then the RHS node is on the
420 // right boundary
421 if ((Nx == 1) && (l2 == (Np - 1)))
422 {
423 add_boundary_node(1, Node_pt[node_count]);
424 }
425
426 // Increment the node number
427 node_count++;
428 }
429 }
430
431 // Now loop over the rest of the elements in the row, apart from the last
432 for (unsigned j = 1; j < (Nx - 1); j++)
433 {
434 // Allocate memory for new element
435 Element_pt[Nx * i + j] = new ELEMENT;
436 // The first row is copied from the lower element
437 for (unsigned l2 = 0; l2 < Np; l2++)
438 {
439 finite_element_pt(Nx * i + j)->node_pt(l2) =
440 finite_element_pt(Nx * (i - 1) + j)->node_pt((Np - 1) * Np + l2);
441 }
442
443 for (unsigned l1 = 1; l1 < Np; l1++)
444 {
445 // First column is same as neighbouring element
446 finite_element_pt(Nx * i + j)->node_pt(l1 * Np) =
447 finite_element_pt(Nx * i + (j - 1))->node_pt(l1 * Np + (Np - 1));
448 // New nodes for other columns
449 for (unsigned l2 = 1; l2 < Np; l2++)
450 {
451 // Allocate memory for the nodes
452 Node_pt[node_count] =
453 finite_element_pt(Nx * i + j)
454 ->construct_node(l1 * Np + l2, time_stepper_pt);
455
456 // Set the position of the node
457 Node_pt[node_count]->x(0) = x_spacing_function(j, l2, i, l1);
458 Node_pt[node_count]->x(1) = y_spacing_function(j, l2, i, l1);
459
460 // Increment the node number
461 node_count++;
462 }
463 }
464 } // End of loop over elements in row
465
466 // Do final element in row
467 // Only if there is more than one column
468 if (Nx > 1)
469 {
470 // Allocate memory for element
471 Element_pt[Nx * i + Nx - 1] = new ELEMENT;
472 // The first row is copied from the lower element
473 for (unsigned l2 = 0; l2 < Np; l2++)
474 {
475 finite_element_pt(Nx * i + Nx - 1)->node_pt(l2) =
476 finite_element_pt(Nx * (i - 1) + Nx - 1)
477 ->node_pt((Np - 1) * Np + l2);
478 }
479
480 for (unsigned l1 = 1; l1 < Np; l1++)
481 {
482 // First column is same as neighbouring element
483 finite_element_pt(Nx * i + Nx - 1)->node_pt(l1 * Np) =
484 finite_element_pt(Nx * i + Nx - 2)->node_pt(l1 * Np + (Np - 1));
485
486 // Middle nodes
487 for (unsigned l2 = 1; l2 < (Np - 1); l2++)
488 {
489 // Allocate memory for node
490 Node_pt[node_count] =
491 finite_element_pt(Nx * i + Nx - 1)
492 ->construct_node(l1 * Np + l2, time_stepper_pt);
493
494 // Set the position of the node
495 Node_pt[node_count]->x(0) = x_spacing_function(Nx - 1, l2, i, l1);
496 Node_pt[node_count]->x(1) = y_spacing_function(Nx - 1, l2, i, l1);
497
498 // Increment the node number
499 node_count++;
500 }
501
502 // Allocate memory for a new boundary node
503 Node_pt[node_count] =
504 finite_element_pt(Nx * i + Nx - 1)
505 ->construct_boundary_node(l1 * Np + (Np - 1), time_stepper_pt);
506
507 // If required make it periodic from the corresponding node on
508 // the other side
509 if (Xperiodic == true)
510 {
511 Node_pt[node_count]->make_periodic(
512 finite_element_pt(Nx * i)->node_pt(l1 * Np));
513 }
514
515 // Set the position of the node
516 Node_pt[node_count]->x(0) = x_spacing_function(Nx - 1, Np - 1, i, l1);
517 Node_pt[node_count]->x(1) = y_spacing_function(Nx - 1, Np - 1, i, l1);
518
519 // Push the node back onto boundaries
520 add_boundary_node(1, Node_pt[node_count]);
521
522 // Increment the node number
523 node_count++;
524 } // End of loop over rows of nodes in the element
525 } // End of if more than one column
526 } // End of loop over rows of elements
527
528 // END OF LOOP OVER CENTRAL ELEMENT ROWS
529
530 // FINAL ELEMENT ROW
531 // ONLY NECESSARY IF THERE IS MORE THAN ONE ROW
532 if (Ny > 1)
533 {
534 // FIRST ELEMENT IN UPPER ROW (upper left corner)
535 // Allocate memory for element
536 Element_pt[Nx * (Ny - 1)] = new ELEMENT;
537 // The first row of nodes is copied from the element below
538 for (unsigned l2 = 0; l2 < Np; l2++)
539 {
540 finite_element_pt(Nx * (Ny - 1))->node_pt(l2) =
541 finite_element_pt(Nx * (Ny - 2))->node_pt((Np - 1) * Np + l2);
542 }
543
544 // Second row of nodes
545 // First column of nodes
546 for (unsigned l1 = 1; l1 < (Np - 1); l1++)
547 {
548 // Allocate memory for node
549 Node_pt[node_count] =
550 finite_element_pt(Nx * (Ny - 1))
551 ->construct_boundary_node(Np * l1, time_stepper_pt);
552
553 // Set the position of the node
554 Node_pt[node_count]->x(0) = x_spacing_function(0, 0, Ny - 1, l1);
555 Node_pt[node_count]->x(1) = y_spacing_function(0, 0, Ny - 1, l1);
556
557 // Push the node back onto boundaries
558 add_boundary_node(3, Node_pt[node_count]);
559
560 // Increment the node number
561 node_count++;
562
563 // Now do the other columns
564 for (unsigned l2 = 1; l2 < Np; l2++)
565 {
566 // Allocate memory for node
567 if ((Nx == 1) && (l2 == Np - 1))
568 {
569 Node_pt[node_count] =
570 finite_element_pt(Nx * (Ny - 1))
571 ->construct_boundary_node(Np * l1 + l2, time_stepper_pt);
572 }
573 else
574 {
575 Node_pt[node_count] =
576 finite_element_pt(Nx * (Ny - 1))
577 ->construct_node(Np * l1 + l2, time_stepper_pt);
578 }
579
580 // Set the position of the node
581 Node_pt[node_count]->x(0) = x_spacing_function(0, l2, Ny - 1, l1);
582 Node_pt[node_count]->x(1) = y_spacing_function(0, l2, Ny - 1, l1);
583
584 // Push the node back onto boundaries
585 if ((Nx == 1) && (l2 == Np - 1))
586 {
587 add_boundary_node(1, Node_pt[node_count]);
588 }
589
590 // Increment the node number
591 node_count++;
592 }
593 }
594
595 // Final row of nodes
596 // First column of nodes
597 // Top left node
598 // Allocate memory for node
599 Node_pt[node_count] =
600 finite_element_pt(Nx * (Ny - 1))
601 ->construct_boundary_node(Np * (Np - 1), time_stepper_pt);
602
603 // Set the position of the node
604 Node_pt[node_count]->x(0) = x_spacing_function(0, 0, Ny - 1, Np - 1);
605 Node_pt[node_count]->x(1) = y_spacing_function(0, 0, Ny - 1, Np - 1);
606
607 // Push the node back onto boundaries
608 add_boundary_node(2, Node_pt[node_count]);
609 add_boundary_node(3, Node_pt[node_count]);
610
611 // Increment the node number
612 node_count++;
613
614 // Now do the other columns
615 for (unsigned l2 = 1; l2 < Np; l2++)
616 {
617 // Allocate memory for the node
618 Node_pt[node_count] =
619 finite_element_pt(Nx * (Ny - 1))
620 ->construct_boundary_node(Np * (Np - 1) + l2, time_stepper_pt);
621
622 // Set the position of the node
623 Node_pt[node_count]->x(0) = x_spacing_function(0, l2, Ny - 1, Np - 1);
624 Node_pt[node_count]->x(1) = y_spacing_function(0, l2, Ny - 1, Np - 1);
625
626 // Push the node back onto boundaries
627 add_boundary_node(2, Node_pt[node_count]);
628
629 // Push the node back onto boundaries
630 if ((Nx == 1) && (l2 == Np - 1))
631 {
632 add_boundary_node(1, Node_pt[node_count]);
633 }
634
635 // Increment the node number
636 node_count++;
637 }
638
639 // Now loop over the rest of the elements in the row, apart from the last
640 for (unsigned j = 1; j < (Nx - 1); j++)
641 {
642 // Allocate memory for element
643 Element_pt[Nx * (Ny - 1) + j] = new ELEMENT;
644 // The first row is copied from the lower element
645 for (unsigned l2 = 0; l2 < Np; l2++)
646 {
647 finite_element_pt(Nx * (Ny - 1) + j)->node_pt(l2) =
648 finite_element_pt(Nx * (Ny - 2) + j)->node_pt((Np - 1) * Np + l2);
649 }
650
651 // Second rows
652 for (unsigned l1 = 1; l1 < (Np - 1); l1++)
653 {
654 // First column is same as neighbouring element
655 finite_element_pt(Nx * (Ny - 1) + j)->node_pt(Np * l1) =
656 finite_element_pt(Nx * (Ny - 1) + (j - 1))
657 ->node_pt(Np * l1 + (Np - 1));
658
659 // New nodes for other columns
660 for (unsigned l2 = 1; l2 < Np; l2++)
661 {
662 // Allocate memory for the node
663 Node_pt[node_count] =
664 finite_element_pt(Nx * (Ny - 1) + j)
665 ->construct_node(Np * l1 + l2, time_stepper_pt);
666
667 // Set the position of the node
668 Node_pt[node_count]->x(0) = x_spacing_function(j, l2, Ny - 1, l1);
669 Node_pt[node_count]->x(1) = y_spacing_function(j, l2, Ny - 1, l1);
670
671 // Increment the node number
672 node_count++;
673 }
674 }
675
676 // Top row
677 // First column is same as neighbouring element
678 finite_element_pt(Nx * (Ny - 1) + j)->node_pt(Np * (Np - 1)) =
679 finite_element_pt(Nx * (Ny - 1) + (j - 1))
680 ->node_pt(Np * (Np - 1) + (Np - 1));
681 // New nodes for other columns
682 for (unsigned l2 = 1; l2 < Np; l2++)
683 {
684 // Allocate memory for node
685 Node_pt[node_count] =
686 finite_element_pt(Nx * (Ny - 1) + j)
687 ->construct_boundary_node(Np * (Np - 1) + l2, time_stepper_pt);
688
689 // Set the position of the node
690 Node_pt[node_count]->x(0) = x_spacing_function(j, l2, Ny - 1, Np - 1);
691 Node_pt[node_count]->x(1) = y_spacing_function(j, l2, Ny - 1, Np - 1);
692
693 // Push the node back onto boundaries
694 add_boundary_node(2, Node_pt[node_count]);
695
696 // Increment the node number
697 node_count++;
698 }
699 } // End of loop over central elements in row
700
701 // FINAL ELEMENT IN ROW (upper right corner)
702 // Only if there is more than one column
703 if (Nx > 1)
704 {
705 // Allocate memory for element
706 Element_pt[Nx * (Ny - 1) + Nx - 1] = new ELEMENT;
707 // The first row is copied from the lower element
708 for (unsigned l2 = 0; l2 < Np; l2++)
709 {
710 finite_element_pt(Nx * (Ny - 1) + Nx - 1)->node_pt(l2) =
711 finite_element_pt(Nx * (Ny - 2) + Nx - 1)
712 ->node_pt((Np - 1) * Np + l2);
713 }
714
715 // Second rows
716 for (unsigned l1 = 1; l1 < (Np - 1); l1++)
717 {
718 // First column is same as neighbouring element
719 finite_element_pt(Nx * (Ny - 1) + Nx - 1)->node_pt(Np * l1) =
720 finite_element_pt(Nx * (Ny - 1) + Nx - 2)
721 ->node_pt(Np * l1 + (Np - 1));
722
723 // Middle nodes
724 for (unsigned l2 = 1; l2 < (Np - 1); l2++)
725 {
726 // Allocate memory for node
727 Node_pt[node_count] =
728 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
729 ->construct_node(Np * l1 + l2, time_stepper_pt);
730
731 // Set the position of the node
732 Node_pt[node_count]->x(0) =
733 x_spacing_function(Nx - 1, l2, Ny - 1, l1);
734 Node_pt[node_count]->x(1) =
735 y_spacing_function(Nx - 1, l2, Ny - 1, l1);
736
737 // Increment the node number
738 node_count++;
739 }
740
741 // Final node
742 // Allocate new memory for a boundary node
743 Node_pt[node_count] =
744 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
745 ->construct_boundary_node(Np * l1 + (Np - 1), time_stepper_pt);
746
747 // If required make it periodic
748 if (Xperiodic == true)
749 {
750 Node_pt[node_count]->make_periodic(
751 finite_element_pt(Nx * (Ny - 1))->node_pt(Np * l1));
752 }
753
754 // Set the position of the node
755 Node_pt[node_count]->x(0) =
756 x_spacing_function(Nx - 1, Np - 1, Ny - 1, l1);
757 Node_pt[node_count]->x(1) =
758 y_spacing_function(Nx - 1, Np - 1, Ny - 1, l1);
759
760 // Push the node back onto boundaries
761 add_boundary_node(1, Node_pt[node_count]);
762
763 // Increment the node number
764 node_count++;
765
766 } // End of loop over middle rows
767
768 // Final row
769 // First column is same as neighbouring element
770 finite_element_pt(Nx * (Ny - 1) + Nx - 1)->node_pt(Np * (Np - 1)) =
771 finite_element_pt(Nx * (Ny - 1) + Nx - 2)
772 ->node_pt(Np * (Np - 1) + (Np - 1));
773
774 // Middle nodes
775 for (unsigned l2 = 1; l2 < (Np - 1); l2++)
776 {
777 // Allocate memory for node
778 Node_pt[node_count] =
779 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
780 ->construct_boundary_node(Np * (Np - 1) + l2, time_stepper_pt);
781
782 // Set the position of the node
783 Node_pt[node_count]->x(0) =
784 x_spacing_function(Nx - 1, l2, Ny - 1, Np - 1);
785 Node_pt[node_count]->x(1) =
786 y_spacing_function(Nx - 1, l2, Ny - 1, Np - 1);
787
788 // Push the node back onto boundaries
789 add_boundary_node(2, Node_pt[node_count]);
790
791 // Increment the node number
792 node_count++;
793 }
794
795 // Final node
796 // Allocate new memory for a periodic node
797 Node_pt[node_count] = finite_element_pt(Nx * (Ny - 1) + Nx - 1)
798 ->construct_boundary_node(
799 Np * (Np - 1) + (Np - 1), time_stepper_pt);
800
801 // If required make it periodic
802 if (Xperiodic == true)
803 {
804 Node_pt[node_count]->make_periodic(
805 finite_element_pt(Nx * (Ny - 1))->node_pt(Np * (Np - 1)));
806 }
807
808 // Set the position of the node
809 Node_pt[node_count]->x(0) =
810 x_spacing_function(Nx - 1, Np - 1, Ny - 1, Np - 1);
811 Node_pt[node_count]->x(1) =
812 y_spacing_function(Nx - 1, Np - 1, Ny - 1, Np - 1);
813
814 // Push the node back onto boundaries
815 add_boundary_node(1, Node_pt[node_count]);
816 add_boundary_node(2, Node_pt[node_count]);
817
818 // Increment the node number
819 node_count++;
820 }
821 // END OF FINAL ELEMENT IN MESH
822 }
823
824 // Setup boundary element lookup schemes
825 setup_boundary_element_info();
826 }
827
828 //============================================================================
829 /// Reorder the elements so they are listed in vertical slices
830 /// (more efficient during the frontal solution if the domain
831 /// is long in the x-direction.
832 //============================================================================
833 template<class ELEMENT>
834 void RectangularQuadMesh<ELEMENT>::element_reorder()
835 {
836 // Find out how many elements there are
837 unsigned long Nelement = nelement();
838 // Create a dummy array of elements
839 Vector<FiniteElement*> dummy;
840
841 // Loop over the elements in horizontal order
842 for (unsigned long j = 0; j < Nx; j++)
843 {
844 // Loop over the elements vertically
845 for (unsigned long i = 0; i < Ny; i++)
846 {
847 // Push back onto the new stack
848 dummy.push_back(finite_element_pt(Nx * i + j));
849 }
850 }
851
852 // Now copy the reordered elements into the element_pt
853 for (unsigned long e = 0; e < Nelement; e++)
854 {
855 Element_pt[e] = dummy[e];
856 }
857 }
858
859} // namespace oomph
860#endif
oomph::GeompackQuadMesh
Quadrilateral mesh generator; Uses input from Geompack++. See: http://members.shaw....
Definition geompack_mesh.template.h:42
oomph::RectangularQuadMesh::element_reorder
virtual void element_reorder()
Reorder the elements: By default they are ordered in "horizontal" layers (increasing in x,...
Definition rectangular_quadmesh.template.cc:834
oomph::RectangularQuadMesh::build_mesh
void build_mesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
Definition rectangular_quadmesh.template.cc:43
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition annular_domain.h:35
rectangular_quadmesh.template.h