oomph::RefineableLinearisedAxisymmetricNavierStokesEquations Class Referenceabstract

#include <refineable_linearised_axisym_navier_stokes_elements.h>

+ Inheritance diagram for oomph::RefineableLinearisedAxisymmetricNavierStokesEquations:

Public Member Functions

 RefineableLinearisedAxisymmetricNavierStokesEquations ()
 Empty 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)
 
double geometric_jacobian (const Vector< double > &x)
 Fill in the geometric Jacobian, which in this case is r. More...
 
void further_build ()
 Further build: pass the pointers down to the sons. More...
 
 RefineableLinearisedAxisymmetricNavierStokesEquations ()
 Empty 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)
 
double geometric_jacobian (const Vector< double > &x)
 Fill in the geometric Jacobian, which in this case is r. More...
 
void further_build ()
 Further build: pass the pointers down to the sons. More...
 
- Public Member Functions inherited from oomph::LinearisedAxisymmetricNavierStokesEquations
 LinearisedAxisymmetricNavierStokesEquations ()
 
const doublere () const
 Reynolds number. More...
 
const doublere_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
const intazimuthal_mode_number () const
 Azimuthal mode number k in e^ik(theta) decomposition. More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
int *& azimuthal_mode_number_pt ()
 Pointer to azimuthal mode number k in e^ik(theta) decomposition. More...
 
const doubleviscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to the viscosity ratio. More...
 
const doubledensity_ratio () const
 
double *& density_ratio_pt ()
 Pointer to the density ratio. More...
 
virtual unsigned npres_linearised_axi_nst () const =0
 
virtual unsigned u_index_linearised_axi_nst (const unsigned &i) const
 
double du_dt_linearised_axi_nst (const unsigned &n, const unsigned &i) const
 
void disable_ALE ()
 
void enable_ALE ()
 
virtual double p_linearised_axi_nst (const unsigned &n_p, const unsigned &i) const =0
 
virtual int p_index_linearised_axi_nst (const unsigned &i) const
 Which nodal value represents the pressure? More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 
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 output_veloc (std::ostream &outfile, const unsigned &nplot, const unsigned &t)
 
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)
 
double interpolated_u_linearised_axi_nst (const Vector< double > &s, const unsigned &i) const
 
double interpolated_p_linearised_axi_nst (const Vector< double > &s, const unsigned &i) const
 
virtual void get_base_flow_u (const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
 
 LinearisedAxisymmetricNavierStokesEquations ()
 
const doublere () const
 Reynolds number. More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
const doublere_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
const doublere_invfr () const
 Global inverse Froude number. More...
 
double *& re_invfr_pt ()
 Pointer to global inverse Froude number. More...
 
const doublere_invro () const
 Global Reynolds number multiplied by inverse Rossby number. More...
 
double *& re_invro_pt ()
 Pointer to global Reynolds number multiplied by inverse Rossby number. More...
 
const Vector< double > & g () const
 Vector of gravitational components. More...
 
Vector< double > *& g_pt ()
 Pointer to Vector of gravitational components. More...
 
const intazimuthal_mode_number () const
 Azimuthal mode number k in e^ik(theta) decomposition. More...
 
int *& azimuthal_mode_number_pt ()
 Pointer to azimuthal mode number k in e^ik(theta) decomposition. More...
 
const doubleviscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to the viscosity ratio. More...
 
const doubledensity_ratio () const
 
double *& density_ratio_pt ()
 Pointer to the density ratio. More...
 
virtual int get_local_eqn_number_corresponding_to_geometric_dofs (const unsigned &n, const unsigned &i)=0
 
virtual unsigned npres_lin_axi_nst () const =0
 
virtual unsigned xhat_index_lin_axi_nst (const unsigned &i) const
 
virtual unsigned u_index_lin_axi_nst (const unsigned &i) const
 
double du_dt_lin_axi_nst (const unsigned &n, const unsigned &i) const
 
double dnodal_position_perturbation_dt_lin_axi_nst (const unsigned &n, const unsigned &i) const
 
unsigned nexternal_H_data ()
 Return the number of external perturbed spine "heights" data. More...
 
void disable_ALE ()
 
void enable_ALE ()
 
virtual double p_lin_axi_nst (const unsigned &n_p, const unsigned &i) const =0
 
virtual int p_index_lin_axi_nst (const unsigned &i) const
 Which nodal value represents the pressure? More...
 
void dkin_energy_dt (double &dkin_en_dt, double &kin_en) const
 Get integral of kinetic energy over element plus deriv w.r.t. time. More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 
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 output_veloc (std::ostream &outfile, const unsigned &nplot, const unsigned &t)
 
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)
 
double interpolated_u_lin_axi_nst (const Vector< double > &s, const unsigned &i) const
 
double interpolated_p_lin_axi_nst (const Vector< double > &s, const unsigned &i) const
 
double interpolated_nodal_position_perturbation_lin_axi_nst (const Vector< double > &s, const unsigned &i) const
 
void assign_additional_local_eqn_numbers ()
 Define the local equation numbering schemes. More...
 
 LinearisedAxisymmetricNavierStokesEquations ()
 
const doublere () const
 Reynolds number. More...
 
const doublere_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
const intazimuthal_mode_number () const
 Azimuthal mode number k in e^ik(theta) decomposition. More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
int *& azimuthal_mode_number_pt ()
 Pointer to azimuthal mode number k in e^ik(theta) decomposition. More...
 
const doubleviscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to the viscosity ratio. More...
 
const doubledensity_ratio () const
 
double *& density_ratio_pt ()
 Pointer to the density ratio. More...
 
virtual unsigned npres_linearised_axi_nst () const =0
 
virtual unsigned u_index_linearised_axi_nst (const unsigned &i) const
 
double du_dt_linearised_axi_nst (const unsigned &n, const unsigned &i) const
 
void disable_ALE ()
 
void enable_ALE ()
 
virtual double p_linearised_axi_nst (const unsigned &n_p, const unsigned &i) const =0
 
virtual int p_index_linearised_axi_nst (const unsigned &i) const
 Which nodal value represents the pressure? More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 
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 output_veloc (std::ostream &outfile, const unsigned &nplot, const unsigned &t)
 
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)
 
double interpolated_u_linearised_axi_nst (const Vector< double > &s, const unsigned &i) const
 
double interpolated_p_linearised_axi_nst (const Vector< double > &s, const unsigned &i) 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)
 
MacroElementmacro_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 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 Nodeconstruct_node (const unsigned &n)
 
virtual Nodeconstruct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual Nodeconstruct_boundary_node (const unsigned &n)
 
virtual Nodeconstruct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
int get_node_number (Node *const &node_pt) const
 
virtual Nodeget_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
 
Datageom_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, 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...
 
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 ()
 
TimeSteppertime_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...
 
Treetree_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 Nodeinterpolating_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 Nodeget_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 RefineableElementroot_element_pt ()
 
virtual RefineableElementfather_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 *, unsignedshape_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 Nodevertex_node_pt (const unsigned &j) const =0
 
virtual unsigned nrecovery_order ()=0
 Order of recovery shape functions. More...
 

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 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 doublemax_integrity_tolerance ()
 Max. allowed discrepancy in element integrity check. More...
 

Protected Member Functions

virtual Nodepressure_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 Nodepressure_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 ()
 
- Protected Member Functions inherited from oomph::LinearisedAxisymmetricNavierStokesEquations
virtual int p_local_eqn (const unsigned &n, const unsigned &i)=0
 
