simple_cubic_mesh.template.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_SIMPLE_CUBIC_MESH_TEMPLATE_CC
27 #define OOMPH_SIMPLE_CUBIC_MESH_TEMPLATE_CC
28 
29 // OOMPH-LIB Header
30 #include "simple_cubic_mesh.template.h"
31 
32 namespace oomph
33 {
34  //=========================================================================
35  /// Generic mesh construction. This function contains all the details of
36  /// the mesh generation process, including all the tedious loops, counting
37  /// spacing and boundary functions.
38  //========================================================================
39  template<class ELEMENT>
40  void SimpleCubicMesh<ELEMENT>::build_mesh(TimeStepper* time_stepper_pt)
41  {
42  // Mesh can only be built with 3D Qelements.
43  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
44 
45  if ((Nx == 1) || (Ny == 1) || (Nz == 1))
46  {
47  std::ostringstream error_message;
48  error_message << "SimpleCubicMesh needs at least two elements in each,\n"
49  << "coordinate direction. You have specified \n"
50  << "Nx=" << Nx << "; Ny=" << Ny << "; Nz=" << Nz
51  << std::endl;
52  throw OomphLibError(
53  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
54  }
55 
56  // Set the number of boundaries
57  set_nboundary(6);
58 
59  // Allocate the store for the elements
60  Element_pt.resize(Nx * Ny * Nz);
61  // Create first element
62  unsigned element_num = 0;
63  Element_pt[element_num] = new ELEMENT;
64 
65  // Read out the number of linear points in the element
66  unsigned n_p = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
67 
68  // Can now allocate the store for the nodes
69  Node_pt.resize((1 + (n_p - 1) * Nx) * (1 + (n_p - 1) * Ny) *
70  (1 + (n_p - 1) * Nz));
71 
72  // Set up geometrical data
73  //------------------------
74 
75  unsigned long node_count = 0;
76 
77  // Set the length of the elments
78  double el_length[3] = {(Xmax - Xmin) / double(Nx),
79  (Ymax - Ymin) / double(Ny),
80  (Zmax - Zmin) / double(Nz)};
81 
82  // Storage for local coordinate in element
83  Vector<double> s_fraction;
84 
85  // Now assign the topology
86  // Boundaries are numbered:
87  // 0 is at the bottom
88  // 1 2 3 4 from the front proceeding anticlockwise
89  // 5 is at the top
90  // Pinned value are denoted by an integer value 1
91  // Thus if a node is on two boundaries, ORing the values of the
92  // boundary conditions will give the most restrictive case (pinning)
93 
94  // FIRST ELEMENT (lower left corner at z = 0 plane )
95  //----------------------------------
96 
97  // Set the corner node
98  // Storage for the local node number
99  unsigned local_node_num = 0;
100  // Allocate memory for the node
101  Node_pt[node_count] =
102  finite_element_pt(element_num)
103  ->construct_boundary_node(local_node_num, time_stepper_pt);
104  // Set the pointer from the element
105  finite_element_pt(element_num)->node_pt(local_node_num) =
106  Node_pt[node_count];
107 
108  // Set the position of the node
109  Node_pt[node_count]->x(0) = Xmin;
110  Node_pt[node_count]->x(1) = Ymin;
111  Node_pt[node_count]->x(2) = Zmin;
112 
113  // Add the node to boundaries 0,1 and 4
114  add_boundary_node(0, Node_pt[node_count]);
115  add_boundary_node(1, Node_pt[node_count]);
116  add_boundary_node(4, Node_pt[node_count]);
117  // Increment the node number
118  node_count++;
119 
120  // Loop over the other nodes in the first row in the x direction in
121  // the z=0 plane
122  for (unsigned l2 = 1; l2 < n_p; l2++)
123  {
124  // Set the local node number
125  local_node_num = l2;
126 
127  // Allocate memory for the nodes
128  Node_pt[node_count] =
129  finite_element_pt(element_num)
130  ->construct_boundary_node(local_node_num, time_stepper_pt);
131  // Set the pointer from the element
132  finite_element_pt(element_num)->node_pt(local_node_num) =
133  Node_pt[node_count];
134 
135  // Get the local fraction of the node
136  finite_element_pt(element_num)
137  ->local_fraction_of_node(local_node_num, s_fraction);
138 
139  // Set the position of the node
140  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
141  Node_pt[node_count]->x(1) = Ymin;
142  Node_pt[node_count]->x(2) = Zmin;
143 
144  // Add the node to the boundary 0 and 1
145  add_boundary_node(0, Node_pt[node_count]);
146  add_boundary_node(1, Node_pt[node_count]);
147  // Increment the node number
148  node_count++;
149  }
150 
151  // Loop over the other node columns in the z = 0 plane
152  for (unsigned l1 = 1; l1 < n_p; l1++)
153  {
154  // Set the local node number
155  local_node_num = l1 * n_p;
156 
157  // Allocate memory for the nodes
158  Node_pt[node_count] =
159  finite_element_pt(element_num)
160  ->construct_boundary_node(local_node_num, time_stepper_pt);
161  // Set the pointer from the element
162  finite_element_pt(element_num)->node_pt(local_node_num) =
163  Node_pt[node_count];
164 
165  // Get the fractional position
166  finite_element_pt(element_num)
167  ->local_fraction_of_node(local_node_num, s_fraction);
168 
169  // Set the position of the first node of the row
170  //(with boundaries with 0 and 4)
171  Node_pt[node_count]->x(0) = Xmin;
172  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
173  Node_pt[node_count]->x(2) = Zmin;
174 
175  // Add the node to the boundaries 0 and 4
176  add_boundary_node(4, Node_pt[node_count]);
177  add_boundary_node(0, Node_pt[node_count]);
178  // Increment the node number
179  node_count++;
180 
181  // Loop over the other nodes in the row
182  for (unsigned l2 = 1; l2 < n_p; l2++)
183  {
184  // Set the local node number
185  local_node_num = l1 * n_p + l2;
186 
187  // Allocate the memory for the node
188  Node_pt[node_count] =
189  finite_element_pt(element_num)
190  ->construct_boundary_node(local_node_num, time_stepper_pt);
191  // Set the pointer from the element
192  finite_element_pt(element_num)->node_pt(local_node_num) =
193  Node_pt[node_count];
194 
195  // Get the fractional position
196  finite_element_pt(element_num)
197  ->local_fraction_of_node(local_node_num, s_fraction);
198 
199  // Set the position of the node
200  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
201  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
202  Node_pt[node_count]->x(2) = Zmin;
203 
204  // Add the node to the boundary 0
205  add_boundary_node(0, Node_pt[node_count]);
206  // Increment the node number
207  node_count++;
208  }
209  }
210 
211 
212  //---------------------------------------------------------------------
213 
214  // Loop over the other node columns in the z direction
215  for (unsigned l3 = 1; l3 < n_p; l3++)
216  {
217  // Set the local node number
218  local_node_num = n_p * n_p * l3;
219 
220  // Allocate memory for the node
221  Node_pt[node_count] =
222  finite_element_pt(element_num)
223  ->construct_boundary_node(local_node_num, time_stepper_pt);
224  // Set the pointer from the element
225  finite_element_pt(element_num)->node_pt(local_node_num) =
226  Node_pt[node_count];
227 
228  // Get the fractional position of the node
229  finite_element_pt(element_num)
230  ->local_fraction_of_node(local_node_num, s_fraction);
231 
232  // Set the position of the node
233  Node_pt[node_count]->x(0) = Xmin;
234  Node_pt[node_count]->x(1) = Ymin;
235  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
236 
237  // Add the node to boundaries 1 and 4
238  add_boundary_node(1, Node_pt[node_count]);
239  add_boundary_node(4, Node_pt[node_count]);
240  // Increment the node number
241  node_count++;
242 
243  // Loop over the other nodes in the first row in the x direction
244  for (unsigned l2 = 1; l2 < n_p; l2++)
245  {
246  // Set the local node number
247  local_node_num = l2 + n_p * n_p * l3;
248 
249  // Allocate memory for the nodes
250  Node_pt[node_count] =
251  finite_element_pt(element_num)
252  ->construct_boundary_node(local_node_num, time_stepper_pt);
253  // Set the pointer from the element
254  finite_element_pt(element_num)->node_pt(local_node_num) =
255  Node_pt[node_count];
256 
257  // Get the fractional position of the node
258  finite_element_pt(element_num)
259  ->local_fraction_of_node(local_node_num, s_fraction);
260 
261  // Set the position of the node
262  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
263  Node_pt[node_count]->x(1) = Ymin;
264  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
265 
266  // Add the node to the boundary 1
267  add_boundary_node(1, Node_pt[node_count]);
268  // Increment the node number
269  node_count++;
270  }
271 
272  // Loop over the other node columns
273  for (unsigned l1 = 1; l1 < n_p; l1++)
274  {
275  // Set the local node number
276  local_node_num = l1 * n_p + n_p * n_p * l3;
277 
278  // Allocate memory for the nodes
279  Node_pt[node_count] =
280  finite_element_pt(element_num)
281  ->construct_boundary_node(local_node_num, time_stepper_pt);
282  // Set the pointer from the element
283  finite_element_pt(element_num)->node_pt(local_node_num) =
284  Node_pt[node_count];
285 
286  // Get the fractional position of the node
287  finite_element_pt(element_num)
288  ->local_fraction_of_node(local_node_num, s_fraction);
289 
290  // Set the position of the first node of the row (with boundary 4)
291  Node_pt[node_count]->x(0) = Xmin;
292  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
293  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
294 
295  // Add the node to the boundary 4
296  add_boundary_node(4, Node_pt[node_count]);
297  // Increment the node number
298  node_count++;
299 
300  // Loop over the other nodes in the row
301  for (unsigned l2 = 1; l2 < n_p; l2++)
302  {
303  // Set the local node number
304  local_node_num = l2 + l1 * n_p + n_p * n_p * l3;
305 
306  // Allocate the memory for the node
307  Node_pt[node_count] =
308  finite_element_pt(element_num)
309  ->construct_node(local_node_num, time_stepper_pt);
310  // Set the pointer from the element
311  finite_element_pt(element_num)->node_pt(local_node_num) =
312  Node_pt[node_count];
313 
314  // Get the fractional position of the node
315  finite_element_pt(element_num)
316  ->local_fraction_of_node(local_node_num, s_fraction);
317 
318  // Set the position of the node
319  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
320  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
321  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
322 
323  // No boundary
324 
325  // Increment the node number
326  node_count++;
327  }
328  }
329  }
330 
331 
332  //----------------------------------------------------------------------
333 
334  // CENTRE OF FIRST ROW OF ELEMENTS (PLANE Z = 0)
335  //--------------------------------
336 
337  // Now loop over the first row of elements, apart from final element
338  for (unsigned j = 1; j < (Nx - 1); j++)
339  {
340  // Allocate memory for new element
341  element_num = j;
342  Element_pt[element_num] = new ELEMENT;
343 
344  // We will begin with all the nodes at the plane z = 0
345 
346  // Do first row of nodes
347 
348  // First column of nodes is same as neighbouring element
349  finite_element_pt(element_num)->node_pt(0) =
350  finite_element_pt(element_num - 1)->node_pt((n_p - 1));
351 
352  // New nodes for other columns
353  for (unsigned l2 = 1; l2 < n_p; l2++)
354  {
355  // Set the local node number
356  local_node_num = l2;
357 
358  // Allocate memory for the nodes
359  Node_pt[node_count] =
360  finite_element_pt(element_num)
361  ->construct_boundary_node(local_node_num, time_stepper_pt);
362  // Set the pointer from the element
363  finite_element_pt(element_num)->node_pt(local_node_num) =
364  Node_pt[node_count];
365 
366  // Get the fractional position of the node
367  finite_element_pt(element_num)
368  ->local_fraction_of_node(local_node_num, s_fraction);
369 
370  // Set the position of the node
371  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
372  Node_pt[node_count]->x(1) = Ymin;
373  Node_pt[node_count]->x(2) = Zmin;
374 
375  // Add the node to the boundaries 0 an 1
376  add_boundary_node(0, Node_pt[node_count]);
377  add_boundary_node(1, Node_pt[node_count]);
378  // Increment the node number
379  node_count++;
380  }
381 
382  // Do the rest of the nodes at the plane z = 0
383  for (unsigned l1 = 1; l1 < n_p; l1++)
384  {
385  // First column of nodes is same as neighbouring element
386  finite_element_pt(element_num)->node_pt(l1 * n_p) =
387  finite_element_pt(element_num - 1)->node_pt(l1 * n_p + (n_p - 1));
388 
389  // New nodes for other columns
390  for (unsigned l2 = 1; l2 < n_p; l2++)
391  {
392  // Set the local node number
393  local_node_num = l2 + l1 * n_p;
394 
395  // Allocate memory for the nodes
396  Node_pt[node_count] =
397  finite_element_pt(element_num)
398  ->construct_boundary_node(local_node_num, time_stepper_pt);
399  // Set the pointer from the element
400  finite_element_pt(element_num)->node_pt(local_node_num) =
401  Node_pt[node_count];
402 
403  // Get the fractional position of the node
404  finite_element_pt(element_num)
405  ->local_fraction_of_node(local_node_num, s_fraction);
406 
407  // Set the position of the node
408  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
409  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
410  Node_pt[node_count]->x(2) = Zmin;
411 
412  // Add the node to the Boundary
413  add_boundary_node(0, Node_pt[node_count]);
414  // Increment the node number
415  node_count++;
416  }
417  }
418 
419  // Loop over the other node columns in the z direction
420  for (unsigned l3 = 1; l3 < n_p; l3++)
421  {
422  // First column of nodes is same as neighbouring element
423  finite_element_pt(j)->node_pt(l3 * n_p * n_p) =
424  finite_element_pt(j - 1)->node_pt(l3 * n_p * n_p + (n_p - 1));
425 
426  // New nodes for other columns
427  for (unsigned l2 = 1; l2 < n_p; l2++)
428  {
429  // Set the local node number
430  local_node_num = l2 + l3 * n_p * n_p;
431 
432  // Allocate memory for the nodes
433  Node_pt[node_count] =
434  finite_element_pt(element_num)
435  ->construct_boundary_node(local_node_num, time_stepper_pt);
436  // Set the pointer from the element
437  finite_element_pt(element_num)->node_pt(local_node_num) =
438  Node_pt[node_count];
439 
440  // Get the fractional position of the node
441  finite_element_pt(element_num)
442  ->local_fraction_of_node(local_node_num, s_fraction);
443 
444  // Set the position of the node
445  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
446  Node_pt[node_count]->x(1) = Ymin;
447  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
448 
449  // Add the node to the boundary 1
450  add_boundary_node(1, Node_pt[node_count]);
451  // Increment the node number
452  node_count++;
453  }
454 
455  // Do the rest of the nodes
456  for (unsigned l1 = 1; l1 < n_p; l1++)
457  {
458  // First column of nodes is same as neighbouring element
459  finite_element_pt(j)->node_pt(l1 * n_p + l3 * n_p * n_p) =
460  finite_element_pt(j - 1)->node_pt(l1 * n_p + (n_p - 1) +
461  l3 * n_p * n_p);
462 
463  // New nodes for other columns
464  for (unsigned l2 = 1; l2 < n_p; l2++)
465  {
466  // Set the local node number
467  local_node_num = l2 + l1 * n_p + l3 * n_p * n_p;
468 
469  // Allocate memory for the nodes
470  Node_pt[node_count] =
471  finite_element_pt(element_num)
472  ->construct_node(local_node_num, time_stepper_pt);
473  // Set the pointer from the element
474  finite_element_pt(element_num)->node_pt(local_node_num) =
475  Node_pt[node_count];
476 
477  // Get the fractional position of the node
478  finite_element_pt(element_num)
479  ->local_fraction_of_node(local_node_num, s_fraction);
480 
481  // Set the position of the node
482  Node_pt[node_count]->x(0) =
483  Xmin + el_length[0] * (j + s_fraction[0]);
484  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
485  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
486 
487  // No boundaries
488 
489  // Increment the node number
490  node_count++;
491  }
492  }
493  }
494  }
495 
496  // MY FINAL ELEMENT IN FIRST ROW (lower right corner)
497  //-----------------------------------------------
498 
499  // Allocate memory for new element
500  element_num = Nx - 1;
501  Element_pt[element_num] = new ELEMENT;
502 
503  // We will begin with all the nodes at the plane z = 0
504 
505  // Do first row of nodes
506 
507  // First node is same as neighbouring element
508  finite_element_pt(element_num)->node_pt(0) =
509  finite_element_pt(element_num - 1)->node_pt((n_p - 1));
510 
511  // New nodes for other columns
512  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
513  {
514  // Set the local node number
515  local_node_num = l2;
516 
517  // Allocate memory for the nodes
518  Node_pt[node_count] =
519  finite_element_pt(element_num)
520  ->construct_boundary_node(local_node_num, time_stepper_pt);
521  // Set the pointer from the element
522  finite_element_pt(element_num)->node_pt(local_node_num) =
523  Node_pt[node_count];
524 
525  // Get the fractional position of the node
526  finite_element_pt(element_num)
527  ->local_fraction_of_node(local_node_num, s_fraction);
528 
529  // Set the position of the node
530  Node_pt[node_count]->x(0) =
531  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
532  Node_pt[node_count]->x(1) = Ymin;
533  Node_pt[node_count]->x(2) = Zmin;
534 
535  // Add the node to the boundaries 0 an 1
536  add_boundary_node(0, Node_pt[node_count]);
537  add_boundary_node(1, Node_pt[node_count]);
538  // Increment the node number
539  node_count++;
540  }
541 
542  // Last node (corner)
543  // Set the local node number
544  local_node_num = n_p - 1;
545 
546  // Allocate memory for the node
547  Node_pt[node_count] =
548  finite_element_pt(element_num)
549  ->construct_boundary_node(local_node_num, time_stepper_pt);
550  // Set the pointer from the element
551  finite_element_pt(element_num)->node_pt(local_node_num) =
552  Node_pt[node_count];
553 
554  // Get the fractional position of the node
555  finite_element_pt(element_num)
556  ->local_fraction_of_node(local_node_num, s_fraction);
557 
558  // Set the position of the node
559  Node_pt[node_count]->x(0) = Xmax;
560  Node_pt[node_count]->x(1) = Ymin;
561  Node_pt[node_count]->x(2) = Zmin;
562 
563  // Add the node to the boundaries
564  add_boundary_node(0, Node_pt[node_count]);
565  add_boundary_node(1, Node_pt[node_count]);
566  add_boundary_node(2, Node_pt[node_count]);
567  // Increment the node number
568  node_count++;
569 
570  // Do the middle nodes at the plane z = 0
571  for (unsigned l1 = 1; l1 < n_p; l1++)
572  {
573  // First column of nodes is same as neighbouring element
574  finite_element_pt(element_num)->node_pt(l1 * n_p) =
575  finite_element_pt(element_num - 1)->node_pt(l1 * n_p + (n_p - 1));
576 
577  // New nodes for other columns
578  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
579  {
580  // Set the local node number
581  local_node_num = l2 + l1 * n_p;
582 
583  // Allocate memory for the nodes
584  Node_pt[node_count] =
585  finite_element_pt(element_num)
586  ->construct_boundary_node(local_node_num, time_stepper_pt);
587  // Set the pointer from the element
588  finite_element_pt(element_num)->node_pt(local_node_num) =
589  Node_pt[node_count];
590 
591  // Get the fractional position of the node
592  finite_element_pt(element_num)
593  ->local_fraction_of_node(local_node_num, s_fraction);
594 
595  // Set the position of the node
596  Node_pt[node_count]->x(0) =
597  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
598  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
599  Node_pt[node_count]->x(2) = Zmin;
600 
601  // Add the node to the boundary 0
602  add_boundary_node(0, Node_pt[node_count]);
603  // Increment the node number
604  node_count++;
605  }
606 
607  // New node for final column
608  // Set the local node number
609  local_node_num = l1 * n_p + (n_p - 1);
610 
611  // Allocate memory for node
612  Node_pt[node_count] =
613  finite_element_pt(element_num)
614  ->construct_boundary_node(local_node_num, time_stepper_pt);
615  // Set the pointer from the element
616  finite_element_pt(element_num)->node_pt(local_node_num) =
617  Node_pt[node_count];
618 
619  // Get the fractional position of the node
620  finite_element_pt(element_num)
621  ->local_fraction_of_node(local_node_num, s_fraction);
622 
623  // Set the position of the node
624  Node_pt[node_count]->x(0) = Xmax;
625  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
626  Node_pt[node_count]->x(2) = Zmin;
627 
628  // Add the node to the boundaries 0 and 2
629  add_boundary_node(2, Node_pt[node_count]);
630  add_boundary_node(0, Node_pt[node_count]);
631  // Increment the node number
632  node_count++;
633  }
634 
635  // Loop over the other node columns in the z direction
636  for (unsigned l3 = 1; l3 < n_p; l3++)
637  {
638  // First column of nodes is same as neighbouring element
639  finite_element_pt(element_num)->node_pt(l3 * n_p * n_p) =
640  finite_element_pt(element_num - 1)->node_pt(l3 * n_p * n_p + (n_p - 1));
641 
642  // New nodes for other rows (y=0)
643  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
644  {
645  // Set the local node number
646  local_node_num = l2 + l3 * n_p * n_p;
647 
648  // Allocate memory for the nodes
649  Node_pt[node_count] =
650  finite_element_pt(element_num)
651  ->construct_boundary_node(local_node_num, time_stepper_pt);
652  // Set the pointer from the element
653  finite_element_pt(element_num)->node_pt(local_node_num) =
654  Node_pt[node_count];
655 
656  // Get the fractional position of the node
657  finite_element_pt(element_num)
658  ->local_fraction_of_node(local_node_num, s_fraction);
659 
660  // Set the position of the node
661  Node_pt[node_count]->x(0) =
662  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
663  Node_pt[node_count]->x(1) = Ymin;
664  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
665 
666  // Add the node to the boundary 1
667  add_boundary_node(1, Node_pt[node_count]);
668  // Increment the node number
669  node_count++;
670  }
671 
672  // Last node of the row (y=0)
673  // Set the local node number
674  local_node_num = n_p - 1 + l3 * n_p * n_p;
675 
676  // Allocate memory for the nodes
677  Node_pt[node_count] =
678  finite_element_pt(element_num)
679  ->construct_boundary_node(local_node_num, time_stepper_pt);
680  // Set the pointer from the element
681  finite_element_pt(element_num)->node_pt(local_node_num) =
682  Node_pt[node_count];
683 
684  // Get the fractional position of the node
685  finite_element_pt(element_num)
686  ->local_fraction_of_node(local_node_num, s_fraction);
687 
688  // Set the position of the node
689  Node_pt[node_count]->x(0) = Xmax;
690  Node_pt[node_count]->x(1) = Ymin;
691  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
692 
693  // Add the node to the boundary 1 and 2
694  add_boundary_node(1, Node_pt[node_count]);
695  add_boundary_node(2, Node_pt[node_count]);
696  // Increment the node number
697  node_count++;
698 
699  // Do the rest of the nodes
700  for (unsigned l1 = 1; l1 < n_p; l1++)
701  {
702  // First column of nodes is same as neighbouring element
703  finite_element_pt(element_num)->node_pt(l1 * n_p + l3 * n_p * n_p) =
704  finite_element_pt(element_num - 1)
705  ->node_pt(l1 * n_p + (n_p - 1) + l3 * n_p * n_p);
706 
707  // New nodes for other columns
708  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
709  {
710  // Set the local node number
711  local_node_num = l2 + l1 * n_p + l3 * n_p * n_p;
712 
713  // Allocate memory for the nodes
714  Node_pt[node_count] =
715  finite_element_pt(element_num)
716  ->construct_node(local_node_num, time_stepper_pt);
717  // Set the pointer from the element
718  finite_element_pt(element_num)->node_pt(local_node_num) =
719  Node_pt[node_count];
720 
721  // Get the fractional position of the node
722  finite_element_pt(element_num)
723  ->local_fraction_of_node(local_node_num, s_fraction);
724 
725  // Set the position of the node
726  Node_pt[node_count]->x(0) =
727  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
728  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
729  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
730 
731  // No boundaries
732 
733  // Increment the node number
734  node_count++;
735  }
736 
737  // Last nodes (at the surface x = Lx)
738  // Set the local nod number
739  local_node_num = l1 * n_p + (n_p - 1) + l3 * n_p * n_p;
740  // Allocate memory for the nodes
741  Node_pt[node_count] =
742  finite_element_pt(element_num)
743  ->construct_boundary_node(local_node_num, time_stepper_pt);
744  // Set the pointer from the element
745  finite_element_pt(element_num)->node_pt(local_node_num) =
746  Node_pt[node_count];
747 
748  // Get the fractional position of the node
749  finite_element_pt(element_num)
750  ->local_fraction_of_node(local_node_num, s_fraction);
751 
752  // Set the position of the node
753  Node_pt[node_count]->x(0) = Xmax;
754  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
755  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
756 
757  // Add the node to the boundary 2
758  add_boundary_node(2, Node_pt[node_count]);
759  // Increment the node number
760  node_count++;
761  }
762  }
763 
764 
765  // ALL CENTRAL ELEMENT ROWS (WE ARE STILL IN THE LAYER z=0)
766  //------------------------
767 
768  // Loop over remaining element rows
769  for (unsigned i = 1; i < (Ny - 1); i++)
770  {
771  // Set the first element in the row
772 
773  // Allocate memory for element (z = 0) (x =0)
774  element_num = Nx * i;
775  Element_pt[element_num] = new ELEMENT;
776 
777  // The first row of nodes is copied from the element below (at z=0)
778  for (unsigned l2 = 0; l2 < n_p; l2++)
779  {
780  finite_element_pt(element_num)->node_pt(l2) =
781  finite_element_pt(element_num - Nx)->node_pt((n_p - 1) * n_p + l2);
782  }
783 
784  // Other rows are new nodes
785  for (unsigned l1 = 1; l1 < n_p; l1++)
786  {
787  // First column of nodes
788  // Set the local node number
789  local_node_num = l1 * n_p;
790 
791  // Allocate memory for the fist node in the x direction (x=0)
792  Node_pt[node_count] =
793  finite_element_pt(element_num)
794  ->construct_boundary_node(local_node_num, time_stepper_pt);
795  // Set the pointer from the element
796  finite_element_pt(element_num)->node_pt(local_node_num) =
797  Node_pt[node_count];
798 
799  // Get the fractional position of the node
800  finite_element_pt(element_num)
801  ->local_fraction_of_node(local_node_num, s_fraction);
802 
803  // Set the position of the node
804  Node_pt[node_count]->x(0) = Xmin;
805  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
806  Node_pt[node_count]->x(2) = Zmin;
807 
808  // Add the node to the boundaries 4 and 0
809  add_boundary_node(0, Node_pt[node_count]);
810  add_boundary_node(4, Node_pt[node_count]);
811  // Increment the node number
812  node_count++;
813 
814  // Now do the other columns
815  for (unsigned l2 = 1; l2 < n_p; l2++)
816  {
817  // Set the local node number
818  local_node_num = l2 + l1 * n_p;
819 
820  // Allocate memory for node
821  Node_pt[node_count] =
822  finite_element_pt(element_num)
823  ->construct_boundary_node(local_node_num, time_stepper_pt);
824  // Set the pointer from the element
825  finite_element_pt(element_num)->node_pt(local_node_num) =
826  Node_pt[node_count];
827 
828  // Get the fractional position of the node
829  finite_element_pt(element_num)
830  ->local_fraction_of_node(local_node_num, s_fraction);
831 
832 
833  // Set the position of the node
834  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
835  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
836  Node_pt[node_count]->x(2) = Zmin;
837 
838  // Add the node to the boundary and 0
839  add_boundary_node(0, Node_pt[node_count]);
840  // Increment the node number
841  node_count++;
842  }
843  }
844 
845  // As always we extend now this element to the third direction
846  for (unsigned l3 = 1; l3 < n_p; l3++)
847  {
848  // The first row of nodes is copied from the element below (at z=0)
849  for (unsigned l2 = 0; l2 < n_p; l2++)
850  {
851  finite_element_pt(element_num)->node_pt(l2 + l3 * n_p * n_p) =
852  finite_element_pt(element_num - Nx)
853  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
854  }
855 
856  // Other rows are new nodes (first the nodes for which x=0)
857  for (unsigned l1 = 1; l1 < n_p; l1++)
858  {
859  // First column of nodes
860  // Set the local node number
861  local_node_num = l1 * n_p + l3 * n_p * n_p;
862 
863  // Allocate memory for node
864  Node_pt[node_count] =
865  finite_element_pt(element_num)
866  ->construct_boundary_node(local_node_num, time_stepper_pt);
867  // Set the pointer from the element
868  finite_element_pt(element_num)->node_pt(local_node_num) =
869  Node_pt[node_count];
870 
871  // Get the fractional position of the node
872  finite_element_pt(element_num)
873  ->local_fraction_of_node(local_node_num, s_fraction);
874 
875  // Set the position of the node
876  Node_pt[node_count]->x(0) = Xmin;
877  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
878  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
879 
880  // Add the node to the boundary 4
881  add_boundary_node(4, Node_pt[node_count]);
882 
883  // Increment the node number
884  node_count++;
885 
886  // Now do the other columns (we extend to the rest of the nodes)
887  for (unsigned l2 = 1; l2 < n_p; l2++)
888  {
889  // Set the local node number
890  local_node_num = l2 + l1 * n_p + n_p * n_p * l3;
891 
892  // Allocate memory for node
893  Node_pt[node_count] =
894  finite_element_pt(element_num)
895  ->construct_node(local_node_num, time_stepper_pt);
896  // Set the pointer from the element
897  finite_element_pt(element_num)->node_pt(local_node_num) =
898  Node_pt[node_count];
899 
900  // Get the fractional position of the node
901  finite_element_pt(element_num)
902  ->local_fraction_of_node(local_node_num, s_fraction);
903 
904  // Set the position of the node
905  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
906  Node_pt[node_count]->x(1) =
907  Ymin + el_length[1] * (i + s_fraction[1]);
908  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
909 
910  // No boundaries
911 
912  // Increment the node number
913  node_count++;
914  }
915  }
916  }
917 
918  // Now loop over the rest of the elements in the row, apart from the last
919  for (unsigned j = 1; j < (Nx - 1); j++)
920  {
921  // Allocate memory for new element
922  element_num = Nx * i + j;
923  Element_pt[element_num] = new ELEMENT;
924 
925  // The first row is copied from the lower element
926  for (unsigned l2 = 0; l2 < n_p; l2++)
927  {
928  finite_element_pt(element_num)->node_pt(l2) =
929  finite_element_pt(element_num - Nx)->node_pt((n_p - 1) * n_p + l2);
930  }
931 
932  for (unsigned l1 = 1; l1 < n_p; l1++)
933  {
934  // First column is same as neighbouring element
935  finite_element_pt(element_num)->node_pt(l1 * n_p) =
936  finite_element_pt(element_num - 1)->node_pt(l1 * n_p + (n_p - 1));
937 
938  // New nodes for other columns
939  for (unsigned l2 = 1; l2 < n_p; l2++)
940  {
941  // Set local node number
942  local_node_num = l1 * n_p + l2;
943 
944  // Allocate memory for the nodes
945  Node_pt[node_count] =
946  finite_element_pt(element_num)
947  ->construct_boundary_node(local_node_num, time_stepper_pt);
948  // Set the pointer
949  finite_element_pt(element_num)->node_pt(local_node_num) =
950  Node_pt[node_count];
951 
952  // Get the fractional position of the node
953  finite_element_pt(element_num)
954  ->local_fraction_of_node(local_node_num, s_fraction);
955 
956  // Set the position of the node
957  Node_pt[node_count]->x(0) =
958  Xmin + el_length[0] * (j + s_fraction[0]);
959  Node_pt[node_count]->x(1) =
960  Ymin + el_length[1] * (i + s_fraction[1]);
961  Node_pt[node_count]->x(2) = Zmin;
962 
963  // Add the node to the boundary and 0
964  add_boundary_node(0, Node_pt[node_count]);
965  // Increment the node number
966  node_count++;
967  }
968  }
969 
970  // We extend to the third dimension
971  for (unsigned l3 = 1; l3 < n_p; l3++)
972  {
973  // The first row is copied from the lower element
974 
975  for (unsigned l2 = 0; l2 < n_p; l2++)
976  {
977  finite_element_pt(element_num)->node_pt(l2 + l3 * n_p * n_p) =
978  finite_element_pt(element_num - Nx)
979  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
980  }
981 
982  for (unsigned l1 = 1; l1 < n_p; l1++)
983  {
984  // First column is same as neighbouring element
985  finite_element_pt(element_num)->node_pt(l1 * n_p + l3 * n_p * n_p) =
986  finite_element_pt(element_num - 1)
987  ->node_pt(l1 * n_p + l3 * n_p * n_p + (n_p - 1));
988 
989  // New nodes for other columns
990  for (unsigned l2 = 1; l2 < n_p; l2++)
991  {
992  // Set the local node number
993  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
994 
995  // Allocate memory for the nodes
996  Node_pt[node_count] =
997  finite_element_pt(element_num)
998  ->construct_node(local_node_num, time_stepper_pt);
999  // Set the pointer
1000  finite_element_pt(element_num)->node_pt(local_node_num) =
1001  Node_pt[node_count];
1002 
1003  // Get the fractional position of the node
1004  finite_element_pt(element_num)
1005  ->local_fraction_of_node(local_node_num, s_fraction);
1006 
1007  // Set the position of the node
1008  Node_pt[node_count]->x(0) =
1009  Xmin + el_length[0] * (j + s_fraction[0]);
1010  Node_pt[node_count]->x(1) =
1011  Ymin + el_length[1] * (i + s_fraction[1]);
1012  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1013 
1014  // No boundaries
1015 
1016  // Increment the node number
1017  node_count++;
1018  }
1019  }
1020  }
1021 
1022  } // End of loop over elements in row
1023 
1024 
1025  // Do final element in row
1026 
1027  // Allocate memory for element
1028  element_num = Nx * i + Nx - 1;
1029  Element_pt[element_num] = new ELEMENT;
1030 
1031  // We begin with z =0
1032  // The first row is copied from the lower element
1033  for (unsigned l2 = 0; l2 < n_p; l2++)
1034  {
1035  finite_element_pt(element_num)->node_pt(l2) =
1036  finite_element_pt(element_num - Nx)->node_pt((n_p - 1) * n_p + l2);
1037  }
1038 
1039  for (unsigned l1 = 1; l1 < n_p; l1++)
1040  {
1041  // First column is same as neighbouring element
1042  finite_element_pt(element_num)->node_pt(l1 * n_p) =
1043  finite_element_pt(element_num - 1)->node_pt(l1 * n_p + (n_p - 1));
1044 
1045  // Middle nodes
1046  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
1047  {
1048  // Set local node number
1049  local_node_num = l1 * n_p + l2;
1050 
1051  // Allocate memory for node
1052  Node_pt[node_count] =
1053  finite_element_pt(element_num)
1054  ->construct_boundary_node(local_node_num, time_stepper_pt);
1055  // Set the pointer
1056  finite_element_pt(element_num)->node_pt(local_node_num) =
1057  Node_pt[node_count];
1058 
1059  // Get the fractional position of the node
1060  finite_element_pt(element_num)
1061  ->local_fraction_of_node(local_node_num, s_fraction);
1062 
1063  // Set the position of the node
1064  Node_pt[node_count]->x(0) =
1065  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
1066  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
1067  Node_pt[node_count]->x(2) = Zmin;
1068 
1069  // Add the node to the boundary and 0
1070  add_boundary_node(0, Node_pt[node_count]);
1071 
1072  // Increment the node number
1073  node_count++;
1074  }
1075 
1076  // Final node
1077 
1078  // Set the local node number
1079  local_node_num = l1 * n_p + (n_p - 1);
1080 
1081  // Allocate memory for node
1082  Node_pt[node_count] =
1083  finite_element_pt(element_num)
1084  ->construct_boundary_node(local_node_num, time_stepper_pt);
1085  // Set the pointer
1086  finite_element_pt(element_num)->node_pt(local_node_num) =
1087  Node_pt[node_count];
1088 
1089  // Get the fractional position of the node
1090  finite_element_pt(element_num)
1091  ->local_fraction_of_node(local_node_num, s_fraction);
1092 
1093 
1094  // Set the position of the node
1095  Node_pt[node_count]->x(0) = Xmax;
1096  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
1097  Node_pt[node_count]->x(2) = Zmin;
1098 
1099  // Add the node to the boundaries - and 1
1100  add_boundary_node(0, Node_pt[node_count]);
1101  add_boundary_node(2, Node_pt[node_count]);
1102 
1103  // Increment the node number
1104  node_count++;
1105 
1106  } // End of loop over rows of nodes in the element
1107 
1108  // We go to the third dimension
1109  for (unsigned l3 = 1; l3 < n_p; l3++)
1110  {
1111  // The first row is copied from the lower element
1112  for (unsigned l2 = 0; l2 < n_p; l2++)
1113  {
1114  finite_element_pt(element_num)->node_pt(l2 + l3 * n_p * n_p) =
1115  finite_element_pt(element_num - Nx)
1116  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
1117  }
1118 
1119  for (unsigned l1 = 1; l1 < n_p; l1++)
1120  {
1121  // First column is same as neighbouring element
1122  finite_element_pt(element_num)->node_pt(l1 * n_p + l3 * n_p * n_p) =
1123  finite_element_pt(element_num - 1)
1124  ->node_pt(l1 * n_p + (n_p - 1) + l3 * n_p * n_p);
1125 
1126  // Middle nodes
1127  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
1128  {
1129  // Set the local node number
1130  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
1131 
1132  // Allocate memory for node
1133  Node_pt[node_count] =
1134  finite_element_pt(element_num)
1135  ->construct_node(local_node_num, time_stepper_pt);
1136  // Set the pointer
1137  finite_element_pt(element_num)->node_pt(local_node_num) =
1138  Node_pt[node_count];
1139 
1140  // Get the fractional position of the node
1141  finite_element_pt(element_num)
1142  ->local_fraction_of_node(local_node_num, s_fraction);
1143 
1144  // Set the position of the node
1145  Node_pt[node_count]->x(0) =
1146  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
1147  Node_pt[node_count]->x(1) =
1148  Ymin + el_length[1] * (i + s_fraction[1]);
1149  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1150 
1151  // No boundaries
1152 
1153  // Increment the node number
1154  node_count++;
1155  }
1156 
1157  // Final node
1158  // Set the local node number
1159  local_node_num = l1 * n_p + (n_p - 1) + l3 * n_p * n_p;
1160 
1161  // Allocate memory for node
1162  Node_pt[node_count] =
1163  finite_element_pt(element_num)
1164  ->construct_boundary_node(local_node_num, time_stepper_pt);
1165  // Set the pointer
1166  finite_element_pt(element_num)->node_pt(local_node_num) =
1167  Node_pt[node_count];
1168 
1169  // Get the fractional position of the node
1170  finite_element_pt(element_num)
1171  ->local_fraction_of_node(local_node_num, s_fraction);
1172 
1173  // Set the position of the node
1174  Node_pt[node_count]->x(0) = Xmax;
1175  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
1176  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1177 
1178  // Add the node to the boundary 2
1179  add_boundary_node(2, Node_pt[node_count]);
1180  // Increment the node number
1181  node_count++;
1182  } // End of loop over rows of nodes in the element
1183 
1184 
1185  } // End of the 3dimension loop
1186 
1187 
1188  } // End of loop over rows of elements
1189 
1190 
1191  // FINAL ELEMENT ROW (IN THE z=0 LAYER)
1192  //=================
1193 
1194 
1195  // FIRST ELEMENT IN UPPER ROW (upper left corner)
1196  //----------------------------------------------
1197 
1198  // Allocate memory for element
1199  element_num = Nx * (Ny - 1);
1200  Element_pt[element_num] = new ELEMENT;
1201  // We begin with all the nodes for which z=0
1202  // The first row of nodes is copied from the element below
1203  for (unsigned l2 = 0; l2 < n_p; l2++)
1204  {
1205  finite_element_pt(element_num)->node_pt(l2) =
1206  finite_element_pt(element_num - Nx)->node_pt((n_p - 1) * n_p + l2);
1207  }
1208 
1209  // Second row of nodes
1210  // First column of nodes
1211  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
1212  {
1213  // Set the local node number
1214  local_node_num = n_p * l1;
1215 
1216  // Allocate memory for node
1217  Node_pt[node_count] =
1218  finite_element_pt(element_num)
1219  ->construct_boundary_node(local_node_num, time_stepper_pt);
1220  // Set the pointer from the element
1221  finite_element_pt(element_num)->node_pt(local_node_num) =
1222  Node_pt[node_count];
1223 
1224  // Get the fractional position of the node
1225  finite_element_pt(element_num)
1226  ->local_fraction_of_node(local_node_num, s_fraction);
1227 
1228  // Set the position of the node
1229  Node_pt[node_count]->x(0) = Xmin;
1230  Node_pt[node_count]->x(1) =
1231  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1232  Node_pt[node_count]->x(2) = Zmin;
1233 
1234  // Add the node to the boundaries 4 and 0
1235  add_boundary_node(0, Node_pt[node_count]);
1236  add_boundary_node(4, Node_pt[node_count]);
1237  // Increment the node number
1238  node_count++;
1239 
1240  // Now do the other columns
1241  for (unsigned l2 = 1; l2 < n_p; l2++)
1242  {
1243  // Set the local node number
1244  local_node_num = n_p * l1 + l2;
1245 
1246  // Allocate memory for node
1247  Node_pt[node_count] =
1248  finite_element_pt(element_num)
1249  ->construct_boundary_node(local_node_num, time_stepper_pt);
1250  // Set the pointer from the element
1251  finite_element_pt(element_num)->node_pt(local_node_num) =
1252  Node_pt[node_count];
1253 
1254  // Get the fractional position of the node
1255  finite_element_pt(element_num)
1256  ->local_fraction_of_node(local_node_num, s_fraction);
1257 
1258  // Set the position of the node
1259  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
1260  Node_pt[node_count]->x(1) =
1261  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1262  Node_pt[node_count]->x(2) = Zmin;
1263 
1264  // Add the node to the boundary 0
1265  add_boundary_node(0, Node_pt[node_count]);
1266  // Increment the node number
1267  node_count++;
1268  }
1269  }
1270 
1271  // Final row of nodes
1272  // First column of nodes
1273  // Top left node
1274  // Set local node number
1275  local_node_num = n_p * (n_p - 1);
1276  // Allocate memory for node
1277  Node_pt[node_count] =
1278  finite_element_pt(element_num)
1279  ->construct_boundary_node(local_node_num, time_stepper_pt);
1280  // Set the pointer from the element
1281  finite_element_pt(element_num)->node_pt(local_node_num) =
1282  Node_pt[node_count];
1283 
1284  // Get the fractional position of the node
1285  finite_element_pt(element_num)
1286  ->local_fraction_of_node(local_node_num, s_fraction);
1287 
1288  // Set the position of the node
1289  Node_pt[node_count]->x(0) = Xmin;
1290  Node_pt[node_count]->x(1) = Ymax;
1291  Node_pt[node_count]->x(2) = Zmin;
1292 
1293  // Add the node to the boundaries 0,3 and 4
1294  add_boundary_node(0, Node_pt[node_count]);
1295  add_boundary_node(3, Node_pt[node_count]);
1296  add_boundary_node(4, Node_pt[node_count]);
1297 
1298  // Increment the node number
1299  node_count++;
1300 
1301  // Now do the other columns
1302  for (unsigned l2 = 1; l2 < n_p; l2++)
1303  {
1304  // Set the local node number
1305  local_node_num = n_p * (n_p - 1) + l2;
1306  // Allocate memory for the node
1307  Node_pt[node_count] =
1308  finite_element_pt(element_num)
1309  ->construct_boundary_node(local_node_num, time_stepper_pt);
1310  // Set the pointer from the element
1311  finite_element_pt(element_num)->node_pt(local_node_num) =
1312  Node_pt[node_count];
1313 
1314  // Get the fractional position of the node
1315  finite_element_pt(element_num)
1316  ->local_fraction_of_node(local_node_num, s_fraction);
1317 
1318  // Set the position of the node
1319  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
1320  Node_pt[node_count]->x(1) = Ymax;
1321  Node_pt[node_count]->x(2) = Zmin;
1322 
1323  // Add the node to the boundaries
1324  add_boundary_node(0, Node_pt[node_count]);
1325  add_boundary_node(3, Node_pt[node_count]);
1326  // Increment the node number
1327  node_count++;
1328  }
1329 
1330  // We jump to the third dimension
1331  for (unsigned l3 = 1; l3 < n_p; l3++)
1332  {
1333  // The first row of nodes is copied from the element below
1334  for (unsigned l2 = 0; l2 < n_p; l2++)
1335  {
1336  finite_element_pt(element_num)->node_pt(l2 + l3 * n_p * n_p) =
1337  finite_element_pt(element_num - Nx)
1338  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
1339  }
1340 
1341  // Second row of nodes
1342  // First column of nodes
1343  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
1344  {
1345  // Set the local node number
1346  local_node_num = n_p * l1 + l3 * n_p * n_p;
1347 
1348  // Allocate memory for node
1349  Node_pt[node_count] =
1350  finite_element_pt(element_num)
1351  ->construct_boundary_node(local_node_num, time_stepper_pt);
1352  // Set the pointer from the element
1353  finite_element_pt(element_num)->node_pt(local_node_num) =
1354  Node_pt[node_count];
1355 
1356  // Get the fractional position of the node
1357  finite_element_pt(element_num)
1358  ->local_fraction_of_node(local_node_num, s_fraction);
1359 
1360  // Set the position of the node
1361  Node_pt[node_count]->x(0) = Xmin;
1362  Node_pt[node_count]->x(1) =
1363  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1364  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1365 
1366  // Add the node to the boundary 4
1367  add_boundary_node(4, Node_pt[node_count]);
1368  // Increment the node number
1369  node_count++;
1370 
1371  // Now do the other columns
1372  for (unsigned l2 = 1; l2 < n_p; l2++)
1373  {
1374  // Set local node number
1375  local_node_num = n_p * l1 + l2 + l3 * n_p * n_p;
1376 
1377  // Allocate memory for node
1378  Node_pt[node_count] =
1379  finite_element_pt(element_num)
1380  ->construct_node(local_node_num, time_stepper_pt);
1381  // Set the pointer from the element
1382  finite_element_pt(element_num)->node_pt(local_node_num) =
1383  Node_pt[node_count];
1384 
1385  // Get the fractional position of the node
1386  finite_element_pt(element_num)
1387  ->local_fraction_of_node(local_node_num, s_fraction);
1388 
1389  // Set the position of the node
1390  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
1391  Node_pt[node_count]->x(1) =
1392  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1393  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1394 
1395  // No boundaries
1396 
1397  // Increment the node number
1398  node_count++;
1399  }
1400  }
1401 
1402  // Final row of nodes
1403  // First column of nodes
1404  // Top left node
1405  local_node_num = n_p * (n_p - 1) + l3 * n_p * n_p;
1406  // Allocate memory for node
1407  Node_pt[node_count] =
1408  finite_element_pt(element_num)
1409  ->construct_boundary_node(local_node_num, time_stepper_pt);
1410  // Set the pointer from the element
1411  finite_element_pt(element_num)->node_pt(local_node_num) =
1412  Node_pt[node_count];
1413 
1414  // Get the fractional position of the node
1415  finite_element_pt(element_num)
1416  ->local_fraction_of_node(local_node_num, s_fraction);
1417 
1418  // Set the position of the node
1419  Node_pt[node_count]->x(0) = Xmin;
1420  Node_pt[node_count]->x(1) = Ymax;
1421  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1422 
1423  // Add the node to the boundaries
1424  add_boundary_node(3, Node_pt[node_count]);
1425  add_boundary_node(4, Node_pt[node_count]);
1426 
1427  // Increment the node number
1428  node_count++;
1429 
1430  // Now do the other columns
1431  for (unsigned l2 = 1; l2 < n_p; l2++)
1432  {
1433  local_node_num = n_p * (n_p - 1) + l2 + l3 * n_p * n_p;
1434  // Allocate memory for the node
1435  Node_pt[node_count] =
1436  finite_element_pt(element_num)
1437  ->construct_boundary_node(local_node_num, time_stepper_pt);
1438  // Set the pointer from the element
1439  finite_element_pt(element_num)->node_pt(local_node_num) =
1440  Node_pt[node_count];
1441 
1442  // Get the fractional position of the node
1443  finite_element_pt(element_num)
1444  ->local_fraction_of_node(local_node_num, s_fraction);
1445 
1446  // Set the position of the node
1447  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
1448  Node_pt[node_count]->x(1) = Ymax;
1449  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1450 
1451  // Add the node to the boundary 3
1452  add_boundary_node(3, Node_pt[node_count]);
1453 
1454 
1455  // Increment the node number
1456  node_count++;
1457  }
1458  }
1459 
1460 
1461  // Now loop over the rest of the elements in the row, apart from the last
1462  // (first plane z = 0)
1463  for (unsigned j = 1; j < (Nx - 1); j++)
1464  {
1465  // Allocate memory for element
1466  element_num = Nx * (Ny - 1) + j;
1467  Element_pt[element_num] = new ELEMENT;
1468  // The first row is copied from the lower element
1469  for (unsigned l2 = 0; l2 < n_p; l2++)
1470  {
1471  finite_element_pt(element_num)->node_pt(l2) =
1472  finite_element_pt(element_num - Nx)->node_pt((n_p - 1) * n_p + l2);
1473  }
1474 
1475  // Second rows
1476  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
1477  {
1478  // First column is same as neighbouring element
1479  finite_element_pt(element_num)->node_pt(n_p * l1) =
1480  finite_element_pt(element_num - 1)->node_pt(n_p * l1 + (n_p - 1));
1481 
1482  // New nodes for other columns
1483  for (unsigned l2 = 1; l2 < n_p; l2++)
1484  {
1485  local_node_num = n_p * l1 + l2;
1486  // Allocate memory for the node
1487  Node_pt[node_count] =
1488  finite_element_pt(element_num)
1489  ->construct_boundary_node(local_node_num, time_stepper_pt);
1490 
1491  // Set the pointer
1492  finite_element_pt(element_num)->node_pt(local_node_num) =
1493  Node_pt[node_count];
1494 
1495  // Get the fractional position of the node
1496  finite_element_pt(element_num)
1497  ->local_fraction_of_node(local_node_num, s_fraction);
1498 
1499  // Set the position of the node
1500  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
1501  Node_pt[node_count]->x(1) =
1502  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1503  Node_pt[node_count]->x(2) = Zmin;
1504 
1505  // Add the node to the boundary
1506  add_boundary_node(0, Node_pt[node_count]);
1507 
1508  // Increment the node number
1509  node_count++;
1510  }
1511  }
1512 
1513  // Top row
1514  // First column is same as neighbouring element
1515  finite_element_pt(element_num)->node_pt(n_p * (n_p - 1)) =
1516  finite_element_pt(element_num - 1)
1517  ->node_pt(n_p * (n_p - 1) + (n_p - 1));
1518  // New nodes for other columns
1519  for (unsigned l2 = 1; l2 < n_p; l2++)
1520  {
1521  local_node_num = n_p * (n_p - 1) + l2;
1522  // Allocate memory for node
1523  Node_pt[node_count] =
1524  finite_element_pt(element_num)
1525  ->construct_boundary_node(local_node_num, time_stepper_pt);
1526  // Set the pointer
1527  finite_element_pt(element_num)->node_pt(local_node_num) =
1528  Node_pt[node_count];
1529 
1530  // Get the fractional position of the node
1531  finite_element_pt(element_num)
1532  ->local_fraction_of_node(local_node_num, s_fraction);
1533 
1534  // Set the position of the node
1535  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
1536  Node_pt[node_count]->x(1) = Ymax;
1537  Node_pt[node_count]->x(2) = Zmin;
1538 
1539  // Add the node to the boundary
1540  add_boundary_node(3, Node_pt[node_count]);
1541  add_boundary_node(0, Node_pt[node_count]);
1542 
1543  // Increment the node number
1544  node_count++;
1545  }
1546 
1547 
1548  // Jump in the third dimension
1549 
1550  for (unsigned l3 = 1; l3 < n_p; l3++)
1551  {
1552  // The first row is copied from the lower element
1553  for (unsigned l2 = 0; l2 < n_p; l2++)
1554  {
1555  finite_element_pt(element_num)->node_pt(l2 + l3 * n_p * n_p) =
1556  finite_element_pt(element_num - Nx)
1557  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
1558  }
1559 
1560  // Second rows
1561  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
1562  {
1563  // First column is same as neighbouring element
1564  finite_element_pt(element_num)->node_pt(n_p * l1 + l3 * n_p * n_p) =
1565  finite_element_pt(element_num - 1)
1566  ->node_pt(n_p * l1 + (n_p - 1) + l3 * n_p * n_p);
1567 
1568  // New nodes for other columns
1569  for (unsigned l2 = 1; l2 < n_p; l2++)
1570  {
1571  local_node_num = n_p * l1 + l2 + l3 * n_p * n_p;
1572  // Allocate memory for the node
1573  Node_pt[node_count] =
1574  finite_element_pt(element_num)
1575  ->construct_node(local_node_num, time_stepper_pt);
1576 
1577  // Set the pointer
1578  finite_element_pt(element_num)->node_pt(local_node_num) =
1579  Node_pt[node_count];
1580 
1581  // Get the fractional position of the node
1582  finite_element_pt(element_num)
1583  ->local_fraction_of_node(local_node_num, s_fraction);
1584 
1585  // Set the position of the node
1586  Node_pt[node_count]->x(0) =
1587  Xmin + el_length[0] * (j + s_fraction[0]);
1588  Node_pt[node_count]->x(1) =
1589  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1590  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1591 
1592  // No boundaries
1593 
1594  // Increment the node number
1595  // add_boundary_node(0,Node_pt[node_count]);
1596  node_count++;
1597  }
1598  }
1599 
1600  // Top row
1601  // First column is same as neighbouring element
1602  finite_element_pt(element_num)
1603  ->node_pt(n_p * (n_p - 1) + l3 * n_p * n_p) =
1604  finite_element_pt(element_num - 1)
1605  ->node_pt(n_p * (n_p - 1) + (n_p - 1) + l3 * n_p * n_p);
1606  // New nodes for other columns
1607  for (unsigned l2 = 1; l2 < n_p; l2++)
1608  {
1609  local_node_num = n_p * (n_p - 1) + l2 + l3 * n_p * n_p;
1610  // Allocate memory for node
1611  Node_pt[node_count] =
1612  finite_element_pt(element_num)
1613  ->construct_boundary_node(local_node_num, time_stepper_pt);
1614  // Set the pointer
1615  finite_element_pt(element_num)->node_pt(local_node_num) =
1616  Node_pt[node_count];
1617 
1618  // Get the fractional position of the node
1619  finite_element_pt(element_num)
1620  ->local_fraction_of_node(local_node_num, s_fraction);
1621 
1622  // Set the position of the node
1623  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
1624  Node_pt[node_count]->x(1) = Ymax;
1625  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1626 
1627  // Add the node to the boundary
1628  add_boundary_node(3, Node_pt[node_count]);
1629 
1630  // Increment the node number
1631  node_count++;
1632  }
1633  }
1634 
1635  } // End of loop over central elements in row
1636 
1637 
1638  // FINAL ELEMENT IN ROW (upper right corner) IN LAYER z = 0
1639  //--------------------------------------------------------
1640 
1641  // Allocate memory for element
1642  element_num = Nx * (Ny - 1) + Nx - 1;
1643  Element_pt[element_num] = new ELEMENT;
1644 
1645  // We work first in the plane z =0
1646  // The first row is copied from the lower element
1647  for (unsigned l2 = 0; l2 < n_p; l2++)
1648  {
1649  finite_element_pt(element_num)->node_pt(l2) =
1650  finite_element_pt(element_num - Nx)->node_pt((n_p - 1) * n_p + l2);
1651  }
1652 
1653  // Second rows
1654  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
1655  {
1656  // First column is same as neighbouring element
1657  finite_element_pt(element_num)->node_pt(n_p * l1) =
1658  finite_element_pt(element_num - 1)->node_pt(n_p * l1 + (n_p - 1));
1659 
1660  // Middle nodes
1661  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
1662  {
1663  local_node_num = n_p * l1 + l2;
1664  // Allocate memory for node
1665  Node_pt[node_count] =
1666  finite_element_pt(element_num)
1667  ->construct_boundary_node(local_node_num, time_stepper_pt);
1668  // Set the pointer
1669  finite_element_pt(element_num)->node_pt(local_node_num) =
1670  Node_pt[node_count];
1671 
1672  // Get the fractional position of the node
1673  finite_element_pt(element_num)
1674  ->local_fraction_of_node(local_node_num, s_fraction);
1675 
1676  // Set the position of the node
1677  Node_pt[node_count]->x(0) =
1678  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
1679  Node_pt[node_count]->x(1) =
1680  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1681  Node_pt[node_count]->x(2) = Zmin;
1682 
1683  // Add the node to the boundary
1684  add_boundary_node(0, Node_pt[node_count]);
1685 
1686  // Increment the node number
1687  node_count++;
1688  }
1689 
1690  // Final node
1691  local_node_num = n_p * l1 + (n_p - 1);
1692  // Allocate memory for node
1693  Node_pt[node_count] =
1694  finite_element_pt(element_num)
1695  ->construct_boundary_node(local_node_num, time_stepper_pt);
1696  // Set the pointer
1697  finite_element_pt(element_num)->node_pt(local_node_num) =
1698  Node_pt[node_count];
1699 
1700  // Get the fractional position of the node
1701  finite_element_pt(element_num)
1702  ->local_fraction_of_node(local_node_num, s_fraction);
1703 
1704  // Set the position of the node
1705  Node_pt[node_count]->x(0) = Xmax;
1706  Node_pt[node_count]->x(1) =
1707  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1708  Node_pt[node_count]->x(2) = Zmin;
1709 
1710  // Add the node to the boundaries 0 and 2
1711  add_boundary_node(0, Node_pt[node_count]);
1712  add_boundary_node(2, Node_pt[node_count]);
1713 
1714 
1715  // Increment the node number
1716  node_count++;
1717 
1718  } // End of loop over middle rows
1719 
1720  // Final row
1721  // First column is same as neighbouring element
1722  finite_element_pt(element_num)->node_pt(n_p * (n_p - 1)) =
1723  finite_element_pt(element_num - 1)->node_pt(n_p * (n_p - 1) + (n_p - 1));
1724 
1725  // Middle nodes
1726  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
1727  {
1728  local_node_num = n_p * (n_p - 1) + l2;
1729  // Allocate memory for node
1730  Node_pt[node_count] =
1731  finite_element_pt(element_num)
1732  ->construct_boundary_node(local_node_num, time_stepper_pt);
1733  // Set the pointer
1734  finite_element_pt(element_num)->node_pt(local_node_num) =
1735  Node_pt[node_count];
1736 
1737  // Get the fractional position of the node
1738  finite_element_pt(element_num)
1739  ->local_fraction_of_node(local_node_num, s_fraction);
1740 
1741  // Set the position of the node
1742  Node_pt[node_count]->x(0) =
1743  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
1744  Node_pt[node_count]->x(1) = Ymax;
1745  Node_pt[node_count]->x(2) = Zmin;
1746 
1747  // Add the node to the boundaries 0 nd 3
1748  add_boundary_node(0, Node_pt[node_count]);
1749  add_boundary_node(3, Node_pt[node_count]);
1750 
1751 
1752  // Increment the node number
1753  node_count++;
1754  }
1755 
1756  // Final node
1757  // Determine number of values
1758  local_node_num = n_p * (n_p - 1) + (n_p - 1);
1759  // Allocate memory for node
1760  Node_pt[node_count] =
1761  finite_element_pt(element_num)
1762  ->construct_boundary_node(local_node_num, time_stepper_pt);
1763  // Set the pointer
1764  finite_element_pt(element_num)->node_pt(local_node_num) =
1765  Node_pt[node_count];
1766 
1767  // Get the fractional position of the node
1768  finite_element_pt(element_num)
1769  ->local_fraction_of_node(local_node_num, s_fraction);
1770 
1771  // Set the position of the node
1772  Node_pt[node_count]->x(0) = Xmax;
1773  Node_pt[node_count]->x(1) = Ymax;
1774  Node_pt[node_count]->x(2) = Zmin;
1775 
1776  // Add the node to the boundaries 0,2 and 3
1777  add_boundary_node(0, Node_pt[node_count]);
1778  add_boundary_node(2, Node_pt[node_count]);
1779  add_boundary_node(3, Node_pt[node_count]);
1780 
1781  // Increment the node number
1782  node_count++;
1783 
1784  // We jump to the third dimension
1785 
1786  for (unsigned l3 = 1; l3 < n_p; l3++)
1787  {
1788  // The first row is copied from the lower element
1789  for (unsigned l2 = 0; l2 < n_p; l2++)
1790  {
1791  finite_element_pt(element_num)->node_pt(l2 + l3 * n_p * n_p) =
1792  finite_element_pt(element_num - Nx)
1793  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
1794  }
1795 
1796  // Second rows
1797  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
1798  {
1799  // First column is same as neighbouring element
1800  finite_element_pt(element_num)->node_pt(n_p * l1 + l3 * n_p * n_p) =
1801  finite_element_pt(element_num - 1)
1802  ->node_pt(n_p * l1 + (n_p - 1) + l3 * n_p * n_p);
1803 
1804  // Middle nodes
1805  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
1806  {
1807  // Determine number of values
1808  local_node_num = n_p * l1 + l2 + l3 * n_p * n_p;
1809  // Allocate memory for node
1810  Node_pt[node_count] =
1811  finite_element_pt(element_num)
1812  ->construct_node(local_node_num, time_stepper_pt);
1813  // Set the pointer
1814  finite_element_pt(element_num)->node_pt(local_node_num) =
1815  Node_pt[node_count];
1816 
1817  // Get the fractional position of the node
1818  finite_element_pt(element_num)
1819  ->local_fraction_of_node(local_node_num, s_fraction);
1820 
1821  // Set the position of the node
1822  Node_pt[node_count]->x(0) =
1823  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
1824  Node_pt[node_count]->x(1) =
1825  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1826  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1827 
1828  // No boundaries
1829 
1830  // Increment the node number
1831  node_count++;
1832  }
1833 
1834  // Final node
1835  // Determine number of values
1836  local_node_num = n_p * l1 + (n_p - 1) + l3 * n_p * n_p;
1837  // Allocate memory for node
1838  Node_pt[node_count] =
1839  finite_element_pt(element_num)
1840  ->construct_boundary_node(local_node_num, time_stepper_pt);
1841  // Set the pointer
1842  finite_element_pt(element_num)->node_pt(local_node_num) =
1843  Node_pt[node_count];
1844 
1845  // Get the fractional position of the node
1846  finite_element_pt(element_num)
1847  ->local_fraction_of_node(local_node_num, s_fraction);
1848 
1849  // Set the position of the node
1850  Node_pt[node_count]->x(0) = Xmax;
1851  Node_pt[node_count]->x(1) =
1852  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
1853  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1854 
1855  // Add the node to the boundary 2
1856  add_boundary_node(2, Node_pt[node_count]);
1857  // Increment the node number
1858  node_count++;
1859 
1860  } // End of loop over middle rows
1861 
1862  // Final row
1863  // First column is same as neighbouring element
1864  finite_element_pt(element_num)
1865  ->node_pt(n_p * (n_p - 1) + l3 * n_p * n_p) =
1866  finite_element_pt(element_num - 1)
1867  ->node_pt(n_p * (n_p - 1) + (n_p - 1) + l3 * n_p * n_p);
1868 
1869  // Middle nodes
1870  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
1871  {
1872  // Determine number of values
1873  local_node_num = n_p * (n_p - 1) + l2 + l3 * n_p * n_p;
1874  // Allocate memory for node
1875  Node_pt[node_count] =
1876  finite_element_pt(element_num)
1877  ->construct_boundary_node(local_node_num, time_stepper_pt);
1878  // Set the pointer
1879  finite_element_pt(element_num)->node_pt(local_node_num) =
1880  Node_pt[node_count];
1881 
1882  // Get the fractional position of the node
1883  finite_element_pt(element_num)
1884  ->local_fraction_of_node(local_node_num, s_fraction);
1885 
1886  // Set the position of the node
1887  Node_pt[node_count]->x(0) =
1888  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
1889  Node_pt[node_count]->x(1) = Ymax;
1890  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1891 
1892  // Add the node to the boundary 3
1893  add_boundary_node(3, Node_pt[node_count]);
1894  // Increment the node number
1895  node_count++;
1896  }
1897 
1898  // Final node
1899  // Determine number of values
1900  local_node_num = n_p * (n_p - 1) + (n_p - 1) + l3 * n_p * n_p;
1901  // Allocate memory for node
1902  Node_pt[node_count] =
1903  finite_element_pt(element_num)
1904  ->construct_boundary_node(local_node_num, time_stepper_pt);
1905  // Set the pointer
1906  finite_element_pt(element_num)->node_pt(local_node_num) =
1907  Node_pt[node_count];
1908 
1909  // Get the fractional position of the node
1910  finite_element_pt(element_num)
1911  ->local_fraction_of_node(local_node_num, s_fraction);
1912 
1913  // Set the position of the node
1914  Node_pt[node_count]->x(0) = Xmax;
1915  Node_pt[node_count]->x(1) = Ymax;
1916  Node_pt[node_count]->x(2) = Zmin + el_length[2] * s_fraction[2];
1917 
1918  // Add the node to the boundaries 2 and 3
1919  add_boundary_node(2, Node_pt[node_count]);
1920  add_boundary_node(3, Node_pt[node_count]);
1921 
1922  // Increment the node number
1923  node_count++;
1924  }
1925 
1926  // END OF THE FIRST LAYER
1927 
1928 
1929  //----------------------------------------------------------------------------------------------------------------------------
1930  // ***************************************NOW WE MAKE ALL THE INTREMEDIATE
1931  // LAYERS**********************************************
1932  //----------------------------------------------------------------------------------------------------------------------------
1933 
1934 
1935  for (unsigned k = 1; k < (Nz - 1); k++) // good loop for the diferent layers
1936  // for(unsigned k=1;k<Nz;k++) // bad loop for testing the hole mesh but the
1937  // last layer
1938  {
1939  // FIRST ELEMENT OF THE LAYER
1940  //----------------------------------
1941 
1942  element_num = k * Nx * Ny;
1943  Element_pt[element_num] = new ELEMENT;
1944 
1945  // The lowest nodes are copied from the lower element
1946  for (unsigned l1 = 0; l1 < n_p; l1++)
1947  {
1948  for (unsigned l2 = 0; l2 < n_p; l2++)
1949  {
1950  finite_element_pt(element_num)->node_pt(l2 + n_p * l1) =
1951  finite_element_pt(element_num - Nx * Ny)
1952  ->node_pt(l2 + n_p * l1 + n_p * n_p * (n_p - 1));
1953  }
1954  }
1955 
1956 
1957  // Loop over the other node columns in the z direction
1958 
1959  for (unsigned l3 = 1; l3 < n_p; l3++)
1960  {
1961  // Set the corner node
1962  // Determine number of values at this node
1963  local_node_num = n_p * n_p * l3;
1964 
1965  // Allocate memory for the node
1966  Node_pt[node_count] =
1967  finite_element_pt(element_num)
1968  ->construct_boundary_node(local_node_num, time_stepper_pt);
1969 
1970  // Set the pointer from the element
1971  finite_element_pt(element_num)->node_pt(local_node_num) =
1972  Node_pt[node_count];
1973 
1974  // Get the fractional position of the node
1975  finite_element_pt(element_num)
1976  ->local_fraction_of_node(local_node_num, s_fraction);
1977 
1978  // Set the position of the node
1979  Node_pt[node_count]->x(0) = Xmin;
1980  Node_pt[node_count]->x(1) = Ymin;
1981  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
1982 
1983  // Add the node to boundaries 1 and 4
1984  add_boundary_node(1, Node_pt[node_count]);
1985  add_boundary_node(4, Node_pt[node_count]);
1986 
1987  // Increment the node number
1988  node_count++;
1989 
1990  // Loop over the other nodes in the first row in the x direction
1991  for (unsigned l2 = 1; l2 < n_p; l2++)
1992  {
1993  // Determine the number of values at this node
1994  local_node_num = l2 + n_p * n_p * l3;
1995 
1996  // Allocate memory for the nodes
1997  Node_pt[node_count] =
1998  finite_element_pt(element_num)
1999  ->construct_boundary_node(local_node_num, time_stepper_pt);
2000  // Set the pointer from the element
2001  finite_element_pt(element_num)->node_pt(local_node_num) =
2002  Node_pt[node_count];
2003 
2004  // Get the fractional position of the node
2005  finite_element_pt(element_num)
2006  ->local_fraction_of_node(local_node_num, s_fraction);
2007 
2008  // Set the position of the node
2009  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
2010  Node_pt[node_count]->x(1) = Ymin;
2011  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
2012 
2013  // Add the node to the boundary 1
2014  add_boundary_node(1, Node_pt[node_count]);
2015  // Increment the node number
2016  node_count++;
2017  }
2018 
2019  // Loop over the other node columns
2020  for (unsigned l1 = 1; l1 < n_p; l1++)
2021  {
2022  // Determine the number of values
2023  local_node_num = l1 * n_p + n_p * n_p * l3;
2024 
2025  // Allocate memory for the nodes
2026  Node_pt[node_count] =
2027  finite_element_pt(element_num)
2028  ->construct_boundary_node(local_node_num, time_stepper_pt);
2029  // Set the pointer from the element
2030  finite_element_pt(element_num)->node_pt(local_node_num) =
2031  Node_pt[node_count];
2032 
2033  // Get the fractional position of the node
2034  finite_element_pt(element_num)
2035  ->local_fraction_of_node(local_node_num, s_fraction);
2036 
2037  // Set the position of the first node of the row (with boundary 4)
2038  Node_pt[node_count]->x(0) = Xmin;
2039  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
2040  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
2041 
2042  // Add the node to the boundary 4
2043  add_boundary_node(4, Node_pt[node_count]);
2044  // Increment the node number
2045  node_count++;
2046 
2047  // Loop over the other nodes in the row
2048  for (unsigned l2 = 1; l2 < n_p; l2++)
2049  {
2050  // Set the number of values
2051  local_node_num = l1 * n_p + l2 + n_p * n_p * l3;
2052 
2053  // Allocate the memory for the node
2054  Node_pt[node_count] =
2055  finite_element_pt(element_num)
2056  ->construct_node(local_node_num, time_stepper_pt);
2057  // Set the pointer from the element
2058  finite_element_pt(element_num)->node_pt(local_node_num) =
2059  Node_pt[node_count];
2060 
2061  // Get the fractional position of the node
2062  finite_element_pt(element_num)
2063  ->local_fraction_of_node(local_node_num, s_fraction);
2064 
2065  // Set the position of the node
2066  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
2067  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
2068  Node_pt[node_count]->x(2) =
2069  Zmin + el_length[2] * (k + s_fraction[2]);
2070 
2071 
2072  // No boundary
2073 
2074  // Increment the node number
2075  node_count++;
2076  }
2077  }
2078  }
2079 
2080 
2081  //----------------------------------------------------------------------
2082 
2083  // CENTRE OF FIRST ROW OF ELEMENTS
2084  //--------------------------------
2085 
2086  // Now loop over the first row of elements, apart from final element
2087  for (unsigned j = 1; j < (Nx - 1); j++)
2088  {
2089  // Allocate memory for new element
2090  element_num = j + k * Nx * Ny;
2091  Element_pt[element_num] = new ELEMENT;
2092 
2093  // The lowest nodes are copied from the lower element
2094  for (unsigned l1 = 0; l1 < n_p; l1++)
2095  {
2096  for (unsigned l2 = 0; l2 < n_p; l2++)
2097  {
2098  finite_element_pt(j + k * Nx * Ny)->node_pt(l2 + n_p * l1) =
2099  finite_element_pt(j + (k - 1) * Nx * Ny)
2100  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
2101  }
2102  }
2103 
2104  // Loop over the other node columns in the z direction
2105  for (unsigned l3 = 1; l3 < n_p; l3++)
2106  {
2107  // First column of nodes is same as neighbouring element
2108  finite_element_pt(j + k * Nx * Ny)->node_pt(l3 * n_p * n_p) =
2109  finite_element_pt(j - 1 + k * Nx * Ny)
2110  ->node_pt(l3 * n_p * n_p + (n_p - 1));
2111 
2112  // New nodes for other columns
2113  for (unsigned l2 = 1; l2 < n_p; l2++)
2114  {
2115  // Determine number of values
2116  local_node_num = l2 + l3 * n_p * n_p;
2117 
2118  // Allocate memory for the nodes
2119  Node_pt[node_count] =
2120  finite_element_pt(element_num)
2121  ->construct_boundary_node(local_node_num, time_stepper_pt);
2122  // Set the pointer from the element
2123  finite_element_pt(element_num)->node_pt(local_node_num) =
2124  Node_pt[node_count];
2125 
2126  // Get the fractional position of the node
2127  finite_element_pt(element_num)
2128  ->local_fraction_of_node(local_node_num, s_fraction);
2129 
2130  // Set the position of the node
2131  Node_pt[node_count]->x(0) =
2132  Xmin + el_length[0] * (j + s_fraction[0]);
2133  Node_pt[node_count]->x(1) = Ymin;
2134  Node_pt[node_count]->x(2) =
2135  Zmin + el_length[2] * (k + s_fraction[2]);
2136 
2137  // Add the node to the boundary 1
2138  add_boundary_node(1, Node_pt[node_count]);
2139  // Increment the node number
2140  node_count++;
2141  }
2142 
2143  // Do the rest of the nodes
2144  for (unsigned l1 = 1; l1 < n_p; l1++)
2145  {
2146  // First column of nodes is same as neighbouring element
2147  finite_element_pt(j + k * Nx * Ny)
2148  ->node_pt(l1 * n_p + l3 * n_p * n_p) =
2149  finite_element_pt(j - 1 + k * Nx * Ny)
2150  ->node_pt(l1 * n_p + (n_p - 1) + l3 * n_p * n_p);
2151 
2152  // New nodes for other columns
2153  for (unsigned l2 = 1; l2 < n_p; l2++)
2154  {
2155  // Determine number of values
2156  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
2157 
2158  // Allocate memory for the nodes
2159  Node_pt[node_count] =
2160  finite_element_pt(element_num)
2161  ->construct_node(local_node_num, time_stepper_pt);
2162  // Set the pointer from the element
2163  finite_element_pt(element_num)->node_pt(local_node_num) =
2164  Node_pt[node_count];
2165 
2166  // Get the fractional position of the node
2167  finite_element_pt(element_num)
2168  ->local_fraction_of_node(local_node_num, s_fraction);
2169 
2170  // Set the position of the node
2171  Node_pt[node_count]->x(0) =
2172  Xmin + el_length[0] * (j + s_fraction[0]);
2173  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
2174  Node_pt[node_count]->x(2) =
2175  Zmin + el_length[2] * (k + s_fraction[2]);
2176 
2177  // No boundaries
2178 
2179  // Increment the node number
2180  node_count++;
2181  }
2182  }
2183  }
2184  }
2185 
2186  // MY FINAL ELEMENT IN FIRST ROW (right corner)
2187  //-----------------------------------------------
2188 
2189 
2190  // Allocate memory for new element
2191  element_num = Nx - 1 + k * Nx * Ny;
2192  Element_pt[element_num] = new ELEMENT;
2193 
2194  // The lowest nodes are copied from the lower element
2195  for (unsigned l1 = 0; l1 < n_p; l1++)
2196  {
2197  for (unsigned l2 = 0; l2 < n_p; l2++)
2198  {
2199  finite_element_pt(Nx - 1 + k * Nx * Ny)->node_pt(l2 + n_p * l1) =
2200  finite_element_pt(Nx - 1 + (k - 1) * Nx * Ny)
2201  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
2202  }
2203  }
2204 
2205 
2206  // Loop over the other node columns in the z direction
2207  for (unsigned l3 = 1; l3 < n_p; l3++)
2208  {
2209  // First column of nodes is same as neighbouring element
2210  finite_element_pt(Nx - 1 + k * Nx * Ny)->node_pt(l3 * n_p * n_p) =
2211  finite_element_pt(Nx - 2 + k * Nx * Ny)
2212  ->node_pt(l3 * n_p * n_p + (n_p - 1));
2213 
2214  // New nodes for other rows (y=0)
2215  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
2216  {
2217  // Determine number of values
2218  local_node_num = l2 + l3 * n_p * n_p;
2219 
2220  // Allocate memory for the nodes
2221  Node_pt[node_count] =
2222  finite_element_pt(element_num)
2223  ->construct_boundary_node(local_node_num, time_stepper_pt);
2224  // Set the pointer from the element
2225  finite_element_pt(element_num)->node_pt(local_node_num) =
2226  Node_pt[node_count];
2227 
2228  // Get the fractional position of the node
2229  finite_element_pt(element_num)
2230  ->local_fraction_of_node(local_node_num, s_fraction);
2231 
2232  // Set the position of the node
2233  Node_pt[node_count]->x(0) =
2234  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
2235  Node_pt[node_count]->x(1) = Ymin;
2236  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
2237 
2238  // Add the node to the boundary 1
2239  add_boundary_node(1, Node_pt[node_count]);
2240  // Increment the node number
2241  node_count++;
2242  }
2243 
2244  // Last node of the row (y=0)
2245 
2246  // Determine number of values
2247  local_node_num = (n_p - 1) + l3 * n_p * n_p;
2248 
2249  // Allocate memory for the nodes
2250  Node_pt[node_count] =
2251  finite_element_pt(element_num)
2252  ->construct_boundary_node(local_node_num, time_stepper_pt);
2253  // Set the pointer from the element
2254  finite_element_pt(element_num)->node_pt(local_node_num) =
2255  Node_pt[node_count];
2256 
2257  // Get the fractional position of the node
2258  finite_element_pt(element_num)
2259  ->local_fraction_of_node(local_node_num, s_fraction);
2260 
2261  // Set the position of the node
2262  Node_pt[node_count]->x(0) = Xmax;
2263  Node_pt[node_count]->x(1) = Ymin;
2264  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
2265 
2266  // Add the node to the boundary 1 and 2
2267  add_boundary_node(1, Node_pt[node_count]);
2268  add_boundary_node(2, Node_pt[node_count]);
2269  // Increment the node number
2270  node_count++;
2271 
2272  // Do the rest of the nodes
2273  for (unsigned l1 = 1; l1 < n_p; l1++)
2274  {
2275  // First column of nodes is same as neighbouring element
2276  finite_element_pt(Nx - 1 + k * Nx * Ny)
2277  ->node_pt(l1 * n_p + l3 * n_p * n_p) =
2278  finite_element_pt(Nx - 2 + k * Nx * Ny)
2279  ->node_pt(l1 * n_p + (n_p - 1) + l3 * n_p * n_p);
2280 
2281  // New nodes for other columns
2282  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
2283  {
2284  // Determine number of values
2285  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
2286 
2287  // Allocate memory for the nodes
2288  Node_pt[node_count] =
2289  finite_element_pt(element_num)
2290  ->construct_node(local_node_num, time_stepper_pt);
2291  // Set the pointer from the element
2292  finite_element_pt(element_num)->node_pt(local_node_num) =
2293  Node_pt[node_count];
2294 
2295  // Get the fractional position of the node
2296  finite_element_pt(element_num)
2297  ->local_fraction_of_node(local_node_num, s_fraction);
2298 
2299  // Set the position of the node
2300  Node_pt[node_count]->x(0) =
2301  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
2302  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
2303  Node_pt[node_count]->x(2) =
2304  Zmin + el_length[2] * (k + s_fraction[2]);
2305 
2306  // No boundaries
2307 
2308  // Increment the node number
2309  node_count++;
2310  }
2311 
2312  // Last nodes (at the surface x = Lx)
2313  // Determine number of values
2314  local_node_num = l1 * n_p + (n_p - 1) + l3 * n_p * n_p;
2315  // Allocate memory for the nodes
2316  Node_pt[node_count] =
2317  finite_element_pt(element_num)
2318  ->construct_boundary_node(local_node_num, time_stepper_pt);
2319  // Set the pointer from the element
2320  finite_element_pt(element_num)->node_pt(local_node_num) =
2321  Node_pt[node_count];
2322 
2323  // Get the fractional position of the node
2324  finite_element_pt(element_num)
2325  ->local_fraction_of_node(local_node_num, s_fraction);
2326 
2327 
2328  // Set the position of the node
2329  Node_pt[node_count]->x(0) = Xmax;
2330  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
2331  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
2332 
2333  // Add the node to the boundary 2
2334  add_boundary_node(2, Node_pt[node_count]);
2335 
2336  // Increment the node number
2337  node_count++;
2338  }
2339  }
2340 
2341 
2342  // ALL CENTRAL ELEMENT ROWS
2343  //------------------------
2344 
2345  // Loop over remaining element rows
2346  for (unsigned i = 1; i < (Ny - 1); i++)
2347  {
2348  // Set the first element in the row
2349 
2350  // Allocate memory for element (x =0)
2351  element_num = Nx * i + Nx * Ny * k;
2352  Element_pt[element_num] = new ELEMENT;
2353 
2354  // The lowest nodes are copied from the lower element
2355  for (unsigned l1 = 0; l1 < n_p; l1++)
2356  {
2357  for (unsigned l2 = 0; l2 < n_p; l2++)
2358  {
2359  finite_element_pt(Nx * i + k * Nx * Ny)->node_pt(l2 + n_p * l1) =
2360  finite_element_pt(Nx * i + (k - 1) * Nx * Ny)
2361  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
2362  }
2363  }
2364 
2365  // We extend now this element to the third direction
2366 
2367  for (unsigned l3 = 1; l3 < n_p; l3++)
2368  {
2369  // The first row of nodes is copied from the element below (at z=0)
2370  for (unsigned l2 = 0; l2 < n_p; l2++)
2371  {
2372  finite_element_pt(Nx * i + k * Nx * Ny)
2373  ->node_pt(l2 + l3 * n_p * n_p) =
2374  finite_element_pt(Nx * (i - 1) + k * Nx * Ny)
2375  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
2376  }
2377 
2378  // Other rows are new nodes (first the nodes for which x=0)
2379  for (unsigned l1 = 1; l1 < n_p; l1++)
2380  {
2381  // First column of nodes
2382 
2383  // Determine number of values
2384  local_node_num = l1 * n_p + l3 * n_p * n_p;
2385 
2386  // Allocate memory for node
2387  Node_pt[node_count] =
2388  finite_element_pt(element_num)
2389  ->construct_boundary_node(local_node_num, time_stepper_pt);
2390  // Set the pointer from the element
2391  finite_element_pt(element_num)->node_pt(local_node_num) =
2392  Node_pt[node_count];
2393 
2394  // Get the fractional position of the node
2395  finite_element_pt(element_num)
2396  ->local_fraction_of_node(local_node_num, s_fraction);
2397 
2398  // Set the position of the node
2399  Node_pt[node_count]->x(0) = Xmin;
2400  Node_pt[node_count]->x(1) =
2401  Ymin + el_length[1] * (i + s_fraction[1]);
2402  Node_pt[node_count]->x(2) =
2403  Zmin + el_length[2] * (k + s_fraction[2]);
2404 
2405  // Add the node to the boundary 4
2406  add_boundary_node(4, Node_pt[node_count]);
2407 
2408  // Increment the node number
2409  node_count++;
2410 
2411  // Now do the other columns (we extend to the rest of the nodes)
2412  for (unsigned l2 = 1; l2 < n_p; l2++)
2413  {
2414  // Determine number of values
2415  local_node_num = l1 * n_p + l2 + n_p * n_p * l3;
2416 
2417  // Allocate memory for node
2418  Node_pt[node_count] =
2419  finite_element_pt(element_num)
2420  ->construct_node(local_node_num, time_stepper_pt);
2421  // Set the pointer from the element
2422  finite_element_pt(element_num)->node_pt(local_node_num) =
2423  Node_pt[node_count];
2424 
2425  // Get the fractional position of the node
2426  finite_element_pt(element_num)
2427  ->local_fraction_of_node(local_node_num, s_fraction);
2428 
2429  // Set the position of the node
2430  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
2431  Node_pt[node_count]->x(1) =
2432  Ymin + el_length[1] * (i + s_fraction[1]);
2433  Node_pt[node_count]->x(2) =
2434  Zmin + el_length[2] * (k + s_fraction[2]);
2435 
2436  // No boundaries
2437 
2438  // Increment the node number
2439  node_count++;
2440  }
2441  }
2442  }
2443 
2444  // Now loop over the rest of the elements in the row, apart from the
2445  // last
2446  for (unsigned j = 1; j < (Nx - 1); j++)
2447  {
2448  // Allocate memory for new element
2449  element_num = Nx * i + j + k * Nx * Ny;
2450  Element_pt[element_num] = new ELEMENT;
2451 
2452  // The lowest nodes are copied from the lower element
2453  for (unsigned l1 = 0; l1 < n_p; l1++)
2454  {
2455  for (unsigned l2 = 0; l2 < n_p; l2++)
2456  {
2457  finite_element_pt(Nx * i + j + k * Nx * Ny)
2458  ->node_pt(l2 + n_p * l1) =
2459  finite_element_pt(Nx * i + j + (k - 1) * Nx * Ny)
2460  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
2461  }
2462  }
2463 
2464  // We extend to the third dimension
2465 
2466  for (unsigned l3 = 1; l3 < n_p; l3++)
2467  {
2468  // The first row is copied from the lower element
2469 
2470  for (unsigned l2 = 0; l2 < n_p; l2++)
2471  {
2472  finite_element_pt(Nx * i + j + k * Nx * Ny)
2473  ->node_pt(l2 + l3 * n_p * n_p) =
2474  finite_element_pt(Nx * (i - 1) + j + k * Nx * Ny)
2475  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
2476  }
2477 
2478  for (unsigned l1 = 1; l1 < n_p; l1++)
2479  {
2480  // First column is same as neighbouring element
2481  finite_element_pt(Nx * i + j + k * Nx * Ny)
2482  ->node_pt(l1 * n_p + l3 * n_p * n_p) =
2483  finite_element_pt(Nx * i + (j - 1) + k * Nx * Ny)
2484  ->node_pt(l1 * n_p + l3 * n_p * n_p + (n_p - 1));
2485 
2486  // New nodes for other columns
2487  for (unsigned l2 = 1; l2 < n_p; l2++)
2488  {
2489  // Determine number of values
2490  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
2491 
2492  // Allocate memory for the nodes
2493  Node_pt[node_count] =
2494  finite_element_pt(element_num)
2495  ->construct_node(local_node_num, time_stepper_pt);
2496  // Set the pointer
2497  finite_element_pt(element_num)->node_pt(local_node_num) =
2498  Node_pt[node_count];
2499  // Get the fractional position of the node
2500  finite_element_pt(element_num)
2501  ->local_fraction_of_node(local_node_num, s_fraction);
2502 
2503 
2504  // Set the position of the node
2505  Node_pt[node_count]->x(0) =
2506  Xmin + el_length[0] * (j + s_fraction[0]);
2507  Node_pt[node_count]->x(1) =
2508  Ymin + el_length[1] * (i + s_fraction[1]);
2509  Node_pt[node_count]->x(2) =
2510  Zmin + el_length[2] * (k + s_fraction[2]);
2511 
2512 
2513  // No boundaries
2514  // Increment the node number
2515  node_count++;
2516  }
2517  }
2518  }
2519 
2520  } // End of loop over elements in row
2521 
2522 
2523  // Do final element in row
2524 
2525  // Allocate memory for element
2526  element_num = Nx * i + Nx - 1 + k * Nx * Ny;
2527  Element_pt[element_num] = new ELEMENT;
2528 
2529  // The lowest nodes are copied from the lower element
2530  for (unsigned l1 = 0; l1 < n_p; l1++)
2531  {
2532  for (unsigned l2 = 0; l2 < n_p; l2++)
2533  {
2534  finite_element_pt(Nx * i + Nx - 1 + k * Nx * Ny)
2535  ->node_pt(l2 + n_p * l1) =
2536  finite_element_pt(Nx * i + Nx - 1 + (k - 1) * Nx * Ny)
2537  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
2538  }
2539  }
2540 
2541 
2542  // We go to the third dimension
2543  for (unsigned l3 = 1; l3 < n_p; l3++)
2544  {
2545  // The first row is copied from the lower element
2546  for (unsigned l2 = 0; l2 < n_p; l2++)
2547  {
2548  finite_element_pt(Nx * i + Nx - 1 + k * Nx * Ny)
2549  ->node_pt(l2 + l3 * n_p * n_p) =
2550  finite_element_pt(Nx * (i - 1) + Nx - 1 + k * Nx * Ny)
2551  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
2552  }
2553 
2554  for (unsigned l1 = 1; l1 < n_p; l1++)
2555  {
2556  // First column is same as neighbouring element
2557  finite_element_pt(Nx * i + Nx - 1 + k * Nx * Ny)
2558  ->node_pt(l1 * n_p + l3 * n_p * n_p) =
2559  finite_element_pt(Nx * i + Nx - 2 + k * Nx * Ny)
2560  ->node_pt(l1 * n_p + (n_p - 1) + l3 * n_p * n_p);
2561 
2562  // Middle nodes
2563  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
2564  {
2565  // Determine number of values
2566  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
2567 
2568  // Allocate memory for node
2569  Node_pt[node_count] =
2570  finite_element_pt(element_num)
2571  ->construct_node(local_node_num, time_stepper_pt);
2572  // Set the pointer
2573  finite_element_pt(element_num)->node_pt(local_node_num) =
2574  Node_pt[node_count];
2575 
2576  // Get the fractional position of the node
2577  finite_element_pt(element_num)
2578  ->local_fraction_of_node(local_node_num, s_fraction);
2579 
2580  // Set the position of the node
2581  Node_pt[node_count]->x(0) =
2582  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
2583  Node_pt[node_count]->x(1) =
2584  Ymin + el_length[1] * (i + s_fraction[1]);
2585  Node_pt[node_count]->x(2) =
2586  Zmin + el_length[2] * (k + s_fraction[2]);
2587 
2588  // No boundaries
2589 
2590  // Increment the node number
2591  node_count++;
2592  }
2593 
2594  // Final node
2595 
2596  // Determine number of values
2597  local_node_num = l1 * n_p + (n_p - 1) + l3 * n_p * n_p;
2598 
2599  // Allocate memory for node
2600  Node_pt[node_count] =
2601  finite_element_pt(element_num)
2602  ->construct_boundary_node(local_node_num, time_stepper_pt);
2603  // Set the pointer
2604  finite_element_pt(element_num)->node_pt(local_node_num) =
2605  Node_pt[node_count];
2606 
2607  // Get the fractional position of the node
2608  finite_element_pt(element_num)
2609  ->local_fraction_of_node(local_node_num, s_fraction);
2610 
2611  // Set the position of the node
2612  Node_pt[node_count]->x(0) = Xmax;
2613  Node_pt[node_count]->x(1) =
2614  Ymin + el_length[1] * (i + s_fraction[1]);
2615  Node_pt[node_count]->x(2) =
2616  Zmin + el_length[2] * (k + s_fraction[2]);
2617 
2618  // Add the node to the boundary 2
2619  add_boundary_node(2, Node_pt[node_count]);
2620 
2621  // Increment the node number
2622  node_count++;
2623  } // End of loop over rows of nodes in the element
2624 
2625 
2626  } // End of the 3dimension loop
2627 
2628 
2629  } // End of loop over rows of elements
2630 
2631 
2632  // FINAL ELEMENT ROW (IN INTERMIDIATE LAYERS)
2633  //=================
2634 
2635 
2636  // FIRST ELEMENT IN UPPER ROW (upper left corner)
2637  //----------------------------------------------
2638 
2639  // Allocate memory for element
2640  element_num = Nx * (Ny - 1) + k * Nx * Ny;
2641  Element_pt[element_num] = new ELEMENT;
2642 
2643  // The lowest nodes are copied from the lower element
2644  for (unsigned l1 = 0; l1 < n_p; l1++)
2645  {
2646  for (unsigned l2 = 0; l2 < n_p; l2++)
2647  {
2648  finite_element_pt(Nx * (Ny - 1) + k * Nx * Ny)
2649  ->node_pt(l2 + n_p * l1) =
2650  finite_element_pt(Nx * (Ny - 1) + (k - 1) * Nx * Ny)
2651  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
2652  }
2653  }
2654 
2655  // We jump to the third dimension
2656  for (unsigned l3 = 1; l3 < n_p; l3++)
2657  {
2658  // The first row of nodes is copied from the element below
2659  for (unsigned l2 = 0; l2 < n_p; l2++)
2660  {
2661  finite_element_pt(Nx * (Ny - 1) + k * Nx * Ny)
2662  ->node_pt(l2 + l3 * n_p * n_p) =
2663  finite_element_pt(Nx * (Ny - 2) + k * Nx * Ny)
2664  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
2665  }
2666 
2667  // Second row of nodes
2668  // First column of nodes
2669  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
2670  {
2671  // Determine number of values
2672  local_node_num = n_p * l1 + l3 * n_p * n_p;
2673 
2674  // Allocate memory for node
2675  Node_pt[node_count] =
2676  finite_element_pt(element_num)
2677  ->construct_boundary_node(local_node_num, time_stepper_pt);
2678  // Set the pointer from the element
2679  finite_element_pt(element_num)->node_pt(local_node_num) =
2680  Node_pt[node_count];
2681 
2682  // Get the fractional position of the node
2683  finite_element_pt(element_num)
2684  ->local_fraction_of_node(local_node_num, s_fraction);
2685 
2686  // Set the position of the node
2687  Node_pt[node_count]->x(0) = Xmin;
2688  Node_pt[node_count]->x(1) =
2689  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
2690  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
2691 
2692  // Add the node to the boundary 4
2693  add_boundary_node(4, Node_pt[node_count]);
2694 
2695  // Increment the node number
2696  node_count++;
2697 
2698  // Now do the other columns
2699  for (unsigned l2 = 1; l2 < n_p; l2++)
2700  {
2701  // Determine number of values
2702  local_node_num = n_p * l1 + l2 + l3 * n_p * n_p;
2703 
2704  // Allocate memory for node
2705  Node_pt[node_count] =
2706  finite_element_pt(element_num)
2707  ->construct_node(local_node_num, time_stepper_pt);
2708  // Set the pointer from the element
2709  finite_element_pt(element_num)->node_pt(local_node_num) =
2710  Node_pt[node_count];
2711 
2712  // Get the fractional position of the node
2713  finite_element_pt(element_num)
2714  ->local_fraction_of_node(local_node_num, s_fraction);
2715 
2716  // Set the position of the node
2717  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
2718  Node_pt[node_count]->x(1) =
2719  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
2720  Node_pt[node_count]->x(2) =
2721  Zmin + el_length[2] * (k + s_fraction[2]);
2722 
2723  // No boundaries
2724 
2725  // Increment the node number
2726  node_count++;
2727  }
2728  }
2729 
2730  // Final row of nodes
2731  // First column of nodes
2732  // Top left node
2733  // Determine number of values
2734  local_node_num = n_p * (n_p - 1) + l3 * n_p * n_p;
2735  // Allocate memory for node
2736  Node_pt[node_count] =
2737  finite_element_pt(element_num)
2738  ->construct_boundary_node(local_node_num, time_stepper_pt);
2739  // Set the pointer from the element
2740  finite_element_pt(element_num)->node_pt(local_node_num) =
2741  Node_pt[node_count];
2742 
2743  // Get the fractional position of the node
2744  finite_element_pt(element_num)
2745  ->local_fraction_of_node(local_node_num, s_fraction);
2746 
2747  // Set the position of the node
2748  Node_pt[node_count]->x(0) = Xmin;
2749  Node_pt[node_count]->x(1) = Ymax;
2750  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
2751 
2752  // Add the node to the boundaries
2753  add_boundary_node(3, Node_pt[node_count]);
2754  add_boundary_node(4, Node_pt[node_count]);
2755 
2756  // Increment the node number
2757  node_count++;
2758 
2759  // Now do the other columns
2760  for (unsigned l2 = 1; l2 < n_p; l2++)
2761  {
2762  // Determine number of values
2763  local_node_num = n_p * (n_p - 1) + l2 + l3 * n_p * n_p;
2764  // Allocate memory for the node
2765  Node_pt[node_count] =
2766  finite_element_pt(element_num)
2767  ->construct_boundary_node(local_node_num, time_stepper_pt);
2768  // Set the pointer from the element
2769  finite_element_pt(element_num)->node_pt(local_node_num) =
2770  Node_pt[node_count];
2771 
2772  // Get the fractional position of the node
2773  finite_element_pt(element_num)
2774  ->local_fraction_of_node(local_node_num, s_fraction);
2775 
2776  // Set the position of the node
2777  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
2778  Node_pt[node_count]->x(1) = Ymax;
2779  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
2780 
2781  // Add the node to the boundary 3
2782  add_boundary_node(3, Node_pt[node_count]);
2783 
2784  // Increment the node number
2785  node_count++;
2786  }
2787  }
2788 
2789 
2790  // Now loop over the rest of the elements in the row, apart from the last
2791  for (unsigned j = 1; j < (Nx - 1); j++)
2792  {
2793  // Allocate memory for element
2794  element_num = Nx * (Ny - 1) + j + k * Nx * Ny;
2795  Element_pt[element_num] = new ELEMENT;
2796 
2797  // The lowest nodes are copied from the lower element
2798  for (unsigned l1 = 0; l1 < n_p; l1++)
2799  {
2800  for (unsigned l2 = 0; l2 < n_p; l2++)
2801  {
2802  finite_element_pt(Nx * (Ny - 1) + j + k * Nx * Ny)
2803  ->node_pt(l2 + n_p * l1) =
2804  finite_element_pt(Nx * (Ny - 1) + j + (k - 1) * Nx * Ny)
2805  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
2806  }
2807  }
2808 
2809  // Jump in the third dimension
2810 
2811  for (unsigned l3 = 1; l3 < n_p; l3++)
2812  {
2813  // The first row is copied from the lower element
2814  for (unsigned l2 = 0; l2 < n_p; l2++)
2815  {
2816  finite_element_pt(Nx * (Ny - 1) + j + k * Nx * Ny)
2817  ->node_pt(l2 + l3 * n_p * n_p) =
2818  finite_element_pt(Nx * (Ny - 2) + j + k * Nx * Ny)
2819  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
2820  }
2821 
2822  // Second rows
2823  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
2824  {
2825  // First column is same as neighbouring element
2826  finite_element_pt(Nx * (Ny - 1) + j + k * Nx * Ny)
2827  ->node_pt(n_p * l1 + l3 * n_p * n_p) =
2828  finite_element_pt(Nx * (Ny - 1) + (j - 1) + k * Nx * Ny)
2829  ->node_pt(n_p * l1 + (n_p - 1) + l3 * n_p * n_p);
2830 
2831  // New nodes for other columns
2832  for (unsigned l2 = 1; l2 < n_p; l2++)
2833  {
2834  // Determine number of values
2835  local_node_num = n_p * l1 + l2 + l3 * n_p * n_p;
2836  // Allocate memory for the node
2837  Node_pt[node_count] =
2838  finite_element_pt(element_num)
2839  ->construct_node(local_node_num, time_stepper_pt);
2840 
2841  // Set the pointer
2842  finite_element_pt(element_num)->node_pt(local_node_num) =
2843  Node_pt[node_count];
2844 
2845  // Get the fractional position of the node
2846  finite_element_pt(element_num)
2847  ->local_fraction_of_node(local_node_num, s_fraction);
2848 
2849  // Set the position of the node
2850  Node_pt[node_count]->x(0) =
2851  Xmin + el_length[0] * (j + s_fraction[0]);
2852  Node_pt[node_count]->x(1) =
2853  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
2854  Node_pt[node_count]->x(2) =
2855  Zmin + el_length[2] * (k + s_fraction[2]);
2856 
2857  // No boundaries
2858 
2859  // Increment the node number
2860  // add_boundary_node(0,Node_pt[node_count]);
2861  node_count++;
2862  }
2863  }
2864 
2865  // Top row
2866  // First column is same as neighbouring element
2867  finite_element_pt(Nx * (Ny - 1) + j + k * Nx * Ny)
2868  ->node_pt(n_p * (n_p - 1) + l3 * n_p * n_p) =
2869  finite_element_pt(Nx * (Ny - 1) + (j - 1) + k * Nx * Ny)
2870  ->node_pt(n_p * (n_p - 1) + (n_p - 1) + l3 * n_p * n_p);
2871  // New nodes for other columns
2872  for (unsigned l2 = 1; l2 < n_p; l2++)
2873  {
2874  // Determine number of values
2875  local_node_num = n_p * (n_p - 1) + l2 + l3 * n_p * n_p;
2876  // Allocate memory for node
2877  Node_pt[node_count] =
2878  finite_element_pt(element_num)
2879  ->construct_boundary_node(local_node_num, time_stepper_pt);
2880  // Set the pointer
2881  finite_element_pt(element_num)->node_pt(local_node_num) =
2882  Node_pt[node_count];
2883 
2884  // Get the fractional position of the node
2885  finite_element_pt(element_num)
2886  ->local_fraction_of_node(local_node_num, s_fraction);
2887 
2888  // Set the position of the node
2889  Node_pt[node_count]->x(0) =
2890  Xmin + el_length[0] * (j + s_fraction[0]);
2891  Node_pt[node_count]->x(1) = Ymax;
2892  Node_pt[node_count]->x(2) =
2893  Zmin + el_length[2] * (k + s_fraction[2]);
2894 
2895  // Add the node to the boundary
2896  add_boundary_node(3, Node_pt[node_count]);
2897 
2898  // Increment the node number
2899  node_count++;
2900  }
2901  }
2902 
2903  } // End of loop over central elements in row
2904 
2905 
2906  // FINAL ELEMENT IN ROW (upper right corner)
2907  //-----------------------------------------
2908 
2909  // Allocate memory for element
2910  element_num = Nx * (Ny - 1) + Nx - 1 + k * Nx * Ny;
2911  Element_pt[element_num] = new ELEMENT;
2912 
2913  // The lowest nodes are copied from the lower element
2914  for (unsigned l1 = 0; l1 < n_p; l1++)
2915  {
2916  for (unsigned l2 = 0; l2 < n_p; l2++)
2917  {
2918  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + k * Nx * Ny)
2919  ->node_pt(l2 + n_p * l1) =
2920  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + (k - 1) * Nx * Ny)
2921  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
2922  }
2923  }
2924 
2925  // We jump to the third dimension
2926 
2927 
2928  for (unsigned l3 = 1; l3 < n_p; l3++)
2929  {
2930  // The first row is copied from the lower element
2931  for (unsigned l2 = 0; l2 < n_p; l2++)
2932  {
2933  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + k * Nx * Ny)
2934  ->node_pt(l2 + l3 * n_p * n_p) =
2935  finite_element_pt(Nx * (Ny - 2) + Nx - 1 + k * Nx * Ny)
2936  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
2937  }
2938 
2939  // Second rows
2940  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
2941  {
2942  // First column is same as neighbouring element
2943  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + k * Nx * Ny)
2944  ->node_pt(n_p * l1 + l3 * n_p * n_p) =
2945  finite_element_pt(Nx * (Ny - 1) + Nx - 2 + k * Nx * Ny)
2946  ->node_pt(n_p * l1 + (n_p - 1) + l3 * n_p * n_p);
2947 
2948  // Middle nodes
2949  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
2950  {
2951  // Determine number of values
2952  local_node_num = n_p * l1 + l2 + l3 * n_p * n_p;
2953  // Allocate memory for node
2954  Node_pt[node_count] =
2955  finite_element_pt(element_num)
2956  ->construct_node(local_node_num, time_stepper_pt);
2957  // Set the pointer
2958  finite_element_pt(element_num)->node_pt(local_node_num) =
2959  Node_pt[node_count];
2960  // Get the fractional position of the node
2961  finite_element_pt(element_num)
2962  ->local_fraction_of_node(local_node_num, s_fraction);
2963 
2964 
2965  // Set the position of the node
2966  Node_pt[node_count]->x(0) =
2967  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
2968  Node_pt[node_count]->x(1) =
2969  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
2970  Node_pt[node_count]->x(2) =
2971  Zmin + el_length[2] * (k + s_fraction[2]);
2972 
2973  // No boundaries
2974 
2975  // Increment the node number
2976  node_count++;
2977  }
2978 
2979  // Final node
2980  // Determine number of values
2981  local_node_num = n_p * l1 + (n_p - 1) + l3 * n_p * n_p;
2982  // Allocate memory for node
2983  Node_pt[node_count] =
2984  finite_element_pt(element_num)
2985  ->construct_boundary_node(local_node_num, time_stepper_pt);
2986  // Set the pointer
2987  finite_element_pt(element_num)->node_pt(local_node_num) =
2988  Node_pt[node_count];
2989 
2990  // Get the fractional position of the node
2991  finite_element_pt(element_num)
2992  ->local_fraction_of_node(local_node_num, s_fraction);
2993 
2994  // Set the position of the node
2995  Node_pt[node_count]->x(0) = Xmax;
2996  Node_pt[node_count]->x(1) =
2997  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
2998  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
2999 
3000  // Add the node to the boundary 2
3001  add_boundary_node(2, Node_pt[node_count]);
3002 
3003  // Increment the node number
3004  node_count++;
3005 
3006  } // End of loop over middle rows
3007 
3008  // Final row
3009  // First column is same as neighbouring element
3010  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + k * Nx * Ny)
3011  ->node_pt(n_p * (n_p - 1) + l3 * n_p * n_p) =
3012  finite_element_pt(Nx * (Ny - 1) + Nx - 2 + k * Nx * Ny)
3013  ->node_pt(n_p * (n_p - 1) + (n_p - 1) + l3 * n_p * n_p);
3014 
3015  // Middle nodes
3016  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
3017  {
3018  // Determine number of values
3019  local_node_num = n_p * (n_p - 1) + l2 + l3 * n_p * n_p;
3020  // Allocate memory for node
3021  Node_pt[node_count] =
3022  finite_element_pt(element_num)
3023  ->construct_boundary_node(local_node_num, time_stepper_pt);
3024  // Set the pointer
3025  finite_element_pt(element_num)->node_pt(local_node_num) =
3026  Node_pt[node_count];
3027 
3028  // Get the fractional position of the node
3029  finite_element_pt(element_num)
3030  ->local_fraction_of_node(local_node_num, s_fraction);
3031 
3032  // Set the position of the node
3033  Node_pt[node_count]->x(0) =
3034  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
3035  Node_pt[node_count]->x(1) = Ymax;
3036  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
3037 
3038  // Add the node to the boundary 3
3039  add_boundary_node(3, Node_pt[node_count]);
3040 
3041  // Increment the node number
3042  node_count++;
3043  }
3044 
3045  // Final node
3046  // Determine number of values
3047  local_node_num = n_p * (n_p - 1) + (n_p - 1) + l3 * n_p * n_p;
3048  // Allocate memory for node
3049  Node_pt[node_count] =
3050  finite_element_pt(element_num)
3051  ->construct_boundary_node(local_node_num, time_stepper_pt);
3052  // Set the pointer
3053  finite_element_pt(element_num)->node_pt(local_node_num) =
3054  Node_pt[node_count];
3055 
3056  // Get the fractional position of the node
3057  finite_element_pt(element_num)
3058  ->local_fraction_of_node(local_node_num, s_fraction);
3059 
3060  // Set the position of the node
3061  Node_pt[node_count]->x(0) = Xmax;
3062  Node_pt[node_count]->x(1) = Ymax;
3063  Node_pt[node_count]->x(2) = Zmin + el_length[2] * (k + s_fraction[2]);
3064 
3065  // Add the node to the boundaries 2 and 3
3066  add_boundary_node(2, Node_pt[node_count]);
3067  add_boundary_node(3, Node_pt[node_count]);
3068 
3069  // Increment the node number
3070  node_count++;
3071  }
3072 
3073  } // End loop of the elements layer
3074 
3075 
3076  // END OF THE INTERMEDIATE LAYERS
3077 
3078  // ----------------------------------------------------------------------------------
3079  // ****************BEGINNING OF THE LAST
3080  // LAYER**************************************
3081  // ----------------------------------------------------------------------------------
3082 
3083 
3084  // FIRST ELEMENT OF THE UPPER LAYER
3085  //----------------------------------
3086 
3087  element_num = (Nz - 1) * Nx * Ny;
3088  Element_pt[element_num] = new ELEMENT;
3089 
3090  // The lowest nodes are copied from the lower element
3091  for (unsigned l1 = 0; l1 < n_p; l1++)
3092  {
3093  for (unsigned l2 = 0; l2 < n_p; l2++)
3094  {
3095  finite_element_pt((Nz - 1) * Nx * Ny)->node_pt(l2 + n_p * l1) =
3096  finite_element_pt((Nz - 2) * Nx * Ny)
3097  ->node_pt(l2 + n_p * l1 + n_p * n_p * (n_p - 1));
3098  }
3099  }
3100 
3101 
3102  // Loop over the other node columns in the z direction but the last
3103 
3104  for (unsigned l3 = 1; l3 < (n_p - 1); l3++)
3105  {
3106  // Set the corner nodes
3107  // Determine number of values at this node
3108  local_node_num = n_p * n_p * l3;
3109 
3110  // Allocate memory for the node
3111  Node_pt[node_count] =
3112  finite_element_pt(element_num)
3113  ->construct_boundary_node(local_node_num, time_stepper_pt);
3114 
3115  // Set the pointer from the element
3116  finite_element_pt(element_num)->node_pt(local_node_num) =
3117  Node_pt[node_count];
3118  // Get the fractional position of the node
3119  finite_element_pt(element_num)
3120  ->local_fraction_of_node(local_node_num, s_fraction);
3121 
3122 
3123  // Set the position of the node
3124  Node_pt[node_count]->x(0) = Xmin;
3125  Node_pt[node_count]->x(1) = Ymin;
3126  Node_pt[node_count]->x(2) =
3127  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3128 
3129  // Add the node to boundaries 1 and 4
3130  add_boundary_node(1, Node_pt[node_count]);
3131  add_boundary_node(4, Node_pt[node_count]);
3132 
3133  // Increment the node number
3134  node_count++;
3135 
3136  // Loop over the other nodes in the first row in the x direction
3137  for (unsigned l2 = 1; l2 < n_p; l2++)
3138  {
3139  // Determine the number of values at this node
3140  local_node_num = l2 + n_p * n_p * l3;
3141 
3142  // Allocate memory for the nodes
3143  Node_pt[node_count] =
3144  finite_element_pt(element_num)
3145  ->construct_boundary_node(local_node_num, time_stepper_pt);
3146  // Set the pointer from the element
3147  finite_element_pt(element_num)->node_pt(local_node_num) =
3148  Node_pt[node_count];
3149  // Get the fractional position of the node
3150  finite_element_pt(element_num)
3151  ->local_fraction_of_node(local_node_num, s_fraction);
3152 
3153 
3154  // Set the position of the node
3155  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
3156  Node_pt[node_count]->x(1) = Ymin;
3157  Node_pt[node_count]->x(2) =
3158  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3159 
3160  // Add the node to the boundary 1
3161  add_boundary_node(1, Node_pt[node_count]);
3162  // Increment the node number
3163  node_count++;
3164  }
3165 
3166  // Loop over the other node columns
3167  for (unsigned l1 = 1; l1 < n_p; l1++)
3168  {
3169  // Determine the number of values
3170  local_node_num = l1 * n_p + n_p * n_p * l3;
3171 
3172  // Allocate memory for the nodes
3173  Node_pt[node_count] =
3174  finite_element_pt(element_num)
3175  ->construct_boundary_node(local_node_num, time_stepper_pt);
3176  // Set the pointer from the element
3177  finite_element_pt(element_num)->node_pt(local_node_num) =
3178  Node_pt[node_count];
3179 
3180  // Get the fractional position of the node
3181  finite_element_pt(element_num)
3182  ->local_fraction_of_node(local_node_num, s_fraction);
3183 
3184  // Set the position of the first node of the row (with boundary 4)
3185  Node_pt[node_count]->x(0) = Xmin;
3186  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3187  Node_pt[node_count]->x(2) =
3188  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3189 
3190  // Add the node to the boundary 4
3191  add_boundary_node(4, Node_pt[node_count]);
3192  // Increment the node number
3193  node_count++;
3194 
3195  // Loop over the other nodes in the row
3196  for (unsigned l2 = 1; l2 < n_p; l2++)
3197  {
3198  // Set the number of values
3199  local_node_num = l1 * n_p + l2 + n_p * n_p * l3;
3200 
3201  // Allocate the memory for the node
3202  Node_pt[node_count] =
3203  finite_element_pt(element_num)
3204  ->construct_node(local_node_num, time_stepper_pt);
3205  // Set the pointer from the element
3206  finite_element_pt(element_num)->node_pt(local_node_num) =
3207  Node_pt[node_count];
3208 
3209  // Get the fractional position of the node
3210  finite_element_pt(element_num)
3211  ->local_fraction_of_node(local_node_num, s_fraction);
3212 
3213  // Set the position of the node
3214  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
3215  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3216  Node_pt[node_count]->x(2) =
3217  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3218 
3219  // No boundary
3220 
3221  // Increment the node number
3222  node_count++;
3223  }
3224  }
3225  }
3226 
3227  // Make the top nodes
3228 
3229  // Set the corner nodes
3230  // Determine number of values at this node
3231  local_node_num = n_p * n_p * (n_p - 1);
3232 
3233  // Allocate memory for the node
3234  Node_pt[node_count] =
3235  finite_element_pt(element_num)
3236  ->construct_boundary_node(local_node_num, time_stepper_pt);
3237 
3238  // Set the pointer from the element
3239  finite_element_pt(element_num)->node_pt(local_node_num) =
3240  Node_pt[node_count];
3241 
3242  // Get the fractional position of the node
3243  finite_element_pt(element_num)
3244  ->local_fraction_of_node(local_node_num, s_fraction);
3245 
3246  // Set the position of the node
3247  Node_pt[node_count]->x(0) = Xmin;
3248  Node_pt[node_count]->x(1) = Ymin;
3249  Node_pt[node_count]->x(2) = Zmax;
3250 
3251  // Add the node to boundaries 1, 4 and 5
3252  add_boundary_node(1, Node_pt[node_count]);
3253  add_boundary_node(4, Node_pt[node_count]);
3254  add_boundary_node(5, Node_pt[node_count]);
3255 
3256  // Increment the node number
3257  node_count++;
3258 
3259  // Loop over the other nodes in the first row in the x direction
3260  for (unsigned l2 = 1; l2 < n_p; l2++)
3261  {
3262  // Determine the number of values at this node
3263  local_node_num = l2 + n_p * n_p * (n_p - 1);
3264 
3265  // Allocate memory for the nodes
3266  Node_pt[node_count] =
3267  finite_element_pt(element_num)
3268  ->construct_boundary_node(local_node_num, time_stepper_pt);
3269  // Set the pointer from the element
3270  finite_element_pt(element_num)->node_pt(local_node_num) =
3271  Node_pt[node_count];
3272 
3273  // Get the fractional position of the node
3274  finite_element_pt(element_num)
3275  ->local_fraction_of_node(local_node_num, s_fraction);
3276 
3277  // Set the position of the node
3278  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
3279  Node_pt[node_count]->x(1) = Ymin;
3280  Node_pt[node_count]->x(2) = Zmax;
3281 
3282  // Add the node to the boundaries 1 and 5
3283  add_boundary_node(1, Node_pt[node_count]);
3284  add_boundary_node(5, Node_pt[node_count]);
3285  // Increment the node number
3286  node_count++;
3287  }
3288 
3289  // Loop over the other node columns
3290  for (unsigned l1 = 1; l1 < n_p; l1++)
3291  {
3292  // Determine the number of values
3293  local_node_num = l1 * n_p + n_p * n_p * (n_p - 1);
3294 
3295  // Allocate memory for the nodes
3296  Node_pt[node_count] =
3297  finite_element_pt(element_num)
3298  ->construct_boundary_node(local_node_num, time_stepper_pt);
3299  // Set the pointer from the element
3300  finite_element_pt(element_num)->node_pt(local_node_num) =
3301  Node_pt[node_count];
3302 
3303  // Get the fractional position of the node
3304  finite_element_pt(element_num)
3305  ->local_fraction_of_node(local_node_num, s_fraction);
3306 
3307  // Set the position of the first node of the row (with boundary 4)
3308  Node_pt[node_count]->x(0) = Xmin;
3309  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3310  Node_pt[node_count]->x(2) = Zmax;
3311 
3312  // Add the node to the boundary 4
3313  add_boundary_node(4, Node_pt[node_count]);
3314  add_boundary_node(5, Node_pt[node_count]);
3315  // Increment the node number
3316  node_count++;
3317 
3318  // Loop over the other nodes in the row
3319  for (unsigned l2 = 1; l2 < n_p; l2++)
3320  {
3321  // Set the number of values
3322  local_node_num = l1 * n_p + l2 + n_p * n_p * (n_p - 1);
3323 
3324  // Allocate the memory for the node
3325  Node_pt[node_count] =
3326  finite_element_pt(element_num)
3327  ->construct_boundary_node(local_node_num, time_stepper_pt);
3328  // Set the pointer from the element
3329  finite_element_pt(element_num)->node_pt(local_node_num) =
3330  Node_pt[node_count];
3331 
3332  // Get the fractional position of the node
3333  finite_element_pt(element_num)
3334  ->local_fraction_of_node(local_node_num, s_fraction);
3335 
3336  // Set the position of the node
3337  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
3338  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3339  Node_pt[node_count]->x(2) = Zmax;
3340 
3341  // Add the node to the boundary 5
3342  add_boundary_node(5, Node_pt[node_count]);
3343 
3344  // Increment the node number
3345  node_count++;
3346  }
3347  }
3348 
3349  //----------------------------------------------------------------------
3350 
3351  // CENTRE OF FIRST ROW OF ELEMENTS OF THE UPPER LAYER
3352  //--------------------------------------------------
3353 
3354  // Now loop over the first row of elements, apart from final element
3355  for (unsigned j = 1; j < (Nx - 1); j++)
3356  {
3357  // Allocate memory for new element
3358  element_num = j + (Nz - 1) * Nx * Ny;
3359  Element_pt[element_num] = new ELEMENT;
3360 
3361  // The lowest nodes are copied from the lower element
3362  for (unsigned l1 = 0; l1 < n_p; l1++)
3363  {
3364  for (unsigned l2 = 0; l2 < n_p; l2++)
3365  {
3366  finite_element_pt(j + (Nz - 1) * Nx * Ny)->node_pt(l2 + n_p * l1) =
3367  finite_element_pt(j + (Nz - 2) * Nx * Ny)
3368  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
3369  }
3370  }
3371 
3372  // Loop over the other node columns in the z direction but the last
3373  for (unsigned l3 = 1; l3 < (n_p - 1); l3++)
3374  {
3375  // First column of nodes is same as neighbouring element
3376  finite_element_pt(j + (Nz - 1) * Nx * Ny)->node_pt(l3 * n_p * n_p) =
3377  finite_element_pt(j - 1 + (Nz - 1) * Nx * Ny)
3378  ->node_pt(l3 * n_p * n_p + (n_p - 1));
3379 
3380  // New nodes for other columns
3381  for (unsigned l2 = 1; l2 < n_p; l2++)
3382  {
3383  // Determine number of values
3384  local_node_num = l2 + l3 * n_p * n_p;
3385 
3386  // Allocate memory for the nodes
3387  Node_pt[node_count] =
3388  finite_element_pt(element_num)
3389  ->construct_boundary_node(local_node_num, time_stepper_pt);
3390  // Set the pointer from the element
3391  finite_element_pt(element_num)->node_pt(local_node_num) =
3392  Node_pt[node_count];
3393 
3394  // Get the fractional position of the node
3395  finite_element_pt(element_num)
3396  ->local_fraction_of_node(local_node_num, s_fraction);
3397 
3398  // Set the position of the node
3399  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
3400  Node_pt[node_count]->x(1) = Ymin;
3401  Node_pt[node_count]->x(2) =
3402  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3403 
3404  // Add the node to the boundary 1
3405  add_boundary_node(1, Node_pt[node_count]);
3406  // Increment the node number
3407  node_count++;
3408  }
3409 
3410  // Do the rest of the nodes
3411  for (unsigned l1 = 1; l1 < n_p; l1++)
3412  {
3413  // First column of nodes is same as neighbouring element
3414  finite_element_pt(j + (Nz - 1) * Nx * Ny)
3415  ->node_pt(l1 * n_p + l3 * n_p * n_p) =
3416  finite_element_pt(j - 1 + (Nz - 1) * Nx * Ny)
3417  ->node_pt(l1 * n_p + (n_p - 1) + l3 * n_p * n_p);
3418 
3419  // New nodes for other columns
3420  for (unsigned l2 = 1; l2 < n_p; l2++)
3421  {
3422  // Determine number of values
3423  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
3424 
3425  // Allocate memory for the nodes
3426  Node_pt[node_count] =
3427  finite_element_pt(element_num)
3428  ->construct_node(local_node_num, time_stepper_pt);
3429  // Set the pointer from the element
3430  finite_element_pt(element_num)->node_pt(local_node_num) =
3431  Node_pt[node_count];
3432 
3433  // Get the fractional position of the node
3434  finite_element_pt(element_num)
3435  ->local_fraction_of_node(local_node_num, s_fraction);
3436 
3437  // Set the position of the node
3438  Node_pt[node_count]->x(0) =
3439  Xmin + el_length[0] * (j + s_fraction[0]);
3440  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3441  Node_pt[node_count]->x(2) =
3442  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3443 
3444  // No boundaries
3445 
3446  // Increment the node number
3447  node_count++;
3448  }
3449  }
3450  }
3451 
3452 
3453  // Top nodes
3454 
3455  // First node is same as neighbouring element
3456  finite_element_pt(j + (Nz - 1) * Nx * Ny)
3457  ->node_pt((n_p - 1) * n_p * n_p) =
3458  finite_element_pt(j - 1 + (Nz - 1) * Nx * Ny)
3459  ->node_pt((n_p - 1) * n_p * n_p + (n_p - 1));
3460 
3461  // New nodes for other columns
3462  for (unsigned l2 = 1; l2 < n_p; l2++)
3463  {
3464  // Determine number of values
3465  local_node_num = l2 + (n_p - 1) * n_p * n_p;
3466 
3467  // Allocate memory for the nodes
3468  Node_pt[node_count] =
3469  finite_element_pt(element_num)
3470  ->construct_boundary_node(local_node_num, time_stepper_pt);
3471  // Set the pointer from the element
3472  finite_element_pt(element_num)->node_pt(local_node_num) =
3473  Node_pt[node_count];
3474 
3475  // Get the fractional position of the node
3476  finite_element_pt(element_num)
3477  ->local_fraction_of_node(local_node_num, s_fraction);
3478 
3479  // Set the position of the node
3480  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
3481  Node_pt[node_count]->x(1) = Ymin;
3482  Node_pt[node_count]->x(2) = Zmax;
3483 
3484  // Add the node to the boundaries 1 and 5
3485  add_boundary_node(1, Node_pt[node_count]);
3486  add_boundary_node(5, Node_pt[node_count]);
3487  // Increment the node number
3488  node_count++;
3489  }
3490 
3491  // Do the rest of the nodes
3492  for (unsigned l1 = 1; l1 < n_p; l1++)
3493  {
3494  // First column of nodes is same as neighbouring element
3495  finite_element_pt(j + (Nz - 1) * Nx * Ny)
3496  ->node_pt(l1 * n_p + (n_p - 1) * n_p * n_p) =
3497  finite_element_pt(j - 1 + (Nz - 1) * Nx * Ny)
3498  ->node_pt(l1 * n_p + (n_p - 1) + (n_p - 1) * n_p * n_p);
3499 
3500  // New nodes for other columns
3501  for (unsigned l2 = 1; l2 < n_p; l2++)
3502  {
3503  // Determine number of values
3504  local_node_num = l1 * n_p + l2 + (n_p - 1) * n_p * n_p;
3505 
3506  // Allocate memory for the nodes
3507  Node_pt[node_count] =
3508  finite_element_pt(element_num)
3509  ->construct_boundary_node(local_node_num, time_stepper_pt);
3510  // Set the pointer from the element
3511  finite_element_pt(element_num)->node_pt(local_node_num) =
3512  Node_pt[node_count];
3513 
3514  // Get the fractional position of the node
3515  finite_element_pt(element_num)
3516  ->local_fraction_of_node(local_node_num, s_fraction);
3517 
3518  // Set the position of the node
3519  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
3520  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3521  Node_pt[node_count]->x(2) = Zmax;
3522 
3523  // Add to the boundary 5
3524  add_boundary_node(5, Node_pt[node_count]);
3525 
3526  // Increment the node number
3527  node_count++;
3528  }
3529  }
3530  } // End loop of first row of elements
3531 
3532 
3533  // MY FINAL ELEMENT IN THE FIRST ROW OF THE UPPER LAYER (right corner)
3534  //--------------------------------------------------------------------
3535 
3536 
3537  // Allocate memory for new element
3538  element_num = Nx - 1 + (Nz - 1) * Nx * Ny;
3539  Element_pt[element_num] = new ELEMENT;
3540 
3541  // The lowest nodes are copied from the lower element
3542  for (unsigned l1 = 0; l1 < n_p; l1++)
3543  {
3544  for (unsigned l2 = 0; l2 < n_p; l2++)
3545  {
3546  finite_element_pt(Nx - 1 + (Nz - 1) * Nx * Ny)->node_pt(l2 + n_p * l1) =
3547  finite_element_pt(Nx - 1 + (Nz - 2) * Nx * Ny)
3548  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
3549  }
3550  }
3551 
3552 
3553  // Loop over the other node columns in the z direction but the last
3554  for (unsigned l3 = 1; l3 < (n_p - 1); l3++)
3555  {
3556  // First column of nodes is same as neighbouring element
3557  finite_element_pt(Nx - 1 + (Nz - 1) * Nx * Ny)->node_pt(l3 * n_p * n_p) =
3558  finite_element_pt(Nx - 2 + (Nz - 1) * Nx * Ny)
3559  ->node_pt(l3 * n_p * n_p + (n_p - 1));
3560 
3561  // New nodes for other rows (y=0)
3562  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
3563  {
3564  // Determine number of values
3565  local_node_num = l2 + l3 * n_p * n_p;
3566 
3567  // Allocate memory for the nodes
3568  Node_pt[node_count] =
3569  finite_element_pt(element_num)
3570  ->construct_boundary_node(local_node_num, time_stepper_pt);
3571  // Set the pointer from the element
3572  finite_element_pt(element_num)->node_pt(local_node_num) =
3573  Node_pt[node_count];
3574 
3575  // Get the fractional position of the node
3576  finite_element_pt(element_num)
3577  ->local_fraction_of_node(local_node_num, s_fraction);
3578 
3579  // Set the position of the node
3580  Node_pt[node_count]->x(0) =
3581  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
3582  Node_pt[node_count]->x(1) = Ymin;
3583  Node_pt[node_count]->x(2) =
3584  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3585 
3586  // Add the node to the boundary 1
3587  add_boundary_node(1, Node_pt[node_count]);
3588  // Increment the node number
3589  node_count++;
3590  }
3591 
3592  // Last node of the row (y=0)
3593 
3594  // Determine number of values
3595  local_node_num = (n_p - 1) + l3 * n_p * n_p;
3596 
3597  // Allocate memory for the nodes
3598  Node_pt[node_count] =
3599  finite_element_pt(element_num)
3600  ->construct_boundary_node(local_node_num, time_stepper_pt);
3601  // Set the pointer from the element
3602  finite_element_pt(element_num)->node_pt(local_node_num) =
3603  Node_pt[node_count];
3604 
3605  // Get the fractional position of the node
3606  finite_element_pt(element_num)
3607  ->local_fraction_of_node(local_node_num, s_fraction);
3608 
3609  // Set the position of the node
3610  Node_pt[node_count]->x(0) = Xmax;
3611  Node_pt[node_count]->x(1) = Ymin;
3612  Node_pt[node_count]->x(2) =
3613  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3614 
3615  // Add the node to the boundary 1 and 2
3616  add_boundary_node(1, Node_pt[node_count]);
3617  add_boundary_node(2, Node_pt[node_count]);
3618  // Increment the node number
3619  node_count++;
3620 
3621  // Do the rest of the nodes
3622  for (unsigned l1 = 1; l1 < n_p; l1++)
3623  {
3624  // First column of nodes is same as neighbouring element
3625  finite_element_pt(Nx - 1 + (Nz - 1) * Nx * Ny)
3626  ->node_pt(l1 * n_p + l3 * n_p * n_p) =
3627  finite_element_pt(Nx - 2 + (Nz - 1) * Nx * Ny)
3628  ->node_pt(l1 * n_p + (n_p - 1) + l3 * n_p * n_p);
3629 
3630  // New nodes for other columns
3631  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
3632  {
3633  // Determine number of values
3634  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
3635 
3636  // Allocate memory for the nodes
3637  Node_pt[node_count] =
3638  finite_element_pt(element_num)
3639  ->construct_node(local_node_num, time_stepper_pt);
3640  // Set the pointer from the element
3641  finite_element_pt(element_num)->node_pt(local_node_num) =
3642  Node_pt[node_count];
3643 
3644  // Get the fractional position of the node
3645  finite_element_pt(element_num)
3646  ->local_fraction_of_node(local_node_num, s_fraction);
3647 
3648  // Set the position of the node
3649  Node_pt[node_count]->x(0) =
3650  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
3651  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3652  Node_pt[node_count]->x(2) =
3653  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3654 
3655  // No boundaries
3656 
3657  // Increment the node number
3658  node_count++;
3659  }
3660 
3661  // Last nodes (at the surface x = Lx)
3662  // Determine number of values
3663  local_node_num = l1 * n_p + (n_p - 1) + l3 * n_p * n_p;
3664  // Allocate memory for the nodes
3665  Node_pt[node_count] =
3666  finite_element_pt(element_num)
3667  ->construct_boundary_node(local_node_num, time_stepper_pt);
3668  // Set the pointer from the element
3669  finite_element_pt(element_num)->node_pt(local_node_num) =
3670  Node_pt[node_count];
3671 
3672  // Set the position of the node
3673  Node_pt[node_count]->x(0) = Xmax;
3674  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3675  Node_pt[node_count]->x(2) =
3676  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3677 
3678  // Add the node to the boundary 2
3679  add_boundary_node(2, Node_pt[node_count]);
3680 
3681  // Increment the node number
3682  node_count++;
3683  }
3684  }
3685 
3686  // We make the top nodes
3687  // First node is same as neighbouring element
3688  finite_element_pt(Nx - 1 + (Nz - 1) * Nx * Ny)
3689  ->node_pt((n_p - 1) * n_p * n_p) =
3690  finite_element_pt(Nx - 2 + (Nz - 1) * Nx * Ny)
3691  ->node_pt((n_p - 1) * n_p * n_p + (n_p - 1));
3692 
3693  // New nodes for other rows (y=0)
3694  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
3695  {
3696  // Determine number of values
3697  local_node_num = l2 + (n_p - 1) * n_p * n_p;
3698 
3699  // Allocate memory for the nodes
3700  Node_pt[node_count] =
3701  finite_element_pt(element_num)
3702  ->construct_boundary_node(local_node_num, time_stepper_pt);
3703  // Set the pointer from the element
3704  finite_element_pt(element_num)->node_pt(local_node_num) =
3705  Node_pt[node_count];
3706 
3707  // Get the fractional position of the node
3708  finite_element_pt(element_num)
3709  ->local_fraction_of_node(local_node_num, s_fraction);
3710 
3711  // Set the position of the node
3712  Node_pt[node_count]->x(0) =
3713  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
3714  Node_pt[node_count]->x(1) = Ymin;
3715  Node_pt[node_count]->x(2) = Zmax;
3716 
3717  // Add the node to the boundaries 1 and 5
3718  add_boundary_node(1, Node_pt[node_count]);
3719  add_boundary_node(5, Node_pt[node_count]);
3720 
3721  // Increment the node number
3722  node_count++;
3723  }
3724 
3725  // Last node of the row (y=0)
3726 
3727  // Determine number of values
3728  local_node_num = (n_p - 1) + (n_p - 1) * n_p * n_p;
3729 
3730  // Allocate memory for the nodes
3731  Node_pt[node_count] =
3732  finite_element_pt(element_num)
3733  ->construct_boundary_node(local_node_num, time_stepper_pt);
3734  // Set the pointer from the element
3735  finite_element_pt(element_num)->node_pt(local_node_num) =
3736  Node_pt[node_count];
3737 
3738  // Get the fractional position of the node
3739  finite_element_pt(element_num)
3740  ->local_fraction_of_node(local_node_num, s_fraction);
3741 
3742  // Set the position of the node
3743  Node_pt[node_count]->x(0) = Xmax;
3744  Node_pt[node_count]->x(1) = Ymin;
3745  Node_pt[node_count]->x(2) = Zmax;
3746 
3747  // Add the node to the boundary 1 and 2
3748  add_boundary_node(1, Node_pt[node_count]);
3749  add_boundary_node(2, Node_pt[node_count]);
3750  add_boundary_node(5, Node_pt[node_count]);
3751  // Increment the node number
3752  node_count++;
3753 
3754  // Do the rest of the nodes
3755  for (unsigned l1 = 1; l1 < n_p; l1++)
3756  {
3757  // First column of nodes is same as neighbouring element
3758  finite_element_pt(Nx - 1 + (Nz - 1) * Nx * Ny)
3759  ->node_pt(l1 * n_p + (n_p - 1) * n_p * n_p) =
3760  finite_element_pt(Nx - 2 + (Nz - 1) * Nx * Ny)
3761  ->node_pt(l1 * n_p + (n_p - 1) + (n_p - 1) * n_p * n_p);
3762 
3763  // New nodes for other columns
3764  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
3765  {
3766  // Determine number of values
3767  local_node_num = l1 * n_p + l2 + (n_p - 1) * n_p * n_p;
3768 
3769  // Allocate memory for the nodes
3770  Node_pt[node_count] =
3771  finite_element_pt(element_num)
3772  ->construct_boundary_node(local_node_num, time_stepper_pt);
3773  // Set the pointer from the element
3774  finite_element_pt(element_num)->node_pt(local_node_num) =
3775  Node_pt[node_count];
3776 
3777  // Get the fractional position of the node
3778  finite_element_pt(element_num)
3779  ->local_fraction_of_node(local_node_num, s_fraction);
3780 
3781  // Set the position of the node
3782  Node_pt[node_count]->x(0) =
3783  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
3784  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3785  Node_pt[node_count]->x(2) = Zmax;
3786 
3787  // Add node to the boundary 5
3788  add_boundary_node(5, Node_pt[node_count]);
3789 
3790  // Increment the node number
3791  node_count++;
3792  }
3793 
3794  // Last nodes (at the surface x = Lx)
3795  // Determine number of values
3796  local_node_num = l1 * n_p + (n_p - 1) + (n_p - 1) * n_p * n_p;
3797  // Allocate memory for the nodes
3798  Node_pt[node_count] =
3799  finite_element_pt(element_num)
3800  ->construct_boundary_node(local_node_num, time_stepper_pt);
3801  // Set the pointer from the element
3802  finite_element_pt(element_num)->node_pt(local_node_num) =
3803  Node_pt[node_count];
3804 
3805  // Get the fractional position of the node
3806  finite_element_pt(element_num)
3807  ->local_fraction_of_node(local_node_num, s_fraction);
3808 
3809  // Set the position of the node
3810  Node_pt[node_count]->x(0) = Xmax;
3811  Node_pt[node_count]->x(1) = Ymin + el_length[1] * s_fraction[1];
3812  Node_pt[node_count]->x(2) = Zmax;
3813 
3814  // Add the node to the boundaries 2 and 5
3815  add_boundary_node(2, Node_pt[node_count]);
3816  add_boundary_node(5, Node_pt[node_count]);
3817 
3818 
3819  // Increment the node number
3820  node_count++;
3821  }
3822 
3823 
3824  // ALL CENTRAL ELEMENT ROWS OF THE TOP LAYER
3825  //------------------------------------------
3826 
3827  // Loop over remaining element rows
3828  for (unsigned i = 1; i < (Ny - 1); i++)
3829  {
3830  // Set the first element in the row
3831 
3832  // Allocate memory for element (x =0)
3833  element_num = Nx * i + Nx * Ny * (Nz - 1);
3834  Element_pt[element_num] = new ELEMENT;
3835 
3836  // The lowest nodes are copied from the lower element
3837  for (unsigned l1 = 0; l1 < n_p; l1++)
3838  {
3839  for (unsigned l2 = 0; l2 < n_p; l2++)
3840  {
3841  finite_element_pt(Nx * i + (Nz - 1) * Nx * Ny)
3842  ->node_pt(l2 + n_p * l1) =
3843  finite_element_pt(Nx * i + (Nz - 2) * Nx * Ny)
3844  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
3845  }
3846  }
3847 
3848  // We extend now this element to the third dimension
3849 
3850  for (unsigned l3 = 1; l3 < (n_p - 1); l3++)
3851  {
3852  // The first row of nodes is copied from the element below (at z=0)
3853  for (unsigned l2 = 0; l2 < n_p; l2++)
3854  {
3855  finite_element_pt(Nx * i + (Nz - 1) * Nx * Ny)
3856  ->node_pt(l2 + l3 * n_p * n_p) =
3857  finite_element_pt(Nx * (i - 1) + (Nz - 1) * Nx * Ny)
3858  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
3859  }
3860 
3861  // Other rows are new nodes (first the nodes for which x=0)
3862  for (unsigned l1 = 1; l1 < n_p; l1++)
3863  {
3864  // First column of nodes
3865 
3866  // Determine number of values
3867  local_node_num = l1 * n_p + l3 * n_p * n_p;
3868 
3869  // Allocate memory for node
3870  Node_pt[node_count] =
3871  finite_element_pt(element_num)
3872  ->construct_boundary_node(local_node_num, time_stepper_pt);
3873  // Set the pointer from the element
3874  finite_element_pt(element_num)->node_pt(local_node_num) =
3875  Node_pt[node_count];
3876 
3877  // Get the fractional position of the node
3878  finite_element_pt(element_num)
3879  ->local_fraction_of_node(local_node_num, s_fraction);
3880 
3881  // Set the position of the node
3882  Node_pt[node_count]->x(0) = Xmin;
3883  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
3884  Node_pt[node_count]->x(2) =
3885  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3886 
3887  // Add the node to the boundary 4
3888  add_boundary_node(4, Node_pt[node_count]);
3889 
3890  // Increment the node number
3891  node_count++;
3892 
3893  // Now do the other columns (we extend to the rest of the nodes)
3894  for (unsigned l2 = 1; l2 < n_p; l2++)
3895  {
3896  // Determine number of values
3897  local_node_num = l1 * n_p + l2 + n_p * n_p * l3;
3898 
3899  // Allocate memory for node
3900  Node_pt[node_count] =
3901  finite_element_pt(element_num)
3902  ->construct_node(local_node_num, time_stepper_pt);
3903  // Set the pointer from the element
3904  finite_element_pt(element_num)->node_pt(local_node_num) =
3905  Node_pt[node_count];
3906 
3907  // Get the fractional position of the node
3908  finite_element_pt(element_num)
3909  ->local_fraction_of_node(local_node_num, s_fraction);
3910 
3911  // Set the position of the node
3912  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
3913  Node_pt[node_count]->x(1) =
3914  Ymin + el_length[1] * (i + s_fraction[1]);
3915  Node_pt[node_count]->x(2) =
3916  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
3917 
3918  // No boundaries
3919 
3920  // Increment the node number
3921  node_count++;
3922  }
3923  }
3924  }
3925 
3926  // Now the top nodes
3927 
3928  // The first row of nodes is copied from the element below (at z=0)
3929  for (unsigned l2 = 0; l2 < n_p; l2++)
3930  {
3931  finite_element_pt(Nx * i + (Nz - 1) * Nx * Ny)
3932  ->node_pt(l2 + (n_p - 1) * n_p * n_p) =
3933  finite_element_pt(Nx * (i - 1) + (Nz - 1) * Nx * Ny)
3934  ->node_pt((n_p - 1) * n_p + l2 + (n_p - 1) * n_p * n_p);
3935  }
3936 
3937  // Other rows are new nodes (first the nodes for which x=0)
3938  for (unsigned l1 = 1; l1 < n_p; l1++)
3939  {
3940  // First column of nodes
3941 
3942  // Determine number of values
3943  local_node_num = l1 * n_p + (n_p - 1) * n_p * n_p;
3944 
3945  // Allocate memory for node
3946  Node_pt[node_count] =
3947  finite_element_pt(element_num)
3948  ->construct_boundary_node(local_node_num, time_stepper_pt);
3949  // Set the pointer from the element
3950  finite_element_pt(element_num)->node_pt(local_node_num) =
3951  Node_pt[node_count];
3952 
3953  // Get the fractional position of the node
3954  finite_element_pt(element_num)
3955  ->local_fraction_of_node(local_node_num, s_fraction);
3956 
3957  // Set the position of the node
3958  Node_pt[node_count]->x(0) = Xmin;
3959  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
3960  Node_pt[node_count]->x(2) = Zmax;
3961 
3962  // Add the node to the boundaries 4 and 5
3963  add_boundary_node(4, Node_pt[node_count]);
3964  add_boundary_node(5, Node_pt[node_count]);
3965 
3966  // Increment the node number
3967  node_count++;
3968 
3969  // Now do the other columns (we extend to the rest of the nodes)
3970  for (unsigned l2 = 1; l2 < n_p; l2++)
3971  {
3972  // Determine number of values
3973  local_node_num = l1 * n_p + l2 + n_p * n_p * (n_p - 1);
3974 
3975  // Allocate memory for node
3976  Node_pt[node_count] =
3977  finite_element_pt(element_num)
3978  ->construct_boundary_node(local_node_num, time_stepper_pt);
3979  // Set the pointer from the element
3980  finite_element_pt(element_num)->node_pt(local_node_num) =
3981  Node_pt[node_count];
3982 
3983  // Get the fractional position of the node
3984  finite_element_pt(element_num)
3985  ->local_fraction_of_node(local_node_num, s_fraction);
3986 
3987  // Set the position of the node
3988  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
3989  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
3990  Node_pt[node_count]->x(2) = Zmax;
3991 
3992  // Add the node to boundary 5
3993  add_boundary_node(5, Node_pt[node_count]);
3994 
3995  // Increment the node number
3996  node_count++;
3997  }
3998  }
3999 
4000 
4001  // Now loop over the rest of the elements in the row, apart from the last
4002  for (unsigned j = 1; j < (Nx - 1); j++)
4003  {
4004  // Allocate memory for new element
4005  element_num = Nx * i + j + (Nz - 1) * Nx * Ny;
4006  Element_pt[element_num] = new ELEMENT;
4007 
4008  // The lowest nodes are copied from the lower element
4009  for (unsigned l1 = 0; l1 < n_p; l1++)
4010  {
4011  for (unsigned l2 = 0; l2 < n_p; l2++)
4012  {
4013  finite_element_pt(Nx * i + j + (Nz - 1) * Nx * Ny)
4014  ->node_pt(l2 + n_p * l1) =
4015  finite_element_pt(Nx * i + j + (Nz - 2) * Nx * Ny)
4016  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
4017  }
4018  }
4019 
4020  // We extend to the third dimension but the last layer of nodes
4021 
4022  for (unsigned l3 = 1; l3 < (n_p - 1); l3++)
4023  {
4024  // The first row is copied from the lower element
4025 
4026  for (unsigned l2 = 0; l2 < n_p; l2++)
4027  {
4028  finite_element_pt(Nx * i + j + (Nz - 1) * Nx * Ny)
4029  ->node_pt(l2 + l3 * n_p * n_p) =
4030  finite_element_pt(Nx * (i - 1) + j + (Nz - 1) * Nx * Ny)
4031  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
4032  }
4033 
4034  for (unsigned l1 = 1; l1 < n_p; l1++)
4035  {
4036  // First column is same as neighbouring element
4037  finite_element_pt(Nx * i + j + (Nz - 1) * Nx * Ny)
4038  ->node_pt(l1 * n_p + l3 * n_p * n_p) =
4039  finite_element_pt(Nx * i + (j - 1) + (Nz - 1) * Nx * Ny)
4040  ->node_pt(l1 * n_p + l3 * n_p * n_p + (n_p - 1));
4041 
4042  // New nodes for other columns
4043  for (unsigned l2 = 1; l2 < n_p; l2++)
4044  {
4045  // Determine number of values
4046  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
4047 
4048  // Allocate memory for the nodes
4049  Node_pt[node_count] =
4050  finite_element_pt(element_num)
4051  ->construct_node(local_node_num, time_stepper_pt);
4052  // Set the pointer
4053  finite_element_pt(element_num)->node_pt(local_node_num) =
4054  Node_pt[node_count];
4055  // Get the fractional position of the node
4056  finite_element_pt(element_num)
4057  ->local_fraction_of_node(local_node_num, s_fraction);
4058 
4059 
4060  // Set the position of the node
4061  Node_pt[node_count]->x(0) =
4062  Xmin + el_length[0] * (j + s_fraction[0]);
4063  Node_pt[node_count]->x(1) =
4064  Ymin + el_length[1] * (i + s_fraction[1]);
4065  Node_pt[node_count]->x(2) =
4066  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4067  // No boundaries
4068 
4069  // Increment the node number
4070  node_count++;
4071  }
4072  }
4073  }
4074 
4075  // Now the top nodes
4076 
4077  // The first row is copied from the lower element
4078 
4079  for (unsigned l2 = 0; l2 < n_p; l2++)
4080  {
4081  finite_element_pt(Nx * i + j + (Nz - 1) * Nx * Ny)
4082  ->node_pt(l2 + (n_p - 1) * n_p * n_p) =
4083  finite_element_pt(Nx * (i - 1) + j + (Nz - 1) * Nx * Ny)
4084  ->node_pt((n_p - 1) * n_p + l2 + (n_p - 1) * n_p * n_p);
4085  }
4086 
4087  for (unsigned l1 = 1; l1 < n_p; l1++)
4088  {
4089  // First column is same as neighbouring element
4090  finite_element_pt(Nx * i + j + (Nz - 1) * Nx * Ny)
4091  ->node_pt(l1 * n_p + (n_p - 1) * n_p * n_p) =
4092  finite_element_pt(Nx * i + (j - 1) + (Nz - 1) * Nx * Ny)
4093  ->node_pt(l1 * n_p + (n_p - 1) * n_p * n_p + (n_p - 1));
4094 
4095  // New nodes for other columns
4096  for (unsigned l2 = 1; l2 < n_p; l2++)
4097  {
4098  // Determine number of values
4099  local_node_num = l1 * n_p + l2 + (n_p - 1) * n_p * n_p;
4100 
4101  // Allocate memory for the nodes
4102  Node_pt[node_count] =
4103  finite_element_pt(element_num)
4104  ->construct_boundary_node(local_node_num, time_stepper_pt);
4105  // Set the pointer
4106  finite_element_pt(element_num)->node_pt(local_node_num) =
4107  Node_pt[node_count];
4108 
4109  // Get the fractional position of the node
4110  finite_element_pt(element_num)
4111  ->local_fraction_of_node(local_node_num, s_fraction);
4112 
4113  // Set the position of the node
4114  Node_pt[node_count]->x(0) =
4115  Xmin + el_length[0] * (j + s_fraction[0]);
4116  Node_pt[node_count]->x(1) =
4117  Ymin + el_length[1] * (i + s_fraction[1]);
4118  Node_pt[node_count]->x(2) = Zmax;
4119 
4120  // Add to boundary 5
4121  add_boundary_node(5, Node_pt[node_count]);
4122 
4123  // Increment the node number
4124  node_count++;
4125  }
4126  }
4127 
4128 
4129  } // End of loop over elements in row
4130 
4131 
4132  // Do final element in row
4133 
4134  // Allocate memory for element
4135  element_num = Nx * i + Nx - 1 + (Nz - 1) * Nx * Ny;
4136  Element_pt[element_num] = new ELEMENT;
4137 
4138  // The lowest nodes are copied from the lower element
4139  for (unsigned l1 = 0; l1 < n_p; l1++)
4140  {
4141  for (unsigned l2 = 0; l2 < n_p; l2++)
4142  {
4143  finite_element_pt(Nx * i + Nx - 1 + (Nz - 1) * Nx * Ny)
4144  ->node_pt(l2 + n_p * l1) =
4145  finite_element_pt(Nx * i + Nx - 1 + (Nz - 2) * Nx * Ny)
4146  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
4147  }
4148  }
4149 
4150 
4151  // We go to the third dimension but we do not reach the top
4152  for (unsigned l3 = 1; l3 < (n_p - 1); l3++)
4153  {
4154  // The first row is copied from the lower element
4155  for (unsigned l2 = 0; l2 < n_p; l2++)
4156  {
4157  finite_element_pt(Nx * i + Nx - 1 + (Nz - 1) * Nx * Ny)
4158  ->node_pt(l2 + l3 * n_p * n_p) =
4159  finite_element_pt(Nx * (i - 1) + Nx - 1 + (Nz - 1) * Nx * Ny)
4160  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
4161  }
4162 
4163  for (unsigned l1 = 1; l1 < n_p; l1++)
4164  {
4165  // First column is same as neighbouring element
4166  finite_element_pt(Nx * i + Nx - 1 + (Nz - 1) * Nx * Ny)
4167  ->node_pt(l1 * n_p + l3 * n_p * n_p) =
4168  finite_element_pt(Nx * i + Nx - 2 + (Nz - 1) * Nx * Ny)
4169  ->node_pt(l1 * n_p + (n_p - 1) + l3 * n_p * n_p);
4170 
4171  // Middle nodes
4172  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
4173  {
4174  // Determine number of values
4175  local_node_num = l1 * n_p + l2 + l3 * n_p * n_p;
4176 
4177  // Allocate memory for node
4178  Node_pt[node_count] =
4179  finite_element_pt(element_num)
4180  ->construct_node(local_node_num, time_stepper_pt);
4181  // Set the pointer
4182  finite_element_pt(element_num)->node_pt(local_node_num) =
4183  Node_pt[node_count];
4184 
4185  // Get the fractional position of the node
4186  finite_element_pt(element_num)
4187  ->local_fraction_of_node(local_node_num, s_fraction);
4188 
4189  // Set the position of the node
4190  Node_pt[node_count]->x(0) =
4191  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
4192  Node_pt[node_count]->x(1) =
4193  Ymin + el_length[1] * (i + s_fraction[1]);
4194  Node_pt[node_count]->x(2) =
4195  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4196 
4197  // No boundaries
4198 
4199  // Increment the node number
4200  node_count++;
4201  }
4202 
4203  // Final node
4204 
4205  // Determine number of values
4206  local_node_num = l1 * n_p + (n_p - 1) + l3 * n_p * n_p;
4207 
4208  // Allocate memory for node
4209  Node_pt[node_count] =
4210  finite_element_pt(element_num)
4211  ->construct_boundary_node(local_node_num, time_stepper_pt);
4212  // Set the pointer
4213  finite_element_pt(element_num)->node_pt(local_node_num) =
4214  Node_pt[node_count];
4215 
4216  // Get the fractional position of the node
4217  finite_element_pt(element_num)
4218  ->local_fraction_of_node(local_node_num, s_fraction);
4219 
4220  // Set the position of the node
4221  Node_pt[node_count]->x(0) = Xmax;
4222  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
4223  Node_pt[node_count]->x(2) =
4224  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4225 
4226  // Add the node to the boundary 2
4227  add_boundary_node(2, Node_pt[node_count]);
4228 
4229  // Increment the node number
4230  node_count++;
4231  } // End of loop over rows of nodes in the element
4232 
4233 
4234  } // End of the 3dimension loop
4235 
4236  // Now the top nodes
4237 
4238  // The first row is copied from the lower element
4239  for (unsigned l2 = 0; l2 < n_p; l2++)
4240  {
4241  finite_element_pt(Nx * i + Nx - 1 + (Nz - 1) * Nx * Ny)
4242  ->node_pt(l2 + (n_p - 1) * n_p * n_p) =
4243  finite_element_pt(Nx * (i - 1) + Nx - 1 + (Nz - 1) * Nx * Ny)
4244  ->node_pt((n_p - 1) * n_p + l2 + (n_p - 1) * n_p * n_p);
4245  }
4246 
4247  for (unsigned l1 = 1; l1 < n_p; l1++)
4248  {
4249  // First column is same as neighbouring element
4250  finite_element_pt(Nx * i + Nx - 1 + (Nz - 1) * Nx * Ny)
4251  ->node_pt(l1 * n_p + (n_p - 1) * n_p * n_p) =
4252  finite_element_pt(Nx * i + Nx - 2 + (Nz - 1) * Nx * Ny)
4253  ->node_pt(l1 * n_p + (n_p - 1) + (n_p - 1) * n_p * n_p);
4254 
4255  // Middle nodes
4256  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
4257  {
4258  // Determine number of values
4259  local_node_num = l1 * n_p + l2 + (n_p - 1) * n_p * n_p;
4260 
4261  // Allocate memory for node
4262  Node_pt[node_count] =
4263  finite_element_pt(element_num)
4264  ->construct_boundary_node(local_node_num, time_stepper_pt);
4265  // Set the pointer
4266  finite_element_pt(element_num)->node_pt(local_node_num) =
4267  Node_pt[node_count];
4268 
4269  // Get the fractional position of the node
4270  finite_element_pt(element_num)
4271  ->local_fraction_of_node(local_node_num, s_fraction);
4272 
4273  // Set the position of the node
4274  Node_pt[node_count]->x(0) =
4275  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
4276  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
4277  Node_pt[node_count]->x(2) = Zmax;
4278 
4279  // Add to boundary 5
4280  add_boundary_node(5, Node_pt[node_count]);
4281 
4282  // Increment the node number
4283  node_count++;
4284  }
4285 
4286  // Final node
4287 
4288  // Determine number of values
4289  local_node_num = l1 * n_p + (n_p - 1) + (n_p - 1) * n_p * n_p;
4290 
4291  // Allocate memory for node
4292  Node_pt[node_count] =
4293  finite_element_pt(element_num)
4294  ->construct_boundary_node(local_node_num, time_stepper_pt);
4295  // Set the pointer
4296  finite_element_pt(element_num)->node_pt(local_node_num) =
4297  Node_pt[node_count];
4298 
4299  // Get the fractional position of the node
4300  finite_element_pt(element_num)
4301  ->local_fraction_of_node(local_node_num, s_fraction);
4302 
4303  // Set the position of the node
4304  Node_pt[node_count]->x(0) = Xmax;
4305  Node_pt[node_count]->x(1) = Ymin + el_length[1] * (i + s_fraction[1]);
4306  Node_pt[node_count]->x(2) = Zmax;
4307 
4308  // Add the node to the boundaries 2 and 5
4309  add_boundary_node(2, Node_pt[node_count]);
4310  add_boundary_node(5, Node_pt[node_count]);
4311 
4312  // Increment the node number
4313  node_count++;
4314  } // End of loop over rows of nodes in the element
4315 
4316 
4317  } // End of loop over rows of elements
4318 
4319 
4320  // FINAL ELEMENT ROW (IN TOP LAYER)
4321  //===========================================
4322 
4323 
4324  // FIRST ELEMENT IN UPPER ROW (upper left corner)
4325  //----------------------------------------------
4326 
4327  // Allocate memory for element
4328  element_num = Nx * (Ny - 1) + (Nz - 1) * Nx * Ny;
4329  Element_pt[element_num] = new ELEMENT;
4330 
4331  // The lowest nodes are copied from the lower element
4332  for (unsigned l1 = 0; l1 < n_p; l1++)
4333  {
4334  for (unsigned l2 = 0; l2 < n_p; l2++)
4335  {
4336  finite_element_pt(Nx * (Ny - 1) + (Nz - 1) * Nx * Ny)
4337  ->node_pt(l2 + n_p * l1) =
4338  finite_element_pt(Nx * (Ny - 1) + (Nz - 2) * Nx * Ny)
4339  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
4340  }
4341  }
4342 
4343  // We jump to the third dimension but the last layer of nodes
4344  for (unsigned l3 = 1; l3 < (n_p - 1); l3++)
4345  {
4346  // The first row of nodes is copied from the element below
4347  for (unsigned l2 = 0; l2 < n_p; l2++)
4348  {
4349  finite_element_pt(Nx * (Ny - 1) + (Nz - 1) * Nx * Ny)
4350  ->node_pt(l2 + l3 * n_p * n_p) =
4351  finite_element_pt(Nx * (Ny - 2) + (Nz - 1) * Nx * Ny)
4352  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
4353  }
4354 
4355  // Second row of nodes
4356  // First column of nodes
4357  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
4358  {
4359  // Determine number of values
4360  local_node_num = n_p * l1 + l3 * n_p * n_p;
4361 
4362  // Allocate memory for node
4363  Node_pt[node_count] =
4364  finite_element_pt(element_num)
4365  ->construct_boundary_node(local_node_num, time_stepper_pt);
4366  // Set the pointer from the element
4367  finite_element_pt(element_num)->node_pt(local_node_num) =
4368  Node_pt[node_count];
4369 
4370  // Get the fractional position of the node
4371  finite_element_pt(element_num)
4372  ->local_fraction_of_node(local_node_num, s_fraction);
4373 
4374  // Set the position of the node
4375  Node_pt[node_count]->x(0) = Xmin;
4376  Node_pt[node_count]->x(1) =
4377  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
4378  Node_pt[node_count]->x(2) =
4379  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4380 
4381  // Add the node to the boundary 4
4382  add_boundary_node(4, Node_pt[node_count]);
4383 
4384  // Increment the node number
4385  node_count++;
4386 
4387  // Now do the other columns
4388  for (unsigned l2 = 1; l2 < n_p; l2++)
4389  {
4390  // Determine number of values
4391  local_node_num = n_p * l1 + l2 + l3 * n_p * n_p;
4392 
4393  // Allocate memory for node
4394  Node_pt[node_count] =
4395  finite_element_pt(element_num)
4396  ->construct_node(local_node_num, time_stepper_pt);
4397  // Set the pointer from the element
4398  finite_element_pt(element_num)->node_pt(local_node_num) =
4399  Node_pt[node_count];
4400 
4401  // Get the fractional position of the node
4402  finite_element_pt(element_num)
4403  ->local_fraction_of_node(local_node_num, s_fraction);
4404 
4405  // Set the position of the node
4406  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
4407  Node_pt[node_count]->x(1) =
4408  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
4409  Node_pt[node_count]->x(2) =
4410  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4411 
4412  // No boundaries
4413 
4414  // Increment the node number
4415  node_count++;
4416  }
4417  }
4418 
4419  // Final row of nodes
4420  // First column of nodes
4421  // Top left node
4422  // Determine number of values
4423  local_node_num = n_p * (n_p - 1) + l3 * n_p * n_p;
4424  // Allocate memory for node
4425  Node_pt[node_count] =
4426  finite_element_pt(element_num)
4427  ->construct_boundary_node(local_node_num, time_stepper_pt);
4428  // Set the pointer from the element
4429  finite_element_pt(element_num)->node_pt(local_node_num) =
4430  Node_pt[node_count];
4431 
4432  // Get the fractional position of the node
4433  finite_element_pt(element_num)
4434  ->local_fraction_of_node(local_node_num, s_fraction);
4435 
4436  // Set the position of the node
4437  Node_pt[node_count]->x(0) = Xmin;
4438  Node_pt[node_count]->x(1) = Ymax;
4439  Node_pt[node_count]->x(2) =
4440  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4441 
4442  // Add the node to the boundaries
4443  add_boundary_node(3, Node_pt[node_count]);
4444  add_boundary_node(4, Node_pt[node_count]);
4445 
4446  // Increment the node number
4447  node_count++;
4448 
4449  // Now do the other columns
4450  for (unsigned l2 = 1; l2 < n_p; l2++)
4451  {
4452  // Determine number of values
4453  local_node_num = n_p * (n_p - 1) + l2 + l3 * n_p * n_p;
4454  // Allocate memory for the node
4455  Node_pt[node_count] =
4456  finite_element_pt(element_num)
4457  ->construct_boundary_node(local_node_num, time_stepper_pt);
4458  // Set the pointer from the element
4459  finite_element_pt(element_num)->node_pt(local_node_num) =
4460  Node_pt[node_count];
4461 
4462  // Get the fractional position of the node
4463  finite_element_pt(element_num)
4464  ->local_fraction_of_node(local_node_num, s_fraction);
4465 
4466  // Set the position of the node
4467  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
4468  Node_pt[node_count]->x(1) = Ymax;
4469  Node_pt[node_count]->x(2) =
4470  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4471 
4472  // Add the node to the boundary 3
4473  add_boundary_node(3, Node_pt[node_count]);
4474 
4475  // Increment the node number
4476  node_count++;
4477  }
4478  }
4479  // Now the top nodes
4480 
4481  // The first row of nodes is copied from the element below
4482  for (unsigned l2 = 0; l2 < n_p; l2++)
4483  {
4484  finite_element_pt(Nx * (Ny - 1) + (Nz - 1) * Nx * Ny)
4485  ->node_pt(l2 + (n_p - 1) * n_p * n_p) =
4486  finite_element_pt(Nx * (Ny - 2) + (Nz - 1) * Nx * Ny)
4487  ->node_pt((n_p - 1) * n_p + l2 + (n_p - 1) * n_p * n_p);
4488  }
4489 
4490  // Second row of nodes
4491  // First column of nodes
4492  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
4493  {
4494  // Determine number of values
4495  local_node_num = n_p * l1 + (n_p - 1) * n_p * n_p;
4496 
4497  // Allocate memory for node
4498  Node_pt[node_count] =
4499  finite_element_pt(element_num)
4500  ->construct_boundary_node(local_node_num, time_stepper_pt);
4501  // Set the pointer from the element
4502  finite_element_pt(element_num)->node_pt(local_node_num) =
4503  Node_pt[node_count];
4504 
4505  // Get the fractional position of the node
4506  finite_element_pt(element_num)
4507  ->local_fraction_of_node(local_node_num, s_fraction);
4508 
4509  // Set the position of the node
4510  Node_pt[node_count]->x(0) = Xmin;
4511  Node_pt[node_count]->x(1) =
4512  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
4513  Node_pt[node_count]->x(2) = Zmax;
4514 
4515  // Add the node to the boundaries 4 and 5
4516  add_boundary_node(4, Node_pt[node_count]);
4517  add_boundary_node(5, Node_pt[node_count]);
4518 
4519  // Increment the node number
4520  node_count++;
4521 
4522  // Now do the other columns
4523  for (unsigned l2 = 1; l2 < n_p; l2++)
4524  {
4525  // Determine number of values
4526  local_node_num = n_p * l1 + l2 + (n_p - 1) * n_p * n_p;
4527 
4528  // Allocate memory for node
4529  Node_pt[node_count] =
4530  finite_element_pt(element_num)
4531  ->construct_boundary_node(local_node_num, time_stepper_pt);
4532  // Set the pointer from the element
4533  finite_element_pt(element_num)->node_pt(local_node_num) =
4534  Node_pt[node_count];
4535  // Get the fractional position of the node
4536  finite_element_pt(element_num)
4537  ->local_fraction_of_node(local_node_num, s_fraction);
4538 
4539  // Set the position of the node
4540  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
4541  Node_pt[node_count]->x(1) =
4542  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
4543  Node_pt[node_count]->x(2) = Zmax;
4544 
4545  // Add the node to the boundary 5
4546  add_boundary_node(5, Node_pt[node_count]);
4547 
4548  // Increment the node number
4549  node_count++;
4550  }
4551  }
4552 
4553  // Final row of nodes
4554  // First column of nodes
4555  // Top left node
4556  // Determine number of values
4557  local_node_num = n_p * (n_p - 1) + (n_p - 1) * n_p * n_p;
4558  // Allocate memory for node
4559  Node_pt[node_count] =
4560  finite_element_pt(element_num)
4561  ->construct_boundary_node(local_node_num, time_stepper_pt);
4562  // Set the pointer from the element
4563  finite_element_pt(element_num)->node_pt(local_node_num) =
4564  Node_pt[node_count];
4565  // Get the fractional position of the node
4566  finite_element_pt(element_num)
4567  ->local_fraction_of_node(local_node_num, s_fraction);
4568 
4569  // Set the position of the node
4570  Node_pt[node_count]->x(0) = Xmin;
4571  Node_pt[node_count]->x(1) = Ymax;
4572  Node_pt[node_count]->x(2) = Zmax;
4573 
4574  // Add the node to the boundaries
4575  add_boundary_node(3, Node_pt[node_count]);
4576  add_boundary_node(4, Node_pt[node_count]);
4577  add_boundary_node(5, Node_pt[node_count]);
4578 
4579  // Increment the node number
4580  node_count++;
4581 
4582  // Now do the other columns
4583  for (unsigned l2 = 1; l2 < n_p; l2++)
4584  {
4585  // Determine number of values
4586  local_node_num = n_p * (n_p - 1) + l2 + (n_p - 1) * n_p * n_p;
4587  // Allocate memory for the node
4588  Node_pt[node_count] =
4589  finite_element_pt(element_num)
4590  ->construct_boundary_node(local_node_num, time_stepper_pt);
4591  // Set the pointer from the element
4592  finite_element_pt(element_num)->node_pt(local_node_num) =
4593  Node_pt[node_count];
4594 
4595  // Get the fractional position of the node
4596  finite_element_pt(element_num)
4597  ->local_fraction_of_node(local_node_num, s_fraction);
4598 
4599  // Set the position of the node
4600  Node_pt[node_count]->x(0) = Xmin + el_length[0] * s_fraction[0];
4601  Node_pt[node_count]->x(1) = Ymax;
4602  Node_pt[node_count]->x(2) = Zmax;
4603 
4604  // Add the node to the boundaries 3 and 5
4605  add_boundary_node(3, Node_pt[node_count]);
4606  add_boundary_node(5, Node_pt[node_count]);
4607 
4608  // Increment the node number
4609  node_count++;
4610  }
4611 
4612 
4613  // Now loop over the rest of the elements in the row, apart from the last
4614  for (unsigned j = 1; j < (Nx - 1); j++)
4615  {
4616  // Allocate memory for element
4617  element_num = Nx * (Ny - 1) + j + (Nz - 1) * Nx * Ny;
4618  Element_pt[element_num] = new ELEMENT;
4619 
4620  // The lowest nodes are copied from the lower element
4621  for (unsigned l1 = 0; l1 < n_p; l1++)
4622  {
4623  for (unsigned l2 = 0; l2 < n_p; l2++)
4624  {
4625  finite_element_pt(Nx * (Ny - 1) + j + (Nz - 1) * Nx * Ny)
4626  ->node_pt(l2 + n_p * l1) =
4627  finite_element_pt(Nx * (Ny - 1) + j + (Nz - 2) * Nx * Ny)
4628  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
4629  }
4630  }
4631 
4632  // Jump in the third dimension but the top nodes
4633 
4634  for (unsigned l3 = 1; l3 < (n_p - 1); l3++)
4635  {
4636  // The first row is copied from the lower element
4637  for (unsigned l2 = 0; l2 < n_p; l2++)
4638  {
4639  finite_element_pt(Nx * (Ny - 1) + j + (Nz - 1) * Nx * Ny)
4640  ->node_pt(l2 + l3 * n_p * n_p) =
4641  finite_element_pt(Nx * (Ny - 2) + j + (Nz - 1) * Nx * Ny)
4642  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
4643  }
4644 
4645  // Second rows
4646  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
4647  {
4648  // First column is same as neighbouring element
4649  finite_element_pt(Nx * (Ny - 1) + j + (Nz - 1) * Nx * Ny)
4650  ->node_pt(n_p * l1 + l3 * n_p * n_p) =
4651  finite_element_pt(Nx * (Ny - 1) + (j - 1) + (Nz - 1) * Nx * Ny)
4652  ->node_pt(n_p * l1 + (n_p - 1) + l3 * n_p * n_p);
4653 
4654  // New nodes for other columns
4655  for (unsigned l2 = 1; l2 < n_p; l2++)
4656  {
4657  // Determine number of values
4658  local_node_num = n_p * l1 + l2 + l3 * n_p * n_p;
4659  // Allocate memory for the node
4660  Node_pt[node_count] =
4661  finite_element_pt(element_num)
4662  ->construct_node(local_node_num, time_stepper_pt);
4663 
4664  // Set the pointer
4665  finite_element_pt(element_num)->node_pt(local_node_num) =
4666  Node_pt[node_count];
4667 
4668  // Get the fractional position of the node
4669  finite_element_pt(element_num)
4670  ->local_fraction_of_node(local_node_num, s_fraction);
4671 
4672  // Set the position of the node
4673  Node_pt[node_count]->x(0) =
4674  Xmin + el_length[0] * (j + s_fraction[0]);
4675  Node_pt[node_count]->x(1) =
4676  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
4677  Node_pt[node_count]->x(2) =
4678  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4679 
4680  // No boundaries
4681 
4682  // Increment the node number
4683  // add_boundary_node(0,Node_pt[node_count]);
4684  node_count++;
4685  }
4686  }
4687 
4688  // Top row
4689  // First column is same as neighbouring element
4690  finite_element_pt(Nx * (Ny - 1) + j + (Nz - 1) * Nx * Ny)
4691  ->node_pt(n_p * (n_p - 1) + l3 * n_p * n_p) =
4692  finite_element_pt(Nx * (Ny - 1) + (j - 1) + (Nz - 1) * Nx * Ny)
4693  ->node_pt(n_p * (n_p - 1) + (n_p - 1) + l3 * n_p * n_p);
4694  // New nodes for other columns
4695  for (unsigned l2 = 1; l2 < n_p; l2++)
4696  {
4697  // Determine number of values
4698  local_node_num = n_p * (n_p - 1) + l2 + l3 * n_p * n_p;
4699  // Allocate memory for node
4700  Node_pt[node_count] =
4701  finite_element_pt(element_num)
4702  ->construct_boundary_node(local_node_num, time_stepper_pt);
4703  // Set the pointer
4704  finite_element_pt(element_num)->node_pt(local_node_num) =
4705  Node_pt[node_count];
4706 
4707  // Get the fractional position of the node
4708  finite_element_pt(element_num)
4709  ->local_fraction_of_node(local_node_num, s_fraction);
4710 
4711  // Set the position of the node
4712  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
4713  Node_pt[node_count]->x(1) = Ymax;
4714  Node_pt[node_count]->x(2) =
4715  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4716 
4717  // Add the node to the boundary
4718  add_boundary_node(3, Node_pt[node_count]);
4719 
4720  // Increment the node number
4721  node_count++;
4722  }
4723  }
4724 
4725  // Now the top nodes
4726 
4727  // The first row is copied from the lower element
4728  for (unsigned l2 = 0; l2 < n_p; l2++)
4729  {
4730  finite_element_pt(Nx * (Ny - 1) + j + (Nz - 1) * Nx * Ny)
4731  ->node_pt(l2 + (n_p - 1) * n_p * n_p) =
4732  finite_element_pt(Nx * (Ny - 2) + j + (Nz - 1) * Nx * Ny)
4733  ->node_pt((n_p - 1) * n_p + l2 + (n_p - 1) * n_p * n_p);
4734  }
4735 
4736  // Second rows
4737  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
4738  {
4739  // First column is same as neighbouring element
4740  finite_element_pt(Nx * (Ny - 1) + j + (Nz - 1) * Nx * Ny)
4741  ->node_pt(n_p * l1 + (n_p - 1) * n_p * n_p) =
4742  finite_element_pt(Nx * (Ny - 1) + (j - 1) + (Nz - 1) * Nx * Ny)
4743  ->node_pt(n_p * l1 + (n_p - 1) + (n_p - 1) * n_p * n_p);
4744 
4745  // New nodes for other columns
4746  for (unsigned l2 = 1; l2 < n_p; l2++)
4747  {
4748  // Determine number of values
4749  local_node_num = n_p * l1 + l2 + (n_p - 1) * n_p * n_p;
4750  // Allocate memory for the node
4751  Node_pt[node_count] =
4752  finite_element_pt(element_num)
4753  ->construct_boundary_node(local_node_num, time_stepper_pt);
4754 
4755  // Set the pointer
4756  finite_element_pt(element_num)->node_pt(local_node_num) =
4757  Node_pt[node_count];
4758 
4759  // Get the fractional position of the node
4760  finite_element_pt(element_num)
4761  ->local_fraction_of_node(local_node_num, s_fraction);
4762 
4763  // Set the position of the node
4764  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
4765  Node_pt[node_count]->x(1) =
4766  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
4767  Node_pt[node_count]->x(2) = Zmax;
4768 
4769  // Add the node to the boundary 5
4770  add_boundary_node(5, Node_pt[node_count]);
4771 
4772  // Increment the node number add_boundary_node(0,Node_pt[node_count]);
4773  node_count++;
4774  }
4775  }
4776 
4777  // Top row
4778  // First column is same as neighbouring element
4779  finite_element_pt(Nx * (Ny - 1) + j + (Nz - 1) * Nx * Ny)
4780  ->node_pt(n_p * (n_p - 1) + (n_p - 1) * n_p * n_p) =
4781  finite_element_pt(Nx * (Ny - 1) + (j - 1) + (Nz - 1) * Nx * Ny)
4782  ->node_pt(n_p * (n_p - 1) + (n_p - 1) + (n_p - 1) * n_p * n_p);
4783  // New nodes for other columns
4784  for (unsigned l2 = 1; l2 < n_p; l2++)
4785  {
4786  // Determine number of values
4787  local_node_num = n_p * (n_p - 1) + l2 + (n_p - 1) * n_p * n_p;
4788  // Allocate memory for node
4789  Node_pt[node_count] =
4790  finite_element_pt(element_num)
4791  ->construct_boundary_node(local_node_num, time_stepper_pt);
4792  // Set the pointer
4793  finite_element_pt(element_num)->node_pt(local_node_num) =
4794  Node_pt[node_count];
4795 
4796  // Get the fractional position of the node
4797  finite_element_pt(element_num)
4798  ->local_fraction_of_node(local_node_num, s_fraction);
4799 
4800  // Set the position of the node
4801  Node_pt[node_count]->x(0) = Xmin + el_length[0] * (j + s_fraction[0]);
4802  Node_pt[node_count]->x(1) = Ymax;
4803  Node_pt[node_count]->x(2) = Zmax;
4804 
4805  // Add the node to the boundaries 3 and 5
4806  add_boundary_node(3, Node_pt[node_count]);
4807  add_boundary_node(5, Node_pt[node_count]);
4808 
4809  // Increment the node number
4810  node_count++;
4811  }
4812 
4813 
4814  } // End of loop over central elements in row
4815 
4816 
4817  // LAST ELEMENT (Really!!)
4818  //-----------------------------------------
4819 
4820  // Allocate memory for element
4821  element_num = Nx * (Ny - 1) + Nx - 1 + (Nz - 1) * Nx * Ny;
4822  Element_pt[element_num] = new ELEMENT;
4823 
4824  // The lowest nodes are copied from the lower element
4825  for (unsigned l1 = 0; l1 < n_p; l1++)
4826  {
4827  for (unsigned l2 = 0; l2 < n_p; l2++)
4828  {
4829  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + (Nz - 1) * Nx * Ny)
4830  ->node_pt(l2 + n_p * l1) =
4831  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + (Nz - 2) * Nx * Ny)
4832  ->node_pt(l2 + n_p * l1 + (n_p - 1) * n_p * n_p);
4833  }
4834  }
4835 
4836  // We jump to the third dimension but the top nodes
4837 
4838  for (unsigned l3 = 1; l3 < (n_p - 1); l3++)
4839  {
4840  // The first row is copied from the lower element
4841  for (unsigned l2 = 0; l2 < n_p; l2++)
4842  {
4843  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + (Nz - 1) * Nx * Ny)
4844  ->node_pt(l2 + l3 * n_p * n_p) =
4845  finite_element_pt(Nx * (Ny - 2) + Nx - 1 + (Nz - 1) * Nx * Ny)
4846  ->node_pt((n_p - 1) * n_p + l2 + l3 * n_p * n_p);
4847  }
4848 
4849  // Second rows
4850  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
4851  {
4852  // First column is same as neighbouring element
4853  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + (Nz - 1) * Nx * Ny)
4854  ->node_pt(n_p * l1 + l3 * n_p * n_p) =
4855  finite_element_pt(Nx * (Ny - 1) + Nx - 2 + (Nz - 1) * Nx * Ny)
4856  ->node_pt(n_p * l1 + (n_p - 1) + l3 * n_p * n_p);
4857 
4858  // Middle nodes
4859  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
4860  {
4861  // Determine number of values
4862  local_node_num = n_p * l1 + l2 + l3 * n_p * n_p;
4863  // Allocate memory for node
4864  Node_pt[node_count] =
4865  finite_element_pt(element_num)
4866  ->construct_node(local_node_num, time_stepper_pt);
4867  // Set the pointer
4868  finite_element_pt(element_num)->node_pt(local_node_num) =
4869  Node_pt[node_count];
4870 
4871  // Get the fractional position of the node
4872  finite_element_pt(element_num)
4873  ->local_fraction_of_node(local_node_num, s_fraction);
4874 
4875  // Set the position of the node
4876  Node_pt[node_count]->x(0) =
4877  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
4878  Node_pt[node_count]->x(1) =
4879  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
4880  Node_pt[node_count]->x(2) =
4881  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4882 
4883  // No boundaries
4884 
4885  // Increment the node number
4886  node_count++;
4887  }
4888 
4889  // Final node
4890  // Determine number of values
4891  local_node_num = n_p * l1 + (n_p - 1) + l3 * n_p * n_p;
4892  // Allocate memory for node
4893  Node_pt[node_count] =
4894  finite_element_pt(element_num)
4895  ->construct_boundary_node(local_node_num, time_stepper_pt);
4896  // Set the pointer
4897  finite_element_pt(element_num)->node_pt(local_node_num) =
4898  Node_pt[node_count];
4899 
4900  // Get the fractional position of the node
4901  finite_element_pt(element_num)
4902  ->local_fraction_of_node(local_node_num, s_fraction);
4903 
4904  // Set the position of the node
4905  Node_pt[node_count]->x(0) = Xmax;
4906  Node_pt[node_count]->x(1) =
4907  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
4908  Node_pt[node_count]->x(2) =
4909  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4910 
4911  // Add the node to the boundary 2
4912  add_boundary_node(2, Node_pt[node_count]);
4913 
4914 
4915  // Increment the node number
4916  node_count++;
4917 
4918  } // End of loop over middle rows
4919 
4920  // Final row
4921  // First column is same as neighbouring element
4922  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + (Nz - 1) * Nx * Ny)
4923  ->node_pt(n_p * (n_p - 1) + l3 * n_p * n_p) =
4924  finite_element_pt(Nx * (Ny - 1) + Nx - 2 + (Nz - 1) * Nx * Ny)
4925  ->node_pt(n_p * (n_p - 1) + (n_p - 1) + l3 * n_p * n_p);
4926 
4927  // Middle nodes
4928  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
4929  {
4930  // Determine number of values
4931  local_node_num = n_p * (n_p - 1) + l2 + l3 * n_p * n_p;
4932  // Allocate memory for node
4933  Node_pt[node_count] =
4934  finite_element_pt(element_num)
4935  ->construct_boundary_node(local_node_num, time_stepper_pt);
4936  // Set the pointer
4937  finite_element_pt(element_num)->node_pt(local_node_num) =
4938  Node_pt[node_count];
4939 
4940  // Get the fractional position of the node
4941  finite_element_pt(element_num)
4942  ->local_fraction_of_node(local_node_num, s_fraction);
4943 
4944  // Set the position of the node
4945  Node_pt[node_count]->x(0) =
4946  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
4947  Node_pt[node_count]->x(1) = Ymax;
4948  Node_pt[node_count]->x(2) =
4949  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4950 
4951  // Add the node to the boundary 3
4952  add_boundary_node(3, Node_pt[node_count]);
4953 
4954  // Increment the node number
4955  node_count++;
4956  }
4957 
4958  // Final node
4959  // Determine number of values
4960  local_node_num = n_p * (n_p - 1) + (n_p - 1) + l3 * n_p * n_p;
4961  // Allocate memory for node
4962  Node_pt[node_count] =
4963  finite_element_pt(element_num)
4964  ->construct_boundary_node(local_node_num, time_stepper_pt);
4965  // Set the pointer
4966  finite_element_pt(element_num)->node_pt(local_node_num) =
4967  Node_pt[node_count];
4968 
4969  // Get the fractional position of the node
4970  finite_element_pt(element_num)
4971  ->local_fraction_of_node(local_node_num, s_fraction);
4972 
4973  // Set the position of the node
4974  Node_pt[node_count]->x(0) = Xmax;
4975  // In fluid 2
4976  Node_pt[node_count]->x(1) = Ymax;
4977  Node_pt[node_count]->x(2) =
4978  Zmin + el_length[2] * (Nz - 1 + s_fraction[2]);
4979 
4980  // Add the node to the boundaries 2 and 3
4981  add_boundary_node(2, Node_pt[node_count]);
4982  add_boundary_node(3, Node_pt[node_count]);
4983 
4984  // Increment the node number
4985  node_count++;
4986  }
4987 
4988 
4989  // Now the top nodes
4990 
4991 
4992  // The first row is copied from the lower element
4993  for (unsigned l2 = 0; l2 < n_p; l2++)
4994  {
4995  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + (Nz - 1) * Nx * Ny)
4996  ->node_pt(l2 + (n_p - 1) * n_p * n_p) =
4997  finite_element_pt(Nx * (Ny - 2) + Nx - 1 + (Nz - 1) * Nx * Ny)
4998  ->node_pt((n_p - 1) * n_p + l2 + (n_p - 1) * n_p * n_p);
4999  }
5000 
5001  // Second rows
5002  for (unsigned l1 = 1; l1 < (n_p - 1); l1++)
5003  {
5004  // First column is same as neighbouring element
5005  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + (Nz - 1) * Nx * Ny)
5006  ->node_pt(n_p * l1 + (n_p - 1) * n_p * n_p) =
5007  finite_element_pt(Nx * (Ny - 1) + Nx - 2 + (Nz - 1) * Nx * Ny)
5008  ->node_pt(n_p * l1 + (n_p - 1) + (n_p - 1) * n_p * n_p);
5009 
5010  // Middle nodes
5011  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
5012  {
5013  // Determine number of values
5014  local_node_num = n_p * l1 + l2 + (n_p - 1) * n_p * n_p;
5015  // Allocate memory for node
5016  Node_pt[node_count] =
5017  finite_element_pt(element_num)
5018  ->construct_boundary_node(local_node_num, time_stepper_pt);
5019  // Set the pointer
5020  finite_element_pt(element_num)->node_pt(local_node_num) =
5021  Node_pt[node_count];
5022 
5023  // Get the fractional position of the node
5024  finite_element_pt(element_num)
5025  ->local_fraction_of_node(local_node_num, s_fraction);
5026 
5027  // Set the position of the node
5028  Node_pt[node_count]->x(0) =
5029  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
5030  Node_pt[node_count]->x(1) =
5031  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
5032  Node_pt[node_count]->x(2) = Zmax;
5033 
5034  // Add to boundary 5
5035  add_boundary_node(5, Node_pt[node_count]);
5036 
5037  // Increment the node number
5038  node_count++;
5039  }
5040 
5041  // Final node
5042  // Determine number of values
5043  local_node_num = n_p * l1 + (n_p - 1) + (n_p - 1) * n_p * n_p;
5044  // Allocate memory for node
5045  Node_pt[node_count] =
5046  finite_element_pt(element_num)
5047  ->construct_boundary_node(local_node_num, time_stepper_pt);
5048  // Set the pointer
5049  finite_element_pt(element_num)->node_pt(local_node_num) =
5050  Node_pt[node_count];
5051 
5052  // Set the position of the node
5053  Node_pt[node_count]->x(0) = Xmax;
5054  Node_pt[node_count]->x(1) =
5055  Ymin + el_length[1] * (Ny - 1 + s_fraction[1]);
5056  Node_pt[node_count]->x(2) = Zmax;
5057 
5058  // Add the node to the boundaries 2 and 5
5059  add_boundary_node(2, Node_pt[node_count]);
5060  add_boundary_node(5, Node_pt[node_count]);
5061 
5062  // Increment the node number
5063  node_count++;
5064 
5065  } // End of loop over middle rows
5066 
5067  // Final row
5068  // First node is same as neighbouring element
5069  finite_element_pt(Nx * (Ny - 1) + Nx - 1 + (Nz - 1) * Nx * Ny)
5070  ->node_pt(n_p * (n_p - 1) + (n_p - 1) * n_p * n_p) =
5071  finite_element_pt(Nx * (Ny - 1) + Nx - 2 + (Nz - 1) * Nx * Ny)
5072  ->node_pt(n_p * (n_p - 1) + (n_p - 1) + (n_p - 1) * n_p * n_p);
5073 
5074  // Middle nodes
5075  for (unsigned l2 = 1; l2 < (n_p - 1); l2++)
5076  {
5077  // Determine number of values
5078  local_node_num = n_p * (n_p - 1) + l2 + (n_p - 1) * n_p * n_p;
5079  // Allocate memory for node
5080  Node_pt[node_count] =
5081  finite_element_pt(element_num)
5082  ->construct_boundary_node(local_node_num, time_stepper_pt);
5083  // Set the pointer
5084  finite_element_pt(element_num)->node_pt(local_node_num) =
5085  Node_pt[node_count];
5086 
5087  // Get the fractional position of the node
5088  finite_element_pt(element_num)
5089  ->local_fraction_of_node(local_node_num, s_fraction);
5090 
5091  // Set the position of the node
5092  Node_pt[node_count]->x(0) =
5093  Xmin + el_length[0] * (Nx - 1 + s_fraction[0]);
5094  Node_pt[node_count]->x(1) = Ymax;
5095  Node_pt[node_count]->x(2) = Zmax;
5096 
5097  // Add the node to the boundary 3
5098  add_boundary_node(3, Node_pt[node_count]);
5099  add_boundary_node(5, Node_pt[node_count]);
5100 
5101 
5102  // Increment the node number
5103  node_count++;
5104  }
5105 
5106  // Final node (really!!)
5107  // Determine number of values
5108  local_node_num = n_p * (n_p - 1) + (n_p - 1) + (n_p - 1) * n_p * n_p;
5109  // Allocate memory for node
5110  Node_pt[node_count] =
5111  finite_element_pt(element_num)
5112  ->construct_boundary_node(local_node_num, time_stepper_pt);
5113  // Set the pointer
5114  finite_element_pt(element_num)->node_pt(local_node_num) =
5115  Node_pt[node_count];
5116 
5117  // Get the fractional position of the node
5118  finite_element_pt(element_num)
5119  ->local_fraction_of_node(local_node_num, s_fraction);
5120 
5121  // Set the position of the node
5122  Node_pt[node_count]->x(0) = Xmax;
5123  Node_pt[node_count]->x(1) = Ymax;
5124  Node_pt[node_count]->x(2) = Zmax;
5125 
5126  // Add the node to the boundaries 2, 3 and 5
5127  add_boundary_node(2, Node_pt[node_count]);
5128  add_boundary_node(3, Node_pt[node_count]);
5129  add_boundary_node(5, Node_pt[node_count]);
5130 
5131  // Increment the node number
5132  node_count++;
5133 
5134 
5135  // Setup lookup scheme that establishes which elements are located
5136  // on the mesh boundaries
5137  setup_boundary_element_info();
5138  }
5139 
5140 
5141 } // namespace oomph
5142 
5143 #endif
oomph::SimpleCubicMesh::build_mesh
void build_mesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
Definition: simple_cubic_mesh.template.cc:40
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: annular_domain.h:35
simple_cubic_mesh.template.h