Public Member Functions | Static Public Member Functions | Protected Member Functions | List of all members
oomph::RefineablePolarNavierStokesEquations Class Referenceabstract

#include <refineable_polar_navier_stokes_elements.h>

+ Inheritance diagram for oomph::RefineablePolarNavierStokesEquations:

Public Member Functions

 RefineablePolarNavierStokesEquations ()
 Constructor. More...
 
unsigned num_Z2_flux_terms ()
 Number of 'flux' terms for Z2 error estimation. More...
 
void get_Z2_flux (const Vector< double > &s, Vector< double > &flux)
 
void further_build ()
 Further build, pass the pointers down to the sons. More...
 
- Public Member Functions inherited from oomph::PolarNavierStokesEquations
 PolarNavierStokesEquations ()
 Constructor: NULL the body force and source function. More...
 
const double & re () const
 Reynolds number. More...
 
const double & alpha () const
 Alpha. More...
 
const double & re_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& alpha_pt ()
 Pointer to Alpha. More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
const double & viscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to Viscosity Ratio. More...
 
const double & density_ratio () const
 
double *& density_ratio_pt ()
 Pointer to Density ratio. More...
 
const double & re_invfr () const
 Global inverse Froude number. More...
 
double *& re_invfr_pt ()
 Pointer to global inverse Froude number. More...
 
const Vector< double > & g () const
 Vector of gravitational components. More...
 
Vector< double > *& g_pt ()
 Pointer to Vector of gravitational components. More...
 
NavierStokesBodyForceFctPt & body_force_fct_pt ()
 Access function for the body-force pointer. More...
 
NavierStokesBodyForceFctPt body_force_fct_pt () const
 Access function for the body-force pointer. Const version. More...
 
NavierStokesSourceFctPt & source_fct_pt ()
 Access function for the source-function pointer. More...
 
NavierStokesSourceFctPt source_fct_pt () const
 Access function for the source-function pointer. Const version. More...
 
virtual unsigned npres_pnst () const =0
 Function to return number of pressure degrees of freedom. More...
 
virtual double u_pnst (const unsigned &n, const unsigned &i) const =0
 
virtual double u_pnst (const unsigned &t, const unsigned &n, const unsigned &i) const =0
 
virtual unsigned u_index_pnst (const unsigned &i) const
 
double du_dt_pnst (const unsigned &n, const unsigned &i) const
 
virtual double p_pnst (const unsigned &n_p) const =0
 
virtual void fix_pressure (const unsigned &p_dof, const double &p_value)=0
 Pin p_dof-th pressure dof and set it to value specified by p_value. More...
 
virtual int p_nodal_index_pnst ()
 
double pressure_integral () const
 Integral of pressure over element. More...
 
double dissipation () const
 Return integral of dissipation over element. More...
 
double dissipation (const Vector< double > &s) const
 Return dissipation at local coordinate s. More...
 
double kin_energy () const
 Get integral of kinetic energy over element. More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 Strain-rate tensor: Now returns polar strain. More...
 