virtual double dshape_and_dtest_eulerian_linearised_axi_nst (const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst (const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual void pshape_linearised_axi_nst (const Vector< double > &s, Shape &psi) const =0
 Compute the pressure shape functions at local coordinate s. More...
 
virtual void pshape_linearised_axi_nst (const Vector< double > &s, Shape &psi, Shape &test) const =0
 Compute the pressure shape and test functions at local coordinate s. More...
 
virtual void get_base_flow_u (const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
 
virtual void get_base_flow_dudx (const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
 
virtual int p_local_eqn (const unsigned &n, const unsigned &i)=0
 
virtual double dshape_and_dtest_eulerian_lin_axi_nst (const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_lin_axi_nst (const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_and_dnodal_coordinates_at_knot_lin_axi_nst (const unsigned &ipt, Shape &psi, DShape &dpsidx, RankFourTensor< double > &d_dpsidx_dX, Shape &test, DShape &dtestdx, RankFourTensor< double > &d_dtestdx_dX, DenseMatrix< double > &djacobian_dX) const =0
 
virtual void pshape_lin_axi_nst (const Vector< double > &s, Shape &psi) const =0
 Compute the pressure shape functions at local coordinate s. More...
 
virtual void pshape_lin_axi_nst (const Vector< double > &s, Shape &psi, Shape &test) const =0
 Compute the pressure shape and test functions at local coordinate s. More...
 
virtual void get_base_flow_dudx (const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
 
virtual void get_base_flow_dudt (const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
 
virtual void get_base_flow_p (const double &time, const unsigned &ipt, const Vector< double > &x, double &result) const
 
virtual void get_base_flow_d_dudx_dX (const double &time, const unsigned &ipt, const Vector< double > &x, RankFourTensor< double > &result) const
 
void get_body_force_base_flow (const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result)
 
double get_source_base_flow (const double &time, const unsigned &ipt, const Vector< double > &x)
 
void get_body_force_gradient_base_flow (const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result)
 
void get_source_gradient_base_flow (const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result)
 
virtual void get_base_flow_duds (const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
 
virtual void fill_in_generic_residual_contribution_lin_axi_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 
virtual int p_local_eqn (const unsigned &n, const unsigned &i)=0
 
virtual double dshape_and_dtest_eulerian_linearised_axi_nst (const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst (const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual void pshape_linearised_axi_nst (const Vector< double > &s, Shape &psi) const =0
 Compute the pressure shape functions at local coordinate s. More...
 
virtual void pshape_linearised_axi_nst (const Vector< double > &s, Shape &psi, Shape &test) const =0
 Compute the pressure shape and test functions at local coordinate s. More...
 
virtual void get_base_flow_u (const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
 
virtual void get_base_flow_dudx (const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
 
- 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)
 
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)
 

Private Member Functions

void fill_in_generic_residual_contribution_linearised_axi_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 
void fill_in_generic_residual_contribution_linearised_axi_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 

Additional Inherited Members

- 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 Attributes inherited from oomph::LinearisedAxisymmetricNavierStokesEquations
void(*&)(const double &time, const Vector< double > &x, Vector< double > &fbase_flow_u_fct_pt ()
 Access function for the base flow solution pointer. More...
 
void(*&)(const double &time, const Vector< double > &x, DenseMatrix< double > &fbase_flow_dudx_fct_pt ()
 
void(*&)(const double &time, const Vector< double > &x, Vector< double > &fbase_flow_dudt_fct_pt ()
 
void(*&)(const double &time, const Vector< double > &x, double &fbase_flow_p_fct_pt ()
 Access function for the base flow pressure pointer. More...
 
void(*&)(const double &time, const Vector< double > &x, RankFourTensor< double > &fbase_flow_d_dudx_dX_fct_pt ()
 
void(*&)(const double &time, const Vector< double > &x, Vector< double > &fbody_force_fct_pt ()
 Access function for the body-force pointer. More...
 
double(*&)(const double &time, const Vector< double > &xsource_fct_pt ()
 Access function for the source-function pointer. More...
 
void(*&)(const double &time, const Vector< double > &x, DenseMatrix< double > &fbase_flow_duds_fct_pt ()
 
- Static Public Attributes inherited from oomph::LinearisedAxisymmetricNavierStokesEquations
static Vector< doubleGamma
 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::LinearisedAxisymmetricNavierStokesEquations
doubleViscosity_Ratio_pt
 
doubleDensity_Ratio_pt
 
doubleRe_pt
 Pointer to global Reynolds number. More...
 
doubleReSt_pt
 Pointer to global Reynolds number x Strouhal number (=Womersley) More...
 
intAzimuthal_Mode_Number_pt
 Pointer to azimuthal mode number k in e^ik(theta) decomposition. More...
 
void(* Base_flow_u_fct_pt )(const double &time, const Vector< double > &x, Vector< double > &result)
 Pointer to base flow solution (velocity components) function. More...
 
void(* Base_flow_dudx_fct_pt )(const double &time, const Vector< double > &x, DenseMatrix< double > &result)
 
bool ALE_is_disabled
 
doubleReInvFr_pt
 
doubleReInvRo_pt
 
Vector< double > * G_pt
 Pointer to global gravity Vector. More...
 
void(* Base_flow_dudt_fct_pt )(const double &time, const Vector< double > &x, Vector< double > &result)
 
void(* Base_flow_p_fct_pt )(const double &time, const Vector< double > &x, double &result)
 Pointer to base flow solution (pressure) function. More...
 
void(* Base_flow_d_dudx_dX_fct_pt )(const double &time, const Vector< double > &x, RankFourTensor< double > &result)
 
void(* Body_force_fct_pt )(const double &time, const Vector< double > &x, Vector< double > &result)
 Pointer to (base flow) body force function. More...
 
double(* Source_fct_pt )(const double &time, const Vector< double > &x)
 Pointer to (base flow) volumetric source function. More...
 
void(* Base_flow_duds_fct_pt )(const double &time, const Vector< double > &x, DenseMatrix< double > &result)
 
Vector< intExternal_H_local_eqn
 
Vector< unsignedExternal_node
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElementMacro_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...
 
TimeStepperGeom_object_time_stepper_pt
 
- Protected Attributes inherited from oomph::RefineableElement
TreeTree_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< doubleDummy_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 the linearised axisymmetric Navier–Stokes equations

Constructor & Destructor Documentation

◆ RefineableLinearisedAxisymmetricNavierStokesEquations() [1/2]

oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::RefineableLinearisedAxisymmetricNavierStokesEquations ( )
inline

Empty Constructor.

70  :
ElementWithZ2ErrorEstimator()
Default empty constructor.
Definition: error_estimator.h:82
LinearisedAxisymmetricNavierStokesEquations()
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:198
RefineableElement()
Definition: refineable_elements.h:188

◆ RefineableLinearisedAxisymmetricNavierStokesEquations() [2/2]

oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::RefineableLinearisedAxisymmetricNavierStokesEquations ( )
inline

Empty Constructor.

Member Function Documentation

◆ fill_in_generic_residual_contribution_linearised_axi_nst() [1/2]

void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_linearised_axi_nst ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
DenseMatrix< double > &  mass_matrix,
unsigned  flag 
)
privatevirtual

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

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

Compute the residuals for the refineable linearised axisymmetric Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobian as well.

Reimplemented from oomph::LinearisedAxisymmetricNavierStokesEquations.

45  {
46  // The dimension is actually two
47 // const unsigned DIM=2;
48 
49  // Determine number of nodes in the element
50  const unsigned n_node = nnode();
51 
52  // Get continuous time from timestepper of first node
53  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
54 
55  // Determine how many pressure values there are associated with
56  // a single pressure component
57  const unsigned n_pres = npres_linearised_axi_nst();
58 
59  // Get the nodal indices at which the velocity is stored
60  unsigned u_nodal_index[6];
61  for(unsigned i=0;i<6;++i)
62  {
63  u_nodal_index[i] = u_index_linearised_axi_nst(i);
64  }
65 
66  // Which nodal values represent the two pressure components?
67  // (Negative if pressure is not based on nodal interpolation).
68  Vector<int> p_index(2);
69  for(unsigned i=0;i<2;i++)
70  {
71  p_index[i] = this->p_index_linearised_axi_nst(i);
72  }
73 
74  // Local array of booleans that are true if the l-th pressure value is
75  // hanging (avoid repeated virtual function calls)
76  bool pressure_dof_is_hanging[n_pres];
77 
78  // Loop over two pressure components
79 // for(unsigned i=0;i<2;i++)
80 // {
81  // If the pressure is stored at a node
82  if(p_index[0] >= 0)
83  {
84  // Read out whether the pressure is hanging
85  for(unsigned l=0;l<n_pres;++l)
86  {
87  pressure_dof_is_hanging[l] =
88  pressure_node_pt(l)->is_hanging(p_index[0]);
89  }
90  }
91  // Otherwise the pressure is not stored at a node and so cannot hang
92  else
93  {
94  for(unsigned l=0;l<n_pres;++l)
95  { pressure_dof_is_hanging[l] = false; }
96  }
97 // } // End of loop over pressure components
98 
99  // Set up memory for the shape and test functions
100  // Note that there are two dimensions, r and z, in this problem
101  Shape psif(n_node), testf(n_node);
102  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
103 
104  // Set up memory for the pressure shape and test functions
105  Shape psip(n_pres), testp(n_pres);
106 
107  // Determine number of integration points
108  const unsigned n_intpt = integral_pt()->nweight();
109 
110  // Set up memory for the vector to hold local coordinates (two dimensions)
111  Vector<double> s(2);
112 
113  // Get physical variables from the element
114  // (Reynolds number must be multiplied by the density ratio)
115  const double scaled_re = re()*density_ratio();
116  const double scaled_re_st = re_st()*density_ratio();
117  const double visc_ratio = viscosity_ratio();
118  const int k = azimuthal_mode_number();
119 
120  // Integers used to store the local equation and unknown numbers
121  int local_eqn = 0, local_unknown = 0;
122 
123  // Local storage for pointers to hang info objects
124  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
125 
126  // Loop over the integration points
127  for(unsigned ipt=0;ipt<n_intpt;ipt++)
128  {
129  // Assign values of the local coordinates s
130  for(unsigned i=0;i<2;i++) { s[i] = integral_pt()->knot(ipt,i); }
131 
132  // Get the integral weight
133  const double w = integral_pt()->weight(ipt);
134 
135  // Calculate the fluid shape and test functions, and their derivatives
136  // w.r.t. the global coordinates
137  const double J =
139  testf,dtestfdx);
140 
141  // Calculate the pressure shape and test functions
142  pshape_linearised_axi_nst(s,psip,testp);
143 
144  // Premultiply the weights and the Jacobian of the mapping between
145  // local and global coordinates
146  const double W = w*J;
147 
148  // Allocate storage for the position and the derivative of the
149  // mesh positions w.r.t. time
150  Vector<double> interpolated_x(2,0.0);
151  Vector<double> mesh_velocity(2,0.0);
152 
153  // Allocate storage for the velocity components (six of these)
154  // and their derivatives w.r.t. time
155  Vector<double> interpolated_u(6,0.0);
156  Vector<double> dudt(6,0.0);
157 
158  // Allocate storage for the pressure components (two of these)
159  Vector<double> interpolated_p(2,0.0);
160 
161  // Allocate storage for the derivatives of the velocity components
162  // w.r.t. global coordinates (r and z)
163  DenseMatrix<double> interpolated_dudx(6,2,0.0);
164 
165  // Calculate pressure at the integration point
166  // -------------------------------------------
167 
168  // Loop over pressure degrees of freedom (associated with a single
169  // pressure component) in the element
170  for(unsigned l=0;l<n_pres;l++)
171  {
172  // Cache the shape function
173  const double psip_ = psip(l);
174 
175  // Loop over the two pressure components
176  for(unsigned i=0;i<2;i++)
177  {
178  // Get the value
179  const double p_value = this->p_linearised_axi_nst(l,i);
180 
181  // Add contribution
182  interpolated_p[i] += p_value*psip_;
183  }
184  } // End of loop over the pressure degrees of freedom in the element
185 
186  // Calculate velocities and their derivatives at the integration point
187  // -------------------------------------------------------------------
188 
189  // Loop over the element's nodes
190  for(unsigned l=0;l<n_node;l++)
191  {
192  // Cache the shape function
193  const double psif_ = psif(l);
194 
195  // Loop over the two coordinate directions
196  for(unsigned i=0;i<2;i++)
197  {
198  interpolated_x[i] += this->nodal_position(l,i)*psif_;
199  }
200 
201  // Loop over the six velocity components
202  for(unsigned i=0;i<6;i++)
203  {
204  // Get the value
205  const double u_value = this->nodal_value(l,u_nodal_index[i]);
206 
207  // Add contribution
208  interpolated_u[i] += u_value*psif_;
209 
210  // Add contribution to dudt
211  dudt[i] += du_dt_linearised_axi_nst(l,i)*psif_;
212 
213  // Loop over two coordinate directions (for derivatives)
214  for(unsigned j=0;j<2;j++)
215  {
216  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
217  }
218  }
219  } // End of loop over the element's nodes
220 
221  // Get the mesh velocity if ALE is enabled
222  if(!ALE_is_disabled)
223  {
224  // Loop over the element's nodes
225  for(unsigned l=0;l<n_node;l++)
226  {
227  // Loop over the two coordinate directions
228  for(unsigned i=0;i<2;i++)
229  {
230  mesh_velocity[i] += this->dnodal_position_dt(l,i)*psif(l);
231  }
232  }
233  }
234 
235  // Get velocities and their derivatives from base flow problem
236  // -----------------------------------------------------------
237 
238  // Allocate storage for the velocity components of the base state
239  // solution (initialise to zero)
240  Vector<double> base_flow_u(3,0.0);
241 
242  // Get the user-defined base state solution velocity components
243  get_base_flow_u(time,ipt,interpolated_x,base_flow_u);
244 
245  // Allocate storage for the derivatives of the base state solution's
246  // velocity components w.r.t. global coordinate (r and z)
247  // N.B. the derivatives of the base flow components w.r.t. the
248  // azimuthal coordinate direction (theta) are always zero since the
249  // base flow is axisymmetric
250  DenseMatrix<double> base_flow_dudx(3,2,0.0);
251 
252  // Get the derivatives of the user-defined base state solution
253  // velocity components w.r.t. global coordinates
254  get_base_flow_dudx(time,ipt,interpolated_x,base_flow_dudx);
255 
256  // Cache base flow velocities and their derivatives
257  const double interpolated_ur = base_flow_u[0];
258  const double interpolated_uz = base_flow_u[1];
259  const double interpolated_utheta = base_flow_u[2];
260  const double interpolated_durdr = base_flow_dudx(0,0);
261  const double interpolated_durdz = base_flow_dudx(0,1);
262  const double interpolated_duzdr = base_flow_dudx(1,0);
263  const double interpolated_duzdz = base_flow_dudx(1,1);
264  const double interpolated_duthetadr = base_flow_dudx(2,0);
265  const double interpolated_duthetadz = base_flow_dudx(2,1);
266 
267  // Cache r-component of position
268  const double r = interpolated_x[0];
269 
270  // Cache unknowns
271  const double interpolated_UC = interpolated_u[0];
272  const double interpolated_US = interpolated_u[1];
273  const double interpolated_WC = interpolated_u[2];
274  const double interpolated_WS = interpolated_u[3];
275  const double interpolated_VC = interpolated_u[4];
276  const double interpolated_VS = interpolated_u[5];
277  const double interpolated_PC = interpolated_p[0];
278  const double interpolated_PS = interpolated_p[1];
279 
280  // Cache derivatives of the unknowns
281  const double interpolated_dUCdr = interpolated_dudx(0,0);
282  const double interpolated_dUCdz = interpolated_dudx(0,1);
283  const double interpolated_dUSdr = interpolated_dudx(1,0);
284  const double interpolated_dUSdz = interpolated_dudx(1,1);
285  const double interpolated_dWCdr = interpolated_dudx(2,0);
286  const double interpolated_dWCdz = interpolated_dudx(2,1);
287  const double interpolated_dWSdr = interpolated_dudx(3,0);
288  const double interpolated_dWSdz = interpolated_dudx(3,1);
289  const double interpolated_dVCdr = interpolated_dudx(4,0);
290  const double interpolated_dVCdz = interpolated_dudx(4,1);
291  const double interpolated_dVSdr = interpolated_dudx(5,0);
292  const double interpolated_dVSdz = interpolated_dudx(5,1);
293 
294  // ==================
295  // MOMENTUM EQUATIONS
296  // ==================
297 
298  // Number of master nodes
299  unsigned n_master = 1;
300 
301  // Storage for the weight of the shape function
302  double hang_weight = 1.0;
303 
304  // Loop over the nodes for the test functions/equations
305  for(unsigned l=0;l<n_node;l++)
306  {
307  // Cache test functions and their derivatives
308  const double testf_ = testf(l);
309  const double dtestfdr = dtestfdx(l,0);
310  const double dtestfdz = dtestfdx(l,1);
311 
312  // Local boolean that indicates the hanging status of the node
313  bool is_node_hanging = node_pt(l)->is_hanging();
314 
315  // If the node is hanging
316  if(is_node_hanging)
317  {
318  // Get the hanging pointer
319  hang_info_pt = node_pt(l)->hanging_pt();
320 
321  // Read out number of master nodes from hanging data
322  n_master = hang_info_pt->nmaster();
323  }
324  // Otherwise the node is its own master
325  else
326  {
327  n_master = 1;
328  }
329 
330  // Loop over the master nodes
331  for(unsigned m=0;m<n_master;m++)
332  {
333  // Loop over velocity components for equations
334  for(unsigned i=0;i<6;i++)
335  {
336  // Get the equation number
337  // -----------------------
338 
339  // If the node is hanging
340  if(is_node_hanging)
341  {
342  // Get the equation number from the master node
343  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
344  u_nodal_index[i]);
345 
346  // Get the hang weight from the master node
347  hang_weight = hang_info_pt->master_weight(m);
348  }
349  // If the node is not hanging
350  else
351  {
352  // Local equation number
353  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
354 
355  // Node contributes with full weight
356  hang_weight = 1.0;
357  }
358 
359  // If it's not a boundary condition...
360  if(local_eqn>= 0)
361  {
362  // initialise for residual calculation
363  double sum = 0.0;
364 
365  switch(i)
366  {
367 
368  // ---------------------------------------------
369  // FIRST (RADIAL) MOMENTUM EQUATION: COSINE PART
370  // ---------------------------------------------
371 
372  case 0:
373 
374  // Pressure gradient term
375  sum += interpolated_PC*(testf_ + r*dtestfdr)*W*hang_weight;
376 
377  // Stress tensor terms
378  sum -= visc_ratio*
379  r*(1.0+Gamma[0])*interpolated_dUCdr*dtestfdr*W*hang_weight;
380 
381  sum -= visc_ratio*r*
382  (interpolated_dUCdz + Gamma[0]*interpolated_dWCdr)
383  *dtestfdz*W*hang_weight;
384 
385  sum += visc_ratio*(
386  (k*Gamma[0]*interpolated_dVSdr)
387  - (k*(2.0+Gamma[0])*interpolated_VS/r)
388  - ((1.0+Gamma[0]+(k*k))*interpolated_UC/r))
389  *testf_*W*hang_weight;
390 
391  // Inertial terms (du/dt)
392  sum -= scaled_re_st*r*dudt[0]*testf_*W*hang_weight;
393 
394  // Inertial terms (convective)
395  sum -=
396  scaled_re*(r*interpolated_ur*interpolated_dUCdr
397  + r*interpolated_durdr*interpolated_UC
398  + k*interpolated_utheta*interpolated_US
399  - 2*interpolated_utheta*interpolated_VC
400  + r*interpolated_uz*interpolated_dUCdz
401  + r*interpolated_durdz*interpolated_WC)
402  *testf_*W*hang_weight;
403 
404  // Mesh velocity terms
405  if(!ALE_is_disabled)
406  {
407  for(unsigned j=0;j<2;j++)
408  {
409  sum +=
410  scaled_re_st*r*mesh_velocity[j]
411  *interpolated_dudx(0,j)*testf_*W*hang_weight;
412  }
413  }
414 
415  break;
416 //asd
417  // --------------------------------------------
418  // SECOND (RADIAL) MOMENTUM EQUATION: SINE PART
419  // --------------------------------------------
420 
421  case 1:
422 
423  // Pressure gradient term
424  sum += interpolated_PS*(testf_ + r*dtestfdr)*W*hang_weight;
425 
426  // Stress tensor terms
427  sum -= visc_ratio*
428  r*(1.0+Gamma[0])*interpolated_dUSdr*dtestfdr*W*hang_weight;
429 
430  sum -= visc_ratio*r*
431  (interpolated_dUSdz + Gamma[0]*interpolated_dWSdr)*
432  dtestfdz*W*hang_weight;
433 
434  sum -= visc_ratio*(
435  (k*Gamma[0]*interpolated_dVCdr)
436  - (k*(2.0+Gamma[0])*interpolated_VC/r)
437  + ((1.0+Gamma[0]+(k*k))*interpolated_US/r))*testf_*W*hang_weight;
438 
439  // Inertial terms (du/dt)
440  sum -= scaled_re_st*r*dudt[1]*testf_*W*hang_weight;
441 
442  // Inertial terms (convective)
443  sum -=
444  scaled_re*(r*interpolated_ur*interpolated_dUSdr
445  + r*interpolated_durdr*interpolated_US
446  - k*interpolated_utheta*interpolated_UC
447  - 2*interpolated_utheta*interpolated_VS
448  + r*interpolated_uz*interpolated_dUSdz
449  + r*interpolated_durdz*interpolated_WS)*testf_*W*hang_weight;
450 
451  // Mesh velocity terms
452  if(!ALE_is_disabled)
453  {
454  for(unsigned j=0;j<2;j++)
455  {
456  sum +=
457  scaled_re_st*r*mesh_velocity[j]
458  *interpolated_dudx(1,j)*testf_*W*hang_weight;
459  }
460  }
461 
462  break;
463 
464  // --------------------------------------------
465  // THIRD (AXIAL) MOMENTUM EQUATION: COSINE PART
466  // --------------------------------------------
467 
468  case 2:
469 
470  // Pressure gradient term
471  sum += r*interpolated_PC*dtestfdz*W*hang_weight;
472 
473  // Stress tensor terms
474  sum -= visc_ratio*r*
475  (interpolated_dWCdr + Gamma[1]*interpolated_dUCdz)
476  *dtestfdr*W*hang_weight;
477 
478  sum -= visc_ratio*r*
479  (1.0 + Gamma[1])*interpolated_dWCdz*dtestfdz*W*hang_weight;
480 
481  sum += visc_ratio*k*
482  (Gamma[1]*interpolated_dVSdz - k*interpolated_WC/r)*testf_*W*hang_weight;
483 
484  // Inertial terms (du/dt)
485  sum -= scaled_re_st*r*dudt[2]*testf_*W*hang_weight;
486 
487  // Inertial terms (convective)
488  sum -=
489  scaled_re*(r*interpolated_ur*interpolated_dWCdr
490  + r*interpolated_duzdr*interpolated_UC
491  + k*interpolated_utheta*interpolated_WS
492  + r*interpolated_uz*interpolated_dWCdz
493  + r*interpolated_duzdz*interpolated_WC)*testf_*W*hang_weight;
494 
495  // Mesh velocity terms
496  if(!ALE_is_disabled)
497  {
498  for(unsigned j=0;j<2;j++)
499  {
500  sum +=
501  scaled_re_st*r*mesh_velocity[j]
502  *interpolated_dudx(2,j)*testf_*W*hang_weight;
503  }
504  }
505 
506  break;
507 
508  // -------------------------------------------
509  // FOURTH (AXIAL) MOMENTUM EQUATION: SINE PART
510  // -------------------------------------------
511 
512  case 3:
513 
514  // Pressure gradient term
515  sum += r*interpolated_PS*dtestfdz*W*hang_weight;
516 
517  // Stress tensor terms
518  sum -= visc_ratio*r*
519  (interpolated_dWSdr + Gamma[1]*interpolated_dUSdz)
520  *dtestfdr*W*hang_weight;
521 
522  sum -= visc_ratio*r*
523  (1.0+Gamma[1])*interpolated_dWSdz*dtestfdz*W*hang_weight;
524 
525  sum -= visc_ratio*k*
526  (Gamma[1]*interpolated_dVCdz + k*interpolated_WS/r)*testf_*W*hang_weight;
527 
528  // Inertial terms (du/dt)
529  sum -= scaled_re_st*r*dudt[3]*testf_*W*hang_weight;
530 
531  // Inertial terms (convective)
532  sum -=
533  scaled_re*(r*interpolated_ur*interpolated_dWSdr
534  + r*interpolated_duzdr*interpolated_US
535  - k*interpolated_utheta*interpolated_WC
536  + r*interpolated_uz*interpolated_dWSdz
537  + r*interpolated_duzdz*interpolated_WS)*testf_*W*hang_weight;
538 
539  // Mesh velocity terms
540  if(!ALE_is_disabled)
541  {
542  for(unsigned j=0;j<2;j++)
543  {
544  sum +=
545  scaled_re_st*r*mesh_velocity[j]
546  *interpolated_dudx(3,j)*testf_*W*hang_weight;
547  }
548  }
549 
550  break;
551 
552  // ------------------------------------------------
553  // FIFTH (AZIMUTHAL) MOMENTUM EQUATION: COSINE PART
554  // ------------------------------------------------
555 
556  case 4:
557 
558  // Pressure gradient term
559  sum -= k*interpolated_PS*testf_*W*hang_weight;
560 
561  // Stress tensor terms
562  sum += visc_ratio*
563  (-r*interpolated_dVCdr
564  - k*Gamma[0]*interpolated_US
565  + Gamma[0]*interpolated_VC)*dtestfdr*W*hang_weight;
566 
567  sum -= visc_ratio*
568  (k*Gamma[0]*interpolated_WS + r*interpolated_dVCdz)*dtestfdz*W*hang_weight;
569 
570  sum += visc_ratio*
571  (Gamma[0]*interpolated_dVCdr
572  + k*(2.0+Gamma[0])*interpolated_US/r
573  - (1.0+(k*k)+(k*k*Gamma[0]))*interpolated_VC/r)*testf_*W*hang_weight;
574 
575  // Inertial terms (du/dt)
576  sum -= scaled_re_st*r*dudt[4]*testf_*W*hang_weight;
577 
578  // Inertial terms (convective)
579  sum -=
580  scaled_re*(r*interpolated_ur*interpolated_dVCdr
581  + r*interpolated_duthetadr*interpolated_UC
582  + k*interpolated_utheta*interpolated_VS
583  + interpolated_utheta*interpolated_UC
584  + interpolated_ur*interpolated_VC
585  + r*interpolated_uz*interpolated_dVCdz
586  + r*interpolated_duthetadz*interpolated_WC)*testf_*W*hang_weight;
587 
588  // Mesh velocity terms
589  if(!ALE_is_disabled)
590  {
591  for(unsigned j=0;j<2;j++)
592  {
593  sum +=
594  scaled_re_st*r*mesh_velocity[j]
595  *interpolated_dudx(4,j)*testf_*W*hang_weight;
596  }
597  }
598 
599  break;
600 
601  // ----------------------------------------------
602  // SIXTH (AZIMUTHAL) MOMENTUM EQUATION: SINE PART
603  // ----------------------------------------------
604 
605  case 5:
606 
607  // Pressure gradient term
608  sum += k*interpolated_PC*testf_*W*hang_weight;
609 
610  // Stress tensor terms
611  sum += visc_ratio*
612  (-r*interpolated_dVSdr
613  + k*Gamma[0]*interpolated_UC
614  + Gamma[0]*interpolated_VS)*dtestfdr*W*hang_weight;
615 
616  sum += visc_ratio*
617  (k*Gamma[0]*interpolated_WC - r*interpolated_dVSdz)*dtestfdz*W*hang_weight;
618 
619  sum += visc_ratio*
620  (Gamma[0]*interpolated_dVSdr
621  - k*(2.0+Gamma[0])*interpolated_UC/r
622  - (1.0+(k*k)+(k*k*Gamma[0]))*interpolated_VS/r)*testf_*W*hang_weight;
623 
624  // Inertial terms (du/dt)
625  sum -= scaled_re_st*r*dudt[5]*testf_*W*hang_weight;
626 
627  // Inertial terms (convective)
628  sum -=
629  scaled_re*(r*interpolated_ur*interpolated_dVSdr
630  + r*interpolated_duthetadr*interpolated_US
631  - k*interpolated_utheta*interpolated_VC
632  + interpolated_utheta*interpolated_US
633  + interpolated_ur*interpolated_VS
634  + r*interpolated_uz*interpolated_dVSdz
635  + r*interpolated_duthetadz*interpolated_WS)*testf_*W*hang_weight;
636 
637  // Mesh velocity terms
638  if(!ALE_is_disabled)
639  {
640  for(unsigned j=0;j<2;j++)
641  {
642  sum +=
643  scaled_re_st*r*mesh_velocity[j]
644  *interpolated_dudx(5,j)*testf_*W*hang_weight;
645  }
646  }
647 
648  break;
649 
650  } // End of switch statement for momentum equations
651 
652  // Add contribution to elemental residual vector
653  residuals[local_eqn] += sum;
654 
655  // ======================
656  // Calculate the Jacobian
657  // ======================
658  if(flag)
659  {
660  // Number of master nodes
661  unsigned n_master2 = 1;
662 
663  // Storage for the weight of the shape function
664  double hang_weight2 = 1.0;
665 
666  // Loop over the velocity nodes for columns
667  for(unsigned l2=0;l2<n_node;l2++)
668  {
669  // Cache velocity shape functions and their derivatives
670  const double psif_ = psif[l2];
671  const double dpsifdr = dpsifdx(l2,0);
672  const double dpsifdz = dpsifdx(l2,1);
673 
674  // Local boolean that indicates the hanging status of the node
675  bool is_node2_hanging = node_pt(l2)->is_hanging();
676 
677  // If the node is hanging
678  if(is_node2_hanging)
679  {
680  // Get the hanging pointer
681  hang_info2_pt = node_pt(l2)->hanging_pt();
682 
683  // Read out number of master nodes from hanging data
684  n_master2 = hang_info2_pt->nmaster();
685  }
686  // Otherwise the node is its own master
687  else
688  {
689  n_master2 = 1;
690  }
691 
692  // Loop over the master nodes
693  for(unsigned m2=0;m2<n_master2;m2++)
694  {
695  // Loop over the velocity components
696  for(unsigned i2=0;i2<6;i2++)
697  {
698  // Get the number of the unknown
699  // -----------------------------
700 
701  // If the node is hanging
702  if(is_node2_hanging)
703  {
704  // Get the equation number from the master node
705  local_unknown =
706  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
707  u_nodal_index[i2]);
708 
709  // Get the hang weight from the master node
710  hang_weight2 = hang_info2_pt->master_weight(m2);
711  }
712  // If the node is not hanging
713  else
714  {
715  // Local equation number
716  local_unknown =
717  this->nodal_local_eqn(l2,u_nodal_index[i2]);
718 
719  // Node contributes with full weight
720  hang_weight2 = 1.0;
721  }
722 
723  // If the unknown is not pinned
724  if(local_unknown >= 0)
725  {
726  // Different results for i and i2
727  switch(i)
728  {
729 
730  // ---------------------------------------------
731  // FIRST (RADIAL) MOMENTUM EQUATION: COSINE PART
732  // ---------------------------------------------
733 
734  case 0:
735 
736  switch(i2)
737  {
738 
739  // Radial velocity component (cosine part) U_k^C
740 
741  case 0:
742 
743  // Add the mass matrix entries
744  if(flag==2)
745  {
746  mass_matrix(local_eqn,local_unknown) +=
747  scaled_re_st*r*psif_*testf_*W*hang_weight*hang_weight2
748  *hang_weight*hang_weight2;
749  }
750 
751  // Add contributions to the Jacobian matrix
752 
753  // Stress tensor terms
754  jacobian(local_eqn,local_unknown)
755  -= visc_ratio*r*(1.0+Gamma[0])*dpsifdr*dtestfdr*W*hang_weight*hang_weight2;
756 
757  jacobian(local_eqn,local_unknown)
758  -= visc_ratio*r*dpsifdz*dtestfdz*W*hang_weight*hang_weight2;
759 
760  jacobian(local_eqn,local_unknown)
761  -= visc_ratio*(1.0+Gamma[0]+(k*k))*psif_*testf_*W*hang_weight*hang_weight2/r;
762 
763  // Inertial terms (du/dt)
764  jacobian(local_eqn,local_unknown) -=
765  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
766  psif_*testf_*W*hang_weight*hang_weight2;
767 
768  // Inertial terms (convective)
769  jacobian(local_eqn,local_unknown) -=
770  scaled_re*r*(psif_*interpolated_durdr
771  + interpolated_ur*dpsifdr
772  + interpolated_uz*dpsifdz)*testf_*W*hang_weight*hang_weight2;
773 
774  // Mesh velocity terms
775  if(!ALE_is_disabled)
776  {
777  for(unsigned j=0;j<2;j++)
778  {
779  jacobian(local_eqn,local_unknown) +=
780  scaled_re_st*r*mesh_velocity[j]
781  *dpsifdx(l2,j)*testf_*W*hang_weight*hang_weight2;
782  }
783  }
784 
785  break;
786 
787  // Radial velocity component (sine part) U_k^S
788 
789  case 1:
790 
791  // Convective term
792  jacobian(local_eqn,local_unknown) -=
793  scaled_re*k*interpolated_utheta*psif_*testf_*W*hang_weight*hang_weight2;
794 
795  break;
796 
797  // Axial velocity component (cosine part) W_k^C
798 
799  case 2:
800 
801  // Stress tensor term
802  jacobian(local_eqn,local_unknown) -=
803  visc_ratio*Gamma[0]*r*dpsifdr*dtestfdz*W*hang_weight*hang_weight2;
804 
805  // Convective term
806  jacobian(local_eqn,local_unknown) -=
807  scaled_re*r*interpolated_durdz*psif_*testf_*W*hang_weight*hang_weight2;
808 
809  break;
810 
811  // Axial velocity component (sine part) W_k^S
812  // has no contribution
813 
814  case 3:
815  break;
816 
817  // Azimuthal velocity component (cosine part) V_k^C
818 
819  case 4:
820 
821  // Convective term
822  jacobian(local_eqn,local_unknown) +=
823  scaled_re*2*interpolated_utheta*psif_*testf_*W*hang_weight*hang_weight2;
824 
825  break;
826 
827  // Azimuthal velocity component (sine part) V_k^S
828 
829  case 5:
830 
831  // Stress tensor term
832  jacobian(local_eqn,local_unknown) += visc_ratio*(
833  (Gamma[0]*k*dpsifdr) - (k*(2.0+Gamma[0])*psif_/r))*testf_*W*hang_weight*hang_weight2;
834 
835  break;
836 
837  } // End of first (radial) momentum equation
838  break;
839 
840  // --------------------------------------------
841  // SECOND (RADIAL) MOMENTUM EQUATION: SINE PART
842  // --------------------------------------------
843 
844  case 1:
845 
846  switch(i2)
847  {
848 
849  // Radial velocity component (cosine part) U_k^C
850 
851  case 0:
852 
853  // Convective term
854  jacobian(local_eqn,local_unknown) +=
855  scaled_re*k*interpolated_utheta*psif_*testf_*W*hang_weight*hang_weight2;
856 
857  break;
858 
859  // Radial velocity component (sine part) U_k^S
860 
861  case 1:
862 
863  // Add the mass matrix entries
864  if(flag==2)
865  {
866  mass_matrix(local_eqn,local_unknown) +=
867  scaled_re_st*r*psif_*testf_*W*hang_weight*hang_weight2;
868  }
869 
870  // Add contributions to the Jacobian matrix
871 
872  // Stress tensor terms
873  jacobian(local_eqn,local_unknown)
874  -= visc_ratio*r*(1.0+Gamma[0])*dpsifdr*dtestfdr*W*hang_weight*hang_weight2;
875 
876  jacobian(local_eqn,local_unknown)
877  -= visc_ratio*r*dpsifdz*dtestfdz*W*hang_weight*hang_weight2;
878 
879  jacobian(local_eqn,local_unknown)
880  -= visc_ratio*(1.0+Gamma[0]+(k*k))*psif_*testf_*W*hang_weight*hang_weight2/r;
881 
882  // Inertial terms (du/dt)
883  jacobian(local_eqn,local_unknown)
884  -= scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
885  psif_*testf_*W*hang_weight*hang_weight2;
886 
887  // Inertial terms (convective)
888  jacobian(local_eqn,local_unknown) -=
889  scaled_re*r*(psif_*interpolated_durdr
890  + interpolated_ur*dpsifdr
891  + interpolated_uz*dpsifdz)*testf_*W*hang_weight*hang_weight2;
892 
893  // Mesh velocity terms
894  if(!ALE_is_disabled)
895  {
896  for(unsigned j=0;j<2;j++)
897  {
898  jacobian(local_eqn,local_unknown) +=
899  scaled_re_st*r*mesh_velocity[j]
900  *dpsifdx(l2,j)*testf_*W*hang_weight*hang_weight2;
901  }
902  }
903 
904  break;
905 
906  // Axial velocity component (cosine part) W_k^C
907  // has no contribution
908 
909  case 2:
910  break;
911 
912  // Axial velocity component (sine part) W_k^S
913 
914  case 3:
915 
916  // Stress tensor term
917  jacobian(local_eqn,local_unknown) -=
918  visc_ratio*Gamma[0]*r*dpsifdr*dtestfdz*W*hang_weight*hang_weight2;
919 
920  // Convective term
921  jacobian(local_eqn,local_unknown) -=
922  scaled_re*r*interpolated_durdz*psif_*testf_*W*hang_weight*hang_weight2;
923 
924  break;
925 
926  // Azimuthal velocity component (cosine part) V_k^C
927 
928  case 4:
929 
930  // Stress tensor terms
931  jacobian(local_eqn,local_unknown) -= visc_ratio*(
932  (Gamma[0]*k*dpsifdr)
933  - (k*(2.0+Gamma[0])*psif_/r))*testf_*W*hang_weight*hang_weight2;
934 
935  break;
936 
937  // Azimuthal velocity component (sine part) V_k^S
938 
939  case 5:
940 
941  // Convective term
942  jacobian(local_eqn,local_unknown) +=
943  scaled_re*2*interpolated_utheta*psif_*testf_*W*hang_weight*hang_weight2;
944 
945  break;
946 
947  } // End of second (radial) momentum equation
948  break;
949 
950  // --------------------------------------------
951  // THIRD (AXIAL) MOMENTUM EQUATION: COSINE PART
952  // --------------------------------------------
953 
954  case 2:
955 
956  switch(i2)
957  {
958 
959  // Radial velocity component (cosine part) U_k^C
960 
961  case 0:
962 
963  // Stress tensor term
964  jacobian(local_eqn,local_unknown) -=
965  visc_ratio*r*Gamma[1]*dpsifdz*dtestfdr*W*hang_weight*hang_weight2;
966 
967  // Convective term
968  jacobian(local_eqn,local_unknown) -=
969  scaled_re*r*psif_*interpolated_duzdr*testf_*W*hang_weight*hang_weight2;
970 
971  break;
972 
973  // Radial velocity component (sine part) U_k^S
974  // has no contribution
975 
976  case 1:
977  break;
978 
979  // Axial velocity component (cosine part) W_k^C
980 
981  case 2:
982 
983  // Add the mass matrix entries
984  if(flag==2)
985  {
986  mass_matrix(local_eqn,local_unknown) +=
987  scaled_re_st*r*psif_*testf_*W*hang_weight*hang_weight2;
988  }
989 
990  // Add contributions to the Jacobian matrix
991 
992  // Stress tensor terms
993  jacobian(local_eqn,local_unknown) -=
994  visc_ratio*r*dpsifdr*dtestfdr*W*hang_weight*hang_weight2;
995 
996  jacobian(local_eqn,local_unknown) -=
997  visc_ratio*r*(1.0+Gamma[1])*dpsifdz*dtestfdz*W*hang_weight*hang_weight2;
998 
999  jacobian(local_eqn,local_unknown) -=
1000  visc_ratio*k*k*psif_*testf_*W*hang_weight*hang_weight2/r;
1001 
1002  // Inertial terms (du/dt)
1003  jacobian(local_eqn,local_unknown) -=
1004  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1005  psif_*testf_*W*hang_weight*hang_weight2;
1006 
1007  // Inertial terms (convective)
1008  jacobian(local_eqn,local_unknown) -=
1009  scaled_re*r*(interpolated_ur*dpsifdr
1010  + psif_*interpolated_duzdz
1011  + interpolated_uz*dpsifdz)*testf_*W*hang_weight*hang_weight2;
1012 
1013 
1014  // Mesh velocity terms
1015  if(!ALE_is_disabled)
1016  {
1017  for(unsigned j=0;j<2;j++)
1018  {
1019  jacobian(local_eqn,local_unknown) +=
1020  scaled_re_st*r*mesh_velocity[j]
1021  *dpsifdx(l2,j)*testf_*W*hang_weight*hang_weight2;
1022  }
1023  }
1024 
1025  break;
1026 
1027  // Axial velocity component (sine part) W_k^S
1028 
1029  case 3:
1030 
1031  // Convective term
1032  jacobian(local_eqn,local_unknown) -=
1033  scaled_re*k*interpolated_utheta*psif_*testf_*W*hang_weight*hang_weight2;
1034 
1035  break;
1036 
1037  // Azimuthal velocity component (cosine part) V_k^C
1038  // has no contribution
1039 
1040  case 4:
1041  break;
1042 
1043  // Azimuthal velocity component (sine part) V_k^S
1044 
1045  case 5:
1046 
1047  // Stress tensor term
1048  jacobian(local_eqn,local_unknown) +=
1049  visc_ratio*Gamma[1]*k*dpsifdz*testf_*W*hang_weight*hang_weight2;
1050 
1051  break;
1052 
1053  } // End of third (axial) momentum equation
1054  break;
1055 
1056  // -------------------------------------------
1057  // FOURTH (AXIAL) MOMENTUM EQUATION: SINE PART
1058  // -------------------------------------------
1059 
1060  case 3:
1061 
1062  switch(i2)
1063  {
1064 
1065  // Radial velocity component (cosine part) U_k^C
1066  // has no contribution
1067 
1068  case 0:
1069  break;
1070 
1071  // Radial velocity component (sine part) U_k^S
1072 
1073  case 1:
1074 
1075  // Stress tensor term
1076  jacobian(local_eqn,local_unknown) -=
1077  visc_ratio*r*Gamma[1]*dpsifdz*dtestfdr*W*hang_weight*hang_weight2;
1078 
1079  // Convective term
1080  jacobian(local_eqn,local_unknown) -=
1081  scaled_re*r*psif_*interpolated_duzdr*testf_*W*hang_weight*hang_weight2;
1082 
1083  break;
1084 
1085  // Axial velocity component (cosine part) W_k^S
1086 
1087  case 2:
1088 
1089  // Convective term
1090  jacobian(local_eqn,local_unknown) +=
1091  scaled_re*k*interpolated_utheta*psif_*testf_*W*hang_weight*hang_weight2;
1092 
1093  break;
1094 
1095  // Axial velocity component (sine part) W_k^S
1096 
1097  case 3:
1098 
1099  // Add the mass matrix entries
1100  if(flag==2)
1101  {
1102  mass_matrix(local_eqn,local_unknown) +=
1103  scaled_re_st*r*psif_*testf_*W*hang_weight*hang_weight2;
1104  }
1105 
1106  // Add contributions to the Jacobian matrix
1107 
1108  // Stress tensor terms
1109  jacobian(local_eqn,local_unknown) -=
1110  visc_ratio*r*dpsifdr*dtestfdr*W*hang_weight*hang_weight2;
1111 
1112  jacobian(local_eqn,local_unknown) -=
1113  visc_ratio*r*(1.0+Gamma[1])*dpsifdz*dtestfdz*W*hang_weight*hang_weight2;
1114 
1115  jacobian(local_eqn,local_unknown) -=
1116  visc_ratio*k*k*psif_*testf_*W*hang_weight*hang_weight2/r;
1117 
1118  // Inertial terms (du/dt)
1119  jacobian(local_eqn,local_unknown) -=
1120  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1121  psif_*testf_*W*hang_weight*hang_weight2;
1122 
1123  // Inertial terms (convective)
1124  jacobian(local_eqn,local_unknown) -=
1125  scaled_re*r*(interpolated_ur*dpsifdr
1126  + psif_*interpolated_duzdz
1127  + interpolated_uz*dpsifdz)*testf_*W*hang_weight*hang_weight2;
1128 
1129 
1130  // Mesh velocity terms
1131  if(!ALE_is_disabled)
1132  {
1133  for(unsigned j=0;j<2;j++)
1134  {
1135  jacobian(local_eqn,local_unknown) +=
1136  scaled_re_st*r*mesh_velocity[j]
1137  *dpsifdx(l2,j)*testf_*W*hang_weight*hang_weight2;
1138  }
1139  }
1140 
1141  break;
1142 
1143  // Azimuthal velocity component (cosine part) V_k^C
1144 
1145  case 4:
1146 
1147  // Stress tensor term
1148  jacobian(local_eqn,local_unknown) -=
1149  visc_ratio*Gamma[1]*k*dpsifdz*testf_*W*hang_weight*hang_weight2;
1150 
1151  break;
1152 
1153  // Azimuthal velocity component (sine part) V_k^S
1154  // has no contribution
1155 
1156  case 5:
1157  break;
1158 
1159  } // End of fourth (axial) momentum equation
1160  break;
1161 
1162  // ------------------------------------------------
1163  // FIFTH (AZIMUTHAL) MOMENTUM EQUATION: COSINE PART
1164  // ------------------------------------------------
1165 
1166  case 4:
1167 
1168  switch(i2)
1169  {
1170 
1171  // Radial velocity component (cosine part) U_k^C
1172 
1173  case 0:
1174 
1175  // Convective terms
1176  jacobian(local_eqn,local_unknown) -=
1177  scaled_re*(r*interpolated_duthetadr
1178  + interpolated_utheta)*psif_*testf_*W*hang_weight*hang_weight2;
1179 
1180  break;
1181 
1182  // Radial velocity component (sine part) U_k^S
1183 
1184  case 1:
1185 
1186  // Stress tensor terms
1187  jacobian(local_eqn,local_unknown) += visc_ratio*k*psif_*(
1188  ((2.0+Gamma[0])*testf_/r) - (Gamma[0]*dtestfdr))*W*hang_weight*hang_weight2;
1189 
1190  break;
1191 
1192  // Axial velocity component (cosine part) W_k^C
1193 
1194  case 2:
1195 
1196  // Convective term
1197  jacobian(local_eqn,local_unknown) -=
1198  scaled_re*r*psif_*interpolated_duthetadz*testf_*W*hang_weight*hang_weight2;
1199 
1200  break;
1201 
1202  // Axial velocity component (sine part) W_k^S
1203 
1204  case 3:
1205 
1206  // Stress tensor term
1207  jacobian(local_eqn,local_unknown) -=
1208  visc_ratio*k*Gamma[0]*psif_*dtestfdz*W*hang_weight*hang_weight2;
1209 
1210  break;
1211 
1212  // Azimuthal velocity component (cosine part) V_k^C
1213 
1214  case 4:
1215 
1216  // Add the mass matrix entries
1217  if(flag==2)
1218  {
1219  mass_matrix(local_eqn,local_unknown) +=
1220  scaled_re_st*r*psif_*testf_*W*hang_weight*hang_weight2;
1221  }
1222 
1223  // Add contributions to the Jacobian matrix
1224 
1225  // Stress tensor terms
1226  jacobian(local_eqn,local_unknown) +=
1227  visc_ratio*(Gamma[0]*psif_ - r*dpsifdr)*dtestfdr*W*hang_weight*hang_weight2;
1228 
1229  jacobian(local_eqn,local_unknown) -=
1230  visc_ratio*r*dpsifdz*dtestfdz*W*hang_weight*hang_weight2;
1231 
1232  jacobian(local_eqn,local_unknown) +=
1233  visc_ratio*(Gamma[0]*dpsifdr
1234  - (1.0+(k*k)+(k*k*Gamma[0]))*psif_/r)*testf_*W*hang_weight*hang_weight2;
1235 
1236  // Inertial terms (du/dt)
1237  jacobian(local_eqn,local_unknown) -=
1238  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1239  psif_*testf_*W*hang_weight*hang_weight2;
1240 
1241  // Inertial terms (convective)
1242  jacobian(local_eqn,local_unknown) -=
1243  scaled_re*(r*interpolated_ur*dpsifdr
1244  + interpolated_ur*psif_
1245  + r*interpolated_uz*dpsifdz)*testf_*W*hang_weight*hang_weight2;
1246 
1247  // Mesh velocity terms
1248  if(!ALE_is_disabled)
1249  {
1250  for(unsigned j=0;j<2;j++)
1251  {
1252  jacobian(local_eqn,local_unknown) +=
1253  scaled_re_st*r*mesh_velocity[j]
1254  *dpsifdx(l2,j)*testf_*W*hang_weight*hang_weight2;
1255  }
1256  }
1257 
1258  break;
1259 
1260  // Azimuthal velocity component (sine part) V_k^S
1261 
1262  case 5:
1263 
1264  // Convective term
1265  jacobian(local_eqn,local_unknown) -=
1266  scaled_re*k*interpolated_utheta*psif_*testf_*W*hang_weight*hang_weight2;
1267 
1268  break;
1269 
1270  } // End of fifth (azimuthal) momentum equation
1271  break;
1272 
1273  // ----------------------------------------------
1274  // SIXTH (AZIMUTHAL) MOMENTUM EQUATION: SINE PART
1275  // ----------------------------------------------
1276 
1277  case 5:
1278 
1279  switch(i2)
1280  {
1281 
1282  // Radial velocity component (cosine part) U_k^C
1283 
1284  case 0:
1285 
1286  // Stress tensor terms
1287  jacobian(local_eqn,local_unknown) += visc_ratio*k*psif_*(
1288  (Gamma[0]*dtestfdr) - ((2.0+Gamma[0])*testf_/r))*W*hang_weight*hang_weight2;
1289 
1290  break;
1291 
1292  // Radial velocity component (sine part) U_k^S
1293 
1294  case 1:
1295 
1296  // Convective terms
1297  jacobian(local_eqn,local_unknown) -=
1298  scaled_re*(r*interpolated_duthetadr
1299  + interpolated_utheta)*psif_*testf_*W*hang_weight*hang_weight2;
1300 
1301  break;
1302 
1303  // Axial velocity component (cosine part) W_k^C
1304 
1305  case 2:
1306 
1307  // Stress tensor term
1308  jacobian(local_eqn,local_unknown) +=
1309  visc_ratio*k*Gamma[0]*psif_*dtestfdz*W*hang_weight*hang_weight2;
1310 
1311  break;
1312 
1313  // Axial velocity component (sine part) W_k^S
1314 
1315  case 3:
1316 
1317  // Convective term
1318  jacobian(local_eqn,local_unknown) -=
1319  scaled_re*r*psif_*interpolated_duthetadz*testf_*W*hang_weight*hang_weight2;
1320 
1321  break;
1322 
1323  // Azimuthal velocity component (cosine part) V_k^C
1324 
1325  case 4:
1326 
1327  // Convective term
1328  jacobian(local_eqn,local_unknown) +=
1329  scaled_re*k*interpolated_utheta*psif_*testf_*W*hang_weight*hang_weight2;
1330 
1331  break;
1332 
1333  // Azimuthal velocity component (sine part) V_k^S
1334 
1335  case 5:
1336 
1337  // Add the mass matrix entries
1338  if(flag==2)
1339  {
1340  mass_matrix(local_eqn,local_unknown) +=
1341  scaled_re_st*r*psif_*testf_*W*hang_weight*hang_weight2;
1342  }
1343 
1344  // Add contributions to the Jacobian matrix
1345 
1346  // Stress tensor terms
1347  jacobian(local_eqn,local_unknown) +=
1348  visc_ratio*(Gamma[0]*psif_ - r*dpsifdr)*dtestfdr*W*hang_weight*hang_weight2;
1349 
1350  jacobian(local_eqn,local_unknown) -=
1351  visc_ratio*r*dpsifdz*dtestfdz*W*hang_weight*hang_weight2;
1352 
1353  jacobian(local_eqn,local_unknown) +=
1354  visc_ratio*(Gamma[0]*dpsifdr
1355  - (1.0+(k*k)+(k*k*Gamma[0]))*psif_/r)*testf_*W*hang_weight*hang_weight2;
1356 
1357  // Inertial terms (du/dt)
1358  jacobian(local_eqn,local_unknown) -=
1359  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1360  psif_*testf_*W*hang_weight*hang_weight2;
1361 
1362  // Convective terms
1363  jacobian(local_eqn,local_unknown) -=
1364  scaled_re*(r*interpolated_ur*dpsifdr
1365  + interpolated_ur*psif_
1366  + r*interpolated_uz*dpsifdz)*testf_*W*hang_weight*hang_weight2;
1367 
1368  // Mesh velocity terms
1369  if(!ALE_is_disabled)
1370  {
1371  for(unsigned j=0;j<2;j++)
1372  {
1373  jacobian(local_eqn,local_unknown) +=
1374  scaled_re_st*r*mesh_velocity[j]
1375  *dpsifdx(l2,j)*testf_*W*hang_weight*hang_weight2;
1376  }
1377  }
1378 
1379  break;
1380 
1381  } // End of sixth (azimuthal) momentum equation
1382  break;
1383  }
1384  } // End of if not boundary condition statement
1385  } // End of loop over velocity components
1386  } // End of loop over master (m2) nodes
1387  } // End of loop over the velocity nodes (l2)
1388 
1389 
1390  // Loop over the pressure shape functions
1391  for(unsigned l2=0;l2<n_pres;l2++)
1392  {
1393  // If the pressure dof is hanging
1394  if(pressure_dof_is_hanging[l2])
1395  {
1396  // Pressure dof is hanging so it must be nodal-based
1397  hang_info2_pt =
1398  pressure_node_pt(l2)->hanging_pt(p_index[0]);
1399 
1400  // Get the number of master nodes from the pressure node
1401  n_master2 = hang_info2_pt->nmaster();
1402  }
1403  // Otherwise the node is its own master
1404  else
1405  {
1406  n_master2 = 1;
1407  }
1408 
1409  // Loop over the master nodes
1410  for(unsigned m2=0;m2<n_master2;m2++)
1411  {
1412  // Loop over the two pressure components
1413  for(unsigned i2=0;i2<2;i2++)
1414  {
1415  // Get the number of the unknown
1416  // -----------------------------
1417 
1418  // If the pressure dof is hanging
1419  if(pressure_dof_is_hanging[l2])
1420  {
1421  // Get the unknown from the node
1422  local_unknown =
1423  local_hang_eqn(hang_info2_pt->master_node_pt(m2),
1424  p_index[i2]);
1425 
1426  // Get the weight from the hanging object
1427  hang_weight2 = hang_info2_pt->master_weight(m2);
1428  }
1429  else
1430  {
1431  local_unknown = p_local_eqn(l2,i2);
1432  hang_weight2 = 1.0;
1433  }
1434 
1435  // If the unknown is not pinned
1436  if(local_unknown >= 0)
1437  {
1438  // Add in contributions to momentum equations
1439  switch(i)
1440  {
1441 
1442  // ---------------------------------------------
1443  // FIRST (RADIAL) MOMENTUM EQUATION: COSINE PART
1444  // ---------------------------------------------
1445 
1446  case 0:
1447 
1448  switch(i2)
1449  {
1450  // Cosine part P_k^C
1451  case 0:
1452 
1453  jacobian(local_eqn,local_unknown) +=
1454  psip[l2]*(testf_ + r*dtestfdr)*W*hang_weight*hang_weight2;
1455 
1456  break;
1457 
1458  // Sine part P_k^S has no contribution
1459  case 1:
1460  break;
1461 
1462  } // End of first (radial) momentum equation
1463  break;
1464 
1465  // --------------------------------------------
1466  // SECOND (RADIAL) MOMENTUM EQUATION: SINE PART
1467  // --------------------------------------------
1468 
1469  case 1:
1470 
1471  switch(i2)
1472  {
1473  // Cosine part P_k^C has no contribution
1474  case 0:
1475  break;
1476 
1477  // Sine part P_k^S
1478  case 1:
1479 
1480  jacobian(local_eqn,local_unknown)
1481  += psip[l2]*(testf_ + r*dtestfdr)*W*hang_weight*hang_weight2;
1482 
1483  break;
1484 
1485  } // End of second (radial) momentum equation
1486  break;
1487 
1488  // --------------------------------------------
1489  // THIRD (AXIAL) MOMENTUM EQUATION: COSINE PART
1490  // --------------------------------------------
1491 
1492  case 2:
1493 
1494  switch(i2)
1495  {
1496  // Cosine part P_k^C
1497  case 0:
1498 
1499  jacobian(local_eqn,local_unknown) +=
1500  r*psip[l2]*dtestfdz*W*hang_weight*hang_weight2;
1501 
1502  break;
1503 
1504  // Sine part P_k^S has no contribution
1505  case 1:
1506  break;
1507 
1508  } // End of third (axial) momentum equation
1509  break;
1510 
1511  // -------------------------------------------
1512  // FOURTH (AXIAL) MOMENTUM EQUATION: SINE PART
1513  // -------------------------------------------
1514 
1515  case 3:
1516 
1517  switch(i2)
1518  {
1519  // Cosine part P_k^C has no contribution
1520  case 0:
1521  break;
1522 
1523  // Sine part P_k^S
1524  case 1:
1525 
1526  jacobian(local_eqn,local_unknown) +=
1527  r*psip[l2]*dtestfdz*W*hang_weight*hang_weight2;
1528 
1529  break;
1530 
1531  } // End of fourth (axial) momentum equation
1532  break;
1533 
1534  // ------------------------------------------------
1535  // FIFTH (AZIMUTHAL) MOMENTUM EQUATION: COSINE PART
1536  // ------------------------------------------------
1537 
1538  case 4:
1539 
1540  switch(i2)
1541  {
1542  // Cosine part P_k^C has no contribution
1543  case 0:
1544  break;
1545 
1546  // Sine part P_k^S
1547  case 1:
1548 
1549  jacobian(local_eqn,local_unknown) -=
1550  k*psip[l2]*testf_*W*hang_weight*hang_weight2;
1551 
1552  break;
1553 
1554  } // End of fifth (azimuthal) momentum equation
1555  break;
1556 
1557  // ----------------------------------------------
1558  // SIXTH (AZIMUTHAL) MOMENTUM EQUATION: SINE PART
1559  // ----------------------------------------------
1560 
1561  case 5:
1562 
1563  switch(i2)
1564  {
1565  // Cosine part P_k^C
1566  case 0:
1567 
1568  jacobian(local_eqn,local_unknown) +=
1569  k*psip[l2]*testf_*W*hang_weight*hang_weight2;
1570 
1571  break;
1572 
1573  // Sine part P_k^S has no contribution
1574  case 1:
1575  break;
1576  } // End of sixth (azimuthal) momentum equation
1577  } // End of add in contributions to different equations
1578  } // End of if unknown is not pinned statement
1579  } // End of loop over pressure components
1580  } // End of loop over master (m2) nodes
1581  } // End of loop over pressure "nodes"
1582  }// End of Jacobian calculation
1583  } // End of if not boundary condition statement
1584  } // End of loop over velocity components
1585  } // End of loop over master (m) nodes
1586  } // End of loop over nodes
1587 
1588  // ====================
1589  // CONTINUITY EQUATIONS
1590  // ====================
1591 
1592  // Loop over the pressure shape functions
1593  for(unsigned l=0;l<n_pres;l++)
1594  {
1595  // Cache test function
1596  const double testp_ = testp[l];
1597 
1598  // If the pressure dof is hanging
1599  if(pressure_dof_is_hanging[l])
1600  {
1601  // Pressure dof is hanging so it must be nodal-based
1602  hang_info_pt = pressure_node_pt(l)->hanging_pt(p_index[0]);
1603 
1604  //Get the number of master nodes from the pressure node
1605  n_master = hang_info_pt->nmaster();
1606  }
1607  // Otherwise the node is its own master
1608  else
1609  {
1610  n_master = 1;
1611  }
1612 
1613  // Loop over the master nodes
1614  for(unsigned m=0;m<n_master;m++)
1615  {
1616  // Loop over the two pressure components
1617  for(unsigned i=0;i<2;i++)
1618  {
1619  // Get the equation number
1620  // -----------------------
1621 
1622  // If the pressure dof is hanging
1623  if(pressure_dof_is_hanging[l])
1624  {
1625  // Get the equation number from the master node
1626  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1627  p_index[i]);
1628 
1629  // Get the hang weight from the master node
1630  hang_weight = hang_info_pt->master_weight(m);
1631  }
1632  else
1633  {
1634  // Local equation number
1635  local_eqn = this->p_local_eqn(l,i);
1636 
1637  // Node contributes with full weight
1638  hang_weight = 1.0;
1639  }
1640 
1641  // If it's not a boundary condition...
1642  if(local_eqn >= 0)
1643  {
1644 
1645  switch(i)
1646  {
1647 
1648  // --------------------------------------
1649  // FIRST CONTINUITY EQUATION: COSINE PART
1650  // --------------------------------------
1651 
1652  case 0:
1653 
1654  // Gradient terms
1655  residuals[local_eqn] +=
1656  (interpolated_UC + r*interpolated_dUCdr
1657  + k*interpolated_VS + r*interpolated_dWCdz)*testp_*W*hang_weight;
1658 
1659  break;
1660 
1661  // -------------------------------------
1662  // SECOND CONTINUITY EQUATION: SINE PART
1663  // -------------------------------------
1664 
1665  case 1:
1666 
1667  // Gradient terms
1668  residuals[local_eqn] +=
1669  (interpolated_US + r*interpolated_dUSdr
1670  - k*interpolated_VC + r*interpolated_dWSdz)*testp_*W*hang_weight;
1671 
1672  break;
1673  }
1674 
1675  // ======================
1676  // Calculate the Jacobian
1677  // ======================
1678  if(flag)
1679  {
1680  // Number of master nodes
1681  unsigned n_master2 = 1;
1682 
1683  // Storage for the weight of the shape function
1684  double hang_weight2 = 1.0;
1685 
1686  // Loop over the velocity nodes for columns
1687  for(unsigned l2=0;l2<n_node;l2++)
1688  {
1689  // Cache velocity shape functions and their derivatives
1690  const double psif_ = psif[l2];
1691  const double dpsifdr = dpsifdx(l2,0);
1692  const double dpsifdz = dpsifdx(l2,1);
1693 
1694  // Local boolean to indicate whether the node is hanging
1695  bool is_node2_hanging = node_pt(l2)->is_hanging();
1696 
1697  // If the node is hanging
1698  if(is_node2_hanging)
1699  {
1700  // Get the hanging pointer
1701  hang_info2_pt = node_pt(l2)->hanging_pt();
1702 
1703  // Read out number of master nodes from hanging data
1704  n_master2 = hang_info2_pt->nmaster();
1705  }
1706  // Otherwise the node is its own master
1707  else
1708  {
1709  n_master2 = 1;
1710  }
1711 
1712  // Loop over the master nodes
1713  for(unsigned m2=0;m2<n_master2;m2++)
1714  {
1715  // Loop over the velocity components
1716  for(unsigned i2=0;i2<6;i2++)
1717  {
1718  // Get the number of the unknown
1719  // -----------------------------
1720 
1721  // If the node is hanging
1722  if(is_node2_hanging)
1723  {
1724  // Get the equation number from the master node
1725  local_unknown =
1726  local_hang_eqn(hang_info2_pt->master_node_pt(m2),
1727  u_nodal_index[i2]);
1728 
1729  // Get the hang weight from the master node
1730  hang_weight2 = hang_info2_pt->master_weight(m2);
1731  }
1732  // If the node is not hanging
1733  else
1734  {
1735  // Local equation number
1736  local_unknown =
1737  this->nodal_local_eqn(l2,u_nodal_index[i2]);
1738 
1739  // Node contributes with full weight
1740  hang_weight2 = 1.0;
1741  }
1742 
1743  // If the unknown is not pinned
1744  if(local_unknown >= 0)
1745  {
1746  // Different results for i and i2
1747  switch(i)
1748  {
1749 
1750  // --------------------------------------
1751  // FIRST CONTINUITY EQUATION: COSINE PART
1752  // --------------------------------------
1753 
1754  case 0:
1755 
1756  switch(i2)
1757  {
1758  // Radial velocity component (cosine part) U_k^C
1759 
1760  case 0:
1761 
1762  jacobian(local_eqn,local_unknown) +=
1763  (psif_ + r*dpsifdr)*testp_*W*hang_weight*hang_weight2;
1764 
1765  break;
1766 
1767  // Radial velocity component (sine part) U_k^S
1768  // has no contribution
1769 
1770  case 1:
1771  break;
1772 
1773  // Axial velocity component (cosine part) W_k^C
1774 
1775  case 2:
1776 
1777  jacobian(local_eqn,local_unknown) +=
1778  r*dpsifdz*testp_*W*hang_weight*hang_weight2;
1779 
1780  break;
1781 
1782  // Axial velocity component (sine part) W_k^S
1783  // has no contribution
1784 
1785  case 3:
1786  break;
1787 
1788  // Azimuthal velocity component (cosine part) V_k^C
1789  // has no contribution
1790 
1791  case 4:
1792  break;
1793 
1794  // Azimuthal velocity component (sine part) V_k^S
1795 
1796  case 5:
1797 
1798  jacobian(local_eqn,local_unknown) +=
1799  k*psif_*testp_*W*hang_weight*hang_weight2;
1800 
1801  break;
1802 
1803  } // End of first continuity equation
1804  break;
1805 
1806  // -------------------------------------
1807  // SECOND CONTINUITY EQUATION: SINE PART
1808  // -------------------------------------
1809 
1810  case 1:
1811 
1812  switch(i2)
1813  {
1814  // Radial velocity component (cosine part) U_k^C
1815  // has no contribution
1816 
1817  case 0:
1818  break;
1819 
1820  // Radial velocity component (sine part) U_k^S
1821 
1822  case 1:
1823 
1824  jacobian(local_eqn,local_unknown) +=
1825  (psif_ + r*dpsifdr)*testp_*W*hang_weight*hang_weight2;
1826 
1827  break;
1828 
1829  // Axial velocity component (cosine part) W_k^C
1830  // has no contribution
1831 
1832  case 2:
1833  break;
1834 
1835  // Axial velocity component (sine part) W_k^S
1836 
1837  case 3:
1838 
1839  jacobian(local_eqn,local_unknown) +=
1840  r*dpsifdz*testp_*W*hang_weight*hang_weight2;
1841 
1842  break;
1843 
1844  // Azimuthal velocity component (cosine part) V_k^C
1845 
1846  case 4:
1847 
1848  jacobian(local_eqn,local_unknown) -=
1849  k*psif_*testp_*W*hang_weight*hang_weight2;
1850 
1851  break;
1852 
1853  // Azimuthal velocity component (sine part) V_k^S
1854  // has no contribution
1855 
1856  case 5:
1857  break;
1858 
1859  } // End of second continuity equation
1860  break;
1861  }
1862  } // End of if unknown is not pinned statement
1863  } // End of loop over velocity components
1864  } // End of loop over master (m2) nodes
1865  } // End of loop over velocity nodes
1866 
1867  // Real and imaginary pressure components, P_k^C and P_k^S,
1868  // have no contribution to Jacobian
1869 
1870  } // End of Jacobian calculation
1871  } // End of if not boundary condition statement
1872  } // End of loop over the two pressure components
1873  } // End of loop over master nodes
1874  } // End of loop over pressure nodes
1875 
1876  } // End of loop over integration points
1877 
1878  } // End of fill_in_generic_residual_contribution_linearised_axi_nst(...)
int i
Definition: BiCGSTAB_step_by_step.cpp:9
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
RowVector3d w
Definition: Matrix_resize_int.cpp:3
MatrixType m2(n_dims)
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
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
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317
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.
Definition: elements.h:2333
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
const double & viscosity_ratio() const
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:245
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:226
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:217
double du_dt_linearised_axi_nst(const unsigned &n, const unsigned &i) const
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:289
const double & re() const
Reynolds number.
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:223
virtual void get_base_flow_u(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:143
virtual int p_local_eqn(const unsigned &n, const unsigned &i)=0
bool ALE_is_disabled
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:110
virtual void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
virtual void get_base_flow_dudx(const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:164
virtual double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
virtual unsigned npres_linearised_axi_nst() const =0
const double & density_ratio() const
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:252
virtual unsigned u_index_linearised_axi_nst(const unsigned &i) const
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:284
const int & azimuthal_mode_number() const
Azimuthal mode number k in e^ik(theta) decomposition.
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:229
virtual double p_linearised_axi_nst(const unsigned &n_p, const unsigned &i) const =0
virtual int p_index_linearised_axi_nst(const unsigned &i) const
Which nodal value represents the pressure?
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:335
HangInfo *const & hanging_pt() const
Definition: nodes.h:1228
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Definition: refineable_elements.h:278
virtual Node * pressure_node_pt(const unsigned &n_p)
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/refineable_linearised_axisym_navier_stokes_elements.h:58
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
RealScalar s
Definition: level1_cplx_impl.h:130
int * m
Definition: level2_cplx_impl.h:294
char char char int int * k
Definition: level2_impl.h:374
r
Definition: UniformPSDSelfTest.py:20
@ W
Definition: quadtree.h:63
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::LinearisedAxisymmetricNavierStokesEquations::ALE_is_disabled, oomph::LinearisedAxisymmetricNavierStokesEquations::azimuthal_mode_number(), oomph::LinearisedAxisymmetricNavierStokesEquations::density_ratio(), oomph::FiniteElement::dnodal_position_dt(), oomph::LinearisedAxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(), oomph::LinearisedAxisymmetricNavierStokesEquations::du_dt_linearised_axi_nst(), oomph::LinearisedAxisymmetricNavierStokesEquations::Gamma, oomph::LinearisedAxisymmetricNavierStokesEquations::get_base_flow_dudx(), oomph::LinearisedAxisymmetricNavierStokesEquations::get_base_flow_u(), oomph::Node::hanging_pt(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), oomph::Node::is_hanging(), J, j, k, 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::LinearisedAxisymmetricNavierStokesEquations::npres_linearised_axi_nst(), oomph::Integral::nweight(), oomph::LinearisedAxisymmetricNavierStokesEquations::p_index_linearised_axi_nst(), oomph::LinearisedAxisymmetricNavierStokesEquations::p_linearised_axi_nst(), oomph::LinearisedAxisymmetricNavierStokesEquations::p_local_eqn(), pressure_node_pt(), oomph::LinearisedAxisymmetricNavierStokesEquations::pshape_linearised_axi_nst(), UniformPSDSelfTest::r, oomph::LinearisedAxisymmetricNavierStokesEquations::re(), oomph::LinearisedAxisymmetricNavierStokesEquations::re_st(), s, oomph::Time::time(), oomph::TimeStepper::time_pt(), oomph::Data::time_stepper_pt(), oomph::LinearisedAxisymmetricNavierStokesEquations::u_index_linearised_axi_nst(), oomph::LinearisedAxisymmetricNavierStokesEquations::viscosity_ratio(), w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and oomph::TimeStepper::weight().

◆ fill_in_generic_residual_contribution_linearised_axi_nst() [2/2]

void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_linearised_axi_nst ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
DenseMatrix< double > &  mass_matrix,
unsigned  flag 
)
privatevirtual

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

Reimplemented from oomph::LinearisedAxisymmetricNavierStokesEquations.

◆ further_build() [1/2]

void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::further_build ( )
inlinevirtual

Further build: pass the pointers down to the sons.

Reimplemented from oomph::RefineableElement.

Reimplemented in oomph::RefineableLinearisedAxisymmetricQCrouzeixRaviartElement, and oomph::RefineableLinearisedAxisymmetricQCrouzeixRaviartElement.

133  {
134  // Find the father element
136  cast_father_element_pt
138  (this->father_element_pt());
139 
140  // Set the viscosity ratio pointer
141  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
142  // Set the density ratio pointer
143  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
144  // Set pointer to global Reynolds number
145  this->Re_pt = cast_father_element_pt->re_pt();
146  // Set pointer to global Reynolds number x Strouhal number (=Womersley)
147  this->ReSt_pt = cast_father_element_pt->re_st_pt();
148  // Set pointer to azimuthal mode number
150  cast_father_element_pt->azimuthal_mode_number_pt();
151  // Set pointer to global Reynolds number x inverse Froude number
152 // this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
153  //Set pointer to the global Reynolds number x inverse Rossby number
154 // this->ReInvRo_pt = cast_father_element_pt->re_invro_pt();
155  //Set pointer to global gravity Vector
156 // this->G_pt = cast_father_element_pt->g_pt();
157 
158  //Set pointer to body force function
159 // this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
160 
161  //Set pointer to volumetric source function
162 // this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
163 
164  // Set the ALE_is_disabled flag
165  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
166  }
double * Re_pt
Pointer to global Reynolds number.
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:88
double * Density_Ratio_pt
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:85
int * Azimuthal_Mode_Number_pt
Pointer to azimuthal mode number k in e^ik(theta) decomposition.
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:94
double * Viscosity_Ratio_pt
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:81
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.h:91
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Definition: refineable_elements.h:539
RefineableLinearisedAxisymmetricNavierStokesEquations()
Empty Constructor.
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/refineable_linearised_axisym_navier_stokes_elements.h:70

References oomph::LinearisedAxisymmetricNavierStokesEquations::ALE_is_disabled, oomph::LinearisedAxisymmetricNavierStokesEquations::Azimuthal_Mode_Number_pt, oomph::LinearisedAxisymmetricNavierStokesEquations::azimuthal_mode_number_pt(), oomph::LinearisedAxisymmetricNavierStokesEquations::Density_Ratio_pt, oomph::LinearisedAxisymmetricNavierStokesEquations::density_ratio_pt(), oomph::RefineableElement::father_element_pt(), oomph::LinearisedAxisymmetricNavierStokesEquations::Re_pt, oomph::LinearisedAxisymmetricNavierStokesEquations::re_pt(), oomph::LinearisedAxisymmetricNavierStokesEquations::re_st_pt(), oomph::LinearisedAxisymmetricNavierStokesEquations::ReSt_pt, oomph::LinearisedAxisymmetricNavierStokesEquations::Viscosity_Ratio_pt, and oomph::LinearisedAxisymmetricNavierStokesEquations::viscosity_ratio_pt().

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

◆ further_build() [2/2]

void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::further_build ( )
inlinevirtual

Further build: pass the pointers down to the sons.

Reimplemented from oomph::RefineableElement.

Reimplemented in oomph::RefineableLinearisedAxisymmetricQCrouzeixRaviartElement.

138  {
139  // Find the father element
141  cast_father_element_pt =
143  this->father_element_pt());
144 
145  // Set the viscosity ratio pointer
146  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
147 
148  // Set the density ratio pointer
149  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
150 
151  // Set pointer to global Reynolds number
152  this->Re_pt = cast_father_element_pt->re_pt();
153 
154  // Set pointer to global Reynolds number x Strouhal number (=Womersley)
155  this->ReSt_pt = cast_father_element_pt->re_st_pt();
156 
157  // Set pointer to azimuthal mode number
159  cast_father_element_pt->azimuthal_mode_number_pt();
160 
161  // Set the ALE_is_disabled flag
162  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
163  }

References oomph::LinearisedAxisymmetricNavierStokesEquations::ALE_is_disabled, oomph::LinearisedAxisymmetricNavierStokesEquations::Azimuthal_Mode_Number_pt, oomph::LinearisedAxisymmetricNavierStokesEquations::Density_Ratio_pt, oomph::RefineableElement::father_element_pt(), oomph::LinearisedAxisymmetricNavierStokesEquations::Re_pt, oomph::LinearisedAxisymmetricNavierStokesEquations::ReSt_pt, and oomph::LinearisedAxisymmetricNavierStokesEquations::Viscosity_Ratio_pt.

◆ geometric_jacobian() [1/2]

double oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::geometric_jacobian ( const Vector< double > &  x)
inlinevirtual

Fill in the geometric Jacobian, which in this case is r.

Reimplemented from oomph::ElementWithZ2ErrorEstimator.

128 { return x[0]; }
list x
Definition: plotDoE.py:28

References plotDoE::x.

◆ geometric_jacobian() [2/2]

double oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::geometric_jacobian ( const Vector< double > &  x)
inlinevirtual

Fill in the geometric Jacobian, which in this case is r.

Reimplemented from oomph::ElementWithZ2ErrorEstimator.

132  {
133  return x[0];
134  }

References plotDoE::x.

◆ get_Z2_flux() [1/2]

void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::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.

85  {
86  // Specify the number of velocity dimensions
87  unsigned DIM=3;
88 #ifdef PARANOID
89  unsigned num_entries=DIM+((DIM*DIM)-DIM)/2;
90  if (flux.size()!=num_entries)
91  {
92  std::ostringstream error_message;
93  error_message << "The flux vector has the wrong number of entries, "
94  << flux.size() << ", whereas it should be "
95  << num_entries << std::endl;
96  throw OomphLibError(
97  error_message.str(),
100  }
101 #endif
102  // Get strain rate matrix
103  DenseMatrix<double> strainrate(DIM);
104  this->strain_rate(s,strainrate);
105 
106  // Pack into flux Vector
107  unsigned icount=0;
108 
109  // Start with diagonal terms
110  for(unsigned i=0;i<DIM;i++)
111  {
112  flux[icount]=strainrate(i,i);
113  icount++;
114  }
115 
116  // Off diagonals row by row
117  for(unsigned i=0;i<DIM;i++)
118  {
119  for(unsigned j=i+1;j<DIM;j++)
120  {
121  flux[icount]=strainrate(i,j);
122  icount++;
123  }
124  }
125  }
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/linearised_axisym_navier_stokes_elements.cc:333
#define DIM
Definition: linearised_navier_stokes_elements.h:44
void flux(const double &time, const Vector< double > &x, double &flux)
Get flux applied along boundary x=0.
Definition: pretend_melt.cc:59
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References DIM, ProblemParameters::flux(), i, j, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::LinearisedAxisymmetricNavierStokesEquations::strain_rate().

◆ get_Z2_flux() [2/2]

void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::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.

86  {
87  // Specify the number of velocity dimensions
88  unsigned DIM = 3;
89 
90 #ifdef PARANOID
91  unsigned num_entries = DIM + ((DIM * DIM) - DIM) / 2;
92  if (flux.size() != num_entries)
93  {
94  std::ostringstream error_message;
95  error_message << "The flux vector has the wrong number of entries, "
96  << flux.size() << ", whereas it should be " << num_entries
97  << std::endl;
98  throw OomphLibError(error_message.str(),
99  "RefineableLinearisedAxisymmetricNavierStokesEquati"
100  "ons::get_Z2_flux()",
102  }
103 #endif
104 
105  // Get strain rate matrix
106  DenseMatrix<double> strainrate(DIM);
107  this->strain_rate(s, strainrate);
108 
109  // Pack into flux Vector
110  unsigned icount = 0;
111 
112  // Start with diagonal terms
113  for (unsigned i = 0; i < DIM; i++)
114  {
115  flux[icount] = strainrate(i, i);
116  icount++;
117  }
118 
119  // Off diagonals row by row
120  for (unsigned i = 0; i < DIM; i++)
121  {
122  for (unsigned j = i + 1; j < DIM; j++)
123  {
124  flux[icount] = strainrate(i, j);
125  icount++;
126  }
127  }
128  }

References DIM, ProblemParameters::flux(), i, j, OOMPH_EXCEPTION_LOCATION, and oomph::LinearisedAxisymmetricNavierStokesEquations::strain_rate().

◆ num_Z2_flux_terms() [1/2]

unsigned oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::num_Z2_flux_terms ( )
inlinevirtual

Number of 'flux' terms for Z2 error estimation.

Implements oomph::ElementWithZ2ErrorEstimator.

77  {
78  // 3 diagonal strain rates, 3 off diagonal
79  return 6;
80  }

◆ num_Z2_flux_terms() [2/2]

unsigned oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::num_Z2_flux_terms ( )
inlinevirtual

Number of 'flux' terms for Z2 error estimation.

Implements oomph::ElementWithZ2ErrorEstimator.

78  {
79  // 3 diagonal strain rates, 3 off diagonal
80  return 6;
81  }

◆ pin_elemental_redundant_nodal_pressure_dofs() [1/2]

virtual void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::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::RefineableLinearisedAxisymmetricQTaylorHoodElement, and oomph::RefineableLinearisedAxisymmetricQTaylorHoodElement.

65 {}

Referenced by pin_redundant_nodal_pressures().

◆ pin_elemental_redundant_nodal_pressure_dofs() [2/2]

virtual void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::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::RefineableLinearisedAxisymmetricQTaylorHoodElement, and oomph::RefineableLinearisedAxisymmetricQTaylorHoodElement.

65 {}

◆ pin_redundant_nodal_pressures() [1/2]

static void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::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

  • RefineableLinearisedAxisymmetricNavierStokesEquations:: pin_all_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.)

180  {
181  // Loop over all elements to brutally pin all nodal pressure degrees of
182  // freedom
183  const unsigned n_element=element_pt.size();
184  for(unsigned e=0;e<n_element;e++)
185  {
188  }
189  }
Array< double, 1, 3 > e(1./3., 0.5, 2.)
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Definition: demo_drivers/axisym_navier_stokes/counter_rotating_disks/refineable_linearised_axisym_navier_stokes_elements.h:65

References e(), and pin_elemental_redundant_nodal_pressure_dofs().

◆ pin_redundant_nodal_pressures() [2/2]

static void oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::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

  • RefineableLinearisedAxisymmetricNavierStokesEquations:: pin_all_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.)

177  {
178  // Loop over all elements to brutally pin all nodal pressure degrees of
179  // freedom
180  const unsigned n_element = element_pt.size();
181  for (unsigned e = 0; e < n_element; e++)
182  {
184  element_pt[e])
186  }
187  }

References e(), and pin_elemental_redundant_nodal_pressure_dofs().

◆ pressure_node_pt() [1/2]

virtual Node* oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::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 in oomph::RefineableLinearisedAxisymmetricQTaylorHoodElement, and oomph::RefineableLinearisedAxisymmetricQTaylorHoodElement.

58 {return NULL;}

Referenced by fill_in_generic_residual_contribution_linearised_axi_nst().

◆ pressure_node_pt() [2/2]

virtual Node* oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::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 in oomph::RefineableLinearisedAxisymmetricQTaylorHoodElement, and oomph::RefineableLinearisedAxisymmetricQTaylorHoodElement.

56  {
57  return NULL;
58  }

◆ unpin_all_pressure_dofs() [1/2]

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

Unpin all pressure dofs in elements listed in Vector.

194  {
195  // Loop over all elements to brutally unpin all nodal pressure degrees of
196  // freedom and internal pressure degrees of freedom
197  const unsigned n_element = element_pt.size();
198  for(unsigned e=0;e<n_element;e++)
199  {
201  (element_pt[e])->unpin_elemental_pressure_dofs();
202  }
203  }
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.

References e(), and unpin_elemental_pressure_dofs().

◆ unpin_all_pressure_dofs() [2/2]

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

Unpin all pressure dofs in elements listed in Vector.

192  {
193  // Loop over all elements to brutally unpin all nodal pressure degrees of
194  // freedom and internal pressure degrees of freedom
195  const unsigned n_element = element_pt.size();
196  for (unsigned e = 0; e < n_element; e++)
197  {
199  element_pt[e])
201  }
202  }

References e(), and unpin_elemental_pressure_dofs().

◆ unpin_elemental_pressure_dofs() [1/2]

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

◆ unpin_elemental_pressure_dofs() [2/2]

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

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