Public Member Functions | Protected Member Functions | List of all members
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D > Class Template Reference

p-refineable version of RefineableQElement<2,INITIAL_NNODE_1D>. More...

#include <hp_refineable_elements.h>

+ Inheritance diagram for oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >:

Public Member Functions

 PRefineableQElement ()
 Constructor. More...
 
virtual ~PRefineableQElement ()
 Destructor. More...
 
void initial_setup (Tree *const &adopted_father_pt=0, const unsigned &initial_p_order=0)
 
void pre_build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt)
 
void p_refine (const int &inc, Mesh *const &mesh_pt, GeneralisedElement *const &clone_pt)
 p-refine the element (refine if inc>0, unrefine if inc<0). More...
 
void shape (const Vector< double > &s, Shape &psi) const
 Overload the shape functions. More...
 
void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsi) const
 Derivatives of shape functions for PRefineableQElement<DIM> More...
 
void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
void further_setup_hanging_nodes ()
 
unsigned nnode_1d () const
 
unsigned initial_p_order () const
 Get the initial P_order. More...
 
Node * get_node_at_local_coordinate (const Vector< double > &s) const
 Return the node at the specified local coordinate. More...
 
Node * node_created_by_neighbour (const Vector< double > &s_fraction, bool &is_periodic)
 
Node * node_created_by_son_of_neighbour (const Vector< double > &s_fraction, bool &is_periodic)
 
void local_coordinate_of_node (const unsigned &n, Vector< double > &s) const
 Get local coordinates of node j in the element; vector sets its own size. More...
 
void local_fraction_of_node (const unsigned &n, Vector< double > &s_fraction)
 Get the local fractino of node j in the element. More...
 
double local_one_d_fraction_of_node (const unsigned &n1d, const unsigned &i)
 The local one-d fraction is the same. More...
 
void rebuild_from_sons (Mesh *&mesh_pt)
 
void check_integrity (double &max_error)
 
- Public Member Functions inherited from oomph::RefineableQElement< 2 >
 RefineableQElement ()
 Constructor: Pass refinement level (default 0 = root) More...
 
 RefineableQElement (const RefineableQElement< 2 > &dummy)=delete
 Broken copy constructor. More...
 
virtual ~RefineableQElement ()
 Broken assignment operator. More...
 
unsigned required_nsons () const
 A refineable quad element has four sons. More...
 
virtual void build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)
 
void check_integrity (double &max_error)
 
void output_corners (std::ostream &outfile, const std::string &colour) const
 Print corner nodes, use colour. More...
 
QuadTree * quadtree_pt ()
 Pointer to quadtree representation of this element. More...
 
QuadTree * quadtree_pt () const
 Pointer to quadtree representation of this element. More...
 
void setup_hanging_nodes (Vector< std::ofstream * > &output_stream)
 
void get_boundaries (const int &edge, std::set< unsigned > &boundaries) const
 
void get_bcs (int bound, Vector< int > &bound_cons) const
 
void interpolated_zeta_on_edge (const unsigned &boundary, const int &edge, const Vector< double > &s, Vector< double > &zeta)
 
- Public Member Functions inherited from oomph::RefineableElement
 RefineableElement ()
 
virtual ~RefineableElement ()
 Destructor, delete the allocated storage for the hanging equations. More...
 
 RefineableElement (const RefineableElement &)=delete
 Broken copy constructor. More...
 
void operator= (const RefineableElement &)=delete
 Broken assignment operator. More...
 
Tree * tree_pt ()
 Access function: Pointer to quadtree representation of this element. More...
 
void set_tree_pt (Tree *my_tree_pt)
 Set pointer to quadtree representation of this element. More...
 
bool refinement_is_enabled ()
 Flag to indicate suppression of any refinement. More...
 
void disable_refinement ()
 Suppress of any refinement for this element. More...
 
void enable_refinement ()
 Emnable refinement for this element. More...
 
template<class ELEMENT >
void split (Vector< ELEMENT * > &son_pt) const
 
int local_hang_eqn (Node *const &node_pt, const unsigned &i)
 
void set_refinement_level (const int &refine_level)
 Set the refinement level. More...
 
unsigned refinement_level () const
 Return the Refinement level. More...
 
void select_for_refinement ()
 Select the element for refinement. More...
 
void deselect_for_refinement ()
 Deselect the element for refinement. More...
 
void select_sons_for_unrefinement ()
 Unrefinement will be performed by merging the four sons of this element. More...
 
void deselect_sons_for_unrefinement ()
 
bool to_be_refined ()
 Has the element been selected for refinement? More...
 
bool sons_to_be_unrefined ()
 Has the element been selected for unrefinement? More...
 
virtual void unbuild ()
 
virtual void deactivate_element ()
 
long number () const
 Element number (for debugging/plotting) More...
 
void set_number (const long &mynumber)
 Set element number (for debugging/plotting) More...
 
virtual unsigned ncont_interpolated_values () const =0
 
virtual void get_interpolated_values (const Vector< double > &s, Vector< double > &values)
 
virtual void get_interpolated_values (const unsigned &t, const Vector< double > &s, Vector< double > &values)=0
 
virtual Node * interpolating_node_pt (const unsigned &n, const int &value_id)
 
virtual double local_one_d_fraction_of_interpolating_node (const unsigned &n1d, const unsigned &i, const int &value_id)
 
virtual Node * get_interpolating_node_at_local_coordinate (const Vector< double > &s, const int &value_id)
 
virtual unsigned ninterpolating_node (const int &value_id)
 
virtual unsigned ninterpolating_node_1d (const int &value_id)
 
virtual void interpolating_basis (const Vector< double > &s, Shape &psi, const int &value_id) const
 
void identify_field_data_for_interactions (std::set< std::pair< Data *, unsigned >> &paired_field_data)
 
void assign_nodal_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual RefineableElement * root_element_pt ()
 
virtual RefineableElement * father_element_pt () const
 Return a pointer to the father element. More...
 
void get_father_at_refinement_level (unsigned &refinement_level, RefineableElement *&father_at_reflevel_pt)
 
virtual void further_build ()
 Further build: e.g. deal with interpolation of internal values. More...
 
void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
unsigned nshape_controlling_nodes ()
 
std::map< Node *, unsigned > shape_controlling_node_lookup ()
 
- Public Member Functions inherited from oomph::FiniteElement
void set_dimension (const unsigned &dim)
 
void set_nodal_dimension (const unsigned &nodal_dim)
 
void set_nnodal_position_type (const unsigned &nposition_type)
 Set the number of types required to interpolate the coordinate. More...
 
void set_n_node (const unsigned &n)
 
int nodal_local_eqn (const unsigned &n, const unsigned &i) const
 
double dJ_eulerian_at_knot (const unsigned &ipt, Shape &psi, DenseMatrix< double > &djacobian_dX) const
 
 FiniteElement ()
 Constructor. More...
 
virtual ~FiniteElement ()
 
 FiniteElement (const FiniteElement &)=delete
 Broken copy constructor. More...
 
void get_centre_of_gravity_and_max_radius_in_terms_of_zeta (Vector< double > &cog, double &max_radius) const
 
MacroElement * macro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void describe_nodal_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void assign_all_generic_local_eqn_numbers (const bool &store_local_dof_pt)
 
Node *& node_pt (const unsigned &n)
 Return a pointer to the local node n. More...
 
Node *const & node_pt (const unsigned &n) const
 Return a pointer to the local node n (const version) More...
 
unsigned nnode () const
 Return the number of nodes. More...
 
double raw_nodal_position (const unsigned &n, const unsigned &i) const
 
double raw_nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position (const unsigned &n, const unsigned &i) const
 
double nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double dnodal_position_dt (const unsigned &n, const unsigned &i) const
 Return the i-th component of nodal velocity: dx/dt at local node n. More...
 
double dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual void disable_ALE ()
 
virtual void enable_ALE ()
 
virtual unsigned required_nvalue (const unsigned &n) const
 
unsigned nnodal_position_type () const
 
bool has_hanging_nodes () const
 
unsigned nodal_dimension () const
 Return the required Eulerian dimension of the nodes in this element. More...
 
virtual Node * construct_node (const unsigned &n)
 
virtual Node * construct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual Node * construct_boundary_node (const unsigned &n)
 
virtual Node * construct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
int get_node_number (Node *const &node_pt) const
 
double raw_nodal_value (const unsigned &n, const unsigned &i) const
 
double raw_nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
unsigned dim () const
 