void strain_rate_by_r (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 Function to return polar strain multiplied by r. More...
 
void get_traction (const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
 
void get_load (const Vector< double > &s, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
 
void output (std::ostream &outfile)
 
void output (std::ostream &outfile, const unsigned &nplot)
 
void output (FILE *file_pt)
 
void output (FILE *file_pt, const unsigned &nplot)
 
void full_output (std::ostream &outfile)
 
void full_output (std::ostream &outfile, const unsigned &nplot)
 
void output_veloc (std::ostream &outfile, const unsigned &nplot, const unsigned &t)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 
void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
void fill_in_contribution_to_residuals (Vector< double > &residuals)
 Compute the element's residual Vector. More...
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void fill_in_contribution_to_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
void interpolated_u_pnst (const Vector< double > &s, Vector< double > &veloc) const
 Compute vector of FE interpolated velocity u at local coordinate s. More...
 
double interpolated_u_pnst (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated velocity u[i] at local coordinate s. More...
 
double interpolated_p_pnst (const Vector< double > &s) const
 Return FE interpolated pressure at local coordinate s. More...
 
double interpolated_dudx_pnst (const Vector< double > &s, const unsigned &i, const unsigned &j) const
 
- 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...
 
virtual bool local_coord_is_valid (const Vector< double > &s)
 Broken assignment operator. More...
 
virtual void move_local_coord_back_into_element (Vector< double > &s) const
 
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)
 
virtual void set_macro_elem_pt (MacroElement *macro_elem_pt)
 
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 get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (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 ElementGeometry::ElementGeometry element_geometry () const
 Return the geometry type of the element (either Q or T usually). More...
 
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 (const unsigned &t, std::ostream &outfile, const unsigned &n_plot) const
 
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, 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...
 
virtual unsigned nnode_on_face () const
 
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::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...
 
virtual unsigned required_nsons () const
 
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)
 
virtual void build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)=0
 
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
 
virtual void check_integrity (double &max_error)=0
 
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 setup_hanging_nodes (Vector< std::ofstream * > &output_stream)
 
virtual void further_setup_hanging_nodes ()
 
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::ElementWithZ2ErrorEstimator
 ElementWithZ2ErrorEstimator ()
 Default empty constructor. More...
 
 ElementWithZ2ErrorEstimator (const ElementWithZ2ErrorEstimator &)=delete
 Broken copy constructor. More...
 
void operator= (const ElementWithZ2ErrorEstimator &)=delete
 Broken assignment operator. More...
 
virtual unsigned ncompound_fluxes ()
 
virtual void compute_exact_Z2_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
 
virtual void get_Z2_compound_flux_indices (Vector< unsigned > &flux_index)
 
virtual unsigned nvertex_node () const =0
 Number of vertex nodes in the element. More...
 
virtual Node * vertex_node_pt (const unsigned &j) const =0
 
virtual unsigned nrecovery_order ()=0
 Order of recovery shape functions. More...
 
virtual double geometric_jacobian (const Vector< double > &x)
 

Static Public Member Functions

static void pin_redundant_nodal_pressures (const Vector< GeneralisedElement * > &element_pt)
 
static void unpin_all_pressure_dofs (const Vector< GeneralisedElement * > &element_pt)
 Unpin all pressure dofs in elements listed in vector. More...
 
- Static Public Member Functions inherited from oomph::RefineableElement
static double & max_integrity_tolerance ()
 Max. allowed discrepancy in element integrity check. More...
 

Protected Member Functions

virtual Node * pressure_node_pt (const unsigned &n_p)
 
virtual void unpin_elemental_pressure_dofs ()=0
 Unpin all pressure dofs in the element. More...
 
virtual void pin_elemental_redundant_nodal_pressure_dofs ()
 
virtual void fill_in_generic_residual_contribution (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 Start of what would've been refineable_navier_stokes_elements.cc //. More...
 
- Protected Member Functions inherited from oomph::PolarNavierStokesEquations
virtual int p_local_eqn (const unsigned &n)=0
 
virtual double dshape_and_dtest_eulerian_pnst (const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_pnst (const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual void pshape_pnst (const Vector< double > &s, Shape &psi) const =0
 Compute the pressure shape functions at local coordinate s. More...
 
virtual void pshape_pnst (const Vector< double > &s, Shape &psi, Shape &test) const =0
 
void get_body_force (double time, const Vector< double > &x, Vector< double > &result)
 Calculate the body force at a given time and Eulerian position. More...
 
double get_source_fct (double time, const Vector< double > &x)
 
- 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)
 
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)
 
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_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)
 
- 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)
 

Additional Inherited Members

- Public Types inherited from oomph::PolarNavierStokesEquations
typedef void(* NavierStokesBodyForceFctPt) (const double &time, const Vector< double > &x, Vector< double > &body_force)
 
typedef double(* NavierStokesSourceFctPt) (const double &time, const Vector< double > &x)
 
- Public Types inherited from oomph::FiniteElement
typedef void(* SteadyExactSolutionFctPt) (const Vector< double > &, Vector< double > &)
 
typedef void(* UnsteadyExactSolutionFctPt) (const double &, const Vector< double > &, Vector< double > &)
 
- Static Public Attributes inherited from oomph::PolarNavierStokesEquations
static Vector< double > Gamma
 Vector to decide whether the stress-divergence form is used or not. 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::PolarNavierStokesEquations
double * Viscosity_Ratio_pt
 
double * Density_Ratio_pt
 
double * Alpha_pt
 Pointer to the angle alpha. More...
 
double * Re_pt
 Pointer to global Reynolds number. More...
 
double * ReSt_pt
 Pointer to global Reynolds number x Strouhal number (=Womersley) More...
 
double * ReInvFr_pt
 
Vector< double > * G_pt
 Pointer to global gravity Vector. More...
 
NavierStokesBodyForceFctPt Body_force_fct_pt
 Pointer to body force function. More...
 
NavierStokesSourceFctPt Source_fct_pt
 Pointer to volumetric source function. More...
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElement * Macro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepper * Geom_object_time_stepper_pt
 
- Protected Attributes inherited from oomph::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...
 
- 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
 
- Static Protected Attributes inherited from oomph::RefineableElement
static double Max_integrity_tolerance = 1.0e-8
 Max. allowed discrepancy in element integrity check. More...
 

Detailed Description

Refineable version of my Polar Navier–Stokes equations

Constructor & Destructor Documentation

◆ RefineablePolarNavierStokesEquations()

oomph::RefineablePolarNavierStokesEquations::RefineablePolarNavierStokesEquations ( )
inline

Constructor.

77  : PolarNavierStokesEquations(),
78  RefineableElement(),
79  ElementWithZ2ErrorEstimator()
80  {
81  }
oomph::ElementWithZ2ErrorEstimator::ElementWithZ2ErrorEstimator
ElementWithZ2ErrorEstimator()
Default empty constructor.
Definition: error_estimator.h:82
oomph::PolarNavierStokesEquations::PolarNavierStokesEquations
PolarNavierStokesEquations()
Constructor: NULL the body force and source function.
Definition: polar_navier_stokes_elements.h:245
oomph::RefineableElement::RefineableElement
RefineableElement()
Definition: refineable_elements.h:188

Member Function Documentation

◆ fill_in_generic_residual_contribution()

void oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
DenseMatrix< double > &  mass_matrix,
unsigned  flag 
)
protectedvirtual

Start of what would've been refineable_navier_stokes_elements.cc //.

Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both flag=0: compute only residual vector

Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobian as well. flag=2 for Residuals, Jacobian and mass_matrix

This is now my new version with Jacobian and dimensionless phi

This version supports hanging nodes

Reimplemented from oomph::PolarNavierStokesEquations.

