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

#include <axisym_navier_stokes_elements.h>

+ Inheritance diagram for oomph::AxisymmetricNavierStokesEquations:

Public Member Functions

 AxisymmetricNavierStokesEquations ()
 Constructor: NULL the body force and source function. More...
 
const double & re () const
 Reynolds number. More...
 
const double & re_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
const double & re_invfr () const
 Global inverse Froude number. More...
 
double *& re_invfr_pt ()
 Pointer to global inverse Froude number. More...
 
const double & re_invro () const
 Global Reynolds number multiplied by inverse Rossby number. More...
 
double *& re_invro_pt ()
 Pointer to global inverse Froude number. More...
 
const Vector< double > & g () const
 Vector of gravitational components. More...
 
Vector< double > *& g_pt ()
 Pointer to Vector of gravitational components. More...
 
const double & density_ratio () const
 
double *& density_ratio_pt ()
 Pointer to Density ratio. More...
 
const double & viscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to Viscosity Ratio. More...
 
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...
 
double dissipation () const
 Return integral of dissipation over element. More...
 
double dissipation (const Vector< double > &s) const
 Return dissipation at local coordinate s. More...
 
double kin_energy () const
 Get integral of kinetic energy over element. More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 Strain-rate tensor: \( e_{ij} \) where \( i,j = r,z,\theta \) (in that order) More...
 
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)
 
MacroElement * macro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void 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 Node * vertex_node_pt (const unsigned &j) const
 
virtual Node * construct_node (const unsigned &n)
 
virtual Node * construct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual Node * construct_boundary_node (const unsigned &n)
 
virtual Node * construct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
int get_node_number (Node *const &node_pt) const
 
virtual Node * get_node_at_local_coordinate (const Vector< double > &s) const
 
double raw_nodal_value (const unsigned &n, const unsigned &i) const
 
double raw_nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
unsigned dim () const
 
virtual ElementGeometry::ElementGeometry element_geometry () const
 Return the geometry type of the element (either Q or T usually). More...
 