virtual double interpolated_x (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated coordinate x[i] at local coordinate s. More...
 
virtual double interpolated_x (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_x (const Vector< double > &s, Vector< double > &x) const
 Return FE interpolated position x[] at local coordinate s as Vector. More...
 
virtual void interpolated_x (const unsigned &t, const Vector< double > &s, Vector< double > &x) const
 
virtual double interpolated_dxdt (const Vector< double > &s, const unsigned &i, const unsigned &t)
 
virtual void interpolated_dxdt (const Vector< double > &s, const unsigned &t, Vector< double > &dxdt)
 
unsigned ngeom_data () const
 
Data * geom_data_pt (const unsigned &j)
 
void position (const Vector< double > &zeta, Vector< double > &r) const
 
void position (const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
 
void dposition_dt (const Vector< double > &zeta, const unsigned &t, Vector< double > &drdt)
 
virtual double zeta_nodal (const unsigned &n, const unsigned &k, const unsigned &i) const
 
void interpolated_zeta (const Vector< double > &s, Vector< double > &zeta) const
 
void locate_zeta (const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
 
virtual void node_update ()
 
virtual void identify_geometric_data (std::set< Data * > &geometric_data_pt)
 
virtual double s_min () const
 Min value of local coordinate. More...
 
virtual double s_max () const
 Max. value of local coordinate. More...
 
double size () const
 
virtual double compute_physical_size () const
 
virtual void point_output_data (const Vector< double > &s, Vector< double > &data)
 
void point_output (std::ostream &outfile, const Vector< double > &s)
 
virtual unsigned nplot_points_paraview (const unsigned &nplot) const
 
virtual unsigned nsub_elements_paraview (const unsigned &nplot) const
 
void output_paraview (std::ofstream &file_out, const unsigned &nplot) const
 
virtual void write_paraview_output_offset_information (std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
 
virtual void write_paraview_type (std::ofstream &file_out, const unsigned &nplot) const
 
virtual void write_paraview_offsets (std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
 
virtual unsigned nscalar_paraview () const
 
virtual void scalar_value_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
 
virtual void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
 
virtual void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
virtual std::string scalar_name_paraview (const unsigned &i) const
 
virtual void output (std::ostream &outfile)
 
virtual void output (std::ostream &outfile, const unsigned &n_plot)
 
virtual void output (const unsigned &t, std::ostream &outfile, const unsigned &n_plot) const
 
virtual void output (FILE *file_pt)
 
virtual void output (FILE *file_pt, const unsigned &n_plot)
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 Output an exact solution over the element. More...
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 Output a time-dependent exact solution over the element. More...
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, const SolutionFunctorBase &exact_soln) const
 Output a time-dependent exact solution over the element. More...
 
virtual void get_s_plot (const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
 
virtual std::string tecplot_zone_string (const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (std::ostream &outfile, const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (FILE *file_pt, const unsigned &nplot) const
 
virtual unsigned nplot_points (const unsigned &nplot) const
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (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, double &error, 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, 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_abs_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error)
 
void integrate_fct (FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt, Vector< double > &integral)
 Integrate Vector-valued function over element. More...
 
void integrate_fct (FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt, const double &time, Vector< double > &integral)
 Integrate Vector-valued time-dep function over element. More...
 
virtual void build_face_element (const int &face_index, FaceElement *face_element_pt)
 
virtual unsigned self_test ()
 
virtual unsigned get_bulk_node_number (const int &face_index, const unsigned &i) const
 
virtual int face_outer_unit_normal_sign (const int &face_index) const
 Get the sign of the outer unit normal on the face given by face_index. More...
 
void face_node_number_error_check (const unsigned &i) const
 Range check for face node numbers. More...
 
virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt (const int &face_index) const
 
virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt (const int &face_index) const
 
- Public Member Functions inherited from oomph::GeneralisedElement
 GeneralisedElement ()
 Constructor: Initialise all pointers and all values to zero. More...
 
virtual ~GeneralisedElement ()
 Virtual destructor to clean up any memory allocated by the object. More...
 
 GeneralisedElement (const GeneralisedElement &)=delete
 Broken copy constructor. More...
 
void operator= (const GeneralisedElement &)=delete
 Broken assignment operator. More...
 
Data *& internal_data_pt (const unsigned &i)
 Return a pointer to i-th internal data object. More...
 
Data *const & internal_data_pt (const unsigned &i) const
 Return a pointer to i-th internal data object (const version) More...
 
Data *& external_data_pt (const unsigned &i)
 Return a pointer to i-th external data object. More...
 
Data *const & external_data_pt (const unsigned &i) const
 Return a pointer to i-th external data object (const version) More...
 
unsigned long eqn_number (const unsigned &ieqn_local) const
 
int local_eqn_number (const unsigned long &ieqn_global) const
 
unsigned add_external_data (Data *const &data_pt, const bool &fd=true)
 
bool external_data_fd (const unsigned &i) const
 
void exclude_external_data_fd (const unsigned &i)
 
void include_external_data_fd (const unsigned &i)
 
void flush_external_data ()
 Flush all external data. More...
 
void flush_external_data (Data *const &data_pt)
 Flush the object addressed by data_pt from the external data array. More...
 
unsigned ninternal_data () const
 Return the number of internal data objects. More...
 
unsigned nexternal_data () const
 Return the number of external data objects. More...
 
unsigned ndof () const
 Return the number of equations/dofs in the element. More...
 
void dof_vector (const unsigned &t, Vector< double > &dof)
 Return the vector of dof values at time level t. More...
 
void dof_pt_vector (Vector< double * > &dof_pt)
 Return the vector of pointers to dof values. More...
 
void set_internal_data_time_stepper (const unsigned &i, TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void assign_internal_eqn_numbers (unsigned long &global_number, Vector< double * > &Dof_pt)
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void add_internal_value_pt_to_map (std::map< unsigned, double * > &map_of_value_pt)
 
virtual void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void complete_setup_of_dependencies ()
 
virtual void get_residuals (Vector< double > &residuals)
 
virtual void get_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
virtual void get_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void get_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void get_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void get_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void get_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void get_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void get_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void get_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_norm (double &norm)
 
virtual unsigned ndof_types () const
 
virtual void get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
 
- Public Member Functions inherited from oomph::GeomObject
 GeomObject ()
 Default constructor. More...
 
 GeomObject (const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)
 
 GeomObject (const GeomObject &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const GeomObject &)=delete
 Broken assignment operator. More...
 
virtual ~GeomObject ()
 (Empty) destructor More...
 
unsigned nlagrangian () const
 Access function to # of Lagrangian coordinates. More...
 
unsigned ndim () const
 Access function to # of Eulerian coordinates. More...
 
void set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
 Set # of Lagrangian and Eulerian coordinates. More...
 
TimeStepper *& time_stepper_pt ()
 
TimeStepper * time_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 
- Public Member Functions inherited from oomph::QuadElementBase
 QuadElementBase ()
 Constructor. Empty. More...
 
virtual unsigned nvertex_node () const =0
 Number of vertex nodes in the element. More...
 
virtual Node * vertex_node_pt (const unsigned &j) const =0
 Pointer to the j-th vertex node in the element. More...
 
- Public Member Functions inherited from oomph::QElementBase
 QElementBase ()
 Constructor: Initialise pointers to macro element reference coords. More...
 
 QElementBase (const QElementBase &)=delete
 Broken copy constructor. More...
 
virtual ~QElementBase ()
 Broken assignment operator. More...
 
bool local_coord_is_valid (const Vector< double > &s)
 Check whether the local coordinate are valid or not. More...
 
void move_local_coord_back_into_element (Vector< double > &s) const
 
virtual void set_macro_elem_pt (MacroElement *macro_elem_pt)
 
double & s_macro_ll (const unsigned &i)
 
double & s_macro_ur (const unsigned &i)
 
double s_macro_ll (const unsigned &i) const
 
double s_macro_ur (const unsigned &i) const
 
void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
unsigned nnode_on_face () const
 
ElementGeometry::ElementGeometry element_geometry () const
 It's a Q element! More...
 
- Public Member Functions inherited from oomph::QElementGeometricBase
 QElementGeometricBase ()
 Empty default constructor. More...
 
 QElementGeometricBase (const QElementGeometricBase &)=delete
 Broken copy constructor. More...
 
- Public Member Functions inherited from oomph::PRefineableElement
 PRefineableElement ()
 Constructor, calls the RefineableElement constructor. More...
 
virtual ~PRefineableElement ()
 Destructor, empty. More...
 
 PRefineableElement (const PRefineableElement &)=delete
 Broken copy constructor. More...
 
void operator= (const PRefineableElement &)=delete
 Broken assignment operator. More...
 
bool p_refinement_is_enabled ()
 Flag to indicate suppression of any refinement. More...
 
void disable_p_refinement ()
 Suppress of any refinement for this element. More...
 
void enable_p_refinement ()
 Emnable refinement for this element. More...
 
unsigned & p_order ()
 Access function to P_order. More...
 
unsigned p_order () const
 Access function to P_order (const version) More...
 
void select_for_p_refinement ()
 Select the element for p-refinement. More...
 
void deselect_for_p_refinement ()
 Deselect the element for p-refinement. More...
 
void select_for_p_unrefinement ()
 Select the element for p-unrefinement. More...
 
void deselect_for_p_unrefinement ()
 Deselect the element for p-unrefinement. More...
 
bool to_be_p_refined ()
 Has the element been selected for refinement? More...
 
bool to_be_p_unrefined ()
 Has the element been selected for p-unrefinement? More...
 
bool nodes_built ()
 Return true if all the nodes have been built, false if not. More...
 

Protected Member Functions

void quad_hang_helper (const int &value_id, const int &my_edge, std::ofstream &output_hangfile)
 
- Protected Member Functions inherited from oomph::RefineableQElement< 2 >
void setup_father_bounds ()
 
void get_edge_bcs (const int &edge, Vector< int > &bound_cons) const
 Determine Vector of boundary conditions along edge (N/S/W/E) More...
 
void setup_hang_for_value (const int &value_id)
 
- Protected Member Functions inherited from oomph::RefineableElement
void assemble_local_to_eulerian_jacobian (const DShape &dpsids, DenseMatrix< double > &jacobian) const
 
void assemble_local_to_eulerian_jacobian2 (const DShape &d2psids, DenseMatrix< double > &jacobian2) const
 
void assemble_eulerian_base_vectors (const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
 
double local_to_eulerian_mapping_diagonal (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
void assign_hanging_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign the local equation numbers for hanging node variables. More...
 
virtual void fill_in_jacobian_from_nodal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
- Protected Member Functions inherited from oomph::FiniteElement
template<unsigned DIM>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double invert_jacobian_mapping (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &inverse_jacobian) const
 
virtual void dJ_eulerian_dnodal_coordinates (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<unsigned DIM>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
virtual void d_dshape_eulerian_dnodal_coordinates (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<unsigned DIM>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
virtual void transform_derivatives (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
void transform_derivatives_diagonal (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
virtual void transform_second_derivatives (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
void fill_in_jacobian_from_nodal_by_fd (DenseMatrix< double > &jacobian)
 
virtual void update_before_nodal_fd ()
 
virtual void reset_after_nodal_fd ()
 
virtual void update_in_nodal_fd (const unsigned &i)
 
virtual void reset_in_nodal_fd (const unsigned &i)
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Zero-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 One-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Two-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
- Protected Member Functions inherited from oomph::GeneralisedElement
unsigned add_internal_data (Data *const &data_pt, const bool &fd=true)
 
bool internal_data_fd (const unsigned &i) const
 
void exclude_internal_data_fd (const unsigned &i)
 
void include_internal_data_fd (const unsigned &i)
 
void clear_global_eqn_numbers ()
 
void add_global_eqn_numbers (std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)
 
virtual void assign_internal_and_external_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void assign_additional_local_eqn_numbers ()
 
int internal_local_eqn (const unsigned &i, const unsigned &j) const
 
int external_local_eqn (const unsigned &i, const unsigned &j)
 
virtual void fill_in_contribution_to_residuals (Vector< double > &residuals)
 
void fill_in_jacobian_from_internal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_internal_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
virtual void update_before_internal_fd ()
 
virtual void reset_after_internal_fd ()
 
virtual void update_in_internal_fd (const unsigned &i)
 
virtual void reset_in_internal_fd (const unsigned &i)
 
virtual void update_before_external_fd ()
 
virtual void reset_after_external_fd ()
 
virtual void update_in_external_fd (const unsigned &i)
 
virtual void reset_in_external_fd (const unsigned &i)
 
virtual void fill_in_contribution_to_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void fill_in_contribution_to_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void fill_in_contribution_to_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void fill_in_contribution_to_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void fill_in_contribution_to_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void fill_in_contribution_to_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void fill_in_contribution_to_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 

Additional Inherited Members

- Public Types inherited from oomph::RefineableQElement< 2 >
typedef void(RefineableQElement< 2 >::* VoidMemberFctPt) ()
 
- Public Types inherited from oomph::FiniteElement
typedef void(* SteadyExactSolutionFctPt) (const Vector< double > &, Vector< double > &)
 
typedef void(* UnsteadyExactSolutionFctPt) (const double &, const Vector< double > &, Vector< double > &)
 
- Static Public Member Functions inherited from oomph::RefineableElement
static double & max_integrity_tolerance ()
 Max. allowed discrepancy in element integrity check. More...
 
- Static Public Attributes inherited from oomph::FiniteElement
static double Tolerance_for_singular_jacobian = 1.0e-16
 Tolerance below which the jacobian is considered singular. More...
 
static bool Accept_negative_jacobian = false
 
static bool Suppress_output_while_checking_for_inverted_elements
 
- Static Public Attributes inherited from oomph::GeneralisedElement
static bool Suppress_warning_about_repeated_internal_data
 
static bool Suppress_warning_about_repeated_external_data = true
 
static double Default_fd_jacobian_step = 1.0e-8
 
- Static Protected Member Functions inherited from oomph::RefineableElement
static void check_value_id (const int &n_continuously_interpolated_values, const int &value_id)
 
- Protected Attributes inherited from oomph::RefineableElement
Tree * Tree_pt
 A pointer to a general tree object. More...
 
unsigned Refine_level
 Refinement level. More...
 
bool To_be_refined
 Flag for refinement. More...
 
bool Refinement_is_enabled
 Flag to indicate suppression of any refinement. More...
 
bool Sons_to_be_unrefined
 Flag for unrefinement. More...
 
long Number
 Global element number – for plotting/validation purposes. More...
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElement * Macro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepper * Geom_object_time_stepper_pt
 
- Protected Attributes inherited from oomph::PRefineableElement
unsigned P_order
 The polynomial expansion order of the elemental basis functions. More...
 
bool To_be_p_refined
 Flag for p-refinement. More...
 
bool P_refinement_is_enabled
 Flag to indicate suppression of any refinement. More...
 
bool To_be_p_unrefined
 Flag for unrefinement. More...
 
- Static Protected Attributes inherited from oomph::RefineableQElement< 2 >
static std::map< unsigned, DenseMatrix< int > > Father_bound
 
- Static Protected Attributes inherited from oomph::RefineableElement
static double Max_integrity_tolerance = 1.0e-8
 Max. allowed discrepancy in element integrity check. More...
 
- Static Protected Attributes inherited from oomph::FiniteElement
static const unsigned Default_Initial_Nvalue = 0
 Default value for the number of values at a node. More...
 
static const double Node_location_tolerance = 1.0e-14
 
static const unsigned N2deriv [] = {0, 1, 3, 6}
 
- Static Protected Attributes inherited from oomph::GeneralisedElement
static DenseMatrix< double > Dummy_matrix
 
static std::deque< double * > Dof_pt_deque
 

Detailed Description

template<unsigned INITIAL_NNODE_1D>
class oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >

p-refineable version of RefineableQElement<2,INITIAL_NNODE_1D>.

Constructor & Destructor Documentation

◆ PRefineableQElement()

template<unsigned INITIAL_NNODE_1D>
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::PRefineableQElement ( )
inline

Constructor.

149 : PRefineableElement(), RefineableQElement<2>() {}
oomph::PRefineableElement::PRefineableElement
PRefineableElement()
Constructor, calls the RefineableElement constructor.
Definition: refineable_elements.h:668

◆ ~PRefineableQElement()

template<unsigned INITIAL_NNODE_1D>
virtual oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::~PRefineableQElement ( )
inlinevirtual

Destructor.

152 {}

Member Function Documentation

◆ check_integrity()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::check_integrity ( double &  max_error)
virtual

Check the integrity of interpolated values across element boundaries. Note: with the mortar method, continuity is enforced weakly across non- conforming element boundaries, so it makes no sense to check the continuity of interpolated values across these boundaries.

Check inter-element continuity of

  • nodal positions
  • (nodally) interpolated function values Overloaded to not check differences in the value. Mortaring doesn't enforce strong continuity between elements.

Implements oomph::RefineableElement.

3750  {
3751  // Overloaded to *not* check for continuity in value of interpolated
3752  // variables. This is necessary because mortaring does not ensure continuity
3753  // across element boundaries. It therefore makes no sense to test for this.
3754 
3755  // Dummy set max_error to 0
3756  max_error = 0.0;
3757 
3758  // With macro-elements, (strong) continuity in position is nolonger
3759  // guaranteed either, so we don't check for this either. In fact, we do
3760  // nothing at all.
3761  if (this->macro_elem_pt() != 0)
3762  {
3763  // We have a macro element, so do nothing!
3764  return;
3765  }
3766 
3767  using namespace QuadTreeNames;
3768 
3769  // Number of nodes along edge
3770  unsigned n_p = nnode_1d();
3771 
3772  // Number of timesteps (incl. present) for which continuity is
3773  // to be checked.
3774  unsigned n_time = 1;
3775 
3776  // Initialise errors
3777  max_error = 0.0;
3778  Vector<double> max_error_x(2, 0.0);
3779  double max_error_val = 0.0;
3780 
3781  Vector<int> edges(4);
3782  edges[0] = S;
3783  edges[1] = N;
3784  edges[2] = W;
3785  edges[3] = E;
3786 
3787  // Loop over the edges
3788  for (unsigned edge_counter = 0; edge_counter < 4; edge_counter++)
3789  {
3790  Vector<unsigned> translate_s(2);
3791  Vector<double> s(2), s_lo_neigh(2), s_hi_neigh(2), s_fraction(2);
3792  int neigh_edge, diff_level;
3793  bool in_neighbouring_tree;
3794 
3795  // Find pointer to neighbour in this direction
3796  QuadTree* neigh_pt;
3797  neigh_pt = quadtree_pt()->gteq_edge_neighbour(edges[edge_counter],
3798  translate_s,
3799  s_lo_neigh,
3800  s_hi_neigh,
3801  neigh_edge,
3802  diff_level,
3803  in_neighbouring_tree);
3804 
3805  // Neighbour exists and has existing nodes
3806  if ((neigh_pt != 0) && (neigh_pt->object_pt()->nodes_built()))
3807  {
3808  // Need to exclude periodic nodes from this check
3809  // There are only periodic nodes if we are in a neighbouring tree
3810  bool is_periodic = false;
3811  if (in_neighbouring_tree)
3812  {
3813  // Is it periodic
3814  is_periodic = this->tree_pt()->root_pt()->is_neighbour_periodic(
3815  edges[edge_counter]);
3816  }
3817 
3818  // We also need to exclude edges which may have hanging nodes
3819  // because mortaring does not guarantee (strong) continuity
3820  // in position or in value at nonconforming element boundaries
3821  bool exclude_this_edge = false;
3822  if (diff_level != 0)
3823  {
3824  // h-type nonconformity (dependent)
3825  exclude_this_edge = true;
3826  }
3827  else if (neigh_pt->nsons() != 0)
3828  {
3829  // h-type nonconformity (master)
3830  exclude_this_edge = true;
3831  }
3832  else
3833  {
3834  unsigned my_p_order = this->p_order();
3835  unsigned neigh_p_order =
3836  dynamic_cast<PRefineableQElement*>(neigh_pt->object_pt())
3837  ->p_order();
3838  if (my_p_order != neigh_p_order)
3839  {
3840  // p-type nonconformity
3841  exclude_this_edge = true;
3842  }
3843  }
3844 
3845  // With macro-elements, (strong) continuity in position is nolonger
3846  // guaranteed either, so we don't check for this either. In fact, we do
3847  // nothing at all.
3848  if (dynamic_cast<FiniteElement*>(neigh_pt->object_pt())
3849  ->macro_elem_pt() != 0)
3850  {
3851  // We have a macro element, so do nothing!
3852  break;
3853  }
3854 
3855  // Check conforming edges
3856  if (!exclude_this_edge)
3857  {
3858  // Loop over nodes along the edge
3859  for (unsigned i0 = 0; i0 < n_p; i0++)
3860  {
3861  // Storage for pointer to the local node
3862  Node* local_node_pt = 0;
3863 
3864  switch (edge_counter)
3865  {
3866  case 0:
3867  // Local fraction of node
3868  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3869  s_fraction[1] = 0.0;
3870  // Get pointer to local node
3871  local_node_pt = this->node_pt(i0);
3872  break;
3873 
3874  case 1:
3875  // Local fraction of node
3876  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3877  s_fraction[1] = 1.0;
3878  // Get pointer to local node
3879  local_node_pt = this->node_pt(i0 + n_p * (n_p - 1));
3880  break;
3881 
3882  case 2:
3883  // Local fraction of node
3884  s_fraction[0] = 0.0;
3885  s_fraction[1] = this->local_one_d_fraction_of_node(i0, 1);
3886  // Get pointer to local node
3887  local_node_pt = this->node_pt(n_p * i0);
3888  break;
3889 
3890  case 3:
3891  // Local fraction of node
3892  s_fraction[0] = 1.0;
3893  s_fraction[1] = this->local_one_d_fraction_of_node(i0, 1);
3894  // Get pointer to local node
3895  local_node_pt = this->node_pt(n_p - 1 + n_p * i0);
3896  break;
3897  }
3898 
3899  // Calculate the local coordinate and the local coordinate as viewed
3900  // from the neighbour
3901  Vector<double> s_in_neighb(2);
3902  for (unsigned i = 0; i < 2; i++)
3903  {
3904  // Local coordinate in this element
3905  s[i] = -1.0 + 2.0 * s_fraction[i];
3906  // Local coordinate in the neighbour
3907  s_in_neighb[i] =
3908  s_lo_neigh[i] +
3909  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
3910  }
3911 
3912  // Loop over timesteps
3913  for (unsigned t = 0; t < n_time; t++)
3914  {
3915  // Get the nodal position from neighbour element
3916  Vector<double> x_in_neighb(2);
3917  neigh_pt->object_pt()->interpolated_x(
3918  t, s_in_neighb, x_in_neighb);
3919 
3920  // Check error only if the node is NOT periodic
3921  if (is_periodic == false)
3922  {
3923  for (int i = 0; i < 2; i++)
3924  {
3925  // Find the spatial error
3926  double err =
3927  std::fabs(local_node_pt->x(t, i) - x_in_neighb[i]);
3928 
3929  // If it's bigger than our tolerance, say so
3930  if (err > 1e-9)
3931  {
3932  oomph_info << "errx " << err << " " << t << " "
3933  << local_node_pt->x(t, i) << " "
3934  << x_in_neighb[i] << std::endl;
3935 
3936  oomph_info << "at " << local_node_pt->x(0) << " "
3937  << local_node_pt->x(1) << std::endl;
3938  }
3939 
3940  // If it's bigger than the previous max error, it is the
3941  // new max error!
3942  if (err > max_error_x[i])
3943  {
3944  max_error_x[i] = err;
3945  }
3946  }
3947  }
3948 
3949  // Get the values from neighbour element. Note: # of values
3950  // gets set by routine (because in general we don't know
3951  // how many interpolated values a certain element has
3952  Vector<double> values_in_neighb;
3953  neigh_pt->object_pt()->get_interpolated_values(
3954  t, s_in_neighb, values_in_neighb);
3955 
3956  // Get the values in current element.
3957  Vector<double> values;
3958  this->get_interpolated_values(t, s, values);
3959 
3960  // Now figure out how many continuously interpolated values there
3961  // are
3962  unsigned num_val =
3963  neigh_pt->object_pt()->ncont_interpolated_values();
3964 
3965  // Check error
3966  for (unsigned ival = 0; ival < num_val; ival++)
3967  {
3968  double err = std::fabs(values[ival] - values_in_neighb[ival]);
3969 
3970  if (err > 1.0e-10)
3971  {
3972  oomph_info << local_node_pt->x(0) << " "
3973  << local_node_pt->x(1) << " \n# "
3974  << "erru (S)" << err << " " << ival << " "
3975  << this->get_node_number(local_node_pt) << " "
3976  << values[ival] << " " << values_in_neighb[ival]
3977  << std::endl;
3978  }
3979 
3980  if (err > max_error_val)
3981  {
3982  max_error_val = err;
3983  }
3984  }
3985  }
3986  }
3987  }
3988  }
3989  }
3990 
3991  max_error = max_error_x[0];
3992  if (max_error_x[1] > max_error) max_error = max_error_x[1];
3993  if (max_error_val > max_error) max_error = max_error_val;
3994 
3995  if (max_error > 1e-9)
3996  {
3997  oomph_info << "\n#------------------------------------ \n#Max error ";
3998  oomph_info << max_error_x[0] << " " << max_error_x[1] << " "
3999  << max_error_val << std::endl;
4000  oomph_info << "#------------------------------------ \n " << std::endl;
4001  }
4002  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::get_node_number
int get_node_number(Node *const &node_pt) const
Definition: elements.cc:3814
oomph::FiniteElement::macro_elem_pt
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Definition: elements.h:1878
oomph::FiniteElement::FiniteElement
FiniteElement()
Constructor.
Definition: elements.h:1782
oomph::PRefineableElement::p_order
unsigned & p_order()
Access function to P_order.
Definition: refineable_elements.h:705
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::PRefineableQElement
PRefineableQElement()
Constructor.
Definition: hp_refineable_elements.h:149
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::nnode_1d
unsigned nnode_1d() const
Definition: hp_refineable_elements.h:187
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::local_one_d_fraction_of_node
double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
The local one-d fraction is the same.
Definition: hp_refineable_elements.cc:1272
oomph::QuadTree::gteq_edge_neighbour
QuadTree * gteq_edge_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, int &edge, int &diff_level, bool &in_neighbouring_tree) const
Definition: quadtree.cc:413
oomph::RefineableElement::tree_pt
Tree * tree_pt()
Access function: Pointer to quadtree representation of this element.
Definition: refineable_elements.h:211
oomph::RefineableElement::get_interpolated_values
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Definition: refineable_elements.h:417
oomph::RefineableQElement< 2 >::quadtree_pt
QuadTree * quadtree_pt()
Pointer to quadtree representation of this element.
Definition: refineable_quad_element.h:152
oomph::TreeRoot::is_neighbour_periodic
bool is_neighbour_periodic(const int &direction)
Definition: tree.h:364
oomph::Tree::root_pt
TreeRoot *& root_pt()
Return pointer to root of the tree.
Definition: tree.h:141
N
@ N
Definition: constructor.cpp:22
s
RealScalar s
Definition: level1_cplx_impl.h:130
Global_Physical_Variables::E
double E
Elastic modulus.
Definition: TwenteMeshGluing.cpp:68
MeshRefinement::max_error
double max_error
Definition: MortaringCantileverCompareToNonMortaring.cpp:188
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
oomph::QuadTreeNames::S
@ S
Definition: quadtree.h:62
oomph::QuadTreeNames::W
@ W
Definition: quadtree.h:63
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319
plotPSD.t
t
Definition: plotPSD.py:36

References Global_Physical_Variables::E, e(), boost::multiprecision::fabs(), oomph::RefineableElement::get_interpolated_values(), oomph::QuadTree::gteq_edge_neighbour(), i, oomph::FiniteElement::interpolated_x(), oomph::FiniteElement::macro_elem_pt(), MeshRefinement::max_error, N, oomph::RefineableElement::ncont_interpolated_values(), oomph::RefineableElement::nodes_built(), oomph::Tree::nsons(), oomph::Tree::object_pt(), oomph::oomph_info, s, oomph::QuadTreeNames::S, plotPSD::t, oomph::QuadTreeNames::W, and oomph::Node::x().

◆ d2shape_local()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::d2shape_local ( const Vector< double > &  s,
Shape &  psi,
DShape &  dpsids,
DShape &  d2psids 
) const
virtual

Second derivatives of shape functions for PRefineableQElement<DIM> d2psids(i,0) = \( d^2 \psi_j / d s^2 \)

Reimplemented from oomph::FiniteElement.

3194  {
3195  std::ostringstream error_message;
3196  error_message
3197  << "\nd2shape_local currently not implemented for this element\n";
3198  throw OomphLibError(
3199  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3200  }
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_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ dshape_local()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::dshape_local ( const Vector< double > &  s,
Shape &  psi,
DShape &  dpsi 
) const
virtual

Derivatives of shape functions for PRefineableQElement<DIM>

Reimplemented from oomph::FiniteElement.

3030  {
3031  switch (p_order())
3032  {
3033  case 2:
3034  {
3035  // Call the shape functions and derivatives
3036  OneDimensionalLegendreShape<2>::calculate_nodal_positions();
3037  OneDimensionalLegendreShape<2> psi1(s[0]), psi2(s[1]);
3038  OneDimensionalLegendreDShape<2> dpsi1ds(s[0]), dpsi2ds(s[1]);
3039 
3040  // Index for the shape functions
3041  unsigned index = 0;
3042  // Loop over shape functions in element
3043  for (unsigned i = 0; i < 2; i++)
3044  {
3045  for (unsigned j = 0; j < 2; j++)
3046  {
3047  // Assign the values
3048  dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3049  dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3050  psi[index] = psi2[i] * psi1[j];
3051  // Increment the index
3052  ++index;
3053  }
3054  }
3055  break;
3056  }
3057  case 3:
3058  {
3059  // Call the shape functions and derivatives
3060  OneDimensionalLegendreShape<3>::calculate_nodal_positions();
3061  OneDimensionalLegendreShape<3> psi1(s[0]), psi2(s[1]);
3062  OneDimensionalLegendreDShape<3> dpsi1ds(s[0]), dpsi2ds(s[1]);
3063 
3064  // Index for the shape functions
3065  unsigned index = 0;
3066  // Loop over shape functions in element
3067  for (unsigned i = 0; i < 3; i++)
3068  {
3069  for (unsigned j = 0; j < 3; j++)
3070  {
3071  // Assign the values
3072  dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3073  dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3074  psi[index] = psi2[i] * psi1[j];
3075  // Increment the index
3076  ++index;
3077  }
3078  }
3079  break;
3080  }
3081  case 4:
3082  {
3083  // Call the shape functions and derivatives
3084  OneDimensionalLegendreShape<4>::calculate_nodal_positions();
3085  OneDimensionalLegendreShape<4> psi1(s[0]), psi2(s[1]);
3086  OneDimensionalLegendreDShape<4> dpsi1ds(s[0]), dpsi2ds(s[1]);
3087 
3088  // Index for the shape functions
3089  unsigned index = 0;
3090  // Loop over shape functions in element
3091  for (unsigned i = 0; i < 4; i++)
3092  {
3093  for (unsigned j = 0; j < 4; j++)
3094  {
3095  // Assign the values
3096  dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3097  dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3098  psi[index] = psi2[i] * psi1[j];
3099  // Increment the index
3100  ++index;
3101  }
3102  }
3103  break;
3104  }
3105  case 5:
3106  {
3107  // Call the shape functions and derivatives
3108  OneDimensionalLegendreShape<5>::calculate_nodal_positions();
3109  OneDimensionalLegendreShape<5> psi1(s[0]), psi2(s[1]);
3110  OneDimensionalLegendreDShape<5> dpsi1ds(s[0]), dpsi2ds(s[1]);
3111 
3112  // Index for the shape functions
3113  unsigned index = 0;
3114  // Loop over shape functions in element
3115  for (unsigned i = 0; i < 5; i++)
3116  {
3117  for (unsigned j = 0; j < 5; j++)
3118  {
3119  // Assign the values
3120  dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3121  dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3122  psi[index] = psi2[i] * psi1[j];
3123  // Increment the index
3124  ++index;
3125  }
3126  }
3127  break;
3128  }
3129  case 6:
3130  {
3131  // Call the shape functions and derivatives
3132  OneDimensionalLegendreShape<6>::calculate_nodal_positions();
3133  OneDimensionalLegendreShape<6> psi1(s[0]), psi2(s[1]);
3134  OneDimensionalLegendreDShape<6> dpsi1ds(s[0]), dpsi2ds(s[1]);
3135 
3136  // Index for the shape functions
3137  unsigned index = 0;
3138  // Loop over shape functions in element
3139  for (unsigned i = 0; i < 6; i++)
3140  {
3141  for (unsigned j = 0; j < 6; j++)
3142  {
3143  // Assign the values
3144  dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3145  dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3146  psi[index] = psi2[i] * psi1[j];
3147  // Increment the index
3148  ++index;
3149  }
3150  }
3151  break;
3152  }
3153  case 7:
3154  {
3155  // Call the shape functions and derivatives
3156  OneDimensionalLegendreShape<7>::calculate_nodal_positions();
3157  OneDimensionalLegendreShape<7> psi1(s[0]), psi2(s[1]);
3158  OneDimensionalLegendreDShape<7> dpsi1ds(s[0]), dpsi2ds(s[1]);
3159 
3160  // Index for the shape functions
3161  unsigned index = 0;
3162  // Loop over shape functions in element
3163  for (unsigned i = 0; i < 7; i++)
3164  {
3165  for (unsigned j = 0; j < 7; j++)
3166  {
3167  // Assign the values
3168  dpsids(index, 0) = psi2[i] * dpsi1ds[j];
3169  dpsids(index, 1) = dpsi2ds[i] * psi1[j];
3170  psi[index] = psi2[i] * psi1[j];
3171  // Increment the index
3172  ++index;
3173  }
3174  }
3175  break;
3176  }
3177  default:
3178  std::ostringstream error_message;
3179  error_message << "\nERROR: Exceeded maximum polynomial order for";
3180  error_message << "\n polynomial order for shape functions.\n";
3181  throw OomphLibError(error_message.str(),
3182  OOMPH_CURRENT_FUNCTION,
3183  OOMPH_EXCEPTION_LOCATION);
3184  }
3185  }
oomph::OneDimensionalLegendreShape::calculate_nodal_positions
static void calculate_nodal_positions()
Static function used to populate the stored positions.
Definition: shape.h:1241
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::OneDimensionalLegendreShape< NNODE_1D >::calculate_nodal_positions(), i, j, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

◆ further_setup_hanging_nodes()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::further_setup_hanging_nodes ( )
inlinevirtual

Perform additional hanging node procedures for variables that are not interpolated by all nodes (e.g. lower order interpolations for the pressure in Taylor Hood).

Implements oomph::RefineableQElement< 2 >.

182 {}

◆ get_node_at_local_coordinate()

template<unsigned INITIAL_NNODE_1D>
Node * oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::get_node_at_local_coordinate ( const Vector< double > &  s) const
virtual

Return the node at the specified local coordinate.

Reimplemented from oomph::FiniteElement.

1319  {
1320  unsigned Nnode_1d = this->nnode_1d();
1321  // Load the tolerance into a local variable
1322  double tol = FiniteElement::Node_location_tolerance;
1323  // There are two possible indices.
1324  Vector<int> index(2);
1325 
1326  // Loop over indices
1327  for (unsigned i = 0; i < 2; i++)
1328  {
1329  // Determine the index
1330  // -------------------
1331 
1332  bool is_found = false;
1333 
1334  // If we are at the lower limit, the index is zero
1335  if (std::fabs(s[i] + 1.0) < tol)
1336  {
1337  index[i] = 0;
1338  is_found = true;
1339  }
1340  // If we are at the upper limit, the index is the number of nodes minus 1
1341  else if (std::fabs(s[i] - 1.0) < tol)
1342  {
1343  index[i] = Nnode_1d - 1;
1344  is_found = true;
1345  }
1346  // Otherwise, we have to calculate the index in general
1347  else
1348  {
1349  // Compute Gauss-Lobatto-Legendre node positions
1350  Vector<double> z;
1351  Orthpoly::gll_nodes(Nnode_1d, z);
1352  // Loop over possible internal nodal positions
1353  for (unsigned n = 1; n < Nnode_1d - 1; n++)
1354  {
1355  if (std::fabs(z[n] - s[i]) < tol)
1356  {
1357  index[i] = n;
1358  is_found = true;
1359  break;
1360  }
1361  }
1362  }
1363 
1364  if (!is_found)
1365  {
1366  // No matching nodes
1367  return 0;
1368  }
1369  }
1370  // If we've got here we have a node, so let's return a pointer to it
1371  return this->node_pt(index[0] + Nnode_1d * index[1]);
1372  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
oomph::FiniteElement::Node_location_tolerance
static const double Node_location_tolerance
Definition: elements.h:1374
GlobalParameters::Nnode_1d
unsigned Nnode_1d
The number of nodes in one direction (default=2)
Definition: structured_cubic_point_source.cc:56
oomph::Orthpoly::gll_nodes
void gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
Definition: orthpoly.cc:33

References boost::multiprecision::fabs(), oomph::Orthpoly::gll_nodes(), i, n, GlobalParameters::Nnode_1d, oomph::FiniteElement::Node_location_tolerance, and s.

Referenced by p_refine().

◆ initial_p_order()

template<unsigned INITIAL_NNODE_1D>
unsigned oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::initial_p_order ( ) const
inlinevirtual

Get the initial P_order.

Implements oomph::PRefineableElement.

194  {
195  return INITIAL_NNODE_1D;
196  }

◆ initial_setup()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::initial_setup ( Tree *const &  adopted_father_pt = 0,
const unsigned &  initial_p_order = 0 
)
virtual

Initial setup of element (set the correct p-order and integration scheme) If an adopted father is specified, information from this is used instead of using the father found from the tree.

Set the correct p-order of the element based on that of its father. Then construct an integration scheme of the correct order. If an adopted father is specified, information from this is used instead of using the father found from the tree.

Reimplemented from oomph::RefineableElement.

1669  {
1670  // Storage for pointer to my father (in quadtree impersonation)
1671  QuadTree* father_pt = 0;
1672 
1673  // Check if an adopted father has been specified
1674  if (adopted_father_pt != 0)
1675  {
1676  // Get pointer to my father (in quadtree impersonation)
1677  father_pt = dynamic_cast<QuadTree*>(adopted_father_pt);
1678  }
1679  // Check if element is in a tree
1680  else if (Tree_pt != 0)
1681  {
1682  // Get pointer to my father (in quadtree impersonation)
1683  father_pt = dynamic_cast<QuadTree*>(quadtree_pt()->father_pt());
1684  }
1685  // else
1686  // {
1687  // throw OomphLibError(
1688  // "Element not in a tree, and no adopted father has been
1689  // specified!", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1690  // }
1691 
1692  // Check if we found a father
1693  if (father_pt != 0 || initial_p_order != 0)
1694  {
1695  if (father_pt != 0)
1696  {
1697  PRefineableQElement<2, INITIAL_NNODE_1D>* father_el_pt =
1698  dynamic_cast<PRefineableQElement<2, INITIAL_NNODE_1D>*>(
1699  father_pt->object_pt());
1700  if (father_el_pt != 0)
1701  {
1702  unsigned father_p_order = father_el_pt->p_order();
1703  // Set p-order to that of father
1704  P_order = father_p_order;
1705  }
1706  }
1707  else
1708  {
1709  P_order = initial_p_order;
1710  }
1711 
1712  // Now sort out the element...
1713  // (has p^2 nodes)
1714  unsigned new_n_node = P_order * P_order;
1715 
1716  // Allocate new space for Nodes (at the element level)
1717  this->set_n_node(new_n_node);
1718 
1719  // Set integration scheme
1720  delete this->integral_pt();
1721  switch (P_order)
1722  {
1723  case 2:
1724  this->set_integration_scheme(new GaussLobattoLegendre<2, 2>);
1725  break;
1726  case 3:
1727  this->set_integration_scheme(new GaussLobattoLegendre<2, 3>);
1728  break;
1729  case 4:
1730  this->set_integration_scheme(new GaussLobattoLegendre<2, 4>);
1731  break;
1732  case 5:
1733  this->set_integration_scheme(new GaussLobattoLegendre<2, 5>);
1734  break;
1735  case 6:
1736  this->set_integration_scheme(new GaussLobattoLegendre<2, 6>);
1737  break;
1738  case 7:
1739  this->set_integration_scheme(new GaussLobattoLegendre<2, 7>);
1740  break;
1741  default:
1742  std::ostringstream error_message;
1743  error_message << "\nERROR: Exceeded maximum polynomial order for";
1744  error_message << "\n integration scheme.\n";
1745  throw OomphLibError(error_message.str(),
1746  OOMPH_CURRENT_FUNCTION,
1747  OOMPH_EXCEPTION_LOCATION);
1748  }
1749  }
1750  }
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
oomph::FiniteElement::set_n_node
void set_n_node(const unsigned &n)
Definition: elements.h:1404
oomph::FiniteElement::set_integration_scheme
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3210
oomph::PRefineableElement::P_order
unsigned P_order
The polynomial expansion order of the elemental basis functions.
Definition: refineable_elements.h:655
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::initial_p_order
unsigned initial_p_order() const
Get the initial P_order.
Definition: hp_refineable_elements.h:193
oomph::RefineableElement::Tree_pt
Tree * Tree_pt
A pointer to a general tree object.
Definition: refineable_elements.h:100
oomph::Tree::father_pt
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:235

References oomph::Tree::object_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::PRefineableElement::p_order().

◆ local_coordinate_of_node()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::local_coordinate_of_node ( const unsigned &  n,
Vector< double > &  s 
) const
virtual

Get local coordinates of node j in the element; vector sets its own size.

Reimplemented from oomph::FiniteElement.

1212  {
1213  s.resize(2);
1214  unsigned Nnode_1d = this->nnode_1d();
1215  unsigned n0 = n % Nnode_1d;
1216  unsigned n1 = n / Nnode_1d;
1217 
1218  switch (Nnode_1d)
1219  {
1220  case 2:
1221  OneDimensionalLegendreShape<2>::calculate_nodal_positions();
1222  s[0] = OneDimensionalLegendreShape<2>::nodal_position(n0);
1223  s[1] = OneDimensionalLegendreShape<2>::nodal_position(n1);
1224  break;
1225  case 3:
1226  OneDimensionalLegendreShape<3>::calculate_nodal_positions();
1227  s[0] = OneDimensionalLegendreShape<3>::nodal_position(n0);
1228  s[1] = OneDimensionalLegendreShape<3>::nodal_position(n1);
1229  break;
1230  case 4:
1231  OneDimensionalLegendreShape<4>::calculate_nodal_positions();
1232  s[0] = OneDimensionalLegendreShape<4>::nodal_position(n0);
1233  s[1] = OneDimensionalLegendreShape<4>::nodal_position(n1);
1234  break;
1235  case 5:
1236  OneDimensionalLegendreShape<5>::calculate_nodal_positions();
1237  s[0] = OneDimensionalLegendreShape<5>::nodal_position(n0);
1238  s[1] = OneDimensionalLegendreShape<5>::nodal_position(n1);
1239  break;
1240  case 6:
1241  OneDimensionalLegendreShape<6>::calculate_nodal_positions();
1242  s[0] = OneDimensionalLegendreShape<6>::nodal_position(n0);
1243  s[1] = OneDimensionalLegendreShape<6>::nodal_position(n1);
1244  break;
1245  case 7:
1246  OneDimensionalLegendreShape<7>::calculate_nodal_positions();
1247  s[0] = OneDimensionalLegendreShape<7>::nodal_position(n0);
1248  s[1] = OneDimensionalLegendreShape<7>::nodal_position(n1);
1249  break;
1250  default:
1251  std::ostringstream error_message;
1252  error_message << "\nERROR: Exceeded maximum polynomial order for";
1253  error_message << "\n shape functions.\n";
1254  throw OomphLibError(error_message.str(),
1255  OOMPH_CURRENT_FUNCTION,
1256  OOMPH_EXCEPTION_LOCATION);
1257  break;
1258  }
1259  }
oomph::OneDimensionalLegendreShape::nodal_position
static double nodal_position(const unsigned &n)
Definition: shape.h:1250

References oomph::OneDimensionalLegendreShape< NNODE_1D >::calculate_nodal_positions(), n, GlobalParameters::Nnode_1d, oomph::OneDimensionalLegendreShape< NNODE_1D >::nodal_position(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

◆ local_fraction_of_node()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::local_fraction_of_node ( const unsigned &  n,
Vector< double > &  s_fraction 
)
virtual

Get the local fractino of node j in the element.

Reimplemented from oomph::FiniteElement.

1265  {
1266  this->local_coordinate_of_node(n, s_fraction);
1267  s_fraction[0] = 0.5 * (s_fraction[0] + 1.0);
1268  s_fraction[1] = 0.5 * (s_fraction[1] + 1.0);
1269  }
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::local_coordinate_of_node
void local_coordinate_of_node(const unsigned &n, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size.
Definition: hp_refineable_elements.cc:1210

References n.

◆ local_one_d_fraction_of_node()

template<unsigned INITIAL_NNODE_1D>
double oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::local_one_d_fraction_of_node ( const unsigned &  n1d,
const unsigned &  i 
)
virtual

The local one-d fraction is the same.

Reimplemented from oomph::FiniteElement.

1274  {
1275  switch (this->nnode_1d())
1276  {
1277  case 2:
1278  OneDimensionalLegendreShape<2>::calculate_nodal_positions();
1279  return 0.5 *
1280  (OneDimensionalLegendreShape<2>::nodal_position(n1d) + 1.0);
1281  case 3:
1282  OneDimensionalLegendreShape<3>::calculate_nodal_positions();
1283  return 0.5 *
1284  (OneDimensionalLegendreShape<3>::nodal_position(n1d) + 1.0);
1285  case 4:
1286  OneDimensionalLegendreShape<4>::calculate_nodal_positions();
1287  return 0.5 *
1288  (OneDimensionalLegendreShape<4>::nodal_position(n1d) + 1.0);
1289  case 5:
1290  OneDimensionalLegendreShape<5>::calculate_nodal_positions();
1291  return 0.5 *
1292  (OneDimensionalLegendreShape<5>::nodal_position(n1d) + 1.0);
1293  case 6:
1294  OneDimensionalLegendreShape<6>::calculate_nodal_positions();
1295  return 0.5 *
1296  (OneDimensionalLegendreShape<6>::nodal_position(n1d) + 1.0);
1297  case 7:
1298  OneDimensionalLegendreShape<7>::calculate_nodal_positions();
1299  return 0.5 *
1300  (OneDimensionalLegendreShape<7>::nodal_position(n1d) + 1.0);
1301  default:
1302  std::ostringstream error_message;
1303  error_message << "\nERROR: Exceeded maximum polynomial order for";
1304  error_message << "\n shape functions.\n";
1305  throw OomphLibError(error_message.str(),
1306  OOMPH_CURRENT_FUNCTION,
1307  OOMPH_EXCEPTION_LOCATION);
1308  return 0.0;
1309  }
1310  }

References oomph::OneDimensionalLegendreShape< NNODE_1D >::calculate_nodal_positions(), oomph::OneDimensionalLegendreShape< NNODE_1D >::nodal_position(), OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ nnode_1d()

template<unsigned INITIAL_NNODE_1D>
unsigned oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::nnode_1d ( ) const
inlinevirtual

Returns the number of nodes along each edge of the element. Overloaded to return the (variable) p-order rather than the template argument.

Reimplemented from oomph::FiniteElement.

188  {
189  return this->p_order();
190  }

◆ node_created_by_neighbour()

template<unsigned INITIAL_NNODE_1D>
Node * oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::node_created_by_neighbour ( const Vector< double > &  s_fraction,
bool &  is_periodic 
)
virtual

If a neighbouring element has already created a node at a position corresponding to the local fractional position within the present element, s_fraction, return a pointer to that node. If not, return NULL (0). If the node is periodic the flag is_periodic will be true

Reimplemented from oomph::RefineableQElement< 2 >.

1384  {
1385  using namespace QuadTreeNames;
1386 
1387  // Calculate the edges on which the node lies
1388  Vector<int> edges;
1389  if (s_fraction[0] == 0.0)
1390  {
1391  edges.push_back(W);
1392  }
1393  if (s_fraction[0] == 1.0)
1394  {
1395  edges.push_back(E);
1396  }
1397  if (s_fraction[1] == 0.0)
1398  {
1399  edges.push_back(S);
1400  }
1401  if (s_fraction[1] == 1.0)
1402  {
1403  edges.push_back(N);
1404  }
1405 
1406  // Find the number of edges
1407  unsigned n_size = edges.size();
1408  // If there are no edges, then there is no neighbour, return 0
1409  if (n_size == 0)
1410  {
1411  return 0;
1412  }
1413 
1414  Vector<unsigned> translate_s(2);
1415  Vector<double> s_lo_neigh(2);
1416  Vector<double> s_hi_neigh(2);
1417  Vector<double> s(2);
1418 
1419  int neigh_edge, diff_level;
1420  bool in_neighbouring_tree;
1421 
1422  // Loop over the edges
1423  for (unsigned j = 0; j < n_size; j++)
1424  {
1425  // Find pointer to neighbouring element along edge
1426  QuadTree* neigh_pt;
1427  neigh_pt = quadtree_pt()->gteq_edge_neighbour(edges[j],
1428  translate_s,
1429  s_lo_neigh,
1430  s_hi_neigh,
1431  neigh_edge,
1432  diff_level,
1433  in_neighbouring_tree);
1434 
1435  // Neighbour exists
1436  if (neigh_pt != 0)
1437  {
1438  // Have any of its vertex nodes been created yet?
1439  // (Must look in incomplete neighbours because after the
1440  // pre-build they may contain pointers to the required nodes. e.g.
1441  // h-refinement of neighbouring linear and quadratic elements)
1442  bool a_vertex_node_is_built = false;
1443  QElement<2, INITIAL_NNODE_1D>* neigh_obj_pt =
1444  dynamic_cast<QElement<2, INITIAL_NNODE_1D>*>(neigh_pt->object_pt());
1445  if (neigh_obj_pt == 0)
1446  {
1447  throw OomphLibError("Not a quad element!",
1448  "PRefineableQElement<2,INITIAL_NNODE_1D>::node_"
1449  "created_by_neighbour()",
1450  OOMPH_EXCEPTION_LOCATION);
1451  }
1452  for (unsigned vnode = 0; vnode < neigh_obj_pt->nvertex_node(); vnode++)
1453  {
1454  if (neigh_obj_pt->vertex_node_pt(vnode) != 0)
1455  a_vertex_node_is_built = true;
1456  break;
1457  }
1458  if (a_vertex_node_is_built)
1459  {
1460  // We now need to translate the nodal location
1461  // as defined in terms of the fractional coordinates of the present
1462  // element into those of its neighbour
1463 
1464  // Calculate the local coordinate in the neighbour
1465  // Note that we need to use the translation scheme in case
1466  // the local coordinates are swapped in the neighbour.
1467  for (unsigned i = 0; i < 2; i++)
1468  {
1469  s[i] = s_lo_neigh[i] +
1470  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
1471  }
1472 
1473  // Find the node in the neighbour
1474  Node* neighbour_node_pt =
1475  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
1476 
1477  // If there is a node, return it
1478  if (neighbour_node_pt != 0)
1479  {
1480  // Now work out whether it's a periodic boundary
1481  // only possible if we have moved into a neighbouring tree
1482  if (in_neighbouring_tree)
1483  {
1484  // Return whether the neighbour is periodic
1485  is_periodic =
1486  quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1487  }
1488  // Return the pointer to the neighbouring node
1489  return neighbour_node_pt;
1490  }
1491  }
1492  }
1493  }
1494  // Node not found, return null
1495  return 0;
1496  }

References Global_Physical_Variables::E, oomph::FiniteElement::get_node_at_local_coordinate(), oomph::QuadTree::gteq_edge_neighbour(), i, j, N, oomph::Tree::object_pt(), OOMPH_EXCEPTION_LOCATION, s, oomph::QuadTreeNames::S, and oomph::QuadTreeNames::W.

◆ node_created_by_son_of_neighbour()

template<unsigned INITIAL_NNODE_1D>
Node * oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::node_created_by_son_of_neighbour ( const Vector< double > &  s_fraction,
bool &  is_periodic 
)
virtual

If a neighbouring element's son has already created a node at a position corresponding to the local fractional position within the present element, s_fraction, return a pointer to that node. If not, return NULL (0). If the node is periodic the flag is_periodic will be true

Reimplemented from oomph::RefineableQElement< 2 >.

1509  {
1510  using namespace QuadTreeNames;
1511 
1512  // Calculate the edges on which the node lies
1513  Vector<int> edges;
1514  if (s_fraction[0] == 0.0)
1515  {
1516  edges.push_back(W);
1517  }
1518  if (s_fraction[0] == 1.0)
1519  {
1520  edges.push_back(E);
1521  }
1522  if (s_fraction[1] == 0.0)
1523  {
1524  edges.push_back(S);
1525  }
1526  if (s_fraction[1] == 1.0)
1527  {
1528  edges.push_back(N);
1529  }
1530 
1531  // Find the number of edges
1532  unsigned n_size = edges.size();
1533  // If there are no edges, then there is no neighbour, return 0
1534  if (n_size == 0)
1535  {
1536  return 0;
1537  }
1538 
1539  Vector<unsigned> translate_s(2);
1540  Vector<double> s_lo_neigh(2);
1541  Vector<double> s_hi_neigh(2);
1542  Vector<double> s(2);
1543 
1544  int neigh_edge, diff_level;
1545  bool in_neighbouring_tree;
1546 
1547  // Loop over the edges
1548  for (unsigned j = 0; j < n_size; j++)
1549  {
1550  // Find pointer to neighbouring element along edge
1551  QuadTree* neigh_pt;
1552  neigh_pt = quadtree_pt()->gteq_edge_neighbour(edges[j],
1553  translate_s,
1554  s_lo_neigh,
1555  s_hi_neigh,
1556  neigh_edge,
1557  diff_level,
1558  in_neighbouring_tree);
1559 
1560  // Neighbour exists
1561  if (neigh_pt != 0)
1562  {
1563  // Have its nodes been created yet?
1564  // (Must look in sons of unfinished neighbours too!!!)
1565  if (1)
1566  {
1567  // We now need to translate the nodal location
1568  // as defined in terms of the fractional coordinates of the present
1569  // element into those of its neighbour
1570 
1571  // Calculate the local coordinate in the neighbour
1572  // Note that we need to use the translation scheme in case
1573  // the local coordinates are swapped in the neighbour.
1574  for (unsigned i = 0; i < 2; i++)
1575  {
1576  s[i] = s_lo_neigh[i] +
1577  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
1578  }
1579 
1580  // Check if the element has sons
1581  if (neigh_pt->nsons() != 0)
1582  {
1583  // First, find the son element in which the node should live
1584 
1585  // Find coordinates in the sons
1586  Vector<double> s_in_son(2);
1587  int son = -10;
1588  // If negative on the west side
1589  if (s[0] < 0.0)
1590  {
1591  // On the south side
1592  if (s[1] < 0.0)
1593  {
1594  // It's the southwest son
1595  son = SW;
1596  s_in_son[0] = 1.0 + 2.0 * s[0];
1597  s_in_son[1] = 1.0 + 2.0 * s[1];
1598  }
1599  // On the north side
1600  else
1601  {
1602  // It's the northwest son
1603  son = NW;
1604  s_in_son[0] = 1.0 + 2.0 * s[0];
1605  s_in_son[1] = -1.0 + 2.0 * s[1];
1606  }
1607  }
1608  else
1609  {
1610  // On the south side
1611  if (s[1] < 0.0)
1612  {
1613  // It's the southeast son
1614  son = SE;
1615  s_in_son[0] = -1.0 + 2.0 * s[0];
1616  s_in_son[1] = 1.0 + 2.0 * s[1];
1617  }
1618  // On the north side
1619  else
1620  {
1621  // It's the northeast son
1622  son = NE;
1623  s_in_son[0] = -1.0 + 2.0 * s[0];
1624  s_in_son[1] = -1.0 + 2.0 * s[1];
1625  }
1626  }
1627 
1628  // Find the node in the neighbour's son
1629  Node* neighbour_son_node_pt =
1630  neigh_pt->son_pt(son)->object_pt()->get_node_at_local_coordinate(
1631  s_in_son);
1632 
1633  // If there is a node, return it
1634  if (neighbour_son_node_pt != 0)
1635  {
1636  // Now work out whether it's a periodic boundary
1637  // only possible if we have moved into a neighbouring tree
1638  if (in_neighbouring_tree)
1639  {
1640  // Return whether the neighbour is periodic
1641  is_periodic =
1642  quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
1643  }
1644  // Return the pointer to the neighbouring node
1645  return neighbour_son_node_pt;
1646  }
1647  }
1648  else
1649  {
1650  // No sons to search in, so no node can be found
1651  return 0;
1652  }
1653  }
1654  }
1655  }
1656  // Node not found, return null
1657  return 0;
1658  }
oomph::QuadTreeNames::SE
@ SE
Definition: quadtree.h:57
oomph::QuadTreeNames::NW
@ NW
Definition: quadtree.h:58
oomph::QuadTreeNames::NE
@ NE
Definition: quadtree.h:59
oomph::QuadTreeNames::SW
@ SW
Definition: quadtree.h:56

References Global_Physical_Variables::E, oomph::FiniteElement::get_node_at_local_coordinate(), oomph::QuadTree::gteq_edge_neighbour(), i, j, N, oomph::QuadTreeNames::NE, oomph::Tree::nsons(), oomph::QuadTreeNames::NW, oomph::Tree::object_pt(), s, oomph::QuadTreeNames::S, oomph::QuadTreeNames::SE, oomph::Tree::son_pt(), oomph::QuadTreeNames::SW, and oomph::QuadTreeNames::W.

◆ p_refine()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::p_refine ( const int &  inc,
Mesh *const &  mesh_pt,
GeneralisedElement *const &  clone_pt 
)
virtual

p-refine the element (refine if inc>0, unrefine if inc<0).

p-refine the element inc times. (If inc<0 then p-unrefinement is performed.)

/ Now loop over nodes in element

/BENFLAG:

/stdcout << "this = " << this << std::endl;

Implements oomph::PRefineableElement.

1992  {
1993  // Cast clone to correct type
1994  PRefineableQElement<2, INITIAL_NNODE_1D>* clone_el_pt =
1995  dynamic_cast<PRefineableQElement<2, INITIAL_NNODE_1D>*>(clone_pt);
1996 
1997  // If it is a MacroElement then we need to copy the geometric data too.
1998  MacroElementNodeUpdateElementBase* m_el_pt =
1999  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
2000  if (m_el_pt != 0)
2001  {
2002  MacroElementNodeUpdateElementBase* clone_m_el_pt =
2003  dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_pt);
2004 
2005  // Store local copy of geom object vector, so it can be passed on
2006  // to son elements (and their nodes) during refinement
2007  unsigned ngeom_object = m_el_pt->ngeom_object();
2008  clone_m_el_pt->geom_object_pt().resize(ngeom_object);
2009  for (unsigned i = 0; i < ngeom_object; i++)
2010  {
2011  clone_m_el_pt->geom_object_pt()[i] = m_el_pt->geom_object_pt(i);
2012  }
2013 
2014  // clone_m_el_pt->set_node_update_info(m_el_pt->geom_object_pt());
2015  }
2016 
2017  // Check if we can use it
2018  if (clone_el_pt == 0)
2019  {
2020  throw OomphLibError(
2021  "Cloned copy must be a PRefineableQElement<2,INITIAL_NNODE_1D>!",
2022  OOMPH_CURRENT_FUNCTION,
2023  OOMPH_EXCEPTION_LOCATION);
2024  }
2025 #ifdef PARANOID
2026  // Clone exists, so check that it is infact a copy of me
2027  else
2028  {
2029  // Flag to keep track of check
2030  bool clone_is_ok = true;
2031 
2032  // Does the clone have the correct p-order?
2033  clone_is_ok = clone_is_ok && (clone_el_pt->p_order() == this->p_order());
2034 
2035  if (!clone_is_ok)
2036  {
2037  std::ostringstream err_stream;
2038  err_stream << "Clone element has a different p-order from me,"
2039  << std::endl
2040  << "but it is supposed to be a copy!" << std::endl;
2041 
2042  throw OomphLibError(
2043  err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2044  }
2045 
2046  // Does the clone have the same number of nodes as me?
2047  clone_is_ok = clone_is_ok && (clone_el_pt->nnode() == this->nnode());
2048 
2049  if (!clone_is_ok)
2050  {
2051  std::ostringstream err_stream;
2052  err_stream << "Clone element has a different number of nodes from me,"
2053  << std::endl
2054  << "but it is supposed to be a copy!" << std::endl;
2055 
2056  throw OomphLibError(
2057  err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2058  }
2059 
2060  // Does the clone have the same nodes as me?
2061  for (unsigned n = 0; n < this->nnode(); n++)
2062  {
2063  clone_is_ok =
2064  clone_is_ok && (clone_el_pt->node_pt(n) == this->node_pt(n));
2065  }
2066 
2067  if (!clone_is_ok)
2068  {
2069  std::ostringstream err_stream;
2070  err_stream << "The nodes belonging to the clone element are different"
2071  << std::endl
2072  << "from mine, but it is supposed to be a copy!"
2073  << std::endl;
2074 
2075  throw OomphLibError(
2076  err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2077  }
2078 
2079  // Is this a macro element?
2080  MacroElementNodeUpdateElementBase* m_el_pt =
2081  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
2082  if (m_el_pt != 0)
2083  {
2084  // Get macro element version of clone
2085  MacroElementNodeUpdateElementBase* clone_m_el_pt =
2086  dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_el_pt);
2087 
2088  // Does the clone have the same geometric objects?
2089  Vector<GeomObject*> clone_geom_object_pt(
2090  clone_m_el_pt->geom_object_pt());
2091  Vector<GeomObject*> geom_object_pt(m_el_pt->geom_object_pt());
2092 
2093  clone_is_ok =
2094  clone_is_ok && (geom_object_pt.size() == clone_geom_object_pt.size());
2095  for (unsigned n = 0;
2096  n < std::min(geom_object_pt.size(), clone_geom_object_pt.size());
2097  n++)
2098  {
2099  clone_is_ok =
2100  clone_is_ok && (clone_geom_object_pt[n] == geom_object_pt[n]);
2101  }
2102 
2103  if (!clone_is_ok)
2104  {
2105  std::ostringstream err_stream;
2106  err_stream << "The clone element has different geometric objects"
2107  << std::endl
2108  << "from mine, but it is supposed to be a copy!"
2109  << std::endl;
2110 
2111  throw OomphLibError(
2112  err_stream.str(),
2113  "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2114  OOMPH_EXCEPTION_LOCATION);
2115  }
2116  }
2117 
2118  // If we get to here then the clone has all the information we require
2119  }
2120 #endif
2121 
2122  // Currently we can't handle the case of generalised coordinates
2123  // since we haven't established how they should be interpolated.
2124  // Buffer this case:
2125  if (clone_el_pt->node_pt(0)->nposition_type() != 1)
2126  {
2127  throw OomphLibError("Can't handle generalised nodal positions (yet).",
2128  OOMPH_CURRENT_FUNCTION,
2129  OOMPH_EXCEPTION_LOCATION);
2130  }
2131 
2132  // Timestepper should be the same for all nodes -- use it
2133  // to create timesteppers for new nodes
2134  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
2135 
2136  // Get number of history values (incl. present)
2137  unsigned ntstorage = time_stepper_pt->ntstorage();
2138 
2139  // Increment p-order of the element
2140  P_order += inc;
2141 
2142  // Change integration scheme
2143  delete this->integral_pt();
2144  switch (P_order)
2145  {
2146  case 2:
2147  this->set_integration_scheme(new GaussLobattoLegendre<2, 2>);
2148  break;
2149  case 3:
2150  this->set_integration_scheme(new GaussLobattoLegendre<2, 3>);
2151  break;
2152  case 4:
2153  this->set_integration_scheme(new GaussLobattoLegendre<2, 4>);
2154  break;
2155  case 5:
2156  this->set_integration_scheme(new GaussLobattoLegendre<2, 5>);
2157  break;
2158  case 6:
2159  this->set_integration_scheme(new GaussLobattoLegendre<2, 6>);
2160  break;
2161  case 7:
2162  this->set_integration_scheme(new GaussLobattoLegendre<2, 7>);
2163  break;
2164  default:
2165  std::ostringstream error_message;
2166  error_message << "\nERROR: Exceeded maximum polynomial order for";
2167  error_message << "\n integration scheme.\n";
2168  throw OomphLibError(error_message.str(),
2169  OOMPH_CURRENT_FUNCTION,
2170  OOMPH_EXCEPTION_LOCATION);
2171  }
2172 
2173  // Allocate new space for Nodes (at the element level)
2174  this->set_n_node(P_order * P_order);
2175 
2176  // Copy vertex nodes and create new edge and internal nodes
2177  //---------------------------------------------------------
2178 
2179  // Setup vertex coordinates in element:
2180  //-------------------------------------
2181  Vector<double> s_lo(2);
2182  Vector<double> s_hi(2);
2183  s_lo[0] = -1.0;
2184  s_hi[0] = 1.0;
2185  s_lo[1] = -1.0;
2186  s_hi[1] = 1.0;
2187 
2188  // Local coordinate in element
2189  Vector<double> s(2);
2190 
2191  unsigned jnod = 0;
2192 
2193  Vector<double> s_fraction(2);
2194  // Loop over nodes in element
2195  for (unsigned i0 = 0; i0 < P_order; i0++)
2196  {
2197  // Get the fractional position of the node in the direction of s[0]
2198  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
2199  // Local coordinate
2200  s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
2201 
2202  for (unsigned i1 = 0; i1 < P_order; i1++)
2203  {
2204  // Get the fractional position of the node in the direction of s[1]
2205  s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
2206  // Local coordinate
2207  s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
2208 
2209  // Local node number
2210  jnod = i0 + P_order * i1;
2211 
2212  // Initialise flag: So far, this node hasn't been built
2213  // or copied yet
2214  bool node_done = false;
2215 
2216  // Get the pointer to the node in this element (or rather, its clone),
2217  // returns NULL if there is not node
2218  Node* created_node_pt = clone_el_pt->get_node_at_local_coordinate(s);
2219 
2220  // Does this node already exist in this element?
2221  //----------------------------------------------
2222  if (created_node_pt != 0)
2223  {
2224  // Copy node across
2225  this->node_pt(jnod) = created_node_pt;
2226 
2227  // Make sure that we update the values at the node so that
2228  // they are consistent with the present representation.
2229  // This is only need for mixed interpolation where the value
2230  // at the father could now become active.
2231 
2232  // Loop over all history values
2233  for (unsigned t = 0; t < ntstorage; t++)
2234  {
2235  // Get values from father element
2236  // Note: get_interpolated_values() sets Vector size itself.
2237  Vector<double> prev_values;
2238  clone_el_pt->get_interpolated_values(t, s, prev_values);
2239  // Find the minimum number of values
2240  //(either those stored at the node, or those returned by
2241  // the function)
2242  unsigned n_val_at_node = created_node_pt->nvalue();
2243  unsigned n_val_from_function = prev_values.size();
2244  // Use the ternary conditional operator here
2245  unsigned n_var = n_val_at_node < n_val_from_function ?
2246  n_val_at_node :
2247  n_val_from_function;
2248  // Assign the values that we can
2249  for (unsigned k = 0; k < n_var; k++)
2250  {
2251  created_node_pt->set_value(t, k, prev_values[k]);
2252  }
2253  }
2254 
2255  // Node has been created by copy
2256  node_done = true;
2257  }
2258  // Node does not exist in this element but might already
2259  //------------------------------------------------------
2260  // have been created by neighbouring elements
2261  //-------------------------------------------
2262  else
2263  {
2264  // Was the node created by one of its neighbours
2265  // Whether or not the node lies on an edge can be calculated
2266  // by from the fractional position
2267  bool is_periodic = false;
2268  created_node_pt = node_created_by_neighbour(s_fraction, is_periodic);
2269 
2270  // If the node was so created, assign the pointers
2271  if (created_node_pt != 0)
2272  {
2273  // If the node is periodic
2274  if (is_periodic)
2275  {
2276  // Now the node must be on a boundary, but we don't know which
2277  // one
2278  // The returned created_node_pt is actually the neighbouring
2279  // periodic node
2280  Node* neighbour_node_pt = created_node_pt;
2281 
2282  // Determine the edge on which the new node will live
2283  //(cannot be a vertex node)
2284  int my_bound = Tree::OMEGA;
2285  if (s_fraction[0] == 0.0) my_bound = QuadTreeNames::W;
2286  else if (s_fraction[0] == 1.0)
2287  my_bound = QuadTreeNames::E;
2288  else if (s_fraction[1] == 0.0)
2289  my_bound = QuadTreeNames::S;
2290  else if (s_fraction[1] == 1.0)
2291  my_bound = QuadTreeNames::N;
2292 
2293  // Storage for the set of Mesh boundaries on which the
2294  // appropriate edge lives.
2295  // [New nodes should always be mid-edge nodes and therefore
2296  // only live on one boundary but just to play it safe...]
2297  std::set<unsigned> boundaries;
2298  // Only get the boundaries if we are at the edge of
2299  // an element. Nodes in the centre of an element cannot be
2300  // on Mesh boundaries
2301  if (my_bound != Tree::OMEGA)
2302  {
2303  clone_el_pt->get_boundaries(my_bound, boundaries);
2304  }
2305 
2306 #ifdef PARANOID
2307  // Case where a new node lives on more than one boundary
2308  // seems fishy enough to flag
2309  if (boundaries.size() > 1)
2310  {
2311  throw OomphLibError(
2312  "boundaries.size()!=1 seems a bit strange..\n",
2313  OOMPH_CURRENT_FUNCTION,
2314  OOMPH_EXCEPTION_LOCATION);
2315  }
2316 
2317  // Case when there are no boundaries, we are in big trouble
2318  if (boundaries.size() == 0)
2319  {
2320  std::ostringstream error_stream;
2321  error_stream << "Periodic node is not on a boundary...\n"
2322  << "Coordinates: " << created_node_pt->x(0) << " "
2323  << created_node_pt->x(1) << "\n";
2324  throw OomphLibError(error_stream.str(),
2325  OOMPH_CURRENT_FUNCTION,
2326  OOMPH_EXCEPTION_LOCATION);
2327  }
2328 #endif
2329 
2330  // Create node and set the pointer to it from the element
2331  created_node_pt =
2332  this->construct_boundary_node(jnod, time_stepper_pt);
2333  // Make the node periodic from the neighbour
2334  created_node_pt->make_periodic(neighbour_node_pt);
2335 
2336  // Loop over # of history values
2337  for (unsigned t = 0; t < ntstorage; t++)
2338  {
2339  // Get position from father element -- this uses the macro
2340  // element representation if appropriate. If the node
2341  // turns out to be a hanging node later on, then
2342  // its position gets adjusted in line with its
2343  // hanging node interpolation.
2344  Vector<double> x_prev(2);
2345  clone_el_pt->get_x(t, s, x_prev);
2346  // Set previous positions of the new node
2347  for (unsigned i = 0; i < 2; i++)
2348  {
2349  created_node_pt->x(t, i) = x_prev[i];
2350  }
2351  }
2352 
2353  // Check if we need to add nodes to the mesh
2354  if (mesh_pt != 0)
2355  {
2356  // Next, we Update the boundary lookup schemes
2357  // Loop over the boundaries stored in the set
2358  for (std::set<unsigned>::iterator it = boundaries.begin();
2359  it != boundaries.end();
2360  ++it)
2361  {
2362  // Add the node to the boundary
2363  mesh_pt->add_boundary_node(*it, created_node_pt);
2364 
2365  // If we have set an intrinsic coordinate on this
2366  // mesh boundary then it must also be interpolated on
2367  // the new node
2368  // Now interpolate the intrinsic boundary coordinate
2369  if (mesh_pt->boundary_coordinate_exists(*it) == true)
2370  {
2371  Vector<double> zeta(1);
2372  clone_el_pt->interpolated_zeta_on_edge(
2373  *it, my_bound, s, zeta);
2374 
2375  created_node_pt->set_coordinates_on_boundary(*it, zeta);
2376  }
2377  }
2378 
2379  // Make sure that we add the node to the mesh
2380  mesh_pt->add_node_pt(created_node_pt);
2381  }
2382  } // End of periodic case
2383  // Otherwise the node is not periodic, so just set the
2384  // pointer to the neighbours node
2385  else
2386  {
2387  this->node_pt(jnod) = created_node_pt;
2388  }
2389  // Node has been created
2390  node_done = true;
2391  }
2392  // Node does not exist in neighbour element but might already
2393  //-----------------------------------------------------------
2394  // have been created by a son of a neighbouring element
2395  //-----------------------------------------------------
2396  else
2397  {
2398  // Was the node created by one of its neighbours' sons
2399  // Whether or not the node lies on an edge can be calculated
2400  // by from the fractional position
2401  bool is_periodic = false;
2402  created_node_pt =
2403  node_created_by_son_of_neighbour(s_fraction, is_periodic);
2404 
2405  // If the node was so created, assign the pointers
2406  if (created_node_pt != 0)
2407  {
2408  // If the node is periodic
2409  if (is_periodic)
2410  {
2411  // Now the node must be on a boundary, but we don't know which
2412  // one
2413  // The returned created_node_pt is actually the neighbouring
2414  // periodic node
2415  Node* neighbour_node_pt = created_node_pt;
2416 
2417  // Determine the edge on which the new node will live
2418  //(cannot be a vertex node)
2419  int my_bound = Tree::OMEGA;
2420  if (s_fraction[0] == 0.0) my_bound = QuadTreeNames::W;
2421  else if (s_fraction[0] == 1.0)
2422  my_bound = QuadTreeNames::E;
2423  else if (s_fraction[1] == 0.0)
2424  my_bound = QuadTreeNames::S;
2425  else if (s_fraction[1] == 1.0)
2426  my_bound = QuadTreeNames::N;
2427 
2428  // Storage for the set of Mesh boundaries on which the
2429  // appropriate edge lives.
2430  // [New nodes should always be mid-edge nodes and therefore
2431  // only live on one boundary but just to play it safe...]
2432  std::set<unsigned> boundaries;
2433  // Only get the boundaries if we are at the edge of
2434  // an element. Nodes in the centre of an element cannot be
2435  // on Mesh boundaries
2436  if (my_bound != Tree::OMEGA)
2437  {
2438  clone_el_pt->get_boundaries(my_bound, boundaries);
2439  }
2440 
2441 #ifdef PARANOID
2442  // Case where a new node lives on more than one boundary
2443  // seems fishy enough to flag
2444  if (boundaries.size() > 1)
2445  {
2446  throw OomphLibError(
2447  "boundaries.size()!=1 seems a bit strange..\n",
2448  OOMPH_CURRENT_FUNCTION,
2449  OOMPH_EXCEPTION_LOCATION);
2450  }
2451 
2452  // Case when there are no boundaries, we are in big trouble
2453  if (boundaries.size() == 0)
2454  {
2455  std::ostringstream error_stream;
2456  error_stream << "Periodic node is not on a boundary...\n"
2457  << "Coordinates: " << created_node_pt->x(0)
2458  << " " << created_node_pt->x(1) << "\n";
2459  throw OomphLibError(error_stream.str(),
2460  OOMPH_CURRENT_FUNCTION,
2461  OOMPH_EXCEPTION_LOCATION);
2462  }
2463 #endif
2464 
2465  // Create node and set the pointer to it from the element
2466  created_node_pt =
2467  this->construct_boundary_node(jnod, time_stepper_pt);
2468  // Make the node periodic from the neighbour
2469  created_node_pt->make_periodic(neighbour_node_pt);
2470 
2471  // Loop over # of history values
2472  for (unsigned t = 0; t < ntstorage; t++)
2473  {
2474  // Get position from father element -- this uses the macro
2475  // element representation if appropriate. If the node
2476  // turns out to be a hanging node later on, then
2477  // its position gets adjusted in line with its
2478  // hanging node interpolation.
2479  Vector<double> x_prev(2);
2480  clone_el_pt->get_x(t, s, x_prev);
2481  // Set previous positions of the new node
2482  for (unsigned i = 0; i < 2; i++)
2483  {
2484  created_node_pt->x(t, i) = x_prev[i];
2485  }
2486  }
2487 
2488  // Check if we need to add nodes to the mesh
2489  if (mesh_pt != 0)
2490  {
2491  // Next, we Update the boundary lookup schemes
2492  // Loop over the boundaries stored in the set
2493  for (std::set<unsigned>::iterator it = boundaries.begin();
2494  it != boundaries.end();
2495  ++it)
2496  {
2497  // Add the node to the boundary
2498  mesh_pt->add_boundary_node(*it, created_node_pt);
2499 
2500  // If we have set an intrinsic coordinate on this
2501  // mesh boundary then it must also be interpolated on
2502  // the new node
2503  // Now interpolate the intrinsic boundary coordinate
2504  if (mesh_pt->boundary_coordinate_exists(*it) == true)
2505  {
2506  Vector<double> zeta(1);
2507  clone_el_pt->interpolated_zeta_on_edge(
2508  *it, my_bound, s, zeta);
2509 
2510  created_node_pt->set_coordinates_on_boundary(*it, zeta);
2511  }
2512  }
2513 
2514  // Make sure that we add the node to the mesh
2515  mesh_pt->add_node_pt(created_node_pt);
2516  }
2517  } // End of periodic case
2518  // Otherwise the node is not periodic, so just set the
2519  // pointer to the neighbours node
2520  else
2521  {
2522  this->node_pt(jnod) = created_node_pt;
2523  }
2524  // Node has been created
2525  node_done = true;
2526  } // Node does not exist in son of neighbouring element
2527  } // Node does not exist in neighbouring element
2528  } // Node does not exist in this element
2529 
2530  // Node has not been built anywhere ---> build it here
2531  if (!node_done)
2532  {
2533  // Firstly, we need to determine whether or not a node lies
2534  // on the boundary before building it, because
2535  // we actually assign a different type of node on boundaries.
2536 
2537  // Determine the edge on which the new node will live
2538  //(cannot be a vertex node)
2539  int my_bound = Tree::OMEGA;
2540  if (s_fraction[0] == 0.0) my_bound = QuadTreeNames::W;
2541  else if (s_fraction[0] == 1.0)
2542  my_bound = QuadTreeNames::E;
2543  else if (s_fraction[1] == 0.0)
2544  my_bound = QuadTreeNames::S;
2545  else if (s_fraction[1] == 1.0)
2546  my_bound = QuadTreeNames::N;
2547 
2548  // Storage for the set of Mesh boundaries on which the
2549  // appropriate edge lives.
2550  // [New nodes should always be mid-edge nodes and therefore
2551  // only live on one boundary but just to play it safe...]
2552  std::set<unsigned> boundaries;
2553  // Only get the boundaries if we are at the edge of
2554  // an element. Nodes in the centre of an element cannot be
2555  // on Mesh boundaries
2556  if (my_bound != Tree::OMEGA)
2557  {
2558  clone_el_pt->get_boundaries(my_bound, boundaries);
2559  }
2560 
2561 #ifdef PARANOID
2562  // Case where a new node lives on more than one boundary
2563  // seems fishy enough to flag
2564  if (boundaries.size() > 1)
2565  {
2566  throw OomphLibError("boundaries.size()!=1 seems a bit strange..\n",
2567  OOMPH_CURRENT_FUNCTION,
2568  OOMPH_EXCEPTION_LOCATION);
2569  }
2570 #endif
2571 
2572  // If the node lives on a mesh boundary,
2573  // then we need to create a boundary node
2574  if (boundaries.size() > 0)
2575  {
2576  // Create node and set the pointer to it from the element
2577  created_node_pt =
2578  this->construct_boundary_node(jnod, time_stepper_pt);
2579 
2580  // Now we need to work out whether to pin the values at
2581  // the new node based on the boundary conditions applied at
2582  // its Mesh boundary
2583 
2584  // Get the boundary conditions from the father
2585  Vector<int> bound_cons(this->ncont_interpolated_values());
2586  clone_el_pt->get_bcs(my_bound, bound_cons);
2587 
2588  // Loop over the values and pin, if necessary
2589  unsigned n_value = created_node_pt->nvalue();
2590  for (unsigned k = 0; k < n_value; k++)
2591  {
2592  if (bound_cons[k])
2593  {
2594  created_node_pt->pin(k);
2595  }
2596  }
2597 
2598  // Solid node? If so, deal with the positional boundary
2599  // conditions:
2600  SolidNode* solid_node_pt =
2601  dynamic_cast<SolidNode*>(created_node_pt);
2602  if (solid_node_pt != 0)
2603  {
2604  // Get the positional boundary conditions from the father:
2605  unsigned n_dim = created_node_pt->ndim();
2606  Vector<int> solid_bound_cons(n_dim);
2607  RefineableSolidQElement<2>* clone_solid_el_pt =
2608  dynamic_cast<RefineableSolidQElement<2>*>(clone_el_pt);
2609 #ifdef PARANOID
2610  if (clone_solid_el_pt == 0)
2611  {
2612  std::string error_message =
2613  "We have a SolidNode outside a refineable SolidElement\n";
2614  error_message +=
2615  "during mesh refinement -- this doesn't make sense";
2616 
2617  throw OomphLibError(error_message,
2618  OOMPH_CURRENT_FUNCTION,
2619  OOMPH_EXCEPTION_LOCATION);
2620  }
2621 #endif
2622  clone_solid_el_pt->get_solid_bcs(my_bound, solid_bound_cons);
2623 
2624  // Loop over the positions and pin, if necessary
2625  for (unsigned k = 0; k < n_dim; k++)
2626  {
2627  if (solid_bound_cons[k])
2628  {
2629  solid_node_pt->pin_position(k);
2630  }
2631  }
2632  } // End of if solid_node_pt
2633 
2634 
2635  // Check if we need to add nodes to the mesh
2636  if (mesh_pt != 0)
2637  {
2638  // Next, we Update the boundary lookup schemes
2639  // Loop over the boundaries stored in the set
2640  for (std::set<unsigned>::iterator it = boundaries.begin();
2641  it != boundaries.end();
2642  ++it)
2643  {
2644  // Add the node to the boundary
2645  mesh_pt->add_boundary_node(*it, created_node_pt);
2646 
2647  // If we have set an intrinsic coordinate on this
2648  // mesh boundary then it must also be interpolated on
2649  // the new node
2650  // Now interpolate the intrinsic boundary coordinate
2651  if (mesh_pt->boundary_coordinate_exists(*it) == true)
2652  {
2653  Vector<double> zeta(1);
2654  clone_el_pt->interpolated_zeta_on_edge(
2655  *it, my_bound, s, zeta);
2656 
2657  created_node_pt->set_coordinates_on_boundary(*it, zeta);
2658  }
2659  }
2660  }
2661  }
2662  // Otherwise the node is not on a Mesh boundary and
2663  // we create a normal "bulk" node
2664  else
2665  {
2666  // Create node and set the pointer to it from the element
2667  created_node_pt = this->construct_node(jnod, time_stepper_pt);
2668  }
2669 
2670  // Now we set the position and values at the newly created node
2671 
2672  // In the first instance use macro element or FE representation
2673  // to create past and present nodal positions.
2674  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
2675  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
2676  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
2677  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
2678  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
2679  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
2680  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
2681 
2682  // Loop over # of history values
2683  for (unsigned t = 0; t < ntstorage; t++)
2684  {
2685  // Get position from father element -- this uses the macro
2686  // element representation if appropriate. If the node
2687  // turns out to be a hanging node later on, then
2688  // its position gets adjusted in line with its
2689  // hanging node interpolation.
2690  Vector<double> x_prev(2);
2691  clone_el_pt->get_x(t, s, x_prev);
2692 
2693  // Set previous positions of the new node
2694  for (unsigned i = 0; i < 2; i++)
2695  {
2696  created_node_pt->x(t, i) = x_prev[i];
2697  }
2698  }
2699 
2700  // Loop over all history values
2701  for (unsigned t = 0; t < ntstorage; t++)
2702  {
2703  // Get values from father element
2704  // Note: get_interpolated_values() sets Vector size itself.
2705  Vector<double> prev_values;
2706  clone_el_pt->get_interpolated_values(t, s, prev_values);
2707  // Initialise the values at the new node
2708  unsigned n_value = created_node_pt->nvalue();
2709  for (unsigned k = 0; k < n_value; k++)
2710  {
2711  created_node_pt->set_value(t, k, prev_values[k]);
2712  }
2713  }
2714 
2715  // Add new node to mesh (if requested)
2716  if (mesh_pt != 0)
2717  {
2718  mesh_pt->add_node_pt(created_node_pt);
2719  }
2720 
2721  AlgebraicElementBase* alg_el_pt =
2722  dynamic_cast<AlgebraicElementBase*>(this);
2723 
2724  // If we do have an algebraic element
2725  if (alg_el_pt != 0)
2726  {
2727  std::string error_message = "Have not implemented p-refinement for";
2728  error_message += "Algebraic p-refineable elements yet\n";
2729 
2730  throw OomphLibError(
2731  error_message,
2732  "PRefineableQElement<2,INITIAL_NNODE_1D>::p_refine()",
2733  OOMPH_EXCEPTION_LOCATION);
2734  }
2735 
2736  } // End of case when we build the node ourselves
2737 
2738  // Check if the element is an algebraic element
2739  AlgebraicElementBase* alg_el_pt =
2740  dynamic_cast<AlgebraicElementBase*>(this);
2741 
2742  // If the element is an algebraic element, setup
2743  // node position (past and present) from algebraic node update
2744  // function. This over-writes previous assingments that
2745  // were made based on the macro-element/FE representation.
2746  // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
2747  // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
2748  // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
2749  // INFO FOR *ALL* ROOT ELEMENTS!
2750  if (alg_el_pt != 0)
2751  {
2752  // Build algebraic node update info for new node
2753  // This sets up the node update data for all node update
2754  // functions that are shared by all nodes in the father
2755  // element
2756  alg_el_pt->setup_algebraic_node_update(
2757  this->node_pt(jnod), s, clone_el_pt);
2758  }
2759 
2760  } // End of vertical loop over nodes in element
2761 
2762  } // End of horizontal loop over nodes in element
2763 
2764 
2765  // If the element is a MacroElementNodeUpdateElement, set
2766  // the update parameters for the current element's nodes --
2767  // all this needs is the vector of (pointers to the)
2768  // geometric objects that affect the MacroElement-based
2769  // node update -- this needs to be done to set the node
2770  // update info for newly created nodes
2771  MacroElementNodeUpdateElementBase* clone_m_el_pt =
2772  dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_el_pt);
2773  if (clone_m_el_pt != 0)
2774  {
2775  // Get vector of geometric objects from father (construct vector
2776  // via copy operation)
2777  Vector<GeomObject*> geom_object_pt(clone_m_el_pt->geom_object_pt());
2778 
2779  // Cast current element to MacroElementNodeUpdateElement:
2780  MacroElementNodeUpdateElementBase* m_el_pt =
2781  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
2782 
2783 #ifdef PARANOID
2784  if (m_el_pt == 0)
2785  {
2786  std::string error_message =
2787  "Failed to cast to MacroElementNodeUpdateElementBase*\n";
2788  error_message +=
2789  "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
2790  error_message += "then I should be too....\n";
2791 
2792  throw OomphLibError(
2793  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2794  }
2795 #endif
2796  // Build update info by passing vector of geometric objects:
2797  // This sets the current element to be the update element
2798  // for all of the element's nodes -- this is reversed
2799  // if the element is ever un-refined in the father element's
2800  // rebuild_from_sons() function which overwrites this
2801  // assignment to avoid nasty segmentation faults that occur
2802  // when a node tries to update itself via an element that no
2803  // longer exists...
2804  m_el_pt->set_node_update_info(geom_object_pt);
2805 
2807  // unsigned n_node = this->nnode();
2808  // for (unsigned j=0;j<n_node;j++)
2809  // {
2810  // // Get local coordinate in element (Vector sets its own size)
2811  // Vector<double> s_in_node_update_element;
2812  // this->local_coordinate_of_node(j,s_in_node_update_element);
2813  //
2814  // // Pass the lot to the node
2815  // static_cast<MacroElementNodeUpdateNode*>(this->node_pt(j))->
2816  // set_node_update_info(this,s_in_node_update_element,m_el_pt->geom_object_pt());
2817  // }
2818 
2820  // std::cout << "Checking that all the nodes have this as their update
2821  // element..." << std::endl;
2823  // for(unsigned j=0; j<this->nnode(); j++)
2824  // {
2825  // //std::cout << this->node_pt(j) << ": [" << this->node_pt(j)->x(0)
2826  // << "," << this->node_pt(j)->x(1) << "] update element: " <<
2827  // dynamic_cast<MacroElementNodeUpdateNode*>(this->node_pt(j))->node_update_element_pt()
2828  // << std::endl; MacroElementNodeUpdateNode* mac_nod_pt =
2829  // dynamic_cast<MacroElementNodeUpdateNode*>(this->node_pt(j));
2830  // if(mac_nod_pt->node_update_element_pt()!=this)
2831  // {
2832  // std::cout << "Something's not right! Update element is wrong..." <<
2833  // std::endl;
2834  // }
2835  // FiniteElement* up_el_pt =
2836  // dynamic_cast<FiniteElement*>(mac_nod_pt->node_update_element_pt());
2837  // bool not_good = true;
2838  // for(unsigned l=0; l<up_el_pt->nnode(); l++)
2839  // {
2840  // if(up_el_pt->node_pt(l)==mac_nod_pt)
2841  // {
2842  // not_good = false;
2843  // break;
2844  // }
2845  // }
2846  // if(not_good==true)
2847  // {
2848  // std::cout << "Macro node doesn't belong to its update element!" <<
2849  // std::endl;
2850  // }
2851  // }
2852 
2853 
2854  // Loop over all nodes in element again, to re-set the positions
2855  // This must be done using the new element's macro-element
2856  // representation, rather than the old version which may be
2857  // of a different p-order!
2858  for (unsigned i0 = 0; i0 < P_order; i0++)
2859  {
2860  // Get the fractional position of the node in the direction of s[0]
2861  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
2862  // Local coordinate
2863  s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
2864 
2865  for (unsigned i1 = 0; i1 < P_order; i1++)
2866  {
2867  // Get the fractional position of the node in the direction of s[1]
2868  s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
2869  // Local coordinate
2870  s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
2871 
2872  // Local node number
2873  jnod = i0 + P_order * i1;
2874 
2875  // Loop over # of history values
2876  for (unsigned t = 0; t < ntstorage; t++)
2877  {
2878  // Get position from father element -- this uses the macro
2879  // element representation if appropriate. If the node
2880  // turns out to be a hanging node later on, then
2881  // its position gets adjusted in line with its
2882  // hanging node interpolation.
2883  Vector<double> x_prev(2);
2884  this->get_x(t, s, x_prev);
2885 
2886  // Set previous positions of the new node
2887  for (unsigned i = 0; i < 2; i++)
2888  {
2889  this->node_pt(jnod)->x(t, i) = x_prev[i];
2890  }
2891  }
2892  }
2893  }
2894  }
2895 
2896  // Not necessary to delete the old nodes since all original nodes are in the
2897  // current mesh and so will be pruned as part of the mesh adaption process.
2898 
2899  // Do any further-build required
2900  this->further_build();
2901  }
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::Data::nvalue
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
oomph::FiniteElement::construct_node
virtual Node * construct_node(const unsigned &n)
Definition: elements.h:2509
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::FiniteElement::get_x
void get_x(const Vector< double > &s, Vector< double > &x) const
Definition: elements.h:1885
oomph::FiniteElement::construct_boundary_node
virtual Node * construct_boundary_node(const unsigned &n)
Definition: elements.h:2538
oomph::GeomObject::time_stepper_pt
TimeStepper *& time_stepper_pt()
Definition: geom_objects.h:192
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::node_created_by_neighbour
Node * node_created_by_neighbour(const Vector< double > &s_fraction, bool &is_periodic)
Definition: hp_refineable_elements.cc:1382
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::node_created_by_son_of_neighbour
Node * node_created_by_son_of_neighbour(const Vector< double > &s_fraction, bool &is_periodic)
Definition: hp_refineable_elements.cc:1507
oomph::RefineableElement::further_build
virtual void further_build()
Further build: e.g. deal with interpolation of internal values.
Definition: refineable_elements.h:599
oomph::RefineableElement::ncont_interpolated_values
virtual unsigned ncont_interpolated_values() const =0
oomph::TimeStepper::ntstorage
unsigned ntstorage() const
Definition: timesteppers.h:601
oomph::Tree::OMEGA
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:262
min
#define min(a, b)
Definition: datatypes.h:22
k
char char char int int * k
Definition: level2_impl.h:374
Eigen::zeta
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
oomph::QuadTreeNames::E
@ E
Definition: quadtree.h:61
oomph::QuadTreeNames::N
@ N
Definition: quadtree.h:60

References oomph::Mesh::add_boundary_node(), oomph::Mesh::add_node_pt(), oomph::Mesh::boundary_coordinate_exists(), oomph::QuadTreeNames::E, oomph::MacroElementNodeUpdateElementBase::geom_object_pt(), oomph::RefineableQElement< 2 >::get_bcs(), oomph::RefineableQElement< 2 >::get_boundaries(), oomph::RefineableElement::get_interpolated_values(), get_node_at_local_coordinate(), oomph::RefineableSolidQElement< 2 >::get_solid_bcs(), oomph::FiniteElement::get_x(), i, oomph::RefineableQElement< 2 >::interpolated_zeta_on_edge(), k, oomph::Node::make_periodic(), min, n, oomph::QuadTreeNames::N, oomph::Node::ndim(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::Node::nposition_type(), oomph::TimeStepper::ntstorage(), oomph::Data::nvalue(), oomph::Tree::OMEGA, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::PRefineableElement::p_order(), oomph::Data::pin(), oomph::SolidNode::pin_position(), s, oomph::QuadTreeNames::S, oomph::Node::set_coordinates_on_boundary(), oomph::MacroElementNodeUpdateElementBase::set_node_update_info(), oomph::Data::set_value(), oomph::AlgebraicElementBase::setup_algebraic_node_update(), oomph::Global_string_for_annotation::string(), plotPSD::t, oomph::QuadTreeNames::W, oomph::Node::x(), and Eigen::zeta().

◆ pre_build()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::pre_build ( Mesh *&  mesh_pt,
Vector< Node * > &  new_node_pt 
)
virtual

Pre-build (search father for required nodes which may already exist)

Check the father element for any required nodes which already exist

/ Build update info by passing vector of geometric objects: / This sets the current element to be the update element / for all of the element's nodes – this is reversed / if the element is ever un-refined in the father element's / rebuild_from_sons() function which overwrites this / assignment to avoid nasty segmentation faults that occur / when a node tries to update itself via an element that no / longer exists...

Reimplemented from oomph::RefineableElement.

1759  {
1760  using namespace QuadTreeNames;
1761 
1762  // Get the number of 1d nodes
1763  unsigned n_p = nnode_1d();
1764 
1765  // Check whether static father_bound needs to be created
1766  if (Father_bound[n_p].nrow() == 0)
1767  {
1768  setup_father_bounds();
1769  }
1770 
1771  // Pointer to my father (in quadtree impersonation)
1772  QuadTree* father_pt = dynamic_cast<QuadTree*>(quadtree_pt()->father_pt());
1773 
1774  // What type of son am I? Ask my quadtree representation...
1775  int son_type = Tree_pt->son_type();
1776 
1777  // Has somebody build me already? (If any nodes have not been built)
1778  if (!nodes_built())
1779  {
1780 #ifdef PARANOID
1781  if (father_pt == 0)
1782  {
1783  std::string error_message =
1784  "Something fishy here: I have no father and yet \n";
1785  error_message += "I have no nodes. Who has created me then?!\n";
1786 
1787  throw OomphLibError(
1788  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1789  }
1790 #endif
1791 
1792  // Return pointer to father element
1793  RefineableQElement<2>* father_el_pt =
1794  dynamic_cast<RefineableQElement<2>*>(father_pt->object_pt());
1795 
1796  // Timestepper should be the same for all nodes in father
1797  // element -- use it create timesteppers for new nodes
1798  TimeStepper* time_stepper_pt =
1799  father_el_pt->node_pt(0)->time_stepper_pt();
1800 
1801  // Number of history values (incl. present)
1802  unsigned ntstorage = time_stepper_pt->ntstorage();
1803 
1804  Vector<double> s_lo(2);
1805  Vector<double> s_hi(2);
1806  Vector<double> s(2);
1807  Vector<double> x(2);
1808 
1809  // Setup vertex coordinates in father element:
1810  //--------------------------------------------
1811  switch (son_type)
1812  {
1813  case SW:
1814  s_lo[0] = -1.0;
1815  s_hi[0] = 0.0;
1816  s_lo[1] = -1.0;
1817  s_hi[1] = 0.0;
1818  break;
1819 
1820  case SE:
1821  s_lo[0] = 0.0;
1822  s_hi[0] = 1.0;
1823  s_lo[1] = -1.0;
1824  s_hi[1] = 0.0;
1825  break;
1826 
1827  case NE:
1828  s_lo[0] = 0.0;
1829  s_hi[0] = 1.0;
1830  s_lo[1] = 0.0;
1831  s_hi[1] = 1.0;
1832  break;
1833 
1834  case NW:
1835  s_lo[0] = -1.0;
1836  s_hi[0] = 0.0;
1837  s_lo[1] = 0.0;
1838  s_hi[1] = 1.0;
1839  break;
1840  }
1841 
1842  // Pass macro element pointer on to sons and
1843  // set coordinates in macro element
1844  // hierher why can I see this?
1845  if (father_el_pt->macro_elem_pt() != 0)
1846  {
1847  set_macro_elem_pt(father_el_pt->macro_elem_pt());
1848  for (unsigned i = 0; i < 2; i++)
1849  {
1850  s_macro_ll(i) =
1851  father_el_pt->s_macro_ll(i) +
1852  0.5 * (s_lo[i] + 1.0) *
1853  (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1854  s_macro_ur(i) =
1855  father_el_pt->s_macro_ll(i) +
1856  0.5 * (s_hi[i] + 1.0) *
1857  (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1858  }
1859  }
1860 
1861 
1862  // If the father element hasn't been generated yet, we're stuck...
1863  if (father_el_pt->node_pt(0) == 0)
1864  {
1865  throw OomphLibError(
1866  "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
1867  OOMPH_CURRENT_FUNCTION,
1868  OOMPH_EXCEPTION_LOCATION);
1869  }
1870  else
1871  {
1872  unsigned jnod = 0;
1873  Vector<double> x_small(2);
1874  Vector<double> x_large(2);
1875 
1876  Vector<double> s_fraction(2);
1877  // Loop over nodes in element
1878  for (unsigned i0 = 0; i0 < n_p; i0++)
1879  {
1880  // Get the fractional position of the node in the direction of s[0]
1881  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
1882  // Local coordinate in father element
1883  s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
1884 
1885  for (unsigned i1 = 0; i1 < n_p; i1++)
1886  {
1887  // Get the fractional position of the node in the direction of s[1]
1888  s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
1889  // Local coordinate in father element
1890  s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
1891 
1892  // Local node number
1893  jnod = i0 + n_p * i1;
1894 
1895  // Check whether the father's node is periodic if so, complain
1896  /* {
1897  Node* father_node_pt = father_el_pt->node_pt(jnod);
1898  if((father_node_pt->is_a_copy()) ||
1899  (father_node_pt->position_is_a_copy()))
1900  {
1901  throw OomphLibError(
1902  "Can't handle periodic nodes (yet).",
1903  OOMPH_CURRENT_FUNCTION,
1904  OOMPH_EXCEPTION_LOCATION);
1905  }
1906  }*/
1907 
1908  // Initialise flag: So far, this node hasn't been built
1909  // or copied yet
1910  // bool node_done=false;
1911 
1912  // Get the pointer to the node in the father, returns NULL
1913  // if there is not node
1914  Node* created_node_pt =
1915  father_el_pt->get_node_at_local_coordinate(s);
1916 
1917  // Does this node already exist in father element?
1918  //------------------------------------------------
1919  if (created_node_pt != 0)
1920  {
1921  // Copy node across
1922  this->node_pt(jnod) = created_node_pt;
1923 
1924  // Make sure that we update the values at the node so that
1925  // they are consistent with the present representation.
1926  // This is only need for mixed interpolation where the value
1927  // at the father could now become active.
1928 
1929  // Loop over all history values
1930  for (unsigned t = 0; t < ntstorage; t++)
1931  {
1932  // Get values from father element
1933  // Note: get_interpolated_values() sets Vector size itself.
1934  Vector<double> prev_values;
1935  father_el_pt->get_interpolated_values(t, s, prev_values);
1936  // Find the minimum number of values
1937  //(either those stored at the node, or those returned by
1938  // the function)
1939  unsigned n_val_at_node = created_node_pt->nvalue();
1940  unsigned n_val_from_function = prev_values.size();
1941  // Use the ternary conditional operator here
1942  unsigned n_var = n_val_at_node < n_val_from_function ?
1943  n_val_at_node :
1944  n_val_from_function;
1945  // Assign the values that we can
1946  for (unsigned k = 0; k < n_var; k++)
1947  {
1948  created_node_pt->set_value(t, k, prev_values[k]);
1949  }
1950  }
1951 
1952  // Node has been created by copy
1953  // node_done=true;
1954  }
1955  } // End of vertical loop over nodes in element
1956  } // End of horizontal loop over nodes in element
1957  } // Sanity check: Father element has been generated
1958 
1959  } // End of nodes not built
1960  else
1961  {
1962  // If this is element updates by macro element node update then we need
1963  // to set the correct info in the nodes here since this won't get done
1964  // later in build() because we already have all our nodes created.
1965  MacroElementNodeUpdateElementBase* m_el_pt =
1966  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
1967  if (m_el_pt != 0)
1968  {
1969  // Get vector of geometric objects
1970  Vector<GeomObject*> geom_object_pt = m_el_pt->geom_object_pt();
1971 
1980  m_el_pt->set_node_update_info(geom_object_pt);
1981  }
1982  }
1983  }
oomph::FiniteElement::set_macro_elem_pt
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Definition: elements.h:1872
oomph::PRefineableElement::nodes_built
bool nodes_built()
Return true if all the nodes have been built, false if not.
Definition: refineable_elements.h:766
oomph::QElementBase::s_macro_ll
double & s_macro_ll(const unsigned &i)
Definition: Qelements.h:186
oomph::QElementBase::s_macro_ur
double & s_macro_ur(const unsigned &i)
Definition: Qelements.h:202
oomph::RefineableQElement< 2 >::Father_bound
static std::map< unsigned, DenseMatrix< int > > Father_bound
Definition: refineable_quad_element.h:176
oomph::RefineableQElement< 2 >::setup_father_bounds
void setup_father_bounds()
Definition: refineable_quad_element.cc:46
oomph::Tree::son_type
int son_type() const
Return son type.
Definition: tree.h:214
plotDoE.x
list x
Definition: plotDoE.py:28

References oomph::MacroElementNodeUpdateElementBase::geom_object_pt(), oomph::RefineableElement::get_interpolated_values(), oomph::FiniteElement::get_node_at_local_coordinate(), i, k, oomph::FiniteElement::macro_elem_pt(), oomph::QuadTreeNames::NE, oomph::FiniteElement::node_pt(), oomph::Data::nvalue(), oomph::QuadTreeNames::NW, oomph::Tree::object_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, oomph::QElementBase::s_macro_ll(), oomph::QElementBase::s_macro_ur(), oomph::QuadTreeNames::SE, oomph::MacroElementNodeUpdateElementBase::set_node_update_info(), oomph::Data::set_value(), oomph::Tree::son_type(), oomph::Global_string_for_annotation::string(), oomph::QuadTreeNames::SW, plotPSD::t, oomph::Data::time_stepper_pt(), and plotDoE::x.

◆ quad_hang_helper()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::quad_hang_helper ( const int &  value_id,
const int &  my_edge,
std::ofstream &  output_hangfile 
)
protectedvirtual

Set up hanging node information. Overloaded to implement the mortar method rather than constrained approximation. This enforces continuity weakly via an integral matching condition at non-conforming element boundaries.

Internal function to set up the hanging nodes on a particular edge of the element. Implements the mortarting method to enforce continuity weakly across non-conforming element boundaries \(\Gamma\) using an integral matching condition

\[ \int_\Gamma (u_{\mbox{S}} - u_{\mbox{M}}) \psi \mbox{d} s = 0 \]

for all polynomials \(\psi\) on \(\Gamma\) of degree at most p-2 (where p is the spectral-order of the dependent element) and a vertex matching condition

\[ (u_{\mbox{S}} - u_{\mbox{M}})\big\vert_{\partial\Gamma} = 0.\]

The algorithm works as follows:

  • First the element determines if its edge my_edge is on the master or dependent side of the non-conformity. At h-type non-conformities we choose long edges to be masters, and at p-type nonconformities the edge with lower p-order is the master.
  • Mortaring is performed by the dependent side.
  • If a vertex node of the mortar is shared between master and dependent element then the vertex matching condition is enforced automatically. Otherwise it must be imposed by constraining its value to that at on the master side.
  • The integral matching condition is discretised and the mortar test functions \( \psi \) are chosen to be derivatives of Legendre polynomials of degree p-1.
  • The mortar mass matrix M is constructed. Its entries are the mortar test functions evaluated at the dependent nodal positions, so it is diagonal.
  • The local projection matrix is constructed for the master element by applying the discretised integral matching condition along the mortar using the appropriate quadrature order.
  • The global projection matrix is then assembled by subtracting contributions from the mortar vertex nodes.
  • The mortar system \( M\xi^s = P\hat{\xi^m} \) is constructed, where \( \xi^m \) and \( \xi^s \) are the nodal values at the master and dependent nodes respectively.
  • The conformity matrix \( C = M^{-1}P \) is computed. This is straightforward since the mass matrix is diagonal.
  • Finally, the master nodes and weights for each dependent node are read from the conformity matrix and stored in the dependents' hanging schemes.

The positions of the dependent nodes are set to be consistent with their hanging schemes.

Reimplemented from oomph::RefineableQElement< 2 >.

4053  {
4054  using namespace QuadTreeNames;
4055 
4056  Vector<unsigned> translate_s(2);
4057  Vector<double> s_lo_neigh(2);
4058  Vector<double> s_hi_neigh(2);
4059  int neigh_edge, diff_level;
4060  bool in_neighbouring_tree;
4061 
4062  // Find pointer to neighbour in this direction
4063  QuadTree* neigh_pt;
4064  neigh_pt = this->quadtree_pt()->gteq_edge_neighbour(my_edge,
4065  translate_s,
4066  s_lo_neigh,
4067  s_hi_neigh,
4068  neigh_edge,
4069  diff_level,
4070  in_neighbouring_tree);
4071 
4072  // Work out master/dependent edges
4073  //----------------------------
4074 
4075  // Set up booleans
4076  // bool h_type_master = false;
4077  bool h_type_dependent = false;
4078  // bool p_type_master = false;
4079  bool p_type_dependent = false;
4080 
4081  // Neighbour exists and all nodes have been created
4082  if (neigh_pt != 0)
4083  {
4084  // Check if neighbour is bigger than me
4085  if (diff_level != 0)
4086  {
4087  // Dependent at h-type non-conformity
4088  h_type_dependent = true;
4089  }
4090  // Check if neighbour is the same size as me
4091  else if (neigh_pt->nsons() == 0)
4092  {
4093  // Neighbour is same size as me
4094  // Find p-orders of each element
4095  unsigned my_p_order =
4096  dynamic_cast<PRefineableQElement<2, INITIAL_NNODE_1D>*>(this)
4097  ->p_order();
4098  unsigned neigh_p_order =
4099  dynamic_cast<PRefineableQElement<2, INITIAL_NNODE_1D>*>(
4100  neigh_pt->object_pt())
4101  ->p_order();
4102 
4103  // Check for p-type non-conformity
4104  if (neigh_p_order == my_p_order)
4105  {
4106  // At a conforming interface
4107  }
4108  else if (neigh_p_order < my_p_order)
4109  {
4110  // Dependent at p-type non-conformity
4111  p_type_dependent = true;
4112  }
4113  else
4114  {
4115  // Master at p-type non-conformity
4116  // p_type_master = true;
4117  }
4118  }
4119  // Neighbour must be smaller than me
4120  else
4121  {
4122  // Master at h-type non-conformity
4123  // h_type_master = true;
4124  }
4125  }
4126  else
4127  {
4128  // Edge is on a boundary
4129  }
4130 
4131  // Now do the mortaring
4132  //---------------------
4133  if (h_type_dependent || p_type_dependent)
4134  {
4135  // Compute the active coordinate index along the this side of mortar
4136  unsigned active_coord_index;
4137  if (my_edge == N || my_edge == S) active_coord_index = 0;
4138  else if (my_edge == E || my_edge == W)
4139  active_coord_index = 1;
4140  else
4141  {
4142  throw OomphLibError("Cannot transform coordinates",
4143  OOMPH_CURRENT_FUNCTION,
4144  OOMPH_EXCEPTION_LOCATION);
4145  }
4146 
4147  // Get pointer to neighbouring master element (in p-refineable form)
4148  PRefineableQElement<2, INITIAL_NNODE_1D>* neigh_obj_pt;
4149  neigh_obj_pt = dynamic_cast<PRefineableQElement<2, INITIAL_NNODE_1D>*>(
4150  neigh_pt->object_pt());
4151 
4152  // Create vector of master and dependent nodes
4153  //----------------------------------------
4154  Vector<Node*> master_node_pt, dependent_node_pt;
4155  Vector<unsigned> master_node_number, dependent_node_number;
4156  Vector<Vector<double>> dependent_node_s_fraction;
4157 
4158  // Number of nodes in one dimension
4159  const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
4160  const unsigned neigh_n_p = neigh_obj_pt->ninterpolating_node_1d(value_id);
4161 
4162  // Test for the periodic node case
4163  // Are we crossing a periodic boundary
4164  bool is_periodic = false;
4165  if (in_neighbouring_tree)
4166  {
4167  is_periodic =
4168  this->tree_pt()->root_pt()->is_neighbour_periodic(my_edge);
4169  }
4170 
4171  // If it is periodic we actually need to get the node in
4172  // the neighbour of the neighbour (which will be a parent of
4173  // the present element) so that the "fixed" coordinate is
4174  // correctly calculated.
4175  // The idea is to replace the neigh_pt and associated data
4176  // with those of the neighbour of the neighbour
4177  if (is_periodic)
4178  {
4179  throw OomphLibError(
4180  "Cannot do mortaring with periodic hanging nodes yet! (Fix me!)",
4181  "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4182  OOMPH_EXCEPTION_LOCATION);
4183  } // End of special treatment for periodic hanging nodes
4184 
4185  // Storage for pointers to the nodes and their numbers along the master
4186  // edge
4187  unsigned neighbour_node_number = 0;
4188  Node* neighbour_node_pt = 0;
4189 
4190  // Loop over nodes along the edge
4191  for (unsigned i0 = 0; i0 < neigh_n_p; i0++)
4192  {
4193  // Find the neighbour's node
4194  switch (neigh_edge)
4195  {
4196  case N:
4197  neighbour_node_number = i0 + neigh_n_p * (neigh_n_p - 1);
4198  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
4199  neighbour_node_number, value_id);
4200  break;
4201 
4202  case S:
4203  neighbour_node_number = i0;
4204  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
4205  neighbour_node_number, value_id);
4206  break;
4207 
4208  case E:
4209  neighbour_node_number = neigh_n_p - 1 + neigh_n_p * i0;
4210  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
4211  neighbour_node_number, value_id);
4212  break;
4213 
4214  case W:
4215  neighbour_node_number = neigh_n_p * i0;
4216  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
4217  neighbour_node_number, value_id);
4218  break;
4219 
4220  default:
4221  throw OomphLibError("my_edge not N, S, W, E\n",
4222  OOMPH_CURRENT_FUNCTION,
4223  OOMPH_EXCEPTION_LOCATION);
4224  }
4225 
4226  // Set node as master node
4227  master_node_number.push_back(neighbour_node_number);
4228  master_node_pt.push_back(neighbour_node_pt);
4229  }
4230 
4231  // Storage for pointers to the local nodes and their numbers along my edge
4232  unsigned local_node_number = 0;
4233  Node* local_node_pt = 0;
4234 
4235  // Loop over the nodes along my edge
4236  for (unsigned i0 = 0; i0 < my_n_p; i0++)
4237  {
4238  // Storage for the fractional position of the node in the element
4239  Vector<double> s_fraction(2);
4240 
4241  // Find the local node and the fractional position of the node
4242  // which depends on the edge, of course
4243  switch (my_edge)
4244  {
4245  case N:
4246  s_fraction[0] =
4247  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
4248  s_fraction[1] = 1.0;
4249  local_node_number = i0 + my_n_p * (my_n_p - 1);
4250  local_node_pt =
4251  this->interpolating_node_pt(local_node_number, value_id);
4252  break;
4253 
4254  case S:
4255  s_fraction[0] =
4256  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
4257  s_fraction[1] = 0.0;
4258  local_node_number = i0;
4259  local_node_pt =
4260  this->interpolating_node_pt(local_node_number, value_id);
4261  break;
4262 
4263  case E:
4264  s_fraction[0] = 1.0;
4265  s_fraction[1] =
4266  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
4267  local_node_number = my_n_p - 1 + my_n_p * i0;
4268  local_node_pt =
4269  this->interpolating_node_pt(local_node_number, value_id);
4270  break;
4271 
4272  case W:
4273  s_fraction[0] = 0.0;
4274  s_fraction[1] =
4275  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
4276  local_node_number = my_n_p * i0;
4277  local_node_pt =
4278  this->interpolating_node_pt(local_node_number, value_id);
4279  break;
4280 
4281  default:
4282  throw OomphLibError("my_edge not N, S, W, E\n",
4283  OOMPH_CURRENT_FUNCTION,
4284  OOMPH_EXCEPTION_LOCATION);
4285  }
4286 
4287  // Add node to vector of dependent element nodes
4288  dependent_node_number.push_back(local_node_number);
4289  dependent_node_pt.push_back(local_node_pt);
4290 
4291  // Store node's local fraction
4292  dependent_node_s_fraction.push_back(s_fraction);
4293  }
4294 
4295  // Store the number of dependent and master nodes for use later
4296  const unsigned n_dependent_nodes = dependent_node_pt.size();
4297  const unsigned n_master_nodes = master_node_pt.size();
4298  const unsigned dependent_element_nnode_1d = my_n_p;
4299  const unsigned master_element_nnode_1d = neigh_n_p;
4300 
4301  // Storage for master shapes
4302  Shape master_shapes(neigh_obj_pt->ninterpolating_node(value_id));
4303 
4304  // Get master and dependent nodal positions
4305  //-------------------------------------
4306  Vector<double> dependent_nodal_position;
4307  Vector<double> dependent_weight(dependent_element_nnode_1d);
4308  Orthpoly::gll_nodes(
4309  dependent_element_nnode_1d, dependent_nodal_position, dependent_weight);
4310  Vector<double> master_nodal_position;
4311  Vector<double> master_weight(master_element_nnode_1d);
4312  Orthpoly::gll_nodes(
4313  master_element_nnode_1d, master_nodal_position, master_weight);
4314 
4315  // Apply the vertex matching condition
4316  //------------------------------------
4317  // Vertiex matching is ensured automatically in cases where there is a
4318  // node at each end of the mortar that is shared between the master and
4319  // dependent elements. Where this is not the case, the vertex matching
4320  // condition must be enforced by constraining the value of the unknown at
4321  // the node on the dependent side to be the same as the value at that
4322  // location in the master.
4323 
4324  // Store positions of the mortar vertex/non-vertex nodes in the dependent
4325  // element
4326  const unsigned n_mortar_vertices = 2;
4327  Vector<unsigned> vertex_pos(n_mortar_vertices);
4328  vertex_pos[0] = 0;
4329  vertex_pos[1] = this->ninterpolating_node_1d(value_id) - 1;
4330  Vector<unsigned> non_vertex_pos(my_n_p - n_mortar_vertices);
4331  for (unsigned i = 0; i < my_n_p - n_mortar_vertices; i++)
4332  {
4333  non_vertex_pos[i] = i + 1;
4334  }
4335 
4336  // Check if the mortar vertices are shared
4337  for (unsigned v = 0; v < n_mortar_vertices; v++)
4338  {
4339  // Search master node storage for the node
4340  bool node_is_shared = false;
4341  for (unsigned i = 0; i < master_node_pt.size(); i++)
4342  {
4343  if (dependent_node_pt[vertex_pos[v]] == master_node_pt[i])
4344  {
4345  node_is_shared = true;
4346  break;
4347  }
4348  }
4349 
4350  // If the node is not shared then we must constrain its value by setting
4351  // up a hanging scheme
4352  if (!node_is_shared)
4353  {
4354  // Calculate weights. These are just the master shapes evaluated at
4355  // this dependent node's position
4356 
4357  // Work out this node's location in the master
4358  Vector<double> s_in_neigh(2);
4359  for (unsigned i = 0; i < 2; i++)
4360  {
4361  s_in_neigh[i] =
4362  s_lo_neigh[i] +
4363  dependent_node_s_fraction[vertex_pos[v]][translate_s[i]] *
4364  (s_hi_neigh[i] - s_lo_neigh[i]);
4365  }
4366 
4367  // Get master shapes at dependent nodal positions
4368  neigh_obj_pt->interpolating_basis(
4369  s_in_neigh, master_shapes, value_id);
4370 
4371  // Remove any existing hanging node info
4372  // (This may be a bit wasteful, but guarantees correctness)
4373  dependent_node_pt[vertex_pos[v]]->set_nonhanging();
4374 
4375  // Set up hanging scheme for this node
4376  HangInfo* explicit_hang_pt = new HangInfo(n_master_nodes);
4377  for (unsigned m = 0; m < n_master_nodes; m++)
4378  {
4379  explicit_hang_pt->set_master_node_pt(
4380  m, master_node_pt[m], master_shapes[master_node_number[m]]);
4381  }
4382 
4383  // Make node hang
4384  dependent_node_pt[vertex_pos[v]]->set_hanging_pt(explicit_hang_pt,
4385  -1);
4386  }
4387  }
4388 
4389  // Check that there are actually dependent nodes for which we still need
4390  // to construct a hanging scheme. If not then there is nothing more to do.
4391  if (n_dependent_nodes - n_mortar_vertices > 0)
4392  {
4393  // Assemble mass matrix for mortar
4394  //--------------------------------
4395  Vector<double> psi(n_dependent_nodes - n_mortar_vertices);
4396  Vector<double> diag_M(n_dependent_nodes - n_mortar_vertices);
4397  Vector<Vector<double>> shared_node_M(n_mortar_vertices);
4398  for (unsigned i = 0; i < shared_node_M.size(); i++)
4399  {
4400  shared_node_M[i].resize(n_dependent_nodes - n_mortar_vertices);
4401  }
4402 
4403  // Fill in part corresponding to dependent nodal positions (unknown)
4404  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4405  {
4406  // Use L'Hosptal's rule:
4407  psi[i] =
4408  pow(-1.0, int((dependent_element_nnode_1d - 1) - i - 1)) *
4409  -Orthpoly::ddlegendre(dependent_element_nnode_1d - 1,
4410  dependent_nodal_position[non_vertex_pos[i]]);
4411  // Put in contribution
4412  diag_M[i] = psi[i] * dependent_weight[non_vertex_pos[i]];
4413  }
4414 
4415  // Fill in part corresponding to dependent element vertices (known)
4416  for (unsigned v = 0; v < shared_node_M.size(); v++)
4417  {
4418  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4419  {
4420  // Check if denominator is zero
4421  if (std::fabs(dependent_nodal_position[non_vertex_pos[i]] -
4422  dependent_nodal_position[vertex_pos[v]]) >= 1.0e-8)
4423  {
4424  // We're ok
4425  psi[i] =
4426  pow(-1.0, int((dependent_element_nnode_1d - 1) - i - 1)) *
4427  Orthpoly::dlegendre(dependent_element_nnode_1d - 1,
4428  dependent_nodal_position[vertex_pos[v]]) /
4429  (dependent_nodal_position[non_vertex_pos[i]] -
4430  dependent_nodal_position[vertex_pos[v]]);
4431  }
4432  // Check if numerator is zero
4433  else if (std::fabs(Orthpoly::dlegendre(
4434  dependent_element_nnode_1d - 1,
4435  dependent_nodal_position[vertex_pos[v]])) < 1.0e-8)
4436  {
4437  // We can use l'hopital's rule
4438  psi[i] =
4439  pow(-1.0, int((dependent_element_nnode_1d - 1) - i - 1)) *
4440  -Orthpoly::ddlegendre(
4441  dependent_element_nnode_1d - 1,
4442  dependent_nodal_position[non_vertex_pos[i]]);
4443  }
4444  else
4445  {
4446  // We can't use l'hopital's rule
4447  throw OomphLibError(
4448  "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4449  "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4450  OOMPH_EXCEPTION_LOCATION);
4451  }
4452  // Put in contribution
4453  shared_node_M[v][i] = psi[i] * dependent_weight[vertex_pos[v]];
4454  }
4455  }
4456 
4457  // Assemble local projection matrix for mortar
4458  //--------------------------------------------
4459 
4460  // Have only one local projection matrix because there is just one
4461  // master
4462  Vector<Vector<double>> P(n_dependent_nodes - n_mortar_vertices);
4463  for (unsigned i = 0; i < P.size(); i++)
4464  {
4465  P[i].resize(n_master_nodes, 0.0);
4466  }
4467 
4468  // Storage for local coordinate
4469  Vector<double> s(2);
4470 
4471  // Sum contributions from master element shapes (quadrature).
4472  // The order of quadrature must be high enough to evaluate a polynomial
4473  // of order N_s + N_m - 3 exactly, where N_s = n_dependent_nodes, N_m =
4474  // n_master_nodes.
4475  // (Use pointers for the quadrature knots and weights so that
4476  // data is not unnecessarily copied)
4477  // unsigned quadrature_order =
4478  // std::max(dependent_element_nnode_1d,master_element_nnode_1d);
4479  Vector<double>*quadrature_knot, *quadrature_weight;
4480  if (dependent_element_nnode_1d >= master_element_nnode_1d)
4481  {
4482  // Use the same quadrature order as the dependent element (me)
4483  quadrature_knot = &dependent_nodal_position;
4484  quadrature_weight = &dependent_weight;
4485  }
4486  else
4487  {
4488  // Use the same quadrature order as the master element (neighbour)
4489  quadrature_knot = &master_nodal_position;
4490  quadrature_weight = &master_weight;
4491  }
4492 
4493  // Quadrature loop
4494  for (unsigned q = 0; q < (*quadrature_weight).size(); q++)
4495  {
4496  // Evaluate mortar test functions at the quadrature knots in the
4497  // dependent
4498  // s[active_coord_index] = (*quadrature_knot)[q];
4499  Vector<double> s_on_mortar(1);
4500  s_on_mortar[0] = (*quadrature_knot)[q];
4501 
4502  // Get psi
4503  for (unsigned k = 0; k < n_dependent_nodes - n_mortar_vertices; k++)
4504  {
4505  // Check if denominator is zero
4506  if (std::fabs(dependent_nodal_position[non_vertex_pos[k]] -
4507  s_on_mortar[0]) >= 1.0e-08)
4508  {
4509  // We're ok
4510  psi[k] =
4511  pow(-1.0, int((dependent_element_nnode_1d - 1) - k - 1)) *
4512  Orthpoly::dlegendre(dependent_element_nnode_1d - 1,
4513  s_on_mortar[0]) /
4514  (dependent_nodal_position[non_vertex_pos[k]] - s_on_mortar[0]);
4515  }
4516  // Check if numerator is zero
4517  else if (std::fabs(Orthpoly::dlegendre(
4518  dependent_element_nnode_1d - 1, s_on_mortar[0])) <
4519  1.0e-8)
4520  {
4521  // We can use l'Hopital's rule
4522  psi[k] =
4523  pow(-1.0, int((dependent_element_nnode_1d - 1) - k - 1)) *
4524  -Orthpoly::ddlegendre(dependent_element_nnode_1d - 1,
4525  s_on_mortar[0]);
4526  }
4527  else
4528  {
4529  // We can't use l'hopital's rule
4530  throw OomphLibError(
4531  "Cannot use l'Hopital's rule. Dividing by zero is not allowed!",
4532  "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4533  OOMPH_EXCEPTION_LOCATION);
4534  }
4535  }
4536 
4537  // Convert coordinate on mortar to local fraction in dependent element
4538  Vector<double> s_fraction(2);
4539  for (unsigned i = 0; i < 2; i++)
4540  {
4541  s_fraction[i] = (i == active_coord_index) ?
4542  0.5 * (s_on_mortar[0] + 1.0) :
4543  dependent_node_s_fraction[vertex_pos[0]][i];
4544  }
4545 
4546  // Project active coordinate into master element
4547  Vector<double> s_in_neigh(2);
4548  for (unsigned i = 0; i < 2; i++)
4549  {
4550  s_in_neigh[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
4551  (s_hi_neigh[i] - s_lo_neigh[i]);
4552  }
4553 
4554  // Evaluate master shapes at projections of local quadrature knots
4555  neigh_obj_pt->interpolating_basis(
4556  s_in_neigh, master_shapes, value_id);
4557 
4558  // Populate local projection matrix
4559  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4560  {
4561  for (unsigned j = 0; j < n_master_nodes; j++)
4562  {
4563  P[i][j] += master_shapes[master_node_number[j]] * psi[i] *
4564  (*quadrature_weight)[q];
4565  }
4566  }
4567  }
4568 
4569  // Assemble global projection matrices for mortar
4570  //-----------------------------------------------
4571  // Need to subtract contributions from the "known unknowns"
4572  // corresponding to the nodes at the vertices of the mortar
4573 
4574  // Assemble contributions from mortar vertex nodes
4575  for (unsigned v = 0; v < n_mortar_vertices; v++)
4576  {
4577  // Convert coordinate on mortar to local fraction in dependent element
4578  Vector<double> s_fraction(2);
4579  for (unsigned i = 0; i < 2; i++)
4580  {
4581  s_fraction[i] =
4582  (i == active_coord_index) ?
4583  0.5 * (dependent_nodal_position[vertex_pos[v]] + 1.0) :
4584  dependent_node_s_fraction[vertex_pos[0]][i];
4585  }
4586 
4587  // Project active coordinate into master element
4588  Vector<double> s_in_neigh(2);
4589  for (unsigned i = 0; i < 2; i++)
4590  {
4591  s_in_neigh[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
4592  (s_hi_neigh[i] - s_lo_neigh[i]);
4593  }
4594 
4595  // Get master shapes at dependent nodal positions
4596  neigh_obj_pt->interpolating_basis(
4597  s_in_neigh, master_shapes, value_id);
4598 
4599  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4600  {
4601  for (unsigned k = 0; k < n_master_nodes; k++)
4602  {
4603  P[i][k] -=
4604  master_shapes[master_node_number[k]] * shared_node_M[v][i];
4605  }
4606  }
4607  }
4608 
4609  // Solve mortar system
4610  //--------------------
4611  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4612  {
4613  for (unsigned j = 0; j < n_master_nodes; j++)
4614  {
4615  P[i][j] /= diag_M[i];
4616  }
4617  }
4618 
4619  // Create structures to hold the hanging info
4620  //-------------------------------------------
4621  Vector<HangInfo*> hang_info_pt(n_dependent_nodes);
4622  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4623  {
4624  hang_info_pt[i] = new HangInfo(n_master_nodes);
4625  }
4626 
4627  // Copy information to hanging nodes
4628  //----------------------------------
4629  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4630  {
4631  for (unsigned j = 0; j < n_master_nodes; j++)
4632  {
4633  hang_info_pt[i]->set_master_node_pt(j, master_node_pt[j], P[i][j]);
4634  }
4635  }
4636 
4637  // Set pointers to hanging info
4638  //-----------------------------
4639  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
4640  {
4641  // Check that the node shoule actually hang.
4642  // That is, if the polynomial orders of the elements at a p-type
4643  // non-conormity are both odd then the middle node on the edge is
4644  // shared but a hanging scheme has been constructed for it.
4645  bool node_is_really_shared = false;
4646  for (unsigned m = 0; m < hang_info_pt[i]->nmaster(); m++)
4647  {
4648  // We can simply check if the hanging scheme lists itself as a
4649  // master
4650  if (hang_info_pt[i]->master_node_pt(m) ==
4651  dependent_node_pt[non_vertex_pos[i]])
4652  {
4653  node_is_really_shared = true;
4654 
4655 #ifdef PARANOID
4656  // Also check the corresponding weight: it should be 1
4657  if (std::fabs(hang_info_pt[i]->master_weight(m) - 1.0) > 1.0e-06)
4658  {
4659  throw OomphLibError(
4660  "Something fishy here -- with shared node at a mortar vertex",
4661  "PRefineableQElemen<2,INITIAL_NNODE_1D>t::quad_hang_helper()",
4662  OOMPH_EXCEPTION_LOCATION);
4663  }
4664 #endif
4665  }
4666  }
4667 
4668  // Now we can make the node hang, if it isn't a shared node
4669  if (!node_is_really_shared)
4670  {
4671  dependent_node_pt[non_vertex_pos[i]]->set_hanging_pt(
4672  hang_info_pt[i], -1);
4673  }
4674  }
4675 
4676  } // End of case where there are still dependent nodes
4677 
4678 #ifdef PARANOID
4679  // Check all dependent nodes, hanging or otherwise
4680  for (unsigned i = 0; i < n_dependent_nodes; i++)
4681  {
4682  // Check that weights sum to 1 for those that hang
4683  if (dependent_node_pt[i]->is_hanging())
4684  {
4685  // Check that weights sum to 1
4686  double weight_sum = 0.0;
4687  for (unsigned m = 0;
4688  m < dependent_node_pt[i]->hanging_pt()->nmaster();
4689  m++)
4690  {
4691  weight_sum += dependent_node_pt[i]->hanging_pt()->master_weight(m);
4692  }
4693 
4694  // Warn if not
4695  if (std::fabs(weight_sum - 1.0) > 1.0e-08)
4696  {
4697  oomph_info << "Sum of master weights: " << weight_sum << std::endl;
4698  OomphLibWarning(
4699  "Weights in hanging scheme do not sum to 1",
4700  "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4701  OOMPH_EXCEPTION_LOCATION);
4702  }
4703  }
4704  else
4705  {
4706  // Check that this node is shared with the master element if it
4707  // isn't hanging
4708  bool is_master = false;
4709  for (unsigned n = 0; n < n_master_nodes; n++)
4710  {
4711  if (dependent_node_pt[i] == master_node_pt[n])
4712  {
4713  // Node is a master
4714  is_master = true;
4715  break;
4716  }
4717  }
4718 
4719  if (!is_master)
4720  {
4721  // Throw error
4722  std::ostringstream error_string;
4723  error_string << "This node in the dependent element is neither"
4724  << std::endl
4725  << "hanging or shared with a master element."
4726  << std::endl;
4727 
4728  throw OomphLibError(
4729  error_string.str(),
4730  "PRefineableQElement<2,INITIAL_NNODE_1D>::quad_hang_helper()",
4731  OOMPH_EXCEPTION_LOCATION);
4732  }
4733  }
4734  }
4735 #endif
4736 
4737  // Finally, Loop over all dependent nodes and fine-tune their positions
4738  //-----------------------------------------------------------------
4739  // Here we simply set the node's positions to be consistent
4740  // with the hanging scheme. This is not strictly necessary
4741  // because it is done in the mesh adaptation before the node
4742  // becomes non-hanging later on. We make no attempt to ensure
4743  // (strong) continuity in the position across the mortar.
4744  for (unsigned i = 0; i < n_dependent_nodes; i++)
4745  {
4746  // Only fine-tune hanging nodes
4747  if (dependent_node_pt[i]->is_hanging())
4748  {
4749  // If we are doing the position, then
4750  if (value_id == -1)
4751  {
4752  // Get the local coordinate of this dependent node
4753  Vector<double> s_local(2);
4754  this->local_coordinate_of_node(dependent_node_number[i], s_local);
4755 
4756  // Get the position from interpolation in this element via
4757  // the hanging scheme
4758  Vector<double> x_in_neighb(2);
4759  this->interpolated_x(s_local, x_in_neighb);
4760 
4761  // Fine adjust the coordinates (macro map will pick up boundary
4762  // accurately but will lead to different element edges)
4763  dependent_node_pt[i]->x(0) = x_in_neighb[0];
4764  dependent_node_pt[i]->x(1) = x_in_neighb[1];
4765  }
4766  }
4767  }
4768  } // End of case where this interface is to be mortared
4769  }
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
oomph::RefineableElement::ninterpolating_node_1d
virtual unsigned ninterpolating_node_1d(const int &value_id)
Definition: refineable_elements.h:474
oomph::RefineableElement::interpolating_node_pt
virtual Node * interpolating_node_pt(const unsigned &n, const int &value_id)
Definition: refineable_elements.h:437
oomph::RefineableElement::local_one_d_fraction_of_interpolating_node
virtual double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &value_id)
Definition: refineable_elements.h:447
m
int * m
Definition: level2_cplx_impl.h:294
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019
Global_Physical_Variables::P
double P
Uniform pressure.
Definition: TwenteMeshGluing.cpp:77
oomph::Orthpoly::dlegendre
double dlegendre(const unsigned &p, const double &x)
Definition: orthpoly.h:121
oomph::Orthpoly::ddlegendre
double ddlegendre(const unsigned &p, const double &x)
Definition: orthpoly.h:144

References oomph::Orthpoly::ddlegendre(), oomph::Orthpoly::dlegendre(), Global_Physical_Variables::E, boost::multiprecision::fabs(), oomph::Orthpoly::gll_nodes(), oomph::QuadTree::gteq_edge_neighbour(), i, oomph::RefineableElement::interpolating_basis(), oomph::RefineableElement::interpolating_node_pt(), j, k, m, n, N, oomph::RefineableElement::ninterpolating_node(), oomph::RefineableElement::ninterpolating_node_1d(), oomph::Tree::nsons(), oomph::Tree::object_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Global_Physical_Variables::P, Eigen::bfloat16_impl::pow(), Eigen::numext::q, s, oomph::QuadTreeNames::S, oomph::HangInfo::set_master_node_pt(), v, and oomph::QuadTreeNames::W.

◆ rebuild_from_sons()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::rebuild_from_sons ( Mesh *&  mesh_pt)
virtual

Rebuild the element. This needs to find any nodes in the sons which are still required.

Rebuild the element from nodes found in its sons Adjusts its p-order to be the maximum of its sons' p-orders

Implements oomph::RefineableElement.

3209  {
3210  using namespace QuadTreeNames;
3211 
3212  // Get p-orders of sons
3213  unsigned n_sons = this->tree_pt()->nsons();
3214  Vector<unsigned> son_p_order(n_sons);
3215  unsigned max_son_p_order = 0;
3216  for (unsigned ison = 0; ison < n_sons; ison++)
3217  {
3218  PRefineableElement* el_pt = dynamic_cast<PRefineableElement*>(
3219  this->tree_pt()->son_pt(ison)->object_pt());
3220  son_p_order[ison] = el_pt->p_order();
3221  if (son_p_order[ison] > max_son_p_order)
3222  max_son_p_order = son_p_order[ison];
3223  }
3224 
3225  unsigned old_Nnode = this->nnode();
3226  unsigned old_P_order = this->p_order();
3227  // Set p-order of the element
3228  this->p_order() = max_son_p_order;
3229 
3230  // Change integration scheme
3231  delete this->integral_pt();
3232  switch (this->p_order())
3233  {
3234  case 2:
3235  this->set_integration_scheme(new GaussLobattoLegendre<2, 2>);
3236  break;
3237  case 3:
3238  this->set_integration_scheme(new GaussLobattoLegendre<2, 3>);
3239  break;
3240  case 4:
3241  this->set_integration_scheme(new GaussLobattoLegendre<2, 4>);
3242  break;
3243  case 5:
3244  this->set_integration_scheme(new GaussLobattoLegendre<2, 5>);
3245  break;
3246  case 6:
3247  this->set_integration_scheme(new GaussLobattoLegendre<2, 6>);
3248  break;
3249  case 7:
3250  this->set_integration_scheme(new GaussLobattoLegendre<2, 7>);
3251  break;
3252  default:
3253  std::ostringstream error_message;
3254  error_message << "\nERROR: Exceeded maximum polynomial order for";
3255  error_message << "\n integration scheme.\n";
3256  throw OomphLibError(error_message.str(),
3257  OOMPH_CURRENT_FUNCTION,
3258  OOMPH_EXCEPTION_LOCATION);
3259  }
3260 
3261  // Back up pointers to old nodes before they are lost
3262  Vector<Node*> old_node_pt(old_Nnode);
3263  for (unsigned n = 0; n < old_Nnode; n++)
3264  {
3265  old_node_pt[n] = this->node_pt(n);
3266  }
3267 
3268  // Allocate new space for Nodes (at the element level)
3269  this->set_n_node(this->p_order() * this->p_order());
3270 
3271  // Copy vertex nodes which were populated in the pre-build
3272  this->node_pt(0) = old_node_pt[0];
3273  this->node_pt(this->p_order() - 1) = old_node_pt[old_P_order - 1];
3274  this->node_pt(this->p_order() * (this->p_order() - 1)) =
3275  old_node_pt[(old_P_order) * (old_P_order - 1)];
3276  this->node_pt(this->p_order() * this->p_order() - 1) =
3277  old_node_pt[(old_P_order) * (old_P_order)-1];
3278 
3279  // Copy midpoint nodes from sons if new p-order is odd
3280  if (this->p_order() % 2 == 1)
3281  {
3282  // Work out which is midpoint node
3283  unsigned midpoint = (this->p_order() - 1) / 2;
3284 
3285  // Bottom edge
3286  this->node_pt(midpoint) = dynamic_cast<RefineableQElement<2>*>(
3287  quadtree_pt()->son_pt(SW)->object_pt())
3288  ->vertex_node_pt(1);
3289  // Left edge
3290  this->node_pt(midpoint * this->p_order()) =
3291  dynamic_cast<RefineableQElement<2>*>(
3292  quadtree_pt()->son_pt(SW)->object_pt())
3293  ->vertex_node_pt(2);
3294  // Top edge
3295  this->node_pt((this->p_order() - 1) * this->p_order() + midpoint) =
3296  dynamic_cast<RefineableQElement<2>*>(
3297  quadtree_pt()->son_pt(NE)->object_pt())
3298  ->vertex_node_pt(2);
3299  // Right edge
3300  this->node_pt((midpoint + 1) * this->p_order() - 1) =
3301  dynamic_cast<RefineableQElement<2>*>(
3302  quadtree_pt()->son_pt(NE)->object_pt())
3303  ->vertex_node_pt(1);
3304  }
3305 
3306 
3307  // The timestepper should be the same for all nodes and node 0 should
3308  // never be deleted.
3309  if (this->node_pt(0) == 0)
3310  {
3311  throw OomphLibError("The Corner node (0) does not exist",
3312  OOMPH_CURRENT_FUNCTION,
3313  OOMPH_EXCEPTION_LOCATION);
3314  }
3315 
3316  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
3317  unsigned ntstorage = time_stepper_pt->ntstorage();
3318 
3319  unsigned jnod = 0;
3320  Vector<double> s_fraction(2), s(2);
3321  // Loop over the nodes in the element
3322  unsigned n_p = this->nnode_1d();
3323  for (unsigned i0 = 0; i0 < n_p; i0++)
3324  {
3325  // Get the fractional position of the node
3326  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3327  // Local coordinate
3328  s[0] = -1.0 + 2.0 * s_fraction[0];
3329 
3330  for (unsigned i1 = 0; i1 < n_p; i1++)
3331  {
3332  // Get the fractional position of the node in the direction of s[1]
3333  s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
3334  // Local coordinate in father element
3335  s[1] = -1.0 + 2.0 * s_fraction[1];
3336 
3337  // Set the local node number
3338  jnod = i0 + n_p * i1;
3339 
3340  // Initialise flag: So far, this node hasn't been built
3341  // or copied yet
3342  bool node_done = false;
3343 
3344  // Get the pointer to the node in this element, returns NULL
3345  // if there is not node
3346  Node* created_node_pt = this->get_node_at_local_coordinate(s);
3347 
3348  // Does this node already exist in this element?
3349  //----------------------------------------------
3350  if (created_node_pt != 0)
3351  {
3352  // Copy node across
3353  this->node_pt(jnod) = created_node_pt;
3354 
3355  // Node has been created by copy
3356  node_done = true;
3357  }
3358  // Node does not exist in this element but might already
3359  //------------------------------------------------------
3360  // have been created by neighbouring elements
3361  //-------------------------------------------
3362  else
3363  {
3364  // Was the node created by one of its neighbours
3365  // Whether or not the node lies on an edge can be calculated
3366  // by from the fractional position
3367  bool is_periodic = false;
3368  created_node_pt = node_created_by_neighbour(s_fraction, is_periodic);
3369 
3370  // If the node was so created, assign the pointers
3371  if (created_node_pt != 0)
3372  {
3373  // If the node is periodic
3374  if (is_periodic)
3375  {
3376  throw OomphLibError("Cannot handle periodic nodes yet",
3377  OOMPH_CURRENT_FUNCTION,
3378  OOMPH_EXCEPTION_LOCATION);
3379  }
3380  // Non-periodic case, just set the pointer
3381  else
3382  {
3383  this->node_pt(jnod) = created_node_pt;
3384  }
3385  // Node has been created
3386  node_done = true;
3387  }
3388  } // Node does not exist in this element
3389 
3390  // Node has not been built anywhere ---> build it here
3391  if (!node_done)
3392  {
3393  // First, find the son element in which the node should live
3394 
3395  // Find coordinates in the sons
3396  Vector<double> s_in_son(2);
3397  using namespace QuadTreeNames;
3398  int son = -10;
3399  // If negative on the west side
3400  if (s_fraction[0] < 0.5)
3401  {
3402  // On the south side
3403  if (s_fraction[1] < 0.5)
3404  {
3405  // It's the southwest son
3406  son = SW;
3407  s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
3408  s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
3409  }
3410  // On the north side
3411  else
3412  {
3413  // It's the northwest son
3414  son = NW;
3415  s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
3416  s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
3417  }
3418  }
3419  else
3420  {
3421  // On the south side
3422  if (s_fraction[1] < 0.5)
3423  {
3424  // It's the southeast son
3425  son = SE;
3426  s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
3427  s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
3428  }
3429  // On the north side
3430  else
3431  {
3432  // It's the northeast son
3433  son = NE;
3434  s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
3435  s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
3436  }
3437  }
3438 
3439  // Get the pointer to the son element in which the new node
3440  // would live
3441  PRefineableQElement<2, INITIAL_NNODE_1D>* son_el_pt =
3442  dynamic_cast<PRefineableQElement<2, INITIAL_NNODE_1D>*>(
3443  this->tree_pt()->son_pt(son)->object_pt());
3444 
3445  // If we are rebuilding, then worry about the boundary conditions
3446  // Find the boundary of the node
3447  // Initially none
3448  int boundary = Tree::OMEGA;
3449  // If we are on the western boundary
3450  if (i0 == 0)
3451  {
3452  boundary = W;
3453  }
3454  // If we are on the eastern boundary
3455  else if (i0 == n_p - 1)
3456  {
3457  boundary = E;
3458  }
3459 
3460  // If we are on the southern boundary
3461  if (i1 == 0)
3462  {
3463  // If we already have already set the boundary, we're on a corner
3464  switch (boundary)
3465  {
3466  case W:
3467  boundary = SW;
3468  break;
3469  case E:
3470  boundary = SE;
3471  break;
3472  // Boundary not set
3473  default:
3474  boundary = S;
3475  break;
3476  }
3477  }
3478  // If we are the northern bounadry
3479  else if (i1 == n_p - 1)
3480  {
3481  // If we already have a boundary
3482  switch (boundary)
3483  {
3484  case W:
3485  boundary = NW;
3486  break;
3487  case E:
3488  boundary = NE;
3489  break;
3490  default:
3491  boundary = N;
3492  break;
3493  }
3494  }
3495 
3496  // set of boundaries that this edge in the son lives on
3497  std::set<unsigned> boundaries;
3498 
3499  // Now get the boundary conditions from the son
3500  // The boundaries will be common to the son because there can be
3501  // no rotations here
3502  if (boundary != Tree::OMEGA)
3503  {
3504  son_el_pt->get_boundaries(boundary, boundaries);
3505  }
3506 
3507  // If the node lives on a boundary:
3508  // Construct a boundary node,
3509  // Get boundary conditions and
3510  // update all lookup schemes
3511  if (boundaries.size() > 0)
3512  {
3513  // Construct the new node
3514  created_node_pt =
3515  this->construct_boundary_node(jnod, time_stepper_pt);
3516 
3517  // Get the boundary conditions from the son
3518  Vector<int> bound_cons(this->ncont_interpolated_values());
3519  son_el_pt->get_bcs(boundary, bound_cons);
3520 
3521  // Loop over the values and pin if necessary
3522  unsigned nval = created_node_pt->nvalue();
3523  for (unsigned k = 0; k < nval; k++)
3524  {
3525  if (bound_cons[k])
3526  {
3527  created_node_pt->pin(k);
3528  }
3529  }
3530 
3531  // Solid node? If so, deal with the positional boundary
3532  // conditions:
3533  SolidNode* solid_node_pt =
3534  dynamic_cast<SolidNode*>(created_node_pt);
3535  if (solid_node_pt != 0)
3536  {
3537  // Get the positional boundary conditions from the father:
3538  unsigned n_dim = created_node_pt->ndim();
3539  Vector<int> solid_bound_cons(n_dim);
3540  RefineableSolidQElement<2>* son_solid_el_pt =
3541  dynamic_cast<RefineableSolidQElement<2>*>(son_el_pt);
3542 #ifdef PARANOID
3543  if (son_solid_el_pt == 0)
3544  {
3545  std::string error_message =
3546  "We have a SolidNode outside a refineable SolidElement\n";
3547  error_message +=
3548  "during mesh refinement -- this doesn't make sense\n";
3549 
3550  throw OomphLibError(error_message,
3551  OOMPH_CURRENT_FUNCTION,
3552  OOMPH_EXCEPTION_LOCATION);
3553  }
3554 #endif
3555  son_solid_el_pt->get_solid_bcs(boundary, solid_bound_cons);
3556 
3557  // Loop over the positions and pin, if necessary
3558  for (unsigned k = 0; k < n_dim; k++)
3559  {
3560  if (solid_bound_cons[k])
3561  {
3562  solid_node_pt->pin_position(k);
3563  }
3564  }
3565  } // End of if solid_node_pt
3566 
3567 
3568  // Next we update the boundary look-up schemes
3569  // Loop over the boundaries stored in the set
3570  for (std::set<unsigned>::iterator it = boundaries.begin();
3571  it != boundaries.end();
3572  ++it)
3573  {
3574  // Add the node to the boundary
3575  mesh_pt->add_boundary_node(*it, created_node_pt);
3576 
3577  // If we have set an intrinsic coordinate on this
3578  // mesh boundary then it must also be interpolated on
3579  // the new node
3580  // Now interpolate the intrinsic boundary coordinate
3581  if (mesh_pt->boundary_coordinate_exists(*it) == true)
3582  {
3583  Vector<double> zeta(1);
3584  son_el_pt->interpolated_zeta_on_edge(
3585  *it, boundary, s_in_son, zeta);
3586 
3587  created_node_pt->set_coordinates_on_boundary(*it, zeta);
3588  }
3589  }
3590  }
3591  // Otherwise the node is not on a Mesh boundary
3592  // and we create a normal "bulk" node
3593  else
3594  {
3595  // Construct the new node
3596  created_node_pt = this->construct_node(jnod, time_stepper_pt);
3597  }
3598 
3599  // Now we set the position and values at the newly created node
3600 
3601  // In the first instance use macro element or FE representation
3602  // to create past and present nodal positions.
3603  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
3604  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
3605  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
3606  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
3607  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
3608  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
3609  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
3610 
3611  // Loop over # of history values
3612  // Loop over # of history values
3613  for (unsigned t = 0; t < ntstorage; t++)
3614  {
3615  using namespace QuadTreeNames;
3616  // Get the position from the son
3617  Vector<double> x_prev(2);
3618 
3619  // Now let's fill in the value
3620  son_el_pt->get_x(t, s_in_son, x_prev);
3621  for (unsigned i = 0; i < 2; i++)
3622  {
3623  created_node_pt->x(t, i) = x_prev[i];
3624  }
3625  }
3626 
3627  // Now set up the values
3628  // Loop over all history values
3629  for (unsigned t = 0; t < ntstorage; t++)
3630  {
3631  // Get values from father element
3632  // Note: get_interpolated_values() sets Vector size itself.
3633  Vector<double> prev_values;
3634  son_el_pt->get_interpolated_values(t, s_in_son, prev_values);
3635 
3636  // Initialise the values at the new node
3637  for (unsigned k = 0; k < created_node_pt->nvalue(); k++)
3638  {
3639  created_node_pt->set_value(t, k, prev_values[k]);
3640  }
3641  }
3642 
3643  // Add the node to the mesh
3644  mesh_pt->add_node_pt(created_node_pt);
3645 
3646  // Check if the element is an algebraic element
3647  AlgebraicElementBase* alg_el_pt =
3648  dynamic_cast<AlgebraicElementBase*>(this);
3649 
3650  // If we do have an algebraic element
3651  if (alg_el_pt != 0)
3652  {
3653  std::string error_message =
3654  "Have not implemented rebuilding from sons for";
3655  error_message += "Algebraic p-refineable elements yet\n";
3656 
3657  throw OomphLibError(
3658  error_message,
3659  "PRefineableQElement<2,INITIAL_NNODE_1D>::rebuild_from_sons()",
3660  OOMPH_EXCEPTION_LOCATION);
3661  }
3662 
3663  } // End of the case when we build the node ourselves
3664  }
3665  } // End of loop over all nodes in element
3666 
3667 
3668  // If the element is a MacroElementNodeUpdateElement, set the update
3669  // parameters for the current element's nodes. These need to be reset
3670  // (as in MacroElementNodeUpdateElement<ELEMENT>::rebuild_from_sons())
3671  // because the nodes in this element have changed
3672  MacroElementNodeUpdateElementBase* m_el_pt =
3673  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
3674  if (m_el_pt != 0)
3675  {
3676  // Loop over the nodes
3677  for (unsigned j = 0; j < this->nnode(); j++)
3678  {
3679  // Get local coordinate in element (Vector sets its own size)
3680  Vector<double> s_in_node_update_element;
3681  this->local_coordinate_of_node(j, s_in_node_update_element);
3682 
3683  // Get vector of geometric objects
3684  Vector<GeomObject*> geom_object_pt(m_el_pt->geom_object_pt());
3685 
3686  // Pass the lot to the node
3687  static_cast<MacroElementNodeUpdateNode*>(this->node_pt(j))
3688  ->set_node_update_info(
3689  this, s_in_node_update_element, geom_object_pt);
3690  }
3691  }
3692 
3693  // MacroElementNodeUpdateElementBase* m_el_pt=dynamic_cast<
3694  // MacroElementNodeUpdateElementBase*>(this);
3695  // if(m_el_pt!=0)
3696  // {
3697  // Loop over all nodes in element again, to re-set the positions
3698  // This must be done using the new element's macro-element
3699  // representation, rather than the old version which may be
3700  // of a different p-order!
3701  for (unsigned i0 = 0; i0 < n_p; i0++)
3702  {
3703  // Get the fractional position of the node
3704  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
3705  // Local coordinate
3706  s[0] = -1.0 + 2.0 * s_fraction[0];
3707 
3708  for (unsigned i1 = 0; i1 < n_p; i1++)
3709  {
3710  // Get the fractional position of the node in the direction of s[1]
3711  s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
3712  // Local coordinate in father element
3713  s[1] = -1.0 + 2.0 * s_fraction[1];
3714 
3715  // Set the local node number
3716  jnod = i0 + n_p * i1;
3717 
3718  // Loop over # of history values
3719  for (unsigned t = 0; t < ntstorage; t++)
3720  {
3721  // Get position from father element -- this uses the macro
3722  // element representation if appropriate. If the node
3723  // turns out to be a hanging node later on, then
3724  // its position gets adjusted in line with its
3725  // hanging node interpolation.
3726  Vector<double> x_prev(2);
3727  this->get_x(t, s, x_prev);
3728 
3729  // Set previous positions of the new node
3730  for (unsigned i = 0; i < 2; i++)
3731  {
3732  this->node_pt(jnod)->x(t, i) = x_prev[i];
3733  }
3734  }
3735  }
3736  }
3737  // }
3738  }
oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::get_node_at_local_coordinate
Node * get_node_at_local_coordinate(const Vector< double > &s) const
Return the node at the specified local coordinate.
Definition: hp_refineable_elements.cc:1317
oomph::QuadElementBase::vertex_node_pt
virtual Node * vertex_node_pt(const unsigned &j) const =0
Pointer to the j-th vertex node in the element.
oomph::Tree::object_pt
RefineableElement * object_pt() const
Definition: tree.h:88
oomph::Tree::nsons
unsigned nsons() const
Return number of sons (zero if it's a leaf node)
Definition: tree.h:129
oomph::Tree::son_pt
Tree * son_pt(const int &son_index) const
Definition: tree.h:103

References oomph::Mesh::add_boundary_node(), oomph::Mesh::add_node_pt(), oomph::Mesh::boundary_coordinate_exists(), Global_Physical_Variables::E, oomph::MacroElementNodeUpdateElementBase::geom_object_pt(), oomph::RefineableQElement< 2 >::get_bcs(), oomph::RefineableQElement< 2 >::get_boundaries(), oomph::RefineableElement::get_interpolated_values(), oomph::RefineableSolidQElement< 2 >::get_solid_bcs(), oomph::FiniteElement::get_x(), i, oomph::RefineableQElement< 2 >::interpolated_zeta_on_edge(), j, k, n, N, oomph::Node::ndim(), oomph::QuadTreeNames::NE, oomph::TimeStepper::ntstorage(), oomph::Data::nvalue(), oomph::QuadTreeNames::NW, oomph::Tree::OMEGA, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::PRefineableElement::p_order(), oomph::Data::pin(), oomph::SolidNode::pin_position(), s, oomph::QuadTreeNames::S, oomph::QuadTreeNames::SE, oomph::Node::set_coordinates_on_boundary(), oomph::Data::set_value(), oomph::Global_string_for_annotation::string(), oomph::QuadTreeNames::SW, plotPSD::t, oomph::QuadTreeNames::W, oomph::Node::x(), and Eigen::zeta().

◆ shape()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::shape ( const Vector< double > &  s,
Shape &  psi 
) const
virtual

Overload the shape functions.

Shape functions for PRefineableQElement<DIM>

Implements oomph::FiniteElement.

2909  {
2910  switch (p_order())
2911  {
2912  case 2:
2913  {
2914  // Call the OneDimensional Shape functions
2915  OneDimensionalLegendreShape<2>::calculate_nodal_positions();
2916  OneDimensionalLegendreShape<2> psi1(s[0]), psi2(s[1]);
2917 
2918  // Now let's loop over the nodal points in the element
2919  // and copy the values back in
2920  for (unsigned i = 0; i < 2; i++)
2921  {
2922  for (unsigned j = 0; j < 2; j++)
2923  {
2924  psi(2 * i + j) = psi2[i] * psi1[j];
2925  }
2926  }
2927  break;
2928  }
2929  case 3:
2930  {
2931  // Call the OneDimensional Shape functions
2932  OneDimensionalLegendreShape<3>::calculate_nodal_positions();
2933  OneDimensionalLegendreShape<3> psi1(s[0]), psi2(s[1]);
2934 
2935  // Now let's loop over the nodal points in the element
2936  // and copy the values back in
2937  for (unsigned i = 0; i < 3; i++)
2938  {
2939  for (unsigned j = 0; j < 3; j++)
2940  {
2941  psi(3 * i + j) = psi2[i] * psi1[j];
2942  }
2943  }
2944  break;
2945  }
2946  case 4:
2947  {
2948  // Call the OneDimensional Shape functions
2949  OneDimensionalLegendreShape<4>::calculate_nodal_positions();
2950  OneDimensionalLegendreShape<4> psi1(s[0]), psi2(s[1]);
2951 
2952  // Now let's loop over the nodal points in the element
2953  // and copy the values back in
2954  for (unsigned i = 0; i < 4; i++)
2955  {
2956  for (unsigned j = 0; j < 4; j++)
2957  {
2958  psi(4 * i + j) = psi2[i] * psi1[j];
2959  }
2960  }
2961  break;
2962  }
2963  case 5:
2964  {
2965  // Call the OneDimensional Shape functions
2966  OneDimensionalLegendreShape<5>::calculate_nodal_positions();
2967  OneDimensionalLegendreShape<5> psi1(s[0]), psi2(s[1]);
2968 
2969  // Now let's loop over the nodal points in the element
2970  // and copy the values back in
2971  for (unsigned i = 0; i < 5; i++)
2972  {
2973  for (unsigned j = 0; j < 5; j++)
2974  {
2975  psi(5 * i + j) = psi2[i] * psi1[j];
2976  }
2977  }
2978  break;
2979  }
2980  case 6:
2981  {
2982  // Call the OneDimensional Shape functions
2983  OneDimensionalLegendreShape<6>::calculate_nodal_positions();
2984  OneDimensionalLegendreShape<6> psi1(s[0]), psi2(s[1]);
2985 
2986  // Now let's loop over the nodal points in the element
2987  // and copy the values back in
2988  for (unsigned i = 0; i < 6; i++)
2989  {
2990  for (unsigned j = 0; j < 6; j++)
2991  {
2992  psi(6 * i + j) = psi2[i] * psi1[j];
2993  }
2994  }
2995  break;
2996  }
2997  case 7:
2998  {
2999  // Call the OneDimensional Shape functions
3000  OneDimensionalLegendreShape<7>::calculate_nodal_positions();
3001  OneDimensionalLegendreShape<7> psi1(s[0]), psi2(s[1]);
3002 
3003  // Now let's loop over the nodal points in the element
3004  // and copy the values back in
3005  for (unsigned i = 0; i < 7; i++)
3006  {
3007  for (unsigned j = 0; j < 7; j++)
3008  {
3009  psi(7 * i + j) = psi2[i] * psi1[j];
3010  }
3011  }
3012  break;
3013  }
3014  default:
3015  std::ostringstream error_message;
3016  error_message << "\nERROR: Exceeded maximum polynomial order for";
3017  error_message << "\n polynomial order for shape functions.\n";
3018  throw OomphLibError(error_message.str(),
3019  OOMPH_CURRENT_FUNCTION,
3020  OOMPH_EXCEPTION_LOCATION);
3021  }
3022  }

References oomph::OneDimensionalLegendreShape< NNODE_1D >::calculate_nodal_positions(), i, j, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

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