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

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

#include <hp_refineable_elements.h>

+ Inheritance diagram for oomph::PRefineableQElement< 3, 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< 3 >
 RefineableQElement ()
 Constructor: Pass refinement level (default 0 = root) More...
 
 RefineableQElement (const RefineableQElement< 3 > &dummy)=delete
 Broken copy constructor. More...
 
virtual ~RefineableQElement ()
 Broken assignment operator. More...
 
unsigned required_nsons () const
 A refineable brick element has eight 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...
 
OcTree * octree_pt ()
 Pointer to octree representation of this element. More...
 
OcTree * octree_pt () const
 Pointer to octree representation of this element. More...
 
void setup_hanging_nodes (Vector< std::ofstream * > &output_stream)
 
- 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::BrickElementBase
 BrickElementBase ()
 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 oc_hang_helper (const int &value_id, const int &my_face, std::ofstream &output_hangfile)
 
- Protected Member Functions inherited from oomph::RefineableQElement< 3 >
void setup_father_bounds ()
 
void get_face_bcs (const int &edge, Vector< int > &bound_cons) const
 
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_face (const unsigned &boundary, const int &face, const Vector< double > &s, Vector< double > &zeta)
 
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< 3 >
typedef void(RefineableQElement< 3 >::* VoidMemFctPt) ()
 
- 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< 3 >
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< 3, INITIAL_NNODE_1D >

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

Constructor & Destructor Documentation

◆ PRefineableQElement()

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

Constructor.

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

◆ ~PRefineableQElement()

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

Destructor.

253 {}

Member Function Documentation

◆ check_integrity()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 3, 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.

8338  {
8339  // Not yet implemented -- doing nothing for now.
8340 
8341  // Dummy set max_error to 0
8342  max_error = 0.0;
8343  return;
8344  }
MeshRefinement::max_error
double max_error
Definition: MortaringCantileverCompareToNonMortaring.cpp:188

References MeshRefinement::max_error.

◆ d2shape_local()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 3, 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.

7620  {
7621  std::ostringstream error_message;
7622  error_message
7623  << "\nd2shape_local currently not implemented for this element\n";
7624  throw OomphLibError(
7625  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
7626  }
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< 3, 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.

7426  {
7427  switch (P_order)
7428  {
7429  case 2:
7430  {
7431  // Call the shape functions and derivatives
7432  OneDimensionalLegendreShape<2>::calculate_nodal_positions();
7433  OneDimensionalLegendreShape<2> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7434  OneDimensionalLegendreDShape<2> dpsi1ds(s[0]), dpsi2ds(s[1]),
7435  dpsi3ds(s[2]);
7436 
7437  // Index for the shape functions
7438  unsigned index = 0;
7439  // Loop over shape functions in element
7440  for (unsigned i = 0; i < P_order; i++)
7441  {
7442  for (unsigned j = 0; j < P_order; j++)
7443  {
7444  for (unsigned k = 0; k < P_order; k++)
7445  {
7446  // Assign the values
7447  dpsids(index, 0) = psi3[i] * psi2[j] * dpsi1ds[k];
7448  dpsids(index, 1) = psi3[i] * dpsi2ds[j] * psi1[k];
7449  dpsids(index, 2) = dpsi3ds[i] * psi2[j] * psi1[k];
7450  psi[index] = psi3[i] * psi2[j] * psi1[k];
7451  // Increment the index
7452  ++index;
7453  }
7454  }
7455  }
7456  break;
7457  }
7458  case 3:
7459  {
7460  // Call the shape functions and derivatives
7461  OneDimensionalLegendreShape<3>::calculate_nodal_positions();
7462  OneDimensionalLegendreShape<3> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7463  OneDimensionalLegendreDShape<3> dpsi1ds(s[0]), dpsi2ds(s[1]),
7464  dpsi3ds(s[2]);
7465 
7466  // Index for the shape functions
7467  unsigned index = 0;
7468  // Loop over shape functions in element
7469  for (unsigned i = 0; i < P_order; i++)
7470  {
7471  for (unsigned j = 0; j < P_order; j++)
7472  {
7473  for (unsigned k = 0; k < P_order; k++)
7474  {
7475  // Assign the values
7476  dpsids(index, 0) = psi3[i] * psi2[j] * dpsi1ds[k];
7477  dpsids(index, 1) = psi3[i] * dpsi2ds[j] * psi1[k];
7478  dpsids(index, 2) = dpsi3ds[i] * psi2[j] * psi1[k];
7479  psi[index] = psi3[i] * psi2[j] * psi1[k];
7480  // Increment the index
7481  ++index;
7482  }
7483  }
7484  }
7485  break;
7486  }
7487  case 4:
7488  {
7489  // Call the shape functions and derivatives
7490  OneDimensionalLegendreShape<4>::calculate_nodal_positions();
7491  OneDimensionalLegendreShape<4> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7492  OneDimensionalLegendreDShape<4> dpsi1ds(s[0]), dpsi2ds(s[1]),
7493  dpsi3ds(s[2]);
7494 
7495  // Index for the shape functions
7496  unsigned index = 0;
7497  // Loop over shape functions in element
7498  for (unsigned i = 0; i < P_order; i++)
7499  {
7500  for (unsigned j = 0; j < P_order; j++)
7501  {
7502  for (unsigned k = 0; k < P_order; k++)
7503  {
7504  // Assign the values
7505  dpsids(index, 0) = psi3[i] * psi2[j] * dpsi1ds[k];
7506  dpsids(index, 1) = psi3[i] * dpsi2ds[j] * psi1[k];
7507  dpsids(index, 2) = dpsi3ds[i] * psi2[j] * psi1[k];
7508  psi[index] = psi3[i] * psi2[j] * psi1[k];
7509  // Increment the index
7510  ++index;
7511  }
7512  }
7513  }
7514  break;
7515  }
7516  case 5:
7517  {
7518  // Call the shape functions and derivatives
7519  OneDimensionalLegendreShape<5>::calculate_nodal_positions();
7520  OneDimensionalLegendreShape<5> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7521  OneDimensionalLegendreDShape<5> dpsi1ds(s[0]), dpsi2ds(s[1]),
7522  dpsi3ds(s[2]);
7523 
7524  // Index for the shape functions
7525  unsigned index = 0;
7526  // Loop over shape functions in element
7527  for (unsigned i = 0; i < P_order; i++)
7528  {
7529  for (unsigned j = 0; j < P_order; j++)
7530  {
7531  for (unsigned k = 0; k < P_order; k++)
7532  {
7533  // Assign the values
7534  dpsids(index, 0) = psi3[i] * psi2[j] * dpsi1ds[k];
7535  dpsids(index, 1) = psi3[i] * dpsi2ds[j] * psi1[k];
7536  dpsids(index, 2) = dpsi3ds[i] * psi2[j] * psi1[k];
7537  psi[index] = psi3[i] * psi2[j] * psi1[k];
7538  // Increment the index
7539  ++index;
7540  }
7541  }
7542  }
7543  break;
7544  }
7545  case 6:
7546  {
7547  // Call the shape functions and derivatives
7548  OneDimensionalLegendreShape<6>::calculate_nodal_positions();
7549  OneDimensionalLegendreShape<6> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7550  OneDimensionalLegendreDShape<6> dpsi1ds(s[0]), dpsi2ds(s[1]),
7551  dpsi3ds(s[2]);
7552 
7553  // Index for the shape functions
7554  unsigned index = 0;
7555  // Loop over shape functions in element
7556  for (unsigned i = 0; i < P_order; i++)
7557  {
7558  for (unsigned j = 0; j < P_order; j++)
7559  {
7560  for (unsigned k = 0; k < P_order; k++)
7561  {
7562  // Assign the values
7563  dpsids(index, 0) = psi3[i] * psi2[j] * dpsi1ds[k];
7564  dpsids(index, 1) = psi3[i] * dpsi2ds[j] * psi1[k];
7565  dpsids(index, 2) = dpsi3ds[i] * psi2[j] * psi1[k];
7566  psi[index] = psi3[i] * psi2[j] * psi1[k];
7567  // Increment the index
7568  ++index;
7569  }
7570  }
7571  }
7572  break;
7573  }
7574  case 7:
7575  {
7576  // Call the shape functions and derivatives
7577  OneDimensionalLegendreShape<7>::calculate_nodal_positions();
7578  OneDimensionalLegendreShape<7> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7579  OneDimensionalLegendreDShape<7> dpsi1ds(s[0]), dpsi2ds(s[1]),
7580  dpsi3ds(s[2]);
7581 
7582  // Index for the shape functions
7583  unsigned index = 0;
7584  // Loop over shape functions in element
7585  for (unsigned i = 0; i < P_order; i++)
7586  {
7587  for (unsigned j = 0; j < P_order; j++)
7588  {
7589  for (unsigned k = 0; k < P_order; k++)
7590  {
7591  // Assign the values
7592  dpsids(index, 0) = psi3[i] * psi2[j] * dpsi1ds[k];
7593  dpsids(index, 1) = psi3[i] * dpsi2ds[j] * psi1[k];
7594  dpsids(index, 2) = dpsi3ds[i] * psi2[j] * psi1[k];
7595  psi[index] = psi3[i] * psi2[j] * psi1[k];
7596  // Increment the index
7597  ++index;
7598  }
7599  }
7600  }
7601  break;
7602  }
7603  default:
7604  std::ostringstream error_message;
7605  error_message << "\nERROR: Exceeded maximum polynomial order for";
7606  error_message << "\n polynomial order for shape functions.\n";
7607  throw OomphLibError(error_message.str(),
7608  OOMPH_CURRENT_FUNCTION,
7609  OOMPH_EXCEPTION_LOCATION);
7610  }
7611  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
oomph::OneDimensionalLegendreShape::calculate_nodal_positions
static void calculate_nodal_positions()
Static function used to populate the stored positions.
Definition: shape.h:1241
oomph::PRefineableElement::P_order
unsigned P_order
The polynomial expansion order of the elemental basis functions.
Definition: refineable_elements.h:655
s
RealScalar s
Definition: level1_cplx_impl.h:130
k
char char char int int * k
Definition: level2_impl.h:374
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

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

◆ further_setup_hanging_nodes()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 3, 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< 3 >.

283 {}

◆ get_node_at_local_coordinate()

template<unsigned INITIAL_NNODE_1D>
Node * oomph::PRefineableQElement< 3, 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.

4893  {
4894  unsigned Nnode_1d = this->nnode_1d();
4895  // Load the tolerance into a local variable
4896  double tol = FiniteElement::Node_location_tolerance;
4897  // There are two possible indices.
4898  Vector<int> index(3);
4899 
4900  // Loop over indices
4901  for (unsigned i = 0; i < 3; i++)
4902  {
4903  // Determine the index
4904  // -------------------
4905 
4906  bool is_found = false;
4907 
4908  // If we are at the lower limit, the index is zero
4909  if (std::fabs(s[i] + 1.0) < tol)
4910  {
4911  index[i] = 0;
4912  is_found = true;
4913  }
4914  // If we are at the upper limit, the index is the number of nodes minus 1
4915  else if (std::fabs(s[i] - 1.0) < tol)
4916  {
4917  index[i] = Nnode_1d - 1;
4918  is_found = true;
4919  }
4920  // Otherwise, we have to calculate the index in general
4921  else
4922  {
4923  // Compute Gauss-Lobatto-Legendre node positions
4924  Vector<double> z;
4925  Orthpoly::gll_nodes(Nnode_1d, z);
4926  // Loop over possible internal nodal positions
4927  for (unsigned n = 1; n < Nnode_1d - 1; n++)
4928  {
4929  if (std::fabs(z[n] - s[i]) < tol)
4930  {
4931  index[i] = n;
4932  is_found = true;
4933  break;
4934  }
4935  }
4936  }
4937 
4938  if (!is_found)
4939  {
4940  // No matching nodes
4941  return 0;
4942  }
4943  }
4944  // If we've got here we have a node, so let's return a pointer to it
4945  return this->node_pt(index[0] + Nnode_1d * index[1] +
4946  Nnode_1d * Nnode_1d * index[2]);
4947  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::Node_location_tolerance
static const double Node_location_tolerance
Definition: elements.h:1374
oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::nnode_1d
unsigned nnode_1d() const
Definition: hp_refineable_elements.h:288
GlobalParameters::Nnode_1d
unsigned Nnode_1d
The number of nodes in one direction (default=2)
Definition: structured_cubic_point_source.cc:56
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
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< 3, INITIAL_NNODE_1D >::initial_p_order ( ) const
inlinevirtual

Get the initial P_order.

Implements oomph::PRefineableElement.

295  {
296  return INITIAL_NNODE_1D;
297  }

◆ initial_setup()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 3, 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.

5770  {
5771  // Storage for pointer to my father (in octree impersonation)
5772  OcTree* father_pt = 0;
5773 
5774  // Check if an adopted father has been specified
5775  if (adopted_father_pt != 0)
5776  {
5777  // Get pointer to my father (in octree impersonation)
5778  father_pt = dynamic_cast<OcTree*>(adopted_father_pt);
5779  }
5780  // Check if element is in a tree
5781  else if (Tree_pt != 0)
5782  {
5783  // Get pointer to my father (in octree impersonation)
5784  father_pt = dynamic_cast<OcTree*>(octree_pt()->father_pt());
5785  }
5786  else
5787  {
5788  // throw OomphLibError(
5789  // "Element not in a tree, and no adopted father has been
5790  // specified!",
5791  // "PRefineableQElement<2,INITIAL_NNODE_1D>::initial_setup()",
5792  // OOMPH_EXCEPTION_LOCATION);
5793  }
5794 
5795  // Check if element has father
5796  if (father_pt != 0 || initial_p_order != 0)
5797  {
5798  if (father_pt != 0)
5799  {
5800  PRefineableQElement<3, INITIAL_NNODE_1D>* father_el_pt =
5801  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(
5802  father_pt->object_pt());
5803  if (father_el_pt != 0)
5804  {
5805  unsigned father_p_order = father_el_pt->p_order();
5806  // Set p-order to that of father
5807  P_order = father_p_order;
5808  }
5809  }
5810  else
5811  {
5812  P_order = initial_p_order;
5813  }
5814 
5815  // Now sort out the element...
5816  // (has p^3 nodes)
5817  unsigned new_n_node = P_order * P_order * P_order;
5818 
5819  // Allocate new space for Nodes (at the element level)
5820  this->set_n_node(new_n_node);
5821 
5822  // Set integration scheme
5823  delete this->integral_pt();
5824  switch (P_order)
5825  {
5826  case 2:
5827  this->set_integration_scheme(new GaussLobattoLegendre<3, 2>);
5828  break;
5829  case 3:
5830  this->set_integration_scheme(new GaussLobattoLegendre<3, 3>);
5831  break;
5832  case 4:
5833  this->set_integration_scheme(new GaussLobattoLegendre<3, 4>);
5834  break;
5835  case 5:
5836  this->set_integration_scheme(new GaussLobattoLegendre<3, 5>);
5837  break;
5838  case 6:
5839  this->set_integration_scheme(new GaussLobattoLegendre<3, 6>);
5840  break;
5841  case 7:
5842  this->set_integration_scheme(new GaussLobattoLegendre<3, 7>);
5843  break;
5844  default:
5845  std::ostringstream error_message;
5846  error_message << "\nERROR: Exceeded maximum polynomial order for";
5847  error_message << "\n integration scheme.\n";
5848  throw OomphLibError(error_message.str(),
5849  OOMPH_CURRENT_FUNCTION,
5850  OOMPH_EXCEPTION_LOCATION);
5851  }
5852  }
5853  }
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::PRefineableQElement< 3, INITIAL_NNODE_1D >::initial_p_order
unsigned initial_p_order() const
Get the initial P_order.
Definition: hp_refineable_elements.h:294
oomph::RefineableElement::Tree_pt
Tree * Tree_pt
A pointer to a general tree object.
Definition: refineable_elements.h:100
oomph::RefineableQElement< 3 >::octree_pt
OcTree * octree_pt()
Pointer to octree representation of this element.
Definition: refineable_brick_element.h:144
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< 3, 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.

4779  {
4780  s.resize(3);
4781  unsigned Nnode_1d = this->nnode_1d();
4782  unsigned n0 = n % Nnode_1d;
4783  unsigned n1 = unsigned(double(n) / double(Nnode_1d)) % Nnode_1d;
4784  unsigned n2 = unsigned(double(n) / double(Nnode_1d * Nnode_1d));
4785 
4786  switch (Nnode_1d)
4787  {
4788  case 2:
4789  OneDimensionalLegendreShape<2>::calculate_nodal_positions();
4790  s[0] = OneDimensionalLegendreShape<2>::nodal_position(n0);
4791  s[1] = OneDimensionalLegendreShape<2>::nodal_position(n1);
4792  s[2] = OneDimensionalLegendreShape<2>::nodal_position(n2);
4793  break;
4794  case 3:
4795  OneDimensionalLegendreShape<3>::calculate_nodal_positions();
4796  s[0] = OneDimensionalLegendreShape<3>::nodal_position(n0);
4797  s[1] = OneDimensionalLegendreShape<3>::nodal_position(n1);
4798  s[2] = OneDimensionalLegendreShape<3>::nodal_position(n2);
4799  break;
4800  case 4:
4801  OneDimensionalLegendreShape<4>::calculate_nodal_positions();
4802  s[0] = OneDimensionalLegendreShape<4>::nodal_position(n0);
4803  s[1] = OneDimensionalLegendreShape<4>::nodal_position(n1);
4804  s[2] = OneDimensionalLegendreShape<4>::nodal_position(n2);
4805  break;
4806  case 5:
4807  OneDimensionalLegendreShape<5>::calculate_nodal_positions();
4808  s[0] = OneDimensionalLegendreShape<5>::nodal_position(n0);
4809  s[1] = OneDimensionalLegendreShape<5>::nodal_position(n1);
4810  s[2] = OneDimensionalLegendreShape<5>::nodal_position(n2);
4811  break;
4812  case 6:
4813  OneDimensionalLegendreShape<6>::calculate_nodal_positions();
4814  s[0] = OneDimensionalLegendreShape<6>::nodal_position(n0);
4815  s[1] = OneDimensionalLegendreShape<6>::nodal_position(n1);
4816  s[2] = OneDimensionalLegendreShape<6>::nodal_position(n2);
4817  break;
4818  case 7:
4819  OneDimensionalLegendreShape<7>::calculate_nodal_positions();
4820  s[0] = OneDimensionalLegendreShape<7>::nodal_position(n0);
4821  s[1] = OneDimensionalLegendreShape<7>::nodal_position(n1);
4822  s[2] = OneDimensionalLegendreShape<7>::nodal_position(n2);
4823  break;
4824  default:
4825  std::ostringstream error_message;
4826  error_message << "\nERROR: Exceeded maximum polynomial order for";
4827  error_message << "\n shape functions.\n";
4828  throw OomphLibError(error_message.str(),
4829  OOMPH_CURRENT_FUNCTION,
4830  OOMPH_EXCEPTION_LOCATION);
4831  break;
4832  }
4833  }
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< 3, 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.

4839  {
4840  this->local_coordinate_of_node(n, s_fraction);
4841  s_fraction[0] = 0.5 * (s_fraction[0] + 1.0);
4842  s_fraction[1] = 0.5 * (s_fraction[1] + 1.0);
4843  s_fraction[2] = 0.5 * (s_fraction[2] + 1.0);
4844  }
oomph::PRefineableQElement< 3, 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:4777

References n.

◆ local_one_d_fraction_of_node()

template<unsigned INITIAL_NNODE_1D>
double oomph::PRefineableQElement< 3, 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.

4849  {
4850  switch (this->nnode_1d())
4851  {
4852  case 2:
4853  OneDimensionalLegendreShape<2>::calculate_nodal_positions();
4854  return 0.5 *
4855  (OneDimensionalLegendreShape<2>::nodal_position(n1d) + 1.0);
4856  case 3:
4857  OneDimensionalLegendreShape<3>::calculate_nodal_positions();
4858  return 0.5 *
4859  (OneDimensionalLegendreShape<3>::nodal_position(n1d) + 1.0);
4860  case 4:
4861  OneDimensionalLegendreShape<4>::calculate_nodal_positions();
4862  return 0.5 *
4863  (OneDimensionalLegendreShape<4>::nodal_position(n1d) + 1.0);
4864  case 5:
4865  OneDimensionalLegendreShape<5>::calculate_nodal_positions();
4866  return 0.5 *
4867  (OneDimensionalLegendreShape<5>::nodal_position(n1d) + 1.0);
4868  case 6:
4869  OneDimensionalLegendreShape<6>::calculate_nodal_positions();
4870  return 0.5 *
4871  (OneDimensionalLegendreShape<6>::nodal_position(n1d) + 1.0);
4872  case 7:
4873  OneDimensionalLegendreShape<7>::calculate_nodal_positions();
4874  return 0.5 *
4875  (OneDimensionalLegendreShape<7>::nodal_position(n1d) + 1.0);
4876  default:
4877  std::ostringstream error_message;
4878  error_message << "\nERROR: Exceeded maximum polynomial order for";
4879  error_message << "\n shape functions.\n";
4880  throw OomphLibError(error_message.str(),
4881  OOMPH_CURRENT_FUNCTION,
4882  OOMPH_EXCEPTION_LOCATION);
4883  return 0.0;
4884  }
4885  }

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

289  {
290  return this->p_order();
291  }
oomph::PRefineableElement::p_order
unsigned & p_order()
Access function to P_order.
Definition: refineable_elements.h:705

◆ node_created_by_neighbour()

template<unsigned INITIAL_NNODE_1D>
Node * oomph::PRefineableQElement< 3, 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).

Reimplemented from oomph::RefineableQElement< 3 >.

4958  {
4959  using namespace OcTreeNames;
4960 
4961  // Calculate the faces/edges on which the node lies
4962  Vector<int> faces;
4963  Vector<int> edges;
4964 
4965  if (s_fraction[0] == 0.0)
4966  {
4967  faces.push_back(L);
4968  if (s_fraction[1] == 0.0)
4969  {
4970  edges.push_back(LD);
4971  }
4972  if (s_fraction[2] == 0.0)
4973  {
4974  edges.push_back(LB);
4975  }
4976  if (s_fraction[1] == 1.0)
4977  {
4978  edges.push_back(LU);
4979  }
4980  if (s_fraction[2] == 1.0)
4981  {
4982  edges.push_back(LF);
4983  }
4984  }
4985 
4986  if (s_fraction[0] == 1.0)
4987  {
4988  faces.push_back(R);
4989  if (s_fraction[1] == 0.0)
4990  {
4991  edges.push_back(RD);
4992  }
4993  if (s_fraction[2] == 0.0)
4994  {
4995  edges.push_back(RB);
4996  }
4997  if (s_fraction[1] == 1.0)
4998  {
4999  edges.push_back(RU);
5000  }
5001  if (s_fraction[2] == 1.0)
5002  {
5003  edges.push_back(RF);
5004  }
5005  }
5006 
5007  if (s_fraction[1] == 0.0)
5008  {
5009  faces.push_back(D);
5010  if (s_fraction[2] == 0.0)
5011  {
5012  edges.push_back(DB);
5013  }
5014  if (s_fraction[2] == 1.0)
5015  {
5016  edges.push_back(DF);
5017  }
5018  }
5019 
5020  if (s_fraction[1] == 1.0)
5021  {
5022  faces.push_back(U);
5023  if (s_fraction[2] == 0.0)
5024  {
5025  edges.push_back(UB);
5026  }
5027  if (s_fraction[2] == 1.0)
5028  {
5029  edges.push_back(UF);
5030  }
5031  }
5032 
5033  if (s_fraction[2] == 0.0)
5034  {
5035  faces.push_back(B);
5036  }
5037 
5038  if (s_fraction[2] == 1.0)
5039  {
5040  faces.push_back(F);
5041  }
5042 
5043  // Find the number of faces
5044  unsigned n_face = faces.size();
5045 
5046  // Find the number of edges
5047  unsigned n_edge = edges.size();
5048 
5049  Vector<unsigned> translate_s(3);
5050  Vector<double> s_lo_neigh(3);
5051  Vector<double> s_hi_neigh(3);
5052  Vector<double> s(3);
5053 
5054  int neigh_face, diff_level;
5055  bool in_neighbouring_tree;
5056 
5057  // Loop over the faces on which the node lies
5058  //------------------------------------------
5059  for (unsigned j = 0; j < n_face; j++)
5060  {
5061  // Find pointer to neighbouring element along face
5062  OcTree* neigh_pt;
5063  neigh_pt = octree_pt()->gteq_face_neighbour(faces[j],
5064  translate_s,
5065  s_lo_neigh,
5066  s_hi_neigh,
5067  neigh_face,
5068  diff_level,
5069  in_neighbouring_tree);
5070 
5071  // Neighbour exists
5072  if (neigh_pt != 0)
5073  {
5074  // Have any of its vertex nodes been created yet?
5075  // (Must look in incomplete neighbours because after the
5076  // pre-build they may contain pointers to the required nodes. e.g.
5077  // h-refinement of neighbouring linear and quadratic elements)
5078  bool a_vertex_node_is_built = false;
5079  QElement<3, INITIAL_NNODE_1D>* neigh_obj_pt =
5080  dynamic_cast<QElement<3, INITIAL_NNODE_1D>*>(neigh_pt->object_pt());
5081  if (neigh_obj_pt == 0)
5082  {
5083  throw OomphLibError("Not a quad element!",
5084  OOMPH_CURRENT_FUNCTION,
5085  OOMPH_EXCEPTION_LOCATION);
5086  }
5087  for (unsigned vnode = 0; vnode < neigh_obj_pt->nvertex_node(); vnode++)
5088  {
5089  if (neigh_obj_pt->vertex_node_pt(vnode) != 0)
5090  a_vertex_node_is_built = true;
5091  break;
5092  }
5093  if (a_vertex_node_is_built)
5094  {
5095  // We now need to translate the nodal location, defined in terms
5096  // of the fractional coordinates of the present element into
5097  // those of its neighbour. For this we use the information returned
5098  // to use from the octree function.
5099 
5100  // Calculate the local coordinate in the neighbour
5101  // Note that we need to use the translation scheme in case
5102  // the local coordinates are swapped in the neighbour.
5103  for (unsigned i = 0; i < 3; i++)
5104  {
5105  s[i] = s_lo_neigh[i] +
5106  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
5107  }
5108 
5109  // Find the node in the neighbour
5110  Node* neighbour_node_pt =
5111  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
5112 
5113  // If there is a node, return it
5114  if (neighbour_node_pt != 0)
5115  {
5116  // Now work out whether it's a periodic boundary
5117  // only possible if we have moved into a neighbouring tree
5118  if (in_neighbouring_tree)
5119  {
5120  // Return whether the neighbour is periodic
5121  is_periodic =
5122  octree_pt()->root_pt()->is_neighbour_periodic(faces[j]);
5123  }
5124 
5125  return neighbour_node_pt;
5126  }
5127  }
5128  }
5129  } // End of loop over faces
5130 
5131 
5132  // Loop over the edges on which the node lies
5133  //------------------------------------------
5134  for (unsigned j = 0; j < n_edge; j++)
5135  {
5136  // Even if we restrict ourselves to true edge neighbours (i.e.
5137  // elements that are not also face neighbours) there may be multiple
5138  // edge neighbours across the edges between multiple root octrees.
5139  // When making the first call to OcTree::gteq_true_edge_neighbour(...)
5140  // we simply return the first one of these multiple edge neighbours
5141  // (if there are any at all, of course) and also return the total number
5142  // of true edge neighbours. If the node in question already exists
5143  // on the first edge neighbour we're done. If it doesn't it may exist
5144  // on other edge neighbours so we repeat the process over all
5145  // other edge neighbours (bailing out if a node is found, of course).
5146 
5147  // Initially return the zero-th true edge neighbour
5148  unsigned i_root_edge_neighbour = 0;
5149 
5150  // Initialise the total number of true edge neighbours
5151  unsigned nroot_edge_neighbour = 0;
5152 
5153  // Keep searching until we've found the node or until we've checked
5154  // all available edge neighbours
5155  bool keep_searching = true;
5156  while (keep_searching)
5157  {
5158  // Find pointer to neighbouring element along edge
5159  OcTree* neigh_pt;
5160  neigh_pt = octree_pt()->gteq_true_edge_neighbour(edges[j],
5161  i_root_edge_neighbour,
5162  nroot_edge_neighbour,
5163  translate_s,
5164  s_lo_neigh,
5165  s_hi_neigh,
5166  neigh_face,
5167  diff_level);
5168 
5169  // Neighbour exists
5170  if (neigh_pt != 0)
5171  {
5172  // Have any of its vertex nodes been created yet?
5173  // (Must look in incomplete neighbours because after the
5174  // pre-build they may contain pointers to the required nodes. e.g.
5175  // h-refinement of neighbouring linear and quadratic elements)
5176  bool a_vertex_node_is_built = false;
5177  QElement<3, INITIAL_NNODE_1D>* neigh_obj_pt =
5178  dynamic_cast<QElement<3, INITIAL_NNODE_1D>*>(neigh_pt->object_pt());
5179  if (neigh_obj_pt == 0)
5180  {
5181  throw OomphLibError("Not a quad element!",
5182  OOMPH_CURRENT_FUNCTION,
5183  OOMPH_EXCEPTION_LOCATION);
5184  }
5185  for (unsigned vnode = 0; vnode < neigh_obj_pt->nvertex_node();
5186  vnode++)
5187  {
5188  if (neigh_obj_pt->vertex_node_pt(vnode) != 0)
5189  a_vertex_node_is_built = true;
5190  break;
5191  }
5192  if (a_vertex_node_is_built)
5193  {
5194  // We now need to translate the nodal location, defined in terms
5195  // of the fractional coordinates of the present element into
5196  // those of its neighbour. For this we use the information returned
5197  // to use from the octree function.
5198 
5199  // Calculate the local coordinate in the neighbour
5200  // Note that we need to use the translation scheme in case
5201  // the local coordinates are swapped in the neighbour.
5202  for (unsigned i = 0; i < 3; i++)
5203  {
5204  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
5205  (s_hi_neigh[i] - s_lo_neigh[i]);
5206  }
5207 
5208  // Find the node in the neighbour
5209  Node* neighbour_node_pt =
5210  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
5211 
5212  // If there is a node, return it
5213  if (neighbour_node_pt != 0)
5214  {
5215  // Get the faces on which the edge lies
5216  Vector<int> faces_attached_to_edge =
5217  OcTree::faces_of_common_edge(edges[j]);
5218 
5219  // Get the number of entries in the vector
5220  unsigned n_faces_attached_to_edge = faces_attached_to_edge.size();
5221 
5222  // Loop over the faces
5223  for (unsigned i_face = 0; i_face < n_faces_attached_to_edge;
5224  i_face++)
5225  {
5226  // Is the node periodic in the face direction?
5227  is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
5228  faces_attached_to_edge[i_face]);
5229 
5230  // Check if the edge is periodic in the i_face-th face direction
5231  if (is_periodic)
5232  {
5233  // We're done!
5234  break;
5235  }
5236  } // for (unsigned
5237  // i_face=0;i_face<n_faces_attached_to_edge;i_face++)
5238 
5239  return neighbour_node_pt;
5240  }
5241  }
5242  }
5243 
5244  // Keep searching, but only if there are further edge neighbours
5245  // Try next root edge neighbour
5246  i_root_edge_neighbour++;
5247 
5248  // Have we exhausted the search
5249  if (i_root_edge_neighbour >= nroot_edge_neighbour)
5250  {
5251  keep_searching = false;
5252  }
5253 
5254  } // End of while keep searching over all true edge neighbours
5255 
5256  } // End of loop over edges
5257 
5258  // Node not found, return null
5259  return 0;
5260  }
D
dominoes D
Definition: Domino.cpp:55
L
MatrixXd L
Definition: LLT_example.cpp:6
R
@ R
Definition: StatisticsVector.h:21
oomph::Matrix
Definition: matrices.h:74
oomph::OcTree::faces_of_common_edge
static Vector< int > faces_of_common_edge(const int &edge)
Function that, given an edge, returns the two faces on which it.
Definition: octree.cc:268
oomph::OcTree::gteq_face_neighbour
OcTree * gteq_face_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_sw, Vector< double > &s_ne, int &face, int &diff_level, bool &in_neighbouring_tree) const
Definition: octree.cc:3373
oomph::OcTree::gteq_true_edge_neighbour
OcTree * gteq_true_edge_neighbour(const int &direction, const unsigned &i_root_edge_neighbour, unsigned &nroot_edge_neighbour, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, int &edge, int &diff_level) const
Definition: octree.cc:3618
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
RachelsAdvectionDiffusion::U
double U
Swimming speed.
Definition: two_d_variable_diff_adapt.cc:53
oomph::OcTreeNames::DF
@ DF
Definition: octree.h:67
oomph::OcTreeNames::RD
@ RD
Definition: octree.h:62
oomph::OcTreeNames::RB
@ RB
Definition: octree.h:58
oomph::OcTreeNames::UF
@ UF
Definition: octree.h:68
oomph::OcTreeNames::RU
@ RU
Definition: octree.h:64
oomph::OcTreeNames::F
@ F
Definition: octree.h:74
oomph::OcTreeNames::LF
@ LF
Definition: octree.h:65
oomph::OcTreeNames::LD
@ LD
Definition: octree.h:61
oomph::OcTreeNames::RF
@ RF
Definition: octree.h:66
oomph::OcTreeNames::LU
@ LU
Definition: octree.h:63
oomph::OcTreeNames::DB
@ DB
Definition: octree.h:59
oomph::OcTreeNames::UB
@ UB
Definition: octree.h:60
oomph::OcTreeNames::LB
@ LB
Definition: octree.h:57

