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