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