Public Types | Public Member Functions | Protected Member Functions | Static Protected Attributes | List of all members
oomph::RefineableQElement< 3 > Class Referenceabstract

#include <refineable_brick_element.h>

+ Inheritance diagram for oomph::RefineableQElement< 3 >:

Public Types

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 > &)
 

Public Member Functions

 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 Node * node_created_by_neighbour (const Vector< double > &s_fraction, bool &is_periodic)
 
virtual Node * node_created_by_son_of_neighbour (const Vector< double > &s_fraction, bool &is_periodic)
 
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)
 
virtual void further_setup_hanging_nodes ()=0
 
- 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 rebuild_from_sons (Mesh *&mesh_pt)=0
 
virtual void unbuild ()
 
virtual void deactivate_element ()
 
virtual bool nodes_built ()
 Return true if all the nodes have been built, false if not. More...
 
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 initial_setup (Tree *const &adopted_father_pt=0, const unsigned &initial_p_order=0)
 
virtual void pre_build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt)
 Pre-build the element. More...
 
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
 
virtual void local_coordinate_of_node (const unsigned &j, Vector< double > &s) const
 
virtual void local_fraction_of_node (const unsigned &j, Vector< double > &s_fraction)
 
virtual double local_one_d_fraction_of_node (const unsigned &n1d, const unsigned &i)
 
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 (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) 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...
 
virtual unsigned nnode_1d () const
 
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
 
virtual Node * get_node_at_local_coordinate (const Vector< double > &s) 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...
 

Protected Member Functions

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)
 
virtual void oc_hang_helper (const int &value_id, const int &my_edge, std::ofstream &output_hangfile)
 
- 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)
 

Static Protected Attributes

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
 

Additional Inherited Members

- 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
 

Detailed Description

Refineable version of QElement<3,NNODE_1D>.

Refinement is performed by octree procedures. When the element is subdivided, the geometry of its sons is established by calls to their father's

get_x(...)
oomph::FiniteElement::get_x
void get_x(const Vector< double > &s, Vector< double > &x) const
Definition: elements.h:1885

function which refers to

  • the father element's geometric FE mapping (this is the default)

or

  • to a MacroElement 's MacroElement::macro_map (if the pointer to the macro element is non-NULL)

The class provides a generic RefineableQElement<3>::build() function which deals with generic isoparametric QElements in which all values are associated with nodes. The RefineableQElement<3>::further_build() function provides an interface for any element-specific non-generic build operations.

Member Typedef Documentation

◆ VoidMemFctPt

typedef void(RefineableQElement<3>::* oomph::RefineableQElement< 3 >::VoidMemFctPt) ()

Shorthand for pointer to an argument-free void member function of the refineable element

Constructor & Destructor Documentation

◆ RefineableQElement() [1/2]

oomph::RefineableQElement< 3 >::RefineableQElement ( )
inline

Constructor: Pass refinement level (default 0 = root)

75  : RefineableElement()
76  {
77 #ifdef LEAK_CHECK
78  LeakCheckNames::RefineableQElement<3> _build += 1;
79 #endif
80  }
oomph::RefineableElement::RefineableElement
RefineableElement()
Definition: refineable_elements.h:188

◆ RefineableQElement() [2/2]

oomph::RefineableQElement< 3 >::RefineableQElement ( const RefineableQElement< 3 > &  dummy)
delete

Broken copy constructor.

◆ ~RefineableQElement()

virtual oomph::RefineableQElement< 3 >::~RefineableQElement ( )
inlinevirtual

Broken assignment operator.

Destructor

95  {
96 #ifdef LEAK_CHECK
97  LeakCheckNames::RefineableQElement<3> _build -= 1;
98 #endif
99  }

Member Function Documentation

◆ build()

void oomph::RefineableQElement< 3 >::build ( Mesh *&  mesh_pt,
Vector< Node * > &  new_node_pt,
bool &  was_already_built,
std::ofstream &  new_nodes_file 
)
virtual

Build the element: i.e. give it nodal positions, apply BCs, etc. Pointers to any new nodes will be returned in new_node_pt. If it is open, the positions of the new nodes will be written to the file stream new_nodes_file