References D, oomph::OcTreeNames::DB, oomph::OcTreeNames::DF, oomph::OcTreeNames::F, oomph::OcTree::faces_of_common_edge(), oomph::FiniteElement::get_node_at_local_coordinate(), oomph::OcTree::gteq_face_neighbour(), oomph::OcTree::gteq_true_edge_neighbour(), i, j, L, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, oomph::Tree::object_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, oomph::OcTreeNames::RF, oomph::OcTreeNames::RU, s, RachelsAdvectionDiffusion::U, oomph::OcTreeNames::UB, and oomph::OcTreeNames::UF.

◆ node_created_by_son_of_neighbour()

template<unsigned INITIAL_NNODE_1D>
Node * oomph::PRefineableQElement< 3, 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< 3 >.

5273  {
5274  using namespace OcTreeNames;
5275 
5276  // Calculate the faces/edges on which the node lies
5277  Vector<int> faces;
5278  Vector<int> edges;
5279 
5280  if (s_fraction[0] == 0.0)
5281  {
5282  faces.push_back(L);
5283  if (s_fraction[1] == 0.0)
5284  {
5285  edges.push_back(LD);
5286  }
5287  if (s_fraction[2] == 0.0)
5288  {
5289  edges.push_back(LB);
5290  }
5291  if (s_fraction[1] == 1.0)
5292  {
5293  edges.push_back(LU);
5294  }
5295  if (s_fraction[2] == 1.0)
5296  {
5297  edges.push_back(LF);
5298  }
5299  }
5300 
5301  if (s_fraction[0] == 1.0)
5302  {
5303  faces.push_back(R);
5304  if (s_fraction[1] == 0.0)
5305  {
5306  edges.push_back(RD);
5307  }
5308  if (s_fraction[2] == 0.0)
5309  {
5310  edges.push_back(RB);
5311  }
5312  if (s_fraction[1] == 1.0)
5313  {
5314  edges.push_back(RU);
5315  }
5316  if (s_fraction[2] == 1.0)
5317  {
5318  edges.push_back(RF);
5319  }
5320  }
5321 
5322  if (s_fraction[1] == 0.0)
5323  {
5324  faces.push_back(D);
5325  if (s_fraction[2] == 0.0)
5326  {
5327  edges.push_back(DB);
5328  }
5329  if (s_fraction[2] == 1.0)
5330  {
5331  edges.push_back(DF);
5332  }
5333  }
5334 
5335  if (s_fraction[1] == 1.0)
5336  {
5337  faces.push_back(U);
5338  if (s_fraction[2] == 0.0)
5339  {
5340  edges.push_back(UB);
5341  }
5342  if (s_fraction[2] == 1.0)
5343  {
5344  edges.push_back(UF);
5345  }
5346  }
5347 
5348  if (s_fraction[2] == 0.0)
5349  {
5350  faces.push_back(B);
5351  }
5352 
5353  if (s_fraction[2] == 1.0)
5354  {
5355  faces.push_back(F);
5356  }
5357 
5358  // Find the number of faces
5359  unsigned n_face = faces.size();
5360 
5361  // Find the number of edges
5362  unsigned n_edge = edges.size();
5363 
5364  Vector<unsigned> translate_s(3);
5365  Vector<double> s_lo_neigh(3);
5366  Vector<double> s_hi_neigh(3);
5367  Vector<double> s(3);
5368 
5369  int neigh_face, diff_level;
5370  bool in_neighbouring_tree;
5371 
5372  // Loop over the faces on which the node lies
5373  //------------------------------------------
5374  for (unsigned j = 0; j < n_face; j++)
5375  {
5376  // Find pointer to neighbouring element along face
5377  OcTree* neigh_pt;
5378  neigh_pt = octree_pt()->gteq_face_neighbour(faces[j],
5379  translate_s,
5380  s_lo_neigh,
5381  s_hi_neigh,
5382  neigh_face,
5383  diff_level,
5384  in_neighbouring_tree);
5385 
5386  // Neighbour exists
5387  if (neigh_pt != 0)
5388  {
5389  // Have its nodes been created yet?
5390  // (Must look in sons of unfinished neighbours too!!!)
5391  if (true)
5392  {
5393  // We now need to translate the nodal location, defined in terms
5394  // of the fractional coordinates of the present element into
5395  // those of its neighbour. For this we use the information returned
5396  // to use from the octree function.
5397 
5398  // Calculate the local coordinate in the neighbour
5399  // Note that we need to use the translation scheme in case
5400  // the local coordinates are swapped in the neighbour.
5401  for (unsigned i = 0; i < 3; i++)
5402  {
5403  s[i] = s_lo_neigh[i] +
5404  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
5405  }
5406 
5407  // Check if the element has sons
5408  if (neigh_pt->nsons() != 0)
5409  {
5410  // First, find the son element in which the node should live
5411 
5412  // Find coordinates in the sons
5413  Vector<double> s_in_son(3);
5414  int son = -10;
5415  // On the left
5416  if (s[0] < 0.0)
5417  {
5418  // On the bottom
5419  if (s[1] < 0.0)
5420  {
5421  // On the back
5422  if (s[2] < 0.0)
5423  {
5424  // It's the LDB son
5425  son = OcTreeNames::LDB;
5426  s_in_son[0] = 1.0 + 2.0 * s[0];
5427  s_in_son[1] = 1.0 + 2.0 * s[1];
5428  s_in_son[2] = 1.0 + 2.0 * s[2];
5429  }
5430  // On the front
5431  else
5432  {
5433  // It's the LDF son
5434  son = OcTreeNames::LDF;
5435  s_in_son[0] = 1.0 + 2.0 * s[0];
5436  s_in_son[1] = 1.0 + 2.0 * s[1];
5437  s_in_son[2] = -1.0 + 2.0 * s[2];
5438  }
5439  }
5440  // On the top
5441  else
5442  {
5443  // On the back
5444  if (s[2] < 0.0)
5445  {
5446  // It's the LUB son
5447  son = OcTreeNames::LUB;
5448  s_in_son[0] = 1.0 + 2.0 * s[0];
5449  s_in_son[1] = -1.0 + 2.0 * s[1];
5450  s_in_son[2] = 1.0 + 2.0 * s[2];
5451  }
5452  // On the front
5453  else
5454  {
5455  // It's the LUF son
5456  son = OcTreeNames::LUF;
5457  s_in_son[0] = 1.0 + 2.0 * s[0];
5458  s_in_son[1] = -1.0 + 2.0 * s[1];
5459  s_in_son[2] = -1.0 + 2.0 * s[2];
5460  }
5461  }
5462  }
5463  // On the right
5464  else
5465  {
5466  // On the bottom
5467  if (s[1] < 0.0)
5468  {
5469  // On the back
5470  if (s[2] < 0.0)
5471  {
5472  // It's the RDB son
5473  son = OcTreeNames::RDB;
5474  s_in_son[0] = -1.0 + 2.0 * s[0];
5475  s_in_son[1] = 1.0 + 2.0 * s[1];
5476  s_in_son[2] = 1.0 + 2.0 * s[2];
5477  }
5478  // On the front
5479  else
5480  {
5481  // It's the RDF son
5482  son = OcTreeNames::RDF;
5483  s_in_son[0] = -1.0 + 2.0 * s[0];
5484  s_in_son[1] = 1.0 + 2.0 * s[1];
5485  s_in_son[2] = -1.0 + 2.0 * s[2];
5486  }
5487  }
5488  // On the top
5489  else
5490  {
5491  // On the back
5492  if (s[2] < 0.0)
5493  {
5494  // It's the RUB son
5495  son = OcTreeNames::RUB;
5496  s_in_son[0] = -1.0 + 2.0 * s[0];
5497  s_in_son[1] = -1.0 + 2.0 * s[1];
5498  s_in_son[2] = 1.0 + 2.0 * s[2];
5499  }
5500  // On the front
5501  else
5502  {
5503  // It's the RUF son
5504  son = OcTreeNames::RUF;
5505  s_in_son[0] = -1.0 + 2.0 * s[0];
5506  s_in_son[1] = -1.0 + 2.0 * s[1];
5507  s_in_son[2] = -1.0 + 2.0 * s[2];
5508  }
5509  }
5510  }
5511 
5512  // Find the node in the neighbour's son
5513  Node* neighbour_son_node_pt =
5514  neigh_pt->son_pt(son)->object_pt()->get_node_at_local_coordinate(
5515  s_in_son);
5516 
5517  // If there is a node, return it
5518  if (neighbour_son_node_pt != 0)
5519  {
5520  // Now work out whether it's a periodic boundary
5521  // only possible if we have moved into a neighbouring tree
5522  if (in_neighbouring_tree)
5523  {
5524  // Return whether the neighbour is periodic
5525  is_periodic =
5526  octree_pt()->root_pt()->is_neighbour_periodic(faces[j]);
5527  }
5528 
5529  // Return the pointer to the neighbouring node
5530  return neighbour_son_node_pt;
5531  }
5532  }
5533  }
5534  }
5535  } // End of loop over faces
5536 
5537 
5538  // Loop over the edges on which the node lies
5539  //------------------------------------------
5540  for (unsigned j = 0; j < n_edge; j++)
5541  {
5542  // Even if we restrict ourselves to true edge neighbours (i.e.
5543  // elements that are not also face neighbours) there may be multiple
5544  // edge neighbours across the edges between multiple root octrees.
5545  // When making the first call to OcTree::gteq_true_edge_neighbour(...)
5546  // we simply return the first one of these multiple edge neighbours
5547  // (if there are any at all, of course) and also return the total number
5548  // of true edge neighbours. If the node in question already exists
5549  // on the first edge neighbour we're done. If it doesn't it may exist
5550  // on other edge neighbours so we repeat the process over all
5551  // other edge neighbours (bailing out if a node is found, of course).
5552 
5553  // Initially return the zero-th true edge neighbour
5554  unsigned i_root_edge_neighbour = 0;
5555 
5556  // Initialise the total number of true edge neighbours
5557  unsigned nroot_edge_neighbour = 0;
5558 
5559  // Keep searching until we've found the node or until we've checked
5560  // all available edge neighbours
5561  bool keep_searching = true;
5562  while (keep_searching)
5563  {
5564  // Find pointer to neighbouring element along edge
5565  OcTree* neigh_pt;
5566  neigh_pt = octree_pt()->gteq_true_edge_neighbour(edges[j],
5567  i_root_edge_neighbour,
5568  nroot_edge_neighbour,
5569  translate_s,
5570  s_lo_neigh,
5571  s_hi_neigh,
5572  neigh_face,
5573  diff_level);
5574 
5575  // Neighbour exists
5576  if (neigh_pt != 0)
5577  {
5578  // Have its nodes been created yet?
5579  // (Must look in sons of unfinished neighbours too!!!)
5580  if (true)
5581  {
5582  // We now need to translate the nodal location, defined in terms
5583  // of the fractional coordinates of the present element into
5584  // those of its neighbour. For this we use the information returned
5585  // to use from the octree function.
5586 
5587  // Calculate the local coordinate in the neighbour
5588  // Note that we need to use the translation scheme in case
5589  // the local coordinates are swapped in the neighbour.
5590  for (unsigned i = 0; i < 3; i++)
5591  {
5592  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
5593  (s_hi_neigh[i] - s_lo_neigh[i]);
5594  }
5595 
5596  // Check if the element has sons
5597  if (neigh_pt->nsons() != 0)
5598  {
5599  // First, find the son element in which the node should live
5600 
5601  // Find coordinates in the sons
5602  Vector<double> s_in_son(3);
5603  int son = -10;
5604  // On the left
5605  if (s[0] < 0.0)
5606  {
5607  // On the bottom
5608  if (s[1] < 0.0)
5609  {
5610  // On the back
5611  if (s[2] < 0.0)
5612  {
5613  // It's the LDB son
5614  son = OcTreeNames::LDB;
5615  s_in_son[0] = 1.0 + 2.0 * s[0];
5616  s_in_son[1] = 1.0 + 2.0 * s[1];
5617  s_in_son[2] = 1.0 + 2.0 * s[2];
5618  }
5619  // On the front
5620  else
5621  {
5622  // It's the LDF son
5623  son = OcTreeNames::LDF;
5624  s_in_son[0] = 1.0 + 2.0 * s[0];
5625  s_in_son[1] = 1.0 + 2.0 * s[1];
5626  s_in_son[2] = -1.0 + 2.0 * s[2];
5627  }
5628  }
5629  // On the top
5630  else
5631  {
5632  // On the back
5633  if (s[2] < 0.0)
5634  {
5635  // It's the LUB son
5636  son = OcTreeNames::LUB;
5637  s_in_son[0] = 1.0 + 2.0 * s[0];
5638  s_in_son[1] = -1.0 + 2.0 * s[1];
5639  s_in_son[2] = 1.0 + 2.0 * s[2];
5640  }
5641  // On the front
5642  else
5643  {
5644  // It's the LUF son
5645  son = OcTreeNames::LUF;
5646  s_in_son[0] = 1.0 + 2.0 * s[0];
5647  s_in_son[1] = -1.0 + 2.0 * s[1];
5648  s_in_son[2] = -1.0 + 2.0 * s[2];
5649  }
5650  }
5651  }
5652  // On the right
5653  else
5654  {
5655  // On the bottom
5656  if (s[1] < 0.0)
5657  {
5658  // On the back
5659  if (s[2] < 0.0)
5660  {
5661  // It's the RDB son
5662  son = OcTreeNames::RDB;
5663  s_in_son[0] = -1.0 + 2.0 * s[0];
5664  s_in_son[1] = 1.0 + 2.0 * s[1];
5665  s_in_son[2] = 1.0 + 2.0 * s[2];
5666  }
5667  // On the front
5668  else
5669  {
5670  // It's the RDF son
5671  son = OcTreeNames::RDF;
5672  s_in_son[0] = -1.0 + 2.0 * s[0];
5673  s_in_son[1] = 1.0 + 2.0 * s[1];
5674  s_in_son[2] = -1.0 + 2.0 * s[2];
5675  }
5676  }
5677  // On the top
5678  else
5679  {
5680  // On the back
5681  if (s[2] < 0.0)
5682  {
5683  // It's the RUB son
5684  son = OcTreeNames::RUB;
5685  s_in_son[0] = -1.0 + 2.0 * s[0];
5686  s_in_son[1] = -1.0 + 2.0 * s[1];
5687  s_in_son[2] = 1.0 + 2.0 * s[2];
5688  }
5689  // On the front
5690  else
5691  {
5692  // It's the RUF son
5693  son = OcTreeNames::RUF;
5694  s_in_son[0] = -1.0 + 2.0 * s[0];
5695  s_in_son[1] = -1.0 + 2.0 * s[1];
5696  s_in_son[2] = -1.0 + 2.0 * s[2];
5697  }
5698  }
5699  }
5700 
5701  // Find the node in the neighbour's son
5702  Node* neighbour_son_node_pt =
5703  neigh_pt->son_pt(son)
5704  ->object_pt()
5705  ->get_node_at_local_coordinate(s_in_son);
5706 
5707  // If there is a node, return it
5708  if (neighbour_son_node_pt != 0)
5709  {
5710  // Get the faces on which the edge lies
5711  Vector<int> faces_attached_to_edge =
5712  OcTree::faces_of_common_edge(edges[j]);
5713 
5714  // Get the number of entries in the vector
5715  unsigned n_faces_attached_to_edge =
5716  faces_attached_to_edge.size();
5717 
5718  // Loop over the faces
5719  for (unsigned i_face = 0; i_face < n_faces_attached_to_edge;
5720  i_face++)
5721  {
5722  // Is the node periodic in the face direction?
5723  is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
5724  faces_attached_to_edge[i_face]);
5725 
5726  // Check if the edge is periodic in the i_face-th face
5727  // direction
5728  if (is_periodic)
5729  {
5730  // We're done!
5731  break;
5732  }
5733  } // for (unsigned
5734  // i_face=0;i_face<n_faces_attached_to_edge;i_face++)
5735 
5736  // Return the pointer to the neighbouring node
5737  return neighbour_son_node_pt;
5738  }
5739  }
5740  }
5741  }
5742 
5743  // Keep searching, but only if there are further edge neighbours
5744  // Try next root edge neighbour
5745  i_root_edge_neighbour++;
5746 
5747  // Have we exhausted the search
5748  if (i_root_edge_neighbour >= nroot_edge_neighbour)
5749  {
5750  keep_searching = false;
5751  }
5752 
5753  } // End of while keep searching over all true edge neighbours
5754 
5755  } // End of loop over edges
5756 
5757  // Node not found, return null
5758  return 0;
5759  }
oomph::OcTreeNames::RDF
@ RDF
Definition: octree.h:54
oomph::OcTreeNames::RUB
@ RUB
Definition: octree.h:52
oomph::OcTreeNames::LUF
@ LUF
Definition: octree.h:55
oomph::OcTreeNames::LDF
@ LDF
Definition: octree.h:53
oomph::OcTreeNames::RDB
@ RDB
Definition: octree.h:50
oomph::OcTreeNames::LUB
@ LUB
Definition: octree.h:51
oomph::OcTreeNames::RUF
@ RUF
Definition: octree.h:56
oomph::OcTreeNames::LDB
@ LDB
Definition: octree.h:49

References D, oomph::OcTreeNames::DB, oomph::OcTreeNames::DF, oomph::OcTreeNames::F, oomph::OcTree::faces_of_common_edge(), oomph::FiniteElement::get_node_at_local_coordinate(), oomph::OcTree::gteq_face_neighbour(), oomph::OcTree::gteq_true_edge_neighbour(), i, j, L, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::OcTreeNames::LDB, oomph::OcTreeNames::LDF, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, oomph::OcTreeNames::LUB, oomph::OcTreeNames::LUF, oomph::Tree::nsons(), oomph::Tree::object_pt(), R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, oomph::OcTreeNames::RDB, oomph::OcTreeNames::RDF, oomph::OcTreeNames::RF, oomph::OcTreeNames::RU, oomph::OcTreeNames::RUB, oomph::OcTreeNames::RUF, s, oomph::Tree::son_pt(), RachelsAdvectionDiffusion::U, oomph::OcTreeNames::UB, and oomph::OcTreeNames::UF.

