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-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_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 
38 namespace 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
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.
////////////////////////////////////////////////////////////////////// //////////////////////////////...