Public Types | Public Member Functions | Private Member Functions | Private Attributes | List of all members
oomph::HermiteQuadMesh< ELEMENT > Class Template Reference

#include <hermite_element_quad_mesh.template.h>

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

Public Types

typedef void(* MeshSpacingFnPtr) (const Vector< double > &m_uniform_spacing, Vector< double > &m_non_uniform_spacing)
 
- 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)
 

Public Member Functions

 HermiteQuadMesh (const unsigned &nx, const unsigned &ny, TopologicallyRectangularDomain *domain, const bool &periodic_in_x=false, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
 HermiteQuadMesh (const unsigned &nx, const unsigned &ny, TopologicallyRectangularDomain *domain, const MeshSpacingFnPtr spacing_fn, const bool &periodic_in_x=false, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
 ~HermiteQuadMesh ()
 Destructor - does nothing - handled in mesh base class. More...
 
unsigned & nelement_in_dim (const unsigned &d)
 Access function for number of elements in mesh in each dimension. More...
 
- 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...
 

Private Member Functions

void macro_coordinate_position (const unsigned &node_num_x, const unsigned &node_num_y, Vector< double > &macro_element_position)
 
void set_position_of_node (const unsigned &node_num_x, const unsigned &node_num_y, Node *node_pt)
 
void set_position_of_boundary_node (const unsigned &node_num_x, const unsigned &node_num_y, BoundaryNode< Node > *node_pt)
 
void generalised_macro_element_position_of_node (const unsigned &node_num_x, const unsigned &node_num_y, DenseMatrix< double > &m_gen)
 
virtual void build_mesh (TimeStepper *time_stepper_pt)
 Generic mesh construction function to build the mesh. More...
 
virtual void setup_boundary_element_info ()
 
virtual void setup_boundary_element_info (std::ostream &outfile)
 

Private Attributes

Vector< unsigned > Nelement
 number of elements in each coordinate direction More...
 
bool Xperiodic
 
TopologicallyRectangularDomain * Domain_pt
 
MeshSpacingFnPtr Node_spacing_fn
 non uniform mesh spacing function pointer More...
 

Additional Inherited Members

- 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...
 
- 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 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
 

Detailed Description

template<class ELEMENT>
class oomph::HermiteQuadMesh< ELEMENT >

A two dimensional Hermite bicubic element quadrilateral mesh for a topologically rectangular domain. The geometry of the problem must be prescribed using the TopologicallyRectangularDomain. Non uniform node spacing can be prescribed using a function pointer.

Member Typedef Documentation

◆ MeshSpacingFnPtr

template<class ELEMENT >
typedef void(* oomph::HermiteQuadMesh< ELEMENT >::MeshSpacingFnPtr) (const Vector< double > &m_uniform_spacing, Vector< double > &m_non_uniform_spacing)

Mesh Spacing Function Pointer - an optional function pointer to prescibe the node spacing in a non-uniformly spaced mesh - takes the position of a node (in macro element coordinates) in the uniformly spaced mesh and return the position in the non-uniformly spaced mesh

Constructor & Destructor Documentation

◆ HermiteQuadMesh() [1/2]

template<class ELEMENT >
oomph::HermiteQuadMesh< ELEMENT >::HermiteQuadMesh ( const unsigned &  nx,
const unsigned &  ny,
TopologicallyRectangularDomain *  domain,
const bool &  periodic_in_x = false,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline

Mesh Constructor (for a uniformly spaced mesh). Takes the following arguments : nx : number of elements in x direction; ny : number of elements in y direction; domain : topologically rectangular domain; periodic_in_x : flag specifying if the mesh is periodic in the x direction (default = false); time_stepper_pt : pointer to the time stepper (default = no timestepper);

78  {
79  // Mesh can only be built with 2D QHermiteElements.
80  MeshChecker::assert_geometric_element<QHermiteElementBase, ELEMENT>(2);
81 
82  // set number of elements in each coordinate direction
83  Nelement.resize(2);
84  Nelement[0] = nx;
85  Nelement[1] = ny;
86 
87  // set x periodicity
88  Xperiodic = periodic_in_x;
89 
90  // set the domain pointer
91  Domain_pt = domain;
92 
93  // set the node spacing function to zero
94  Node_spacing_fn = 0;
95 
96  // builds the mesh
97  build_mesh(time_stepper_pt);
98  }
oomph::HermiteQuadMesh::Nelement
Vector< unsigned > Nelement
number of elements in each coordinate direction
Definition: hermite_element_quad_mesh.template.h:233
oomph::HermiteQuadMesh::Node_spacing_fn
MeshSpacingFnPtr Node_spacing_fn
non uniform mesh spacing function pointer
Definition: hermite_element_quad_mesh.template.h:244
oomph::HermiteQuadMesh::Xperiodic
bool Xperiodic
Definition: hermite_element_quad_mesh.template.h:237
oomph::HermiteQuadMesh::Domain_pt
TopologicallyRectangularDomain * Domain_pt
Definition: hermite_element_quad_mesh.template.h:241
oomph::HermiteQuadMesh::build_mesh
virtual void build_mesh(TimeStepper *time_stepper_pt)
Generic mesh construction function to build the mesh.
Definition: hermite_element_quad_mesh.template.cc:603
Mesh_Parameters::nx
const unsigned nx
Definition: ConstraintElementsUnitTest.cpp:30
Mesh_Parameters::ny
const unsigned ny
Definition: ConstraintElementsUnitTest.cpp:31
domain
SteadyAxisymAdvectionDiffusion problem on rectangular domain
Definition: steady_axisym_advection_diffusion.cc:151

References oomph::HermiteQuadMesh< ELEMENT >::build_mesh(), domain, oomph::HermiteQuadMesh< ELEMENT >::Domain_pt, oomph::HermiteQuadMesh< ELEMENT >::Nelement, oomph::HermiteQuadMesh< ELEMENT >::Node_spacing_fn, Mesh_Parameters::nx, Mesh_Parameters::ny, and oomph::HermiteQuadMesh< ELEMENT >::Xperiodic.

◆ HermiteQuadMesh() [2/2]

template<class ELEMENT >
oomph::HermiteQuadMesh< ELEMENT >::HermiteQuadMesh ( const unsigned &  nx,
const unsigned &  ny,
TopologicallyRectangularDomain *  domain,
const MeshSpacingFnPtr  spacing_fn,
const bool &  periodic_in_x = false,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline

Mesh Constructor (for a non-uniformly spaced mesh). Takes the following arguments : nx : number of elements in x direction; ny : number of elements in y direction; domain : topologically rectangular domain; spacing_fn : spacing function prescribing a non-uniformly spaced mesh periodic_in_x : flag specifying if the mesh is periodic in the x direction (default = false); time_stepper_pt : pointer to the time stepper (default = notimestepper);

build the mesh

120  {
121  // Mesh can only be built with 2D QHermiteElements.
122  MeshChecker::assert_geometric_element<QHermiteElementBase, ELEMENT>(2);
123 
124  // set number of elements in each coordinate direction
125  Nelement.resize(2);
126  Nelement[0] = nx;
127  Nelement[1] = ny;
128 
129  // set x periodicity
130  Xperiodic = periodic_in_x;
131 
132  // set the domain pointer
133  Domain_pt = domain;
134 
135  // set the node spacing function to zero
136  Node_spacing_fn = spacing_fn;
137 
139  build_mesh(time_stepper_pt);
140  }

References oomph::HermiteQuadMesh< ELEMENT >::build_mesh(), domain, oomph::HermiteQuadMesh< ELEMENT >::Domain_pt, oomph::HermiteQuadMesh< ELEMENT >::Nelement, oomph::HermiteQuadMesh< ELEMENT >::Node_spacing_fn, Mesh_Parameters::nx, Mesh_Parameters::ny, and oomph::HermiteQuadMesh< ELEMENT >::Xperiodic.

◆ ~HermiteQuadMesh()

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

Destructor - does nothing - handled in mesh base class.

144 {};

Member Function Documentation

◆ build_mesh()

template<class ELEMENT >
void oomph::HermiteQuadMesh< ELEMENT >::build_mesh ( TimeStepper *  time_stepper_pt)
privatevirtual

Generic mesh construction function to build the mesh.

604  {
605  // Mesh can only be built with 2D QHermiteElements.
606  MeshChecker::assert_geometric_element<QHermiteElementBase, ELEMENT>(2);
607 
608  // Set the number of boundaries
609  set_nboundary(4);
610 
611  // Allocate the store for the elements
612  Element_pt.resize(Nelement[0] * Nelement[1]);
613 
614  // Allocate the memory for the first element
615  Element_pt[0] = new ELEMENT;
616  finite_element_pt(0)->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
617 
618  // Can now allocate the store for the nodes
619  Node_pt.resize((1 + Nelement[0]) * (1 + Nelement[1]));
620 
621  // Set up geometrical data
622  unsigned long node_count = 0;
623 
624  // Now assign the topology
625  // Boundaries are numbered 0 1 2 3 from the bottom proceeding anticlockwise
626  // Pinned value are denoted by an integer value 1
627  // Thus if a node is on two boundaries, ORing the values of the
628  // boundary conditions will give the most restrictive case (pinning)
629 
630 
631  // we first create the lowest row of elements
632  // ##########################################
633  // ##########################################
634 
635 
636  // FIRST ELEMENT - lower left hand corner
637  // **************************************
638 
639 
640  // LOWER LEFT HAND NODE
641 
642  // Set the corner node
643  // Allocate memory for the node
644  Node_pt[node_count] =
645  finite_element_pt(0)->construct_boundary_node(0, time_stepper_pt);
646 
647  // Push the node back onto boundaries
648  add_boundary_node(0, Node_pt[node_count]);
649  add_boundary_node(3, Node_pt[node_count]);
650 
651  // set the position of the boundary node
652  set_position_of_boundary_node(
653  0, 0, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
654 
655  // Increment the node number
656  node_count++;
657 
658  // LOWER RIGHT HAND NODE
659 
660  // Allocate memory for the node
661  Node_pt[node_count] =
662  finite_element_pt(0)->construct_boundary_node(1, time_stepper_pt);
663 
664  // Push the node back onto boundaries
665  add_boundary_node(0, Node_pt[node_count]);
666 
667  // If we only have one column then the RHS node is on the right boundary
668  if (Nelement[0] == 1)
669  {
670  add_boundary_node(1, Node_pt[node_count]);
671  }
672 
673  // set the position of the node
674  set_position_of_boundary_node(
675  1, 0, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
676 
677  // Increment the node number
678  node_count++;
679 
680  // UPPER LEFT HAND NODE
681 
682  // Allocate memory for the node
683  Node_pt[node_count] =
684  finite_element_pt(0)->construct_boundary_node(2, time_stepper_pt);
685 
686  // Push the node back onto boundaries
687  add_boundary_node(3, Node_pt[node_count]);
688 
689  // If we only have one row, then the top node is on the top boundary
690  if (Nelement[1] == 1)
691  {
692  add_boundary_node(2, Node_pt[node_count]);
693  }
694 
695  // set the position of the boundary node
696  set_position_of_boundary_node(
697  0, 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
698 
699  // Increment the node number
700  node_count++;
701 
702  // UPPER RIGHT NODE
703 
704  // Allocate the memory for the node
705  Node_pt[node_count] =
706  finite_element_pt(0)->construct_node(3, time_stepper_pt);
707 
708  // If we only have one column then the RHS node is on the right boundary
709  if (Nelement[0] == 1)
710  {
711  add_boundary_node(1, Node_pt[node_count]);
712  }
713  // If we only have one row, then the top node is on the top boundary
714  if (Nelement[1] == 1)
715  {
716  add_boundary_node(2, Node_pt[node_count]);
717  }
718 
719  // if the node is a boundary node
720  if (Nelement[0] == 1 || Nelement[1] == 1)
721  {
722  // set the position of the boundary node
723  set_position_of_boundary_node(
724  1, 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
725  }
726  else
727  {
728  // set the position of the node
729  set_position_of_node(1, 1, Node_pt[node_count]);
730  }
731 
732  // Increment the node number
733  node_count++;
734 
735  // END OF FIRST ELEMENT
736 
737 
738  // CENTRE OF FIRST ROW OF ELEMENTS
739  // *******************************
740 
741 
742  // Now loop over the first row of elements, apart from final element
743  for (unsigned j = 1; j < (Nelement[0] - 1); j++)
744  {
745  // Allocate memory for new element
746  Element_pt[j] = new ELEMENT;
747  finite_element_pt(j)->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
748 
749  // LOWER LEFT NODE
750 
751  // lower left hand node column of nodes is same as lower right hand node
752  // in neighbouring element
753  finite_element_pt(j)->node_pt(0) = finite_element_pt(j - 1)->node_pt(1);
754 
755  // LOWER RIGHT NODE
756 
757  // Allocate memory for the nodes
758  Node_pt[node_count] =
759  finite_element_pt(j)->construct_boundary_node(1, time_stepper_pt);
760 
761  // Push the node back onto boundaries
762  add_boundary_node(0, Node_pt[node_count]);
763 
764  // set the position of the boundary node
765  set_position_of_boundary_node(
766  j + 1, 0, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
767 
768  // Increment the node number
769  node_count++;
770 
771  // UPPER LEFT NODE
772 
773  // First column of nodes is same as neighbouring element
774  finite_element_pt(j)->node_pt(2) = finite_element_pt(j - 1)->node_pt(3);
775 
776  // UPPER RIGHT NODE
777 
778  // Allocate memory for the nodes
779  Node_pt[node_count] =
780  finite_element_pt(j)->construct_node(3, time_stepper_pt);
781 
782  // If we only have one row, then the top node is on the top boundary
783  if (Nelement[0] == 1)
784  {
785  add_boundary_node(2, Node_pt[node_count]);
786  }
787 
788  // set the position of the (boundary) node
789  if (Nelement[0] == 1)
790  set_position_of_boundary_node(
791  j + 1, 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
792  else
793  set_position_of_node(j + 1, 1, Node_pt[node_count]);
794 
795  // Increment the node number
796  node_count++;
797  }
798 
799  // FINAL ELEMENT IN FIRST ROW (lower right corner element)
800  // **************************
801 
802 
803  // Only allocate if there is more than one element in the row
804  if (Nelement[0] > 1)
805  {
806  // Allocate memory for element
807  Element_pt[Nelement[0] - 1] = new ELEMENT;
808  finite_element_pt(Nelement[0] - 1)
809  ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
810 
811  // LOWER LEFT NODE
812 
813  // First column of nodes is same as neighbouring element
814  finite_element_pt(Nelement[0] - 1)->node_pt(0) =
815  finite_element_pt(Nelement[0] - 2)->node_pt(1);
816 
817  // LOWER RIGHT NODE
818 
819  // If we have an Xperiodic mesh then the final node is the first node in
820  // the row
821  if (Xperiodic == true)
822  {
823  // Note that this is periodic in the x-direction
824  finite_element_pt(Nelement[0] - 1)->node_pt(1) =
825  finite_element_pt(0)->node_pt(0);
826  }
827  // Otherwise it's a new final node
828  else
829  {
830  // Allocate memory for the node
831  Node_pt[node_count] = finite_element_pt(Nelement[0] - 1)
832  ->construct_boundary_node(1, time_stepper_pt);
833  }
834 
835  // Push the node back onto boundaries
836  add_boundary_node(0, Node_pt[node_count]);
837  add_boundary_node(1, Node_pt[node_count]);
838 
839  // set the position of the boundary node
840  set_position_of_boundary_node(
841  Nelement[0], 0, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
842 
843  // if not periodic mesh
844  if (Xperiodic == false)
845  {
846  node_count++;
847  }
848 
849  // UPPER LEFT NODE
850 
851  // same as upper right node in element to left
852  finite_element_pt(Nelement[0] - 1)->node_pt(2) =
853  finite_element_pt(Nelement[0] - 2)->node_pt(3);
854 
855  // If we only have one row, then the top node is on the top boundary
856  if (Nelement[1] == 1)
857  {
858  add_boundary_node(2, Node_pt[node_count]);
859  }
860 
861  // set the position of the (boundary) node
862  if (Nelement[1] == 1)
863  set_position_of_boundary_node(
864  Nelement[0] - 1,
865  1,
866  dynamic_cast<BoundaryNode<Node>*>(
867  finite_element_pt(Nelement[0] - 2)->node_pt(3)));
868 
869  // UPPER RIGHT NODE
870 
871  // If we have an Xperiodic mesh then the nodes in the final column are
872  // those in the first column
873  if (Xperiodic == true)
874  {
875  // Note that these are periodic in the x-direction
876  finite_element_pt(Nelement[0] - 1)->node_pt(3) =
877  finite_element_pt(0)->node_pt(2);
878  }
879  // Otherwise we need new nodes for the final column
880  else
881  {
882  // Allocate memory for node
883  Node_pt[node_count] = finite_element_pt(Nelement[0] - 1)
884  ->construct_boundary_node(3, time_stepper_pt);
885  }
886 
887  // If we only have one row, then the top node is on the top boundary
888  if (Nelement[1] == 1)
889  {
890  add_boundary_node(2, Node_pt[node_count]);
891  }
892 
893  // boundary 1 node
894  add_boundary_node(1, Node_pt[node_count]);
895 
896  // set the position of the boundary node
897  set_position_of_boundary_node(
898  Nelement[0], 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
899 
900  // Increment the node number
901  if (Xperiodic == false)
902  {
903  // Increment the node number
904  node_count++;
905  }
906  }
907 
908  // now create all remaining central rows
909  // #####################################
910  // #####################################
911 
912 
913  // Loop over remaining element rows in the mesh
914  for (unsigned i = 1; i < (Nelement[1] - 1); i++)
915  {
916  // set the first element in the row
917  // ********************************
918 
919 
920  // Allocate memory for element
921  Element_pt[Nelement[0] * i] = new ELEMENT;
922  finite_element_pt(Nelement[0] * i)
923  ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
924 
925  // The first row of nodes is copied from the element below
926  for (unsigned l = 0; l < 2; l++)
927  {
928  finite_element_pt(Nelement[0] * i)->node_pt(l) =
929  finite_element_pt(Nelement[0] * (i - 1))->node_pt(2 + l);
930  }
931 
932  // UPPER LEFT HAND NODE
933 
934  // Allocate memory for node
935  Node_pt[node_count] = finite_element_pt(Nelement[0] * i)
936  ->construct_boundary_node(2, time_stepper_pt);
937 
938  // Push the node back onto boundaries
939  add_boundary_node(3, Node_pt[node_count]);
940 
941  // set the position of the boundary node
942  set_position_of_boundary_node(
943  0, i + 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
944 
945  // Increment the node number
946  node_count++;
947 
948  // UPPER RIGHT HAND NODE
949 
950  // If we only have one column, the node could be on the boundary
951  if (Nelement[0] == 1)
952  {
953  Node_pt[node_count] = finite_element_pt(Nelement[0] * i)
954  ->construct_boundary_node(3, time_stepper_pt);
955  }
956  else
957  {
958  Node_pt[node_count] = finite_element_pt(Nelement[0] * i)
959  ->construct_node(3, time_stepper_pt);
960  }
961 
962  // If we only have one column then the RHS node is on the
963  // right boundary
964  if (Nelement[0] == 1)
965  {
966  add_boundary_node(1, Node_pt[node_count]);
967  }
968 
969  // set the position of the (boundary) node
970  if (Nelement[0] == 1)
971  set_position_of_boundary_node(
972  1, i + 1, dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
973  else
974  set_position_of_node(1, i + 1, Node_pt[node_count]);
975 
976  // Increment the node number
977  node_count++;
978 
979 
980  // loop over central elements in row
981  // *********************************
982 
983  // Now loop over the rest of the elements in the row, apart from the last
984  for (unsigned j = 1; j < (Nelement[0] - 1); j++)
985  {
986  // Allocate memory for new element
987  Element_pt[Nelement[0] * i + j] = new ELEMENT;
988  finite_element_pt(Nelement[0] * i + j)
989  ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
990 
991  // LOWER LEFT AND RIGHT NODES
992 
993  // The first row is copied from the lower element
994  for (unsigned l = 0; l < 2; l++)
995  {
996  finite_element_pt(Nelement[0] * i + j)->node_pt(l) =
997  finite_element_pt(Nelement[0] * (i - 1) + j)->node_pt(2 + l);
998  }
999 
1000  // UPPER LEFT NODE
1001 
1002  // First column is same as neighbouring element
1003  finite_element_pt(Nelement[0] * i + j)->node_pt(2) =
1004  finite_element_pt(Nelement[0] * i + (j - 1))->node_pt(3);
1005 
1006  // UPPER RIGHT NODE
1007 
1008  // Allocate memory for the nodes
1009  Node_pt[node_count] = finite_element_pt(Nelement[0] * i + j)
1010  ->construct_node(3, time_stepper_pt);
1011 
1012  // set position of node
1013  set_position_of_node(j + 1, i + 1, Node_pt[node_count]);
1014 
1015  // Increment the node number
1016  node_count++;
1017  }
1018 
1019 
1020  // final element in row
1021  // ********************
1022 
1023 
1024  // Only if there is more than one column
1025  if (Nelement[0] > 1)
1026  {
1027  // Allocate memory for element
1028  Element_pt[Nelement[0] * i + Nelement[0] - 1] = new ELEMENT;
1029  finite_element_pt(Nelement[0] * i + Nelement[0] - 1)
1030  ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1031 
1032  // LOWER LEFT AND RIGHT NODES
1033 
1034  // The first row is copied from the lower element
1035  for (unsigned l = 0; l < 2; l++)
1036  {
1037  finite_element_pt(Nelement[0] * i + Nelement[0] - 1)->node_pt(l) =
1038  finite_element_pt(Nelement[0] * (i - 1) + Nelement[0] - 1)
1039  ->node_pt(2 + l);
1040  }
1041 
1042  // UPPER LEFT NODE
1043 
1044  // First node is same as neighbouring element
1045  finite_element_pt(Nelement[0] * i + Nelement[0] - 1)->node_pt(2) =
1046  finite_element_pt(Nelement[0] * i + Nelement[0] - 2)->node_pt(3);
1047 
1048  // UPPER RIGHT NODE
1049 
1050  // If we have an Xperiodic mesh then this node is the same
1051  // as the first node
1052  if (Xperiodic == true)
1053  {
1054  // Allocate memory for node, periodic in x-direction
1055  finite_element_pt(Nelement[0] * i + Nelement[0] - 1)->node_pt(3) =
1056  finite_element_pt(Nelement[0] * i)->node_pt(2);
1057  }
1058  // Otherwise allocate memory for a new node
1059  else
1060  {
1061  Node_pt[node_count] =
1062  finite_element_pt(Nelement[0] * i + Nelement[0] - 1)
1063  ->construct_boundary_node(3, time_stepper_pt);
1064  }
1065 
1066  // Push the node back onto boundaries
1067  add_boundary_node(1, Node_pt[node_count]);
1068 
1069  // set position of boundary node
1070  set_position_of_boundary_node(
1071  Nelement[0],
1072  i + 1,
1073  dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1074 
1075  // Increment the node number
1076  node_count++;
1077  }
1078  }
1079 
1080 
1081  // final element row
1082  // #################
1083  // #################
1084 
1085 
1086  // only if there is more than one row
1087  if (Nelement[1] > 1)
1088  {
1089  // first element in upper row (upper left corner)
1090  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1091 
1092  // Allocate memory for element
1093  Element_pt[Nelement[0] * (Nelement[1] - 1)] = new ELEMENT;
1094  finite_element_pt(Nelement[0] * (Nelement[1] - 1))
1095  ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1096 
1097  // LOWER LEFT AND LOWER RIGHT NODES
1098 
1099  // The first row of nodes is copied from the element below
1100  for (unsigned l = 0; l < 2; l++)
1101  {
1102  finite_element_pt(Nelement[0] * (Nelement[1] - 1))->node_pt(l) =
1103  finite_element_pt(Nelement[0] * (Nelement[1] - 2))->node_pt(2 + l);
1104  }
1105 
1106  // UPPER LEFT NODE
1107 
1108  // Allocate memory for node
1109  Node_pt[node_count] = finite_element_pt(Nelement[0] * (Nelement[1] - 1))
1110  ->construct_boundary_node(2, time_stepper_pt);
1111 
1112  // Push the node back onto boundaries
1113  add_boundary_node(2, Node_pt[node_count]);
1114  add_boundary_node(3, Node_pt[node_count]);
1115 
1116  // set position of boundary node
1117  set_position_of_boundary_node(
1118  0, Nelement[1], dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1119 
1120  // Increment the node number
1121  node_count++;
1122 
1123  // UPPER RIGHT NODE
1124 
1125  // Allocate memory for the node
1126  Node_pt[node_count] = finite_element_pt(Nelement[0] * (Nelement[1] - 1))
1127  ->construct_boundary_node(3, time_stepper_pt);
1128 
1129  // Push the node back onto boundaries
1130  add_boundary_node(2, Node_pt[node_count]);
1131 
1132  // set position of boundary node
1133  set_position_of_boundary_node(
1134  1, Nelement[1], dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1135 
1136  // Increment the node number
1137  node_count++;
1138 
1139 
1140  // loop over central elements of upper element row
1141  // ***********************************************
1142 
1143 
1144  // Now loop over the rest of the elements in the row, apart from the last
1145  for (unsigned j = 1; j < (Nelement[0] - 1); j++)
1146  {
1147  // Allocate memory for element
1148  Element_pt[Nelement[0] * (Nelement[1] - 1) + j] = new ELEMENT;
1149  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + j)
1150  ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1151 
1152  // LOWER LEFT AND LOWER RIGHT NODES
1153 
1154  // The first row is copied from the lower element
1155  for (unsigned l = 0; l < 2; l++)
1156  {
1157  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + j)->node_pt(l) =
1158  finite_element_pt(Nelement[0] * (Nelement[1] - 2) + j)
1159  ->node_pt(2 + l);
1160  }
1161 
1162  // UPPER LEFT NODE
1163 
1164  // First column is same as neighbouring element
1165  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + j)->node_pt(2) =
1166  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + (j - 1))
1167  ->node_pt(3);
1168 
1169  // UPPER RIGHT NODE
1170 
1171  // Allocate memory for node
1172  Node_pt[node_count] =
1173  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + j)
1174  ->construct_boundary_node(3, time_stepper_pt);
1175 
1176  // Push the node back onto boundaries
1177  add_boundary_node(2, Node_pt[node_count]);
1178 
1179  // set position of boundary node
1180  set_position_of_boundary_node(
1181  j + 1,
1182  Nelement[1],
1183  dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1184 
1185  // Increment the node number
1186  node_count++;
1187 
1188  } // End of loop over central elements in row
1189 
1190 
1191  // final element in row (upper right corner)
1192  // *****************************************
1193 
1194  // Only if there is more than one column
1195  if (Nelement[0] > 1)
1196  {
1197  // Allocate memory for element
1198  Element_pt[Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1] =
1199  new ELEMENT;
1200  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1201  ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1202 
1203  // LOWER LEFT AND LOWER RIGHT NODES
1204 
1205  // The first row is copied from the lower element
1206  for (unsigned l = 0; l < 2; l++)
1207  {
1208  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1209  ->node_pt(l) =
1210  finite_element_pt(Nelement[0] * (Nelement[1] - 2) + Nelement[0] - 1)
1211  ->node_pt(2 + l);
1212  }
1213 
1214  // UPPER LEFT NODE
1215 
1216  // First column is same as neighbouring element
1217  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1218  ->node_pt(2) =
1219  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 2)
1220  ->node_pt(3);
1221 
1222  // UPPER RIGHT NODE
1223 
1224  // If we have an Xperiodic mesh, the node must be copied from
1225  // the first column
1226  if (Xperiodic == true)
1227  {
1228  // Allocate memory for node, periodic in x-direction
1229  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1230  ->node_pt(3) =
1231  finite_element_pt(Nelement[0] * (Nelement[1] - 1))->node_pt(2);
1232  }
1233  // Otherwise, allocate new memory for node
1234  else
1235  {
1236  Node_pt[node_count] =
1237  finite_element_pt(Nelement[0] * (Nelement[1] - 1) + Nelement[0] - 1)
1238  ->construct_boundary_node(3, time_stepper_pt);
1239  }
1240 
1241  // Push the node back onto boundaries
1242  add_boundary_node(1, Node_pt[node_count]);
1243  add_boundary_node(2, Node_pt[node_count]);
1244 
1245  // set position of boundary node
1246  set_position_of_boundary_node(
1247  Nelement[0],
1248  Nelement[1],
1249  dynamic_cast<BoundaryNode<Node>*>(Node_pt[node_count]));
1250 
1251  // Increment the node number
1252  node_count++;
1253  }
1254  }
1255 
1256  // Setup boundary element lookup schemes
1257  setup_boundary_element_info();
1258  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
oomph::Domain::macro_element_pt
MacroElement * macro_element_pt(const unsigned &i)
Access to i-th macro element.
Definition: domain.h:116
oomph::FiniteElement::set_macro_elem_pt
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Definition: elements.h:1872
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::HermiteQuadMesh::setup_boundary_element_info
virtual void setup_boundary_element_info()
Definition: hermite_element_quad_mesh.template.h:215
oomph::HermiteQuadMesh::set_position_of_boundary_node
void set_position_of_boundary_node(const unsigned &node_num_x, const unsigned &node_num_y, BoundaryNode< Node > *node_pt)
Definition: hermite_element_quad_mesh.template.cc:116
oomph::HermiteQuadMesh::set_position_of_node
void set_position_of_node(const unsigned &node_num_x, const unsigned &node_num_y, Node *node_pt)
Definition: hermite_element_quad_mesh.template.cc:48
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
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References i, and j.

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

◆ generalised_macro_element_position_of_node()

template<class ELEMENT >
void oomph::HermiteQuadMesh< ELEMENT >::generalised_macro_element_position_of_node ( const unsigned &  node_num_x,
const unsigned &  node_num_y,
DenseMatrix< double > &  m_gen 
)
private

computes the generalised position of the node at position (node_num_x, node_num_y) in the macro element coordinate scheme. index 0 of m_gen : 0 - m_i 1 - dm_i/ds_0 2 - dm_i/ds_1 3 - d2m_i/ds_0ds_1 (where i is index 1 of m_gen)

204  {
205  // obtain position of node
206  Vector<double> s_macro_N(2);
207  macro_coordinate_position(node_num_x, node_num_y, s_macro_N);
208 
209  // set position of node in macro element coordinates
210  m_gen(0, 0) = s_macro_N[0];
211  m_gen(0, 1) = s_macro_N[1];
212 
213  // if on left hand edge
214  if (node_num_x == 0)
215  {
216  // first dm0ds0 and dm1ds0
217 
218  // get position of right node
219  Vector<double> s_macro_R(2);
220  macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
221 
222  // compute dm0ds0
223  m_gen(1, 0) = 0.5 * (s_macro_R[0] - s_macro_N[0]);
224 
225  // compute dm1ds0
226  m_gen(1, 1) = 0.5 * (s_macro_R[1] - s_macro_N[1]);
227 
228  // now d2m0/ds0ds1 and d2m1/ds0ds1
229 
230  // if lower left hand corner
231  if (node_num_y == 0)
232  {
233  // get position of upper right node
234  Vector<double> s_macro_UR(2);
235  macro_coordinate_position(node_num_x + 1, node_num_y + 1, s_macro_UR);
236 
237  // get position of upper node
238  Vector<double> s_macro_U(2);
239  macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
240 
241  // compute d2m0/ds0ds1
242  m_gen(3, 0) =
243  0.5 * (0.5 * (s_macro_UR[0] - s_macro_U[0]) - m_gen(1, 0));
244 
245  // compute d2m1/ds0ds1
246  m_gen(3, 1) =
247  0.5 * (0.5 * (s_macro_UR[1] - s_macro_U[1]) - m_gen(1, 1));
248  }
249 
250  // else if upper left hand corner
251  else if (node_num_y == Nelement[1])
252  {
253  // get position of lower right node
254  Vector<double> s_macro_DR(2);
255  macro_coordinate_position(node_num_x + 1, node_num_y - 1, s_macro_DR);
256 
257  // get position of lower node
258  Vector<double> s_macro_D(2);
259  macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
260 
261  // compute d2m0/ds0ds1
262  m_gen(3, 0) =
263  0.5 * (m_gen(1, 0) - 0.5 * (s_macro_DR[0] - s_macro_D[0]));
264 
265  // compute d2m0/ds0ds1
266  m_gen(3, 1) =
267  0.5 * (m_gen(1, 1) - 0.5 * (s_macro_DR[1] - s_macro_D[1]));
268  }
269 
270  // else left hand edge (not corner)
271  else
272  {
273  // get position of lower right node
274  Vector<double> s_macro_DR(2);
275  macro_coordinate_position(node_num_x + 1, node_num_y - 1, s_macro_DR);
276 
277  // get position of lower node
278  Vector<double> s_macro_D(2);
279  macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
280 
281  // get position of upper right node
282  Vector<double> s_macro_UR(2);
283  macro_coordinate_position(node_num_x + 1, node_num_y + 1, s_macro_UR);
284 
285  // get position of upper node
286  Vector<double> s_macro_U(2);
287  macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
288 
289  // compute d2m0/ds0ds1
290  m_gen(3, 0) =
291  0.5 * std::min(m_gen(1, 0) - 0.5 * (s_macro_DR[0] - s_macro_D[0]),
292  0.5 * (s_macro_UR[0] - s_macro_U[0]) - m_gen(1, 0));
293 
294  // compute d2m1/ds0ds1
295  m_gen(3, 1) =
296  0.5 * std::min(m_gen(1, 1) - 0.5 * (s_macro_DR[1] - s_macro_D[1]),
297  0.5 * (s_macro_UR[1] - s_macro_U[1]) - m_gen(1, 1));
298  }
299  }
300 
301  // if on right hand edge
302  else if (node_num_x == Nelement[0])
303  {
304  // first dm0ds0 and dm1ds0
305 
306  // get position of left node
307  Vector<double> s_macro_L(2);
308  macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
309 
310  // compute dm0ds0
311  m_gen(1, 0) = 0.5 * (s_macro_N[0] - s_macro_L[0]);
312 
313  // compute dm1ds0
314  m_gen(1, 1) = 0.5 * (s_macro_N[1] - s_macro_L[1]);
315 
316  // now d2m0/ds0ds1 and d2m1/ds0ds1
317 
318  // if lower right hand corner
319  if (node_num_y == 0)
320  {
321  // get position of upper left node
322  Vector<double> s_macro_UL(2);
323  macro_coordinate_position(node_num_x - 1, node_num_y + 1, s_macro_UL);
324 
325  // get position of upper node
326  Vector<double> s_macro_U(2);
327  macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
328 
329  // compute d2m0/ds0ds1
330  m_gen(3, 0) =
331  0.5 * (0.5 * (s_macro_U[0] - s_macro_UL[0]) - m_gen(1, 0));
332 
333  // compute d2m1/ds0ds1
334  m_gen(3, 1) =
335  0.5 * (0.5 * (s_macro_U[1] - s_macro_UL[1]) - m_gen(1, 1));
336  }
337 
338  // else if upper right hand corner
339  else if (node_num_y == Nelement[1])
340  {
341  // get position of lower left node
342  Vector<double> s_macro_DL(2);
343  macro_coordinate_position(node_num_x - 1, node_num_y - 1, s_macro_DL);
344 
345  // get position of lower node
346  Vector<double> s_macro_D(2);
347  macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
348 
349  // compute d2m0/ds0ds1
350  m_gen(3, 0) =
351  0.5 * (m_gen(1, 0) - 0.5 * (s_macro_D[0] - s_macro_DL[0]));
352 
353  // compute d2m0/ds0ds1
354  m_gen(3, 1) =
355  0.5 * (m_gen(1, 1) - 0.5 * (s_macro_D[1] - s_macro_DL[1]));
356  }
357 
358  // else left hand edge (not corner)
359  else
360  {
361  // get position of lower left node
362  Vector<double> s_macro_DL(2);
363  macro_coordinate_position(node_num_x - 1, node_num_y - 1, s_macro_DL);
364 
365  // get position of lower node
366  Vector<double> s_macro_D(2);
367  macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
368 
369  // get position of upper left node
370  Vector<double> s_macro_UL(2);
371  macro_coordinate_position(node_num_x - 1, node_num_y + 1, s_macro_UL);
372 
373  // get position of upper node
374  Vector<double> s_macro_U(2);
375  macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
376 
377  // compute d2m0/ds0ds1
378  m_gen(3, 0) =
379  0.5 * std::min(m_gen(1, 0) - 0.5 * (s_macro_D[0] - s_macro_DL[0]),
380  0.5 * (s_macro_U[0] - s_macro_UL[0]) - m_gen(1, 0));
381 
382  // compute d2m0/ds0ds1
383  m_gen(3, 1) =
384  0.5 * std::min(m_gen(1, 1) - 0.5 * (s_macro_D[1] - s_macro_DL[1]),
385  0.5 * (s_macro_U[1] - s_macro_UL[1]) - m_gen(1, 1));
386  }
387  }
388  //
389  else
390  {
391  // get position of left node
392  Vector<double> s_macro_L(2);
393  macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
394 
395  // get position of right node
396  Vector<double> s_macro_R(2);
397  macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
398 
399  // compute dm0/ds0
400  m_gen(1, 0) = 0.5 * std::min(s_macro_N[0] - s_macro_L[0],
401  s_macro_R[0] - s_macro_N[0]);
402 
403  // compute dm1/ds0
404  Vector<double> node_space(2);
405  node_space[0] = s_macro_N[1] - s_macro_L[1];
406  node_space[1] = s_macro_R[1] - s_macro_N[1];
407  // set nodal dof
408  if (node_space[0] > 0 && node_space[1] > 0)
409  m_gen(1, 1) = 0.5 * std::min(node_space[0], node_space[1]);
410  else if (node_space[0] < 0 && node_space[1] < 0)
411  m_gen(1, 1) = 0.5 * std::max(node_space[0], node_space[1]);
412  else
413  m_gen(1, 1) = 0;
414  }
415 
416  // if node in lower row
417  if (node_num_y == 0)
418  {
419  // now dm0ds1 and dm1ds1
420 
421  // get position of upper node
422  Vector<double> s_macro_U(2);
423  macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
424 
425  // compute dm0ds1
426  m_gen(2, 0) = 0.5 * (s_macro_U[0] - s_macro_N[0]);
427 
428  // compute dm1ds1
429  m_gen(2, 1) = 0.5 * (s_macro_U[1] - s_macro_N[1]);
430 
431  // and if not corner node
432  if (node_num_x > 0 && node_num_x < Nelement[0])
433  {
434  // now dm0/ds0ds1 and dm1/ds0ds1
435 
436  // get position of upper left node
437  Vector<double> s_macro_UL(2);
438  macro_coordinate_position(node_num_x - 1, node_num_y + 1, s_macro_UL);
439 
440  // get position of left node
441  Vector<double> s_macro_L(2);
442  macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
443 
444  // get position of right node
445  Vector<double> s_macro_R(2);
446  macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
447 
448  // get position of upper right node
449  Vector<double> s_macro_UR(2);
450  macro_coordinate_position(node_num_x + 1, node_num_y + 1, s_macro_UR);
451 
452  // set dm0/ds0ds1
453  m_gen(3, 0) =
454  0.5 * std::min(m_gen(2, 0) - 0.5 * (s_macro_UL[0] - s_macro_L[0]),
455  0.5 * (s_macro_UR[0] - s_macro_R[0]) - m_gen(2, 0));
456 
457  // set dm1/ds0ds1
458  m_gen(3, 1) =
459  0.5 * std::min(m_gen(2, 1) - 0.5 * (s_macro_UL[1] - s_macro_L[1]),
460  0.5 * (s_macro_UR[1] - s_macro_R[1]) - m_gen(2, 1));
461  }
462  }
463 
464  // else if node in upper row
465  else if (node_num_y == Nelement[1])
466  {
467  // now dm0ds1 and dm1ds1
468 
469  // get position of lower node
470  Vector<double> s_macro_D(2);
471  macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
472 
473  // compute dm0ds1
474  m_gen(2, 0) = 0.5 * (s_macro_N[0] - s_macro_D[0]);
475 
476  // compute dm1ds1
477  m_gen(2, 1) = 0.5 * (s_macro_N[1] - s_macro_D[1]);
478 
479  // and if not corner node
480  if (node_num_x > 0 && node_num_x < Nelement[0])
481  {
482  // now dm0/ds0ds1 and dm1/ds0ds1
483 
484  // get position of lower left node
485  Vector<double> s_macro_DL(2);
486  macro_coordinate_position(node_num_x - 1, node_num_y - 1, s_macro_DL);
487 
488  // get position of left node
489  Vector<double> s_macro_L(2);
490  macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
491 
492  // get position of right node
493  Vector<double> s_macro_R(2);
494  macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
495 
496  // get position of lower right node
497  Vector<double> s_macro_DR(2);
498  macro_coordinate_position(node_num_x + 1, node_num_y - 1, s_macro_DR);
499 
500  // set dm0/ds0ds1
501  m_gen(3, 0) =
502  0.5 * std::min(m_gen(2, 0) - 0.5 * (s_macro_L[0] - s_macro_DL[0]),
503  0.5 * (s_macro_R[0] - s_macro_DR[0]) - m_gen(2, 0));
504 
505  // set dm1/ds0ds1
506  m_gen(3, 1) =
507  0.5 * std::min(m_gen(2, 1) - 0.5 * (s_macro_L[1] - s_macro_DL[1]),
508  0.5 * (s_macro_R[1] - s_macro_DR[1]) - m_gen(2, 1));
509  }
510  }
511  else
512  {
513  // get position of upper node
514  Vector<double> s_macro_U(2);
515  macro_coordinate_position(node_num_x, node_num_y + 1, s_macro_U);
516 
517  // get position of lower node
518  Vector<double> s_macro_D(2);
519  macro_coordinate_position(node_num_x, node_num_y - 1, s_macro_D);
520 
521  // compute dm0/ds1
522  Vector<double> node_space(2);
523  node_space[0] = s_macro_N[0] - s_macro_D[0];
524  node_space[1] = s_macro_U[0] - s_macro_N[0];
525  // set nodal coordinate
526  if (node_space[0] > 0 && node_space[1] > 0)
527  m_gen(2, 0) = 0.5 * std::min(node_space[0], node_space[1]);
528  else if (node_space[0] < 0 && node_space[1] < 0)
529  m_gen(2, 0) = 0.5 * std::max(node_space[0], node_space[1]);
530  else
531  m_gen(2, 0) = 0;
532 
533  // compute dm1/ds1
534  m_gen(2, 1) = 0.5 * std::min(s_macro_N[1] - s_macro_D[1],
535  s_macro_U[1] - s_macro_N[1]);
536 
537  // for interior nodes
538  if (node_num_x > 0 && node_num_x < Nelement[0])
539  {
540  // get position of left node
541  Vector<double> s_macro_L(2);
542  macro_coordinate_position(node_num_x - 1, node_num_y, s_macro_L);
543 
544  // get position of upper left node
545  Vector<double> s_macro_UL(2);
546  macro_coordinate_position(node_num_x - 1, node_num_y + 1, s_macro_UL);
547 
548  // get position of upper right node
549  Vector<double> s_macro_UR(2);
550  macro_coordinate_position(node_num_x + 1, node_num_y + 1, s_macro_UR);
551 
552  // get position of right node
553  Vector<double> s_macro_R(2);
554  macro_coordinate_position(node_num_x + 1, node_num_y, s_macro_R);
555 
556  // get position of lower right node
557  Vector<double> s_macro_DR(2);
558  macro_coordinate_position(node_num_x + 1, node_num_y - 1, s_macro_DR);
559 
560  // get position of lower left node
561  Vector<double> s_macro_DL(2);
562  macro_coordinate_position(node_num_x - 1, node_num_y - 1, s_macro_DL);
563 
564  Vector<double> node_space(2);
565  // comute dm0/ds0ds1 wrt to node above and below
566  node_space[0] =
567  m_gen(1, 0) - 0.5 * std::min(s_macro_D[0] - s_macro_DL[0],
568  s_macro_DR[0] - s_macro_D[0]);
569  node_space[1] = 0.5 * std::min(s_macro_U[0] - s_macro_UL[0],
570  s_macro_UR[0] - s_macro_U[0]) -
571  m_gen(1, 0);
572  // set nodal dof
573  if (node_space[0] > 0 && node_space[1] > 0)
574  m_gen(3, 0) = 0.5 * std::min(node_space[0], node_space[1]);
575  else if (node_space[0] < 0 && node_space[1] < 0)
576  m_gen(3, 0) = 0.5 * std::max(node_space[0], node_space[1]);
577  else
578  m_gen(3, 0) = 0;
579 
580  // comute dm1/ds0ds1 wrt node left and right
581  node_space[0] =
582  m_gen(2, 1) - 0.5 * std::min(s_macro_L[0] - s_macro_DL[0],
583  s_macro_UL[0] - s_macro_L[0]);
584  node_space[1] = 0.5 * std::min(s_macro_R[0] - s_macro_DR[0],
585  s_macro_UR[0] - s_macro_R[0]) -
586  m_gen(2, 1);
587  // set nodal dof
588  if (node_space[0] > 0 && node_space[1] > 0)
589  m_gen(3, 1) = 0.5 * std::min(node_space[0], node_space[1]);
590  else if (node_space[0] < 0 && node_space[1] < 0)
591  m_gen(3, 1) = 0.5 * std::max(node_space[0], node_space[1]);
592  else
593  m_gen(3, 1) = 0;
594  }
595  }
596  }
oomph::HermiteQuadMesh::macro_coordinate_position
void macro_coordinate_position(const unsigned &node_num_x, const unsigned &node_num_y, Vector< double > &macro_element_position)
Definition: hermite_element_quad_mesh.template.h:157
min
#define min(a, b)
Definition: datatypes.h:22
max
#define max(a, b)
Definition: datatypes.h:23

References max, and min.

◆ macro_coordinate_position()

template<class ELEMENT >
void oomph::HermiteQuadMesh< ELEMENT >::macro_coordinate_position ( const unsigned &  node_num_x,
const unsigned &  node_num_y,
Vector< double > &  macro_element_position 
)
inlineprivate

returns the macro element position of the node that is the x-th node along from the LHS and the y-th node up from the lower edge

160  {
161  // compute macro element position in uniformly spaced mesh
162  macro_element_position[0] = 2 * node_num_x / double(Nelement[0]) - 1;
163  macro_element_position[1] = 2 * node_num_y / double(Nelement[1]) - 1;
164 
165  // if a non unform spacing function is provided
166  if (Node_spacing_fn != 0)
167  {
168  Vector<double> temp(macro_element_position);
169  (*Node_spacing_fn)(temp, macro_element_position);
170  }
171  }

References oomph::HermiteQuadMesh< ELEMENT >::Nelement, and oomph::HermiteQuadMesh< ELEMENT >::Node_spacing_fn.

◆ nelement_in_dim()

template<class ELEMENT >
unsigned& oomph::HermiteQuadMesh< ELEMENT >::nelement_in_dim ( const unsigned &  d)
inline

Access function for number of elements in mesh in each dimension.

149  {
150  return Nelement[d];
151  }

References oomph::HermiteQuadMesh< ELEMENT >::Nelement.

◆ set_position_of_boundary_node()

template<class ELEMENT >
void oomph::HermiteQuadMesh< ELEMENT >::set_position_of_boundary_node ( const unsigned &  node_num_x,
const unsigned &  node_num_y,
BoundaryNode< Node > *  node_pt 
)
private

sets the generalised position of the node (i.e. - x_i, dx_i/ds_0, dx_i/ds_1 & d2x_i/ds_0ds_1 for i = 1,2). Takes the x,y coordinates of the node from which its position can be determined. Also sets coordinates on boundary vector for the node to be the generalised position of the node in macro element coordinates

120  {
121  // get the generalised macro element position of the node
122  DenseMatrix<double> m_gen(4, 2, 0.0);
123  generalised_macro_element_position_of_node(node_num_x, node_num_y, m_gen);
124 
125  // copy the macro element coordinates
126  Vector<double> node_macro_position(2);
127  node_macro_position[0] = m_gen(0, 0);
128  node_macro_position[1] = m_gen(0, 1);
129 
130  // get the global position of the nodes
131  Vector<double> x_node(2);
132  Domain_pt->macro_element_pt(0)->macro_map(node_macro_position, x_node);
133 
134  // get the jacobian of the mapping from macro to eulerian coordinates
135  DenseMatrix<double> jacobian_MX(2, 2);
136  Domain_pt->macro_element_pt(0)->assemble_macro_to_eulerian_jacobian(
137  node_macro_position, jacobian_MX);
138 
139  // get the jacobian2 of the mapping from macro to eulerian coordinates
140  DenseMatrix<double> jacobian2_MX(3, 2);
141  Domain_pt->macro_element_pt(0)->assemble_macro_to_eulerian_jacobian2(
142  node_macro_position, jacobian2_MX);
143 
144  // set x_0
145  node_pt->x_gen(0, 0) = x_node[0];
146  // set x_1
147  node_pt->x_gen(0, 1) = x_node[1];
148 
149  // set dxi/dsj for i = 0,1 and j = 0,1
150  // computation : dxi/dsj = dm0/dsj*dxi/dm0 + dm1/dsj*dxi/dm1
151  for (unsigned i = 0; i < 2; i++)
152  {
153  for (unsigned j = 0; j < 2; j++)
154  {
155  node_pt->x_gen(j + 1, i) = m_gen(j + 1, 0) * jacobian_MX(0, i) +
156  m_gen(j + 1, 1) * jacobian_MX(1, i);
157  }
158  }
159 
160  // set d2xi/ds0ds1 for i = 0,1
161  // d2xi/ds0ds1 = d2m0/ds0ds1*dxi/dm0
162  // + d2m1/ds0ds1*dxi/dm1
163  // + (dm0/ds1*d2xi/dm0^2 + dm1/ds1*d2xi/dm0dm1)*dm0/ds0
164  // + (dm0/ds1*d2xi/dm0dm1 + dm1/ds1*d2xi/dm1^2)*dm1/ds0
165  for (unsigned i = 0; i < 2; i++)
166  {
167  node_pt->x_gen(3, i) =
168  m_gen(3, 0) * jacobian_MX(0, i) + m_gen(3, 1) * jacobian_MX(1, i) +
169  m_gen(1, 0) * (m_gen(2, 0) * jacobian2_MX(0, i) +
170  m_gen(2, 1) * jacobian2_MX(2, i)) +
171  m_gen(1, 1) *
172  (m_gen(2, 0) * jacobian2_MX(2, i) + m_gen(2, 1) * jacobian2_MX(1, i));
173  }
174 
175 
176  // pass generalised macro element position to vector to pass to boundary
177  // node
178  for (unsigned b = 0; b < 4; b++)
179  {
180  if (node_pt->is_on_boundary(b))
181  {
182  Vector<double> boundary_macro_position(2, 0);
183  boundary_macro_position[0] = m_gen(0, b % 2);
184  boundary_macro_position[1] = m_gen(1 + b % 2, b % 2);
185  node_pt->set_coordinates_on_boundary(b, boundary_macro_position);
186  }
187  }
188  }
b
Scalar * b
Definition: benchVecAdd.cpp:17
oomph::HermiteQuadMesh::generalised_macro_element_position_of_node
void generalised_macro_element_position_of_node(const unsigned &node_num_x, const unsigned &node_num_y, DenseMatrix< double > &m_gen)
Definition: hermite_element_quad_mesh.template.cc:200
oomph::MacroElement::assemble_macro_to_eulerian_jacobian2
virtual void assemble_macro_to_eulerian_jacobian2(const unsigned &t, const Vector< double > &s, DenseMatrix< double > &jacobian2)
Definition: macro_element.h:186
oomph::MacroElement::macro_map
void macro_map(const Vector< double > &s, Vector< double > &r)
The mapping from local to global coordinates at the current time : r(s)
Definition: macro_element.h:126
oomph::MacroElement::assemble_macro_to_eulerian_jacobian
virtual void assemble_macro_to_eulerian_jacobian(const unsigned &t, const Vector< double > &s, DenseMatrix< double > &jacobian)
Definition: macro_element.h:170
oomph::Node::x_gen
double & x_gen(const unsigned &k, const unsigned &i)
Definition: nodes.h:1126
oomph::Node::is_on_boundary
virtual bool is_on_boundary() const
Definition: nodes.h:1373
oomph::Node::set_coordinates_on_boundary
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Definition: nodes.cc:2394

References b, i, oomph::BoundaryNode< NODE_TYPE >::is_on_boundary(), j, and oomph::BoundaryNode< NODE_TYPE >::set_coordinates_on_boundary().

◆ set_position_of_node()

template<class ELEMENT >
void oomph::HermiteQuadMesh< ELEMENT >::set_position_of_node ( const unsigned &  node_num_x,
const unsigned &  node_num_y,
Node *  node_pt 
)
private

sets the generalised position of the node (i.e. - x_i, dx_i/ds_0, dx_i/ds_1 & d2x_i/ds_0ds_1 for i = 1,2). Takes the x,y coordinates of the node from which its position can be determined.

50  {
51  // get the generalised macro element position of the node
52  DenseMatrix<double> m_gen(4, 2);
53  generalised_macro_element_position_of_node(node_num_x, node_num_y, m_gen);
54 
55  // copy the macro element coordinates
56  Vector<double> node_macro_position(2);
57  node_macro_position[0] = m_gen(0, 0);
58  node_macro_position[1] = m_gen(0, 1);
59 
60  // get the global position of the nodes
61  Vector<double> x_node(2);
62  Domain_pt->macro_element_pt(0)->macro_map(node_macro_position, x_node);
63 
64  // get the jacobian of the mapping from macro to eulerian coordinates
65  DenseMatrix<double> jacobian_MX(2, 2);
66  Domain_pt->macro_element_pt(0)->assemble_macro_to_eulerian_jacobian(
67  node_macro_position, jacobian_MX);
68 
69 
70  // get the jacobian2 of the mapping from macro to eulerian coordinates
71  DenseMatrix<double> jacobian2_MX(3, 2);
72  Domain_pt->macro_element_pt(0)->assemble_macro_to_eulerian_jacobian2(
73  node_macro_position, jacobian2_MX);
74 
75  // set x_0
76  node_pt->x_gen(0, 0) = x_node[0];
77  // set x_1
78  node_pt->x_gen(0, 1) = x_node[1];
79 
80  // set dxi/dsj for i = 0,1 and j = 0,1
81  // computation : dxi/dsj = dm0/dsj*dxi/dm0 + dm1/dsj*dxi/dm1
82  for (unsigned i = 0; i < 2; i++)
83  {
84  for (unsigned j = 0; j < 2; j++)
85  {
86  node_pt->x_gen(j + 1, i) = m_gen(j + 1, 0) * jacobian_MX(0, i) +
87  m_gen(j + 1, 1) * jacobian_MX(1, i);
88  }
89  }
90 
91  // set d2xi/ds0ds1 for i = 0,1
92  // d2xi/ds0ds1 = d2m0/ds0ds1*dxi/dm0
93  // + d2m1/ds0ds1*dxi/dm1
94  // + (dm0/ds1*d2xi/dm0^2 + dm1/ds1*d2xi/dm0dm1)*dm0/ds0
95  // + (dm0/ds1*d2xi/dm0dm1 + dm1/ds1*d2xi/dm1^2)*dm1/ds0
96  for (unsigned i = 0; i < 2; i++)
97  {
98  node_pt->x_gen(3, i) =
99  m_gen(3, 0) * jacobian_MX(0, i) + m_gen(3, 1) * jacobian_MX(1, i) +
100  m_gen(1, 0) * (m_gen(2, 0) * jacobian2_MX(0, i) +
101  m_gen(2, 1) * jacobian2_MX(2, i)) +
102  m_gen(1, 1) *
103  (m_gen(2, 0) * jacobian2_MX(2, i) + m_gen(2, 1) * jacobian2_MX(1, i));
104  }
105  }

References i, j, and oomph::Node::x_gen().

◆ setup_boundary_element_info() [1/2]

template<class ELEMENT >
virtual void oomph::HermiteQuadMesh< ELEMENT >::setup_boundary_element_info ( )
inlineprivatevirtual

Setup lookup schemes which establish whic elements are located next to mesh's boundaries (wrapper to suppress doc). Specific version for HermiteQuadMesh to ensure that the order of the elements in Boundary_element_pt matches the actual order along the boundary. This is required when hijacking the BiharmonicElement to apply the BiharmonicFluidBoundaryElement in BiharmonicFluidProblem::impose_traction_free_edge(...)

Reimplemented from oomph::Mesh.

216  {
217  std::ofstream outfile;
218  setup_boundary_element_info(outfile);
219  }

◆ setup_boundary_element_info() [2/2]

template<class ELEMENT >
void oomph::HermiteQuadMesh< ELEMENT >::setup_boundary_element_info ( std::ostream &  outfile)
privatevirtual

Setup lookup schemes which establish which elements are located next to which boundaries (Doc to outfile if it's open). Specific version for HermiteQuadMesh to ensure that the order of the elements in Boundary_element_pt matches the actual order along the boundary. This is required when hijacking the BiharmonicElement to apply the BiharmonicFluidBoundaryElement in BiharmonicFluidProblem::impose_traction_free_edge(...)

Reimplemented from oomph::Mesh.

1273  {
1274  bool doc = false;
1275  if (outfile) doc = true;
1276 
1277  // Number of boundaries
1278  unsigned nbound = nboundary();
1279 
1280  // Wipe/allocate storage for arrays
1281  Boundary_element_pt.clear();
1282  Face_index_at_boundary.clear();
1283  Boundary_element_pt.resize(nbound);
1284  Face_index_at_boundary.resize(nbound);
1285 
1286  // Temporary vector of sets of pointers to elements on the boundaries:
1287  Vector<Vector<FiniteElement*>> vector_of_boundary_element_pt;
1288  vector_of_boundary_element_pt.resize(nbound);
1289 
1290  // Matrix map for working out the fixed local coord for elements on boundary
1291  MapMatrixMixed<unsigned, FiniteElement*, Vector<int>*> boundary_identifier;
1292 
1293 
1294  // Loop over elements
1295  //-------------------
1296  unsigned nel = nelement();
1297  for (unsigned e = 0; e < nel; e++)
1298  {
1299  // Get pointer to element
1300  FiniteElement* fe_pt = finite_element_pt(e);
1301 
1302  if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
1303 
1304  // Only include 2D elements! Some meshes contain interface elements too.
1305  if (fe_pt->dim() == 2)
1306  {
1307  // Loop over the element's nodes and find out which boundaries they're
1308  // on
1309  // ----------------------------------------------------------------------
1310  unsigned nnode_1d = fe_pt->nnode_1d();
1311 
1312  // Loop over nodes in order
1313  for (unsigned i0 = 0; i0 < nnode_1d; i0++)
1314  {
1315  for (unsigned i1 = 0; i1 < nnode_1d; i1++)
1316  {
1317  // Local node number
1318  unsigned j = i0 + i1 * nnode_1d;
1319 
1320  // Get pointer to vector of boundaries that this
1321  // node lives on
1322  std::set<unsigned>* boundaries_pt = 0;
1323  fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);
1324 
1325  // If the node lives on some boundaries....
1326  if (boundaries_pt != 0)
1327  {
1328  // Loop over boundaries
1329  // unsigned nbound=(*boundaries_pt).size();
1330  for (std::set<unsigned>::iterator it = boundaries_pt->begin();
1331  it != boundaries_pt->end();
1332  ++it)
1333  {
1334  // Add pointer to finite element to set for the appropriate
1335  // boundary -- storage in set makes sure we don't count elements
1336  // multiple times
1337  unsigned temp_size = vector_of_boundary_element_pt[*it].size();
1338  bool temp_flag = true;
1339  for (unsigned temp_i = 0; temp_i < temp_size; temp_i++)
1340  {
1341  if (vector_of_boundary_element_pt[*it][temp_i] == fe_pt)
1342  temp_flag = false;
1343  }
1344  if (temp_flag)
1345  vector_of_boundary_element_pt[*it].push_back(fe_pt);
1346 
1347  // For the current element/boundary combination, create
1348  // a vector that stores an indicator which element boundaries
1349  // the node is located (boundary_identifier=-/+1 for nodes
1350  // on the left/right boundary; boundary_identifier=-/+2 for
1351  // nodes on the lower/upper boundary. We determine these indices
1352  // for all corner nodes of the element and add them to a vector
1353  // to a vector. This allows us to decide which face of the
1354  // element coincides with the boundary since the (quad!) element
1355  // must have exactly two corner nodes on the boundary.
1356  if (boundary_identifier(*it, fe_pt) == 0)
1357  {
1358  boundary_identifier(*it, fe_pt) = new Vector<int>;
1359  }
1360 
1361  // Are we at a corner node?
1362  if (((i0 == 0) || (i0 == nnode_1d - 1)) &&
1363  ((i1 == 0) || (i1 == nnode_1d - 1)))
1364  {
1365  // Create index to represent position relative to s_0
1366  (*boundary_identifier(*it, fe_pt))
1367  .push_back(1 * (2 * i0 / (nnode_1d - 1) - 1));
1368 
1369  // Create index to represent position relative to s_1
1370  (*boundary_identifier(*it, fe_pt))
1371  .push_back(2 * (2 * i1 / (nnode_1d - 1) - 1));
1372  }
1373  }
1374  }
1375  // else
1376  // {
1377  // oomph_info << "...does not live on any boundaries " <<
1378  // std::endl;
1379  // }
1380  }
1381  }
1382  }
1383  // else
1384  //{
1385  // oomph_info << "Element " << e << " does not qualify" << std::endl;
1386  //}
1387  }
1388 
1389 
1390  // Now copy everything across into permanent arrays
1391  //-------------------------------------------------
1392 
1393  // Note: vector_of_boundary_element_pt contains all elements
1394  // that have (at least) one corner node on a boundary -- can't copy
1395  // them across into Boundary_element_pt
1396  // yet because some of them might have only one node on the
1397  // the boundary in which case they don't qualify as
1398  // boundary elements!
1399 
1400  // Loop over boundaries
1401  //---------------------
1402  for (unsigned i = 0; i < nbound; i++)
1403  {
1404  // Number of elements on this boundary (currently stored in a set)
1405  unsigned nel = vector_of_boundary_element_pt[i].size();
1406 
1407  // Allocate storage for the coordinate identifiers
1408  Face_index_at_boundary[i].resize(nel);
1409 
1410  // Loop over elements that have at least one corner node on this boundary
1411  //-----------------------------------------------------------------------
1412  unsigned e_count = 0;
1413  typedef Vector<FiniteElement*>::iterator IT;
1414  for (IT it = vector_of_boundary_element_pt[i].begin();
1415  it != vector_of_boundary_element_pt[i].end();
1416  it++)
1417  {
1418  // Recover pointer to element
1419  FiniteElement* fe_pt = *it;
1420 
1421  // Initialise count for boundary identiers (-2,-1,1,2)
1422  std::map<int, unsigned> count;
1423 
1424  // Loop over coordinates
1425  for (unsigned ii = 0; ii < 2; ii++)
1426  {
1427  // Loop over upper/lower end of coordinates
1428  for (int sign = -1; sign < 3; sign += 2)
1429  {
1430  count[(ii + 1) * sign] = 0;
1431  }
1432  }
1433 
1434  // Loop over boundary indicators for this element/boundary
1435  unsigned n_indicators = (*boundary_identifier(i, fe_pt)).size();
1436  for (unsigned j = 0; j < n_indicators; j++)
1437  {
1438  count[(*boundary_identifier(i, fe_pt))[j]]++;
1439  }
1440  delete boundary_identifier(i, fe_pt);
1441 
1442  // Determine the correct boundary indicator by checking that it
1443  // occurs twice (since two corner nodes of the element's boundary
1444  // need to be located on the domain boundary
1445  int indicator = -10;
1446 
1447  // Check that we're finding exactly one boundary indicator
1448  bool found = false;
1449 
1450  // Loop over coordinates
1451  for (unsigned ii = 0; ii < 2; ii++)
1452  {
1453  // Loop over upper/lower end of coordinates
1454  for (int sign = -1; sign < 3; sign += 2)
1455  {
1456  if (count[(ii + 1) * sign] == 2)
1457  {
1458  // Check that we haven't found multiple boundaries
1459  if (found)
1460  {
1461  std::string error_message =
1462  "Trouble: Multiple boundary identifiers!\n";
1463  error_message += "Elements should only have at most 2 corner ";
1464  error_message += "nodes on any one boundary.\n";
1465 
1466  throw OomphLibError(error_message,
1467  OOMPH_CURRENT_FUNCTION,
1468  OOMPH_EXCEPTION_LOCATION);
1469  }
1470  found = true;
1471  indicator = (ii + 1) * sign;
1472  }
1473  }
1474  }
1475 
1476  // Element has exactly two corner nodes on boundary
1477  if (found)
1478  {
1479  // Add to permanent storage
1480  Boundary_element_pt[i].push_back(fe_pt);
1481 
1482  // Now convert boundary indicator into information required
1483  // for FaceElements
1484  switch (indicator)
1485  {
1486  case -2:
1487 
1488  // s_1 is fixed at -1.0:
1489  Face_index_at_boundary[i][e_count] = -2;
1490  break;
1491 
1492  case -1:
1493 
1494  // s_0 is fixed at -1.0:
1495  Face_index_at_boundary[i][e_count] = -1;
1496  break;
1497 
1498 
1499  case 1:
1500 
1501  // s_0 is fixed at 1.0:
1502  Face_index_at_boundary[i][e_count] = 1;
1503  break;
1504 
1505  case 2:
1506 
1507  // s_1 is fixed at 1.0:
1508  Face_index_at_boundary[i][e_count] = 2;
1509  break;
1510 
1511  default:
1512 
1513  throw OomphLibError("Never get here",
1514  OOMPH_CURRENT_FUNCTION,
1515  OOMPH_EXCEPTION_LOCATION);
1516  }
1517 
1518  // Increment counter
1519  e_count++;
1520  }
1521  }
1522  }
1523 
1524 
1525  // Doc?
1526  //-----
1527  if (doc)
1528  {
1529  // Loop over boundaries
1530  for (unsigned i = 0; i < nbound; i++)
1531  {
1532  unsigned nel = Boundary_element_pt[i].size();
1533  outfile << "Boundary: " << i << " is adjacent to " << nel << " elements"
1534  << std::endl;
1535 
1536  // Loop over elements on given boundary
1537  for (unsigned e = 0; e < nel; e++)
1538  {
1539  FiniteElement* fe_pt = Boundary_element_pt[i][e];
1540  outfile << "Boundary element:" << fe_pt
1541  << " Face index of element along boundary is "
1542  << Face_index_at_boundary[i][e] << std::endl;
1543  }
1544  }
1545  }
1546 
1547 
1548  // Lookup scheme has now been setup yet
1549  Lookup_for_elements_next_boundary_is_setup = true;
1550  }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
oomph::Mesh::Boundary_element_pt
Vector< Vector< FiniteElement * > > Boundary_element_pt
Definition: mesh.h:91
oomph::Mesh::Lookup_for_elements_next_boundary_is_setup
bool Lookup_for_elements_next_boundary_is_setup
Definition: mesh.h:87
oomph::Mesh::Face_index_at_boundary
Vector< Vector< int > > Face_index_at_boundary
Definition: mesh.h:95
oomph::Mesh::nboundary
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
oomph::Mesh::nelement
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
MergeRestartFiles.found
bool found
Definition: MergeRestartFiles.py:24
SYCL::sign
T sign(T x)
Definition: cxx11_tensor_builtins_sycl.cpp:172
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References oomph::FiniteElement::dim(), e(), MergeRestartFiles::found, oomph::Node::get_boundaries_pt(), i, j, oomph::FiniteElement::nnode_1d(), oomph::FiniteElement::node_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, SYCL::sign(), size, and oomph::Global_string_for_annotation::string().

Member Data Documentation

◆ Domain_pt

template<class ELEMENT >
TopologicallyRectangularDomain* oomph::HermiteQuadMesh< ELEMENT >::Domain_pt
private

Pointer to the topologically rectangular domain which prescribes the problem domain

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

◆ Nelement

template<class ELEMENT >
Vector<unsigned> oomph::HermiteQuadMesh< ELEMENT >::Nelement
private

number of elements in each coordinate direction

Referenced by oomph::HermiteQuadMesh< ELEMENT >::HermiteQuadMesh(), oomph::HermiteQuadMesh< ELEMENT >::macro_coordinate_position(), and oomph::HermiteQuadMesh< ELEMENT >::nelement_in_dim().

◆ Node_spacing_fn

template<class ELEMENT >
MeshSpacingFnPtr oomph::HermiteQuadMesh< ELEMENT >::Node_spacing_fn
private

non uniform mesh spacing function pointer

Referenced by oomph::HermiteQuadMesh< ELEMENT >::HermiteQuadMesh(), and oomph::HermiteQuadMesh< ELEMENT >::macro_coordinate_position().

◆ Xperiodic

template<class ELEMENT >
bool oomph::HermiteQuadMesh< ELEMENT >::Xperiodic
private

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

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

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