simple_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-2025 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_SIMPLE_RECTANGULAR_QUADMESH_TEMPLATE_HEADER
27#define OOMPH_SIMPLE_RECTANGULAR_QUADMESH_TEMPLATE_HEADER
28
29#ifndef OOMPH_SIMPLE_RECTANGULAR_QUADMESH_HEADER
30#error __FILE__ should only be included from simple_rectangular_quadmesh.h.
31#endif // OOMPH_SIMPLE_RECTANGULAR_QUADMESH_HEADER
32
33#include "generic/Qelements.h"
34
35namespace oomph
36{
37 //=======================================================================
38 /// Constructor for 2D Quad mesh class:
39 ///
40 /// Nx : number of elements in the x direction
41 ///
42 /// Ny : number of elements in the y direction
43 ///
44 /// Lx : length in the x direction
45 ///
46 /// Ly : length in the y direction
47 ///
48 /// Ordering of elements: 'Lower left' to 'lower right' then 'upwards'
49 ///
50 /// Timestepper defaults to Steady.
51 //=======================================================================
52 template<class ELEMENT>
54 const unsigned& Nx,
55 const unsigned& Ny,
56 const double& Lx,
57 const double& Ly,
58 TimeStepper* time_stepper_pt)
59 {
60 // Mesh can only be built with 2D Qelements.
61 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
62
63 // Set the internal values
64 NX = Nx;
65 NY = Ny;
66
67 // Set the number of boundaries
68 set_nboundary(4);
69
70 // Allocate the store for the elements
71 Element_pt.resize(Nx * Ny);
72
73 // Create first element
74 Element_pt[0] = new ELEMENT;
75
76 // Read out the number of linear points in the element
77 unsigned n_p = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
78
79 // Can now allocate the store for the nodes
80 Node_pt.resize((1 + (n_p - 1) * Nx) * (1 + (n_p - 1) * Ny));
81
82 // Set up geometrical data
83 //------------------------
84
85 unsigned long node_count = 0;
86 double xinit = 0.0, yinit = 0.0;
87
88 // Set the values of the increments
89 // double xstep = Lx/((n_p-1)*Nx);
90 // double ystep = Ly/((n_p-1)*Ny);
91
92 // Set the length of the element
93 double el_length_x = Lx / double(Nx);
94 double el_length_y = Ly / double(Ny);
95
96 // Storage for local coordinate in element
98
99 // Now assign the topology
100 // Boundaries are numbered 0 1 2 3 from the bottom proceeding anticlockwise
101 // Pinned value are denoted by an integer value 1
102 // Thus if a node is on two boundaries, ORing the values of the
103 // boundary conditions will give the most restrictive case (pinning)
104
105 // FIRST ELEMENT (lower left corner)
106 //----------------------------------
107
108 // Set the corner node
109
110 // Allocate memory for the node
112 finite_element_pt(0)->construct_boundary_node(0, time_stepper_pt);
113
114 // Set the pointer from the element
115 finite_element_pt(0)->node_pt(0) = Node_pt[node_count];
116
117 // Set the position of the node
118 Node_pt[node_count]->x(0) = xinit;
119 Node_pt[node_count]->x(1) = yinit;
120
121 // Add the node to boundaries 0 and 3
122 add_boundary_node(0, Node_pt[node_count]);
123 add_boundary_node(3, Node_pt[node_count]);
124
125 // Increment the node number
126 node_count++;
127
128 // Loop over the other nodes in the first row
129 for (unsigned l2 = 1; l2 < n_p; l2++)
130 {
131 // Allocate memory for the nodes
133 finite_element_pt(0)->construct_boundary_node(l2, time_stepper_pt);
134
135 // Set the pointer from the element
136 finite_element_pt(0)->node_pt(l2) = Node_pt[node_count];
137
138 // Get the local fraction of the node
139 finite_element_pt(0)->local_fraction_of_node(l2, s_fraction);
140
141 // Set the position of the node
143 Node_pt[node_count]->x(1) = yinit;
144
145 // Add the node to the boundary
146 add_boundary_node(0, Node_pt[node_count]);
147
148 // Increment the node number
149 node_count++;
150 }
151
152 // Loop over the other node columns
153 for (unsigned l1 = 1; l1 < n_p; l1++)
154 {
155 // Allocate memory for the nodes
156 Node_pt[node_count] = finite_element_pt(0)->construct_boundary_node(
158
159 // Set the pointer from the element
160 finite_element_pt(0)->node_pt(l1 * n_p) = Node_pt[node_count];
161
162 // Get the fractional position
163 finite_element_pt(0)->local_fraction_of_node(l1 * n_p, s_fraction);
164
165 // Set the position of the node
166 Node_pt[node_count]->x(0) = xinit;
168
169 // Add the node to the boundary
170 add_boundary_node(3, Node_pt[node_count]);
171
172 // Increment the node number
173 node_count++;
174
175 // Loop over the other nodes in the row
176 for (unsigned l2 = 1; l2 < n_p; l2++)
177 {
178 // Allocate the memory for the node
180 finite_element_pt(0)->construct_node(l1 * n_p + l2, time_stepper_pt);
181
182 // Set the pointer from the element
183 finite_element_pt(0)->node_pt(l1 * n_p + l2) = Node_pt[node_count];
184
185 // Get the fractional position
186 finite_element_pt(0)->local_fraction_of_node(l1 * n_p + l2, s_fraction);
187
188 // Set the position of the node
191
192 // Increment the node number
193 node_count++;
194 }
195 }
196
197 // CENTRE OF FIRST ROW OF ELEMENTS
198 //--------------------------------
199 // Now loop over the first row of elements, apart from final element
200 for (unsigned j = 1; j < (Nx - 1); j++)
201 {
202 // Allocate memory for new element
203 Element_pt[j] = new ELEMENT;
204
205 // Do first row of nodes
206
207 // First column of nodes is same as neighbouring element
208 finite_element_pt(j)->node_pt(0) =
209 finite_element_pt(j - 1)->node_pt((n_p - 1));
210
211 // New nodes for other columns
212 for (unsigned l2 = 1; l2 < n_p; l2++)
213 {
214 // Allocate memory for the nodes
216 finite_element_pt(j)->construct_boundary_node(l2, time_stepper_pt);
217
218 // Set the pointer from the element
219 finite_element_pt(j)->node_pt(l2) = Node_pt[node_count];
220
221 // Get the fractional position of the node
222 finite_element_pt(j)->local_fraction_of_node(l2, s_fraction);
223
224 // Set the position of the node
225 Node_pt[node_count]->x(0) = xinit + el_length_x * (j + s_fraction[0]);
226 Node_pt[node_count]->x(1) = yinit;
227
228 // Add the node to the boundary
229 add_boundary_node(0, Node_pt[node_count]);
230
231 // Increment the node number
232 node_count++;
233 }
234
235 // Do the rest of the nodes
236 for (unsigned l1 = 1; l1 < n_p; l1++)
237 {
238 // First column of nodes is same as neighbouring element
239 finite_element_pt(j)->node_pt(l1 * n_p) =
240 finite_element_pt(j - 1)->node_pt(l1 * n_p + (n_p - 1));
241
242 // New nodes for other columns
243 for (unsigned l2 = 1; l2 < n_p; l2++)
244 {
245 // Allocate memory for the nodes
246 Node_pt[node_count] = finite_element_pt(j)->construct_node(
247 l1 * n_p + l2, time_stepper_pt);
248
249 // Set the pointer from the element
250 finite_element_pt(j)->node_pt(l1 * n_p + l2) = Node_pt[node_count];
251
252 // Get the fractional position of the node
253 finite_element_pt(j)->local_fraction_of_node(l1 * n_p + l2,
254 s_fraction);
255
256 // Set the position of the node
257 Node_pt[node_count]->x(0) = xinit + el_length_x * (j + s_fraction[0]);
259
260 // Increment the node number
261 node_count++;
262 }
263 }
264 }
265
266 // FINAL ELEMENT IN FIRST ROW (lower right corner)
267 //-----------------------------------------------
268
269 // Allocate memory for element
270 Element_pt[Nx - 1] = new ELEMENT;
271 // First column of nodes is same as neighbouring element
272 finite_element_pt(Nx - 1)->node_pt(0) =
273 finite_element_pt(Nx - 2)->node_pt(n_p - 1);
274
275 // New middle nodes
276 for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
277 {
278 // Allocate memory for node
280 finite_element_pt(Nx - 1)->construct_boundary_node(l2, time_stepper_pt);
281
282 // Set the pointer from the element
283 finite_element_pt(Nx - 1)->node_pt(l2) = Node_pt[node_count];
284
285 // Get the fractional position of the node
286 finite_element_pt(Nx - 1)->local_fraction_of_node(l2, s_fraction);
287
288 // Set the position of the node
289 Node_pt[node_count]->x(0) =
290 xinit + el_length_x * (Nx - 1 + s_fraction[0]);
291 Node_pt[node_count]->x(1) = yinit;
292
293 // Add the node to the boundary
294 add_boundary_node(0, Node_pt[node_count]);
295
296 // Increment the node number
297 node_count++;
298 }
299
300 // New final node
301
302 // Allocate memory for the node
303 Node_pt[node_count] = finite_element_pt(Nx - 1)->construct_boundary_node(
304 n_p - 1, time_stepper_pt);
305
306 // Set the pointer from the element
307 finite_element_pt(Nx - 1)->node_pt(n_p - 1) = Node_pt[node_count];
308
309 // Get the fractional position of the node
310 finite_element_pt(Nx - 1)->local_fraction_of_node(n_p - 1, s_fraction);
311
312 // Set the position of the node
313 Node_pt[node_count]->x(0) = xinit + el_length_x * (Nx - 1 + s_fraction[0]);
314 Node_pt[node_count]->x(1) = yinit;
315
316 // Add the node to the boundaries
317 add_boundary_node(0, Node_pt[node_count]);
318 add_boundary_node(1, Node_pt[node_count]);
319
320 // Increment the node number
321 node_count++;
322
323 // Do the rest of the nodes
324 for (unsigned l1 = 1; l1 < n_p; l1++)
325 {
326 // First column of nodes is same as neighbouring element
327 finite_element_pt(Nx - 1)->node_pt(l1 * n_p) =
328 finite_element_pt(Nx - 2)->node_pt(l1 * n_p + (n_p - 1));
329
330 // New node for middle column
331 for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
332 {
333 // Allocate memory for node
334 Node_pt[node_count] = finite_element_pt(Nx - 1)->construct_node(
335 l1 * n_p + l2, time_stepper_pt);
336
337 // Set the pointer from the element
338 finite_element_pt(Nx - 1)->node_pt(l1 * n_p + l2) = Node_pt[node_count];
339
340 // Get the fractional position
341 finite_element_pt(Nx - 1)->local_fraction_of_node(l1 * n_p + l2,
342 s_fraction);
343
344 // Set the position of the node
345 Node_pt[node_count]->x(0) =
346 xinit + el_length_x * (Nx - 1 + s_fraction[0]);
348
349 // Increment the node number
350 node_count++;
351 }
352
353 // New node for final column
354
355 // Allocate memory for node
356 Node_pt[node_count] = finite_element_pt(Nx - 1)->construct_boundary_node(
357 l1 * n_p + (n_p - 1), time_stepper_pt);
358
359 // Set the pointer from the element
360 finite_element_pt(Nx - 1)->node_pt(l1 * n_p + (n_p - 1)) =
362
363 // Get the fractional position
364 finite_element_pt(Nx - 1)->local_fraction_of_node(l1 * n_p + (n_p - 1),
365 s_fraction);
366
367 // Set the position of the node
368 Node_pt[node_count]->x(0) =
369 xinit + el_length_x * (Nx - 1 + s_fraction[0]);
371
372 // Add the node to the boundary
373 add_boundary_node(1, Node_pt[node_count]);
374
375 // Increment the node number
376 node_count++;
377 }
378
379 // ALL CENTRAL ELEMENT ROWS
380 //------------------------
381
382 // Loop over remaining element rows
383 for (unsigned i = 1; i < (Ny - 1); i++)
384 {
385 // Set the first element in the row
386
387 // Allocate memory for element
388 Element_pt[Nx * i] = new ELEMENT;
389
390 // The first row of nodes is copied from the element below
391 for (unsigned l2 = 0; l2 < n_p; l2++)
392 {
393 finite_element_pt(Nx * i)->node_pt(l2) =
394 finite_element_pt(Nx * (i - 1))->node_pt((n_p - 1) * n_p + l2);
395 }
396
397 // Other rows are new nodes
398 for (unsigned l1 = 1; l1 < n_p; l1++)
399 {
400 // First column of nodes
401
402 // Allocate memory for node
404 finite_element_pt(Nx * i)->construct_boundary_node(l1 * n_p,
406
407 // Set the pointer from the element
408 finite_element_pt(Nx * i)->node_pt(l1 * n_p) = Node_pt[node_count];
409
410 // Get the fractional position of the node
411 finite_element_pt(Nx * i)->local_fraction_of_node(l1 * n_p, s_fraction);
412
413 // Set the position of the node
414 Node_pt[node_count]->x(0) = xinit;
415 Node_pt[node_count]->x(1) = yinit + el_length_y * (i + s_fraction[1]);
416
417 // Add the node to the boundary
418 add_boundary_node(3, Node_pt[node_count]);
419
420 // Increment the node number
421 node_count++;
422
423 // Now do the other columns
424 for (unsigned l2 = 1; l2 < n_p; l2++)
425 {
426 // Allocate memory for node
427 Node_pt[node_count] = finite_element_pt(Nx * i)->construct_node(
428 l1 * n_p + l2, time_stepper_pt);
429
430 // Set the pointer from the element
431 finite_element_pt(Nx * i)->node_pt(l1 * n_p + l2) =
433
434 // Get the fractional position of the node
435 finite_element_pt(Nx * i)->local_fraction_of_node(l1 * n_p + l2,
436 s_fraction);
437
438 // Set the position of the node
440 Node_pt[node_count]->x(1) = yinit + el_length_y * (i + s_fraction[1]);
441
442 // Increment the node number
443 node_count++;
444 }
445 }
446
447 // Now loop over the rest of the elements in the row, apart from the last
448 for (unsigned j = 1; j < (Nx - 1); j++)
449 {
450 // Allocate memory for new element
451 Element_pt[Nx * i + j] = new ELEMENT;
452
453 // The first row is copied from the lower element
454 for (unsigned l2 = 0; l2 < n_p; l2++)
455 {
456 finite_element_pt(Nx * i + j)->node_pt(l2) =
457 finite_element_pt(Nx * (i - 1) + j)->node_pt((n_p - 1) * n_p + l2);
458 }
459
460 for (unsigned l1 = 1; l1 < n_p; l1++)
461 {
462 // First column is same as neighbouring element
463 finite_element_pt(Nx * i + j)->node_pt(l1 * n_p) =
464 finite_element_pt(Nx * i + (j - 1))->node_pt(l1 * n_p + (n_p - 1));
465
466 // New nodes for other columns
467 for (unsigned l2 = 1; l2 < n_p; l2++)
468 {
469 // Allocate memory for the nodes
471 finite_element_pt(Nx * i + j)
472 ->construct_node(l1 * n_p + l2, time_stepper_pt);
473
474 // Set the pointer
475 finite_element_pt(Nx * i + j)->node_pt(l1 * n_p + l2) =
477
478 // Get the fractional position of the node
479 finite_element_pt(Nx * i + j)
480 ->local_fraction_of_node(l1 * n_p + l2, s_fraction);
481
482 // Set the position of the node
483 Node_pt[node_count]->x(0) =
484 xinit + el_length_x * (j + s_fraction[0]);
485 Node_pt[node_count]->x(1) =
486 yinit + el_length_y * (i + s_fraction[1]);
487
488 // Increment the node number
489 node_count++;
490 }
491 }
492 } // End of loop over elements in row
493
494 // Do final element in row
495
496 // Allocate memory for element
497 Element_pt[Nx * i + Nx - 1] = new ELEMENT;
498
499 // The first row is copied from the lower element
500 for (unsigned l2 = 0; l2 < n_p; l2++)
501 {
502 finite_element_pt(Nx * i + Nx - 1)->node_pt(l2) =
503 finite_element_pt(Nx * (i - 1) + Nx - 1)
504 ->node_pt((n_p - 1) * n_p + l2);
505 }
506
507 for (unsigned l1 = 1; l1 < n_p; l1++)
508 {
509 // First column is same as neighbouring element
510 finite_element_pt(Nx * i + Nx - 1)->node_pt(l1 * n_p) =
511 finite_element_pt(Nx * i + Nx - 2)->node_pt(l1 * n_p + (n_p - 1));
512
513 // Middle nodes
514 for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
515 {
516 // Allocate memory for node
518 finite_element_pt(Nx * i + Nx - 1)
519 ->construct_node(l1 * n_p + l2, time_stepper_pt);
520
521 // Set the pointer
522 finite_element_pt(Nx * i + Nx - 1)->node_pt(l1 * n_p + l2) =
524
525 // Get the fractional position of the node
526 finite_element_pt(Nx * i + Nx - 1)
527 ->local_fraction_of_node(l1 * n_p + l2, s_fraction);
528
529 // Set the position of the node
530 Node_pt[node_count]->x(0) =
531 xinit + el_length_x * (Nx - 1 + s_fraction[0]);
532 Node_pt[node_count]->x(1) = yinit + el_length_y * (i + s_fraction[1]);
533
534 // Increment the node number
535 node_count++;
536 }
537
538 // Final node
539
540 // Allocate memory for node
542 finite_element_pt(Nx * i + Nx - 1)
543 ->construct_boundary_node(l1 * n_p + (n_p - 1), time_stepper_pt);
544
545 // Set the pointer
546 finite_element_pt(Nx * i + Nx - 1)->node_pt(l1 * n_p + (n_p - 1)) =
548
549 // Get the fractional position of the node
550 finite_element_pt(Nx * i + Nx - 1)
551 ->local_fraction_of_node(l1 * n_p + (n_p - 1), s_fraction);
552
553 // Set the position of the node
554 Node_pt[node_count]->x(0) =
555 xinit + el_length_x * (Nx - 1 + s_fraction[0]);
556 Node_pt[node_count]->x(1) = yinit + el_length_y * (i + s_fraction[1]);
557
558 // Add the node to the boundary
559 add_boundary_node(1, Node_pt[node_count]);
560
561 // Increment the node number
562 node_count++;
563
564 } // End of loop over rows of nodes in the element
565 } // End of loop over rows of elements
566
567 // FINAL ELEMENT ROW
568 //=================
569
570 // FIRST ELEMENT IN UPPER ROW (upper left corner)
571 //----------------------------------------------
572
573 // Allocate memory for element
574 Element_pt[Nx * (Ny - 1)] = new ELEMENT;
575
576 // The first row of nodes is copied from the element below
577 for (unsigned l2 = 0; l2 < n_p; l2++)
578 {
579 finite_element_pt(Nx * (Ny - 1))->node_pt(l2) =
580 finite_element_pt(Nx * (Ny - 2))->node_pt((n_p - 1) * n_p + l2);
581 }
582
583 // Second row of nodes
584 // First column of nodes
585 for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
586 {
587 // Allocate memory for node
589 finite_element_pt(Nx * (Ny - 1))
590 ->construct_boundary_node(n_p * l1, time_stepper_pt);
591
592 // Set the pointer from the element
593 finite_element_pt(Nx * (Ny - 1))->node_pt(n_p * l1) = Node_pt[node_count];
594
595 // Get the fractional position of the element
596 finite_element_pt(Nx * (Ny - 1))
597 ->local_fraction_of_node(n_p * l1, s_fraction);
598
599 // Set the position of the node
600 Node_pt[node_count]->x(0) = xinit;
601 Node_pt[node_count]->x(1) =
602 yinit + el_length_y * (Ny - 1 + s_fraction[1]);
603
604 // Add the node to the boundary
605 add_boundary_node(3, Node_pt[node_count]);
606
607 // Increment the node number
608 node_count++;
609
610 // Now do the other columns
611 for (unsigned l2 = 1; l2 < n_p; l2++)
612 {
613 // Allocate memory for node
615 finite_element_pt(Nx * (Ny - 1))
616 ->construct_node(n_p * l1 + l2, time_stepper_pt);
617
618 // Set the pointer from the element
619 finite_element_pt(Nx * (Ny - 1))->node_pt(n_p * l1 + l2) =
621
622 // Get the fractional position of the node
623 finite_element_pt(Nx * (Ny - 1))
624 ->local_fraction_of_node(n_p * l1 + l2, s_fraction);
625
626 // Set the position of the node
628 Node_pt[node_count]->x(1) =
629 yinit + el_length_y * (Ny - 1 + s_fraction[1]);
630
631 // Increment the node number
632 node_count++;
633 }
634 }
635
636 // Final row of nodes
637 // First column of nodes
638 // Top left node
639
640 // Allocate memory for node
642 finite_element_pt(Nx * (Ny - 1))
643 ->construct_boundary_node(n_p * (n_p - 1), time_stepper_pt);
644 // Set the pointer from the element
645 finite_element_pt(Nx * (Ny - 1))->node_pt(n_p * (n_p - 1)) =
647
648 // Get the fractional position of the node
649 finite_element_pt(Nx * (Ny - 1))
650 ->local_fraction_of_node(n_p * (n_p - 1), s_fraction);
651
652 // Set the position of the node
653 Node_pt[node_count]->x(0) = xinit;
654 Node_pt[node_count]->x(1) = yinit + el_length_y * Ny;
655
656 // Add the node to the boundaries
657 add_boundary_node(2, Node_pt[node_count]);
658 add_boundary_node(3, Node_pt[node_count]);
659
660 // Increment the node number
661 node_count++;
662
663 // Now do the other columns
664 for (unsigned l2 = 1; l2 < n_p; l2++)
665 {
666 // Allocate memory for the node
668 finite_element_pt(Nx * (Ny - 1))
669 ->construct_boundary_node(n_p * (n_p - 1) + l2, time_stepper_pt);
670
671 // Set the pointer from the element
672 finite_element_pt(Nx * (Ny - 1))->node_pt(n_p * (n_p - 1) + l2) =
674
675 // Get the fractional position of the node
676 finite_element_pt(Nx * (Ny - 1))
677 ->local_fraction_of_node(n_p * (n_p - 1) + l2, s_fraction);
678
679 // Set the position of the node
681 Node_pt[node_count]->x(1) = yinit + el_length_y * Ny;
682
683 // Add the node to the boundary
684 add_boundary_node(2, Node_pt[node_count]);
685
686 // Increment the node number
687 node_count++;
688 }
689
690 // Now loop over the rest of the elements in the row, apart from the last
691 for (unsigned j = 1; j < (Nx - 1); j++)
692 {
693 // Allocate memory for element
694 Element_pt[Nx * (Ny - 1) + j] = new ELEMENT;
695 // The first row is copied from the lower element
696 for (unsigned l2 = 0; l2 < n_p; l2++)
697 {
698 finite_element_pt(Nx * (Ny - 1) + j)->node_pt(l2) =
699 finite_element_pt(Nx * (Ny - 2) + j)->node_pt((n_p - 1) * n_p + l2);
700 }
701
702 // Second rows
703 for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
704 {
705 // First column is same as neighbouring element
706 finite_element_pt(Nx * (Ny - 1) + j)->node_pt(n_p * l1) =
707 finite_element_pt(Nx * (Ny - 1) + (j - 1))
708 ->node_pt(n_p * l1 + (n_p - 1));
709
710 // New nodes for other columns
711 for (unsigned l2 = 1; l2 < n_p; l2++)
712 {
713 // Allocate memory for the node
715 finite_element_pt(Nx * (Ny - 1) + j)
716 ->construct_node(n_p * l1 + l2, time_stepper_pt);
717 // Set the pointer
718 finite_element_pt(Nx * (Ny - 1) + j)->node_pt(n_p * l1 + l2) =
720
721 // Get the fractional position of the node
722 finite_element_pt(Nx * (Ny - 1) + j)
723 ->local_fraction_of_node(n_p * l1 + l2, s_fraction);
724
725 // Set the position of the node
726 Node_pt[node_count]->x(0) = xinit + el_length_x * (j + s_fraction[0]);
727 Node_pt[node_count]->x(1) =
728 yinit + el_length_y * (Ny - 1 + s_fraction[1]);
729
730 // Increment the node number
731 node_count++;
732 }
733 }
734
735 // Top row
736 // First column is same as neighbouring element
737 finite_element_pt(Nx * (Ny - 1) + j)->node_pt(n_p * (n_p - 1)) =
738 finite_element_pt(Nx * (Ny - 1) + (j - 1))
739 ->node_pt(n_p * (n_p - 1) + (n_p - 1));
740 // New nodes for other columns
741 for (unsigned l2 = 1; l2 < n_p; l2++)
742 {
743 // Allocate memory for node
745 finite_element_pt(Nx * (Ny - 1) + j)
746 ->construct_boundary_node(n_p * (n_p - 1) + l2, time_stepper_pt);
747 // Set the pointer
748 finite_element_pt(Nx * (Ny - 1) + j)->node_pt(n_p * (n_p - 1) + l2) =
750
751 // Get the fractional position of the node
752 finite_element_pt(Nx * (Ny - 1) + j)
753 ->local_fraction_of_node(n_p * (n_p - 1) + l2, s_fraction);
754
755 // Set the position of the node
756 Node_pt[node_count]->x(0) = xinit + el_length_x * (j + s_fraction[0]);
757 Node_pt[node_count]->x(1) = yinit + el_length_y * Ny;
758
759 // Add the node to the boundary
760 add_boundary_node(2, Node_pt[node_count]);
761
762 // Increment the node number
763 node_count++;
764 }
765 } // End of loop over central elements in row
766
767 // FINAL ELEMENT IN ROW (upper right corner)
768 //-----------------------------------------
769
770 // Allocate memory for element
771 Element_pt[Nx * (Ny - 1) + Nx - 1] = new ELEMENT;
772 // The first row is copied from the lower element
773 for (unsigned l2 = 0; l2 < n_p; l2++)
774 {
775 finite_element_pt(Nx * (Ny - 1) + Nx - 1)->node_pt(l2) =
776 finite_element_pt(Nx * (Ny - 2) + Nx - 1)
777 ->node_pt((n_p - 1) * n_p + l2);
778 }
779
780 // Second rows
781 for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
782 {
783 // First column is same as neighbouring element
784 finite_element_pt(Nx * (Ny - 1) + Nx - 1)->node_pt(n_p * l1) =
785 finite_element_pt(Nx * (Ny - 1) + Nx - 2)
786 ->node_pt(n_p * l1 + (n_p - 1));
787
788 // Middle nodes
789 for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
790 {
791 // Allocate memory for node
793 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
794 ->construct_node(n_p * l1 + l2, time_stepper_pt);
795 // Set the pointer
796 finite_element_pt(Nx * (Ny - 1) + Nx - 1)->node_pt(n_p * l1 + l2) =
798
799 // Get the fractional position of the node
800 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
801 ->local_fraction_of_node(n_p * l1 + l2, s_fraction);
802
803 // Set the position of the node
804 Node_pt[node_count]->x(0) =
805 xinit + el_length_x * (Nx - 1 + s_fraction[0]);
806 Node_pt[node_count]->x(1) =
807 yinit + el_length_y * (Ny - 1 + s_fraction[1]);
808
809 // Increment the node number
810 node_count++;
811 }
812
813 // Final node
814
815 // Allocate memory for node
817 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
818 ->construct_boundary_node(n_p * l1 + (n_p - 1), time_stepper_pt);
819 // Set the pointer
820 finite_element_pt(Nx * (Ny - 1) + Nx - 1)->node_pt(n_p * l1 + (n_p - 1)) =
822
823 // Get the fractional position
824 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
825 ->local_fraction_of_node(n_p * l1 + (n_p - 1), s_fraction);
826
827 // Set the position of the node
828 Node_pt[node_count]->x(0) = xinit + el_length_x * Nx;
829 Node_pt[node_count]->x(1) =
830 yinit + el_length_y * (Ny - 1 + s_fraction[1]);
831
832 // Add the node to the boundary
833 add_boundary_node(1, Node_pt[node_count]);
834
835 // Increment the node number
836 node_count++;
837
838 } // End of loop over middle rows
839
840 // Final row
841 // First column is same as neighbouring element
842 finite_element_pt(Nx * (Ny - 1) + Nx - 1)->node_pt(n_p * (n_p - 1)) =
843 finite_element_pt(Nx * (Ny - 1) + Nx - 2)
844 ->node_pt(n_p * (n_p - 1) + (n_p - 1));
845
846 // Middle nodes
847 for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
848 {
849 // Allocate memory for node
851 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
852 ->construct_boundary_node(n_p * (n_p - 1) + l2, time_stepper_pt);
853 // Set the pointer
854 finite_element_pt(Nx * (Ny - 1) + Nx - 1)->node_pt(n_p * (n_p - 1) + l2) =
856
857 // Get the fractional position of the node
858 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
859 ->local_fraction_of_node(n_p * (n_p - 1) + l2, s_fraction);
860
861 // Set the position of the node
862 Node_pt[node_count]->x(0) =
863 xinit + el_length_x * (Nx - 1 + s_fraction[0]);
864 // In fluid 2
865 Node_pt[node_count]->x(1) = yinit + el_length_y * Ny;
866
867 // Add the node to the boundary
868 add_boundary_node(2, Node_pt[node_count]);
869
870 // Increment the node number
871 node_count++;
872 }
873
874 // Final node
875
876 // Allocate memory for node
878 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
879 ->construct_boundary_node(n_p * (n_p - 1) + (n_p - 1), time_stepper_pt);
880 // Set the pointer
881 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
882 ->node_pt(n_p * (n_p - 1) + (n_p - 1)) = Node_pt[node_count];
883
884 // Get the fractional position of the node
885 finite_element_pt(Nx * (Ny - 1) + Nx - 1)
886 ->local_fraction_of_node(n_p * (n_p - 1) + (n_p - 1), s_fraction);
887
888 // Set the position of the node
889 Node_pt[node_count]->x(0) = xinit + el_length_x * Nx;
890 Node_pt[node_count]->x(1) = yinit + el_length_y * Ny;
891
892 // Add the node to the boundaries
893 add_boundary_node(1, Node_pt[node_count]);
894 add_boundary_node(2, Node_pt[node_count]);
895
896 // Increment the node number
897 node_count++;
898
899 // Setup lookup scheme that establishes which elements are located
900 // on the mesh boundaries
901 setup_boundary_element_info();
902 }
903
904} // namespace oomph
905#endif
cstr elem_len * i
Definition cfortran.h:603
Node ** Node_pt
Storage for pointers to the nodes in the element.
Definition elements.h:1323
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition elements.h:2222
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition nodes.h:1060
void resize(const unsigned &n_value)
Resize the number of equations.
Definition nodes.cc:2167
SimpleRectangularQuadMesh(const unsigned &Nx, const unsigned &Ny, const double &Lx, const double &Ly, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in the horizontal and vertical directions, and the corresponding...
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).