Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
oomph::SimpleCubicMesh< ELEMENT > Class Template Reference

Simple cubic 3D Brick mesh class. More...

#include <simple_cubic_mesh.template.h>

+ Inheritance diagram for oomph::SimpleCubicMesh< ELEMENT >:

Public Member Functions

 SimpleCubicMesh (const unsigned &nx, const unsigned &ny, const unsigned &nz, const double &lx, const double &ly, const double &lz, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
 SimpleCubicMesh (const unsigned &nx, const unsigned &ny, const unsigned &nz, const double &xmin, const double &xmax, const double &ymin, const double &ymax, const double &zmin, const double &zmax, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
const unsigned & nx () const
 Access function for number of elements in x directions. More...
 
const unsigned & ny () const
 Access function for number of elements in y directions. More...
 
const unsigned & nz () const
 Access function for number of elements in y directions. More...
 
- Public Member Functions inherited from oomph::BrickMeshBase
 BrickMeshBase ()
 Constructor (empty) More...
 
 BrickMeshBase (const BrickMeshBase &)=delete
 Broken copy constructor. More...
 
virtual ~BrickMeshBase ()
 Broken assignment operator. More...
 
void setup_boundary_element_info ()
 
void setup_boundary_element_info (std::ostream &outfile)
 
- Public Member Functions inherited from oomph::Mesh
 Mesh ()
 Default constructor. More...
 
 Mesh (const Vector< Mesh * > &sub_mesh_pt)
 
void merge_meshes (const Vector< Mesh * > &sub_mesh_pt)
 
virtual void reset_boundary_element_info (Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned >> &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
 Virtual function to perform the reset boundary elements info rutines. More...
 
template<class BULK_ELEMENT >
void doc_boundary_coordinates (const unsigned &b, std::ofstream &the_file)
 
virtual void scale_mesh (const double &factor)
 
 Mesh (const Mesh &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const Mesh &)=delete
 Broken assignment operator. More...
 
virtual ~Mesh ()
 Virtual Destructor to clean up all memory. More...
 
void flush_element_and_node_storage ()
 
void flush_element_storage ()
 
void flush_node_storage ()
 
Node *& node_pt (const unsigned long &n)
 Return pointer to global node n. More...
 
Node * node_pt (const unsigned long &n) const
 Return pointer to global node n (const version) More...
 
GeneralisedElement *& element_pt (const unsigned long &e)
 Return pointer to element e. More...
 
GeneralisedElement * element_pt (const unsigned long &e) const
 Return pointer to element e (const version) More...
 
const Vector< GeneralisedElement * > & element_pt () const
 Return reference to the Vector of elements. More...
 
Vector< GeneralisedElement * > & element_pt ()
 Return reference to the Vector of elements. More...
 
FiniteElement * finite_element_pt (const unsigned &e) const
 
Node *& boundary_node_pt (const unsigned &b, const unsigned &n)
 Return pointer to node n on boundary b. More...
 
Node * boundary_node_pt (const unsigned &b, const unsigned &n) const
 Return pointer to node n on boundary b. More...
 
void set_nboundary (const unsigned &nbound)
 Set the number of boundaries in the mesh. More...
 
void remove_boundary_nodes ()
 Clear all pointers to boundary nodes. More...
 
void remove_boundary_nodes (const unsigned &b)
 
void remove_boundary_node (const unsigned &b, Node *const &node_pt)
 Remove a node from the boundary b. More...
 
void add_boundary_node (const unsigned &b, Node *const &node_pt)
 Add a (pointer to) a node to the b-th boundary. More...
 
void copy_boundary_node_data_from_nodes ()
 
bool boundary_coordinate_exists (const unsigned &i) const
 Indicate whether the i-th boundary has an intrinsic coordinate. More...
 
unsigned long nelement () const
 Return number of elements in the mesh. More...
 
unsigned long nnode () const
 Return number of nodes in the mesh. More...
 
unsigned ndof_types () const
 Return number of dof types in mesh. More...
 
unsigned elemental_dimension () const
 Return number of elemental dimension in mesh. More...
 
unsigned nodal_dimension () const
 Return number of nodal dimension in mesh. More...
 
void add_node_pt (Node *const &node_pt)
 Add a (pointer to a) node to the mesh. More...
 
void add_element_pt (GeneralisedElement *const &element_pt)
 Add a (pointer to) an element to the mesh. More...
 
virtual void node_update (const bool &update_all_solid_nodes=false)
 
virtual void reorder_nodes (const bool &use_old_ordering=true)
 
virtual void get_node_reordering (Vector< Node * > &reordering, const bool &use_old_ordering=true) const
 
template<class BULK_ELEMENT , template< class > class FACE_ELEMENT>
void build_face_mesh (const unsigned &b, Mesh *const &face_mesh_pt)
 
unsigned self_test ()
 Self-test: Check elements and nodes. Return 0 for OK. More...
 
void max_and_min_element_size (double &max_size, double &min_size)
 
double total_size ()
 
void check_inverted_elements (bool &mesh_has_inverted_elements, std::ofstream &inverted_element_file)
 
void check_inverted_elements (bool &mesh_has_inverted_elements)
 
unsigned check_for_repeated_nodes (const double &epsilon=1.0e-12)
 
Vector< Node * > prune_dead_nodes ()
 
unsigned nboundary () const
 Return number of boundaries. More...
 
unsigned long nboundary_node (const unsigned &ibound) const
 Return number of nodes on a particular boundary. More...
 
FiniteElement * boundary_element_pt (const unsigned &b, const unsigned &e) const
 Return pointer to e-th finite element on boundary b. More...
 
Node * get_some_non_boundary_node () const
 
unsigned nboundary_element (const unsigned &b) const
 Return number of finite elements that are adjacent to boundary b. More...
 
int face_index_at_boundary (const unsigned &b, const unsigned &e) const
 
virtual void dump (std::ofstream &dump_file, const bool &use_old_ordering=true) const
 Dump the data in the mesh into a file for restart. More...
 
void dump (const std::string &dump_file_name, const bool &use_old_ordering=true) const
 Dump the data in the mesh into a file for restart. More...
 
virtual void read (std::ifstream &restart_file)
 Read solution from restart file. More...
 
void output_paraview (std::ofstream &file_out, const unsigned &nplot) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
void output (std::ostream &outfile)
 Output for all elements. More...
 
void output (std::ostream &outfile, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
void output (FILE *file_pt)
 Output for all elements (C-style output) More...
 
void output (FILE *file_pt, const unsigned &nplot)
 Output at f(n_plot) points in each element (C-style output) More...
 
void output (const std::string &output_filename)
 Output for all elements. More...
 
void output (const std::string &output_filename, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
void output_fct (std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt)
 Output a given Vector function at f(n_plot) points in each element. More...
 
void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt)
 
void output_boundaries (std::ostream &outfile)
 Output the nodes on the boundaries (into separate tecplot zones) More...
 
void output_boundaries (const std::string &output_filename)
 
void assign_initial_values_impulsive ()
 Assign initial values for an impulsive start. More...
 
void shift_time_values ()
 
void calculate_predictions ()
 
void set_nodal_and_elemental_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void set_mesh_level_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void set_consistent_pinned_values_for_continuation (ContinuationStorageScheme *const &continuation_stepper_pt)
 Set consistent values for pinned data in continuation. More...
 
bool does_pointer_correspond_to_mesh_data (double *const &parameter_pt)
 Does the double pointer correspond to any mesh data. More...
 
void set_nodal_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 Set the timestepper associated with the nodal data in the mesh. More...
 
void set_elemental_internal_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void compute_norm (double &norm)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Returns the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
bool is_mesh_distributed () const
 Boolean to indicate if Mesh has been distributed. More...
 
OomphCommunicator * communicator_pt () const
 
void delete_all_external_storage ()
 Wipe the storage for all externally-based elements. More...
 

Protected Member Functions

void build_mesh (TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 Generic mesh construction function: contains all the hard work. More...
 
- Protected Member Functions inherited from oomph::Mesh
unsigned long assign_global_eqn_numbers (Vector< double * > &Dof_pt)
 Assign (global) equation numbers to the nodes. More...
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign local equation numbers in all elements. More...
 
void convert_to_boundary_node (Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
 
void convert_to_boundary_node (Node *&node_pt)
 

Protected Attributes

unsigned Nx
 Number of elements in x direction. More...
 
unsigned Ny
 Number of elements in y direction. More...
 
unsigned Nz
 Number of elements in y direction. More...
 
double Xmin
 Minimum value of x coordinate. More...
 
double Xmax
 Maximum value of x coordinate. More...
 
double Ymin
 Minimum value of y coordinate. More...
 
double Ymax
 Minimum value of y coordinate. More...
 
double Zmin
 Minimum value of z coordinate. More...
 
double Zmax
 Maximum value of z coordinate. More...
 
- Protected Attributes inherited from oomph::Mesh
Vector< Vector< Node * > > Boundary_node_pt
 
bool Lookup_for_elements_next_boundary_is_setup
 
Vector< Vector< FiniteElement * > > Boundary_element_pt
 
Vector< Vector< int > > Face_index_at_boundary
 
Vector< Node * > Node_pt
 Vector of pointers to nodes. More...
 
Vector< GeneralisedElement * > Element_pt
 Vector of pointers to generalised elements. More...
 
std::vector< bool > Boundary_coordinate_exists
 

Additional Inherited Members

- Public Types inherited from oomph::Mesh
typedef void(FiniteElement::* SteadyExactSolutionFctPt) (const Vector< double > &x, Vector< double > &soln)
 
typedef void(FiniteElement::* UnsteadyExactSolutionFctPt) (const double &time, const Vector< double > &x, Vector< double > &soln)
 
- Static Public Attributes inherited from oomph::Mesh
static Steady< 0 > Default_TimeStepper
 The Steady Timestepper. More...
 
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
 Static boolean flag to control warning about mesh level timesteppers. More...
 

Detailed Description

template<class ELEMENT>
class oomph::SimpleCubicMesh< ELEMENT >

Simple cubic 3D Brick mesh class.

Constructor & Destructor Documentation

◆ SimpleCubicMesh() [1/2]

template<class ELEMENT >
oomph::SimpleCubicMesh< ELEMENT >::SimpleCubicMesh ( const unsigned &  nx,
const unsigned &  ny,
const unsigned &  nz,
const double &  lx,
const double &  ly,
const double &  lz,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline

Constructor: Pass number of elements in the x, y, and z directions, and the corresponding dimensions. Assume that the back lower left corner is located at (0,0,0) Timestepper defaults to Steady.

59  : Nx(nx),
60  Ny(ny),
61  Nz(nz),
62  Xmin(0.0),
63  Xmax(lx),
64  Ymin(0.0),
65  Ymax(ly),
66  Zmin(0.0),
67  Zmax(lz)
68  {
69  // Mesh can only be built with 3D Qelements.
70  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
71 
72  // Call the generic build function
73  build_mesh(time_stepper_pt);
74  }
oomph::SimpleCubicMesh::Ny
unsigned Ny
Number of elements in y direction.
Definition: simple_cubic_mesh.template.h:130
oomph::SimpleCubicMesh::ny
const unsigned & ny() const
Access function for number of elements in y directions.
Definition: simple_cubic_mesh.template.h:114
oomph::SimpleCubicMesh::Nz
unsigned Nz
Number of elements in y direction.
Definition: simple_cubic_mesh.template.h:133
oomph::SimpleCubicMesh::Ymax
double Ymax
Minimum value of y coordinate.
Definition: simple_cubic_mesh.template.h:145
oomph::SimpleCubicMesh::Ymin
double Ymin
Minimum value of y coordinate.
Definition: simple_cubic_mesh.template.h:142
oomph::SimpleCubicMesh::Nx
unsigned Nx
Number of elements in x direction.
Definition: simple_cubic_mesh.template.h:127
oomph::SimpleCubicMesh::Xmin
double Xmin
Minimum value of x coordinate.
Definition: simple_cubic_mesh.template.h:136
oomph::SimpleCubicMesh::nx
const unsigned & nx() const
Access function for number of elements in x directions.
Definition: simple_cubic_mesh.template.h:108
oomph::SimpleCubicMesh::nz
const unsigned & nz() const
Access function for number of elements in y directions.
Definition: simple_cubic_mesh.template.h:120
oomph::SimpleCubicMesh::build_mesh
void build_mesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
Definition: simple_cubic_mesh.template.cc:40
oomph::SimpleCubicMesh::Zmax
double Zmax
Maximum value of z coordinate.
Definition: simple_cubic_mesh.template.h:151
oomph::SimpleCubicMesh::Zmin
double Zmin
Minimum value of z coordinate.
Definition: simple_cubic_mesh.template.h:148
oomph::SimpleCubicMesh::Xmax
double Xmax
Maximum value of x coordinate.
Definition: simple_cubic_mesh.template.h:139
Mesh_Parameters::ly
const double ly
Definition: ConstraintElementsUnitTest.cpp:34
Mesh_Parameters::lx
const double lx
Definition: ConstraintElementsUnitTest.cpp:33
Mesh_Parameters::lz
const double lz
Definition: ConstraintElementsUnitTest.cpp:35

References oomph::SimpleCubicMesh< ELEMENT >::build_mesh().

◆ SimpleCubicMesh() [2/2]

template<class ELEMENT >
oomph::SimpleCubicMesh< ELEMENT >::SimpleCubicMesh ( const unsigned &  nx,
const unsigned &  ny,
const unsigned &  nz,
const double &  xmin,
const double &  xmax,
const double &  ymin,
const double &  ymax,
const double &  zmin,
const double &  zmax,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline

Constructor: Pass the number of elements in the x,y and z directions and the correspoding minimum and maximum values of the coordinates in each direction

89  : Nx(nx),
90  Ny(ny),
91  Nz(nz),
92  Xmin(xmin),
93  Xmax(xmax),
94  Ymin(ymin),
95  Ymax(ymax),
96  Zmin(zmin),
97  Zmax(zmax)
98  {
99  // Mesh can only be built with 3D Qelements.
100  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
101 
102  // Call the generic mesh constructor
103  build_mesh(time_stepper_pt);
104  }
Global_Parameters::zmax
double zmax
Definition: time_harmonic_fourier_decomposed_linear_elasticity/cylinder/cylinder.cc:71
Global_Parameters::zmin
double zmin
Definition: time_harmonic_fourier_decomposed_linear_elasticity/cylinder/cylinder.cc:69

References oomph::SimpleCubicMesh< ELEMENT >::build_mesh().

Member Function Documentation

◆ build_mesh()

template<class ELEMENT >
void oomph::SimpleCubicMesh< ELEMENT >::build_mesh ( TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper)
protected

Generic mesh construction function: contains all the hard work.

Generic mesh construction. This function contains all the details of the mesh generation process, including all the tedious loops, counting spacing and boundary functions.

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

References i, j, k, GlobalParameters::Nx, GlobalParameters::Ny, Global_Parameters::Nz, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Global_Parameters::Zmax, and Global_Parameters::Zmin.

Referenced by oomph::SimpleCubicMesh< ELEMENT >::SimpleCubicMesh().

◆ nx()

template<class ELEMENT >
const unsigned& oomph::SimpleCubicMesh< ELEMENT >::nx ( ) const
inline

Access function for number of elements in x directions.

109  {
110  return Nx;
111  }

References oomph::SimpleCubicMesh< ELEMENT >::Nx.

Referenced by CombCanSpineMesh< ELEMENT, INTERFACE_ELEMENT >::add_line_interface_elements(), CombCanSpineMesh< ELEMENT, INTERFACE_ELEMENT >::add_outlet_bulk_elements(), and SolidProblem< ELEMENT_TYPE >::saveSolidMesh().

◆ ny()

template<class ELEMENT >
const unsigned& oomph::SimpleCubicMesh< ELEMENT >::ny ( ) const
inline

Access function for number of elements in y directions.

115  {
116  return Ny;
117  }

References oomph::SimpleCubicMesh< ELEMENT >::Ny.

Referenced by CombCanSpineMesh< ELEMENT, INTERFACE_ELEMENT >::add_line_interface_elements(), and CombCanSpineMesh< ELEMENT, INTERFACE_ELEMENT >::add_outlet_bulk_elements().

◆ nz()

template<class ELEMENT >
const unsigned& oomph::SimpleCubicMesh< ELEMENT >::nz ( ) const
inline

Access function for number of elements in y directions.

121  {
122  return Nz;
123  }

References oomph::SimpleCubicMesh< ELEMENT >::Nz.

Referenced by CombCanSpineMesh< ELEMENT, INTERFACE_ELEMENT >::add_outlet_bulk_elements(), oomph::ElasticRefineableQuarterPipeMesh< ELEMENT >::ElasticRefineableQuarterPipeMesh(), oomph::QuarterPipeMesh< ELEMENT >::QuarterPipeMesh(), and oomph::RefineableQuarterPipeMesh< ELEMENT >::RefineableQuarterPipeMesh().

Member Data Documentation

◆ Nx

template<class ELEMENT >
unsigned oomph::SimpleCubicMesh< ELEMENT >::Nx
protected

Number of elements in x direction.

Referenced by oomph::SimpleCubicMesh< ELEMENT >::nx().

◆ Ny

template<class ELEMENT >
unsigned oomph::SimpleCubicMesh< ELEMENT >::Ny
protected

Number of elements in y direction.

Referenced by oomph::SimpleCubicMesh< ELEMENT >::ny().

◆ Nz

template<class ELEMENT >
unsigned oomph::SimpleCubicMesh< ELEMENT >::Nz
protected

Number of elements in y direction.

Referenced by oomph::SimpleCubicMesh< ELEMENT >::nz().

◆ Xmax

template<class ELEMENT >
double oomph::SimpleCubicMesh< ELEMENT >::Xmax
protected

Maximum value of x coordinate.

◆ Xmin

template<class ELEMENT >
double oomph::SimpleCubicMesh< ELEMENT >::Xmin
protected

Minimum value of x coordinate.

◆ Ymax

template<class ELEMENT >
double oomph::SimpleCubicMesh< ELEMENT >::Ymax
protected

Minimum value of y coordinate.

◆ Ymin

template<class ELEMENT >
double oomph::SimpleCubicMesh< ELEMENT >::Ymin
protected

Minimum value of y coordinate.

◆ Zmax

template<class ELEMENT >
double oomph::SimpleCubicMesh< ELEMENT >::Zmax
protected

Maximum value of z coordinate.

◆ Zmin

template<class ELEMENT >
double oomph::SimpleCubicMesh< ELEMENT >::Zmin
protected

Minimum value of z coordinate.

Referenced by oomph::SingleLayerCubicSpineMesh< ELEMENT >::spine_node_update().

The documentation for this class was generated from the following files: