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