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