oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations Class Referenceabstract

#include <generalised_newtonian_axisym_navier_stokes_elements.h>

+ Inheritance diagram for oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations:

Public Member Functions

 GeneralisedNewtonianAxisymmetricNavierStokesEquations ()
 Constructor: NULL the body force and source function. More...
 
const doublere () const
 Reynolds number. More...
 
const doublere_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 doublere_invfr () const
 Global inverse Froude number. More...
 
double *& re_invfr_pt ()
 Pointer to global inverse Froude number. More...
 
const doublere_invro () const
 Global Reynolds number multiplied by inverse Rossby number. More...
 
double *& re_invro_pt ()
 Pointer to global 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 doubledensity_ratio () const
 
double *& density_ratio_pt ()
 Pointer to Density ratio. More...
 
const doubleviscosity_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)
 
virtual void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
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...
 
virtual void dinterpolated_u_axi_nst_ddata (const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
 
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)
 
MacroElementmacro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void assign_nodal_local_eqn_numbers (const bool &store_local_dof_pt)
 
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 unsigned nvertex_node () const
 
virtual Nodevertex_node_pt (const unsigned &j) const
 
virtual Nodeconstruct_node (const unsigned &n)
 
virtual Nodeconstruct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual Nodeconstruct_boundary_node (const unsigned &n)
 
virtual Nodeconstruct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
int get_node_number (Node *const &node_pt) const
 
virtual Nodeget_node_at_local_coordinate (const Vector< double > &s) const
 
double raw_nodal_value (const unsigned &n, const unsigned &i) const
 
double raw_nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
unsigned dim () const
 
virtual ElementGeometry::ElementGeometry element_geometry () const
 Return the geometry type of the element (either Q or T usually). More...
 
