oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT > Class Template Reference

#include <extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h>

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

Public Member Functions

 ExtrudedCubeMeshFromQuadMesh (QuadMeshBase *quad_mesh_pt, const unsigned &nz, const double &lz, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
 ExtrudedCubeMeshFromQuadMesh (QuadMeshBase *quad_mesh_pt, const unsigned &nz, const double &zmin, const double &zmax, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
 ~ExtrudedCubeMeshFromQuadMesh ()
 
virtual double z_spacing_function (const unsigned &z_element, const unsigned &z_node) const
 
Vector< std::pair< unsigned, int > > get_element_boundary_information (QuadMeshBase *quad_mesh_pt, FiniteElement *quad_el_pt)
 
const unsignednz () const
 Access function for number of elements in z-direction (const version) 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...
 
Nodenode_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...
 
GeneralisedElementelement_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...
 
FiniteElementfinite_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...
 
Nodeboundary_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...
 
FiniteElementboundary_element_pt (const unsigned &b, const unsigned &e) const
 Return pointer to e-th finite element on boundary b. More...
 
Nodeget_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...
 
OomphCommunicatorcommunicator_pt () const
 
void delete_all_external_storage ()
 Wipe the storage for all externally-based elements. More...
 

Protected Member Functions

void build_mesh (QuadMeshBase *quad_mesh_pt, 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 N_node_1d
 The number of nodes in each direction. More...
 
unsigned Nz
 Number of elements in the z-direction. More...
 
double Zmin
 Minimum value of z coordinate. More...
 
double Zmax
 Maximum value of z coordinate. More...
 
Vector< ExtrudedDomain * > Extruded_domain_pt
 Vector of pointers to extruded domain objects. 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< boolBoundary_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::ExtrudedCubeMeshFromQuadMesh< ELEMENT >

Mesh class that takes a 2D mesh consisting of quadrilateral elements and "extrudes" it in the z-direction.

Constructor & Destructor Documentation

◆ ExtrudedCubeMeshFromQuadMesh() [1/2]

template<class ELEMENT >
oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::ExtrudedCubeMeshFromQuadMesh ( QuadMeshBase quad_mesh_pt,
const unsigned nz,
const double lz,
TimeStepper time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline

Constructor: Pass a mesh consisting of quad elements, specify the number of elements in the z direction, and the corresponding length in this direction. Assumes that the back lower left corner is located at (0,0,0). Timestepper defaults to the Steady timestepper.

111  : Nz(nz), Zmin(0.0), Zmax(lz)
112  {
113  // Call the generic build function
114  build_mesh(quad_mesh_pt, time_stepper_pt);
115  }
double Zmin
Minimum value of z coordinate.
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h:198
void build_mesh(QuadMeshBase *quad_mesh_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.cc:99
unsigned Nz
Number of elements in the z-direction.
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h:195
const unsigned & nz() const
Access function for number of elements in z-direction (const version)
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h:184
double Zmax
Maximum value of z coordinate.
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h:201
const double lz
Definition: ConstraintElementsUnitTest.cpp:35

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

◆ ExtrudedCubeMeshFromQuadMesh() [2/2]

template<class ELEMENT >
oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::ExtrudedCubeMeshFromQuadMesh ( QuadMeshBase quad_mesh_pt,
const unsigned nz,
const double zmin,
const double zmax,
TimeStepper time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline

Constructor: Pass a mesh consisting of quad elements, specify the number of elements in the z direction, and the corresponding minimum and maximum z-value of the mesh. Again, timestepper defaults to Steady.

128  : Nz(nz), Zmin(zmin), Zmax(zmax)
129  {
130  // Call the generic mesh constructor
131  build_mesh(quad_mesh_pt, time_stepper_pt);
132  }
double zmax
Definition: time_harmonic_fourier_decomposed_linear_elasticity/cylinder/cylinder.cc:71
double zmin
Definition: time_harmonic_fourier_decomposed_linear_elasticity/cylinder/cylinder.cc:69

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

◆ ~ExtrudedCubeMeshFromQuadMesh()

template<class ELEMENT >
oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::~ExtrudedCubeMeshFromQuadMesh ( )
inline

Destructor: If the underlying spatial domain was made up of any macro elements then we will have created an extruded domain which in turn creates extruded macro elements. As we're responsible for creating the domain, we need to kill it and it'll kill the extruded macro elements...

141  {
142  // Are there any extruded domains?
143  if (Extruded_domain_pt.size() > 0)
144  {
145  // How many extruded domains did we create?
146  unsigned n_extruded_domain = Extruded_domain_pt.size();
147 
148  // Loop over the extruded domains
149  for (unsigned i = 0; i < n_extruded_domain; i++)
150  {
151  // Delete the i-th extruded domain
152  delete Extruded_domain_pt[i];
153 
154  // Make it a null pointer
155  Extruded_domain_pt[i] = 0;
156  }
157  } // if (Extruded_domain_pt.size()>0)
158  } // End of ExtrudedCubeMeshFromQuadMesh
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Vector< ExtrudedDomain * > Extruded_domain_pt
Vector of pointers to extruded domain objects.
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h:204

References oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::Extruded_domain_pt, and i.

Member Function Documentation

◆ build_mesh()

template<class ELEMENT >
void oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::build_mesh ( QuadMeshBase quad_mesh_pt,
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. NOTE: The boundary number of the extruded mesh will follow the same numbering as the input quad mesh. The newly created "front" and "back" face of the 3D mesh will added after those boundaries. For example, if the input mesh has 4 boundaries; b_0, b_1, b_2 & b_3 then the extruded mesh will have 6 boundaries; eb_0, eb_1, eb_2, eb_3, eb_4 & eb_5 where eb_4 is the "front" face and eb5 is the "back" face. The boundaries eb_i here satisfy eb_i = (b_i) x [Zmin,Zmax].

/ Check if the node that we're about to construct lies on the

101  {
102  // Mesh can only be built with 3D Qelements.
103  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
104 
105  // Need at least 2 elements in the z-direction
106  if (Nz == 0)
107  {
108  // Create an ostringstream object to create an error message
109  std::ostringstream error_message;
110 
111  // Create the error message
112  error_message << "Extrude2DQuadMeshTo3DCubeMesh needs at least two "
113  << "elements in each,\ncoordinate direction. You have "
114  << "specified Nz=" << Nz << std::endl;
115 
116  // Throw an error
117  throw OomphLibError(
118  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
119  }
120 
121  // Get the current (i.e. start) time
122  double t_mesh_setup_start = TimingHelpers::timer();
123 
124  // Get the current (i.e. start) time
125  double t_mesh_construct_start = TimingHelpers::timer();
126 
127  // Get the number of boundaries in the input mesh
128  unsigned n_boundary_in_quad_mesh = quad_mesh_pt->nboundary();
129 
130  // The number of boundaries in the extruded mesh
131  unsigned n_boundary_in_extruded_mesh = n_boundary_in_quad_mesh + 2;
132 
133  // Set the number of boundaries in the extruded mesh (just need
134  // to add the front and back face to the number of boundaries)
135  set_nboundary(n_boundary_in_extruded_mesh);
136 
137  // Get the number of elements in the input mesh
138  unsigned n_element_in_quad_mesh = quad_mesh_pt->nelement();
139 
140  // Allocate storage for the boundary information of each quad element.
141  // Takes up more memory but avoids recomputing it many, many times
142  Vector<Vector<std::pair<unsigned, int>>> boundary_information(
143  n_element_in_quad_mesh);
144 
145  // Loop over all the elements in the quad mesh
146  for (unsigned e = 0; e < n_element_in_quad_mesh; e++)
147  {
148  // Get a pointer to the j-th element in the mesh
149  FiniteElement* quad_el_pt = quad_mesh_pt->finite_element_pt(e);
150 
151  // Get the boundary information
152  boundary_information[e] =
153  get_element_boundary_information(quad_mesh_pt, quad_el_pt);
154  }
155 
156  // Allocate storage for the number of elements in the extruded mesh
157  Element_pt.resize(n_element_in_quad_mesh * Nz);
158 
159  // Create the first element
160  ELEMENT* temp_el_pt = new ELEMENT;
161 
162  // Read out the number of linear points in the element (in one direction)
163  unsigned n_node_1d = temp_el_pt->nnode_1d();
164 
165  // Store the variable n_node_1d as the private variable, Np
166  N_node_1d = n_node_1d;
167 
168  // Delete the element
169  delete temp_el_pt;
170 
171  // Make it a null pointer
172  temp_el_pt = 0;
173 
174  // Need the same number of nodes in one direction in the 2D and 3D element
175  if ((n_node_1d != quad_mesh_pt->finite_element_pt(0)->nnode_1d()) ||
176  (n_node_1d * n_node_1d != quad_mesh_pt->finite_element_pt(0)->nnode()))
177  {
178  // Create an ostringstream object to create an error message
179  std::ostringstream error_message;
180 
181  // Create the error message
182  error_message
183  << "Extrude2DQuadMeshTo3DCubeMesh needs the number of "
184  << "nodes (in the 2D element) in one direction\n to match "
185  << "the number of nodes in one direction in the 3D element!";
186 
187  // Throw an error
188  throw OomphLibError(
189  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
190  }
191 
192  // Get the total number of nodes in the input mesh
193  unsigned long n_node_in_quad_mesh = quad_mesh_pt->nnode();
194 
195  // Loop over the nodes in the spatial mesh
196  for (unsigned i = 0; i < n_node_in_quad_mesh; i++)
197  {
198  // Check to see if the i-th node is hanging
199  if (quad_mesh_pt->node_pt(i)->is_hanging())
200  {
201  // Create an output stream
202  std::ostringstream warning_message_stream;
203 
204  // Create an error message
205  warning_message_stream << "Extrusion machinery is not currently able "
206  << "to deal with hanging nodes!" << std::endl;
207 
208  // Throw an error
209  throw OomphLibError(warning_message_stream.str(),
212  }
213  } // for (unsigned i=0;i<n_node_in_quad_mesh;i++)
214 
215  // Number of nodes in the z-direction
216  unsigned n_node_in_z_direction = (1 + (n_node_1d - 1) * Nz);
217 
218  // Get the total number of nodes in the input mesh
219  unsigned long n_node_in_cube_mesh =
220  n_node_in_quad_mesh * n_node_in_z_direction;
221 
222  // Allocate storage for the number of nodes in the extruded mesh
223  Node_pt.resize(n_node_in_cube_mesh);
224 
225  // Counter for the number of nodes in the extruded mesh
226  unsigned long node_count = 0;
227 
228  // Loop over each element slice in the z-direction
229  for (unsigned i = 0; i < Nz; i++)
230  {
231  // Need a map at each element slice to tell us whether or not we've
232  // already created an edge (and thus the nodes on that edge). If the
233  // nodes have already been created then all we have to do is copy them
234  // over in reverse order
235  std::map<std::pair<Edge, double>, Vector<Node*>> edge_to_nodes_map;
236 
237  // As several edges can attach to a single (corner) node we can
238  // accidentally create duplicate nodes (i.e. nodes that lie at
239  // the same spatial position) because the edges aren't covered in
240  // the required order (something we shouldn't rely on anyway!). To
241  // get around this we store the corner nodes of each quad element
242  // with its corresponding z-value. This will in turn return the
243  // corresponding node in the extruded mesh
244  std::map<std::pair<Node*, double>, Node*>
245  quad_corner_node_to_cube_node_map;
246 
247  // Loop over all the elements in one slice
248  for (unsigned e = 0; e < n_element_in_quad_mesh; e++)
249  {
250  // Get a pointer to the e-th element in the input mesh
251  FiniteElement* quad_el_pt = quad_mesh_pt->finite_element_pt(e);
252 
253  // Get the index of the cube element in the extruded mesh
254  unsigned element_index = i * n_element_in_quad_mesh + e;
255 
256  // Create the e-th element in the i-th slice of the extruded mesh
257  ELEMENT* cube_el_pt = new ELEMENT;
258 
259  // Put it in global storage!
260  Element_pt[element_index] = cube_el_pt;
261 
262  // Index of the local node number (in the quad element)
263  unsigned quad_local_node_number = 0;
264 
265  // Index of the local node number (in the cube element)
266  unsigned cube_local_node_number = 0;
267 
268  // Cache the boundary information associated with this element
269  Vector<std::pair<unsigned, int>> el_boundary_information =
270  boundary_information[e];
271 
272  // Loop over the nodes in the third local coordinate direction
273  for (unsigned n_2 = 0; n_2 < n_node_1d; n_2++)
274  {
275  // Calculate the z-value at the n_2-th (z-)node in the i-th element
276  // slice
277  double z_value = z_spacing_function(i, n_2);
278 
279  //---------------------------------------------------------------
280  // At a given node slice, there are three options:
281  // (1) If i>0 and n_2=0 then we already created the nodes
282  // (in the current slice) when we created the nodes in
283  // the previous slice;
284  // (2) If (1) doesn't hold we may have still have created
285  // certain nodes along a face/edge of an element that
286  // we previously created (in this node slice);
287  // (3) No nodes in this node slice have been created yet so
288  // set them all up!
289  // NOTE: if (1) occurs, we can skip to the next node slice but
290  // we can't if (2) occurs (or (3), obviously!).
291  //---------------------------------------------------------------
292  //---------------------------------------------------------------
293  // Case (1): If we're past the first element slice and we're on
294  // the first node slice (in the element) then copy the nodes from
295  // the appropriate element.
296  //---------------------------------------------------------------
297  if ((i > 0) && (n_2 == 0))
298  {
299  // Get the index of the cube element (to copy) in the extruded mesh
300  unsigned copy_element_index = (i - 1) * n_element_in_quad_mesh + e;
301 
302  // Get a pointer to the element in the previous slice
303  // containing the node we want
304  FiniteElement* copy_cube_el_pt =
305  finite_element_pt(copy_element_index);
306 
307  // Storage for the local node number in the element that we're
308  // copying
309  unsigned copy_cube_local_node_number = 0;
310 
311  // Loop over the nodes in the second local coordinate direction
312  for (unsigned n_1 = 0; n_1 < n_node_1d; n_1++)
313  {
314  // Loop over the nodes in the first local coordinate direction
315  for (unsigned n_0 = 0; n_0 < n_node_1d; n_0++)
316  {
317  // Calculate the local node number (in the extruded mesh)
318  cube_local_node_number = n_0 + (n_node_1d * n_1);
319 
320  // Calculate the local node number (in the input mesh)
321  quad_local_node_number = n_0 + (n_node_1d * n_1);
322 
323  // Calculate the local node number (in the element we're
324  // copying it from)
325  copy_cube_local_node_number =
326  cube_local_node_number +
327  n_node_1d * n_node_1d * (n_node_1d - 1);
328 
329  // Copy the node over
330  cube_el_pt->node_pt(cube_local_node_number) =
331  copy_cube_el_pt->node_pt(copy_cube_local_node_number);
332  } // for (unsigned n_0=0;n_0<n_node_1d;n_0++)
333  } // for (unsigned n_1=0;n_1<n_node_1d;n_1++)
334  }
335  //---------------------------------------------------------------
336  // Case (2) & (3): Loop over the edges in the current node slice
337  // and check if any of them have already been set up (and thus
338  // stored in edge_to_nodes_map). If they have then copy them over
339  // and construct all the other nodes. If they haven't then simply
340  // construct all of the nodes.
341  // NOTE: if an edge has been created already, we can't make the
342  // nodes first, assign the coordinates THEN check as it would
343  // be a wasteful operation. Instead, we get the nodes describing
344  // an edge from the input mesh (quad_mesh_pt) and pair it with
345  // the current z-value.
346  //---------------------------------------------------------------
347  else
348  {
349  // Create a vector of bools to tell us if we've set the nodes up
350  // in the current node slice
351  std::vector<bool> has_node_been_setup(n_node_1d * n_node_1d, false);
352 
353  // Number of edges in a slice (N/E/S/W)
354  unsigned n_edge = 4;
355 
356  // Storage for the edges
357  Vector<std::pair<Edge, double>> edges_and_z_value;
358 
359  // Vector to tell us if we are copying any edges (to make sure
360  // we don't store the edge again later!)
361  std::vector<bool> has_edge_been_done(n_edge, false);
362 
363  // What is the last edge? Convert this to an int so the compiler
364  // doesn't complain about comparisons between unsigned and signed
365  // integers...
366  int last_edge = int(QuadTreeNames::N + n_edge);
367 
368  // Loop over the edges
369  for (int i_edge = QuadTreeNames::N; i_edge < last_edge; i_edge++)
370  {
371  // Pointer for the first node associated with this edge
372  Node* node1_pt = 0;
373 
374  // Pointer for the second node associated with this edge
375  Node* node2_pt = 0;
376 
377  // Find the corner nodes associated with this edge
378  switch (i_edge)
379  {
380  case QuadTreeNames::N:
381  // The first node
382  node1_pt = quad_el_pt->node_pt((n_node_1d * n_node_1d) - 1);
383 
384  // The second node
385  node2_pt = quad_el_pt->node_pt(n_node_1d * (n_node_1d - 1));
386 
387  // Break
388  break;
389  case QuadTreeNames::E:
390  // The first node
391  node1_pt = quad_el_pt->node_pt(n_node_1d - 1);
392 
393  // The second node
394  node2_pt = quad_el_pt->node_pt((n_node_1d * n_node_1d) - 1);
395 
396  // Break
397  break;
398  case QuadTreeNames::S:
399  // The first node
400  node1_pt = quad_el_pt->node_pt(0);
401 
402  // The second node
403  node2_pt = quad_el_pt->node_pt(n_node_1d - 1);
404 
405  // Break
406  break;
407  case QuadTreeNames::W:
408  // The first node
409  node1_pt = quad_el_pt->node_pt(n_node_1d * (n_node_1d - 1));
410 
411  // The second node
412  node2_pt = quad_el_pt->node_pt(0);
413 
414  // Break
415  break;
416  default:
417  // Create an ostringstream object to create an error message
418  std::ostringstream error_message_stream;
419 
420  // Create the error message
421  error_message_stream
422  << "Input is " << i_edge << " but it can only\n"
423  << "be either N/E/S/W, i.e. " << QuadTreeNames::N << ","
424  << QuadTreeNames::E << "," << QuadTreeNames::S << " or "
425  << QuadTreeNames::W << std::endl;
426 
427  // Throw an error
428  throw OomphLibError(error_message_stream.str(),
431  }
432 
433  // Create the edge associated with the i_edge-th edge of this
434  // element
435  Edge edge(node1_pt, node2_pt);
436 
437  // Pair it up with the z-value at the current node slice
438  std::pair<Edge, double> edge_and_z_value_pair(edge, z_value);
439 
440  // Store the edge (with its corresponding z-value)
441  edges_and_z_value.push_back(edge_and_z_value_pair);
442 
443  // Is there is a matching entry in the map?
444  std::map<std::pair<Edge, double>, Vector<Node*>>::iterator it =
445  edge_to_nodes_map.find(edge_and_z_value_pair);
446 
447  // If we're not at the end then it has already been created
448  if (it != edge_to_nodes_map.end())
449  {
450  // Get the vector of node pointers
451  Vector<Node*> edge_nodes_pt = it->second;
452 
453 #ifdef PARANOID
454  // Sanity check; make sure enough nodes have been stored
455  if (edge_nodes_pt.size() != n_node_1d)
456  {
457  // Throw an error
458  throw OomphLibError("Not enough nodes in edge_nodes_pt!",
461  }
462 #endif
463 
464  // Remember that the i_edge-th edge has already been created!
465  switch (i_edge)
466  {
467  case QuadTreeNames::N:
468  // Indicate that the Northern edge has been done
469  has_edge_been_done[0] = true;
470 
471  // Break
472  break;
473  case QuadTreeNames::E:
474  // Indicate that the Eastern edge has been done
475  has_edge_been_done[1] = true;
476 
477  // Break
478  break;
479  case QuadTreeNames::S:
480  // Indicate that the Southern edge has been done
481  has_edge_been_done[2] = true;
482 
483  // Break
484  break;
485  case QuadTreeNames::W:
486  // Indicate that the Western edge has been done
487  has_edge_been_done[3] = true;
488 
489  // Break
490  break;
491  default:
492  // Throw an error
493  throw OomphLibError("Incorrect edge type!",
496  }
497 
498  //---------------------------------------------------------------
499  // Assign the nodes appropriately depending on which edge has
500  // already been created.
501  // NOTE: the nodes will have been created in anti-clockwise
502  // order which means we will have to assign them in clockwise
503  // order (from the perspective of the current element).
504  //---------------------------------------------------------------
505  // Loop over the nodes
506  for (unsigned i_node = 0; i_node < n_node_1d; i_node++)
507  {
508  // Calculate the value of quad_local_node_number depending on
509  // which edge the node lies on
510  switch (i_edge)
511  {
512  case QuadTreeNames::N:
513  // Calculate this Node's quad local node number
514  quad_local_node_number =
515  ((n_node_1d * n_node_1d) - 1) - i_node;
516 
517  // Break
518  break;
519  case QuadTreeNames::E:
520  // Calculate this Node's quad local node number
521  quad_local_node_number =
522  (n_node_1d - 1) + (i_node * n_node_1d);
523 
524  // Break
525  break;
526  case QuadTreeNames::S:
527  // Calculate this Node's quad local node number
528  quad_local_node_number = i_node;
529 
530  // Break
531  break;
532  case QuadTreeNames::W:
533  // Calculate this Node's quad local node number
534  quad_local_node_number =
535  (n_node_1d * (n_node_1d - 1)) - (i_node * n_node_1d);
536 
537  // Break
538  break;
539  }
540 
541  // Calculate this Node's cube local node number
542  cube_local_node_number =
543  quad_local_node_number + (n_node_1d * n_node_1d * n_2);
544 
545  // Assign the i_node-th node
546  cube_el_pt->node_pt(cube_local_node_number) =
547  edge_nodes_pt[(n_node_1d - 1) - i_node];
548 
549  // Indicate that the appropriate node has been set up
550  has_node_been_setup[quad_local_node_number] = true;
551  } // for (unsigned i_node=0;i_node<n_node_1d;i_node++)
552  } // if (it!=edge_to_nodes_map.end())
553  } // for (unsigned i_edge=0;i<n_edge;i_edge++)
554 
555  // The number of corners
556  unsigned n_corner = 4;
557 
558  // Loop over the corner nodes
559  for (unsigned i_corner = 0; i_corner < n_corner; i_corner++)
560  {
561  // Calculate the value of quad_local_node_number depending on
562  // which corner the node lies on
563  switch (i_corner)
564  {
565  case 0:
566  // Bottom-left corner
567  quad_local_node_number = 0;
568 
569  // Break
570  break;
571  case 1:
572  // Bottom-right corner
573  quad_local_node_number = n_node_1d - 1;
574 
575  // Break
576  break;
577  case 2:
578  // Top-left corner
579  quad_local_node_number = n_node_1d * (n_node_1d - 1);
580 
581  // Break
582  break;
583  case 3:
584  // Top-right corner
585  quad_local_node_number = (n_node_1d * n_node_1d) - 1;
586 
587  // Break
588  break;
589  }
590 
591  // If the node hasn't already been copied from some other edge
592  if (!has_node_been_setup[quad_local_node_number])
593  {
594  // Get a pointer to the first corner node (in the quad mesh)
595  Node* quad_corner_node_pt =
596  quad_el_pt->node_pt(quad_local_node_number);
597 
598  // Pair up the corner node and z-value
599  std::pair<Node*, double> quad_corner_node_and_z_value(
600  quad_corner_node_pt, z_value);
601 
602  // Have we made this Node yet?
603  std::map<std::pair<Node*, double>, Node*>::iterator it =
604  quad_corner_node_to_cube_node_map.find(
605  quad_corner_node_and_z_value);
606 
607  // If we found a corresponding entry
608  if (it != quad_corner_node_to_cube_node_map.end())
609  {
610  // Index of the local node number
611  cube_local_node_number =
612  quad_local_node_number + (n_node_1d * n_node_1d * n_2);
613 
614  // Get the node and store it in the element
615  cube_el_pt->node_pt(cube_local_node_number) = it->second;
616 
617  // Indicate that the node has now been created
618  has_node_been_setup[quad_local_node_number] = true;
619  }
620  } // if (!has_node_been_setup[quad_local_node_number])
621  } // for (unsigned i_corner=0;i_corner<n_corner;i_corner++)
622 
623  //------------------------------------------------------------------
624  // By this point, all nodes which have already been created have
625  // been copied over. All the other nodes now need to be constructed.
626  //------------------------------------------------------------------
627  // Loop over the nodes in the second local coordinate direction
628  for (unsigned n_1 = 0; n_1 < n_node_1d; n_1++)
629  {
630  // Loop over the nodes in the first local coordinate direction
631  for (unsigned n_0 = 0; n_0 < n_node_1d; n_0++)
632  {
633  // Calculate the local node number (in the input mesh)
634  quad_local_node_number = n_0 + (n_node_1d * n_1);
635 
636  // Calculate the local node number (in the extruded mesh)
637  cube_local_node_number =
638  quad_local_node_number + (n_node_1d * n_node_1d * n_2);
639 
640  // Check if the node has already been created (i.e. by copying)
641  if (has_node_been_setup[quad_local_node_number])
642  {
643  // If this node has already been set up then skip it
644  continue;
645  }
646 
647  // Pointer to point to the node we are about to create
648  Node* cube_node_pt = 0;
649 
650  //--------------------------------------------------------------
653  // Front or back face of the mesh. The node can only lie on the
654  // front face if we're on the first element slice and first node
655  // slice. In contrast, the node can only lie on the back face if
656  // we're on the final element slice and final node slice.
657  //--------------------------------------------------------------
658  // The node lies on the front face or the back face
659  if (((i == 0) && (n_2 == 0)) ||
660  ((i == Nz - 1) && (n_2 == n_node_1d - 1)))
661  {
662  // The node lies on the front face
663  if ((i == 0) && (n_2 == 0))
664  {
665  // Create the boundary node
666  cube_node_pt = cube_el_pt->construct_boundary_node(
667  cube_local_node_number, time_stepper_pt);
668 
669  // Indicate that the node has been created now
670  has_node_been_setup[quad_local_node_number] = true;
671 
672  // Add the node to the appropriate boundary
673  add_boundary_node(n_boundary_in_quad_mesh, cube_node_pt);
674  }
675  // The node lies on the back face
676  else if ((i == Nz - 1) && (n_2 == n_node_1d - 1))
677  {
678  // Create the boundary node
679  cube_node_pt = cube_el_pt->construct_boundary_node(
680  cube_local_node_number, time_stepper_pt);
681 
682  // Indicate that the node has been created now
683  has_node_been_setup[quad_local_node_number] = true;
684 
685  // Add the node to the appropriate boundary
686  add_boundary_node(n_boundary_in_quad_mesh + 1,
687  cube_node_pt);
688  } // if ((i==0)&&(n_2==0))
689 
690  // Get a pointer to the matching quad mesh node
691  Node* quad_node_pt =
692  quad_el_pt->node_pt(quad_local_node_number);
693 
694  // If this Node lies on a boundary
695  if (quad_node_pt->is_on_boundary())
696  {
697  // Storage for the set of boundaries that this Node lies on
698  std::set<unsigned>* boundaries_pt;
699 
700  // Get the boundaries that this Node lies on
701  quad_node_pt->get_boundaries_pt(boundaries_pt);
702 
703  // Create an iterator to loop over the entries of the set
704  std::set<unsigned>::const_iterator const_it;
705 
706  // Iterate over the boundaries
707  for (const_it = boundaries_pt->begin();
708  const_it != boundaries_pt->end();
709  const_it++)
710  {
711  // Add the node to the appropriate boundary (if it hasn't
712  // already been added)
713  add_boundary_node(*const_it, cube_node_pt);
714  }
715  } // quad_node_pt->is_on_boundary()
716  } // if (((i==0)&&(n_2==0))||((i==Nz-1)&&(n_2==n_node_1d-1)))
717 
718 
719  // Check if there is any other boundary information
720  if (el_boundary_information.size() > 0)
721  {
722  // How many boundaries are there to cover?
723  unsigned n_boundaries = el_boundary_information.size();
724 
725  // Loop over the boundaries
726  for (unsigned i_bound = 0; i_bound < n_boundaries; i_bound++)
727  {
728  // If the node still hasn't been created yet
729  if (!has_node_been_setup[quad_local_node_number])
730  {
731  // Create the boundary node
732  cube_node_pt = cube_el_pt->construct_boundary_node(
733  cube_local_node_number, time_stepper_pt);
734 
735  // Indicate that the node has been created now
736  has_node_been_setup[quad_local_node_number] = true;
737  }
738 
739  // Get the boundary information associated with the
740  // i_bound-th entry
741  std::pair<unsigned, int> boundary_and_face_index =
742  el_boundary_information[i_bound];
743 
744  // Decide whether or not the node lies on the boundary. If
745  // it does then add it to the boundary
746  switch (boundary_and_face_index.second)
747  {
748  case -1:
749  // If the s[0]=-1 boundary of the element lies on the
750  // boundary the node can only lie on the boundary if
751  // n_0=0.
752  if (n_0 == 0)
753  {
754  // Now add the node to the appropriate boundary
755  add_boundary_node(boundary_and_face_index.first,
756  cube_node_pt);
757  }
758 
759  // Break
760  break;
761  case 1:
762  // If the s[0]=1 boundary of the element lies on the
763  // boundary the node can only lie on the boundary if
764  // n_0=n_node_1d-1
765  if (n_0 == n_node_1d - 1)
766  {
767  // Now add the node to the appropriate boundary
768  add_boundary_node(boundary_and_face_index.first,
769  cube_node_pt);
770  }
771 
772  // Break
773  break;
774  case -2:
775  // If the s[1]=-1 boundary of the element lies on the
776  // boundary the node can only lie on the boundary if
777  // n_1=0.
778  if (n_1 == 0)
779  {
780  // Now add the node to the appropriate boundary
781  add_boundary_node(boundary_and_face_index.first,
782  cube_node_pt);
783  }
784 
785  // Break
786  break;
787  case 2:
788  // If the s[1]=1 boundary of the element lies on the
789  // boundary the node can only lie on the boundary if
790  // n_1=n_node_1d-1.
791  if (n_1 == n_node_1d - 1)
792  {
793  // Now add the node to the appropriate boundary
794  add_boundary_node(boundary_and_face_index.first,
795  cube_node_pt);
796  }
797 
798  // Break
799  break;
800  } // switch (boundary_and_face_index.second)
801  } // for (unsigned i_bound=0;i_bound<n_boundaries;i_bound++)
802  } // if (el_boundary_information.size()>0)
803 
804  // Okay, so this Node does not lie on a face of a quad element
805  // that lies on a boundary. However(!), there is still another
806  // case in which it might lie on a mesh boundary; when the Node
807  // is the ONLY Node in the element that lies on the boundary.
808  // It's a slightly odd scenario but can easily occur on meshes
809  // created using the QuadFromTriangleMesh class.
810  Node* quad_node_pt =
811  quad_el_pt->node_pt(quad_local_node_number);
812 
813  // If this Node lies on a boundary
814  if (quad_node_pt->is_on_boundary())
815  {
816  // If we haven't created the Node yet AND the corresponding
817  // Node in the quad mesh definitely lies on the boundary
818  if (cube_node_pt == 0)
819  {
820  // Create the boundary node
821  cube_node_pt = cube_el_pt->construct_boundary_node(
822  cube_local_node_number, time_stepper_pt);
823 
824  // Indicate that the node has been created now
825  has_node_been_setup[quad_local_node_number] = true;
826  }
827 
828  // Storage for the set of boundaries that this Node lies on
829  std::set<unsigned>* boundaries_pt;
830 
831  // Get the boundaries that this Node lies on
832  quad_node_pt->get_boundaries_pt(boundaries_pt);
833 
834  // Create an iterator to loop over the entries of the set
835  std::set<unsigned>::const_iterator const_it;
836 
837  // Iterate over the boundaries
838  for (const_it = boundaries_pt->begin();
839  const_it != boundaries_pt->end();
840  const_it++)
841  {
842  // Add the node to the appropriate boundary (if it hasn't
843  // already been added)
844  add_boundary_node(*const_it, cube_node_pt);
845  }
846  } // quad_node_pt->is_on_boundary()
847 
848  // If the node still hasn't been created yet
849  if (!has_node_been_setup[quad_local_node_number])
850  {
851  // The only other possibility is that it is a regular node
852  cube_node_pt = cube_el_pt->construct_node(
853  cube_local_node_number, time_stepper_pt);
854 
855  // Indicate that the node has now been set up
856  has_node_been_setup[quad_local_node_number] = true;
857  }
858 
859  // Set the x-coordinate of the node
860  cube_node_pt->x(0) = quad_node_pt->x(0);
861 
862  // Set the y-coordinate of the node
863  cube_node_pt->x(1) = quad_node_pt->x(1);
864 
865  // Set the z-coordinate of the node
866  cube_node_pt->x(2) = z_value;
867 
868  // Add it to global storage
869  Node_pt[node_count] = cube_node_pt;
870 
871  // Increment the counter
872  node_count++;
873  } // for (unsigned n_0=0;n_0<n_node_1d;n_0++)
874  } // for (unsigned n_1=0;n_1<n_node_1d;n_1++)
875 
876  //------------------------------------------------------------------
877  // We've set up the nodes in one slice, now set up the edge-to-node
878  // information and store it in the map edge_to_nodes_map.
879  //------------------------------------------------------------------
880  for (unsigned i_edge = 0; i_edge < n_edge; i_edge++)
881  {
882  // If it hasn't already been done
883  if (!has_edge_been_done[i_edge])
884  {
885  // Create a vector to store the nodes on this edge
886  Vector<Node*> edge_nodes_pt(n_node_1d, 0);
887 
888  // Loop over the nodes and add them
889  for (unsigned i_node = 0; i_node < n_node_1d; i_node++)
890  {
891  // Calculate the value of quad_local_node_number depending on
892  // which corner the node lies on
893  switch (i_edge)
894  {
895  case 0:
896  // Northern edge
897  quad_local_node_number =
898  ((n_node_1d * n_node_1d) - 1) - i_node;
899 
900  // Break
901  break;
902  case 1:
903  // Eastern edge
904  quad_local_node_number =
905  (n_node_1d - 1) + (i_node * n_node_1d);
906 
907  // Break
908  break;
909  case 2:
910  // Southern edge
911  quad_local_node_number = i_node;
912 
913  // Break
914  break;
915  case 3:
916  // Western edge
917  quad_local_node_number =
918  (n_node_1d * (n_node_1d - 1)) - (i_node * n_node_1d);
919 
920  // Break
921  break;
922  }
923 
924  // Index of the local node number
925  cube_local_node_number =
926  quad_local_node_number + (n_node_1d * n_node_1d * n_2);
927 
928  // Store the i_node-th node along this edge
929  edge_nodes_pt[i_node] =
930  cube_el_pt->node_pt(cube_local_node_number);
931  } // for (unsigned i_node=0;i_node<n_node_1d;i_node++)
932 
933  // Store the nodes along the north edge of the element
934  edge_to_nodes_map[edges_and_z_value[i_edge]] = edge_nodes_pt;
935  } // if (!has_edge_been_done[i_edge])
936  } // for (unsigned i_edge=0;i_edge<n_edge;i_edge++)
937 
938  //--------------------------------------------------------------------
939  // Now store the corner nodes of the element at the current node
940  // slice in the map quad_corner_node_to_cube_node_map.
941  //--------------------------------------------------------------------
942  // Loop over the corner nodes
943  for (unsigned i_corner = 0; i_corner < n_corner; i_corner++)
944  {
945  // Calculate the value of quad_local_node_number depending on
946  // which corner the node lies on
947  switch (i_corner)
948  {
949  case 0:
950  // Bottom-left corner
951  quad_local_node_number = 0;
952 
953  // Break
954  break;
955  case 1:
956  // Bottom-right corner
957  quad_local_node_number = n_node_1d - 1;
958 
959  // Break
960  break;
961  case 2:
962  // Top-left corner
963  quad_local_node_number = n_node_1d * (n_node_1d - 1);
964 
965  // Break
966  break;
967  case 3:
968  // Top-right corner
969  quad_local_node_number = (n_node_1d * n_node_1d) - 1;
970 
971  // Break
972  break;
973  }
974 
975  // Index of the corresponding corner node (in the cube element)
976  cube_local_node_number =
977  quad_local_node_number + (n_node_1d * n_node_1d * n_2);
978 
979  // Get a pointer to the first corner node (in the quad mesh)
980  Node* quad_corner_node_pt =
981  quad_el_pt->node_pt(quad_local_node_number);
982 
983  // Get a pointer to the first corner node (in the cube mesh)
984  Node* cube_corner_node_pt =
985  cube_el_pt->node_pt(cube_local_node_number);
986 
987  // Pair up the corner node and z-value
988  std::pair<Node*, double> quad_corner_node_and_z_value(
989  quad_corner_node_pt, z_value);
990 
991  // Store it in the map
992  quad_corner_node_to_cube_node_map[quad_corner_node_and_z_value] =
993  cube_corner_node_pt;
994  } // for (unsigned i_corner=0;i_corner<n_corner;i_corner++)
995 
996 #ifdef PARANOID
997  // Loop over all the nodes in the current node slice
998  for (unsigned i_node = 0; i_node < (n_node_1d * n_node_1d);
999  i_node++)
1000  {
1001  // Make sure every node in the current slice has been set up
1002  if (!has_node_been_setup[i_node])
1003  {
1004  // Create an ostringstream object to create an error message
1005  std::ostringstream error_message_stream;
1006 
1007  // Create the error message
1008  error_message_stream << "There are nodes in element " << e
1009  << " which have not been constructed!"
1010  << std::endl;
1011 
1012  // Throw an error
1013  throw OomphLibError(error_message_stream.str(),
1016  }
1017  } // for (unsigned i_node=0;i_node<(n_node_1d*n_node_1d);i_node++)
1018 #endif
1019  } // if ((i>0)&&(n_2==0))
1020  } // for (unsigned n_2=0;n_2<n_node_1d;n_2++)
1021  } // for (unsigned e=0;e<n_element_in_quad_mesh;e++)
1022 
1023  // Sanity check; make sure enough nodes have been stored
1024  if (node_count != (1 + (n_node_1d - 1) * (i + 1)) * n_node_in_quad_mesh)
1025  {
1026  // Create an ostringstream object to create an error message
1027  std::ostringstream error_message_stream;
1028 
1029  // The number of nodes we expect there to be in the mesh
1030  int expected_n_node =
1031  (1 + (n_node_1d - 1) * (i + 1)) * n_node_in_quad_mesh;
1032 
1033  // Difference in node count
1034  int node_diff = expected_n_node - node_count;
1035 
1036  // Create the error message
1037  error_message_stream
1038  << "There are meant to be " << expected_n_node
1039  << " nodes in the extruded mesh by the\nend of slice " << i
1040  << " but you have " << node_count << ". "
1041  << "That's a difference of " << node_diff << "!";
1042 
1043  // Throw an error
1044  throw OomphLibError(error_message_stream.str(),
1047  }
1048  } // for (unsigned i=0;i<Nz;i++)
1049 
1050  //------------------------------------------------------------------------
1051  // Count the number of nodes on the boundaries of the quad mesh; there
1052  // should really be a helper function that does this (hint hint)
1053  //------------------------------------------------------------------------
1054  // Use a set to make sure we don't count edge/corner nodes several times
1055  std::set<const Node*> all_quad_boundary_nodes_pt;
1056 
1057  // The number of boundary nodes in our extruded mesh
1058  int n_quad_mesh_boundary_node = 0;
1059 
1060  // The number of boundaries in the quad mesh
1061  for (unsigned b = 0; b < n_boundary_in_quad_mesh; b++)
1062  {
1063  // Add on the number of nodes on the b-th boundary
1064  unsigned n_boundary_node = quad_mesh_pt->nboundary_node(b);
1065 
1066  // Loop over the nodes on this boundary
1067  for (unsigned n = 0; n < n_boundary_node; n++)
1068  {
1069  // We haven't come across this boundary node yet
1070  if (all_quad_boundary_nodes_pt.find(quad_mesh_pt->boundary_node_pt(
1071  b, n)) == all_quad_boundary_nodes_pt.end())
1072  {
1073  // Count this node
1074  n_quad_mesh_boundary_node++;
1075 
1076  // Now remember it so we don't count it twice
1077  all_quad_boundary_nodes_pt.insert(
1078  quad_mesh_pt->boundary_node_pt(b, n));
1079  }
1080  } // for (unsigned n=0; n<n_boundary_node; n++)
1081  } // for (unsigned b=0; b<n_boundary_in_quad_mesh; b++)
1082 
1083  //------------------------------------------------------------------------
1084  // Count the number of nodes on the boundaries of the extruded cube mesh;
1085  // (recall hint hint from above...)
1086  //------------------------------------------------------------------------
1087  // Number of boundary nodes expected in the extruded mesh; initialise to
1088  // the number of boundary nodes on the front and back faces
1089  int n_expected_boundary_node = 2 * quad_mesh_pt->nnode();
1090 
1091  // Add on the nodes expected on the extruded mesh boundaries (but ignoring
1092  // the boundary nodes on the front and back faces)
1093  n_expected_boundary_node +=
1094  (((n_node_1d - 1) * Nz) - 1) * n_quad_mesh_boundary_node;
1095 
1096  // The number of boundary nodes in our extruded mesh
1097  int n_extruded_mesh_boundary_node = 0;
1098 
1099  // Use a set to make sure we don't count edge/corner nodes several times
1100  std::set<const Node*> all_extruded_boundary_nodes_pt;
1101 
1102  // Loop over the boundaries of the mesh
1103  for (unsigned b = 0; b < n_boundary_in_extruded_mesh; b++)
1104  {
1105  // Add on the number of nodes on the b-th boundary
1106  unsigned n_boundary_node = nboundary_node(b);
1107 
1108  // Loop over the nodes on this boundary
1109  for (unsigned n = 0; n < n_boundary_node; n++)
1110  {
1111  // We haven't come across this boundary node yet
1112  if (all_extruded_boundary_nodes_pt.find(this->boundary_node_pt(b, n)) ==
1113  all_extruded_boundary_nodes_pt.end())
1114  {
1115  // Count this node
1116  n_extruded_mesh_boundary_node++;
1117 
1118  // Now remember it so we don't count it twice
1119  all_extruded_boundary_nodes_pt.insert(this->boundary_node_pt(b, n));
1120  }
1121  } // for (unsigned n=0; n<n_boundary_node; n++)
1122  } // for (unsigned b=0; b<n_boundary_in_extruded_mesh; b++)
1123 
1124  //------------------------------------------------------------------------
1125  // Check that there is the correct number of boundary nodes
1126  //------------------------------------------------------------------------
1127  // Error: we don't have the correct number of boundary nodes
1128  if (n_extruded_mesh_boundary_node != n_expected_boundary_node)
1129  {
1130  // Create an ostringstream object to create an error message
1131  std::ostringstream error_message_stream;
1132 
1133  // Difference in node count
1134  int node_diff = n_expected_boundary_node - n_extruded_mesh_boundary_node;
1135 
1136  // Create the error message
1137  error_message_stream << "There should be " << n_expected_boundary_node
1138  << " boundary nodes in the extruded mesh but there"
1139  << "\nare only " << n_extruded_mesh_boundary_node
1140  << " boundary nodes. That's a difference of "
1141  << node_diff << "!" << std::endl;
1142 
1143  // Throw an error
1144  throw OomphLibError(error_message_stream.str(),
1147  }
1148 
1149  // If the user wishes the mesh setup time to be doc-ed
1150  if (MeshExtrusionHelpers::Mesh_extrusion_helper.doc_mesh_setup_time())
1151  {
1152  // Tell the user
1153  oomph_info << "Time taken to extrude mesh [sec]: "
1154  << TimingHelpers::timer() - t_mesh_construct_start
1155  << std::endl;
1156  }
1157 
1158  // Get the current (i.e. start) time
1159  double t_mesh_macro_start = TimingHelpers::timer();
1160 
1161  //----------------------------------------------------------------
1162  // If an element has a macro-element representation in the quad
1163  // mesh then it need an extruded macro-element representation in
1164  // the extruded mesh so set that up here. Then, use the extruded
1165  // macro-element representation of the elements to move all the
1166  // nodes that need to be moved.
1167  //----------------------------------------------------------------
1168  // Create a map that takes a Domain and returns the corresponding
1169  // ExtrudedDomain object
1170  std::map<Domain*, ExtrudedDomain*> domain_to_extruded_domain_map;
1171 
1172  // Loop over the elements in the input mesh
1173  for (unsigned e = 0; e < n_element_in_quad_mesh; e++)
1174  {
1175  // Get a pointer to the e-th element in the input mesh.
1176  // NOTE: The element is upcast to the QElementBase class so we
1177  // can access the s_macro_ll() and s_macro_ur() functions.
1178  QElementBase* quad_el_pt =
1179  dynamic_cast<QElementBase*>(quad_mesh_pt->finite_element_pt(e));
1180 
1181  // Check if the element has a macro-element representation
1182  if (quad_el_pt->macro_elem_pt() != 0)
1183  {
1184  // Get a pointer to the corresponding Domain object
1185  Domain* domain_pt = quad_el_pt->macro_elem_pt()->domain_pt();
1186 
1187  // How many macro elements are there in this Domain?
1188  unsigned n_macro_element = domain_pt->nmacro_element();
1189 
1190  // Create a pointer to point to the extruded version of the above Domain
1191  ExtrudedDomain* extruded_domain_pt = 0;
1192 
1193  // Create an iterator for the map
1194  std::map<Domain*, ExtrudedDomain*>::iterator it;
1195 
1196  // Check if there is a corresponding entry in the map
1197  it = domain_to_extruded_domain_map.find(domain_pt);
1198 
1199  // Check if the extruded version of this domain has already been created
1200  if (it != domain_to_extruded_domain_map.end())
1201  {
1202  // Get the associated ExtrudedDomain pointer
1203  extruded_domain_pt = it->second;
1204  }
1205  // If there is no entry then we need to create the extruded domain
1206  else
1207  {
1208  // Create an ExtrudedDomain object associated with the mesh
1209  extruded_domain_pt = new ExtrudedDomain(domain_pt, Nz, Zmin, Zmax);
1210 
1211  // Now store it in the map
1212  domain_to_extruded_domain_map[domain_pt] = extruded_domain_pt;
1213 
1214  // Add it to the mesh's list of extruded domains
1215  Extruded_domain_pt.push_back(extruded_domain_pt);
1216  }
1217 
1218  // Loop over each element slice in the z-direction
1219  for (unsigned i = 0; i < Nz; i++)
1220  {
1221  // What is the ID number of the macro element associated with
1222  // this quad element?
1223  unsigned macro_id =
1224  (i * n_macro_element +
1225  quad_el_pt->macro_elem_pt()->macro_element_number());
1226 
1227  // Get a pointer to the extruded macro element
1228  ExtrudedMacroElement* extruded_macro_element_pt =
1229  extruded_domain_pt->macro_element_pt(macro_id);
1230 
1231  // Get the index of the cube element in the extruded mesh
1232  unsigned element_index = i * n_element_in_quad_mesh + e;
1233 
1234  // Get the corresponding extruded element from the global storage
1235  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Element_pt[element_index]);
1236 
1237  // Set its macro element pointer
1238  el_pt->set_macro_elem_pt(extruded_macro_element_pt);
1239 
1240  // The number of dimensions in the extruded mesh
1241  unsigned n_dim = 3;
1242 
1243  // Loop over the coordinate directions
1244  for (unsigned i = 0; i < n_dim - 1; i++)
1245  {
1246  // The i-th local coordinate value of the "lower-left" corner of
1247  // the current element
1248  el_pt->s_macro_ll(i) = quad_el_pt->s_macro_ll(i);
1249 
1250  // The i-th local coordinate value of the "upper-right" corner of
1251  // the current element
1252  el_pt->s_macro_ur(i) = quad_el_pt->s_macro_ur(i);
1253  }
1254 
1255  // The local coordinate value of the "lower-left" corner which
1256  // corresponds to the time direction
1257  el_pt->s_macro_ll(n_dim - 1) = -1.0;
1258 
1259  // The local coordinate value of the "upper-right" corner which
1260  // corresponds to the time direction
1261  el_pt->s_macro_ur(n_dim - 1) = 1.0;
1262  } // for (unsigned i=0;i<Nz;i++)
1263  } // if (quad_el_pt->macro_element_pt()!=0)
1264  } // for (unsigned e=0;e<n_element_in_quad_mesh;e++)
1265 
1266  // If the user wishes the mesh setup time to be doc-ed
1267  if (MeshExtrusionHelpers::Mesh_extrusion_helper.doc_mesh_setup_time())
1268  {
1269  // Tell the user
1270  oomph_info << "Time taken to set up extruded macro-elements [sec]: "
1271  << TimingHelpers::timer() - t_mesh_macro_start << std::endl;
1272  }
1273 
1274  // Call the node update function to do all the relevant node moving
1275  node_update();
1276 
1277 #ifdef PARANOID
1278  // Make sure there are no duplicate nodes (i.e. two or more nodes that
1279  // occupy the same Eulerian position)
1281  {
1282  // Throw a warning
1283  OomphLibWarning("The mesh contains repeated nodes!",
1286  }
1287 
1288  // Also make sure enough nodes have been been created
1289  if (node_count != n_node_in_cube_mesh)
1290  {
1291  // Create an ostringstream object to create an error message
1292  std::ostringstream error_message_stream;
1293 
1294  // Create the error message
1295  error_message_stream << "An incorrect number of nodes have been created! "
1296  << "There are\nmeant to be " << n_node_in_cube_mesh
1297  << " nodes in the extruded mesh but\nonly "
1298  << node_count << " have actually been constructed!"
1299  << std::endl;
1300 
1301  // Throw a warning
1302  OomphLibWarning(error_message_stream.str(),
1305  }
1306 #endif
1307 
1308  // Setup lookup scheme that establishes which elements are located
1309  // on the mesh boundaries
1311 
1312  // If the user wishes the mesh setup time to be doc-ed
1313  if (MeshExtrusionHelpers::Mesh_extrusion_helper.doc_mesh_setup_time())
1314  {
1315  // Tell the user
1316  oomph_info << "Total time for extruded mesh creation/setup [sec]: "
1317  << TimingHelpers::timer() - t_mesh_setup_start << std::endl;
1318  }
1319  } // End of build_mesh
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Scalar * b
Definition: benchVecAdd.cpp:17
The simulation can be subdivided into Domain's used in parallel code.
Definition: Domain.h:42
void setup_boundary_element_info()
Definition: brick_mesh.h:195
virtual double z_spacing_function(const unsigned &z_element, const unsigned &z_node) const
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h:163
Vector< std::pair< unsigned, int > > get_element_boundary_information(QuadMeshBase *quad_mesh_pt, FiniteElement *quad_el_pt)
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.cc:45
unsigned N_node_1d
The number of nodes in each direction.
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h:192
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
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
virtual void node_update(const bool &update_all_solid_nodes=false)
Definition: mesh.cc:287
unsigned check_for_repeated_nodes(const double &epsilon=1.0e-12)
Definition: mesh.h:752
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
return int(ret)+1
class oomph::MeshExtrusionHelpers::ExtrusionHelper Mesh_extrusion_helper
@ E
Definition: quadtree.h:61
@ S
Definition: quadtree.h:62
@ W
Definition: quadtree.h:63
@ N
Definition: quadtree.h:60
double timer()
returns the time in seconds after some point in past
Definition: oomph_utilities.cc:1295
OomphInfo oomph_info
Definition: oomph_definitions.cc:319
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References b, oomph::Mesh::boundary_node_pt(), oomph::MacroElement::domain_pt(), e(), oomph::QuadTreeNames::E, oomph::Mesh::finite_element_pt(), oomph::Node::get_boundaries_pt(), i, int(), oomph::Node::is_hanging(), oomph::Node::is_on_boundary(), oomph::FiniteElement::macro_elem_pt(), oomph::MacroElement::macro_element_number(), oomph::ExtrudedDomain::macro_element_pt(), oomph::MeshExtrusionHelpers::Mesh_extrusion_helper, n, oomph::QuadTreeNames::N, oomph::Mesh::nboundary(), oomph::Mesh::nboundary_node(), oomph::Mesh::nelement(), oomph::Domain::nmacro_element(), oomph::FiniteElement::nnode(), oomph::Mesh::nnode(), oomph::FiniteElement::nnode_1d(), oomph::FiniteElement::node_pt(), oomph::Mesh::node_pt(), Global_Parameters::Nz, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::QuadTreeNames::S, oomph::QElementBase::s_macro_ll(), oomph::QElementBase::s_macro_ur(), oomph::TimingHelpers::timer(), oomph::QuadTreeNames::W, oomph::Node::x(), Global_Parameters::Zmax, and Global_Parameters::Zmin.

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

◆ get_element_boundary_information()

template<class ELEMENT >
Vector< std::pair< unsigned, int > > oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::get_element_boundary_information ( QuadMeshBase quad_mesh_pt,
FiniteElement quad_el_pt 
)

Get all the boundary information of an element using the input (quad_mesh_pt) mesh. If the element lies on a boundary then the user will be given the corresponding boundary index and the index of the face of quad_el_pt attached to the boundary. If the element does NOT lie on any boundaries, this function simply returns a vector of size zero.

47  {
48  // Allocate space for the boundary info (initialised to empty)
49  Vector<std::pair<unsigned, int>> boundary_info;
50 
51  // Find the number of boundaries in the input mesh
52  unsigned n_boundary = quad_mesh_pt->nboundary();
53 
54  // Loop over the boundaries of the mesh
55  for (unsigned b = 0; b < n_boundary; b++)
56  {
57  // How many elements are there on the b-th boundary?
58  unsigned n_boundary_element = quad_mesh_pt->nboundary_element(b);
59 
60  // Loop over the elements on the b-th boundary
61  for (unsigned e = 0; e < n_boundary_element; e++)
62  {
63  // Is the e-th element on the b-th boundary the input element?
64  if (quad_el_pt == quad_mesh_pt->boundary_element_pt(b, e))
65  {
66  // Create a pair to hold the boundary index and element face index
67  std::pair<unsigned, int> boundary_and_face_index;
68 
69  // Set the first entry (boundary index)
70  boundary_and_face_index.first = b;
71 
72  // Set the second entry (face index)
73  boundary_and_face_index.second =
74  quad_mesh_pt->face_index_at_boundary(b, e);
75 
76  // Add the boundary index to the boundary_info vector
77  boundary_info.push_back(boundary_and_face_index);
78  }
79  } // for (unsigned e=0;e<n_boundary_element;e++)
80  } // for (unsigned b=0;b<n_boundary;b++)
81 
82  // Return the result
83  return boundary_info;
84  } // End of get_element_boundary_information

References b, oomph::Mesh::boundary_element_pt(), e(), oomph::Mesh::face_index_at_boundary(), oomph::Mesh::nboundary(), and oomph::Mesh::nboundary_element().

◆ nz()

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

Access function for number of elements in z-direction (const version)

185  {
186  // Return the value of the appropriate private variable
187  return Nz;
188  } // End of nz

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

◆ z_spacing_function()

template<class ELEMENT >
virtual double oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::z_spacing_function ( const unsigned z_element,
const unsigned z_node 
) const
inlinevirtual

Return the value of the z-coordinate at the node given by the local node number, znode.

165  {
166  // Calculate the values of equal increments in nodal values
167  double z_step = (Zmax - Zmin) / ((N_node_1d - 1) * Nz);
168 
169  // Return the appropriate value
170  return (Zmin + z_step * ((N_node_1d - 1) * z_element + z_node));
171  } // End of z_spacing_function

References oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::N_node_1d, oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::Nz, oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::Zmax, and oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::Zmin.

Member Data Documentation

◆ Extruded_domain_pt

template<class ELEMENT >
Vector<ExtrudedDomain*> oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::Extruded_domain_pt
protected

Vector of pointers to extruded domain objects.

Referenced by oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::~ExtrudedCubeMeshFromQuadMesh().

◆ N_node_1d

template<class ELEMENT >
unsigned oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::N_node_1d
protected

The number of nodes in each direction.

Referenced by oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::z_spacing_function().

◆ Nz

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

◆ Zmax

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

Maximum value of z coordinate.

Referenced by oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::z_spacing_function().

◆ Zmin

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

Minimum value of z coordinate.

Referenced by oomph::ExtrudedCubeMeshFromQuadMesh< ELEMENT >::z_spacing_function().


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