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

#include <rectangular_quadmesh.template.h>

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

Public Member Functions

 RectangularQuadMesh (const unsigned &nx, const unsigned &ny, const double &lx, const double &ly, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
 RectangularQuadMesh (const unsigned &nx, const unsigned &ny, const double &xmin, const double &xmax, const double &ymin, const double &ymax, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
 RectangularQuadMesh (const unsigned &nx, const unsigned &ny, const double &lx, const double &ly, const bool &periodic_in_x, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
 RectangularQuadMesh (const unsigned &nx, const unsigned &ny, const double &xmin, const double &xmax, const double &ymin, const double &ymax, const bool &periodic_in_x, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
const unsigned & nx () const
 Return number of elements in x direction. More...
 
const unsigned & ny () const
 Return number of elements in y direction. More...
 
const double x_min () const
 Return the minimum value of x coordinate. More...
 
const double x_max () const
 Return the maximum value of x coordinate. More...
 
const double y_min () const
 Return the minimum value of y coordinate. More...
 
const double y_max () const
 Return the maximum value of y coordinate. More...
 
virtual void element_reorder ()
 
virtual double x_spacing_function (unsigned xelement, unsigned xnode, unsigned yelement, unsigned ynode)
 
virtual double y_spacing_function (unsigned xelement, unsigned xnode, unsigned yelement, unsigned ynode)
 
- Public Member Functions inherited from oomph::QuadMeshBase
 QuadMeshBase ()
 Constructor (empty) More...
 
 QuadMeshBase (const QuadMeshBase &node)=delete
 Broken copy constructor. More...
 
void operator= (const QuadMeshBase &)=delete
 Broken assignment operator. More...
 
virtual ~QuadMeshBase ()
 Destructor (empty) 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...
 
 RectangularQuadMesh (const unsigned &nx, const unsigned &ny, const double &xmin, const double &xmax, const double &ymin, const double &ymax, const bool &periodic_in_x, const bool &build, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
- 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
 Nx: number of elements in x-direction. More...
 
unsigned Ny
 Ny: number of elements in y-direction. More...
 
unsigned Np
 Np: number of (linear) points in the element. 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
 Maximum value of y coordinate. More...
 
bool Xperiodic
 
- 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::RectangularQuadMesh< ELEMENT >

RectangularQuadMesh is a two-dimensional mesh of Quad elements with Nx elements in the "x" (horizonal) direction and Ny elements in the "y" (vertical) direction. Two Constructors are provided. The basic constructor assumes that the lower-left-hand corner of the mesh is (0,0) and takes only the arguments, Nx, Ny, Xmax and Ymax. The more complex constructor takes the additional arguments Xmin and Ymin.

This class is designed to be used as a Base class for more complex two dimensional meshes. The virtual functions x_spacing_function() and y_spacing_function() may be overloaded to provide arbitrary node spacing. The default is uniformly spaced nodes in each direction.

It is also possible to make the solution periodic in the x direction.

Constructor & Destructor Documentation

◆ RectangularQuadMesh() [1/5]

template<class ELEMENT >
oomph::RectangularQuadMesh< ELEMENT >::RectangularQuadMesh ( const unsigned &  nx,
const unsigned &  ny,
const double &  xmin,
const double &  xmax,
const double &  ymin,
const double &  ymax,
const bool &  periodic_in_x,
const bool &  build,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inlineprotected

Constructor that allows the specification of minimum and maximum values of x and y coordinates and does not build the mesh This is intend to be used in derived classes that overload the spacing functions. THis is scheduled to be changed, however. The reason why this MUST be done is because the virtual spacing functions cannot be called in the base constructur, because they will not have been overloaded yet!!

103  : Nx(nx),
104  Ny(ny),
105  Xmin(xmin),
106  Xmax(xmax),
107  Ymin(ymin),
108  Ymax(ymax),
109  Xperiodic(periodic_in_x)
110  {
111  if (build)
112  {
113  // Call the generic mesh constructor
114  build_mesh(time_stepper_pt);
115  }
116  }
oomph::RectangularQuadMesh::Nx
unsigned Nx
Nx: number of elements in x-direction.
Definition: rectangular_quadmesh.template.h:63
oomph::RectangularQuadMesh::Ny
unsigned Ny
Ny: number of elements in y-direction.
Definition: rectangular_quadmesh.template.h:65
oomph::RectangularQuadMesh::Xperiodic
bool Xperiodic
Definition: rectangular_quadmesh.template.h:81
oomph::RectangularQuadMesh::ny
const unsigned & ny() const
Return number of elements in y direction.
Definition: rectangular_quadmesh.template.h:231
oomph::RectangularQuadMesh::Ymax
double Ymax
Maximum value of y coordinate.
Definition: rectangular_quadmesh.template.h:77
oomph::RectangularQuadMesh::Xmax
double Xmax
Maximum value of x coordinate.
Definition: rectangular_quadmesh.template.h:72
oomph::RectangularQuadMesh::nx
const unsigned & nx() const
Return number of elements in x direction.
Definition: rectangular_quadmesh.template.h:224
oomph::RectangularQuadMesh::Ymin
double Ymin
Minimum value of y coordinate.
Definition: rectangular_quadmesh.template.h:75
oomph::RectangularQuadMesh::Xmin
double Xmin
Minimum value of x coordinate.
Definition: rectangular_quadmesh.template.h:70
oomph::RectangularQuadMesh::build_mesh
void build_mesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
Definition: rectangular_quadmesh.template.cc:43

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

◆ RectangularQuadMesh() [2/5]

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

Simple constructor: nx: number of elements in x direction; ny: number of elements in y direction; lx, length of domain in x direction (0,lx); ly, length of domain in y direction (0,ly) Also pass pointer to timestepper (defaults to Steady)

129  : Nx(nx),
130  Ny(ny),
131  Xmin(0.0),
132  Xmax(lx),
133  Ymin(0.0),
134  Ymax(ly),
135  Xperiodic(false)
136  {
137  // Mesh can only be built with 2D Qelements.
138  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
139 
140  // Call the generic mesh constructor
141  build_mesh(time_stepper_pt);
142  }
Mesh_Parameters::ly
const double ly
Definition: ConstraintElementsUnitTest.cpp:34
Mesh_Parameters::lx
const double lx
Definition: ConstraintElementsUnitTest.cpp:33

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

◆ RectangularQuadMesh() [3/5]

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

Constructor that allows the specification of minimum and maximum values of x and y coordinates

154  : Nx(nx),
155  Ny(ny),
156  Xmin(xmin),
157  Xmax(xmax),
158  Ymin(ymin),
159  Ymax(ymax),
160  Xperiodic(false)
161  {
162  // Mesh can only be built with 2D Qelements.
163  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
164 
165  // Call the generic mesh constructor
166  build_mesh(time_stepper_pt);
167  }

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

◆ RectangularQuadMesh() [4/5]

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

Simple constructor: nx: number of elements in x direction; ny: number of elements in y direction; lx, length of domain in x direction (0,lx); ly, length of domain in y direction (0,ly) Boolean flag specifies if the mesh is periodic in the x-direction. Also pass pointer to timestepper (defaults to Steady)

181  : Nx(nx),
182  Ny(ny),
183  Xmin(0.0),
184  Xmax(lx),
185  Ymin(0.0),
186  Ymax(ly),
187  Xperiodic(periodic_in_x)
188  {
189  // Mesh can only be built with 2D Qelements.
190  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
191 
192  // Call the generic mesh constructor
193  build_mesh(time_stepper_pt);
194  }

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

◆ RectangularQuadMesh() [5/5]

template<class ELEMENT >
oomph::RectangularQuadMesh< ELEMENT >::RectangularQuadMesh ( const unsigned &  nx,
const unsigned &  ny,
const double &  xmin,
const double &  xmax,
const double &  ymin,
const double &  ymax,
const bool &  periodic_in_x,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline

Constructor that allows the specification of minimum and maximum values of x and y coordinates. Boolean flag specifies if the mesh is periodic in the x-direction.

208  : Nx(nx),
209  Ny(ny),
210  Xmin(xmin),
211  Xmax(xmax),
212  Ymin(ymin),
213  Ymax(ymax),
214  Xperiodic(periodic_in_x)
215  {
216  // Mesh can only be built with 2D Qelements.
217  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
218 
219  // Call the generic mesh constructor
220  build_mesh(time_stepper_pt);
221  }

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

Member Function Documentation

◆ build_mesh()

template<class ELEMENT >
void oomph::RectangularQuadMesh< 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 the "guts" of the mesh generation process, including all the tedious loops, counting and spacing functions. The function should be called in all constuctors of any derived classes.

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

References GlobalParameters::Nx, and GlobalParameters::Ny.

Referenced by oomph::ChannelSpineMesh< ELEMENT >::build_channel_spine_mesh(), oomph::HorizontalSingleLayerSpineMesh< ELEMENT >::build_horizontal_single_layer_mesh(), oomph::TwoLayerPerturbedSpineMesh< ELEMENT >::build_two_layer_mesh(), oomph::TwoLayerSpineMesh< BASE_ELEMENT >::build_two_layer_mesh(), oomph::RectangularQuadMesh< ELEMENT >::RectangularQuadMesh(), oomph::TwoLayerPerturbedSpineMesh< ELEMENT >::TwoLayerPerturbedSpineMesh(), and oomph::TwoLayerSpineMesh< ELEMENT >::TwoLayerSpineMesh().

◆ element_reorder()

template<class ELEMENT >
void oomph::RectangularQuadMesh< ELEMENT >::element_reorder
virtual

Reorder the elements: By default they are ordered in "horizontal" layers (increasing in x, then in y). This function changes this to an ordering in the vertical direction (y first, then x). This is more efficient if a frontal solver is used and the mesh has more elements in the x than the y direction. Can be overloaded in specific derived meshes.

Reorder the elements so they are listed in vertical slices (more efficient during the frontal solution if the domain is long in the x-direction.

Reimplemented in oomph::ChannelSpineMesh< ELEMENT >, oomph::TwoLayerPerturbedSpineMesh< ELEMENT >, and oomph::TwoLayerPerturbedSpineMesh< PERTURBED_ELEMENT >.

835  {
836  // Find out how many elements there are
837  unsigned long Nelement = nelement();
838  // Create a dummy array of elements
839  Vector<FiniteElement*> dummy;
840 
841  // Loop over the elements in horizontal order
842  for (unsigned long j = 0; j < Nx; j++)
843  {
844  // Loop over the elements vertically
845  for (unsigned long i = 0; i < Ny; i++)
846  {
847  // Push back onto the new stack
848  dummy.push_back(finite_element_pt(Nx * i + j));
849  }
850  }
851 
852  // Now copy the reordered elements into the element_pt
853  for (unsigned long e = 0; e < Nelement; e++)
854  {
855  Element_pt[e] = dummy[e];
856  }
857  }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::Mesh::nelement
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590

References e(), i, j, GlobalParameters::Nx, and GlobalParameters::Ny.

◆ nx()

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

Return number of elements in x direction.

225  {
226  // Return the value of Nx
227  return Nx;
228  }

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

Referenced by oomph::BackwardStepQuadMesh< ELEMENT >::BackwardStepQuadMesh(), oomph::CircularCylindricalShellMesh< ELEMENT >::CircularCylindricalShellMesh(), ElasticTwoLayerMesh< ELEMENT >::ElasticTwoLayerMesh(), and oomph::PMLQuadMeshBase< ELEMENT >::pml_locate_zeta().

◆ ny()

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

Return number of elements in y direction.

232  {
233  // Return the value of Ny
234  return Ny;
235  }

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

Referenced by oomph::BackwardStepQuadMesh< ELEMENT >::BackwardStepQuadMesh(), oomph::CircularCylindricalShellMesh< ELEMENT >::CircularCylindricalShellMesh(), and oomph::PMLQuadMeshBase< ELEMENT >::pml_locate_zeta().

◆ x_max()

template<class ELEMENT >
const double oomph::RectangularQuadMesh< ELEMENT >::x_max ( ) const
inline

Return the maximum value of x coordinate.

246  {
247  // Return the value of Xmax
248  return Xmax;
249  }

References oomph::RectangularQuadMesh< ELEMENT >::Xmax.

Referenced by oomph::PMLQuadMeshBase< ELEMENT >::pml_locate_zeta().

◆ x_min()

template<class ELEMENT >
const double oomph::RectangularQuadMesh< ELEMENT >::x_min ( ) const
inline

Return the minimum value of x coordinate.

239  {
240  // Return the value of Xmin
241  return Xmin;
242  }

References oomph::RectangularQuadMesh< ELEMENT >::Xmin.

Referenced by oomph::PMLQuadMeshBase< ELEMENT >::pml_locate_zeta().

◆ x_spacing_function()

template<class ELEMENT >
virtual double oomph::RectangularQuadMesh< ELEMENT >::x_spacing_function ( unsigned  xelement,
unsigned  xnode,
unsigned  yelement,
unsigned  ynode 
)
inlinevirtual

Return the value of the x-coordinate at the node given by the local node number (xnode, ynode) in the element (xelement,yelement). The description is in a "psudeo" two-dimensional coordinate system, so the range of xelement is [0,Nx-1], yelement is [0,Ny-1], and that of xnode and ynode is [0,Np-1]. The default is to return nodes that are equally spaced in the x coodinate.

Reimplemented in oomph::TwoLayerSpineMesh< ELEMENT >, oomph::TwoLayerSpineMesh< BASE_ELEMENT >, oomph::TwoLayerSpineMesh< SpineElement< ELEMENT > >, oomph::ChannelSpineMesh< ELEMENT >, CylinderMesh< ELEMENT >, oomph::TwoLayerPerturbedSpineMesh< ELEMENT >, oomph::TwoLayerPerturbedSpineMesh< PERTURBED_ELEMENT >, CylinderMesh< ELEMENT >, and CylinderMesh< ELEMENT >.

283  {
284  // Calculate the values of equal increments in nodal values
285  double xstep = (Xmax - Xmin) / ((Np - 1) * Nx);
286  // Return the appropriate value
287  return (Xmin + xstep * ((Np - 1) * xelement + xnode));
288  }

References oomph::RectangularQuadMesh< ELEMENT >::Np, oomph::RectangularQuadMesh< ELEMENT >::Nx, oomph::RectangularQuadMesh< ELEMENT >::Xmax, and oomph::RectangularQuadMesh< ELEMENT >::Xmin.

◆ y_max()

template<class ELEMENT >
const double oomph::RectangularQuadMesh< ELEMENT >::y_max ( ) const
inline

Return the maximum value of y coordinate.

260  {
261  // Return the value of Ymax
262  return Ymax;
263  }

References oomph::RectangularQuadMesh< ELEMENT >::Ymax.

Referenced by oomph::PMLQuadMeshBase< ELEMENT >::pml_locate_zeta().

◆ y_min()

template<class ELEMENT >
const double oomph::RectangularQuadMesh< ELEMENT >::y_min ( ) const
inline

Return the minimum value of y coordinate.

253  {
254  // Return the value of Ymin
255  return Ymin;
256  }

References oomph::RectangularQuadMesh< ELEMENT >::Ymin.

Referenced by oomph::PMLQuadMeshBase< ELEMENT >::pml_locate_zeta().

◆ y_spacing_function()

template<class ELEMENT >
virtual double oomph::RectangularQuadMesh< ELEMENT >::y_spacing_function ( unsigned  xelement,
unsigned  xnode,
unsigned  yelement,
unsigned  ynode 
)
inlinevirtual

Return the value of the y-coordinate at the node given by the local node number (xnode, ynode) in the element (xelement,yelement). The description is in a "psudeo" two-dimensional coordinate system, so the range of xelement is [0,Nx-1], yelement is [0,Ny-1], and that of xnode and ynode is [0,Np-1]. The default it to return nodes that are equally spaced in the y coordinate.

Reimplemented in oomph::TwoLayerSpineMesh< ELEMENT >, oomph::TwoLayerSpineMesh< BASE_ELEMENT >, oomph::TwoLayerSpineMesh< SpineElement< ELEMENT > >, CylinderMesh< ELEMENT >, oomph::TwoLayerPerturbedSpineMesh< ELEMENT >, oomph::TwoLayerPerturbedSpineMesh< PERTURBED_ELEMENT >, CylinderMesh< ELEMENT >, and CylinderMesh< ELEMENT >.

300  {
301  double ystep = (Ymax - Ymin) / ((Np - 1) * Ny);
302  // Return the appropriate value
303  return (Ymin + ystep * ((Np - 1) * yelement + ynode));
304  }

References oomph::RectangularQuadMesh< ELEMENT >::Np, oomph::RectangularQuadMesh< ELEMENT >::Ny, oomph::RectangularQuadMesh< ELEMENT >::Ymax, and oomph::RectangularQuadMesh< ELEMENT >::Ymin.

Member Data Documentation

◆ Np

template<class ELEMENT >
unsigned oomph::RectangularQuadMesh< ELEMENT >::Np
protected

Np: number of (linear) points in the element.

Referenced by FlatPlateMesh< ELEMENT >::FlatPlateMesh(), ShellMesh< ELEMENT >::ShellMesh(), oomph::ChannelSpineMesh< ELEMENT >::x_spacing_function(), oomph::RectangularQuadMesh< ELEMENT >::x_spacing_function(), and oomph::RectangularQuadMesh< ELEMENT >::y_spacing_function().

◆ Nx

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

Nx: number of elements in x-direction.

Referenced by FlatPlateMesh< ELEMENT >::FlatPlateMesh(), oomph::RectangularQuadMesh< ELEMENT >::nx(), oomph::CircularCylindricalShellMesh< ELEMENT >::piecewise_linear_axial_bl_stretching_fct(), ShellMesh< ELEMENT >::ShellMesh(), oomph::Refineable_r_mesh< ELEMENT >::stretch_mesh(), and oomph::RectangularQuadMesh< ELEMENT >::x_spacing_function().

◆ Ny

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

Ny: number of elements in y-direction.

Referenced by FlatPlateMesh< ELEMENT >::FlatPlateMesh(), oomph::RectangularQuadMesh< ELEMENT >::ny(), ShellMesh< ELEMENT >::ShellMesh(), and oomph::RectangularQuadMesh< ELEMENT >::y_spacing_function().

◆ Xmax

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

Maximum value of x coordinate.

Referenced by FlatPlateMesh< ELEMENT >::FlatPlateMesh(), oomph::CircularCylindricalShellMesh< ELEMENT >::piecewise_linear_axial_bl_stretching_fct(), ShellMesh< ELEMENT >::ShellMesh(), oomph::Refineable_r_mesh< ELEMENT >::stretch_mesh(), oomph::RectangularQuadMesh< ELEMENT >::x_max(), and oomph::RectangularQuadMesh< ELEMENT >::x_spacing_function().

◆ Xmin

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

Minimum value of x coordinate.

Referenced by FlatPlateMesh< ELEMENT >::FlatPlateMesh(), oomph::CircularCylindricalShellMesh< ELEMENT >::piecewise_linear_axial_bl_stretching_fct(), ShellMesh< ELEMENT >::ShellMesh(), oomph::Refineable_r_mesh< ELEMENT >::stretch_mesh(), oomph::RectangularQuadMesh< ELEMENT >::x_min(), oomph::ChannelSpineMesh< ELEMENT >::x_spacing_function(), and oomph::RectangularQuadMesh< ELEMENT >::x_spacing_function().

◆ Xperiodic

template<class ELEMENT >
bool oomph::RectangularQuadMesh< ELEMENT >::Xperiodic
protected

Boolean variable used to determine whether the mesh is periodic in the x-direction

◆ Ymax

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

Maximum value of y coordinate.

Referenced by oomph::ElasticRefineableRectangularQuadMesh< ELEMENT >::ElasticRefineableRectangularQuadMesh(), FlatPlateMesh< ELEMENT >::FlatPlateMesh(), ShellMesh< ELEMENT >::ShellMesh(), oomph::TwoLayerPerturbedSpineMesh< ELEMENT >::spine_node_update_upper(), oomph::TwoLayerSpineMesh< ELEMENT >::spine_node_update_upper(), oomph::RectangularQuadMesh< ELEMENT >::y_max(), and oomph::RectangularQuadMesh< ELEMENT >::y_spacing_function().

◆ Ymin

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

Minimum value of y coordinate.

Referenced by oomph::ElasticRefineableRectangularQuadMesh< ELEMENT >::ElasticRefineableRectangularQuadMesh(), ElasticTwoLayerMesh< ELEMENT >::ElasticTwoLayerMesh(), FlatPlateMesh< ELEMENT >::FlatPlateMesh(), ShellMesh< ELEMENT >::ShellMesh(), oomph::ChannelSpineMesh< ELEMENT >::spine_node_update(), oomph::SingleLayerSpineMesh< ELEMENT >::spine_node_update(), oomph::TwoLayerPerturbedSpineMesh< ELEMENT >::spine_node_update_lower(), oomph::TwoLayerSpineMesh< ELEMENT >::spine_node_update_lower(), oomph::TwoLayerPerturbedSpineMesh< ELEMENT >::spine_node_update_upper(), oomph::TwoLayerSpineMesh< ELEMENT >::spine_node_update_upper(), oomph::RectangularQuadMesh< ELEMENT >::y_min(), and oomph::RectangularQuadMesh< ELEMENT >::y_spacing_function().

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