virtual double interpolated_x (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated coordinate x[i] at local coordinate s. More...
 
virtual double interpolated_x (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_x (const Vector< double > &s, Vector< double > &x) const
 Return FE interpolated position x[] at local coordinate s as Vector. More...
 
virtual void interpolated_x (const unsigned &t, const Vector< double > &s, Vector< double > &x) const
 
virtual double interpolated_dxdt (const Vector< double > &s, const unsigned &i, const unsigned &t)
 
virtual void interpolated_dxdt (const Vector< double > &s, const unsigned &t, Vector< double > &dxdt)
 
unsigned ngeom_data () const
 
Datageom_data_pt (const unsigned &j)
 
void position (const Vector< double > &zeta, Vector< double > &r) const
 
void position (const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
 
void dposition_dt (const Vector< double > &zeta, const unsigned &t, Vector< double > &drdt)
 
virtual double zeta_nodal (const unsigned &n, const unsigned &k, const unsigned &i) const
 
void interpolated_zeta (const Vector< double > &s, Vector< double > &zeta) const
 
void locate_zeta (const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
 
virtual void node_update ()
 
virtual void identify_field_data_for_interactions (std::set< std::pair< Data *, unsigned >> &paired_field_data)
 
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 ()
 
TimeSteppertime_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 
- Public Member Functions inherited from oomph::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 Attributes

void(*&)(const double &time, const Vector< double > &x, Vector< double > &faxi_nst_body_force_fct_pt ()
 Access function for the body-force pointer. More...
 
double(*&)(const double &time, const Vector< double > &xsource_fct_pt ()
 Access function for the source-function pointer. More...
 

Static Public Attributes

static Vector< doubleGamma
 Vector to decide whether the stress-divergence form is used or not. More...
 
- Static Public Attributes inherited from oomph::FiniteElement
static double Tolerance_for_singular_jacobian = 1.0e-16
 Tolerance below which the jacobian is considered singular. More...
 
static bool Accept_negative_jacobian = false
 
static bool Suppress_output_while_checking_for_inverted_elements
 
- Static Public Attributes inherited from oomph::GeneralisedElement
static bool Suppress_warning_about_repeated_internal_data
 
static bool Suppress_warning_about_repeated_external_data = true
 
static double Default_fd_jacobian_step = 1.0e-8
 

Protected Member Functions

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)
 
virtual void fill_in_generic_residual_contribution_axi_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 
virtual 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)
 
- Protected Member Functions inherited from oomph::FiniteElement
virtual void assemble_local_to_eulerian_jacobian (const DShape &dpsids, DenseMatrix< double > &jacobian) const
 
virtual void assemble_local_to_eulerian_jacobian2 (const DShape &d2psids, DenseMatrix< double > &jacobian2) const
 
virtual void assemble_eulerian_base_vectors (const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
 
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 double local_to_eulerian_mapping_diagonal (const DShape &dpsids, DenseMatrix< double > &jacobian, 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
 
virtual void fill_in_jacobian_from_nodal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
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 Attributes

doubleViscosity_Ratio_pt
 
doubleDensity_Ratio_pt
 
doubleRe_pt
 Pointer to global Reynolds number. More...
 
doubleReSt_pt
 Pointer to global Reynolds number x Strouhal number (=Womersley) More...
 
doubleReInvFr_pt
 
doubleReInvRo_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
MacroElementMacro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepperGeom_object_time_stepper_pt
 

Static Protected Attributes

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< doubleDummy_matrix
 
static std::deque< double * > Dof_pt_deque
 

Static Private Attributes

static int Pressure_not_stored_at_node
 
static double Default_Physical_Constant_Value
 
static double Default_Physical_Ratio_Value
 
static Vector< doubleDefault_Gravity_vector
 Static default value for the gravity vector. More...
 

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

Detailed Description

A class for elements that solve the unsteady axisymmetric Navier–Stokes equations in cylindrical polar coordinates, \( x_0^* = r^*\) and \( x_1^* = z^* \) with \( \partial / \partial \theta = 0 \). We're solving for the radial, axial and azimuthal (swirl) velocities, \( u_0^* = u_r^*(r^*,z^*,t^*) = u^*(r^*,z^*,t^*), \ u_1^* = u_z^*(r^*,z^*,t^*) = w^*(r^*,z^*,t^*)\) and \( u_2^* = u_\theta^*(r^*,z^*,t^*) = v^*(r^*,z^*,t^*) \), respectively, and the pressure \( p(r^*,z^*,t^*) \). This class contains the generic maths – any concrete implementation must be derived from this.

In dimensional form the axisymmetric Navier-Stokes equations are given by the momentum equations (for the \( r^* \), \( z^* \) and \( \theta \) directions, respectively)

\[ \rho\left(\frac{\partial u^*}{\partial t^*} + {u^*}\frac{\partial u^*}{\partial r^*} - \frac{{v^*}^2}{r^*} + {w^*}\frac{\partial u^*}{\partial z^*} \right) = B_r^*\left(r^*,z^*,t^*\right)+ \rho G_r^*+ \frac{1}{r^*} \frac{\partial\left({r^*}\sigma_{rr}^*\right)}{\partial r^*} - \frac{\sigma_{\theta\theta}^*}{r^*} + \frac{\partial\sigma_{rz}^*}{\partial z^*}, \]

\[ \rho\left(\frac{\partial w^*}{\partial t^*} + {u^*}\frac{\partial w^*}{\partial r^*} + {w^*}\frac{\partial w^*}{\partial z^*} \right) = B_z^*\left(r^*,z^*,t^*\right)+\rho G_z^*+ \frac{1}{r^*}\frac{\partial\left({r^*}\sigma_{zr}^*\right)}{\partial r^*} + \frac{\partial\sigma_{zz}^*}{\partial z^*}, \]

\[ \rho\left(\frac{\partial v^*}{\partial t^*} + {u^*}\frac{\partial v^*}{\partial r^*} + \frac{u^* v^*}{r^*} +{w^*}\frac{\partial v^*}{\partial z^*} \right)= B_\theta^*\left(r^*,z^*,t^*\right)+ \rho G_\theta^*+ \frac{1}{r^*}\frac{\partial\left({r^*}\sigma_{\theta r}^*\right)}{\partial r^*} + \frac{\sigma_{r\theta}^*}{r^*} + \frac{\partial\sigma_{\theta z}^*}{\partial z^*}, \]

and

\[ \frac{1}{r^*}\frac{\partial\left(r^*u^*\right)}{\partial r^*} + \frac{\partial w^*}{\partial z^*} = Q^*. \]

The dimensional, symmetric stress tensor is defined as:

\[ \sigma_{rr}^* = -p^* + 2\mu\frac{\partial u^*}{\partial r^*}, \qquad \sigma_{\theta\theta}^* = -p^* +2\mu\frac{u^*}{r^*}, \]

\[ \sigma_{zz}^* = -p^* + 2\mu\frac{\partial w^*}{\partial z^*}, \qquad \sigma_{rz}^* = \mu\left(\frac{\partial u^*}{\partial z^*} + \frac{\partial w^*}{\partial r^*}\right), \]

\[ \sigma_{\theta r}^* = \mu r^*\frac{\partial}{\partial r^*} \left(\frac{v^*}{r^*}\right), \qquad \sigma_{\theta z}^* = \mu\frac{\partial v^*}{\partial z^*}. \]

Here, the (dimensional) velocity components are denoted by \( u^* \), \( w^* \) and \( v^* \) for the radial, axial and azimuthal velocities, respectively, and we have split the body force into two components: A constant vector \( \rho \ G_i^* \) which typically represents gravitational forces; and a variable body force, \( B_i^*(r^*,z^*,t^*) \). \( Q^*(r^*,z^*,t^*) \) is a volumetric source term for the continuity equation and is typically equal to zero.

We non-dimensionalise the equations, using problem-specific reference quantities for the velocity, \( U \), length, \( L \), and time, \( T \), and scale the constant body force vector on the gravitational acceleration, \( g \), so that

\[ u^* = U\, u, \qquad w^* = U\, w, \qquad v^* = U\, v, \]

\[ r^* = L\, r, \qquad z^* = L\, z, \qquad t^* = T\, t, \]

\[ G_i^* = g\, G_i, \qquad B_i^* = \frac{U\mu_{ref}}{L^2}\, B_i, \qquad p^* = \frac{\mu_{ref} U}{L}\, p, \qquad Q^* = \frac{U}{L}\, Q. \]

where we note that the pressure and the variable body force have been non-dimensionalised on the viscous scale. \( \mu_{ref} \) and \( \rho_{ref} \) (used below) are reference values for the fluid viscosity and density, respectively. In single-fluid problems, they are identical to the viscosity \( \mu \) and density \( \rho \) of the (one and only) fluid in the problem.

The non-dimensional form of the axisymmetric Navier-Stokes equations is then given by

\[ R_{\rho} Re\left(St\frac{\partial u}{\partial t} + {u}\frac{\partial u}{\partial r} - \frac{{v}^2}{r} + {w}\frac{\partial u}{\partial z} \right) = B_r\left(r,z,t\right)+ R_\rho \frac{Re}{Fr} G_r + \frac{1}{r} \frac{\partial\left({r}\sigma_{rr}\right)}{\partial r} - \frac{\sigma_{\theta\theta}}{r} + \frac{\partial\sigma_{rz}}{\partial z}, \]

\[ R_{\rho} Re\left(St\frac{\partial w}{\partial t} + {u}\frac{\partial w}{\partial r} + {w}\frac{\partial w}{\partial z} \right) = B_z\left(r,z,t\right)+ R_\rho \frac{Re}{Fr} G_z+ \frac{1}{r}\frac{\partial\left({r}\sigma_{zr}\right)}{\partial r} + \frac{\partial\sigma_{zz}}{\partial z}, \]

\[ R_{\rho} Re\left(St\frac{\partial v}{\partial t} + {u}\frac{\partial v}{\partial r} + \frac{u v}{r} +{w}\frac{\partial v}{\partial z} \right)= B_\theta\left(r,z,t\right)+ R_\rho \frac{Re}{Fr} G_\theta+ \frac{1}{r}\frac{\partial\left({r}\sigma_{\theta r}\right)}{\partial r} + \frac{\sigma_{r\theta}}{r} + \frac{\partial\sigma_{\theta z}}{\partial z}, \]

and

\[ \frac{1}{r}\frac{\partial\left(ru\right)}{\partial r} + \frac{\partial w}{\partial z} = Q. \]

Here the non-dimensional, symmetric stress tensor is defined as:

\[ \sigma_{rr} = -p + 2R_\mu \frac{\partial u}{\partial r}, \qquad \sigma_{\theta\theta} = -p +2R_\mu \frac{u}{r}, \]

\[ \sigma_{zz} = -p + 2R_\mu \frac{\partial w}{\partial z}, \qquad \sigma_{rz} = R_\mu \left(\frac{\partial u}{\partial z} + \frac{\partial w}{\partial r}\right), \]

\[ \sigma_{\theta r} = R_\mu r \frac{\partial}{\partial r}\left(\frac{v}{r}\right), \qquad \sigma_{\theta z} = R_\mu \frac{\partial v}{\partial z}. \]

and the dimensionless parameters

\[ Re = \frac{UL\rho_{ref}}{\mu_{ref}}, \qquad St = \frac{L}{UT}, \qquad Fr = \frac{U^2}{gL}, \]

are the Reynolds number, Strouhal number and Froude number respectively. \( R_\rho=\rho/\rho_{ref} \) and \( R_\mu(T) =\mu(T)/\mu_{ref}\) represent the ratios of the fluid's density and its dynamic viscosity, relative to the density and viscosity values used to form the non-dimensional parameters (By default, \( R_\rho = R_\mu = 1 \); other values tend to be used in problems involving multiple fluids).

Constructor & Destructor Documentation

◆ GeneralisedNewtonianAxisymmetricNavierStokesEquations()

oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::GeneralisedNewtonianAxisymmetricNavierStokesEquations ( )
inline

Constructor: NULL the body force and source function.

394  : Body_force_fct_pt(0),
395  Source_fct_pt(0),
396  Constitutive_eqn_pt(new NewtonianConstitutiveEquation<3>),
398  ALE_is_disabled(false)
399  {
400  // Set all the Physical parameter pointers to the default value zero
406  // Set the Physical ratios to the default value of 1
409  }
double * ReInvRo_pt
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:165
bool Use_extrapolated_strainrate_to_compute_second_invariant
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:183
GeneralisedNewtonianConstitutiveEquation< 3 > * Constitutive_eqn_pt
Pointer to the generalised Newtonian constitutive equation.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:179
double * Viscosity_Ratio_pt
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:145
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:134
double * ReInvFr_pt
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:161
double * Re_pt
Pointer to global Reynolds number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:154
double * Density_Ratio_pt
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:149
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
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:157
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
static double Default_Physical_Constant_Value
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:127
static double Default_Physical_Ratio_Value
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:131
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:168
bool ALE_is_disabled
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:188

References Default_Gravity_vector, Default_Physical_Constant_Value, Default_Physical_Ratio_Value, Density_Ratio_pt, G_pt, Re_pt, ReInvFr_pt, ReInvRo_pt, ReSt_pt, and Viscosity_Ratio_pt.

Member Function Documentation

◆ compute_error() [1/2]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::compute_error ( std::ostream &  outfile,
FiniteElement::SteadyExactSolutionFctPt  exact_soln_pt,
double error,
double norm 
)
virtual

Validate against exact solution. Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element

Validate against exact velocity solution Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element.

Reimplemented from oomph::FiniteElement.

246  {
247  error = 0.0;
248  norm = 0.0;
249 
250  // Vector of local coordinates
251  Vector<double> s(2);
252 
253  // Vector for coordintes
254  Vector<double> x(2);
255 
256  // Set the value of Nintpt
257  unsigned Nintpt = integral_pt()->nweight();
258 
259  outfile << "ZONE" << std::endl;
260 
261  // Exact solution Vector (u,v,w,p)
262  Vector<double> exact_soln(4);
263 
264  // Loop over the integration points
265  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
266  {
267  // Assign values of s
268  for (unsigned i = 0; i < 2; i++)
269  {
270  s[i] = integral_pt()->knot(ipt, i);
271  }
272 
273  // Get the integral weight
274  double w = integral_pt()->weight(ipt);
275 
276  // Get jacobian of mapping
277  double J = J_eulerian(s);
278 
279  // Get x position as Vector
280  interpolated_x(s, x);
281 
282  // Premultiply the weights and the Jacobian and r
283  double W = w * J * x[0];
284 
285  // Get exact solution at this point
286  (*exact_soln_pt)(x, exact_soln);
287 
288  // Velocity error
289  for (unsigned i = 0; i < 3; i++)
290  {
291  norm += exact_soln[i] * exact_soln[i] * W;
294  }
295 
296  // Output x,y,...,u_exact
297  for (unsigned i = 0; i < 2; i++)
298  {
299  outfile << x[i] << " ";
300  }
301 
302  // Output x,y,u_error,v_error,w_error
303  for (unsigned i = 0; i < 3; i++)
304  {
305  outfile << exact_soln[i] - interpolated_u_axi_nst(s, i) << " ";
306  }
307 
308  outfile << std::endl;
309  }
310  }
int i
Definition: BiCGSTAB_step_by_step.cpp:9
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
RowVector3d w
Definition: Matrix_resize_int.cpp:3
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual double J_eulerian(const Vector< double > &s) const
Definition: elements.cc:4103
void interpolated_u_axi_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:944
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
RealScalar s
Definition: level1_cplx_impl.h:130
unsigned Nintpt
Number of integration points for new integration scheme (if used)
Definition: stefan_boltzmann.cc:112
void exact_soln(const double &time, const Vector< double > &x, Vector< double > &soln)
Definition: unstructured_two_d_curved.cc:301
int error
Definition: calibrate.py:297
@ W
Definition: quadtree.h:63
list x
Definition: plotDoE.py:28

References calibrate::error, ProblemParameters::exact_soln(), i, oomph::FiniteElement::integral_pt(), interpolated_u_axi_nst(), oomph::FiniteElement::interpolated_x(), J, oomph::FiniteElement::J_eulerian(), oomph::Integral::knot(), GlobalParameters::Nintpt, oomph::Integral::nweight(), s, w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and plotDoE::x.

◆ compute_error() [2/2]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::compute_error ( std::ostream &  outfile,
FiniteElement::UnsteadyExactSolutionFctPt  exact_soln_pt,
const double time,
double error,
double norm 
)
virtual

Validate against exact solution at given time Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element

Validate against exact velocity solution at given time. Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element.

Reimplemented from oomph::FiniteElement.

169  {
170  error = 0.0;
171  norm = 0.0;
172 
173  // Vector of local coordinates
174  Vector<double> s(2);
175 
176  // Vector for coordintes
177  Vector<double> x(2);
178 
179  // Set the value of Nintpt
180  unsigned Nintpt = integral_pt()->nweight();
181 
182  outfile << "ZONE" << std::endl;
183 
184  // Exact solution Vector (u,v,w,p)
185  Vector<double> exact_soln(4);
186 
187  // Loop over the integration points
188  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
189  {
190  // Assign values of s
191  for (unsigned i = 0; i < 2; i++)
192  {
193  s[i] = integral_pt()->knot(ipt, i);
194  }
195 
196  // Get the integral weight
197  double w = integral_pt()->weight(ipt);
198 
199  // Get jacobian of mapping
200  double J = J_eulerian(s);
201 
202  // Get x position as Vector
203  interpolated_x(s, x);
204 
205  // Premultiply the weights and the Jacobian and r
206  double W = w * J * x[0];
207 
208  // Get exact solution at this point
209  (*exact_soln_pt)(time, x, exact_soln);
210 
211  // Velocity error
212  for (unsigned i = 0; i < 3; i++)
213  {
214  norm += exact_soln[i] * exact_soln[i] * W;
217  }
218 
219  // Output x,y,...,u_exact
220  for (unsigned i = 0; i < 2; i++)
221  {
222  outfile << x[i] << " ";
223  }
224 
225  // Output x,y,[z],u_error,v_error,[w_error]
226  for (unsigned i = 0; i < 3; i++)
227  {
228  outfile << exact_soln[i] - interpolated_u_axi_nst(s, i) << " ";
229  }
230 
231  outfile << std::endl;
232  }
233  }

References calibrate::error, ProblemParameters::exact_soln(), i, oomph::FiniteElement::integral_pt(), interpolated_u_axi_nst(), oomph::FiniteElement::interpolated_x(), J, oomph::FiniteElement::J_eulerian(), oomph::Integral::knot(), GlobalParameters::Nintpt, oomph::Integral::nweight(), s, w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and plotDoE::x.

◆ compute_physical_size()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::compute_physical_size ( ) const
inlinevirtual

Compute the volume of the element.

Reimplemented from oomph::FiniteElement.

1252  {
1253  // Initialise result
1254  double result = 0.0;
1255 
1256  // Figure out the global (Eulerian) spatial dimension of the
1257  // element by checking the Eulerian dimension of the nodes
1258  const unsigned n_dim_eulerian = nodal_dimension();
1259 
1260  // Allocate Vector of global Eulerian coordinates
1261  Vector<double> x(n_dim_eulerian);
1262 
1263  // Set the value of n_intpt
1264  const unsigned n_intpt = integral_pt()->nweight();
1265 
1266  // Vector of local coordinates
1267  const unsigned n_dim = dim();
1268  Vector<double> s(n_dim);
1269 
1270  // Loop over the integration points
1271  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1272  {
1273  // Assign the values of s
1274  for (unsigned i = 0; i < n_dim; i++)
1275  {
1276  s[i] = integral_pt()->knot(ipt, i);
1277  }
1278 
1279  // Assign the values of the global Eulerian coordinates
1280  for (unsigned i = 0; i < n_dim_eulerian; i++)
1281  {
1282  x[i] = interpolated_x(s, i);
1283  }
1284 
1285  // Get the integral weight
1286  const double w = integral_pt()->weight(ipt);
1287 
1288  // Get Jacobian of mapping
1289  const double J = J_eulerian(s);
1290 
1291  // The integrand at the current integration point is r
1292  result += x[0] * w * J;
1293  }
1294 
1295  // Multiply by 2pi (integrating in azimuthal direction)
1296  return (2.0 * MathematicalConstants::Pi * result);
1297  }
unsigned dim() const
Definition: elements.h:2611
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
const double Pi
50 digits from maple
Definition: oomph_utilities.h:157

References oomph::FiniteElement::dim(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), J, oomph::FiniteElement::J_eulerian(), oomph::Integral::knot(), oomph::FiniteElement::nodal_dimension(), oomph::Integral::nweight(), oomph::MathematicalConstants::Pi, s, w, oomph::Integral::weight(), and plotDoE::x.

◆ constitutive_eqn_pt()

GeneralisedNewtonianConstitutiveEquation<3>*& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::constitutive_eqn_pt ( )
inline

Access function for the constitutive equation pointer.

521  {
522  return Constitutive_eqn_pt;
523  }

References Constitutive_eqn_pt.

◆ density_ratio()

const double& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::density_ratio ( ) const
inline

◆ density_ratio_pt()

double*& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::density_ratio_pt ( )
inline

Pointer to Density ratio.

488  {
489  return Density_Ratio_pt;
490  }

References Density_Ratio_pt.

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

◆ dinterpolated_u_axi_nst_ddata()

virtual void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::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. The function is virtual so that it can be overloaded in the refineable version

Reimplemented in oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations.

1032  {
1033  // Find number of nodes
1034  unsigned n_node = nnode();
1035  // Local shape function
1036  Shape psi(n_node);
1037  // Find values of shape function
1038  shape(s, psi);
1039 
1040  // Find the index at which the velocity component is stored
1041  const unsigned u_nodal_index = u_index_axi_nst(i);
1042 
1043  // Find the number of dofs associated with interpolated u
1044  unsigned n_u_dof = 0;
1045  for (unsigned l = 0; l < n_node; l++)
1046  {
1047  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
1048  // If it's positive add to the count
1049  if (global_eqn >= 0)
1050  {
1051  ++n_u_dof;
1052  }
1053  }
1054 
1055  // Now resize the storage schemes
1056  du_ddata.resize(n_u_dof, 0.0);
1057  global_eqn_number.resize(n_u_dof, 0);
1058 
1059  // Loop over th nodes again and set the derivatives
1060  unsigned count = 0;
1061  // Loop over the local nodes and sum
1062  for (unsigned l = 0; l < n_node; l++)
1063  {
1064  // Get the global equation number
1065  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
1066  if (global_eqn >= 0)
1067  {
1068  // Set the global equation number
1069  global_eqn_number[count] = global_eqn;
1070  // Set the derivative with respect to the unknown
1071  du_ddata[count] = psi[l];
1072  // Increase the counter
1073  ++count;
1074  }
1075  }
1076  }
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void shape(const Vector< double > &s, Shape &psi) const =0
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual unsigned u_index_axi_nst(const unsigned &i) const
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:546

References oomph::Data::eqn_number(), i, oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), s, oomph::FiniteElement::shape(), and u_index_axi_nst().

◆ disable_ALE()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::disable_ALE ( )
inlinevirtual

Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!

Reimplemented from oomph::FiniteElement.

600  {
601  ALE_is_disabled = true;
602  }

References ALE_is_disabled.

◆ dissipation() [1/2]

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::dissipation ( ) const

Return integral of dissipation over element.

603  {
604  throw OomphLibError(
605  "Check the dissipation calculation for axisymmetric NSt",
608 
609  // Initialise
610  double diss = 0.0;
611 
612  // Set the value of Nintpt
613  unsigned Nintpt = integral_pt()->nweight();
614 
615  // Set the Vector to hold local coordinates
616  Vector<double> s(2);
617 
618  // Loop over the integration points
619  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
620  {
621  // Assign values of s
622  for (unsigned i = 0; i < 2; i++)
623  {
624  s[i] = integral_pt()->knot(ipt, i);
625  }
626 
627  // Get the integral weight
628  double w = integral_pt()->weight(ipt);
629 
630  // Get Jacobian of mapping
631  double J = J_eulerian(s);
632 
633  // Get strain rate matrix
634  DenseMatrix<double> strainrate(3, 3);
635  strain_rate(s, strainrate);
636 
637  // Initialise
638  double local_diss = 0.0;
639  for (unsigned i = 0; i < 3; i++)
640  {
641  for (unsigned j = 0; j < 3; j++)
642  {
643  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
644  }
645  }
646 
647  diss += local_diss * w * J;
648  }
649 
650  return diss;
651  }
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Definition: generalised_newtonian_axisym_navier_stokes_elements.cc:749
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References i, oomph::FiniteElement::integral_pt(), J, j, oomph::FiniteElement::J_eulerian(), oomph::Integral::knot(), GlobalParameters::Nintpt, oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, strain_rate(), w, and oomph::Integral::weight().

◆ dissipation() [2/2]

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::dissipation ( const Vector< double > &  s) const

Return dissipation at local coordinate s.

722  {
723  throw OomphLibError(
724  "Check the dissipation calculation for axisymmetric NSt",
727 
728  // Get strain rate matrix
729  DenseMatrix<double> strainrate(3, 3);
730  strain_rate(s, strainrate);
731 
732  // Initialise
733  double local_diss = 0.0;
734  for (unsigned i = 0; i < 3; i++)
735  {
736  for (unsigned j = 0; j < 3; j++)
737  {
738  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
739  }
740  }
741 
742  return local_diss;
743  }

References i, j, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, and strain_rate().

◆ dshape_and_dtest_eulerian_at_knot_axi_nst() [1/2]

virtual double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::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
protectedpure virtual

Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of mapping (J). Also compute derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.

Implemented in oomph::GeneralisedNewtonianAxisymmetricTTaylorHoodElement, oomph::GeneralisedNewtonianAxisymmetricTCrouzeixRaviartElement, oomph::GeneralisedNewtonianAxisymmetricQTaylorHoodElement, and oomph::GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement.

◆ dshape_and_dtest_eulerian_at_knot_axi_nst() [2/2]

◆ dshape_and_dtest_eulerian_axi_nst()

virtual double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_axi_nst ( const Vector< double > &  s,
Shape psi,
DShape dpsidx,
Shape test,
DShape dtestdx 
) const
protectedpure virtual

Compute the shape functions and derivatives w.r.t. global coords at local coordinate s. Return Jacobian of mapping between local and global coordinates.

Implemented in oomph::GeneralisedNewtonianAxisymmetricTTaylorHoodElement, oomph::GeneralisedNewtonianAxisymmetricTCrouzeixRaviartElement, oomph::GeneralisedNewtonianAxisymmetricQTaylorHoodElement, and oomph::GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement.

◆ du_dt_axi_nst()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::du_dt_axi_nst ( const unsigned n,
const unsigned i 
) const
inline

i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes.

571  {
572  // Get the data's timestepper
573  TimeStepper* time_stepper_pt = this->node_pt(n)->time_stepper_pt();
574 
575  // Initialise dudt
576  double dudt = 0.0;
577  // Loop over the timesteps, if there is a non Steady timestepper
578  if (!time_stepper_pt->is_steady())
579  {
580  // Get the index at which the velocity is stored
581  const unsigned u_nodal_index = u_index_axi_nst(i);
582 
583  // Number of timsteps (past & present)
584  const unsigned n_time = time_stepper_pt->ntstorage();
585 
586  // Add the contributions to the time derivative
587  for (unsigned t = 0; t < n_time; t++)
588  {
589  dudt +=
590  time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
591  }
592  }
593 
594  return dudt;
595  }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
TimeStepper *& time_stepper_pt()
Definition: geom_objects.h:192
unsigned ntstorage() const
Definition: timesteppers.h:601
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
bool is_steady() const
Definition: timesteppers.h:389
t
Definition: plotPSD.py:36

References i, oomph::TimeStepper::is_steady(), n, oomph::FiniteElement::nodal_value(), oomph::FiniteElement::node_pt(), oomph::TimeStepper::ntstorage(), plotPSD::t, oomph::GeomObject::time_stepper_pt(), oomph::Data::time_stepper_pt(), u_index_axi_nst(), and oomph::TimeStepper::weight().

Referenced by fill_in_generic_dresidual_contribution_axi_nst(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates(), and interpolated_dudt_axi_nst().

◆ enable_ALE()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::enable_ALE ( )
inlinevirtual

(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative. Note: By default, ALE is enabled, at the expense of possibly creating unnecessary work in problems where the mesh is, in fact, stationary.

Reimplemented from oomph::FiniteElement.

609  {
610  ALE_is_disabled = false;
611  }

References ALE_is_disabled.

◆ extrapolated_strain_rate() [1/2]

virtual void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::extrapolated_strain_rate ( const unsigned ipt,
DenseMatrix< double > &  strain_rate 
) const
inlinevirtual

Extrapolated strain-rate tensor: \( e_{ij} \) where \( i,j = r,z,\theta \) (in that order) based on the previous time steps evaluated at ipt-th integration point

667  {
668  // Set the Vector to hold local coordinates (two dimensions)
669  Vector<double> s(2);
670  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
672  }
virtual void extrapolated_strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Definition: generalised_newtonian_axisym_navier_stokes_elements.cc:977

References extrapolated_strain_rate(), i, oomph::FiniteElement::integral_pt(), oomph::Integral::knot(), s, and strain_rate().

◆ extrapolated_strain_rate() [2/2]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::extrapolated_strain_rate ( const Vector< double > &  s,
DenseMatrix< double > &  strainrate 
) const
virtual

Extrapolated strain-rate tensor: \( e_{ij} \) where \( i,j = r,z,\theta \) (in that order) based on the previous time steps evaluated at local coordinate s

Get strain-rate tensor: \( e_{ij} \) where \( i,j = r,z,\theta \) (in that order). Extrapolated from history values.

979  {
980 #ifdef PARANOID
981  if ((strainrate.ncol() != 3) || (strainrate.nrow() != 3))
982  {
983  std::ostringstream error_message;
984  error_message << "The strain rate has incorrect dimensions "
985  << strainrate.ncol() << " x " << strainrate.nrow()
986  << " Not 3" << std::endl;
987 
988  throw OomphLibError(
989  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
990  }
991 #endif
992 
993  // Get strain rates at two previous times
994  DenseMatrix<double> strain_rate_minus_two(3, 3);
995  strain_rate(2, s, strain_rate_minus_two);
996 
997  DenseMatrix<double> strain_rate_minus_one(3, 3);
998  strain_rate(1, s, strain_rate_minus_one);
999 
1000  // Get timestepper from first node
1001  TimeStepper* time_stepper_pt = node_pt(0)->time_stepper_pt();
1002 
1003  // Current and previous timesteps
1004  double dt_current = time_stepper_pt->time_pt()->dt(0);
1005  double dt_prev = time_stepper_pt->time_pt()->dt(1);
1006 
1007  // Extrapolate
1008  for (unsigned i = 0; i < 3; i++)
1009  {
1010  for (unsigned j = 0; j < 3; j++)
1011  {
1012  // Rate of changed based on previous two solutions
1013  double slope =
1014  (strain_rate_minus_one(i, j) - strain_rate_minus_two(i, j)) / dt_prev;
1015 
1016  // Extrapolated value from previous computed one to current one
1017  strainrate(i, j) = strain_rate_minus_one(i, j) + slope * dt_current;
1018  }
1019  }
1020  }
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:485
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:491
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
double & dt(const unsigned &t=0)
Definition: timesteppers.h:136

References oomph::Time::dt(), i, j, oomph::DenseMatrix< T >::ncol(), oomph::FiniteElement::node_pt(), oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, strain_rate(), oomph::TimeStepper::time_pt(), oomph::GeomObject::time_stepper_pt(), and oomph::Data::time_stepper_pt().

Referenced by extrapolated_strain_rate(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), and traction().

◆ fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::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 
)
inlinevirtual

Add the element's contribution to its residuals vector, jacobian matrix and mass matrix

Reimplemented from oomph::GeneralisedElement.

936  {
937  // Call the generic routine with the flag set to 2
939  parameter_pt, dres_dparam, djac_dparam, dmass_matrix_dparam, 2);
940  }
virtual 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)
Definition: generalised_newtonian_axisym_navier_stokes_elements.cc:3369

References fill_in_generic_dresidual_contribution_axi_nst().

◆ fill_in_contribution_to_djacobian_dparameter()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_contribution_to_djacobian_dparameter ( double *const &  parameter_pt,
Vector< double > &  dres_dparam,
DenseMatrix< double > &  djac_dparam 
)
inlinevirtual

Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by hanging-node version

Reimplemented from oomph::GeneralisedElement.

919  {
920  // Call the generic routine with the flag set to 1
922  parameter_pt,
923  dres_dparam,
924  djac_dparam,
926  1);
927  }
static DenseMatrix< double > Dummy_matrix
Definition: elements.h:227

References oomph::GeneralisedElement::Dummy_matrix, and fill_in_generic_dresidual_contribution_axi_nst().

◆ fill_in_contribution_to_dresiduals_dparameter()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_contribution_to_dresiduals_dparameter ( double *const &  parameter_pt,
Vector< double > &  dres_dparam 
)
inlinevirtual

Compute the element's residual Vector.

Reimplemented from oomph::GeneralisedElement.

902  {
903  // Call the generic residuals function with flag set to 0
904  // and using a dummy matrix argument
906  parameter_pt,
907  dres_dparam,
910  0);
911  }

References oomph::GeneralisedElement::Dummy_matrix, and fill_in_generic_dresidual_contribution_axi_nst().

◆ fill_in_contribution_to_hessian_vector_products()

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

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

Reimplemented from oomph::GeneralisedElement.

Reimplemented in oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations.

4102  {
4103  // Find out how many nodes there are
4104  unsigned n_node = nnode();
4105 
4106  // Get the nodal indices at which the velocity is stored
4107  unsigned u_nodal_index[3];
4108  for (unsigned i = 0; i < 3; ++i)
4109  {
4110  u_nodal_index[i] = u_index_axi_nst(i);
4111  }
4112 
4113  // Set up memory for the shape and test functions
4114  // Note that there are only two dimensions, r and z in this problem
4115  Shape psif(n_node), testf(n_node);
4116  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
4117 
4118  // Number of integration points
4119  unsigned Nintpt = integral_pt()->nweight();
4120 
4121  // Set the Vector to hold local coordinates (two dimensions)
4122  Vector<double> s(2);
4123 
4124  // Get Physical Variables from Element
4125  // Reynolds number must be multiplied by the density ratio
4126  double scaled_re = re() * density_ratio();
4127  // double visc_ratio = viscosity_ratio();
4128  Vector<double> G = g();
4129 
4130  // Integers used to store the local equation and unknown numbers
4131  int local_eqn = 0, local_unknown = 0, local_freedom = 0;
4132 
4133  // How may dofs are there
4134  const unsigned n_dof = this->ndof();
4135 
4136  // Create a local matrix eigenvector product contribution
4137  DenseMatrix<double> jac_y(n_dof, n_dof, 0.0);
4138 
4139  // Loop over the integration points
4140  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
4141  {
4142  // Assign values of s
4143  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
4144  // Get the integral weight
4145  double w = integral_pt()->weight(ipt);
4146 
4147  // Call the derivatives of the shape and test functions
4149  ipt, psif, dpsifdx, testf, dtestfdx);
4150 
4151  // Premultiply the weights and the Jacobian
4152  double W = w * J;
4153 
4154  // Allocate storage for the position and the derivative of the
4155  // mesh positions wrt time
4156  Vector<double> interpolated_x(2, 0.0);
4157  // Vector<double> mesh_velocity(2,0.0);
4158  // Allocate storage for the pressure, velocity components and their
4159  // time and spatial derivatives
4160  Vector<double> interpolated_u(3, 0.0);
4161  // Vector<double> dudt(3,0.0);
4162  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
4163 
4164  // Calculate velocities and derivatives at integration point
4165 
4166  // Loop over nodes
4167  for (unsigned l = 0; l < n_node; l++)
4168  {
4169  // Cache the shape function
4170  const double psif_ = psif(l);
4171  // Loop over the two coordinate directions
4172  for (unsigned i = 0; i < 2; i++)
4173  {
4174  interpolated_x[i] += this->raw_nodal_position(l, i) * psif_;
4175  }
4176 
4177  // Loop over the three velocity directions
4178  for (unsigned i = 0; i < 3; i++)
4179  {
4180  // Get the u_value
4181  const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
4182  interpolated_u[i] += u_value * psif_;
4183  // dudt[i]+= du_dt_axi_nst(l,i)*psif_;
4184  // Loop over derivative directions
4185  for (unsigned j = 0; j < 2; j++)
4186  {
4187  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
4188  }
4189  }
4190  }
4191 
4192  // Get the mesh velocity if ALE is enabled
4193  if (!ALE_is_disabled)
4194  {
4195  throw OomphLibError("Moving nodes not implemented\n",
4198  }
4199 
4200  // r is the first position component
4201  double r = interpolated_x[0];
4202 
4203 
4204  // MOMENTUM EQUATIONS
4205  //------------------
4206 
4207  // Loop over the test functions
4208  for (unsigned l = 0; l < n_node; l++)
4209  {
4210  // FIRST (RADIAL) MOMENTUM EQUATION
4211  local_eqn = nodal_local_eqn(l, u_nodal_index[0]);
4212  // If it's not a boundary condition
4213  if (local_eqn >= 0)
4214  {
4215  // Loop over the velocity shape functions yet again
4216  for (unsigned l3 = 0; l3 < n_node; l3++)
4217  {
4218  // Derivative of jacobian terms with respect to radial velocity
4219  local_freedom = nodal_local_eqn(l3, u_nodal_index[0]);
4220  if (local_freedom >= 0)
4221  {
4222  // Storage for the sums
4223  double temp = 0.0;
4224 
4225  // Loop over the velocity shape functions again
4226  for (unsigned l2 = 0; l2 < n_node; l2++)
4227  {
4228  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
4229  // Radial velocity component
4230  if (local_unknown >= 0)
4231  {
4232  // Remains of convective terms
4233  temp -= scaled_re *
4234  (r * psif[l2] * dpsifdx(l3, 0) +
4235  r * psif[l3] * dpsifdx(l2, 0)) *
4236  Y[local_unknown] * testf[l] * W;
4237  }
4238 
4239  // Axial velocity component
4240  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
4241  if (local_unknown >= 0)
4242  {
4243  // Convective terms
4244  temp -= scaled_re * r * psif[l2] * dpsifdx(l3, 1) *
4245  Y[local_unknown] * testf[l] * W;
4246  }
4247  }
4248  // Add the summed contribution to the product matrix
4249  jac_y(local_eqn, local_freedom) += temp;
4250  } // End of derivative wrt radial coordinate
4251 
4252 
4253  // Derivative of jacobian terms with respect to axial velocity
4254  local_freedom = nodal_local_eqn(l3, u_nodal_index[1]);
4255  if (local_freedom >= 0)
4256  {
4257  double temp = 0.0;
4258 
4259  // Loop over the velocity shape functions again
4260  for (unsigned l2 = 0; l2 < n_node; l2++)
4261  {
4262  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
4263  // Radial velocity component
4264  if (local_unknown >= 0)
4265  {
4266  // Remains of convective terms
4267  temp -= scaled_re * (r * psif[l3] * dpsifdx(l2, 1)) *
4268  Y[local_unknown] * testf[l] * W;
4269  }
4270  }
4271  // Add the summed contribution to the product matrix
4272  jac_y(local_eqn, local_freedom) += temp;
4273  } // End of derivative wrt axial coordinate
4274 
4275  // Derivative of jacobian terms with respect to azimuthal velocity
4276  local_freedom = nodal_local_eqn(l3, u_nodal_index[2]);
4277  if (local_freedom >= 0)
4278  {
4279  double temp = 0.0;
4280 
4281  // Loop over the velocity shape functions again
4282  for (unsigned l2 = 0; l2 < n_node; l2++)
4283  {
4284  // Azimuthal velocity component
4285  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
4286  if (local_unknown >= 0)
4287  {
4288  // Convective terms
4289  temp -= -scaled_re * 2.0 * psif[l3] * psif[l2] *
4290  Y[local_unknown] * testf[l] * W;
4291  }
4292  }
4293  // Add the summed contibution
4294  jac_y(local_eqn, local_freedom) += temp;
4295 
4296  } // End of if not boundary condition statement
4297  } // End of loop over freedoms
4298  } // End of RADIAL MOMENTUM EQUATION
4299 
4300 
4301  // SECOND (AXIAL) MOMENTUM EQUATION
4302  local_eqn = nodal_local_eqn(l, u_nodal_index[1]);
4303  // If it's not a boundary condition
4304  if (local_eqn >= 0)
4305  {
4306  // Loop over the velocity shape functions yet again
4307  for (unsigned l3 = 0; l3 < n_node; l3++)
4308  {
4309  // Derivative of jacobian terms with respect to radial velocity
4310  local_freedom = nodal_local_eqn(l3, u_nodal_index[0]);
4311  if (local_freedom >= 0)
4312  {
4313  double temp = 0.0;
4314 
4315  // Loop over the velocity shape functions again
4316  for (unsigned l2 = 0; l2 < n_node; l2++)
4317  {
4318  // Axial velocity component
4319  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
4320  if (local_unknown >= 0)
4321  {
4322  // Convective terms
4323  temp -= scaled_re * (r * psif[l3] * dpsifdx(l2, 0)) *
4324  Y[local_unknown] * testf[l] * W;
4325  }
4326  }
4327  jac_y(local_eqn, local_freedom) += temp;
4328 
4329  // There are no azimithal terms in the axial momentum equation
4330  } // End of loop over velocity shape functions
4331 
4332 
4333  // Derivative of jacobian terms with respect to axial velocity
4334  local_freedom = nodal_local_eqn(l3, u_nodal_index[1]);
4335  if (local_freedom >= 0)
4336  {
4337  double temp = 0.0;
4338 
4339  // Loop over the velocity shape functions again
4340  for (unsigned l2 = 0; l2 < n_node; l2++)
4341  {
4342  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
4343  // Radial velocity component
4344  if (local_unknown >= 0)
4345  {
4346  // Convective terms
4347  temp -= scaled_re * r * psif[l2] * dpsifdx(l3, 0) *
4348  Y[local_unknown] * testf[l] * W;
4349  }
4350 
4351  // Axial velocity component
4352  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
4353  if (local_unknown >= 0)
4354  {
4355  // Convective terms
4356  temp -= scaled_re *
4357  (r * psif[l2] * dpsifdx(l3, 1) +
4358  r * psif[l3] * dpsifdx(l2, 1)) *
4359  Y[local_unknown] * testf[l] * W;
4360  }
4361 
4362  // There are no azimithal terms in the axial momentum equation
4363  } // End of loop over velocity shape functions
4364 
4365  // Add summed contributiont to jacobian product matrix
4366  jac_y(local_eqn, local_freedom) += temp;
4367  }
4368  } // End of loop over local freedoms
4369 
4370  } // End of AXIAL MOMENTUM EQUATION
4371 
4372  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
4373  local_eqn = nodal_local_eqn(l, u_nodal_index[2]);
4374  if (local_eqn >= 0)
4375  {
4376  // Loop over the velocity shape functions yet again
4377  for (unsigned l3 = 0; l3 < n_node; l3++)
4378  {
4379  // Derivative of jacobian terms with respect to radial velocity
4380  local_freedom = nodal_local_eqn(l3, u_nodal_index[0]);
4381  if (local_freedom >= 0)
4382  {
4383  double temp = 0.0;
4384 
4385  // Loop over the velocity shape functions again
4386  for (unsigned l2 = 0; l2 < n_node; l2++)
4387  {
4388  // Azimuthal velocity component
4389  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
4390  if (local_unknown >= 0)
4391  {
4392  // Convective terms
4393  temp -=
4394  scaled_re *
4395  (r * psif[l3] * dpsifdx(l2, 0) + psif[l3] * psif[l2]) *
4396  Y[local_unknown] * testf[l] * W;
4397  }
4398  }
4399  // Add the summed contribution to the jacobian eigenvector sum
4400  jac_y(local_eqn, local_freedom) += temp;
4401  }
4402 
4403  // Derivative of jacobian terms with respect to axial velocity
4404  local_freedom = nodal_local_eqn(l3, u_nodal_index[1]);
4405  if (local_freedom >= 0)
4406  {
4407  double temp = 0.0;
4408 
4409  // Loop over the velocity shape functions again
4410  for (unsigned l2 = 0; l2 < n_node; l2++)
4411  {
4412  // Azimuthal velocity component
4413  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
4414  if (local_unknown >= 0)
4415  {
4416  // Convective terms
4417  temp -= scaled_re * (r * psif[l3] * dpsifdx(l2, 1)) *
4418  Y[local_unknown] * testf[l] * W;
4419  }
4420  }
4421  // Add the summed contribution to the jacobian eigenvector sum
4422  jac_y(local_eqn, local_freedom) += temp;
4423  }
4424 
4425 
4426  // Derivative of jacobian terms with respect to azimuthal velocity
4427  local_freedom = nodal_local_eqn(l3, u_nodal_index[2]);
4428  if (local_freedom >= 0)
4429  {
4430  double temp = 0.0;
4431 
4432 
4433  // Loop over the velocity shape functions again
4434  for (unsigned l2 = 0; l2 < n_node; l2++)
4435  {
4436  // Radial velocity component
4437  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
4438  if (local_unknown >= 0)
4439  {
4440  // Convective terms
4441  temp -=
4442  scaled_re *
4443  (r * psif[l2] * dpsifdx(l3, 0) + psif[l2] * psif[l3]) *
4444  Y[local_unknown] * testf[l] * W;
4445  }
4446 
4447  // Axial velocity component
4448  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
4449  if (local_unknown >= 0)
4450  {
4451  // Convective terms
4452  temp -= scaled_re * r * psif[l2] * dpsifdx(l3, 1) *
4453  Y[local_unknown] * testf[l] * W;
4454  }
4455  }
4456  // Add the fully summed contribution
4457  jac_y(local_eqn, local_freedom) += temp;
4458  }
4459  } // End of loop over freedoms
4460 
4461  // There are no pressure terms
4462  } // End of AZIMUTHAL EQUATION
4463 
4464  } // End of loop over shape functions
4465  }
4466 
4467  // We have now assembled the matrix (d J_{ij} Y_j)/d u_{k}
4468  // and simply need to sum over the vectors
4469  const unsigned n_vec = C.nrow();
4470  for (unsigned i = 0; i < n_dof; i++)
4471  {
4472  for (unsigned k = 0; k < n_dof; k++)
4473  {
4474  // Cache the value of the hessian y product
4475  const double j_y = jac_y(i, k);
4476  // Loop over the possible vectors
4477  for (unsigned v = 0; v < n_vec; v++)
4478  {
4479  product(v, i) += j_y * C(v, k);
4480  }
4481  }
4482  }
4483  }
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
JacobiRotation< float > G
Definition: Jacobi_makeGivens.cpp:2
Matrix< Scalar, Dynamic, Dynamic > C
Definition: bench_gemm.cpp:49
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2576
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.cc:1686
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
const double & re() const
Reynolds number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:420
const double & density_ratio() const
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:481
virtual double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
const Vector< double > & g() const
Vector of gravitational components.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:468
Definition: matrices.h:74
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
char char char int int * k
Definition: level2_impl.h:374
r
Definition: UniformPSDSelfTest.py:20
void product(const MatrixType &m)
Definition: product.h:42
const char Y
Definition: test/EulerAngles.cpp:32

References ALE_is_disabled, density_ratio(), dshape_and_dtest_eulerian_at_knot_axi_nst(), G, g(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), J, j, k, oomph::Integral::knot(), oomph::GeneralisedElement::ndof(), GlobalParameters::Nintpt, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::Matrix< T, MATRIX_TYPE >::nrow(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, product(), UniformPSDSelfTest::r, oomph::FiniteElement::raw_nodal_position(), oomph::FiniteElement::raw_nodal_value(), re(), s, u_index_axi_nst(), v, w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and Y.

◆ fill_in_contribution_to_jacobian()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_contribution_to_jacobian ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian 
)
inlinevirtual

Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by hanging-node version

Reimplemented from oomph::GeneralisedElement.

Reimplemented in oomph::PseudoSolidNodeUpdateElement< GeneralisedNewtonianAxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >.

869  {
870  // Call the generic routine with the flag set to 1
871 
872  // obacht
873  // Specific element routine is commented out and instead the
874  // default FD version is used
876  residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
877  // FiniteElement::fill_in_contribution_to_jacobian(residuals,jacobian);
878  }
virtual void fill_in_generic_residual_contribution_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Definition: generalised_newtonian_axisym_navier_stokes_elements.cc:1167

References oomph::GeneralisedElement::Dummy_matrix, and fill_in_generic_residual_contribution_axi_nst().

◆ fill_in_contribution_to_jacobian_and_mass_matrix()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_contribution_to_jacobian_and_mass_matrix ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
DenseMatrix< double > &  mass_matrix 
)
inlinevirtual

Add the element's contribution to its residuals vector, jacobian matrix and mass matrix

Reimplemented from oomph::GeneralisedElement.

Reimplemented in oomph::PseudoSolidNodeUpdateElement< GeneralisedNewtonianAxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >.

886  {
887  // Call the generic routine with the flag set to 2
889  residuals, jacobian, mass_matrix, 2);
890  }

References fill_in_generic_residual_contribution_axi_nst().

◆ fill_in_contribution_to_residuals()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_contribution_to_residuals ( Vector< double > &  residuals)
inlinevirtual

Compute the element's residual Vector.

Reimplemented from oomph::GeneralisedElement.

Reimplemented in oomph::PseudoSolidNodeUpdateElement< GeneralisedNewtonianAxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >.

855  {
856  // Call the generic residuals function with flag set to 0
857  // and using a dummy matrix argument
859  residuals,
862  0);
863  }

References oomph::GeneralisedElement::Dummy_matrix, and fill_in_generic_residual_contribution_axi_nst().

◆ fill_in_generic_dresidual_contribution_axi_nst()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::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 
)
protectedvirtual

Compute the derivative of residuals for the Navier–Stokes equations; with respect to a parameeter flag=2 or 1 or 0: compute the Jacobian and/or mass matrix as well

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

Reimplemented in oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations.

3375  {
3376  // Die if the parameter is not the Reynolds number
3377  if (parameter_pt != this->re_pt())
3378  {
3379  std::ostringstream error_stream;
3380  error_stream
3381  << "Cannot compute analytic jacobian for parameter addressed by "
3382  << parameter_pt << "\n";
3383  error_stream << "Can only compute derivatives wrt Re (" << Re_pt << ")\n";
3384  throw OomphLibError(
3385  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3386  }
3387 
3388  // Which parameters are we differentiating with respect to
3389  bool diff_re = false;
3390  bool diff_re_st = false;
3391  bool diff_re_inv_fr = false;
3392  bool diff_re_inv_ro = false;
3393 
3394  // Set the boolean flags according to the parameter pointer
3395  if (parameter_pt == this->re_pt())
3396  {
3397  diff_re = true;
3398  }
3399  if (parameter_pt == this->re_st_pt())
3400  {
3401  diff_re_st = true;
3402  }
3403  if (parameter_pt == this->re_invfr_pt())
3404  {
3405  diff_re_inv_fr = true;
3406  }
3407  if (parameter_pt == this->re_invro_pt())
3408  {
3409  diff_re_inv_ro = true;
3410  }
3411 
3412 
3413  // Find out how many nodes there are
3414  unsigned n_node = nnode();
3415 
3416  // Find out how many pressure dofs there are
3417  unsigned n_pres = npres_axi_nst();
3418 
3419  // Get the nodal indices at which the velocity is stored
3420  unsigned u_nodal_index[3];
3421  for (unsigned i = 0; i < 3; ++i)
3422  {
3423  u_nodal_index[i] = u_index_axi_nst(i);
3424  }
3425 
3426  // Set up memory for the shape and test functions
3427  // Note that there are only two dimensions, r and z in this problem
3428  Shape psif(n_node), testf(n_node);
3429  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
3430 
3431  // Set up memory for pressure shape and test functions
3432  Shape psip(n_pres), testp(n_pres);
3433 
3434  // Number of integration points
3435  unsigned Nintpt = integral_pt()->nweight();
3436 
3437  // Set the Vector to hold local coordinates (two dimensions)
3438  Vector<double> s(2);
3439 
3440  // Get Physical Variables from Element
3441  // Reynolds number must be multiplied by the density ratio
3442  // double scaled_re = re()*density_ratio();
3443  // double scaled_re_st = re_st()*density_ratio();
3444  // double scaled_re_inv_fr = re_invfr()*density_ratio();
3445  // double scaled_re_inv_ro = re_invro()*density_ratio();
3446  double dens_ratio = this->density_ratio();
3447  // double visc_ratio = viscosity_ratio();
3448  Vector<double> G = g();
3449 
3450  // Integers used to store the local equation and unknown numbers
3451  int local_eqn = 0, local_unknown = 0;
3452 
3453  // Loop over the integration points
3454  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
3455  {
3456  // Assign values of s
3457  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
3458  // Get the integral weight
3459  double w = integral_pt()->weight(ipt);
3460 
3461  // Call the derivatives of the shape and test functions
3463  ipt, psif, dpsifdx, testf, dtestfdx);
3464 
3465  // Call the pressure shape and test functions
3466  pshape_axi_nst(s, psip, testp);
3467 
3468  // Premultiply the weights and the Jacobian
3469  double W = w * J;
3470 
3471  // Allocate storage for the position and the derivative of the
3472  // mesh positions wrt time
3473  Vector<double> interpolated_x(2, 0.0);
3474  Vector<double> mesh_velocity(2, 0.0);
3475  // Allocate storage for the pressure, velocity components and their
3476  // time and spatial derivatives
3477  double interpolated_p = 0.0;
3478  Vector<double> interpolated_u(3, 0.0);
3479  Vector<double> dudt(3, 0.0);
3480  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
3481 
3482  // Calculate pressure at integration point
3483  for (unsigned l = 0; l < n_pres; l++)
3484  {
3485  interpolated_p += p_axi_nst(l) * psip[l];
3486  }
3487 
3488  // Calculate velocities and derivatives at integration point
3489 
3490  // Loop over nodes
3491  for (unsigned l = 0; l < n_node; l++)
3492  {
3493  // Cache the shape function
3494  const double psif_ = psif(l);
3495  // Loop over the two coordinate directions
3496  for (unsigned i = 0; i < 2; i++)
3497  {
3498  interpolated_x[i] += this->raw_nodal_position(l, i) * psif_;
3499  }
3500  // mesh_velocity[i] += dnodal_position_dt(l,i)*psif[l];
3501 
3502  // Loop over the three velocity directions
3503  for (unsigned i = 0; i < 3; i++)
3504  {
3505  // Get the u_value
3506  const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
3507  interpolated_u[i] += u_value * psif_;
3508  dudt[i] += du_dt_axi_nst(l, i) * psif_;
3509  // Loop over derivative directions
3510  for (unsigned j = 0; j < 2; j++)
3511  {
3512  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
3513  }
3514  }
3515  }
3516 
3517  // Get the mesh velocity if ALE is enabled
3518  if (!ALE_is_disabled)
3519  {
3520  // Loop over nodes
3521  for (unsigned l = 0; l < n_node; l++)
3522  {
3523  // Loop over the two coordinate directions
3524  for (unsigned i = 0; i < 2; i++)
3525  {
3526  mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif(l);
3527  }
3528  }
3529  }
3530 
3531 
3532  // Get the user-defined body force terms
3533  // Vector<double> body_force(3);
3534  // get_body_force(time(),ipt,interpolated_x,body_force);
3535 
3536  // Get the user-defined source function
3537  // double source = get_source_fct(time(),ipt,interpolated_x);
3538 
3539  // Get the user-defined viscosity function
3540  double visc_ratio;
3541  get_viscosity_ratio_axisym_nst(ipt, s, interpolated_x, visc_ratio);
3542 
3543  // r is the first position component
3544  double r = interpolated_x[0];
3545 
3546 
3547  // MOMENTUM EQUATIONS
3548  //------------------
3549 
3550  // Loop over the test functions
3551  for (unsigned l = 0; l < n_node; l++)
3552  {
3553  // FIRST (RADIAL) MOMENTUM EQUATION
3554  local_eqn = nodal_local_eqn(l, u_nodal_index[0]);
3555  // If it's not a boundary condition
3556  if (local_eqn >= 0)
3557  {
3558  // No user-defined body force terms
3559  // dres_dparam[local_eqn] +=
3560  // r*body_force[0]*testf[l]*W;
3561 
3562  // Add the gravitational body force term if the reynolds number
3563  // is equal to re_inv_fr
3564  if (diff_re_inv_fr)
3565  {
3566  dres_dparam[local_eqn] += r * dens_ratio * testf[l] * G[0] * W;
3567  }
3568 
3569  // No pressure gradient term
3570  // residuals[local_eqn] +=
3571  // interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
3572 
3573  // No in the stress tensor terms
3574  // The viscosity ratio needs to go in here to ensure
3575  // continuity of normal stress is satisfied even in flows
3576  // with zero pressure gradient!
3577  // residuals[local_eqn] -= visc_ratio*
3578  // r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
3579 
3580  // residuals[local_eqn] -= visc_ratio*r*
3581  // (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
3582  // dtestfdx(l,1)*W;
3583 
3584  // residuals[local_eqn] -=
3585  // visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
3586 
3587  // Add in the inertial terms
3588  // du/dt term
3589  if (diff_re_st)
3590  {
3591  dres_dparam[local_eqn] -= dens_ratio * r * dudt[0] * testf[l] * W;
3592  }
3593 
3594  // Convective terms
3595  if (diff_re)
3596  {
3597  dres_dparam[local_eqn] -=
3598  dens_ratio *
3599  (r * interpolated_u[0] * interpolated_dudx(0, 0) -
3600  interpolated_u[2] * interpolated_u[2] +
3601  r * interpolated_u[1] * interpolated_dudx(0, 1)) *
3602  testf[l] * W;
3603  }
3604 
3605  // Mesh velocity terms
3606  if (!ALE_is_disabled)
3607  {
3608  if (diff_re_st)
3609  {
3610  for (unsigned k = 0; k < 2; k++)
3611  {
3612  dres_dparam[local_eqn] += dens_ratio * r * mesh_velocity[k] *
3613  interpolated_dudx(0, k) * testf[l] *
3614  W;
3615  }
3616  }
3617  }
3618 
3619  // Add the Coriolis term
3620  if (diff_re_inv_ro)
3621  {
3622  dres_dparam[local_eqn] +=
3623  2.0 * r * dens_ratio * interpolated_u[2] * testf[l] * W;
3624  }
3625 
3626  // CALCULATE THE JACOBIAN
3627  if (flag)
3628  {
3629  // Loop over the velocity shape functions again
3630  for (unsigned l2 = 0; l2 < n_node; l2++)
3631  {
3632  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
3633  // Radial velocity component
3634  if (local_unknown >= 0)
3635  {
3636  if (flag == 2)
3637  {
3638  if (diff_re_st)
3639  {
3640  // Add the mass matrix
3641  dmass_matrix_dparam(local_eqn, local_unknown) +=
3642  dens_ratio * r * psif[l2] * testf[l] * W;
3643  }
3644  }
3645 
3646  // Add contribution to the Jacobian matrix
3647  // jacobian(local_eqn,local_unknown)
3648  // -= visc_ratio*r*(1.0+Gamma[0])
3649  //*dpsifdx(l2,0)*dtestfdx(l,0)*W;
3650 
3651  // jacobian(local_eqn,local_unknown)
3652  // -= visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
3653 
3654  // jacobian(local_eqn,local_unknown)
3655  // -= visc_ratio*(1.0 + Gamma[0])*psif[l2]*testf[l]*W/r;
3656 
3657  // Add in the inertial terms
3658  // du/dt term
3659  if (diff_re_st)
3660  {
3661  djac_dparam(local_eqn, local_unknown) -=
3662  dens_ratio * r *
3663  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
3664  testf[l] * W;
3665  }
3666 
3667  // Convective terms
3668  if (diff_re)
3669  {
3670  djac_dparam(local_eqn, local_unknown) -=
3671  dens_ratio *
3672  (r * psif[l2] * interpolated_dudx(0, 0) +
3673  r * interpolated_u[0] * dpsifdx(l2, 0) +
3674  r * interpolated_u[1] * dpsifdx(l2, 1)) *
3675  testf[l] * W;
3676  }
3677 
3678  // Mesh velocity terms
3679  if (!ALE_is_disabled)
3680  {
3681  for (unsigned k = 0; k < 2; k++)
3682  {
3683  if (diff_re_st)
3684  {
3685  djac_dparam(local_eqn, local_unknown) +=
3686  dens_ratio * r * mesh_velocity[k] * dpsifdx(l2, k) *
3687  testf[l] * W;
3688  }
3689  }
3690  }
3691  }
3692 
3693  // Axial velocity component
3694  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
3695  if (local_unknown >= 0)
3696  {
3697  // jacobian(local_eqn,local_unknown) -=
3698  // visc_ratio*r*Gamma[0]*dpsifdx(l2,0)*dtestfdx(l,1)*W;
3699 
3700  // Convective terms
3701  if (diff_re)
3702  {
3703  djac_dparam(local_eqn, local_unknown) -=
3704  dens_ratio * r * psif[l2] * interpolated_dudx(0, 1) *
3705  testf[l] * W;
3706  }
3707  }
3708 
3709  // Azimuthal velocity component
3710  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
3711  if (local_unknown >= 0)
3712  {
3713  // Convective terms
3714  if (diff_re)
3715  {
3716  djac_dparam(local_eqn, local_unknown) -=
3717  -dens_ratio * 2.0 * interpolated_u[2] * psif[l2] *
3718  testf[l] * W;
3719  }
3720  // Coriolis terms
3721  if (diff_re_inv_ro)
3722  {
3723  djac_dparam(local_eqn, local_unknown) +=
3724  2.0 * r * dens_ratio * psif[l2] * testf[l] * W;
3725  }
3726  }
3727  }
3728 
3729  /*Now loop over pressure shape functions*/
3730  /*This is the contribution from pressure gradient*/
3731  // for(unsigned l2=0;l2<n_pres;l2++)
3732  // {
3733  // local_unknown = p_local_eqn(l2);
3734  // /*If we are at a non-zero degree of freedom in the entry*/
3735  // if(local_unknown >= 0)
3736  // {
3737  // jacobian(local_eqn,local_unknown)
3738  // += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
3739  // }
3740  // }
3741  } /*End of Jacobian calculation*/
3742 
3743  } // End of if not boundary condition statement
3744 
3745  // SECOND (AXIAL) MOMENTUM EQUATION
3746  local_eqn = nodal_local_eqn(l, u_nodal_index[1]);
3747  // If it's not a boundary condition
3748  if (local_eqn >= 0)
3749  {
3750  // Add the user-defined body force terms
3751  // Remember to multiply by the density ratio!
3752  // residuals[local_eqn] +=
3753  // r*body_force[1]*testf[l]*W;
3754 
3755  // Add the gravitational body force term
3756  if (diff_re_inv_fr)
3757  {
3758  dres_dparam[local_eqn] += r * dens_ratio * testf[l] * G[1] * W;
3759  }
3760 
3761  // Add the pressure gradient term
3762  // residuals[local_eqn] += r*interpolated_p*dtestfdx(l,1)*W;
3763 
3764  // Add in the stress tensor terms
3765  // The viscosity ratio needs to go in here to ensure
3766  // continuity of normal stress is satisfied even in flows
3767  // with zero pressure gradient!
3768  // residuals[local_eqn] -= visc_ratio*
3769  // r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
3770  // *dtestfdx(l,0)*W;
3771 
3772  // residuals[local_eqn] -= visc_ratio*r*
3773  // (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
3774 
3775  // Add in the inertial terms
3776  // du/dt term
3777  if (diff_re_st)
3778  {
3779  dres_dparam[local_eqn] -= dens_ratio * r * dudt[1] * testf[l] * W;
3780  }
3781 
3782  // Convective terms
3783  if (diff_re)
3784  {
3785  dres_dparam[local_eqn] -=
3786  dens_ratio *
3787  (r * interpolated_u[0] * interpolated_dudx(1, 0) +
3788  r * interpolated_u[1] * interpolated_dudx(1, 1)) *
3789  testf[l] * W;
3790  }
3791 
3792  // Mesh velocity terms
3793  if (!ALE_is_disabled)
3794  {
3795  if (diff_re_st)
3796  {
3797  for (unsigned k = 0; k < 2; k++)
3798  {
3799  dres_dparam[local_eqn] += dens_ratio * r * mesh_velocity[k] *
3800  interpolated_dudx(1, k) * testf[l] *
3801  W;
3802  }
3803  }
3804  }
3805 
3806  // CALCULATE THE JACOBIAN
3807  if (flag)
3808  {
3809  // Loop over the velocity shape functions again
3810  for (unsigned l2 = 0; l2 < n_node; l2++)
3811  {
3812  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
3813  // Radial velocity component
3814  if (local_unknown >= 0)
3815  {
3816  // Add in the stress tensor terms
3817  // The viscosity ratio needs to go in here to ensure
3818  // continuity of normal stress is satisfied even in flows
3819  // with zero pressure gradient!
3820  // jacobian(local_eqn,local_unknown) -=
3821  // visc_ratio*r*Gamma[1]*dpsifdx(l2,1)*dtestfdx(l,0)*W;
3822 
3823  // Convective terms
3824  if (diff_re)
3825  {
3826  djac_dparam(local_eqn, local_unknown) -=
3827  dens_ratio * r * psif[l2] * interpolated_dudx(1, 0) *
3828  testf[l] * W;
3829  }
3830  }
3831 
3832  // Axial velocity component
3833  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
3834  if (local_unknown >= 0)
3835  {
3836  if (flag == 2)
3837  {
3838  if (diff_re_st)
3839  {
3840  // Add the mass matrix
3841  dmass_matrix_dparam(local_eqn, local_unknown) +=
3842  dens_ratio * r * psif[l2] * testf[l] * W;
3843  }
3844  }
3845 
3846 
3847  // Add in the stress tensor terms
3848  // The viscosity ratio needs to go in here to ensure
3849  // continuity of normal stress is satisfied even in flows
3850  // with zero pressure gradient!
3851  // jacobian(local_eqn,local_unknown) -=
3852  // visc_ratio*r*dpsifdx(l2,0)*dtestfdx(l,0)*W;
3853 
3854  // jacobian(local_eqn,local_unknown) -=
3855  // visc_ratio*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
3856  // dtestfdx(l,1)*W;
3857 
3858  // Add in the inertial terms
3859  // du/dt term
3860  if (diff_re_st)
3861  {
3862  djac_dparam(local_eqn, local_unknown) -=
3863  dens_ratio * r *
3864  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
3865  testf[l] * W;
3866  }
3867 
3868  // Convective terms
3869  if (diff_re)
3870  {
3871  djac_dparam(local_eqn, local_unknown) -=
3872  dens_ratio *
3873  (r * interpolated_u[0] * dpsifdx(l2, 0) +
3874  r * psif[l2] * interpolated_dudx(1, 1) +
3875  r * interpolated_u[1] * dpsifdx(l2, 1)) *
3876  testf[l] * W;
3877  }
3878 
3879  // Mesh velocity terms
3880  if (!ALE_is_disabled)
3881  {
3882  for (unsigned k = 0; k < 2; k++)
3883  {
3884  if (diff_re_st)
3885  {
3886  djac_dparam(local_eqn, local_unknown) +=
3887  dens_ratio * r * mesh_velocity[k] * dpsifdx(l2, k) *
3888  testf[l] * W;
3889  }
3890  }
3891  }
3892  }
3893 
3894  // There are no azimithal terms in the axial momentum equation
3895  } // End of loop over velocity shape functions
3896 
3897  } /*End of Jacobian calculation*/
3898 
3899  } // End of AXIAL MOMENTUM EQUATION
3900 
3901  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
3902  local_eqn = nodal_local_eqn(l, u_nodal_index[2]);
3903  if (local_eqn >= 0)
3904  {
3905  // Add the user-defined body force terms
3906  // Remember to multiply by the density ratio!
3907  // residuals[local_eqn] +=
3908  // r*body_force[2]*testf[l]*W;
3909 
3910  // Add the gravitational body force term
3911  if (diff_re_inv_fr)
3912  {
3913  dres_dparam[local_eqn] += r * dens_ratio * testf[l] * G[2] * W;
3914  }
3915 
3916  // There is NO pressure gradient term
3917 
3918  // Add in the stress tensor terms
3919  // The viscosity ratio needs to go in here to ensure
3920  // continuity of normal stress is satisfied even in flows
3921  // with zero pressure gradient!
3922  // residuals[local_eqn] -= visc_ratio*
3923  // (r*interpolated_dudx(2,0) -
3924  // Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
3925 
3926  // residuals[local_eqn] -= visc_ratio*r*
3927  // interpolated_dudx(2,1)*dtestfdx(l,1)*W;
3928 
3929  // residuals[local_eqn] -= visc_ratio*
3930  // ((interpolated_u[2]/r) -
3931  // Gamma[0]*interpolated_dudx(2,0))*testf[l]*W;
3932 
3933 
3934  // Add in the inertial terms
3935  // du/dt term
3936  if (diff_re_st)
3937  {
3938  dres_dparam[local_eqn] -= dens_ratio * r * dudt[2] * testf[l] * W;
3939  }
3940 
3941  // Convective terms
3942  if (diff_re)
3943  {
3944  dres_dparam[local_eqn] -=
3945  dens_ratio *
3946  (r * interpolated_u[0] * interpolated_dudx(2, 0) +
3947  interpolated_u[0] * interpolated_u[2] +
3948  r * interpolated_u[1] * interpolated_dudx(2, 1)) *
3949  testf[l] * W;
3950  }
3951 
3952  // Mesh velocity terms
3953  if (!ALE_is_disabled)
3954  {
3955  if (diff_re_st)
3956  {
3957  for (unsigned k = 0; k < 2; k++)
3958  {
3959  dres_dparam[local_eqn] += dens_ratio * r * mesh_velocity[k] *
3960  interpolated_dudx(2, k) * testf[l] *
3961  W;
3962  }
3963  }
3964  }
3965 
3966  // Add the Coriolis term
3967  if (diff_re_inv_ro)
3968  {
3969  dres_dparam[local_eqn] -=
3970  2.0 * r * dens_ratio * interpolated_u[0] * testf[l] * W;
3971  }
3972 
3973  // CALCULATE THE JACOBIAN
3974  if (flag)
3975  {
3976  // Loop over the velocity shape functions again
3977  for (unsigned l2 = 0; l2 < n_node; l2++)
3978  {
3979  // Radial velocity component
3980  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
3981  if (local_unknown >= 0)
3982  {
3983  // Convective terms
3984  if (diff_re)
3985  {
3986  djac_dparam(local_eqn, local_unknown) -=
3987  dens_ratio *
3988  (r * psif[l2] * interpolated_dudx(2, 0) +
3989  psif[l2] * interpolated_u[2]) *
3990  testf[l] * W;
3991  }
3992 
3993  // Coriolis term
3994  if (diff_re_inv_ro)
3995  {
3996  djac_dparam(local_eqn, local_unknown) -=
3997  2.0 * r * dens_ratio * psif[l2] * testf[l] * W;
3998  }
3999  }
4000 
4001  // Axial velocity component
4002  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
4003  if (local_unknown >= 0)
4004  {
4005  // Convective terms
4006  if (diff_re)
4007  {
4008  djac_dparam(local_eqn, local_unknown) -=
4009  dens_ratio * r * psif[l2] * interpolated_dudx(2, 1) *
4010  testf[l] * W;
4011  }
4012  }
4013 
4014  // Azimuthal velocity component
4015  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
4016  if (local_unknown >= 0)
4017  {
4018  if (flag == 2)
4019  {
4020  // Add the mass matrix
4021  if (diff_re_st)
4022  {
4023  dmass_matrix_dparam(local_eqn, local_unknown) +=
4024  dens_ratio * r * psif[l2] * testf[l] * W;
4025  }
4026  }
4027 
4028  // Add in the stress tensor terms
4029  // The viscosity ratio needs to go in here to ensure
4030  // continuity of normal stress is satisfied even in flows
4031  // with zero pressure gradient!
4032  // jacobian(local_eqn,local_unknown) -=
4033  // visc_ratio*(r*dpsifdx(l2,0) -
4034  // Gamma[0]*psif[l2])*dtestfdx(l,0)*W;
4035 
4036  // jacobian(local_eqn,local_unknown) -=
4037  // visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
4038 
4039  // jacobian(local_eqn,local_unknown) -=
4040  // visc_ratio*((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
4041  // *testf[l]*W;
4042 
4043  // Add in the inertial terms
4044  // du/dt term
4045  if (diff_re_st)
4046  {
4047  djac_dparam(local_eqn, local_unknown) -=
4048  dens_ratio * r *
4049  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
4050  testf[l] * W;
4051  }
4052 
4053  // Convective terms
4054  if (diff_re)
4055  {
4056  djac_dparam(local_eqn, local_unknown) -=
4057  dens_ratio *
4058  (r * interpolated_u[0] * dpsifdx(l2, 0) +
4059  interpolated_u[0] * psif[l2] +
4060  r * interpolated_u[1] * dpsifdx(l2, 1)) *
4061  testf[l] * W;
4062  }
4063 
4064  // Mesh velocity terms
4065  if (!ALE_is_disabled)
4066  {
4067  if (diff_re_st)
4068  {
4069  for (unsigned k = 0; k < 2; k++)
4070  {
4071  djac_dparam(local_eqn, local_unknown) +=
4072  dens_ratio * r * mesh_velocity[k] * dpsifdx(l2, k) *
4073  testf[l] * W;
4074  }
4075  }
4076  }
4077  }
4078  }
4079 
4080  // There are no pressure terms
4081  } // End of Jacobian
4082 
4083  } // End of AZIMUTHAL EQUATION
4084 
4085  } // End of loop over shape functions
4086 
4087 
4088  // CONTINUITY EQUATION NO PARAMETERS
4089  //-------------------
4090  }
4091  }
double raw_dnodal_position_dt(const unsigned &n, const unsigned &i) const
Definition: elements.h:2256
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:438
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
double du_dt_axi_nst(const unsigned &n, const unsigned &i) const
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:570
virtual void get_viscosity_ratio_axisym_nst(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, double &visc_ratio)
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:357
double *& re_invro_pt()
Pointer to global inverse Froude number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:462
virtual double p_axi_nst(const unsigned &n_p) const =0
double *& re_invfr_pt()
Pointer to global inverse Froude number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:450
double *& re_pt()
Pointer to Reynolds number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:432

References ALE_is_disabled, density_ratio(), dshape_and_dtest_eulerian_at_knot_axi_nst(), du_dt_axi_nst(), G, g(), get_viscosity_ratio_axisym_nst(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), J, j, k, oomph::Integral::knot(), GlobalParameters::Nintpt, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::node_pt(), npres_axi_nst(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p_axi_nst(), pshape_axi_nst(), UniformPSDSelfTest::r, oomph::FiniteElement::raw_dnodal_position_dt(), oomph::FiniteElement::raw_nodal_position(), oomph::FiniteElement::raw_nodal_value(), re_invfr_pt(), re_invro_pt(), Re_pt, re_pt(), re_st_pt(), s, oomph::Data::time_stepper_pt(), u_index_axi_nst(), w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and oomph::TimeStepper::weight().

Referenced by fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter(), fill_in_contribution_to_djacobian_dparameter(), and fill_in_contribution_to_dresiduals_dparameter().

◆ fill_in_generic_residual_contribution_axi_nst()

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

Compute the residuals for the Navier–Stokes equations; flag=2 or 1 or 0: compute the Jacobian and/or mass matrix as well.

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

Reimplemented in oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations.

1172  {
1173  // Return immediately if there are no dofs
1174  if (ndof() == 0) return;
1175 
1176  // Find out how many nodes there are
1177  unsigned n_node = nnode();
1178 
1179  // Get continuous time from timestepper of first node
1180  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1181 
1182  // Find out how many pressure dofs there are
1183  unsigned n_pres = npres_axi_nst();
1184 
1185  // Get the nodal indices at which the velocity is stored
1186  unsigned u_nodal_index[3];
1187  for (unsigned i = 0; i < 3; ++i)
1188  {
1189  u_nodal_index[i] = u_index_axi_nst(i);
1190  }
1191 
1192  // Set up memory for the shape and test functions
1193  // Note that there are only two dimensions, r and z in this problem
1194  Shape psif(n_node), testf(n_node);
1195  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
1196 
1197  // Set up memory for pressure shape and test functions
1198  Shape psip(n_pres), testp(n_pres);
1199 
1200  // Number of integration points
1201  unsigned Nintpt = integral_pt()->nweight();
1202 
1203  // Set the Vector to hold local coordinates (two dimensions)
1204  Vector<double> s(2);
1205 
1206  // Get Physical Variables from Element
1207  // Reynolds number must be multiplied by the density ratio
1208  double scaled_re = re() * density_ratio();
1209  double scaled_re_st = re_st() * density_ratio();
1210  double scaled_re_inv_fr = re_invfr() * density_ratio();
1211  double scaled_re_inv_ro = re_invro() * density_ratio();
1212  // double visc_ratio = viscosity_ratio();
1213  Vector<double> G = g();
1214 
1215  // Integers used to store the local equation and unknown numbers
1216  int local_eqn = 0, local_unknown = 0;
1217 
1218  // Loop over the integration points
1219  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
1220  {
1221  // Assign values of s
1222  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
1223  // Get the integral weight
1224  double w = integral_pt()->weight(ipt);
1225 
1226  // Call the derivatives of the shape and test functions
1228  ipt, psif, dpsifdx, testf, dtestfdx);
1229 
1230  // Call the pressure shape and test functions
1231  pshape_axi_nst(s, psip, testp);
1232 
1233  // Premultiply the weights and the Jacobian
1234  double W = w * J;
1235 
1236  // Allocate storage for the position and the derivative of the
1237  // mesh positions wrt time
1238  Vector<double> interpolated_x(2, 0.0);
1239  Vector<double> mesh_velocity(2, 0.0);
1240  // Allocate storage for the pressure, velocity components and their
1241  // time and spatial derivatives
1242  double interpolated_p = 0.0;
1243  Vector<double> interpolated_u(3, 0.0);
1244  Vector<double> dudt(3, 0.0);
1245  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
1246 
1247  // Calculate pressure at integration point
1248  for (unsigned l = 0; l < n_pres; l++)
1249  {
1250  interpolated_p += p_axi_nst(l) * psip[l];
1251  }
1252 
1253  // Calculate velocities and derivatives at integration point
1254 
1255  // Loop over nodes
1256  for (unsigned l = 0; l < n_node; l++)
1257  {
1258  // Cache the shape function
1259  const double psif_ = psif(l);
1260  // Loop over the two coordinate directions
1261  for (unsigned i = 0; i < 2; i++)
1262  {
1263  interpolated_x[i] += this->raw_nodal_position(l, i) * psif_;
1264  }
1265 
1266  // Loop over the three velocity directions
1267  for (unsigned i = 0; i < 3; i++)
1268  {
1269  // Get the u_value
1270  const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
1271  interpolated_u[i] += u_value * psif_;
1272  dudt[i] += du_dt_axi_nst(l, i) * psif_;
1273  // Loop over derivative directions
1274  for (unsigned j = 0; j < 2; j++)
1275  {
1276  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1277  }
1278  }
1279  }
1280 
1281  // Get the mesh velocity if ALE is enabled
1282  if (!ALE_is_disabled)
1283  {
1284  // Loop over nodes
1285  for (unsigned l = 0; l < n_node; l++)
1286  {
1287  // Loop over the two coordinate directions
1288  for (unsigned i = 0; i < 2; i++)
1289  {
1290  mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif(l);
1291  }
1292  }
1293  }
1294 
1295  // The strainrate used to compute the second invariant
1296  DenseMatrix<double> strainrate_to_compute_second_invariant(3, 3, 0.0);
1297 
1298  // the strainrate used to calculate the second invariant
1299  // can be either the current one or the one extrapolated from
1300  // previous velocity values
1302  {
1303  strain_rate(s, strainrate_to_compute_second_invariant);
1304  }
1305  else
1306  {
1307  extrapolated_strain_rate(ipt, strainrate_to_compute_second_invariant);
1308  }
1309 
1310  // Calculate the second invariant
1312  strainrate_to_compute_second_invariant);
1313 
1314  // Get the viscosity according to the constitutive equation
1315  double viscosity = Constitutive_eqn_pt->viscosity(second_invariant);
1316 
1317  // Get the user-defined body force terms
1318  Vector<double> body_force(3);
1320 
1321  // Get the user-defined source function
1322  double source = get_source_fct(time, ipt, interpolated_x);
1323 
1324  // Get the user-defined viscosity function
1325  double visc_ratio;
1326  get_viscosity_ratio_axisym_nst(ipt, s, interpolated_x, visc_ratio);
1327 
1328  // r is the first position component
1329  double r = interpolated_x[0];
1330 
1331  // obacht set up vectors of the viscosity differentiated w.r.t.
1332  // the velocity components (radial, axial, azimuthal)
1333  Vector<double> dviscosity_dUr(n_node, 0.0);
1334  Vector<double> dviscosity_dUz(n_node, 0.0);
1335  Vector<double> dviscosity_dUphi(n_node, 0.0);
1336 
1338  {
1339  // Calculate the derivate of the viscosity w.r.t. the second invariant
1340  double dviscosity_dsecond_invariant =
1342 
1343  // FD step
1344  // double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
1345 
1346  // calculate a reference strainrate
1347  DenseMatrix<double> strainrate_ref(3, 3, 0.0);
1348  strain_rate(s, strainrate_ref);
1349 
1350  // pre-compute the derivative of the second invariant w.r.t. the
1351  // entries in the rate of strain tensor
1352  DenseMatrix<double> dinvariant_dstrainrate(3, 3, 0.0);
1353 
1354  // d I_2 / d epsilon_{r,r}
1355  dinvariant_dstrainrate(0, 0) =
1356  strainrate_ref(1, 1) + strainrate_ref(2, 2);
1357  // d I_2 / d epsilon_{z,z}
1358  dinvariant_dstrainrate(1, 1) =
1359  strainrate_ref(0, 0) + strainrate_ref(2, 2);
1360  // d I_2 / d epsilon_{phi,phi}
1361  dinvariant_dstrainrate(2, 2) =
1362  strainrate_ref(0, 0) + strainrate_ref(1, 1);
1363  // d I_2 / d epsilon_{r,z}
1364  dinvariant_dstrainrate(0, 1) = -strainrate_ref(1, 0);
1365  // d I_2 / d epsilon_{z,r}
1366  dinvariant_dstrainrate(1, 0) = -strainrate_ref(0, 1);
1367  // d I_2 / d epsilon_{r,phi}
1368  dinvariant_dstrainrate(0, 2) = -strainrate_ref(2, 0);
1369  // d I_2 / d epsilon_{phi,r}
1370  dinvariant_dstrainrate(2, 0) = -strainrate_ref(0, 2);
1371  // d I_2 / d epsilon_{phi,z}
1372  dinvariant_dstrainrate(2, 1) = -strainrate_ref(1, 2);
1373  // d I_2 / d epsilon_{z,phi}
1374  dinvariant_dstrainrate(1, 2) = -strainrate_ref(2, 1);
1375 
1376  // loop over the nodes
1377  for (unsigned l = 0; l < n_node; l++)
1378  {
1379  // Get pointer to l-th local node
1380  // Node* nod_pt = node_pt(l);
1381 
1382  // loop over the three velocity components
1383  for (unsigned i = 0; i < 3; i++)
1384  {
1385  // back up
1386  // double backup = nod_pt->value(u_nodal_index[i]);
1387 
1388  // do the FD step
1389  // nod_pt->set_value(u_nodal_index[i],
1390  // nod_pt->value(u_nodal_index[i])+eps_fd);
1391 
1392  // calculate updated strainrate
1393  // DenseMatrix<double> strainrate_fd(3,3,0.0);
1394  // strain_rate(s,strainrate_fd);
1395 
1396  // initialise the derivative of the second invariant w.r.t. the
1397  // unknown velocity U_{i,l}
1398  double dinvariant_dunknown = 0.0;
1399 
1400  // loop over the entries of the rate of strain tensor
1401  for (unsigned m = 0; m < 3; m++)
1402  {
1403  for (unsigned n = 0; n < 3; n++)
1404  {
1405  // initialise the derivative of the strainrate w.r.t. the
1406  // unknown velocity U_{i,l}
1407  double dstrainrate_dunknown = 0.0;
1408 
1409  // switch based on first index
1410  switch (m)
1411  {
1412  // epsilon_{r ...}
1413  case 0:
1414 
1415  // switch for second index
1416  switch (n)
1417  {
1418  // epsilon_{r r}
1419  case 0:
1420  if (i == 0)
1421  {
1422  dstrainrate_dunknown = dpsifdx(l, 0);
1423  }
1424  break;
1425 
1426  // epsilon_{r z}
1427  case 1:
1428  if (i == 0)
1429  {
1430  dstrainrate_dunknown = 0.5 * dpsifdx(l, 1);
1431  }
1432  else if (i == 1)
1433  {
1434  dstrainrate_dunknown = 0.5 * dpsifdx(l, 0);
1435  }
1436  break;
1437 
1438  // epsilon_{r phi}
1439  case 2:
1440  if (i == 2)
1441  {
1442  dstrainrate_dunknown =
1443  0.5 * dpsifdx(l, 0) - 0.5 / r * psif[l];
1444  }
1445  break;
1446 
1447  default:
1448  std::ostringstream error_stream;
1449  error_stream << "Should never get here...";
1450  throw OomphLibError(error_stream.str(),
1453  }
1454 
1455  break;
1456 
1457  // epsilon_{z ...}
1458  case 1:
1459 
1460  // switch for second index
1461  switch (n)
1462  {
1463  // epsilon_{z r}
1464  case 0:
1465  if (i == 0)
1466  {
1467  dstrainrate_dunknown = 0.5 * dpsifdx(l, 1);
1468  }
1469  else if (i == 1)
1470  {
1471  dstrainrate_dunknown = 0.5 * dpsifdx(l, 0);
1472  }
1473  break;
1474 
1475  // epsilon_{z z}
1476  case 1:
1477  if (i == 1)
1478  {
1479  dstrainrate_dunknown = dpsifdx(l, 1);
1480  }
1481  else
1482  {
1483  // dstrainrate_dunknown=0.0;
1484  }
1485  break;
1486 
1487  // epsilon_{z phi}
1488  case 2:
1489  if (i == 2)
1490  {
1491  dstrainrate_dunknown = 0.5 * dpsifdx(l, 1);
1492  }
1493  break;
1494 
1495  default:
1496  std::ostringstream error_stream;
1497  error_stream << "Should never get here...";
1498  throw OomphLibError(error_stream.str(),
1501  }
1502 
1503  break;
1504 
1505  // epsilon_{phi ...}
1506  case 2:
1507 
1508  // switch for second index
1509  switch (n)
1510  {
1511  // epsilon_{phi r}
1512  case 0:
1513  if (i == 2)
1514  {
1515  dstrainrate_dunknown =
1516  0.5 * dpsifdx(l, 0) - 0.5 / r * psif[l];
1517  }
1518  break;
1519 
1520  // epsilon_{phi z}
1521  case 1:
1522  if (i == 2)
1523  {
1524  dstrainrate_dunknown = 0.5 * dpsifdx(l, 1);
1525  }
1526  break;
1527 
1528  // epsilon_{phi phi}
1529  case 2:
1530  if (i == 0)
1531  {
1532  dstrainrate_dunknown = 1.0 / r * psif[l];
1533  }
1534  break;
1535 
1536  default:
1537  std::ostringstream error_stream;
1538  error_stream << "Should never get here...";
1539  throw OomphLibError(error_stream.str(),
1542  }
1543 
1544  break;
1545 
1546  default:
1547  std::ostringstream error_stream;
1548  error_stream << "Should never get here...";
1549  throw OomphLibError(error_stream.str(),
1552  }
1553  // calculate the difference
1554  // double dstrainrate_dunknown =
1555  // (strainrate_fd(m,n)-strainrate_ref(m,n))/eps_fd;
1556 
1557  dinvariant_dunknown +=
1558  dinvariant_dstrainrate(m, n) * dstrainrate_dunknown;
1559  }
1560  }
1561 
1562  // // get the invariant of the updated strainrate
1563  // double second_invariant_fd=
1564  // SecondInvariantHelper::second_invariant(strainrate_fd);
1565 
1566  // // calculate the difference
1567  // double dinvariant_dunknown =
1568  // (second_invariant_fd - second_invariant)/eps_fd;
1569 
1570  switch (i)
1571  {
1572  case 0:
1573  dviscosity_dUr[l] =
1574  dviscosity_dsecond_invariant * dinvariant_dunknown;
1575  break;
1576 
1577  case 1:
1578  dviscosity_dUz[l] =
1579  dviscosity_dsecond_invariant * dinvariant_dunknown;
1580  break;
1581 
1582  case 2:
1583  dviscosity_dUphi[l] =
1584  dviscosity_dsecond_invariant * dinvariant_dunknown;
1585  break;
1586 
1587  default:
1588  std::ostringstream error_stream;
1589  error_stream << "Should never get here...";
1590  throw OomphLibError(error_stream.str(),
1593  }
1594 
1595  // Reset
1596  // nod_pt->set_value(u_nodal_index[i],backup);
1597  }
1598  }
1599  }
1600 
1601  // MOMENTUM EQUATIONS
1602  //------------------
1603 
1604  // Loop over the test functions
1605  for (unsigned l = 0; l < n_node; l++)
1606  {
1607  // FIRST (RADIAL) MOMENTUM EQUATION
1608  local_eqn = nodal_local_eqn(l, u_nodal_index[0]);
1609  // If it's not a boundary condition
1610  if (local_eqn >= 0)
1611  {
1612  // Add the user-defined body force terms
1613  residuals[local_eqn] += r * body_force[0] * testf[l] * W;
1614 
1615  // Add the gravitational body force term
1616  residuals[local_eqn] += r * scaled_re_inv_fr * testf[l] * G[0] * W;
1617 
1618  // Add the pressure gradient term
1619  // Potentially pre-multiply by viscosity ratio
1621  {
1622  residuals[local_eqn] += visc_ratio * viscosity * interpolated_p *
1623  (testf[l] + r * dtestfdx(l, 0)) * W;
1624  }
1625  else
1626  {
1627  residuals[local_eqn] +=
1628  interpolated_p * (testf[l] + r * dtestfdx(l, 0)) * W;
1629  }
1630 
1631  // Add in the stress tensor terms
1632  // The viscosity ratio needs to go in here to ensure
1633  // continuity of normal stress is satisfied even in flows
1634  // with zero pressure gradient!
1635  // stress term 1
1636  residuals[local_eqn] -= visc_ratio * viscosity * r *
1637  (1.0 + Gamma[0]) * interpolated_dudx(0, 0) *
1638  dtestfdx(l, 0) * W;
1639 
1640  // stress term 2
1641  residuals[local_eqn] -=
1642  visc_ratio * viscosity * r *
1643  (interpolated_dudx(0, 1) + Gamma[0] * interpolated_dudx(1, 0)) *
1644  dtestfdx(l, 1) * W;
1645 
1646  // stress term 3
1647  residuals[local_eqn] -= visc_ratio * viscosity * (1.0 + Gamma[0]) *
1648  interpolated_u[0] * testf[l] * W / r;
1649 
1650  // Add in the inertial terms
1651  // du/dt term
1652  residuals[local_eqn] -= scaled_re_st * r * dudt[0] * testf[l] * W;
1653 
1654  // Convective terms
1655  residuals[local_eqn] -=
1656  scaled_re *
1657  (r * interpolated_u[0] * interpolated_dudx(0, 0) -
1658  interpolated_u[2] * interpolated_u[2] +
1659  r * interpolated_u[1] * interpolated_dudx(0, 1)) *
1660  testf[l] * W;
1661 
1662  // Mesh velocity terms
1663  if (!ALE_is_disabled)
1664  {
1665  for (unsigned k = 0; k < 2; k++)
1666  {
1667  residuals[local_eqn] += scaled_re_st * r * mesh_velocity[k] *
1668  interpolated_dudx(0, k) * testf[l] * W;
1669  }
1670  }
1671 
1672  // Add the Coriolis term
1673  residuals[local_eqn] +=
1674  2.0 * r * scaled_re_inv_ro * interpolated_u[2] * testf[l] * W;
1675 
1676  // CALCULATE THE JACOBIAN
1677  if (flag)
1678  {
1679  // Loop over the velocity shape functions again
1680  for (unsigned l2 = 0; l2 < n_node; l2++)
1681  {
1682  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
1683  // Radial velocity component
1684  if (local_unknown >= 0)
1685  {
1686  if (flag == 2)
1687  {
1688  // Add the mass matrix
1689  mass_matrix(local_eqn, local_unknown) +=
1690  scaled_re_st * r * psif[l2] * testf[l] * W;
1691  }
1692 
1693  // stress term 1
1694  jacobian(local_eqn, local_unknown) -=
1695  visc_ratio * viscosity * r * (1.0 + Gamma[0]) *
1696  dpsifdx(l2, 0) * dtestfdx(l, 0) * W;
1697 
1699  {
1700  // stress term 1 contribution from generalised Newtonian model
1701  jacobian(local_eqn, local_unknown) -=
1702  visc_ratio * dviscosity_dUr[l2] * r * (1.0 + Gamma[0]) *
1703  interpolated_dudx(0, 0) * dtestfdx(l, 0) * W;
1704  }
1705 
1706  // stress term 2
1707  jacobian(local_eqn, local_unknown) -= visc_ratio * viscosity *
1708  r * dpsifdx(l2, 1) *
1709  dtestfdx(l, 1) * W;
1710 
1712  {
1713  // stress term 2 contribution from generalised Newtonian model
1714  jacobian(local_eqn, local_unknown) -=
1715  visc_ratio * dviscosity_dUr[l2] * r *
1716  (interpolated_dudx(0, 1) +
1717  Gamma[0] * interpolated_dudx(1, 0)) *
1718  dtestfdx(l, 1) * W;
1719  }
1720 
1721  // stress term 3
1722  jacobian(local_eqn, local_unknown) -=
1723  visc_ratio * viscosity * (1.0 + Gamma[0]) * psif[l2] *
1724  testf[l] * W / r;
1725 
1727  {
1728  // stress term 3 contribution from generalised Newtonian model
1729  jacobian(local_eqn, local_unknown) -=
1730  visc_ratio * dviscosity_dUr[l2] * (1.0 + Gamma[0]) *
1731  interpolated_u[0] * testf[l] * W / r;
1732  }
1733 
1734  // Add in the inertial terms
1735  // du/dt term
1736  jacobian(local_eqn, local_unknown) -=
1737  scaled_re_st * r *
1738  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
1739  testf[l] * W;
1740 
1741  // Convective terms
1742  jacobian(local_eqn, local_unknown) -=
1743  scaled_re *
1744  (r * psif[l2] * interpolated_dudx(0, 0) +
1745  r * interpolated_u[0] * dpsifdx(l2, 0) +
1746  r * interpolated_u[1] * dpsifdx(l2, 1)) *
1747  testf[l] * W;
1748 
1749  // Mesh velocity terms
1750  if (!ALE_is_disabled)
1751  {
1752  for (unsigned k = 0; k < 2; k++)
1753  {
1754  jacobian(local_eqn, local_unknown) +=
1755  scaled_re_st * r * mesh_velocity[k] * dpsifdx(l2, k) *
1756  testf[l] * W;
1757  }
1758  }
1759  }
1760 
1761 
1762  // Axial velocity component
1763  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
1764  if (local_unknown >= 0)
1765  {
1766  // no stress term 1
1767 
1769  {
1770  // stress term 1 contribution from generalised Newtonian model
1771  jacobian(local_eqn, local_unknown) -=
1772  visc_ratio * dviscosity_dUz[l2] * r * (1.0 + Gamma[0]) *
1773  interpolated_dudx(0, 0) * dtestfdx(l, 0) * W;
1774  }
1775 
1776  // stress term 2
1777  jacobian(local_eqn, local_unknown) -=
1778  visc_ratio * viscosity * r * Gamma[0] * dpsifdx(l2, 0) *
1779  dtestfdx(l, 1) * W;
1780 
1782  {
1783  // stress term 2 contribution from generalised Newtonian model
1784  jacobian(local_eqn, local_unknown) -=
1785  visc_ratio * dviscosity_dUz[l2] * r *
1786  (interpolated_dudx(0, 1) +
1787  Gamma[0] * interpolated_dudx(1, 0)) *
1788  dtestfdx(l, 1) * W;
1789  }
1790 
1791  // no stress term 3
1792 
1794  {
1795  // stress term 3 contribution from generalised Newtonian model
1796  jacobian(local_eqn, local_unknown) -=
1797  visc_ratio * dviscosity_dUz[l2] * (1.0 + Gamma[0]) *
1798  interpolated_u[0] * testf[l] * W / r;
1799  }
1800 
1801  // Convective terms
1802  jacobian(local_eqn, local_unknown) -= scaled_re * r * psif[l2] *
1803  interpolated_dudx(0, 1) *
1804  testf[l] * W;
1805  }
1806 
1807  // Azimuthal velocity component
1808  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
1809  if (local_unknown >= 0)
1810  {
1811  // no stress term 1
1812 
1814  {
1815  // stress term 1 contribution from generalised Newtonian model
1816  jacobian(local_eqn, local_unknown) -=
1817  visc_ratio * dviscosity_dUphi[l2] * r * (1.0 + Gamma[0]) *
1818  interpolated_dudx(0, 0) * dtestfdx(l, 0) * W;
1819  }
1820 
1821  // no stress term 2
1822 
1824  {
1825  // stress term 2 contribution from generalised Newtonian model
1826  jacobian(local_eqn, local_unknown) -=
1827  visc_ratio * dviscosity_dUphi[l2] * r *
1828  (interpolated_dudx(0, 1) +
1829  Gamma[0] * interpolated_dudx(1, 0)) *
1830  dtestfdx(l, 1) * W;
1831  }
1832 
1833  // no stress term 3
1834 
1836  {
1837  // stress term 3 contribution from generalised Newtonian model
1838  jacobian(local_eqn, local_unknown) -=
1839  visc_ratio * dviscosity_dUphi[l2] * (1.0 + Gamma[0]) *
1840  interpolated_u[0] * testf[l] * W / r;
1841  }
1842 
1843  // Convective terms
1844  jacobian(local_eqn, local_unknown) -= -scaled_re * 2.0 *
1845  interpolated_u[2] *
1846  psif[l2] * testf[l] * W;
1847 
1848  // Coriolis terms
1849  jacobian(local_eqn, local_unknown) +=
1850  2.0 * r * scaled_re_inv_ro * psif[l2] * testf[l] * W;
1851  }
1852  }
1853 
1854  /*Now loop over pressure shape functions*/
1855  /*This is the contribution from pressure gradient*/
1856  // Potentially pre-multiply by viscosity ratio
1857  for (unsigned l2 = 0; l2 < n_pres; l2++)
1858  {
1859  local_unknown = p_local_eqn(l2);
1860  /*If we are at a non-zero degree of freedom in the entry*/
1861  if (local_unknown >= 0)
1862  {
1864  {
1865  jacobian(local_eqn, local_unknown) +=
1866  visc_ratio * viscosity * psip[l2] *
1867  (testf[l] + r * dtestfdx(l, 0)) * W;
1868  }
1869  else
1870  {
1871  jacobian(local_eqn, local_unknown) +=
1872  psip[l2] * (testf[l] + r * dtestfdx(l, 0)) * W;
1873  }
1874  }
1875  }
1876  } /*End of Jacobian calculation*/
1877 
1878  } // End of if not boundary condition statement
1879 
1880  // SECOND (AXIAL) MOMENTUM EQUATION
1881  local_eqn = nodal_local_eqn(l, u_nodal_index[1]);
1882  // If it's not a boundary condition
1883  if (local_eqn >= 0)
1884  {
1885  // Add the user-defined body force terms
1886  // Remember to multiply by the density ratio!
1887  residuals[local_eqn] += r * body_force[1] * testf[l] * W;
1888 
1889  // Add the gravitational body force term
1890  residuals[local_eqn] += r * scaled_re_inv_fr * testf[l] * G[1] * W;
1891 
1892  // Add the pressure gradient term
1893  // Potentially pre-multiply by viscosity ratio
1895  {
1896  residuals[local_eqn] +=
1897  r * visc_ratio * viscosity * interpolated_p * dtestfdx(l, 1) * W;
1898  }
1899  else
1900  {
1901  residuals[local_eqn] += r * interpolated_p * dtestfdx(l, 1) * W;
1902  }
1903 
1904  // Add in the stress tensor terms
1905  // The viscosity ratio needs to go in here to ensure
1906  // continuity of normal stress is satisfied even in flows
1907  // with zero pressure gradient!
1908  // stress term 1
1909  residuals[local_eqn] -=
1910  visc_ratio * viscosity * r *
1911  (interpolated_dudx(1, 0) + Gamma[1] * interpolated_dudx(0, 1)) *
1912  dtestfdx(l, 0) * W;
1913 
1914  // stress term 2
1915  residuals[local_eqn] -= visc_ratio * viscosity * r *
1916  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
1917  dtestfdx(l, 1) * W;
1918 
1919  // Add in the inertial terms
1920  // du/dt term
1921  residuals[local_eqn] -= scaled_re_st * r * dudt[1] * testf[l] * W;
1922 
1923  // Convective terms
1924  residuals[local_eqn] -=
1925  scaled_re *
1926  (r * interpolated_u[0] * interpolated_dudx(1, 0) +
1927  r * interpolated_u[1] * interpolated_dudx(1, 1)) *
1928  testf[l] * W;
1929 
1930  // Mesh velocity terms
1931  if (!ALE_is_disabled)
1932  {
1933  for (unsigned k = 0; k < 2; k++)
1934  {
1935  residuals[local_eqn] += scaled_re_st * r * mesh_velocity[k] *
1936  interpolated_dudx(1, k) * testf[l] * W;
1937  }
1938  }
1939 
1940  // CALCULATE THE JACOBIAN
1941  if (flag)
1942  {
1943  // Loop over the velocity shape functions again
1944  for (unsigned l2 = 0; l2 < n_node; l2++)
1945  {
1946  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
1947  // Radial velocity component
1948  if (local_unknown >= 0)
1949  {
1950  // Add in the stress tensor terms
1951  // The viscosity ratio needs to go in here to ensure
1952  // continuity of normal stress is satisfied even in flows
1953  // with zero pressure gradient!
1954  // stress term 1
1955  jacobian(local_eqn, local_unknown) -=
1956  visc_ratio * viscosity * r * Gamma[1] * dpsifdx(l2, 1) *
1957  dtestfdx(l, 0) * W;
1958 
1960  {
1961  // stress term 1 contribution from generalised Newtonian model
1962  jacobian(local_eqn, local_unknown) -=
1963  visc_ratio * dviscosity_dUr[l2] * r *
1964  (interpolated_dudx(1, 0) +
1965  Gamma[1] * interpolated_dudx(0, 1)) *
1966  dtestfdx(l, 0) * W;
1967  }
1968 
1969  // no stress term 2
1970 
1972  {
1973  // stress term 2 contribution from generalised Newtonian model
1974  jacobian(local_eqn, local_unknown) -=
1975  visc_ratio * dviscosity_dUr[l2] * r * (1.0 + Gamma[1]) *
1976  interpolated_dudx(1, 1) * dtestfdx(l, 1) * W;
1977  }
1978 
1979  // Convective terms
1980  jacobian(local_eqn, local_unknown) -= scaled_re * r * psif[l2] *
1981  interpolated_dudx(1, 0) *
1982  testf[l] * W;
1983  }
1984 
1985  // Axial velocity component
1986  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
1987  if (local_unknown >= 0)
1988  {
1989  if (flag == 2)
1990  {
1991  // Add the mass matrix
1992  mass_matrix(local_eqn, local_unknown) +=
1993  scaled_re_st * r * psif[l2] * testf[l] * W;
1994  }
1995 
1996  // Add in the stress tensor terms
1997  // The viscosity ratio needs to go in here to ensure
1998  // continuity of normal stress is satisfied even in flows
1999  // with zero pressure gradient!
2000  // stress term 1
2001  jacobian(local_eqn, local_unknown) -= visc_ratio * viscosity *
2002  r * dpsifdx(l2, 0) *
2003  dtestfdx(l, 0) * W;
2004 
2006  {
2007  // stress term 1 contribution from generalised Newtonian model
2008  jacobian(local_eqn, local_unknown) -=
2009  visc_ratio * dviscosity_dUz[l2] * r *
2010  (interpolated_dudx(1, 0) +
2011  Gamma[1] * interpolated_dudx(0, 1)) *
2012  dtestfdx(l, 0) * W;
2013  }
2014 
2015  // stress term 2
2016  jacobian(local_eqn, local_unknown) -=
2017  visc_ratio * viscosity * r * (1.0 + Gamma[1]) *
2018  dpsifdx(l2, 1) * dtestfdx(l, 1) * W;
2019 
2021  {
2022  // stress term 2 contribution from generalised Newtonian model
2023  jacobian(local_eqn, local_unknown) -=
2024  visc_ratio * dviscosity_dUz[l2] * r * (1.0 + Gamma[1]) *
2025  interpolated_dudx(1, 1) * dtestfdx(l, 1) * W;
2026  }
2027 
2028  // Add in the inertial terms
2029  // du/dt term
2030  jacobian(local_eqn, local_unknown) -=
2031  scaled_re_st * r *
2032  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
2033  testf[l] * W;
2034 
2035  // Convective terms
2036  jacobian(local_eqn, local_unknown) -=
2037  scaled_re *
2038  (r * interpolated_u[0] * dpsifdx(l2, 0) +
2039  r * psif[l2] * interpolated_dudx(1, 1) +
2040  r * interpolated_u[1] * dpsifdx(l2, 1)) *
2041  testf[l] * W;
2042 
2043 
2044  // Mesh velocity terms
2045  if (!ALE_is_disabled)
2046  {
2047  for (unsigned k = 0; k < 2; k++)
2048  {
2049  jacobian(local_eqn, local_unknown) +=
2050  scaled_re_st * r * mesh_velocity[k] * dpsifdx(l2, k) *
2051  testf[l] * W;
2052  }
2053  }
2054  }
2055 
2056  // There are no azimithal terms in the axial momentum equation
2057  // edit: for the generalised Newtonian elements there are...
2058  // Azimuthal velocity component
2059  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
2060  if (local_unknown >= 0)
2061  {
2063  {
2064  // stress term 1 contribution from generalised Newtonian model
2065  jacobian(local_eqn, local_unknown) -=
2066  visc_ratio * dviscosity_dUphi[l2] * r *
2067  (interpolated_dudx(1, 0) +
2068  Gamma[1] * interpolated_dudx(0, 1)) *
2069  dtestfdx(l, 0) * W;
2070 
2071  // stress term 2 contribution from generalised Newtonian model
2072  jacobian(local_eqn, local_unknown) -=
2073  visc_ratio * dviscosity_dUphi[l2] * r * (1.0 + Gamma[1]) *
2074  interpolated_dudx(1, 1) * dtestfdx(l, 1) * W;
2075  }
2076  }
2077  } // End of loop over velocity shape functions
2078 
2079  /*Now loop over pressure shape functions*/
2080  /*This is the contribution from pressure gradient*/
2081  // Potentially pre-multiply by viscosity ratio
2082  for (unsigned l2 = 0; l2 < n_pres; l2++)
2083  {
2084  local_unknown = p_local_eqn(l2);
2085  /*If we are at a non-zero degree of freedom in the entry*/
2086  if (local_unknown >= 0)
2087  {
2089  {
2090  jacobian(local_eqn, local_unknown) +=
2091  r * visc_ratio * viscosity * psip[l2] * dtestfdx(l, 1) * W;
2092  }
2093  else
2094  {
2095  jacobian(local_eqn, local_unknown) +=
2096  r * psip[l2] * dtestfdx(l, 1) * W;
2097  }
2098  }
2099  }
2100  } /*End of Jacobian calculation*/
2101 
2102  } // End of AXIAL MOMENTUM EQUATION
2103 
2104  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
2105  local_eqn = nodal_local_eqn(l, u_nodal_index[2]);
2106  if (local_eqn >= 0)
2107  {
2108  // Add the user-defined body force terms
2109  // Remember to multiply by the density ratio!
2110  residuals[local_eqn] += r * body_force[2] * testf[l] * W;
2111 
2112  // Add the gravitational body force term
2113  residuals[local_eqn] += r * scaled_re_inv_fr * testf[l] * G[2] * W;
2114 
2115  // There is NO pressure gradient term
2116 
2117  // Add in the stress tensor terms
2118  // The viscosity ratio needs to go in here to ensure
2119  // continuity of normal stress is satisfied even in flows
2120  // with zero pressure gradient!
2121  // stress term 1
2122  residuals[local_eqn] -=
2123  visc_ratio * viscosity *
2124  (r * interpolated_dudx(2, 0) - Gamma[0] * interpolated_u[2]) *
2125  dtestfdx(l, 0) * W;
2126 
2127  // stress term 2
2128  residuals[local_eqn] -= visc_ratio * viscosity * r *
2129  interpolated_dudx(2, 1) * dtestfdx(l, 1) * W;
2130 
2131  // stress term 3
2132  residuals[local_eqn] -=
2133  visc_ratio * viscosity *
2134  ((interpolated_u[2] / r) - Gamma[0] * interpolated_dudx(2, 0)) *
2135  testf[l] * W;
2136 
2137 
2138  // Add in the inertial terms
2139  // du/dt term
2140  residuals[local_eqn] -= scaled_re_st * r * dudt[2] * testf[l] * W;
2141 
2142  // Convective terms
2143  residuals[local_eqn] -=
2144  scaled_re *
2145  (r * interpolated_u[0] * interpolated_dudx(2, 0) +
2146  interpolated_u[0] * interpolated_u[2] +
2147  r * interpolated_u[1] * interpolated_dudx(2, 1)) *
2148  testf[l] * W;
2149 
2150  // Mesh velocity terms
2151  if (!ALE_is_disabled)
2152  {
2153  for (unsigned k = 0; k < 2; k++)
2154  {
2155  residuals[local_eqn] += scaled_re_st * r * mesh_velocity[k] *
2156  interpolated_dudx(2, k) * testf[l] * W;
2157  }
2158  }
2159 
2160  // Add the Coriolis term
2161  residuals[local_eqn] -=
2162  2.0 * r * scaled_re_inv_ro * interpolated_u[0] * testf[l] * W;
2163 
2164  // CALCULATE THE JACOBIAN
2165  if (flag)
2166  {
2167  // Loop over the velocity shape functions again
2168  for (unsigned l2 = 0; l2 < n_node; l2++)
2169  {
2170  // Radial velocity component
2171  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
2172  if (local_unknown >= 0)
2173  {
2175  {
2176  // stress term 1 contribution from generalised Newtonian model
2177  jacobian(local_eqn, local_unknown) -=
2178  visc_ratio * dviscosity_dUr[l2] *
2179  (r * interpolated_dudx(2, 0) -
2180  Gamma[0] * interpolated_u[2]) *
2181  dtestfdx(l, 0) * W;
2182 
2183  // stress term 2 contribution from generalised Newtonian model
2184  jacobian(local_eqn, local_unknown) -=
2185  visc_ratio * dviscosity_dUr[l2] * r *
2186  interpolated_dudx(2, 1) * dtestfdx(l, 1) * W;
2187 
2188  // stress term 3 contribution from generalised Newtonian model
2189  jacobian(local_eqn, local_unknown) -=
2190  visc_ratio * dviscosity_dUr[l2] *
2191  ((interpolated_u[2] / r) -
2192  Gamma[0] * interpolated_dudx(2, 0)) *
2193  testf[l] * W;
2194  }
2195 
2196  // Convective terms
2197  jacobian(local_eqn, local_unknown) -=
2198  scaled_re *
2199  (r * psif[l2] * interpolated_dudx(2, 0) +
2200  psif[l2] * interpolated_u[2]) *
2201  testf[l] * W;
2202 
2203  // Coriolis term
2204  jacobian(local_eqn, local_unknown) -=
2205  2.0 * r * scaled_re_inv_ro * psif[l2] * testf[l] * W;
2206  }
2207 
2208  // Axial velocity component
2209  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
2210  if (local_unknown >= 0)
2211  {
2213  {
2214  // stress term 1 contribution from generalised Newtonian model
2215  jacobian(local_eqn, local_unknown) -=
2216  visc_ratio * dviscosity_dUz[l2] *
2217  (r * interpolated_dudx(2, 0) -
2218  Gamma[0] * interpolated_u[2]) *
2219  dtestfdx(l, 0) * W;
2220 
2221  // stress term 2 contribution from generalised Newtonian model
2222  jacobian(local_eqn, local_unknown) -=
2223  visc_ratio * dviscosity_dUz[l2] * r *
2224  interpolated_dudx(2, 1) * dtestfdx(l, 1) * W;
2225 
2226  // stress term 3 contribution from generalised Newtonian model
2227  jacobian(local_eqn, local_unknown) -=
2228  visc_ratio * dviscosity_dUz[l2] *
2229  ((interpolated_u[2] / r) -
2230  Gamma[0] * interpolated_dudx(2, 0)) *
2231  testf[l] * W;
2232  }
2233 
2234  // Convective terms
2235  jacobian(local_eqn, local_unknown) -= scaled_re * r * psif[l2] *
2236  interpolated_dudx(2, 1) *
2237  testf[l] * W;
2238  }
2239 
2240  // Azimuthal velocity component
2241  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
2242  if (local_unknown >= 0)
2243  {
2244  if (flag == 2)
2245  {
2246  // Add the mass matrix
2247  mass_matrix(local_eqn, local_unknown) +=
2248  scaled_re_st * r * psif[l2] * testf[l] * W;
2249  }
2250 
2251  // Add in 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  // stress term 1
2256  jacobian(local_eqn, local_unknown) -=
2257  visc_ratio * viscosity *
2258  (r * dpsifdx(l2, 0) - Gamma[0] * psif[l2]) * dtestfdx(l, 0) *
2259  W;
2260 
2262  {
2263  // stress term 1 contribution from generalised Newtonian model
2264  jacobian(local_eqn, local_unknown) -=
2265  visc_ratio * dviscosity_dUphi[l2] *
2266  (r * interpolated_dudx(2, 0) -
2267  Gamma[0] * interpolated_u[2]) *
2268  dtestfdx(l, 0) * W;
2269  }
2270 
2271  // stress term 2
2272  jacobian(local_eqn, local_unknown) -= visc_ratio * viscosity *
2273  r * dpsifdx(l2, 1) *
2274  dtestfdx(l, 1) * W;
2275 
2277  {
2278  // stress term 2 contribution from generalised Newtonian model
2279  jacobian(local_eqn, local_unknown) -=
2280  visc_ratio * dviscosity_dUphi[l2] * r *
2281  interpolated_dudx(2, 1) * dtestfdx(l, 1) * W;
2282  }
2283 
2284  // stress term 3
2285  jacobian(local_eqn, local_unknown) -=
2286  visc_ratio * viscosity *
2287  ((psif[l2] / r) - Gamma[0] * dpsifdx(l2, 0)) * testf[l] * W;
2288 
2290  {
2291  // stress term 3 contribution from generalised Newtonian model
2292  jacobian(local_eqn, local_unknown) -=
2293  visc_ratio * dviscosity_dUphi[l2] *
2294  ((interpolated_u[2] / r) -
2295  Gamma[0] * interpolated_dudx(2, 0)) *
2296  testf[l] * W;
2297  }
2298 
2299  // Add in the inertial terms
2300  // du/dt term
2301  jacobian(local_eqn, local_unknown) -=
2302  scaled_re_st * r *
2303  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
2304  testf[l] * W;
2305 
2306  // Convective terms
2307  jacobian(local_eqn, local_unknown) -=
2308  scaled_re *
2309  (r * interpolated_u[0] * dpsifdx(l2, 0) +
2310  interpolated_u[0] * psif[l2] +
2311  r * interpolated_u[1] * dpsifdx(l2, 1)) *
2312  testf[l] * W;
2313 
2314  // Mesh velocity terms
2315  if (!ALE_is_disabled)
2316  {
2317  for (unsigned k = 0; k < 2; k++)
2318  {
2319  jacobian(local_eqn, local_unknown) +=
2320  scaled_re_st * r * mesh_velocity[k] * dpsifdx(l2, k) *
2321  testf[l] * W;
2322  }
2323  }
2324  }
2325  }
2326 
2327  // There are no pressure terms
2328  } // End of Jacobian
2329 
2330  } // End of AZIMUTHAL EQUATION
2331 
2332  } // End of loop over shape functions
2333 
2334 
2335  // CONTINUITY EQUATION
2336  //-------------------
2337 
2338  // Loop over the shape functions
2339  for (unsigned l = 0; l < n_pres; l++)
2340  {
2341  local_eqn = p_local_eqn(l);
2342  // If not a boundary conditions
2343  if (local_eqn >= 0)
2344  {
2345  // Source term
2346  residuals[local_eqn] -= r * source * testp[l] * W;
2347 
2348  // Gradient terms
2349  // Potentially pre-multiply by viscosity ratio
2351  {
2352  residuals[local_eqn] +=
2353  visc_ratio * viscosity *
2354  (interpolated_u[0] + r * interpolated_dudx(0, 0) +
2355  r * interpolated_dudx(1, 1)) *
2356  testp[l] * W;
2357  }
2358  else
2359  {
2360  residuals[local_eqn] +=
2361  (interpolated_u[0] + r * interpolated_dudx(0, 0) +
2362  r * interpolated_dudx(1, 1)) *
2363  testp[l] * W;
2364  }
2365 
2366  /*CALCULATE THE JACOBIAN*/
2367  if (flag)
2368  {
2369  /*Loop over the velocity shape functions*/
2370  for (unsigned l2 = 0; l2 < n_node; l2++)
2371  {
2372  // Radial velocity component
2373  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
2374  if (local_unknown >= 0)
2375  {
2377  {
2378  jacobian(local_eqn, local_unknown) +=
2379  visc_ratio * viscosity * (psif[l2] + r * dpsifdx(l2, 0)) *
2380  testp[l] * W;
2381  }
2382  else
2383  {
2384  jacobian(local_eqn, local_unknown) +=
2385  (psif[l2] + r * dpsifdx(l2, 0)) * testp[l] * W;
2386  }
2387  }
2388 
2389  // Axial velocity component
2390  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
2391  if (local_unknown >= 0)
2392  {
2394  {
2395  jacobian(local_eqn, local_unknown) +=
2396  r * visc_ratio * viscosity * dpsifdx(l2, 1) * testp[l] * W;
2397  }
2398  else
2399  {
2400  jacobian(local_eqn, local_unknown) +=
2401  r * dpsifdx(l2, 1) * testp[l] * W;
2402  }
2403  }
2404  } /*End of loop over l2*/
2405  } /*End of Jacobian calculation*/
2406 
2407  } // End of if not boundary condition
2408 
2409  } // End of loop over l
2410  }
2411  }
static bool Pre_multiply_by_viscosity_ratio
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:139
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:456
virtual int p_local_eqn(const unsigned &n) const =0
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:426
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
const double & re_invfr() const
Global inverse Froude number.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:444
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
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
virtual double viscosity(const double &second_invariant_of_rate_of_strain_tensor)=0
virtual double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)=0
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
int * m
Definition: level2_cplx_impl.h:294
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Definition: axisym_linear_elasticity/cylinder/cylinder.cc:96
void source(const Vector< double > &x, Vector< double > &f)
Source function.
Definition: unstructured_two_d_circle.cc:46
double second_invariant(const DenseMatrix< double > &tensor)
Compute second invariant of tensor.
Definition: oomph_utilities.cc:163

References ALE_is_disabled, Global_Parameters::body_force(), Constitutive_eqn_pt, density_ratio(), dshape_and_dtest_eulerian_at_knot_axi_nst(), du_dt_axi_nst(), oomph::GeneralisedNewtonianConstitutiveEquation< DIM >::dviscosity_dinvariant(), extrapolated_strain_rate(), G, g(), Gamma, get_body_force_axi_nst(), get_source_fct(), get_viscosity_ratio_axisym_nst(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), J, j, k, oomph::Integral::knot(), m, n, oomph::GeneralisedElement::ndof(), GlobalParameters::Nintpt, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::node_pt(), npres_axi_nst(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p_axi_nst(), p_local_eqn(), Pre_multiply_by_viscosity_ratio, pshape_axi_nst(), UniformPSDSelfTest::r, oomph::FiniteElement::raw_dnodal_position_dt(), oomph::FiniteElement::raw_nodal_position(), oomph::FiniteElement::raw_nodal_value(), re(), re_invfr(), re_invro(), re_st(), s, oomph::SecondInvariantHelper::second_invariant(), TestProblem::source(), strain_rate(), oomph::Time::time(), oomph::TimeStepper::time_pt(), oomph::Data::time_stepper_pt(), u_index_axi_nst(), Use_extrapolated_strainrate_to_compute_second_invariant, oomph::GeneralisedNewtonianConstitutiveEquation< DIM >::viscosity(), w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and oomph::TimeStepper::weight().

Referenced by fill_in_contribution_to_jacobian(), fill_in_contribution_to_jacobian_and_mass_matrix(), and fill_in_contribution_to_residuals().

◆ g()

◆ g_pt()

Vector<double>*& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::g_pt ( )
inline

Pointer to Vector of gravitational components.

475  {
476  return G_pt;
477  }

References G_pt.

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

◆ get_body_force_axi_nst()

virtual void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_body_force_axi_nst ( const double time,
const unsigned ipt,
const Vector< double > &  s,
const Vector< double > &  x,
Vector< double > &  result 
)
inlineprotectedvirtual

Calculate the body force fct at a given time and Eulerian position.

242  {
243  // If the function pointer is zero return zero
244  if (Body_force_fct_pt == 0)
245  {
246  // Loop over dimensions and set body forces to zero
247  for (unsigned i = 0; i < 3; i++)
248  {
249  result[i] = 0.0;
250  }
251  }
252  // Otherwise call the function
253  else
254  {
255  (*Body_force_fct_pt)(time, x, result);
256  }
257  }

References Body_force_fct_pt, i, and plotDoE::x.

Referenced by fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_body_force_gradient_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ get_body_force_gradient_axi_nst()

virtual void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::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 
)
inlineprotectedvirtual

Get gradient of body force term at (Eulerian) position x. Computed via function pointer (if set) or by finite differencing (default)

268  {
269  // hierher: Implement function pointer version
270  /* //If no gradient function has been set, FD it */
271  /* if(Body_force_fct_gradient_pt==0) */
272  {
273  // Reference value
274  Vector<double> body_force(3, 0.0);
275  get_body_force_axi_nst(time, ipt, s, x, body_force);
276 
277  // FD it
278  const double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
279  Vector<double> body_force_pls(3, 0.0);
280  Vector<double> x_pls(x);
281  for (unsigned i = 0; i < 2; i++)
282  {
283  x_pls[i] += eps_fd;
284  get_body_force_axi_nst(time, ipt, s, x_pls, body_force_pls);
285  for (unsigned j = 0; j < 3; j++)
286  {
287  d_body_force_dx(j, i) =
288  (body_force_pls[j] - body_force[j]) / eps_fd;
289  }
290  x_pls[i] = x[i];
291  }
292  }
293  /* else */
294  /* { */
295  /* // Get gradient */
296  /* (*Source_fct_gradient_pt)(time,x,gradient); */
297  /* } */
298  }
static double Default_fd_jacobian_step
Definition: elements.h:1198

References Global_Parameters::body_force(), oomph::GeneralisedElement::Default_fd_jacobian_step, get_body_force_axi_nst(), i, j, s, and plotDoE::x.

Referenced by get_dresidual_dnodal_coordinates(), and oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ get_dresidual_dnodal_coordinates()

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

Reimplemented in oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations.

2423  {
2424  throw OomphLibError(
2425  "This has not been checked for generalised Newtonian elements!",
2428 
2429  // Return immediately if there are no dofs
2430  if (ndof() == 0)
2431  {
2432  return;
2433  }
2434 
2435  // Determine number of nodes in element
2436  const unsigned n_node = nnode();
2437 
2438  // Get continuous time from timestepper of first node
2439  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
2440 
2441  // Determine number of pressure dofs in element
2442  const unsigned n_pres = npres_axi_nst();
2443 
2444  // Find the indices at which the local velocities are stored
2445  unsigned u_nodal_index[3];
2446  for (unsigned i = 0; i < 3; i++)
2447  {
2448  u_nodal_index[i] = u_index_axi_nst(i);
2449  }
2450 
2451  // Set up memory for the shape and test functions
2452  // Note that there are only two dimensions, r and z, in this problem
2453  Shape psif(n_node), testf(n_node);
2454  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
2455 
2456  // Set up memory for pressure shape and test functions
2457  Shape psip(n_pres), testp(n_pres);
2458 
2459  // Deriatives of shape fct derivatives w.r.t. nodal coords
2460  RankFourTensor<double> d_dpsifdx_dX(2, n_node, n_node, 2);
2461  RankFourTensor<double> d_dtestfdx_dX(2, n_node, n_node, 2);
2462 
2463  // Derivative of Jacobian of mapping w.r.t. to nodal coords
2464  DenseMatrix<double> dJ_dX(2, n_node);
2465 
2466  // Derivatives of derivative of u w.r.t. nodal coords
2467  // Note that the entry d_dudx_dX(p,q,i,k) contains the derivative w.r.t.
2468  // nodal coordinate X_pq of du_i/dx_k. Since there are three velocity
2469  // components, i can take on the values 0, 1 and 2, while k and p only
2470  // take on the values 0 and 1 since there are only two spatial dimensions.
2471  RankFourTensor<double> d_dudx_dX(2, n_node, 3, 2);
2472 
2473  // Derivatives of nodal velocities w.r.t. nodal coords:
2474  // Assumption: Interaction only local via no-slip so
2475  // X_pq only affects U_iq.
2476  // Note that the entry d_U_dX(p,q,i) contains the derivative w.r.t. nodal
2477  // coordinate X_pq of nodal value U_iq. In other words this entry is the
2478  // derivative of the i-th velocity component w.r.t. the p-th spatial
2479  // coordinate, taken at the q-th local node.
2480  RankThreeTensor<double> d_U_dX(2, n_node, 3, 0.0);
2481 
2482  // Determine the number of integration points
2483  const unsigned n_intpt = integral_pt()->nweight();
2484 
2485  // Vector to hold local coordinates (two dimensions)
2486  Vector<double> s(2);
2487 
2488  // Get physical variables from element
2489  // (Reynolds number must be multiplied by the density ratio)
2490  const double scaled_re = re() * density_ratio();
2491  const double scaled_re_st = re_st() * density_ratio();
2492  const double scaled_re_inv_fr = re_invfr() * density_ratio();
2493  const double scaled_re_inv_ro = re_invro() * density_ratio();
2494  const double visc_ratio = viscosity_ratio();
2495  const Vector<double> G = g();
2496 
2497  // FD step
2499 
2500  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
2501  // Assumption: Interaction only local via no-slip so
2502  // X_ij only affects U_ij.
2503  bool element_has_node_with_aux_node_update_fct = false;
2504  for (unsigned q = 0; q < n_node; q++)
2505  {
2506  // Get pointer to q-th local node
2507  Node* nod_pt = node_pt(q);
2508 
2509  // Only compute if there's a node-update fct involved
2510  if (nod_pt->has_auxiliary_node_update_fct_pt())
2511  {
2512  element_has_node_with_aux_node_update_fct = true;
2513 
2514  // This functionality has not been tested so issue a warning
2515  // to this effect
2516  std::ostringstream warning_stream;
2517  warning_stream << "\nThe functionality to evaluate the additional"
2518  << "\ncontribution to the deriv of the residual eqn"
2519  << "\nw.r.t. the nodal coordinates which comes about"
2520  << "\nif a node's values are updated using an auxiliary"
2521  << "\nnode update function has NOT been tested for"
2522  << "\naxisymmetric Navier-Stokes elements. Use at your"
2523  << "\nown risk" << std::endl;
2524  OomphLibWarning(warning_stream.str(),
2525  "GeneralisedNewtonianAxisymmetricNavierStokesEquations:"
2526  ":get_dresidual_dnodal_coordinates",
2528 
2529  // Current nodal velocity
2530  Vector<double> u_ref(3);
2531  for (unsigned i = 0; i < 3; i++)
2532  {
2533  u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
2534  }
2535 
2536  // FD
2537  for (unsigned p = 0; p < 2; p++)
2538  {
2539  // Make backup
2540  double backup = nod_pt->x(p);
2541 
2542  // Do FD step. No node update required as we're
2543  // attacking the coordinate directly...
2544  nod_pt->x(p) += eps_fd;
2545 
2546  // Do auxiliary node update (to apply no slip)
2547  nod_pt->perform_auxiliary_node_update_fct();
2548 
2549  // Loop over velocity components
2550  for (unsigned i = 0; i < 3; i++)
2551  {
2552  // Evaluate
2553  d_U_dX(p, q, i) =
2554  (*(nod_pt->value_pt(u_nodal_index[i])) - u_ref[i]) / eps_fd;
2555  }
2556 
2557  // Reset
2558  nod_pt->x(p) = backup;
2559 
2560  // Do auxiliary node update (to apply no slip)
2561  nod_pt->perform_auxiliary_node_update_fct();
2562  }
2563  }
2564  }
2565 
2566  // Integer to store the local equation number
2567  int local_eqn = 0;
2568 
2569  // Loop over the integration points
2570  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2571  {
2572  // Assign values of s
2573  for (unsigned i = 0; i < 2; i++)
2574  {
2575  s[i] = integral_pt()->knot(ipt, i);
2576  }
2577 
2578  // Get the integral weight
2579  const double w = integral_pt()->weight(ipt);
2580 
2581  // Call the derivatives of the shape and test functions
2582  const double J = dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,
2583  psif,
2584  dpsifdx,
2585  d_dpsifdx_dX,
2586  testf,
2587  dtestfdx,
2588  d_dtestfdx_dX,
2589  dJ_dX);
2590 
2591  // Call the pressure shape and test functions
2592  pshape_axi_nst(s, psip, testp);
2593 
2594  // Allocate storage for the position and the derivative of the
2595  // mesh positions w.r.t. time
2596  Vector<double> interpolated_x(2, 0.0);
2597  Vector<double> mesh_velocity(2, 0.0);
2598 
2599  // Allocate storage for the pressure, velocity components and their
2600  // time and spatial derivatives
2601  double interpolated_p = 0.0;
2602  Vector<double> interpolated_u(3, 0.0);
2603  Vector<double> dudt(3, 0.0);
2604  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
2605 
2606  // Calculate pressure at integration point
2607  for (unsigned l = 0; l < n_pres; l++)
2608  {
2609  interpolated_p += p_axi_nst(l) * psip[l];
2610  }
2611 
2612  // Calculate velocities and derivatives at integration point
2613  // ---------------------------------------------------------
2614 
2615  // Loop over nodes
2616  for (unsigned l = 0; l < n_node; l++)
2617  {
2618  // Cache the shape function
2619  const double psif_ = psif(l);
2620 
2621  // Loop over the two coordinate directions
2622  for (unsigned i = 0; i < 2; i++)
2623  {
2624  interpolated_x[i] += this->raw_nodal_position(l, i) * psif_;
2625  }
2626 
2627  // Loop over the three velocity directions
2628  for (unsigned i = 0; i < 3; i++)
2629  {
2630  // Get the nodal value
2631  const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
2632  interpolated_u[i] += u_value * psif_;
2633  dudt[i] += du_dt_axi_nst(l, i) * psif_;
2634 
2635  // Loop over derivative directions
2636  for (unsigned j = 0; j < 2; j++)
2637  {
2638  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
2639  }
2640  }
2641  }
2642 
2643  // Get the mesh velocity if ALE is enabled
2644  if (!ALE_is_disabled)
2645  {
2646  // Loop over nodes
2647  for (unsigned l = 0; l < n_node; l++)
2648  {
2649  // Loop over the two coordinate directions
2650  for (unsigned i = 0; i < 2; i++)
2651  {
2652  mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif[l];
2653  }
2654  }
2655  }
2656 
2657  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
2658  for (unsigned q = 0; q < n_node; q++)
2659  {
2660  // Loop over the two coordinate directions
2661  for (unsigned p = 0; p < 2; p++)
2662  {
2663  // Loop over the three velocity components
2664  for (unsigned i = 0; i < 3; i++)
2665  {
2666  // Loop over the two coordinate directions
2667  for (unsigned k = 0; k < 2; k++)
2668  {
2669  double aux = 0.0;
2670 
2671  // Loop over nodes and add contribution
2672  for (unsigned j = 0; j < n_node; j++)
2673  {
2674  aux += this->raw_nodal_value(j, u_nodal_index[i]) *
2675  d_dpsifdx_dX(p, q, j, k);
2676  }
2677  d_dudx_dX(p, q, i, k) = aux;
2678  }
2679  }
2680  }
2681  }
2682 
2683  // Get weight of actual nodal position/value in computation of mesh
2684  // velocity from positional/value time stepper
2685  const double pos_time_weight =
2687  const double val_time_weight =
2688  node_pt(0)->time_stepper_pt()->weight(1, 0);
2689 
2690  // Get the user-defined body force terms
2691  Vector<double> body_force(3);
2693 
2694  // Get the user-defined source function
2695  const double source = get_source_fct(time, ipt, interpolated_x);
2696 
2697  // Note: Can use raw values and avoid bypassing hanging information
2698  // because this is the non-refineable version!
2699 
2700  // Get gradient of body force function
2701  DenseMatrix<double> d_body_force_dx(3, 2, 0.0);
2703  time, ipt, s, interpolated_x, d_body_force_dx);
2704 
2705  // Get gradient of source function
2706  Vector<double> source_gradient(2, 0.0);
2707  get_source_fct_gradient(time, ipt, interpolated_x, source_gradient);
2708 
2709  // Cache r (the first position component)
2710  const double r = interpolated_x[0];
2711 
2712  // Assemble shape derivatives
2713  // --------------------------
2714 
2715  // ==================
2716  // MOMENTUM EQUATIONS
2717  // ==================
2718 
2719  // Loop over the test functions
2720  for (unsigned l = 0; l < n_node; l++)
2721  {
2722  // Cache the test function
2723  const double testf_ = testf[l];
2724 
2725  // --------------------------------
2726  // FIRST (RADIAL) MOMENTUM EQUATION
2727  // --------------------------------
2728  local_eqn = nodal_local_eqn(l, u_nodal_index[0]);
2729  ;
2730 
2731  // If it's not a boundary condition
2732  if (local_eqn >= 0)
2733  {
2734  // Loop over the two coordinate directions
2735  for (unsigned p = 0; p < 2; p++)
2736  {
2737  // Loop over nodes
2738  for (unsigned q = 0; q < n_node; q++)
2739  {
2740  // Residual x deriv of Jacobian
2741  // ----------------------------
2742 
2743  // Add the user-defined body force terms
2744  double sum = r * body_force[0] * testf_;
2745 
2746  // Add the gravitational body force term
2747  sum += r * scaled_re_inv_fr * testf_ * G[0];
2748 
2749  // Add the pressure gradient term
2750  sum += interpolated_p * (testf_ + r * dtestfdx(l, 0));
2751 
2752  // Add the stress tensor terms
2753  // The viscosity ratio needs to go in here to ensure
2754  // continuity of normal stress is satisfied even in flows
2755  // with zero pressure gradient!
2756  sum -= visc_ratio * r * (1.0 + Gamma[0]) *
2757  interpolated_dudx(0, 0) * dtestfdx(l, 0);
2758 
2759  sum -=
2760  visc_ratio * r *
2761  (interpolated_dudx(0, 1) + Gamma[0] * interpolated_dudx(1, 0)) *
2762  dtestfdx(l, 1);
2763 
2764  sum -=
2765  visc_ratio * (1.0 + Gamma[0]) * interpolated_u[0] * testf_ / r;
2766 
2767  // Add the du/dt term
2768  sum -= scaled_re_st * r * dudt[0] * testf_;
2769 
2770  // Add the convective terms
2771  sum -= scaled_re *
2772  (r * interpolated_u[0] * interpolated_dudx(0, 0) -
2773  interpolated_u[2] * interpolated_u[2] +
2774  r * interpolated_u[1] * interpolated_dudx(0, 1)) *
2775  testf_;
2776 
2777  // Add the mesh velocity terms
2778  if (!ALE_is_disabled)
2779  {
2780  for (unsigned k = 0; k < 2; k++)
2781  {
2782  sum += scaled_re_st * r * mesh_velocity[k] *
2783  interpolated_dudx(0, k) * testf_;
2784  }
2785  }
2786 
2787  // Add the Coriolis term
2788  sum += 2.0 * r * scaled_re_inv_ro * interpolated_u[2] * testf_;
2789 
2790  // Multiply through by deriv of Jacobian and integration weight
2791  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2792  sum * dJ_dX(p, q) * w;
2793 
2794  // Derivative of residual x Jacobian
2795  // ---------------------------------
2796 
2797  // Body force
2798  sum = r * d_body_force_dx(0, p) * psif[q] * testf_;
2799  if (p == 0)
2800  {
2801  sum += body_force[0] * psif[q] * testf_;
2802  }
2803 
2804  // Gravitational body force
2805  if (p == 0)
2806  {
2807  sum += scaled_re_inv_fr * G[0] * psif[q] * testf_;
2808  }
2809 
2810  // Pressure gradient term
2811  sum += r * interpolated_p * d_dtestfdx_dX(p, q, l, 0);
2812  if (p == 0)
2813  {
2814  sum += interpolated_p * psif[q] * dtestfdx(l, 0);
2815  }
2816 
2817  // Viscous terms
2818  sum -=
2819  r * visc_ratio *
2820  ((1.0 + Gamma[0]) *
2821  (d_dudx_dX(p, q, 0, 0) * dtestfdx(l, 0) +
2822  interpolated_dudx(0, 0) * d_dtestfdx_dX(p, q, l, 0)) +
2823  (d_dudx_dX(p, q, 0, 1) + Gamma[0] * d_dudx_dX(p, q, 1, 0)) *
2824  dtestfdx(l, 1) +
2825  (interpolated_dudx(0, 1) +
2826  Gamma[0] * interpolated_dudx(1, 0)) *
2827  d_dtestfdx_dX(p, q, l, 1));
2828  if (p == 0)
2829  {
2830  sum -= visc_ratio *
2831  ((1.0 + Gamma[0]) *
2832  (interpolated_dudx(0, 0) * psif[q] * dtestfdx(l, 0) -
2833  interpolated_u[0] * psif[q] * testf_ / (r * r)) +
2834  (interpolated_dudx(0, 1) +
2835  Gamma[0] * interpolated_dudx(1, 0)) *
2836  psif[q] * dtestfdx(l, 1));
2837  }
2838 
2839  // Convective terms, including mesh velocity
2840  for (unsigned k = 0; k < 2; k++)
2841  {
2842  double tmp = scaled_re * interpolated_u[k];
2843  if (!ALE_is_disabled)
2844  {
2845  tmp -= scaled_re_st * mesh_velocity[k];
2846  }
2847  sum -= r * tmp * d_dudx_dX(p, q, 0, k) * testf_;
2848  if (p == 0)
2849  {
2850  sum -= tmp * interpolated_dudx(0, k) * psif[q] * testf_;
2851  }
2852  }
2853  if (!ALE_is_disabled)
2854  {
2855  sum += scaled_re_st * r * pos_time_weight *
2856  interpolated_dudx(0, p) * psif[q] * testf_;
2857  }
2858 
2859  // du/dt term
2860  if (p == 0)
2861  {
2862  sum -= scaled_re_st * dudt[0] * psif[q] * testf_;
2863  }
2864 
2865  // Coriolis term
2866  if (p == 0)
2867  {
2868  sum +=
2869  2.0 * scaled_re_inv_ro * interpolated_u[2] * psif[q] * testf_;
2870  }
2871 
2872  // Multiply through by Jacobian and integration weight
2873  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2874 
2875  // Derivs w.r.t. nodal velocities
2876  // ------------------------------
2877  if (element_has_node_with_aux_node_update_fct)
2878  {
2879  // Contribution from deriv of first velocity component
2880  double tmp =
2881  scaled_re_st * r * val_time_weight * psif[q] * testf_;
2882  tmp +=
2883  r * scaled_re * interpolated_dudx(0, 0) * psif[q] * testf_;
2884  for (unsigned k = 0; k < 2; k++)
2885  {
2886  double tmp2 = scaled_re * interpolated_u[k];
2887  if (!ALE_is_disabled)
2888  {
2889  tmp2 -= scaled_re_st * mesh_velocity[k];
2890  }
2891  tmp += r * tmp2 * dpsifdx(q, k) * testf_;
2892  }
2893 
2894  tmp += r * visc_ratio * (1.0 + Gamma[0]) * dpsifdx(q, 0) *
2895  dtestfdx(l, 0);
2896  tmp += r * visc_ratio * dpsifdx(q, 1) * dtestfdx(l, 1);
2897  tmp += visc_ratio * (1.0 + Gamma[0]) * psif[q] * testf_ / r;
2898 
2899  // Multiply through by dU_0q/dX_pq
2900  sum = -d_U_dX(p, q, 0) * tmp;
2901 
2902  // Contribution from deriv of second velocity component
2903  tmp =
2904  scaled_re * r * interpolated_dudx(0, 1) * psif[q] * testf_;
2905  tmp +=
2906  r * visc_ratio * Gamma[0] * dpsifdx(q, 0) * dtestfdx(l, 1);
2907 
2908  // Multiply through by dU_1q/dX_pq
2909  sum -= d_U_dX(p, q, 1) * tmp;
2910 
2911  // Contribution from deriv of third velocity component
2912  tmp = 2.0 *
2913  (r * scaled_re_inv_ro + scaled_re * interpolated_u[2]) *
2914  psif[q] * testf_;
2915 
2916  // Multiply through by dU_2q/dX_pq
2917  sum += d_U_dX(p, q, 2) * tmp;
2918 
2919  // Multiply through by Jacobian and integration weight
2920  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2921  }
2922  } // End of loop over q
2923  } // End of loop over p
2924  } // End of if it's not a boundary condition
2925 
2926  // --------------------------------
2927  // SECOND (AXIAL) MOMENTUM EQUATION
2928  // --------------------------------
2929  local_eqn = nodal_local_eqn(l, u_nodal_index[1]);
2930  ;
2931 
2932  // If it's not a boundary condition
2933  if (local_eqn >= 0)
2934  {
2935  // Loop over the two coordinate directions
2936  for (unsigned p = 0; p < 2; p++)
2937  {
2938  // Loop over nodes
2939  for (unsigned q = 0; q < n_node; q++)
2940  {
2941  // Residual x deriv of Jacobian
2942  // ----------------------------
2943 
2944  // Add the user-defined body force terms
2945  double sum = r * body_force[1] * testf_;
2946 
2947  // Add the gravitational body force term
2948  sum += r * scaled_re_inv_fr * testf_ * G[1];
2949 
2950  // Add the pressure gradient term
2951  sum += r * interpolated_p * dtestfdx(l, 1);
2952 
2953  // Add the stress tensor terms
2954  // The viscosity ratio needs to go in here to ensure
2955  // continuity of normal stress is satisfied even in flows
2956  // with zero pressure gradient!
2957  sum -=
2958  visc_ratio * r *
2959  (interpolated_dudx(1, 0) + Gamma[1] * interpolated_dudx(0, 1)) *
2960  dtestfdx(l, 0);
2961 
2962  sum -= visc_ratio * r * (1.0 + Gamma[1]) *
2963  interpolated_dudx(1, 1) * dtestfdx(l, 1);
2964 
2965  // Add the du/dt term
2966  sum -= scaled_re_st * r * dudt[1] * testf_;
2967 
2968  // Add the convective terms
2969  sum -= scaled_re *
2970  (r * interpolated_u[0] * interpolated_dudx(1, 0) +
2971  r * interpolated_u[1] * interpolated_dudx(1, 1)) *
2972  testf_;
2973 
2974  // Add the mesh velocity terms
2975  if (!ALE_is_disabled)
2976  {
2977  for (unsigned k = 0; k < 2; k++)
2978  {
2979  sum += scaled_re_st * r * mesh_velocity[k] *
2980  interpolated_dudx(1, k) * testf_;
2981  }
2982  }
2983 
2984  // Multiply through by deriv of Jacobian and integration weight
2985  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2986  sum * dJ_dX(p, q) * w;
2987 
2988  // Derivative of residual x Jacobian
2989  // ---------------------------------
2990 
2991  // Body force
2992  sum = r * d_body_force_dx(1, p) * psif[q] * testf_;
2993  if (p == 0)
2994  {
2995  sum += body_force[1] * psif[q] * testf_;
2996  }
2997 
2998  // Gravitational body force
2999  if (p == 0)
3000  {
3001  sum += scaled_re_inv_fr * G[1] * psif[q] * testf_;
3002  }
3003 
3004  // Pressure gradient term
3005  sum += r * interpolated_p * d_dtestfdx_dX(p, q, l, 1);
3006  if (p == 0)
3007  {
3008  sum += interpolated_p * psif[q] * dtestfdx(l, 1);
3009  }
3010 
3011  // Viscous terms
3012  sum -=
3013  r * visc_ratio *
3014  ((d_dudx_dX(p, q, 1, 0) + Gamma[1] * d_dudx_dX(p, q, 0, 1)) *
3015  dtestfdx(l, 0) +
3016  (interpolated_dudx(1, 0) +
3017  Gamma[1] * interpolated_dudx(0, 1)) *
3018  d_dtestfdx_dX(p, q, l, 0) +
3019  (1.0 + Gamma[1]) * d_dudx_dX(p, q, 1, 1) * dtestfdx(l, 1) +
3020  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
3021  d_dtestfdx_dX(p, q, l, 1));
3022  if (p == 0)
3023  {
3024  sum -=
3025  visc_ratio * ((interpolated_dudx(1, 0) +
3026  Gamma[1] * interpolated_dudx(0, 1)) *
3027  psif[q] * dtestfdx(l, 0) +
3028  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
3029  psif[q] * dtestfdx(l, 1));
3030  }
3031 
3032  // Convective terms, including mesh velocity
3033  for (unsigned k = 0; k < 2; k++)
3034  {
3035  double tmp = scaled_re * interpolated_u[k];
3036  if (!ALE_is_disabled)
3037  {
3038  tmp -= scaled_re_st * mesh_velocity[k];
3039  }
3040  sum -= r * tmp * d_dudx_dX(p, q, 1, k) * testf_;
3041  if (p == 0)
3042  {
3043  sum -= tmp * interpolated_dudx(1, k) * psif[q] * testf_;
3044  }
3045  }
3046  if (!ALE_is_disabled)
3047  {
3048  sum += scaled_re_st * r * pos_time_weight *
3049  interpolated_dudx(1, p) * psif[q] * testf_;
3050  }
3051 
3052  // du/dt term
3053  if (p == 0)
3054  {
3055  sum -= scaled_re_st * dudt[1] * psif[q] * testf_;
3056  }
3057 
3058  // Multiply through by Jacobian and integration weight
3059  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
3060 
3061  // Derivs w.r.t. to nodal velocities
3062  // ---------------------------------
3063  if (element_has_node_with_aux_node_update_fct)
3064  {
3065  // Contribution from deriv of first velocity component
3066  double tmp =
3067  scaled_re * r * interpolated_dudx(1, 0) * psif[q] * testf_;
3068  tmp +=
3069  r * visc_ratio * Gamma[1] * dpsifdx(q, 1) * dtestfdx(l, 0);
3070 
3071  // Multiply through by dU_0q/dX_pq
3072  sum = -d_U_dX(p, q, 0) * tmp;
3073 
3074  // Contribution from deriv of second velocity component
3075  tmp = scaled_re_st * r * val_time_weight * psif[q] * testf_;
3076  tmp +=
3077  r * scaled_re * interpolated_dudx(1, 1) * psif[q] * testf_;
3078  for (unsigned k = 0; k < 2; k++)
3079  {
3080  double tmp2 = scaled_re * interpolated_u[k];
3081  if (!ALE_is_disabled)
3082  {
3083  tmp2 -= scaled_re_st * mesh_velocity[k];
3084  }
3085  tmp += r * tmp2 * dpsifdx(q, k) * testf_;
3086  }
3087  tmp += r * visc_ratio *
3088  (dpsifdx(q, 0) * dtestfdx(l, 0) +
3089  (1.0 + Gamma[1]) * dpsifdx(q, 1) * dtestfdx(l, 1));
3090 
3091  // Multiply through by dU_1q/dX_pq
3092  sum -= d_U_dX(p, q, 1) * tmp;
3093 
3094  // Multiply through by Jacobian and integration weight
3095  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
3096  }
3097  } // End of loop over q
3098  } // End of loop over p
3099  } // End of if it's not a boundary condition
3100 
3101  // -----------------------------------
3102  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
3103  // -----------------------------------
3104  local_eqn = nodal_local_eqn(l, u_nodal_index[2]);
3105  ;
3106 
3107  // If it's not a boundary condition
3108  if (local_eqn >= 0)
3109  {
3110  // Loop over the two coordinate directions
3111  for (unsigned p = 0; p < 2; p++)
3112  {
3113  // Loop over nodes
3114  for (unsigned q = 0; q < n_node; q++)
3115  {
3116  // Residual x deriv of Jacobian
3117  // ----------------------------
3118 
3119  // Add the user-defined body force terms
3120  double sum = r * body_force[2] * testf_;
3121 
3122  // Add the gravitational body force term
3123  sum += r * scaled_re_inv_fr * testf_ * G[2];
3124 
3125  // There is NO pressure gradient term
3126 
3127  // Add in the stress tensor terms
3128  // The viscosity ratio needs to go in here to ensure
3129  // continuity of normal stress is satisfied even in flows
3130  // with zero pressure gradient!
3131  sum -=
3132  visc_ratio *
3133  (r * interpolated_dudx(2, 0) - Gamma[0] * interpolated_u[2]) *
3134  dtestfdx(l, 0);
3135 
3136  sum -= visc_ratio * r * interpolated_dudx(2, 1) * dtestfdx(l, 1);
3137 
3138  sum -=
3139  visc_ratio *
3140  ((interpolated_u[2] / r) - Gamma[0] * interpolated_dudx(2, 0)) *
3141  testf_;
3142 
3143  // Add the du/dt term
3144  sum -= scaled_re_st * r * dudt[2] * testf_;
3145 
3146  // Add the convective terms
3147  sum -= scaled_re *
3148  (r * interpolated_u[0] * interpolated_dudx(2, 0) +
3149  interpolated_u[0] * interpolated_u[2] +
3150  r * interpolated_u[1] * interpolated_dudx(2, 1)) *
3151  testf_;
3152 
3153  // Add the mesh velocity terms
3154  if (!ALE_is_disabled)
3155  {
3156  for (unsigned k = 0; k < 2; k++)
3157  {
3158  sum += scaled_re_st * r * mesh_velocity[k] *
3159  interpolated_dudx(2, k) * testf_;
3160  }
3161  }
3162 
3163  // Add the Coriolis term
3164  sum -= 2.0 * r * scaled_re_inv_ro * interpolated_u[0] * testf_;
3165 
3166  // Multiply through by deriv of Jacobian and integration weight
3167  dresidual_dnodal_coordinates(local_eqn, p, q) +=
3168  sum * dJ_dX(p, q) * w;
3169 
3170  // Derivative of residual x Jacobian
3171  // ---------------------------------
3172 
3173  // Body force
3174  sum = r * d_body_force_dx(2, p) * psif[q] * testf_;
3175  if (p == 0)
3176  {
3177  sum += body_force[2] * psif[q] * testf_;
3178  }
3179 
3180  // Gravitational body force
3181  if (p == 0)
3182  {
3183  sum += scaled_re_inv_fr * G[2] * psif[q] * testf_;
3184  }
3185 
3186  // There is NO pressure gradient term
3187 
3188  // Viscous terms
3189  sum -= visc_ratio * r *
3190  (d_dudx_dX(p, q, 2, 0) * dtestfdx(l, 0) +
3191  interpolated_dudx(2, 0) * d_dtestfdx_dX(p, q, l, 0) +
3192  d_dudx_dX(p, q, 2, 1) * dtestfdx(l, 1) +
3193  interpolated_dudx(2, 1) * d_dtestfdx_dX(p, q, l, 1));
3194 
3195  sum += visc_ratio * Gamma[0] * d_dudx_dX(p, q, 2, 0) * testf_;
3196 
3197  if (p == 0)
3198  {
3199  sum -= visc_ratio *
3200  (interpolated_dudx(2, 0) * psif[q] * dtestfdx(l, 0) +
3201  interpolated_dudx(2, 1) * psif[q] * dtestfdx(l, 1) +
3202  interpolated_u[2] * psif[q] * testf_ / (r * r));
3203  }
3204 
3205  // Convective terms, including mesh velocity
3206  for (unsigned k = 0; k < 2; k++)
3207  {
3208  double tmp = scaled_re * interpolated_u[k];
3209  if (!ALE_is_disabled)
3210  {
3211  tmp -= scaled_re_st * mesh_velocity[k];
3212  }
3213  sum -= r * tmp * d_dudx_dX(p, q, 2, k) * testf_;
3214  if (p == 0)
3215  {
3216  sum -= tmp * interpolated_dudx(2, k) * psif[q] * testf_;
3217  }
3218  }
3219  if (!ALE_is_disabled)
3220  {
3221  sum += scaled_re_st * r * pos_time_weight *
3222  interpolated_dudx(2, p) * psif[q] * testf_;
3223  }
3224 
3225  // du/dt term
3226  if (p == 0)
3227  {
3228  sum -= scaled_re_st * dudt[2] * psif[q] * testf_;
3229  }
3230 
3231  // Coriolis term
3232  if (p == 0)
3233  {
3234  sum -=
3235  2.0 * scaled_re_inv_ro * interpolated_u[0] * psif[q] * testf_;
3236  }
3237 
3238  // Multiply through by Jacobian and integration weight
3239  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
3240 
3241  // Derivs w.r.t. to nodal velocities
3242  // ---------------------------------
3243  if (element_has_node_with_aux_node_update_fct)
3244  {
3245  // Contribution from deriv of first velocity component
3246  double tmp = (2.0 * r * scaled_re_inv_ro +
3247  r * scaled_re * interpolated_dudx(2, 0) -
3248  scaled_re * interpolated_u[2]) *
3249  psif[q] * testf_;
3250 
3251  // Multiply through by dU_0q/dX_pq
3252  sum = -d_U_dX(p, q, 0) * tmp;
3253 
3254  // Contribution from deriv of second velocity component
3255  // Multiply through by dU_1q/dX_pq
3256  sum -= d_U_dX(p, q, 1) * scaled_re * r *
3257  interpolated_dudx(2, 1) * psif[q] * testf_;
3258 
3259  // Contribution from deriv of third velocity component
3260  tmp = scaled_re_st * r * val_time_weight * psif[q] * testf_;
3261  tmp -= scaled_re * interpolated_u[0] * psif[q] * testf_;
3262  for (unsigned k = 0; k < 2; k++)
3263  {
3264  double tmp2 = scaled_re * interpolated_u[k];
3265  if (!ALE_is_disabled)
3266  {
3267  tmp2 -= scaled_re_st * mesh_velocity[k];
3268  }
3269  tmp += r * tmp2 * dpsifdx(q, k) * testf_;
3270  }
3271  tmp += visc_ratio * (r * dpsifdx(q, 0) - Gamma[0] * psif[q]) *
3272  dtestfdx(l, 0);
3273  tmp += r * visc_ratio * dpsifdx(q, 1) * dtestfdx(l, 1);
3274  tmp += visc_ratio * ((psif[q] / r) - Gamma[0] * dpsifdx(q, 0)) *
3275  testf_;
3276 
3277  // Multiply through by dU_2q/dX_pq
3278  sum -= d_U_dX(p, q, 2) * tmp;
3279 
3280  // Multiply through by Jacobian and integration weight
3281  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
3282  }
3283  } // End of loop over q
3284  } // End of loop over p
3285  } // End of if it's not a boundary condition
3286 
3287  } // End of loop over test functions
3288 
3289 
3290  // ===================
3291  // CONTINUITY EQUATION
3292  // ===================
3293 
3294  // Loop over the shape functions
3295  for (unsigned l = 0; l < n_pres; l++)
3296  {
3297  local_eqn = p_local_eqn(l);
3298 
3299  // Cache the test function
3300  const double testp_ = testp[l];
3301 
3302  // If not a boundary conditions
3303  if (local_eqn >= 0)
3304  {
3305  // Loop over the two coordinate directions
3306  for (unsigned p = 0; p < 2; p++)
3307  {
3308  // Loop over nodes
3309  for (unsigned q = 0; q < n_node; q++)
3310  {
3311  // Residual x deriv of Jacobian
3312  //-----------------------------
3313 
3314  // Source term
3315  double aux = -r * source;
3316 
3317  // Gradient terms
3318  aux += (interpolated_u[0] + r * interpolated_dudx(0, 0) +
3319  r * interpolated_dudx(1, 1));
3320 
3321  // Multiply through by deriv of Jacobian and integration weight
3322  dresidual_dnodal_coordinates(local_eqn, p, q) +=
3323  aux * dJ_dX(p, q) * testp_ * w;
3324 
3325  // Derivative of residual x Jacobian
3326  // ---------------------------------
3327 
3328  // Gradient of source function
3329  aux = -r * source_gradient[p] * psif[q];
3330  if (p == 0)
3331  {
3332  aux -= source * psif[q];
3333  }
3334 
3335  // Gradient terms
3336  aux += r * (d_dudx_dX(p, q, 0, 0) + d_dudx_dX(p, q, 1, 1));
3337  if (p == 0)
3338  {
3339  aux +=
3340  (interpolated_dudx(0, 0) + interpolated_dudx(1, 1)) * psif[q];
3341  }
3342 
3343  // Derivs w.r.t. nodal velocities
3344  // ------------------------------
3345  if (element_has_node_with_aux_node_update_fct)
3346  {
3347  aux += d_U_dX(p, q, 0) * (psif[q] + r * dpsifdx(q, 0));
3348  aux += d_U_dX(p, q, 1) * r * dpsifdx(q, 1);
3349  }
3350 
3351  // Multiply through by Jacobian, test fct and integration weight
3352  dresidual_dnodal_coordinates(local_eqn, p, q) +=
3353  aux * testp_ * J * w;
3354  }
3355  }
3356  }
3357  } // End of loop over shape functions for continuity eqn
3358 
3359  } // End of loop over integration points
3360  }
float * p
Definition: Tutorial_Map_using.cpp:9
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
const double & viscosity_ratio() const
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:494
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
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
Eigen::Matrix< Scalar, Dynamic, Dynamic, ColMajor > tmp
Definition: level3_impl.h:365
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019

References ALE_is_disabled, Global_Parameters::body_force(), oomph::GeneralisedElement::Default_fd_jacobian_step, density_ratio(), dshape_and_dtest_eulerian_at_knot_axi_nst(), du_dt_axi_nst(), G, g(), Gamma, get_body_force_axi_nst(), get_body_force_gradient_axi_nst(), get_source_fct(), get_source_fct_gradient(), oomph::Node::has_auxiliary_node_update_fct_pt(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), J, j, k, oomph::Integral::knot(), oomph::GeneralisedElement::ndof(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::node_pt(), npres_axi_nst(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p, p_axi_nst(), p_local_eqn(), oomph::Node::perform_auxiliary_node_update_fct(), oomph::Node::position_time_stepper_pt(), pshape_axi_nst(), Eigen::numext::q, UniformPSDSelfTest::r, oomph::FiniteElement::raw_dnodal_position_dt(), oomph::FiniteElement::raw_nodal_position(), oomph::FiniteElement::raw_nodal_value(), re(), re_invfr(), re_invro(), re_st(), s, TestProblem::source(), oomph::Time::time(), oomph::TimeStepper::time_pt(), oomph::Data::time_stepper_pt(), tmp, u_index_axi_nst(), oomph::Data::value_pt(), viscosity_ratio(), w, oomph::Integral::weight(), oomph::TimeStepper::weight(), and oomph::Node::x().

◆ get_pressure_and_velocity_mass_matrix_diagonal()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_pressure_and_velocity_mass_matrix_diagonal ( Vector< double > &  press_mass_diag,
Vector< double > &  veloc_mass_diag,
const unsigned which_one = 0 
)
virtual

Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed, otherwise only the pressure (which_one=1) or the velocity mass matrix (which_one=2 – the LSC version of the preconditioner only needs that one) NOTE: pressure versions isn't implemented yet because this element has never been tried with Fp preconditoner.

Implements oomph::NavierStokesElementWithDiagonalMassMatrices.

75  {
76 #ifdef PARANOID
77  if ((which_one == 0) || (which_one == 1))
78  {
79  throw OomphLibError("Computation of diagonal of pressure mass matrix is "
80  "not impmented yet.\n",
83  }
84 #endif
85 
86  // Resize and initialise
87  veloc_mass_diag.assign(ndof(), 0.0);
88 
89  // find out how many nodes there are
90  const unsigned n_node = nnode();
91 
92  // find number of coordinates
93  const unsigned n_value = 3;
94 
95  // find the indices at which the local velocities are stored
96  Vector<unsigned> u_nodal_index(n_value);
97  for (unsigned i = 0; i < n_value; i++)
98  {
99  u_nodal_index[i] = this->u_index_axi_nst(i);
100  }
101 
102  // Set up memory for test functions
103  Shape test(n_node);
104 
105  // Number of integration points
106  unsigned n_intpt = integral_pt()->nweight();
107 
108  // Integer to store the local equations no
109  int local_eqn = 0;
110 
111  // Loop over the integration points
112  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
113  {
114  // Get the integral weight
115  double w = integral_pt()->weight(ipt);
116 
117  // Get determinant of Jacobian of the mapping
118  double J = J_eulerian_at_knot(ipt);
119 
120  // Premultiply weights and Jacobian
121  double W = w * J;
122 
123  // Get the velocity test functions - there is no explicit
124  // function to give the test function so use shape
125  shape_at_knot(ipt, test);
126 
127  // Need to get the position to sort out the jacobian properly
128  double r = 0.0;
129  for (unsigned l = 0; l < n_node; l++)
130  {
131  r += this->nodal_position(l, 0) * test(l);
132  }
133 
134  // Multiply by the geometric factor
135  W *= r;
136 
137  // Loop over the veclocity test functions
138  for (unsigned l = 0; l < n_node; l++)
139  {
140  // Loop over the velocity components
141  for (unsigned i = 0; i < n_value; i++)
142  {
143  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
144 
145  // If not a boundary condition
146  if (local_eqn >= 0)
147  {
148  // Add the contribution
149  veloc_mass_diag[local_eqn] += test[l] * test[l] * W;
150  } // End of if not boundary condition statement
151  } // End of loop over dimension
152  } // End of loop over test functions
153  }
154  }
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317
virtual double J_eulerian_at_knot(const unsigned &ipt) const
Definition: elements.cc:4168
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Definition: elements.cc:3220
squared absolute sa ArrayBase::abs2 DOXCOMMA MatrixBase::cwiseAbs2 sa Eigen::abs2 DOXCOMMA Eigen::pow DOXCOMMA ArrayBase::square nearest sa Eigen::floor DOXCOMMA Eigen::ceil DOXCOMMA ArrayBase::round nearest integer not less than the given sa Eigen::floor DOXCOMMA ArrayBase::ceil not a number test
Definition: GlobalFunctions.h:109
Definition: indexed_view.cpp:20

References i, oomph::FiniteElement::integral_pt(), J, oomph::FiniteElement::J_eulerian_at_knot(), oomph::GeneralisedElement::ndof(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::nodal_position(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, UniformPSDSelfTest::r, oomph::FiniteElement::shape_at_knot(), Eigen::test, u_index_axi_nst(), w, oomph::QuadTreeNames::W, and oomph::Integral::weight().

◆ get_source_fct()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_source_fct ( const double time,
const unsigned ipt,
const Vector< double > &  x 
)
inlineprotected

Calculate the source fct at given time and Eulerian position.

304  {
305  // If the function pointer is zero return zero
306  if (Source_fct_pt == 0)
307  {
308  return 0;
309  }
310 
311  // Otherwise call the function
312  else
313  {
314  return (*Source_fct_pt)(time, x);
315  }
316  }

References Source_fct_pt, and plotDoE::x.

Referenced by fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates(), and get_source_fct_gradient().

◆ get_source_fct_gradient()

virtual void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_source_fct_gradient ( const double time,
const unsigned ipt,
const Vector< double > &  x,
Vector< double > &  gradient 
)
inlineprotectedvirtual

Get gradient of source term at (Eulerian) position x. Computed via function pointer (if set) or by finite differencing (default)

324  {
325  // hierher: Implement function pointer version
326  /* //If no gradient function has been set, FD it */
327  /* if(Source_fct_gradient_pt==0) */
328  {
329  // Reference value
330  const double source = get_source_fct(time, ipt, x);
331 
332  // FD it
333  const double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
334  double source_pls = 0.0;
335  Vector<double> x_pls(x);
336  for (unsigned i = 0; i < 2; i++)
337  {
338  x_pls[i] += eps_fd;
339  source_pls = get_source_fct(time, ipt, x_pls);
340  gradient[i] = (source_pls - source) / eps_fd;
341  x_pls[i] = x[i];
342  }
343  }
344  /* else */
345  /* { */
346  /* // Get gradient */
347  /* (*Source_fct_gradient_pt)(time,x,gradient); */
348  /* } */
349  }

References oomph::GeneralisedElement::Default_fd_jacobian_step, get_source_fct(), i, TestProblem::source(), and plotDoE::x.

Referenced by get_dresidual_dnodal_coordinates(), and oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ get_viscosity_ratio_axisym_nst()

virtual void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::get_viscosity_ratio_axisym_nst ( const unsigned ipt,
const Vector< double > &  s,
const Vector< double > &  x,
double visc_ratio 
)
inlineprotectedvirtual

Calculate the viscosity ratio relative to the viscosity used in the definition of the Reynolds number at given time and Eulerian position This is the equivalent of the Newtonian case and basically allows the flow of two Generalised Newtonian fluids with different reference viscosities

361  {
362  visc_ratio = *Viscosity_Ratio_pt;
363  }

References Viscosity_Ratio_pt.

Referenced by fill_in_generic_dresidual_contribution_axi_nst(), and fill_in_generic_residual_contribution_axi_nst().

◆ interpolated_d_dudx_dX_axi_nst()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::interpolated_d_dudx_dX_axi_nst ( const Vector< double > &  s,
const unsigned p,
const unsigned q,
const unsigned i,
const unsigned k 
) const
inline

Return FE interpolated derivatives w.r.t. nodal coordinates X_{pq} of the spatial derivatives of the velocity components du_i/dx_k at local coordinate s

1198  {
1199  // Determine number of nodes
1200  const unsigned n_node = nnode();
1201 
1202  // Allocate storage for local shape function and its derivatives
1203  // with respect to space
1204  Shape psif(n_node);
1205  DShape dpsifds(n_node, 2);
1206 
1207  // Find values of shape function (ignore jacobian)
1208  (void)this->dshape_local(s, psif, dpsifds);
1209 
1210  // Allocate memory for the jacobian and the inverse of the jacobian
1211  DenseMatrix<double> jacobian(2), inverse_jacobian(2);
1212 
1213  // Allocate memory for the derivative of the jacobian w.r.t. nodal coords
1214  DenseMatrix<double> djacobian_dX(2, n_node);
1215 
1216  // Allocate memory for the derivative w.r.t. nodal coords of the
1217  // spatial derivatives of the shape functions
1218  RankFourTensor<double> d_dpsifdx_dX(2, n_node, 3, 2);
1219 
1220  // Now calculate the inverse jacobian
1221  const double det =
1222  local_to_eulerian_mapping(dpsifds, jacobian, inverse_jacobian);
1223 
1224  // Calculate the derivative of the jacobian w.r.t. nodal coordinates
1225  // Note: must call this before "transform_derivatives(...)" since this
1226  // function requires dpsids rather than dpsidx
1227  dJ_eulerian_dnodal_coordinates(jacobian, dpsifds, djacobian_dX);
1228 
1229  // Calculate the derivative of dpsidx w.r.t. nodal coordinates
1230  // Note: this function also requires dpsids rather than dpsidx
1232  det, jacobian, djacobian_dX, inverse_jacobian, dpsifds, d_dpsifdx_dX);
1233 
1234  // Get the index at which the velocity is stored
1235  const unsigned u_nodal_index = u_index_axi_nst(i);
1236 
1237  // Initialise value of dudx
1238  double interpolated_d_dudx_dX = 0.0;
1239 
1240  // Loop over the local nodes and sum
1241  for (unsigned l = 0; l < n_node; l++)
1242  {
1243  interpolated_d_dudx_dX +=
1244  nodal_value(l, u_nodal_index) * d_dpsifdx_dX(p, q, l, k);
1245  }
1246 
1247  return (interpolated_d_dudx_dX);
1248  }
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
Definition: elements.cc:2749
virtual void dJ_eulerian_dnodal_coordinates(const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
Definition: elements.cc:2669
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Definition: elements.h:1508
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Definition: elements.h:1981

References oomph::FiniteElement::d_dshape_eulerian_dnodal_coordinates(), oomph::FiniteElement::dJ_eulerian_dnodal_coordinates(), oomph::FiniteElement::dshape_local(), i, k, oomph::FiniteElement::local_to_eulerian_mapping(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), p, Eigen::numext::q, and u_index_axi_nst().

◆ interpolated_duds_axi_nst()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::interpolated_duds_axi_nst ( const Vector< double > &  s,
const unsigned i,
const unsigned j 
) const
inline

Return FE interpolated derivatives of velocity component u[i] w.r.t spatial local coordinate direction s[j] at local coordinate s

1105  {
1106  // Determine number of nodes
1107  const unsigned n_node = nnode();
1108 
1109  // Allocate storage for local shape function and its derivatives
1110  // with respect to space
1111  Shape psif(n_node);
1112  DShape dpsifds(n_node, 2);
1113 
1114  // Find values of shape function (ignore jacobian)
1115  (void)this->dshape_local(s, psif, dpsifds);
1116 
1117  // Get the index at which the velocity is stored
1118  const unsigned u_nodal_index = u_index_axi_nst(i);
1119 
1120  // Initialise value of duds
1121  double interpolated_duds = 0.0;
1122 
1123  // Loop over the local nodes and sum
1124  for (unsigned l = 0; l < n_node; l++)
1125  {
1126  interpolated_duds += nodal_value(l, u_nodal_index) * dpsifds(l, j);
1127  }
1128 
1129  return (interpolated_duds);
1130  }

References oomph::FiniteElement::dshape_local(), i, j, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), and u_index_axi_nst().

◆ interpolated_dudt_axi_nst()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::interpolated_dudt_axi_nst ( const Vector< double > &  s,
const unsigned i 
) const
inline

Return FE interpolated derivatives of velocity component u[i] w.r.t time at local coordinate s

1168  {
1169  // Determine number of nodes
1170  const unsigned n_node = nnode();
1171 
1172  // Allocate storage for local shape function
1173  Shape psif(n_node);
1174 
1175  // Find values of shape function
1176  shape(s, psif);
1177 
1178  // Initialise value of dudt
1179  double interpolated_dudt = 0.0;
1180 
1181  // Loop over the local nodes and sum
1182  for (unsigned l = 0; l < n_node; l++)
1183  {
1184  interpolated_dudt += du_dt_axi_nst(l, i) * psif[l];
1185  }
1186 
1187  return (interpolated_dudt);
1188  }

References du_dt_axi_nst(), i, oomph::FiniteElement::nnode(), s, and oomph::FiniteElement::shape().

◆ interpolated_dudx_axi_nst()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::interpolated_dudx_axi_nst ( const Vector< double > &  s,
const unsigned i,
const unsigned j 
) const
inline

Return FE interpolated derivatives of velocity component u[i] w.r.t spatial global coordinate direction x[j] at local coordinate s

1137  {
1138  // Determine number of nodes
1139  const unsigned n_node = nnode();
1140 
1141  // Allocate storage for local shape function and its derivatives
1142  // with respect to space
1143  Shape psif(n_node);
1144  DShape dpsifdx(n_node, 2);
1145 
1146  // Find values of shape function (ignore jacobian)
1147  (void)this->dshape_eulerian(s, psif, dpsifdx);
1148 
1149  // Get the index at which the velocity is stored
1150  const unsigned u_nodal_index = u_index_axi_nst(i);
1151 
1152  // Initialise value of dudx
1153  double interpolated_dudx = 0.0;
1154 
1155  // Loop over the local nodes and sum
1156  for (unsigned l = 0; l < n_node; l++)
1157  {
1158  interpolated_dudx += nodal_value(l, u_nodal_index) * dpsifdx(l, j);
1159  }
1160 
1161  return (interpolated_dudx);
1162  }
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Definition: elements.cc:3298

References oomph::FiniteElement::dshape_eulerian(), i, j, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), and u_index_axi_nst().

◆ interpolated_p_axi_nst()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::interpolated_p_axi_nst ( const Vector< double > &  s) const
inline

Return FE interpolated pressure at local coordinate s.

1081  {
1082  // Find number of nodes
1083  unsigned n_pres = npres_axi_nst();
1084  // Local shape function
1085  Shape psi(n_pres);
1086  // Find values of shape function
1087  pshape_axi_nst(s, psi);
1088 
1089  // Initialise value of p
1090  double interpolated_p = 0.0;
1091  // Loop over the local nodes and sum
1092  for (unsigned l = 0; l < n_pres; l++)
1093  {
1094  interpolated_p += p_axi_nst(l) * psi[l];
1095  }
1096 
1097  return (interpolated_p);
1098  }

References npres_axi_nst(), p_axi_nst(), pshape_axi_nst(), and s.

Referenced by oomph::RefineableGeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement::further_build(), oomph::RefineableGeneralisedNewtonianAxisymmetricQTaylorHoodElement::get_interpolated_values(), output(), point_output_data(), pressure_integral(), scalar_value_paraview(), and traction().

◆ interpolated_u_axi_nst() [1/3]

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::interpolated_u_axi_nst ( const unsigned t,
const Vector< double > &  s,
const unsigned i 
) const
inline

Return FE interpolated velocity u[i] at local coordinate s.

999  {
1000  // Find number of nodes
1001  unsigned n_node = nnode();
1002  // Local shape function
1003  Shape psi(n_node);
1004  // Find values of shape function
1005  shape(s, psi);
1006 
1007  // Get the index at which the velocity is stored
1008  unsigned u_nodal_index = u_index_axi_nst(i);
1009 
1010  // Initialise value of u
1011  double interpolated_u = 0.0;
1012  // Loop over the local nodes and sum
1013  for (unsigned l = 0; l < n_node; l++)
1014  {
1015  interpolated_u += nodal_value(t, l, u_nodal_index) * psi[l];
1016  }
1017 
1018  return (interpolated_u);
1019  }

References i, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), s, oomph::FiniteElement::shape(), plotPSD::t, and u_index_axi_nst().

◆ interpolated_u_axi_nst() [2/3]

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::interpolated_u_axi_nst ( const Vector< double > &  s,
const unsigned i 
) const
inline

Return FE interpolated velocity u[i] at local coordinate s.

971  {
972  // Find number of nodes
973  unsigned n_node = nnode();
974  // Local shape function
975  Shape psi(n_node);
976  // Find values of shape function
977  shape(s, psi);
978 
979  // Get the index at which the velocity is stored
980  unsigned u_nodal_index = u_index_axi_nst(i);
981 
982  // Initialise value of u
983  double interpolated_u = 0.0;
984  // Loop over the local nodes and sum
985  for (unsigned l = 0; l < n_node; l++)
986  {
987  interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
988  }
989 
990  return (interpolated_u);
991  }

References i, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), s, oomph::FiniteElement::shape(), and u_index_axi_nst().

◆ interpolated_u_axi_nst() [3/3]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::interpolated_u_axi_nst ( const Vector< double > &  s,
Vector< double > &  veloc 
) const
inline

Compute vector of FE interpolated velocity u at local coordinate s.

946  {
947  // Find number of nodes
948  unsigned n_node = nnode();
949  // Local shape function
950  Shape psi(n_node);
951  // Find values of shape function
952  shape(s, psi);
953 
954  for (unsigned i = 0; i < 3; i++)
955  {
956  // Index at which the nodal value is stored
957  unsigned u_nodal_index = u_index_axi_nst(i);
958  // Initialise value of u
959  veloc[i] = 0.0;
960  // Loop over the local nodes and sum
961  for (unsigned l = 0; l < n_node; l++)
962  {
963  veloc[i] += nodal_value(l, u_nodal_index) * psi[l];
964  }
965  }
966  }

References i, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), s, oomph::FiniteElement::shape(), and u_index_axi_nst().

Referenced by compute_error(), oomph::RefineableGeneralisedNewtonianAxisymmetricQTaylorHoodElement::get_interpolated_values(), oomph::RefineableGeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement::get_interpolated_values(), kin_energy(), output(), point_output_data(), and scalar_value_paraview().

◆ kin_energy()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::kin_energy ( ) const

Get integral of kinetic energy over element.

Get integral of kinetic energy over element:

1027  {
1028  throw OomphLibError(
1029  "Check the kinetic energy calculation for axisymmetric NSt",
1032 
1033  // Initialise
1034  double kin_en = 0.0;
1035 
1036  // Set the value of Nintpt
1037  unsigned Nintpt = integral_pt()->nweight();
1038 
1039  // Set the Vector to hold local coordinates
1040  Vector<double> s(2);
1041 
1042  // Loop over the integration points
1043  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
1044  {
1045  // Assign values of s
1046  for (unsigned i = 0; i < 2; i++)
1047  {
1048  s[i] = integral_pt()->knot(ipt, i);
1049  }
1050 
1051  // Get the integral weight
1052  double w = integral_pt()->weight(ipt);
1053 
1054  // Get Jacobian of mapping
1055  double J = J_eulerian(s);
1056 
1057  // Loop over directions
1058  double veloc_squared = 0.0;
1059  for (unsigned i = 0; i < 3; i++)
1060  {
1061  veloc_squared +=
1063  }
1064 
1065  kin_en += 0.5 * veloc_squared * w * J * interpolated_x(s, 0);
1066  }
1067 
1068  return kin_en;
1069  }

References i, oomph::FiniteElement::integral_pt(), interpolated_u_axi_nst(), oomph::FiniteElement::interpolated_x(), J, oomph::FiniteElement::J_eulerian(), oomph::Integral::knot(), GlobalParameters::Nintpt, oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, w, and oomph::Integral::weight().

◆ max_and_min_invariant_and_viscosity()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::max_and_min_invariant_and_viscosity ( double min_invariant,
double max_invariant,
double min_viscosity,
double max_viscosity 
) const

Get max. and min. strain rate invariant and visocosity over all integration points in element

1123  {
1124  // Initialise
1125  min_invariant = DBL_MAX;
1126  max_invariant = -DBL_MAX;
1127  min_viscosity = DBL_MAX;
1128  max_viscosity = -DBL_MAX;
1129 
1130  // Number of integration points
1131  unsigned Nintpt = integral_pt()->nweight();
1132 
1133  // Set the Vector to hold local coordinates (two dimensions)
1134  Vector<double> s(2);
1135 
1136  // Loop over the integration points
1137  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
1138  {
1139  // Assign values of s
1140  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
1141 
1142  // The strainrate
1143  DenseMatrix<double> strainrate(3, 3, 0.0);
1144  strain_rate(s, strainrate);
1145 
1146  // Calculate the second invariant
1147  double second_invariant =
1149 
1150  // Get the viscosity according to the constitutive equation
1151  double viscosity = Constitutive_eqn_pt->viscosity(second_invariant);
1152 
1153  min_invariant = std::min(min_invariant, second_invariant);
1154  max_invariant = std::max(max_invariant, second_invariant);
1155  min_viscosity = std::min(min_viscosity, viscosity);
1156  max_viscosity = std::max(max_viscosity, viscosity);
1157  }
1158  }
#define min(a, b)
Definition: datatypes.h:22
#define max(a, b)
Definition: datatypes.h:23

References Constitutive_eqn_pt, i, oomph::FiniteElement::integral_pt(), oomph::Integral::knot(), max, min, GlobalParameters::Nintpt, oomph::Integral::nweight(), s, oomph::SecondInvariantHelper::second_invariant(), strain_rate(), and oomph::GeneralisedNewtonianConstitutiveEquation< DIM >::viscosity().

◆ n_u_nst()

unsigned oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::n_u_nst ( ) const
inline

Return the number of velocity components for use in general FluidInterface class

564  {
565  return 3;
566  }

◆ npres_axi_nst()

◆ nscalar_paraview()

unsigned oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::nscalar_paraview ( ) const
inlinevirtual

Number of scalars/fields output by this element. Reimplements broken virtual function in base class.

Reimplemented from oomph::FiniteElement.

694  {
695  return 4;
696  }

◆ output() [1/4]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::output ( FILE *  file_pt)
inlinevirtual

Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::PseudoSolidNodeUpdateElement< GeneralisedNewtonianAxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >, oomph::GeneralisedNewtonianAxisymmetricTTaylorHoodElement, oomph::GeneralisedNewtonianAxisymmetricTCrouzeixRaviartElement, oomph::GeneralisedNewtonianAxisymmetricQTaylorHoodElement, and oomph::GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement.

803  {
804  unsigned nplot = 5;
805  output(file_pt, nplot);
806  }
void output(std::ostream &outfile)
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:789

References output().

◆ output() [2/4]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::output ( FILE *  file_pt,
const unsigned nplot 
)
virtual

Output function: x,y,[z],u,v,[w],p in tecplot format. nplot points in each coordinate direction

Output function: r,z,u,v,w,p in tecplot format. Specified number of plot points in each coordinate direction.

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::GeneralisedNewtonianAxisymmetricTTaylorHoodElement, oomph::GeneralisedNewtonianAxisymmetricTCrouzeixRaviartElement, oomph::GeneralisedNewtonianAxisymmetricQTaylorHoodElement, oomph::GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement, and oomph::PseudoSolidNodeUpdateElement< GeneralisedNewtonianAxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >.

555  {
556  // Vector of local coordinates
557  Vector<double> s(2);
558 
559  // Tecplot header info
560  fprintf(file_pt, "%s ", tecplot_zone_string(nplot).c_str());
561 
562  // Loop over plot points
563  unsigned num_plot_points = nplot_points(nplot);
564  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
565  {
566  // Get local coordinates of plot point
567  get_s_plot(iplot, nplot, s);
568 
569  // Coordinates
570  for (unsigned i = 0; i < 2; i++)
571  {
572  // outfile << interpolated_x(s,i) << " ";
573  fprintf(file_pt, "%g ", interpolated_x(s, i));
574  }
575 
576  // Velocities
577  for (unsigned i = 0; i < 3; i++)
578  {
579  // outfile << interpolated_u(s,i) << " ";
580  fprintf(file_pt, "%g ", interpolated_u_axi_nst(s, i));
581  }
582 
583  // Pressure
584  // outfile << interpolated_p(s) << " ";
585  fprintf(file_pt, "%g ", interpolated_p_axi_nst(s));
586 
587  // outfile << std::endl;
588  fprintf(file_pt, "\n");
589  }
590  // outfile << std::endl;
591  fprintf(file_pt, "\n");
592 
593  // Write tecplot footer (e.g. FE connectivity lists)
594  write_tecplot_zone_footer(file_pt, nplot);
595  }
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Definition: elements.h:3161
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Definition: elements.h:3186
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Definition: elements.h:3174
double interpolated_p_axi_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
Definition: generalised_newtonian_axisym_navier_stokes_elements.h:1080

References oomph::FiniteElement::get_s_plot(), i, interpolated_p_axi_nst(), interpolated_u_axi_nst(), oomph::FiniteElement::interpolated_x(), oomph::FiniteElement::nplot_points(), s, oomph::FiniteElement::tecplot_zone_string(), and oomph::FiniteElement::write_tecplot_zone_footer().

◆ output() [3/4]

◆ output() [4/4]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::output ( std::ostream &  outfile,
const unsigned nplot 
)
virtual

Output function: x,y,[z],u,v,[w],p in tecplot format. nplot points in each coordinate direction

Output function: r,z,u,v,w,p in tecplot format. Specified number of plot points in each coordinate direction.

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::GeneralisedNewtonianAxisymmetricTTaylorHoodElement, oomph::GeneralisedNewtonianAxisymmetricTCrouzeixRaviartElement, oomph::GeneralisedNewtonianAxisymmetricQTaylorHoodElement, oomph::GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement, and oomph::PseudoSolidNodeUpdateElement< GeneralisedNewtonianAxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >.

509  {
510  // Vector of local coordinates
511  Vector<double> s(2);
512 
513  // Tecplot header info
514  outfile << tecplot_zone_string(nplot);
515 
516  // Loop over plot points
517  unsigned num_plot_points = nplot_points(nplot);
518  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
519  {
520  // Get local coordinates of plot point
521  get_s_plot(iplot, nplot, s);
522 
523  // Coordinates
524  for (unsigned i = 0; i < 2; i++)
525  {
526  outfile << interpolated_x(s, i) << " ";
527  }
528 
529  // Velocities
530  for (unsigned i = 0; i < 3; i++)
531  {
532  outfile << interpolated_u_axi_nst(s, i) << " ";
533  }
534 
535  // Pressure
536  outfile << interpolated_p_axi_nst(s) << " ";
537 
538  outfile << std::endl;
539  }
540  outfile << std::endl;
541 
542  // Write tecplot footer (e.g. FE connectivity lists)
543  write_tecplot_zone_footer(outfile, nplot);
544  }

References oomph::FiniteElement::get_s_plot(), i, interpolated_p_axi_nst(), interpolated_u_axi_nst(), oomph::FiniteElement::interpolated_x(), oomph::FiniteElement::nplot_points(), s, oomph::FiniteElement::tecplot_zone_string(), and oomph::FiniteElement::write_tecplot_zone_footer().

◆ output_fct() [1/2]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::output_fct ( std::ostream &  outfile,
const unsigned nplot,
const double time,
FiniteElement::UnsteadyExactSolutionFctPt  exact_soln_pt 
)
virtual

Output exact solution specified via function pointer at a given time and at a given number of plot points. Function prints as many components as are returned in solution Vector.

Output "exact" solution at a given time Solution is provided via function pointer. Plot at a given number of plot points. Function prints as many components as are returned in solution Vector.

Reimplemented from oomph::FiniteElement.

378  {
379  // Vector of local coordinates
380  Vector<double> s(2);
381 
382  // Vector for coordintes
383  Vector<double> x(2);
384 
385  // Tecplot header info
386  outfile << tecplot_zone_string(nplot);
387 
388  // Exact solution Vector
389  Vector<double> exact_soln;
390 
391  // Loop over plot points
392  unsigned num_plot_points = nplot_points(nplot);
393  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
394  {
395  // Get local coordinates of plot point
396  get_s_plot(iplot, nplot, s);
397 
398  // Get x position as Vector
399  interpolated_x(s, x);
400 
401  // Get exact solution at this point
402  (*exact_soln_pt)(time, x, exact_soln);
403 
404  // Output x,y,...
405  for (unsigned i = 0; i < 2; i++)
406  {
407  outfile << x[i] << " ";
408  }
409 
410  // Output "exact solution"
411  for (unsigned i = 0; i < exact_soln.size(); i++)
412  {
413  outfile << exact_soln[i] << " ";
414  }
415 
416  outfile << std::endl;
417  }
418 
419  // Write tecplot footer (e.g. FE connectivity lists)
420  write_tecplot_zone_footer(outfile, nplot);
421  }

References ProblemParameters::exact_soln(), oomph::FiniteElement::get_s_plot(), i, oomph::FiniteElement::interpolated_x(), oomph::FiniteElement::nplot_points(), s, oomph::FiniteElement::tecplot_zone_string(), oomph::FiniteElement::write_tecplot_zone_footer(), and plotDoE::x.

◆ output_fct() [2/2]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::output_fct ( std::ostream &  outfile,
const unsigned nplot,
FiniteElement::SteadyExactSolutionFctPt  exact_soln_pt 
)
virtual

Output exact solution specified via function pointer at a given number of plot points. Function prints as many components as are returned in solution Vector

Output "exact" solution Solution is provided via function pointer. Plot at a given number of plot points. Function prints as many components as are returned in solution Vector.

Reimplemented from oomph::FiniteElement.

322  {
323  // Vector of local coordinates
324  Vector<double> s(2);
325 
326  // Vector for coordintes
327  Vector<double> x(2);
328 
329  // Tecplot header info
330  outfile << tecplot_zone_string(nplot);
331 
332  // Exact solution Vector
333  Vector<double> exact_soln;
334 
335  // Loop over plot points
336  unsigned num_plot_points = nplot_points(nplot);
337  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
338  {
339  // Get local coordinates of plot point
340  get_s_plot(iplot, nplot, s);
341 
342  // Get x position as Vector
343  interpolated_x(s, x);
344 
345  // Get exact solution at this point
346  (*exact_soln_pt)(x, exact_soln);
347 
348  // Output x,y,...
349  for (unsigned i = 0; i < 2; i++)
350  {
351  outfile << x[i] << " ";
352  }
353 
354  // Output "exact solution"
355  for (unsigned i = 0; i < exact_soln.size(); i++)
356  {
357  outfile << exact_soln[i] << " ";
358  }
359 
360  outfile << std::endl;
361  }
362 
363  // Write tecplot footer (e.g. FE connectivity lists)
364  write_tecplot_zone_footer(outfile, nplot);
365  }

References ProblemParameters::exact_soln(), oomph::FiniteElement::get_s_plot(), i, oomph::FiniteElement::interpolated_x(), oomph::FiniteElement::nplot_points(), s, oomph::FiniteElement::tecplot_zone_string(), oomph::FiniteElement::write_tecplot_zone_footer(), and plotDoE::x.

◆ output_veloc()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::output_veloc ( std::ostream &  outfile,
const unsigned nplot,
const unsigned t 
)

Output function: x,y,[z],u,v,[w] in tecplot format. nplot points in each coordinate direction at timestep t (t=0: present; t>0: previous timestep)

Output function: Velocities only x,y,[z],u,v,[w] in tecplot format at specified previous timestep (t=0: present; t>0: previous timestep). Specified number of plot points in each coordinate direction.

432  {
433  // Find number of nodes
434  unsigned n_node = nnode();
435 
436  // Local shape function
437  Shape psi(n_node);
438 
439  // Vectors of local coordinates and coords and velocities
440  Vector<double> s(2);
441  Vector<double> interpolated_x(2);
442  Vector<double> interpolated_u(3);
443 
444 
445  // Tecplot header info
446  outfile << tecplot_zone_string(nplot);
447 
448  // Loop over plot points
449  unsigned num_plot_points = nplot_points(nplot);
450  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
451  {
452  // Get local coordinates of plot point
453  get_s_plot(iplot, nplot, s);
454 
455  // Get shape functions
456  shape(s, psi);
457 
458  // Loop over coordinate directions
459  for (unsigned i = 0; i < 2; i++)
460  {
461  interpolated_x[i] = 0.0;
462  // Loop over the local nodes and sum
463  for (unsigned l = 0; l < n_node; l++)
464  {
465  interpolated_x[i] += nodal_position(t, l, i) * psi[l];
466  }
467  }
468 
469  // Loop over the velocity components
470  for (unsigned i = 0; i < 3; i++)
471  {
472  // Get the index at which the velocity is stored
473  unsigned u_nodal_index = u_index_axi_nst(i);
474  interpolated_u[i] = 0.0;
475  // Loop over the local nodes and sum
476  for (unsigned l = 0; l < n_node; l++)
477  {
478  interpolated_u[i] += nodal_value(t, l, u_nodal_index) * psi[l];
479  }
480  }
481 
482  // Coordinates
483  for (unsigned i = 0; i < 2; i++)
484  {
485  outfile << interpolated_x[i] << " ";
486  }
487 
488  // Velocities
489  for (unsigned i = 0; i < 3; i++)
490  {
491  outfile << interpolated_u[i] << " ";
492  }
493 
494  outfile << std::endl;
495  }
496 
497  // Write tecplot footer (e.g. FE connectivity lists)
498  write_tecplot_zone_footer(outfile, nplot);
499  }

References oomph::FiniteElement::get_s_plot(), i, oomph::FiniteElement::interpolated_x(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::FiniteElement::nplot_points(), s, oomph::FiniteElement::shape(), plotPSD::t, oomph::FiniteElement::tecplot_zone_string(), u_index_axi_nst(), and oomph::FiniteElement::write_tecplot_zone_footer().

◆ p_axi_nst()

◆ p_local_eqn()

◆ p_nodal_index_axi_nst()

virtual int oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::p_nodal_index_axi_nst ( ) const
inlinevirtual

◆ point_output_data()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::point_output_data ( const Vector< double > &  s,
Vector< double > &  data 
)
inlinevirtual

Output solution in data vector at local cordinates s: r,z,u_r,u_z,u_phi,p

Reimplemented from oomph::FiniteElement.

769  {
770  // Output the components of the position
771  for (unsigned i = 0; i < 2; i++)
772  {
773  data.push_back(interpolated_x(s, i));
774  }
775 
776  // Output the components of the FE representation of u at s
777  for (unsigned i = 0; i < 3; i++)
778  {
779  data.push_back(interpolated_u_axi_nst(s, i));
780  }
781 
782  // Output FE representation of p at s
783  data.push_back(interpolated_p_axi_nst(s));
784  }
int data[]
Definition: Map_placement_new.cpp:1

References data, i, interpolated_p_axi_nst(), interpolated_u_axi_nst(), oomph::FiniteElement::interpolated_x(), and s.

◆ pressure_integral()

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::pressure_integral ( ) const

Integral of pressure over element.

Return pressure integrated over the element.

1076  {
1077  // Initialise
1078  double press_int = 0;
1079 
1080  // Set the value of Nintpt
1081  unsigned Nintpt = integral_pt()->nweight();
1082 
1083  // Set the Vector to hold local coordinates
1084  Vector<double> s(2);
1085 
1086  // Loop over the integration points
1087  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
1088  {
1089  // Assign values of s
1090  for (unsigned i = 0; i < 2; i++)
1091  {
1092  s[i] = integral_pt()->knot(ipt, i);
1093  }
1094 
1095  // Get the integral weight
1096  double w = integral_pt()->weight(ipt);
1097 
1098  // Get Jacobian of mapping
1099  double J = J_eulerian(s);
1100 
1101  // Premultiply the weights and the Jacobian
1102  double W = w * J * interpolated_x(s, 0);
1103 
1104  // Get pressure
1105  double press = interpolated_p_axi_nst(s);
1106 
1107  // Add
1108  press_int += press * W;
1109  }
1110 
1111  return press_int;
1112  }

References i, oomph::FiniteElement::integral_pt(), interpolated_p_axi_nst(), oomph::FiniteElement::interpolated_x(), J, oomph::FiniteElement::J_eulerian(), oomph::Integral::knot(), GlobalParameters::Nintpt, oomph::Integral::nweight(), s, w, oomph::QuadTreeNames::W, and oomph::Integral::weight().

◆ pshape_axi_nst() [1/2]

◆ pshape_axi_nst() [2/2]

virtual void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::pshape_axi_nst ( const Vector< double > &  s,
Shape psi,
Shape test 
) const
protectedpure virtual

◆ re()

◆ re_invfr()

◆ re_invfr_pt()

double*& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invfr_pt ( )
inline

Pointer to global inverse Froude number.

451  {
452  return ReInvFr_pt;
453  }

References ReInvFr_pt.

Referenced by fill_in_generic_dresidual_contribution_axi_nst(), and oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::further_build().

◆ re_invro()

const double& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invro ( ) const
inline

◆ re_invro_pt()

double*& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_invro_pt ( )
inline

Pointer to global inverse Froude number.

463  {
464  return ReInvRo_pt;
465  }

References ReInvRo_pt.

Referenced by fill_in_generic_dresidual_contribution_axi_nst(), and oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::further_build().

◆ re_pt()

double*& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_pt ( )
inline

◆ re_st()

const double& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_st ( ) const
inline

◆ re_st_pt()

double*& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::re_st_pt ( )
inline

Pointer to product of Reynolds and Strouhal number (=Womersley number)

439  {
440  return ReSt_pt;
441  }

References ReSt_pt.

Referenced by fill_in_generic_dresidual_contribution_axi_nst(), and oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::further_build().

◆ scalar_name_paraview()

std::string oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::scalar_name_paraview ( const unsigned i) const
inlinevirtual

Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.

Reimplemented from oomph::FiniteElement.

742  {
743  // Veloc
744  if (i < 3)
745  {
746  return "Velocity " + StringConversion::to_string(i);
747  }
748  // Pressure field
749  else if (i == 3)
750  {
751  return "Pressure";
752  }
753  // Never get here
754  else
755  {
756  std::stringstream error_stream;
757  error_stream
758  << "Axisymmetric Navier-Stokes Elements only store 4 fields "
759  << std::endl;
760  throw OomphLibError(
761  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
762  return " ";
763  }
764  }
std::string to_string(T object, unsigned float_precision=8)
Definition: oomph_utilities.h:189

References i, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::StringConversion::to_string().

◆ scalar_value_paraview()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::scalar_value_paraview ( std::ofstream &  file_out,
const unsigned i,
const unsigned nplot 
) const
inlinevirtual

Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specific element type.

Reimplemented from oomph::FiniteElement.

703  {
704  // Vector of local coordinates
705  Vector<double> s(2);
706 
707  // Loop over plot points
708  unsigned num_plot_points = nplot_points_paraview(nplot);
709  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
710  {
711  // Get local coordinates of plot point
712  get_s_plot(iplot, nplot, s);
713 
714  // Velocities
715  if (i < 3)
716  {
717  file_out << interpolated_u_axi_nst(s, i) << std::endl;
718  }
719  // Pressure
720  else if (i == 3)
721  {
722  file_out << interpolated_p_axi_nst(s) << std::endl;
723  }
724  // Never get here
725  else
726  {
727  std::stringstream error_stream;
728  error_stream
729  << "Axisymmetric Navier-Stokes Elements only store 4 fields "
730  << std::endl;
731  throw OomphLibError(error_stream.str(),
734  }
735  }
736  }
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Definition: elements.h:2862

References oomph::FiniteElement::get_s_plot(), i, interpolated_p_axi_nst(), interpolated_u_axi_nst(), oomph::FiniteElement::nplot_points_paraview(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

◆ strain_rate() [1/2]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::strain_rate ( const unsigned t,
const Vector< double > &  s,
DenseMatrix< double > &  strainrate 
) const

Strain-rate tensor: \( e_{ij} \) where \( i,j = r,z,\theta \) (in that order) at a specific time history value

Get strain-rate tensor: \( e_{ij} \) where \( i,j = r,z,\theta \) (in that order) at a specific time

854  {
855 #ifdef PARANOID
856  if ((strainrate.ncol() != 3) || (strainrate.nrow() != 3))
857  {
858  std::ostringstream error_message;
859  error_message << "The strain rate has incorrect dimensions "
860  << strainrate.ncol() << " x " << strainrate.nrow()
861  << " Not 3" << std::endl;
862 
863  throw OomphLibError(
864  error_message.str(),
865  "GeneralisedNewtonianAxisymmetricNavierStokeEquations::strain_rate()",
867  }
868 #endif
869 
870  // Find out how many nodes there are in the element
871  unsigned n_node = nnode();
872 
873  // Set up memory for the shape and test functions
874  Shape psi(n_node);
875  DShape dpsidx(n_node, 2);
876 
877  // Loop over all nodes to back up current positions and over-write them
878  // with the appropriate history values
879  DenseMatrix<double> backed_up_nodal_position(n_node, 2);
880  for (unsigned j = 0; j < n_node; j++)
881  {
882  backed_up_nodal_position(j, 0) = node_pt(j)->x(0);
883  node_pt(j)->x(0) = node_pt(j)->x(t, 0);
884  backed_up_nodal_position(j, 1) = node_pt(j)->x(1);
885  node_pt(j)->x(1) = node_pt(j)->x(t, 1);
886  }
887 
888  // Call the derivatives of the shape functions
889  dshape_eulerian(s, psi, dpsidx);
890 
891  // Radius
892  double interpolated_r = 0.0;
893  double interpolated_z = 0.0;
894 
895  // Velocity components and their derivatives
896  double ur = 0.0;
897  double durdr = 0.0;
898  double durdz = 0.0;
899  double uz = 0.0;
900  double duzdr = 0.0;
901  double duzdz = 0.0;
902  double uphi = 0.0;
903  double duphidr = 0.0;
904  double duphidz = 0.0;
905 
906  // Get the local storage for the indices
907  unsigned u_nodal_index[3];
908  for (unsigned i = 0; i < 3; ++i)
909  {
910  u_nodal_index[i] = u_index_axi_nst(i);
911  }
912 
913  // Loop over nodes to assemble velocities and their derivatives
914  // w.r.t. to r and z (x_0 and x_1)
915  for (unsigned l = 0; l < n_node; l++)
916  {
917  interpolated_r += nodal_position(l, 0) * psi[l];
918  interpolated_z += nodal_position(l, 1) * psi[l];
919 
920  ur += nodal_value(t, l, u_nodal_index[0]) * psi[l];
921  uz += nodal_value(t, l, u_nodal_index[1]) * psi[l];
922  uphi += nodal_value(t, l, u_nodal_index[2]) * psi[l];
923 
924  durdr += nodal_value(t, l, u_nodal_index[0]) * dpsidx(l, 0);
925  durdz += nodal_value(t, l, u_nodal_index[0]) * dpsidx(l, 1);
926 
927  duzdr += nodal_value(t, l, u_nodal_index[1]) * dpsidx(l, 0);
928  duzdz += nodal_value(t, l, u_nodal_index[1]) * dpsidx(l, 1);
929 
930  duphidr += nodal_value(t, l, u_nodal_index[2]) * dpsidx(l, 0);
931  duphidz += nodal_value(t, l, u_nodal_index[2]) * dpsidx(l, 1);
932  }
933 
934 
935  // Assign strain rates without negative powers of the radius
936  // and zero those with:
937  strainrate(0, 0) = durdr;
938  strainrate(0, 1) = 0.5 * (durdz + duzdr);
939  strainrate(1, 0) = strainrate(0, 1);
940  strainrate(0, 2) = 0.0;
941  strainrate(2, 0) = strainrate(0, 2);
942  strainrate(1, 1) = duzdz;
943  strainrate(1, 2) = 0.5 * duphidz;
944  strainrate(2, 1) = strainrate(1, 2);
945  strainrate(2, 2) = 0.0;
946 
947 
948  // Overwrite the strain rates with negative powers of the radius
949  // unless we're at the origin
950  if (std::fabs(interpolated_r) > 1.0e-16)
951  {
952  double inverse_radius = 1.0 / interpolated_r;
953  strainrate(0, 2) = 0.5 * (duphidr - inverse_radius * uphi);
954  strainrate(2, 0) = strainrate(0, 2);
955  strainrate(2, 2) = inverse_radius * ur;
956  }
957  else
958  {
959  // hierher L'Hospital? --- also for strain(0,2)
960  strainrate(2, 2) = durdr;
961  }
962 
963  // Reset current nodal positions
964  for (unsigned j = 0; j < n_node; j++)
965  {
966  node_pt(j)->x(0) = backed_up_nodal_position(j, 0);
967  node_pt(j)->x(1) = backed_up_nodal_position(j, 1);
968  }
969  }
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117

References oomph::FiniteElement::dshape_eulerian(), boost::multiprecision::fabs(), i, j, oomph::DenseMatrix< T >::ncol(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::FiniteElement::node_pt(), oomph::DenseMatrix< T >::nrow(), OOMPH_EXCEPTION_LOCATION, s, plotPSD::t, u_index_axi_nst(), and oomph::Node::x().

◆ strain_rate() [2/2]

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::strain_rate ( const Vector< double > &  s,
DenseMatrix< double > &  strainrate 
) const

Strain-rate tensor: \( e_{ij} \) where \( i,j = r,z,\theta \) (in that order)

Get strain-rate tensor: \( e_{ij} \) where \( i,j = r,z,\theta \) (in that order)

751  {
752 #ifdef PARANOID
753  if ((strainrate.ncol() != 3) || (strainrate.nrow() != 3))
754  {
755  std::ostringstream error_message;
756  error_message << "The strain rate has incorrect dimensions "
757  << strainrate.ncol() << " x " << strainrate.nrow()
758  << " Not 3" << std::endl;
759 
760  throw OomphLibError(
761  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
762  }
763 #endif
764 
765  // Find out how many nodes there are in the element
766  unsigned n_node = nnode();
767 
768  // Set up memory for the shape and test functions
769  Shape psi(n_node);
770  DShape dpsidx(n_node, 2);
771 
772  // Call the derivatives of the shape functions
773  dshape_eulerian(s, psi, dpsidx);
774 
775  // Radius
776  double interpolated_r = 0.0;
777 
778  // Velocity components and their derivatives
779  double ur = 0.0;
780  double durdr = 0.0;
781  double durdz = 0.0;
782  double uz = 0.0;
783  double duzdr = 0.0;
784  double duzdz = 0.0;
785  double uphi = 0.0;
786  double duphidr = 0.0;
787  double duphidz = 0.0;
788 
789  // Get the local storage for the indices
790  unsigned u_nodal_index[3];
791  for (unsigned i = 0; i < 3; ++i)
792  {
793  u_nodal_index[i] = u_index_axi_nst(i);
794  }
795 
796  // Loop over nodes to assemble velocities and their derivatives
797  // w.r.t. to r and z (x_0 and x_1)
798  for (unsigned l = 0; l < n_node; l++)
799  {
800  interpolated_r += nodal_position(l, 0) * psi[l];
801 
802  ur += nodal_value(l, u_nodal_index[0]) * psi[l];
803  uz += nodal_value(l, u_nodal_index[1]) * psi[l];
804  uphi += nodal_value(l, u_nodal_index[2]) * psi[l];
805 
806  durdr += nodal_value(l, u_nodal_index[0]) * dpsidx(l, 0);
807  durdz += nodal_value(l, u_nodal_index[0]) * dpsidx(l, 1);
808 
809  duzdr += nodal_value(l, u_nodal_index[1]) * dpsidx(l, 0);
810  duzdz += nodal_value(l, u_nodal_index[1]) * dpsidx(l, 1);
811 
812  duphidr += nodal_value(l, u_nodal_index[2]) * dpsidx(l, 0);
813  duphidz += nodal_value(l, u_nodal_index[2]) * dpsidx(l, 1);
814  }
815 
816 
817  // Assign strain rates without negative powers of the radius
818  // and zero those with:
819  strainrate(0, 0) = durdr;
820  strainrate(0, 1) = 0.5 * (durdz + duzdr);
821  strainrate(1, 0) = strainrate(0, 1);
822  strainrate(0, 2) = 0.0;
823  strainrate(2, 0) = strainrate(0, 2);
824  strainrate(1, 1) = duzdz;
825  strainrate(1, 2) = 0.5 * duphidz;
826  strainrate(2, 1) = strainrate(1, 2);
827  strainrate(2, 2) = 0.0;
828 
829 
830  // Overwrite the strain rates with negative powers of the radius
831  // unless we're at the origin
832  if (std::fabs(interpolated_r) > 1.0e-16)
833  {
834  double inverse_radius = 1.0 / interpolated_r;
835  strainrate(0, 2) = 0.5 * (duphidr - inverse_radius * uphi);
836  strainrate(2, 0) = strainrate(0, 2);
837  strainrate(2, 2) = inverse_radius * ur;
838  }
839  else
840  {
841  // hierher L'Hospital? -- also for (2,0)
842  strainrate(2, 2) = durdr;
843  }
844  }

References oomph::FiniteElement::dshape_eulerian(), boost::multiprecision::fabs(), i, oomph::DenseMatrix< T >::ncol(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, and u_index_axi_nst().

Referenced by dissipation(), extrapolated_strain_rate(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::get_Z2_flux(), oomph::GeneralisedNewtonianAxisymmetricTCrouzeixRaviartElement::get_Z2_flux(), oomph::GeneralisedNewtonianAxisymmetricTTaylorHoodElement::get_Z2_flux(), max_and_min_invariant_and_viscosity(), and traction().

◆ traction()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::traction ( const Vector< double > &  s,
const Vector< double > &  N,
Vector< double > &  traction 
) const

Compute traction (on the viscous scale) at local coordinate s for outer unit normal N

661  {
662  // throw OomphLibError(
663  // "Check the traction calculation for axisymmetric NSt",
664  // OOMPH_CURRENT_FUNCTION,
665  // OOMPH_EXCEPTION_LOCATION);
666 
667  // Pad out normal vector if required
668  Vector<double> n_local(3, 0.0);
669  n_local[0] = N[0];
670  n_local[1] = N[1];
671 
672 #ifdef PARANOID
673  if ((N.size() == 3) && (N[2] != 0.0))
674  {
675  throw OomphLibError(
676  "Unit normal passed into this fct should either be 2D (r,z) or have a "
677  "zero component in the theta-direction",
680  }
681 #endif
682 
683  // Get velocity gradients
684  DenseMatrix<double> strainrate(3, 3, 0.0);
685 
686  // Do we use the current or extrapolated strainrate to compute
687  // the second invariant?
689  {
690  strain_rate(s, strainrate);
691  }
692  else
693  {
694  extrapolated_strain_rate(s, strainrate);
695  }
696 
697  // Get the second invariant of the rate of strain tensor
698  double second_invariant =
700 
702 
703  // Get pressure
704  double press = interpolated_p_axi_nst(s);
705 
706  // Loop over traction components
707  for (unsigned i = 0; i < 3; i++)
708  {
709  traction[i] = -press * n_local[i];
710  for (unsigned j = 0; j < 3; j++)
711  {
712  traction[i] += visc * 2.0 * strainrate(i, j) * n_local[j];
713  }
714  }
715  }
void traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Definition: generalised_newtonian_axisym_navier_stokes_elements.cc:657
@ N
Definition: constructor.cpp:22

References Constitutive_eqn_pt, extrapolated_strain_rate(), i, interpolated_p_axi_nst(), j, N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, oomph::SecondInvariantHelper::second_invariant(), strain_rate(), Use_extrapolated_strainrate_to_compute_second_invariant, and oomph::GeneralisedNewtonianConstitutiveEquation< DIM >::viscosity().

Referenced by oomph::GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement::get_traction(), and oomph::GeneralisedNewtonianAxisymmetricQTaylorHoodElement::get_traction().

◆ u_index_axi_nst()

virtual unsigned oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::u_index_axi_nst ( const unsigned i) const
inlinevirtual

◆ u_index_nst()

unsigned oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::u_index_nst ( const unsigned i) const
inline

Return the index at which the i-th unknown velocity component is stored with a common interface for use in general FluidInterface and similar elements. To do: Merge all common storage etc to a common base class for Navier–Stokes elements in all coordinate systems.

557  {
558  return this->u_index_axi_nst(i);
559  }

References u_index_axi_nst().

◆ use_current_strainrate_to_compute_second_invariant()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::use_current_strainrate_to_compute_second_invariant ( )
inline

Switch to fully implict evaluation of non-Newtonian const eqn.

References Use_extrapolated_strainrate_to_compute_second_invariant.

◆ use_extrapolated_strainrate_to_compute_second_invariant()

void oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::use_extrapolated_strainrate_to_compute_second_invariant ( )
inline

Use extrapolation for non-Newtonian const eqn.

References Use_extrapolated_strainrate_to_compute_second_invariant.

◆ viscosity_ratio()

const double& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::viscosity_ratio ( ) const
inline

◆ viscosity_ratio_pt()

double*& oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::viscosity_ratio_pt ( )
inline

Pointer to Viscosity Ratio.

501  {
502  return Viscosity_Ratio_pt;
503  }

References Viscosity_Ratio_pt.

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

Member Data Documentation

◆ ALE_is_disabled

◆ axi_nst_body_force_fct_pt

void(*&)(const double& time, const Vector<double>& x, Vector<double>& f) oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::axi_nst_body_force_fct_pt()
inline

Access function for the body-force pointer.

509  {
510  return Body_force_fct_pt;
511  }

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

◆ Body_force_fct_pt

void(* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Body_force_fct_pt) (const double &time, const Vector< double > &x, Vector< double > &result)
protected

◆ Constitutive_eqn_pt

GeneralisedNewtonianConstitutiveEquation<3>* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Constitutive_eqn_pt
protected

◆ Default_Gravity_vector

Vector< double > oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Default_Gravity_vector
staticprivate

Static default value for the gravity vector.

Navier-Stokes equations default gravity vector.

Referenced by GeneralisedNewtonianAxisymmetricNavierStokesEquations().

◆ Default_Physical_Constant_Value

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Default_Physical_Constant_Value
staticprivate

Static default value for the physical constants (all initialised to zero)

Referenced by GeneralisedNewtonianAxisymmetricNavierStokesEquations().

◆ Default_Physical_Ratio_Value

double oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Default_Physical_Ratio_Value
staticprivate

Static default value for the physical ratios (all are initialised to one)

Referenced by GeneralisedNewtonianAxisymmetricNavierStokesEquations().

◆ Density_Ratio_pt

double* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Density_Ratio_pt
protected

Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)

Referenced by density_ratio(), density_ratio_pt(), oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::further_build(), and GeneralisedNewtonianAxisymmetricNavierStokesEquations().

◆ G_pt

Vector<double>* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::G_pt
protected

◆ Gamma

Vector< double > oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Gamma
static

◆ Pre_multiply_by_viscosity_ratio

bool oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Pre_multiply_by_viscosity_ratio
staticprotected

Static boolean telling us whether we pre-multiply the pressure and continuity by the viscosity ratio

Referenced by fill_in_generic_residual_contribution_axi_nst().

◆ Pressure_not_stored_at_node

int oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Pressure_not_stored_at_node
staticprivate

Static "magic" number that indicates that the pressure is not stored at a node

Referenced by p_nodal_index_axi_nst().

◆ Re_pt

double* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Re_pt
protected

◆ ReInvFr_pt

double* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ReInvFr_pt
protected

Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)

Referenced by oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::further_build(), GeneralisedNewtonianAxisymmetricNavierStokesEquations(), re_invfr(), and re_invfr_pt().

◆ ReInvRo_pt

double* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ReInvRo_pt
protected

Pointer to global Reynolds number x inverse Rossby number (used when in a rotating frame)

Referenced by oomph::RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::further_build(), GeneralisedNewtonianAxisymmetricNavierStokesEquations(), re_invro(), and re_invro_pt().

◆ ReSt_pt

double* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::ReSt_pt
protected

◆ source_fct_pt

double(*&)(const double& time, const Vector<double>& x) oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::source_fct_pt()
inline

Access function for the source-function pointer.

515  {
516  return Source_fct_pt;
517  }

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

◆ Source_fct_pt

double(* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Source_fct_pt) (const double &time, const Vector< double > &x)
protected

◆ Use_extrapolated_strainrate_to_compute_second_invariant

◆ Viscosity_Ratio_pt

double* oomph::GeneralisedNewtonianAxisymmetricNavierStokesEquations::Viscosity_Ratio_pt
protected

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