52  {
53  // Find out how many nodes there are
54  unsigned n_node = nnode();
55 
56  // Find out how many pressure dofs there are
57  unsigned n_pres = npres_pnst();
58 
59  // Find the indices at which the local velocities are stored
60  unsigned u_nodal_index[2];
61  for (unsigned i = 0; i < 2; i++)
62  {
63  u_nodal_index[i] = u_index_pnst(i);
64  }
65 
66  // Which nodal value represents the pressure? (Negative if pressure
67  // is not based on nodal interpolation).
68  int p_index = this->p_nodal_index_pnst();
69 
70  // Local array of booleans that are true if the l-th pressure value is
71  // hanging (avoid repeated virtual function calls)
72  bool pressure_dof_is_hanging[n_pres];
73  // If the pressure is stored at a node
74  if (p_index >= 0)
75  {
76  // Read out whether the pressure is hanging
77  for (unsigned l = 0; l < n_pres; ++l)
78  {
79  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
80  }
81  }
82  // Otherwise the pressure is not stored at a node and so cannot hang
83  else
84  {
85  for (unsigned l = 0; l < n_pres; ++l)
86  {
87  pressure_dof_is_hanging[l] = false;
88  }
89  }
90 
91  // Set up memory for the shape and test functions
92  Shape psif(n_node), testf(n_node);
93  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
94 
95  // Set up memory for pressure shape and test functions
96  Shape psip(n_pres), testp(n_pres);
97 
98  // Number of integration points
99  unsigned n_intpt = integral_pt()->nweight();
100 
101  // Set the Vector to hold local coordinates
102  Vector<double> s(2);
103 
104  // Get the reynolds number and Alpha
105  const double Re = re();
106  const double Alpha = alpha();
107  const double Re_St = re_st();
108 
109  // Integers to store the local equations and unknowns
110  int local_eqn = 0, local_unknown = 0;
111 
112  // Pointers to hang info objects
113  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
114 
115  // Loop over the integration points
116  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
117  {
118  // Assign values of s
119  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
120  // Get the integral weight
121  double w = integral_pt()->weight(ipt);
122 
123  // Call the derivatives of the shape and test functions
124  double J = dshape_and_dtest_eulerian_at_knot_pnst(
125  ipt, psif, dpsifdx, testf, dtestfdx);
126 
127  // Call the pressure shape and test functions
128  this->pshape_pnst(s, psip, testp);
129 
130  // Premultiply the weights and the Jacobian
131  double W = w * J;
132 
133  // Calculate local values of the pressure and velocity components
134  // Allocate storage initialised to zero
135  double interpolated_p = 0.0;
136  Vector<double> interpolated_u(2, 0.0);
137  Vector<double> interpolated_x(2, 0.0);
138  // Vector<double> dudt(2);
139  DenseMatrix<double> interpolated_dudx(2, 2, 0.0);
140 
141  // Initialise to zero
142  for (unsigned i = 0; i < 2; i++)
143  {
144  // dudt[i] = 0.0;
145  interpolated_u[i] = 0.0;
146  interpolated_x[i] = 0.0;
147  for (unsigned j = 0; j < 2; j++)
148  {
149  interpolated_dudx(i, j) = 0.0;
150  }
151  }
152 
153  // Calculate pressure
154  for (unsigned l = 0; l < n_pres; l++)
155  interpolated_p += this->p_pnst(l) * psip[l];
156 
157  // Calculate velocities and derivatives:
158 
159  // Loop over nodes
160  for (unsigned l = 0; l < n_node; l++)
161  {
162  // Loop over directions
163  for (unsigned i = 0; i < 2; i++)
164  {
165  // Get the nodal value
166  double u_value = this->nodal_value(l, u_nodal_index[i]);
167  interpolated_u[i] += u_value * psif[l];
168  interpolated_x[i] += this->nodal_position(l, i) * psif[l];
169  // dudt[i]+=du_dt_pnst(l,i)*psif[l];
170 
171  // Loop over derivative directions
172  for (unsigned j = 0; j < 2; j++)
173  {
174  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
175  }
176  }
177  }
178 
179  // MOMENTUM EQUATIONS
180  //------------------
181  // Number of master nodes and storage for the weight of the shape function
182  unsigned n_master = 1;
183  double hang_weight = 1.0;
184  // Loop over the nodes for the test functions
185  for (unsigned l = 0; l < n_node; l++)
186  {
187  // Local boolean to indicate whether the node is hanging
188  bool is_node_hanging = node_pt(l)->is_hanging();
189 
190  // If the node is hanging
191  if (is_node_hanging)
192  {
193  hang_info_pt = node_pt(l)->hanging_pt();
194  // Read out number of master nodes from hanging data
195  n_master = hang_info_pt->nmaster();
196  }
197  // Otherwise the node is its own master
198  else
199  {
200  n_master = 1;
201  }
202 
203  // Now add in a new loop over the master nodes
204  for (unsigned m = 0; m < n_master; m++)
205  {
206  // Can't loop over velocity components as don't have identical
207  // contributions Do seperately for i = {0,1} instead
208  unsigned i = 0;
209  {
210  // Get the equation number
211  // If the node is hanging
212  if (is_node_hanging)
213  {
214  // Get the equation number from the master node
215  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
216  u_nodal_index[i]);
217  // Get the hang weight from the master node
218  hang_weight = hang_info_pt->master_weight(m);
219  }
220  // If the node is not hanging
221  else
222  {
223  // Local equation number
224  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
225  // Node contributes with full weight
226  hang_weight = 1.0;
227  }
228 
229  /*IF it's not a boundary condition*/
230  if (local_eqn >= 0)
231  {
232  // Add the testf[l] term of the stress tensor
233  residuals[local_eqn] +=
234  ((interpolated_p / interpolated_x[0]) -
235  ((1. + Gamma[i]) / pow(interpolated_x[0], 2.)) *
236  ((1. / Alpha) * interpolated_dudx(1, 1) +
237  interpolated_u[0])) *
238  testf[l] * interpolated_x[0] * Alpha * W * hang_weight;
239 
240  // Add the dtestfdx(l,0) term of the stress tensor
241  residuals[local_eqn] +=
242  (interpolated_p - (1. + Gamma[i]) * interpolated_dudx(0, 0)) *
243  dtestfdx(l, 0) * interpolated_x[0] * Alpha * W * hang_weight;
244 
245  // Add the dtestfdx(l,1) term of the stress tensor
246  residuals[local_eqn] -=
247  ((1. / (interpolated_x[0] * Alpha)) * interpolated_dudx(0, 1) -
248  (interpolated_u[1] / interpolated_x[0]) +
249  Gamma[i] * interpolated_dudx(1, 0)) *
250  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
251  interpolated_x[0] * Alpha * W * hang_weight;
252 
253  // Convective terms
254  residuals[local_eqn] -=
255  Re *
256  (interpolated_u[0] * interpolated_dudx(0, 0) +
257  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
258  interpolated_dudx(0, 1) -
259  (pow(interpolated_u[1], 2.) / interpolated_x[0])) *
260  testf[l] * interpolated_x[0] * Alpha * W * hang_weight;
261 
262 
263  // CALCULATE THE JACOBIAN
264  if (flag)
265  {
266  // Number of master nodes and weights
267  unsigned n_master2 = 1;
268  double hang_weight2 = 1.0;
269  // Loop over the velocity shape functions again
270  for (unsigned l2 = 0; l2 < n_node; l2++)
271  {
272  // Local boolean to indicate whether the node is hanging
273  bool is_node2_hanging = node_pt(l2)->is_hanging();
274 
275  // If the node is hanging
276  if (is_node2_hanging)
277  {
278  hang_info2_pt = node_pt(l2)->hanging_pt();
279  // Read out number of master nodes from hanging data
280  n_master2 = hang_info2_pt->nmaster();
281  }
282  // Otherwise the node is its own master
283  else
284  {
285  n_master2 = 1;
286  }
287 
288  // Loop over the master nodes
289  for (unsigned m2 = 0; m2 < n_master2; m2++)
290  {
291  // Again can't loop over velocity components due to loss of
292  // symmetry
293  unsigned i2 = 0;
294  {
295  // Get the number of the unknown
296  // If the node is hanging
297  if (is_node2_hanging)
298  {
299  // Get the equation number from the master node
300  local_unknown = this->local_hang_eqn(
301  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
302  // Get the hang weights
303  hang_weight2 = hang_info2_pt->master_weight(m2);
304  }
305  else
306  {
307  local_unknown =
308  this->nodal_local_eqn(l2, u_nodal_index[i2]);
309  hang_weight2 = 1.0;
310  }
311 
312  // If at a non-zero degree of freedom add in the entry
313  if (local_unknown >= 0)
314  {
315  // Add contribution to Elemental Matrix
316  jacobian(local_eqn, local_unknown) -=
317  (1. + Gamma[i]) *
318  (psif[l2] / pow(interpolated_x[0], 2.)) * testf[l] *
319  interpolated_x[0] * Alpha * W * hang_weight *
320  hang_weight2;
321 
322  jacobian(local_eqn, local_unknown) -=
323  (1. + Gamma[i]) * dpsifdx(l2, 0) * dtestfdx(l, 0) *
324  interpolated_x[0] * Alpha * W * hang_weight *
325  hang_weight2;
326 
327  jacobian(local_eqn, local_unknown) -=
328  (1. / (interpolated_x[0] * Alpha)) * dpsifdx(l2, 1) *
329  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
330  interpolated_x[0] * Alpha * W * hang_weight *
331  hang_weight2;
332 
333  // Now add in the inertial terms
334  jacobian(local_eqn, local_unknown) -=
335  Re *
336  (psif[l2] * interpolated_dudx(0, 0) +
337  interpolated_u[0] * dpsifdx(l2, 0) +
338  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
339  dpsifdx(l2, 1)) *
340  testf[l] * interpolated_x[0] * Alpha * W *
341  hang_weight * hang_weight2;
342 
343  // extra bit for mass matrix
344  if (flag == 2)
345  {
346  mass_matrix(local_eqn, local_unknown) +=
347  Re_St * psif[l2] * testf[l] * interpolated_x[0] *
348  Alpha * W * hang_weight * hang_weight2;
349  }
350 
351  } // End of (Jacobian's) if not boundary condition
352  // statement
353  } // End of i2=0 section
354 
355  i2 = 1;
356  {
357  // Get the number of the unknown
358  // If the node is hanging
359  if (is_node2_hanging)
360  {
361  // Get the equation number from the master node
362  local_unknown = this->local_hang_eqn(
363  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
364  // Get the hang weights
365  hang_weight2 = hang_info2_pt->master_weight(m2);
366  }
367  else
368  {
369  local_unknown =
370  this->nodal_local_eqn(l2, u_nodal_index[i2]);
371  hang_weight2 = 1.0;
372  }
373 
374  // If at a non-zero degree of freedom add in the entry
375  if (local_unknown >= 0)
376  {
377  // Add contribution to Elemental Matrix
378  jacobian(local_eqn, local_unknown) -=
379  ((1. + Gamma[i]) /
380  (pow(interpolated_x[0], 2.) * Alpha)) *
381  dpsifdx(l2, 1) * testf[l] * interpolated_x[0] *
382  Alpha * W * hang_weight * hang_weight2;
383 
384  jacobian(local_eqn, local_unknown) -=
385  (-(psif[l2] / interpolated_x[0]) +
386  Gamma[i] * dpsifdx(l2, 0)) *
387  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
388  interpolated_x[0] * Alpha * W * hang_weight *
389  hang_weight2;
390 
391  // Now add in the inertial terms
392  jacobian(local_eqn, local_unknown) -=
393  Re *
394  ((psif[l2] / (interpolated_x[0] * Alpha)) *
395  interpolated_dudx(0, 1) -
396  2 * interpolated_u[1] *
397  (psif[l2] / interpolated_x[0])) *
398  testf[l] * interpolated_x[0] * Alpha * W *
399  hang_weight * hang_weight2;
400 
401  } // End of (Jacobian's) if not boundary condition
402  // statement
403  } // End of i2=1 section
404 
405  } // End of loop over master nodes m2
406  } // End of l2 loop
407 
408  /*Now loop over pressure shape functions*/
409  /*This is the contribution from pressure gradient*/
410  for (unsigned l2 = 0; l2 < n_pres; l2++)
411  {
412  // If the pressure dof is hanging
413  if (pressure_dof_is_hanging[l2])
414  {
415  hang_info2_pt =
416  this->pressure_node_pt(l2)->hanging_pt(p_index);
417  // Pressure dof is hanging so it must be nodal-based
418  // Get the number of master nodes from the pressure node
419  n_master2 = hang_info2_pt->nmaster();
420  }
421  // Otherwise the node is its own master
422  else
423  {
424  n_master2 = 1;
425  }
426 
427  // Loop over the master nodes
428  for (unsigned m2 = 0; m2 < n_master2; m2++)
429  {
430  // Get the number of the unknown
431  // If the pressure dof is hanging
432  if (pressure_dof_is_hanging[l2])
433  {
434  // Get the unknown from the master node
435  local_unknown = this->local_hang_eqn(
436  hang_info2_pt->master_node_pt(m2), p_index);
437  // Get the weight from the hanging object
438  hang_weight2 = hang_info2_pt->master_weight(m2);
439  }
440  else
441  {
442  local_unknown = this->p_local_eqn(l2);
443  hang_weight2 = 1.0;
444  }
445 
446  /*If we are at a non-zero degree of freedom in the entry*/
447  if (local_unknown >= 0)
448  {
449  jacobian(local_eqn, local_unknown) +=
450  (psip[l2] / interpolated_x[0]) * testf[l] *
451  interpolated_x[0] * Alpha * W * hang_weight *
452  hang_weight2;
453 
454  jacobian(local_eqn, local_unknown) +=
455  psip[l2] * dtestfdx(l, 0) * interpolated_x[0] * Alpha *
456  W * hang_weight * hang_weight2;
457 
458  } // End of Jacobian pressure if not a boundary condition
459  // statement
460 
461  } // End of loop over master nodes m2
462  } // End of loop over pressure shape functions l2
463 
464  } /*End of Jacobian calculation*/
465 
466  } // End of if not boundary condition statement
467  } // End of i=0 section
468 
469  i = 1;
470  {
471  // Get the equation number
472  // If the node is hanging
473  if (is_node_hanging)
474  {
475  // Get the equation number from the master node
476  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
477  u_nodal_index[i]);
478  // Get the hang weight from the master node
479  hang_weight = hang_info_pt->master_weight(m);
480  }
481  // If the node is not hanging
482  else
483  {
484  // Local equation number
485  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
486  // Node contributes with full weight
487  hang_weight = 1.0;
488  }
489 
490  /*IF it's not a boundary condition*/
491  if (local_eqn >= 0)
492  {
493  // Add the testf[l] term of the stress tensor
494  residuals[local_eqn] +=
495  ((1. / (pow(interpolated_x[0], 2.) * Alpha)) *
496  interpolated_dudx(0, 1) -
497  (interpolated_u[1] / pow(interpolated_x[0], 2.)) +
498  Gamma[i] * (1. / interpolated_x[0]) *
499  interpolated_dudx(1, 0)) *
500  testf[l] * interpolated_x[0] * Alpha * W * hang_weight;
501 
502  // Add the dtestfdx(l,0) term of the stress tensor
503  residuals[local_eqn] -=
504  (interpolated_dudx(1, 0) +
505  Gamma[i] * ((1. / (interpolated_x[0] * Alpha)) *
506  interpolated_dudx(0, 1) -
507  (interpolated_u[1] / interpolated_x[0]))) *
508  dtestfdx(l, 0) * interpolated_x[0] * Alpha * W * hang_weight;
509 
510  // Add the dtestfdx(l,1) term of the stress tensor
511  residuals[local_eqn] +=
512  (interpolated_p -
513  (1. + Gamma[i]) * ((1. / (interpolated_x[0] * Alpha)) *
514  interpolated_dudx(1, 1) +
515  (interpolated_u[0] / interpolated_x[0]))) *
516  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
517  interpolated_x[0] * Alpha * W * hang_weight;
518 
519  // Convective terms
520  residuals[local_eqn] -=
521  Re *
522  (interpolated_u[0] * interpolated_dudx(1, 0) +
523  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
524  interpolated_dudx(1, 1) +
525  ((interpolated_u[0] * interpolated_u[1]) /
526  interpolated_x[0])) *
527  testf[l] * interpolated_x[0] * Alpha * W * hang_weight;
528 
529  // CALCULATE THE JACOBIAN
530  if (flag)
531  {
532  // Number of master nodes and weights
533  unsigned n_master2 = 1;
534  double hang_weight2 = 1.0;
535 
536  // Loop over the velocity shape functions again
537  for (unsigned l2 = 0; l2 < n_node; l2++)
538  {
539  // Local boolean to indicate whether the node is hanging
540  bool is_node2_hanging = node_pt(l2)->is_hanging();
541 
542  // If the node is hanging
543  if (is_node2_hanging)
544  {
545  hang_info2_pt = node_pt(l2)->hanging_pt();
546  // Read out number of master nodes from hanging data
547  n_master2 = hang_info2_pt->nmaster();
548  }
549  // Otherwise the node is its own master
550  else
551  {
552  n_master2 = 1;
553  }
554 
555  // Loop over the master nodes
556  for (unsigned m2 = 0; m2 < n_master2; m2++)
557  {
558  // Again can't loop over velocity components due to loss of
559  // symmetry
560  unsigned i2 = 0;
561  {
562  // Get the number of the unknown
563  // If the node is hanging
564  if (is_node2_hanging)
565  {
566  // Get the equation number from the master node
567  local_unknown = this->local_hang_eqn(
568  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
569  // Get the hang weights
570  hang_weight2 = hang_info2_pt->master_weight(m2);
571  }
572  else
573  {
574  local_unknown =
575  this->nodal_local_eqn(l2, u_nodal_index[i2]);
576  hang_weight2 = 1.0;
577  }
578 
579  // If at a non-zero degree of freedom add in the entry
580  if (local_unknown >= 0)
581  {
582  // Add contribution to Elemental Matrix
583  jacobian(local_eqn, local_unknown) +=
584  (1. / (pow(interpolated_x[0], 2.) * Alpha)) *
585  dpsifdx(l2, 1) * testf[l] * interpolated_x[0] *
586  Alpha * W * hang_weight * hang_weight2;
587 
588  jacobian(local_eqn, local_unknown) -=
589  Gamma[i] * (1. / (interpolated_x[0] * Alpha)) *
590  dpsifdx(l2, 1) * dtestfdx(l, 0) * interpolated_x[0] *
591  Alpha * W * hang_weight * hang_weight2;
592 
593  jacobian(local_eqn, local_unknown) -=
594  (1 + Gamma[i]) * (psif[l2] / interpolated_x[0]) *
595  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
596  interpolated_x[0] * Alpha * W * hang_weight *
597  hang_weight2;
598 
599  // Now add in the inertial terms
600  jacobian(local_eqn, local_unknown) -=
601  Re *
602  (psif[l2] * interpolated_dudx(1, 0) +
603  (psif[l2] * interpolated_u[1] / interpolated_x[0])) *
604  testf[l] * interpolated_x[0] * Alpha * W *
605  hang_weight * hang_weight2;
606 
607  } // End of (Jacobian's) if not boundary condition
608  // statement
609  } // End of i2=0 section
610 
611  i2 = 1;
612  {
613  // Get the number of the unknown
614  // If the node is hanging
615  if (is_node2_hanging)
616  {
617  // Get the equation number from the master node
618  local_unknown = this->local_hang_eqn(
619  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
620  // Get the hang weights
621  hang_weight2 = hang_info2_pt->master_weight(m2);
622  }
623  else
624  {
625  local_unknown =
626  this->nodal_local_eqn(l2, u_nodal_index[i2]);
627  hang_weight2 = 1.0;
628  }
629 
630  // If at a non-zero degree of freedom add in the entry
631  if (local_unknown >= 0)
632  {
633  // Add contribution to Elemental Matrix
634  jacobian(local_eqn, local_unknown) +=
635  (-(psif[l2] / pow(interpolated_x[0], 2.)) +
636  Gamma[i] * (1. / interpolated_x[0]) *
637  dpsifdx(l2, 0)) *
638  testf[l] * interpolated_x[0] * Alpha * W *
639  hang_weight * hang_weight2;
640 
641  jacobian(local_eqn, local_unknown) -=
642  (dpsifdx(l2, 0) -
643  Gamma[i] * (psif[l2] / interpolated_x[0])) *
644  dtestfdx(l, 0) * interpolated_x[0] * Alpha * W *
645  hang_weight * hang_weight2;
646 
647  jacobian(local_eqn, local_unknown) -=
648  (1. + Gamma[i]) * (1. / (interpolated_x[0] * Alpha)) *
649  dpsifdx(l2, 1) * (1. / (interpolated_x[0] * Alpha)) *
650  dtestfdx(l, 1) * interpolated_x[0] * Alpha * W *
651  hang_weight * hang_weight2;
652 
653  // Now add in the inertial terms
654  jacobian(local_eqn, local_unknown) -=
655  Re *
656  (interpolated_u[0] * dpsifdx(l2, 0) +
657  (psif[l2] / (interpolated_x[0] * Alpha)) *
658  interpolated_dudx(1, 1) +
659  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
660  dpsifdx(l2, 1) +
661  (interpolated_u[0] * psif[l2] / interpolated_x[0])) *
662  testf[l] * interpolated_x[0] * Alpha * W *
663  hang_weight * hang_weight2;
664 
665  // extra bit for mass matrix
666  if (flag == 2)
667  {
668  mass_matrix(local_eqn, local_unknown) +=
669  Re_St * psif[l2] * testf[l] * interpolated_x[0] *
670  Alpha * W * hang_weight * hang_weight2;
671  }
672 
673  } // End of (Jacobian's) if not boundary condition
674  // statement
675  } // End of i2=1 section
676 
677  } // End of loop over master nodes m2
678  } // End of l2 loop
679 
680  /*Now loop over pressure shape functions*/
681  /*This is the contribution from pressure gradient*/
682  for (unsigned l2 = 0; l2 < n_pres; l2++)
683  {
684  // If the pressure dof is hanging
685  if (pressure_dof_is_hanging[l2])
686  {
687  hang_info2_pt =
688  this->pressure_node_pt(l2)->hanging_pt(p_index);
689  // Pressure dof is hanging so it must be nodal-based
690  // Get the number of master nodes from the pressure node
691  n_master2 = hang_info2_pt->nmaster();
692  }
693  // Otherwise the node is its own master
694  else
695  {
696  n_master2 = 1;
697  }
698 
699  // Loop over the master nodes
700  for (unsigned m2 = 0; m2 < n_master2; m2++)
701  {
702  // Get the number of the unknown
703  // If the pressure dof is hanging
704  if (pressure_dof_is_hanging[l2])
705  {
706  // Get the unknown from the master node
707  local_unknown = this->local_hang_eqn(
708  hang_info2_pt->master_node_pt(m2), p_index);
709  // Get the weight from the hanging object
710  hang_weight2 = hang_info2_pt->master_weight(m2);
711  }
712  else
713  {
714  local_unknown = this->p_local_eqn(l2);
715  hang_weight2 = 1.0;
716  }
717 
718 
719  /*If we are at a non-zero degree of freedom in the entry*/
720  if (local_unknown >= 0)
721  {
722  jacobian(local_eqn, local_unknown) +=
723  (psip[l2] / interpolated_x[0]) * (1. / Alpha) *
724  dtestfdx(l, 1) * interpolated_x[0] * Alpha * W *
725  hang_weight * hang_weight2;
726 
727  } // End of if not boundary condition for pressure in
728  // jacobian
729 
730  } // End of loop over master nodes m2
731  } // End of loop over pressure test functions l2
732 
733  } /*End of Jacobian calculation*/
734 
735  } // End of if not boundary condition statement
736  } // End of i=1 section
737 
738  } // End of loop over master nodes m
739  } // End of loop over shape functions
740 
741 
742  // CONTINUITY EQUATION
743  //-------------------
744 
745  // Loop over the shape functions
746  for (unsigned l = 0; l < n_pres; l++)
747  {
748  // If the pressure dof is hanging
749  if (pressure_dof_is_hanging[l])
750  {
751  // Pressure dof is hanging so it must be nodal-based
752  // Get the hang info object
753  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
754  // Get the number of master nodes from the pressure node
755  n_master = hang_info_pt->nmaster();
756  }
757  // Otherwise the node is its own master
758  else
759  {
760  n_master = 1;
761  }
762 
763  // Loop over the master nodes
764  for (unsigned m = 0; m < n_master; m++)
765  {
766  // Get the number of the unknown
767  // If the pressure dof is hanging
768  if (pressure_dof_is_hanging[l])
769  {
770  // Get the local equation from the master node
771  local_eqn =
772  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
773  // Get the weight for the node
774  hang_weight = hang_info_pt->master_weight(m);
775  }
776  else
777  {
778  local_eqn = this->p_local_eqn(l);
779  hang_weight = 1.0;
780  }
781 
782  // If not a boundary conditions
783  if (local_eqn >= 0)
784  {
785  residuals[local_eqn] +=
786  (interpolated_dudx(0, 0) +
787  (interpolated_u[0] / interpolated_x[0]) +
788  (1. / (interpolated_x[0] * Alpha)) * interpolated_dudx(1, 1)) *
789  testp[l] * interpolated_x[0] * Alpha * W * hang_weight;
790 
791  /*CALCULATE THE JACOBIAN*/
792  if (flag)
793  {
794  unsigned n_master2 = 1;
795  double hang_weight2 = 1.0;
796  /*Loop over the velocity shape functions*/
797  for (unsigned l2 = 0; l2 < n_node; l2++)
798  {
799  // Local boolean to indicate whether the node is hanging
800  bool is_node2_hanging = node_pt(l2)->is_hanging();
801 
802  // If the node is hanging
803  if (is_node2_hanging)
804  {
805  hang_info2_pt = node_pt(l2)->hanging_pt();
806  // Read out number of master nodes from hanging data
807  n_master2 = hang_info2_pt->nmaster();
808  }
809  // Otherwise the node is its own master
810  else
811  {
812  n_master2 = 1;
813  }
814 
815  // Loop over the master nodes
816  for (unsigned m2 = 0; m2 < n_master2; m2++)
817  {
818  unsigned i2 = 0;
819  {
820  // Get the number of the unknown
821  // If the node is hanging
822  if (is_node2_hanging)
823  {
824  // Get the equation number from the master node
825  local_unknown = this->local_hang_eqn(
826  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
827  hang_weight2 = hang_info2_pt->master_weight(m2);
828  }
829  else
830  {
831  local_unknown =
832  this->nodal_local_eqn(l2, u_nodal_index[i2]);
833  hang_weight2 = 1.0;
834  }
835  /*If we're at a non-zero degree of freedom add it in*/
836  if (local_unknown >= 0)
837  {
838  jacobian(local_eqn, local_unknown) +=
839  (dpsifdx(l2, 0) + (psif[l2] / interpolated_x[0])) *
840  testp[l] * interpolated_x[0] * Alpha * W * hang_weight *
841  hang_weight2;
842  }
843  } // End of i2=0 section
844 
845  i2 = 1;
846  {
847  // Get the number of the unknown
848  // If the node is hanging
849  if (is_node2_hanging)
850  {
851  // Get the equation number from the master node
852  local_unknown = this->local_hang_eqn(
853  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
854  hang_weight2 = hang_info2_pt->master_weight(m2);
855  }
856  else
857  {
858  local_unknown =
859  this->nodal_local_eqn(l2, u_nodal_index[i2]);
860  hang_weight2 = 1.0;
861  }
862 
863  /*If we're at a non-zero degree of freedom add it in*/
864  if (local_unknown >= 0)
865  {
866  jacobian(local_eqn, local_unknown) +=
867  (1. / (interpolated_x[0] * Alpha)) * dpsifdx(l2, 1) *
868  testp[l] * interpolated_x[0] * Alpha * W * hang_weight *
869  hang_weight2;
870  }
871  } // End of i2=1 section
872 
873  } // End of loop over master nodes m2
874  } /*End of loop over l2*/
875  } /*End of Jacobian calculation*/
876 
877  } // End of if not boundary condition
878  } // End of loop over master nodes m
879  } // End of loop over pressure test functions l
880 
881  } // End of loop over integration points
882 
883  } // End of add_generic_residual_contribution
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
J
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
m2
MatrixType m2(n_dims)
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::nodal_value
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
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::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
oomph::FiniteElement::nodal_position
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317
oomph::HangInfo::nmaster
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::Node::hanging_pt
HangInfo *const & hanging_pt() const
Definition: nodes.h:1228
oomph::Node::is_hanging
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
oomph::PolarNavierStokesEquations::pshape_pnst
virtual void pshape_pnst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
oomph::PolarNavierStokesEquations::Gamma
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition: polar_navier_stokes_elements.h:261
oomph::PolarNavierStokesEquations::p_local_eqn
virtual int p_local_eqn(const unsigned &n)=0
oomph::PolarNavierStokesEquations::alpha
const double & alpha() const
Alpha.
Definition: polar_navier_stokes_elements.h:272
oomph::PolarNavierStokesEquations::p_pnst
virtual double p_pnst(const unsigned &n_p) const =0
oomph::PolarNavierStokesEquations::p_nodal_index_pnst
virtual int p_nodal_index_pnst()
Definition: polar_navier_stokes_elements.h:432
oomph::PolarNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_pnst
virtual double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
oomph::PolarNavierStokesEquations::u_index_pnst
virtual unsigned u_index_pnst(const unsigned &i) const
Definition: polar_navier_stokes_elements.h:394
oomph::PolarNavierStokesEquations::re
const double & re() const
Reynolds number.
Definition: polar_navier_stokes_elements.h:266
oomph::PolarNavierStokesEquations::re_st
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: polar_navier_stokes_elements.h:278
oomph::PolarNavierStokesEquations::npres_pnst
virtual unsigned npres_pnst() const =0
Function to return number of pressure degrees of freedom.
oomph::RefineableElement::local_hang_eqn
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Definition: refineable_elements.h:278
oomph::RefineablePolarNavierStokesEquations::pressure_node_pt
virtual Node * pressure_node_pt(const unsigned &n_p)
Definition: refineable_polar_navier_stokes_elements.h:62
s
RealScalar s
Definition: level1_cplx_impl.h:130
m
int * m
Definition: level2_cplx_impl.h:294
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
GlobalPhysicalVariables::Re
double Re
Reynolds number.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:61
TanhSolnForAdvectionDiffusion::Alpha
double Alpha
Parameter for steepness of step.
Definition: two_d_adv_diff_adapt.cc:53
oomph::Problem_Parameter::Re_St
double Re_St
Reynolds multiplied Strouhal.
Definition: axisym_vibrating_shell_non_newtonian.cc:91
oomph::QuadTreeNames::W
@ W
Definition: quadtree.h:63
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References TanhSolnForAdvectionDiffusion::Alpha, oomph::PolarNavierStokesEquations::alpha(), oomph::PolarNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_pnst(), oomph::PolarNavierStokesEquations::Gamma, oomph::Node::hanging_pt(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), oomph::Node::is_hanging(), J, j, oomph::Integral::knot(), oomph::RefineableElement::local_hang_eqn(), m, m2(), oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::FiniteElement::node_pt(), oomph::PolarNavierStokesEquations::npres_pnst(), oomph::Integral::nweight(), oomph::PolarNavierStokesEquations::p_local_eqn(), oomph::PolarNavierStokesEquations::p_nodal_index_pnst(), oomph::PolarNavierStokesEquations::p_pnst(), Eigen::bfloat16_impl::pow(), pressure_node_pt(), oomph::PolarNavierStokesEquations::pshape_pnst(), GlobalPhysicalVariables::Re, oomph::PolarNavierStokesEquations::re(), oomph::Problem_Parameter::Re_St, oomph::PolarNavierStokesEquations::re_st(), s, oomph::PolarNavierStokesEquations::u_index_pnst(), w, oomph::QuadTreeNames::W, and oomph::Integral::weight().

◆ further_build()

void oomph::RefineablePolarNavierStokesEquations::further_build ( )
inlinevirtual

Further build, pass the pointers down to the sons.

Reimplemented from oomph::RefineableElement.

Reimplemented in oomph::RefineablePolarCrouzeixRaviartElement.

175  {
176  // Find the father element
177  RefineablePolarNavierStokesEquations* cast_father_element_pt =
178  dynamic_cast<RefineablePolarNavierStokesEquations*>(
179  this->father_element_pt());
180 
181  // Set the viscosity ratio pointer
182  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
183  // Set the density ratio pointer
184  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
185  // Set pointer to global Reynolds number
186  this->Re_pt = cast_father_element_pt->re_pt();
187  // Set pointer to global Reynolds number x Strouhal number (=Womersley)
188  this->ReSt_pt = cast_father_element_pt->re_st_pt();
189  // Set pointer to global Reynolds number x inverse Froude number
190  this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
191  // Set pointer to global gravity Vector
192  this->G_pt = cast_father_element_pt->g_pt();
193  // Set pointer to alpha
194  this->Alpha_pt = cast_father_element_pt->alpha_pt();
195 
196  // Set pointer to body force function
197  this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
198 
199  // Set pointer to volumetric source function
200  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
201  }
oomph::PolarNavierStokesEquations::ReSt_pt
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: polar_navier_stokes_elements.h:157
oomph::PolarNavierStokesEquations::Viscosity_Ratio_pt
double * Viscosity_Ratio_pt
Definition: polar_navier_stokes_elements.h:142
oomph::PolarNavierStokesEquations::G_pt
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: polar_navier_stokes_elements.h:164
oomph::PolarNavierStokesEquations::Body_force_fct_pt
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
Definition: polar_navier_stokes_elements.h:167
oomph::PolarNavierStokesEquations::Re_pt
double * Re_pt
Pointer to global Reynolds number.
Definition: polar_navier_stokes_elements.h:154
oomph::PolarNavierStokesEquations::Density_Ratio_pt
double * Density_Ratio_pt
Definition: polar_navier_stokes_elements.h:146
oomph::PolarNavierStokesEquations::Source_fct_pt
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
Definition: polar_navier_stokes_elements.h:170
oomph::PolarNavierStokesEquations::ReInvFr_pt
double * ReInvFr_pt
Definition: polar_navier_stokes_elements.h:161
oomph::PolarNavierStokesEquations::Alpha_pt
double * Alpha_pt
Pointer to the angle alpha.
Definition: polar_navier_stokes_elements.h:151
oomph::RefineableElement::father_element_pt
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Definition: refineable_elements.h:539
oomph::RefineablePolarNavierStokesEquations::RefineablePolarNavierStokesEquations
RefineablePolarNavierStokesEquations()
Constructor.
Definition: refineable_polar_navier_stokes_elements.h:76

References oomph::PolarNavierStokesEquations::Alpha_pt, oomph::PolarNavierStokesEquations::alpha_pt(), oomph::PolarNavierStokesEquations::Body_force_fct_pt, oomph::PolarNavierStokesEquations::body_force_fct_pt(), oomph::PolarNavierStokesEquations::Density_Ratio_pt, oomph::PolarNavierStokesEquations::density_ratio_pt(), oomph::RefineableElement::father_element_pt(), oomph::PolarNavierStokesEquations::G_pt, oomph::PolarNavierStokesEquations::g_pt(), oomph::PolarNavierStokesEquations::re_invfr_pt(), oomph::PolarNavierStokesEquations::Re_pt, oomph::PolarNavierStokesEquations::re_pt(), oomph::PolarNavierStokesEquations::re_st_pt(), oomph::PolarNavierStokesEquations::ReInvFr_pt, oomph::PolarNavierStokesEquations::ReSt_pt, oomph::PolarNavierStokesEquations::Source_fct_pt, oomph::PolarNavierStokesEquations::source_fct_pt(), oomph::PolarNavierStokesEquations::Viscosity_Ratio_pt, and oomph::PolarNavierStokesEquations::viscosity_ratio_pt().

Referenced by oomph::RefineablePolarCrouzeixRaviartElement::further_build().

◆ get_Z2_flux()

void oomph::RefineablePolarNavierStokesEquations::get_Z2_flux ( const Vector< double > &  s,
Vector< double > &  flux 
)
inlinevirtual

Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.

Implements oomph::ElementWithZ2ErrorEstimator.

132  {
133 #ifdef PARANOID
134  unsigned num_entries = 3;
135  if (flux.size() != num_entries)
136  {
137  std::ostringstream error_message;
138  error_message << "The flux vector has the wrong number of entries, "
139  << flux.size() << ", whereas it should be " << num_entries
140  << std::endl;
141  throw OomphLibError(error_message.str(),
142  OOMPH_CURRENT_FUNCTION,
143  OOMPH_EXCEPTION_LOCATION);
144  }
145 #endif
146 
147  // Get strain rate matrix
148  DenseMatrix<double> strainrate(2);
149  // this->strain_rate(s,strainrate);
150  this->strain_rate_by_r(s, strainrate);
151 
152  // Pack into flux Vector
153  unsigned icount = 0;
154 
155  // Start with diagonal terms
156  for (unsigned i = 0; i < 2; i++)
157  {
158  flux[icount] = strainrate(i, i);
159  icount++;
160  }
161 
162  // Off diagonals row by row
163  for (unsigned i = 0; i < 2; i++)
164  {
165  for (unsigned j = i + 1; j < 2; j++)
166  {
167  flux[icount] = strainrate(i, j);
168  icount++;
169  }
170  }
171  }
oomph::PolarNavierStokesEquations::strain_rate_by_r
void strain_rate_by_r(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Function to return polar strain multiplied by r.
Definition: polar_navier_stokes_elements.cc:909
ProblemParameters::flux
void flux(const double &time, const Vector< double > &x, double &flux)
Get flux applied along boundary x=0.
Definition: pretend_melt.cc:59
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 ProblemParameters::flux(), i, j, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::PolarNavierStokesEquations::strain_rate_by_r().

◆ num_Z2_flux_terms()

unsigned oomph::RefineablePolarNavierStokesEquations::num_Z2_flux_terms ( )
inlinevirtual

Number of 'flux' terms for Z2 error estimation.

Implements oomph::ElementWithZ2ErrorEstimator.

124  {
125  // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
126  return 3;
127  }

◆ pin_elemental_redundant_nodal_pressure_dofs()

virtual void oomph::RefineablePolarNavierStokesEquations::pin_elemental_redundant_nodal_pressure_dofs ( )
inlineprotectedvirtual

Pin unused nodal pressure dofs (empty by default, because by default pressure dofs are not associated with nodes)

Reimplemented in oomph::RefineablePolarTaylorHoodElement.

72 {}

Referenced by pin_redundant_nodal_pressures().

◆ pin_redundant_nodal_pressures()

static void oomph::RefineablePolarNavierStokesEquations::pin_redundant_nodal_pressures ( const Vector< GeneralisedElement * > &  element_pt)
inlinestatic

Loop over all elements in Vector (which typically contains all the elements in a fluid mesh) and pin the nodal pressure degrees of freedom that are not being used. Function uses the member function

  • RefineablePolarNavierStokesEquations:: pin_elemental_redundant_nodal_pressure_dofs()

which is empty by default and should be implemented for elements with nodal pressure degrees of freedom (e.g. for refineable Taylor-Hood.)

96  {
97  // Loop over all elements and call the function that pins their
98  // unused nodal pressure data
99  unsigned n_element = element_pt.size();
100  for (unsigned e = 0; e < n_element; e++)
101  {
102  dynamic_cast<RefineablePolarNavierStokesEquations*>(element_pt[e])
103  ->pin_elemental_redundant_nodal_pressure_dofs();
104  }
105  }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::RefineablePolarNavierStokesEquations::pin_elemental_redundant_nodal_pressure_dofs
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Definition: refineable_polar_navier_stokes_elements.h:72

References e(), and pin_elemental_redundant_nodal_pressure_dofs().

◆ pressure_node_pt()

virtual Node* oomph::RefineablePolarNavierStokesEquations::pressure_node_pt ( const unsigned &  n_p)
inlineprotectedvirtual

Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interpolation).

Reimplemented from oomph::PolarNavierStokesEquations.

Reimplemented in oomph::RefineablePolarTaylorHoodElement.

63  {
64  return NULL;
65  }

Referenced by fill_in_generic_residual_contribution().

◆ unpin_all_pressure_dofs()

static void oomph::RefineablePolarNavierStokesEquations::unpin_all_pressure_dofs ( const Vector< GeneralisedElement * > &  element_pt)
inlinestatic

Unpin all pressure dofs in elements listed in vector.

110  {
111  // Loop over all elements to brutally unpin all nodal pressure degrees of
112  // freedom and internal pressure degrees of freedom
113  unsigned n_element = element_pt.size();
114  for (unsigned e = 0; e < n_element; e++)
115  {
116  dynamic_cast<RefineablePolarNavierStokesEquations*>(element_pt[e])
117  ->unpin_elemental_pressure_dofs();
118  }
119  }
oomph::RefineablePolarNavierStokesEquations::unpin_elemental_pressure_dofs
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.

References e(), and unpin_elemental_pressure_dofs().

◆ unpin_elemental_pressure_dofs()

virtual void oomph::RefineablePolarNavierStokesEquations::unpin_elemental_pressure_dofs ( )
protectedpure virtual

Unpin all pressure dofs in the element.

Implemented in oomph::RefineablePolarCrouzeixRaviartElement, and oomph::RefineablePolarTaylorHoodElement.

Referenced by unpin_all_pressure_dofs().

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