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