virtual double interpolated_x (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated coordinate x[i] at local coordinate s. More...
 
virtual double interpolated_x (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_x (const Vector< double > &s, Vector< double > &x) const
 Return FE interpolated position x[] at local coordinate s as Vector. More...
 
virtual void interpolated_x (const unsigned &t, const Vector< double > &s, Vector< double > &x) const
 
virtual double interpolated_dxdt (const Vector< double > &s, const unsigned &i, const unsigned &t)
 
virtual void interpolated_dxdt (const Vector< double > &s, const unsigned &t, Vector< double > &dxdt)
 
unsigned ngeom_data () const
 
Data * geom_data_pt (const unsigned &j)
 
void position (const Vector< double > &zeta, Vector< double > &r) const
 
void position (const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
 
void dposition_dt (const Vector< double > &zeta, const unsigned &t, Vector< double > &drdt)
 
virtual double zeta_nodal (const unsigned &n, const unsigned &k, const unsigned &i) const
 
void interpolated_zeta (const Vector< double > &s, Vector< double > &zeta) const
 
void locate_zeta (const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
 
virtual void node_update ()
 
virtual void identify_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 ()
 
TimeStepper * time_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 
- Public Member Functions inherited from oomph::NavierStokesElementWithDiagonalMassMatrices
 NavierStokesElementWithDiagonalMassMatrices ()
 Empty constructor. More...
 
virtual ~NavierStokesElementWithDiagonalMassMatrices ()
 Virtual destructor. More...
 
 NavierStokesElementWithDiagonalMassMatrices (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken copy constructor. More...
 
void operator= (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken assignment operator. More...
 

Public Attributes

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

Static Public Attributes

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

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

double * Viscosity_Ratio_pt
 
double * Density_Ratio_pt
 
double * Re_pt
 Pointer to global Reynolds number. More...
 
double * ReSt_pt
 Pointer to global Reynolds number x Strouhal number (=Womersley) More...
 
double * ReInvFr_pt
 
double * ReInvRo_pt
 
Vector< double > * G_pt
 Pointer to global gravity Vector. More...
 
void(* Body_force_fct_pt )(const double &time, const Vector< double > &x, Vector< double > &result)
 Pointer to body force function. More...
 
double(* Source_fct_pt )(const double &time, const Vector< double > &x)
 Pointer to volumetric source function. More...
 
bool ALE_is_disabled
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElement * Macro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepper * Geom_object_time_stepper_pt
 

Static Private Attributes

static int Pressure_not_stored_at_node = -100
 
static double Default_Physical_Constant_Value
 Navier–Stokes equations static data. More...
 
static double Default_Physical_Ratio_Value = 1.0
 
static Vector< double > Default_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 > &)
 
- Static Protected Attributes inherited from oomph::FiniteElement
static const unsigned Default_Initial_Nvalue = 0
 Default value for the number of values at a node. More...
 
static const double Node_location_tolerance = 1.0e-14
 
static const unsigned N2deriv [] = {0, 1, 3, 6}
 
- Static Protected Attributes inherited from oomph::GeneralisedElement
static DenseMatrix< double > Dummy_matrix
 
static std::deque< double * > Dof_pt_deque
 

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

◆ AxisymmetricNavierStokesEquations()

oomph::AxisymmetricNavierStokesEquations::AxisymmetricNavierStokesEquations ( )
inline

Constructor: NULL the body force and source function.

376  : Body_force_fct_pt(0), Source_fct_pt(0), ALE_is_disabled(false)
377  {
378  // Set all the Physical parameter pointers to the default value zero
379  Re_pt = &Default_Physical_Constant_Value;
380  ReSt_pt = &Default_Physical_Constant_Value;
381  ReInvFr_pt = &Default_Physical_Constant_Value;
382  ReInvRo_pt = &Default_Physical_Constant_Value;
383  G_pt = &Default_Gravity_vector;
384  // Set the Physical ratios to the default value of 1
385  Viscosity_Ratio_pt = &Default_Physical_Ratio_Value;
386  Density_Ratio_pt = &Default_Physical_Ratio_Value;
387  }
oomph::AxisymmetricNavierStokesEquations::Default_Gravity_vector
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
Definition: axisym_navier_stokes_elements.h:130
oomph::AxisymmetricNavierStokesEquations::Source_fct_pt
double(* Source_fct_pt)(const double &time, const Vector< double > &x)
Pointer to volumetric source function.
Definition: axisym_navier_stokes_elements.h:168
oomph::AxisymmetricNavierStokesEquations::G_pt
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: axisym_navier_stokes_elements.h:160
oomph::AxisymmetricNavierStokesEquations::ReInvFr_pt
double * ReInvFr_pt
Definition: axisym_navier_stokes_elements.h:153
oomph::AxisymmetricNavierStokesEquations::Re_pt
double * Re_pt
Pointer to global Reynolds number.
Definition: axisym_navier_stokes_elements.h:146
oomph::AxisymmetricNavierStokesEquations::ReSt_pt
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: axisym_navier_stokes_elements.h:149
oomph::AxisymmetricNavierStokesEquations::Density_Ratio_pt
double * Density_Ratio_pt
Definition: axisym_navier_stokes_elements.h:141
oomph::AxisymmetricNavierStokesEquations::Default_Physical_Ratio_Value
static double Default_Physical_Ratio_Value
Definition: axisym_navier_stokes_elements.h:127
oomph::AxisymmetricNavierStokesEquations::Body_force_fct_pt
void(* Body_force_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to body force function.
Definition: axisym_navier_stokes_elements.h:163
oomph::AxisymmetricNavierStokesEquations::Default_Physical_Constant_Value
static double Default_Physical_Constant_Value
Navier–Stokes equations static data.
Definition: axisym_navier_stokes_elements.h:123
oomph::AxisymmetricNavierStokesEquations::ALE_is_disabled
bool ALE_is_disabled
Definition: axisym_navier_stokes_elements.h:173
oomph::AxisymmetricNavierStokesEquations::ReInvRo_pt
double * ReInvRo_pt
Definition: axisym_navier_stokes_elements.h:157
oomph::AxisymmetricNavierStokesEquations::Viscosity_Ratio_pt
double * Viscosity_Ratio_pt
Definition: axisym_navier_stokes_elements.h:137

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::AxisymmetricNavierStokesEquations::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.

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

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

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::AxisymmetricNavierStokesEquations::compute_physical_size ( ) const
inlinevirtual

Compute the volume of the element.

Reimplemented from oomph::FiniteElement.

1173  {
1174  // Initialise result
1175  double result = 0.0;
1176 
1177  // Figure out the global (Eulerian) spatial dimension of the
1178  // element by checking the Eulerian dimension of the nodes
1179  const unsigned n_dim_eulerian = nodal_dimension();
1180 
1181  // Allocate Vector of global Eulerian coordinates
1182  Vector<double> x(n_dim_eulerian);
1183 
1184  // Set the value of n_intpt
1185  const unsigned n_intpt = integral_pt()->nweight();
1186 
1187  // Vector of local coordinates
1188  const unsigned n_dim = dim();
1189  Vector<double> s(n_dim);
1190 
1191  // Loop over the integration points
1192  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1193  {
1194  // Assign the values of s
1195  for (unsigned i = 0; i < n_dim; i++)
1196  {
1197  s[i] = integral_pt()->knot(ipt, i);
1198  }
1199 
1200  // Assign the values of the global Eulerian coordinates
1201  for (unsigned i = 0; i < n_dim_eulerian; i++)
1202  {
1203  x[i] = interpolated_x(s, i);
1204  }
1205 
1206  // Get the integral weight
1207  const double w = integral_pt()->weight(ipt);
1208 
1209  // Get Jacobian of mapping
1210  const double J = J_eulerian(s);
1211 
1212  // The integrand at the current integration point is r
1213  result += x[0] * w * J;
1214  }
1215 
1216  // Multiply by 2pi (integrating in azimuthal direction)
1217  return (2.0 * MathematicalConstants::Pi * result);
1218  }
oomph::FiniteElement::dim
unsigned dim() const
Definition: elements.h:2611
oomph::FiniteElement::nodal_dimension
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
oomph::MathematicalConstants::Pi
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.

◆ density_ratio()

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

Density ratio for element: Element's density relative to the viscosity used in the definition of the Reynolds number

460  {
461  return *Density_Ratio_pt;
462  }

References Density_Ratio_pt.

Referenced by fill_in_contribution_to_hessian_vector_products(), fill_in_generic_dresidual_contribution_axi_nst(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ density_ratio_pt()

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

Pointer to Density ratio.

466  {
467  return Density_Ratio_pt;
468  }

References Density_Ratio_pt.

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

◆ dinterpolated_u_axi_nst_ddata()

virtual void oomph::AxisymmetricNavierStokesEquations::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::RefineableAxisymmetricNavierStokesEquations.

953  {
954  // Find number of nodes
955  unsigned n_node = nnode();
956  // Local shape function
957  Shape psi(n_node);
958  // Find values of shape function
959  shape(s, psi);
960 
961  // Find the index at which the velocity component is stored
962  const unsigned u_nodal_index = u_index_axi_nst(i);
963 
964  // Find the number of dofs associated with interpolated u
965  unsigned n_u_dof = 0;
966  for (unsigned l = 0; l < n_node; l++)
967  {
968  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
969  // If it's positive add to the count
970  if (global_eqn >= 0)
971  {
972  ++n_u_dof;
973  }
974  }
975 
976  // Now resize the storage schemes
977  du_ddata.resize(n_u_dof, 0.0);
978  global_eqn_number.resize(n_u_dof, 0);
979 
980  // Loop over th nodes again and set the derivatives
981  unsigned count = 0;
982  // Loop over the local nodes and sum
983  for (unsigned l = 0; l < n_node; l++)
984  {
985  // Get the global equation number
986  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
987  if (global_eqn >= 0)
988  {
989  // Set the global equation number
990  global_eqn_number[count] = global_eqn;
991  // Set the derivative with respect to the unknown
992  du_ddata[count] = psi[l];
993  // Increase the counter
994  ++count;
995  }
996  }
997  }
oomph::AxisymmetricNavierStokesEquations::u_index_axi_nst
virtual unsigned u_index_axi_nst(const unsigned &i) const
Definition: axisym_navier_stokes_elements.h:506
oomph::Data::eqn_number
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::shape
virtual void shape(const Vector< double > &s, Shape &psi) const =0
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210

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

Referenced by QAxisymAdvectionDiffusionElementWithExternalElement::get_dwind_axi_adv_diff_dexternal_element_data().

◆ disable_ALE()

void oomph::AxisymmetricNavierStokesEquations::disable_ALE ( )
inlinevirtual

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

Reimplemented from oomph::FiniteElement.

561  {
562  ALE_is_disabled = true;
563  }

References ALE_is_disabled.

Referenced by oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::disable_ALE(), and oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::disable_ALE().

◆ dissipation() [1/2]

double oomph::AxisymmetricNavierStokesEquations::dissipation ( ) const

Return integral of dissipation over element.

595  {
596  throw OomphLibError(
597  "Check the dissipation calculation for axisymmetric NSt",
598  OOMPH_CURRENT_FUNCTION,
599  OOMPH_EXCEPTION_LOCATION);
600 
601  // Initialise
602  double diss = 0.0;
603 
604  // Set the value of Nintpt
605  unsigned Nintpt = integral_pt()->nweight();
606 
607  // Set the Vector to hold local coordinates
608  Vector<double> s(2);
609 
610  // Loop over the integration points
611  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
612  {
613  // Assign values of s
614  for (unsigned i = 0; i < 2; i++)
615  {
616  s[i] = integral_pt()->knot(ipt, i);
617  }
618 
619  // Get the integral weight
620  double w = integral_pt()->weight(ipt);
621 
622  // Get Jacobian of mapping
623  double J = J_eulerian(s);
624 
625  // Get strain rate matrix
626  DenseMatrix<double> strainrate(3, 3);
627  strain_rate(s, strainrate);
628 
629  // Initialise
630  double local_diss = 0.0;
631  for (unsigned i = 0; i < 3; i++)
632  {
633  for (unsigned j = 0; j < 3; j++)
634  {
635  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
636  }
637  }
638 
639  diss += local_diss * w * J;
640  }
641 
642  return diss;
643  }
oomph::AxisymmetricNavierStokesEquations::strain_rate
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
Definition: axisym_navier_stokes_elements.cc:725
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
j
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::AxisymmetricNavierStokesEquations::dissipation ( const Vector< double > &  s) const

Return dissipation at local coordinate s.

698  {
699  throw OomphLibError(
700  "Check the dissipation calculation for axisymmetric NSt",
701  OOMPH_CURRENT_FUNCTION,
702  OOMPH_EXCEPTION_LOCATION);
703 
704  // Get strain rate matrix
705  DenseMatrix<double> strainrate(3, 3);
706  strain_rate(s, strainrate);
707 
708  // Initialise
709  double local_diss = 0.0;
710  for (unsigned i = 0; i < 3; i++)
711  {
712  for (unsigned j = 0; j < 3; j++)
713  {
714  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
715  }
716  }
717 
718  return local_diss;
719  }

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::AxisymmetricNavierStokesEquations::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::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

◆ dshape_and_dtest_eulerian_at_knot_axi_nst() [2/2]

virtual double oomph::AxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_axi_nst ( const unsigned &  ipt,
Shape &  psi,
DShape &  dpsidx,
Shape &  test,
DShape &  dtestdx 
) const
protectedpure virtual

Compute the shape functions and derivatives w.r.t. global coords at ipt-th integration point Return Jacobian of mapping between local and global coordinates.

Implemented in oomph::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

Referenced by fill_in_contribution_to_hessian_vector_products(), fill_in_generic_dresidual_contribution_axi_nst(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ dshape_and_dtest_eulerian_axi_nst()

virtual double oomph::AxisymmetricNavierStokesEquations::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::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

◆ du_dt_axi_nst()

double oomph::AxisymmetricNavierStokesEquations::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.

532  {
533  // Get the data's timestepper
534  TimeStepper* time_stepper_pt = this->node_pt(n)->time_stepper_pt();
535 
536  // Initialise dudt
537  double dudt = 0.0;
538  // Loop over the timesteps, if there is a non Steady timestepper
539  if (!time_stepper_pt->is_steady())
540  {
541  // Get the index at which the velocity is stored
542  const unsigned u_nodal_index = u_index_axi_nst(i);
543 
544  // Number of timsteps (past & present)
545  const unsigned n_time = time_stepper_pt->ntstorage();
546 
547  // Add the contributions to the time derivative
548  for (unsigned t = 0; t < n_time; t++)
549  {
550  dudt +=
551  time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
552  }
553  }
554 
555  return dudt;
556  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::FiniteElement::nodal_value
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
oomph::GeomObject::time_stepper_pt
TimeStepper *& time_stepper_pt()
Definition: geom_objects.h:192
oomph::TimeStepper::ntstorage
unsigned ntstorage() const
Definition: timesteppers.h:601
oomph::TimeStepper::weight
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
oomph::TimeStepper::is_steady
bool is_steady() const
Definition: timesteppers.h:389
plotPSD.t
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::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates(), and interpolated_dudt_axi_nst().

◆ enable_ALE()

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

570  {
571  ALE_is_disabled = false;
572  }

References ALE_is_disabled.

Referenced by oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::enable_ALE(), and oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::enable_ALE().

◆ fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter()

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

857  {
858  // Call the generic routine with the flag set to 2
859  fill_in_generic_dresidual_contribution_axi_nst(
860  parameter_pt, dres_dparam, djac_dparam, dmass_matrix_dparam, 2);
861  }
oomph::AxisymmetricNavierStokesEquations::fill_in_generic_dresidual_contribution_axi_nst
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: axisym_navier_stokes_elements.cc:2510

References fill_in_generic_dresidual_contribution_axi_nst().

◆ fill_in_contribution_to_djacobian_dparameter()

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

840  {
841  // Call the generic routine with the flag set to 1
842  fill_in_generic_dresidual_contribution_axi_nst(
843  parameter_pt,
844  dres_dparam,
845  djac_dparam,
846  GeneralisedElement::Dummy_matrix,
847  1);
848  }
oomph::GeneralisedElement::Dummy_matrix
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::AxisymmetricNavierStokesEquations::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.

823  {
824  // Call the generic residuals function with flag set to 0
825  // and using a dummy matrix argument
826  fill_in_generic_dresidual_contribution_axi_nst(
827  parameter_pt,
828  dres_dparam,
829  GeneralisedElement::Dummy_matrix,
830  GeneralisedElement::Dummy_matrix,
831  0);
832  }

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

◆ fill_in_contribution_to_hessian_vector_products()

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

3243  {
3244  // Find out how many nodes there are
3245  unsigned n_node = nnode();
3246 
3247  // Get the nodal indices at which the velocity is stored
3248  unsigned u_nodal_index[3];
3249  for (unsigned i = 0; i < 3; ++i)
3250  {
3251  u_nodal_index[i] = u_index_axi_nst(i);
3252  }
3253 
3254  // Set up memory for the shape and test functions
3255  // Note that there are only two dimensions, r and z in this problem
3256  Shape psif(n_node), testf(n_node);
3257  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
3258 
3259  // Number of integration points
3260  unsigned Nintpt = integral_pt()->nweight();
3261 
3262  // Set the Vector to hold local coordinates (two dimensions)
3263  Vector<double> s(2);
3264 
3265  // Get Physical Variables from Element
3266  // Reynolds number must be multiplied by the density ratio
3267  double scaled_re = re() * density_ratio();
3268  // double visc_ratio = viscosity_ratio();
3269  Vector<double> G = g();
3270 
3271  // Integers used to store the local equation and unknown numbers
3272  int local_eqn = 0, local_unknown = 0, local_freedom = 0;
3273 
3274  // How may dofs are there
3275  const unsigned n_dof = this->ndof();
3276 
3277  // Create a local matrix eigenvector product contribution
3278  DenseMatrix<double> jac_y(n_dof, n_dof, 0.0);
3279 
3280  // Loop over the integration points
3281  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
3282  {
3283  // Assign values of s
3284  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
3285  // Get the integral weight
3286  double w = integral_pt()->weight(ipt);
3287 
3288  // Call the derivatives of the shape and test functions
3289  double J = dshape_and_dtest_eulerian_at_knot_axi_nst(
3290  ipt, psif, dpsifdx, testf, dtestfdx);
3291 
3292  // Premultiply the weights and the Jacobian
3293  double W = w * J;
3294 
3295  // Allocate storage for the position and the derivative of the
3296  // mesh positions wrt time
3297  Vector<double> interpolated_x(2, 0.0);
3298  // Vector<double> mesh_velocity(2,0.0);
3299  // Allocate storage for the pressure, velocity components and their
3300  // time and spatial derivatives
3301  Vector<double> interpolated_u(3, 0.0);
3302  // Vector<double> dudt(3,0.0);
3303  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
3304 
3305  // Calculate velocities and derivatives at integration point
3306 
3307  // Loop over nodes
3308  for (unsigned l = 0; l < n_node; l++)
3309  {
3310  // Cache the shape function
3311  const double psif_ = psif(l);
3312  // Loop over the two coordinate directions
3313  for (unsigned i = 0; i < 2; i++)
3314  {
3315  interpolated_x[i] += this->raw_nodal_position(l, i) * psif_;
3316  }
3317 
3318  // Loop over the three velocity directions
3319  for (unsigned i = 0; i < 3; i++)
3320  {
3321  // Get the u_value
3322  const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
3323  interpolated_u[i] += u_value * psif_;
3324  // dudt[i]+= du_dt_axi_nst(l,i)*psif_;
3325  // Loop over derivative directions
3326  for (unsigned j = 0; j < 2; j++)
3327  {
3328  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
3329  }
3330  }
3331  }
3332 
3333  // Get the mesh velocity if ALE is enabled
3334  if (!ALE_is_disabled)
3335  {
3336  throw OomphLibError("Moving nodes not implemented\n",
3337  OOMPH_CURRENT_FUNCTION,
3338  OOMPH_EXCEPTION_LOCATION);
3339  }
3340 
3341  // r is the first position component
3342  double r = interpolated_x[0];
3343 
3344 
3345  // MOMENTUM EQUATIONS
3346  //------------------
3347 
3348  // Loop over the test functions
3349  for (unsigned l = 0; l < n_node; l++)
3350  {
3351  // FIRST (RADIAL) MOMENTUM EQUATION
3352  local_eqn = nodal_local_eqn(l, u_nodal_index[0]);
3353  // If it's not a boundary condition
3354  if (local_eqn >= 0)
3355  {
3356  // Loop over the velocity shape functions yet again
3357  for (unsigned l3 = 0; l3 < n_node; l3++)
3358  {
3359  // Derivative of jacobian terms with respect to radial velocity
3360  local_freedom = nodal_local_eqn(l3, u_nodal_index[0]);
3361  if (local_freedom >= 0)
3362  {
3363  // Storage for the sums
3364  double temp = 0.0;
3365 
3366  // Loop over the velocity shape functions again
3367  for (unsigned l2 = 0; l2 < n_node; l2++)
3368  {
3369  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
3370  // Radial velocity component
3371  if (local_unknown >= 0)
3372  {
3373  // Remains of convective terms
3374  temp -= scaled_re *
3375  (r * psif[l2] * dpsifdx(l3, 0) +
3376  r * psif[l3] * dpsifdx(l2, 0)) *
3377  Y[local_unknown] * testf[l] * W;
3378  }
3379 
3380  // Axial velocity component
3381  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
3382  if (local_unknown >= 0)
3383  {
3384  // Convective terms
3385  temp -= scaled_re * r * psif[l2] * dpsifdx(l3, 1) *
3386  Y[local_unknown] * testf[l] * W;
3387  }
3388  }
3389  // Add the summed contribution to the product matrix
3390  jac_y(local_eqn, local_freedom) += temp;
3391  } // End of derivative wrt radial coordinate
3392 
3393 
3394  // Derivative of jacobian terms with respect to axial velocity
3395  local_freedom = nodal_local_eqn(l3, u_nodal_index[1]);
3396  if (local_freedom >= 0)
3397  {
3398  double temp = 0.0;
3399 
3400  // Loop over the velocity shape functions again
3401  for (unsigned l2 = 0; l2 < n_node; l2++)
3402  {
3403  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
3404  // Radial velocity component
3405  if (local_unknown >= 0)
3406  {
3407  // Remains of convective terms
3408  temp -= scaled_re * (r * psif[l3] * dpsifdx(l2, 1)) *
3409  Y[local_unknown] * testf[l] * W;
3410  }
3411  }
3412  // Add the summed contribution to the product matrix
3413  jac_y(local_eqn, local_freedom) += temp;
3414  } // End of derivative wrt axial coordinate
3415 
3416  // Derivative of jacobian terms with respect to azimuthal velocity
3417  local_freedom = nodal_local_eqn(l3, u_nodal_index[2]);
3418  if (local_freedom >= 0)
3419  {
3420  double temp = 0.0;
3421 
3422  // Loop over the velocity shape functions again
3423  for (unsigned l2 = 0; l2 < n_node; l2++)
3424  {
3425  // Azimuthal velocity component
3426  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
3427  if (local_unknown >= 0)
3428  {
3429  // Convective terms
3430  temp -= -scaled_re * 2.0 * psif[l3] * psif[l2] *
3431  Y[local_unknown] * testf[l] * W;
3432  }
3433  }
3434  // Add the summed contibution
3435  jac_y(local_eqn, local_freedom) += temp;
3436 
3437  } // End of if not boundary condition statement
3438  } // End of loop over freedoms
3439  } // End of RADIAL MOMENTUM EQUATION
3440 
3441 
3442  // SECOND (AXIAL) MOMENTUM EQUATION
3443  local_eqn = nodal_local_eqn(l, u_nodal_index[1]);
3444  // If it's not a boundary condition
3445  if (local_eqn >= 0)
3446  {
3447  // Loop over the velocity shape functions yet again
3448  for (unsigned l3 = 0; l3 < n_node; l3++)
3449  {
3450  // Derivative of jacobian terms with respect to radial velocity
3451  local_freedom = nodal_local_eqn(l3, u_nodal_index[0]);
3452  if (local_freedom >= 0)
3453  {
3454  double temp = 0.0;
3455 
3456  // Loop over the velocity shape functions again
3457  for (unsigned l2 = 0; l2 < n_node; l2++)
3458  {
3459  // Axial velocity component
3460  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
3461  if (local_unknown >= 0)
3462  {
3463  // Convective terms
3464  temp -= scaled_re * (r * psif[l3] * dpsifdx(l2, 0)) *
3465  Y[local_unknown] * testf[l] * W;
3466  }
3467  }
3468  jac_y(local_eqn, local_freedom) += temp;
3469 
3470  // There are no azimithal terms in the axial momentum equation
3471  } // End of loop over velocity shape functions
3472 
3473 
3474  // Derivative of jacobian terms with respect to axial velocity
3475  local_freedom = nodal_local_eqn(l3, u_nodal_index[1]);
3476  if (local_freedom >= 0)
3477  {
3478  double temp = 0.0;
3479 
3480  // Loop over the velocity shape functions again
3481  for (unsigned l2 = 0; l2 < n_node; l2++)
3482  {
3483  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
3484  // Radial velocity component
3485  if (local_unknown >= 0)
3486  {
3487  // Convective terms
3488  temp -= scaled_re * r * psif[l2] * dpsifdx(l3, 0) *
3489  Y[local_unknown] * testf[l] * W;
3490  }
3491 
3492  // Axial velocity component
3493  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
3494  if (local_unknown >= 0)
3495  {
3496  // Convective terms
3497  temp -= scaled_re *
3498  (r * psif[l2] * dpsifdx(l3, 1) +
3499  r * psif[l3] * dpsifdx(l2, 1)) *
3500  Y[local_unknown] * testf[l] * W;
3501  }
3502 
3503  // There are no azimithal terms in the axial momentum equation
3504  } // End of loop over velocity shape functions
3505 
3506  // Add summed contributiont to jacobian product matrix
3507  jac_y(local_eqn, local_freedom) += temp;
3508  }
3509  } // End of loop over local freedoms
3510 
3511  } // End of AXIAL MOMENTUM EQUATION
3512 
3513  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
3514  local_eqn = nodal_local_eqn(l, u_nodal_index[2]);
3515  if (local_eqn >= 0)
3516  {
3517  // Loop over the velocity shape functions yet again
3518  for (unsigned l3 = 0; l3 < n_node; l3++)
3519  {
3520  // Derivative of jacobian terms with respect to radial velocity
3521  local_freedom = nodal_local_eqn(l3, u_nodal_index[0]);
3522  if (local_freedom >= 0)
3523  {
3524  double temp = 0.0;
3525 
3526  // Loop over the velocity shape functions again
3527  for (unsigned l2 = 0; l2 < n_node; l2++)
3528  {
3529  // Azimuthal velocity component
3530  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
3531  if (local_unknown >= 0)
3532  {
3533  // Convective terms
3534  temp -=
3535  scaled_re *
3536  (r * psif[l3] * dpsifdx(l2, 0) + psif[l3] * psif[l2]) *
3537  Y[local_unknown] * testf[l] * W;
3538  }
3539  }
3540  // Add the summed contribution to the jacobian eigenvector sum
3541  jac_y(local_eqn, local_freedom) += temp;
3542  }
3543 
3544  // Derivative of jacobian terms with respect to axial velocity
3545  local_freedom = nodal_local_eqn(l3, u_nodal_index[1]);
3546  if (local_freedom >= 0)
3547  {
3548  double temp = 0.0;
3549 
3550  // Loop over the velocity shape functions again
3551  for (unsigned l2 = 0; l2 < n_node; l2++)
3552  {
3553  // Azimuthal velocity component
3554  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
3555  if (local_unknown >= 0)
3556  {
3557  // Convective terms
3558  temp -= scaled_re * (r * psif[l3] * dpsifdx(l2, 1)) *
3559  Y[local_unknown] * testf[l] * W;
3560  }
3561  }
3562  // Add the summed contribution to the jacobian eigenvector sum
3563  jac_y(local_eqn, local_freedom) += temp;
3564  }
3565 
3566 
3567  // Derivative of jacobian terms with respect to azimuthal velocity
3568  local_freedom = nodal_local_eqn(l3, u_nodal_index[2]);
3569  if (local_freedom >= 0)
3570  {
3571  double temp = 0.0;
3572 
3573 
3574  // Loop over the velocity shape functions again
3575  for (unsigned l2 = 0; l2 < n_node; l2++)
3576  {
3577  // Radial velocity component
3578  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
3579  if (local_unknown >= 0)
3580  {
3581  // Convective terms
3582  temp -=
3583  scaled_re *
3584  (r * psif[l2] * dpsifdx(l3, 0) + psif[l2] * psif[l3]) *
3585  Y[local_unknown] * testf[l] * W;
3586  }
3587 
3588  // Axial velocity component
3589  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
3590  if (local_unknown >= 0)
3591  {
3592  // Convective terms
3593  temp -= scaled_re * r * psif[l2] * dpsifdx(l3, 1) *
3594  Y[local_unknown] * testf[l] * W;
3595  }
3596  }
3597  // Add the fully summed contribution
3598  jac_y(local_eqn, local_freedom) += temp;
3599  }
3600  } // End of loop over freedoms
3601 
3602  // There are no pressure terms
3603  } // End of AZIMUTHAL EQUATION
3604 
3605  } // End of loop over shape functions
3606  }
3607 
3608  // We have now assembled the matrix (d J_{ij} Y_j)/d u_{k}
3609  // and simply need to sum over the vectors
3610  const unsigned n_vec = C.nrow();
3611  for (unsigned i = 0; i < n_dof; i++)
3612  {
3613  for (unsigned k = 0; k < n_dof; k++)
3614  {
3615  // Cache the value of the hessian y product
3616  const double j_y = jac_y(i, k);
3617  // Loop over the possible vectors
3618  for (unsigned v = 0; v < n_vec; v++)
3619  {
3620  product(v, i) += j_y * C(v, k);
3621  }
3622  }
3623  }
3624  }
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
G
JacobiRotation< float > G
Definition: Jacobi_makeGivens.cpp:2
C
Matrix< Scalar, Dynamic, Dynamic > C
Definition: bench_gemm.cpp:49
oomph::AxisymmetricNavierStokesEquations::density_ratio
const double & density_ratio() const
Definition: axisym_navier_stokes_elements.h:459
oomph::AxisymmetricNavierStokesEquations::g
const Vector< double > & g() const
Vector of gravitational components.
Definition: axisym_navier_stokes_elements.h:446
oomph::AxisymmetricNavierStokesEquations::re
const double & re() const
Reynolds number.
Definition: axisym_navier_stokes_elements.h:398
oomph::AxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_axi_nst
virtual double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
oomph::FiniteElement::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
oomph::FiniteElement::raw_nodal_value
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2576
oomph::FiniteElement::raw_nodal_position
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.cc:1686
oomph::GeneralisedElement::ndof
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
oomph::Matrix
Definition: matrices.h:74
oomph::Matrix::nrow
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
k
char char char int int * k
Definition: level2_impl.h:374
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20
product
void product(const MatrixType &m)
Definition: product.h:42
Y
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::AxisymmetricNavierStokesEquations::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< AxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >.

795  {
796  // Call the generic routine with the flag set to 1
797  fill_in_generic_residual_contribution_axi_nst(
798  residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
799  }
oomph::AxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst
virtual void fill_in_generic_residual_contribution_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Definition: axisym_navier_stokes_elements.cc:914

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

Referenced by oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::fill_in_contribution_to_jacobian(), RefineableQAxisymCrouzeixRaviartBoussinesqElement::fill_in_contribution_to_jacobian(), QAxisymCrouzeixRaviartElementWithExternalElement::fill_in_contribution_to_jacobian(), and oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::fill_in_contribution_to_jacobian().

◆ fill_in_contribution_to_jacobian_and_mass_matrix()

void oomph::AxisymmetricNavierStokesEquations::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< AxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >.

807  {
808  // Call the generic routine with the flag set to 2
809  fill_in_generic_residual_contribution_axi_nst(
810  residuals, jacobian, mass_matrix, 2);
811  }

References fill_in_generic_residual_contribution_axi_nst().

Referenced by oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::fill_in_contribution_to_jacobian_and_mass_matrix().

◆ fill_in_contribution_to_residuals()

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

Compute the element's residual Vector.

Reimplemented from oomph::GeneralisedElement.

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

781  {
782  // Call the generic residuals function with flag set to 0
783  // and using a dummy matrix argument
784  fill_in_generic_residual_contribution_axi_nst(
785  residuals,
786  GeneralisedElement::Dummy_matrix,
787  GeneralisedElement::Dummy_matrix,
788  0);
789  }

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

Referenced by oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::fill_in_contribution_to_residuals(), and oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::fill_in_contribution_to_residuals().

◆ fill_in_generic_dresidual_contribution_axi_nst()

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

2516  {
2517  // Die if the parameter is not the Reynolds number
2518  if (parameter_pt != this->re_pt())
2519  {
2520  std::ostringstream error_stream;
2521  error_stream
2522  << "Cannot compute analytic jacobian for parameter addressed by "
2523  << parameter_pt << "\n";
2524  error_stream << "Can only compute derivatives wrt Re (" << Re_pt << ")\n";
2525  throw OomphLibError(
2526  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2527  }
2528 
2529  // Which parameters are we differentiating with respect to
2530  bool diff_re = false;
2531  bool diff_re_st = false;
2532  bool diff_re_inv_fr = false;
2533  bool diff_re_inv_ro = false;
2534 
2535  // Set the boolean flags according to the parameter pointer
2536  if (parameter_pt == this->re_pt())
2537  {
2538  diff_re = true;
2539  }
2540  if (parameter_pt == this->re_st_pt())
2541  {
2542  diff_re_st = true;
2543  }
2544  if (parameter_pt == this->re_invfr_pt())
2545  {
2546  diff_re_inv_fr = true;
2547  }
2548  if (parameter_pt == this->re_invro_pt())
2549  {
2550  diff_re_inv_ro = true;
2551  }
2552 
2553 
2554  // Find out how many nodes there are
2555  unsigned n_node = nnode();
2556 
2557  // Find out how many pressure dofs there are
2558  unsigned n_pres = npres_axi_nst();
2559 
2560  // Get the nodal indices at which the velocity is stored
2561  unsigned u_nodal_index[3];
2562  for (unsigned i = 0; i < 3; ++i)
2563  {
2564  u_nodal_index[i] = u_index_axi_nst(i);
2565  }
2566 
2567  // Set up memory for the shape and test functions
2568  // Note that there are only two dimensions, r and z in this problem
2569  Shape psif(n_node), testf(n_node);
2570  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
2571 
2572  // Set up memory for pressure shape and test functions
2573  Shape psip(n_pres), testp(n_pres);
2574 
2575  // Number of integration points
2576  unsigned Nintpt = integral_pt()->nweight();
2577 
2578  // Set the Vector to hold local coordinates (two dimensions)
2579  Vector<double> s(2);
2580 
2581  // Get Physical Variables from Element
2582  // Reynolds number must be multiplied by the density ratio
2583  // double scaled_re = re()*density_ratio();
2584  // double scaled_re_st = re_st()*density_ratio();
2585  // double scaled_re_inv_fr = re_invfr()*density_ratio();
2586  // double scaled_re_inv_ro = re_invro()*density_ratio();
2587  double dens_ratio = this->density_ratio();
2588  // double visc_ratio = viscosity_ratio();
2589  Vector<double> G = g();
2590 
2591  // Integers used to store the local equation and unknown numbers
2592  int local_eqn = 0, local_unknown = 0;
2593 
2594  // Loop over the integration points
2595  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
2596  {
2597  // Assign values of s
2598  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
2599  // Get the integral weight
2600  double w = integral_pt()->weight(ipt);
2601 
2602  // Call the derivatives of the shape and test functions
2603  double J = dshape_and_dtest_eulerian_at_knot_axi_nst(
2604  ipt, psif, dpsifdx, testf, dtestfdx);
2605 
2606  // Call the pressure shape and test functions
2607  pshape_axi_nst(s, psip, testp);
2608 
2609  // Premultiply the weights and the Jacobian
2610  double W = w * J;
2611 
2612  // Allocate storage for the position and the derivative of the
2613  // mesh positions wrt time
2614  Vector<double> interpolated_x(2, 0.0);
2615  Vector<double> mesh_velocity(2, 0.0);
2616  // Allocate storage for the pressure, velocity components and their
2617  // time and spatial derivatives
2618  double interpolated_p = 0.0;
2619  Vector<double> interpolated_u(3, 0.0);
2620  Vector<double> dudt(3, 0.0);
2621  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
2622 
2623  // Calculate pressure at integration point
2624  for (unsigned l = 0; l < n_pres; l++)
2625  {
2626  interpolated_p += p_axi_nst(l) * psip[l];
2627  }
2628 
2629  // Calculate velocities and derivatives at integration point
2630 
2631  // Loop over nodes
2632  for (unsigned l = 0; l < n_node; l++)
2633  {
2634  // Cache the shape function
2635  const double psif_ = psif(l);
2636  // Loop over the two coordinate directions
2637  for (unsigned i = 0; i < 2; i++)
2638  {
2639  interpolated_x[i] += this->raw_nodal_position(l, i) * psif_;
2640  }
2641  // mesh_velocity[i] += dnodal_position_dt(l,i)*psif[l];
2642 
2643  // Loop over the three velocity directions
2644  for (unsigned i = 0; i < 3; i++)
2645  {
2646  // Get the u_value
2647  const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
2648  interpolated_u[i] += u_value * psif_;
2649  dudt[i] += du_dt_axi_nst(l, i) * psif_;
2650  // Loop over derivative directions
2651  for (unsigned j = 0; j < 2; j++)
2652  {
2653  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
2654  }
2655  }
2656  }
2657 
2658  // Get the mesh velocity if ALE is enabled
2659  if (!ALE_is_disabled)
2660  {
2661  // Loop over nodes
2662  for (unsigned l = 0; l < n_node; l++)
2663  {
2664  // Loop over the two coordinate directions
2665  for (unsigned i = 0; i < 2; i++)
2666  {
2667  mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif(l);
2668  }
2669  }
2670  }
2671 
2672 
2673  // Get the user-defined body force terms
2674  // Vector<double> body_force(3);
2675  // get_body_force(time(),ipt,interpolated_x,body_force);
2676 
2677  // Get the user-defined source function
2678  // double source = get_source_fct(time(),ipt,interpolated_x);
2679 
2680  // Get the user-defined viscosity function
2681  double visc_ratio;
2682  get_viscosity_ratio_axisym_nst(ipt, s, interpolated_x, visc_ratio);
2683 
2684  // r is the first position component
2685  double r = interpolated_x[0];
2686 
2687 
2688  // MOMENTUM EQUATIONS
2689  //------------------
2690 
2691  // Loop over the test functions
2692  for (unsigned l = 0; l < n_node; l++)
2693  {
2694  // FIRST (RADIAL) MOMENTUM EQUATION
2695  local_eqn = nodal_local_eqn(l, u_nodal_index[0]);
2696  // If it's not a boundary condition
2697  if (local_eqn >= 0)
2698  {
2699  // No user-defined body force terms
2700  // dres_dparam[local_eqn] +=
2701  // r*body_force[0]*testf[l]*W;
2702 
2703  // Add the gravitational body force term if the reynolds number
2704  // is equal to re_inv_fr
2705  if (diff_re_inv_fr)
2706  {
2707  dres_dparam[local_eqn] += r * dens_ratio * testf[l] * G[0] * W;
2708  }
2709 
2710  // No pressure gradient term
2711  // residuals[local_eqn] +=
2712  // interpolated_p*(testf[l] + r*dtestfdx(l,0))*W;
2713 
2714  // No in the stress tensor terms
2715  // The viscosity ratio needs to go in here to ensure
2716  // continuity of normal stress is satisfied even in flows
2717  // with zero pressure gradient!
2718  // residuals[local_eqn] -= visc_ratio*
2719  // r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W;
2720 
2721  // residuals[local_eqn] -= visc_ratio*r*
2722  // (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
2723  // dtestfdx(l,1)*W;
2724 
2725  // residuals[local_eqn] -=
2726  // visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W/r;
2727 
2728  // Add in the inertial terms
2729  // du/dt term
2730  if (diff_re_st)
2731  {
2732  dres_dparam[local_eqn] -= dens_ratio * r * dudt[0] * testf[l] * W;
2733  }
2734 
2735  // Convective terms
2736  if (diff_re)
2737  {
2738  dres_dparam[local_eqn] -=
2739  dens_ratio *
2740  (r * interpolated_u[0] * interpolated_dudx(0, 0) -
2741  interpolated_u[2] * interpolated_u[2] +
2742  r * interpolated_u[1] * interpolated_dudx(0, 1)) *
2743  testf[l] * W;
2744  }
2745 
2746  // Mesh velocity terms
2747  if (!ALE_is_disabled)
2748  {
2749  if (diff_re_st)
2750  {
2751  for (unsigned k = 0; k < 2; k++)
2752  {
2753  dres_dparam[local_eqn] += dens_ratio * r * mesh_velocity[k] *
2754  interpolated_dudx(0, k) * testf[l] *
2755  W;
2756  }
2757  }
2758  }
2759 
2760  // Add the Coriolis term
2761  if (diff_re_inv_ro)
2762  {
2763  dres_dparam[local_eqn] +=
2764  2.0 * r * dens_ratio * interpolated_u[2] * testf[l] * W;
2765  }
2766 
2767  // CALCULATE THE JACOBIAN
2768  if (flag)
2769  {
2770  // Loop over the velocity shape functions again
2771  for (unsigned l2 = 0; l2 < n_node; l2++)
2772  {
2773  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
2774  // Radial velocity component
2775  if (local_unknown >= 0)
2776  {
2777  if (flag == 2)
2778  {
2779  if (diff_re_st)
2780  {
2781  // Add the mass matrix
2782  dmass_matrix_dparam(local_eqn, local_unknown) +=
2783  dens_ratio * r * psif[l2] * testf[l] * W;
2784  }
2785  }
2786 
2787  // Add contribution to the Jacobian matrix
2788  // jacobian(local_eqn,local_unknown)
2789  // -= visc_ratio*r*(1.0+Gamma[0])
2790  //*dpsifdx(l2,0)*dtestfdx(l,0)*W;
2791 
2792  // jacobian(local_eqn,local_unknown)
2793  // -= visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
2794 
2795  // jacobian(local_eqn,local_unknown)
2796  // -= visc_ratio*(1.0 + Gamma[0])*psif[l2]*testf[l]*W/r;
2797 
2798  // Add in the inertial terms
2799  // du/dt term
2800  if (diff_re_st)
2801  {
2802  djac_dparam(local_eqn, local_unknown) -=
2803  dens_ratio * r *
2804  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
2805  testf[l] * W;
2806  }
2807 
2808  // Convective terms
2809  if (diff_re)
2810  {
2811  djac_dparam(local_eqn, local_unknown) -=
2812  dens_ratio *
2813  (r * psif[l2] * interpolated_dudx(0, 0) +
2814  r * interpolated_u[0] * dpsifdx(l2, 0) +
2815  r * interpolated_u[1] * dpsifdx(l2, 1)) *
2816  testf[l] * W;
2817  }
2818 
2819  // Mesh velocity terms
2820  if (!ALE_is_disabled)
2821  {
2822  for (unsigned k = 0; k < 2; k++)
2823  {
2824  if (diff_re_st)
2825  {
2826  djac_dparam(local_eqn, local_unknown) +=
2827  dens_ratio * r * mesh_velocity[k] * dpsifdx(l2, k) *
2828  testf[l] * W;
2829  }
2830  }
2831  }
2832  }
2833 
2834  // Axial velocity component
2835  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
2836  if (local_unknown >= 0)
2837  {
2838  // jacobian(local_eqn,local_unknown) -=
2839  // visc_ratio*r*Gamma[0]*dpsifdx(l2,0)*dtestfdx(l,1)*W;
2840 
2841  // Convective terms
2842  if (diff_re)
2843  {
2844  djac_dparam(local_eqn, local_unknown) -=
2845  dens_ratio * r * psif[l2] * interpolated_dudx(0, 1) *
2846  testf[l] * W;
2847  }
2848  }
2849 
2850  // Azimuthal velocity component
2851  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
2852  if (local_unknown >= 0)
2853  {
2854  // Convective terms
2855  if (diff_re)
2856  {
2857  djac_dparam(local_eqn, local_unknown) -=
2858  -dens_ratio * 2.0 * interpolated_u[2] * psif[l2] *
2859  testf[l] * W;
2860  }
2861  // Coriolis terms
2862  if (diff_re_inv_ro)
2863  {
2864  djac_dparam(local_eqn, local_unknown) +=
2865  2.0 * r * dens_ratio * psif[l2] * testf[l] * W;
2866  }
2867  }
2868  }
2869 
2870  /*Now loop over pressure shape functions*/
2871  /*This is the contribution from pressure gradient*/
2872  // for(unsigned l2=0;l2<n_pres;l2++)
2873  // {
2874  // local_unknown = p_local_eqn(l2);
2875  // /*If we are at a non-zero degree of freedom in the entry*/
2876  // if(local_unknown >= 0)
2877  // {
2878  // jacobian(local_eqn,local_unknown)
2879  // += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W;
2880  // }
2881  // }
2882  } /*End of Jacobian calculation*/
2883 
2884  } // End of if not boundary condition statement
2885 
2886  // SECOND (AXIAL) MOMENTUM EQUATION
2887  local_eqn = nodal_local_eqn(l, u_nodal_index[1]);
2888  // If it's not a boundary condition
2889  if (local_eqn >= 0)
2890  {
2891  // Add the user-defined body force terms
2892  // Remember to multiply by the density ratio!
2893  // residuals[local_eqn] +=
2894  // r*body_force[1]*testf[l]*W;
2895 
2896  // Add the gravitational body force term
2897  if (diff_re_inv_fr)
2898  {
2899  dres_dparam[local_eqn] += r * dens_ratio * testf[l] * G[1] * W;
2900  }
2901 
2902  // Add the pressure gradient term
2903  // residuals[local_eqn] += r*interpolated_p*dtestfdx(l,1)*W;
2904 
2905  // Add in the stress tensor terms
2906  // The viscosity ratio needs to go in here to ensure
2907  // continuity of normal stress is satisfied even in flows
2908  // with zero pressure gradient!
2909  // residuals[local_eqn] -= visc_ratio*
2910  // r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
2911  // *dtestfdx(l,0)*W;
2912 
2913  // residuals[local_eqn] -= visc_ratio*r*
2914  // (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W;
2915 
2916  // Add in the inertial terms
2917  // du/dt term
2918  if (diff_re_st)
2919  {
2920  dres_dparam[local_eqn] -= dens_ratio * r * dudt[1] * testf[l] * W;
2921  }
2922 
2923  // Convective terms
2924  if (diff_re)
2925  {
2926  dres_dparam[local_eqn] -=
2927  dens_ratio *
2928  (r * interpolated_u[0] * interpolated_dudx(1, 0) +
2929  r * interpolated_u[1] * interpolated_dudx(1, 1)) *
2930  testf[l] * W;
2931  }
2932 
2933  // Mesh velocity terms
2934  if (!ALE_is_disabled)
2935  {
2936  if (diff_re_st)
2937  {
2938  for (unsigned k = 0; k < 2; k++)
2939  {
2940  dres_dparam[local_eqn] += dens_ratio * r * mesh_velocity[k] *
2941  interpolated_dudx(1, k) * testf[l] *
2942  W;
2943  }
2944  }
2945  }
2946 
2947  // CALCULATE THE JACOBIAN
2948  if (flag)
2949  {
2950  // Loop over the velocity shape functions again
2951  for (unsigned l2 = 0; l2 < n_node; l2++)
2952  {
2953  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
2954  // Radial velocity component
2955  if (local_unknown >= 0)
2956  {
2957  // Add in the stress tensor terms
2958  // The viscosity ratio needs to go in here to ensure
2959  // continuity of normal stress is satisfied even in flows
2960  // with zero pressure gradient!
2961  // jacobian(local_eqn,local_unknown) -=
2962  // visc_ratio*r*Gamma[1]*dpsifdx(l2,1)*dtestfdx(l,0)*W;
2963 
2964  // Convective terms
2965  if (diff_re)
2966  {
2967  djac_dparam(local_eqn, local_unknown) -=
2968  dens_ratio * r * psif[l2] * interpolated_dudx(1, 0) *
2969  testf[l] * W;
2970  }
2971  }
2972 
2973  // Axial velocity component
2974  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
2975  if (local_unknown >= 0)
2976  {
2977  if (flag == 2)
2978  {
2979  if (diff_re_st)
2980  {
2981  // Add the mass matrix
2982  dmass_matrix_dparam(local_eqn, local_unknown) +=
2983  dens_ratio * r * psif[l2] * testf[l] * W;
2984  }
2985  }
2986 
2987 
2988  // Add in the stress tensor terms
2989  // The viscosity ratio needs to go in here to ensure
2990  // continuity of normal stress is satisfied even in flows
2991  // with zero pressure gradient!
2992  // jacobian(local_eqn,local_unknown) -=
2993  // visc_ratio*r*dpsifdx(l2,0)*dtestfdx(l,0)*W;
2994 
2995  // jacobian(local_eqn,local_unknown) -=
2996  // visc_ratio*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
2997  // dtestfdx(l,1)*W;
2998 
2999  // Add in the inertial terms
3000  // du/dt term
3001  if (diff_re_st)
3002  {
3003  djac_dparam(local_eqn, local_unknown) -=
3004  dens_ratio * r *
3005  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
3006  testf[l] * W;
3007  }
3008 
3009  // Convective terms
3010  if (diff_re)
3011  {
3012  djac_dparam(local_eqn, local_unknown) -=
3013  dens_ratio *
3014  (r * interpolated_u[0] * dpsifdx(l2, 0) +
3015  r * psif[l2] * interpolated_dudx(1, 1) +
3016  r * interpolated_u[1] * dpsifdx(l2, 1)) *
3017  testf[l] * W;
3018  }
3019 
3020  // Mesh velocity terms
3021  if (!ALE_is_disabled)
3022  {
3023  for (unsigned k = 0; k < 2; k++)
3024  {
3025  if (diff_re_st)
3026  {
3027  djac_dparam(local_eqn, local_unknown) +=
3028  dens_ratio * r * mesh_velocity[k] * dpsifdx(l2, k) *
3029  testf[l] * W;
3030  }
3031  }
3032  }
3033  }
3034 
3035  // There are no azimithal terms in the axial momentum equation
3036  } // End of loop over velocity shape functions
3037 
3038  } /*End of Jacobian calculation*/
3039 
3040  } // End of AXIAL MOMENTUM EQUATION
3041 
3042  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
3043  local_eqn = nodal_local_eqn(l, u_nodal_index[2]);
3044  if (local_eqn >= 0)
3045  {
3046  // Add the user-defined body force terms
3047  // Remember to multiply by the density ratio!
3048  // residuals[local_eqn] +=
3049  // r*body_force[2]*testf[l]*W;
3050 
3051  // Add the gravitational body force term
3052  if (diff_re_inv_fr)
3053  {
3054  dres_dparam[local_eqn] += r * dens_ratio * testf[l] * G[2] * W;
3055  }
3056 
3057  // There is NO pressure gradient term
3058 
3059  // Add in the stress tensor terms
3060  // The viscosity ratio needs to go in here to ensure
3061  // continuity of normal stress is satisfied even in flows
3062  // with zero pressure gradient!
3063  // residuals[local_eqn] -= visc_ratio*
3064  // (r*interpolated_dudx(2,0) -
3065  // Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W;
3066 
3067  // residuals[local_eqn] -= visc_ratio*r*
3068  // interpolated_dudx(2,1)*dtestfdx(l,1)*W;
3069 
3070  // residuals[local_eqn] -= visc_ratio*
3071  // ((interpolated_u[2]/r) -
3072  // Gamma[0]*interpolated_dudx(2,0))*testf[l]*W;
3073 
3074 
3075  // Add in the inertial terms
3076  // du/dt term
3077  if (diff_re_st)
3078  {
3079  dres_dparam[local_eqn] -= dens_ratio * r * dudt[2] * testf[l] * W;
3080  }
3081 
3082  // Convective terms
3083  if (diff_re)
3084  {
3085  dres_dparam[local_eqn] -=
3086  dens_ratio *
3087  (r * interpolated_u[0] * interpolated_dudx(2, 0) +
3088  interpolated_u[0] * interpolated_u[2] +
3089  r * interpolated_u[1] * interpolated_dudx(2, 1)) *
3090  testf[l] * W;
3091  }
3092 
3093  // Mesh velocity terms
3094  if (!ALE_is_disabled)
3095  {
3096  if (diff_re_st)
3097  {
3098  for (unsigned k = 0; k < 2; k++)
3099  {
3100  dres_dparam[local_eqn] += dens_ratio * r * mesh_velocity[k] *
3101  interpolated_dudx(2, k) * testf[l] *
3102  W;
3103  }
3104  }
3105  }
3106 
3107  // Add the Coriolis term
3108  if (diff_re_inv_ro)
3109  {
3110  dres_dparam[local_eqn] -=
3111  2.0 * r * dens_ratio * interpolated_u[0] * testf[l] * W;
3112  }
3113 
3114  // CALCULATE THE JACOBIAN
3115  if (flag)
3116  {
3117  // Loop over the velocity shape functions again
3118  for (unsigned l2 = 0; l2 < n_node; l2++)
3119  {
3120  // Radial velocity component
3121  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
3122  if (local_unknown >= 0)
3123  {
3124  // Convective terms
3125  if (diff_re)
3126  {
3127  djac_dparam(local_eqn, local_unknown) -=
3128  dens_ratio *
3129  (r * psif[l2] * interpolated_dudx(2, 0) +
3130  psif[l2] * interpolated_u[2]) *
3131  testf[l] * W;
3132  }
3133 
3134  // Coriolis term
3135  if (diff_re_inv_ro)
3136  {
3137  djac_dparam(local_eqn, local_unknown) -=
3138  2.0 * r * dens_ratio * psif[l2] * testf[l] * W;
3139  }
3140  }
3141 
3142  // Axial velocity component
3143  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
3144  if (local_unknown >= 0)
3145  {
3146  // Convective terms
3147  if (diff_re)
3148  {
3149  djac_dparam(local_eqn, local_unknown) -=
3150  dens_ratio * r * psif[l2] * interpolated_dudx(2, 1) *
3151  testf[l] * W;
3152  }
3153  }
3154 
3155  // Azimuthal velocity component
3156  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
3157  if (local_unknown >= 0)
3158  {
3159  if (flag == 2)
3160  {
3161  // Add the mass matrix
3162  if (diff_re_st)
3163  {
3164  dmass_matrix_dparam(local_eqn, local_unknown) +=
3165  dens_ratio * r * psif[l2] * testf[l] * W;
3166  }
3167  }
3168 
3169  // Add in the stress tensor terms
3170  // The viscosity ratio needs to go in here to ensure
3171  // continuity of normal stress is satisfied even in flows
3172  // with zero pressure gradient!
3173  // jacobian(local_eqn,local_unknown) -=
3174  // visc_ratio*(r*dpsifdx(l2,0) -
3175  // Gamma[0]*psif[l2])*dtestfdx(l,0)*W;
3176 
3177  // jacobian(local_eqn,local_unknown) -=
3178  // visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)*W;
3179 
3180  // jacobian(local_eqn,local_unknown) -=
3181  // visc_ratio*((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
3182  // *testf[l]*W;
3183 
3184  // Add in the inertial terms
3185  // du/dt term
3186  if (diff_re_st)
3187  {
3188  djac_dparam(local_eqn, local_unknown) -=
3189  dens_ratio * r *
3190  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
3191  testf[l] * W;
3192  }
3193 
3194  // Convective terms
3195  if (diff_re)
3196  {
3197  djac_dparam(local_eqn, local_unknown) -=
3198  dens_ratio *
3199  (r * interpolated_u[0] * dpsifdx(l2, 0) +
3200  interpolated_u[0] * psif[l2] +
3201  r * interpolated_u[1] * dpsifdx(l2, 1)) *
3202  testf[l] * W;
3203  }
3204 
3205  // Mesh velocity terms
3206  if (!ALE_is_disabled)
3207  {
3208  if (diff_re_st)
3209  {
3210  for (unsigned k = 0; k < 2; k++)
3211  {
3212  djac_dparam(local_eqn, local_unknown) +=
3213  dens_ratio * r * mesh_velocity[k] * dpsifdx(l2, k) *
3214  testf[l] * W;
3215  }
3216  }
3217  }
3218  }
3219  }
3220 
3221  // There are no pressure terms
3222  } // End of Jacobian
3223 
3224  } // End of AZIMUTHAL EQUATION
3225 
3226  } // End of loop over shape functions
3227 
3228 
3229  // CONTINUITY EQUATION NO PARAMETERS
3230  //-------------------
3231  }
3232  }
oomph::AxisymmetricNavierStokesEquations::p_axi_nst
virtual double p_axi_nst(const unsigned &n_p) const =0
oomph::AxisymmetricNavierStokesEquations::re_invro_pt
double *& re_invro_pt()
Pointer to global inverse Froude number.
Definition: axisym_navier_stokes_elements.h:440
oomph::AxisymmetricNavierStokesEquations::pshape_axi_nst
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
oomph::AxisymmetricNavierStokesEquations::re_pt
double *& re_pt()
Pointer to Reynolds number.
Definition: axisym_navier_stokes_elements.h:410
oomph::AxisymmetricNavierStokesEquations::get_viscosity_ratio_axisym_nst
virtual void get_viscosity_ratio_axisym_nst(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, double &visc_ratio)
Definition: axisym_navier_stokes_elements.h:339
oomph::AxisymmetricNavierStokesEquations::re_st_pt
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
Definition: axisym_navier_stokes_elements.h:416
oomph::AxisymmetricNavierStokesEquations::du_dt_axi_nst
double du_dt_axi_nst(const unsigned &n, const unsigned &i) const
Definition: axisym_navier_stokes_elements.h:531
oomph::AxisymmetricNavierStokesEquations::re_invfr_pt
double *& re_invfr_pt()
Pointer to global inverse Froude number.
Definition: axisym_navier_stokes_elements.h:428
oomph::AxisymmetricNavierStokesEquations::npres_axi_nst
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
oomph::FiniteElement::raw_dnodal_position_dt
double raw_dnodal_position_dt(const unsigned &n, const unsigned &i) const
Definition: elements.h:2256

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::AxisymmetricNavierStokesEquations::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::RefineableAxisymmetricNavierStokesEquations.

919  {
920  // Find out how many nodes there are
921  unsigned n_node = nnode();
922 
923  // Get continuous time from timestepper of first node
924  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
925 
926  // Find out how many pressure dofs there are
927  unsigned n_pres = npres_axi_nst();
928 
929  // Get the nodal indices at which the velocity is stored
930  unsigned u_nodal_index[3];
931  for (unsigned i = 0; i < 3; ++i)
932  {
933  u_nodal_index[i] = u_index_axi_nst(i);
934  }
935 
936  // Set up memory for the shape and test functions
937  // Note that there are only two dimensions, r and z in this problem
938  Shape psif(n_node), testf(n_node);
939  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
940 
941  // Set up memory for pressure shape and test functions
942  Shape psip(n_pres), testp(n_pres);
943 
944  // Number of integration points
945  unsigned Nintpt = integral_pt()->nweight();
946 
947  // Set the Vector to hold local coordinates (two dimensions)
948  Vector<double> s(2);
949 
950  // Get Physical Variables from Element
951  // Reynolds number must be multiplied by the density ratio
952  double scaled_re = re() * density_ratio();
953  double scaled_re_st = re_st() * density_ratio();
954  double scaled_re_inv_fr = re_invfr() * density_ratio();
955  double scaled_re_inv_ro = re_invro() * density_ratio();
956  // double visc_ratio = viscosity_ratio();
957  Vector<double> G = g();
958 
959  // Integers used to store the local equation and unknown numbers
960  int local_eqn = 0, local_unknown = 0;
961 
962  // Loop over the integration points
963  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
964  {
965  // Assign values of s
966  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
967  // Get the integral weight
968  double w = integral_pt()->weight(ipt);
969 
970  // Call the derivatives of the shape and test functions
971  double J = dshape_and_dtest_eulerian_at_knot_axi_nst(
972  ipt, psif, dpsifdx, testf, dtestfdx);
973 
974  // Call the pressure shape and test functions
975  pshape_axi_nst(s, psip, testp);
976 
977  // Premultiply the weights and the Jacobian
978  double W = w * J;
979 
980  // Allocate storage for the position and the derivative of the
981  // mesh positions wrt time
982  Vector<double> interpolated_x(2, 0.0);
983  Vector<double> mesh_velocity(2, 0.0);
984  // Allocate storage for the pressure, velocity components and their
985  // time and spatial derivatives
986  double interpolated_p = 0.0;
987  Vector<double> interpolated_u(3, 0.0);
988  Vector<double> dudt(3, 0.0);
989  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
990 
991  // Calculate pressure at integration point
992  for (unsigned l = 0; l < n_pres; l++)
993  {
994  interpolated_p += p_axi_nst(l) * psip[l];
995  }
996 
997  // Calculate velocities and derivatives at integration point
998 
999  // Loop over nodes
1000  for (unsigned l = 0; l < n_node; l++)
1001  {
1002  // Cache the shape function
1003  const double psif_ = psif(l);
1004  // Loop over the two coordinate directions
1005  for (unsigned i = 0; i < 2; i++)
1006  {
1007  interpolated_x[i] += this->raw_nodal_position(l, i) * psif_;
1008  }
1009  // mesh_velocity[i] += dnodal_position_dt(l,i)*psif[l];
1010 
1011  // Loop over the three velocity directions
1012  for (unsigned i = 0; i < 3; i++)
1013  {
1014  // Get the u_value
1015  const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
1016  interpolated_u[i] += u_value * psif_;
1017  dudt[i] += du_dt_axi_nst(l, i) * psif_;
1018  // Loop over derivative directions
1019  for (unsigned j = 0; j < 2; j++)
1020  {
1021  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1022  }
1023  }
1024  }
1025 
1026  // Get the mesh velocity if ALE is enabled
1027  if (!ALE_is_disabled)
1028  {
1029  // Loop over nodes
1030  for (unsigned l = 0; l < n_node; l++)
1031  {
1032  // Loop over the two coordinate directions
1033  for (unsigned i = 0; i < 2; i++)
1034  {
1035  mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif(l);
1036  }
1037  }
1038  }
1039 
1040 
1041  // Get the user-defined body force terms
1042  Vector<double> body_force(3);
1043  get_body_force_axi_nst(time, ipt, s, interpolated_x, body_force);
1044 
1045  // Get the user-defined source function
1046  double source = get_source_fct(time, ipt, interpolated_x);
1047 
1048  // Get the user-defined viscosity function
1049  double visc_ratio;
1050  get_viscosity_ratio_axisym_nst(ipt, s, interpolated_x, visc_ratio);
1051 
1052  // r is the first position component
1053  double r = interpolated_x[0];
1054 
1055  // MOMENTUM EQUATIONS
1056  //------------------
1057 
1058  // Loop over the test functions
1059  for (unsigned l = 0; l < n_node; l++)
1060  {
1061  // FIRST (RADIAL) MOMENTUM EQUATION
1062  local_eqn = nodal_local_eqn(l, u_nodal_index[0]);
1063  // If it's not a boundary condition
1064  if (local_eqn >= 0)
1065  {
1066  // Add the user-defined body force terms
1067  residuals[local_eqn] += r * body_force[0] * testf[l] * W;
1068 
1069  // Add the gravitational body force term
1070  residuals[local_eqn] += r * scaled_re_inv_fr * testf[l] * G[0] * W;
1071 
1072  // Add the pressure gradient term
1073  residuals[local_eqn] +=
1074  interpolated_p * (testf[l] + r * dtestfdx(l, 0)) * W;
1075 
1076  // Add in the stress tensor terms
1077  // The viscosity ratio needs to go in here to ensure
1078  // continuity of normal stress is satisfied even in flows
1079  // with zero pressure gradient!
1080  residuals[local_eqn] -= visc_ratio * r * (1.0 + Gamma[0]) *
1081  interpolated_dudx(0, 0) * dtestfdx(l, 0) * W;
1082 
1083  residuals[local_eqn] -=
1084  visc_ratio * r *
1085  (interpolated_dudx(0, 1) + Gamma[0] * interpolated_dudx(1, 0)) *
1086  dtestfdx(l, 1) * W;
1087 
1088  residuals[local_eqn] -= visc_ratio * (1.0 + Gamma[0]) *
1089  interpolated_u[0] * testf[l] * W / r;
1090 
1091  // Add in the inertial terms
1092  // du/dt term
1093  residuals[local_eqn] -= scaled_re_st * r * dudt[0] * testf[l] * W;
1094 
1095  // Convective terms
1096  residuals[local_eqn] -=
1097  scaled_re *
1098  (r * interpolated_u[0] * interpolated_dudx(0, 0) -
1099  interpolated_u[2] * interpolated_u[2] +
1100  r * interpolated_u[1] * interpolated_dudx(0, 1)) *
1101  testf[l] * W;
1102 
1103  // Mesh velocity terms
1104  if (!ALE_is_disabled)
1105  {
1106  for (unsigned k = 0; k < 2; k++)
1107  {
1108  residuals[local_eqn] += scaled_re_st * r * mesh_velocity[k] *
1109  interpolated_dudx(0, k) * testf[l] * W;
1110  }
1111  }
1112 
1113  // Add the Coriolis term
1114  residuals[local_eqn] +=
1115  2.0 * r * scaled_re_inv_ro * interpolated_u[2] * testf[l] * W;
1116 
1117  // CALCULATE THE JACOBIAN
1118  if (flag)
1119  {
1120  // Loop over the velocity shape functions again
1121  for (unsigned l2 = 0; l2 < n_node; l2++)
1122  {
1123  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
1124  // Radial velocity component
1125  if (local_unknown >= 0)
1126  {
1127  if (flag == 2)
1128  {
1129  // Add the mass matrix
1130  mass_matrix(local_eqn, local_unknown) +=
1131  scaled_re_st * r * psif[l2] * testf[l] * W;
1132  }
1133 
1134  // Add contribution to the Jacobian matrix
1135  jacobian(local_eqn, local_unknown) -=
1136  visc_ratio * r * (1.0 + Gamma[0]) * dpsifdx(l2, 0) *
1137  dtestfdx(l, 0) * W;
1138 
1139  jacobian(local_eqn, local_unknown) -=
1140  visc_ratio * r * dpsifdx(l2, 1) * dtestfdx(l, 1) * W;
1141 
1142  jacobian(local_eqn, local_unknown) -=
1143  visc_ratio * (1.0 + Gamma[0]) * psif[l2] * testf[l] * W / r;
1144 
1145  // Add in the inertial terms
1146  // du/dt term
1147  jacobian(local_eqn, local_unknown) -=
1148  scaled_re_st * r *
1149  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
1150  testf[l] * W;
1151 
1152  // Convective terms
1153  jacobian(local_eqn, local_unknown) -=
1154  scaled_re *
1155  (r * psif[l2] * interpolated_dudx(0, 0) +
1156  r * interpolated_u[0] * dpsifdx(l2, 0) +
1157  r * interpolated_u[1] * dpsifdx(l2, 1)) *
1158  testf[l] * W;
1159 
1160  // Mesh velocity terms
1161  if (!ALE_is_disabled)
1162  {
1163  for (unsigned k = 0; k < 2; k++)
1164  {
1165  jacobian(local_eqn, local_unknown) +=
1166  scaled_re_st * r * mesh_velocity[k] * dpsifdx(l2, k) *
1167  testf[l] * W;
1168  }
1169  }
1170  }
1171 
1172 
1173  // Axial velocity component
1174  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
1175  if (local_unknown >= 0)
1176  {
1177  jacobian(local_eqn, local_unknown) -=
1178  visc_ratio * r * Gamma[0] * dpsifdx(l2, 0) * dtestfdx(l, 1) *
1179  W;
1180 
1181  // Convective terms
1182  jacobian(local_eqn, local_unknown) -= scaled_re * r * psif[l2] *
1183  interpolated_dudx(0, 1) *
1184  testf[l] * W;
1185  }
1186 
1187  // Azimuthal velocity component
1188  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
1189  if (local_unknown >= 0)
1190  {
1191  // Convective terms
1192  jacobian(local_eqn, local_unknown) -= -scaled_re * 2.0 *
1193  interpolated_u[2] *
1194  psif[l2] * testf[l] * W;
1195 
1196  // Coriolis terms
1197  jacobian(local_eqn, local_unknown) +=
1198  2.0 * r * scaled_re_inv_ro * psif[l2] * testf[l] * W;
1199  }
1200  }
1201 
1202  /*Now loop over pressure shape functions*/
1203  /*This is the contribution from pressure gradient*/
1204  for (unsigned l2 = 0; l2 < n_pres; l2++)
1205  {
1206  local_unknown = p_local_eqn(l2);
1207  /*If we are at a non-zero degree of freedom in the entry*/
1208  if (local_unknown >= 0)
1209  {
1210  jacobian(local_eqn, local_unknown) +=
1211  psip[l2] * (testf[l] + r * dtestfdx(l, 0)) * W;
1212  }
1213  }
1214  } /*End of Jacobian calculation*/
1215 
1216  } // End of if not boundary condition statement
1217 
1218  // SECOND (AXIAL) MOMENTUM EQUATION
1219  local_eqn = nodal_local_eqn(l, u_nodal_index[1]);
1220  // If it's not a boundary condition
1221  if (local_eqn >= 0)
1222  {
1223  // Add the user-defined body force terms
1224  // Remember to multiply by the density ratio!
1225  residuals[local_eqn] += r * body_force[1] * testf[l] * W;
1226 
1227  // Add the gravitational body force term
1228  residuals[local_eqn] += r * scaled_re_inv_fr * testf[l] * G[1] * W;
1229 
1230  // Add the pressure gradient term
1231  residuals[local_eqn] += r * interpolated_p * dtestfdx(l, 1) * W;
1232 
1233  // Add in the stress tensor terms
1234  // The viscosity ratio needs to go in here to ensure
1235  // continuity of normal stress is satisfied even in flows
1236  // with zero pressure gradient!
1237  residuals[local_eqn] -=
1238  visc_ratio * r *
1239  (interpolated_dudx(1, 0) + Gamma[1] * interpolated_dudx(0, 1)) *
1240  dtestfdx(l, 0) * W;
1241 
1242  residuals[local_eqn] -= visc_ratio * r * (1.0 + Gamma[1]) *
1243  interpolated_dudx(1, 1) * dtestfdx(l, 1) * W;
1244 
1245  // Add in the inertial terms
1246  // du/dt term
1247  residuals[local_eqn] -= scaled_re_st * r * dudt[1] * testf[l] * W;
1248 
1249  // Convective terms
1250  residuals[local_eqn] -=
1251  scaled_re *
1252  (r * interpolated_u[0] * interpolated_dudx(1, 0) +
1253  r * interpolated_u[1] * interpolated_dudx(1, 1)) *
1254  testf[l] * W;
1255 
1256  // Mesh velocity terms
1257  if (!ALE_is_disabled)
1258  {
1259  for (unsigned k = 0; k < 2; k++)
1260  {
1261  residuals[local_eqn] += scaled_re_st * r * mesh_velocity[k] *
1262  interpolated_dudx(1, k) * testf[l] * W;
1263  }
1264  }
1265 
1266  // CALCULATE THE JACOBIAN
1267  if (flag)
1268  {
1269  // Loop over the velocity shape functions again
1270  for (unsigned l2 = 0; l2 < n_node; l2++)
1271  {
1272  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
1273  // Radial velocity component
1274  if (local_unknown >= 0)
1275  {
1276  // Add in the stress tensor terms
1277  // The viscosity ratio needs to go in here to ensure
1278  // continuity of normal stress is satisfied even in flows
1279  // with zero pressure gradient!
1280  jacobian(local_eqn, local_unknown) -=
1281  visc_ratio * r * Gamma[1] * dpsifdx(l2, 1) * dtestfdx(l, 0) *
1282  W;
1283 
1284  // Convective terms
1285  jacobian(local_eqn, local_unknown) -= scaled_re * r * psif[l2] *
1286  interpolated_dudx(1, 0) *
1287  testf[l] * W;
1288  }
1289 
1290  // Axial velocity component
1291  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
1292  if (local_unknown >= 0)
1293  {
1294  if (flag == 2)
1295  {
1296  // Add the mass matrix
1297  mass_matrix(local_eqn, local_unknown) +=
1298  scaled_re_st * r * psif[l2] * testf[l] * W;
1299  }
1300 
1301 
1302  // Add in the stress tensor terms
1303  // The viscosity ratio needs to go in here to ensure
1304  // continuity of normal stress is satisfied even in flows
1305  // with zero pressure gradient!
1306  jacobian(local_eqn, local_unknown) -=
1307  visc_ratio * r * dpsifdx(l2, 0) * dtestfdx(l, 0) * W;
1308 
1309  jacobian(local_eqn, local_unknown) -=
1310  visc_ratio * r * (1.0 + Gamma[1]) * dpsifdx(l2, 1) *
1311  dtestfdx(l, 1) * W;
1312 
1313  // Add in the inertial terms
1314  // du/dt term
1315  jacobian(local_eqn, local_unknown) -=
1316  scaled_re_st * r *
1317  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
1318  testf[l] * W;
1319 
1320  // Convective terms
1321  jacobian(local_eqn, local_unknown) -=
1322  scaled_re *
1323  (r * interpolated_u[0] * dpsifdx(l2, 0) +
1324  r * psif[l2] * interpolated_dudx(1, 1) +
1325  r * interpolated_u[1] * dpsifdx(l2, 1)) *
1326  testf[l] * W;
1327 
1328 
1329  // Mesh velocity terms
1330  if (!ALE_is_disabled)
1331  {
1332  for (unsigned k = 0; k < 2; k++)
1333  {
1334  jacobian(local_eqn, local_unknown) +=
1335  scaled_re_st * r * mesh_velocity[k] * dpsifdx(l2, k) *
1336  testf[l] * W;
1337  }
1338  }
1339  }
1340 
1341  // There are no azimithal terms in the axial momentum equation
1342  } // End of loop over velocity shape functions
1343 
1344  /*Now loop over pressure shape functions*/
1345  /*This is the contribution from pressure gradient*/
1346  for (unsigned l2 = 0; l2 < n_pres; l2++)
1347  {
1348  local_unknown = p_local_eqn(l2);
1349  /*If we are at a non-zero degree of freedom in the entry*/
1350  if (local_unknown >= 0)
1351  {
1352  jacobian(local_eqn, local_unknown) +=
1353  r * psip[l2] * dtestfdx(l, 1) * W;
1354  }
1355  }
1356  } /*End of Jacobian calculation*/
1357 
1358  } // End of AXIAL MOMENTUM EQUATION
1359 
1360  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
1361  local_eqn = nodal_local_eqn(l, u_nodal_index[2]);
1362  if (local_eqn >= 0)
1363  {
1364  // Add the user-defined body force terms
1365  // Remember to multiply by the density ratio!
1366  residuals[local_eqn] += r * body_force[2] * testf[l] * W;
1367 
1368  // Add the gravitational body force term
1369  residuals[local_eqn] += r * scaled_re_inv_fr * testf[l] * G[2] * W;
1370 
1371  // There is NO pressure gradient term
1372 
1373  // Add in the stress tensor terms
1374  // The viscosity ratio needs to go in here to ensure
1375  // continuity of normal stress is satisfied even in flows
1376  // with zero pressure gradient!
1377  residuals[local_eqn] -=
1378  visc_ratio *
1379  (r * interpolated_dudx(2, 0) - Gamma[0] * interpolated_u[2]) *
1380  dtestfdx(l, 0) * W;
1381 
1382  residuals[local_eqn] -=
1383  visc_ratio * r * interpolated_dudx(2, 1) * dtestfdx(l, 1) * W;
1384 
1385  residuals[local_eqn] -=
1386  visc_ratio *
1387  ((interpolated_u[2] / r) - Gamma[0] * interpolated_dudx(2, 0)) *
1388  testf[l] * W;
1389 
1390 
1391  // Add in the inertial terms
1392  // du/dt term
1393  residuals[local_eqn] -= scaled_re_st * r * dudt[2] * testf[l] * W;
1394 
1395  // Convective terms
1396  residuals[local_eqn] -=
1397  scaled_re *
1398  (r * interpolated_u[0] * interpolated_dudx(2, 0) +
1399  interpolated_u[0] * interpolated_u[2] +
1400  r * interpolated_u[1] * interpolated_dudx(2, 1)) *
1401  testf[l] * W;
1402 
1403  // Mesh velocity terms
1404  if (!ALE_is_disabled)
1405  {
1406  for (unsigned k = 0; k < 2; k++)
1407  {
1408  residuals[local_eqn] += scaled_re_st * r * mesh_velocity[k] *
1409  interpolated_dudx(2, k) * testf[l] * W;
1410  }
1411  }
1412 
1413  // Add the Coriolis term
1414  residuals[local_eqn] -=
1415  2.0 * r * scaled_re_inv_ro * interpolated_u[0] * testf[l] * W;
1416 
1417  // CALCULATE THE JACOBIAN
1418  if (flag)
1419  {
1420  // Loop over the velocity shape functions again
1421  for (unsigned l2 = 0; l2 < n_node; l2++)
1422  {
1423  // Radial velocity component
1424  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
1425  if (local_unknown >= 0)
1426  {
1427  // Convective terms
1428  jacobian(local_eqn, local_unknown) -=
1429  scaled_re *
1430  (r * psif[l2] * interpolated_dudx(2, 0) +
1431  psif[l2] * interpolated_u[2]) *
1432  testf[l] * W;
1433 
1434  // Coriolis term
1435  jacobian(local_eqn, local_unknown) -=
1436  2.0 * r * scaled_re_inv_ro * psif[l2] * testf[l] * W;
1437  }
1438 
1439  // Axial velocity component
1440  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
1441  if (local_unknown >= 0)
1442  {
1443  // Convective terms
1444  jacobian(local_eqn, local_unknown) -= scaled_re * r * psif[l2] *
1445  interpolated_dudx(2, 1) *
1446  testf[l] * W;
1447  }
1448 
1449  // Azimuthal velocity component
1450  local_unknown = nodal_local_eqn(l2, u_nodal_index[2]);
1451  if (local_unknown >= 0)
1452  {
1453  if (flag == 2)
1454  {
1455  // Add the mass matrix
1456  mass_matrix(local_eqn, local_unknown) +=
1457  scaled_re_st * r * psif[l2] * testf[l] * W;
1458  }
1459 
1460  // Add in the stress tensor terms
1461  // The viscosity ratio needs to go in here to ensure
1462  // continuity of normal stress is satisfied even in flows
1463  // with zero pressure gradient!
1464  jacobian(local_eqn, local_unknown) -=
1465  visc_ratio * (r * dpsifdx(l2, 0) - Gamma[0] * psif[l2]) *
1466  dtestfdx(l, 0) * W;
1467 
1468  jacobian(local_eqn, local_unknown) -=
1469  visc_ratio * r * dpsifdx(l2, 1) * dtestfdx(l, 1) * W;
1470 
1471  jacobian(local_eqn, local_unknown) -=
1472  visc_ratio * ((psif[l2] / r) - Gamma[0] * dpsifdx(l2, 0)) *
1473  testf[l] * W;
1474 
1475  // Add in the inertial terms
1476  // du/dt term
1477  jacobian(local_eqn, local_unknown) -=
1478  scaled_re_st * r *
1479  node_pt(l2)->time_stepper_pt()->weight(1, 0) * psif[l2] *
1480  testf[l] * W;
1481 
1482  // Convective terms
1483  jacobian(local_eqn, local_unknown) -=
1484  scaled_re *
1485  (r * interpolated_u[0] * dpsifdx(l2, 0) +
1486  interpolated_u[0] * psif[l2] +
1487  r * interpolated_u[1] * dpsifdx(l2, 1)) *
1488  testf[l] * W;
1489 
1490  // Mesh velocity terms
1491  if (!ALE_is_disabled)
1492  {
1493  for (unsigned k = 0; k < 2; k++)
1494  {
1495  jacobian(local_eqn, local_unknown) +=
1496  scaled_re_st * r * mesh_velocity[k] * dpsifdx(l2, k) *
1497  testf[l] * W;
1498  }
1499  }
1500  }
1501  }
1502 
1503  // There are no pressure terms
1504  } // End of Jacobian
1505 
1506  } // End of AZIMUTHAL EQUATION
1507 
1508  } // End of loop over shape functions
1509 
1510 
1511  // CONTINUITY EQUATION
1512  //-------------------
1513 
1514  // Loop over the shape functions
1515  for (unsigned l = 0; l < n_pres; l++)
1516  {
1517  local_eqn = p_local_eqn(l);
1518  // If not a boundary conditions
1519  if (local_eqn >= 0)
1520  {
1521  // Source term
1522  residuals[local_eqn] -= r * source * testp[l] * W;
1523 
1524  // Gradient terms
1525  residuals[local_eqn] +=
1526  (interpolated_u[0] + r * interpolated_dudx(0, 0) +
1527  r * interpolated_dudx(1, 1)) *
1528  testp[l] * W;
1529 
1530  /*CALCULATE THE JACOBIAN*/
1531  if (flag)
1532  {
1533  /*Loop over the velocity shape functions*/
1534  for (unsigned l2 = 0; l2 < n_node; l2++)
1535  {
1536  // Radial velocity component
1537  local_unknown = nodal_local_eqn(l2, u_nodal_index[0]);
1538  if (local_unknown >= 0)
1539  {
1540  jacobian(local_eqn, local_unknown) +=
1541  (psif[l2] + r * dpsifdx(l2, 0)) * testp[l] * W;
1542  }
1543 
1544  // Axial velocity component
1545  local_unknown = nodal_local_eqn(l2, u_nodal_index[1]);
1546  if (local_unknown >= 0)
1547  {
1548  jacobian(local_eqn, local_unknown) +=
1549  r * dpsifdx(l2, 1) * testp[l] * W;
1550  }
1551  } /*End of loop over l2*/
1552  } /*End of Jacobian calculation*/
1553 
1554  } // End of if not boundary condition
1555 
1556  } // End of loop over l
1557  }
1558  }
oomph::AxisymmetricNavierStokesEquations::re_invro
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
Definition: axisym_navier_stokes_elements.h:434
oomph::AxisymmetricNavierStokesEquations::get_source_fct
double get_source_fct(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
Definition: axisym_navier_stokes_elements.h:286
oomph::AxisymmetricNavierStokesEquations::re_st
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: axisym_navier_stokes_elements.h:404
oomph::AxisymmetricNavierStokesEquations::re_invfr
const double & re_invfr() const
Global inverse Froude number.
Definition: axisym_navier_stokes_elements.h:422
oomph::AxisymmetricNavierStokesEquations::p_local_eqn
virtual int p_local_eqn(const unsigned &n) const =0
oomph::AxisymmetricNavierStokesEquations::get_body_force_axi_nst
virtual void get_body_force_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force fct at a given time and Eulerian position.
Definition: axisym_navier_stokes_elements.h:222
oomph::AxisymmetricNavierStokesEquations::Gamma
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition: axisym_navier_stokes_elements.h:393
oomph::TimeStepper::time_pt
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
oomph::Time::time
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
Global_Parameters::body_force
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Definition: axisym_linear_elasticity/cylinder/cylinder.cc:96
TestProblem::source
void source(const Vector< double > &x, Vector< double > &f)
Source function.
Definition: unstructured_two_d_circle.cc:46

References ALE_is_disabled, Global_Parameters::body_force(), density_ratio(), dshape_and_dtest_eulerian_at_knot_axi_nst(), du_dt_axi_nst(), 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(), GlobalParameters::Nintpt, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::node_pt(), npres_axi_nst(), oomph::Integral::nweight(), p_axi_nst(), p_local_eqn(), 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, TestProblem::source(), oomph::Time::time(), oomph::TimeStepper::time_pt(), 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_jacobian(), fill_in_contribution_to_jacobian_and_mass_matrix(), and fill_in_contribution_to_residuals().

◆ g()

const Vector<double>& oomph::AxisymmetricNavierStokesEquations::g ( ) const
inline

Vector of gravitational components.

447  {
448  return *G_pt;
449  }

References G_pt.

Referenced by fill_in_contribution_to_hessian_vector_products(), fill_in_generic_dresidual_contribution_axi_nst(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::fill_in_off_diagonal_jacobian_blocks_analytic(), RefineableQAxisymCrouzeixRaviartBoussinesqElement::get_body_force_axi_nst(), oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::get_body_force_axi_nst(), QAxisymCrouzeixRaviartElementWithExternalElement::get_body_force_axi_nst(), RefineableQAxisymCrouzeixRaviartBoussinesqElement::get_dbody_force_axi_nst_dexternal_element_data(), QAxisymCrouzeixRaviartElementWithExternalElement::get_dbody_force_axi_nst_dexternal_element_data(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ g_pt()

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

Pointer to Vector of gravitational components.

453  {
454  return G_pt;
455  }

References G_pt.

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

◆ get_body_force_axi_nst()

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

Reimplemented in QAxisymCrouzeixRaviartElementWithExternalElement, oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement, QAxisymCrouzeixRaviartElementWithExternalElement, oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement, RefineableQAxisymCrouzeixRaviartBoussinesqElement, and RefineableQAxisymCrouzeixRaviartBoussinesqElement.

227  {
228  // If the function pointer is zero return zero
229  if (Body_force_fct_pt == 0)
230  {
231  // Loop over dimensions and set body forces to zero
232  for (unsigned i = 0; i < 3; i++)
233  {
234  result[i] = 0.0;
235  }
236  }
237  // Otherwise call the function
238  else
239  {
240  (*Body_force_fct_pt)(time, x, result);
241  }
242  }

References Body_force_fct_pt, i, and plotDoE::x.

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

◆ get_body_force_gradient_axi_nst()

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

253  {
254  // hierher: Implement function pointer version
255  /* //If no gradient function has been set, FD it */
256  /* if(Body_force_fct_gradient_pt==0) */
257  {
258  // Reference value
259  Vector<double> body_force(3, 0.0);
260  get_body_force_axi_nst(time, ipt, s, x, body_force);
261 
262  // FD it
263  const double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
264  Vector<double> body_force_pls(3, 0.0);
265  Vector<double> x_pls(x);
266  for (unsigned i = 0; i < 2; i++)
267  {
268  x_pls[i] += eps_fd;
269  get_body_force_axi_nst(time, ipt, s, x_pls, body_force_pls);
270  for (unsigned j = 0; j < 3; j++)
271  {
272  d_body_force_dx(j, i) =
273  (body_force_pls[j] - body_force[j]) / eps_fd;
274  }
275  x_pls[i] = x[i];
276  }
277  }
278  /* else */
279  /* { */
280  /* // Get gradient */
281  /* (*Source_fct_gradient_pt)(time,x,gradient); */
282  /* } */
283  }
oomph::GeneralisedElement::Default_fd_jacobian_step
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::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ get_dresidual_dnodal_coordinates()

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

1569  {
1570  // Return immediately if there are no dofs
1571  if (ndof() == 0)
1572  {
1573  return;
1574  }
1575 
1576  // Get continuous time from timestepper of first node
1577  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1578 
1579  // Determine number of nodes in element
1580  const unsigned n_node = nnode();
1581 
1582  // Determine number of pressure dofs in element
1583  const unsigned n_pres = npres_axi_nst();
1584 
1585  // Find the indices at which the local velocities are stored
1586  unsigned u_nodal_index[3];
1587  for (unsigned i = 0; i < 3; i++)
1588  {
1589  u_nodal_index[i] = u_index_axi_nst(i);
1590  }
1591 
1592  // Set up memory for the shape and test functions
1593  // Note that there are only two dimensions, r and z, in this problem
1594  Shape psif(n_node), testf(n_node);
1595  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
1596 
1597  // Set up memory for pressure shape and test functions
1598  Shape psip(n_pres), testp(n_pres);
1599 
1600  // Deriatives of shape fct derivatives w.r.t. nodal coords
1601  RankFourTensor<double> d_dpsifdx_dX(2, n_node, n_node, 2);
1602  RankFourTensor<double> d_dtestfdx_dX(2, n_node, n_node, 2);
1603 
1604  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1605  DenseMatrix<double> dJ_dX(2, n_node);
1606 
1607  // Derivatives of derivative of u w.r.t. nodal coords
1608  // Note that the entry d_dudx_dX(p,q,i,k) contains the derivative w.r.t.
1609  // nodal coordinate X_pq of du_i/dx_k. Since there are three velocity
1610  // components, i can take on the values 0, 1 and 2, while k and p only
1611  // take on the values 0 and 1 since there are only two spatial dimensions.
1612  RankFourTensor<double> d_dudx_dX(2, n_node, 3, 2);
1613 
1614  // Derivatives of nodal velocities w.r.t. nodal coords:
1615  // Assumption: Interaction only local via no-slip so
1616  // X_pq only affects U_iq.
1617  // Note that the entry d_U_dX(p,q,i) contains the derivative w.r.t. nodal
1618  // coordinate X_pq of nodal value U_iq. In other words this entry is the
1619  // derivative of the i-th velocity component w.r.t. the p-th spatial
1620  // coordinate, taken at the q-th local node.
1621  RankThreeTensor<double> d_U_dX(2, n_node, 3, 0.0);
1622 
1623  // Determine the number of integration points
1624  const unsigned n_intpt = integral_pt()->nweight();
1625 
1626  // Vector to hold local coordinates (two dimensions)
1627  Vector<double> s(2);
1628 
1629  // Get physical variables from element
1630  // (Reynolds number must be multiplied by the density ratio)
1631  const double scaled_re = re() * density_ratio();
1632  const double scaled_re_st = re_st() * density_ratio();
1633  const double scaled_re_inv_fr = re_invfr() * density_ratio();
1634  const double scaled_re_inv_ro = re_invro() * density_ratio();
1635  const double visc_ratio = viscosity_ratio();
1636  const Vector<double> G = g();
1637 
1638  // FD step
1639  double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
1640 
1641  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1642  // Assumption: Interaction only local via no-slip so
1643  // X_ij only affects U_ij.
1644  bool element_has_node_with_aux_node_update_fct = false;
1645  for (unsigned q = 0; q < n_node; q++)
1646  {
1647  // Get pointer to q-th local node
1648  Node* nod_pt = node_pt(q);
1649 
1650  // Only compute if there's a node-update fct involved
1651  if (nod_pt->has_auxiliary_node_update_fct_pt())
1652  {
1653  element_has_node_with_aux_node_update_fct = true;
1654 
1655  // This functionality has not been tested so issue a warning
1656  // to this effect
1657  std::ostringstream warning_stream;
1658  warning_stream << "\nThe functionality to evaluate the additional"
1659  << "\ncontribution to the deriv of the residual eqn"
1660  << "\nw.r.t. the nodal coordinates which comes about"
1661  << "\nif a node's values are updated using an auxiliary"
1662  << "\nnode update function has NOT been tested for"
1663  << "\naxisymmetric Navier-Stokes elements. Use at your"
1664  << "\nown risk" << std::endl;
1665  OomphLibWarning(
1666  warning_stream.str(),
1667  "AxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates",
1668  OOMPH_EXCEPTION_LOCATION);
1669 
1670  // Current nodal velocity
1671  Vector<double> u_ref(3);
1672  for (unsigned i = 0; i < 3; i++)
1673  {
1674  u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
1675  }
1676 
1677  // FD
1678  for (unsigned p = 0; p < 2; p++)
1679  {
1680  // Make backup
1681  double backup = nod_pt->x(p);
1682 
1683  // Do FD step. No node update required as we're
1684  // attacking the coordinate directly...
1685  nod_pt->x(p) += eps_fd;
1686 
1687  // Do auxiliary node update (to apply no slip)
1688  nod_pt->perform_auxiliary_node_update_fct();
1689 
1690  // Loop over velocity components
1691  for (unsigned i = 0; i < 3; i++)
1692  {
1693  // Evaluate
1694  d_U_dX(p, q, i) =
1695  (*(nod_pt->value_pt(u_nodal_index[i])) - u_ref[i]) / eps_fd;
1696  }
1697 
1698  // Reset
1699  nod_pt->x(p) = backup;
1700 
1701  // Do auxiliary node update (to apply no slip)
1702  nod_pt->perform_auxiliary_node_update_fct();
1703  }
1704  }
1705  }
1706 
1707  // Integer to store the local equation number
1708  int local_eqn = 0;
1709 
1710  // Loop over the integration points
1711  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1712  {
1713  // Assign values of s
1714  for (unsigned i = 0; i < 2; i++)
1715  {
1716  s[i] = integral_pt()->knot(ipt, i);
1717  }
1718 
1719  // Get the integral weight
1720  const double w = integral_pt()->weight(ipt);
1721 
1722  // Call the derivatives of the shape and test functions
1723  const double J = dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,
1724  psif,
1725  dpsifdx,
1726  d_dpsifdx_dX,
1727  testf,
1728  dtestfdx,
1729  d_dtestfdx_dX,
1730  dJ_dX);
1731 
1732  // Call the pressure shape and test functions
1733  pshape_axi_nst(s, psip, testp);
1734 
1735  // Allocate storage for the position and the derivative of the
1736  // mesh positions w.r.t. time
1737  Vector<double> interpolated_x(2, 0.0);
1738  Vector<double> mesh_velocity(2, 0.0);
1739 
1740  // Allocate storage for the pressure, velocity components and their
1741  // time and spatial derivatives
1742  double interpolated_p = 0.0;
1743  Vector<double> interpolated_u(3, 0.0);
1744  Vector<double> dudt(3, 0.0);
1745  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
1746 
1747  // Calculate pressure at integration point
1748  for (unsigned l = 0; l < n_pres; l++)
1749  {
1750  interpolated_p += p_axi_nst(l) * psip[l];
1751  }
1752 
1753  // Calculate velocities and derivatives at integration point
1754  // ---------------------------------------------------------
1755 
1756  // Loop over nodes
1757  for (unsigned l = 0; l < n_node; l++)
1758  {
1759  // Cache the shape function
1760  const double psif_ = psif(l);
1761 
1762  // Loop over the two coordinate directions
1763  for (unsigned i = 0; i < 2; i++)
1764  {
1765  interpolated_x[i] += this->raw_nodal_position(l, i) * psif_;
1766  }
1767 
1768  // Loop over the three velocity directions
1769  for (unsigned i = 0; i < 3; i++)
1770  {
1771  // Get the nodal value
1772  const double u_value = this->raw_nodal_value(l, u_nodal_index[i]);
1773  interpolated_u[i] += u_value * psif_;
1774  dudt[i] += du_dt_axi_nst(l, i) * psif_;
1775 
1776  // Loop over derivative directions
1777  for (unsigned j = 0; j < 2; j++)
1778  {
1779  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1780  }
1781  }
1782  }
1783 
1784  // Get the mesh velocity if ALE is enabled
1785  if (!ALE_is_disabled)
1786  {
1787  // Loop over nodes
1788  for (unsigned l = 0; l < n_node; l++)
1789  {
1790  // Loop over the two coordinate directions
1791  for (unsigned i = 0; i < 2; i++)
1792  {
1793  mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif[l];
1794  }
1795  }
1796  }
1797 
1798  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1799  for (unsigned q = 0; q < n_node; q++)
1800  {
1801  // Loop over the two coordinate directions
1802  for (unsigned p = 0; p < 2; p++)
1803  {
1804  // Loop over the three velocity components
1805  for (unsigned i = 0; i < 3; i++)
1806  {
1807  // Loop over the two coordinate directions
1808  for (unsigned k = 0; k < 2; k++)
1809  {
1810  double aux = 0.0;
1811 
1812  // Loop over nodes and add contribution
1813  for (unsigned j = 0; j < n_node; j++)
1814  {
1815  aux += this->raw_nodal_value(j, u_nodal_index[i]) *
1816  d_dpsifdx_dX(p, q, j, k);
1817  }
1818  d_dudx_dX(p, q, i, k) = aux;
1819  }
1820  }
1821  }
1822  }
1823 
1824  // Get weight of actual nodal position/value in computation of mesh
1825  // velocity from positional/value time stepper
1826  const double pos_time_weight =
1827  node_pt(0)->position_time_stepper_pt()->weight(1, 0);
1828  const double val_time_weight =
1829  node_pt(0)->time_stepper_pt()->weight(1, 0);
1830 
1831  // Get the user-defined body force terms
1832  Vector<double> body_force(3);
1833  get_body_force_axi_nst(time, ipt, s, interpolated_x, body_force);
1834 
1835  // Get the user-defined source function
1836  const double source = get_source_fct(time, ipt, interpolated_x);
1837 
1838  // Note: Can use raw values and avoid bypassing hanging information
1839  // because this is the non-refineable version!
1840 
1841  // Get gradient of body force function
1842  DenseMatrix<double> d_body_force_dx(3, 2, 0.0);
1843  get_body_force_gradient_axi_nst(
1844  time, ipt, s, interpolated_x, d_body_force_dx);
1845 
1846  // Get gradient of source function
1847  Vector<double> source_gradient(2, 0.0);
1848  get_source_fct_gradient(time, ipt, interpolated_x, source_gradient);
1849 
1850  // Cache r (the first position component)
1851  const double r = interpolated_x[0];
1852 
1853  // Assemble shape derivatives
1854  // --------------------------
1855 
1856  // ==================
1857  // MOMENTUM EQUATIONS
1858  // ==================
1859 
1860  // Loop over the test functions
1861  for (unsigned l = 0; l < n_node; l++)
1862  {
1863  // Cache the test function
1864  const double testf_ = testf[l];
1865 
1866  // --------------------------------
1867  // FIRST (RADIAL) MOMENTUM EQUATION
1868  // --------------------------------
1869  local_eqn = nodal_local_eqn(l, u_nodal_index[0]);
1870  ;
1871 
1872  // If it's not a boundary condition
1873  if (local_eqn >= 0)
1874  {
1875  // Loop over the two coordinate directions
1876  for (unsigned p = 0; p < 2; p++)
1877  {
1878  // Loop over nodes
1879  for (unsigned q = 0; q < n_node; q++)
1880  {
1881  // Residual x deriv of Jacobian
1882  // ----------------------------
1883 
1884  // Add the user-defined body force terms
1885  double sum = r * body_force[0] * testf_;
1886 
1887  // Add the gravitational body force term
1888  sum += r * scaled_re_inv_fr * testf_ * G[0];
1889 
1890  // Add the pressure gradient term
1891  sum += interpolated_p * (testf_ + r * dtestfdx(l, 0));
1892 
1893  // Add the stress tensor terms
1894  // The viscosity ratio needs to go in here to ensure
1895  // continuity of normal stress is satisfied even in flows
1896  // with zero pressure gradient!
1897  sum -= visc_ratio * r * (1.0 + Gamma[0]) *
1898  interpolated_dudx(0, 0) * dtestfdx(l, 0);
1899 
1900  sum -=
1901  visc_ratio * r *
1902  (interpolated_dudx(0, 1) + Gamma[0] * interpolated_dudx(1, 0)) *
1903  dtestfdx(l, 1);
1904 
1905  sum -=
1906  visc_ratio * (1.0 + Gamma[0]) * interpolated_u[0] * testf_ / r;
1907 
1908  // Add the du/dt term
1909  sum -= scaled_re_st * r * dudt[0] * testf_;
1910 
1911  // Add the convective terms
1912  sum -= scaled_re *
1913  (r * interpolated_u[0] * interpolated_dudx(0, 0) -
1914  interpolated_u[2] * interpolated_u[2] +
1915  r * interpolated_u[1] * interpolated_dudx(0, 1)) *
1916  testf_;
1917 
1918  // Add the mesh velocity terms
1919  if (!ALE_is_disabled)
1920  {
1921  for (unsigned k = 0; k < 2; k++)
1922  {
1923  sum += scaled_re_st * r * mesh_velocity[k] *
1924  interpolated_dudx(0, k) * testf_;
1925  }
1926  }
1927 
1928  // Add the Coriolis term
1929  sum += 2.0 * r * scaled_re_inv_ro * interpolated_u[2] * testf_;
1930 
1931  // Multiply through by deriv of Jacobian and integration weight
1932  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1933  sum * dJ_dX(p, q) * w;
1934 
1935  // Derivative of residual x Jacobian
1936  // ---------------------------------
1937 
1938  // Body force
1939  sum = r * d_body_force_dx(0, p) * psif[q] * testf_;
1940  if (p == 0)
1941  {
1942  sum += body_force[0] * psif[q] * testf_;
1943  }
1944 
1945  // Gravitational body force
1946  if (p == 0)
1947  {
1948  sum += scaled_re_inv_fr * G[0] * psif[q] * testf_;
1949  }
1950 
1951  // Pressure gradient term
1952  sum += r * interpolated_p * d_dtestfdx_dX(p, q, l, 0);
1953  if (p == 0)
1954  {
1955  sum += interpolated_p * psif[q] * dtestfdx(l, 0);
1956  }
1957 
1958  // Viscous terms
1959  sum -=
1960  r * visc_ratio *
1961  ((1.0 + Gamma[0]) *
1962  (d_dudx_dX(p, q, 0, 0) * dtestfdx(l, 0) +
1963  interpolated_dudx(0, 0) * d_dtestfdx_dX(p, q, l, 0)) +
1964  (d_dudx_dX(p, q, 0, 1) + Gamma[0] * d_dudx_dX(p, q, 1, 0)) *
1965  dtestfdx(l, 1) +
1966  (interpolated_dudx(0, 1) +
1967  Gamma[0] * interpolated_dudx(1, 0)) *
1968  d_dtestfdx_dX(p, q, l, 1));
1969  if (p == 0)
1970  {
1971  sum -= visc_ratio *
1972  ((1.0 + Gamma[0]) *
1973  (interpolated_dudx(0, 0) * psif[q] * dtestfdx(l, 0) -
1974  interpolated_u[0] * psif[q] * testf_ / (r * r)) +
1975  (interpolated_dudx(0, 1) +
1976  Gamma[0] * interpolated_dudx(1, 0)) *
1977  psif[q] * dtestfdx(l, 1));
1978  }
1979 
1980  // Convective terms, including mesh velocity
1981  for (unsigned k = 0; k < 2; k++)
1982  {
1983  double tmp = scaled_re * interpolated_u[k];
1984  if (!ALE_is_disabled)
1985  {
1986  tmp -= scaled_re_st * mesh_velocity[k];
1987  }
1988  sum -= r * tmp * d_dudx_dX(p, q, 0, k) * testf_;
1989  if (p == 0)
1990  {
1991  sum -= tmp * interpolated_dudx(0, k) * psif[q] * testf_;
1992  }
1993  }
1994  if (!ALE_is_disabled)
1995  {
1996  sum += scaled_re_st * r * pos_time_weight *
1997  interpolated_dudx(0, p) * psif[q] * testf_;
1998  }
1999 
2000  // du/dt term
2001  if (p == 0)
2002  {
2003  sum -= scaled_re_st * dudt[0] * psif[q] * testf_;
2004  }
2005 
2006  // Coriolis term
2007  if (p == 0)
2008  {
2009  sum +=
2010  2.0 * scaled_re_inv_ro * interpolated_u[2] * psif[q] * testf_;
2011  }
2012 
2013  // Multiply through by Jacobian and integration weight
2014  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2015 
2016  // Derivs w.r.t. nodal velocities
2017  // ------------------------------
2018  if (element_has_node_with_aux_node_update_fct)
2019  {
2020  // Contribution from deriv of first velocity component
2021  double tmp =
2022  scaled_re_st * r * val_time_weight * psif[q] * testf_;
2023  tmp +=
2024  r * scaled_re * interpolated_dudx(0, 0) * psif[q] * testf_;
2025  for (unsigned k = 0; k < 2; k++)
2026  {
2027  double tmp2 = scaled_re * interpolated_u[k];
2028  if (!ALE_is_disabled)
2029  {
2030  tmp2 -= scaled_re_st * mesh_velocity[k];
2031  }
2032  tmp += r * tmp2 * dpsifdx(q, k) * testf_;
2033  }
2034 
2035  tmp += r * visc_ratio * (1.0 + Gamma[0]) * dpsifdx(q, 0) *
2036  dtestfdx(l, 0);
2037  tmp += r * visc_ratio * dpsifdx(q, 1) * dtestfdx(l, 1);
2038  tmp += visc_ratio * (1.0 + Gamma[0]) * psif[q] * testf_ / r;
2039 
2040  // Multiply through by dU_0q/dX_pq
2041  sum = -d_U_dX(p, q, 0) * tmp;
2042 
2043  // Contribution from deriv of second velocity component
2044  tmp =
2045  scaled_re * r * interpolated_dudx(0, 1) * psif[q] * testf_;
2046  tmp +=
2047  r * visc_ratio * Gamma[0] * dpsifdx(q, 0) * dtestfdx(l, 1);
2048 
2049  // Multiply through by dU_1q/dX_pq
2050  sum -= d_U_dX(p, q, 1) * tmp;
2051 
2052  // Contribution from deriv of third velocity component
2053  tmp = 2.0 *
2054  (r * scaled_re_inv_ro + scaled_re * interpolated_u[2]) *
2055  psif[q] * testf_;
2056 
2057  // Multiply through by dU_2q/dX_pq
2058  sum += d_U_dX(p, q, 2) * tmp;
2059 
2060  // Multiply through by Jacobian and integration weight
2061  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2062  }
2063  } // End of loop over q
2064  } // End of loop over p
2065  } // End of if it's not a boundary condition
2066 
2067  // --------------------------------
2068  // SECOND (AXIAL) MOMENTUM EQUATION
2069  // --------------------------------
2070  local_eqn = nodal_local_eqn(l, u_nodal_index[1]);
2071  ;
2072 
2073  // If it's not a boundary condition
2074  if (local_eqn >= 0)
2075  {
2076  // Loop over the two coordinate directions
2077  for (unsigned p = 0; p < 2; p++)
2078  {
2079  // Loop over nodes
2080  for (unsigned q = 0; q < n_node; q++)
2081  {
2082  // Residual x deriv of Jacobian
2083  // ----------------------------
2084 
2085  // Add the user-defined body force terms
2086  double sum = r * body_force[1] * testf_;
2087 
2088  // Add the gravitational body force term
2089  sum += r * scaled_re_inv_fr * testf_ * G[1];
2090 
2091  // Add the pressure gradient term
2092  sum += r * interpolated_p * dtestfdx(l, 1);
2093 
2094  // Add the stress tensor terms
2095  // The viscosity ratio needs to go in here to ensure
2096  // continuity of normal stress is satisfied even in flows
2097  // with zero pressure gradient!
2098  sum -=
2099  visc_ratio * r *
2100  (interpolated_dudx(1, 0) + Gamma[1] * interpolated_dudx(0, 1)) *
2101  dtestfdx(l, 0);
2102 
2103  sum -= visc_ratio * r * (1.0 + Gamma[1]) *
2104  interpolated_dudx(1, 1) * dtestfdx(l, 1);
2105 
2106  // Add the du/dt term
2107  sum -= scaled_re_st * r * dudt[1] * testf_;
2108 
2109  // Add the convective terms
2110  sum -= scaled_re *
2111  (r * interpolated_u[0] * interpolated_dudx(1, 0) +
2112  r * interpolated_u[1] * interpolated_dudx(1, 1)) *
2113  testf_;
2114 
2115  // Add the mesh velocity terms
2116  if (!ALE_is_disabled)
2117  {
2118  for (unsigned k = 0; k < 2; k++)
2119  {
2120  sum += scaled_re_st * r * mesh_velocity[k] *
2121  interpolated_dudx(1, k) * testf_;
2122  }
2123  }
2124 
2125  // Multiply through by deriv of Jacobian and integration weight
2126  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2127  sum * dJ_dX(p, q) * w;
2128 
2129  // Derivative of residual x Jacobian
2130  // ---------------------------------
2131 
2132  // Body force
2133  sum = r * d_body_force_dx(1, p) * psif[q] * testf_;
2134  if (p == 0)
2135  {
2136  sum += body_force[1] * psif[q] * testf_;
2137  }
2138 
2139  // Gravitational body force
2140  if (p == 0)
2141  {
2142  sum += scaled_re_inv_fr * G[1] * psif[q] * testf_;
2143  }
2144 
2145  // Pressure gradient term
2146  sum += r * interpolated_p * d_dtestfdx_dX(p, q, l, 1);
2147  if (p == 0)
2148  {
2149  sum += interpolated_p * psif[q] * dtestfdx(l, 1);
2150  }
2151 
2152  // Viscous terms
2153  sum -=
2154  r * visc_ratio *
2155  ((d_dudx_dX(p, q, 1, 0) + Gamma[1] * d_dudx_dX(p, q, 0, 1)) *
2156  dtestfdx(l, 0) +
2157  (interpolated_dudx(1, 0) +
2158  Gamma[1] * interpolated_dudx(0, 1)) *
2159  d_dtestfdx_dX(p, q, l, 0) +
2160  (1.0 + Gamma[1]) * d_dudx_dX(p, q, 1, 1) * dtestfdx(l, 1) +
2161  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
2162  d_dtestfdx_dX(p, q, l, 1));
2163  if (p == 0)
2164  {
2165  sum -=
2166  visc_ratio * ((interpolated_dudx(1, 0) +
2167  Gamma[1] * interpolated_dudx(0, 1)) *
2168  psif[q] * dtestfdx(l, 0) +
2169  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
2170  psif[q] * dtestfdx(l, 1));
2171  }
2172 
2173  // Convective terms, including mesh velocity
2174  for (unsigned k = 0; k < 2; k++)
2175  {
2176  double tmp = scaled_re * interpolated_u[k];
2177  if (!ALE_is_disabled)
2178  {
2179  tmp -= scaled_re_st * mesh_velocity[k];
2180  }
2181  sum -= r * tmp * d_dudx_dX(p, q, 1, k) * testf_;
2182  if (p == 0)
2183  {
2184  sum -= tmp * interpolated_dudx(1, k) * psif[q] * testf_;
2185  }
2186  }
2187  if (!ALE_is_disabled)
2188  {
2189  sum += scaled_re_st * r * pos_time_weight *
2190  interpolated_dudx(1, p) * psif[q] * testf_;
2191  }
2192 
2193  // du/dt term
2194  if (p == 0)
2195  {
2196  sum -= scaled_re_st * dudt[1] * psif[q] * testf_;
2197  }
2198 
2199  // Multiply through by Jacobian and integration weight
2200  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2201 
2202  // Derivs w.r.t. to nodal velocities
2203  // ---------------------------------
2204  if (element_has_node_with_aux_node_update_fct)
2205  {
2206  // Contribution from deriv of first velocity component
2207  double tmp =
2208  scaled_re * r * interpolated_dudx(1, 0) * psif[q] * testf_;
2209  tmp +=
2210  r * visc_ratio * Gamma[1] * dpsifdx(q, 1) * dtestfdx(l, 0);
2211 
2212  // Multiply through by dU_0q/dX_pq
2213  sum = -d_U_dX(p, q, 0) * tmp;
2214 
2215  // Contribution from deriv of second velocity component
2216  tmp = scaled_re_st * r * val_time_weight * psif[q] * testf_;
2217  tmp +=
2218  r * scaled_re * interpolated_dudx(1, 1) * psif[q] * testf_;
2219  for (unsigned k = 0; k < 2; k++)
2220  {
2221  double tmp2 = scaled_re * interpolated_u[k];
2222  if (!ALE_is_disabled)
2223  {
2224  tmp2 -= scaled_re_st * mesh_velocity[k];
2225  }
2226  tmp += r * tmp2 * dpsifdx(q, k) * testf_;
2227  }
2228  tmp += r * visc_ratio *
2229  (dpsifdx(q, 0) * dtestfdx(l, 0) +
2230  (1.0 + Gamma[1]) * dpsifdx(q, 1) * dtestfdx(l, 1));
2231 
2232  // Multiply through by dU_1q/dX_pq
2233  sum -= d_U_dX(p, q, 1) * tmp;
2234 
2235  // Multiply through by Jacobian and integration weight
2236  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2237  }
2238  } // End of loop over q
2239  } // End of loop over p
2240  } // End of if it's not a boundary condition
2241 
2242  // -----------------------------------
2243  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
2244  // -----------------------------------
2245  local_eqn = nodal_local_eqn(l, u_nodal_index[2]);
2246  ;
2247 
2248  // If it's not a boundary condition
2249  if (local_eqn >= 0)
2250  {
2251  // Loop over the two coordinate directions
2252  for (unsigned p = 0; p < 2; p++)
2253  {
2254  // Loop over nodes
2255  for (unsigned q = 0; q < n_node; q++)
2256  {
2257  // Residual x deriv of Jacobian
2258  // ----------------------------
2259 
2260  // Add the user-defined body force terms
2261  double sum = r * body_force[2] * testf_;
2262 
2263  // Add the gravitational body force term
2264  sum += r * scaled_re_inv_fr * testf_ * G[2];
2265 
2266  // There is NO pressure gradient term
2267 
2268  // Add in the stress tensor terms
2269  // The viscosity ratio needs to go in here to ensure
2270  // continuity of normal stress is satisfied even in flows
2271  // with zero pressure gradient!
2272  sum -=
2273  visc_ratio *
2274  (r * interpolated_dudx(2, 0) - Gamma[0] * interpolated_u[2]) *
2275  dtestfdx(l, 0);
2276 
2277  sum -= visc_ratio * r * interpolated_dudx(2, 1) * dtestfdx(l, 1);
2278 
2279  sum -=
2280  visc_ratio *
2281  ((interpolated_u[2] / r) - Gamma[0] * interpolated_dudx(2, 0)) *
2282  testf_;
2283 
2284  // Add the du/dt term
2285  sum -= scaled_re_st * r * dudt[2] * testf_;
2286 
2287  // Add the convective terms
2288  sum -= scaled_re *
2289  (r * interpolated_u[0] * interpolated_dudx(2, 0) +
2290  interpolated_u[0] * interpolated_u[2] +
2291  r * interpolated_u[1] * interpolated_dudx(2, 1)) *
2292  testf_;
2293 
2294  // Add the mesh velocity terms
2295  if (!ALE_is_disabled)
2296  {
2297  for (unsigned k = 0; k < 2; k++)
2298  {
2299  sum += scaled_re_st * r * mesh_velocity[k] *
2300  interpolated_dudx(2, k) * testf_;
2301  }
2302  }
2303 
2304  // Add the Coriolis term
2305  sum -= 2.0 * r * scaled_re_inv_ro * interpolated_u[0] * testf_;
2306 
2307  // Multiply through by deriv of Jacobian and integration weight
2308  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2309  sum * dJ_dX(p, q) * w;
2310 
2311  // Derivative of residual x Jacobian
2312  // ---------------------------------
2313 
2314  // Body force
2315  sum = r * d_body_force_dx(2, p) * psif[q] * testf_;
2316  if (p == 0)
2317  {
2318  sum += body_force[2] * psif[q] * testf_;
2319  }
2320 
2321  // Gravitational body force
2322  if (p == 0)
2323  {
2324  sum += scaled_re_inv_fr * G[2] * psif[q] * testf_;
2325  }
2326 
2327  // There is NO pressure gradient term
2328 
2329  // Viscous terms
2330  sum -= visc_ratio * r *
2331  (d_dudx_dX(p, q, 2, 0) * dtestfdx(l, 0) +
2332  interpolated_dudx(2, 0) * d_dtestfdx_dX(p, q, l, 0) +
2333  d_dudx_dX(p, q, 2, 1) * dtestfdx(l, 1) +
2334  interpolated_dudx(2, 1) * d_dtestfdx_dX(p, q, l, 1));
2335 
2336  sum += visc_ratio * Gamma[0] * d_dudx_dX(p, q, 2, 0) * testf_;
2337 
2338  if (p == 0)
2339  {
2340  sum -= visc_ratio *
2341  (interpolated_dudx(2, 0) * psif[q] * dtestfdx(l, 0) +
2342  interpolated_dudx(2, 1) * psif[q] * dtestfdx(l, 1) +
2343  interpolated_u[2] * psif[q] * testf_ / (r * r));
2344  }
2345 
2346  // Convective terms, including mesh velocity
2347  for (unsigned k = 0; k < 2; k++)
2348  {
2349  double tmp = scaled_re * interpolated_u[k];
2350  if (!ALE_is_disabled)
2351  {
2352  tmp -= scaled_re_st * mesh_velocity[k];
2353  }
2354  sum -= r * tmp * d_dudx_dX(p, q, 2, k) * testf_;
2355  if (p == 0)
2356  {
2357  sum -= tmp * interpolated_dudx(2, k) * psif[q] * testf_;
2358  }
2359  }
2360  if (!ALE_is_disabled)
2361  {
2362  sum += scaled_re_st * r * pos_time_weight *
2363  interpolated_dudx(2, p) * psif[q] * testf_;
2364  }
2365 
2366  // du/dt term
2367  if (p == 0)
2368  {
2369  sum -= scaled_re_st * dudt[2] * psif[q] * testf_;
2370  }
2371 
2372  // Coriolis term
2373  if (p == 0)
2374  {
2375  sum -=
2376  2.0 * scaled_re_inv_ro * interpolated_u[0] * psif[q] * testf_;
2377  }
2378 
2379  // Multiply through by Jacobian and integration weight
2380  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2381 
2382  // Derivs w.r.t. to nodal velocities
2383  // ---------------------------------
2384  if (element_has_node_with_aux_node_update_fct)
2385  {
2386  // Contribution from deriv of first velocity component
2387  double tmp = (2.0 * r * scaled_re_inv_ro +
2388  r * scaled_re * interpolated_dudx(2, 0) -
2389  scaled_re * interpolated_u[2]) *
2390  psif[q] * testf_;
2391 
2392  // Multiply through by dU_0q/dX_pq
2393  sum = -d_U_dX(p, q, 0) * tmp;
2394 
2395  // Contribution from deriv of second velocity component
2396  // Multiply through by dU_1q/dX_pq
2397  sum -= d_U_dX(p, q, 1) * scaled_re * r *
2398  interpolated_dudx(2, 1) * psif[q] * testf_;
2399 
2400  // Contribution from deriv of third velocity component
2401  tmp = scaled_re_st * r * val_time_weight * psif[q] * testf_;
2402  tmp -= scaled_re * interpolated_u[0] * psif[q] * testf_;
2403  for (unsigned k = 0; k < 2; k++)
2404  {
2405  double tmp2 = scaled_re * interpolated_u[k];
2406  if (!ALE_is_disabled)
2407  {
2408  tmp2 -= scaled_re_st * mesh_velocity[k];
2409  }
2410  tmp += r * tmp2 * dpsifdx(q, k) * testf_;
2411  }
2412  tmp += visc_ratio * (r * dpsifdx(q, 0) - Gamma[0] * psif[q]) *
2413  dtestfdx(l, 0);
2414  tmp += r * visc_ratio * dpsifdx(q, 1) * dtestfdx(l, 1);
2415  tmp += visc_ratio * ((psif[q] / r) - Gamma[0] * dpsifdx(q, 0)) *
2416  testf_;
2417 
2418  // Multiply through by dU_2q/dX_pq
2419  sum -= d_U_dX(p, q, 2) * tmp;
2420 
2421  // Multiply through by Jacobian and integration weight
2422  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2423  }
2424  } // End of loop over q
2425  } // End of loop over p
2426  } // End of if it's not a boundary condition
2427 
2428  } // End of loop over test functions
2429 
2430 
2431  // ===================
2432  // CONTINUITY EQUATION
2433  // ===================
2434 
2435  // Loop over the shape functions
2436  for (unsigned l = 0; l < n_pres; l++)
2437  {
2438  local_eqn = p_local_eqn(l);
2439 
2440  // Cache the test function
2441  const double testp_ = testp[l];
2442 
2443  // If not a boundary conditions
2444  if (local_eqn >= 0)
2445  {
2446  // Loop over the two coordinate directions
2447  for (unsigned p = 0; p < 2; p++)
2448  {
2449  // Loop over nodes
2450  for (unsigned q = 0; q < n_node; q++)
2451  {
2452  // Residual x deriv of Jacobian
2453  //-----------------------------
2454 
2455  // Source term
2456  double aux = -r * source;
2457 
2458  // Gradient terms
2459  aux += (interpolated_u[0] + r * interpolated_dudx(0, 0) +
2460  r * interpolated_dudx(1, 1));
2461 
2462  // Multiply through by deriv of Jacobian and integration weight
2463  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2464  aux * dJ_dX(p, q) * testp_ * w;
2465 
2466  // Derivative of residual x Jacobian
2467  // ---------------------------------
2468 
2469  // Gradient of source function
2470  aux = -r * source_gradient[p] * psif[q];
2471  if (p == 0)
2472  {
2473  aux -= source * psif[q];
2474  }
2475 
2476  // Gradient terms
2477  aux += r * (d_dudx_dX(p, q, 0, 0) + d_dudx_dX(p, q, 1, 1));
2478  if (p == 0)
2479  {
2480  aux +=
2481  (interpolated_dudx(0, 0) + interpolated_dudx(1, 1)) * psif[q];
2482  }
2483 
2484  // Derivs w.r.t. nodal velocities
2485  // ------------------------------
2486  if (element_has_node_with_aux_node_update_fct)
2487  {
2488  aux += d_U_dX(p, q, 0) * (psif[q] + r * dpsifdx(q, 0));
2489  aux += d_U_dX(p, q, 1) * r * dpsifdx(q, 1);
2490  }
2491 
2492  // Multiply through by Jacobian, test fct and integration weight
2493  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2494  aux * testp_ * J * w;
2495  }
2496  }
2497  }
2498  } // End of loop over shape functions for continuity eqn
2499 
2500  } // End of loop over integration points
2501  }
p
float * p
Definition: Tutorial_Map_using.cpp:9
oomph::AxisymmetricNavierStokesEquations::viscosity_ratio
const double & viscosity_ratio() const
Definition: axisym_navier_stokes_elements.h:472
oomph::AxisymmetricNavierStokesEquations::get_source_fct_gradient
virtual void get_source_fct_gradient(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
Definition: axisym_navier_stokes_elements.h:305
oomph::AxisymmetricNavierStokesEquations::get_body_force_gradient_axi_nst
virtual void get_body_force_gradient_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
Definition: axisym_navier_stokes_elements.h:247
oomph::Node::position_time_stepper_pt
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
tmp
Eigen::Matrix< Scalar, Dynamic, Dynamic, ColMajor > tmp
Definition: level3_impl.h:365
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019

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_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::AxisymmetricNavierStokesEquations::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.

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

Calculate the source fct at given time and Eulerian position.

289  {
290  // If the function pointer is zero return zero
291  if (Source_fct_pt == 0)
292  {
293  return 0;
294  }
295 
296  // Otherwise call the function
297  else
298  {
299  return (*Source_fct_pt)(time, x);
300  }
301  }

References Source_fct_pt, and plotDoE::x.

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

◆ get_source_fct_gradient()

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

309  {
310  // hierher: Implement function pointer version
311  /* //If no gradient function has been set, FD it */
312  /* if(Source_fct_gradient_pt==0) */
313  {
314  // Reference value
315  const double source = get_source_fct(time, ipt, x);
316 
317  // FD it
318  const double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
319  double source_pls = 0.0;
320  Vector<double> x_pls(x);
321  for (unsigned i = 0; i < 2; i++)
322  {
323  x_pls[i] += eps_fd;
324  source_pls = get_source_fct(time, ipt, x_pls);
325  gradient[i] = (source_pls - source) / eps_fd;
326  x_pls[i] = x[i];
327  }
328  }
329  /* else */
330  /* { */
331  /* // Get gradient */
332  /* (*Source_fct_gradient_pt)(time,x,gradient); */
333  /* } */
334  }

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

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

◆ get_viscosity_ratio_axisym_nst()

virtual void oomph::AxisymmetricNavierStokesEquations::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

343  {
344  visc_ratio = *Viscosity_Ratio_pt;
345  }

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::AxisymmetricNavierStokesEquations::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

1119  {
1120  // Determine number of nodes
1121  const unsigned n_node = nnode();
1122 
1123  // Allocate storage for local shape function and its derivatives
1124  // with respect to space
1125  Shape psif(n_node);
1126  DShape dpsifds(n_node, 2);
1127 
1128  // Find values of shape function (ignore jacobian)
1129  (void)this->dshape_local(s, psif, dpsifds);
1130 
1131  // Allocate memory for the jacobian and the inverse of the jacobian
1132  DenseMatrix<double> jacobian(2), inverse_jacobian(2);
1133 
1134  // Allocate memory for the derivative of the jacobian w.r.t. nodal coords
1135  DenseMatrix<double> djacobian_dX(2, n_node);
1136 
1137  // Allocate memory for the derivative w.r.t. nodal coords of the
1138  // spatial derivatives of the shape functions
1139  RankFourTensor<double> d_dpsifdx_dX(2, n_node, 3, 2);
1140 
1141  // Now calculate the inverse jacobian
1142  const double det =
1143  local_to_eulerian_mapping(dpsifds, jacobian, inverse_jacobian);
1144 
1145  // Calculate the derivative of the jacobian w.r.t. nodal coordinates
1146  // Note: must call this before "transform_derivatives(...)" since this
1147  // function requires dpsids rather than dpsidx
1148  dJ_eulerian_dnodal_coordinates(jacobian, dpsifds, djacobian_dX);
1149 
1150  // Calculate the derivative of dpsidx w.r.t. nodal coordinates
1151  // Note: this function also requires dpsids rather than dpsidx
1152  d_dshape_eulerian_dnodal_coordinates(
1153  det, jacobian, djacobian_dX, inverse_jacobian, dpsifds, d_dpsifdx_dX);
1154 
1155  // Get the index at which the velocity is stored
1156  const unsigned u_nodal_index = u_index_axi_nst(i);
1157 
1158  // Initialise value of dudx
1159  double interpolated_d_dudx_dX = 0.0;
1160 
1161  // Loop over the local nodes and sum
1162  for (unsigned l = 0; l < n_node; l++)
1163  {
1164  interpolated_d_dudx_dX +=
1165  nodal_value(l, u_nodal_index) * d_dpsifdx_dX(p, q, l, k);
1166  }
1167 
1168  return (interpolated_d_dudx_dX);
1169  }
oomph::FiniteElement::d_dshape_eulerian_dnodal_coordinates
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
oomph::FiniteElement::dJ_eulerian_dnodal_coordinates
virtual void dJ_eulerian_dnodal_coordinates(const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
Definition: elements.cc:2669
oomph::FiniteElement::local_to_eulerian_mapping
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Definition: elements.h:1508
oomph::FiniteElement::dshape_local
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::AxisymmetricNavierStokesEquations::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

1026  {
1027  // Determine number of nodes
1028  const unsigned n_node = nnode();
1029 
1030  // Allocate storage for local shape function and its derivatives
1031  // with respect to space
1032  Shape psif(n_node);
1033  DShape dpsifds(n_node, 2);
1034 
1035  // Find values of shape function (ignore jacobian)
1036  (void)this->dshape_local(s, psif, dpsifds);
1037 
1038  // Get the index at which the velocity is stored
1039  const unsigned u_nodal_index = u_index_axi_nst(i);
1040 
1041  // Initialise value of duds
1042  double interpolated_duds = 0.0;
1043 
1044  // Loop over the local nodes and sum
1045  for (unsigned l = 0; l < n_node; l++)
1046  {
1047  interpolated_duds += nodal_value(l, u_nodal_index) * dpsifds(l, j);
1048  }
1049 
1050  return (interpolated_duds);
1051  }

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

Referenced by LinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement::get_base_flow_duds().

◆ interpolated_dudt_axi_nst()

double oomph::AxisymmetricNavierStokesEquations::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

1089  {
1090  // Determine number of nodes
1091  const unsigned n_node = nnode();
1092 
1093  // Allocate storage for local shape function
1094  Shape psif(n_node);
1095 
1096  // Find values of shape function
1097  shape(s, psif);
1098 
1099  // Initialise value of dudt
1100  double interpolated_dudt = 0.0;
1101 
1102  // Loop over the local nodes and sum
1103  for (unsigned l = 0; l < n_node; l++)
1104  {
1105  interpolated_dudt += du_dt_axi_nst(l, i) * psif[l];
1106  }
1107 
1108  return (interpolated_dudt);
1109  }

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

Referenced by LinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement::get_base_flow_dudt().

◆ interpolated_dudx_axi_nst()

double oomph::AxisymmetricNavierStokesEquations::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

1058  {
1059  // Determine number of nodes
1060  const unsigned n_node = nnode();
1061 
1062  // Allocate storage for local shape function and its derivatives
1063  // with respect to space
1064  Shape psif(n_node);
1065  DShape dpsifdx(n_node, 2);
1066 
1067  // Find values of shape function (ignore jacobian)
1068  (void)this->dshape_eulerian(s, psif, dpsifdx);
1069 
1070  // Get the index at which the velocity is stored
1071  const unsigned u_nodal_index = u_index_axi_nst(i);
1072 
1073  // Initialise value of dudx
1074  double interpolated_dudx = 0.0;
1075 
1076  // Loop over the local nodes and sum
1077  for (unsigned l = 0; l < n_node; l++)
1078  {
1079  interpolated_dudx += nodal_value(l, u_nodal_index) * dpsifdx(l, j);
1080  }
1081 
1082  return (interpolated_dudx);
1083  }
oomph::FiniteElement::dshape_eulerian
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().

Referenced by LinearisedAxisymmetricQTaylorHoodMultiDomainElement::get_base_flow_dudx(), LinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement::get_base_flow_dudx(), RefineableLinearisedAxisymmetricQTaylorHoodMultiDomainElement::get_base_flow_dudx(), and RefineableLinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement::get_base_flow_dudx().

◆ interpolated_p_axi_nst()

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

Return FE interpolated pressure at local coordinate s.

1002  {
1003  // Find number of nodes
1004  unsigned n_pres = npres_axi_nst();
1005  // Local shape function
1006  Shape psi(n_pres);
1007  // Find values of shape function
1008  pshape_axi_nst(s, psi);
1009 
1010  // Initialise value of p
1011  double interpolated_p = 0.0;
1012  // Loop over the local nodes and sum
1013  for (unsigned l = 0; l < n_pres; l++)
1014  {
1015  interpolated_p += p_axi_nst(l) * psi[l];
1016  }
1017 
1018  return (interpolated_p);
1019  }

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

Referenced by oomph::RefineableAxisymmetricQCrouzeixRaviartElement::further_build(), LinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement::get_base_flow_p(), oomph::RefineableAxisymmetricQTaylorHoodElement::get_interpolated_values(), oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::output(), oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::output(), output(), point_output_data(), pressure_integral(), scalar_value_paraview(), and traction().

◆ interpolated_u_axi_nst() [1/3]

double oomph::AxisymmetricNavierStokesEquations::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.

920  {
921  // Find number of nodes
922  unsigned n_node = nnode();
923  // Local shape function
924  Shape psi(n_node);
925  // Find values of shape function
926  shape(s, psi);
927 
928  // Get the index at which the velocity is stored
929  unsigned u_nodal_index = u_index_axi_nst(i);
930 
931  // Initialise value of u
932  double interpolated_u = 0.0;
933  // Loop over the local nodes and sum
934  for (unsigned l = 0; l < n_node; l++)
935  {
936  interpolated_u += nodal_value(t, l, u_nodal_index) * psi[l];
937  }
938 
939  return (interpolated_u);
940  }

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::AxisymmetricNavierStokesEquations::interpolated_u_axi_nst ( const Vector< double > &  s,
const unsigned &  i 
) const
inline

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

892  {
893  // Find number of nodes
894  unsigned n_node = nnode();
895  // Local shape function
896  Shape psi(n_node);
897  // Find values of shape function
898  shape(s, psi);
899 
900  // Get the index at which the velocity is stored
901  unsigned u_nodal_index = u_index_axi_nst(i);
902 
903  // Initialise value of u
904  double interpolated_u = 0.0;
905  // Loop over the local nodes and sum
906  for (unsigned l = 0; l < n_node; l++)
907  {
908  interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
909  }
910 
911  return (interpolated_u);
912  }

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::AxisymmetricNavierStokesEquations::interpolated_u_axi_nst ( const Vector< double > &  s,
Vector< double > &  veloc 
) const
inline

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

867  {
868  // Find number of nodes
869  unsigned n_node = nnode();
870  // Local shape function
871  Shape psi(n_node);
872  // Find values of shape function
873  shape(s, psi);
874 
875  for (unsigned i = 0; i < 3; i++)
876  {
877  // Index at which the nodal value is stored
878  unsigned u_nodal_index = u_index_axi_nst(i);
879  // Initialise value of u
880  veloc[i] = 0.0;
881  // Loop over the local nodes and sum
882  for (unsigned l = 0; l < n_node; l++)
883  {
884  veloc[i] += nodal_value(l, u_nodal_index) * psi[l];
885  }
886  }
887  }

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

Referenced by compute_error(), LinearisedAxisymmetricQTaylorHoodMultiDomainElement::get_base_flow_u(), LinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement::get_base_flow_u(), RefineableLinearisedAxisymmetricQTaylorHoodMultiDomainElement::get_base_flow_u(), RefineableLinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement::get_base_flow_u(), oomph::RefineableAxisymmetricQTaylorHoodElement::get_interpolated_values(), oomph::RefineableAxisymmetricQCrouzeixRaviartElement::get_interpolated_values(), oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::get_wind_axi_adv_diff(), RefineableQAxisymAdvectionDiffusionBoussinesqElement::get_wind_axi_adv_diff(), QAxisymAdvectionDiffusionElementWithExternalElement::get_wind_axi_adv_diff(), oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::get_wind_axi_adv_diff(), kin_energy(), oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::output(), oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::output(), output(), point_output_data(), and scalar_value_paraview().

◆ kin_energy()

double oomph::AxisymmetricNavierStokesEquations::kin_energy ( ) const

Get integral of kinetic energy over element.

Get integral of kinetic energy over element:

822  {
823  throw OomphLibError(
824  "Check the kinetic energy calculation for axisymmetric NSt",
825  OOMPH_CURRENT_FUNCTION,
826  OOMPH_EXCEPTION_LOCATION);
827 
828  // Initialise
829  double kin_en = 0.0;
830 
831  // Set the value of Nintpt
832  unsigned Nintpt = integral_pt()->nweight();
833 
834  // Set the Vector to hold local coordinates
835  Vector<double> s(2);
836 
837  // Loop over the integration points
838  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
839  {
840  // Assign values of s
841  for (unsigned i = 0; i < 2; i++)
842  {
843  s[i] = integral_pt()->knot(ipt, i);
844  }
845 
846  // Get the integral weight
847  double w = integral_pt()->weight(ipt);
848 
849  // Get Jacobian of mapping
850  double J = J_eulerian(s);
851 
852  // Loop over directions
853  double veloc_squared = 0.0;
854  for (unsigned i = 0; i < 3; i++)
855  {
856  veloc_squared +=
857  interpolated_u_axi_nst(s, i) * interpolated_u_axi_nst(s, i);
858  }
859 
860  kin_en += 0.5 * veloc_squared * w * J * interpolated_x(s, 0);
861  }
862 
863  return kin_en;
864  }

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().

◆ n_u_nst()

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

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

524  {
525  return 3;
526  }

◆ npres_axi_nst()

virtual unsigned oomph::AxisymmetricNavierStokesEquations::npres_axi_nst ( ) const
pure virtual

Function to return number of pressure degrees of freedom.

Implemented in oomph::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

Referenced by fill_in_generic_dresidual_contribution_axi_nst(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates(), and interpolated_p_axi_nst().

◆ nscalar_paraview()

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

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

Reimplemented from oomph::FiniteElement.

620  {
621  return 4;
622  }

Referenced by oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::nscalar_paraview().

◆ output() [1/4]

void oomph::AxisymmetricNavierStokesEquations::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< AxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >, oomph::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

729  {
730  unsigned nplot = 5;
731  output(file_pt, nplot);
732  }
oomph::AxisymmetricNavierStokesEquations::output
void output(std::ostream &outfile)
Definition: axisym_navier_stokes_elements.h:715

References output().

◆ output() [2/4]

void oomph::AxisymmetricNavierStokesEquations::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::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, oomph::AxisymmetricQCrouzeixRaviartElement, and oomph::PseudoSolidNodeUpdateElement< AxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >.

548  {
549  // Vector of local coordinates
550  Vector<double> s(2);
551 
552  // Tecplot header info
553  fprintf(file_pt, "%s ", tecplot_zone_string(nplot).c_str());
554 
555  // Loop over plot points
556  unsigned num_plot_points = nplot_points(nplot);
557  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
558  {
559  // Get local coordinates of plot point
560  get_s_plot(iplot, nplot, s);
561 
562  // Coordinates
563  for (unsigned i = 0; i < 2; i++)
564  {
565  // outfile << interpolated_x(s,i) << " ";
566  fprintf(file_pt, "%g ", interpolated_x(s, i));
567  }
568 
569  // Velocities
570  for (unsigned i = 0; i < 3; i++)
571  {
572  // outfile << interpolated_u(s,i) << " ";
573  fprintf(file_pt, "%g ", interpolated_u_axi_nst(s, i));
574  }
575 
576  // Pressure
577  // outfile << interpolated_p(s) << " ";
578  fprintf(file_pt, "%g ", interpolated_p_axi_nst(s));
579 
580  // outfile << std::endl;
581  fprintf(file_pt, "\n");
582  }
583  // outfile << std::endl;
584  fprintf(file_pt, "\n");
585 
586  // Write tecplot footer (e.g. FE connectivity lists)
587  write_tecplot_zone_footer(file_pt, nplot);
588  }
oomph::AxisymmetricNavierStokesEquations::interpolated_p_axi_nst
double interpolated_p_axi_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
Definition: axisym_navier_stokes_elements.h:1001
oomph::FiniteElement::tecplot_zone_string
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Definition: elements.h:3161
oomph::FiniteElement::get_s_plot
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
oomph::FiniteElement::nplot_points
virtual unsigned nplot_points(const unsigned &nplot) const
Definition: elements.h:3186
oomph::FiniteElement::write_tecplot_zone_footer
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Definition: elements.h:3174

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]

void oomph::AxisymmetricNavierStokesEquations::output ( std::ostream &  outfile)
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< AxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >, oomph::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

716  {
717  unsigned nplot = 5;
718  output(outfile, nplot);
719  }

Referenced by output(), oomph::AxisymmetricQCrouzeixRaviartElement::output(), oomph::AxisymmetricQTaylorHoodElement::output(), oomph::AxisymmetricTCrouzeixRaviartElement::output(), and oomph::AxisymmetricTTaylorHoodElement::output().

◆ output() [4/4]

void oomph::AxisymmetricNavierStokesEquations::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::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, oomph::AxisymmetricQCrouzeixRaviartElement, and oomph::PseudoSolidNodeUpdateElement< AxisymmetricTTaylorHoodElement, TPVDElement< 2, 3 > >.

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

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::AxisymmetricNavierStokesEquations::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.

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

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::AxisymmetricNavierStokesEquations::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.

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

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::AxisymmetricNavierStokesEquations::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.

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

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

virtual double oomph::AxisymmetricNavierStokesEquations::p_axi_nst ( const unsigned &  n_p) const
pure virtual

Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes.

Implemented in oomph::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricQTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

Referenced by fill_in_generic_dresidual_contribution_axi_nst(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates(), and interpolated_p_axi_nst().

◆ p_local_eqn()

virtual int oomph::AxisymmetricNavierStokesEquations::p_local_eqn ( const unsigned &  n) const
protectedpure virtual

Access function for the local equation number information for the pressure. p_local_eqn[n] = local equation number or < 0 if pinned

Implemented in oomph::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

Referenced by fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ p_nodal_index_axi_nst()

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

Which nodal value represents the pressure?

Reimplemented in oomph::AxisymmetricTTaylorHoodElement, and oomph::AxisymmetricQTaylorHoodElement.

580  {
581  return Pressure_not_stored_at_node;
582  }
oomph::AxisymmetricNavierStokesEquations::Pressure_not_stored_at_node
static int Pressure_not_stored_at_node
Definition: axisym_navier_stokes_elements.h:119

References Pressure_not_stored_at_node.

Referenced by oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ point_output_data()

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

695  {
696  // Output the components of the position
697  for (unsigned i = 0; i < 2; i++)
698  {
699  data.push_back(interpolated_x(s, i));
700  }
701 
702  // Output the components of the FE representation of u at s
703  for (unsigned i = 0; i < 3; i++)
704  {
705  data.push_back(interpolated_u_axi_nst(s, i));
706  }
707 
708  // Output FE representation of p at s
709  data.push_back(interpolated_p_axi_nst(s));
710  }
data
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::AxisymmetricNavierStokesEquations::pressure_integral ( ) const

Integral of pressure over element.

Return pressure integrated over the element.

870  {
871  // Initialise
872  double press_int = 0;
873 
874  // Set the value of Nintpt
875  unsigned Nintpt = integral_pt()->nweight();
876 
877  // Set the Vector to hold local coordinates
878  Vector<double> s(2);
879 
880  // Loop over the integration points
881  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
882  {
883  // Assign values of s
884  for (unsigned i = 0; i < 2; i++)
885  {
886  s[i] = integral_pt()->knot(ipt, i);
887  }
888 
889  // Get the integral weight
890  double w = integral_pt()->weight(ipt);
891 
892  // Get Jacobian of mapping
893  double J = J_eulerian(s);
894 
895  // Premultiply the weights and the Jacobian
896  double W = w * J * interpolated_x(s, 0);
897 
898  // Get pressure
899  double press = interpolated_p_axi_nst(s);
900 
901  // Add
902  press_int += press * W;
903  }
904 
905  return press_int;
906  }

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]

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

Compute the pressure shape functions at local coordinate s.

Implemented in oomph::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

Referenced by fill_in_generic_dresidual_contribution_axi_nst(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates(), and interpolated_p_axi_nst().

◆ pshape_axi_nst() [2/2]

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

Compute the pressure shape and test functions at local coordinate s

Implemented in oomph::AxisymmetricTTaylorHoodElement, oomph::AxisymmetricTCrouzeixRaviartElement, oomph::AxisymmetricQTaylorHoodElement, and oomph::AxisymmetricQCrouzeixRaviartElement.

◆ re()

const double& oomph::AxisymmetricNavierStokesEquations::re ( ) const
inline

Reynolds number.

399  {
400  return *Re_pt;
401  }

References Re_pt.

Referenced by fill_in_contribution_to_hessian_vector_products(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ re_invfr()

const double& oomph::AxisymmetricNavierStokesEquations::re_invfr ( ) const
inline

Global inverse Froude number.

423  {
424  return *ReInvFr_pt;
425  }

References ReInvFr_pt.

Referenced by fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ re_invfr_pt()

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

Pointer to global inverse Froude number.

429  {
430  return ReInvFr_pt;
431  }

References ReInvFr_pt.

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

◆ re_invro()

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

Global Reynolds number multiplied by inverse Rossby number.

435  {
436  return *ReInvRo_pt;
437  }

References ReInvRo_pt.

Referenced by fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ re_invro_pt()

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

Pointer to global inverse Froude number.

441  {
442  return ReInvRo_pt;
443  }

References ReInvRo_pt.

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

◆ re_pt()

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

Pointer to Reynolds number.

411  {
412  return Re_pt;
413  }

References Re_pt.

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

◆ re_st()

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

Product of Reynolds and Strouhal number (=Womersley number)

405  {
406  return *ReSt_pt;
407  }

References ReSt_pt.

Referenced by fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ re_st_pt()

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

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

417  {
418  return ReSt_pt;
419  }

References ReSt_pt.

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

◆ scalar_name_paraview()

std::string oomph::AxisymmetricNavierStokesEquations::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.

668  {
669  // Veloc
670  if (i < 3)
671  {
672  return "Velocity " + StringConversion::to_string(i);
673  }
674  // Pressure field
675  else if (i == 3)
676  {
677  return "Pressure";
678  }
679  // Never get here
680  else
681  {
682  std::stringstream error_stream;
683  error_stream
684  << "Axisymmetric Navier-Stokes Elements only store 4 fields "
685  << std::endl;
686  throw OomphLibError(
687  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
688  return " ";
689  }
690  }
oomph::StringConversion::to_string
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().

Referenced by oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::scalar_name_paraview().

◆ scalar_value_paraview()

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

629  {
630  // Vector of local coordinates
631  Vector<double> s(2);
632 
633  // Loop over plot points
634  unsigned num_plot_points = nplot_points_paraview(nplot);
635  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
636  {
637  // Get local coordinates of plot point
638  get_s_plot(iplot, nplot, s);
639 
640  // Velocities
641  if (i < 3)
642  {
643  file_out << interpolated_u_axi_nst(s, i) << std::endl;
644  }
645  // Pressure
646  else if (i == 3)
647  {
648  file_out << interpolated_p_axi_nst(s) << std::endl;
649  }
650  // Never get here
651  else
652  {
653  std::stringstream error_stream;
654  error_stream
655  << "Axisymmetric Navier-Stokes Elements only store 4 fields "
656  << std::endl;
657  throw OomphLibError(error_stream.str(),
658  OOMPH_CURRENT_FUNCTION,
659  OOMPH_EXCEPTION_LOCATION);
660  }
661  }
662  }
oomph::FiniteElement::nplot_points_paraview
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.

Referenced by oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::scalar_value_paraview().

◆ strain_rate()

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

727  {
728 #ifdef PARANOID
729  if ((strainrate.ncol() != 3) || (strainrate.nrow() != 3))
730  {
731  std::ostringstream error_message;
732  error_message << "The strain rate has incorrect dimensions "
733  << strainrate.ncol() << " x " << strainrate.nrow()
734  << " Not 3" << std::endl;
735 
736  throw OomphLibError(
737  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
738  }
739 #endif
740 
741  // Find out how many nodes there are in the element
742  unsigned n_node = nnode();
743 
744  // Set up memory for the shape and test functions
745  Shape psi(n_node);
746  DShape dpsidx(n_node, 2);
747 
748  // Call the derivatives of the shape functions
749  dshape_eulerian(s, psi, dpsidx);
750 
751  // Radius
752  double interpolated_r = 0.0;
753 
754  // Velocity components and their derivatives
755  double ur = 0.0;
756  double durdr = 0.0;
757  double durdz = 0.0;
758  double uz = 0.0;
759  double duzdr = 0.0;
760  double duzdz = 0.0;
761  double uphi = 0.0;
762  double duphidr = 0.0;
763  double duphidz = 0.0;
764 
765  // Get the local storage for the indices
766  unsigned u_nodal_index[3];
767  for (unsigned i = 0; i < 3; ++i)
768  {
769  u_nodal_index[i] = u_index_axi_nst(i);
770  }
771 
772  // Loop over nodes to assemble velocities and their derivatives
773  // w.r.t. to r and z (x_0 and x_1)
774  for (unsigned l = 0; l < n_node; l++)
775  {
776  interpolated_r += nodal_position(l, 0) * psi[l];
777 
778  ur += nodal_value(l, u_nodal_index[0]) * psi[l];
779  uz += nodal_value(l, u_nodal_index[1]) * psi[l];
780  uphi += nodal_value(l, u_nodal_index[2]) * psi[l];
781 
782  durdr += nodal_value(l, u_nodal_index[0]) * dpsidx(l, 0);
783  durdz += nodal_value(l, u_nodal_index[0]) * dpsidx(l, 1);
784 
785  duzdr += nodal_value(l, u_nodal_index[1]) * dpsidx(l, 0);
786  duzdz += nodal_value(l, u_nodal_index[1]) * dpsidx(l, 1);
787 
788  duphidr += nodal_value(l, u_nodal_index[2]) * dpsidx(l, 0);
789  duphidz += nodal_value(l, u_nodal_index[2]) * dpsidx(l, 1);
790  }
791 
792 
793  // Assign strain rates without negative powers of the radius
794  // and zero those with:
795  strainrate(0, 0) = durdr;
796  strainrate(0, 1) = 0.5 * (durdz + duzdr);
797  strainrate(1, 0) = strainrate(0, 1);
798  strainrate(0, 2) = 0.0;
799  strainrate(2, 0) = strainrate(0, 2);
800  strainrate(1, 1) = duzdz;
801  strainrate(1, 2) = 0.5 * duphidz;
802  strainrate(2, 1) = strainrate(1, 2);
803  strainrate(2, 2) = 0.0;
804 
805 
806  // Overwrite the strain rates with negative powers of the radius
807  // unless we're at the origin
808  if (std::fabs(interpolated_r) > 1.0e-16)
809  {
810  double inverse_radius = 1.0 / interpolated_r;
811  strainrate(0, 2) = 0.5 * (duphidr - inverse_radius * uphi);
812  strainrate(2, 0) = strainrate(0, 2);
813  strainrate(2, 2) = inverse_radius * ur;
814  }
815  }
oomph::DenseMatrix::nrow
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:485
oomph::DenseMatrix::ncol
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:491
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117

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(), oomph::RefineableAxisymmetricNavierStokesEquations::get_Z2_flux(), oomph::AxisymmetricTCrouzeixRaviartElement::get_Z2_flux(), oomph::AxisymmetricTTaylorHoodElement::get_Z2_flux(), and traction().

◆ traction()

void oomph::AxisymmetricNavierStokesEquations::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

653  {
654  // throw OomphLibError(
655  // "Check the traction calculation for axisymmetric NSt",
656  // OOMPH_CURRENT_FUNCTION,
657  // OOMPH_EXCEPTION_LOCATION);
658 
659  // Pad out normal vector if required
660  Vector<double> n_local(3, 0.0);
661  n_local[0] = N[0];
662  n_local[1] = N[1];
663 
664 #ifdef PARANOID
665  if ((N.size() == 3) && (N[2] != 0.0))
666  {
667  throw OomphLibError(
668  "Unit normal passed into this fct should either be 2D (r,z) or have a "
669  "zero component in the theta-direction",
670  OOMPH_CURRENT_FUNCTION,
671  OOMPH_EXCEPTION_LOCATION);
672  }
673 #endif
674 
675  // Get velocity gradients
676  DenseMatrix<double> strainrate(3, 3);
677  strain_rate(s, strainrate);
678 
679  // Get pressure
680  double press = interpolated_p_axi_nst(s);
681 
682  // Loop over traction components
683  for (unsigned i = 0; i < 3; i++)
684  {
685  traction[i] = -press * n_local[i];
686  for (unsigned j = 0; j < 3; j++)
687  {
688  traction[i] += 2.0 * strainrate(i, j) * n_local[j];
689  }
690  }
691  }
oomph::AxisymmetricNavierStokesEquations::traction
void traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Definition: axisym_navier_stokes_elements.cc:649
N
@ N
Definition: constructor.cpp:22

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

Referenced by oomph::AxisymmetricQCrouzeixRaviartElement::get_traction(), and oomph::AxisymmetricQTaylorHoodElement::get_traction().

◆ u_index_axi_nst()

virtual unsigned oomph::AxisymmetricNavierStokesEquations::u_index_axi_nst ( const unsigned &  i) const
inlinevirtual

Return the index at which the i-th unknown velocity component is stored. The default value, i, is appropriate for single-physics problems. In derived multi-physics elements, this function should be overloaded to reflect the chosen storage scheme. Note that these equations require that the unknowns are always stored at the same indices at each node.

507  {
508  return i;
509  }

References i.

Referenced by dinterpolated_u_axi_nst_ddata(), oomph::RefineableAxisymmetricNavierStokesEquations::dinterpolated_u_axi_nst_ddata(), du_dt_axi_nst(), fill_in_contribution_to_hessian_vector_products(), fill_in_generic_dresidual_contribution_axi_nst(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::fill_in_off_diagonal_jacobian_blocks_analytic(), get_dresidual_dnodal_coordinates(), oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates(), oomph::RefineableAxisymmetricQTaylorHoodElement::get_interpolated_values(), oomph::RefineableAxisymmetricQCrouzeixRaviartElement::get_interpolated_values(), get_pressure_and_velocity_mass_matrix_diagonal(), RefineableQAxisymAdvectionDiffusionBoussinesqElement::identify_all_field_data_for_external_interaction(), oomph::AxisymmetricTCrouzeixRaviartElement::identify_load_data(), interpolated_d_dudx_dX_axi_nst(), interpolated_duds_axi_nst(), interpolated_dudx_axi_nst(), interpolated_u_axi_nst(), output_veloc(), strain_rate(), and u_index_nst().

◆ u_index_nst()

unsigned oomph::AxisymmetricNavierStokesEquations::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.

517  {
518  return this->u_index_axi_nst(i);
519  }

References u_index_axi_nst().

Referenced by oomph::AxisymmetricQAdvectionCrouzeixRaviartElement::fill_in_off_diagonal_jacobian_blocks_by_fd(), and oomph::FSIAxisymmetricQTaylorHoodElement::identify_load_data().

◆ viscosity_ratio()

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

Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of the Reynolds number

473  {
474  return *Viscosity_Ratio_pt;
475  }

References Viscosity_Ratio_pt.

Referenced by oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ viscosity_ratio_pt()

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

Pointer to Viscosity Ratio.

479  {
480  return Viscosity_Ratio_pt;
481  }

References Viscosity_Ratio_pt.

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

Member Data Documentation

◆ ALE_is_disabled

bool oomph::AxisymmetricNavierStokesEquations::ALE_is_disabled
protected

Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed. Only set to true if you're sure that the mesh is stationary

Referenced by disable_ALE(), enable_ALE(), fill_in_contribution_to_hessian_vector_products(), fill_in_generic_dresidual_contribution_axi_nst(), fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::further_build(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ axi_nst_body_force_fct_pt

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

Access function for the body-force pointer.

487  {
488  return Body_force_fct_pt;
489  }

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

◆ Body_force_fct_pt

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

Pointer to body force function.

Referenced by oomph::RefineableAxisymmetricNavierStokesEquations::further_build(), and get_body_force_axi_nst().

◆ Default_Gravity_vector

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

Static default value for the gravity vector.

Navier-Stokes equations default gravity vector.

Referenced by AxisymmetricNavierStokesEquations().

◆ Default_Physical_Constant_Value

double oomph::AxisymmetricNavierStokesEquations::Default_Physical_Constant_Value
staticprivate
Initial value:
=
0.0

Navier–Stokes equations static data.

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

Referenced by AxisymmetricNavierStokesEquations().

◆ Default_Physical_Ratio_Value

double oomph::AxisymmetricNavierStokesEquations::Default_Physical_Ratio_Value = 1.0
staticprivate

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

Referenced by AxisymmetricNavierStokesEquations().

◆ Density_Ratio_pt

double* oomph::AxisymmetricNavierStokesEquations::Density_Ratio_pt
protected

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

Referenced by AxisymmetricNavierStokesEquations(), density_ratio(), density_ratio_pt(), and oomph::RefineableAxisymmetricNavierStokesEquations::further_build().

◆ G_pt

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

Pointer to global gravity Vector.

Referenced by AxisymmetricNavierStokesEquations(), oomph::RefineableAxisymmetricNavierStokesEquations::further_build(), g(), and g_pt().

◆ Gamma

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

Vector to decide whether the stress-divergence form is used or not.

Navier–Stokes equations static data.

Referenced by fill_in_generic_residual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_axi_nst(), get_dresidual_dnodal_coordinates(), and oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates().

◆ Pressure_not_stored_at_node

int oomph::AxisymmetricNavierStokesEquations::Pressure_not_stored_at_node = -100
staticprivate

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

"Magic" negative number that indicates that the pressure is not stored at a node. This cannot be -1 because that represents the positional hanging scheme in the hanging_pt object of nodes

Referenced by p_nodal_index_axi_nst().

◆ Re_pt

double* oomph::AxisymmetricNavierStokesEquations::Re_pt
protected

Pointer to global Reynolds number.

Referenced by AxisymmetricNavierStokesEquations(), fill_in_generic_dresidual_contribution_axi_nst(), oomph::RefineableAxisymmetricNavierStokesEquations::further_build(), re(), and re_pt().

◆ ReInvFr_pt

double* oomph::AxisymmetricNavierStokesEquations::ReInvFr_pt
protected

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

Referenced by AxisymmetricNavierStokesEquations(), oomph::RefineableAxisymmetricNavierStokesEquations::further_build(), re_invfr(), and re_invfr_pt().

◆ ReInvRo_pt

double* oomph::AxisymmetricNavierStokesEquations::ReInvRo_pt
protected

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

Referenced by AxisymmetricNavierStokesEquations(), oomph::RefineableAxisymmetricNavierStokesEquations::further_build(), re_invro(), and re_invro_pt().

◆ ReSt_pt

double* oomph::AxisymmetricNavierStokesEquations::ReSt_pt
protected

Pointer to global Reynolds number x Strouhal number (=Womersley)

Referenced by AxisymmetricNavierStokesEquations(), oomph::RefineableAxisymmetricNavierStokesEquations::further_build(), re_st(), and re_st_pt().

◆ source_fct_pt

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

Access function for the source-function pointer.

493  {
494  return Source_fct_pt;
495  }

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

◆ Source_fct_pt

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

Pointer to volumetric source function.

Referenced by oomph::RefineableAxisymmetricNavierStokesEquations::further_build(), and get_source_fct().

◆ Viscosity_Ratio_pt

double* oomph::AxisymmetricNavierStokesEquations::Viscosity_Ratio_pt
protected

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

Referenced by AxisymmetricNavierStokesEquations(), oomph::RefineableAxisymmetricNavierStokesEquations::further_build(), get_viscosity_ratio_axisym_nst(), viscosity_ratio(), and viscosity_ratio_pt().

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