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$
7 //LIC//
8 //LIC// $LastChangedDate$
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_RECTANGULAR_QUADMESH_TEMPLATE_CC
31 #define OOMPH_RECTANGULAR_QUADMESH_TEMPLATE_CC
32 
33 
34 // OOMPH-LIB Headers
35 #include "rectangular_quadmesh.template.h"
36 
37 
38 namespace oomph
39 {
40 
41 //===========================================================================
42 /// Generic mesh construction. This function contains the "guts" of the
43 /// mesh generation process, including all the tedious loops, counting
44 /// and spacing functions. The function should be called in all constuctors
45 /// of any derived classes.
46 //===========================================================================
47 template <class ELEMENT>
48 void RectangularQuadMesh<ELEMENT>::build_mesh(TimeStepper* time_stepper_pt)
49 {
50 
51  // Mesh can only be built with 2D Qelements.
52  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
53 
54  //Set the number of boundaries
55  set_nboundary(4);
56 
57  //Allocate the store for the elements
58  Element_pt.resize(Nx*Ny);
59  //Allocate the memory for the first element
60  Element_pt[0] = new ELEMENT;
61 
62  //Read out the number of linear points in the element
63  Np = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
64 
65  //Can now allocate the store for the nodes
66  Node_pt.resize((1 + (Np-1)*Nx)*(1 + (Np-1)*Ny));
67 
68  //Set up geometrical data
69  unsigned long node_count=0;
70 
71  //Now assign the topology
72  //Boundaries are numbered 0 1 2 3 from the bottom proceeding anticlockwise
73  //Pinned value are denoted by an integer value 1
74  //Thus if a node is on two boundaries, ORing the values of the
75  //boundary conditions will give the most restrictive case (pinning)
76 
77  //FIRST ELEMENT (lower left corner)
78 
79  //Set the corner node
80  //Allocate memory for the node
81  Node_pt[node_count] =
82  finite_element_pt(0)->construct_boundary_node(0,time_stepper_pt);
83 
84  //Set the position of the node
85  Node_pt[node_count]->x(0) = x_spacing_function(0,0,0,0);
86  Node_pt[node_count]->x(1) = y_spacing_function(0,0,0,0);
87 
88  //Push the node back onto boundaries
89  add_boundary_node(0,Node_pt[node_count]);
90  add_boundary_node(3,Node_pt[node_count]);
91 
92  //Increment the node number
93  node_count++;
94 
95  //Loop over the other nodes in the first row
96  for(unsigned l2=1;l2<Np;l2++)
97  {
98  //Allocate memory for the nodes
99  Node_pt[node_count] =
100  finite_element_pt(0)->construct_boundary_node(l2,time_stepper_pt);
101 
102  //Set the position of the node
103  Node_pt[node_count]->x(0) = x_spacing_function(0,l2,0,0);
104  Node_pt[node_count]->x(1) = y_spacing_function(0,l2,0,0);
105 
106  //Push the node back onto boundaries
107  add_boundary_node(0,Node_pt[node_count]);
108 
109  //If we only have one column then the RHS node is on the right boundary
110  if((Nx==1) && (l2==(Np-1)))
111  {add_boundary_node(1,Node_pt[node_count]);}
112 
113  //Increment the node number
114  node_count++;
115  }
116 
117  //Loop over the other node columns
118  for(unsigned l1=1;l1<Np;l1++)
119  {
120  //Allocate memory for the nodes
121  Node_pt[node_count] =
122  finite_element_pt(0)->construct_boundary_node(l1*Np,time_stepper_pt);
123 
124  //Set the position of the node
125  Node_pt[node_count]->x(0) = x_spacing_function(0,0,0,l1);
126  Node_pt[node_count]->x(1) = y_spacing_function(0,0,0,l1);
127 
128  //Push the node back onto boundaries
129  add_boundary_node(3,Node_pt[node_count]);
130 
131  //If we only have one row, then the top node is on the top boundary
132  if((Ny==1) && (l1==(Np-1)))
133  {add_boundary_node(2,Node_pt[node_count]);}
134 
135  //Increment the node number
136  node_count++;
137 
138  //Loop over the other nodes in the row
139  for(unsigned l2=1;l2<Np;l2++)
140  {
141  //Allocate the memory for the node
142  //If it lies on a boundary make a boundary node
143  if(((Nx==1) && (l2==(Np-1))) || ((Ny==1) && (l1==(Np-1))))
144  {
145  Node_pt[node_count] = finite_element_pt(0)->
146  construct_boundary_node(l1*Np+l2,time_stepper_pt);
147  }
148  //Otherwise its a normal node
149  else
150  {
151  Node_pt[node_count] =
152  finite_element_pt(0)->construct_node(l1*Np+l2,time_stepper_pt);
153  }
154 
155  //Set the position of the node
156  Node_pt[node_count]->x(0) = x_spacing_function(0,l2,0,l1);
157  Node_pt[node_count]->x(1) = y_spacing_function(0,l2,0,l1);
158 
159  //If we only have one column then the RHS node is on the right boundary
160  if((Nx==1) && l2==(Np-1))
161  {add_boundary_node(1,Node_pt[node_count]);}
162  //If we only have one row, then the top node is on the top boundary
163  if((Ny==1) && (l1==(Np-1)))
164  {add_boundary_node(2,Node_pt[node_count]);}
165 
166  //Increment the node number
167  node_count++;
168  }
169  }
170  //END OF FIRST ELEMENT
171 
172  //CENTRE OF FIRST ROW OF ELEMENTS
173  //Now loop over the first row of elements, apart from final element
174  for(unsigned j=1;j<(Nx-1);j++)
175  {
176  //Allocate memory for new element
177  Element_pt[j] = new ELEMENT;
178  //Do first row of nodes
179  //First column of nodes is same as neighbouring element
180  finite_element_pt(j)->node_pt(0) =
181  finite_element_pt(j-1)->node_pt((Np-1));
182  //New nodes for other columns
183  for(unsigned l2=1;l2<Np;l2++)
184  {
185  //Allocate memory for the nodes
186  Node_pt[node_count] =
187  finite_element_pt(j)->construct_boundary_node(l2,time_stepper_pt);
188 
189  //Set the position of the node
190  Node_pt[node_count]->x(0) = x_spacing_function(j,l2,0,0);
191  Node_pt[node_count]->x(1) = y_spacing_function(j,l2,0,0);
192 
193  //Push the node back onto boundaries
194  add_boundary_node(0,Node_pt[node_count]);
195 
196  //Increment the node number
197  node_count++;
198  }
199 
200  //Do the rest of the nodes
201  for(unsigned l1=1;l1<Np;l1++)
202  {
203  //First column of nodes is same as neighbouring element
204  finite_element_pt(j)->node_pt(l1*Np) =
205  finite_element_pt(j-1)->node_pt(l1*Np+(Np-1));
206  //New nodes for other columns
207  for(unsigned l2=1;l2<Np;l2++)
208  {
209  //Allocate memory for the nodes
210  //If we only have one row, this node could be on the boundary
211  if((Ny==1) && (l1==(Np-1)))
212  {
213  Node_pt[node_count] = finite_element_pt(j)->
214  construct_boundary_node(l1*Np+l2,time_stepper_pt);
215  }
216  //Otherwise create a normal node
217  else
218  {
219  Node_pt[node_count] =
220  finite_element_pt(j)->construct_node(l1*Np+l2,time_stepper_pt);
221  }
222 
223  //Set the position of the node
224  Node_pt[node_count]->x(0) = x_spacing_function(j,l2,0,l1);
225  Node_pt[node_count]->x(1) = y_spacing_function(j,l2,0,l1);
226 
227  //If we only have one row, then the top node is on the top boundary
228  if((Ny==1) && (l1==(Np-1)))
229  {add_boundary_node(2,Node_pt[node_count]);}
230 
231  //Increment the node number
232  node_count++;
233  }
234  }
235  }
236  //END OF CENTRE OF FIRST ROW OF ELEMENTS
237 
238  //FINAL ELEMENT IN FIRST ROW (lower right corner)
239  //Only allocate if there is more than one element in the row
240  if(Nx > 1)
241  {
242  //Allocate memory for element
243  Element_pt[Nx-1] = new ELEMENT;
244 
245  //First column of nodes is same as neighbouring element
246  finite_element_pt(Nx-1)->node_pt(0) =
247  finite_element_pt(Nx-2)->node_pt(Np-1);
248 
249  //New middle nodes
250  for(unsigned l2=1;l2<(Np-1);l2++)
251  {
252  //Allocate memory for node
253  Node_pt[node_count] =
254  finite_element_pt(Nx-1)->construct_boundary_node(l2,time_stepper_pt);
255 
256  //Set the position of the node
257  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,l2,0,0);
258  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,l2,0,0);
259 
260  //Push the node back onto boundaries
261  add_boundary_node(0,Node_pt[node_count]);
262 
263  //Increment the node number
264  node_count++;
265  }
266 
267  //Allocate memory for a new boundary node
268  Node_pt[node_count] =
269  finite_element_pt(Nx-1)->construct_boundary_node(Np-1,time_stepper_pt);
270 
271  //If required make it periodic from the node on the other side
272  if(Xperiodic==true)
273  {
274  Node_pt[node_count]->
275  make_periodic(finite_element_pt(0)->node_pt(0));
276  }
277 
278  //Set the position of the node
279  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,Np-1,0,0);
280  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,Np-1,0,0);
281 
282  //Push the node back onto boundaries
283  add_boundary_node(0,Node_pt[node_count]);
284  add_boundary_node(1,Node_pt[node_count]);
285 
286  //Increment the node number
287  node_count++;
288 
289  //Do the rest of the nodes
290  for(unsigned l1=1;l1<Np;l1++)
291  {
292  //First column of nodes is same as neighbouring element
293  finite_element_pt(Nx-1)->node_pt(l1*Np) =
294  finite_element_pt(Nx-2)->node_pt(l1*Np+(Np-1));
295 
296  //New node for middle column
297  for(unsigned l2=1;l2<(Np-1);l2++)
298  {
299  //Allocate memory for node
300  //If it lies on the boundary, create a boundary node
301  if((Ny==1) && (l1==(Np-1)))
302  {
303  Node_pt[node_count] = finite_element_pt(Nx-1)->
304  construct_boundary_node(l1*Np+l2,time_stepper_pt);
305  }
306  //Otherwise create a normal node
307  else
308  {
309  Node_pt[node_count] =
310  finite_element_pt(Nx-1)->construct_node(l1*Np+l2,time_stepper_pt);
311  }
312 
313  //Set the position of the node
314  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,l2,0,l1);
315  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,l2,0,l1);
316 
317  //If we only have one row, then the top node is on the top boundary
318  if((Ny==1) && (l1==(Np-1)))
319  {add_boundary_node(2,Node_pt[node_count]);}
320 
321  //Increment the node number
322  node_count++;
323  }
324 
325  //Allocate memory for a new boundary node
326  Node_pt[node_count] = finite_element_pt(Nx-1)->
327  construct_boundary_node(l1*Np+(Np-1),time_stepper_pt);
328 
329  //If required make it periodic from the corresponding node on the other
330  //side
331  if(Xperiodic==true)
332  {
333  Node_pt[node_count]->
334  make_periodic(finite_element_pt(0)->node_pt(l1*Np));
335  }
336 
337  //Set the position of the node
338  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,Np-1,0,l1);
339  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,Np-1,0,l1);
340 
341  //Push the node back onto boundaries
342  add_boundary_node(1,Node_pt[node_count]);
343 
344  //If we only have one row, then the top node is on the top boundary
345  if((Ny==1) && (l1==(Np-1)))
346  {add_boundary_node(2,Node_pt[node_count]);}
347 
348  //Increment the node number
349  node_count++;
350  }
351  }
352  //END OF FIRST ROW OF ELEMENTS
353 
354  //ALL CENTRAL ELEMENT ROWS
355  //Loop over remaining element rows
356  for(unsigned i=1;i<(Ny-1);i++)
357  {
358  //Set the first element in the row
359  //Allocate memory for element
360  Element_pt[Nx*i] = new ELEMENT;
361 
362  //The first row of nodes is copied from the element below
363  for(unsigned l2=0;l2<Np;l2++)
364  {
365  finite_element_pt(Nx*i)->node_pt(l2) =
366  finite_element_pt(Nx*(i-1))->node_pt((Np-1)*Np + l2);
367  }
368 
369  //Other rows are new nodes
370  for(unsigned l1=1;l1<Np;l1++)
371  {
372  //First column of nodes
373  //Allocate memory for node
374  Node_pt[node_count] = finite_element_pt(Nx*i)->
375  construct_boundary_node(l1*Np,time_stepper_pt);
376 
377  //Set the position of the node
378  Node_pt[node_count]->x(0) = x_spacing_function(0,0,i,l1);
379  Node_pt[node_count]->x(1) = y_spacing_function(0,0,i,l1);
380 
381  //Push the node back onto boundaries
382  add_boundary_node(3,Node_pt[node_count]);
383 
384  //Increment the node number
385  node_count++;
386 
387  //Now do the other columns
388  for(unsigned l2=1;l2<Np;l2++)
389  {
390  //Allocate memory for node
391  //If we only have one column, the node could be on the boundary
392  if((Nx==1) && (l2==(Np-1)))
393  {
394  Node_pt[node_count] = finite_element_pt(Nx*i)->
395  construct_boundary_node(l1*Np+l2,time_stepper_pt);
396  }
397  else
398  {
399  Node_pt[node_count] =
400  finite_element_pt(Nx*i)->construct_node(l1*Np+l2,time_stepper_pt);
401  }
402 
403  //Set the position of the node
404  Node_pt[node_count]->x(0) = x_spacing_function(0,l2,i,l1);
405  Node_pt[node_count]->x(1) = y_spacing_function(0,l2,i,l1);
406 
407  //If we only have one column then the RHS node is on the
408  //right boundary
409  if((Nx==1) && (l2==(Np-1)))
410  {add_boundary_node(1,Node_pt[node_count]);}
411 
412  //Increment the node number
413  node_count++;
414  }
415  }
416 
417  //Now loop over the rest of the elements in the row, apart from the last
418  for(unsigned j=1;j<(Nx-1);j++)
419  {
420  //Allocate memory for new element
421  Element_pt[Nx*i+j] = new ELEMENT;
422  //The first row is copied from the lower element
423  for(unsigned l2=0;l2<Np;l2++)
424  {
425  finite_element_pt(Nx*i+j)->node_pt(l2) =
426  finite_element_pt(Nx*(i-1)+j)->node_pt((Np-1)*Np + l2);
427  }
428 
429  for(unsigned l1=1;l1<Np;l1++)
430  {
431  //First column is same as neighbouring element
432  finite_element_pt(Nx*i+j)->node_pt(l1*Np) =
433  finite_element_pt(Nx*i+(j-1))->node_pt(l1*Np+(Np-1));
434  //New nodes for other columns
435  for(unsigned l2=1;l2<Np;l2++)
436  {
437  //Allocate memory for the nodes
438  Node_pt[node_count] =
439  finite_element_pt(Nx*i+j)->
440  construct_node(l1*Np+l2,time_stepper_pt);
441 
442  //Set the position of the node
443  Node_pt[node_count]->x(0) = x_spacing_function(j,l2,i,l1);
444  Node_pt[node_count]->x(1) = y_spacing_function(j,l2,i,l1);
445 
446  //Increment the node number
447  node_count++;
448  }
449  }
450  } //End of loop over elements in row
451 
452  //Do final element in row
453  //Only if there is more than one column
454  if(Nx > 1)
455  {
456  //Allocate memory for element
457  Element_pt[Nx*i+Nx-1] = new ELEMENT;
458  //The first row is copied from the lower element
459  for(unsigned l2=0;l2<Np;l2++)
460  {
461  finite_element_pt(Nx*i+Nx-1)->node_pt(l2) =
462  finite_element_pt(Nx*(i-1)+Nx-1)->node_pt((Np-1)*Np + l2);
463  }
464 
465  for(unsigned l1=1;l1<Np;l1++)
466  {
467  //First column is same as neighbouring element
468  finite_element_pt(Nx*i+Nx-1)->node_pt(l1*Np) =
469  finite_element_pt(Nx*i+Nx-2)->node_pt(l1*Np+(Np-1));
470 
471  //Middle nodes
472  for(unsigned l2=1;l2<(Np-1);l2++)
473  {
474  //Allocate memory for node
475  Node_pt[node_count] =
476  finite_element_pt(Nx*i+Nx-1)->
477  construct_node(l1*Np+l2,time_stepper_pt);
478 
479  //Set the position of the node
480  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,l2,i,l1);
481  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,l2,i,l1);
482 
483  //Increment the node number
484  node_count++;
485  }
486 
487  //Allocate memory for a new boundary node
488  Node_pt[node_count] =
489  finite_element_pt(Nx*i+Nx-1)->
490  construct_boundary_node(l1*Np+(Np-1),time_stepper_pt);
491 
492  //If required make it periodic from the corresponding node on
493  //the other side
494  if(Xperiodic==true)
495  {
496  Node_pt[node_count]->
497  make_periodic(finite_element_pt(Nx*i)->node_pt(l1*Np));
498  }
499 
500  //Set the position of the node
501  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,Np-1,i,l1);
502  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,Np-1,i,l1);
503 
504  //Push the node back onto boundaries
505  add_boundary_node(1,Node_pt[node_count]);
506 
507  //Increment the node number
508  node_count++;
509  } //End of loop over rows of nodes in the element
510  } //End of if more than one column
511  } //End of loop over rows of elements
512 
513  //END OF LOOP OVER CENTRAL ELEMENT ROWS
514 
515  //FINAL ELEMENT ROW
516  //ONLY NECESSARY IF THERE IS MORE THAN ONE ROW
517  if(Ny > 1)
518  {
519  //FIRST ELEMENT IN UPPER ROW (upper left corner)
520  //Allocate memory for element
521  Element_pt[Nx*(Ny-1)] = new ELEMENT;
522  //The first row of nodes is copied from the element below
523  for(unsigned l2=0;l2<Np;l2++)
524  {
525  finite_element_pt(Nx*(Ny-1))->node_pt(l2)
526  = finite_element_pt(Nx*(Ny-2))->node_pt((Np-1)*Np + l2);
527  }
528 
529  //Second row of nodes
530  //First column of nodes
531  for(unsigned l1=1;l1<(Np-1);l1++)
532  {
533  //Allocate memory for node
534  Node_pt[node_count] = finite_element_pt(Nx*(Ny-1))->
535  construct_boundary_node(Np*l1,time_stepper_pt);
536 
537  //Set the position of the node
538  Node_pt[node_count]->x(0) = x_spacing_function(0,0,Ny-1,l1);
539  Node_pt[node_count]->x(1) = y_spacing_function(0,0,Ny-1,l1);
540 
541  //Push the node back onto boundaries
542  add_boundary_node(3,Node_pt[node_count]);
543 
544  //Increment the node number
545  node_count++;
546 
547  //Now do the other columns
548  for(unsigned l2=1;l2<Np;l2++)
549  {
550 
551  //Allocate memory for node
552  if ((Nx==1)&&(l2==Np-1))
553  {
554  Node_pt[node_count] =
555  finite_element_pt(Nx*(Ny-1))->
556  construct_boundary_node(Np*l1+l2,time_stepper_pt);
557  }
558  else
559  {
560  Node_pt[node_count] =
561  finite_element_pt(Nx*(Ny-1))->
562  construct_node(Np*l1+l2,time_stepper_pt);
563  }
564 
565  //Set the position of the node
566  Node_pt[node_count]->x(0) = x_spacing_function(0,l2,Ny-1,l1);
567  Node_pt[node_count]->x(1) = y_spacing_function(0,l2,Ny-1,l1);
568 
569  //Push the node back onto boundaries
570  if ((Nx==1)&&(l2==Np-1))
571  {
572  add_boundary_node(1,Node_pt[node_count]);
573  }
574 
575  //Increment the node number
576  node_count++;
577  }
578  }
579 
580  //Final row of nodes
581  //First column of nodes
582  //Top left node
583  //Allocate memory for node
584  Node_pt[node_count] = finite_element_pt(Nx*(Ny-1))->
585  construct_boundary_node(Np*(Np-1),time_stepper_pt);
586 
587  //Set the position of the node
588  Node_pt[node_count]->x(0) = x_spacing_function(0,0,Ny-1,Np-1);
589  Node_pt[node_count]->x(1) = y_spacing_function(0,0,Ny-1,Np-1);
590 
591  //Push the node back onto boundaries
592  add_boundary_node(2,Node_pt[node_count]);
593  add_boundary_node(3,Node_pt[node_count]);
594 
595  //Increment the node number
596  node_count++;
597 
598  //Now do the other columns
599  for(unsigned l2=1;l2<Np;l2++)
600  {
601  //Allocate memory for the node
602  Node_pt[node_count] =
603  finite_element_pt(Nx*(Ny-1))->
604  construct_boundary_node(Np*(Np-1)+l2,time_stepper_pt);
605 
606  //Set the position of the node
607  Node_pt[node_count]->x(0) = x_spacing_function(0,l2,Ny-1,Np-1);
608  Node_pt[node_count]->x(1) = y_spacing_function(0,l2,Ny-1,Np-1);
609 
610  //Push the node back onto boundaries
611  add_boundary_node(2,Node_pt[node_count]);
612 
613  //Push the node back onto boundaries
614  if ((Nx==1)&&(l2==Np-1))
615  {
616  add_boundary_node(1,Node_pt[node_count]);
617  }
618 
619  //Increment the node number
620  node_count++;
621  }
622 
623  //Now loop over the rest of the elements in the row, apart from the last
624  for(unsigned j=1;j<(Nx-1);j++)
625  {
626  //Allocate memory for element
627  Element_pt[Nx*(Ny-1)+j] = new ELEMENT;
628  //The first row is copied from the lower element
629  for(unsigned l2=0;l2<Np;l2++)
630  {
631  finite_element_pt(Nx*(Ny-1)+j)->node_pt(l2) =
632  finite_element_pt(Nx*(Ny-2)+j)->node_pt((Np-1)*Np + l2);
633  }
634 
635  //Second rows
636  for(unsigned l1=1;l1<(Np-1);l1++)
637  {
638  //First column is same as neighbouring element
639  finite_element_pt(Nx*(Ny-1)+j)->node_pt(Np*l1)
640  = finite_element_pt(Nx*(Ny-1)+(j-1))->node_pt(Np*l1+(Np-1));
641 
642  //New nodes for other columns
643  for(unsigned l2=1;l2<Np;l2++)
644  {
645  //Allocate memory for the node
646  Node_pt[node_count] =
647  finite_element_pt(Nx*(Ny-1)+j)->
648  construct_node(Np*l1+l2,time_stepper_pt);
649 
650  //Set the position of the node
651  Node_pt[node_count]->x(0) = x_spacing_function(j,l2,Ny-1,l1);
652  Node_pt[node_count]->x(1) = y_spacing_function(j,l2,Ny-1,l1);
653 
654  //Increment the node number
655  node_count++;
656  }
657  }
658 
659  //Top row
660  //First column is same as neighbouring element
661  finite_element_pt(Nx*(Ny-1)+j)->node_pt(Np*(Np-1))
662  = finite_element_pt(Nx*(Ny-1)+(j-1))->node_pt(Np*(Np-1)+(Np-1));
663  //New nodes for other columns
664  for(unsigned l2=1;l2<Np;l2++)
665  {
666  //Allocate memory for node
667  Node_pt[node_count] =
668  finite_element_pt(Nx*(Ny-1)+j)->
669  construct_boundary_node(Np*(Np-1)+l2,time_stepper_pt);
670 
671  //Set the position of the node
672  Node_pt[node_count]->x(0) = x_spacing_function(j,l2,Ny-1,Np-1);
673  Node_pt[node_count]->x(1) = y_spacing_function(j,l2,Ny-1,Np-1);
674 
675  //Push the node back onto boundaries
676  add_boundary_node(2,Node_pt[node_count]);
677 
678  //Increment the node number
679  node_count++;
680  }
681  } //End of loop over central elements in row
682 
683  //FINAL ELEMENT IN ROW (upper right corner)
684  //Only if there is more than one column
685  if(Nx > 1)
686  {
687  //Allocate memory for element
688  Element_pt[Nx*(Ny-1)+Nx-1] = new ELEMENT;
689  //The first row is copied from the lower element
690  for(unsigned l2=0;l2<Np;l2++)
691  {
692  finite_element_pt(Nx*(Ny-1)+Nx-1)->node_pt(l2) =
693  finite_element_pt(Nx*(Ny-2)+Nx-1)->node_pt((Np-1)*Np + l2);
694  }
695 
696  //Second rows
697  for(unsigned l1=1;l1<(Np-1);l1++)
698  {
699  //First column is same as neighbouring element
700  finite_element_pt(Nx*(Ny-1)+Nx-1)->node_pt(Np*l1)
701  = finite_element_pt(Nx*(Ny-1)+Nx-2)->node_pt(Np*l1+(Np-1));
702 
703  //Middle nodes
704  for(unsigned l2=1;l2<(Np-1);l2++)
705  {
706  //Allocate memory for node
707  Node_pt[node_count] =
708  finite_element_pt(Nx*(Ny-1)+Nx-1)->
709  construct_node(Np*l1+l2,time_stepper_pt);
710 
711  //Set the position of the node
712  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,l2,Ny-1,l1);
713  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,l2,Ny-1,l1);
714 
715  //Increment the node number
716  node_count++;
717  }
718 
719  //Final node
720  //Allocate new memory for a boundary node
721  Node_pt[node_count] =
722  finite_element_pt(Nx*(Ny-1)+Nx-1)->
723  construct_boundary_node(Np*l1+(Np-1),time_stepper_pt);
724 
725  //If required make it periodic
726  if(Xperiodic==true)
727  {
728  Node_pt[node_count]->
729  make_periodic(finite_element_pt(Nx*(Ny-1))->node_pt(Np*l1));
730  }
731 
732  //Set the position of the node
733  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,Np-1,Ny-1,l1);
734  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,Np-1,Ny-1,l1);
735 
736  //Push the node back onto boundaries
737  add_boundary_node(1,Node_pt[node_count]);
738 
739  //Increment the node number
740  node_count++;
741 
742  } //End of loop over middle rows
743 
744  //Final row
745  //First column is same as neighbouring element
746  finite_element_pt(Nx*(Ny-1)+Nx-1)->node_pt(Np*(Np-1))
747  = finite_element_pt(Nx*(Ny-1)+Nx-2)->node_pt(Np*(Np-1)+(Np-1));
748 
749  //Middle nodes
750  for(unsigned l2=1;l2<(Np-1);l2++)
751  {
752  //Allocate memory for node
753  Node_pt[node_count] =
754  finite_element_pt(Nx*(Ny-1)+Nx-1)->
755  construct_boundary_node(Np*(Np-1)+l2,time_stepper_pt);
756 
757  //Set the position of the node
758  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,l2,Ny-1,Np-1);
759  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,l2,Ny-1,Np-1);
760 
761  //Push the node back onto boundaries
762  add_boundary_node(2,Node_pt[node_count]);
763 
764  //Increment the node number
765  node_count++;
766  }
767 
768  //Final node
769  //Allocate new memory for a periodic node
770  Node_pt[node_count] = finite_element_pt(Nx*(Ny-1)+Nx-1)->
771  construct_boundary_node(Np*(Np-1)+(Np-1),time_stepper_pt);
772 
773  //If required make it periodic
774  if(Xperiodic==true)
775  {
776  Node_pt[node_count]->
777  make_periodic(finite_element_pt(Nx*(Ny-1))->node_pt(Np*(Np-1)));
778  }
779 
780  //Set the position of the node
781  Node_pt[node_count]->x(0) = x_spacing_function(Nx-1,Np-1,Ny-1,Np-1);
782  Node_pt[node_count]->x(1) = y_spacing_function(Nx-1,Np-1,Ny-1,Np-1);
783 
784  //Push the node back onto boundaries
785  add_boundary_node(1,Node_pt[node_count]);
786  add_boundary_node(2,Node_pt[node_count]);
787 
788  //Increment the node number
789  node_count++;
790  }
791  //END OF FINAL ELEMENT IN MESH
792  }
793 
794  // Setup boundary element lookup schemes
795  setup_boundary_element_info();
796 
797 }
798 
799 //============================================================================
800 /// Reorder the elements so they are listed in vertical slices
801 /// (more efficient during the frontal solution if the domain
802 /// is long in the x-direction.
803 //============================================================================
804 template<class ELEMENT>
805 void RectangularQuadMesh<ELEMENT>::element_reorder()
806  {
807  //Find out how many elements there are
808  unsigned long Nelement = nelement();
809  //Create a dummy array of elements
810  Vector<FiniteElement *> dummy;
811 
812  //Loop over the elements in horizontal order
813  for(unsigned long j=0;j<Nx;j++)
814  {
815  //Loop over the elements vertically
816  for(unsigned long i=0;i<Ny;i++)
817  {
818  //Push back onto the new stack
819  dummy.push_back(finite_element_pt(Nx*i + j));
820  }
821  }
822 
823  //Now copy the reordered elements into the element_pt
824  for(unsigned long e=0;e<Nelement;e++) {Element_pt[e] = dummy[e];}
825  }
826 
827 }
828 #endif
oomph::RectangularQuadMesh::element_reorder
virtual void element_reorder()
Reorder the elements: By default they are ordered in "horizontal" layers (increasing in x...
Definition: rectangular_quadmesh.template.cc:805
oomph::RectangularQuadMesh::build_mesh
void build_mesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
Definition: rectangular_quadmesh.template.cc:48
rectangular_quadmesh.template.h