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

Refineable version of the Axisymmetric Navier–Stokes equations. More...

#include <generalised_newtonian_refineable_axisym_navier_stokes_elements.h>

+ Inheritance diagram for oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations:

Public Member Functions

 RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations ()
 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...
 
void dinterpolated_u_axi_nst_ddata (const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
 
virtual void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
- Public Member Functions inherited from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations
 GeneralisedNewtonianAxisymmetricNavierStokesEquations ()
 Constructor: NULL the body force and source function. More...
 
const double & re () const
 Reynolds number. More...
 
const double & re_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
const double & re_invfr () const
 Global inverse Froude number. More...
 
double *& re_invfr_pt ()
 Pointer to global inverse Froude number. More...
 
const double & re_invro () const
 Global Reynolds number multiplied by inverse Rossby number. More...
 
double *& re_invro_pt ()
 Pointer to global inverse Froude number. More...
 
const Vector< double > & g () const
 Vector of gravitational components. More...
 
Vector< double > *& g_pt ()
 Pointer to Vector of gravitational components. More...
 
const double & density_ratio () const
 
double *& density_ratio_pt ()
 Pointer to Density ratio. More...
 
const double & viscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to Viscosity Ratio. More...
 
GeneralisedNewtonianConstitutiveEquation< 3 > *& constitutive_eqn_pt ()
 Access function for the constitutive equation pointer. More...
 
void use_current_strainrate_to_compute_second_invariant ()
 Switch to fully implict evaluation of non-Newtonian const eqn. More...
 
void use_extrapolated_strainrate_to_compute_second_invariant ()
 Use extrapolation for non-Newtonian const eqn. More...
 
virtual unsigned npres_axi_nst () const =0
 Function to return number of pressure degrees of freedom. More...
 
virtual unsigned u_index_axi_nst (const unsigned &i) const
 
unsigned u_index_nst (const unsigned &i) const
 
unsigned n_u_nst () const
 
double du_dt_axi_nst (const unsigned &n, const unsigned &i) const
 
void disable_ALE ()
 
void enable_ALE ()
 
virtual double p_axi_nst (const unsigned &n_p) const =0
 
virtual int p_nodal_index_axi_nst () const
 Which nodal value represents the pressure? More...
 
double pressure_integral () const
 Integral of pressure over element. More...
 
void max_and_min_invariant_and_viscosity (double &min_invariant, double &max_invariant, double &min_viscosity, double &max_viscosity) const
 
double dissipation () const
 Return integral of dissipation over element. More...
 
double dissipation (const Vector< double > &s) const
 Return dissipation at local coordinate s. More...
 
double kin_energy () const
 Get integral of kinetic energy over element. More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 
void strain_rate (const unsigned &t, const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 
virtual void extrapolated_strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 
virtual void extrapolated_strain_rate (const unsigned &ipt, DenseMatrix< double > &strain_rate) const
 
void traction (const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
 
void get_pressure_and_velocity_mass_matrix_diagonal (Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
 
unsigned nscalar_paraview () const
 
void scalar_value_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
 
std::string scalar_name_paraview (const unsigned &i) const
 
void point_output_data (const Vector< double > &s, Vector< double > &data)
 
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 output_fct (std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 
void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
void fill_in_contribution_to_residuals (Vector< double > &residuals)
 Compute the element's residual Vector. More...
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void fill_in_contribution_to_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
void fill_in_contribution_to_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 Compute the element's residual Vector. More...
 
void fill_in_contribution_to_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
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)
 
void interpolated_u_axi_nst (const Vector< double > &s, Vector< double > &veloc) const
 Compute vector of FE interpolated velocity u at local coordinate s. More...
 
double interpolated_u_axi_nst (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated velocity u[i] at local coordinate s. More...
 
double interpolated_u_axi_nst (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 Return FE interpolated velocity u[i] at local coordinate s. More...
 
double interpolated_p_axi_nst (const Vector< double > &s) const
 Return FE interpolated pressure at local coordinate s. More...
 
double interpolated_duds_axi_nst (const Vector< double > &s, const unsigned &i, const unsigned &j) const
 
double interpolated_dudx_axi_nst (const Vector< double > &s, const unsigned &i, const unsigned &j) const
 
double interpolated_dudt_axi_nst (const Vector< double > &s, const unsigned &i) const
 
double interpolated_d_dudx_dX_axi_nst (const Vector< double > &s, const unsigned &p, const unsigned &q, const unsigned &i, const unsigned &k) const
 
double compute_physical_size () const
 Compute the volume of the element. More...
 
- Public Member Functions inherited from oomph::FiniteElement
void set_dimension (const unsigned &dim)
 
void set_nodal_dimension (const unsigned &nodal_dim)
 
void set_nnodal_position_type (const unsigned &nposition_type)
 Set the number of types required to interpolate the coordinate. More...
 
void set_n_node (const unsigned &n)
 
int nodal_local_eqn (const unsigned &n, const unsigned &i) const
 
double dJ_eulerian_at_knot (const unsigned &ipt, Shape &psi, DenseMatrix< double > &djacobian_dX) const
 
 FiniteElement ()
 Constructor. More...
 
virtual ~FiniteElement ()
 
 FiniteElement (const FiniteElement &)=delete
 Broken copy constructor. More...
 
virtual bool local_coord_is_valid (const Vector< double > &s)
 Broken assignment operator. More...
 
virtual void move_local_coord_back_into_element (Vector< double > &s) const
 
void get_centre_of_gravity_and_max_radius_in_terms_of_zeta (Vector< double > &cog, double &max_radius) const
 
virtual void local_coordinate_of_node (const unsigned &j, Vector< double > &s) const
 
virtual void local_fraction_of_node (const unsigned &j, Vector< double > &s_fraction)
 
virtual double local_one_d_fraction_of_node (const unsigned &n1d, const unsigned &i)
 
virtual void set_macro_elem_pt (MacroElement *macro_elem_pt)
 
MacroElement * macro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void describe_nodal_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void assign_all_generic_local_eqn_numbers (const bool &store_local_dof_pt)
 
Node *& node_pt (const unsigned &n)
 Return a pointer to the local node n. More...
 
Node *const & node_pt (const unsigned &n) const
 Return a pointer to the local node n (const version) More...
 
unsigned nnode () const
 Return the number of nodes. More...
 
virtual unsigned nnode_1d () const
 
double raw_nodal_position (const unsigned &n, const unsigned &i) const
 
double raw_nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position (const unsigned &n, const unsigned &i) const
 
double nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double dnodal_position_dt (const unsigned &n, const unsigned &i) const
 Return the i-th component of nodal velocity: dx/dt at local node n. More...
 
double dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual unsigned required_nvalue (const unsigned &n) const
 
unsigned nnodal_position_type () const
 
bool has_hanging_nodes () const
 
unsigned nodal_dimension () const
 Return the required Eulerian dimension of the nodes in this element. More...
 
virtual Node * construct_node (const unsigned &n)
 
virtual Node * construct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual Node * construct_boundary_node (const unsigned &n)
 
virtual Node * construct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
int get_node_number (Node *const &node_pt) const
 
virtual Node * get_node_at_local_coordinate (const Vector< double > &s) const
 
double raw_nodal_value (const unsigned &n, const unsigned &i) const
 
double raw_nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
unsigned dim () const
 
virtual ElementGeometry::ElementGeometry element_geometry () const
 Return the geometry type of the element (either Q or T usually). More...
 
virtual double interpolated_x (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated coordinate x[i] at local coordinate s. More...
 
virtual double interpolated_x (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_x (const Vector< double > &s, Vector< double > &x) const
 Return FE interpolated position x[] at local coordinate s as Vector. More...
 
virtual void interpolated_x (const unsigned &t, const Vector< double > &s, Vector< double > &x) const
 
virtual double interpolated_dxdt (const Vector< double > &s, const unsigned &i, const unsigned &t)
 
virtual void interpolated_dxdt (const Vector< double > &s, const unsigned &t, Vector< double > &dxdt)
 
unsigned ngeom_data () const
 
Data * geom_data_pt (const unsigned &j)
 
void position (const Vector< double > &zeta, Vector< double > &r) const
 
void position (const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
 
void dposition_dt (const Vector< double > &zeta, const unsigned &t, Vector< double > &drdt)
 
virtual double zeta_nodal (const unsigned &n, const unsigned &k, const unsigned &i) const
 
void interpolated_zeta (const Vector< double > &s, Vector< double > &zeta) const
 
void locate_zeta (const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
 
virtual void node_update ()
 
virtual void identify_geometric_data (std::set< Data * > &geometric_data_pt)
 
virtual double s_min () const
 Min value of local coordinate. More...
 
virtual double s_max () const
 Max. value of local coordinate. More...
 
double size () const
 
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 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 void output (const unsigned &t, std::ostream &outfile, const unsigned &n_plot) const
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, const SolutionFunctorBase &exact_soln) const
 Output a time-dependent exact solution over the element. More...
 
virtual void get_s_plot (const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
 
virtual std::string tecplot_zone_string (const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (std::ostream &outfile, const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (FILE *file_pt, const unsigned &nplot) const
 
virtual unsigned nplot_points (const unsigned &nplot) const
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_abs_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error)
 
void integrate_fct (FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt, Vector< double > &integral)
 Integrate Vector-valued function over element. More...
 
void integrate_fct (FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt, const double &time, Vector< double > &integral)
 Integrate Vector-valued time-dep function over element. More...
 
virtual void build_face_element (const int &face_index, FaceElement *face_element_pt)
 
virtual unsigned self_test ()
 
virtual unsigned get_bulk_node_number (const int &face_index, const unsigned &i) const
 
virtual int face_outer_unit_normal_sign (const int &face_index) const
 Get the sign of the outer unit normal on the face given by face_index. More...
 
virtual unsigned nnode_on_face () const
 
void face_node_number_error_check (const unsigned &i) const
 Range check for face node numbers. More...
 
virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt (const int &face_index) const
 
virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt (const int &face_index) const
 
- Public Member Functions inherited from oomph::GeneralisedElement
 GeneralisedElement ()
 Constructor: Initialise all pointers and all values to zero. More...
 
virtual ~GeneralisedElement ()
 Virtual destructor to clean up any memory allocated by the object. More...
 
 GeneralisedElement (const GeneralisedElement &)=delete
 Broken copy constructor. More...
 
void operator= (const GeneralisedElement &)=delete
 Broken assignment operator. More...
 
Data *& internal_data_pt (const unsigned &i)
 Return a pointer to i-th internal data object. More...
 
Data *const & internal_data_pt (const unsigned &i) const
 Return a pointer to i-th internal data object (const version) More...
 
Data *& external_data_pt (const unsigned &i)
 Return a pointer to i-th external data object. More...
 
Data *const & external_data_pt (const unsigned &i) const
 Return a pointer to i-th external data object (const version) More...
 
unsigned long eqn_number (const unsigned &ieqn_local) const
 
int local_eqn_number (const unsigned long &ieqn_global) const
 
unsigned add_external_data (Data *const &data_pt, const bool &fd=true)
 
bool external_data_fd (const unsigned &i) const
 
void exclude_external_data_fd (const unsigned &i)
 
void include_external_data_fd (const unsigned &i)
 
void flush_external_data ()
 Flush all external data. More...
 
void flush_external_data (Data *const &data_pt)
 Flush the object addressed by data_pt from the external data array. More...
 
unsigned ninternal_data () const
 Return the number of internal data objects. More...
 
unsigned nexternal_data () const
 Return the number of external data objects. More...
 
unsigned ndof () const
 Return the number of equations/dofs in the element. More...
 
void dof_vector (const unsigned &t, Vector< double > &dof)
 Return the vector of dof values at time level t. More...
 
void dof_pt_vector (Vector< double * > &dof_pt)
 Return the vector of pointers to dof values. More...
 
void set_internal_data_time_stepper (const unsigned &i, TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void assign_internal_eqn_numbers (unsigned long &global_number, Vector< double * > &Dof_pt)
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void add_internal_value_pt_to_map (std::map< unsigned, double * > &map_of_value_pt)
 
virtual void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void complete_setup_of_dependencies ()
 
virtual void get_residuals (Vector< double > &residuals)
 
virtual void get_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
virtual void get_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void get_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void get_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void get_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void get_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void get_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void get_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void get_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_norm (double &norm)
 
virtual unsigned ndof_types () const
 
virtual void get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
 
- Public Member Functions inherited from oomph::GeomObject
 GeomObject ()
 Default constructor. More...
 
 GeomObject (const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)
 
 GeomObject (const GeomObject &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const GeomObject &)=delete
 Broken assignment operator. More...
 
virtual ~GeomObject ()
 (Empty) destructor More...
 
unsigned nlagrangian () const
 Access function to # of Lagrangian coordinates. More...
 
unsigned ndim () const
 Access function to # of Eulerian coordinates. More...
 
void set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
 Set # of Lagrangian and Eulerian coordinates. More...
 
TimeStepper *& time_stepper_pt ()
 
TimeStepper * time_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 
- Public Member Functions inherited from oomph::NavierStokesElementWithDiagonalMassMatrices
 NavierStokesElementWithDiagonalMassMatrices ()
 Empty constructor. More...
 
virtual ~NavierStokesElementWithDiagonalMassMatrices ()
 Virtual destructor. More...
 
 NavierStokesElementWithDiagonalMassMatrices (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken copy constructor. More...
 
void operator= (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken assignment operator. More...
 
- Public Member Functions inherited from oomph::RefineableElement
 RefineableElement ()
 
virtual ~RefineableElement ()
 Destructor, delete the allocated storage for the hanging equations. More...
 
 RefineableElement (const RefineableElement &)=delete
 Broken copy constructor. More...
 
void operator= (const RefineableElement &)=delete
 Broken assignment operator. More...
 
Tree * tree_pt ()
 Access function: Pointer to quadtree representation of this element. More...
 
void set_tree_pt (Tree *my_tree_pt)
 Set pointer to quadtree representation of this element. More...
 
virtual unsigned required_nsons () const
 
bool refinement_is_enabled ()
 Flag to indicate suppression of any refinement. More...
 
void disable_refinement ()
 Suppress of any refinement for this element. More...
 
void enable_refinement ()
 Emnable refinement for this element. More...
 
template<class ELEMENT >
void split (Vector< ELEMENT * > &son_pt) const
 
int local_hang_eqn (Node *const &node_pt, const unsigned &i)
 
virtual void build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)=0
 
void set_refinement_level (const int &refine_level)
 Set the refinement level. More...
 
unsigned refinement_level () const
 Return the Refinement level. More...
 
void select_for_refinement ()
 Select the element for refinement. More...
 
void deselect_for_refinement ()
 Deselect the element for refinement. More...
 
void select_sons_for_unrefinement ()
 Unrefinement will be performed by merging the four sons of this element. More...
 
void deselect_sons_for_unrefinement ()
 
bool to_be_refined ()
 Has the element been selected for refinement? More...
 
bool sons_to_be_unrefined ()
 Has the element been selected for unrefinement? More...
 
virtual void rebuild_from_sons (Mesh *&mesh_pt)=0
 
virtual void unbuild ()
 
virtual void deactivate_element ()
 
virtual bool nodes_built ()
 Return true if all the nodes have been built, false if not. More...
 
long number () const
 Element number (for debugging/plotting) More...
 
void set_number (const long &mynumber)
 Set element number (for debugging/plotting) More...
 
virtual unsigned ncont_interpolated_values () const =0
 
virtual void get_interpolated_values (const Vector< double > &s, Vector< double > &values)
 
virtual void get_interpolated_values (const unsigned &t, const Vector< double > &s, Vector< double > &values)=0
 
virtual Node * interpolating_node_pt (const unsigned &n, const int &value_id)
 
virtual double local_one_d_fraction_of_interpolating_node (const unsigned &n1d, const unsigned &i, const int &value_id)
 
virtual Node * get_interpolating_node_at_local_coordinate (const Vector< double > &s, const int &value_id)
 
virtual unsigned ninterpolating_node (const int &value_id)
 
virtual unsigned ninterpolating_node_1d (const int &value_id)
 
virtual void interpolating_basis (const Vector< double > &s, Shape &psi, const int &value_id) const
 
virtual void check_integrity (double &max_error)=0
 
void identify_field_data_for_interactions (std::set< std::pair< Data *, unsigned >> &paired_field_data)
 
void assign_nodal_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual RefineableElement * root_element_pt ()
 
virtual RefineableElement * father_element_pt () const
 Return a pointer to the father element. More...
 
void get_father_at_refinement_level (unsigned &refinement_level, RefineableElement *&father_at_reflevel_pt)
 
virtual void initial_setup (Tree *const &adopted_father_pt=0, const unsigned &initial_p_order=0)
 
virtual void pre_build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt)
 Pre-build the element. More...
 
virtual void setup_hanging_nodes (Vector< std::ofstream * > &output_stream)
 
virtual void further_setup_hanging_nodes ()
 
void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
unsigned nshape_controlling_nodes ()
 
std::map< Node *, unsigned > shape_controlling_node_lookup ()
 
- Public Member Functions inherited from oomph::ElementWithZ2ErrorEstimator
 ElementWithZ2ErrorEstimator ()
 Default empty constructor. More...
 
 ElementWithZ2ErrorEstimator (const ElementWithZ2ErrorEstimator &)=delete
 Broken copy constructor. More...
 
void operator= (const ElementWithZ2ErrorEstimator &)=delete
 Broken assignment operator. More...
 
virtual unsigned ncompound_fluxes ()
 
virtual void compute_exact_Z2_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
 
virtual void get_Z2_compound_flux_indices (Vector< unsigned > &flux_index)
 
virtual unsigned nvertex_node () const =0
 Number of vertex nodes in the element. More...
 
virtual Node * vertex_node_pt (const unsigned &j) const =0
 
virtual unsigned nrecovery_order ()=0
 Order of recovery shape functions. More...
 

Static Public Member Functions

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

Protected Member Functions

virtual Node * pressure_node_pt (const unsigned &n_p)
 
virtual void unpin_elemental_pressure_dofs ()=0
 Unpin all pressure dofs in the element. More...
 
virtual void pin_elemental_redundant_nodal_pressure_dofs ()
 
- Protected Member Functions inherited from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations
virtual int p_local_eqn (const unsigned &n) const =0
 
virtual double dshape_and_dtest_eulerian_axi_nst (const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_axi_nst (const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_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_axi_nst (const Vector< double > &s, Shape &psi) const =0
 Compute the pressure shape functions at local coordinate s. More...
 
virtual void pshape_axi_nst (const Vector< double > &s, Shape &psi, Shape &test) const =0
 
virtual void get_body_force_axi_nst (const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
 Calculate the body force fct at a given time and Eulerian position. More...
 
virtual void get_body_force_gradient_axi_nst (const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
 
double get_source_fct (const double &time, const unsigned &ipt, const Vector< double > &x)
 Calculate the source fct at given time and Eulerian position. More...
 
virtual void get_source_fct_gradient (const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
 
virtual void get_viscosity_ratio_axisym_nst (const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, double &visc_ratio)
 
- Protected Member Functions inherited from oomph::FiniteElement
template<unsigned DIM>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double invert_jacobian_mapping (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &inverse_jacobian) const
 
virtual void dJ_eulerian_dnodal_coordinates (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<unsigned DIM>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
virtual void d_dshape_eulerian_dnodal_coordinates (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<unsigned DIM>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
virtual void transform_derivatives (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
void transform_derivatives_diagonal (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
virtual void transform_second_derivatives (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
void fill_in_jacobian_from_nodal_by_fd (DenseMatrix< double > &jacobian)
 
virtual void update_before_nodal_fd ()
 
virtual void reset_after_nodal_fd ()
 
virtual void update_in_nodal_fd (const unsigned &i)
 
virtual void reset_in_nodal_fd (const unsigned &i)
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Zero-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 One-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Two-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
- Protected Member Functions inherited from oomph::GeneralisedElement
unsigned add_internal_data (Data *const &data_pt, const bool &fd=true)
 
bool internal_data_fd (const unsigned &i) const
 
void exclude_internal_data_fd (const unsigned &i)
 
void include_internal_data_fd (const unsigned &i)
 
void clear_global_eqn_numbers ()
 
void add_global_eqn_numbers (std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)
 
virtual void assign_internal_and_external_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void assign_additional_local_eqn_numbers ()
 
int internal_local_eqn (const unsigned &i, const unsigned &j) const
 
int external_local_eqn (const unsigned &i, const unsigned &j)
 
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_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_axi_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 
void fill_in_generic_dresidual_contribution_axi_nst (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
 
void fill_in_contribution_to_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 

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::GeneralisedNewtonianAxisymmetricNavierStokesEquations
void(*&)(const double &time, const Vector< double > &x, Vector< double > &f) axi_nst_body_force_fct_pt ()
 Access function for the body-force pointer. More...
 
double(*&)(const double &time, const Vector< double > &x) source_fct_pt ()
 Access function for the source-function pointer. More...
 
- Static Public Attributes inherited from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations
static Vector< double > Gamma
 Vector to decide whether the stress-divergence form is used or not. More...
 
- Static Public Attributes inherited from oomph::FiniteElement
static double Tolerance_for_singular_jacobian = 1.0e-16
 Tolerance below which the jacobian is considered singular. More...
 
static bool Accept_negative_jacobian = false
 
static bool Suppress_output_while_checking_for_inverted_elements
 
- Static Public Attributes inherited from oomph::GeneralisedElement
static bool Suppress_warning_about_repeated_internal_data
 
static bool Suppress_warning_about_repeated_external_data = true
 
static double Default_fd_jacobian_step = 1.0e-8
 
- Static Protected Member Functions inherited from oomph::RefineableElement
static void check_value_id (const int &n_continuously_interpolated_values, const int &value_id)
 
- Protected Attributes inherited from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations
double * Viscosity_Ratio_pt
 
double * Density_Ratio_pt
 
double * Re_pt
 Pointer to global Reynolds number. More...
 
double * ReSt_pt
 Pointer to global Reynolds number x Strouhal number (=Womersley) More...
 
double * ReInvFr_pt
 
double * ReInvRo_pt
 
Vector< double > * G_pt
 Pointer to global gravity Vector. More...
 
void(* Body_force_fct_pt )(const double &time, const Vector< double > &x, Vector< double > &result)
 Pointer to body force function. More...
 
double(* Source_fct_pt )(const double &time, const Vector< double > &x)
 Pointer to volumetric source function. More...
 
GeneralisedNewtonianConstitutiveEquation< 3 > * Constitutive_eqn_pt
 Pointer to the generalised Newtonian constitutive equation. More...
 
bool Use_extrapolated_strainrate_to_compute_second_invariant
 
bool ALE_is_disabled
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElement * Macro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepper * Geom_object_time_stepper_pt
 
- Protected Attributes inherited from oomph::RefineableElement
Tree * Tree_pt
 A pointer to a general tree object. More...
 
unsigned Refine_level
 Refinement level. More...
 
bool To_be_refined
 Flag for refinement. More...
 
bool Refinement_is_enabled
 Flag to indicate suppression of any refinement. More...
 
bool Sons_to_be_unrefined
 Flag for unrefinement. More...
 
long Number
 Global element number – for plotting/validation purposes. More...
 
- Static Protected Attributes inherited from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations
static bool Pre_multiply_by_viscosity_ratio
 
- Static Protected Attributes inherited from oomph::FiniteElement
static const unsigned Default_Initial_Nvalue = 0
 Default value for the number of values at a node. More...
 
static const double Node_location_tolerance = 1.0e-14
 
static const unsigned N2deriv [] = {0, 1, 3, 6}
 
- Static Protected Attributes inherited from oomph::GeneralisedElement
static DenseMatrix< double > Dummy_matrix
 
static std::deque< double * > Dof_pt_deque
 
- Static Protected Attributes inherited from oomph::RefineableElement
static double Max_integrity_tolerance = 1.0e-8
 Max. allowed discrepancy in element integrity check. More...
 

Detailed Description

Refineable version of the Axisymmetric Navier–Stokes equations.

Constructor & Destructor Documentation

◆ RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations()

oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations ( )
inline

Empty Constructor.

68  : GeneralisedNewtonianAxisymmetricNavierStokesEquations(),
69  RefineableElement(),
70  ElementWithZ2ErrorEstimator()
71  {
72  }
oomph::ElementWithZ2ErrorEstimator::ElementWithZ2ErrorEstimator
ElementWithZ2ErrorEstimator()
Default empty constructor.
Definition: error_estimator.h:82
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::GeneralisedNewtonianAxisymmetricNavierStokesEquations
GeneralisedNewtonianAxisymmetricNavierStokesEquations()
Constructor: NULL the body force and source function.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:393
oomph::RefineableElement::RefineableElement
RefineableElement()
Definition: refineable_elements.h:188

Member Function Documentation

◆ dinterpolated_u_axi_nst_ddata()

void oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::dinterpolated_u_axi_nst_ddata ( const Vector< double > &  s,
const unsigned &  i,
Vector< double > &  du_ddata,
Vector< unsigned > &  global_eqn_number 
)
inlinevirtual

Compute the derivatives of the i-th component of velocity at point s with respect to all data that can affect its value. In addition, return the global equation numbers corresponding to the data. Overload the non-refineable version to take account of hanging node information

Reimplemented from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations.

180  {
181  // Find number of nodes
182  unsigned n_node = this->nnode();
183  // Local shape function
184  Shape psi(n_node);
185  // Find values of shape function at the given local coordinate
186  this->shape(s, psi);
187 
188  // Find the index at which the velocity component is stored
189  const unsigned u_nodal_index = this->u_index_axi_nst(i);
190 
191  // Storage for hang info pointer
192  HangInfo* hang_info_pt = 0;
193  // Storage for global equation
194  int global_eqn = 0;
195 
196  // Find the number of dofs associated with interpolated u
197  unsigned n_u_dof = 0;
198  for (unsigned l = 0; l < n_node; l++)
199  {
200  unsigned n_master = 1;
201 
202  // Local bool (is the node hanging)
203  bool is_node_hanging = this->node_pt(l)->is_hanging();
204 
205  // If the node is hanging, get the number of master nodes
206  if (is_node_hanging)
207  {
208  hang_info_pt = this->node_pt(l)->hanging_pt();
209  n_master = hang_info_pt->nmaster();
210  }
211  // Otherwise there is just one master node, the node itself
212  else
213  {
214  n_master = 1;
215  }
216 
217  // Loop over the master nodes
218  for (unsigned m = 0; m < n_master; m++)
219  {
220  // Get the equation number
221  if (is_node_hanging)
222  {
223  // Get the equation number from the master node
224  global_eqn =
225  hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
226  }
227  else
228  {
229  // Global equation number
230  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
231  }
232 
233  // If it's positive add to the count
234  if (global_eqn >= 0)
235  {
236  ++n_u_dof;
237  }
238  }
239  }
240 
241  // Now resize the storage schemes
242  du_ddata.resize(n_u_dof, 0.0);
243  global_eqn_number.resize(n_u_dof, 0);
244 
245  // Loop over th nodes again and set the derivatives
246  unsigned count = 0;
247  // Loop over the local nodes and sum
248  for (unsigned l = 0; l < n_node; l++)
249  {
250  unsigned n_master = 1;
251  double hang_weight = 1.0;
252 
253  // Local bool (is the node hanging)
254  bool is_node_hanging = this->node_pt(l)->is_hanging();
255 
256  // If the node is hanging, get the number of master nodes
257  if (is_node_hanging)
258  {
259  hang_info_pt = this->node_pt(l)->hanging_pt();
260  n_master = hang_info_pt->nmaster();
261  }
262  // Otherwise there is just one master node, the node itself
263  else
264  {
265  n_master = 1;
266  }
267 
268  // Loop over the master nodes
269  for (unsigned m = 0; m < n_master; m++)
270  {
271  // If the node is hanging get weight from master node
272  if (is_node_hanging)
273  {
274  // Get the hang weight from the master node
275  hang_weight = hang_info_pt->master_weight(m);
276  }
277  else
278  {
279  // Node contributes with full weight
280  hang_weight = 1.0;
281  }
282 
283  // Get the equation number
284  if (is_node_hanging)
285  {
286  // Get the equation number from the master node
287  global_eqn =
288  hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
289  }
290  else
291  {
292  // Local equation number
293  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
294  }
295 
296  if (global_eqn >= 0)
297  {
298  // Set the global equation number
299  global_eqn_number[count] = global_eqn;
300  // Set the derivative with respect to the unknown
301  du_ddata[count] = psi[l] * hang_weight;
302  // Increase the counter
303  ++count;
304  }
305  }
306  }
307  }
oomph::Data::eqn_number
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::shape
virtual void shape(const Vector< double > &s, Shape &psi) const =0
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::u_index_axi_nst
virtual unsigned u_index_axi_nst(const unsigned &i) const
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:546
oomph::Node::hanging_pt
HangInfo *const & hanging_pt() const
Definition: nodes.h:1228
oomph::Node::is_hanging
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
m
int * m
Definition: level2_cplx_impl.h:294

References oomph::Data::eqn_number(), oomph::Node::hanging_pt(), oomph::Node::is_hanging(), m, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::FiniteElement::shape(), and oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::u_index_axi_nst().

◆ fill_in_contribution_to_hessian_vector_products()

void oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_contribution_to_hessian_vector_products ( Vector< double > const &  Y,
DenseMatrix< double > const &  C,
DenseMatrix< double > &  product 
)
inlineprivatevirtual

Compute the hessian tensor vector products required to perform continuation of bifurcations analytically

Reimplemented from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations.

394  {
395  throw OomphLibError("Not yet implemented\n",
396  OOMPH_CURRENT_FUNCTION,
397  OOMPH_EXCEPTION_LOCATION);
398  }
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ fill_in_generic_dresidual_contribution_axi_nst()

void oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_generic_dresidual_contribution_axi_nst ( double *const &  parameter_pt,
Vector< double > &  dres_dparam,
DenseMatrix< double > &  djac_dparam,
DenseMatrix< double > &  dmass_matrix_dparam,
unsigned  flag 
)
inlineprivatevirtual

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

Reimplemented from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations.

382  {
383  throw OomphLibError("Not yet implemented\n",
384  OOMPH_CURRENT_FUNCTION,
385  OOMPH_EXCEPTION_LOCATION);
386  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ fill_in_generic_residual_contribution_axi_nst()

void oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_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 and mass matrix flag=2: compute all flag=1: compute both residual and Jacobian 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

Reimplemented from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations.

43  {
44  // The dimension is actually two
45  unsigned DIM = 2;
46 
47  // Find out how many nodes there are
48  unsigned n_node = nnode();
49 
50  // Get continuous time from timestepper of first node
51  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
52 
53  // Find out how many pressure dofs there are
54  unsigned n_pres = npres_axi_nst();
55 
56  // Get the local indices of the nodal coordinates
57  unsigned u_nodal_index[3];
58  for (unsigned i = 0; i < 3; ++i)
59  {
60  u_nodal_index[i] = u_index_axi_nst(i);
61  }
62 
63  // Which nodal value represents the pressure? (Negative if pressure
64  // is not based on nodal interpolation).
65  int p_index = this->p_nodal_index_axi_nst();
66 
67  // Local array of booleans that are true if the l-th pressure value is
68  // hanging (avoid repeated virtual function calls)
69  bool pressure_dof_is_hanging[n_pres];
70  // If the pressure is stored at a node
71  if (p_index >= 0)
72  {
73  // Read out whether the pressure is hanging
74  for (unsigned l = 0; l < n_pres; ++l)
75  {
76  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
77  }
78  }
79  // Otherwise the pressure is not stored at a node and so cannot hang
80  else
81  {
82  for (unsigned l = 0; l < n_pres; ++l)
83  {
84  pressure_dof_is_hanging[l] = false;
85  }
86  }
87 
88 
89  // Set up memory for the shape and test functions
90  Shape psif(n_node), testf(n_node);
91  DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
92 
93 
94  // Set up memory for pressure shape and test functions
95  Shape psip(n_pres), testp(n_pres);
96 
97  // Set the value of Nintpt
98  unsigned Nintpt = integral_pt()->nweight();
99 
100  // Set the Vector to hold local coordinates
101  Vector<double> s(DIM);
102 
103  // Get Physical Variables from Element
104  // Reynolds number must be multiplied by the density ratio
105  double scaled_re = re() * density_ratio();
106  double scaled_re_st = re_st() * density_ratio();
107  double scaled_re_inv_fr = re_invfr() * density_ratio();
108  double scaled_re_inv_ro = re_invro() * density_ratio();
109  double visc_ratio = viscosity_ratio(); // hierher -- rewrite and
110  // make consistent with
111  // non-refineable version
112  Vector<double> G = g();
113 
114  // Integers to store the local equation and unknown numbers
115  int local_eqn = 0, local_unknown = 0;
116 
117  // Local storage for pointers to hang info objects
118  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
119 
120  // Loop over the integration points
121  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
122  {
123  // Assign values of s
124  for (unsigned i = 0; i < DIM; i++)
125  {
126  s[i] = integral_pt()->knot(ipt, i);
127  }
128 
129  // Get the integral weight
130  double w = integral_pt()->weight(ipt);
131 
132  // Call the derivatives of the shape and test functions
133  double J = dshape_and_dtest_eulerian_at_knot_axi_nst(
134  ipt, psif, dpsifdx, testf, dtestfdx);
135 
136  // Call the pressure shape and test functions
137  pshape_axi_nst(s, psip, testp);
138 
139  // Premultiply the weights and the Jacobian
140  double W = w * J;
141 
142  // Calculate local values of the pressure and velocity components
143  //--------------------------------------------------------------
144  double interpolated_p = 0.0;
145  Vector<double> interpolated_u(DIM + 1, 0.0);
146  Vector<double> interpolated_x(DIM, 0.0);
147  Vector<double> mesh_veloc(DIM, 0.0);
148  Vector<double> dudt(DIM + 1, 0.0);
149  DenseMatrix<double> interpolated_dudx(DIM + 1, DIM, 0.0);
150 
151  // Calculate pressure
152  for (unsigned l = 0; l < n_pres; l++)
153  {
154  interpolated_p += p_axi_nst(l) * psip[l];
155  }
156 
157 
158  // Calculate velocities and derivatives
159 
160  // Loop over nodes
161  for (unsigned l = 0; l < n_node; l++)
162  {
163  // Cache the shape function
164  const double psif_ = psif(l);
165  // Loop over directions
166  for (unsigned i = 0; i < DIM; i++)
167  {
168  interpolated_x[i] += nodal_position(l, i) * psif_;
169  // mesh_veloc[i] +=dnodal_position_dt(l,i)*psif(l);
170  }
171 
172  for (unsigned i = 0; i < DIM + 1; i++)
173  {
174  const double u_value = nodal_value(l, u_nodal_index[i]);
175  interpolated_u[i] += u_value * psif_;
176  dudt[i] += du_dt_axi_nst(l, i) * psif_;
177  // Loop over derivative directions for gradients
178  for (unsigned j = 0; j < DIM; j++)
179  {
180  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
181  }
182  }
183  }
184 
185  // Get the mesh velocity if ALE is enabled
186  if (!ALE_is_disabled)
187  {
188  // Loop over nodes
189  for (unsigned l = 0; l < n_node; l++)
190  {
191  // Loop over the two coordinate directions
192  for (unsigned i = 0; i < 2; i++)
193  {
194  mesh_veloc[i] += this->dnodal_position_dt(l, i) * psif(l);
195  }
196  }
197  }
198 
199  // The strainrate used to compute the second invariant
200  DenseMatrix<double> strainrate_to_compute_second_invariant(3, 3, 0.0);
201 
202  // the strainrate used to calculate the second invariant
203  // can be either the current one or the one extrapolated from
204  // previous velocity values
205  if (!Use_extrapolated_strainrate_to_compute_second_invariant)
206  {
207  strain_rate(s, strainrate_to_compute_second_invariant);
208  }
209  else
210  {
211  extrapolated_strain_rate(ipt, strainrate_to_compute_second_invariant);
212  }
213 
214  // Calculate the second invariant
215  double second_invariant = SecondInvariantHelper::second_invariant(
216  strainrate_to_compute_second_invariant);
217 
218  // Get the viscosity according to the constitutive equation
219  double viscosity = Constitutive_eqn_pt->viscosity(second_invariant);
220 
221  // Get the user-defined body force terms
222  Vector<double> body_force(DIM + 1);
223  get_body_force_axi_nst(time, ipt, s, interpolated_x, body_force);
224 
225  // Get the user-defined source function
226  double source = get_source_fct(time, ipt, interpolated_x);
227 
228  // r is the first postition component
229  double r = interpolated_x[0];
230 
231  // obacht set up vectors of the viscosity differentiated w.r.t.
232  // the velocity components (radial, axial, azimuthal)
233  Vector<double> dviscosity_dUr(n_node, 0.0);
234  Vector<double> dviscosity_dUz(n_node, 0.0);
235  Vector<double> dviscosity_dUphi(n_node, 0.0);
236 
237  if (!Use_extrapolated_strainrate_to_compute_second_invariant)
238  {
239  // Calculate the derivate of the viscosity w.r.t. the second invariant
240  double dviscosity_dsecond_invariant =
241  Constitutive_eqn_pt->dviscosity_dinvariant(second_invariant);
242 
243  // calculate a reference strainrate
244  DenseMatrix<double> strainrate_ref(3, 3, 0.0);
245  strain_rate(s, strainrate_ref);
246 
247  // pre-compute the derivative of the second invariant w.r.t. the
248  // entries in the rate of strain tensor
249  DenseMatrix<double> dinvariant_dstrainrate(3, 3, 0.0);
250 
251  // d I_2 / d epsilon_{r,r}
252  dinvariant_dstrainrate(0, 0) =
253  strainrate_ref(1, 1) + strainrate_ref(2, 2);
254  // d I_2 / d epsilon_{z,z}
255  dinvariant_dstrainrate(1, 1) =
256  strainrate_ref(0, 0) + strainrate_ref(2, 2);
257  // d I_2 / d epsilon_{phi,phi}
258  dinvariant_dstrainrate(2, 2) =
259  strainrate_ref(0, 0) + strainrate_ref(1, 1);
260  // d I_2 / d epsilon_{r,z}
261  dinvariant_dstrainrate(0, 1) = -strainrate_ref(1, 0);
262  // d I_2 / d epsilon_{z,r}
263  dinvariant_dstrainrate(1, 0) = -strainrate_ref(0, 1);
264  // d I_2 / d epsilon_{r,phi}
265  dinvariant_dstrainrate(0, 2) = -strainrate_ref(2, 0);
266  // d I_2 / d epsilon_{phi,r}
267  dinvariant_dstrainrate(2, 0) = -strainrate_ref(0, 2);
268  // d I_2 / d epsilon_{phi,z}
269  dinvariant_dstrainrate(2, 1) = -strainrate_ref(1, 2);
270  // d I_2 / d epsilon_{z,phi}
271  dinvariant_dstrainrate(1, 2) = -strainrate_ref(2, 1);
272 
273  // loop over the nodes
274  for (unsigned l = 0; l < n_node; l++)
275  {
276  // Get pointer to l-th local node
277  // Node* nod_pt = node_pt(l);
278 
279  // loop over the three velocity components
280  for (unsigned i = 0; i < 3; i++)
281  {
282  // initialise the derivative of the second invariant w.r.t. the
283  // unknown velocity U_{i,l}
284  double dinvariant_dunknown = 0.0;
285 
286  // loop over the entries of the rate of strain tensor
287  for (unsigned m = 0; m < 3; m++)
288  {
289  for (unsigned n = 0; n < 3; n++)
290  {
291  // initialise the derivative of the strainrate w.r.t. the
292  // unknown velocity U_{i,l}
293  double dstrainrate_dunknown = 0.0;
294 
295  // switch based on first index
296  switch (m)
297  {
298  // epsilon_{r ...}
299  case 0:
300 
301  // switch for second index
302  switch (n)
303  {
304  // epsilon_{r r}
305  case 0:
306  if (i == 0)
307  {
308  dstrainrate_dunknown = dpsifdx(l, 0);
309  }
310  break;
311 
312  // epsilon_{r z}
313  case 1:
314  if (i == 0)
315  {
316  dstrainrate_dunknown = 0.5 * dpsifdx(l, 1);
317  }
318  else if (i == 1)
319  {
320  dstrainrate_dunknown = 0.5 * dpsifdx(l, 0);
321  }
322  break;
323 
324  // epsilon_{r phi}
325  case 2:
326  if (i == 2)
327  {
328  dstrainrate_dunknown =
329  0.5 * dpsifdx(l, 0) - 0.5 / r * psif[l];
330  }
331  break;
332 
333  default:
334  std::ostringstream error_stream;
335  error_stream << "Should never get here...";
336  throw OomphLibError(error_stream.str(),
337  OOMPH_CURRENT_FUNCTION,
338  OOMPH_EXCEPTION_LOCATION);
339  }
340 
341  break;
342 
343  // epsilon_{z ...}
344  case 1:
345 
346  // switch for second index
347  switch (n)
348  {
349  // epsilon_{z r}
350  case 0:
351  if (i == 0)
352  {
353  dstrainrate_dunknown = 0.5 * dpsifdx(l, 1);
354  }
355  else if (i == 1)
356  {
357  dstrainrate_dunknown = 0.5 * dpsifdx(l, 0);
358  }
359  break;
360 
361  // epsilon_{z z}
362  case 1:
363  if (i == 1)
364  {
365  dstrainrate_dunknown = dpsifdx(l, 1);
366  }
367  else
368  {
369  // dstrainrate_dunknown=0.0;
370  }
371  break;
372 
373  // epsilon_{z phi}
374  case 2:
375  if (i == 2)
376  {
377  dstrainrate_dunknown = 0.5 * dpsifdx(l, 1);
378  }
379  break;
380 
381  default:
382  std::ostringstream error_stream;
383  error_stream << "Should never get here...";
384  throw OomphLibError(error_stream.str(),
385  OOMPH_CURRENT_FUNCTION,
386  OOMPH_EXCEPTION_LOCATION);
387  }
388 
389  break;
390 
391  // epsilon_{phi ...}
392  case 2:
393 
394  // switch for second index
395  switch (n)
396  {
397  // epsilon_{phi r}
398  case 0:
399  if (i == 2)
400  {
401  dstrainrate_dunknown =
402  0.5 * dpsifdx(l, 0) - 0.5 / r * psif[l];
403  }
404  break;
405 
406  // epsilon_{phi z}
407  case 1:
408  if (i == 2)
409  {
410  dstrainrate_dunknown = 0.5 * dpsifdx(l, 1);
411  }
412  break;
413 
414  // epsilon_{phi phi}
415  case 2:
416  if (i == 0)
417  {
418  dstrainrate_dunknown = 1.0 / r * psif[l];
419  }
420  break;
421 
422  default:
423  std::ostringstream error_stream;
424  error_stream << "Should never get here...";
425  throw OomphLibError(error_stream.str(),
426  OOMPH_CURRENT_FUNCTION,
427  OOMPH_EXCEPTION_LOCATION);
428  }
429 
430  break;
431 
432  default:
433  std::ostringstream error_stream;
434  error_stream << "Should never get here...";
435  throw OomphLibError(error_stream.str(),
436  OOMPH_CURRENT_FUNCTION,
437  OOMPH_EXCEPTION_LOCATION);
438  }
439 
440  dinvariant_dunknown +=
441  dinvariant_dstrainrate(m, n) * dstrainrate_dunknown;
442  }
443  }
444 
445  switch (i)
446  {
447  case 0:
448  dviscosity_dUr[l] =
449  dviscosity_dsecond_invariant * dinvariant_dunknown;
450  break;
451 
452  case 1:
453  dviscosity_dUz[l] =
454  dviscosity_dsecond_invariant * dinvariant_dunknown;
455  break;
456 
457  case 2:
458  dviscosity_dUphi[l] =
459  dviscosity_dsecond_invariant * dinvariant_dunknown;
460  break;
461 
462  default:
463  std::ostringstream error_stream;
464  error_stream << "Should never get here...";
465  throw OomphLibError(error_stream.str(),
466  OOMPH_CURRENT_FUNCTION,
467  OOMPH_EXCEPTION_LOCATION);
468  }
469  }
470  }
471  }
472 
473 
474  // MOMENTUM EQUATIONS
475  //==================
476  // Number of master nodes and storage for the weight of the shape function
477  unsigned n_master = 1;
478  double hang_weight = 1.0;
479 
480  // Loop over the nodes for the test functions/equations
481  //----------------------------------------------------
482  for (unsigned l = 0; l < n_node; l++)
483  {
484  // Local boolean that indicates the hanging status of the node
485  bool is_node_hanging = node_pt(l)->is_hanging();
486 
487  // If the node is hanging
488  if (is_node_hanging)
489  {
490  // Get the hanging pointer
491  hang_info_pt = node_pt(l)->hanging_pt();
492  // Read out number of master nodes from hanging data
493  n_master = hang_info_pt->nmaster();
494  }
495  // Otherwise the node is its own master
496  else
497  {
498  n_master = 1;
499  }
500 
501  // Loop over the master nodes
502  for (unsigned m = 0; m < n_master; m++)
503  {
504  // Loop over velocity components for equations
505  for (unsigned i = 0; i < DIM + 1; i++)
506  {
507  // Get the equation number
508  // If the node is hanging
509  if (is_node_hanging)
510  {
511  // Get the equation number from the master node
512  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
513  u_nodal_index[i]);
514  // Get the hang weight from the master node
515  hang_weight = hang_info_pt->master_weight(m);
516  }
517  // If the node is not hanging
518  else
519  {
520  // Local equation number
521  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
522  // Node contributes with full weight
523  hang_weight = 1.0;
524  }
525 
526  // If it's not a boundary condition...
527  if (local_eqn >= 0)
528  {
529  // initialise for residual calculation
530  double sum = 0.0;
531 
532  switch (i)
533  {
534  // RADIAL MOMENTUM EQUATION
535  case 0:
536  // Add the user-defined body force terms
537  sum += r * body_force[0] * testf[l] * W * hang_weight;
538 
539  // Add the gravitational body force term
540  sum +=
541  r * scaled_re_inv_fr * testf[l] * G[0] * W * hang_weight;
542 
543  // Add the pressure gradient term
544  sum += interpolated_p * (testf[l] + r * dtestfdx(l, 0)) * W *
545  hang_weight;
546 
547  // Add in the stress tensor terms
548  // The viscosity ratio needs to go in here to ensure
549  // continuity of normal stress is satisfied even in flows
550  // with zero pressure gradient!
551  // stress term 1
552  sum -= visc_ratio * viscosity * r * (1.0 + Gamma[0]) *
553  interpolated_dudx(0, 0) * dtestfdx(l, 0) * W *
554  hang_weight;
555 
556  // stress term 2
557  sum -= visc_ratio * viscosity * r *
558  (interpolated_dudx(0, 1) +
559  Gamma[0] * interpolated_dudx(1, 0)) *
560  dtestfdx(l, 1) * W * hang_weight;
561 
562  // stress term 3
563  sum -= visc_ratio * viscosity * (1.0 + Gamma[0]) *
564  interpolated_u[0] * testf[l] * W * hang_weight / r;
565 
566  // Add in the inertial terms
567  // du/dt term
568  sum -=
569  scaled_re_st * r * dudt[0] * testf[l] * W * hang_weight;
570 
571  // Convective terms
572  sum -= scaled_re *
573  (r * interpolated_u[0] * interpolated_dudx(0, 0) -
574  interpolated_u[2] * interpolated_u[2] +
575  r * interpolated_u[1] * interpolated_dudx(0, 1)) *
576  testf[l] * W * hang_weight;
577 
578  // Mesh velocity terms
579  if (!ALE_is_disabled)
580  {
581  for (unsigned k = 0; k < 2; k++)
582  {
583  sum += scaled_re_st * r * mesh_veloc[k] *
584  interpolated_dudx(0, k) * testf[l] * W *
585  hang_weight;
586  }
587  }
588 
589  // Coriolis term
590  sum += 2.0 * r * scaled_re_inv_ro * interpolated_u[2] *
591  testf[l] * W * hang_weight;
592 
593  break;
594 
595  // AXIAL MOMENTUM EQUATION
596  case 1:
597  // If it's not a boundary condition
598  // Add the user-defined body force terms
599  // Remember to multiply by the density ratio!
600  sum += r * body_force[1] * testf[l] * W * hang_weight;
601 
602  // Add the gravitational body force term
603  sum +=
604  r * scaled_re_inv_fr * testf[l] * G[1] * W * hang_weight;
605 
606  // Add the pressure gradient term
607  sum += r * interpolated_p * dtestfdx(l, 1) * W * hang_weight;
608 
609  // Add in the stress tensor terms
610  // The viscosity ratio needs to go in here to ensure
611  // continuity of normal stress is satisfied even in flows
612  // with zero pressure gradient!
613  // stress term 1
614  sum -= visc_ratio * viscosity * r *
615  (interpolated_dudx(1, 0) +
616  Gamma[1] * interpolated_dudx(0, 1)) *
617  dtestfdx(l, 0) * W * hang_weight;
618 
619  // stress term 2
620  sum -= visc_ratio * viscosity * r * (1.0 + Gamma[1]) *
621  interpolated_dudx(1, 1) * dtestfdx(l, 1) * W *
622  hang_weight;
623 
624  // Add in the inertial terms
625  // du/dt term
626  sum -=
627  scaled_re_st * r * dudt[1] * testf[l] * W * hang_weight;
628 
629  // Convective terms
630  sum -= scaled_re *
631  (r * interpolated_u[0] * interpolated_dudx(1, 0) +
632  r * interpolated_u[1] * interpolated_dudx(1, 1)) *
633  testf[l] * W * hang_weight;
634 
635  // Mesh velocity terms
636  if (!ALE_is_disabled)
637  {
638  for (unsigned k = 0; k < 2; k++)
639  {
640  sum += scaled_re_st * r * mesh_veloc[k] *
641  interpolated_dudx(1, k) * testf[l] * W *
642  hang_weight;
643  }
644  }
645  break;
646 
647  // AZIMUTHAL MOMENTUM EQUATION
648  case 2:
649  // Add the user-defined body force terms
650  // Remember to multiply by the density ratio!
651  sum += r * body_force[2] * testf[l] * W * hang_weight;
652 
653  // Add the gravitational body force term
654  sum +=
655  r * scaled_re_inv_fr * testf[l] * G[2] * W * hang_weight;
656 
657  // There is NO pressure gradient term
658 
659  // Add in the stress tensor terms
660  // The viscosity ratio needs to go in here to ensure
661  // continuity of normal stress is satisfied even in flows
662  // with zero pressure gradient!
663  // stress term 1
664  sum -= visc_ratio * viscosity *
665  (r * interpolated_dudx(2, 0) -
666  Gamma[0] * interpolated_u[2]) *
667  dtestfdx(l, 0) * W * hang_weight;
668 
669  // stress term 2
670  sum -= visc_ratio * viscosity * r * interpolated_dudx(2, 1) *
671  dtestfdx(l, 1) * W * hang_weight;
672 
673  // stress term 3
674  sum -= visc_ratio * viscosity *
675  ((interpolated_u[2] / r) -
676  Gamma[0] * interpolated_dudx(2, 0)) *
677  testf[l] * W * hang_weight;
678 
679 
680  // Add in the inertial terms
681  // du/dt term
682  sum -=
683  scaled_re_st * r * dudt[2] * testf[l] * W * hang_weight;
684 
685  // Convective terms
686  sum -= scaled_re *
687  (r * interpolated_u[0] * interpolated_dudx(2, 0) +
688  interpolated_u[0] * interpolated_u[2] +
689  r * interpolated_u[1] * interpolated_dudx(2, 1)) *
690  testf[l] * W * hang_weight;
691 
692  // Mesh velocity terms
693  if (!ALE_is_disabled)
694  {
695  for (unsigned k = 0; k < 2; k++)
696  {
697  sum += scaled_re_st * r * mesh_veloc[k] *
698  interpolated_dudx(2, k) * testf[l] * W *
699  hang_weight;
700  }
701  }
702 
703  // Coriolis term
704  sum -= 2.0 * r * scaled_re_inv_ro * interpolated_u[0] *
705  testf[l] * W * hang_weight;
706 
707  break;
708  }
709 
710  // Add contribution
711  residuals[local_eqn] += sum;
712 
713  // CALCULATE THE JACOBIAN
714  if (flag)
715  {
716  // Number of master nodes and weights
717  unsigned n_master2 = 1;
718  double hang_weight2 = 1.0;
719  // Loop over the velocity nodes for columns
720  for (unsigned l2 = 0; l2 < n_node; l2++)
721  {
722  // Local boolean for hanging status
723  bool is_node2_hanging = node_pt(l2)->is_hanging();
724 
725  // If the node is hanging
726  if (is_node2_hanging)
727  {
728  hang_info2_pt = node_pt(l2)->hanging_pt();
729  // Read out number of master nodes from hanging data
730  n_master2 = hang_info2_pt->nmaster();
731  }
732  // Otherwise the node is its own master
733  else
734  {
735  n_master2 = 1;
736  }
737 
738  // Loop over the master nodes
739  for (unsigned m2 = 0; m2 < n_master2; m2++)
740  {
741  // Loop over the velocity components
742  for (unsigned i2 = 0; i2 < DIM + 1; i2++)
743  {
744  // Get the number of the unknown
745  // If the node is hanging
746  if (is_node2_hanging)
747  {
748  // Get the equation number from the master node
749  local_unknown = this->local_hang_eqn(
750  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
751  hang_weight2 = hang_info2_pt->master_weight(m2);
752  }
753  else
754  {
755  local_unknown =
756  this->nodal_local_eqn(l2, u_nodal_index[i2]);
757  hang_weight2 = 1.0;
758  }
759 
760  // If the unknown is non-zero
761  if (local_unknown >= 0)
762  {
763  // Different results for i and i2
764  switch (i)
765  {
766  // RADIAL MOMENTUM EQUATION
767  case 0:
768  switch (i2)
769  {
770  // radial component
771  case 0:
772 
773  // Add the mass matrix entries
774  if (flag == 2)
775  {
776  mass_matrix(local_eqn, local_unknown) +=
777  scaled_re_st * r * psif[l2] * testf[l] * W *
778  hang_weight * hang_weight2;
779  }
780 
781  // Add contribution to the Jacobian matrix
782  // stress term 1
783  jacobian(local_eqn, local_unknown) -=
784  visc_ratio * viscosity * r *
785  (1.0 + Gamma[0]) * dpsifdx(l2, 0) *
786  dtestfdx(l, 0) * W * hang_weight *
787  hang_weight2;
788 
789  if (
790  !Use_extrapolated_strainrate_to_compute_second_invariant)
791  {
792  // stress term 1 contribution from
793  // generalised Newtonian model
794  jacobian(local_eqn, local_unknown) -=
795  visc_ratio * dviscosity_dUr[l2] * r *
796  (1.0 + Gamma[0]) * interpolated_dudx(0, 0) *
797  dtestfdx(l, 0) * W * hang_weight *
798  hang_weight2;
799  }
800 
801  // stress term 2
802  jacobian(local_eqn, local_unknown) -=
803  visc_ratio * viscosity * r * dpsifdx(l2, 1) *
804  dtestfdx(l, 1) * W * hang_weight *
805  hang_weight2;
806 
807  if (
808  !Use_extrapolated_strainrate_to_compute_second_invariant)
809  {
810  // stress term 2 contribution from
811  // generalised Newtonian model
812  jacobian(local_eqn, local_unknown) -=
813  visc_ratio * dviscosity_dUr[l2] * r *
814  (interpolated_dudx(0, 1) +
815  Gamma[0] * interpolated_dudx(1, 0)) *
816  dtestfdx(l, 1) * W * hang_weight *
817  hang_weight2;
818  }
819 
820  // stress term 3
821  jacobian(local_eqn, local_unknown) -=
822  visc_ratio * viscosity * (1.0 + Gamma[0]) *
823  psif[l2] * testf[l] * W * hang_weight *
824  hang_weight2 / r;
825 
826  if (
827  !Use_extrapolated_strainrate_to_compute_second_invariant)
828  {
829  // stress term 3 contribution from
830  // generalised Newtonian model
831  jacobian(local_eqn, local_unknown) -=
832  visc_ratio * dviscosity_dUr[l2] *
833  (1.0 + Gamma[0]) * interpolated_u[0] *
834  testf[l] * W * hang_weight * hang_weight2 /
835  r;
836  }
837 
838  // Add in the inertial terms
839  // du/dt term
840  jacobian(local_eqn, local_unknown) -=
841  scaled_re_st * r *
842  node_pt(l2)->time_stepper_pt()->weight(1, 0) *
843  psif[l2] * testf[l] * W * hang_weight *
844  hang_weight2;
845 
846  // Convective terms
847  jacobian(local_eqn, local_unknown) -=
848  scaled_re *
849  (r * psif[l2] * interpolated_dudx(0, 0) +
850  r * interpolated_u[0] * dpsifdx(l2, 0) +
851  r * interpolated_u[1] * dpsifdx(l2, 1)) *
852  testf[l] * W * hang_weight * hang_weight2;
853 
854  // Mesh velocity terms
855  if (!ALE_is_disabled)
856  {
857  for (unsigned k = 0; k < 2; k++)
858  {
859  jacobian(local_eqn, local_unknown) +=
860  scaled_re_st * r * mesh_veloc[k] *
861  dpsifdx(l2, k) * testf[l] * W *
862  hang_weight * hang_weight2;
863  }
864  }
865  break;
866 
867  // axial component
868  case 1:
869 
870  // no stress term 1
871 
872  if (
873  !Use_extrapolated_strainrate_to_compute_second_invariant)
874  {
875  // stress term 1 contribution from
876  // generalised Newtonian model
877  jacobian(local_eqn, local_unknown) -=
878  visc_ratio * dviscosity_dUz[l2] * r *
879  (1.0 + Gamma[0]) * interpolated_dudx(0, 0) *
880  dtestfdx(l, 0) * W * hang_weight *
881  hang_weight2;
882  }
883 
884  // stress term 2
885  jacobian(local_eqn, local_unknown) -=
886  visc_ratio * viscosity * r * Gamma[0] *
887  dpsifdx(l2, 0) * dtestfdx(l, 1) * W *
888  hang_weight * hang_weight2;
889 
890  if (
891  !Use_extrapolated_strainrate_to_compute_second_invariant)
892  {
893  // stress term 2 contribution from
894  // generalised Newtonian model
895  jacobian(local_eqn, local_unknown) -=
896  visc_ratio * dviscosity_dUz[l2] * r *
897  (interpolated_dudx(0, 1) +
898  Gamma[0] * interpolated_dudx(1, 0)) *
899  dtestfdx(l, 1) * W * hang_weight *
900  hang_weight2;
901  }
902 
903  // no stress term 3
904 
905  if (
906  !Use_extrapolated_strainrate_to_compute_second_invariant)
907  {
908  // stress term 3 contribution from
909  // generalised Newtonian model
910  jacobian(local_eqn, local_unknown) -=
911  visc_ratio * dviscosity_dUz[l2] *
912  (1.0 + Gamma[0]) * interpolated_u[0] *
913  testf[l] * W * hang_weight * hang_weight2 /
914  r;
915  }
916 
917  // Convective terms
918  jacobian(local_eqn, local_unknown) -=
919  scaled_re * r * psif[l2] *
920  interpolated_dudx(0, 1) * testf[l] * W *
921  hang_weight * hang_weight2;
922  break;
923 
924  // azimuthal component
925  case 2:
926 
927  // no stress term 1
928 
929  if (
930  !Use_extrapolated_strainrate_to_compute_second_invariant)
931  {
932  // stress term 1 contribution from
933  // generalised Newtonian model
934  jacobian(local_eqn, local_unknown) -=
935  visc_ratio * dviscosity_dUphi[l2] * r *
936  (1.0 + Gamma[0]) * interpolated_dudx(0, 0) *
937  dtestfdx(l, 0) * W * hang_weight *
938  hang_weight2;
939  }
940 
941  // no stress term 2
942 
943  if (
944  !Use_extrapolated_strainrate_to_compute_second_invariant)
945  {
946  // stress term 2 contribution from
947  // generalised Newtonian model
948  jacobian(local_eqn, local_unknown) -=
949  visc_ratio * dviscosity_dUphi[l2] * r *
950  (interpolated_dudx(0, 1) +
951  Gamma[0] * interpolated_dudx(1, 0)) *
952  dtestfdx(l, 1) * W * hang_weight *
953  hang_weight2;
954  }
955 
956  // no stress term 3
957 
958  if (
959  !Use_extrapolated_strainrate_to_compute_second_invariant)
960  {
961  // stress term 3 contribution from
962  // generalised Newtonian model
963  jacobian(local_eqn, local_unknown) -=
964  visc_ratio * dviscosity_dUphi[l2] *
965  (1.0 + Gamma[0]) * interpolated_u[0] *
966  testf[l] * W * hang_weight * hang_weight2 /
967  r;
968  }
969 
970  // Convective terms
971  jacobian(local_eqn, local_unknown) -=
972  -scaled_re * 2.0 * interpolated_u[2] *
973  psif[l2] * testf[l] * W * hang_weight *
974  hang_weight2;
975 
976  // Coriolis terms
977  jacobian(local_eqn, local_unknown) +=
978  2.0 * r * scaled_re_inv_ro * psif[l2] *
979  testf[l] * W * hang_weight * hang_weight2;
980 
981  break;
982  } /*End of contribution radial momentum eqn*/
983  break;
984 
985  // AXIAL MOMENTUM EQUATION
986  case 1:
987  switch (i2)
988  {
989  // radial component
990  case 0:
991  // Add in the stress tensor terms
992  // The viscosity ratio needs to go in here to
993  // ensure continuity of normal stress is
994  // satisfied even in flows with zero pressure
995  // gradient!
996  // stress term 1
997  jacobian(local_eqn, local_unknown) -=
998  visc_ratio * viscosity * r * Gamma[1] *
999  dpsifdx(l2, 1) * dtestfdx(l, 0) * W *
1000  hang_weight * hang_weight2;
1001 
1002  if (
1003  !Use_extrapolated_strainrate_to_compute_second_invariant)
1004  {
1005  // stress term 1 contribution from
1006  // generalised Newtonian model
1007  jacobian(local_eqn, local_unknown) -=
1008  visc_ratio * dviscosity_dUr[l2] * r *
1009  (interpolated_dudx(1, 0) +
1010  Gamma[1] * interpolated_dudx(0, 1)) *
1011  dtestfdx(l, 0) * W * hang_weight *
1012  hang_weight2;
1013  }
1014 
1015  // no stress term 2
1016 
1017  if (
1018  !Use_extrapolated_strainrate_to_compute_second_invariant)
1019  {
1020  // stress term 2 contribution from
1021  // generalised Newtonian model
1022  jacobian(local_eqn, local_unknown) -=
1023  visc_ratio * dviscosity_dUr[l2] * r *
1024  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
1025  dtestfdx(l, 1) * W * hang_weight *
1026  hang_weight2;
1027  }
1028 
1029  // Convective terms
1030  jacobian(local_eqn, local_unknown) -=
1031  scaled_re * r * psif[l2] *
1032  interpolated_dudx(1, 0) * testf[l] * W *
1033  hang_weight * hang_weight2;
1034  break;
1035 
1036  // axial component
1037  case 1:
1038 
1039  // Add the mass matrix terms
1040  if (flag == 2)
1041  {
1042  mass_matrix(local_eqn, local_unknown) +=
1043  scaled_re_st * r * psif[l2] * testf[l] * W *
1044  hang_weight * hang_weight2;
1045  }
1046 
1047  // Add in the stress tensor terms
1048  // The viscosity ratio needs to go in here to
1049  // ensure continuity of normal stress is
1050  // satisfied even in flows with zero pressure
1051  // gradient!
1052  // stress term 1
1053  jacobian(local_eqn, local_unknown) -=
1054  visc_ratio * viscosity * r * dpsifdx(l2, 0) *
1055  dtestfdx(l, 0) * W * hang_weight *
1056  hang_weight2;
1057 
1058  if (
1059  !Use_extrapolated_strainrate_to_compute_second_invariant)
1060  {
1061  // stress term 1 contribution from
1062  // generalised Newtonian model
1063  jacobian(local_eqn, local_unknown) -=
1064  visc_ratio * dviscosity_dUz[l2] * r *
1065  (interpolated_dudx(1, 0) +
1066  Gamma[1] * interpolated_dudx(0, 1)) *
1067  dtestfdx(l, 0) * W * hang_weight *
1068  hang_weight2;
1069  }
1070 
1071  // stress term 2
1072  jacobian(local_eqn, local_unknown) -=
1073  visc_ratio * viscosity * r *
1074  (1.0 + Gamma[1]) * dpsifdx(l2, 1) *
1075  dtestfdx(l, 1) * W * hang_weight *
1076  hang_weight2;
1077 
1078  if (
1079  !Use_extrapolated_strainrate_to_compute_second_invariant)
1080  {
1081  // stress term 2 contribution from
1082  // generalised Newtonian model
1083  jacobian(local_eqn, local_unknown) -=
1084  visc_ratio * dviscosity_dUz[l2] * r *
1085  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
1086  dtestfdx(l, 1) * W * hang_weight *
1087  hang_weight2;
1088  }
1089 
1090  // Add in the inertial terms
1091  // du/dt term
1092  jacobian(local_eqn, local_unknown) -=
1093  scaled_re_st * r *
1094  node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1095  psif[l2] * testf[l] * W * hang_weight *
1096  hang_weight2;
1097 
1098  // Convective terms
1099  jacobian(local_eqn, local_unknown) -=
1100  scaled_re *
1101  (r * interpolated_u[0] * dpsifdx(l2, 0) +
1102  r * psif[l2] * interpolated_dudx(1, 1) +
1103  r * interpolated_u[1] * dpsifdx(l2, 1)) *
1104  testf[l] * W * hang_weight * hang_weight2;
1105 
1106  // Mesh velocity terms
1107  if (!ALE_is_disabled)
1108  {
1109  for (unsigned k = 0; k < 2; k++)
1110  {
1111  jacobian(local_eqn, local_unknown) +=
1112  scaled_re_st * r * mesh_veloc[k] *
1113  dpsifdx(l2, k) * testf[l] * W *
1114  hang_weight * hang_weight2;
1115  }
1116  }
1117  break;
1118 
1119  // azimuthal component
1120  case 2:
1121  // There are no azimithal terms in the axial
1122  // momentum equation
1123  // edit: for the generalised Newtonian elements
1124  // there are...
1125 
1126  if (
1127  !Use_extrapolated_strainrate_to_compute_second_invariant)
1128  {
1129  // stress term 1 contribution from
1130  // generalised Newtonian model
1131  jacobian(local_eqn, local_unknown) -=
1132  visc_ratio * dviscosity_dUphi[l2] * r *
1133  (interpolated_dudx(1, 0) +
1134  Gamma[1] * interpolated_dudx(0, 1)) *
1135  dtestfdx(l, 0) * W * hang_weight *
1136  hang_weight2;
1137 
1138  // stress term 2 contribution from
1139  // generalised Newtonian model
1140  jacobian(local_eqn, local_unknown) -=
1141  visc_ratio * dviscosity_dUphi[l2] * r *
1142  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
1143  dtestfdx(l, 1) * W * hang_weight *
1144  hang_weight2;
1145  }
1146  break;
1147  }
1148  break;
1149 
1150  // AZIMUTHAL MOMENTUM EQUATION
1151  case 2:
1152  switch (i2)
1153  {
1154  // radial component
1155  case 0:
1156  if (
1157  !Use_extrapolated_strainrate_to_compute_second_invariant)
1158  {
1159  // stress term 1 contribution from
1160  // generalised Newtonian model
1161  jacobian(local_eqn, local_unknown) -=
1162  visc_ratio * dviscosity_dUr[l2] *
1163  (r * interpolated_dudx(2, 0) -
1164  Gamma[0] * interpolated_u[2]) *
1165  dtestfdx(l, 0) * W * hang_weight *
1166  hang_weight2;
1167 
1168  // stress term 2 contribution from
1169  // generalised Newtonian model
1170  jacobian(local_eqn, local_unknown) -=
1171  visc_ratio * dviscosity_dUr[l2] * r *
1172  interpolated_dudx(2, 1) * dtestfdx(l, 1) *
1173  W * hang_weight * hang_weight2;
1174 
1175  // stress term 3 contribution from
1176  // generalised Newtonian model
1177  jacobian(local_eqn, local_unknown) -=
1178  visc_ratio * dviscosity_dUr[l2] *
1179  ((interpolated_u[2] / r) -
1180  Gamma[0] * interpolated_dudx(2, 0)) *
1181  testf[l] * W * hang_weight * hang_weight2;
1182  }
1183 
1184  // Convective terms
1185  jacobian(local_eqn, local_unknown) -=
1186  scaled_re *
1187  (r * psif[l2] * interpolated_dudx(2, 0) +
1188  psif[l2] * interpolated_u[2]) *
1189  testf[l] * W * hang_weight * hang_weight2;
1190 
1191  // Coriolis term
1192  jacobian(local_eqn, local_unknown) -=
1193  2.0 * r * scaled_re_inv_ro * psif[l2] *
1194  testf[l] * W * hang_weight * hang_weight2;
1195 
1196  break;
1197 
1198  // axial component
1199  case 1:
1200  if (
1201  !Use_extrapolated_strainrate_to_compute_second_invariant)
1202  {
1203  // stress term 1 contribution from
1204  // generalised Newtonian model
1205  jacobian(local_eqn, local_unknown) -=
1206  visc_ratio * dviscosity_dUz[l2] *
1207  (r * interpolated_dudx(2, 0) -
1208  Gamma[0] * interpolated_u[2]) *
1209  dtestfdx(l, 0) * W * hang_weight *
1210  hang_weight2;
1211 
1212  // stress term 2 contribution from
1213  // generalised Newtonian model
1214  jacobian(local_eqn, local_unknown) -=
1215  visc_ratio * dviscosity_dUz[l2] * r *
1216  interpolated_dudx(2, 1) * dtestfdx(l, 1) *
1217  W * hang_weight * hang_weight2;
1218 
1219  // stress term 3 contribution from
1220  // generalised Newtonian model
1221  jacobian(local_eqn, local_unknown) -=
1222  visc_ratio * dviscosity_dUz[l2] *
1223  ((interpolated_u[2] / r) -
1224  Gamma[0] * interpolated_dudx(2, 0)) *
1225  testf[l] * W * hang_weight * hang_weight2;
1226  }
1227 
1228  // Convective terms
1229  jacobian(local_eqn, local_unknown) -=
1230  scaled_re * r * psif[l2] *
1231  interpolated_dudx(2, 1) * testf[l] * W *
1232  hang_weight * hang_weight2;
1233  break;
1234 
1235  // azimuthal component
1236  case 2:
1237 
1238  // Add the mass matrix terms
1239  if (flag == 2)
1240  {
1241  mass_matrix(local_eqn, local_unknown) +=
1242  scaled_re_st * r * psif[l2] * testf[l] * W *
1243  hang_weight * hang_weight2;
1244  }
1245 
1246  // Add in the stress tensor terms
1247  // The viscosity ratio needs to go in here to
1248  // ensure continuity of normal stress is
1249  // satisfied even in flows with zero pressure
1250  // gradient!
1251  // stress term 1
1252  jacobian(local_eqn, local_unknown) -=
1253  visc_ratio * viscosity *
1254  (r * dpsifdx(l2, 0) - Gamma[0] * psif[l2]) *
1255  dtestfdx(l, 0) * W * hang_weight *
1256  hang_weight2;
1257 
1258  if (
1259  !Use_extrapolated_strainrate_to_compute_second_invariant)
1260  {
1261  // stress term 1 contribution from
1262  // generalised Newtonian model
1263  jacobian(local_eqn, local_unknown) -=
1264  visc_ratio * dviscosity_dUphi[l2] *
1265  (r * interpolated_dudx(2, 0) -
1266  Gamma[0] * interpolated_u[2]) *
1267  dtestfdx(l, 0) * W * hang_weight *
1268  hang_weight2;
1269  }
1270 
1271  // stress term 2
1272  jacobian(local_eqn, local_unknown) -=
1273  visc_ratio * viscosity * r * dpsifdx(l2, 1) *
1274  dtestfdx(l, 1) * W * hang_weight *
1275  hang_weight2;
1276 
1277  if (
1278  !Use_extrapolated_strainrate_to_compute_second_invariant)
1279  {
1280  // stress term 2 contribution from
1281  // generalised Newtonian model
1282  jacobian(local_eqn, local_unknown) -=
1283  visc_ratio * dviscosity_dUphi[l2] * r *
1284  interpolated_dudx(2, 1) * dtestfdx(l, 1) *
1285  W * hang_weight * hang_weight2;
1286  }
1287 
1288  // stress term 3
1289  jacobian(local_eqn, local_unknown) -=
1290  visc_ratio * viscosity *
1291  ((psif[l2] / r) - Gamma[0] * dpsifdx(l2, 0)) *
1292  testf[l] * W * hang_weight * hang_weight2;
1293 
1294  if (
1295  !Use_extrapolated_strainrate_to_compute_second_invariant)
1296  {
1297  // stress term 3 contribution from
1298  // generalised Newtonian model
1299  jacobian(local_eqn, local_unknown) -=
1300  visc_ratio * dviscosity_dUphi[l2] *
1301  ((interpolated_u[2] / r) -
1302  Gamma[0] * interpolated_dudx(2, 0)) *
1303  testf[l] * W * hang_weight * hang_weight2;
1304  }
1305 
1306  // Add in the inertial terms
1307  // du/dt term
1308  jacobian(local_eqn, local_unknown) -=
1309  scaled_re_st * r *
1310  node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1311  psif[l2] * testf[l] * W * hang_weight *
1312  hang_weight2;
1313 
1314  // Convective terms
1315  jacobian(local_eqn, local_unknown) -=
1316  scaled_re *
1317  (r * interpolated_u[0] * dpsifdx(l2, 0) +
1318  interpolated_u[0] * psif[l2] +
1319  r * interpolated_u[1] * dpsifdx(l2, 1)) *
1320  testf[l] * W * hang_weight * hang_weight2;
1321 
1322  // Mesh velocity terms
1323  if (!ALE_is_disabled)
1324  {
1325  for (unsigned k = 0; k < 2; k++)
1326  {
1327  jacobian(local_eqn, local_unknown) +=
1328  scaled_re_st * r * mesh_veloc[k] *
1329  dpsifdx(l2, k) * testf[l] * W *
1330  hang_weight * hang_weight2;
1331  }
1332  }
1333  break;
1334  }
1335  break;
1336  }
1337  }
1338  }
1339  }
1340  } // End of loop over the nodes
1341 
1342 
1343  // Loop over the pressure shape functions
1344  for (unsigned l2 = 0; l2 < n_pres; l2++)
1345  {
1346  // If the pressure dof is hanging
1347  if (pressure_dof_is_hanging[l2])
1348  {
1349  // Pressure dof is hanging so it must be nodal-based
1350  hang_info2_pt = pressure_node_pt(l2)->hanging_pt(p_index);
1351 
1352  // Get the number of master nodes from the pressure node
1353  n_master2 = hang_info2_pt->nmaster();
1354  }
1355  // Otherwise the node is its own master
1356  else
1357  {
1358  n_master2 = 1;
1359  }
1360 
1361  // Loop over the master nodes
1362  for (unsigned m2 = 0; m2 < n_master2; m2++)
1363  {
1364  // Get the number of the unknown
1365  // If the pressure dof is hanging
1366  if (pressure_dof_is_hanging[l2])
1367  {
1368  // Get the unknown from the node
1369  local_unknown = local_hang_eqn(
1370  hang_info2_pt->master_node_pt(m2), p_index);
1371  // Get the weight from the hanging object
1372  hang_weight2 = hang_info2_pt->master_weight(m2);
1373  }
1374  else
1375  {
1376  local_unknown = p_local_eqn(l2);
1377  hang_weight2 = 1.0;
1378  }
1379 
1380  // If the unknown is not pinned
1381  if (local_unknown >= 0)
1382  {
1383  // Add in contributions to different equations
1384  switch (i)
1385  {
1386  // RADIAL MOMENTUM EQUATION
1387  case 0:
1388  jacobian(local_eqn, local_unknown) +=
1389  psip[l2] * (testf[l] + r * dtestfdx(l, 0)) * W *
1390  hang_weight * hang_weight2;
1391  break;
1392 
1393  // AXIAL MOMENTUM EQUATION
1394  case 1:
1395  jacobian(local_eqn, local_unknown) +=
1396  r * psip[l2] * dtestfdx(l, 1) * W * hang_weight *
1397  hang_weight2;
1398  break;
1399 
1400  // AZIMUTHAL MOMENTUM EQUATION
1401  case 2:
1402  break;
1403  }
1404  }
1405  }
1406  } // End of loop over pressure dofs
1407  } // End of Jacobian calculation
1408  } // End of if not boundary condition statement
1409  } // End of loop over velocity components
1410  } // End of loop over master nodes
1411  } // End of loop over nodes
1412 
1413 
1414  // CONTINUITY EQUATION
1415  //===================
1416 
1417  // Loop over the pressure shape functions
1418  for (unsigned l = 0; l < n_pres; l++)
1419  {
1420  // If the pressure dof is hanging
1421  if (pressure_dof_is_hanging[l])
1422  {
1423  // Pressure dof is hanging so it must be nodal-based
1424  hang_info_pt = pressure_node_pt(l)->hanging_pt(p_index);
1425  // Get the number of master nodes from the pressure node
1426  n_master = hang_info_pt->nmaster();
1427  }
1428  // Otherwise the node is its own master
1429  else
1430  {
1431  n_master = 1;
1432  }
1433 
1434  // Loop over the master nodes
1435  for (unsigned m = 0; m < n_master; m++)
1436  {
1437  // Get the number of the unknown
1438  // If the pressure dof is hanging
1439  if (pressure_dof_is_hanging[l])
1440  {
1441  local_eqn =
1442  local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1443  hang_weight = hang_info_pt->master_weight(m);
1444  }
1445  else
1446  {
1447  local_eqn = p_local_eqn(l);
1448  hang_weight = 1.0;
1449  }
1450 
1451  // If the equation is not pinned
1452  if (local_eqn >= 0)
1453  {
1454  // Source term
1455  residuals[local_eqn] -= r * source * testp[l] * W * hang_weight;
1456 
1457  // Gradient terms
1458  residuals[local_eqn] +=
1459  (interpolated_u[0] + r * interpolated_dudx(0, 0) +
1460  r * interpolated_dudx(1, 1)) *
1461  testp[l] * W * hang_weight;
1462 
1463  // CALCULATE THE JACOBIAN
1464  //======================
1465  if (flag)
1466  {
1467  unsigned n_master2 = 1;
1468  double hang_weight2 = 1.0;
1469  // Loop over the velocity nodes for columns
1470  for (unsigned l2 = 0; l2 < n_node; l2++)
1471  {
1472  // Local boolean to indicate whether the node is hanging
1473  bool is_node2_hanging = node_pt(l2)->is_hanging();
1474 
1475  // If the node is hanging
1476  if (is_node2_hanging)
1477  {
1478  hang_info2_pt = node_pt(l2)->hanging_pt();
1479  // Read out number of master nodes from hanging data
1480  n_master2 = hang_info2_pt->nmaster();
1481  }
1482  // Otherwise the node is its own master
1483  else
1484  {
1485  n_master2 = 1;
1486  }
1487 
1488  // Loop over the master nodes
1489  for (unsigned m2 = 0; m2 < n_master2; m2++)
1490  {
1491  // Loop over the velocity components
1492  for (unsigned i2 = 0; i2 < DIM + 1; i2++)
1493  {
1494  // Get the number of the unknown
1495  // If the node is hanging
1496  if (is_node2_hanging)
1497  {
1498  // Get the equation number from the master node
1499  local_unknown = local_hang_eqn(
1500  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
1501  hang_weight2 = hang_info2_pt->master_weight(m2);
1502  }
1503  else
1504  {
1505  local_unknown =
1506  this->nodal_local_eqn(l2, u_nodal_index[i2]);
1507  hang_weight2 = 1.0;
1508  }
1509 
1510  // If the unknown is not pinned
1511  if (local_unknown >= 0)
1512  {
1513  switch (i2)
1514  {
1515  // radial component
1516  case 0:
1517  jacobian(local_eqn, local_unknown) +=
1518  (psif[l2] + r * dpsifdx(l2, 0)) * testp[l] * W *
1519  hang_weight * hang_weight2;
1520  break;
1521 
1522  // axial component
1523  case 1:
1524  jacobian(local_eqn, local_unknown) +=
1525  r * dpsifdx(l2, 1) * testp[l] * W * hang_weight *
1526  hang_weight2;
1527  break;
1528 
1529  // azimuthal component
1530  case 2:
1531  break;
1532  }
1533  }
1534  }
1535  }
1536  } // End of loop over nodes
1537 
1538  // NO PRESSURE CONTRIBUTIONS TO CONTINUITY EQUATION
1539  } // End of Jacobian calculation
1540  }
1541  }
1542  } // End of loop over pressure nodes
1543 
1544  } // End of loop over integration points
1545 
1546 
1547  } // End of fill_in_generic_residual_contribution_axi_nst(...)
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
G
JacobiRotation< float > G
Definition: Jacobi_makeGivens.cpp:2
J
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
m2
MatrixType m2(n_dims)
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::FiniteElement::nodal_value
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
oomph::FiniteElement::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
oomph::FiniteElement::nodal_position
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317
oomph::FiniteElement::dnodal_position_dt
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
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invro
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:456
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_local_eqn
virtual int p_local_eqn(const unsigned &n) const =0
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::viscosity_ratio
const double & viscosity_ratio() const
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:494
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::extrapolated_strain_rate
virtual void extrapolated_strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Definition: generalised_newtonian_axisym_navier_stokes_elements.cc:977
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::pshape_axi_nst
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_st
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:426
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Gamma
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:415
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::npres_axi_nst
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invfr
const double & re_invfr() const
Global inverse Froude number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:444
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Use_extrapolated_strainrate_to_compute_second_invariant
bool Use_extrapolated_strainrate_to_compute_second_invariant
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:183
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Constitutive_eqn_pt
GeneralisedNewtonianConstitutiveEquation< 3 > * Constitutive_eqn_pt
Pointer to the generalised Newtonian constitutive equation.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:179
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re
const double & re() const
Reynolds number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:420
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::du_dt_axi_nst
double du_dt_axi_nst(const unsigned &n, const unsigned &i) const
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:570
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_source_fct
double get_source_fct(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:301
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_body_force_axi_nst
virtual void get_body_force_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force fct at a given time and Eulerian position.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:237
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_nodal_index_axi_nst
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:618
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_axi_nst
virtual double p_axi_nst(const unsigned &n_p) const =0
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::density_ratio
const double & density_ratio() const
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:481
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_axi_nst
virtual double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::g
const Vector< double > & g() const
Vector of gravitational components.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:468
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ALE_is_disabled
bool ALE_is_disabled
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:188
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::strain_rate
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Definition: generalised_newtonian_axisym_navier_stokes_elements.cc:749
oomph::GeneralisedNewtonianConstitutiveEquation::viscosity
virtual double viscosity(const double &second_invariant_of_rate_of_strain_tensor)=0
oomph::GeneralisedNewtonianConstitutiveEquation::dviscosity_dinvariant
virtual double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)=0
oomph::HangInfo::nmaster
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::RefineableElement::local_hang_eqn
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Definition: refineable_elements.h:278
oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::pressure_node_pt
virtual Node * pressure_node_pt(const unsigned &n_p)
Definition: generalised_newtonian_refineable_axisym_navier_stokes_elements.h:53
oomph::TimeStepper::weight
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
oomph::TimeStepper::time_pt
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
oomph::Time::time
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
s
RealScalar s
Definition: level1_cplx_impl.h:130
k
char char char int int * k
Definition: level2_impl.h:374
DIM
#define DIM
Definition: linearised_navier_stokes_elements.h:44
GlobalParameters::Nintpt
unsigned Nintpt
Number of integration points for new integration scheme (if used)
Definition: stefan_boltzmann.cc:112
Global_Parameters::body_force
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Definition: axisym_linear_elasticity/cylinder/cylinder.cc:96
TestProblem::source
void source(const Vector< double > &x, Vector< double > &f)
Source function.
Definition: unstructured_two_d_circle.cc:46
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20
oomph::QuadTreeNames::W
@ W
Definition: quadtree.h:63
oomph::SecondInvariantHelper::second_invariant
double second_invariant(const DenseMatrix< double > &tensor)
Compute second invariant of tensor.
Definition: oomph_utilities.cc:163
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ALE_is_disabled, Global_Parameters::body_force(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Constitutive_eqn_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::density_ratio(), DIM, oomph::FiniteElement::dnodal_position_dt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_axi_nst(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::du_dt_axi_nst(), oomph::GeneralisedNewtonianConstitutiveEquation< DIM >::dviscosity_dinvariant(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::extrapolated_strain_rate(), G, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::g(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Gamma, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_body_force_axi_nst(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_source_fct(), 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(), n, GlobalParameters::Nintpt, oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::FiniteElement::node_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::npres_axi_nst(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_axi_nst(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_local_eqn(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_nodal_index_axi_nst(), pressure_node_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::pshape_axi_nst(), UniformPSDSelfTest::r, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invfr(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invro(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_st(), s, oomph::SecondInvariantHelper::second_invariant(), TestProblem::source(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::strain_rate(), oomph::Time::time(), oomph::TimeStepper::time_pt(), oomph::Data::time_stepper_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::u_index_axi_nst(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Use_extrapolated_strainrate_to_compute_second_invariant, oomph::GeneralisedNewtonianConstitutiveEquation< DIM >::viscosity(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::viscosity_ratio(), w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and oomph::TimeStepper::weight().

◆ further_build()

void oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::further_build ( )
inlinevirtual

Further build: pass the pointers down to the sons.

Reimplemented from oomph::RefineableElement.

Reimplemented in oomph::RefineableGeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement.

137  {
138  // Find the father element
139  RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations*
140  cast_father_element_pt = dynamic_cast<
141  RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations*>(
142  this->father_element_pt());
143 
144  // Set the viscosity ratio pointer
145  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
146  // Set the density ratio pointer
147  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
148  // Set pointer to global Reynolds number
149  this->Re_pt = cast_father_element_pt->re_pt();
150  // Set pointer to global Reynolds number x Strouhal number (=Womersley)
151  this->ReSt_pt = cast_father_element_pt->re_st_pt();
152  // Set pointer to global Reynolds number x inverse Froude number
153  this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
154  // Set pointer to the global Reynolds number x inverse Rossby number
155  this->ReInvRo_pt = cast_father_element_pt->re_invro_pt();
156  // Set pointer to global gravity Vector
157  this->G_pt = cast_father_element_pt->g_pt();
158 
159  // Set pointer to body force function
160  this->Body_force_fct_pt =
161  cast_father_element_pt->axi_nst_body_force_fct_pt();
162 
163  // Set pointer to volumetric source function
164  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
165 
166  // Set the ALE_is_disabled flag
167  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
168  }
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ReInvRo_pt
double * ReInvRo_pt
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:165
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Viscosity_Ratio_pt
double * Viscosity_Ratio_pt
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:145
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ReInvFr_pt
double * ReInvFr_pt
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:161
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Re_pt
double * Re_pt
Pointer to global Reynolds number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:154
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Density_Ratio_pt
double * Density_Ratio_pt
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:149
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Body_force_fct_pt
void(* Body_force_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to body force function.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:171
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ReSt_pt
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:157
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Source_fct_pt
double(* Source_fct_pt)(const double &time, const Vector< double > &x)
Pointer to volumetric source function.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:176
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::G_pt
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:168
oomph::RefineableElement::father_element_pt
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Definition: refineable_elements.h:539
oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations
RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations()
Empty Constructor.
Definition: generalised_newtonian_refineable_axisym_navier_stokes_elements.h:67

References oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ALE_is_disabled, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::axi_nst_body_force_fct_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Body_force_fct_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Density_Ratio_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::density_ratio_pt(), oomph::RefineableElement::father_element_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::G_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::g_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invfr_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invro_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Re_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_st_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ReInvFr_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ReInvRo_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ReSt_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::source_fct_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Source_fct_pt, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Viscosity_Ratio_pt, and oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::viscosity_ratio_pt().

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

◆ geometric_jacobian()

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

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

Reimplemented from oomph::ElementWithZ2ErrorEstimator.

130  {
131  return x[0];
132  }
plotDoE.x
list x
Definition: plotDoE.py:28

References plotDoE::x.

◆ get_dresidual_dnodal_coordinates()

void oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates ( RankThreeTensor< double > &  dresidual_dnodal_coordinates)
virtual

Compute derivatives of elemental residual vector with respect to nodal coordinates. This function computes these terms analytically and overwrites the default implementation in the FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}

Compute derivatives of elemental residual vector with respect to nodal coordinates. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij} Overloads the FD-based version in the FiniteElement base class.

Reimplemented from oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations.

1559  {
1560  // Create an Oomph Lib warning
1561  std::string warning_message = "This function has not been tested.\n";
1562  std::string function =
1563  "RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::\n";
1564  function += "get_dresidual_dnodal_coordinates(...)";
1565  OomphLibWarning(warning_message, function, OOMPH_EXCEPTION_LOCATION);
1566 
1567  // Return immediately if there are no dofs
1568  if (ndof() == 0)
1569  {
1570  return;
1571  }
1572 
1573  // Determine number of nodes in element
1574  const unsigned n_node = nnode();
1575 
1576  // Get continuous time from timestepper of first node
1577  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1578 
1579  // Determine number of pressure dofs in element
1580  const unsigned n_pres = this->npres_axi_nst();
1581 
1582  // Find the indices at which the local velocities are stored
1583  unsigned u_nodal_index[3];
1584  for (unsigned i = 0; i < 3; i++)
1585  {
1586  u_nodal_index[i] = this->u_index_axi_nst(i);
1587  }
1588 
1589  // Which nodal value represents the pressure? (Negative if pressure
1590  // is not based on nodal interpolation).
1591  const int p_index = this->p_nodal_index_axi_nst();
1592 
1593  // Local array of booleans that are true if the l-th pressure value is
1594  // hanging (avoid repeated virtual function calls)
1595  bool pressure_dof_is_hanging[n_pres];
1596 
1597  // If the pressure is stored at a node
1598  if (p_index >= 0)
1599  {
1600  // Read out whether the pressure is hanging
1601  for (unsigned l = 0; l < n_pres; ++l)
1602  {
1603  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
1604  }
1605  }
1606  // Otherwise the pressure is not stored at a node and so cannot hang
1607  else
1608  {
1609  for (unsigned l = 0; l < n_pres; ++l)
1610  {
1611  pressure_dof_is_hanging[l] = false;
1612  }
1613  }
1614 
1615  // Set up memory for the shape and test functions
1616  // Note that there are only two dimensions, r and z, in this problem
1617  Shape psif(n_node), testf(n_node);
1618  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
1619 
1620  // Set up memory for pressure shape and test functions
1621  Shape psip(n_pres), testp(n_pres);
1622 
1623  // Determine number of shape controlling nodes
1624  const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1625 
1626  // Deriatives of shape fct derivatives w.r.t. nodal coords
1627  RankFourTensor<double> d_dpsifdx_dX(2, n_shape_controlling_node, n_node, 2);
1628  RankFourTensor<double> d_dtestfdx_dX(
1629  2, n_shape_controlling_node, n_node, 2);
1630 
1631  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1632  DenseMatrix<double> dJ_dX(2, n_shape_controlling_node);
1633 
1634  // Derivatives of derivative of u w.r.t. nodal coords
1635  // Note that the entry d_dudx_dX(p,q,i,k) contains the derivative w.r.t.
1636  // nodal coordinate X_pq of du_i/dx_k. Since there are three velocity
1637  // components, i can take on the values 0, 1 and 2, while k and p only
1638  // take on the values 0 and 1 since there are only two spatial dimensions.
1639  RankFourTensor<double> d_dudx_dX(2, n_shape_controlling_node, 3, 2);
1640 
1641  // Derivatives of nodal velocities w.r.t. nodal coords:
1642  // Assumption: Interaction only local via no-slip so
1643  // X_pq only affects U_iq.
1644  // Note that the entry d_U_dX(p,q,i) contains the derivative w.r.t. nodal
1645  // coordinate X_pq of nodal value U_iq. In other words this entry is the
1646  // derivative of the i-th velocity component w.r.t. the p-th spatial
1647  // coordinate, taken at the q-th local node.
1648  RankThreeTensor<double> d_U_dX(2, n_shape_controlling_node, 3, 0.0);
1649 
1650  // Determine the number of integration points
1651  const unsigned n_intpt = integral_pt()->nweight();
1652 
1653  // Vector to hold local coordinates (two dimensions)
1654  Vector<double> s(2);
1655 
1656  // Get physical variables from element
1657  // (Reynolds number must be multiplied by the density ratio)
1658  const double scaled_re = this->re() * this->density_ratio();
1659  const double scaled_re_st = this->re_st() * this->density_ratio();
1660  const double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
1661  const double scaled_re_inv_ro = this->re_invro() * this->density_ratio();
1662  const double visc_ratio = this->viscosity_ratio();
1663  const Vector<double> G = this->g();
1664 
1665  // FD step
1666  double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
1667 
1668  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1669  // Assumption: Interaction only local via no-slip so
1670  // X_ij only affects U_ij.
1671  bool element_has_node_with_aux_node_update_fct = false;
1672 
1673  std::map<Node*, unsigned> local_shape_controlling_node_lookup =
1674  shape_controlling_node_lookup();
1675 
1676  // FD loop over shape-controlling nodes
1677  for (std::map<Node*, unsigned>::iterator it =
1678  local_shape_controlling_node_lookup.begin();
1679  it != local_shape_controlling_node_lookup.end();
1680  it++)
1681  {
1682  // Get pointer to q-th local shape-controlling node
1683  Node* nod_pt = it->first;
1684 
1685  // Get its number
1686  unsigned q = it->second;
1687 
1688  // Only compute if there's a node-update fct involved
1689  if (nod_pt->has_auxiliary_node_update_fct_pt())
1690  {
1691  element_has_node_with_aux_node_update_fct = true;
1692 
1693  // This functionality has not been tested so issue a warning
1694  // to this effect
1695  std::ostringstream warning_stream;
1696  warning_stream << "\nThe functionality to evaluate the additional"
1697  << "\ncontribution to the deriv of the residual eqn"
1698  << "\nw.r.t. the nodal coordinates which comes about"
1699  << "\nif a node's values are updated using an auxiliary"
1700  << "\nnode update function has NOT been tested for"
1701  << "\nrefineable axisymmetric Navier-Stokes elements."
1702  << "\nUse at your own risk" << std::endl;
1703  OomphLibWarning(warning_stream.str(),
1704  "RefineableGeneralisedNewtonianAxisymmetricNavierStokes"
1705  "Equations::get_dresidual_dnodal_coordinates",
1706  OOMPH_EXCEPTION_LOCATION);
1707 
1708  // Current nodal velocity
1709  Vector<double> u_ref(3);
1710  for (unsigned i = 0; i < 3; i++)
1711  {
1712  u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
1713  }
1714 
1715  // FD
1716  for (unsigned p = 0; p < 2; p++)
1717  {
1718  // Make backup
1719  double backup = nod_pt->x(p);
1720 
1721  // Do FD step. No node update required as we're
1722  // attacking the coordinate directly...
1723  nod_pt->x(p) += eps_fd;
1724 
1725  // Do auxiliary node update (to apply no slip)
1726  nod_pt->perform_auxiliary_node_update_fct();
1727 
1728  // Loop over velocity components
1729  for (unsigned i = 0; i < 3; i++)
1730  {
1731  // Evaluate
1732  d_U_dX(p, q, i) =
1733  (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
1734  }
1735 
1736  // Reset
1737  nod_pt->x(p) = backup;
1738 
1739  // Do auxiliary node update (to apply no slip)
1740  nod_pt->perform_auxiliary_node_update_fct();
1741  }
1742  }
1743  }
1744 
1745  // Integer to store the local equation number
1746  int local_eqn = 0;
1747 
1748  // Pointers to hang info object
1749  HangInfo* hang_info_pt = 0;
1750 
1751  // Loop over the integration points
1752  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1753  {
1754  // Assign values of s
1755  for (unsigned i = 0; i < 2; i++)
1756  {
1757  s[i] = integral_pt()->knot(ipt, i);
1758  }
1759 
1760  // Get the integral weight
1761  const double w = integral_pt()->weight(ipt);
1762 
1763  // Call the derivatives of the shape and test functions
1764  const double J =
1765  this->dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,
1766  psif,
1767  dpsifdx,
1768  d_dpsifdx_dX,
1769  testf,
1770  dtestfdx,
1771  d_dtestfdx_dX,
1772  dJ_dX);
1773 
1774  // Call the pressure shape and test functions
1775  this->pshape_axi_nst(s, psip, testp);
1776 
1777  // Allocate storage for the position and the derivative of the
1778  // mesh positions w.r.t. time
1779  Vector<double> interpolated_x(2, 0.0);
1780  Vector<double> mesh_velocity(2, 0.0);
1781 
1782  // Allocate storage for the pressure, velocity components and their
1783  // time and spatial derivatives
1784  double interpolated_p = 0.0;
1785  Vector<double> interpolated_u(3, 0.0);
1786  Vector<double> dudt(3, 0.0);
1787  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
1788 
1789  // Calculate pressure at integration point
1790  for (unsigned l = 0; l < n_pres; l++)
1791  {
1792  interpolated_p += this->p_axi_nst(l) * psip[l];
1793  }
1794 
1795  // Calculate velocities and derivatives at integration point
1796  // ---------------------------------------------------------
1797 
1798  // Loop over nodes
1799  for (unsigned l = 0; l < n_node; l++)
1800  {
1801  // Cache the shape function
1802  const double psif_ = psif(l);
1803 
1804  // Loop over the two coordinate directions
1805  for (unsigned i = 0; i < 2; i++)
1806  {
1807  interpolated_x[i] += nodal_position(l, i) * psif_;
1808  }
1809 
1810  // Loop over the three velocity directions
1811  for (unsigned i = 0; i < 3; i++)
1812  {
1813  // Get the nodal value
1814  const double u_value = nodal_value(l, u_nodal_index[i]);
1815  interpolated_u[i] += u_value * psif_;
1816  dudt[i] += this->du_dt_axi_nst(l, i) * psif_;
1817 
1818  // Loop over derivative directions
1819  for (unsigned j = 0; j < 2; j++)
1820  {
1821  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1822  }
1823  }
1824  }
1825 
1826  // Get the mesh velocity if ALE is enabled
1827  if (!this->ALE_is_disabled)
1828  {
1829  // Loop over nodes
1830  for (unsigned l = 0; l < n_node; l++)
1831  {
1832  // Loop over the two coordinate directions
1833  for (unsigned i = 0; i < 2; i++)
1834  {
1835  mesh_velocity[i] += this->dnodal_position_dt(l, i) * psif[l];
1836  }
1837  }
1838  }
1839 
1840  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1841 
1842  // Loop over shape-controlling nodes
1843  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1844  {
1845  // Loop over the two coordinate directions
1846  for (unsigned p = 0; p < 2; p++)
1847  {
1848  // Loop over the three velocity components
1849  for (unsigned i = 0; i < 3; i++)
1850  {
1851  // Loop over the two coordinate directions
1852  for (unsigned k = 0; k < 2; k++)
1853  {
1854  double aux = 0.0;
1855 
1856  // Loop over nodes and add contribution
1857  for (unsigned j = 0; j < n_node; j++)
1858  {
1859  aux +=
1860  nodal_value(j, u_nodal_index[i]) * d_dpsifdx_dX(p, q, j, k);
1861  }
1862  d_dudx_dX(p, q, i, k) = aux;
1863  }
1864  }
1865  }
1866  }
1867 
1868  // Get weight of actual nodal position/value in computation of mesh
1869  // velocity from positional/value time stepper
1870  const double pos_time_weight =
1871  node_pt(0)->position_time_stepper_pt()->weight(1, 0);
1872  const double val_time_weight =
1873  node_pt(0)->time_stepper_pt()->weight(1, 0);
1874 
1875  // Get the user-defined body force terms
1876  Vector<double> body_force(3);
1877  this->get_body_force_axi_nst(time, ipt, s, interpolated_x, body_force);
1878 
1879  // Get the user-defined source function
1880  const double source = this->get_source_fct(time, ipt, interpolated_x);
1881 
1882  // Get gradient of body force function
1883  DenseMatrix<double> d_body_force_dx(3, 2, 0.0);
1884  this->get_body_force_gradient_axi_nst(
1885  time, ipt, s, interpolated_x, d_body_force_dx);
1886 
1887  // Get gradient of source function
1888  Vector<double> source_gradient(2, 0.0);
1889  this->get_source_fct_gradient(time, ipt, interpolated_x, source_gradient);
1890 
1891  // Cache r (the first position component)
1892  const double r = interpolated_x[0];
1893 
1894  // Assemble shape derivatives
1895  // --------------------------
1896 
1897  // ==================
1898  // MOMENTUM EQUATIONS
1899  // ==================
1900 
1901  // Number of master nodes and storage for the weight of the shape function
1902  unsigned n_master = 1;
1903  double hang_weight = 1.0;
1904 
1905  // Loop over the test functions
1906  for (unsigned l = 0; l < n_node; l++)
1907  {
1908  // Cache the test function
1909  const double testf_ = testf[l];
1910 
1911  // Local boolean to indicate whether the node is hanging
1912  bool is_node_hanging = node_pt(l)->is_hanging();
1913 
1914  // If the node is hanging
1915  if (is_node_hanging)
1916  {
1917  hang_info_pt = node_pt(l)->hanging_pt();
1918 
1919  // Read out number of master nodes from hanging data
1920  n_master = hang_info_pt->nmaster();
1921  }
1922  // Otherwise the node is its own master
1923  else
1924  {
1925  n_master = 1;
1926  }
1927 
1928  // Loop over the master nodes
1929  for (unsigned m = 0; m < n_master; m++)
1930  {
1931  // --------------------------------
1932  // FIRST (RADIAL) MOMENTUM EQUATION
1933  // --------------------------------
1934 
1935  // Get the equation number
1936  // If the node is hanging
1937  if (is_node_hanging)
1938  {
1939  // Get the equation number from the master node
1940  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1941  u_nodal_index[0]);
1942  // Get the hang weight from the master node
1943  hang_weight = hang_info_pt->master_weight(m);
1944  }
1945  // If the node is not hanging
1946  else
1947  {
1948  // Local equation number
1949  local_eqn = this->nodal_local_eqn(l, u_nodal_index[0]);
1950 
1951  // Node contributes with full weight
1952  hang_weight = 1.0;
1953  }
1954 
1955  // If it's not a boundary condition
1956  if (local_eqn >= 0)
1957  {
1958  // Loop over the two coordinate directions
1959  for (unsigned p = 0; p < 2; p++)
1960  {
1961  // Loop over shape controlling nodes
1962  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1963  {
1964  // Residual x deriv of Jacobian
1965  // ----------------------------
1966 
1967  // Add the user-defined body force terms
1968  double sum = r * body_force[0] * testf_;
1969 
1970  // Add the gravitational body force term
1971  sum += r * scaled_re_inv_fr * testf_ * G[0];
1972 
1973  // Add the pressure gradient term
1974  sum += interpolated_p * (testf_ + r * dtestfdx(l, 0));
1975 
1976  // Add the stress tensor terms
1977  // The viscosity ratio needs to go in here to ensure
1978  // continuity of normal stress is satisfied even in flows
1979  // with zero pressure gradient!
1980  sum -= visc_ratio * r * (1.0 + Gamma[0]) *
1981  interpolated_dudx(0, 0) * dtestfdx(l, 0);
1982 
1983  sum -= visc_ratio * r *
1984  (interpolated_dudx(0, 1) +
1985  Gamma[0] * interpolated_dudx(1, 0)) *
1986  dtestfdx(l, 1);
1987 
1988  sum -= visc_ratio * (1.0 + Gamma[0]) * interpolated_u[0] *
1989  testf_ / r;
1990 
1991  // Add the du/dt term
1992  sum -= scaled_re_st * r * dudt[0] * testf_;
1993 
1994  // Add the convective terms
1995  sum -= scaled_re *
1996  (r * interpolated_u[0] * interpolated_dudx(0, 0) -
1997  interpolated_u[2] * interpolated_u[2] +
1998  r * interpolated_u[1] * interpolated_dudx(0, 1)) *
1999  testf_;
2000 
2001  // Add the mesh velocity terms
2002  if (!ALE_is_disabled)
2003  {
2004  for (unsigned k = 0; k < 2; k++)
2005  {
2006  sum += scaled_re_st * r * mesh_velocity[k] *
2007  interpolated_dudx(0, k) * testf_;
2008  }
2009  }
2010 
2011  // Add the Coriolis term
2012  sum += 2.0 * r * scaled_re_inv_ro * interpolated_u[2] * testf_;
2013 
2014  // Multiply through by deriv of Jacobian and integration weight
2015  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2016  sum * dJ_dX(p, q) * w * hang_weight;
2017 
2018  // Derivative of residual x Jacobian
2019  // ---------------------------------
2020 
2021  // Body force
2022  sum = r * d_body_force_dx(0, p) * psif[q] * testf_;
2023  if (p == 0)
2024  {
2025  sum += body_force[0] * psif[q] * testf_;
2026  }
2027 
2028  // Gravitational body force
2029  if (p == 0)
2030  {
2031  sum += scaled_re_inv_fr * G[0] * psif[q] * testf_;
2032  }
2033 
2034  // Pressure gradient term
2035  sum += r * interpolated_p * d_dtestfdx_dX(p, q, l, 0);
2036  if (p == 0)
2037  {
2038  sum += interpolated_p * psif[q] * dtestfdx(l, 0);
2039  }
2040 
2041  // Viscous terms
2042  sum -=
2043  r * visc_ratio *
2044  ((1.0 + Gamma[0]) *
2045  (d_dudx_dX(p, q, 0, 0) * dtestfdx(l, 0) +
2046  interpolated_dudx(0, 0) * d_dtestfdx_dX(p, q, l, 0)) +
2047  (d_dudx_dX(p, q, 0, 1) + Gamma[0] * d_dudx_dX(p, q, 1, 0)) *
2048  dtestfdx(l, 1) +
2049  (interpolated_dudx(0, 1) +
2050  Gamma[0] * interpolated_dudx(1, 0)) *
2051  d_dtestfdx_dX(p, q, l, 1));
2052  if (p == 0)
2053  {
2054  sum -=
2055  visc_ratio *
2056  ((1.0 + Gamma[0]) *
2057  (interpolated_dudx(0, 0) * psif[q] * dtestfdx(l, 0) -
2058  interpolated_u[0] * psif[q] * testf_ / (r * r)) +
2059  (interpolated_dudx(0, 1) +
2060  Gamma[0] * interpolated_dudx(1, 0)) *
2061  psif[q] * dtestfdx(l, 1));
2062  }
2063 
2064  // Convective terms, including mesh velocity
2065  for (unsigned k = 0; k < 2; k++)
2066  {
2067  double tmp = scaled_re * interpolated_u[k];
2068  if (!ALE_is_disabled)
2069  {
2070  tmp -= scaled_re_st * mesh_velocity[k];
2071  }
2072  sum -= r * tmp * d_dudx_dX(p, q, 0, k) * testf_;
2073  if (p == 0)
2074  {
2075  sum -= tmp * interpolated_dudx(0, k) * psif[q] * testf_;
2076  }
2077  }
2078  if (!ALE_is_disabled)
2079  {
2080  sum += scaled_re_st * r * pos_time_weight *
2081  interpolated_dudx(0, p) * psif[q] * testf_;
2082  }
2083 
2084  // du/dt term
2085  if (p == 0)
2086  {
2087  sum -= scaled_re_st * dudt[0] * psif[q] * testf_;
2088  }
2089 
2090  // Coriolis term
2091  if (p == 0)
2092  {
2093  sum += 2.0 * scaled_re_inv_ro * interpolated_u[2] * psif[q] *
2094  testf_;
2095  }
2096 
2097  // Multiply through by Jacobian and integration weight
2098  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2099  sum * J * w * hang_weight;
2100 
2101  } // End of loop over shape controlling nodes q
2102  } // End of loop over coordinate directions p
2103 
2104  // Derivs w.r.t. to nodal velocities
2105  // ---------------------------------
2106  if (element_has_node_with_aux_node_update_fct)
2107  {
2108  // Loop over local nodes
2109  for (unsigned q_local = 0; q_local < n_node; q_local++)
2110  {
2111  // Number of master nodes and storage for the weight of
2112  // the shape function
2113  unsigned n_master2 = 1;
2114  double hang_weight2 = 1.0;
2115  HangInfo* hang_info2_pt = 0;
2116 
2117  // Local boolean to indicate whether the node is hanging
2118  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2119 
2120  Node* actual_shape_controlling_node_pt = node_pt(q_local);
2121 
2122  // If the node is hanging
2123  if (is_node_hanging2)
2124  {
2125  hang_info2_pt = node_pt(q_local)->hanging_pt();
2126 
2127  // Read out number of master nodes from hanging data
2128  n_master2 = hang_info2_pt->nmaster();
2129  }
2130  // Otherwise the node is its own master
2131  else
2132  {
2133  n_master2 = 1;
2134  }
2135 
2136  // Loop over the master nodes
2137  for (unsigned mm = 0; mm < n_master2; mm++)
2138  {
2139  if (is_node_hanging2)
2140  {
2141  actual_shape_controlling_node_pt =
2142  hang_info2_pt->master_node_pt(mm);
2143  hang_weight2 = hang_info2_pt->master_weight(mm);
2144  }
2145 
2146  // Find the corresponding number
2147  unsigned q = local_shape_controlling_node_lookup
2148  [actual_shape_controlling_node_pt];
2149 
2150  // Loop over the two coordinate directions
2151  for (unsigned p = 0; p < 2; p++)
2152  {
2153  // Contribution from deriv of first velocity component
2154  double tmp = scaled_re_st * r * val_time_weight *
2155  psif[q_local] * testf_;
2156  tmp += r * scaled_re * interpolated_dudx(0, 0) *
2157  psif[q_local] * testf_;
2158 
2159  for (unsigned k = 0; k < 2; k++)
2160  {
2161  double tmp2 = scaled_re * interpolated_u[k];
2162  if (!ALE_is_disabled)
2163  {
2164  tmp2 -= scaled_re_st * mesh_velocity[k];
2165  }
2166  tmp += r * tmp2 * dpsifdx(q_local, k) * testf_;
2167  }
2168 
2169  tmp += r * visc_ratio * (1.0 + Gamma[0]) *
2170  dpsifdx(q_local, 0) * dtestfdx(l, 0);
2171  tmp +=
2172  r * visc_ratio * dpsifdx(q_local, 1) * dtestfdx(l, 1);
2173  tmp += visc_ratio * (1.0 + Gamma[0]) * psif[q_local] *
2174  testf_ / r;
2175 
2176  // Multiply through by dU_0q/dX_pq
2177  double sum = -d_U_dX(p, q_local, 0) * tmp;
2178 
2179  // Contribution from deriv of second velocity component
2180  tmp = scaled_re * r * interpolated_dudx(0, 1) *
2181  psif[q_local] * testf_;
2182  tmp += r * visc_ratio * Gamma[0] * dpsifdx(q_local, 0) *
2183  dtestfdx(l, 1);
2184 
2185  // Multiply through by dU_1q/dX_pq
2186  sum -= d_U_dX(p, q, 1) * tmp;
2187 
2188  // Contribution from deriv of third velocity component
2189  tmp =
2190  2.0 *
2191  (r * scaled_re_inv_ro + scaled_re * interpolated_u[2]) *
2192  psif[q_local] * testf_;
2193 
2194  // Multiply through by dU_2q/dX_pq
2195  sum += d_U_dX(p, q, 2) * tmp;
2196 
2197  // Multiply through by Jacobian and integration weight
2198  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2199  sum * J * w * hang_weight * hang_weight2;
2200  }
2201  } // End of loop over master nodes
2202  } // End of loop over local nodes
2203  } // End of if(element_has_node_with_aux_node_update_fct)
2204  } // End of if it's not a boundary condition
2205 
2206  // --------------------------------
2207  // SECOND (AXIAL) MOMENTUM EQUATION
2208  // --------------------------------
2209 
2210  // Get the equation number
2211  // If the node is hanging
2212  if (is_node_hanging)
2213  {
2214  // Get the equation number from the master node
2215  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
2216  u_nodal_index[1]);
2217  // Get the hang weight from the master node
2218  hang_weight = hang_info_pt->master_weight(m);
2219  }
2220  // If the node is not hanging
2221  else
2222  {
2223  // Local equation number
2224  local_eqn = this->nodal_local_eqn(l, u_nodal_index[1]);
2225 
2226  // Node contributes with full weight
2227  hang_weight = 1.0;
2228  }
2229 
2230  // If it's not a boundary condition
2231  if (local_eqn >= 0)
2232  {
2233  // Loop over the two coordinate directions
2234  for (unsigned p = 0; p < 2; p++)
2235  {
2236  // Loop over shape controlling nodes
2237  for (unsigned q = 0; q < n_shape_controlling_node; q++)
2238  {
2239  // Residual x deriv of Jacobian
2240  // ----------------------------
2241 
2242  // Add the user-defined body force terms
2243  double sum = r * body_force[1] * testf_;
2244 
2245  // Add the gravitational body force term
2246  sum += r * scaled_re_inv_fr * testf_ * G[1];
2247 
2248  // Add the pressure gradient term
2249  sum += r * interpolated_p * dtestfdx(l, 1);
2250 
2251  // Add the stress tensor terms
2252  // The viscosity ratio needs to go in here to ensure
2253  // continuity of normal stress is satisfied even in flows
2254  // with zero pressure gradient!
2255  sum -= visc_ratio * r *
2256  (interpolated_dudx(1, 0) +
2257  Gamma[1] * interpolated_dudx(0, 1)) *
2258  dtestfdx(l, 0);
2259 
2260  sum -= visc_ratio * r * (1.0 + Gamma[1]) *
2261  interpolated_dudx(1, 1) * dtestfdx(l, 1);
2262 
2263  // Add the du/dt term
2264  sum -= scaled_re_st * r * dudt[1] * testf_;
2265 
2266  // Add the convective terms
2267  sum -= scaled_re *
2268  (r * interpolated_u[0] * interpolated_dudx(1, 0) +
2269  r * interpolated_u[1] * interpolated_dudx(1, 1)) *
2270  testf_;
2271 
2272  // Add the mesh velocity terms
2273  if (!ALE_is_disabled)
2274  {
2275  for (unsigned k = 0; k < 2; k++)
2276  {
2277  sum += scaled_re_st * r * mesh_velocity[k] *
2278  interpolated_dudx(1, k) * testf_;
2279  }
2280  }
2281 
2282  // Multiply through by deriv of Jacobian and integration weight
2283  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2284  sum * dJ_dX(p, q) * w * hang_weight;
2285 
2286  // Derivative of residual x Jacobian
2287  // ---------------------------------
2288 
2289  // Body force
2290  sum = r * d_body_force_dx(1, p) * psif[q] * testf_;
2291  if (p == 0)
2292  {
2293  sum += body_force[1] * psif[q] * testf_;
2294  }
2295 
2296  // Gravitational body force
2297  if (p == 0)
2298  {
2299  sum += scaled_re_inv_fr * G[1] * psif[q] * testf_;
2300  }
2301 
2302  // Pressure gradient term
2303  sum += r * interpolated_p * d_dtestfdx_dX(p, q, l, 1);
2304  if (p == 0)
2305  {
2306  sum += interpolated_p * psif[q] * dtestfdx(l, 1);
2307  }
2308 
2309  // Viscous terms
2310  sum -=
2311  r * visc_ratio *
2312  ((d_dudx_dX(p, q, 1, 0) + Gamma[1] * d_dudx_dX(p, q, 0, 1)) *
2313  dtestfdx(l, 0) +
2314  (interpolated_dudx(1, 0) +
2315  Gamma[1] * interpolated_dudx(0, 1)) *
2316  d_dtestfdx_dX(p, q, l, 0) +
2317  (1.0 + Gamma[1]) * d_dudx_dX(p, q, 1, 1) * dtestfdx(l, 1) +
2318  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
2319  d_dtestfdx_dX(p, q, l, 1));
2320  if (p == 0)
2321  {
2322  sum -=
2323  visc_ratio * ((interpolated_dudx(1, 0) +
2324  Gamma[1] * interpolated_dudx(0, 1)) *
2325  psif[q] * dtestfdx(l, 0) +
2326  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
2327  psif[q] * dtestfdx(l, 1));
2328  }
2329 
2330  // Convective terms, including mesh velocity
2331  for (unsigned k = 0; k < 2; k++)
2332  {
2333  double tmp = scaled_re * interpolated_u[k];
2334  if (!ALE_is_disabled)
2335  {
2336  tmp -= scaled_re_st * mesh_velocity[k];
2337  }
2338  sum -= r * tmp * d_dudx_dX(p, q, 1, k) * testf_;
2339  if (p == 0)
2340  {
2341  sum -= tmp * interpolated_dudx(1, k) * psif[q] * testf_;
2342  }
2343  }
2344  if (!ALE_is_disabled)
2345  {
2346  sum += scaled_re_st * r * pos_time_weight *
2347  interpolated_dudx(1, p) * psif[q] * testf_;
2348  }
2349 
2350  // du/dt term
2351  if (p == 0)
2352  {
2353  sum -= scaled_re_st * dudt[1] * psif[q] * testf_;
2354  }
2355 
2356  // Multiply through by Jacobian and integration weight
2357  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2358  sum * J * w * hang_weight;
2359 
2360  } // End of loop over shape controlling nodes q
2361  } // End of loop over coordinate directions p
2362 
2363  // Derivs w.r.t. to nodal velocities
2364  // ---------------------------------
2365  if (element_has_node_with_aux_node_update_fct)
2366  {
2367  // Loop over local nodes
2368  for (unsigned q_local = 0; q_local < n_node; q_local++)
2369  {
2370  // Number of master nodes and storage for the weight of
2371  // the shape function
2372  unsigned n_master2 = 1;
2373  double hang_weight2 = 1.0;
2374  HangInfo* hang_info2_pt = 0;
2375 
2376  // Local boolean to indicate whether the node is hanging
2377  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2378 
2379  Node* actual_shape_controlling_node_pt = node_pt(q_local);
2380 
2381  // If the node is hanging
2382  if (is_node_hanging2)
2383  {
2384  hang_info2_pt = node_pt(q_local)->hanging_pt();
2385 
2386  // Read out number of master nodes from hanging data
2387  n_master2 = hang_info2_pt->nmaster();
2388  }
2389  // Otherwise the node is its own master
2390  else
2391  {
2392  n_master2 = 1;
2393  }
2394 
2395  // Loop over the master nodes
2396  for (unsigned mm = 0; mm < n_master2; mm++)
2397  {
2398  if (is_node_hanging2)
2399  {
2400  actual_shape_controlling_node_pt =
2401  hang_info2_pt->master_node_pt(mm);
2402  hang_weight2 = hang_info2_pt->master_weight(mm);
2403  }
2404 
2405  // Find the corresponding number
2406  unsigned q = local_shape_controlling_node_lookup
2407  [actual_shape_controlling_node_pt];
2408 
2409  // Loop over the two coordinate directions
2410  for (unsigned p = 0; p < 2; p++)
2411  {
2412  // Contribution from deriv of first velocity component
2413  double tmp = scaled_re * r * interpolated_dudx(1, 0) *
2414  psif[q_local] * testf_;
2415  tmp += r * visc_ratio * Gamma[1] * dpsifdx(q_local, 1) *
2416  dtestfdx(l, 0);
2417 
2418  // Multiply through by dU_0q/dX_pq
2419  double sum = -d_U_dX(p, q, 0) * tmp;
2420 
2421  // Contribution from deriv of second velocity component
2422  tmp = scaled_re_st * r * val_time_weight * psif[q_local] *
2423  testf_;
2424  tmp += r * scaled_re * interpolated_dudx(1, 1) *
2425  psif[q_local] * testf_;
2426 
2427  for (unsigned k = 0; k < 2; k++)
2428  {
2429  double tmp2 = scaled_re * interpolated_u[k];
2430  if (!ALE_is_disabled)
2431  {
2432  tmp2 -= scaled_re_st * mesh_velocity[k];
2433  }
2434  tmp += r * tmp2 * dpsifdx(q_local, k) * testf_;
2435  }
2436  tmp +=
2437  r * visc_ratio *
2438  (dpsifdx(q_local, 0) * dtestfdx(l, 0) +
2439  (1.0 + Gamma[1]) * dpsifdx(q_local, 1) * dtestfdx(l, 1));
2440 
2441  // Multiply through by dU_1q/dX_pq
2442  sum -= d_U_dX(p, q, 1) * tmp;
2443 
2444  // Multiply through by Jacobian and integration weight
2445  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2446  sum * J * w * hang_weight * hang_weight2;
2447  }
2448  } // End of loop over master nodes
2449  } // End of loop over local nodes
2450  } // End of if(element_has_node_with_aux_node_update_fct)
2451  } // End of if it's not a boundary condition
2452 
2453  // -----------------------------------
2454  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
2455  // -----------------------------------
2456 
2457  // Get the equation number
2458  // If the node is hanging
2459  if (is_node_hanging)
2460  {
2461  // Get the equation number from the master node
2462  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
2463  u_nodal_index[2]);
2464  // Get the hang weight from the master node
2465  hang_weight = hang_info_pt->master_weight(m);
2466  }
2467  // If the node is not hanging
2468  else
2469  {
2470  // Local equation number
2471  local_eqn = this->nodal_local_eqn(l, u_nodal_index[2]);
2472 
2473  // Node contributes with full weight
2474  hang_weight = 1.0;
2475  }
2476 
2477  // If it's not a boundary condition
2478  if (local_eqn >= 0)
2479  {
2480  // Loop over the two coordinate directions
2481  for (unsigned p = 0; p < 2; p++)
2482  {
2483  // Loop over shape controlling nodes
2484  for (unsigned q = 0; q < n_shape_controlling_node; q++)
2485  {
2486  // Residual x deriv of Jacobian
2487  // ----------------------------
2488 
2489  // Add the user-defined body force terms
2490  double sum = r * body_force[2] * testf_;
2491 
2492  // Add the gravitational body force term
2493  sum += r * scaled_re_inv_fr * testf_ * G[2];
2494 
2495  // There is NO pressure gradient term
2496 
2497  // Add in the stress tensor terms
2498  // The viscosity ratio needs to go in here to ensure
2499  // continuity of normal stress is satisfied even in flows
2500  // with zero pressure gradient!
2501  sum -=
2502  visc_ratio *
2503  (r * interpolated_dudx(2, 0) - Gamma[0] * interpolated_u[2]) *
2504  dtestfdx(l, 0);
2505 
2506  sum -=
2507  visc_ratio * r * interpolated_dudx(2, 1) * dtestfdx(l, 1);
2508 
2509  sum -= visc_ratio *
2510  ((interpolated_u[2] / r) -
2511  Gamma[0] * interpolated_dudx(2, 0)) *
2512  testf_;
2513 
2514  // Add the du/dt term
2515  sum -= scaled_re_st * r * dudt[2] * testf_;
2516 
2517  // Add the convective terms
2518  sum -= scaled_re *
2519  (r * interpolated_u[0] * interpolated_dudx(2, 0) +
2520  interpolated_u[0] * interpolated_u[2] +
2521  r * interpolated_u[1] * interpolated_dudx(2, 1)) *
2522  testf_;
2523 
2524  // Add the mesh velocity terms
2525  if (!ALE_is_disabled)
2526  {
2527  for (unsigned k = 0; k < 2; k++)
2528  {
2529  sum += scaled_re_st * r * mesh_velocity[k] *
2530  interpolated_dudx(2, k) * testf_;
2531  }
2532  }
2533 
2534  // Add the Coriolis term
2535  sum -= 2.0 * r * scaled_re_inv_ro * interpolated_u[0] * testf_;
2536 
2537  // Multiply through by deriv of Jacobian and integration weight
2538  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2539  sum * dJ_dX(p, q) * w * hang_weight;
2540 
2541  // Derivative of residual x Jacobian
2542  // ---------------------------------
2543 
2544  // Body force
2545  sum = r * d_body_force_dx(2, p) * psif[q] * testf_;
2546  if (p == 0)
2547  {
2548  sum += body_force[2] * psif[q] * testf_;
2549  }
2550 
2551  // Gravitational body force
2552  if (p == 0)
2553  {
2554  sum += scaled_re_inv_fr * G[2] * psif[q] * testf_;
2555  }
2556 
2557  // There is NO pressure gradient term
2558 
2559  // Viscous terms
2560  sum -= visc_ratio * r *
2561  (d_dudx_dX(p, q, 2, 0) * dtestfdx(l, 0) +
2562  interpolated_dudx(2, 0) * d_dtestfdx_dX(p, q, l, 0) +
2563  d_dudx_dX(p, q, 2, 1) * dtestfdx(l, 1) +
2564  interpolated_dudx(2, 1) * d_dtestfdx_dX(p, q, l, 1));
2565 
2566  sum += visc_ratio * Gamma[0] * d_dudx_dX(p, q, 2, 0) * testf_;
2567 
2568  if (p == 0)
2569  {
2570  sum -= visc_ratio *
2571  (interpolated_dudx(2, 0) * psif[q] * dtestfdx(l, 0) +
2572  interpolated_dudx(2, 1) * psif[q] * dtestfdx(l, 1) +
2573  interpolated_u[2] * psif[q] * testf_ / (r * r));
2574  }
2575 
2576  // Convective terms, including mesh velocity
2577  for (unsigned k = 0; k < 2; k++)
2578  {
2579  double tmp = scaled_re * interpolated_u[k];
2580  if (!ALE_is_disabled)
2581  {
2582  tmp -= scaled_re_st * mesh_velocity[k];
2583  }
2584  sum -= r * tmp * d_dudx_dX(p, q, 2, k) * testf_;
2585  if (p == 0)
2586  {
2587  sum -= tmp * interpolated_dudx(2, k) * psif[q] * testf_;
2588  }
2589  }
2590  if (!ALE_is_disabled)
2591  {
2592  sum += scaled_re_st * r * pos_time_weight *
2593  interpolated_dudx(2, p) * psif[q] * testf_;
2594  }
2595 
2596  // du/dt term
2597  if (p == 0)
2598  {
2599  sum -= scaled_re_st * dudt[2] * psif[q] * testf_;
2600  }
2601 
2602  // Coriolis term
2603  if (p == 0)
2604  {
2605  sum -= 2.0 * scaled_re_inv_ro * interpolated_u[0] * psif[q] *
2606  testf_;
2607  }
2608 
2609  // Multiply through by Jacobian and integration weight
2610  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2611  sum * J * w * hang_weight;
2612 
2613  } // End of loop over shape controlling nodes q
2614  } // End of loop over coordinate directions p
2615 
2616  // Derivs w.r.t. to nodal velocities
2617  // ---------------------------------
2618  if (element_has_node_with_aux_node_update_fct)
2619  {
2620  // Loop over local nodes
2621  for (unsigned q_local = 0; q_local < n_node; q_local++)
2622  {
2623  // Number of master nodes and storage for the weight of
2624  // the shape function
2625  unsigned n_master2 = 1;
2626  double hang_weight2 = 1.0;
2627  HangInfo* hang_info2_pt = 0;
2628 
2629  // Local boolean to indicate whether the node is hanging
2630  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2631 
2632  Node* actual_shape_controlling_node_pt = node_pt(q_local);
2633 
2634  // If the node is hanging
2635  if (is_node_hanging2)
2636  {
2637  hang_info2_pt = node_pt(q_local)->hanging_pt();
2638 
2639  // Read out number of master nodes from hanging data
2640  n_master2 = hang_info2_pt->nmaster();
2641  }
2642  // Otherwise the node is its own master
2643  else
2644  {
2645  n_master2 = 1;
2646  }
2647 
2648  // Loop over the master nodes
2649  for (unsigned mm = 0; mm < n_master2; mm++)
2650  {
2651  if (is_node_hanging2)
2652  {
2653  actual_shape_controlling_node_pt =
2654  hang_info2_pt->master_node_pt(mm);
2655  hang_weight2 = hang_info2_pt->master_weight(mm);
2656  }
2657 
2658  // Find the corresponding number
2659  unsigned q = local_shape_controlling_node_lookup
2660  [actual_shape_controlling_node_pt];
2661 
2662  // Loop over the two coordinate directions
2663  for (unsigned p = 0; p < 2; p++)
2664  {
2665  // Contribution from deriv of first velocity component
2666  double tmp = (2.0 * r * scaled_re_inv_ro +
2667  r * scaled_re * interpolated_dudx(2, 0) -
2668  scaled_re * interpolated_u[2]) *
2669  psif[q_local] * testf_;
2670 
2671  // Multiply through by dU_0q/dX_pq
2672  double sum = -d_U_dX(p, q, 0) * tmp;
2673 
2674  // Contribution from deriv of second velocity component
2675  // Multiply through by dU_1q/dX_pq
2676  sum -= d_U_dX(p, q, 1) * scaled_re * r *
2677  interpolated_dudx(2, 1) * psif[q_local] * testf_;
2678 
2679  // Contribution from deriv of third velocity component
2680  tmp = scaled_re_st * r * val_time_weight * psif[q_local] *
2681  testf_;
2682  tmp -=
2683  scaled_re * interpolated_u[0] * psif[q_local] * testf_;
2684 
2685  for (unsigned k = 0; k < 2; k++)
2686  {
2687  double tmp2 = scaled_re * interpolated_u[k];
2688  if (!ALE_is_disabled)
2689  {
2690  tmp2 -= scaled_re_st * mesh_velocity[k];
2691  }
2692  tmp += r * tmp2 * dpsifdx(q_local, k) * testf_;
2693  }
2694  tmp +=
2695  visc_ratio *
2696  (r * dpsifdx(q_local, 0) - Gamma[0] * psif[q_local]) *
2697  dtestfdx(l, 0);
2698  tmp +=
2699  r * visc_ratio * dpsifdx(q_local, 1) * dtestfdx(l, 1);
2700  tmp +=
2701  visc_ratio *
2702  ((psif[q_local] / r) - Gamma[0] * dpsifdx(q_local, 0)) *
2703  testf_;
2704 
2705  // Multiply through by dU_2q/dX_pq
2706  sum -= d_U_dX(p, q, 2) * tmp;
2707 
2708  // Multiply through by Jacobian and integration weight
2709  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2710  sum * J * w * hang_weight * hang_weight2;
2711  }
2712  } // End of loop over master nodes
2713  } // End of loop over local nodes
2714  } // End of if(element_has_node_with_aux_node_update_fct)
2715  } // End of if it's not a boundary condition
2716  } // End of loop over master nodes
2717  } // End of loop over test functions
2718 
2719 
2720  // ===================
2721  // CONTINUITY EQUATION
2722  // ===================
2723 
2724  // Loop over the shape functions
2725  for (unsigned l = 0; l < n_pres; l++)
2726  {
2727  // If the pressure dof is hanging
2728  if (pressure_dof_is_hanging[l])
2729  {
2730  // Pressure dof is hanging so it must be nodal-based
2731  // Get the hang info object
2732  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
2733 
2734  // Get the number of master nodes from the pressure node
2735  n_master = hang_info_pt->nmaster();
2736  }
2737  // Otherwise the node is its own master
2738  else
2739  {
2740  n_master = 1;
2741  }
2742 
2743  // Loop over the master nodes
2744  for (unsigned m = 0; m < n_master; m++)
2745  {
2746  // Get the number of the unknown
2747  // If the pressure dof is hanging
2748  if (pressure_dof_is_hanging[l])
2749  {
2750  // Get the local equation from the master node
2751  local_eqn =
2752  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
2753  // Get the weight for the node
2754  hang_weight = hang_info_pt->master_weight(m);
2755  }
2756  else
2757  {
2758  local_eqn = this->p_local_eqn(l);
2759  hang_weight = 1.0;
2760  }
2761 
2762  // Cache the test function
2763  const double testp_ = testp[l];
2764 
2765  // If not a boundary conditions
2766  if (local_eqn >= 0)
2767  {
2768  // Loop over the two coordinate directions
2769  for (unsigned p = 0; p < 2; p++)
2770  {
2771  // Loop over shape controlling nodes
2772  for (unsigned q = 0; q < n_shape_controlling_node; q++)
2773  {
2774  // Residual x deriv of Jacobian
2775  //-----------------------------
2776 
2777  // Source term
2778  double aux = -r * source;
2779 
2780  // Gradient terms
2781  aux += (interpolated_u[0] + r * interpolated_dudx(0, 0) +
2782  r * interpolated_dudx(1, 1));
2783 
2784  // Multiply through by deriv of Jacobian and integration weight
2785  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2786  aux * dJ_dX(p, q) * testp_ * w * hang_weight;
2787 
2788  // Derivative of residual x Jacobian
2789  // ---------------------------------
2790 
2791  // Gradient of source function
2792  aux = -r * source_gradient[p] * psif[q];
2793  if (p == 0)
2794  {
2795  aux -= source * psif[q];
2796  }
2797 
2798  // Gradient terms
2799  aux += r * (d_dudx_dX(p, q, 0, 0) + d_dudx_dX(p, q, 1, 1));
2800  if (p == 0)
2801  {
2802  aux += (interpolated_dudx(0, 0) + interpolated_dudx(1, 1)) *
2803  psif[q];
2804  }
2805 
2806  // Multiply through by Jacobian and integration weight
2807  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2808  aux * testp_ * J * w * hang_weight;
2809  }
2810  }
2811 
2812  // Derivs w.r.t. to nodal velocities
2813  // ---------------------------------
2814  if (element_has_node_with_aux_node_update_fct)
2815  {
2816  // Loop over local nodes
2817  for (unsigned q_local = 0; q_local < n_node; q_local++)
2818  {
2819  // Number of master nodes and storage for the weight of
2820  // the shape function
2821  unsigned n_master2 = 1;
2822  double hang_weight2 = 1.0;
2823  HangInfo* hang_info2_pt = 0;
2824 
2825  // Local boolean to indicate whether the node is hanging
2826  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2827 
2828  Node* actual_shape_controlling_node_pt = node_pt(q_local);
2829 
2830  // If the node is hanging
2831  if (is_node_hanging2)
2832  {
2833  hang_info2_pt = node_pt(q_local)->hanging_pt();
2834 
2835  // Read out number of master nodes from hanging data
2836  n_master2 = hang_info2_pt->nmaster();
2837  }
2838  // Otherwise the node is its own master
2839  else
2840  {
2841  n_master2 = 1;
2842  }
2843 
2844  // Loop over the master nodes
2845  for (unsigned mm = 0; mm < n_master2; mm++)
2846  {
2847  if (is_node_hanging2)
2848  {
2849  actual_shape_controlling_node_pt =
2850  hang_info2_pt->master_node_pt(mm);
2851  hang_weight2 = hang_info2_pt->master_weight(mm);
2852  }
2853 
2854  // Find the corresponding number
2855  unsigned q = local_shape_controlling_node_lookup
2856  [actual_shape_controlling_node_pt];
2857 
2858  // Loop over the two coordinate directions
2859  for (unsigned p = 0; p < 2; p++)
2860  {
2861  double aux = d_U_dX(p, q, 0) *
2862  (psif[q_local] + r * dpsifdx(q_local, 0)) +
2863  d_U_dX(p, q, 1) * r * dpsifdx(q_local, 1);
2864 
2865  // Multiply through by Jacobian, test function and
2866  // integration weight
2867  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2868  aux * testp_ * J * w * hang_weight * hang_weight2;
2869  }
2870  } // End of loop over (mm) master nodes
2871  } // End of loop over local nodes q_local
2872  } // End of if(element_has_node_with_aux_node_update_fct)
2873  } // End of if it's not a boundary condition
2874  } // End of loop over (m) master nodes
2875  } // End of loop over shape functions for continuity eqn
2876 
2877  } // End of loop over integration points
2878 
2879  } // End of get_dresidual_dnodal_coordinates(...)
p
float * p
Definition: Tutorial_Map_using.cpp:9
oomph::GeneralisedElement::Default_fd_jacobian_step
static double Default_fd_jacobian_step
Definition: elements.h:1198
oomph::GeneralisedElement::ndof
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_source_fct_gradient
virtual void get_source_fct_gradient(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:320
oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_body_force_gradient_axi_nst
virtual void get_body_force_gradient_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:262
oomph::Node::position_time_stepper_pt
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
oomph::RefineableElement::nshape_controlling_nodes
unsigned nshape_controlling_nodes()
Definition: refineable_elements.h:627
oomph::RefineableElement::shape_controlling_node_lookup
std::map< Node *, unsigned > shape_controlling_node_lookup()
Definition: refineable_elements.h:636
tmp
Eigen::Matrix< Scalar, Dynamic, Dynamic, ColMajor > tmp
Definition: level3_impl.h:365
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286

References oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ALE_is_disabled, Global_Parameters::body_force(), oomph::GeneralisedElement::Default_fd_jacobian_step, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::density_ratio(), oomph::FiniteElement::dnodal_position_dt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_axi_nst(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::du_dt_axi_nst(), G, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::g(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Gamma, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_body_force_axi_nst(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_body_force_gradient_axi_nst(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_source_fct(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_source_fct_gradient(), oomph::Node::hanging_pt(), oomph::Node::has_auxiliary_node_update_fct_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, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::GeneralisedElement::ndof(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::FiniteElement::node_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::npres_axi_nst(), oomph::RefineableElement::nshape_controlling_nodes(), oomph::Integral::nweight(), OOMPH_EXCEPTION_LOCATION, p, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_axi_nst(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_local_eqn(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_nodal_index_axi_nst(), oomph::Node::perform_auxiliary_node_update_fct(), oomph::Node::position_time_stepper_pt(), pressure_node_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::pshape_axi_nst(), Eigen::numext::q, UniformPSDSelfTest::r, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invfr(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invro(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_st(), s, oomph::RefineableElement::shape_controlling_node_lookup(), TestProblem::source(), oomph::Global_string_for_annotation::string(), oomph::Time::time(), oomph::TimeStepper::time_pt(), oomph::Data::time_stepper_pt(), tmp, oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::u_index_axi_nst(), oomph::Data::value_pt(), oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::viscosity_ratio(), w, oomph::Integral::weight(), oomph::TimeStepper::weight(), and oomph::Node::x().

◆ get_Z2_flux()

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

84  {
85  // Specify the number of velocity dimensions
86  unsigned DIM = 3;
87 
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 " << num_entries
95  << std::endl;
96  throw OomphLibError(error_message.str(),
97  OOMPH_CURRENT_FUNCTION,
98  OOMPH_EXCEPTION_LOCATION);
99  }
100 #endif
101 
102 
103  // Get strain rate matrix
104  DenseMatrix<double> strainrate(DIM);
105  this->strain_rate(s, strainrate);
106 
107  // Pack into flux Vector
108  unsigned icount = 0;
109 
110  // Start with diagonal terms
111  for (unsigned i = 0; i < DIM; i++)
112  {
113  flux[icount] = strainrate(i, i);
114  icount++;
115  }
116 
117  // Off diagonals row by row
118  for (unsigned i = 0; i < DIM; i++)
119  {
120  for (unsigned j = i + 1; j < DIM; j++)
121  {
122  flux[icount] = strainrate(i, j);
123  icount++;
124  }
125  }
126  }
ProblemParameters::flux
void flux(const double &time, const Vector< double > &x, double &flux)
Get flux applied along boundary x=0.
Definition: pretend_melt.cc:59

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

◆ num_Z2_flux_terms()

unsigned oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::num_Z2_flux_terms ( )
inlinevirtual

Number of 'flux' terms for Z2 error estimation.

Implements oomph::ElementWithZ2ErrorEstimator.

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

◆ pin_elemental_redundant_nodal_pressure_dofs()

virtual void oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::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::RefineableGeneralisedNewtonianAxisymmetricQTaylorHoodElement.

63 {}

Referenced by pin_redundant_nodal_pressures().

◆ pin_redundant_nodal_pressures()

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

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

322  {
323  // Loop over all elements to brutally pin all nodal pressure degrees of
324  // freedom
325  unsigned n_element = element_pt.size();
326  for (unsigned e = 0; e < n_element; e++)
327  {
328  dynamic_cast<
329  RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations*>(
330  element_pt[e])
331  ->pin_elemental_redundant_nodal_pressure_dofs();
332  }
333  }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::pin_elemental_redundant_nodal_pressure_dofs
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Definition: generalised_newtonian_refineable_axisym_navier_stokes_elements.h:63

References e(), and pin_elemental_redundant_nodal_pressure_dofs().

◆ pressure_node_pt()

virtual Node* oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::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::RefineableGeneralisedNewtonianAxisymmetricQTaylorHoodElement.

54  {
55  return NULL;
56  }

Referenced by fill_in_generic_residual_contribution_axi_nst(), and get_dresidual_dnodal_coordinates().

◆ unpin_all_pressure_dofs()

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

Unpin all pressure dofs in elements listed in vector.

338  {
339  // Loop over all elements to brutally unpin all nodal pressure degrees of
340  // freedom and internal pressure degrees of freedom
341  unsigned n_element = element_pt.size();
342  for (unsigned e = 0; e < n_element; e++)
343  {
344  dynamic_cast<
345  RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations*>(
346  element_pt[e])
347  ->unpin_elemental_pressure_dofs();
348  }
349  }
oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::unpin_elemental_pressure_dofs
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.

References e(), and unpin_elemental_pressure_dofs().

◆ unpin_elemental_pressure_dofs()

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

Unpin all pressure dofs in the element.

Implemented in oomph::RefineableGeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement, and oomph::RefineableGeneralisedNewtonianAxisymmetricQTaylorHoodElement.

Referenced by unpin_all_pressure_dofs().

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