Build the element by doing the following:

  • Give it nodal positions (by establishing the pointers to its nodes)
  • In the process create new nodes where required (i.e. if they don't exist in father element or have already been created while building new neighbour elements). Node building involves the following steps:
    • Get nodal position from father element.
    • Establish the time-history of the newly created nodal point (its coordinates and the previous values) consistent with the father's history.
    • Determine the boundary conditions of the nodes (newly created nodes can only lie on the interior of any edges of the father element – this makes it possible to to figure out what their bc should be...)
    • Add node to the mesh's stoarge scheme for the boundary nodes.
    • Add the new node to the mesh itself
    • Doc newly created nodes in file "new_nodes.dat" in the directory / of the DocInfo object (only if it's open!)
  • Finally, excute the element-specific further_build() (empty by default – must be overloaded for specific elements). This deals with any build operations that are not included in the generic process outlined above. For instance, in Crouzeix Raviart elements we need to initialise the internal pressure values in manner consistent with the pressure distribution in the father element.

Implements oomph::RefineableElement.

Reimplemented in oomph::RefineableSolidQElement< 3 >.

1297  {
1298  using namespace OcTreeNames;
1299 
1300  // Number of dimensions
1301  unsigned n_dim = 3;
1302 
1303  // Get the number of 1d nodes
1304  unsigned n_p = nnode_1d();
1305 
1306  // Check whether static father_bound needs to be created
1307  if (Father_bound[n_p].nrow() == 0)
1308  {
1309  setup_father_bounds();
1310  }
1311 
1312  // Pointer to my father (in octree impersonation)
1313  OcTree* father_pt = dynamic_cast<OcTree*>(octree_pt()->father_pt());
1314 
1315  // What type of son am I? Ask my octree representation...
1316  int son_type = octree_pt()->son_type();
1317 
1318  // Has somebody build me already? (If any nodes have not been built)
1319  if (!nodes_built())
1320  {
1321 #ifdef PARANOID
1322  if (father_pt == 0)
1323  {
1324  std::string error_message =
1325  "Something fishy here: I have no father and yet \n";
1326  error_message += "I have no nodes. Who has created me then?!\n";
1327 
1328  throw OomphLibError(
1329  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1330  }
1331 #endif
1332 
1333  // Indicate status:
1334  was_already_built = false;
1335 
1336  // Return pointer to father element
1337  RefineableQElement<3>* father_el_pt =
1338  dynamic_cast<RefineableQElement<3>*>(father_pt->object_pt());
1339 
1340  // Timestepper should be the same for all nodes in father
1341  // element -- use it create timesteppers for new nodes
1342  TimeStepper* time_stepper_pt =
1343  father_el_pt->node_pt(0)->time_stepper_pt();
1344 
1345  // Number of history values (incl. present)
1346  unsigned ntstorage = time_stepper_pt->ntstorage();
1347 
1348  // Currently we can't handle the case of generalised coordinates
1349  // since we haven't established how they should be interpolated
1350  // Buffer this case:
1351  if (father_el_pt->node_pt(0)->nposition_type() != 1)
1352  {
1353  throw OomphLibError("Can't handle generalised nodal positions (yet).",
1354  OOMPH_CURRENT_FUNCTION,
1355  OOMPH_EXCEPTION_LOCATION);
1356  }
1357 
1358  Vector<int> s_lo(n_dim);
1359  Vector<int> s_hi(n_dim);
1360  Vector<double> s(n_dim);
1361  Vector<double> x(n_dim);
1362 
1363  // Setup vertex coordinates in father element:
1364  //--------------------------------------------
1365  // find the s_lo coordinates
1366  s_lo = octree_pt()->Direction_to_vector[son_type];
1367 
1368  // Just scale them, because the Direction_to_vector
1369  // doesn't really gives s_lo;
1370  for (unsigned i = 0; i < n_dim; i++)
1371  {
1372  s_lo[i] = (s_lo[i] + 1) / 2 - 1;
1373  }
1374 
1375  // setup s_hi (Actually s_hi[i]=s_lo[i]+1)
1376  for (unsigned i = 0; i < n_dim; i++)
1377  {
1378  s_hi[i] = s_lo[i] + 1;
1379  }
1380 
1381  // Pass macro element pointer on to sons and
1382  // set coordinates in macro element
1383  if (father_el_pt->macro_elem_pt() != 0)
1384  {
1385  set_macro_elem_pt(father_el_pt->macro_elem_pt());
1386  for (unsigned i = 0; i < n_dim; i++)
1387  {
1388  s_macro_ll(i) =
1389  father_el_pt->s_macro_ll(i) +
1390  0.5 * (s_lo[i] + 1.0) *
1391  (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1392  s_macro_ur(i) =
1393  father_el_pt->s_macro_ll(i) +
1394  0.5 * (s_hi[i] + 1.0) *
1395  (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1396  }
1397  }
1398 
1399 
1400  // If the father element hasn't been generated yet, we're stuck...
1401  if (father_el_pt->node_pt(0) == 0)
1402  {
1403  throw OomphLibError(
1404  "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
1405  OOMPH_CURRENT_FUNCTION,
1406  OOMPH_EXCEPTION_LOCATION);
1407  }
1408  else
1409  {
1410  unsigned jnod = 0;
1411  Vector<double> s_fraction(n_dim);
1412 
1413  // Loop over nodes in element
1414  for (unsigned i0 = 0; i0 < n_p; i0++)
1415  {
1416  // Get the fractional position of the node in the direction of s[0]
1417  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
1418  // Local coordinate in father element
1419  s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
1420 
1421  for (unsigned i1 = 0; i1 < n_p; i1++)
1422  {
1423  // Get the fractional position of the node in the direction of s[1]
1424  s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
1425  // Local coordinate in father element
1426  s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
1427 
1428  for (unsigned i2 = 0; i2 < n_p; i2++)
1429  {
1430  // Get the fractional position of the node in the direction of
1431  // s[2]
1432  s_fraction[2] = local_one_d_fraction_of_node(i2, 2);
1433 
1434  // Local coordinate in father element
1435  s[2] = s_lo[2] + (s_hi[2] - s_lo[2]) * s_fraction[2];
1436 
1437  // Local node number
1438  jnod = i0 + n_p * i1 + n_p * n_p * i2;
1439 
1440  // Initialise flag: So far, this node hasn't been built
1441  // or copied yet
1442  bool node_done = false;
1443 
1444  // Get the pointer to the node in the father; returns NULL
1445  // if there is not a node
1446  Node* created_node_pt =
1447  father_el_pt->get_node_at_local_coordinate(s);
1448 
1449  // Does this node already exist in father element?
1450  //------------------------------------------------
1451  bool node_exists_in_father = false;
1452  if (created_node_pt != 0)
1453  {
1454  // Remember this!
1455  node_exists_in_father = true;
1456 
1457  // Copy node across
1458  node_pt(jnod) = created_node_pt;
1459 
1460  // Make sure that we update the values at the node so that
1461  // they are consistent with the present representation.
1462  // This is only need for mixed interpolation where the value
1463  // at the father could now become active.
1464 
1465  // Loop over all history values
1466  for (unsigned t = 0; t < ntstorage; t++)
1467  {
1468  // Get values from father element
1469  // Note: get_interpolated_values() sets Vector size itself.
1470  Vector<double> prev_values;
1471  father_el_pt->get_interpolated_values(t, s, prev_values);
1472  // Find the minimum number of values
1473  //(either those stored at the node, or those returned by
1474  // the function)
1475  unsigned n_val_at_node = created_node_pt->nvalue();
1476  unsigned n_val_from_function = prev_values.size();
1477  // Use the ternary conditional operator here
1478  unsigned n_var = n_val_at_node < n_val_from_function ?
1479  n_val_at_node :
1480  n_val_from_function;
1481  // Assign the values that we can
1482  for (unsigned k = 0; k < n_var; k++)
1483  {
1484  created_node_pt->set_value(t, k, prev_values[k]);
1485  }
1486  }
1487 
1488  // Node has been created by copy
1489  node_done = true;
1490  }
1491  // Node does not exist in father element but might already
1492  //--------------------------------------------------------
1493  // have been created by neighbouring elements
1494  //-------------------------------------------
1495  else
1496  {
1497  // Boolean to check if the node is periodic
1498  bool is_periodic = false;
1499 
1500  // Was the node created by one of its neighbours
1501  // Whether or not the node lies on an edge can be determined
1502  // from the fractional position
1503  created_node_pt =
1504  node_created_by_neighbour(s_fraction, is_periodic);
1505 
1506  // If so, then copy the pointer across
1507  if (created_node_pt != 0)
1508  {
1509  // Now the node must be on a boundary, but we don't know which
1510  // one. The returned created_node_pt is actually the
1511  // neighbouring periodic node
1512  Node* neighbour_node_pt = created_node_pt;
1513 
1514  // Determine the edge on which the new node will live
1515  int father_bound = Father_bound[n_p](jnod, son_type);
1516 
1517  // Storage for the set of Mesh boundaries on which the
1518  // appropriate father edge lives.
1519  std::set<unsigned> boundaries;
1520 
1521  // Only get the boundaries if we are at the edge of
1522  // an element. Nodes in the centre of an element cannot be
1523  // on Mesh boundaries
1524  if (father_bound != Tree::OMEGA)
1525  {
1526  father_el_pt->get_boundaries(father_bound, boundaries);
1527  }
1528 
1529 #ifdef PARANOID
1530  // Case where a new node lives on more than one boundary
1531  // seems fishy enough to flag
1532  if (boundaries.size() > 2)
1533  {
1534  throw OomphLibError(
1535  "boundaries.size()>2 seems a bit strange..\n",
1536  OOMPH_CURRENT_FUNCTION,
1537  OOMPH_EXCEPTION_LOCATION);
1538  }
1539 #endif
1540 
1541  // If the node is periodic and is definitely a boundary node.
1542  // NOTE: the reason for this is that confusion can arise when
1543  // a node is created on an edge that joins a periodic face.
1544  if ((is_periodic) && (boundaries.size() > 0))
1545  {
1546  // Create node and set the pointer to it from the element
1547  created_node_pt =
1548  construct_boundary_node(jnod, time_stepper_pt);
1549 
1550  // Make the node periodic from the neighbour
1551  created_node_pt->make_periodic(neighbour_node_pt);
1552 
1553  // Add to vector of new nodes
1554  new_node_pt.push_back(created_node_pt);
1555 
1556  // Loop over # of history values
1557  for (unsigned t = 0; t < ntstorage; t++)
1558  {
1559  // Get position from father element -- this uses the macro
1560  // element representation if appropriate. If the node
1561  // turns out to be a hanging node later on, then
1562  // its position gets adjusted in line with its
1563  // hanging node interpolation.
1564  Vector<double> x_prev(n_dim);
1565  father_el_pt->get_x(t, s, x_prev);
1566  // Set previous positions of the new node
1567  for (unsigned i = 0; i < n_dim; i++)
1568  {
1569  created_node_pt->x(t, i) = x_prev[i];
1570  }
1571  }
1572 
1573  // Next, we Update the boundary lookup schemes
1574  // Loop over the boundaries stored in the set
1575  for (std::set<unsigned>::iterator it = boundaries.begin();
1576  it != boundaries.end();
1577  ++it)
1578  {
1579  // Add the node to the boundary
1580  mesh_pt->add_boundary_node(*it, created_node_pt);
1581 
1582  // If we have set an intrinsic coordinate on this
1583  // mesh boundary then it must also be interpolated on
1584  // the new node
1585  // Now interpolate the intrinsic boundary coordinate
1586  if (mesh_pt->boundary_coordinate_exists(*it) == true)
1587  {
1588  Vector<double> zeta(2, 0.0);
1589  father_el_pt->interpolated_zeta_on_face(
1590  *it, father_bound, s, zeta);
1591  created_node_pt->set_coordinates_on_boundary(*it, zeta);
1592  }
1593  }
1594 
1595  // Make sure that we add the node to the mesh
1596  mesh_pt->add_node_pt(created_node_pt);
1597  } // End of periodic case
1598  // Otherwise the node is not periodic, so just set the
1599  // pointer to the neighbours node
1600  else
1601  {
1602  node_pt(jnod) = created_node_pt;
1603  }
1604  node_done = true;
1605  }
1606  // Node does not exist in neighbour element but might already
1607  //-----------------------------------------------------------
1608  // have been created by a son of a neighbouring element
1609  //-----------------------------------------------------
1610  else
1611  {
1612  // Was the node created by one of its neighbours' sons
1613  // Whether or not the node lies on an edge can be calculated
1614  // by from the fractional position
1615  bool is_periodic = false;
1616  created_node_pt =
1617  node_created_by_son_of_neighbour(s_fraction, is_periodic);
1618 
1619  // If the node was so created, assign the pointers
1620  if (created_node_pt != 0)
1621  {
1622  // If the node is periodic
1623  if (is_periodic)
1624  {
1625  // Now the node must be on a boundary, but we don't know
1626  // which one The returned created_node_pt is actually the
1627  // neighbouring periodic node
1628  Node* neighbour_node_pt = created_node_pt;
1629 
1630  // Determine the edge on which the new node will live
1631  int father_bound = Father_bound[n_p](jnod, son_type);
1632 
1633  // Storage for the set of Mesh boundaries on which the
1634  // appropriate father edge lives.
1635  std::set<unsigned> boundaries;
1636 
1637  // Only get the boundaries if we are at the edge of
1638  // an element. Nodes in the centre of an element cannot be
1639  // on Mesh boundaries
1640  if (father_bound != Tree::OMEGA)
1641  {
1642  father_el_pt->get_boundaries(father_bound, boundaries);
1643  }
1644 
1645 #ifdef PARANOID
1646  // Case where a new node lives on more than one boundary
1647  // seems fishy enough to flag
1648  if (boundaries.size() > 2)
1649  {
1650  throw OomphLibError(
1651  "boundaries.size()>2 seems a bit strange..\n",
1652  OOMPH_CURRENT_FUNCTION,
1653  OOMPH_EXCEPTION_LOCATION);
1654  }
1655 
1656  // Case when there are no boundaries, we are in big
1657  // trouble
1658  if (boundaries.size() == 0)
1659  {
1660  std::ostringstream error_stream;
1661  error_stream
1662  << "Periodic node is not on a boundary...\n"
1663  << "Coordinates: " << created_node_pt->x(0) << " "
1664  << created_node_pt->x(1) << "\n";
1665  throw OomphLibError(error_stream.str(),
1666  OOMPH_CURRENT_FUNCTION,
1667  OOMPH_EXCEPTION_LOCATION);
1668  }
1669 #endif
1670 
1671  // Create node and set the pointer to it from the element
1672  created_node_pt =
1673  construct_boundary_node(jnod, time_stepper_pt);
1674 
1675  // Make the node periodic from the neighbour
1676  created_node_pt->make_periodic(neighbour_node_pt);
1677 
1678  // Add to vector of new nodes
1679  new_node_pt.push_back(created_node_pt);
1680 
1681  // Loop over # of history values
1682  for (unsigned t = 0; t < ntstorage; t++)
1683  {
1684  // Get position from father element -- this uses the
1685  // macro element representation if appropriate. If the
1686  // node turns out to be a hanging node later on, then
1687  // its position gets adjusted in line with its
1688  // hanging node interpolation.
1689  Vector<double> x_prev(n_dim, 0.0);
1690  father_el_pt->get_x(t, s, x_prev);
1691  // Set previous positions of the new node
1692  for (unsigned i = 0; i < n_dim; i++)
1693  {
1694  created_node_pt->x(t, i) = x_prev[i];
1695  }
1696  }
1697 
1698  // Next, we Update the boundary lookup schemes
1699  // Loop over the boundaries stored in the set
1700  for (std::set<unsigned>::iterator it = boundaries.begin();
1701  it != boundaries.end();
1702  ++it)
1703  {
1704  // Add the node to the boundary
1705  mesh_pt->add_boundary_node(*it, created_node_pt);
1706 
1707  // If we have set an intrinsic coordinate on this
1708  // mesh boundary then it must also be interpolated on
1709  // the new node
1710  // Now interpolate the intrinsic boundary coordinate
1711  if (mesh_pt->boundary_coordinate_exists(*it) == true)
1712  {
1713  Vector<double> zeta(2, 0.0);
1714  father_el_pt->interpolated_zeta_on_face(
1715  *it, father_bound, s, zeta);
1716  created_node_pt->set_coordinates_on_boundary(*it,
1717  zeta);
1718  }
1719  }
1720 
1721  // Make sure that we add the node to the mesh
1722  mesh_pt->add_node_pt(created_node_pt);
1723  } // End of periodic case
1724  // Otherwise the node is not periodic, so just set the
1725  // pointer to the neighbours node
1726  else
1727  {
1728  node_pt(jnod) = created_node_pt;
1729  }
1730  // Node has been created
1731  node_done = true;
1732  } // Node does not exist in son of neighbouring element
1733  } // Node does not exist in neighbouring element
1734  } // Node does not exist in father element
1735 
1736  // Node already exists in father: No need to do anything else!
1737  // otherwise deal with boundary information etc.
1738  if (!node_exists_in_father)
1739  {
1740  // Check boundary status
1741  //----------------------
1742 
1743  // Firstly, we need to determine whether or not a node lies
1744  // on the boundary before building it, because
1745  // we actually assign a different type of node on boundaries.
1746 
1747  // If the new node lives on a face that is
1748  // shared with a face of its father element,
1749  // it needs to inherit the bounday conditions
1750  // from the father face
1751  int father_bound = Father_bound[n_p](jnod, son_type);
1752 
1753  // Storage for the set of Mesh boundaries on which the
1754  // appropriate father face lives.
1755  // [New nodes should always be mid-edge nodes in father
1756  // and therefore only live on one boundary but just to
1757  // play it safe...]
1758  std::set<unsigned> boundaries;
1759 
1760  // Only get the boundaries if we are at the edge of
1761  // an element. Nodes in the centre of an element cannot be
1762  // on Mesh boundaries
1763  if (father_bound != Tree::OMEGA)
1764  {
1765  father_el_pt->get_boundaries(father_bound, boundaries);
1766  }
1767 
1768 #ifdef PARANOID
1769  // Case where a new node lives on more than two boundaries
1770  // seems fishy enough to flag
1771  if (boundaries.size() > 2)
1772  {
1773  throw OomphLibError(
1774  "boundaries.size()>2 seems a bit strange..\n",
1775  OOMPH_CURRENT_FUNCTION,
1776  OOMPH_EXCEPTION_LOCATION);
1777  }
1778 #endif
1779 
1780  // If the node lives on a mesh boundary,
1781  // then we need to create a boundary node
1782  if (boundaries.size() > 0)
1783  {
1784  // Do we need a new node?
1785  if (!node_done)
1786  {
1787  // Create node and set the internal pointer
1788  created_node_pt =
1789  construct_boundary_node(jnod, time_stepper_pt);
1790  // Add to vector of new nodes
1791  new_node_pt.push_back(created_node_pt);
1792  }
1793 
1794  // Now we need to work out whether to pin the values at
1795  // the new node based on the boundary conditions applied at
1796  // its Mesh boundary
1797 
1798  // Get the boundary conditions from the father.
1799  // Note: We can only deal with the values that are
1800  // continuously interpolated in the bulk element.
1801  unsigned n_cont = ncont_interpolated_values();
1802  Vector<int> bound_cons(n_cont);
1803  father_el_pt->get_bcs(father_bound, bound_cons);
1804 
1805  // Loop over the continuously interpolated values and pin,
1806  // if necessary
1807  for (unsigned k = 0; k < n_cont; k++)
1808  {
1809  if (bound_cons[k])
1810  {
1811  created_node_pt->pin(k);
1812  }
1813  }
1814 
1815  // Solid node? If so, deal with the positional boundary
1816  // conditions:
1817  SolidNode* solid_node_pt =
1818  dynamic_cast<SolidNode*>(created_node_pt);
1819  if (solid_node_pt != 0)
1820  {
1821  // Get the positional boundary conditions from the father:
1822  unsigned n_dim = created_node_pt->ndim();
1823  Vector<int> solid_bound_cons(n_dim);
1824  RefineableSolidQElement<3>* father_solid_el_pt =
1825  dynamic_cast<RefineableSolidQElement<3>*>(father_el_pt);
1826 #ifdef PARANOID
1827  if (father_solid_el_pt == 0)
1828  {
1829  std::string error_message = "We have a SolidNode outside "
1830  "a refineable SolidElement\n";
1831  error_message +=
1832  "during mesh refinement -- this doesn't make sense";
1833 
1834  throw OomphLibError(error_message,
1835  OOMPH_CURRENT_FUNCTION,
1836  OOMPH_EXCEPTION_LOCATION);
1837  }
1838 #endif
1839  father_solid_el_pt->get_solid_bcs(father_bound,
1840  solid_bound_cons);
1841 
1842  // Loop over the positions and pin, if necessary
1843  for (unsigned k = 0; k < n_dim; k++)
1844  {
1845  if (solid_bound_cons[k])
1846  {
1847  solid_node_pt->pin_position(k);
1848  }
1849  }
1850  } // End of if solid_node_pt
1851 
1852  // Next update the boundary look-up schemes
1853 
1854  // Loop over the boundaries in the set
1855  for (std::set<unsigned>::iterator it = boundaries.begin();
1856  it != boundaries.end();
1857  ++it)
1858  {
1859  // Add the node to the bounadry
1860  mesh_pt->add_boundary_node(*it, created_node_pt);
1861 
1862  // If we have set an intrinsic coordinate on this
1863  // mesh boundary then it must also be interpolated on
1864  // the new node
1865  // Now interpolate the intrinsic boundary coordinate
1866  if (mesh_pt->boundary_coordinate_exists(*it) == true)
1867  {
1868  // Usually there will be two coordinates
1869  Vector<double> zeta(2);
1870  father_el_pt->interpolated_zeta_on_face(
1871  *it, father_bound, s, zeta);
1872 
1873  created_node_pt->set_coordinates_on_boundary(*it, zeta);
1874  }
1875  }
1876  }
1877  // Otherwise the node is not on a Mesh boundary and
1878  // we create a normal "bulk" node
1879  else
1880  {
1881  // Do we need a new node?
1882  if (!node_done)
1883  {
1884  // Create node and set the pointer to it from the element
1885  created_node_pt = construct_node(jnod, time_stepper_pt);
1886  // Add to vector of new nodes
1887  new_node_pt.push_back(created_node_pt);
1888  }
1889  }
1890 
1891  // In the first instance use macro element or FE representation
1892  // to create past and present nodal positions.
1893  // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
1894  // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
1895  // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
1896  // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
1897  // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
1898  // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
1899  // NOT ASSIGN SENSIBLE INITIAL POSITONS!
1900 
1901 
1902  // Have we created a new node?
1903  if (!node_done)
1904  {
1905  // Loop over # of history values
1906  for (unsigned t = 0; t < ntstorage; t++)
1907  {
1908  // Get position from father element -- this uses the macro
1909  // element representation if appropriate. If the node
1910  // turns out to be a hanging node later on, then
1911  // its position gets adjusted in line with its
1912  // hanging node interpolation.
1913  Vector<double> x_prev(n_dim);
1914  father_el_pt->get_x(t, s, x_prev);
1915 
1916  // Set previous positions of the new node
1917  for (unsigned i = 0; i < n_dim; i++)
1918  {
1919  created_node_pt->x(t, i) = x_prev[i];
1920  }
1921  }
1922 
1923  // Now set the values
1924  // Loop over all history values
1925  for (unsigned t = 0; t < ntstorage; t++)
1926  {
1927  // Get values from father element
1928  // Note: get_interpolated_values() sets Vector size itself.
1929  Vector<double> prev_values;
1930  father_el_pt->get_interpolated_values(t, s, prev_values);
1931 
1932  // Initialise the values at the new node
1933  unsigned n_value = created_node_pt->nvalue();
1934  for (unsigned k = 0; k < n_value; k++)
1935  {
1936  created_node_pt->set_value(t, k, prev_values[k]);
1937  }
1938  }
1939 
1940  // Add new node to mesh
1941  mesh_pt->add_node_pt(created_node_pt);
1942 
1943  } // End of whether the node has been created by us or not
1944 
1945  } // End of if node already existed in father in which case
1946  // everything above gets bypassed.
1947 
1948  // Check if the element is an algebraic element
1949  AlgebraicElementBase* alg_el_pt =
1950  dynamic_cast<AlgebraicElementBase*>(this);
1951 
1952  // If the element is an algebraic element, setup
1953  // node position (past and present) from algebraic update
1954  // function.
1955  // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
1956  // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
1957  // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
1958  // INFO FOR *ALL* ROOT ELEMENTS!
1959  if (alg_el_pt != 0)
1960  {
1961  // Build algebraic node update info for new node
1962  // This sets up the node update data for all node update
1963  // functions that are shared by all nodes in the father
1964  // element
1965  alg_el_pt->setup_algebraic_node_update(
1966  node_pt(jnod), s, father_el_pt);
1967  }
1968 
1969  // We have built the node and we are documenting
1970  if ((!node_done) && (new_nodes_file.is_open()))
1971  {
1972  new_nodes_file << node_pt(jnod)->x(0) << " "
1973  << node_pt(jnod)->x(1) << " "
1974  << node_pt(jnod)->x(2) << std::endl;
1975  }
1976 
1977  } // End of Z loop over nodes in element
1978 
1979  } // End of vertical loop over nodes in element
1980 
1981  } // End of horizontal loop over nodes in element
1982 
1983  // If the element is a MacroElementNodeUpdateElement, set
1984  // the update parameters for the current element's nodes --
1985  // all this needs is the vector of (pointers to the)
1986  // geometric objects that affect the MacroElement-based
1987  // node update -- this is the same as that in the father element
1988  MacroElementNodeUpdateElementBase* father_m_el_pt =
1989  dynamic_cast<MacroElementNodeUpdateElementBase*>(father_el_pt);
1990  if (father_m_el_pt != 0)
1991  {
1992  // Get vector of geometric objects from father (construct vector
1993  // via copy operation)
1994  Vector<GeomObject*> geom_object_pt(father_m_el_pt->geom_object_pt());
1995 
1996  // Cast current element to MacroElementNodeUpdateElement:
1997  MacroElementNodeUpdateElementBase* m_el_pt =
1998  dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
1999 
2000 #ifdef PARANOID
2001  if (m_el_pt == 0)
2002  {
2003  std::string error_message =
2004  "Failed to cast to MacroElementNodeUpdateElementBase*\n";
2005  error_message +=
2006  "Strange -- if the father is a MacroElementNodeUpdateElement\n";
2007  error_message += "the son should be too....\n";
2008 
2009  throw OomphLibError(
2010  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2011  }
2012 #endif
2013  // Build update info by passing vector of geometric objects:
2014  // This sets the current element to be the update element
2015  // for all of the element's nodes -- this is reversed
2016  // if the element is ever un-refined in the father element's
2017  // rebuild_from_sons() function which overwrites this
2018  // assignment to avoid nasty segmentation faults that occur
2019  // when a node tries to update itself via an element that no
2020  // longer exists...
2021  m_el_pt->set_node_update_info(geom_object_pt);
2022  }
2023 
2024 #ifdef OOMPH_HAS_MPI
2025  // Is the new element a halo element?
2026  if (tree_pt()->father_pt()->object_pt()->is_halo())
2027  {
2028  Non_halo_proc_ID =
2029  tree_pt()->father_pt()->object_pt()->non_halo_proc_ID();
2030  }
2031 #endif
2032 
2033  // Is it an ElementWithMovingNodes?
2034  ElementWithMovingNodes* aux_el_pt =
2035  dynamic_cast<ElementWithMovingNodes*>(this);
2036 
2037  // Pass down the information re the method for the evaluation
2038  // of the shape derivatives
2039  if (aux_el_pt != 0)
2040  {
2041  ElementWithMovingNodes* aux_father_el_pt =
2042  dynamic_cast<ElementWithMovingNodes*>(father_el_pt);
2043 
2044 #ifdef PARANOID
2045  if (aux_father_el_pt == 0)
2046  {
2047  std::string error_message =
2048  "Failed to cast to ElementWithMovingNodes*\n";
2049  error_message +=
2050  "Strange -- if the son is a ElementWithMovingNodes\n";
2051  error_message += "the father should be too....\n";
2052 
2053  throw OomphLibError(
2054  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2055  }
2056 #endif
2057 
2058  // If evaluating the residuals by finite differences in the father
2059  // continue to do so in the child
2060  if (aux_father_el_pt
2061  ->are_dresidual_dnodal_coordinates_always_evaluated_by_fd())
2062  {
2063  aux_el_pt
2064  ->enable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
2065  }
2066 
2067 
2068  aux_el_pt->method_for_shape_derivs() =
2069  aux_father_el_pt->method_for_shape_derivs();
2070 
2071  // If bypassing the evaluation of fill_in_jacobian_from_geometric_data
2072  // continue to do so
2073  if (aux_father_el_pt
2074  ->is_fill_in_jacobian_from_geometric_data_bypassed())
2075  {
2076  aux_el_pt->enable_bypass_fill_in_jacobian_from_geometric_data();
2077  }
2078  }
2079 
2080  // Now do further build (if any)
2081  further_build();
2082 
2083  } // Sanity check: Father element has been generated
2084 
2085  } // End for element has not been built yet
2086  else
2087  {
2088  was_already_built = true;
2089  }
2090  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
oomph::Data::nvalue
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
oomph::FiniteElement::set_macro_elem_pt
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Definition: elements.h:1872
oomph::FiniteElement::construct_node
virtual Node * construct_node(const unsigned &n)
Definition: elements.h:2509
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::nnode_1d
virtual unsigned nnode_1d() const
Definition: elements.h:2218
oomph::FiniteElement::construct_boundary_node
virtual Node * construct_boundary_node(const unsigned &n)
Definition: elements.h:2538
oomph::FiniteElement::local_one_d_fraction_of_node
virtual double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
Definition: elements.h:1858
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::OcTree::Direction_to_vector
static Vector< Vector< int > > Direction_to_vector
Definition: octree.h:353
oomph::QElementBase::s_macro_ll
double & s_macro_ll(const unsigned &i)
Definition: Qelements.h:186
oomph::QElementBase::s_macro_ur
double & s_macro_ur(const unsigned &i)
Definition: Qelements.h:202
oomph::RefineableElement::tree_pt
Tree * tree_pt()
Access function: Pointer to quadtree representation of this element.
Definition: refineable_elements.h:211
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::nodes_built
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
Definition: refineable_elements.h:389
oomph::RefineableElement::ncont_interpolated_values
virtual unsigned ncont_interpolated_values() const =0
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::RefineableQElement< 3 >::node_created_by_son_of_neighbour
virtual Node * node_created_by_son_of_neighbour(const Vector< double > &s_fraction, bool &is_periodic)
Definition: refineable_brick_element.h:118
oomph::RefineableQElement< 3 >::node_created_by_neighbour
virtual Node * node_created_by_neighbour(const Vector< double > &s_fraction, bool &is_periodic)
Definition: refineable_brick_element.cc:992
oomph::RefineableQElement< 3 >::octree_pt
OcTree * octree_pt()
Pointer to octree representation of this element.
Definition: refineable_brick_element.h:144
oomph::TimeStepper::ntstorage
unsigned ntstorage() const
Definition: timesteppers.h:601
oomph::Tree::object_pt
RefineableElement * object_pt() const
Definition: tree.h:88
oomph::Tree::father_pt
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:235
oomph::Tree::son_type
int son_type() const
Return son type.
Definition: tree.h:214
oomph::Tree::OMEGA
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:262
s
RealScalar s
Definition: level1_cplx_impl.h:130
k
char char char int int * k
Definition: level2_impl.h:374
Eigen::zeta
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
plotDoE.x
list x
Definition: plotDoE.py:28
plotPSD.t
t
Definition: plotPSD.py:36
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::Mesh::add_boundary_node(), oomph::Mesh::add_node_pt(), oomph::Mesh::boundary_coordinate_exists(), oomph::ElementWithMovingNodes::enable_always_evaluate_dresidual_dnodal_coordinates_by_fd(), oomph::ElementWithMovingNodes::enable_bypass_fill_in_jacobian_from_geometric_data(), oomph::MacroElementNodeUpdateElementBase::geom_object_pt(), get_bcs(), get_boundaries(), oomph::RefineableElement::get_interpolated_values(), oomph::FiniteElement::get_node_at_local_coordinate(), oomph::RefineableSolidQElement< 3 >::get_solid_bcs(), oomph::FiniteElement::get_x(), i, interpolated_zeta_on_face(), k, oomph::FiniteElement::macro_elem_pt(), oomph::Node::make_periodic(), oomph::ElementWithMovingNodes::method_for_shape_derivs(), oomph::Node::ndim(), oomph::FiniteElement::node_pt(), oomph::Node::nposition_type(), oomph::Data::nvalue(), oomph::Tree::object_pt(), oomph::Tree::OMEGA, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::Data::pin(), oomph::SolidNode::pin_position(), s, oomph::QElementBase::s_macro_ll(), oomph::QElementBase::s_macro_ur(), oomph::Node::set_coordinates_on_boundary(), oomph::MacroElementNodeUpdateElementBase::set_node_update_info(), oomph::Data::set_value(), oomph::AlgebraicElementBase::setup_algebraic_node_update(), oomph::Tree::son_type(), oomph::Global_string_for_annotation::string(), plotPSD::t, oomph::Data::time_stepper_pt(), plotDoE::x, oomph::Node::x(), and Eigen::zeta().

◆ check_integrity()

void oomph::RefineableQElement< 3 >::check_integrity ( double &  max_error)
virtual

Check the integrity of the element: ensure that the position and values are continuous across the element faces

Check inter-element continuity of

  • nodal positions
  • (nodally) interpolated function values

Implements oomph::RefineableElement.

2510  {
2511  using namespace OcTreeNames;
2512 
2513  // std::ofstream error_file("errors.dat");
2514 
2515  // Number of nodes along edge
2516  unsigned n_p = nnode_1d();
2517 
2518  // Number of timesteps (incl. present) for which continuity is
2519  // to be checked.
2520  unsigned n_time = 1;
2521 
2522  // Initialise errors
2523  max_error = 0.0;
2524  Vector<double> max_error_x(3, 0.0);
2525  double max_error_val = 0.0;
2526 
2527  // Set the faces
2528  Vector<int> faces(6);
2529  faces[0] = D;
2530  faces[1] = U;
2531  faces[2] = L;
2532  faces[3] = R;
2533  faces[4] = B;
2534  faces[5] = F;
2535 
2536  // Loop over the edges
2537  for (unsigned face_counter = 0; face_counter < 6; face_counter++)
2538  {
2539  Vector<double> s(3), s_lo_neigh(3), s_hi_neigh(3);
2540  Vector<double> s_fraction(3);
2541  Vector<unsigned> translate_s(3);
2542  int neigh_face, diff_level;
2543  int my_face;
2544  bool in_neighbouring_tree;
2545 
2546  my_face = faces[face_counter];
2547 
2548  // Find pointer to neighbour in this direction
2549  OcTree* neigh_pt;
2550  neigh_pt = octree_pt()->gteq_face_neighbour(my_face,
2551  translate_s,
2552  s_lo_neigh,
2553  s_hi_neigh,
2554  neigh_face,
2555  diff_level,
2556  in_neighbouring_tree);
2557 
2558  // Neighbour exists and has existing nodes
2559  if ((neigh_pt != 0) && (neigh_pt->object_pt()->nodes_built()))
2560  {
2561  // if (diff_level!=0)
2562  {
2563  // Need to exclude periodic nodes from this check
2564  // There are only periodic nodes if we are in a neighbouring tree
2565  bool is_periodic = false;
2566  if (in_neighbouring_tree)
2567  {
2568  // Is it periodic
2569  is_periodic =
2570  tree_pt()->root_pt()->is_neighbour_periodic(faces[face_counter]);
2571  }
2572 
2573  // Loop over nodes along the edge
2574  for (unsigned i0 = 0; i0 < n_p; i0++)
2575  {
2576  for (unsigned i1 = 0; i1 < n_p; i1++)
2577  {
2578  // Local node number
2579  unsigned n = 0;
2580  switch (face_counter)
2581  {
2582  case 0:
2583  // Fractions
2584  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2585  s_fraction[1] = 0.0;
2586  s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2587  // Set local node number
2588  n = i0 + i1 * n_p * n_p;
2589  break;
2590 
2591  case 1:
2592  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2593  s_fraction[1] = 1.0;
2594  s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2595  // Set local node number
2596  n = i0 + n_p * (n_p - 1) + i1 * n_p * n_p;
2597  break;
2598 
2599  case 2:
2600  s_fraction[0] = 0.0;
2601  s_fraction[1] = local_one_d_fraction_of_node(i0, 1);
2602  s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2603  // Set local node number
2604  n = n_p * i0 + i1 * n_p * n_p;
2605  break;
2606 
2607  case 3:
2608  s_fraction[0] = 1.0;
2609  s_fraction[1] = local_one_d_fraction_of_node(i0, 1);
2610  s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2611  // Set local node number
2612  n = n_p - 1 + n_p * i0 + n_p * n_p * i1;
2613  break;
2614 
2615  case 4:
2616  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2617  s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
2618  s_fraction[2] = 0.0;
2619  // Set local node number
2620  n = i0 + n_p * i1;
2621  break;
2622 
2623  case 5:
2624  s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2625  s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
2626  s_fraction[2] = 1.0;
2627  // Set local node number
2628  n = i0 + n_p * i1 + n_p * n_p * (n_p - 1);
2629  break;
2630  }
2631 
2632 
2633  // We have to check if the hi and lo directions along the
2634  // face are inverted or not
2635  Vector<double> s_in_neighb(3);
2636  for (unsigned i = 0; i < 3; i++)
2637  {
2638  s[i] = -1.0 + 2.0 * s_fraction[i];
2639  s_in_neighb[i] =
2640  s_lo_neigh[i] +
2641  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
2642  }
2643 
2644 
2645  // Loop over timesteps
2646  for (unsigned t = 0; t < n_time; t++)
2647  {
2648  // Get the nodal position from neighbour element
2649  Vector<double> x_in_neighb(3);
2650  neigh_pt->object_pt()->interpolated_x(
2651  t, s_in_neighb, x_in_neighb);
2652 
2653  // Check error only if the node is NOT periodic
2654  if (is_periodic == false)
2655  {
2656  // Check error
2657  for (unsigned i = 0; i < 3; i++)
2658  {
2659  // Find the spatial error
2660  double err =
2661  std::fabs(node_pt(n)->x(t, i) - x_in_neighb[i]);
2662 
2663  // If it's bigger than our tolerance, say so
2664  if (err > 1e-9)
2665  {
2666  oomph_info << "errx[" << i << "], t x, x_neigh: " << err
2667  << " " << t << " " << node_pt(n)->x(t, i)
2668  << " " << x_in_neighb[i] << std::endl;
2669  oomph_info << "at " << node_pt(n)->x(0) << " "
2670  << node_pt(n)->x(1) << " " << node_pt(n)->x(2)
2671  << " " << std::endl;
2672  }
2673 
2674  // If it's bigger than the previous max error, it is the
2675  // new max error!
2676  if (err > max_error_x[i])
2677  {
2678  max_error_x[i] = err;
2679  }
2680  }
2681  }
2682 
2683  // Get the values from neighbour element. Note: # of values
2684  // gets set by routine (because in general we don't know
2685  // how many interpolated values a certain element has
2686  Vector<double> values_in_neighb;
2687  neigh_pt->object_pt()->get_interpolated_values(
2688  t, s_in_neighb, values_in_neighb);
2689 
2690  // Get the values in current element.
2691  Vector<double> values;
2692  get_interpolated_values(t, s, values);
2693 
2694  // Now figure out how many continuously interpolated
2695  // values there are
2696  unsigned num_val =
2697  neigh_pt->object_pt()->ncont_interpolated_values();
2698 
2699  // Check error
2700  for (unsigned ival = 0; ival < num_val; ival++)
2701  {
2702  double err = std::fabs(values[ival] - values_in_neighb[ival]);
2703 
2704  if (err > 1.0e-10)
2705  {
2706  oomph_info << node_pt(n)->x(0) << " " << node_pt(n)->x(1)
2707  << " " << node_pt(n)->x(2) << " \n# "
2708  << "erru (S)" << err << " " << ival << " " << n
2709  << " " << values[ival] << " "
2710  << values_in_neighb[ival] << std::endl;
2711  // error_file<<"ZONE"<<std::endl
2712  // << node_pt(n)->x(0) << " "
2713  // << node_pt(n)->x(1) << " "
2714  // << node_pt(n)->x(2) << std::endl;
2715  }
2716 
2717  if (err > max_error_val)
2718  {
2719  max_error_val = err;
2720  }
2721  }
2722  }
2723  }
2724  }
2725  }
2726  }
2727  }
2728 
2729  max_error = max_error_x[0];
2730  if (max_error_x[1] > max_error) max_error = max_error_x[1];
2731  if (max_error_x[2] > max_error) max_error = max_error_x[2];
2732  if (max_error_val > max_error) max_error = max_error_val;
2733 
2734  if (max_error > 1e-9)
2735  {
2736  oomph_info << "\n#------------------------------------ \n#Max error ";
2737  oomph_info << max_error_x[0] << " " << max_error_x[1] << " "
2738  << max_error_x[2] << " " << max_error_val << std::endl;
2739  oomph_info << "#------------------------------------ \n " << std::endl;
2740  }
2741 
2742  // error_file.close();
2743  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
D
dominoes D
Definition: Domino.cpp:55
L
MatrixXd L
Definition: LLT_example.cpp:6
R
@ R
Definition: StatisticsVector.h:21
B
Matrix< SCALARB, Dynamic, Dynamic, opt_B > B
Definition: bench_gemm.cpp:48
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::RefineableElement::get_interpolated_values
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Definition: refineable_elements.h:417
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
MeshRefinement::max_error
double max_error
Definition: MortaringCantileverCompareToNonMortaring.cpp:188
RachelsAdvectionDiffusion::U
double U
Swimming speed.
Definition: two_d_variable_diff_adapt.cc:53
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
oomph::OcTreeNames::F
@ F
Definition: octree.h:74
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

References D, e(), oomph::OcTreeNames::F, boost::multiprecision::fabs(), oomph::RefineableElement::get_interpolated_values(), oomph::OcTree::gteq_face_neighbour(), i, oomph::FiniteElement::interpolated_x(), L, MeshRefinement::max_error, n, oomph::RefineableElement::ncont_interpolated_values(), oomph::RefineableElement::nodes_built(), oomph::Tree::object_pt(), oomph::oomph_info, R, s, plotPSD::t, RachelsAdvectionDiffusion::U, and plotDoE::x.

◆ further_setup_hanging_nodes()

virtual void oomph::RefineableQElement< 3 >::further_setup_hanging_nodes ( )
pure virtual

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

Reimplemented from oomph::RefineableElement.

Implemented in oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >.

◆ get_bcs()

void oomph::RefineableQElement< 3 >::get_bcs ( int  bound,
Vector< int > &  bound_cons 
) const
protected

Determine Vector of boundary conditions along the element's boundary (or vertex) bound (S/W/N/E/SW/SE/NW/NE).

This function assumes that the same boundary condition is applied along the entire area of an element's face (of course, the vertices combine the boundary conditions of their two adjacent edges in the most restrictive combination. Hence, if we're at a vertex, we apply the most restrictive boundary condition of the two adjacent edges. If we're on an edge (in its proper interior), we apply the least restrictive boundary condition of all nodes along the edge.

Usual convention:

  • bound_cons[ival]=0 if value ival on this boundary is free
  • bound_cons[ival]=1 if value ival on this boundary is pinned

Determine Vector of boundary conditions along the element's boundary bound.

This function assumes that the same boundary condition is applied along the entire area of an element's face (of course, the vertices combine the boundary conditions of their two adjacent edges in the most restrictive combination. Hence, if we're at a vertex, we apply the most restrictive boundary condition of the two adjacent edges. If we're on an edge (in its proper interior), we apply the least restrictive boundary condition of all nodes along the edge.

Usual convention:

  • bound_cons[ival]=0 if value ival on this boundary is free
  • bound_cons[ival]=1 if value ival on this boundary is pinned
239  {
240  using namespace OcTreeNames;
241 
242  // Max. number of nodal data values in element
243  unsigned nvalue = ncont_interpolated_values();
244  // Set up temporary vectors to hold edge boundary conditions
245  Vector<int> bound_cons1(nvalue), bound_cons2(nvalue);
246  Vector<int> bound_cons3(nvalue);
247 
248  Vector<int> vect1(3), vect2(3), vect3(3);
249  Vector<int> vect_elem;
250  Vector<int> notzero;
251  int n = 0;
252 
253  vect_elem = OcTree::Direction_to_vector[bound];
254 
255  // Just to see if bound is a face, an edge, or a vertex, n stores
256  // the number of non-zero values in the vector reprensenting the bound
257  // and the vector notzero stores the position of these values
258  for (int i = 0; i < 3; i++)
259  {
260  if (vect_elem[i] != 0)
261  {
262  n++;
263  notzero.push_back(i);
264  }
265  }
266 
267  switch (n)
268  {
269  // If there is only one non-zero value, bound is a face
270  case 1:
271  get_face_bcs(bound, bound_cons);
272  break;
273 
274  // If there are two non-zero values, bound is an edge
275  case 2:
276 
277  for (unsigned i = 0; i < 3; i++)
278  {
279  vect1[i] = 0;
280  vect2[i] = 0;
281  }
282  // vect1 and vect2 are the vector of the two faces adjacent to bound
283  vect1[notzero[0]] = vect_elem[notzero[0]];
284  vect2[notzero[1]] = vect_elem[notzero[1]];
285 
286  get_face_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
287  get_face_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
288  // get the most restrictive bc
289  for (unsigned k = 0; k < nvalue; k++)
290  {
291  bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
292  }
293  break;
294 
295  // If there are three non-zero value, bound is a vertex
296  case 3:
297 
298  for (unsigned i = 0; i < 3; i++)
299  {
300  vect1[i] = 0;
301  vect2[i] = 0;
302  vect3[i] = 0;
303  }
304  // vectors to the three adjacent faces of the vertex
305  vect1[0] = vect_elem[0];
306  vect2[1] = vect_elem[1];
307  vect3[2] = vect_elem[2];
308 
309  get_face_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
310  get_face_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
311  get_face_bcs(OcTree::Vector_to_direction[vect3], bound_cons3);
312 
313 
314  // set the bcs to the most restrictive ones
315  for (unsigned k = 0; k < nvalue; k++)
316  {
317  bound_cons[k] = (bound_cons1[k] || bound_cons2[k] || bound_cons3[k]);
318  }
319  break;
320 
321  default:
322  throw OomphLibError("Make sure you are not giving OMEGA as bound",
323  OOMPH_CURRENT_FUNCTION,
324  OOMPH_EXCEPTION_LOCATION);
325  }
326  }
oomph::OcTree::Vector_to_direction
static std::map< Vector< int >, int > Vector_to_direction
Definition: octree.h:348
oomph::RefineableQElement< 3 >::get_face_bcs
void get_face_bcs(const int &edge, Vector< int > &bound_cons) const
Definition: refineable_brick_element.cc:337

References oomph::OcTree::Direction_to_vector, i, k, n, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::OcTree::Vector_to_direction.

Referenced by build(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::p_refine(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::rebuild_from_sons(), and oomph::RefineableQSpectralElement< 3 >::rebuild_from_sons().

◆ get_boundaries()

void oomph::RefineableQElement< 3 >::get_boundaries ( const int &  element,
std::set< unsigned > &  boundary 
) const
protected

Given an element edge/vertex, return a Vector which contains all the (mesh-)boundary numbers that this element edge/vertex lives on.

For proper edges, the boundary is the one (if any) that is shared by both vertex nodes). For vertex nodes, we just return their boundaries.

426  {
427  using namespace OcTreeNames;
428 
429  // Number of 1d nodes along an edge
430  unsigned n_p = nnode_1d();
431  // Left and right-hand nodes
432  int node[4];
433  int a = 0, b = 0;
434  int n = 0;
435  Vector<int> vect_face(3);
436  vect_face = OcTree::Direction_to_vector[element];
437  // Set the left (lower) and right (upper) nodes for the edge
438 
439  // this is to know what is the type of element (face, edge, vertex)
440  // we just need to count the number of values equal to 0 in the
441  // vector representing this element
442 
443  // Local storage for the node-numbers in given directions
444  // Initialise to zero (assume at LH end of each domain)
445  int one_d_node_number[3] = {0, 0, 0};
446  // n stores the number of values equal to 0, a is the position of the
447  // last 0-value ;b is the position of the last non0-value
448  for (int i = 0; i < 3; i++)
449  {
450  if (vect_face[i] == 0)
451  {
452  a = i;
453  n++;
454  }
455  else
456  {
457  b = i;
458  // If we are at the extreme (RH) end of the face,
459  // set the node number accordingly
460  if (vect_face[i] == 1)
461  {
462  one_d_node_number[i] = n_p - 1;
463  }
464  }
465  }
466 
467  switch (n)
468  {
469  // if n=0 element is a vertex, and need to look at only one node
470  case 0:
471  node[0] = one_d_node_number[0] + n_p * one_d_node_number[1] +
472  n_p * n_p * one_d_node_number[2];
473  node[1] = node[0];
474  node[2] = node[0];
475  node[3] = node[0];
476  break;
477 
478  // if n=1 element is an edge, and we need to look at two nodes
479  case 1:
480  if (a == 0)
481  {
482  node[0] = (n_p - 1) + n_p * one_d_node_number[1] +
483  n_p * n_p * one_d_node_number[2];
484  node[1] =
485  n_p * one_d_node_number[1] + n_p * n_p * one_d_node_number[2];
486  }
487  else if (a == 1)
488  {
489  node[0] = n_p * (n_p - 1) + one_d_node_number[0] +
490  n_p * n_p * one_d_node_number[2];
491  node[1] = one_d_node_number[0] + n_p * n_p * one_d_node_number[2];
492  }
493  else if (a == 2)
494  {
495  node[0] = one_d_node_number[0] + n_p * one_d_node_number[1] +
496  n_p * n_p * (n_p - 1);
497  node[1] = one_d_node_number[0] + n_p * one_d_node_number[1];
498  }
499  node[2] = node[1];
500  node[3] = node[1];
501  break;
502 
503  // if n=2 element is a face, and we need to look at its 4 nodes
504  case 2:
505  if (b == 0)
506  {
507  node[0] =
508  one_d_node_number[0] + n_p * n_p * (n_p - 1) + n_p * (n_p - 1);
509  node[1] = one_d_node_number[0] + n_p * (n_p - 1);
510  node[2] = one_d_node_number[0] + n_p * n_p * (n_p - 1);
511  node[3] = one_d_node_number[0];
512  }
513  else if (b == 1)
514  {
515  node[0] =
516  n_p * one_d_node_number[1] + n_p * n_p * (n_p - 1) + (n_p - 1);
517  node[1] = n_p * one_d_node_number[1] + (n_p - 1);
518  node[2] = n_p * one_d_node_number[1] + n_p * n_p * (n_p - 1);
519  node[3] = n_p * one_d_node_number[1];
520  }
521  else if (b == 2)
522  {
523  node[0] =
524  n_p * n_p * one_d_node_number[2] + n_p * (n_p - 1) + (n_p - 1);
525  node[1] = n_p * n_p * one_d_node_number[2] + (n_p - 1);
526  node[2] = n_p * n_p * one_d_node_number[2] + n_p * (n_p - 1);
527  node[3] = n_p * n_p * one_d_node_number[2];
528  }
529  break;
530  default:
531  throw OomphLibError("Make sure you are not giving OMEGA as boundary",
532  OOMPH_CURRENT_FUNCTION,
533  OOMPH_EXCEPTION_LOCATION);
534  }
535 
536 
537  // Empty boundary set: Edge does not live on any boundary
538  boundary.clear();
539 
540  // Storage for the boundaries at the four nodes
541  Vector<std::set<unsigned>*> node_boundaries_pt(4, 0);
542 
543  // Loop over the four nodes and get the boundary information
544  for (unsigned i = 0; i < 4; i++)
545  {
546  node_pt(node[i])->get_boundaries_pt(node_boundaries_pt[i]);
547  }
548 
549 
550  // Now work out the intersections
551  Vector<std::set<unsigned>> boundary_aux(2);
552 
553  for (unsigned i = 0; i < 2; i++)
554  {
555  // If the two nodes both lie on boundaries
556  if ((node_boundaries_pt[2 * i] != 0) &&
557  (node_boundaries_pt[2 * i + 1] != 0))
558  {
559  // Find the intersection (the common boundaries) of these nodes
560  std::set_intersection(
561  node_boundaries_pt[2 * i]->begin(),
562  node_boundaries_pt[2 * i]->end(),
563  node_boundaries_pt[2 * i + 1]->begin(),
564  node_boundaries_pt[2 * i + 1]->end(),
565  inserter(boundary_aux[i], boundary_aux[i].begin()));
566  }
567  }
568 
569  // Now calculate the total intersection
570  set_intersection(boundary_aux[0].begin(),
571  boundary_aux[0].end(),
572  boundary_aux[1].begin(),
573  boundary_aux[1].end(),
574  inserter(boundary, boundary.begin()));
575  }
b
Scalar * b
Definition: benchVecAdd.cpp:17
oomph::Node::get_boundaries_pt
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Definition: nodes.h:1365
Eigen::placeholders::end
static constexpr lastp1_t end
Definition: IndexedViewHelper.h:79
a
const Scalar * a
Definition: level2_cplx_impl.h:32

References a, b, oomph::OcTree::Direction_to_vector, Eigen::placeholders::end, i, n, OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

Referenced by build(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::p_refine(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::rebuild_from_sons(), and oomph::RefineableQSpectralElement< 3 >::rebuild_from_sons().

◆ get_face_bcs()

void oomph::RefineableQElement< 3 >::get_face_bcs ( const int &  face,
Vector< int > &  bound_cons 
) const
protected

Determine Vector of boundary conditions along the element's face (R/L/U/D/B/F) – BC is the least restrictive combination of all the nodes on this face.

Usual convention:

  • bound_cons[ival]=0 if value ival on this boundary is free
  • bound_cons[ival]=1 if value ival on this boundary is pinned

Determine Vector of boundary conditions along the element's face (R/L/U/D/B/F) – BC is the least restrictive combination of all the nodes on this face

Usual convention:

  • bound_cons[ival]=0 if value ival on this boundary is free
  • bound_cons[ival]=1 if value ival on this boundary is pinned
339  {
340  using namespace OcTreeNames;
341 
342  // Number of nodes along 1D edge
343  unsigned n_p = nnode_1d();
344  // the four corner nodes on the boundary
345  unsigned node1, node2, node3, node4;
346 
347  // Set the four corner nodes for the face
348  switch (face)
349  {
350  case U:
351  node1 = n_p * n_p * n_p - 1;
352  node2 = n_p * n_p - 1;
353  node3 = n_p * (n_p - 1);
354  node4 = n_p * (n_p * n_p - 1);
355  break;
356 
357  case D:
358  node1 = 0;
359  node2 = n_p - 1;
360  node3 = (n_p * n_p + 1) * (n_p - 1);
361  node4 = n_p * n_p * (n_p - 1);
362  break;
363 
364  case R:
365  node1 = n_p - 1;
366  node2 = (n_p * n_p + 1) * (n_p - 1);
367  node3 = n_p * n_p * n_p - 1;
368  node4 = n_p * n_p - 1;
369  break;
370 
371  case L:
372  node1 = 0;
373  node2 = n_p * (n_p - 1);
374  node3 = n_p * (n_p * n_p - 1);
375  node4 = n_p * n_p * (n_p - 1);
376  break;
377 
378  case B:
379  node1 = 0;
380  node2 = n_p - 1;
381  node3 = n_p * n_p - 1;
382  node4 = n_p * (n_p - 1);
383  break;
384 
385  case F:
386  node1 = n_p * n_p * n_p - 1;
387  node2 = n_p * (n_p * n_p - 1);
388  node3 = n_p * n_p * (n_p - 1);
389  node4 = (n_p - 1) * (n_p * n_p + 1);
390  break;
391 
392  default:
393  std::ostringstream error_stream;
394  error_stream << "Wrong edge " << face << " passed\n";
395 
396  throw OomphLibError(
397  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
398  }
399 
400  // Max. number of nodal data values in element
401  unsigned maxnvalue = ncont_interpolated_values();
402 
403  // Loop over all values, the least restrictive value is
404  // the multiplication of the boundary conditions at the 4 nodes
405  // Assuming that free is always zero and pinned is one
406  for (unsigned k = 0; k < maxnvalue; k++)
407  {
408  bound_cons[k] =
409  node_pt(node1)->is_pinned(k) * node_pt(node2)->is_pinned(k) *
410  node_pt(node3)->is_pinned(k) * node_pt(node4)->is_pinned(k);
411  }
412  }
oomph::Data::is_pinned
bool is_pinned(const unsigned &i) const
Test whether the i-th variable is pinned (1: true; 0: false).
Definition: nodes.h:417

References D, oomph::OcTreeNames::F, k, L, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, and RachelsAdvectionDiffusion::U.

◆ interpolated_zeta_on_face()

void oomph::RefineableQElement< 3 >::interpolated_zeta_on_face ( const unsigned &  boundary,
const int &  face,
const Vector< double > &  s,
Vector< double > &  zeta 
)
protected

Return the value of the intrinsic boundary coordinate interpolated along the face

587  {
588  using namespace OcTreeNames;
589 
590  // Number of nodes along an edge
591  unsigned nnodes_1d = nnode_1d();
592 
593  // Number of nodes on a face
594  unsigned nnodes_2d = nnodes_1d * nnodes_1d;
595 
596  // Total number of nodes
597  unsigned nnodes_3d = nnode();
598 
599  // Storage for the shape functions
600  Shape psi(nnodes_3d);
601 
602  // Get the shape functions at the passed position
603  this->shape(s, psi);
604 
605  // Unsigned data that give starts and increments for the loop
606  // over the nodes on the faces.
607  unsigned start = 0, increment1 = 1, increment2 = 1;
608 
609  // Flag to record if actually on a face or an edge
610  bool on_edge = true;
611 
612  // Which face?
613  switch (face)
614  {
615  case L:
616 #ifdef PARANOID
617  if (s[0] != -1.0)
618  {
619  std::ostringstream error_stream;
620  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
621  << " is not on Left face\n";
622 
623  throw OomphLibError(error_stream.str(),
624  OOMPH_CURRENT_FUNCTION,
625  OOMPH_EXCEPTION_LOCATION);
626  }
627 #endif
628  // Start is zero (bottom-left-back corner)
629  increment1 = nnodes_1d;
630  increment2 = 0;
631  on_edge = false;
632  break;
633 
634  case R:
635 #ifdef PARANOID
636  if (s[0] != 1.0)
637  {
638  std::ostringstream error_stream;
639  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
640  << " is not on Right face\n";
641 
642  throw OomphLibError(error_stream.str(),
643  OOMPH_CURRENT_FUNCTION,
644  OOMPH_EXCEPTION_LOCATION);
645  }
646 #endif
647  // Start is bottom-right-back corner
648  start = nnodes_1d - 1;
649  increment1 = nnodes_1d;
650  increment2 = 0;
651  on_edge = false;
652  break;
653 
654  case D:
655 #ifdef PARANOID
656  if (s[1] != -1.0)
657  {
658  std::ostringstream error_stream;
659  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
660  << " is not on Bottom face\n";
661 
662  throw OomphLibError(error_stream.str(),
663  OOMPH_CURRENT_FUNCTION,
664  OOMPH_EXCEPTION_LOCATION);
665  }
666 #endif
667  // Start is zero and increments2 is nnode_2d-nnode_1d
668  increment2 = nnodes_2d - nnodes_1d;
669  on_edge = false;
670  break;
671 
672  case U:
673 #ifdef PARANOID
674  if (s[1] != 1.0)
675  {
676  std::ostringstream error_stream;
677  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
678  << " is not on Upper face\n";
679 
680  throw OomphLibError(error_stream.str(),
681  OOMPH_CURRENT_FUNCTION,
682  OOMPH_EXCEPTION_LOCATION);
683  }
684 #endif
685  // Start is top-left-back corner and increments2 is nnode_2d-nnode_1d
686  start = nnodes_2d - nnodes_1d;
687  increment2 = start;
688  on_edge = false;
689  break;
690 
691  case B:
692 #ifdef PARANOID
693  if (s[2] != -1.0)
694  {
695  std::ostringstream error_stream;
696  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
697  << " is not on Back face\n";
698 
699  throw OomphLibError(error_stream.str(),
700  OOMPH_CURRENT_FUNCTION,
701  OOMPH_EXCEPTION_LOCATION);
702  }
703 #endif
704  // Start is zero and increments are 1 and 0
705  increment2 = 0;
706  on_edge = false;
707  break;
708 
709  case F:
710 #ifdef PARANOID
711  if (s[2] != 1.0)
712  {
713  std::ostringstream error_stream;
714  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
715  << " is not on Front face\n";
716 
717  throw OomphLibError(error_stream.str(),
718  OOMPH_CURRENT_FUNCTION,
719  OOMPH_EXCEPTION_LOCATION);
720  }
721 #endif
722  // Start is bottom-left-front corner
723  start = nnodes_3d - nnodes_2d;
724  increment2 = 0;
725  on_edge = false;
726  break;
727 
728  case LF:
729 #ifdef PARANOID
730  if ((s[0] != -1.0) || (s[2] != 1.0))
731  {
732  std::ostringstream error_stream;
733  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
734  << " is not on Front-Left edge\n";
735 
736  throw OomphLibError(error_stream.str(),
737  OOMPH_CURRENT_FUNCTION,
738  OOMPH_EXCEPTION_LOCATION);
739  }
740 #endif
741  start = nnodes_3d - nnodes_2d;
742  increment1 = nnodes_1d;
743  break;
744 
745  case LD:
746 #ifdef PARANOID
747  if ((s[0] != -1.0) || (s[1] != -1.0))
748  {
749  std::ostringstream error_stream;
750  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
751  << " is not on Bottom-Left edge\n";
752 
753  throw OomphLibError(error_stream.str(),
754  OOMPH_CURRENT_FUNCTION,
755  OOMPH_EXCEPTION_LOCATION);
756  }
757 #endif
758  increment1 = nnodes_2d;
759  break;
760 
761  case LU:
762 #ifdef PARANOID
763  if ((s[0] != -1.0) || (s[1] != 1.0))
764  {
765  std::ostringstream error_stream;
766  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
767  << " is not on Upper-Left edge\n";
768 
769  throw OomphLibError(error_stream.str(),
770  OOMPH_CURRENT_FUNCTION,
771  OOMPH_EXCEPTION_LOCATION);
772  }
773 #endif
774  start = nnodes_2d - nnodes_1d;
775  increment1 = nnodes_2d;
776  break;
777 
778  case LB:
779 #ifdef PARANOID
780  if ((s[0] != -1.0) || (s[2] != -1.0))
781  {
782  std::ostringstream error_stream;
783  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
784  << " is not on Back-Left edge\n";
785 
786  throw OomphLibError(error_stream.str(),
787  OOMPH_CURRENT_FUNCTION,
788  OOMPH_EXCEPTION_LOCATION);
789  }
790 #endif
791  increment1 = nnodes_1d;
792  break;
793 
794  case RF:
795 #ifdef PARANOID
796  if ((s[0] != 1.0) || (s[2] != 1.0))
797  {
798  std::ostringstream error_stream;
799  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
800  << " is not on Front-Right edge\n";
801 
802  throw OomphLibError(error_stream.str(),
803  OOMPH_CURRENT_FUNCTION,
804  OOMPH_EXCEPTION_LOCATION);
805  }
806 #endif
807  start = nnodes_3d - 1;
808  increment1 = -nnodes_1d;
809  break;
810 
811  case RD:
812 #ifdef PARANOID
813  if ((s[0] != 1.0) || (s[1] != -1.0))
814  {
815  std::ostringstream error_stream;
816  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
817  << " is not on Bottom-Right edge\n";
818 
819  throw OomphLibError(error_stream.str(),
820  OOMPH_CURRENT_FUNCTION,
821  OOMPH_EXCEPTION_LOCATION);
822  }
823 #endif
824  start = nnodes_1d - 1;
825  increment1 = nnodes_2d;
826  break;
827 
828  case RU:
829 #ifdef PARANOID
830  if ((s[0] != 1.0) || (s[1] != 1.0))
831  {
832  std::ostringstream error_stream;
833  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
834  << " is not on Upper-Right edge\n";
835 
836  throw OomphLibError(error_stream.str(),
837  OOMPH_CURRENT_FUNCTION,
838  OOMPH_EXCEPTION_LOCATION);
839  }
840 #endif
841  start = nnodes_2d - 1;
842  increment1 = nnodes_2d;
843  break;
844 
845  case RB:
846 #ifdef PARANOID
847  if ((s[0] != 1.0) || (s[2] != -1.0))
848  {
849  std::ostringstream error_stream;
850  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
851  << " is not on Back-Right edge\n";
852 
853  throw OomphLibError(error_stream.str(),
854  OOMPH_CURRENT_FUNCTION,
855  OOMPH_EXCEPTION_LOCATION);
856  }
857 #endif
858  start = nnodes_1d - 1;
859  increment1 = nnodes_1d;
860  break;
861 
862  case DB:
863 #ifdef PARANOID
864  if ((s[1] != -1.0) || (s[2] != -1.0))
865  {
866  std::ostringstream error_stream;
867  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
868  << " is not on Back-Bottom edge\n";
869 
870  throw OomphLibError(error_stream.str(),
871  OOMPH_CURRENT_FUNCTION,
872  OOMPH_EXCEPTION_LOCATION);
873  }
874 #endif
875  break;
876 
877  case DF:
878 #ifdef PARANOID
879  if ((s[1] != -1.0) || (s[2] != 1.0))
880  {
881  std::ostringstream error_stream;
882  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
883  << " is not on Front-Bottom edge\n";
884 
885  throw OomphLibError(error_stream.str(),
886  OOMPH_CURRENT_FUNCTION,
887  OOMPH_EXCEPTION_LOCATION);
888  }
889 #endif
890  start = nnodes_3d - nnodes_2d;
891  break;
892 
893  case UB:
894 #ifdef PARANOID
895  if ((s[1] != 1.0) || (s[2] != -1.0))
896  {
897  std::ostringstream error_stream;
898  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
899  << " is not on Back-Upper edge\n";
900 
901  throw OomphLibError(error_stream.str(),
902  OOMPH_CURRENT_FUNCTION,
903  OOMPH_EXCEPTION_LOCATION);
904  }
905 #endif
906  start = nnodes_2d - nnodes_1d;
907  break;
908 
909  case UF:
910 #ifdef PARANOID
911  if ((s[1] != 1.0) || (s[2] != 1.0))
912  {
913  std::ostringstream error_stream;
914  error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
915  << " is not on Upper-Front edge\n";
916 
917  throw OomphLibError(error_stream.str(),
918  OOMPH_CURRENT_FUNCTION,
919  OOMPH_EXCEPTION_LOCATION);
920  }
921 #endif
922  start = nnodes_3d - nnodes_1d;
923  break;
924 
925  default:
926 
927  std::ostringstream error_stream;
928  error_stream << "Wrong face " << OcTree::Direct_string[face]
929  << " passed" << std::endl;
930  error_stream << "Trouble at : s= [" << s[0] << " " << s[1] << " "
931  << s[2] << "]\n";
932  Vector<double> x(3);
933  this->interpolated_x(s, x);
934  error_stream << "corresponding to : x= [" << x[0] << " " << x[1] << " "
935  << x[2] << "]\n";
936  throw OomphLibError(
937  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
938  }
939 
940  // Initialise the intrinsic coordinate
941  zeta[0] = 0.0;
942  zeta[1] = 0.0;
943 
944  // Set the start node number
945  unsigned node = start;
946 
947  if (on_edge)
948  {
949  for (unsigned l1 = 0; l1 < nnodes_1d; l1++)
950  {
951  // Get the intrinsic coordinate
952  Vector<double> node_zeta(2);
953  node_pt(node)->get_coordinates_on_boundary(boundary, node_zeta);
954 
955  // Now multiply by the shape function
956  zeta[0] += node_zeta[0] * psi(node);
957  zeta[1] += node_zeta[1] * psi(node);
958 
959  // Update node
960  node += increment1;
961  }
962  }
963  else
964  {
965  for (unsigned l2 = 0; l2 < nnodes_1d; l2++)
966  {
967  for (unsigned l1 = 0; l1 < nnodes_1d; l1++)
968  {
969  // Get the intrinsic coordinate
970  Vector<double> node_zeta(2);
971  node_pt(node)->get_coordinates_on_boundary(boundary, node_zeta);
972 
973  // Now multiply by the shape function
974  zeta[0] += node_zeta[0] * psi(node);
975  zeta[1] += node_zeta[1] * psi(node);
976 
977  // Update node
978  node += increment1;
979  }
980  // Update node
981  node += increment2;
982  }
983  }
984  }
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::FiniteElement::shape
virtual void shape(const Vector< double > &s, Shape &psi) const =0
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::Node::get_coordinates_on_boundary
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Definition: nodes.cc:2379
oomph::OcTree::Direct_string
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition: octree.h:329
oomph::CumulativeTimings::start
void start(const unsigned &i)
(Re-)start i-th timer
Definition: oomph_utilities.cc:243
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::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::OcTree::Direct_string, oomph::OcTreeNames::F, L, oomph::OcTreeNames::LB, oomph::OcTreeNames::LD, oomph::OcTreeNames::LF, oomph::OcTreeNames::LU, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, oomph::OcTreeNames::RB, oomph::OcTreeNames::RD, oomph::OcTreeNames::RF, oomph::OcTreeNames::RU, s, oomph::OneDimLagrange::shape(), oomph::CumulativeTimings::start(), RachelsAdvectionDiffusion::U, oomph::OcTreeNames::UB, oomph::OcTreeNames::UF, plotDoE::x, and Eigen::zeta().

Referenced by build(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::p_refine(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::rebuild_from_sons(), and oomph::RefineableQSpectralElement< 3 >::rebuild_from_sons().

◆ node_created_by_neighbour()

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

994  {
995  using namespace OcTreeNames;
996 
997  // Calculate the faces/edges on which the node lies
998  Vector<int> faces;
999  Vector<int> edges;
1000 
1001  if (s_fraction[0] == 0.0)
1002  {
1003  faces.push_back(L);
1004  if (s_fraction[1] == 0.0)
1005  {
1006  edges.push_back(LD);
1007  }
1008  if (s_fraction[2] == 0.0)
1009  {
1010  edges.push_back(LB);
1011  }
1012  if (s_fraction[1] == 1.0)
1013  {
1014  edges.push_back(LU);
1015  }
1016  if (s_fraction[2] == 1.0)
1017  {
1018  edges.push_back(LF);
1019  }
1020  }
1021 
1022  if (s_fraction[0] == 1.0)
1023  {
1024  faces.push_back(R);
1025  if (s_fraction[1] == 0.0)
1026  {
1027  edges.push_back(RD);
1028  }
1029  if (s_fraction[2] == 0.0)
1030  {
1031  edges.push_back(RB);
1032  }
1033  if (s_fraction[1] == 1.0)
1034  {
1035  edges.push_back(RU);
1036  }
1037  if (s_fraction[2] == 1.0)
1038  {
1039  edges.push_back(RF);
1040  }
1041  }
1042 
1043  if (s_fraction[1] == 0.0)
1044  {
1045  faces.push_back(D);
1046  if (s_fraction[2] == 0.0)
1047  {
1048  edges.push_back(DB);
1049  }
1050  if (s_fraction[2] == 1.0)
1051  {
1052  edges.push_back(DF);
1053  }
1054  }
1055 
1056  if (s_fraction[1] == 1.0)
1057  {
1058  faces.push_back(U);
1059  if (s_fraction[2] == 0.0)
1060  {
1061  edges.push_back(UB);
1062  }
1063  if (s_fraction[2] == 1.0)
1064  {
1065  edges.push_back(UF);
1066  }
1067  }
1068 
1069  if (s_fraction[2] == 0.0)
1070  {
1071  faces.push_back(B);
1072  }
1073 
1074  if (s_fraction[2] == 1.0)
1075  {
1076  faces.push_back(F);
1077  }
1078 
1079  // Find the number of faces
1080  unsigned n_face = faces.size();
1081 
1082  // Find the number of edges
1083  unsigned n_edge = edges.size();
1084 
1085  Vector<unsigned> translate_s(3);
1086  Vector<double> s_lo_neigh(3);
1087  Vector<double> s_hi_neigh(3);
1088  Vector<double> s(3);
1089 
1090  int neigh_face, diff_level;
1091 
1092  // Loop over the faces on which the node lies
1093  //-------------------------------------------
1094  for (unsigned j = 0; j < n_face; j++)
1095  {
1096  // Boolean to indicate whether or not the neighbour has a
1097  // different Tree root
1098  bool in_neighbouring_tree;
1099 
1100  // Find pointer to neighbouring element along face
1101  OcTree* neigh_pt;
1102  neigh_pt = octree_pt()->gteq_face_neighbour(faces[j],
1103  translate_s,
1104  s_lo_neigh,
1105  s_hi_neigh,
1106  neigh_face,
1107  diff_level,
1108  in_neighbouring_tree);
1109 
1110  // Neighbour exists
1111  if (neigh_pt != 0)
1112  {
1113  // Have its nodes been created yet?
1114  if (neigh_pt->object_pt()->nodes_built())
1115  {
1116  // We now need to translate the nodal location, defined in terms
1117  // of the fractional coordinates of the present element into
1118  // those of its neighbour. For this we use the information returned
1119  // to use from the octree function.
1120 
1121  // Calculate the local coordinate in the neighbour
1122  // Note that we need to use the translation scheme in case
1123  // the local coordinates are swapped in the neighbour.
1124  for (unsigned i = 0; i < 3; i++)
1125  {
1126  s[i] = s_lo_neigh[i] +
1127  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
1128  }
1129 
1130  // Find the node in the neighbour
1131  Node* neighbour_node_pt =
1132  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
1133 
1134  // If there is a node, return it
1135  if (neighbour_node_pt != 0)
1136  {
1137  // Now work out whether it's a periodic boundary. This is
1138  // only possible if we have moved into a neighbouring tree
1139  if (in_neighbouring_tree)
1140  {
1141  // Return whether the neighbour is periodic
1142  is_periodic =
1143  octree_pt()->root_pt()->is_neighbour_periodic(faces[j]);
1144  }
1145 
1146  // Return the neighbour node pointer
1147  return neighbour_node_pt;
1148  }
1149  } // if (neigh_pt->object_pt()->nodes_built())
1150  } // if (neigh_pt!=0)
1151  } // for (unsigned j=0;j<n_face;j++)
1152 
1153  // Loop over the edges on which the node lies
1154  //------------------------------------------
1155  for (unsigned j = 0; j < n_edge; j++)
1156  {
1157  // Even if we restrict ourselves to true edge neighbours (i.e.
1158  // elements that are not also face neighbours) there may be multiple
1159  // edge neighbours across the edges between multiple root octrees.
1160  // When making the first call to OcTree::gteq_true_edge_neighbour(...)
1161  // we simply return the first one of these multiple edge neighbours
1162  // (if there are any at all, of course) and also return the total number
1163  // of true edge neighbours. If the node in question already exists
1164  // on the first edge neighbour we're done. If it doesn't it may exist
1165  // on other edge neighbours so we repeat the process over all
1166  // other edge neighbours (bailing out if a node is found, of course).
1167 
1168  // Initially return the zero-th true edge neighbour
1169  unsigned i_root_edge_neighbour = 0;
1170 
1171  // Initialise the total number of true edge neighbours
1172  unsigned nroot_edge_neighbour = 0;
1173 
1174  // Keep searching until we've found the node or until we've checked
1175  // all available edge neighbours
1176  bool keep_searching = true;
1177  while (keep_searching)
1178  {
1179  // Find pointer to neighbouring element along edge
1180  OcTree* neigh_pt;
1181  neigh_pt = octree_pt()->gteq_true_edge_neighbour(edges[j],
1182  i_root_edge_neighbour,
1183  nroot_edge_neighbour,
1184  translate_s,
1185  s_lo_neigh,
1186  s_hi_neigh,
1187  neigh_face,
1188  diff_level);
1189 
1190  // Neighbour exists
1191  if (neigh_pt != 0)
1192  {
1193  // Have its nodes been created yet?
1194  if (neigh_pt->object_pt()->nodes_built())
1195  {
1196  // We now need to translate the nodal location, defined in terms
1197  // of the fractional coordinates of the present element into
1198  // those of its neighbour. For this we use the information returned
1199  // to use from the octree function.
1200 
1201  // Calculate the local coordinate in the neighbour
1202  // Note that we need to use the translation scheme in case
1203  // the local coordinates are swapped in the neighbour.
1204  for (unsigned i = 0; i < 3; i++)
1205  {
1206  s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
1207  (s_hi_neigh[i] - s_lo_neigh[i]);
1208  }
1209 
1210  // Find the node in the neighbour
1211  Node* neighbour_node_pt =
1212  neigh_pt->object_pt()->get_node_at_local_coordinate(s);
1213 
1214  // If there is a node, return it
1215  if (neighbour_node_pt != 0)
1216  {
1217  // Get the faces on which the edge lies
1218  Vector<int> faces_attached_to_edge =
1219  OcTree::faces_of_common_edge(edges[j]);
1220 
1221  // Get the number of entries in the vector
1222  unsigned n_faces_attached_to_edge = faces_attached_to_edge.size();
1223 
1224  // Loop over the faces
1225  for (unsigned i_face = 0; i_face < n_faces_attached_to_edge;
1226  i_face++)
1227  {
1228  // Is the node periodic in the face direction?
1229  is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
1230  faces_attached_to_edge[i_face]);
1231 
1232  // Check if the edge is periodic in the i_face-th face direction
1233  if (is_periodic)
1234  {
1235  // We're done!
1236  break;
1237  }
1238  } // for (unsigned
1239  // i_face=0;i_face<n_faces_attached_to_edge;i_face++)
1240 
1241  // Return the neighbour node pointer
1242  return neighbour_node_pt;
1243  } // if (neighbour_node_pt!=0)
1244  } // if (neigh_pt->object_pt()->nodes_built())
1245  } // if (neigh_pt!=0)
1246 
1247  // Keep searching, but only if there are further edge neighbours
1248  // Try next root edge neighbour
1249  i_root_edge_neighbour++;
1250 
1251  // Have we exhausted the search?
1252  if (i_root_edge_neighbour >= nroot_edge_neighbour)
1253  {
1254  // Stop searching
1255  keep_searching = false;
1256  }
1257  } // End of while keep searching over all true edge neighbours
1258  } // End of loop over edges
1259 
1260  // Node not found, return null
1261  return 0;
1262  }
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_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
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

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::RefineableElement::nodes_built(), oomph::Tree::object_pt(), 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()

virtual Node* oomph::RefineableQElement< 3 >::node_created_by_son_of_neighbour ( const Vector< double > &  s_fraction,
bool &  is_periodic 
)
inlinevirtual

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 in oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >.

120  {
121  // It is impossible for this situation to arise in meshes
122  // containing elements of uniform p-order. This is here so
123  // that it can be overloaded for p-refineable elements.
124  return 0;
125  }

◆ oc_hang_helper()

void oomph::RefineableQElement< 3 >::oc_hang_helper ( const int &  value_id,
const int &  my_face,
std::ofstream &  output_hangfile 
)
protectedvirtual

Internal helper function that is used to construct the hanging node schemes for the positions.

Internal function to set up the hanging nodes on a particular face of the element

Reimplemented in oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >.

2145  {
2146  using namespace OcTreeNames;
2147 
2148  // Number of dimensions
2149  unsigned n_dim = 3;
2150 
2151  Vector<unsigned> translate_s(n_dim);
2152  Vector<double> s_lo_neigh(n_dim);
2153  Vector<double> s_hi_neigh(n_dim);
2154  int neigh_face;
2155  int diff_level;
2156  bool in_neighbouring_tree;
2157 
2158  // Find pointer to neighbour in this direction
2159  OcTree* neigh_pt;
2160  neigh_pt = octree_pt()->gteq_face_neighbour(my_face,
2161  translate_s,
2162  s_lo_neigh,
2163  s_hi_neigh,
2164  neigh_face,
2165  diff_level,
2166  in_neighbouring_tree);
2167 
2168  // Neighbour exists
2169  if (neigh_pt != 0)
2170  {
2171  // Different sized element?
2172  if (diff_level != 0)
2173  {
2174  // Test for the periodic node case
2175  // Are we crossing a periodic boundary
2176  bool is_periodic = false;
2177  if (in_neighbouring_tree)
2178  {
2179  is_periodic = tree_pt()->root_pt()->is_neighbour_periodic(my_face);
2180  }
2181 
2182  // If it is periodic we actually need to get the node in
2183  // the neighbour of the neighbour (which will be a parent of
2184  // the present element) so that the "fixed" coordinate is
2185  // correctly calculated.
2186  // The idea is to replace the neigh_pt and associated data
2187  // with those of the neighbour of the neighbour
2188  if (is_periodic)
2189  {
2190  // Required data for the neighbour finding routine
2191  Vector<unsigned> translate_s_in_neigh(n_dim, 0.0);
2192  Vector<double> s_lo_neigh_of_neigh(n_dim, 0.0);
2193  Vector<double> s_hi_neigh_of_neigh(n_dim, 0.0);
2194  int neigh_face_of_neigh;
2195  int diff_level_of_neigh;
2196  bool in_neighbouring_tree_of_neigh;
2197 
2198  // Find pointer to neighbour of the neighbour on the edge
2199  // that we are currently considering
2200  OcTree* neigh_of_neigh_pt;
2201  neigh_of_neigh_pt =
2202  neigh_pt->gteq_face_neighbour(neigh_face,
2203  translate_s_in_neigh,
2204  s_lo_neigh_of_neigh,
2205  s_hi_neigh_of_neigh,
2206  neigh_face_of_neigh,
2207  diff_level_of_neigh,
2208  in_neighbouring_tree_of_neigh);
2209 
2210  // Set the value of the NEW neighbour and edge
2211  neigh_pt = neigh_of_neigh_pt;
2212  neigh_face = neigh_face_of_neigh;
2213 
2214  // Set the values of the translation schemes
2215  // Need to find the values of s_lo and s_hi
2216  // in the neighbour of the neighbour
2217 
2218  // Get the minimum and maximum values of the coordinate
2219  // in the neighbour (don't like this, but I think it's
2220  // necessary) Note that these values are hardcoded
2221  // in the quadtrees at some point!!
2222  double s_min = neigh_pt->object_pt()->s_min();
2223  double s_max = neigh_pt->object_pt()->s_max();
2224  Vector<double> s_lo_frac(n_dim), s_hi_frac(n_dim);
2225  // Work out the fractional position of the low and high points
2226  // of the original element
2227  for (unsigned i = 0; i < n_dim; i++)
2228  {
2229  s_lo_frac[i] = (s_lo_neigh[i] - s_min) / (s_max - s_min);
2230  s_hi_frac[i] = (s_hi_neigh[i] - s_min) / (s_max - s_min);
2231  }
2232 
2233  // We should now be able to construct the low and high points in
2234  // the neighbour of the neighbour
2235  for (unsigned i = 0; i < n_dim; i++)
2236  {
2237  s_lo_neigh[i] = s_lo_neigh_of_neigh[i] +
2238  s_lo_frac[translate_s_in_neigh[i]] *
2239  (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
2240  s_hi_neigh[i] = s_lo_neigh_of_neigh[i] +
2241  s_hi_frac[translate_s_in_neigh[i]] *
2242  (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
2243  }
2244 
2245  // Finally we must sort out the translation scheme
2246  Vector<unsigned> temp_translate(n_dim, 0.0);
2247  for (unsigned i = 0; i < n_dim; i++)
2248  {
2249  temp_translate[i] = translate_s_in_neigh[translate_s[i]];
2250  }
2251  for (unsigned i = 0; i < n_dim; i++)
2252  {
2253  translate_s[i] = temp_translate[i];
2254  }
2255  } // End of special treatment for periodic hanging nodes
2256 
2257  // Number of nodes in one dimension
2258  unsigned n_p = ninterpolating_node_1d(value_id);
2259  // Storage for the local nodes along the face of the element
2260  Node* local_node_pt = 0;
2261 
2262  // Loop over nodes along the face
2263  for (unsigned i0 = 0; i0 < n_p; i0++)
2264  {
2265  for (unsigned i1 = 0; i1 < n_p; i1++)
2266  {
2267  // Storage for the fractional position of the node in the element
2268  Vector<double> s_fraction(n_dim);
2269 
2270  // Local node number
2271  switch (my_face)
2272  {
2273  case U:
2274  s_fraction[0] =
2275  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2276  s_fraction[1] = 1.0;
2277  s_fraction[2] =
2278  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2279  local_node_pt = interpolating_node_pt(
2280  i0 + (n_p - 1) * n_p + n_p * n_p * i1, value_id);
2281  break;
2282 
2283  case D:
2284  s_fraction[0] =
2285  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2286  s_fraction[1] = 0.0;
2287  s_fraction[2] =
2288  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2289  local_node_pt =
2290  interpolating_node_pt(i0 + n_p * n_p * i1, value_id);
2291  break;
2292 
2293  case R:
2294  s_fraction[0] = 1.0;
2295  s_fraction[1] =
2296  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
2297  s_fraction[2] =
2298  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2299  local_node_pt = interpolating_node_pt(
2300  n_p - 1 + i0 * n_p + i1 * n_p * n_p, value_id);
2301  break;
2302 
2303  case L:
2304  s_fraction[0] = 0.0;
2305  s_fraction[1] =
2306  local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
2307  s_fraction[2] =
2308  local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2309  local_node_pt =
2310  interpolating_node_pt(n_p * i0 + i1 * n_p * n_p, value_id);
2311  break;
2312 
2313  case B:
2314  s_fraction[0] =
2315  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2316  s_fraction[1] =
2317  local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
2318  s_fraction[2] = 0.0;
2319  local_node_pt = interpolating_node_pt(i0 + i1 * n_p, value_id);
2320  break;
2321 
2322  case F:
2323  s_fraction[0] =
2324  local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2325  s_fraction[1] =
2326  local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
2327  s_fraction[2] = 1.0;
2328  local_node_pt = interpolating_node_pt(
2329  i0 + i1 * n_p + (n_p - 1) * n_p * n_p, value_id);
2330  break;
2331 
2332  default:
2333  throw OomphLibError("my_face not U, D, L, R, B, F\n",
2334  OOMPH_CURRENT_FUNCTION,
2335  OOMPH_EXCEPTION_LOCATION);
2336  }
2337 
2338  // Set up the coordinate in the neighbour
2339  Vector<double> s_in_neighb(n_dim);
2340  for (unsigned i = 0; i < n_dim; i++)
2341  {
2342  s_in_neighb[i] =
2343  s_lo_neigh[i] +
2344  s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
2345  }
2346 
2347  // Find the Node in the neighbouring element
2348  Node* neighbouring_node_pt =
2349  neigh_pt->object_pt()->get_interpolating_node_at_local_coordinate(
2350  s_in_neighb, value_id);
2351 
2352  // If the neighbour does not have a node at this point
2353  if (0 == neighbouring_node_pt)
2354  {
2355  // Do we need to make a hanging node, we assume that we don't
2356  // initially
2357  bool make_hanging_node = false;
2358 
2359  // If the node is not hanging geometrically, then we must make it
2360  // hang
2361  if (!local_node_pt->is_hanging())
2362  {
2363  make_hanging_node = true;
2364  }
2365  // Otherwise is could be hanging geometrically, but still require
2366  // a different hanging scheme for this data value
2367  else
2368  {
2369  if ((value_id != -1) && (local_node_pt->hanging_pt(value_id) ==
2370  local_node_pt->hanging_pt()))
2371  {
2372  make_hanging_node = true;
2373  }
2374  }
2375 
2376  if (make_hanging_node == true)
2377  {
2378  // Cache refineable element used here
2379  RefineableElement* const obj_pt = neigh_pt->object_pt();
2380 
2381  // Get shape functions in neighbour element
2382  Shape psi(obj_pt->ninterpolating_node(value_id));
2383  obj_pt->interpolating_basis(s_in_neighb, psi, value_id);
2384 
2385  // Allocate the storage for the Hang pointer
2386  // We know that it will hold n_p*n_p nodes
2387  HangInfo* hang_pt = new HangInfo(n_p * n_p);
2388 
2389  // Loop over nodes on edge in neighbour and mark them as nodes
2390  // that this node depends on
2391  unsigned n_neighbour;
2392 
2393  // Number of nodes along edge in neighbour element
2394  for (unsigned ii0 = 0; ii0 < n_p; ii0++)
2395  {
2396  for (unsigned ii1 = 0; ii1 < n_p; ii1++)
2397  {
2398  switch (neigh_face)
2399  {
2400  case U:
2401  n_neighbour = ii0 + n_p * (n_p - 1) + ii1 * n_p * n_p;
2402  break;
2403 
2404  case D:
2405  n_neighbour = ii0 + ii1 * n_p * n_p;
2406  break;
2407 
2408  case L:
2409  n_neighbour = ii0 * n_p + ii1 * n_p * n_p;
2410  break;
2411 
2412  case R:
2413  n_neighbour = (n_p - 1) + ii0 * n_p + ii1 * n_p * n_p;
2414  break;
2415 
2416  case B:
2417  n_neighbour = ii0 + ii1 * n_p;
2418  break;
2419 
2420  case F:
2421  n_neighbour = ii0 + ii1 * n_p + n_p * n_p * (n_p - 1);
2422  break;
2423  default:
2424  throw OomphLibError("neigh_face not U, L, R, B, F\n",
2425  OOMPH_CURRENT_FUNCTION,
2426  OOMPH_EXCEPTION_LOCATION);
2427  }
2428 
2429  // Push back neighbouring node and weight into
2430  // Vector of (pointers to)
2431  // master nodes and weights
2432  // The weight is merely the value of the shape function
2433  // corresponding to the node in the neighbour
2434  hang_pt->set_master_node_pt(
2435  ii0 * n_p + ii1,
2436  obj_pt->interpolating_node_pt(n_neighbour, value_id),
2437  psi[n_neighbour]);
2438  }
2439  }
2440  // Now set the hanging data for the position
2441  // This also constrains the data values associated with the
2442  // value id
2443  local_node_pt->set_hanging_pt(hang_pt, value_id);
2444  }
2445 
2446  if (output_hangfile.is_open())
2447  {
2448  // output_hangfile
2449  output_hangfile << local_node_pt->x(0) << " "
2450  << local_node_pt->x(1) << " "
2451  << local_node_pt->x(2) << std::endl;
2452  }
2453  }
2454  else
2455  {
2456 #ifdef PARANOID
2457  if (local_node_pt != neighbouring_node_pt)
2458  {
2459  std::ofstream reportage("dodgy.dat", std::ios_base::app);
2460  reportage << local_node_pt->x(0) << " " << local_node_pt->x(1)
2461  << " " << local_node_pt->x(2) << std::endl;
2462  reportage.close();
2463 
2464  std::ostringstream warning_stream;
2465  warning_stream
2466  << "SANITY CHECK in oc_hang_helper \n"
2467  << "Current node " << local_node_pt << " at "
2468  << "(" << local_node_pt->x(0) << ", " << local_node_pt->x(1)
2469  << ", " << local_node_pt->x(2) << ")" << std::endl
2470  << " is not hanging and has " << std::endl
2471  << "Neighbour's node " << neighbouring_node_pt << " at "
2472  << "(" << neighbouring_node_pt->x(0) << ", "
2473  << neighbouring_node_pt->x(1) << ", "
2474  << neighbouring_node_pt->x(2) << ")" << std::endl
2475  << "even though the two should be "
2476  << "identical" << std::endl;
2477  OomphLibWarning(warning_stream.str(),
2478  OOMPH_CURRENT_FUNCTION,
2479  OOMPH_EXCEPTION_LOCATION);
2480  }
2481 #endif
2482  }
2483 
2484  // If we are doing the position then
2485  if (value_id == -1)
2486  {
2487  // Get the nodal position from neighbour element
2488  Vector<double> x_in_neighb(n_dim);
2489  neigh_pt->object_pt()->interpolated_x(s_in_neighb, x_in_neighb);
2490 
2491  // Fine adjust the coordinates (macro map will pick up boundary
2492  // accurately but will lead to different element edges)
2493  local_node_pt->x(0) = x_in_neighb[0];
2494  local_node_pt->x(1) = x_in_neighb[1];
2495  local_node_pt->x(2) = x_in_neighb[2];
2496  }
2497  }
2498  }
2499  }
2500  }
2501  }
oomph::FiniteElement::s_min
virtual double s_min() const
Min value of local coordinate.
Definition: elements.h:2793
oomph::FiniteElement::s_max
virtual double s_max() const
Max. value of local coordinate.
Definition: elements.h:2803
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

References D, oomph::OcTreeNames::F, oomph::RefineableElement::get_interpolating_node_at_local_coordinate(), oomph::OcTree::gteq_face_neighbour(), oomph::Node::hanging_pt(), i, oomph::FiniteElement::interpolated_x(), oomph::RefineableElement::interpolating_basis(), oomph::RefineableElement::interpolating_node_pt(), oomph::Node::is_hanging(), L, oomph::RefineableElement::ninterpolating_node(), oomph::Tree::object_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, oomph::FiniteElement::s_max(), oomph::FiniteElement::s_min(), oomph::Node::set_hanging_pt(), oomph::HangInfo::set_master_node_pt(), RachelsAdvectionDiffusion::U, and oomph::Node::x().

◆ octree_pt() [1/2]

OcTree* oomph::RefineableQElement< 3 >::octree_pt ( )
inline

Pointer to octree representation of this element.

145  {
146  return dynamic_cast<OcTree*>(Tree_pt);
147  }
oomph::RefineableElement::Tree_pt
Tree * Tree_pt
A pointer to a general tree object.
Definition: refineable_elements.h:100

◆ octree_pt() [2/2]

OcTree* oomph::RefineableQElement< 3 >::octree_pt ( ) const
inline

Pointer to octree representation of this element.

151  {
152  return dynamic_cast<OcTree*>(Tree_pt);
153  }

◆ output_corners()

void oomph::RefineableQElement< 3 >::output_corners ( std::ostream &  outfile,
const std::string &  colour 
) const

Print corner nodes, use colour.

Print corner nodes, use colour (default "BLACK") in right order so that tecplot can draw a cube without crossed lines

42  {
43  Vector<double> s(3);
44  Vector<double> corner(3);
45 
46  outfile << "ZONE I=2, J=2, K=2 C=" << colour << std::endl;
47 
48  s[0] = -1.0;
49  s[1] = -1.0;
50  s[2] = -1.0;
51  get_x(s, corner);
52  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
53  << Number << std::endl;
54 
55  s[0] = -1.0;
56  s[1] = -1.0;
57  s[2] = 1.0;
58  get_x(s, corner);
59  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
60  << Number << std::endl;
61 
62  s[0] = -1.0;
63  s[1] = 1.0;
64  s[2] = -1.0;
65  get_x(s, corner);
66  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
67  << Number << std::endl;
68 
69  s[0] = -1.0;
70  s[1] = 1.0;
71  s[2] = 1.0;
72  get_x(s, corner);
73  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
74  << Number << std::endl;
75 
76  // next face
77 
78 
79  s[0] = 1.0;
80  s[1] = -1.0;
81  s[2] = -1.0;
82  get_x(s, corner);
83  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
84  << Number << std::endl;
85 
86  s[0] = 1.0;
87  s[1] = -1.0;
88  s[2] = 1.0;
89  get_x(s, corner);
90  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
91  << Number << std::endl;
92 
93  s[0] = 1.0;
94  s[1] = 1.0;
95  s[2] = -1.0;
96  get_x(s, corner);
97  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
98  << Number << std::endl;
99 
100  s[0] = 1.0;
101  s[1] = 1.0;
102  s[2] = 1.0;
103  get_x(s, corner);
104  outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
105  << Number << std::endl;
106 
107 
108  // outfile << "TEXT CS = GRID, X = " << corner[0] <<
109  // ",Y = " << corner[1] << ",Z = " << corner[2] <<
110  // ", HU = GRID, H = 0.01, AN = MIDCENTER, T =\""
111  // << Number << "\"" << std::endl;
112  }
oomph::RefineableElement::Number
long Number
Global element number – for plotting/validation purposes.
Definition: refineable_elements.h:115
oomph::TecplotNames::colour
Vector< std::string > colour
Tecplot colours.
Definition: oomph_utilities.cc:1159

References oomph::TecplotNames::colour, oomph::Global_unsigned::Number, and s.

Referenced by TestPoissonProblem< ELEMENT >::doc_solution().

◆ required_nsons()

unsigned oomph::RefineableQElement< 3 >::required_nsons ( ) const
inlinevirtual

A refineable brick element has eight sons.

Reimplemented from oomph::RefineableElement.

103  {
104  return 8;
105  }

◆ setup_father_bounds()

void oomph::RefineableQElement< 3 >::setup_father_bounds ( )
protected

Setup static matrix for coincidence between son nodal points and father boundaries

Setup static matrix for coincidence between son nodal points and father boundaries:

Father_bound[nnode_1d](nnode_son,son_type)={RU/RF/RD/RB/.../ OMEGA}

so that node nnode_son in element of type son_type lies on boundary/vertex Father_bound[nnode_1d](nnode_son,son_type) in its father element. If the node doesn't lie on a boundary the value is OMEGA.

127  {
128  using namespace OcTreeNames;
129 
130  // Find the number of nodes along a 1D edge
131  unsigned n_p = nnode_1d();
132 
133  // Allocate space for the boundary information
134  Father_bound[n_p].resize(n_p * n_p * n_p, 8);
135 
136  // Initialise: By default points are not on the boundary
137  for (unsigned n = 0; n < n_p * n_p * n_p; n++)
138  {
139  for (unsigned ison = 0; ison < 8; ison++)
140  {
141  Father_bound[n_p](n, ison) = Tree::OMEGA;
142  }
143  }
144 
145  for (int i_son = LDB; i_son <= RUF; i_son++)
146  {
147  // vector representing the son
148  Vector<int> vect_son(3);
149  // vector representing (at the end) the boundaries
150  Vector<int> vect_bound(3);
151  vect_son = OcTree::Direction_to_vector[i_son];
152  for (unsigned i0 = 0; i0 < n_p; i0++)
153  {
154  for (unsigned i1 = 0; i1 < n_p; i1++)
155  {
156  for (unsigned i2 = 0; i2 < n_p; i2++)
157  {
158  // Initialisation to make it work
159  for (unsigned i = 0; i < 3; i++)
160  {
161  vect_bound[i] = -vect_son[i];
162  }
163  // Seting up the boundaries coordinates as if the coordinates
164  // of vect_bound had been initialised to 0 if it were so,
165  // vect_bound would be the vector of the boundaries in the son
166  // itself.
167 
168  if (i0 == 0)
169  {
170  vect_bound[0] = -1;
171  }
172  if (i0 == n_p - 1)
173  {
174  vect_bound[0] = 1;
175  }
176  if (i1 == 0)
177  {
178  vect_bound[1] = -1;
179  }
180  if (i1 == n_p - 1)
181  {
182  vect_bound[1] = 1;
183  }
184  if (i2 == 0)
185  {
186  vect_bound[2] = -1;
187  }
188  if (i2 == n_p - 1)
189  {
190  vect_bound[2] = 1;
191  }
192 
193  // The effect of this is to filter the boundaries to keep only the
194  // ones which are effectively father boundaries.
195  // -- if the node is not on a "i0 boundary", we still
196  // have vect_bound[0]=-vect_son[0]
197  // and the result is vect_bound[0]=0
198  // (he is not on this boundary for his father)
199  // -- if he is on a son's boundary which is not one of
200  // the father -> same thing
201  // -- if he is on a boundary which is one of his father's,
202  // vect_bound[i]=vect_son[i]
203  // and the new vect_bound[i] is the same as the old one
204  for (int i = 0; i < 3; i++)
205  {
206  vect_bound[i] = (vect_bound[i] + vect_son[i]) / 2;
207  }
208 
209  // Return the result as {U,R,D,...RDB,LUF,OMEGA}
210  Father_bound[n_p](i0 + n_p * i1 + n_p * n_p * i2, i_son) =
211  OcTree::Vector_to_direction[vect_bound];
212 
213  } // Loop over i2
214  } // Loop over i1
215  } // Loop over i0
216  } // Loop over i_son
217 
218  } // setup_father_bounds()
oomph::OcTreeNames::RUF
@ RUF
Definition: octree.h:56
oomph::OcTreeNames::LDB
@ LDB
Definition: octree.h:49

References oomph::OcTree::Direction_to_vector, i, oomph::OcTreeNames::LDB, n, oomph::Tree::OMEGA, oomph::OcTreeNames::RUF, and oomph::OcTree::Vector_to_direction.

◆ setup_hang_for_value()

void oomph::RefineableQElement< 3 >::setup_hang_for_value ( const int &  value_id)
protected

Internal helper function that is used to construct the hanging node schemes for the value_id-th interpolated value

Internal function that sets up the hanging node scheme for a particular value

2125  {
2126  using namespace OcTreeNames;
2127 
2128  std::ofstream dummy_hangfile;
2129  // Setup hanging nodes on each edge of the element
2130  oc_hang_helper(value_id, U, dummy_hangfile);
2131  oc_hang_helper(value_id, D, dummy_hangfile);
2132  oc_hang_helper(value_id, R, dummy_hangfile);
2133  oc_hang_helper(value_id, L, dummy_hangfile);
2134  oc_hang_helper(value_id, B, dummy_hangfile);
2135  oc_hang_helper(value_id, F, dummy_hangfile);
2136  }
oomph::RefineableQElement< 3 >::oc_hang_helper
virtual void oc_hang_helper(const int &value_id, const int &my_edge, std::ofstream &output_hangfile)
Definition: refineable_brick_element.cc:2142

References D, oomph::OcTreeNames::F, L, R, and RachelsAdvectionDiffusion::U.

◆ setup_hanging_nodes()

void oomph::RefineableQElement< 3 >::setup_hanging_nodes ( Vector< std::ofstream * > &  output_stream)
virtual

Markup all hanging nodes & document the results in the output streams contained in the vector output_stream, if they are open.

Set up all hanging nodes.

(Wrapper to avoid specification of the unwanted output file).

Reimplemented from oomph::RefineableElement.

2099  {
2100 #ifdef PARANOID
2101  if (output_stream.size() != 6)
2102  {
2103  throw OomphLibError("There must be six output streams",
2104  OOMPH_CURRENT_FUNCTION,
2105  OOMPH_EXCEPTION_LOCATION);
2106  }
2107 #endif
2108 
2109  using namespace OcTreeNames;
2110 
2111  // Setup hanging nodes on each edge of the element
2112  oc_hang_helper(-1, U, *(output_stream[0]));
2113  oc_hang_helper(-1, D, *(output_stream[1]));
2114  oc_hang_helper(-1, L, *(output_stream[2]));
2115  oc_hang_helper(-1, R, *(output_stream[3]));
2116  oc_hang_helper(-1, B, *(output_stream[4]));
2117  oc_hang_helper(-1, F, *(output_stream[5]));
2118  }

References D, oomph::OcTreeNames::F, L, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, R, and RachelsAdvectionDiffusion::U.

Member Data Documentation

◆ Father_bound

std::map< unsigned, DenseMatrix< int > > oomph::RefineableQElement< 3 >::Father_bound
staticprotected

Coincidence between son nodal points and father boundaries: Father_bound[nnode_1d](nnode_son,son_type)={RU/RF/RD/RB/.../ OMEGA} so that node nnode_son in element of type son_type lies on boundary/vertex Father_bound[nnode_1d](nnode_son,son_type) in its father element. If the node doesn't lie on a boundary the value is OMEGA.

Static matrix for coincidence between son nodal points and father boundaries

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