simple_rectangular_tri_mesh.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_TRIMESH_TEMPLATE_CC
27#define OOMPH_SIMPLE_RECTANGULAR_TRIMESH_TEMPLATE_CC
28
29
30#include <algorithm>
31
32#include "../generic/Telements.h"
33
34// Simple 2D triangle mesh class
36
37
38namespace oomph
39{
40 //====================================================================
41 /// Constructor for simple 2D triangular mesh class:
42 ///
43 /// n_x : number of elements in the x direction
44 ///
45 /// n_y : number of elements in the y direction
46 ///
47 /// l_x : length in the x direction
48 ///
49 /// l_y : length in the y direction
50 ///
51 /// Ordering of elements: 'lower left' to 'lower right' then 'upwards'
52 //====================================================================
53 template<class ELEMENT>
55 const unsigned& n_x,
56 const unsigned& n_y,
57 const double& l_x,
58 const double& l_y,
59 TimeStepper* time_stepper_pt)
60 : Nx(n_x), Ny(n_y), Lx(l_x), Ly(l_y)
61 {
62 using namespace MathematicalConstants;
63
64 // Mesh can only be built with 2D Telements.
65 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
66
67 // Set number of boundaries
68 this->set_nboundary(4);
69
70 // Allocate the store for the elements
71 unsigned n_element = (Nx) * (Ny)*2;
72 Element_pt.resize(n_element, 0);
73
74 // Create first element
75 Element_pt[0] = new ELEMENT;
76
77 // Currently this mesh only works for 3 and 6 noded triangles
78 if ((finite_element_pt(0)->nnode_1d() != 2) &&
79 (finite_element_pt(0)->nnode_1d() != 3) &&
80 (finite_element_pt(0)->nnode_1d() != 4))
81 {
82 throw OomphLibError(
83 "Currently this mesh only works for 3, 6 & 10-noded triangles",
84 OOMPH_CURRENT_FUNCTION,
85 OOMPH_EXCEPTION_LOCATION);
86 }
87
88 unsigned n_node = 0;
89 // Unless nnode_1d returned as !=2 default no extra nodes
90 unsigned additional_nnode = 0;
91
92 // Allocate storage if no extra nodes
93 if (finite_element_pt(0)->nnode_1d() == 2)
94 {
95 n_node = (Nx + 1) * (Ny + 1);
96 }
97
98 if (finite_element_pt(0)->nnode_1d() == 3)
99 {
100 additional_nnode = 3;
101 // Allocate storage for nodes (including extra)
102 n_node = (2 * Nx + 1) * (2 * Ny + 1);
103 }
104
105 if (finite_element_pt(0)->nnode_1d() == 4)
106 {
107 additional_nnode = 7;
108 // Allocate storage for notes (including extra)
109 n_node = (3 * Nx + 1) * (3 * Ny + 1);
110 }
111
112 Node_pt.resize(n_node);
113
114 // Set up geometrical data
115 //------------------------
116 unsigned long node_count = 0;
117 unsigned long element_count = 0;
118 double xinit = 0.0, yinit = 0.0;
119 // Set the values of the increments
120 double xstep = Lx / (Nx);
121 double ystep = Ly / (Ny);
122
123
124 // FIRST NODE (lower left corner)
125 //-------------------------------
126
127 // Set the corner node
128
129 // Allocate memory for the corner node of the first element
130 //(which is not created loop as all subsequent bottom corners already exist)
131 Node_pt[node_count] =
132 finite_element_pt(0)->construct_node(1, time_stepper_pt);
133
134 // Set the pointer from the element
135 finite_element_pt(0)->node_pt(1) = Node_pt[node_count];
136
137 // Don't increment node_count yet as we still need its containing box
138 //(position and boundaries for first node are allocated in the main loop)
139
140 // CREATE THE ELEMENTS (each has 3 local nodes)
141 //--------------------------------------------------
142 // Elements are created two at a time, the first being lower triangular
143 // and the second upper triangular, so that they form a containing box.
144 // Local nodes are numbered anti-clockwise with node_pt(1) being the
145 // right-angle corner of the triangle.
146 // The global node number, node_count, is considered to define
147 // a containing box, with node_count defined as the node number
148 // of the bottom left corner of the box.
149 for (unsigned j = 0; j < Ny + 1; ++j)
150 {
151 for (unsigned i = 0; i < Nx + 1; ++i)
152 {
153 // CURRENT BOX
154 //(nodes on RHS or top edge of domain do not define a box)
155 if (j < Ny && i < Nx)
156 {
157 // Create two elements
158 //------------------------------
159 if (element_count != 0) // 0th element already exists
160 {
161 Element_pt[element_count] = new ELEMENT;
162 }
163
164 Element_pt[element_count + 1] = new ELEMENT;
165
166
167 // Allocate memory for nodes in the current box
168 //--------------------------------------------
169 // If node on LHS of domain, allocate the top left corner node
170 if (i == 0)
171 {
172 Node_pt[node_count + Nx + 1] =
173 finite_element_pt(element_count)
174 ->construct_node(0, time_stepper_pt);
175 }
176
177 // If on bottom row, allocate the bottom right corner node
178 if (j == 0)
179 {
180 Node_pt[node_count + 1] = finite_element_pt(element_count)
181 ->construct_node(2, time_stepper_pt);
182 }
183
184 // Always need to allocate top right corner node
185 Node_pt[node_count + Nx + 2] = finite_element_pt(element_count + 1)
186 ->construct_node(1, time_stepper_pt);
187
188 // Bottom left corner node is already allocated
189
190
191 // Set the pointers from the elements
192 //----------------------------------
193 finite_element_pt(element_count)->node_pt(0) =
194 Node_pt[node_count + Nx + 1];
195 finite_element_pt(element_count)->node_pt(1) = Node_pt[node_count];
196 finite_element_pt(element_count)->node_pt(2) =
197 Node_pt[node_count + 1];
198 finite_element_pt(element_count + 1)->node_pt(0) =
199 Node_pt[node_count + 1];
200 finite_element_pt(element_count + 1)->node_pt(1) =
201 Node_pt[node_count + Nx + 2];
202 finite_element_pt(element_count + 1)->node_pt(2) =
203 Node_pt[node_count + Nx + 1];
204
205 element_count += 2;
206 }
207
208 // CURRENT (GLOBAL) NODE
209 // Set the position
210 Node_pt[node_count]->x(0) = xinit + i * xstep;
211 Node_pt[node_count]->x(1) = yinit + j * ystep;
212
213 // Add node to any relevant boundaries
214 if (j == 0)
215 {
216 this->convert_to_boundary_node(Node_pt[node_count]);
217 add_boundary_node(0, Node_pt[node_count]);
218 }
219 if (j == Ny)
220 {
221 this->convert_to_boundary_node(Node_pt[node_count]);
222 add_boundary_node(2, Node_pt[node_count]);
223 }
224 if (i == 0)
225 {
226 this->convert_to_boundary_node(Node_pt[node_count]);
227 add_boundary_node(3, Node_pt[node_count]);
228 }
229 if (i == Nx)
230 {
231 this->convert_to_boundary_node(Node_pt[node_count]);
232 add_boundary_node(1, Node_pt[node_count]);
233 }
234
235 // Increment counter
236 node_count++;
237 }
238 }
239
240
241 if (additional_nnode == 3)
242 {
243 // Reset element counter to allow looping over existing elements
244 // to add extra nodes.
245 // Note that the node_count is not reset so no entries are overwritten
246 element_count = 0;
247 for (unsigned j = 0; j < Ny + 1; ++j)
248 {
249 // Note: i counter reduced by 1 since i axis runs through middle of
250 // elements on x-axis
251 for (unsigned i = 0; i < Nx; ++i)
252 {
253 // The additional nodes will be added in stages for ease of
254 // application. Note: local numbering follows the anti-clockwise
255 // pattern, node 3 halfway between nodes 0-1, 4 between 1-2 and the
256 // 5th local node between 2-0.
257
258 // Restricted to stop creation of additional nodes outside the mesh
259 if (j < Ny)
260 {
261 // If on the bottom row middle-node at bottom needs allocating
262 if (j == 0)
263 {
264 Node_pt[node_count] = finite_element_pt(element_count)
265 ->construct_node(4, time_stepper_pt);
266 }
267
268 // Due to directions of iteration node at middle of top box edge
269 // (in next element) always needs allocating
270 Node_pt[node_count + Nx] = finite_element_pt(element_count + 1)
271 ->construct_node(4, time_stepper_pt);
272
273 // Set pointers to the top/bottom middle nodes
274 // Bottom node
275 finite_element_pt(element_count)->node_pt(4) = Node_pt[node_count];
276 // Top node
277 finite_element_pt(element_count + 1)->node_pt(4) =
278 Node_pt[node_count + Nx];
279
280 // Increase the element counter to add top/bottom nodes to
281 // next pair of element on next pass
282 element_count += 2;
283 } // End extra top/bottom node construction
284
285 // Set position of current (Global) Node
286 Node_pt[node_count]->x(0) = xinit + double(i + 0.5) * xstep;
287 Node_pt[node_count]->x(1) = yinit + j * ystep;
288
289 // Add node to any applicable boundaries (node 4's can only be top
290 // or bottom boundaries)
291 if (j == 0)
292 {
293 this->convert_to_boundary_node(Node_pt[node_count]);
294 add_boundary_node(0, Node_pt[node_count]);
295 }
296 if (j == Ny)
297 {
298 this->convert_to_boundary_node(Node_pt[node_count]);
299 add_boundary_node(2, Node_pt[node_count]);
300 }
301
302 // Update node_count
303 node_count++;
304 }
305 }
306
307 // Next stage of additional node implementation involes the middle left
308 // and right nodes (local number 3 on each triangle)
309
310 // Element counter reset for second loop over existing elements
311 element_count = 0;
312 // Note: j counter reduced by 1 since j axis runs through middle of
313 // elements on y-axis
314 for (unsigned j = 0; j < Ny; ++j)
315 {
316 for (unsigned i = 0; i < Nx + 1; ++i)
317 {
318 if (j < Ny && i < Nx)
319 {
320 // If on the left hand boundary the node at middle of the left
321 // side needs allocating
322 if (i == 0)
323 {
324 Node_pt[node_count] = finite_element_pt(element_count)
325 ->construct_node(3, time_stepper_pt);
326 }
327
328 // The mid node on the right hand side always needs constructing
329 // within the bounds of the mesh
330 Node_pt[node_count + 1] = finite_element_pt(element_count + 1)
331 ->construct_node(3, time_stepper_pt);
332
333 // Set pointers from the elements to new nodes
334 finite_element_pt(element_count)->node_pt(3) = Node_pt[node_count];
335 finite_element_pt(element_count + 1)->node_pt(3) =
336 Node_pt[node_count + 1];
337
338 // Increase element_count by 2
339 element_count += 2;
340 } // End extra left/right node construction
341
342 // Set position of current (Global) Node
343 Node_pt[node_count]->x(0) = xinit + i * xstep;
344 Node_pt[node_count]->x(1) = yinit + double(j + 0.5) * ystep;
345
346 // Add node to any applicable boundaries again - only be left/right
347 if (i == 0)
348 {
349 this->convert_to_boundary_node(Node_pt[node_count]);
350 add_boundary_node(3, Node_pt[node_count]);
351 }
352 if (i == Nx)
353 {
354 this->convert_to_boundary_node(Node_pt[node_count]);
355 add_boundary_node(1, Node_pt[node_count]);
356 }
357
358 // Update node_count
359 node_count++;
360 }
361 }
362
363 // Final stage of inserting extra nodes - inclusion of the local
364 // number 5 (middle of hypot. edge)
365
366 element_count = 0;
367 // Note: both i,j counters reduced by 1 since j axis runs through middle
368 // of elements in both x,y
369 for (unsigned j = 0; j < Ny; ++j)
370 {
371 for (unsigned i = 0; i < Nx; ++i)
372 {
373 // The middle node always needs constructing
374 Node_pt[node_count] = finite_element_pt(element_count)
375 ->construct_node(5, time_stepper_pt);
376
377 // Set pointers from the elements to new nodes
378 finite_element_pt(element_count)->node_pt(5) = Node_pt[node_count];
379 finite_element_pt(element_count + 1)->node_pt(5) =
380 Node_pt[node_count];
381
382 // Increase element_count by 2
383 element_count += 2;
384 // End extra left/right node construction
385
386 // Set position of current (Global) Node
387 Node_pt[node_count]->x(0) = xinit + double(i + 0.5) * xstep;
388 Node_pt[node_count]->x(1) = yinit + double(j + 0.5) * ystep;
389
390 // All nodes are internal in this structure so no boundaries
391 // applicable
392
393 // Update node_count
394 node_count++;
395 }
396 }
397 }
398 // End of extra nodes for 6 noded trianglur elements
399
400
401 if (additional_nnode == 7)
402 {
403 // Reset element counter to allow looping over existing elements
404 // to add extra nodes.
405 // Note that the node_count is not reset so no entries are overwritten
406 element_count = 0;
407 for (unsigned j = 0; j < Ny + 1; ++j)
408 {
409 // Note: i counter reduced by 1 for implementation of lower element
410 // node 5 and upper element node 6's.
411 for (unsigned i = 0; i < Nx; ++i)
412 {
413 // The additional nodes will be added in stages for ease of
414 // application. Note: local numbering follows the anti-clockwise
415 // pattern, nodes 3 and 4 equidistant between nodes 0-1, 5 and 6
416 // between 1-2, 7 and 8 between 2-0 and last node 9 located
417 // (internally) at the centroid of the triangle.
418
419 // Restricted to stop creation of additional nodes outside the mesh
420 if (j < Ny)
421 {
422 // If on the bottom row middle-left-node at bottom needs allocating
423 if (j == 0)
424 {
425 Node_pt[node_count] = finite_element_pt(element_count)
426 ->construct_node(5, time_stepper_pt);
427 }
428
429 // Due to directions of iteration node at middle left of top box
430 // edge (in next element) always needs allocating
431 Node_pt[node_count + Nx] = finite_element_pt(element_count + 1)
432 ->construct_node(6, time_stepper_pt);
433
434 // Set pointers to the top/bottom middle nodes
435 // Bottom node
436 finite_element_pt(element_count)->node_pt(5) = Node_pt[node_count];
437 // Top node
438 finite_element_pt(element_count + 1)->node_pt(6) =
439 Node_pt[node_count + Nx];
440
441 // Increase the element counter to add top/bottom nodes to
442 // next pair of element on next pass
443 element_count += 2;
444 } // End extra top/bottom node construction
445
446 // Set position of current (Global) Node
447 Node_pt[node_count]->x(0) = xinit + double(i + 1.0 / 3.0) * xstep;
448 Node_pt[node_count]->x(1) = yinit + j * ystep;
449
450 // Add node to any applicable boundaries (exactly as previous)
451 if (j == 0)
452 {
453 this->convert_to_boundary_node(Node_pt[node_count]);
454 add_boundary_node(0, Node_pt[node_count]);
455 }
456 if (j == Ny)
457 {
458 this->convert_to_boundary_node(Node_pt[node_count]);
459 add_boundary_node(2, Node_pt[node_count]);
460 }
461
462 // Update node_count
463 node_count++;
464 }
465 }
466
467
468 // Next stage: bottom-middle-right node (node 6 in lower tri.el.)
469 // and top-middle-right node (node 5 in upper tri.el.)
470
471 element_count = 0;
472 for (unsigned j = 0; j < Ny + 1; ++j)
473 {
474 // Note: i counter as for above 5/6's
475 for (unsigned i = 0; i < Nx; ++i)
476 {
477 // Restricted to stop creation of additional nodes outside the mesh
478 if (j < Ny)
479 {
480 // If on the bottom row middle-right-node at bottom needs allocating
481 if (j == 0)
482 {
483 Node_pt[node_count] = finite_element_pt(element_count)
484 ->construct_node(6, time_stepper_pt);
485 }
486
487 // Due to directions of iteration node at middle left of top box
488 // edge (in next element) always needs allocating
489 Node_pt[node_count + Nx] = finite_element_pt(element_count + 1)
490 ->construct_node(5, time_stepper_pt);
491
492 // Set pointers to the top/bottom middle nodes
493 // Bottom node
494 finite_element_pt(element_count)->node_pt(6) = Node_pt[node_count];
495 // Top node
496 finite_element_pt(element_count + 1)->node_pt(5) =
497 Node_pt[node_count + Nx];
498
499 // Increase the element counter to add top/bottom nodes to
500 // next pair of element on next pass
501 element_count += 2;
502 } // End extra top/bottom node construction
503
504 // Set position of current (Global) Node
505 Node_pt[node_count]->x(0) = xinit + double(i + 2.0 / 3.0) * xstep;
506 Node_pt[node_count]->x(1) = yinit + j * ystep;
507
508 // Add node to any applicable boundaries (exactly as previous)
509 if (j == 0)
510 {
511 this->convert_to_boundary_node(Node_pt[node_count]);
512 add_boundary_node(0, Node_pt[node_count]);
513 }
514 if (j == Ny)
515 {
516 this->convert_to_boundary_node(Node_pt[node_count]);
517 add_boundary_node(2, Node_pt[node_count]);
518 }
519
520 // Update node_count
521 node_count++;
522 }
523 }
524
525
526 // Next stage of additional node implementation involes node 4 on lower
527 // tri. el. and node 3 on upper tri. el.
528
529 // Element counter reset for next loop over existing elements
530 element_count = 0;
531 // Note: j counter reduced by 1 similar to others above
532 for (unsigned j = 0; j < Ny; ++j)
533 {
534 for (unsigned i = 0; i < Nx + 1; ++i)
535 {
536 if (j < Ny && i < Nx)
537 {
538 // If on the left hand boundary the lower middle node of the left
539 // side needs allocating
540 if (i == 0)
541 {
542 Node_pt[node_count] = finite_element_pt(element_count)
543 ->construct_node(4, time_stepper_pt);
544 }
545
546 // The mid node on the right hand side always needs constructing
547 // within the bounds of the mesh
548 Node_pt[node_count + 1] = finite_element_pt(element_count + 1)
549 ->construct_node(3, time_stepper_pt);
550
551 // Set pointers from the elements to new nodes
552 finite_element_pt(element_count)->node_pt(4) = Node_pt[node_count];
553 finite_element_pt(element_count + 1)->node_pt(3) =
554 Node_pt[node_count + 1];
555
556 // Increase element_count by 2
557 element_count += 2;
558 } // End extra left/right node construction
559
560 // Set position of current (Global) Node
561 Node_pt[node_count]->x(0) = xinit + i * xstep;
562 Node_pt[node_count]->x(1) = yinit + double(j + 1.0 / 3.0) * ystep;
563
564 // Add node to any applicable boundaries again - only be left/right
565 if (i == 0)
566 {
567 this->convert_to_boundary_node(Node_pt[node_count]);
568 add_boundary_node(3, Node_pt[node_count]);
569 }
570 if (i == Nx)
571 {
572 this->convert_to_boundary_node(Node_pt[node_count]);
573 add_boundary_node(1, Node_pt[node_count]);
574 }
575
576 // Update node_count
577 node_count++;
578 }
579 }
580
581
582 // Next stage of additional node implementation involes node 3 on lower
583 // tri. el. and node 4 on upper tri. el.
584
585 // Element counter reset for next loop over existing elements
586 element_count = 0;
587 // Note: j counter reduced by 1 similar to others above
588 for (unsigned j = 0; j < Ny; ++j)
589 {
590 for (unsigned i = 0; i < Nx + 1; ++i)
591 {
592 if (j < Ny && i < Nx)
593 {
594 // If on the left hand boundary the lower middle node of the left
595 // side needs allocating
596 if (i == 0)
597 {
598 Node_pt[node_count] = finite_element_pt(element_count)
599 ->construct_node(3, time_stepper_pt);
600 }
601
602 // The mid node on the right hand side always needs constructing
603 // within the bounds of the mesh
604 Node_pt[node_count + 1] = finite_element_pt(element_count + 1)
605 ->construct_node(4, time_stepper_pt);
606
607 // Set pointers from the elements to new nodes
608 finite_element_pt(element_count)->node_pt(3) = Node_pt[node_count];
609 finite_element_pt(element_count + 1)->node_pt(4) =
610 Node_pt[node_count + 1];
611
612 // Increase element_count by 2
613 element_count += 2;
614 } // End extra left/right node construction
615
616 // Set position of current (Global) Node
617 Node_pt[node_count]->x(0) = xinit + i * xstep;
618 Node_pt[node_count]->x(1) = yinit + double(j + 2.0 / 3.0) * ystep;
619
620 // Add node to any applicable boundaries again - only be left/right
621 if (i == 0)
622 {
623 this->convert_to_boundary_node(Node_pt[node_count]);
624 add_boundary_node(3, Node_pt[node_count]);
625 }
626 if (i == Nx)
627 {
628 this->convert_to_boundary_node(Node_pt[node_count]);
629 add_boundary_node(1, Node_pt[node_count]);
630 }
631
632 // Update node_count
633 node_count++;
634 }
635 }
636
637
638 // Next is the lower tri. el. totally internal node with local number 9
639 element_count = 0;
640 // Note: both i,j counters reduced by 1
641 for (unsigned j = 0; j < Ny; ++j)
642 {
643 for (unsigned i = 0; i < Nx; ++i)
644 {
645 // The middle node always needs constructing
646 Node_pt[node_count] = finite_element_pt(element_count)
647 ->construct_node(9, time_stepper_pt);
648
649 // Set pointers from the elements to new nodes
650 finite_element_pt(element_count)->node_pt(9) = Node_pt[node_count];
651
652 // Increase element_count by 2
653 element_count += 2;
654
655 // Set position of current (Global) Node
656 Node_pt[node_count]->x(0) = xinit + double(i + 1.0 / 3.0) * xstep;
657 Node_pt[node_count]->x(1) = yinit + double(j + 1.0 / 3.0) * ystep;
658
659 // All nodes are internal in this structure so no boundaries
660 // applicable
661
662 // Update node_count
663 node_count++;
664 }
665 }
666
667
668 // Next is the bottom hyp. node -
669 // lower tri. el. number 7, upper tri. el. number 8
670 element_count = 0;
671 // Note: both i,j counters reduced by 1
672 for (unsigned j = 0; j < Ny; ++j)
673 {
674 for (unsigned i = 0; i < Nx; ++i)
675 {
676 // The node always needs constructing
677 Node_pt[node_count] = finite_element_pt(element_count)
678 ->construct_node(7, time_stepper_pt);
679
680 // Set pointers from the elements to new nodes
681 finite_element_pt(element_count)->node_pt(7) = Node_pt[node_count];
682 finite_element_pt(element_count + 1)->node_pt(8) =
683 Node_pt[node_count];
684
685 // Increase element_count by 2
686 element_count += 2;
687
688 // Set position of current (Global) Node
689 Node_pt[node_count]->x(0) = xinit + double(i + 2.0 / 3.0) * xstep;
690 Node_pt[node_count]->x(1) = yinit + double(j + 1.0 / 3.0) * ystep;
691
692 // All nodes are internal in this structure so no boundaries
693 // applicable
694
695 // Update node_count
696 node_count++;
697 }
698 }
699
700
701 // Next is the upper hyp. node -
702 // lower tri. el. number 8, upper tri. el. number 7
703 element_count = 0;
704 // Note: both i,j counters reduced by 1
705 for (unsigned j = 0; j < Ny; ++j)
706 {
707 for (unsigned i = 0; i < Nx; ++i)
708 {
709 // The node always needs constructing
710 Node_pt[node_count] = finite_element_pt(element_count)
711 ->construct_node(8, time_stepper_pt);
712
713 // Set pointers from the elements to new nodes
714 finite_element_pt(element_count)->node_pt(8) = Node_pt[node_count];
715 finite_element_pt(element_count + 1)->node_pt(7) =
716 Node_pt[node_count];
717
718 // Increase element_count by 2
719 element_count += 2;
720
721 // Set position of current (Global) Node
722 Node_pt[node_count]->x(0) = xinit + double(i + 1.0 / 3.0) * xstep;
723 Node_pt[node_count]->x(1) = yinit + double(j + 2.0 / 3.0) * ystep;
724
725 // All nodes are internal in this structure so no boundaries
726 // applicable
727
728 // Update node_count
729 node_count++;
730 }
731 }
732
733
734 // Next is the upper tri. el. totally internal node with local number 9
735 element_count = 0;
736 // Note: both i,j counters reduced by 1
737 for (unsigned j = 0; j < Ny; ++j)
738 {
739 for (unsigned i = 0; i < Nx; ++i)
740 {
741 // The middle node always needs constructing
742 Node_pt[node_count] = finite_element_pt(element_count + 1)
743 ->construct_node(9, time_stepper_pt);
744
745 // Set pointers from the elements to new nodes
746 finite_element_pt(element_count + 1)->node_pt(9) =
747 Node_pt[node_count];
748
749 // Increase element_count by 2
750 element_count += 2;
751
752 // Set position of current (Global) Node
753 Node_pt[node_count]->x(0) = xinit + double(i + 2.0 / 3.0) * xstep;
754 Node_pt[node_count]->x(1) = yinit + double(j + 2.0 / 3.0) * ystep;
755
756 // All nodes are internal in this structure so no boundaries
757 // applicable
758
759 // Update node_count
760 node_count++;
761 }
762 }
763 }
764 }
765
766} // namespace oomph
767#endif
cstr elem_len * i
Definition: cfortran.h:603
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
Definition: elements.h:2509
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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:2218
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
void convert_to_boundary_node(Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
A function that upgrades an ordinary node to a boundary node We shouldn't ever really use this,...
Definition: mesh.cc:2590
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
An OomphLibError object which should be thrown when an run-time error is encountered....
SimpleRectangularTriMesh(const unsigned &n_x, const unsigned &n_y, const double &l_x, const double &l_y, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor n_x : number of elements in the x direction; n_y : number of elements in the y direction;...
double Lx
Length of mesh in x-direction.
unsigned Ny
Number of elements in y directions.
unsigned Nx
Number of elements in x direction.
double Ly
Length of mesh in y-direction.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
//////////////////////////////////////////////////////////////////// ////////////////////////////////...