◆ oc_hang_helper()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::oc_hang_helper ( const int &  value_id,
const int &  my_face,
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 mortar method.

/ Set up hanging scheme for this node

/ Print out hanging scheme

/ Print out local projection matrix

/ Print out global projection matrix

/ Print out solved global projection matrix

/ Create structures to hold the hanging info /----------------------------------------—

/ Copy information to hanging nodes /-------------------------------—

/ Print out hanging scheme

/ Print out dependent and master node coords

/ Print out diag(M)

/ Print out diag(M)

/ Print out local projection matrix

/ Print out global projection matrix

/ Print out solved matrix

Reimplemented from oomph::RefineableQElement< 3 >.

8353  {
8354  using namespace OcTreeNames;
8355 
8356  Vector<unsigned> translate_s(3);
8357  Vector<double> s_lo_neigh(3);
8358  Vector<double> s_hi_neigh(3);
8359  int neigh_face, diff_level;
8360  bool in_neighbouring_tree;
8361 
8362  // Find pointer to neighbour in this direction
8363  OcTree* neigh_pt;
8364  neigh_pt = this->octree_pt()->gteq_face_neighbour(my_face,
8365  translate_s,
8366  s_lo_neigh,
8367  s_hi_neigh,
8368  neigh_face,
8369  diff_level,
8370  in_neighbouring_tree);
8371 
8372  // Work out master/dependent faces
8373  //----------------------------
8374 
8375  // Set up booleans
8376  // bool h_type_master = false;
8377  bool h_type_dependent = false;
8378  // bool p_type_master = false;
8379  bool p_type_dependent = false;
8380 
8381  // Neighbour exists
8382  if (neigh_pt != 0)
8383  {
8384  // Check if neighbour is bigger than me
8385  if (diff_level != 0)
8386  {
8387  // Dependent at h-type non-conformity
8388  h_type_dependent = true;
8389  }
8390  // Check if neighbour is the same size as me
8391  else if (neigh_pt->nsons() == 0)
8392  {
8393  // Neighbour is same size as me
8394  // Find p-orders of each element
8395  unsigned my_p_order =
8396  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(this)
8397  ->p_order();
8398  unsigned neigh_p_order =
8399  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(
8400  neigh_pt->object_pt())
8401  ->p_order();
8402 
8403  // Check for p-type non-conformity
8404  if (neigh_p_order == my_p_order)
8405  {
8406  // At a conforming interface
8407  }
8408  else if (neigh_p_order < my_p_order)
8409  {
8410  // Dependent at p-type non-conformity
8411  p_type_dependent = true;
8412  }
8413  else
8414  {
8415  // Master at p-type non-conformity
8416  // p_type_master = true;
8417  }
8418  }
8419  // Neighbour must be smaller than me
8420  else
8421  {
8422  // Master at h-type non-conformity
8423  // h_type_master = true;
8424  }
8425  }
8426  else
8427  {
8428  // Face is on a boundary
8429  }
8430 
8431  // Work out master/dependent edges
8432  //----------------------------
8433  if (h_type_dependent || p_type_dependent || true)
8434  {
8435  // Get edges of face and the face that shares that edge
8436  Vector<unsigned> face_edge(4), face_edge_other_face(4);
8437  switch (my_face)
8438  {
8439  case L:
8440  face_edge[0] = LB;
8441  face_edge_other_face[0] = B;
8442  face_edge[1] = LF;
8443  face_edge_other_face[1] = F;
8444  face_edge[2] = LD;
8445  face_edge_other_face[2] = D;
8446  face_edge[3] = LU;
8447  face_edge_other_face[3] = U;
8448  break;
8449  case R:
8450  face_edge[0] = RB;
8451  face_edge_other_face[0] = B;
8452  face_edge[1] = RF;
8453  face_edge_other_face[1] = F;
8454  face_edge[2] = RD;
8455  face_edge_other_face[2] = D;
8456  face_edge[3] = RU;
8457  face_edge_other_face[3] = U;
8458  break;
8459  case U:
8460  face_edge[0] = UB;
8461  face_edge_other_face[0] = B;
8462  face_edge[1] = UF;
8463  face_edge_other_face[1] = F;
8464  face_edge[2] = LU;
8465  face_edge_other_face[2] = L;
8466  face_edge[3] = RU;
8467  face_edge_other_face[3] = R;
8468  break;
8469  case D:
8470  face_edge[0] = DB;
8471  face_edge_other_face[0] = B;
8472  face_edge[1] = DF;
8473  face_edge_other_face[1] = F;
8474  face_edge[2] = LD;
8475  face_edge_other_face[2] = L;
8476  face_edge[3] = RD;
8477  face_edge_other_face[3] = R;
8478  break;
8479  case B:
8480  face_edge[0] = DB;
8481  face_edge_other_face[0] = D;
8482  face_edge[1] = UB;
8483  face_edge_other_face[1] = U;
8484  face_edge[2] = LB;
8485  face_edge_other_face[2] = L;
8486  face_edge[3] = RB;
8487  face_edge_other_face[3] = R;
8488  break;
8489  case F:
8490  face_edge[0] = DF;
8491  face_edge_other_face[0] = D;
8492  face_edge[1] = UF;
8493  face_edge_other_face[1] = U;
8494  face_edge[2] = LF;
8495  face_edge_other_face[2] = L;
8496  face_edge[3] = RF;
8497  face_edge_other_face[3] = R;
8498  break;
8499  default:
8500  throw OomphLibError(
8501  "my_face not L, R, D, U, B, F\n",
8502  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8503  OOMPH_EXCEPTION_LOCATION);
8504  }
8505 
8506  // Loop over edges of face
8507  for (unsigned i = 0; i < 4; i++)
8508  {
8509  // Get edge
8510  unsigned my_edge = face_edge[i];
8511 
8512  // Separate storage for edge mortaring
8513  OcTree* edge_neigh_pt = 0;
8514  Vector<unsigned> edge_translate_s(3);
8515  Vector<double> edge_s_lo_neigh(3);
8516  Vector<double> edge_s_hi_neigh(3);
8517  int neigh_edge = 0, edge_diff_level = 0;
8518  unsigned edge_p_order =
8519  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(this)
8520  ->p_order();
8521 
8522  // Temporary storage to keep track of master edge
8523  Vector<unsigned> tmp_edge_translate_s(3);
8524  Vector<double> tmp_edge_s_lo_neigh(3);
8525  Vector<double> tmp_edge_s_hi_neigh(3);
8526  int tmp_neigh_edge, tmp_edge_diff_level;
8527 
8528  // Initially return the zero-th true edge neighbour
8529  unsigned i_root_edge_neighbour = 0;
8530 
8531  // Initialise the total number of true edge neighbours
8532  unsigned nroot_edge_neighbour = 0;
8533 
8534  // Flag to keep track of master status of my_edge
8535  bool my_edge_is_master = true;
8536  // unsigned master_edge_index=0;
8537 
8538  // Keep searching until we've found the node or until we've checked
8539  // all available edge neighbours
8540  bool keep_searching = true;
8541  bool search_faces = false;
8542  bool first_face_searched = false;
8543  // unsigned index=0;
8544  while (keep_searching)
8545  {
8546  // Pointer to edge neighbouring OcTree
8547  OcTree* tmp_edge_neigh_pt;
8548 
8549  // Looking in edge neighbours that are also face neighbours
8550  // if(index==0 || index==1)
8551  if (search_faces)
8552  {
8553  if (!first_face_searched)
8554  {
8555  // Find pointer to neighbour in this direction
8556  tmp_edge_neigh_pt =
8557  this->octree_pt()->gteq_face_neighbour(my_face,
8558  tmp_edge_translate_s,
8559  tmp_edge_s_lo_neigh,
8560  tmp_edge_s_hi_neigh,
8561  tmp_neigh_edge,
8562  tmp_edge_diff_level,
8563  in_neighbouring_tree);
8564 
8565  // Mark first face as searched
8566  first_face_searched = true;
8567  }
8568  else
8569  {
8570  // Find pointer to neighbour in this direction
8571  tmp_edge_neigh_pt =
8572  this->octree_pt()->gteq_face_neighbour(face_edge_other_face[i],
8573  tmp_edge_translate_s,
8574  tmp_edge_s_lo_neigh,
8575  tmp_edge_s_hi_neigh,
8576  tmp_neigh_edge,
8577  tmp_edge_diff_level,
8578  in_neighbouring_tree);
8579 
8580  // Search is finally exhausted
8581  keep_searching = false;
8582  }
8583  }
8584  // Looking in a true edge neighbour
8585  else
8586  {
8587  // Find pointer to neighbour in this direction
8588  tmp_edge_neigh_pt =
8589  this->octree_pt()->gteq_true_edge_neighbour(my_edge,
8590  i_root_edge_neighbour,
8591  nroot_edge_neighbour,
8592  tmp_edge_translate_s,
8593  tmp_edge_s_lo_neigh,
8594  tmp_edge_s_hi_neigh,
8595  tmp_neigh_edge,
8596  tmp_edge_diff_level);
8597  }
8598 
8599  // Set up booleans
8600  // bool h_type_edge_master = false;
8601  // bool h_type_edge_dependent = false;
8602  // bool p_type_edge_master = false;
8603  // bool p_type_edge_dependent = false;
8604 
8605  // Flag to check if we have a new edge master
8606  bool new_edge_master = false;
8607 
8608  // Edge neighbour exists
8609  if (tmp_edge_neigh_pt != 0)
8610  {
8611  // Check if neighbour is bigger than my biggest neighbour
8612  if (tmp_edge_diff_level < edge_diff_level)
8613  {
8614  // Dependent at h-type non-conformity
8615  // h_type_edge_dependent = true;
8616  new_edge_master = true;
8617  // Update edge_diff_level and p-order
8618  edge_diff_level = tmp_edge_diff_level;
8619  edge_p_order =
8620  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(
8621  tmp_edge_neigh_pt->object_pt())
8622  ->p_order();
8623  }
8624  // Check if neighbour is the same size as my biggest neighbour
8625  else if (tmp_edge_diff_level == edge_diff_level &&
8626  tmp_edge_neigh_pt->nsons() == 0)
8627  {
8628  // Neighbour is same size as me
8629  // Find p-orders of each element
8630  // unsigned my_p_order =
8631  // dynamic_cast<PRefineableQElement<3,INITIAL_NNODE_1D>*>
8632  // (this)->p_order();
8633  unsigned tmp_edge_neigh_p_order =
8634  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(
8635  tmp_edge_neigh_pt->object_pt())
8636  ->p_order();
8637 
8638  // Check for p-type non-conformity
8639  if (tmp_edge_neigh_p_order == edge_p_order)
8640  {
8641  // At a conforming interface
8642  }
8643  else if (tmp_edge_neigh_p_order < edge_p_order)
8644  {
8645  // Dependent at p-type non-conformity
8646  // p_type_edge_dependent = true;
8647  new_edge_master = true;
8648  // Update edge_diff_level and p-order
8649  edge_diff_level = tmp_edge_diff_level;
8650  edge_p_order =
8651  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(
8652  tmp_edge_neigh_pt->object_pt())
8653  ->p_order();
8654  }
8655  else
8656  {
8657  // Master at p-type non-conformity
8658  // p_type_edge_master = true;
8659  }
8660  }
8661  // Neighbour must be smaller than me
8662  else
8663  {
8664  // Master at h-type non-conformity
8665  // h_type_edge_master = true;
8666  }
8667  }
8668  else
8669  {
8670  // Edge is on a boundary
8671  }
8672 
8673  // Update master neighbour information
8674  if (new_edge_master)
8675  {
8676  // Store new master edge information
8677  edge_neigh_pt = tmp_edge_neigh_pt;
8678  edge_translate_s = tmp_edge_translate_s;
8679  edge_s_lo_neigh = tmp_edge_s_lo_neigh;
8680  edge_s_hi_neigh = tmp_edge_s_hi_neigh;
8681  neigh_edge = tmp_neigh_edge;
8682  edge_diff_level = tmp_edge_diff_level;
8683 
8684  // Set this edge as dependent
8685  my_edge_is_master = false;
8686  }
8687 
8688  // Keep searching, but only if there are further edge neighbours
8689  // Try next root edge neighbour
8690  i_root_edge_neighbour++;
8691 
8692  // Increment counter
8693  // index++;
8694 
8695  // Have we exhausted the search over true edge neighbours
8696  // if (index>2 && i_root_edge_neighbour>=nroot_edge_neighbour)
8697  if (i_root_edge_neighbour >= nroot_edge_neighbour)
8698  {
8699  // keep_searching = false;
8700  // Extend search to face neighours (these are edge neighbours too)
8701  search_faces = true;
8702  }
8703 
8704  } // End of while keep searching over all face and true edge neighbours
8705 
8706  // Now do edge mortaring to enforce the mortar vertex matching condition
8707  if (!my_edge_is_master)
8708  {
8709  // Compute the active coordinate index along the this side of mortar
8710  unsigned active_coord_index;
8711  switch (my_edge)
8712  {
8713  case DB:
8714  case DF:
8715  case UB:
8716  case UF:
8717  active_coord_index = 0;
8718  break;
8719  case LB:
8720  case RB:
8721  case LF:
8722  case RF:
8723  active_coord_index = 1;
8724  break;
8725  case LD:
8726  case RD:
8727  case LU:
8728  case RU:
8729  active_coord_index = 2;
8730  break;
8731  default:
8732  throw OomphLibError(
8733  "Cannot transform coordinates",
8734  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
8735  OOMPH_EXCEPTION_LOCATION);
8736  }
8737 
8738  // Get pointer to neighbouring master element (in p-refineable form)
8739  PRefineableQElement<3, INITIAL_NNODE_1D>* edge_neigh_obj_pt;
8740  edge_neigh_obj_pt =
8741  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(
8742  edge_neigh_pt->object_pt());
8743 
8744  // Create vector of master and dependent nodes
8745  //----------------------------------------
8746  Vector<Node*> master_node_pt, dependent_node_pt;
8747  Vector<unsigned> master_node_number, dependent_node_number;
8748  Vector<Vector<double>> dependent_node_s_fraction;
8749 
8750  // Number of nodes in one dimension
8751  const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
8752  const unsigned neigh_n_p =
8753  edge_neigh_obj_pt->ninterpolating_node_1d(value_id);
8754 
8755  // Storage for pointers to the nodes and their numbers along the
8756  // master edge
8757  unsigned neighbour_node_number = 0;
8758  Node* neighbour_node_pt = 0;
8759 
8760  // Loop over nodes along the edge
8761  bool master_is_not_edge = false;
8762  for (unsigned i0 = 0; i0 < neigh_n_p; i0++)
8763  {
8764  const unsigned s0space = 1;
8765  const unsigned s1space = neigh_n_p;
8766  const unsigned s2space = neigh_n_p * neigh_n_p;
8767 
8768  // Find the neighbour's node
8769  switch (neigh_edge)
8770  {
8771  case DB:
8772  neighbour_node_number = i0 * s0space;
8773  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8774  neighbour_node_number, value_id);
8775  break;
8776  case UB:
8777  neighbour_node_number =
8778  (neigh_n_p - 1) * s1space + i0 * s0space;
8779  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8780  neighbour_node_number, value_id);
8781  break;
8782  case DF:
8783  neighbour_node_number =
8784  (neigh_n_p - 1) * s2space + i0 * s0space;
8785  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8786  neighbour_node_number, value_id);
8787  break;
8788  case UF:
8789  neighbour_node_number = (neigh_n_p - 1) * s1space +
8790  (neigh_n_p - 1) * s2space +
8791  i0 * s0space;
8792  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8793  neighbour_node_number, value_id);
8794  break;
8795 
8796  case LB:
8797  neighbour_node_number = i0 * s1space;
8798  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8799  neighbour_node_number, value_id);
8800  break;
8801  case RB:
8802  neighbour_node_number =
8803  (neigh_n_p - 1) * s0space + i0 * s1space;
8804  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8805  neighbour_node_number, value_id);
8806  break;
8807  case LF:
8808  neighbour_node_number =
8809  (neigh_n_p - 1) * s2space + i0 * s1space;
8810  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8811  neighbour_node_number, value_id);
8812  break;
8813  case RF:
8814  neighbour_node_number = (neigh_n_p - 1) * s0space +
8815  (neigh_n_p - 1) * s2space +
8816  i0 * s1space;
8817  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8818  neighbour_node_number, value_id);
8819  break;
8820 
8821  case LD:
8822  neighbour_node_number = i0 * s2space;
8823  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8824  neighbour_node_number, value_id);
8825  break;
8826  case RD:
8827  neighbour_node_number =
8828  (neigh_n_p - 1) * s0space + i0 * s2space;
8829  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8830  neighbour_node_number, value_id);
8831  break;
8832  case LU:
8833  neighbour_node_number =
8834  (neigh_n_p - 1) * s1space + i0 * s2space;
8835  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8836  neighbour_node_number, value_id);
8837  break;
8838  case RU:
8839  neighbour_node_number = (neigh_n_p - 1) * s0space +
8840  (neigh_n_p - 1) * s1space +
8841  i0 * s2space;
8842  neighbour_node_pt = edge_neigh_obj_pt->interpolating_node_pt(
8843  neighbour_node_number, value_id);
8844  break;
8845 
8846  default:
8847  // Master 'edge' may be a face instead, so no need to throw an
8848  // error
8849  master_is_not_edge = true;
8850  }
8851 
8852  if (master_is_not_edge) break;
8853 
8854  // Set node as master node
8855  master_node_number.push_back(neighbour_node_number);
8856  master_node_pt.push_back(neighbour_node_pt);
8857  }
8858 
8859  // Now if edge is really a face (from an edge neighbour that isn't
8860  // a true edge neighbour) each node on the face is (potentially) a
8861  // master
8862  if (master_is_not_edge)
8863  {
8864  // Loop over nodes along the face
8865  for (unsigned i0 = 0; i0 < neigh_n_p; i0++)
8866  {
8867  // Loop over nodes along the face
8868  for (unsigned i1 = 0; i1 < neigh_n_p; i1++)
8869  {
8870  const unsigned s0space = 1;
8871  const unsigned s1space = neigh_n_p;
8872  const unsigned s2space = neigh_n_p * neigh_n_p;
8873 
8874  // Find the neighbour's node
8875  switch (neigh_edge)
8876  {
8877  case B:
8878  neighbour_node_number = i0 * s0space + i1 * s1space;
8879  neighbour_node_pt =
8880  edge_neigh_obj_pt->interpolating_node_pt(
8881  neighbour_node_number, value_id);
8882  break;
8883  case F:
8884  neighbour_node_number =
8885  (neigh_n_p - 1) * s2space + i0 * s0space + i1 * s1space;
8886  neighbour_node_pt =
8887  edge_neigh_obj_pt->interpolating_node_pt(
8888  neighbour_node_number, value_id);
8889  break;
8890  case D:
8891  neighbour_node_number = i0 * s0space + i1 * s2space;
8892  neighbour_node_pt =
8893  edge_neigh_obj_pt->interpolating_node_pt(
8894  neighbour_node_number, value_id);
8895  break;
8896  case U:
8897  neighbour_node_number =
8898  (neigh_n_p - 1) * s1space + i0 * s0space + i1 * s2space;
8899  neighbour_node_pt =
8900  edge_neigh_obj_pt->interpolating_node_pt(
8901  neighbour_node_number, value_id);
8902  break;
8903  case L:
8904  neighbour_node_number = i0 * s1space + i1 * s2space;
8905  neighbour_node_pt =
8906  edge_neigh_obj_pt->interpolating_node_pt(
8907  neighbour_node_number, value_id);
8908  break;
8909  case R:
8910  neighbour_node_number =
8911  (neigh_n_p - 1) * s0space + i0 * s1space + i1 * s2space;
8912  neighbour_node_pt =
8913  edge_neigh_obj_pt->interpolating_node_pt(
8914  neighbour_node_number, value_id);
8915  break;
8916 
8917  default:
8918  throw OomphLibError("neigh_edge not recognised\n",
8919  "PRefineableQElement<3,INITIAL_NNODE_"
8920  "1D>::oc_hang_helper()",
8921  OOMPH_EXCEPTION_LOCATION);
8922  }
8923 
8924  // Set node as master node
8925  master_node_number.push_back(neighbour_node_number);
8926  master_node_pt.push_back(neighbour_node_pt);
8927  }
8928  }
8929  }
8930 
8931  // Storage for pointers to the local nodes and their numbers along my
8932  // edge
8933  unsigned local_node_number = 0;
8934  Node* local_node_pt = 0;
8935 
8936  // Loop over the nodes along my edge
8937  for (unsigned i0 = 0; i0 < my_n_p; i0++)
8938  {
8939  // Storage for the fractional position of the node in the element
8940  Vector<double> s_fraction(3);
8941 
8942  const unsigned s0space = 1;
8943  const unsigned s1space = my_n_p;
8944  const unsigned s2space = my_n_p * my_n_p;
8945 
8946  // Find the local node and the fractional position of the node
8947  // which depends on the edge, of course
8948  switch (my_edge)
8949  {
8950  case DB:
8951  s_fraction[0] =
8952  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
8953  s_fraction[1] = 0.0;
8954  s_fraction[2] = 0.0;
8955  local_node_number = i0 * s0space;
8956  local_node_pt =
8957  this->interpolating_node_pt(local_node_number, value_id);
8958  break;
8959  case UB:
8960  s_fraction[0] =
8961  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
8962  s_fraction[1] = 1.0;
8963  s_fraction[2] = 0.0;
8964  local_node_number = (my_n_p - 1) * s1space + i0 * s0space;
8965  local_node_pt =
8966  this->interpolating_node_pt(local_node_number, value_id);
8967  break;
8968  case DF:
8969  s_fraction[0] =
8970  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
8971  s_fraction[1] = 0.0;
8972  s_fraction[2] = 1.0;
8973  local_node_number = (my_n_p - 1) * s2space + i0 * s0space;
8974  local_node_pt =
8975  this->interpolating_node_pt(local_node_number, value_id);
8976  break;
8977  case UF:
8978  s_fraction[0] =
8979  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
8980  s_fraction[1] = 1.0;
8981  s_fraction[2] = 1.0;
8982  local_node_number = (my_n_p - 1) * s1space +
8983  (my_n_p - 1) * s2space + i0 * s0space;
8984  local_node_pt =
8985  this->interpolating_node_pt(local_node_number, value_id);
8986  break;
8987 
8988  case LB:
8989  s_fraction[0] = 0.0;
8990  s_fraction[1] =
8991  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
8992  s_fraction[2] = 0.0;
8993  local_node_number = i0 * s1space;
8994  local_node_pt =
8995  this->interpolating_node_pt(local_node_number, value_id);
8996  break;
8997  case RB:
8998  s_fraction[0] = 1.0;
8999  s_fraction[1] =
9000  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9001  s_fraction[2] = 0.0;
9002  local_node_number = (my_n_p - 1) * s0space + i0 * s1space;
9003  local_node_pt =
9004  this->interpolating_node_pt(local_node_number, value_id);
9005  break;
9006  case LF:
9007  s_fraction[0] = 0.0;
9008  s_fraction[1] =
9009  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9010  s_fraction[2] = 1.0;
9011  local_node_number = (my_n_p - 1) * s2space + i0 * s1space;
9012  local_node_pt =
9013  this->interpolating_node_pt(local_node_number, value_id);
9014  break;
9015  case RF:
9016  s_fraction[0] = 1.0;
9017  s_fraction[1] =
9018  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9019  s_fraction[2] = 1.0;
9020  local_node_number = (my_n_p - 1) * s0space +
9021  (my_n_p - 1) * s2space + i0 * s1space;
9022  local_node_pt =
9023  this->interpolating_node_pt(local_node_number, value_id);
9024  break;
9025 
9026  case LD:
9027  s_fraction[0] = 0.0;
9028  s_fraction[1] = 0.0;
9029  s_fraction[2] =
9030  local_one_d_fraction_of_interpolating_node(i0, 2, value_id);
9031  local_node_number = i0 * s2space;
9032  local_node_pt =
9033  this->interpolating_node_pt(local_node_number, value_id);
9034  break;
9035  case RD:
9036  s_fraction[0] = 1.0;
9037  s_fraction[1] = 0.0;
9038  s_fraction[2] =
9039  local_one_d_fraction_of_interpolating_node(i0, 2, value_id);
9040  local_node_number = (my_n_p - 1) * s0space + i0 * s2space;
9041  local_node_pt =
9042  this->interpolating_node_pt(local_node_number, value_id);
9043  break;
9044  case LU:
9045  s_fraction[0] = 0.0;
9046  s_fraction[1] = 1.0;
9047  s_fraction[2] =
9048  local_one_d_fraction_of_interpolating_node(i0, 2, value_id);
9049  local_node_number = (my_n_p - 1) * s1space + i0 * s2space;
9050  local_node_pt =
9051  this->interpolating_node_pt(local_node_number, value_id);
9052  break;
9053  case RU:
9054  s_fraction[0] = 1.0;
9055  s_fraction[1] = 1.0;
9056  s_fraction[2] =
9057  local_one_d_fraction_of_interpolating_node(i0, 2, value_id);
9058  local_node_number = (my_n_p - 1) * s0space +
9059  (my_n_p - 1) * s1space + i0 * s2space;
9060  local_node_pt =
9061  this->interpolating_node_pt(local_node_number, value_id);
9062  break;
9063 
9064  default:
9065  throw OomphLibError(
9066  "my_edge not recognised\n",
9067  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9068  OOMPH_EXCEPTION_LOCATION);
9069  }
9070 
9071  // Add node to vector of dependent element nodes
9072  dependent_node_number.push_back(local_node_number);
9073  dependent_node_pt.push_back(local_node_pt);
9074 
9075  // Store node's local fraction
9076  dependent_node_s_fraction.push_back(s_fraction);
9077  }
9078 
9079  // Store the number of dependent and master nodes for use later
9080  const unsigned n_dependent_nodes = dependent_node_pt.size();
9081  const unsigned n_master_nodes = master_node_pt.size();
9082  const unsigned dependent_element_nnode_1d = my_n_p;
9083  const unsigned master_element_nnode_1d = neigh_n_p;
9084 
9085  // Storage for master shapes
9086  Shape master_shapes(edge_neigh_obj_pt->ninterpolating_node(value_id));
9087 
9088  // Get master and dependent nodal positions
9089  //-------------------------------------
9090  Vector<double> dependent_nodal_position;
9091  Vector<double> dependent_weight(dependent_element_nnode_1d);
9092  Orthpoly::gll_nodes(dependent_element_nnode_1d,
9093  dependent_nodal_position,
9094  dependent_weight);
9095  Vector<double> master_nodal_position;
9096  Vector<double> master_weight(master_element_nnode_1d);
9097  Orthpoly::gll_nodes(
9098  master_element_nnode_1d, master_nodal_position, master_weight);
9099 
9100  // Apply the (1D) vertex matching condition
9101  //-----------------------------------------
9102  // Vertiex matching is ensured automatically in cases where there is a
9103  // node at each end of the mortar that is shared between the master
9104  // and dependent elements. Where this is not the case, the vertex
9105  // matching condition must be enforced by constraining the value of
9106  // the unknown at the node on the dependent side to be the same as the
9107  // value at that location in the master.
9108 
9109  // Store positions of the mortar vertex/non-vertex nodes in the
9110  // dependent element
9111  const unsigned n_mortar_vertices = 2;
9112  Vector<unsigned> vertex_pos(n_mortar_vertices);
9113  vertex_pos[0] = 0;
9114  vertex_pos[1] = this->ninterpolating_node_1d(value_id) - 1;
9115  Vector<unsigned> non_vertex_pos(my_n_p - n_mortar_vertices);
9116  for (unsigned i = 0; i < my_n_p - n_mortar_vertices; i++)
9117  {
9118  non_vertex_pos[i] = i + 1;
9119  }
9120 
9121  // If the node is on a master edge, we may be setting the
9122  // hanging info incorrectly. Hanging schemes (if they are
9123  // required) for such nodes are instead computed by the
9124  // dependent edge. We dont't want to overwrite them here!
9125  std::vector<bool> mortar_vertex_on_master_edge(n_mortar_vertices);
9126  for (unsigned v = 0; v < n_mortar_vertices; v++)
9127  {
9128  // Check if each node is on the master edge
9129  mortar_vertex_on_master_edge[v] = true;
9130  unsigned non_extreme_coordinate = 0;
9131  Vector<double> s_in_neigh(3);
9132  for (unsigned i = 0; i < 3; i++)
9133  {
9134  // Work out this node's location in the master
9135  s_in_neigh[i] =
9136  edge_s_lo_neigh[i] +
9137  dependent_node_s_fraction[vertex_pos[v]][edge_translate_s[i]] *
9138  (edge_s_hi_neigh[i] - edge_s_lo_neigh[i]);
9139 
9140  // Check if local coordinate in master element takes non-extreme
9141  // value
9142  if (std::fabs(std::fabs(s_in_neigh[i]) - 1.0) > 1.0e-14)
9143  {
9144  non_extreme_coordinate++;
9145  if (non_extreme_coordinate > 1)
9146  {
9147  mortar_vertex_on_master_edge[v] = false;
9148  break;
9149  }
9150  }
9151  }
9152  }
9153 
9154  // Now work out if my edge coincides with the master edge
9155  bool my_edge_coincides_with_master = true;
9156  for (unsigned v = 0; v < n_mortar_vertices; v++)
9157  {
9158  my_edge_coincides_with_master =
9159  my_edge_coincides_with_master && mortar_vertex_on_master_edge[v];
9160  }
9161 
9162  // Check if we need to apply the (1D) vertex matching condition at the
9163  // mortar vertices. This is trivially satisfied if the node is shared
9164  // with the master, but if not then we need to constrain it.
9165  for (unsigned v = 0; v < n_mortar_vertices; v++)
9166  {
9167  // Don't make hanging node if my edge doesn't coincide with
9168  // the master edge *and* this node is on the master edge!
9169  // (We are not in a position to determine its hanging status.)
9170  if ((!my_edge_coincides_with_master) &&
9171  mortar_vertex_on_master_edge[v])
9172  continue;
9173 
9174  // Search master node storage for the node
9175  bool node_is_shared = false;
9176  for (unsigned i = 0; i < master_node_pt.size(); i++)
9177  {
9178  if (dependent_node_pt[vertex_pos[v]] == master_node_pt[i])
9179  {
9180  node_is_shared = true;
9181  break;
9182  }
9183  }
9184 
9185  // If the node is not shared then we must constrain its value by
9186  // setting up a hanging scheme
9187  if (!node_is_shared)
9188  {
9189  // Calculate weights. These are just the master shapes evaluated
9190  // at this dependent node's position
9191 
9192  // Work out this node's location in the master
9193  Vector<double> s_in_neigh(3);
9194  for (unsigned i = 0; i < 3; i++)
9195  {
9196  s_in_neigh[i] = edge_s_lo_neigh[i] +
9197  dependent_node_s_fraction[vertex_pos[v]]
9198  [edge_translate_s[i]] *
9199  (edge_s_hi_neigh[i] - edge_s_lo_neigh[i]);
9200  }
9201 
9202  // Get master shapes at dependent nodal positions
9203  edge_neigh_obj_pt->interpolating_basis(
9204  s_in_neigh, master_shapes, value_id);
9205 
9206  // Remove any existing hanging node info
9207  // (This may be a bit wasteful, but guarantees correctness)
9208  dependent_node_pt[vertex_pos[v]]->set_nonhanging();
9209 
9210  // Don't include master nodes with zero weights
9211  Vector<unsigned> master_node_zero_weight;
9212  for (unsigned m = 0; m < n_master_nodes; m++)
9213  {
9214  // Compare weights to some (small) tolerance
9215  if (std::fabs(master_shapes[master_node_number[m]]) < 1.0e-14)
9216  {
9217  // Store
9218  master_node_zero_weight.push_back(m);
9219  }
9220  }
9221 
9222  // Set up hanging scheme for this node
9223  HangInfo* explicit_hang_pt =
9224  new HangInfo(n_master_nodes - master_node_zero_weight.size());
9225  unsigned mindex = 0;
9226  for (unsigned m = 0; m < n_master_nodes; m++)
9227  {
9228  // Check that master doesn't have zero weight
9229  bool skip = false;
9230  for (unsigned i = 0; i < master_node_zero_weight.size(); i++)
9231  {
9232  if (m == master_node_zero_weight[i]) skip = true;
9233  }
9234 
9235  // Add pointer and weight to hang info
9236  if (!skip)
9237  explicit_hang_pt->set_master_node_pt(
9238  mindex++,
9239  master_node_pt[m],
9240  master_shapes[master_node_number[m]]);
9241  }
9242 
9244  // HangInfo* explicit_hang_pt = new HangInfo(n_master_nodes);
9245  // for(unsigned m=0; m<n_master_nodes; m++)
9246  // {
9247  // explicit_hang_pt->set_master_node_pt(m,master_node_pt[m],master_shapes[master_node_number[m]]);
9248  // }
9249 
9250  // Make node hang
9251  dependent_node_pt[vertex_pos[v]]->set_hanging_pt(explicit_hang_pt,
9252  -1);
9253 
9255  // std::cout << "Hanging node: "
9256  // << dependent_node_pt[vertex_pos[v]]->x(0) << " "
9257  // << dependent_node_pt[vertex_pos[v]]->x(1) << " "
9258  // << dependent_node_pt[vertex_pos[v]]->x(2) << " "
9259  // << std::endl;
9260  // for(unsigned m=0; m<explicit_hang_pt->nmaster(); m++)
9261  // {
9262  // std::cout << " m = " << m << ": "
9263  // << explicit_hang_pt->master_node_pt(m)->x(0) << " "
9264  // << explicit_hang_pt->master_node_pt(m)->x(1) << " "
9265  // << explicit_hang_pt->master_node_pt(m)->x(2) << " "
9266  // << "w = " << explicit_hang_pt->master_weight(m) << "
9267  // "
9268  // << std::endl;
9269  // }
9270 
9271  // Check there are no zero weights at this stage
9272  for (unsigned m = 0; m < explicit_hang_pt->nmaster(); m++)
9273  {
9274  if (std::fabs(explicit_hang_pt->master_weight(m)) < 1.0e-14)
9275  {
9276  throw OomphLibError(
9277  "Master has zero weight!",
9278  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9279  OOMPH_EXCEPTION_LOCATION);
9280  }
9281  }
9282  }
9283  }
9284 
9285  // Check that there are actually dependent nodes for which we still
9286  // need to construct a hanging scheme. If not then there is nothing
9287  // more to do.
9288  if (n_dependent_nodes - n_mortar_vertices > 0)
9289  {
9290  // Assemble mass matrix for mortar
9291  //--------------------------------
9292  Vector<double> psi(n_dependent_nodes - n_mortar_vertices);
9293  Vector<double> diag_M(n_dependent_nodes - n_mortar_vertices);
9294  Vector<Vector<double>> shared_node_M(n_mortar_vertices);
9295  for (unsigned i = 0; i < shared_node_M.size(); i++)
9296  {
9297  shared_node_M[i].resize(n_dependent_nodes - n_mortar_vertices);
9298  }
9299 
9300  // Fill in part corresponding to dependent nodal positions (unknown)
9301  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
9302  {
9303  // Use L'Hosptal's rule:
9304  psi[i] =
9305  pow(-1.0, int((dependent_element_nnode_1d - 1) - i - 1)) *
9306  -Orthpoly::ddlegendre(
9307  dependent_element_nnode_1d - 1,
9308  dependent_nodal_position[non_vertex_pos[i]]);
9309  // Put in contribution
9310  diag_M[i] = psi[i] * dependent_weight[non_vertex_pos[i]];
9311  }
9312 
9313  // Fill in part corresponding to dependent element vertices (known)
9314  for (unsigned v = 0; v < shared_node_M.size(); v++)
9315  {
9316  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices;
9317  i++)
9318  {
9319  // Check if denominator is zero
9320  if (std::fabs(dependent_nodal_position[non_vertex_pos[i]] -
9321  dependent_nodal_position[vertex_pos[v]]) >=
9322  1.0e-8)
9323  {
9324  // We're ok
9325  psi[i] =
9326  pow(-1.0, int((dependent_element_nnode_1d - 1) - i - 1)) *
9327  Orthpoly::dlegendre(
9328  dependent_element_nnode_1d - 1,
9329  dependent_nodal_position[vertex_pos[v]]) /
9330  (dependent_nodal_position[non_vertex_pos[i]] -
9331  dependent_nodal_position[vertex_pos[v]]);
9332  }
9333  // Check if numerator is zero
9334  else if (std::fabs(Orthpoly::dlegendre(
9335  dependent_element_nnode_1d - 1,
9336  dependent_nodal_position[vertex_pos[v]])) < 1.0e-8)
9337  {
9338  // We can use l'hopital's rule
9339  psi[i] =
9340  pow(-1.0, int((dependent_element_nnode_1d - 1) - i - 1)) *
9341  -Orthpoly::ddlegendre(
9342  dependent_element_nnode_1d - 1,
9343  dependent_nodal_position[non_vertex_pos[i]]);
9344  }
9345  else
9346  {
9347  // We can't use l'hopital's rule
9348  throw OomphLibError(
9349  "Cannot use l'Hopital's rule. Dividing by zero is not "
9350  "allowed!",
9351  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9352  OOMPH_EXCEPTION_LOCATION);
9353  }
9354  // Put in contribution
9355  shared_node_M[v][i] = psi[i] * dependent_weight[vertex_pos[v]];
9356  }
9357  }
9358 
9359  // Assemble local projection matrix for mortar
9360  //--------------------------------------------
9361 
9362  // Have only one local projection matrix because there is just one
9363  // master
9364  Vector<Vector<double>> P(n_dependent_nodes - n_mortar_vertices);
9365  for (unsigned i = 0; i < P.size(); i++)
9366  {
9367  P[i].resize(n_master_nodes, 0.0);
9368  }
9369 
9370  // Storage for local coordinate
9371  Vector<double> s(3);
9372 
9373  // Sum contributions from master element shapes (quadrature).
9374  // The order of quadrature must be high enough to evaluate a
9375  // polynomial of order N_s + N_m - 3 exactly, where N_s =
9376  // n_dependent_nodes, N_m = n_master_nodes. (Use pointers for the
9377  // quadrature knots and weights so that data is not unnecessarily
9378  // copied)
9379  // unsigned quadrature_order =
9380  // std::max(dependent_element_nnode_1d,master_element_nnode_1d);
9381  Vector<double>*quadrature_knot, *quadrature_weight;
9382  if (dependent_element_nnode_1d >= master_element_nnode_1d)
9383  {
9384  // Use the same quadrature order as the dependent element (me)
9385  quadrature_knot = &dependent_nodal_position;
9386  quadrature_weight = &dependent_weight;
9387  }
9388  else
9389  {
9390  // Use the same quadrature order as the master element (neighbour)
9391  quadrature_knot = &master_nodal_position;
9392  quadrature_weight = &master_weight;
9393  }
9394 
9395  // Quadrature loop
9396  for (unsigned q = 0; q < (*quadrature_weight).size(); q++)
9397  {
9398  // Evaluate mortar test functions at the quadrature knots in the
9399  // dependent
9400  // s[active_coord_index] = (*quadrature_knot)[q];
9401  Vector<double> s_on_mortar(1);
9402  s_on_mortar[0] = (*quadrature_knot)[q];
9403 
9404  // Get psi
9405  for (unsigned k = 0; k < n_dependent_nodes - n_mortar_vertices;
9406  k++)
9407  {
9408  // Check if denominator is zero
9409  if (std::fabs(dependent_nodal_position[non_vertex_pos[k]] -
9410  s_on_mortar[0]) >= 1.0e-08)
9411  {
9412  // We're ok
9413  psi[k] =
9414  pow(-1.0, int((dependent_element_nnode_1d - 1) - k - 1)) *
9415  Orthpoly::dlegendre(dependent_element_nnode_1d - 1,
9416  s_on_mortar[0]) /
9417  (dependent_nodal_position[non_vertex_pos[k]] -
9418  s_on_mortar[0]);
9419  }
9420  // Check if numerator is zero
9421  else if (std::fabs(Orthpoly::dlegendre(
9422  dependent_element_nnode_1d - 1, s_on_mortar[0])) <
9423  1.0e-8)
9424  {
9425  // We can use l'Hopital's rule
9426  psi[k] =
9427  pow(-1.0, int((dependent_element_nnode_1d - 1) - k - 1)) *
9428  -Orthpoly::ddlegendre(dependent_element_nnode_1d - 1,
9429  s_on_mortar[0]);
9430  }
9431  else
9432  {
9433  // We can't use l'hopital's rule
9434  throw OomphLibError(
9435  "Cannot use l'Hopital's rule. Dividing by zero is not "
9436  "allowed!",
9437  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9438  OOMPH_EXCEPTION_LOCATION);
9439  }
9440  }
9441 
9442  // Convert coordinate on mortar to local fraction in dependent
9443  // element
9444  Vector<double> s_fraction(3);
9445  for (unsigned i = 0; i < 3; i++)
9446  {
9447  s_fraction[i] = (i == active_coord_index) ?
9448  0.5 * (s_on_mortar[0] + 1.0) :
9449  dependent_node_s_fraction[vertex_pos[0]][i];
9450  }
9451 
9452  // Project active coordinate into master element
9453  Vector<double> s_in_neigh(3);
9454  for (unsigned i = 0; i < 3; i++)
9455  {
9456  s_in_neigh[i] = edge_s_lo_neigh[i] +
9457  s_fraction[edge_translate_s[i]] *
9458  (edge_s_hi_neigh[i] - edge_s_lo_neigh[i]);
9459  }
9460 
9461  // Evaluate master shapes at projections of local quadrature knots
9462  edge_neigh_obj_pt->interpolating_basis(
9463  s_in_neigh, master_shapes, value_id);
9464 
9465  // Populate local projection matrix
9466  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices;
9467  i++)
9468  {
9469  for (unsigned j = 0; j < n_master_nodes; j++)
9470  {
9471  P[i][j] += master_shapes[master_node_number[j]] * psi[i] *
9472  (*quadrature_weight)[q];
9473  }
9474  }
9475  }
9476 
9478  // std::cout << "P_e:" << std::endl;
9479  // for(unsigned i=0; i<P.size(); i++)
9480  // {
9481  // for(unsigned j=0; j<P[i].size(); j++)
9482  // {
9483  // std::cout << " " << P[i][j];
9484  // }
9485  // }
9486  // std::cout << std::endl;
9487 
9488  // Assemble global projection matrices for mortar
9489  //-----------------------------------------------
9490  // Need to subtract contributions from the "known unknowns"
9491  // corresponding to the nodes at the vertices of the mortar
9492 
9493  // Assemble contributions from mortar vertex nodes
9494  for (unsigned v = 0; v < n_mortar_vertices; v++)
9495  {
9496  // Convert coordinate on mortar to local fraction in dependent
9497  // element
9498  Vector<double> s_fraction(3);
9499  for (unsigned i = 0; i < 3; i++)
9500  {
9501  s_fraction[i] =
9502  (i == active_coord_index) ?
9503  0.5 * (dependent_nodal_position[vertex_pos[v]] + 1.0) :
9504  dependent_node_s_fraction[vertex_pos[0]][i];
9505  }
9506 
9507  // Project active coordinate into master element
9508  Vector<double> s_in_neigh(3);
9509  for (unsigned i = 0; i < 3; i++)
9510  {
9511  s_in_neigh[i] = edge_s_lo_neigh[i] +
9512  s_fraction[edge_translate_s[i]] *
9513  (edge_s_hi_neigh[i] - edge_s_lo_neigh[i]);
9514  }
9515 
9516  // Get master shapes at dependent nodal positions
9517  edge_neigh_obj_pt->interpolating_basis(
9518  s_in_neigh, master_shapes, value_id);
9519 
9520  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices;
9521  i++)
9522  {
9523  for (unsigned k = 0; k < n_master_nodes; k++)
9524  {
9525  P[i][k] -=
9526  master_shapes[master_node_number[k]] * shared_node_M[v][i];
9527  }
9528  }
9529  }
9530 
9532  // std::cout << "P:" << std::endl;
9533  // for(unsigned i=0; i<P.size(); i++)
9534  // {
9535  // for(unsigned j=0; j<P[i].size(); j++)
9536  // {
9537  // std::cout << " " << P[i][j];
9538  // }
9539  // std::cout << std::endl;
9540  // }
9541  // std::cout << std::endl;
9542 
9543  // Solve mortar system
9544  //--------------------
9545  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
9546  {
9547  for (unsigned j = 0; j < n_master_nodes; j++)
9548  {
9549  P[i][j] /= diag_M[i];
9550  }
9551  }
9552 
9554  // std::cout << "solved P:" << std::endl;
9555  // for(unsigned i=0; i<P.size(); i++)
9556  // {
9557  // for(unsigned j=0; j<P[i].size(); j++)
9558  // {
9559  // std::cout << " " << P[i][j];
9560  // }
9561  // std::cout << std::endl;
9562  // }
9563  // std::cout << std::endl;
9564 
9565  // Create and populate structures to hold the hanging info
9566  //--------------------------------------------------------
9567  Vector<HangInfo*> hang_info_pt(n_dependent_nodes);
9568  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
9569  {
9570  // Don't include master nodes with zero weights
9571  Vector<unsigned> master_node_zero_weight;
9572  for (unsigned m = 0; m < n_master_nodes; m++)
9573  {
9574  // Compare weights to some (small) tolerance
9575  if (std::fabs(P[i][m]) < 1.0e-14)
9576  {
9577  // Store
9578  master_node_zero_weight.push_back(m);
9579  }
9580  }
9581 
9582  // Set up hanging scheme for this node
9583  hang_info_pt[i] =
9584  new HangInfo(n_master_nodes - master_node_zero_weight.size());
9585  unsigned mindex = 0;
9586  for (unsigned m = 0; m < n_master_nodes; m++)
9587  {
9588  // Check that master doesn't have zero weight
9589  bool skip = false;
9590  for (unsigned k = 0; k < master_node_zero_weight.size(); k++)
9591  {
9592  if (m == master_node_zero_weight[k]) skip = true;
9593  }
9594 
9595  // Add pointer and weight to hang info
9596  if (!skip)
9597  hang_info_pt[i]->set_master_node_pt(
9598  mindex++, master_node_pt[m], P[i][m]);
9599  }
9600  }
9601 
9604  // Vector<HangInfo*> hang_info_pt(n_dependent_nodes);
9605  // for (unsigned i=0; i<n_dependent_nodes-n_mortar_vertices; i++)
9606  // {
9607  // hang_info_pt[i] = new HangInfo(n_master_nodes);
9608  // }
9609  //
9612  // for(unsigned i=0; i<n_dependent_nodes-n_mortar_vertices; i++)
9613  // {
9614  // for(unsigned j=0; j<n_master_nodes; j++)
9615  // {
9616  // hang_info_pt[i]->set_master_node_pt(j,master_node_pt[j],P[i][j]);
9617  // }
9618  // }
9619 
9620  // Set pointers to hanging info
9621  //-----------------------------
9622  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
9623  {
9624  // Check that the node shoule actually hang.
9625  // That is, if the polynomial orders of the elements at a p-type
9626  // non-conormity are both odd then the middle node on the edge is
9627  // shared but a hanging scheme has been constructed for it.
9628  bool node_is_really_shared = false;
9629  for (unsigned m = 0; m < hang_info_pt[i]->nmaster(); m++)
9630  {
9631  // We can simply check if the hanging scheme lists itself as a
9632  // master
9633  if (hang_info_pt[i]->master_node_pt(m) ==
9634  dependent_node_pt[non_vertex_pos[i]])
9635  {
9636  node_is_really_shared = true;
9637 
9638 #ifdef PARANOID
9639  // Also check the corresponding weight: it should be 1
9640  if (std::fabs(hang_info_pt[i]->master_weight(m) - 1.0) >
9641  1.0e-06)
9642  {
9643  throw OomphLibError("Something fishy here -- with shared "
9644  "node at a mortar vertex",
9645  "PRefineableQElemen<2,INITIAL_NNODE_1D>"
9646  "t::quad_hang_helper()",
9647  OOMPH_EXCEPTION_LOCATION);
9648  }
9649 #endif
9650  }
9651  }
9652 
9653  // Now we can make the node hang, if it isn't a shared node
9654  if (!node_is_really_shared)
9655  {
9656  dependent_node_pt[non_vertex_pos[i]]->set_hanging_pt(
9657  hang_info_pt[i], -1);
9658 
9660  // std::cout << "Hanging node: "
9661  // << dependent_node_pt[non_vertex_pos[i]]->x(0) << " "
9662  // << dependent_node_pt[non_vertex_pos[i]]->x(1) << " "
9663  // << dependent_node_pt[non_vertex_pos[i]]->x(2) << " "
9664  // << std::endl;
9665  // for(unsigned m=0; m<hang_info_pt[i]->nmaster(); m++)
9666  // {
9667  // std::cout << " m = " << m << ": "
9668  // << hang_info_pt[i]->master_node_pt(m)->x(0) << " "
9669  // << hang_info_pt[i]->master_node_pt(m)->x(1) << " "
9670  // << hang_info_pt[i]->master_node_pt(m)->x(2) << " "
9671  // << "w = " << hang_info_pt[i]->master_weight(m) <<
9672  // " "
9673  // << std::endl;
9674  // }
9675  }
9676  }
9677 
9678  } // End of case where there are still dependent nodes
9679 
9680  } // End of edge is dependent
9681 
9682  } // End of loop over face edges
9683 
9684  } // End of if face is dependent
9685 
9686  // Now do the mortaring
9687  //---------------------
9688  if (h_type_dependent || p_type_dependent)
9689  {
9690  // Compute the active coordinate indices along the this side of mortar
9691  Vector<unsigned> active_coord_index(2);
9692  if (my_face == B || my_face == F)
9693  {
9694  active_coord_index[0] = 0;
9695  active_coord_index[1] = 1;
9696  }
9697  else if (my_face == D || my_face == U)
9698  {
9699  active_coord_index[0] = 0;
9700  active_coord_index[1] = 2;
9701  }
9702  else if (my_face == L || my_face == R)
9703  {
9704  active_coord_index[0] = 1;
9705  active_coord_index[1] = 2;
9706  }
9707  else
9708  {
9709  throw OomphLibError(
9710  "Cannot transform coordinates",
9711  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9712  OOMPH_EXCEPTION_LOCATION);
9713  }
9714 
9715  // Get pointer to neighbouring master element (in p-refineable form)
9716  PRefineableQElement<3, INITIAL_NNODE_1D>* neigh_obj_pt;
9717  neigh_obj_pt = dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(
9718  neigh_pt->object_pt());
9719 
9720  // Create vector of master and dependent nodes
9721  //----------------------------------------
9722  Vector<Node*> master_node_pt, dependent_node_pt;
9723  Vector<unsigned> master_node_number, dependent_node_number;
9724  Vector<Vector<double>> dependent_node_s_fraction;
9725 
9726  // Number of nodes in one dimension
9727  const unsigned my_n_p = this->ninterpolating_node_1d(value_id);
9728  const unsigned neigh_n_p = neigh_obj_pt->ninterpolating_node_1d(value_id);
9729 
9730  // Storage for pointers to the nodes and their numbers along the master
9731  // edge
9732  unsigned neighbour_node_number = 0;
9733  Node* neighbour_node_pt = 0;
9734 
9735  // Loop over nodes on the face
9736  for (unsigned i0 = 0; i0 < neigh_n_p; i0++)
9737  {
9738  for (unsigned i1 = 0; i1 < neigh_n_p; i1++)
9739  {
9740  const unsigned s0space = 1;
9741  const unsigned s1space = neigh_n_p;
9742  const unsigned s2space = neigh_n_p * neigh_n_p;
9743 
9744  // Find the neighbour's node
9745  switch (neigh_face)
9746  {
9747  case B:
9748  neighbour_node_number = i0 * s0space + i1 * s1space;
9749  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
9750  neighbour_node_number, value_id);
9751  break;
9752 
9753  case F:
9754  neighbour_node_number =
9755  (neigh_n_p - 1) * s2space + i0 * s0space + i1 * s1space;
9756  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
9757  neighbour_node_number, value_id);
9758  break;
9759 
9760  case D:
9761  neighbour_node_number = i0 * s0space + i1 * s2space;
9762  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
9763  neighbour_node_number, value_id);
9764  break;
9765 
9766  case U:
9767  neighbour_node_number =
9768  (neigh_n_p - 1) * s1space + i0 * s0space + i1 * s2space;
9769  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
9770  neighbour_node_number, value_id);
9771  break;
9772 
9773  case L:
9774  neighbour_node_number = i0 * s1space + i1 * s2space;
9775  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
9776  neighbour_node_number, value_id);
9777  break;
9778 
9779  case R:
9780  neighbour_node_number =
9781  (neigh_n_p - 1) * s0space + i0 * s1space + i1 * s2space;
9782  neighbour_node_pt = neigh_obj_pt->interpolating_node_pt(
9783  neighbour_node_number, value_id);
9784  break;
9785 
9786  default:
9787  throw OomphLibError(
9788  "my_face not L, R, D, U, B, F\n",
9789  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9790  OOMPH_EXCEPTION_LOCATION);
9791  }
9792 
9793  // Set node as master node
9794  master_node_number.push_back(neighbour_node_number);
9795  master_node_pt.push_back(neighbour_node_pt);
9796  }
9797  }
9798 
9799  // Storage for pointers to the local nodes and their numbers along my edge
9800  unsigned local_node_number = 0;
9801  Node* local_node_pt = 0;
9802 
9803  // Loop over the nodes along my edge
9804  for (unsigned i0 = 0; i0 < my_n_p; i0++)
9805  {
9806  for (unsigned i1 = 0; i1 < my_n_p; i1++)
9807  {
9808  // Storage for the fractional position of the node in the element
9809  Vector<double> s_fraction(3);
9810 
9811  const unsigned s0space = 1;
9812  const unsigned s1space = my_n_p;
9813  const unsigned s2space = my_n_p * my_n_p;
9814 
9815  // Find the local node and the fractional position of the node
9816  // which depends on the edge, of course
9817  switch (my_face)
9818  {
9819  case B:
9820  s_fraction[0] =
9821  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
9822  s_fraction[1] =
9823  local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
9824  s_fraction[2] = 0.0;
9825  local_node_number = i0 * s0space + i1 * s1space;
9826  local_node_pt =
9827  this->interpolating_node_pt(local_node_number, value_id);
9828  break;
9829 
9830  case F:
9831  s_fraction[0] =
9832  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
9833  s_fraction[1] =
9834  local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
9835  s_fraction[2] = 1.0;
9836  local_node_number =
9837  (my_n_p - 1) * s2space + i0 * s0space + i1 * s1space;
9838  local_node_pt =
9839  this->interpolating_node_pt(local_node_number, value_id);
9840  break;
9841 
9842  case D:
9843  s_fraction[0] =
9844  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
9845  s_fraction[1] = 0.0;
9846  s_fraction[2] =
9847  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
9848  local_node_number = i0 * s0space + i1 * s2space;
9849  local_node_pt =
9850  this->interpolating_node_pt(local_node_number, value_id);
9851  break;
9852 
9853  case U:
9854  s_fraction[0] =
9855  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
9856  s_fraction[1] = 1.0;
9857  s_fraction[2] =
9858  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
9859  local_node_number =
9860  (my_n_p - 1) * s1space + i0 * s0space + i1 * s2space;
9861  local_node_pt =
9862  this->interpolating_node_pt(local_node_number, value_id);
9863  break;
9864 
9865  case L:
9866  s_fraction[0] = 0.0;
9867  s_fraction[1] =
9868  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9869  s_fraction[2] =
9870  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
9871  local_node_number = i0 * s1space + i1 * s2space;
9872  local_node_pt =
9873  this->interpolating_node_pt(local_node_number, value_id);
9874  break;
9875 
9876  case R:
9877  s_fraction[0] = 1.0;
9878  s_fraction[1] =
9879  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
9880  s_fraction[2] =
9881  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
9882  local_node_number =
9883  (my_n_p - 1) * s0space + i0 * s1space + i1 * s2space;
9884  local_node_pt =
9885  this->interpolating_node_pt(local_node_number, value_id);
9886  break;
9887 
9888  default:
9889  throw OomphLibError(
9890  "my_face not L, R, D, U, B, F\n",
9891  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
9892  OOMPH_EXCEPTION_LOCATION);
9893  }
9894 
9895  // Add node to vector of dependent element nodes
9896  dependent_node_number.push_back(local_node_number);
9897  dependent_node_pt.push_back(local_node_pt);
9898 
9899  // Store node's local fraction
9900  dependent_node_s_fraction.push_back(s_fraction);
9901  }
9902  }
9903 
9904  // Store the number of dependent and master nodes for use later
9905  const unsigned n_dependent_nodes = dependent_node_pt.size();
9906  const unsigned n_master_nodes = master_node_pt.size();
9907  const unsigned dependent_element_nnode_1d = my_n_p;
9908  const unsigned master_element_nnode_1d = neigh_n_p;
9909 
9911  // std::cout << "Dependent nodes on face: " <<
9912  // OcTree::Direct_string[my_face] << std::endl; for(unsigned i=0;
9913  // i<dependent_node_pt.size(); i++)
9914  // {
9915  // std::cout << i << ": "
9916  // << dependent_node_pt[i]->x(0) << " "
9917  // << dependent_node_pt[i]->x(1) << " "
9918  // << dependent_node_pt[i]->x(2) << " "
9919  // << std::endl;
9920  // }
9921  // std::cout << "Master nodes on face: " << OcTree::Direct_string[my_face]
9922  // << std::endl; for(unsigned i=0; i<master_node_pt.size(); i++)
9923  // {
9924  // std::cout << i << ": "
9925  // << master_node_pt[i]->x(0) << " "
9926  // << master_node_pt[i]->x(1) << " "
9927  // << master_node_pt[i]->x(2) << " "
9928  // << std::endl;
9929  // }
9930 
9931  // Storage for master shapes
9932  Shape master_shapes(neigh_obj_pt->ninterpolating_node(value_id));
9933 
9934  // Get master and dependent nodal positions
9935  //-------------------------------------
9936  // Get in 1D
9937  Vector<double> dependent_nodal_position_1d;
9938  Vector<double> dependent_weight_1d(dependent_element_nnode_1d);
9939  Orthpoly::gll_nodes(dependent_element_nnode_1d,
9940  dependent_nodal_position_1d,
9941  dependent_weight_1d);
9942  Vector<double> master_nodal_position_1d;
9943  Vector<double> master_weight_1d(master_element_nnode_1d);
9944  Orthpoly::gll_nodes(
9945  master_element_nnode_1d, master_nodal_position_1d, master_weight_1d);
9946 
9947  // Storage for 2D
9948  Vector<Vector<double>> dependent_nodal_position(
9949  dependent_element_nnode_1d * dependent_element_nnode_1d);
9950  for (unsigned i = 0; i < dependent_nodal_position.size(); i++)
9951  {
9952  dependent_nodal_position[i].resize(2);
9953  }
9954  Vector<double> dependent_weight(dependent_element_nnode_1d *
9955  dependent_element_nnode_1d);
9956  Vector<Vector<double>> master_nodal_position(master_element_nnode_1d *
9957  master_element_nnode_1d);
9958  for (unsigned i = 0; i < master_nodal_position.size(); i++)
9959  {
9960  master_nodal_position[i].resize(2);
9961  }
9962  Vector<double> master_weight(master_element_nnode_1d *
9963  master_element_nnode_1d);
9964 
9965  // Fill in coordinates and weights in 2D
9966  unsigned dependent_index = 0;
9967  for (unsigned i = 0; i < dependent_element_nnode_1d; i++)
9968  {
9969  for (unsigned j = 0; j < dependent_element_nnode_1d; j++)
9970  {
9971  dependent_nodal_position[dependent_index][0] =
9972  dependent_nodal_position_1d[i];
9973  dependent_nodal_position[dependent_index][1] =
9974  dependent_nodal_position_1d[j];
9975  dependent_weight[dependent_index] =
9976  dependent_weight_1d[i] * dependent_weight_1d[j];
9977  dependent_index++;
9978  }
9979  }
9980  unsigned master_index = 0;
9981  for (unsigned i = 0; i < master_element_nnode_1d; i++)
9982  {
9983  for (unsigned j = 0; j < master_element_nnode_1d; j++)
9984  {
9985  master_nodal_position[master_index][0] = master_nodal_position_1d[i];
9986  master_nodal_position[master_index][1] = master_nodal_position_1d[j];
9987  master_weight[master_index] =
9988  master_weight_1d[i] * master_weight_1d[j];
9989  master_index++;
9990  }
9991  }
9992 
9993  // Apply the vertex matching condition
9994  //------------------------------------
9995  // Vertiex matching is ensured automatically in cases where there is a
9996  // node at each end of the mortar that is shared between the master and
9997  // dependent elements. Where this is not the case, the vertex matching
9998  // condition must be enforced by constraining the value of the unknown at
9999  // the node on the dependent side to be the same as the value at that
10000  // location in the master.
10001 
10002  // Store positions of the mortar vertex/non-vertex nodes in the dependent
10003  // element
10004  // const unsigned n_mortar_vertices = 4;
10005  // Vector<unsigned> vertex_pos(n_mortar_vertices);
10006  // vertex_pos[0] = 0;
10007  // vertex_pos[1] = my_n_p-1;
10008  // vertex_pos[2] = my_n_p*(my_n_p-1);
10009  // vertex_pos[3] = my_n_p*my_n_p-1;
10010  Vector<unsigned> non_vertex_pos((dependent_element_nnode_1d - 2) *
10011  (dependent_element_nnode_1d - 2));
10012  unsigned nvindex = 0;
10013  for (unsigned i = 1; i < dependent_element_nnode_1d - 1; i++)
10014  {
10015  for (unsigned j = 1; j < dependent_element_nnode_1d - 1; j++)
10016  {
10017  non_vertex_pos[nvindex++] = i * dependent_element_nnode_1d + j;
10018  }
10019  }
10020  Vector<unsigned> vertex_pos;
10021  for (unsigned i = 0; i < n_dependent_nodes; i++)
10022  {
10023  // Check if node number is in the non-vertex storage
10024  bool node_is_vertex = true;
10025  for (unsigned j = 0; j < non_vertex_pos.size(); j++)
10026  {
10027  if (i == non_vertex_pos[j])
10028  {
10029  // Node is not a vertex
10030  node_is_vertex = false;
10031  break;
10032  }
10033  }
10034  // If we get here and the node is a vertex then add it's index
10035  if (node_is_vertex)
10036  {
10037  vertex_pos.push_back(i);
10038  }
10039  }
10040  // Store number of mortar vertices
10041  const unsigned n_mortar_vertices = vertex_pos.size();
10042 
10043  // Check that there are actually dependent nodes for which we still need
10044  // to construct a hanging scheme. If not then there is nothing more to do.
10045  if (n_dependent_nodes - n_mortar_vertices > 0)
10046  {
10047  // Assemble mass matrix for mortar
10048  //--------------------------------
10049  Vector<double> psi(n_dependent_nodes - n_mortar_vertices);
10050  Vector<double> diag_M(n_dependent_nodes - n_mortar_vertices);
10051  Vector<Vector<double>> shared_node_M(n_mortar_vertices);
10052  for (unsigned i = 0; i < shared_node_M.size(); i++)
10053  {
10054  shared_node_M[i].resize(n_dependent_nodes - n_mortar_vertices);
10055  }
10056 
10057  // Fill in part corresponding to dependent nodal positions (unknown)
10058  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
10059  {
10060  // Mortar test functions in 2D are just the cross product of the 1D
10061  // test functions Initialise to 1
10062  psi[i] = 1.0;
10063  // Take product in each direction
10064  for (unsigned dir = 0; dir < 2; dir++)
10065  {
10066  unsigned index1d =
10067  (dir == 0) ? i : i % (dependent_element_nnode_1d - 2);
10068  // Use L'Hosptal's rule:
10069  psi[i] *=
10070  pow(-1.0, int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10071  -Orthpoly::ddlegendre(
10072  dependent_element_nnode_1d - 1,
10073  dependent_nodal_position[non_vertex_pos[i]][dir]);
10074  }
10075  // Put in contribution
10076  diag_M[i] = psi[i] * dependent_weight[non_vertex_pos[i]];
10077  }
10078 
10080  // std::cout << "diag(M):" << std::endl;
10081  // for(unsigned i=0; i<diag_M.size(); i++)
10082  // {
10083  // std::cout << " " << diag_M[i];
10084  // }
10085  // std::cout << std::endl;
10086 
10087  // Fill in part corresponding to dependent element vertices (known)
10088  for (unsigned v = 0; v < shared_node_M.size(); v++)
10089  {
10090  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
10091  {
10092  // Mortar test functions in 2D are just the cross product of the 1D
10093  // test functions Initialise to 1
10094  psi[i] = 1.0;
10095  // Take product in each direction
10096  for (unsigned dir = 0; dir < 2; dir++)
10097  {
10098  unsigned index1d =
10099  (dir == 0) ? i : i % (dependent_element_nnode_1d - 2);
10100  // Check if denominator is zero
10101  if (std::fabs(dependent_nodal_position[non_vertex_pos[i]][dir] -
10102  dependent_nodal_position[vertex_pos[v]][dir]) >=
10103  1.0e-8)
10104  {
10105  // We're ok
10106  psi[i] *=
10107  pow(-1.0,
10108  int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10109  Orthpoly::dlegendre(
10110  dependent_element_nnode_1d - 1,
10111  dependent_nodal_position[vertex_pos[v]][dir]) /
10112  (dependent_nodal_position[non_vertex_pos[i]][dir] -
10113  dependent_nodal_position[vertex_pos[v]][dir]);
10114  }
10115  // Check if numerator is zero
10116  else if (std::fabs(Orthpoly::dlegendre(
10117  dependent_element_nnode_1d - 1,
10118  dependent_nodal_position[vertex_pos[v]][dir])) <
10119  1.0e-8)
10120  {
10121  // We can use l'hopital's rule
10122  psi[i] *=
10123  pow(-1.0,
10124  int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10125  -Orthpoly::ddlegendre(
10126  dependent_element_nnode_1d - 1,
10127  dependent_nodal_position[non_vertex_pos[i]][dir]);
10128  }
10129  else
10130  {
10131  // We can't use l'hopital's rule
10132  throw OomphLibError(
10133  "Cannot use l'Hopital's rule. Dividing by zero is not "
10134  "allowed!",
10135  "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
10136  OOMPH_EXCEPTION_LOCATION);
10137  }
10138  }
10139  // Put in contribution
10140  shared_node_M[v][i] = psi[i] * dependent_weight[vertex_pos[v]];
10141  }
10142  }
10143 
10145  // std::cout << "shared node M:" << std::endl;
10146  // for(unsigned i=0; i<shared_node_M.size(); i++)
10147  // {
10148  // for(unsigned j=0; j<shared_node_M[i].size(); j++)
10149  // {
10150  // std::cout << " " << shared_node_M[i][j];
10151  // }
10152  // }
10153  // std::cout << std::endl;
10154 
10155  // Assemble local projection matrix for mortar
10156  //--------------------------------------------
10157 
10158  // Have only one local projection matrix because there is just one
10159  // master
10160  Vector<Vector<double>> P(n_dependent_nodes - n_mortar_vertices);
10161  for (unsigned i = 0; i < P.size(); i++)
10162  {
10163  P[i].resize(n_master_nodes, 0.0);
10164  }
10165 
10166  // Storage for local coordinate
10167  Vector<double> s(3);
10168 
10169  // Sum contributions from master element shapes (quadrature).
10170  // The order of quadrature must be high enough to evaluate a polynomial
10171  // of order N_s + N_m - 3 exactly, where N_s = n_dependent_nodes, N_m =
10172  // n_master_nodes.
10173  // (Use pointers for the quadrature knots and weights so that
10174  // data is not unnecessarily copied)
10175  // unsigned quadrature_order =
10176  // std::max(dependent_element_nnode_1d,master_element_nnode_1d);
10177  Vector<Vector<double>>* quadrature_knot;
10178  Vector<double>* quadrature_weight;
10179  if (dependent_element_nnode_1d >= master_element_nnode_1d)
10180  {
10181  // Use the same quadrature order as the dependent element (me)
10182  quadrature_knot = &dependent_nodal_position;
10183  quadrature_weight = &dependent_weight;
10184  }
10185  else
10186  {
10187  // Use the same quadrature order as the master element (neighbour)
10188  quadrature_knot = &master_nodal_position;
10189  quadrature_weight = &master_weight;
10190  }
10191 
10192  // Quadrature loop
10193  for (unsigned q = 0; q < (*quadrature_weight).size(); q++)
10194  {
10195  // Evaluate mortar test functions at the quadrature knots in the
10196  // dependent
10197  Vector<double> s_on_mortar(2);
10198  for (unsigned i = 0; i < 2; i++)
10199  {
10200  s_on_mortar[i] = (*quadrature_knot)[q][i];
10201  }
10202 
10203  // Get psi
10204  for (unsigned k = 0; k < n_dependent_nodes - n_mortar_vertices; k++)
10205  {
10206  // Mortar test functions in 2D are just the cross product of the 1D
10207  // test functions Initialise to 1
10208  psi[k] = 1.0;
10209  // Take product in each direction
10210  for (unsigned dir = 0; dir < 2; dir++)
10211  {
10212  unsigned index1d =
10213  (dir == 0) ? k : k % (dependent_element_nnode_1d - 2);
10214  // Check if denominator is zero
10215  if (std::fabs(dependent_nodal_position[non_vertex_pos[k]][dir] -
10216  s_on_mortar[dir]) >= 1.0e-08)
10217  {
10218  // We're ok
10219  psi[k] *=
10220  pow(-1.0,
10221  int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10222  Orthpoly::dlegendre(dependent_element_nnode_1d - 1,
10223  s_on_mortar[dir]) /
10224  (dependent_nodal_position[non_vertex_pos[k]][dir] -
10225  s_on_mortar[dir]);
10226  }
10227  // Check if numerator is zero
10228  else if (std::fabs(Orthpoly::dlegendre(
10229  dependent_element_nnode_1d - 1, s_on_mortar[dir])) <
10230  1.0e-8)
10231  {
10232  // We can use l'Hopital's rule
10233  psi[k] *=
10234  pow(-1.0,
10235  int((dependent_element_nnode_1d - 1) - index1d - 1)) *
10236  -Orthpoly::ddlegendre(dependent_element_nnode_1d - 1,
10237  s_on_mortar[dir]);
10238  }
10239  else
10240  {
10241  // We can't use l'hopital's rule
10242  throw OomphLibError(
10243  "Cannot use l'Hopital's rule. Dividing by zero is not "
10244  "allowed!",
10245  "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
10246  OOMPH_EXCEPTION_LOCATION);
10247  }
10248  }
10249  }
10250 
10251  // Convert coordinate on mortar to local fraction in dependent element
10252  Vector<double> s_fraction(3);
10253  for (unsigned i = 0; i < 3; i++)
10254  {
10255  if (i == active_coord_index[0])
10256  {
10257  s_fraction[i] = 0.5 * (s_on_mortar[0] + 1.0);
10258  }
10259  else if (i == active_coord_index[1])
10260  {
10261  s_fraction[i] = 0.5 * (s_on_mortar[1] + 1.0);
10262  }
10263  else
10264  {
10265  s_fraction[i] = dependent_node_s_fraction[vertex_pos[0]][i];
10266  }
10267  }
10268 
10269  // Project active coordinate into master element
10270  Vector<double> s_in_neigh(3);
10271  for (unsigned i = 0; i < 3; i++)
10272  {
10273  s_in_neigh[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
10274  (s_hi_neigh[i] - s_lo_neigh[i]);
10275  }
10276 
10277  // Evaluate master shapes at projections of local quadrature knots
10278  neigh_obj_pt->interpolating_basis(
10279  s_in_neigh, master_shapes, value_id);
10280 
10281  // Populate local projection matrix
10282  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
10283  {
10284  for (unsigned j = 0; j < n_master_nodes; j++)
10285  {
10286  P[i][j] += master_shapes[master_node_number[j]] * psi[i] *
10287  (*quadrature_weight)[q];
10288  }
10289  }
10290  }
10291 
10293  // std::cout << "P_e:" << std::endl;
10294  // for(unsigned i=0; i<P.size(); i++)
10295  // {
10296  // for(unsigned j=0; j<P[i].size(); j++)
10297  // {
10298  // std::cout << " " << P[i][j];
10299  // }
10300  // }
10301  // std::cout << std::endl;
10302 
10303  // Assemble global projection matrices for mortar
10304  //-----------------------------------------------
10305  // Need to subtract contributions from the "known unknowns"
10306  // corresponding to the nodes at the vertices of the mortar
10307 
10308  // Assemble contributions from mortar vertex nodes
10309  for (unsigned v = 0; v < n_mortar_vertices; v++)
10310  {
10311  // Convert coordinate on mortar to local fraction in dependent element
10312  Vector<double> s_fraction(3);
10313  for (unsigned i = 0; i < 3; i++)
10314  {
10315  if (i == active_coord_index[0])
10316  {
10317  s_fraction[i] =
10318  0.5 * (dependent_nodal_position[vertex_pos[v]][0] + 1.0);
10319  }
10320  else if (i == active_coord_index[1])
10321  {
10322  s_fraction[i] =
10323  0.5 * (dependent_nodal_position[vertex_pos[v]][1] + 1.0);
10324  }
10325  else
10326  {
10327  s_fraction[i] = dependent_node_s_fraction[vertex_pos[0]][i];
10328  }
10329  }
10330 
10331  // Project active coordinate into master element
10332  Vector<double> s_in_neigh(3);
10333  for (unsigned i = 0; i < 3; i++)
10334  {
10335  s_in_neigh[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
10336  (s_hi_neigh[i] - s_lo_neigh[i]);
10337  }
10338 
10339  // Get master shapes at dependent nodal positions
10340  neigh_obj_pt->interpolating_basis(
10341  s_in_neigh, master_shapes, value_id);
10342 
10343  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
10344  {
10345  for (unsigned k = 0; k < n_master_nodes; k++)
10346  {
10347  P[i][k] -=
10348  master_shapes[master_node_number[k]] * shared_node_M[v][i];
10349  }
10350  }
10351  }
10352 
10354  // std::cout << "P:" << std::endl;
10355  // for(unsigned i=0; i<P.size(); i++)
10356  // {
10357  // for(unsigned j=0; j<P[i].size(); j++)
10358  // {
10359  // std::cout << " " << P[i][j];
10360  // }
10361  // }
10362  // std::cout << std::endl;
10363 
10364  // Solve mortar system
10365  //--------------------
10366  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
10367  {
10368  for (unsigned j = 0; j < n_master_nodes; j++)
10369  {
10370  P[i][j] /= diag_M[i];
10371  }
10372  }
10373 
10375  // std::cout << "solved P:" << std::endl;
10376  // for(unsigned i=0; i<P.size(); i++)
10377  // {
10378  // for(unsigned j=0; j<P[i].size(); j++)
10379  // {
10380  // std::cout << " " << P[i][j];
10381  // }
10382  // }
10383  // std::cout << std::endl;
10384 
10385  // Create structures to hold the hanging info
10386  //-------------------------------------------
10387  Vector<HangInfo*> hang_info_pt(n_dependent_nodes);
10388  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
10389  {
10390  hang_info_pt[i] = new HangInfo(n_master_nodes);
10391  }
10392 
10393  // Copy information to hanging nodes
10394  //----------------------------------
10395  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
10396  {
10397  for (unsigned j = 0; j < n_master_nodes; j++)
10398  {
10399  hang_info_pt[i]->set_master_node_pt(j, master_node_pt[j], P[i][j]);
10400  }
10401  }
10402 
10403  // Set pointers to hanging info
10404  //-----------------------------
10405  for (unsigned i = 0; i < n_dependent_nodes - n_mortar_vertices; i++)
10406  {
10407  // Check that the node shoule actually hang.
10408  // That is, if the polynomial orders of the elements at a p-type
10409  // non-conormity are both odd then the middle node on the edge is
10410  // shared but a hanging scheme has been constructed for it.
10411  bool node_is_really_shared = false;
10412  for (unsigned m = 0; m < hang_info_pt[i]->nmaster(); m++)
10413  {
10414  // We can simply check if the hanging scheme lists itself as a
10415  // master
10416  if (hang_info_pt[i]->master_node_pt(m) ==
10417  dependent_node_pt[non_vertex_pos[i]])
10418  {
10419  node_is_really_shared = true;
10420 
10421 #ifdef PARANOID
10422  // Also check the corresponding weight: it should be 1
10423  if (std::fabs(hang_info_pt[i]->master_weight(m) - 1.0) > 1.0e-06)
10424  {
10425  throw OomphLibError(
10426  "Something fishy here -- with shared node at a mortar vertex",
10427  "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
10428  OOMPH_EXCEPTION_LOCATION);
10429  }
10430 #endif
10431  }
10432  }
10433 
10434  // Now we can make the node hang, if it isn't a shared node
10435  if (!node_is_really_shared)
10436  {
10437  dependent_node_pt[non_vertex_pos[i]]->set_hanging_pt(
10438  hang_info_pt[i], -1);
10439  }
10440  }
10441 
10442  } // End of case where there are still dependent nodes
10443 
10444 #ifdef PARANOID
10445  // Check all dependent nodes, hanging or otherwise
10446  for (unsigned i = 0; i < n_dependent_nodes; i++)
10447  {
10448  // Check that weights sum to 1 for those that hang
10449  if (dependent_node_pt[i]->is_hanging())
10450  {
10451  // Check that weights sum to 1 and no zero weights
10452  double weight_sum = 0.0;
10453  bool zero_weight = false;
10454  for (unsigned m = 0;
10455  m < dependent_node_pt[i]->hanging_pt()->nmaster();
10456  m++)
10457  {
10458  weight_sum += dependent_node_pt[i]->hanging_pt()->master_weight(m);
10459  if (std::fabs(dependent_node_pt[i]->hanging_pt()->master_weight(
10460  m)) < 1.0e-14)
10461  {
10462  zero_weight = true;
10463  oomph_info << "In the hanging scheme for dependent node " << i
10464  << ", master node " << m << " has weight "
10465  << dependent_node_pt[i]->hanging_pt()->master_weight(m)
10466  << " < 1.0e-14" << std::endl;
10467  }
10468  }
10469 
10470  // Warn if not
10471  if (std::fabs(weight_sum - 1.0) > 1.0e-08)
10472  {
10473  oomph_info << "Sum of master weights: " << weight_sum << std::endl;
10474  OomphLibWarning(
10475  "Weights in hanging scheme do not sum to 1",
10476  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10477  OOMPH_EXCEPTION_LOCATION);
10478  }
10479  if (zero_weight)
10480  {
10481  OomphLibWarning(
10482  "Zero weights present in hanging schemes",
10483  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10484  OOMPH_EXCEPTION_LOCATION);
10485  }
10486 
10487  // Also check that dependent nodes do not have themselves as masters
10488  // bool dependent_should_not_be_hanging = false;
10489  for (unsigned m = 0;
10490  m < dependent_node_pt[i]->hanging_pt()->nmaster();
10491  m++)
10492  {
10493  if (dependent_node_pt[i] ==
10494  dependent_node_pt[i]->hanging_pt()->master_node_pt(m))
10495  {
10496  // This shouldn't happen!
10497  throw OomphLibError(
10498  "Dependent node has itself as a master!",
10499  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10500  OOMPH_EXCEPTION_LOCATION);
10501  }
10502  if (dependent_node_pt[i]
10503  ->hanging_pt()
10504  ->master_node_pt(m)
10505  ->is_hanging())
10506  {
10507  // Check if node is master of master
10508  Node* new_nod_pt =
10509  dependent_node_pt[i]->hanging_pt()->master_node_pt(m);
10510  for (unsigned mm = 0; mm < new_nod_pt->hanging_pt()->nmaster();
10511  mm++)
10512  {
10513  if (dependent_node_pt[i] ==
10514  new_nod_pt->hanging_pt()->master_node_pt(mm))
10515  {
10516  std::cout << "++++++++++++++++++++++++++++++++++++++++"
10517  << std::endl;
10518  std::cout
10519  << " Dependent node: " << dependent_node_pt[i]
10520  << " = " << dependent_node_pt[i]->x(0) << " "
10521  << dependent_node_pt[i]->x(1) << " "
10522  << dependent_node_pt[i]->x(2) << " " << std::endl;
10523  std::cout
10524  << " Master node: "
10525  << dependent_node_pt[i]->hanging_pt()->master_node_pt(m)
10526  << " = "
10527  << dependent_node_pt[i]->hanging_pt()->master_node_pt(m)->x(
10528  0)
10529  << " "
10530  << dependent_node_pt[i]->hanging_pt()->master_node_pt(m)->x(
10531  1)
10532  << " "
10533  << dependent_node_pt[i]->hanging_pt()->master_node_pt(m)->x(
10534  2)
10535  << " " << std::endl;
10536  std::cout << "Master's master node: "
10537  << dependent_node_pt[i]
10538  ->hanging_pt()
10539  ->master_node_pt(m)
10540  ->hanging_pt()
10541  ->master_node_pt(mm)
10542  << " = "
10543  << dependent_node_pt[i]
10544  ->hanging_pt()
10545  ->master_node_pt(m)
10546  ->hanging_pt()
10547  ->master_node_pt(mm)
10548  ->x(0)
10549  << " "
10550  << dependent_node_pt[i]
10551  ->hanging_pt()
10552  ->master_node_pt(m)
10553  ->hanging_pt()
10554  ->master_node_pt(mm)
10555  ->x(1)
10556  << " "
10557  << dependent_node_pt[i]
10558  ->hanging_pt()
10559  ->master_node_pt(m)
10560  ->hanging_pt()
10561  ->master_node_pt(mm)
10562  ->x(2)
10563  << " " << std::endl;
10564 
10565 
10566  // Print out hanging scheme
10567  std::cout << "Hanging node: " << dependent_node_pt[i]->x(0)
10568  << " " << dependent_node_pt[i]->x(1) << " "
10569  << dependent_node_pt[i]->x(2) << " " << std::endl;
10570  for (unsigned m_tmp = 0;
10571  m_tmp < dependent_node_pt[i]->hanging_pt()->nmaster();
10572  m_tmp++)
10573  {
10574  std::cout
10575  << " m = " << m_tmp << ": "
10576  << dependent_node_pt[i]
10577  ->hanging_pt()
10578  ->master_node_pt(m_tmp)
10579  ->x(0)
10580  << " "
10581  << dependent_node_pt[i]
10582  ->hanging_pt()
10583  ->master_node_pt(m_tmp)
10584  ->x(1)
10585  << " "
10586  << dependent_node_pt[i]
10587  ->hanging_pt()
10588  ->master_node_pt(m_tmp)
10589  ->x(2)
10590  << " "
10591  << "w = "
10592  << dependent_node_pt[i]->hanging_pt()->master_weight(
10593  m_tmp)
10594  << " " << std::endl;
10595  }
10596 
10597  // Print out hanging scheme
10598  std::cout << "Master node " << m
10599  << " of Hanging node: " << new_nod_pt->x(0) << " "
10600  << new_nod_pt->x(1) << " " << new_nod_pt->x(2)
10601  << " " << std::endl;
10602  for (unsigned mm_tmp = 0;
10603  mm_tmp < new_nod_pt->hanging_pt()->nmaster();
10604  mm_tmp++)
10605  {
10606  std::cout
10607  << " mm = " << mm_tmp << ": "
10608  << new_nod_pt->hanging_pt()->master_node_pt(mm_tmp)->x(0)
10609  << " "
10610  << new_nod_pt->hanging_pt()->master_node_pt(mm_tmp)->x(1)
10611  << " "
10612  << new_nod_pt->hanging_pt()->master_node_pt(mm_tmp)->x(2)
10613  << " "
10614  << "w = "
10615  << new_nod_pt->hanging_pt()->master_weight(mm_tmp) << " "
10616  << std::endl;
10617  }
10618 
10619  // This shouldn't happen!
10620  throw OomphLibError(
10621  "Dependent node has itself as a master of a master!",
10622  "PRefineableQElement<3,INITIAL_NNODE_1D>::oc_hang_helper()",
10623  OOMPH_EXCEPTION_LOCATION);
10624  }
10625  }
10626  }
10627  }
10628  }
10629  else
10630  {
10631  // Check that this node is shared with the master element if it
10632  // isn't hanging
10633  bool is_master = false;
10634  for (unsigned n = 0; n < n_master_nodes; n++)
10635  {
10636  if (dependent_node_pt[i] == master_node_pt[n])
10637  {
10638  // Node is a master
10639  is_master = true;
10640  break;
10641  }
10642  }
10643 
10644  if (!is_master)
10645  {
10647  // std::ostringstream error_string;
10648  // error_string
10649  // << "This node in the dependent element is neither" << std::endl
10650  // << "hanging or shared with a master element." << std::endl;
10651  //
10652  // throw OomphLibError(
10653  // error_string.str(),
10654  // "PRefineableQElement<3,INITIAL_NNODE_1D>::quad_hang_helper()",
10655  // OOMPH_EXCEPTION_LOCATION);
10656  }
10657  }
10658  }
10659 #endif
10660 
10661  // Finally, Loop over all dependent nodes and fine-tune their positions
10662  //-----------------------------------------------------------------
10663  // Here we simply set the node's positions to be consistent
10664  // with the hanging scheme. This is not strictly necessary
10665  // because it is done in the mesh adaptation before the node
10666  // becomes non-hanging later on. We make no attempt to ensure
10667  // (strong) continuity in the position across the mortar.
10668  for (unsigned i = 0; i < n_dependent_nodes; i++)
10669  {
10670  // Only fine-tune hanging nodes
10671  if (dependent_node_pt[i]->is_hanging())
10672  {
10673  // If we are doing the position, then
10674  if (value_id == -1)
10675  {
10676  // Get the local coordinate of this dependent node
10677  Vector<double> s_local(3);
10678  this->local_coordinate_of_node(dependent_node_number[i], s_local);
10679 
10680  // Get the position from interpolation in this element via
10681  // the hanging scheme
10682  Vector<double> x_in_neighb(3);
10683  this->interpolated_x(s_local, x_in_neighb);
10684 
10685  // Fine adjust the coordinates (macro map will pick up boundary
10686  // accurately but will lead to different element edges)
10687  dependent_node_pt[i]->x(0) = x_in_neighb[0];
10688  dependent_node_pt[i]->x(1) = x_in_neighb[1];
10689  dependent_node_pt[i]->x(2) = x_in_neighb[2];
10690  }
10691  }
10692  }
10693  } // End of case where this interface is to be mortared
10694  }
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
B
Matrix< SCALARB, Dynamic, Dynamic, opt_B > B
Definition: bench_gemm.cpp:48
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
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

References D, oomph::OcTreeNames::DB, oomph::Orthpoly::ddlegendre(), oomph::OcTreeNames::DF, oomph::Orthpoly::dlegendre(), oomph::OcTreeNames::F, boost::multiprecision::fabs(), oomph::Orthpoly::gll_nodes(), oomph::OcTree::gteq_face_neighbour(), oomph::OcTree::gteq_true_edge_neighbour(), oomph::Node::hanging_pt(), i, oomph::RefineableElement::interpolating_basis(), oomph::RefineableElement::interpolating_node_pt(), j, k, L, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, m, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), n, oomph::RefineableElement::ninterpolating_node(), oomph::RefineableElement::ninterpolating_node_1d(), oomph::HangInfo::nmaster(), oomph::Tree::nsons(), oomph::Tree::object_pt(), OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Global_Physical_Variables::P, Eigen::bfloat16_impl::pow(), Eigen::numext::q, R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, oomph::OcTreeNames::RF, oomph::OcTreeNames::RU, s, oomph::HangInfo::set_master_node_pt(), RachelsAdvectionDiffusion::U, oomph::OcTreeNames::UB, oomph::OcTreeNames::UF, v, and oomph::Node::x().

◆ p_refine()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 3, 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.)

Implements oomph::PRefineableElement.

6106  {
6107  // Number of dimensions
6108  unsigned n_dim = 3;
6109 
6110  // Cast clone to correct type
6111  PRefineableQElement<3, INITIAL_NNODE_1D>* clone_el_pt =
6112  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(clone_pt);
6113 
6114  // Check if we can use it
6115  if (clone_el_pt == 0)
6116  {
6117  throw OomphLibError(
6118  "Cloned copy must be a PRefineableQElement<3,INITIAL_NNODE_1D>!",
6119  OOMPH_CURRENT_FUNCTION,
6120  OOMPH_EXCEPTION_LOCATION);
6121  }
6122 #ifdef PARANOID
6123  // Clone exists, so check that it is infact a copy of me
6124  else
6125  {
6126  // Flag to keep track of check
6127  bool clone_is_ok = true;
6128 
6129  // Does the clone have the correct p-order?
6130  clone_is_ok = clone_is_ok && (clone_el_pt->p_order() == this->p_order());
6131 
6132  if (!clone_is_ok)
6133  {
6134  std::ostringstream err_stream;
6135  err_stream << "Clone element has a different p-order from me,"
6136  << std::endl
6137  << "but it is supposed to be a copy!" << std::endl;
6138 
6139  throw OomphLibError(
6140  err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6141  }
6142 
6143  // Does the clone have the same number of nodes as me?
6144  clone_is_ok = clone_is_ok && (clone_el_pt->nnode() == this->nnode());
6145 
6146  if (!clone_is_ok)
6147  {
6148  std::ostringstream err_stream;
6149  err_stream << "Clone element has a different number of nodes from me,"
6150  << std::endl
6151  << "but it is supposed to be a copy!" << std::endl;
6152 
6153  throw OomphLibError(
6154  err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6155  }
6156 
6157  // Does the clone have the same nodes as me?
6158  for (unsigned n = 0; n < this->nnode(); n++)
6159  {
6160  clone_is_ok =
6161  clone_is_ok && (clone_el_pt->node_pt(n) == this->node_pt(n));
6162  }
6163 
6164  if (!clone_is_ok)
6165  {
6166  std::ostringstream err_stream;
6167  err_stream << "The nodes belonging to the clone element are different"
6168  << std::endl
6169  << "from mine, but it is supposed to be a copy!"
6170  << std::endl;
6171 
6172  throw OomphLibError(
6173  err_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6174  }
6175 
6176  // If we get to here then the clone has all the information we require
6177  }
6178 #endif
6179 
6180  // Currently we can't handle the case of generalised coordinates
6181  // since we haven't established how they should be interpolated.
6182  // Buffer this case:
6183  if (clone_el_pt->node_pt(0)->nposition_type() != 1)
6184  {
6185  throw OomphLibError("Can't handle generalised nodal positions (yet).",
6186  OOMPH_CURRENT_FUNCTION,
6187  OOMPH_EXCEPTION_LOCATION);
6188  }
6189 
6190  // Timestepper should be the same for all nodes -- use it
6191  // to create timesteppers for new nodes
6192  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
6193 
6194  // Get number of history values (incl. present)
6195  unsigned ntstorage = time_stepper_pt->ntstorage();
6196 
6197  // Increment p-order of the element
6198  this->P_order += inc;
6199 
6200  // Change integration scheme
6201  delete this->integral_pt();
6202  switch (this->P_order)
6203  {
6204  case 2:
6205  this->set_integration_scheme(new GaussLobattoLegendre<3, 2>);
6206  break;
6207  case 3:
6208  this->set_integration_scheme(new GaussLobattoLegendre<3, 3>);
6209  break;
6210  case 4:
6211  this->set_integration_scheme(new GaussLobattoLegendre<3, 4>);
6212  break;
6213  case 5:
6214  this->set_integration_scheme(new GaussLobattoLegendre<3, 5>);
6215  break;
6216  case 6:
6217  this->set_integration_scheme(new GaussLobattoLegendre<3, 6>);
6218  break;
6219  case 7:
6220  this->set_integration_scheme(new GaussLobattoLegendre<3, 7>);
6221  break;
6222  default:
6223  std::ostringstream error_message;
6224  error_message << "\nERROR: Exceeded maximum polynomial order for";
6225  error_message << "\n integration scheme.\n";
6226  throw OomphLibError(error_message.str(),
6227  OOMPH_CURRENT_FUNCTION,
6228  OOMPH_EXCEPTION_LOCATION);
6229  }
6230 
6231  // Allocate new space for Nodes (at the element level)
6232  this->set_n_node(this->P_order * this->P_order * this->P_order);
6233 
6234  // Copy vertex nodes and create new edge and internal nodes
6235  //---------------------------------------------------------
6236 
6237  // Setup vertex coordinates in element:
6238  //-------------------------------------
6239  Vector<double> s_lo(n_dim, 0.0);
6240  Vector<double> s_hi(n_dim, 0.0);
6241  s_lo[0] = -1.0;
6242  s_hi[0] = 1.0;
6243  s_lo[1] = -1.0;
6244  s_hi[1] = 1.0;
6245  s_lo[2] = -1.0;
6246  s_hi[2] = 1.0;
6247 
6248  // Local coordinate in element
6249  Vector<double> s(n_dim, 0.0);
6250 
6251  unsigned jnod = 0;
6252 
6253  Vector<double> s_fraction(n_dim, 0.0);
6254  // Loop over nodes in element
6255  for (unsigned i0 = 0; i0 < this->P_order; i0++)
6256  {
6257  // Get the fractional position of the node in the direction of s[0]
6258  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
6259  // Local coordinate
6260  s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
6261 
6262  for (unsigned i1 = 0; i1 < this->P_order; i1++)
6263  {
6264  // Get the fractional position of the node in the direction of s[1]
6265  s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
6266  // Local coordinate
6267  s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
6268 
6269  for (unsigned i2 = 0; i2 < this->P_order; i2++)
6270  {
6271  // Get the fractional position of the node in the direction of s[2]
6272  s_fraction[2] = this->local_one_d_fraction_of_node(i2, 2);
6273  // Local coordinate
6274  s[2] = s_lo[2] + (s_hi[2] - s_lo[2]) * s_fraction[2];
6275 
6276  // Local node number
6277  jnod = i0 + this->P_order * i1 + this->P_order * this->P_order * i2;
6278 
6279  // Initialise flag: So far, this node hasn't been built
6280  // or copied yet
6281  bool node_done = false;
6282 
6283  // Get the pointer to the node in this element (or rather, its clone),
6284  // returns NULL if there is not node
6285  Node* created_node_pt = clone_el_pt->get_node_at_local_coordinate(s);
6286  // created_node_pt = 0;
6287 
6288  // Does this node already exist in this element?
6289  //----------------------------------------------
6290  if (created_node_pt != 0)
6291  {
6292  // Copy node across
6293  this->node_pt(jnod) = created_node_pt;
6294 
6295  // Make sure that we update the values at the node so that
6296  // they are consistent with the present representation.
6297  // This is only need for mixed interpolation where the value
6298  // at the father could now become active.
6299 
6300  // Loop over all history values
6301  for (unsigned t = 0; t < ntstorage; t++)
6302  {
6303  // Get values from father element
6304  // Note: get_interpolated_values() sets Vector size itself.
6305  Vector<double> prev_values;
6306  clone_el_pt->get_interpolated_values(t, s, prev_values);
6307  // Find the minimum number of values
6308  //(either those stored at the node, or those returned by
6309  // the function)
6310  unsigned n_val_at_node = created_node_pt->nvalue();
6311  unsigned n_val_from_function = prev_values.size();
6312  // Use the ternary conditional operator here
6313  unsigned n_var = n_val_at_node < n_val_from_function ?
6314  n_val_at_node :
6315  n_val_from_function;
6316  // Assign the values that we can
6317  for (unsigned k = 0; k < n_var; k++)
6318  {
6319  created_node_pt->set_value(t, k, prev_values[k]);
6320  }
6321  }
6322 
6323  // Node has been created by copy
6324  node_done = true;
6325  }
6326  // Node does not exist in this element but might already
6327  //------------------------------------------------------
6328  // have been created by neighbouring elements
6329  //-------------------------------------------
6330  else
6331  {
6332  // Was the node created by one of its neighbours
6333  // Whether or not the node lies on an edge can be calculated
6334  // by from the fractional position
6335  bool is_periodic = false;
6336  created_node_pt =
6337  this->node_created_by_neighbour(s_fraction, is_periodic);
6338 
6339  // If the node was so created, assign the pointers
6340  if (created_node_pt != 0)
6341  {
6342  // If the node is periodic
6343  if (is_periodic)
6344  {
6345  // Now the node must be on a boundary, but we don't know which
6346  // one
6347  // The returned created_node_pt is actually the neighbouring
6348  // periodic node
6349  Node* neighbour_node_pt = created_node_pt;
6350 
6351  // Determine the edge on which the new node will live
6352  //(cannot be a vertex node)
6353  int my_bound = Tree::OMEGA;
6354  // If we are on the left face
6355  if (i0 == 0)
6356  {
6357  my_bound = OcTreeNames::L;
6358  }
6359  // If we are on the right face
6360  else if (i0 == this->P_order - 1)
6361  {
6362  my_bound = OcTreeNames::R;
6363  }
6364 
6365  // If we are on the bottom face
6366  if (i1 == 0)
6367  {
6368  // If we already have already set the boundary, we're on an
6369  // edge
6370  switch (my_bound)
6371  {
6372  case OcTreeNames::L:
6373  my_bound = OcTreeNames::LD;
6374  break;
6375  case OcTreeNames::R:
6376  my_bound = OcTreeNames::RD;
6377  break;
6378  // Boundary not set
6379  default:
6380  my_bound = OcTreeNames::D;
6381  break;
6382  }
6383  }
6384  // If we are the top face
6385  else if (i1 == this->P_order - 1)
6386  {
6387  // If we already have a boundary
6388  switch (my_bound)
6389  {
6390  case OcTreeNames::L:
6391  my_bound = OcTreeNames::LU;
6392  break;
6393  case OcTreeNames::R:
6394  my_bound = OcTreeNames::RU;
6395  break;
6396  default:
6397  my_bound = OcTreeNames::U;
6398  break;
6399  }
6400  }
6401 
6402  // If we are on the back face
6403  if (i2 == 0)
6404  {
6405  // If we already have already set the boundary, we're on an
6406  // edge
6407  switch (my_bound)
6408  {
6409  case OcTreeNames::L:
6410  my_bound = OcTreeNames::LB;
6411  break;
6412  case OcTreeNames::LD:
6413  my_bound = OcTreeNames::LDB;
6414  break;
6415  case OcTreeNames::LU:
6416  my_bound = OcTreeNames::LUB;
6417  break;
6418  case OcTreeNames::R:
6419  my_bound = OcTreeNames::RB;
6420  break;
6421  case OcTreeNames::RD:
6422  my_bound = OcTreeNames::RDB;
6423  break;
6424  case OcTreeNames::RU:
6425  my_bound = OcTreeNames::RUB;
6426  break;
6427  case OcTreeNames::D:
6428  my_bound = OcTreeNames::DB;
6429  break;
6430  case OcTreeNames::U:
6431  my_bound = OcTreeNames::UB;
6432  break;
6433  // Boundary not set
6434  default:
6435  my_bound = OcTreeNames::B;
6436  break;
6437  }
6438  }
6439  // If we are the front face
6440  else if (i2 == this->P_order - 1)
6441  {
6442  // If we already have a boundary
6443  switch (my_bound)
6444  {
6445  case OcTreeNames::L:
6446  my_bound = OcTreeNames::LF;
6447  break;
6448  case OcTreeNames::LD:
6449  my_bound = OcTreeNames::LDF;
6450  break;
6451  case OcTreeNames::LU:
6452  my_bound = OcTreeNames::LUF;
6453  break;
6454  case OcTreeNames::R:
6455  my_bound = OcTreeNames::RF;
6456  break;
6457  case OcTreeNames::RD:
6458  my_bound = OcTreeNames::RDF;
6459  break;
6460  case OcTreeNames::RU:
6461  my_bound = OcTreeNames::RUF;
6462  break;
6463  case OcTreeNames::D:
6464  my_bound = OcTreeNames::DF;
6465  break;
6466  case OcTreeNames::U:
6467  my_bound = OcTreeNames::UF;
6468  break;
6469  default:
6470  my_bound = OcTreeNames::F;
6471  break;
6472  }
6473  }
6474 
6475  // Storage for the set of Mesh boundaries on which the
6476  // appropriate edge lives.
6477  // [New nodes should always be mid-edge nodes and therefore
6478  // only live on one boundary but just to play it safe...]
6479  std::set<unsigned> boundaries;
6480  // Only get the boundaries if we are at the edge of
6481  // an element. Nodes in the centre of an element cannot be
6482  // on Mesh boundaries
6483  if (my_bound != Tree::OMEGA)
6484  {
6485  clone_el_pt->get_boundaries(my_bound, boundaries);
6486  }
6487 
6488 #ifdef PARANOID
6489  // Case where a new node lives on more than one boundary
6490  // seems fishy enough to flag
6491  if (boundaries.size() > 2)
6492  {
6493  throw OomphLibError(
6494  "boundaries.size()>2 seems a bit strange..\n",
6495  OOMPH_CURRENT_FUNCTION,
6496  OOMPH_EXCEPTION_LOCATION);
6497  }
6498 
6499  // Case when there are no boundaries, we are in big trouble
6500  if (boundaries.size() == 0)
6501  {
6502  std::ostringstream error_stream;
6503  error_stream << "Periodic node is not on a boundary...\n"
6504  << "Coordinates: " << created_node_pt->x(0)
6505  << " " << created_node_pt->x(1) << "\n";
6506  throw OomphLibError(error_stream.str(),
6507  OOMPH_CURRENT_FUNCTION,
6508  OOMPH_EXCEPTION_LOCATION);
6509  }
6510 #endif
6511 
6512  // Create node and set the pointer to it from the element
6513  created_node_pt =
6514  this->construct_boundary_node(jnod, time_stepper_pt);
6515  // Make the node periodic from the neighbour
6516  created_node_pt->make_periodic(neighbour_node_pt);
6517 
6518  // Loop over # of history values
6519  for (unsigned t = 0; t < ntstorage; t++)
6520  {
6521  // Get position from father element -- this uses the macro
6522  // element representation if appropriate. If the node
6523  // turns out to be a hanging node later on, then
6524  // its position gets adjusted in line with its
6525  // hanging node interpolation.
6526  Vector<double> x_prev(n_dim, 0.0);
6527  clone_el_pt->get_x(t, s, x_prev);
6528  // Set previous positions of the new node
6529  for (unsigned i = 0; i < n_dim; i++)
6530  {
6531  created_node_pt->x(t, i) = x_prev[i];
6532  }
6533  }
6534 
6535  // Check if we need to add nodes to the mesh
6536  if (mesh_pt != 0)
6537  {
6538  // Next, we Update the boundary lookup schemes
6539  // Loop over the boundaries stored in the set
6540  for (std::set<unsigned>::iterator it = boundaries.begin();
6541  it != boundaries.end();
6542  ++it)
6543  {
6544  // Add the node to the boundary
6545  mesh_pt->add_boundary_node(*it, created_node_pt);
6546 
6547  // If we have set an intrinsic coordinate on this
6548  // mesh boundary then it must also be interpolated on
6549  // the new node
6550  // Now interpolate the intrinsic boundary coordinate
6551  if (mesh_pt->boundary_coordinate_exists(*it) == true)
6552  {
6553  Vector<double> zeta(2, 0.0);
6554  clone_el_pt->interpolated_zeta_on_face(
6555  *it, my_bound, s, zeta);
6556 
6557  created_node_pt->set_coordinates_on_boundary(*it, zeta);
6558  }
6559  }
6560 
6561  // Make sure that we add the node to the mesh
6562  mesh_pt->add_node_pt(created_node_pt);
6563  }
6564  } // End of periodic case
6565  // Otherwise the node is not periodic, so just set the
6566  // pointer to the neighbours node
6567  else
6568  {
6569  this->node_pt(jnod) = created_node_pt;
6570  }
6571  // Node has been created
6572  node_done = true;
6573  }
6574  // Node does not exist in neighbour element but might already
6575  //-----------------------------------------------------------
6576  // have been created by a son of a neighbouring element
6577  //-----------------------------------------------------
6578  else
6579  {
6580  // Was the node created by one of its neighbours' sons
6581  // Whether or not the node lies on an edge can be calculated
6582  // by from the fractional position
6583  bool is_periodic = false;
6584  created_node_pt =
6585  this->node_created_by_son_of_neighbour(s_fraction, is_periodic);
6586 
6587  // If the node was so created, assign the pointers
6588  if (created_node_pt != 0)
6589  {
6590  // If the node is periodic
6591  if (is_periodic)
6592  {
6593  // Now the node must be on a boundary, but we don't know which
6594  // one
6595  // The returned created_node_pt is actually the neighbouring
6596  // periodic node
6597  Node* neighbour_node_pt = created_node_pt;
6598 
6599  // Determine the edge on which the new node will live
6600  //(cannot be a vertex node)
6601  int my_bound = Tree::OMEGA;
6602  // If we are on the left face
6603  if (i0 == 0)
6604  {
6605  my_bound = OcTreeNames::L;
6606  }
6607  // If we are on the right face
6608  else if (i0 == this->P_order - 1)
6609  {
6610  my_bound = OcTreeNames::R;
6611  }
6612 
6613  // If we are on the bottom face
6614  if (i1 == 0)
6615  {
6616  // If we already have already set the boundary, we're on an
6617  // edge
6618  switch (my_bound)
6619  {
6620  case OcTreeNames::L:
6621  my_bound = OcTreeNames::LD;
6622  break;
6623  case OcTreeNames::R:
6624  my_bound = OcTreeNames::RD;
6625  break;
6626  // Boundary not set
6627  default:
6628  my_bound = OcTreeNames::D;
6629  break;
6630  }
6631  }
6632  // If we are the top face
6633  else if (i1 == this->P_order - 1)
6634  {
6635  // If we already have a boundary
6636  switch (my_bound)
6637  {
6638  case OcTreeNames::L:
6639  my_bound = OcTreeNames::LU;
6640  break;
6641  case OcTreeNames::R:
6642  my_bound = OcTreeNames::RU;
6643  break;
6644  default:
6645  my_bound = OcTreeNames::U;
6646  break;
6647  }
6648  }
6649 
6650  // If we are on the back face
6651  if (i2 == 0)
6652  {
6653  // If we already have already set the boundary, we're on an
6654  // edge
6655  switch (my_bound)
6656  {
6657  case OcTreeNames::L:
6658  my_bound = OcTreeNames::LB;
6659  break;
6660  case OcTreeNames::LD:
6661  my_bound = OcTreeNames::LDB;
6662  break;
6663  case OcTreeNames::LU:
6664  my_bound = OcTreeNames::LUB;
6665  break;
6666  case OcTreeNames::R:
6667  my_bound = OcTreeNames::RB;
6668  break;
6669  case OcTreeNames::RD:
6670  my_bound = OcTreeNames::RDB;
6671  break;
6672  case OcTreeNames::RU:
6673  my_bound = OcTreeNames::RUB;
6674  break;
6675  case OcTreeNames::D:
6676  my_bound = OcTreeNames::DB;
6677  break;
6678  case OcTreeNames::U:
6679  my_bound = OcTreeNames::UB;
6680  break;
6681  // Boundary not set
6682  default:
6683  my_bound = OcTreeNames::B;
6684  break;
6685  }
6686  }
6687  // If we are the front face
6688  else if (i2 == this->P_order - 1)
6689  {
6690  // If we already have a boundary
6691  switch (my_bound)
6692  {
6693  case OcTreeNames::L:
6694  my_bound = OcTreeNames::LF;
6695  break;
6696  case OcTreeNames::LD:
6697  my_bound = OcTreeNames::LDF;
6698  break;
6699  case OcTreeNames::LU:
6700  my_bound = OcTreeNames::LUF;
6701  break;
6702  case OcTreeNames::R:
6703  my_bound = OcTreeNames::RF;
6704  break;
6705  case OcTreeNames::RD:
6706  my_bound = OcTreeNames::RDF;
6707  break;
6708  case OcTreeNames::RU:
6709  my_bound = OcTreeNames::RUF;
6710  break;
6711  case OcTreeNames::D:
6712  my_bound = OcTreeNames::DF;
6713  break;
6714  case OcTreeNames::U:
6715  my_bound = OcTreeNames::UF;
6716  break;
6717  default:
6718  my_bound = OcTreeNames::F;
6719  break;
6720  }
6721  }
6722 
6723  // Storage for the set of Mesh boundaries on which the
6724  // appropriate edge lives.
6725  // [New nodes should always be mid-edge nodes and therefore
6726  // only live on one boundary but just to play it safe...]
6727  std::set<unsigned> boundaries;
6728  // Only get the boundaries if we are at the edge of
6729  // an element. Nodes in the centre of an element cannot be
6730  // on Mesh boundaries
6731  if (my_bound != Tree::OMEGA)
6732  {
6733  clone_el_pt->get_boundaries(my_bound, boundaries);
6734  }
6735 
6736 #ifdef PARANOID
6737  // Case where a new node lives on more than one boundary
6738  // seems fishy enough to flag
6739  if (boundaries.size() > 2)
6740  {
6741  throw OomphLibError(
6742  "boundaries.size()>2 seems a bit strange..\n",
6743  OOMPH_CURRENT_FUNCTION,
6744  OOMPH_EXCEPTION_LOCATION);
6745  }
6746 
6747  // Case when there are no boundaries, we are in big trouble
6748  if (boundaries.size() == 0)
6749  {
6750  std::ostringstream error_stream;
6751  error_stream << "Periodic node is not on a boundary...\n"
6752  << "Coordinates: " << created_node_pt->x(0)
6753  << " " << created_node_pt->x(1) << "\n";
6754  throw OomphLibError(error_stream.str(),
6755  OOMPH_CURRENT_FUNCTION,
6756  OOMPH_EXCEPTION_LOCATION);
6757  }
6758 #endif
6759 
6760  // Create node and set the pointer to it from the element
6761  created_node_pt =
6762  this->construct_boundary_node(jnod, time_stepper_pt);
6763  // Make the node periodic from the neighbour
6764  created_node_pt->make_periodic(neighbour_node_pt);
6765 
6766  // Loop over # of history values
6767  for (unsigned t = 0; t < ntstorage; t++)
6768  {
6769  // Get position from father element -- this uses the macro
6770  // element representation if appropriate. If the node
6771  // turns out to be a hanging node later on, then
6772  // its position gets adjusted in line with its
6773  // hanging node interpolation.
6774  Vector<double> x_prev(n_dim);
6775  clone_el_pt->get_x(t, s, x_prev);
6776  // Set previous positions of the new node
6777  for (unsigned i = 0; i < n_dim; i++)
6778  {
6779  created_node_pt->x(t, i) = x_prev[i];
6780  }
6781  }
6782 
6783  // Check if we need to add nodes to the mesh
6784  if (mesh_pt != 0)
6785  {
6786  // Next, we Update the boundary lookup schemes
6787  // Loop over the boundaries stored in the set
6788  for (std::set<unsigned>::iterator it = boundaries.begin();
6789  it != boundaries.end();
6790  ++it)
6791  {
6792  // Add the node to the boundary
6793  mesh_pt->add_boundary_node(*it, created_node_pt);
6794 
6795  // If we have set an intrinsic coordinate on this
6796  // mesh boundary then it must also be interpolated on
6797  // the new node
6798  // Now interpolate the intrinsic boundary coordinate
6799  if (mesh_pt->boundary_coordinate_exists(*it) == true)
6800  {
6801  Vector<double> zeta(2, 0.0);
6802  clone_el_pt->interpolated_zeta_on_face(
6803  *it, my_bound, s, zeta);
6804 
6805  created_node_pt->set_coordinates_on_boundary(*it, zeta);
6806  }
6807  }
6808 
6809  // Make sure that we add the node to the mesh
6810  mesh_pt->add_node_pt(created_node_pt);
6811  }
6812  } // End of periodic case
6813  // Otherwise the node is not periodic, so just set the
6814  // pointer to the neighbours node
6815  else
6816  {
6817  this->node_pt(jnod) = created_node_pt;
6818  }
6819  // Node has been created
6820  node_done = true;
6821  } // Node does not exist in son of neighbouring element
6822  } // Node does not exist in neighbouring element
6823  } // Node does not exist in this element
6824 
6825  // Node has not been built anywhere ---> build it here
6826  if (!node_done)
6827  {
6828  // Firstly, we need to determine whether or not a node lies
6829  // on the boundary before building it, because
6830  // we actually assign a different type of node on boundaries.
6831 
6832  // Initially none
6833  int my_bound = Tree::OMEGA;
6834  // If we are on the left face
6835  if (i0 == 0)
6836  {
6837  my_bound = OcTreeNames::L;
6838  }
6839  // If we are on the right face
6840  else if (i0 == this->P_order - 1)
6841  {
6842  my_bound = OcTreeNames::R;
6843  }
6844 
6845  // If we are on the bottom face
6846  if (i1 == 0)
6847  {
6848  // If we already have already set the boundary, we're on an edge
6849  switch (my_bound)
6850  {
6851  case OcTreeNames::L:
6852  my_bound = OcTreeNames::LD;
6853  break;
6854  case OcTreeNames::R:
6855  my_bound = OcTreeNames::RD;
6856  break;
6857  // Boundary not set
6858  default:
6859  my_bound = OcTreeNames::D;
6860  break;
6861  }
6862  }
6863  // If we are the top face
6864  else if (i1 == this->P_order - 1)
6865  {
6866  // If we already have a boundary
6867  switch (my_bound)
6868  {
6869  case OcTreeNames::L:
6870  my_bound = OcTreeNames::LU;
6871  break;
6872  case OcTreeNames::R:
6873  my_bound = OcTreeNames::RU;
6874  break;
6875  default:
6876  my_bound = OcTreeNames::U;
6877  break;
6878  }
6879  }
6880 
6881  // If we are on the back face
6882  if (i2 == 0)
6883  {
6884  // If we already have already set the boundary, we're on an edge
6885  switch (my_bound)
6886  {
6887  case OcTreeNames::L:
6888  my_bound = OcTreeNames::LB;
6889  break;
6890  case OcTreeNames::LD:
6891  my_bound = OcTreeNames::LDB;
6892  break;
6893  case OcTreeNames::LU:
6894  my_bound = OcTreeNames::LUB;
6895  break;
6896  case OcTreeNames::R:
6897  my_bound = OcTreeNames::RB;
6898  break;
6899  case OcTreeNames::RD:
6900  my_bound = OcTreeNames::RDB;
6901  break;
6902  case OcTreeNames::RU:
6903  my_bound = OcTreeNames::RUB;
6904  break;
6905  case OcTreeNames::D:
6906  my_bound = OcTreeNames::DB;
6907  break;
6908  case OcTreeNames::U:
6909  my_bound = OcTreeNames::UB;
6910  break;
6911  // Boundary not set
6912  default:
6913  my_bound = OcTreeNames::B;
6914  break;
6915  }
6916  }
6917  // If we are the front face
6918  else if (i2 == this->P_order - 1)
6919  {
6920  // If we already have a boundary
6921  switch (my_bound)
6922  {
6923  case OcTreeNames::L:
6924  my_bound = OcTreeNames::LF;
6925  break;
6926  case OcTreeNames::LD:
6927  my_bound = OcTreeNames::LDF;
6928  break;
6929  case OcTreeNames::LU:
6930  my_bound = OcTreeNames::LUF;
6931  break;
6932  case OcTreeNames::R:
6933  my_bound = OcTreeNames::RF;
6934  break;
6935  case OcTreeNames::RD:
6936  my_bound = OcTreeNames::RDF;
6937  break;
6938  case OcTreeNames::RU:
6939  my_bound = OcTreeNames::RUF;
6940  break;
6941  case OcTreeNames::D:
6942  my_bound = OcTreeNames::DF;
6943  break;
6944  case OcTreeNames::U:
6945  my_bound = OcTreeNames::UF;
6946  break;
6947  default:
6948  my_bound = OcTreeNames::F;
6949  break;
6950  }
6951  }
6952 
6953  // Storage for the set of Mesh boundaries on which the
6954  // appropriate edge lives.
6955  // [New nodes should always be mid-edge nodes and therefore
6956  // only live on one boundary but just to play it safe...]
6957  std::set<unsigned> boundaries;
6958  // Only get the boundaries if we are at the edge of
6959  // an element. Nodes in the centre of an element cannot be
6960  // on Mesh boundaries
6961  if (my_bound != Tree::OMEGA)
6962  clone_el_pt->get_boundaries(my_bound, boundaries);
6963 
6964 #ifdef PARANOID
6965  // Case where a new node lives on more than two boundaries
6966  // seems fishy enough to flag
6967  if (boundaries.size() > 2)
6968  {
6969  throw OomphLibError(
6970  "boundaries.size()>2 seems a bit strange..\n",
6971  "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
6972  OOMPH_EXCEPTION_LOCATION);
6973  }
6974 #endif
6975 
6976  // If the node lives on a mesh boundary,
6977  // then we need to create a boundary node
6978  if (boundaries.size() > 0)
6979  {
6980  // Create node and set the pointer to it from the element
6981  created_node_pt =
6982  this->construct_boundary_node(jnod, time_stepper_pt);
6983 
6984  // Now we need to work out whether to pin the values at
6985  // the new node based on the boundary conditions applied at
6986  // its Mesh boundary
6987 
6988  // Get the boundary conditions from the father
6989  Vector<int> bound_cons(this->ncont_interpolated_values());
6990  clone_el_pt->get_bcs(my_bound, bound_cons);
6991 
6992  // Loop over the values and pin, if necessary
6993  unsigned n_value = created_node_pt->nvalue();
6994  for (unsigned k = 0; k < n_value; k++)
6995  {
6996  if (bound_cons[k])
6997  {
6998  created_node_pt->pin(k);
6999  }
7000  }
7001 
7002  // Solid node? If so, deal with the positional boundary
7003  // conditions:
7004  SolidNode* solid_node_pt =
7005  dynamic_cast<SolidNode*>(created_node_pt);
7006  if (solid_node_pt != 0)
7007  {
7008  // Get the positional boundary conditions from the father:
7009  unsigned n_dim = created_node_pt->ndim();
7010  Vector<int> solid_bound_cons(n_dim);
7011  RefineableSolidQElement<3>* clone_solid_el_pt =
7012  dynamic_cast<RefineableSolidQElement<3>*>(clone_el_pt);
7013 #ifdef PARANOID
7014  if (clone_solid_el_pt == 0)
7015  {
7016  std::string error_message =
7017  "We have a SolidNode outside a refineable SolidElement\n";
7018  error_message +=
7019  "during mesh refinement -- this doesn't make sense";
7020 
7021  throw OomphLibError(
7022  error_message,
7023  "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
7024  OOMPH_EXCEPTION_LOCATION);
7025  }
7026 #endif
7027  clone_solid_el_pt->get_solid_bcs(my_bound, solid_bound_cons);
7028 
7029  // Loop over the positions and pin, if necessary
7030  for (unsigned k = 0; k < n_dim; k++)
7031  {
7032  if (solid_bound_cons[k])
7033  {
7034  solid_node_pt->pin_position(k);
7035  }
7036  }
7037  } // End of if solid_node_pt
7038 
7039  // Next, we Update the boundary lookup schemes
7040 
7041  // Check if we need to add nodes to the mesh
7042  if (mesh_pt != 0)
7043  {
7044  // Loop over the boundaries stored in the set
7045  for (std::set<unsigned>::iterator it = boundaries.begin();
7046  it != boundaries.end();
7047  ++it)
7048  {
7049  // Add the node to the boundary
7050  mesh_pt->add_boundary_node(*it, created_node_pt);
7051 
7052  // If we have set an intrinsic coordinate on this
7053  // mesh boundary then it must also be interpolated on
7054  // the new node
7055  // Now interpolate the intrinsic boundary coordinate
7056  if (mesh_pt->boundary_coordinate_exists(*it) == true)
7057  {
7058  Vector<double> zeta(2);
7059  clone_el_pt->interpolated_zeta_on_face(
7060  *it, my_bound, s, zeta);
7061 
7062  created_node_pt->set_coordinates_on_boundary(*it, zeta);
7063  }
7064  }
7065  }
7066  }
7067  // Otherwise the node is not on a Mesh boundary and
7068  // we create a normal "bulk" node
7069  else
7070  {
7071  // Create node and set the pointer to it from the element
7072  created_node_pt = this->construct_node(jnod, time_stepper_pt);
7073  }
7074 
7075  // Now we set the position and values at the newly created node
7076 
7077  // In the first instance use macro element or FE representation
7078  // to create past and present nodal positions.
7079  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
7080  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
7081  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
7082  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
7083  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
7084  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
7085  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
7086 
7087  // Loop over # of history values
7088  for (unsigned t = 0; t < ntstorage; t++)
7089  {
7090  // Get position from father element -- this uses the macro
7091  // element representation if appropriate. If the node
7092  // turns out to be a hanging node later on, then
7093  // its position gets adjusted in line with its
7094  // hanging node interpolation.
7095  Vector<double> x_prev(3);
7096  clone_el_pt->get_x(t, s, x_prev);
7097 
7098  // Set previous positions of the new node
7099  for (unsigned i = 0; i < 3; i++)
7100  {
7101  created_node_pt->x(t, i) = x_prev[i];
7102  }
7103  }
7104 
7105  // Loop over all history values
7106  for (unsigned t = 0; t < ntstorage; t++)
7107  {
7108  // Get values from father element
7109  // Note: get_interpolated_values() sets Vector size itself.
7110  Vector<double> prev_values;
7111  clone_el_pt->get_interpolated_values(t, s, prev_values);
7112 
7113  // Initialise the values at the new node
7114  unsigned n_value = created_node_pt->nvalue();
7115  for (unsigned k = 0; k < n_value; k++)
7116  {
7117  created_node_pt->set_value(t, k, prev_values[k]);
7118  }
7119  }
7120 
7121  // Add new node to mesh (if requested)
7122  if (mesh_pt != 0)
7123  {
7124  mesh_pt->add_node_pt(created_node_pt);
7125  }
7126 
7127  AlgebraicElementBase* alg_el_pt =
7128  dynamic_cast<AlgebraicElementBase*>(this);
7129 
7130  // If we do have an algebraic element
7131  if (alg_el_pt != 0)
7132  {
7133  std::string error_message =
7134  "Have not implemented p-refinement for";
7135  error_message += "Algebraic p-refineable elements yet\n";
7136 
7137  throw OomphLibError(
7138  error_message,
7139  "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
7140  OOMPH_EXCEPTION_LOCATION);
7141  }
7142 
7143  } // End of case when we build the node ourselves
7144 
7145  // Check if the element is an algebraic element
7146  AlgebraicElementBase* alg_el_pt =
7147  dynamic_cast<AlgebraicElementBase*>(this);
7148 
7149  // If the element is an algebraic element, setup
7150  // node position (past and present) from algebraic node update
7151  // function. This over-writes previous assingments that
7152  // were made based on the macro-element/FE representation.
7153  // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
7154  // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
7155  // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
7156  // INFO FOR *ALL* ROOT ELEMENTS!
7157  if (alg_el_pt != 0)
7158  {
7159  // Build algebraic node update info for new node
7160  // This sets up the node update data for all node update
7161  // functions that are shared by all nodes in the father
7162  // element
7163  alg_el_pt->setup_algebraic_node_update(
7164  this->node_pt(jnod), s, clone_el_pt);
7165  }
7166 
7167  } // End of vertical loop over nodes in element (i2)
7168 
7169  } // End of horizontal loop over nodes in element (i1)
7170 
7171  } // End of horizontal loop over nodes in element (i0)
7172 
7173 
7174  // Loop over all nodes in element again, to re-set the positions
7175  // This must be done using the new element's macro-element
7176  // representation, rather than the old version which may be
7177  // of a different p-order!
7178  for (unsigned i0 = 0; i0 < this->P_order; i0++)
7179  {
7180  // Get the fractional position of the node in the direction of s[0]
7181  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
7182  // Local coordinate
7183  s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
7184 
7185  for (unsigned i1 = 0; i1 < this->P_order; i1++)
7186  {
7187  // Get the fractional position of the node in the direction of s[1]
7188  s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
7189  // Local coordinate
7190  s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
7191 
7192  for (unsigned i2 = 0; i2 < this->P_order; i2++)
7193  {
7194  // Get the fractional position of the node in the direction of s[2]
7195  s_fraction[2] = this->local_one_d_fraction_of_node(i2, 2);
7196  // Local coordinate
7197  s[2] = s_lo[2] + (s_hi[2] - s_lo[2]) * s_fraction[2];
7198 
7199  // Local node number
7200  jnod = i0 + this->P_order * i1 + this->P_order * this->P_order * i2;
7201 
7202  // Loop over # of history values
7203  for (unsigned t = 0; t < ntstorage; t++)
7204  {
7205  // Get position from father element -- this uses the macro
7206  // element representation if appropriate. If the node
7207  // turns out to be a hanging node later on, then
7208  // its position gets adjusted in line with its
7209  // hanging node interpolation.
7210  Vector<double> x_prev(3);
7211  this->get_x(t, s, x_prev);
7212 
7213  // Set previous positions of the new node
7214  for (unsigned i = 0; i < 3; i++)
7215  {
7216  this->node_pt(jnod)->x(t, i) = x_prev[i];
7217  }
7218  }
7219  }
7220  }
7221  }
7222 
7223 
7224  // If the element is a MacroElementNodeUpdateElement, set
7225  // the update parameters for the current element's nodes --
7226  // all this needs is the vector of (pointers to the)
7227  // geometric objects that affect the MacroElement-based
7228  // node update -- this needs to be done to set the node
7229  // update info for newly created nodes
7230  MacroElementNodeUpdateElementBase* clone_m_el_pt =
7231  dynamic_cast<MacroElementNodeUpdateElementBase*>(clone_el_pt);
7232  if (clone_m_el_pt != 0)
7233  {
7234  // Get vector of geometric objects from father (construct vector
7235  // via copy operation)
7236  Vector<GeomObject*> geom_object_pt(clone_m_el_pt->geom_object_pt());
7237 
7238  // Cast current element to MacroElementNodeUpdateElement:
7239  MacroElementNodeUpdateElementBase* m_el_pt =
7240  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
7241 
7242 #ifdef PARANOID
7243  if (m_el_pt == 0)
7244  {
7245  std::string error_message =
7246  "Failed to cast to MacroElementNodeUpdateElementBase*\n";
7247  error_message +=
7248  "Strange -- if my clone is a MacroElementNodeUpdateElement\n";
7249  error_message += "then I should be too....\n";
7250 
7251  throw OomphLibError(
7252  error_message,
7253  "PRefineableQElement<3,INITIAL_NNODE_1D>::p_refine()",
7254  OOMPH_EXCEPTION_LOCATION);
7255  }
7256 #endif
7257  // Build update info by passing vector of geometric objects:
7258  // This sets the current element to be the update element
7259  // for all of the element's nodes -- this is reversed
7260  // if the element is ever un-refined in the father element's
7261  // rebuild_from_sons() function which overwrites this
7262  // assignment to avoid nasty segmentation faults that occur
7263  // when a node tries to update itself via an element that no
7264  // longer exists...
7265  m_el_pt->set_node_update_info(geom_object_pt);
7266  }
7267 
7268  // Not necessary to delete the old nodes since all original nodes are in the
7269  // current mesh and so will be pruned as part of the mesh adaption process.
7270 
7271  // Do any further-build required
7272  this->further_build();
7273  }
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< 3, 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:4956
oomph::PRefineableQElement< 3, 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:4847
oomph::PRefineableQElement< 3, 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:5271
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
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::OcTreeNames::B
@ B
Definition: octree.h:73
oomph::OcTreeNames::D
@ D
Definition: octree.h:71
oomph::OcTreeNames::L
@ L
Definition: octree.h:69
oomph::OcTreeNames::R
@ R
Definition: octree.h:70
oomph::OcTreeNames::U
@ U
Definition: octree.h:72
plotPSD.t
t
Definition: plotPSD.py:36

References oomph::Mesh::add_boundary_node(), oomph::Mesh::add_node_pt(), oomph::OcTreeNames::B, oomph::Mesh::boundary_coordinate_exists(), oomph::OcTreeNames::D, oomph::OcTreeNames::DB, oomph::OcTreeNames::DF, oomph::OcTreeNames::F, oomph::MacroElementNodeUpdateElementBase::geom_object_pt(), oomph::RefineableQElement< 3 >::get_bcs(), oomph::RefineableQElement< 3 >::get_boundaries(), oomph::RefineableElement::get_interpolated_values(), get_node_at_local_coordinate(), oomph::RefineableSolidQElement< 3 >::get_solid_bcs(), oomph::FiniteElement::get_x(), i, oomph::RefineableQElement< 3 >::interpolated_zeta_on_face(), k, oomph::OcTreeNames::L, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::OcTreeNames::LDB, oomph::OcTreeNames::LDF, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, oomph::OcTreeNames::LUB, oomph::OcTreeNames::LUF, oomph::Node::make_periodic(), 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(), oomph::OcTreeNames::R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, oomph::OcTreeNames::RDB, oomph::OcTreeNames::RDF, oomph::OcTreeNames::RF, oomph::OcTreeNames::RU, oomph::OcTreeNames::RUB, oomph::OcTreeNames::RUF, 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::OcTreeNames::U, oomph::OcTreeNames::UB, oomph::OcTreeNames::UF, oomph::Node::x(), and Eigen::zeta().

◆ pre_build()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 3, 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

/ Pass pointer to time object:

/ Pass macro element pointer on to sons and / set coordinates in macro element / hierher why can I see this?

Reimplemented from oomph::RefineableElement.

5862  {
5863  using namespace OcTreeNames;
5864 
5865  // Get the number of 1d nodes
5866  unsigned n_p = nnode_1d();
5867 
5868  // Check whether static father_bound needs to be created
5869  if (Father_bound[n_p].nrow() == 0)
5870  {
5871  setup_father_bounds();
5872  }
5873 
5874  // Pointer to my father (in quadtree impersonation)
5875  OcTree* father_pt = dynamic_cast<OcTree*>(octree_pt()->father_pt());
5876 
5877  // What type of son am I? Ask my quadtree representation...
5878  int son_type = Tree_pt->son_type();
5879 
5880  // Has somebody build me already? (If any nodes have not been built)
5881  if (!nodes_built())
5882  {
5883 #ifdef PARANOID
5884  if (father_pt == 0)
5885  {
5886  std::string error_message =
5887  "Something fishy here: I have no father and yet \n";
5888  error_message += "I have no nodes. Who has created me then?!\n";
5889 
5890  throw OomphLibError(
5891  error_message,
5892  "PRefineableQElement<3,INITIAL_NNODE_1D>::pre_build()",
5893  OOMPH_EXCEPTION_LOCATION);
5894  }
5895 #endif
5896 
5897  // Return pointer to father element
5898  RefineableQElement<3>* father_el_pt =
5899  dynamic_cast<RefineableQElement<3>*>(father_pt->object_pt());
5900 
5901  // Timestepper should be the same for all nodes in father
5902  // element -- use it create timesteppers for new nodes
5903  TimeStepper* time_stepper_pt =
5904  father_el_pt->node_pt(0)->time_stepper_pt();
5905 
5906  // Number of history values (incl. present)
5907  unsigned ntstorage = time_stepper_pt->ntstorage();
5908 
5910  // this->time_pt()=father_el_pt->time_pt();
5911 
5912  Vector<double> s_lo(3);
5913  Vector<double> s_hi(3);
5914  Vector<double> s(3);
5915  Vector<double> x(3);
5916 
5917  // Setup vertex coordinates in father element:
5918  //--------------------------------------------
5919  switch (son_type)
5920  {
5921  case LDB:
5922  s_lo[0] = -1.0;
5923  s_hi[0] = 0.0;
5924  s_lo[1] = -1.0;
5925  s_hi[1] = 0.0;
5926  s_lo[2] = -1.0;
5927  s_hi[2] = 0.0;
5928  break;
5929 
5930  case LDF:
5931  s_lo[0] = -1.0;
5932  s_hi[0] = 0.0;
5933  s_lo[1] = -1.0;
5934  s_hi[1] = 0.0;
5935  s_lo[2] = 0.0;
5936  s_hi[2] = 1.0;
5937  break;
5938 
5939  case LUB:
5940  s_lo[0] = -1.0;
5941  s_hi[0] = 0.0;
5942  s_lo[1] = 0.0;
5943  s_hi[1] = 1.0;
5944  s_lo[2] = -1.0;
5945  s_hi[2] = 0.0;
5946  break;
5947 
5948  case LUF:
5949  s_lo[0] = -1.0;
5950  s_hi[0] = 0.0;
5951  s_lo[1] = 0.0;
5952  s_hi[1] = 1.0;
5953  s_lo[2] = 0.0;
5954  s_hi[2] = 1.0;
5955  break;
5956 
5957  case RDB:
5958  s_lo[0] = 0.0;
5959  s_hi[0] = 1.0;
5960  s_lo[1] = -1.0;
5961  s_hi[1] = 0.0;
5962  s_lo[2] = -1.0;
5963  s_hi[2] = 0.0;
5964  break;
5965 
5966  case RDF:
5967  s_lo[0] = 0.0;
5968  s_hi[0] = 1.0;
5969  s_lo[1] = -1.0;
5970  s_hi[1] = 0.0;
5971  s_lo[2] = 0.0;
5972  s_hi[2] = 1.0;
5973  break;
5974 
5975  case RUB:
5976  s_lo[0] = 0.0;
5977  s_hi[0] = 1.0;
5978  s_lo[1] = 0.0;
5979  s_hi[1] = 1.0;
5980  s_lo[2] = -1.0;
5981  s_hi[2] = 0.0;
5982  break;
5983 
5984  case RUF:
5985  s_lo[0] = 0.0;
5986  s_hi[0] = 1.0;
5987  s_lo[1] = 0.0;
5988  s_hi[1] = 1.0;
5989  s_lo[2] = 0.0;
5990  s_hi[2] = 1.0;
5991  break;
5992  }
5993 
5997  // if(father_el_pt->macro_elem_pt()!=0)
5998  // {
5999  // set_macro_elem_pt(father_el_pt->macro_elem_pt());
6000  // for(unsigned i=0;i<2;i++)
6001  // {
6002  // s_macro_ll(i)= father_el_pt->s_macro_ll(i)+
6003  // 0.5*(s_lo[i]+1.0)*(father_el_pt->s_macro_ur(i)-
6004  // father_el_pt->s_macro_ll(i));
6005  // s_macro_ur(i)= father_el_pt->s_macro_ll(i)+
6006  // 0.5*(s_hi[i]+1.0)*(father_el_pt->s_macro_ur(i)-
6007  // father_el_pt->s_macro_ll(i));
6008  // }
6009  // }
6010 
6011 
6012  // If the father element hasn't been generated yet, we're stuck...
6013  if (father_el_pt->node_pt(0) == 0)
6014  {
6015  throw OomphLibError(
6016  "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
6017  "PRefineableQElement<3,INITIAL_NNODE_1D>::pre_build()",
6018  OOMPH_EXCEPTION_LOCATION);
6019  }
6020  else
6021  {
6022  unsigned jnod = 0;
6023 
6024  Vector<double> s_fraction(3);
6025  // Loop over nodes in element
6026  for (unsigned i0 = 0; i0 < n_p; i0++)
6027  {
6028  // Get the fractional position of the node in the direction of s[0]
6029  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
6030  // Local coordinate in father element
6031  s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
6032 
6033  for (unsigned i1 = 0; i1 < n_p; i1++)
6034  {
6035  // Get the fractional position of the node in the direction of s[1]
6036  s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
6037  // Local coordinate in father element
6038  s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
6039 
6040  for (unsigned i2 = 0; i2 < n_p; i2++)
6041  {
6042  // Get the fractional position of the node in the direction of
6043  // s[2]
6044  s_fraction[2] = this->local_one_d_fraction_of_node(i2, 2);
6045  // Local coordinate in father element
6046  s[2] = s_lo[2] + (s_hi[2] - s_lo[2]) * s_fraction[2];
6047 
6048  // Local node number
6049  jnod = i0 + n_p * i1 + n_p * n_p * i2;
6050 
6051  // Get the pointer to the node in the father, returns NULL
6052  // if there is not node
6053  Node* created_node_pt =
6054  father_el_pt->get_node_at_local_coordinate(s);
6055 
6056  // Does this node already exist in father element?
6057  //------------------------------------------------
6058  if (created_node_pt != 0)
6059  {
6060  // Copy node across
6061  this->node_pt(jnod) = created_node_pt;
6062 
6063  // Make sure that we update the values at the node so that
6064  // they are consistent with the present representation.
6065  // This is only need for mixed interpolation where the value
6066  // at the father could now become active.
6067 
6068  // Loop over all history values
6069  for (unsigned t = 0; t < ntstorage; t++)
6070  {
6071  // Get values from father element
6072  // Note: get_interpolated_values() sets Vector size itself.
6073  Vector<double> prev_values;
6074  father_el_pt->get_interpolated_values(t, s, prev_values);
6075  // Find the minimum number of values
6076  //(either those stored at the node, or those returned by
6077  // the function)
6078  unsigned n_val_at_node = created_node_pt->nvalue();
6079  unsigned n_val_from_function = prev_values.size();
6080  // Use the ternary conditional operator here
6081  unsigned n_var = n_val_at_node < n_val_from_function ?
6082  n_val_at_node :
6083  n_val_from_function;
6084  // Assign the values that we can
6085  for (unsigned k = 0; k < n_var; k++)
6086  {
6087  created_node_pt->set_value(t, k, prev_values[k]);
6088  }
6089  }
6090  }
6091  } // End of loop i2
6092  } // End of loop i1
6093  } // End of loop i0
6094  } // Sanity check: Father element has been generated
6095 
6096  } // End of nodes not built
6097  }
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::RefineableQElement< 3 >::Father_bound
static std::map< unsigned, DenseMatrix< int > > Father_bound
Definition: refineable_brick_element.h:172
oomph::RefineableQElement< 3 >::setup_father_bounds
void setup_father_bounds()
Definition: refineable_brick_element.cc:126
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::RefineableElement::get_interpolated_values(), oomph::FiniteElement::get_node_at_local_coordinate(), k, oomph::OcTreeNames::LDB, oomph::OcTreeNames::LDF, oomph::OcTreeNames::LUB, oomph::OcTreeNames::LUF, oomph::FiniteElement::node_pt(), oomph::Data::nvalue(), oomph::Tree::object_pt(), OOMPH_EXCEPTION_LOCATION, oomph::OcTreeNames::RDB, oomph::OcTreeNames::RDF, oomph::OcTreeNames::RUB, oomph::OcTreeNames::RUF, s, oomph::Data::set_value(), oomph::Tree::son_type(), oomph::Global_string_for_annotation::string(), plotPSD::t, oomph::Data::time_stepper_pt(), and plotDoE::x.

◆ rebuild_from_sons()

template<unsigned INITIAL_NNODE_1D>
void oomph::PRefineableQElement< 3, 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.

7635  {
7636  // Get p-orders of sons
7637  unsigned n_sons = this->tree_pt()->nsons();
7638  Vector<unsigned> son_p_order(n_sons);
7639  unsigned max_son_p_order = 0;
7640  for (unsigned ison = 0; ison < n_sons; ison++)
7641  {
7642  PRefineableElement* el_pt = dynamic_cast<PRefineableElement*>(
7643  this->tree_pt()->son_pt(ison)->object_pt());
7644  son_p_order[ison] = el_pt->p_order();
7645  if (son_p_order[ison] > max_son_p_order)
7646  max_son_p_order = son_p_order[ison];
7647  }
7648 
7649  unsigned old_Nnode = this->nnode();
7650  unsigned old_P_order = this->p_order();
7651  // Set p-order of the element
7652  this->p_order() = max_son_p_order;
7653 
7654  // Change integration scheme
7655  delete this->integral_pt();
7656  switch (this->p_order())
7657  {
7658  case 2:
7659  this->set_integration_scheme(new GaussLobattoLegendre<3, 2>);
7660  break;
7661  case 3:
7662  this->set_integration_scheme(new GaussLobattoLegendre<3, 3>);
7663  break;
7664  case 4:
7665  this->set_integration_scheme(new GaussLobattoLegendre<3, 4>);
7666  break;
7667  case 5:
7668  this->set_integration_scheme(new GaussLobattoLegendre<3, 5>);
7669  break;
7670  case 6:
7671  this->set_integration_scheme(new GaussLobattoLegendre<3, 6>);
7672  break;
7673  case 7:
7674  this->set_integration_scheme(new GaussLobattoLegendre<3, 7>);
7675  break;
7676  default:
7677  std::ostringstream error_message;
7678  error_message << "\nERROR: Exceeded maximum polynomial order for";
7679  error_message << "\n integration scheme.\n";
7680  throw OomphLibError(
7681  error_message.str(),
7682  "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
7683  OOMPH_EXCEPTION_LOCATION);
7684  }
7685 
7686  // Back up pointers to old nodes before they are lost
7687  Vector<Node*> old_node_pt(old_Nnode);
7688  for (unsigned n = 0; n < old_Nnode; n++)
7689  {
7690  old_node_pt[n] = this->node_pt(n);
7691  }
7692 
7693  // Allocate new space for Nodes (at the element level)
7694  this->set_n_node(this->p_order() * this->p_order() * this->p_order());
7695 
7696  // Copy vertex nodes which were populated in the pre-build
7697  this->node_pt(0) = old_node_pt[0];
7698  this->node_pt(this->p_order() - 1) = old_node_pt[old_P_order - 1];
7699  this->node_pt(this->p_order() * (this->p_order() - 1)) =
7700  old_node_pt[(old_P_order) * (old_P_order - 1)];
7701  this->node_pt(this->p_order() * this->p_order() - 1) =
7702  old_node_pt[(old_P_order) * (old_P_order)-1];
7703  this->node_pt(this->p_order() * this->p_order() * (this->p_order() - 1)) =
7704  old_node_pt[old_P_order * old_P_order * (old_P_order - 1)];
7705  this->node_pt((this->p_order() * this->p_order() + 1) *
7706  (this->p_order() - 1)) =
7707  old_node_pt[(old_P_order * old_P_order + 1) * (old_P_order - 1)];
7708  this->node_pt(this->p_order() * (this->p_order() + 1) *
7709  (this->p_order() - 1)) =
7710  old_node_pt[old_P_order * (old_P_order + 1) * (old_P_order - 1)];
7711  this->node_pt(this->p_order() * this->p_order() * this->p_order() - 1) =
7712  old_node_pt[old_P_order * old_P_order * old_P_order - 1];
7713 
7714  // Copy midpoint nodes from sons if new p-order is odd
7715  if (this->p_order() % 2 == 1)
7716  {
7717  // Work out which is midpoint node
7718  unsigned n_p = this->p_order();
7719  unsigned m = (n_p - 1) / 2;
7720 
7721  unsigned s0space = m;
7722  unsigned s1space = m * n_p;
7723  unsigned s2space = m * n_p * n_p;
7724 
7725  // Back face
7726  // DB edge
7727  this->node_pt(1 * s0space + 0 * s1space + 0 * s2space) =
7728  dynamic_cast<RefineableQElement<3>*>(
7729  octree_pt()->son_pt(OcTreeNames::RDB)->object_pt())
7730  ->vertex_node_pt(0);
7731  // LB edge
7732  this->node_pt(0 * s0space + 1 * s1space + 0 * s2space) =
7733  dynamic_cast<RefineableQElement<3>*>(
7734  octree_pt()->son_pt(OcTreeNames::LUB)->object_pt())
7735  ->vertex_node_pt(0);
7736  // B centre
7737  this->node_pt(1 * s0space + 1 * s1space + 0 * s2space) =
7738  dynamic_cast<RefineableQElement<3>*>(
7739  octree_pt()->son_pt(OcTreeNames::RUB)->object_pt())
7740  ->vertex_node_pt(0);
7741  // RB edge
7742  this->node_pt(2 * s0space + 1 * s1space + 0 * s2space) =
7743  dynamic_cast<RefineableQElement<3>*>(
7744  octree_pt()->son_pt(OcTreeNames::RUB)->object_pt())
7745  ->vertex_node_pt(1);
7746  // UB edge
7747  this->node_pt(1 * s0space + 2 * s1space + 0 * s2space) =
7748  dynamic_cast<RefineableQElement<3>*>(
7749  octree_pt()->son_pt(OcTreeNames::RUB)->object_pt())
7750  ->vertex_node_pt(2);
7751 
7752  // Mid-way between between back and front faces
7753  // LD edge
7754  this->node_pt(0 * s0space + 0 * s1space + 1 * s2space) =
7755  dynamic_cast<RefineableQElement<3>*>(
7756  octree_pt()->son_pt(OcTreeNames::LDF)->object_pt())
7757  ->vertex_node_pt(0);
7758  // D centre
7759  this->node_pt(1 * s0space + 0 * s1space + 1 * s2space) =
7760  dynamic_cast<RefineableQElement<3>*>(
7761  octree_pt()->son_pt(OcTreeNames::LDF)->object_pt())
7762  ->vertex_node_pt(1);
7763  // RD edge
7764  this->node_pt(2 * s0space + 0 * s1space + 1 * s2space) =
7765  dynamic_cast<RefineableQElement<3>*>(
7766  octree_pt()->son_pt(OcTreeNames::RDF)->object_pt())
7767  ->vertex_node_pt(1);
7768  // L centre
7769  this->node_pt(0 * s0space + 1 * s1space + 1 * s2space) =
7770  dynamic_cast<RefineableQElement<3>*>(
7771  octree_pt()->son_pt(OcTreeNames::LUF)->object_pt())
7772  ->vertex_node_pt(0);
7773  // Centre
7774  this->node_pt(1 * s0space + 1 * s1space + 1 * s2space) =
7775  dynamic_cast<RefineableQElement<3>*>(
7776  octree_pt()->son_pt(OcTreeNames::RUF)->object_pt())
7777  ->vertex_node_pt(0);
7778  // R centre
7779  this->node_pt(2 * s0space + 1 * s1space + 1 * s2space) =
7780  dynamic_cast<RefineableQElement<3>*>(
7781  octree_pt()->son_pt(OcTreeNames::RUF)->object_pt())
7782  ->vertex_node_pt(1);
7783  // LU edge
7784  this->node_pt(0 * s0space + 2 * s1space + 1 * s2space) =
7785  dynamic_cast<RefineableQElement<3>*>(
7786  octree_pt()->son_pt(OcTreeNames::LUF)->object_pt())
7787  ->vertex_node_pt(2);
7788  // U center
7789  this->node_pt(1 * s0space + 2 * s1space + 1 * s2space) =
7790  dynamic_cast<RefineableQElement<3>*>(
7791  octree_pt()->son_pt(OcTreeNames::RUF)->object_pt())
7792  ->vertex_node_pt(2);
7793  // RU edge
7794  this->node_pt(2 * s0space + 2 * s1space + 1 * s2space) =
7795  dynamic_cast<RefineableQElement<3>*>(
7796  octree_pt()->son_pt(OcTreeNames::RUF)->object_pt())
7797  ->vertex_node_pt(3);
7798 
7799  // Front face
7800  // DF edge
7801  this->node_pt(1 * s0space + 0 * s1space + 2 * s2space) =
7802  dynamic_cast<RefineableQElement<3>*>(
7803  octree_pt()->son_pt(OcTreeNames::LDF)->object_pt())
7804  ->vertex_node_pt(5);
7805  // LF edge
7806  this->node_pt(0 * s0space + 1 * s1space + 2 * s2space) =
7807  dynamic_cast<RefineableQElement<3>*>(
7808  octree_pt()->son_pt(OcTreeNames::LUF)->object_pt())
7809  ->vertex_node_pt(4);
7810  // F centre
7811  this->node_pt(1 * s0space + 1 * s1space + 2 * s2space) =
7812  dynamic_cast<RefineableQElement<3>*>(
7813  octree_pt()->son_pt(OcTreeNames::RUF)->object_pt())
7814  ->vertex_node_pt(4);
7815  // RF edge
7816  this->node_pt(2 * s0space + 1 * s1space + 2 * s2space) =
7817  dynamic_cast<RefineableQElement<3>*>(
7818  octree_pt()->son_pt(OcTreeNames::RUF)->object_pt())
7819  ->vertex_node_pt(5);
7820  // UF edge
7821  this->node_pt(1 * s0space + 2 * s1space + 2 * s2space) =
7822  dynamic_cast<RefineableQElement<3>*>(
7823  octree_pt()->son_pt(OcTreeNames::RUF)->object_pt())
7824  ->vertex_node_pt(6);
7825  }
7826 
7827  // The timestepper should be the same for all nodes and node 0 should
7828  // never be deleted.
7829  if (this->node_pt(0) == 0)
7830  {
7831  throw OomphLibError(
7832  "The Corner node (0) does not exist",
7833  "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
7834  OOMPH_EXCEPTION_LOCATION);
7835  }
7836 
7837  TimeStepper* time_stepper_pt = this->node_pt(0)->time_stepper_pt();
7838  unsigned ntstorage = time_stepper_pt->ntstorage();
7839 
7840  unsigned jnod = 0;
7841  Vector<double> s_fraction(3), s(3);
7842  // Loop over the nodes in the element
7843  unsigned n_p = this->nnode_1d();
7844  for (unsigned i0 = 0; i0 < n_p; i0++)
7845  {
7846  // Get the fractional position of the node
7847  s_fraction[0] = this->local_one_d_fraction_of_node(i0, 0);
7848  // Local coordinate
7849  s[0] = -1.0 + 2.0 * s_fraction[0];
7850 
7851  for (unsigned i1 = 0; i1 < n_p; i1++)
7852  {
7853  // Get the fractional position of the node in the direction of s[1]
7854  s_fraction[1] = this->local_one_d_fraction_of_node(i1, 1);
7855  // Local coordinate in father element
7856  s[1] = -1.0 + 2.0 * s_fraction[1];
7857 
7858  for (unsigned i2 = 0; i2 < n_p; i2++)
7859  {
7860  // Get the fractional position of the node in the direction of s[2]
7861  s_fraction[2] = this->local_one_d_fraction_of_node(i2, 2);
7862  // Local coordinate in father element
7863  s[2] = -1.0 + 2.0 * s_fraction[2];
7864 
7865  // Set the local node number
7866  jnod = i0 + n_p * i1 + n_p * n_p * i2;
7867 
7868  // Initialise flag: So far, this node hasn't been built
7869  // or copied yet
7870  bool node_done = false;
7871 
7872  // Get the pointer to the node in this element, returns NULL
7873  // if there is not node
7874  Node* created_node_pt = this->get_node_at_local_coordinate(s);
7875 
7876  // Does this node already exist in this element?
7877  //----------------------------------------------
7878  if (created_node_pt != 0)
7879  {
7880  // Copy node across
7881  this->node_pt(jnod) = created_node_pt;
7882 
7883  // Node has been created by copy
7884  node_done = true;
7885  }
7886  // Node does not exist in this element but might already
7887  //------------------------------------------------------
7888  // have been created by neighbouring elements
7889  //-------------------------------------------
7890  else
7891  {
7892  // Was the node created by one of its neighbours
7893  // Whether or not the node lies on an edge can be calculated
7894  // by from the fractional position
7895  bool is_periodic = false;
7896  created_node_pt =
7897  node_created_by_neighbour(s_fraction, is_periodic);
7898 
7899  // If the node was so created, assign the pointers
7900  if (created_node_pt != 0)
7901  {
7902  // If the node is periodic
7903  if (is_periodic)
7904  {
7905  throw OomphLibError("Cannot handle periodic nodes yet",
7906  OOMPH_CURRENT_FUNCTION,
7907  OOMPH_EXCEPTION_LOCATION);
7908  }
7909  // Non-periodic case, just set the pointer
7910  else
7911  {
7912  this->node_pt(jnod) = created_node_pt;
7913  }
7914  // Node has been created
7915  node_done = true;
7916  }
7917  } // Node does not exist in this element
7918 
7919  // Node has not been built anywhere ---> build it here
7920  if (!node_done)
7921  {
7922  // First, find the son element in which the node should live
7923 
7924  // Find coordinates in the sons
7925  Vector<double> s_in_son(3);
7926  // using namespace OcTreeNames;
7927  int son = -10;
7928  // On the left
7929  if (s_fraction[0] < 0.5)
7930  {
7931  // On the bottom
7932  if (s_fraction[1] < 0.5)
7933  {
7934  // On the back
7935  if (s_fraction[2] < 0.5)
7936  {
7937  // It's the LDB son
7938  son = OcTreeNames::LDB;
7939  s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
7940  s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
7941  s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
7942  }
7943  // On the front
7944  else
7945  {
7946  // It's the LDF son
7947  son = OcTreeNames::LDF;
7948  s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
7949  s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
7950  s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
7951  }
7952  }
7953  // On the top
7954  else
7955  {
7956  // On the back
7957  if (s[2] < 0.0)
7958  {
7959  // It's the LUB son
7960  son = OcTreeNames::LUB;
7961  s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
7962  s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
7963  s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
7964  }
7965  // On the front
7966  else
7967  {
7968  // It's the LUF son
7969  son = OcTreeNames::LUF;
7970  s_in_son[0] = -1.0 + 4.0 * s_fraction[0];
7971  s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
7972  s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
7973  }
7974  }
7975  }
7976  // On the right
7977  else
7978  {
7979  // On the bottom
7980  if (s[1] < 0.0)
7981  {
7982  // On the back
7983  if (s[2] < 0.0)
7984  {
7985  // It's the RDB son
7986  son = OcTreeNames::RDB;
7987  s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
7988  s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
7989  s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
7990  }
7991  // On the front
7992  else
7993  {
7994  // It's the RDF son
7995  son = OcTreeNames::RDF;
7996  s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
7997  s_in_son[1] = -1.0 + 4.0 * s_fraction[1];
7998  s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
7999  }
8000  }
8001  // On the top
8002  else
8003  {
8004  // On the back
8005  if (s[2] < 0.0)
8006  {
8007  // It's the RUB son
8008  son = OcTreeNames::RUB;
8009  s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
8010  s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
8011  s_in_son[2] = -1.0 + 4.0 * s_fraction[2];
8012  }
8013  // On the front
8014  else
8015  {
8016  // It's the RUF son
8017  son = OcTreeNames::RUF;
8018  s_in_son[0] = -1.0 + 4.0 * (s_fraction[0] - 0.5);
8019  s_in_son[1] = -1.0 + 4.0 * (s_fraction[1] - 0.5);
8020  s_in_son[2] = -1.0 + 4.0 * (s_fraction[2] - 0.5);
8021  }
8022  }
8023  }
8024 
8025  // Get the pointer to the son element in which the new node
8026  // would live
8027  PRefineableQElement<3, INITIAL_NNODE_1D>* son_el_pt =
8028  dynamic_cast<PRefineableQElement<3, INITIAL_NNODE_1D>*>(
8029  this->tree_pt()->son_pt(son)->object_pt());
8030 
8031  // If we are rebuilding, then worry about the boundary conditions
8032  // Find the boundary of the node
8033  // Initially none
8034  int boundary = Tree::OMEGA;
8035  // If we are on the left face
8036  if (i0 == 0)
8037  {
8038  boundary = OcTreeNames::L;
8039  }
8040  // If we are on the right face
8041  else if (i0 == n_p - 1)
8042  {
8043  boundary = OcTreeNames::R;
8044  }
8045 
8046  // If we are on the bottom face
8047  if (i1 == 0)
8048  {
8049  // If we already have already set the boundary, we're on an edge
8050  switch (boundary)
8051  {
8052  case OcTreeNames::L:
8053  boundary = OcTreeNames::LD;
8054  break;
8055  case OcTreeNames::R:
8056  boundary = OcTreeNames::RD;
8057  break;
8058  // Boundary not set
8059  default:
8060  boundary = OcTreeNames::D;
8061  break;
8062  }
8063  }
8064  // If we are the top face
8065  else if (i1 == n_p - 1)
8066  {
8067  // If we already have a boundary
8068  switch (boundary)
8069  {
8070  case OcTreeNames::L:
8071  boundary = OcTreeNames::LU;
8072  break;
8073  case OcTreeNames::R:
8074  boundary = OcTreeNames::RU;
8075  break;
8076  default:
8077  boundary = OcTreeNames::U;
8078  break;
8079  }
8080  }
8081 
8082  // If we are on the back face
8083  if (i2 == 0)
8084  {
8085  // If we already have already set the boundary, we're on an edge
8086  switch (boundary)
8087  {
8088  case OcTreeNames::L:
8089  boundary = OcTreeNames::LB;
8090  break;
8091  case OcTreeNames::LD:
8092  boundary = OcTreeNames::LDB;
8093  break;
8094  case OcTreeNames::LU:
8095  boundary = OcTreeNames::LUB;
8096  break;
8097  case OcTreeNames::R:
8098  boundary = OcTreeNames::RB;
8099  break;
8100  case OcTreeNames::RD:
8101  boundary = OcTreeNames::RDB;
8102  break;
8103  case OcTreeNames::RU:
8104  boundary = OcTreeNames::RUB;
8105  break;
8106  case OcTreeNames::D:
8107  boundary = OcTreeNames::DB;
8108  break;
8109  case OcTreeNames::U:
8110  boundary = OcTreeNames::UB;
8111  break;
8112  // Boundary not set
8113  default:
8114  boundary = OcTreeNames::B;
8115  break;
8116  }
8117  }
8118  // If we are the front face
8119  else if (i2 == n_p - 1)
8120  {
8121  // If we already have a boundary
8122  switch (boundary)
8123  {
8124  case OcTreeNames::L:
8125  boundary = OcTreeNames::LF;
8126  break;
8127  case OcTreeNames::LD:
8128  boundary = OcTreeNames::LDF;
8129  break;
8130  case OcTreeNames::LU:
8131  boundary = OcTreeNames::LUF;
8132  break;
8133  case OcTreeNames::R:
8134  boundary = OcTreeNames::RF;
8135  break;
8136  case OcTreeNames::RD:
8137  boundary = OcTreeNames::RDF;
8138  break;
8139  case OcTreeNames::RU:
8140  boundary = OcTreeNames::RUF;
8141  break;
8142  case OcTreeNames::D:
8143  boundary = OcTreeNames::DF;
8144  break;
8145  case OcTreeNames::U:
8146  boundary = OcTreeNames::UF;
8147  break;
8148  default:
8149  boundary = OcTreeNames::F;
8150  break;
8151  }
8152  }
8153 
8154  // set of boundaries that this edge in the son lives on
8155  std::set<unsigned> boundaries;
8156 
8157  // Now get the boundary conditions from the son
8158  // The boundaries will be common to the son because there can be
8159  // no rotations here
8160  if (boundary != Tree::OMEGA)
8161  {
8162  son_el_pt->get_boundaries(boundary, boundaries);
8163  }
8164 
8165  // If the node lives on a boundary:
8166  // Construct a boundary node,
8167  // Get boundary conditions and
8168  // update all lookup schemes
8169  if (boundaries.size() > 0)
8170  {
8171  // Construct the new node
8172  created_node_pt =
8173  this->construct_boundary_node(jnod, time_stepper_pt);
8174 
8175  // Get the boundary conditions from the son
8176  Vector<int> bound_cons(this->ncont_interpolated_values());
8177  son_el_pt->get_bcs(boundary, bound_cons);
8178 
8179  // Loop over the values and pin if necessary
8180  unsigned nval = created_node_pt->nvalue();
8181  for (unsigned k = 0; k < nval; k++)
8182  {
8183  if (bound_cons[k])
8184  {
8185  created_node_pt->pin(k);
8186  }
8187  }
8188 
8189  // Solid node? If so, deal with the positional boundary
8190  // conditions:
8191  SolidNode* solid_node_pt =
8192  dynamic_cast<SolidNode*>(created_node_pt);
8193  if (solid_node_pt != 0)
8194  {
8195  // Get the positional boundary conditions from the father:
8196  unsigned n_dim = created_node_pt->ndim();
8197  Vector<int> solid_bound_cons(n_dim);
8198  RefineableSolidQElement<3>* son_solid_el_pt =
8199  dynamic_cast<RefineableSolidQElement<3>*>(son_el_pt);
8200 #ifdef PARANOID
8201  if (son_solid_el_pt == 0)
8202  {
8203  std::string error_message =
8204  "We have a SolidNode outside a refineable SolidElement\n";
8205  error_message +=
8206  "during mesh refinement -- this doesn't make sense\n";
8207 
8208  throw OomphLibError(error_message,
8209  "PRefineableQElement<3,INITIAL_NNODE_1D>:"
8210  ":rebuild_from_sons()",
8211  OOMPH_EXCEPTION_LOCATION);
8212  }
8213 #endif
8214  son_solid_el_pt->get_solid_bcs(boundary, solid_bound_cons);
8215 
8216  // Loop over the positions and pin, if necessary
8217  for (unsigned k = 0; k < n_dim; k++)
8218  {
8219  if (solid_bound_cons[k])
8220  {
8221  solid_node_pt->pin_position(k);
8222  }
8223  }
8224  } // End of if solid_node_pt
8225 
8226 
8227  // Next we update the boundary look-up schemes
8228  // Loop over the boundaries stored in the set
8229  for (std::set<unsigned>::iterator it = boundaries.begin();
8230  it != boundaries.end();
8231  ++it)
8232  {
8233  // Add the node to the boundary
8234  mesh_pt->add_boundary_node(*it, created_node_pt);
8235 
8236  // If we have set an intrinsic coordinate on this
8237  // mesh boundary then it must also be interpolated on
8238  // the new node
8239  // Now interpolate the intrinsic boundary coordinate
8240  if (mesh_pt->boundary_coordinate_exists(*it) == true)
8241  {
8242  Vector<double> zeta(2);
8243  son_el_pt->interpolated_zeta_on_face(
8244  *it, boundary, s_in_son, zeta);
8245 
8246  created_node_pt->set_coordinates_on_boundary(*it, zeta);
8247  }
8248  }
8249  }
8250  // Otherwise the node is not on a Mesh boundary
8251  // and we create a normal "bulk" node
8252  else
8253  {
8254  // Construct the new node
8255  created_node_pt = this->construct_node(jnod, time_stepper_pt);
8256  }
8257 
8258  // Now we set the position and values at the newly created node
8259 
8260  // In the first instance use macro element or FE representation
8261  // to create past and present nodal positions.
8262  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
8263  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
8264  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
8265  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
8266  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
8267  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
8268  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
8269 
8270  // Loop over # of history values
8271  // Loop over # of history values
8272  for (unsigned t = 0; t < ntstorage; t++)
8273  {
8274  using namespace QuadTreeNames;
8275  // Get the position from the son
8276  Vector<double> x_prev(3);
8277 
8278  // Now let's fill in the value
8279  son_el_pt->get_x(t, s_in_son, x_prev);
8280  for (unsigned i = 0; i < 3; i++)
8281  {
8282  created_node_pt->x(t, i) = x_prev[i];
8283  }
8284  }
8285 
8286  // Now set up the values
8287  // Loop over all history values
8288  for (unsigned t = 0; t < ntstorage; t++)
8289  {
8290  // Get values from father element
8291  // Note: get_interpolated_values() sets Vector size itself.
8292  Vector<double> prev_values;
8293  son_el_pt->get_interpolated_values(t, s_in_son, prev_values);
8294 
8295  // Initialise the values at the new node
8296  for (unsigned k = 0; k < created_node_pt->nvalue(); k++)
8297  {
8298  created_node_pt->set_value(t, k, prev_values[k]);
8299  }
8300  }
8301 
8302  // Add the node to the mesh
8303  mesh_pt->add_node_pt(created_node_pt);
8304 
8305  // Check if the element is an algebraic element
8306  AlgebraicElementBase* alg_el_pt =
8307  dynamic_cast<AlgebraicElementBase*>(this);
8308 
8309  // If we do have an algebraic element
8310  if (alg_el_pt != 0)
8311  {
8312  std::string error_message =
8313  "Have not implemented rebuilding from sons for";
8314  error_message += "Algebraic p-refineable elements yet\n";
8315 
8316  throw OomphLibError(
8317  error_message,
8318  "PRefineableQElement<3,INITIAL_NNODE_1D>::rebuild_from_sons()",
8319  OOMPH_EXCEPTION_LOCATION);
8320  }
8321 
8322  } // End of the case when we build the node ourselves
8323  }
8324  }
8325  }
8326  }
oomph::BrickElementBase::vertex_node_pt
virtual Node * vertex_node_pt(const unsigned &j) const =0
Pointer to the j-th vertex node in the element.
oomph::PRefineableQElement< 3, 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:4891
oomph::RefineableElement::tree_pt
Tree * tree_pt()
Access function: Pointer to quadtree representation of this element.
Definition: refineable_elements.h:211
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::OcTreeNames::B, oomph::Mesh::boundary_coordinate_exists(), oomph::OcTreeNames::D, oomph::OcTreeNames::DB, oomph::OcTreeNames::DF, oomph::OcTreeNames::F, oomph::RefineableQElement< 3 >::get_bcs(), oomph::RefineableQElement< 3 >::get_boundaries(), oomph::RefineableElement::get_interpolated_values(), oomph::RefineableSolidQElement< 3 >::get_solid_bcs(), oomph::FiniteElement::get_x(), i, oomph::RefineableQElement< 3 >::interpolated_zeta_on_face(), k, oomph::OcTreeNames::L, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::OcTreeNames::LDB, oomph::OcTreeNames::LDF, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, oomph::OcTreeNames::LUB, oomph::OcTreeNames::LUF, m, n, oomph::Node::ndim(), 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(), oomph::OcTreeNames::R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, oomph::OcTreeNames::RDB, oomph::OcTreeNames::RDF, oomph::OcTreeNames::RF, oomph::OcTreeNames::RU, oomph::OcTreeNames::RUB, oomph::OcTreeNames::RUF, s, oomph::Node::set_coordinates_on_boundary(), oomph::Data::set_value(), oomph::Global_string_for_annotation::string(), plotPSD::t, oomph::OcTreeNames::U, oomph::OcTreeNames::UB, oomph::OcTreeNames::UF, oomph::Node::x(), and Eigen::zeta().

◆ shape()

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

Overload the shape functions.

Shape functions for PRefineableQElement<DIM>

Implements oomph::FiniteElement.

7281  {
7282  switch (P_order)
7283  {
7284  case 2:
7285  {
7286  // Call the OneDimensional Shape functions
7287  OneDimensionalLegendreShape<2>::calculate_nodal_positions();
7288  OneDimensionalLegendreShape<2> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7289 
7290  // Now let's loop over the nodal points in the element
7291  // and copy the values back in
7292  for (unsigned i = 0; i < P_order; i++)
7293  {
7294  for (unsigned j = 0; j < P_order; j++)
7295  {
7296  for (unsigned k = 0; k < P_order; k++)
7297  {
7298  psi(P_order * P_order * i + P_order * j + k) =
7299  psi3[i] * psi2[j] * psi1[k];
7300  }
7301  }
7302  }
7303  break;
7304  }
7305  case 3:
7306  {
7307  // Call the OneDimensional Shape functions
7308  OneDimensionalLegendreShape<3>::calculate_nodal_positions();
7309  OneDimensionalLegendreShape<3> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7310 
7311  // Now let's loop over the nodal points in the element
7312  // and copy the values back in
7313  for (unsigned i = 0; i < P_order; i++)
7314  {
7315  for (unsigned j = 0; j < P_order; j++)
7316  {
7317  for (unsigned k = 0; k < P_order; k++)
7318  {
7319  psi(P_order * P_order * i + P_order * j + k) =
7320  psi3[i] * psi2[j] * psi1[k];
7321  }
7322  }
7323  }
7324  break;
7325  }
7326  case 4:
7327  {
7328  // Call the OneDimensional Shape functions
7329  OneDimensionalLegendreShape<4>::calculate_nodal_positions();
7330  OneDimensionalLegendreShape<4> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7331 
7332  // Now let's loop over the nodal points in the element
7333  // and copy the values back in
7334  for (unsigned i = 0; i < P_order; i++)
7335  {
7336  for (unsigned j = 0; j < P_order; j++)
7337  {
7338  for (unsigned k = 0; k < P_order; k++)
7339  {
7340  psi(P_order * P_order * i + P_order * j + k) =
7341  psi3[i] * psi2[j] * psi1[k];
7342  }
7343  }
7344  }
7345  break;
7346  }
7347  case 5:
7348  {
7349  // Call the OneDimensional Shape functions
7350  OneDimensionalLegendreShape<5>::calculate_nodal_positions();
7351  OneDimensionalLegendreShape<5> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7352 
7353  // Now let's loop over the nodal points in the element
7354  // and copy the values back in
7355  for (unsigned i = 0; i < P_order; i++)
7356  {
7357  for (unsigned j = 0; j < P_order; j++)
7358  {
7359  for (unsigned k = 0; k < P_order; k++)
7360  {
7361  psi(P_order * P_order * i + P_order * j + k) =
7362  psi3[i] * psi2[j] * psi1[k];
7363  }
7364  }
7365  }
7366  break;
7367  }
7368  case 6:
7369  {
7370  // Call the OneDimensional Shape functions
7371  OneDimensionalLegendreShape<6>::calculate_nodal_positions();
7372  OneDimensionalLegendreShape<6> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7373 
7374  // Now let's loop over the nodal points in the element
7375  // and copy the values back in
7376  for (unsigned i = 0; i < P_order; i++)
7377  {
7378  for (unsigned j = 0; j < P_order; j++)
7379  {
7380  for (unsigned k = 0; k < P_order; k++)
7381  {
7382  psi(P_order * P_order * i + P_order * j + k) =
7383  psi3[i] * psi2[j] * psi1[k];
7384  }
7385  }
7386  }
7387  break;
7388  }
7389  case 7:
7390  {
7391  // Call the OneDimensional Shape functions
7392  OneDimensionalLegendreShape<7>::calculate_nodal_positions();
7393  OneDimensionalLegendreShape<7> psi1(s[0]), psi2(s[1]), psi3(s[2]);
7394 
7395  // Now let's loop over the nodal points in the element
7396  // and copy the values back in
7397  for (unsigned i = 0; i < P_order; i++)
7398  {
7399  for (unsigned j = 0; j < P_order; j++)
7400  {
7401  for (unsigned k = 0; k < P_order; k++)
7402  {
7403  psi(P_order * P_order * i + P_order * j + k) =
7404  psi3[i] * psi2[j] * psi1[k];
7405  }
7406  }
7407  }
7408  break;
7409  }
7410  default:
7411  std::ostringstream error_message;
7412  error_message << "\nERROR: Exceeded maximum polynomial order for";
7413  error_message << "\n polynomial order for shape functions.\n";
7414  throw OomphLibError(error_message.str(),
7415  OOMPH_CURRENT_FUNCTION,
7416  OOMPH_EXCEPTION_LOCATION);
7417  }
7418  }

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

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