#include <spherical_navier_stokes_elements.h>

Inheritance diagram for oomph::SphericalNavierStokesEquations:

Public Types
typedef void(*	SphericalNavierStokesBodyForceFctPt) (const double &time, const Vector< double > &x, Vector< double > &body_force)

typedef double(*	SphericalNavierStokesSourceFctPt) (const double &time, const Vector< double > &x)

Public Types inherited from oomph::FiniteElement
typedef void(*	SteadyExactSolutionFctPt) (const Vector< double > &, Vector< double > &)

typedef void(*	UnsteadyExactSolutionFctPt) (const double &, const Vector< double > &, Vector< double > &)

Public Member Functions
double	cot (const double &th) const

	SphericalNavierStokesEquations ()

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 &	viscosity_ratio () const

double *&	viscosity_ratio_pt ()
	Pointer to Viscosity Ratio. More...

const double &	density_ratio () const

double *&	density_ratio_pt ()
	Pointer to Density ratio. More...

const double &	re_invfr () const
	Global inverse Froude number. More...

double *&	re_invfr_pt ()
	Pointer to global inverse Froude number. More...

const 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...

SphericalNavierStokesBodyForceFctPt &	body_force_fct_pt ()
	Access function for the body-force pointer. More...

SphericalNavierStokesBodyForceFctPt	body_force_fct_pt () const
	Access function for the body-force pointer. Const version. More...

SphericalNavierStokesSourceFctPt &	source_fct_pt ()
	Access function for the source-function pointer. More...

SphericalNavierStokesSourceFctPt	source_fct_pt () const
	Access function for the source-function pointer. Const version. More...

virtual unsigned	npres_spherical_nst () const =0
	Function to return number of pressure degrees of freedom. More...

double	u_spherical_nst (const unsigned &n, const unsigned &i) const

double	u_spherical_nst (const unsigned &t, const unsigned &n, const unsigned &i) const

virtual unsigned	u_index_spherical_nst (const unsigned &i) const

double	du_dt_spherical_nst (const unsigned &n, const unsigned &i) const

void	disable_ALE ()

void	enable_ALE ()

virtual double	p_spherical_nst (const unsigned &n_p) const =0

virtual void	fix_pressure (const unsigned &p_dof, const double &p_value)=0
	Pin p_dof-th pressure dof and set it to value specified by p_value. More...

virtual int	p_nodal_index_spherical_nst () const

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

void	get_vorticity (const Vector< double > &s, Vector< double > &vorticity) const
	Compute the vorticity vector at local coordinate s. More...

double	kin_energy () const
	Get integral of kinetic energy over element. More...

double	d_kin_energy_dt () const
	Get integral of time derivative of kinetic energy over element. More...

void	strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
	Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i) More...

void	get_traction (const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)

void	get_load (const Vector< double > &s, const Vector< double > &N, Vector< double > &load)

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

void	output (std::ostream &outfile)

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

void	output (FILE *file_pt)

void	output (FILE *file_pt, const unsigned &nplot)

void	full_output (std::ostream &outfile)

void	full_output (std::ostream &outfile, const unsigned &nplot)

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

void	output_vorticity (std::ostream &outfile, const unsigned &nplot)

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	compute_error_e (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, FiniteElement::SteadyExactSolutionFctPt exact_soln_dr_pt, FiniteElement::SteadyExactSolutionFctPt exact_soln_dtheta_pt, double &u_error, double &u_norm, double &p_error, double &p_norm)

void	compute_shear_stress (std::ostream &outfile)

void	extract_velocity (std::ostream &outfile)

Vector< double >	actual (const Vector< double > &x)

Vector< double >	actual_dr (const Vector< double > &x)

Vector< double >	actual_dth (const Vector< double > &x)

void	fill_in_contribution_to_residuals (Vector< double > &residuals)
	Compute the element's residual Vector. More...

void	fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)

void	fill_in_contribution_to_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)

void	interpolated_u_spherical_nst (const Vector< double > &s, Vector< double > &veloc) const
	Compute vector of FE interpolated velocity u at local coordinate s. More...

double	interpolated_u_spherical_nst (const Vector< double > &s, const unsigned &i) const
	Return FE interpolated velocity u[i] at local coordinate s. More...

double	interpolated_dudx_spherical_nst (const Vector< double > &s, const unsigned &i, const unsigned &j) const

double	interpolated_p_spherical_nst (const Vector< double > &s) const
	Return FE interpolated pressure at local coordinate s. More...

Public Member Functions inherited from oomph::FSIFluidElement
	FSIFluidElement ()
	Constructor. More...

	FSIFluidElement (const FSIFluidElement &)=delete
	Broken copy constructor. More...

void	operator= (const FSIFluidElement &)=delete
	Broken assignment operator. More...

virtual void	identify_load_data (std::set< std::pair< Data *, unsigned >> &paired_load_data)=0

virtual void	identify_pressure_data (std::set< std::pair< Data *, unsigned >> &paired_pressure_data)=0

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 void	get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)

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

virtual double	compute_physical_size () const

virtual void	point_output_data (const Vector< double > &s, Vector< double > &data)

void	point_output (std::ostream &outfile, const Vector< double > &s)

virtual unsigned	nplot_points_paraview (const unsigned &nplot) const

virtual unsigned	nsub_elements_paraview (const unsigned &nplot) const

void	output_paraview (std::ofstream &file_out, const unsigned &nplot) const

virtual void	write_paraview_output_offset_information (std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const

virtual void	write_paraview_type (std::ofstream &file_out, const unsigned &nplot) const

virtual void	write_paraview_offsets (std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const

virtual unsigned	nscalar_paraview () const

virtual void	scalar_value_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const

virtual void	scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const

virtual void	scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const

virtual std::string	scalar_name_paraview (const unsigned &i) const

virtual void	output (const unsigned &t, std::ostream &outfile, const unsigned &n_plot) const

virtual void	output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, const SolutionFunctorBase &exact_soln) const
	Output a time-dependent exact solution over the element. More...

virtual void	get_s_plot (const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const

virtual std::string	tecplot_zone_string (const unsigned &nplot) const

virtual void	write_tecplot_zone_footer (std::ostream &outfile, const unsigned &nplot) const

virtual void	write_tecplot_zone_footer (FILE *file_pt, const unsigned &nplot) const

virtual unsigned	nplot_points (const unsigned &nplot) const

virtual void	compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
	Calculate the norm of the error and that of the exact solution. More...

virtual void	compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
	Calculate the norm of the error and that of the exact solution. More...

virtual void	compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)

virtual void	compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)

virtual void	compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)

virtual void	compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)

virtual void	compute_abs_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error)

void	integrate_fct (FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt, Vector< double > &integral)
	Integrate Vector-valued function over element. More...

void	integrate_fct (FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt, const double &time, Vector< double > &integral)
	Integrate Vector-valued time-dep function over element. More...

virtual void	build_face_element (const int &face_index, FaceElement *face_element_pt)

virtual unsigned	self_test ()

virtual unsigned	get_bulk_node_number (const int &face_index, const unsigned &i) const

virtual int	face_outer_unit_normal_sign (const int &face_index) const
	Get the sign of the outer unit normal on the face given by face_index. More...

virtual unsigned	nnode_on_face () const

void	face_node_number_error_check (const unsigned &i) const
	Range check for face node numbers. More...

virtual CoordinateMappingFctPt	face_to_bulk_coordinate_fct_pt (const int &face_index) const

virtual BulkCoordinateDerivativesFctPt	bulk_coordinate_derivatives_fct_pt (const int &face_index) const

Public Member Functions inherited from oomph::GeneralisedElement
	GeneralisedElement ()
	Constructor: Initialise all pointers and all values to zero. More...

virtual	~GeneralisedElement ()
	Virtual destructor to clean up any memory allocated by the object. More...

	GeneralisedElement (const GeneralisedElement &)=delete
	Broken copy constructor. More...

void	operator= (const GeneralisedElement &)=delete
	Broken assignment operator. More...

Data *&	internal_data_pt (const unsigned &i)
	Return a pointer to i-th internal data object. More...

Data *const &	internal_data_pt (const unsigned &i) const
	Return a pointer to i-th internal data object (const version) More...

Data *&	external_data_pt (const unsigned &i)
	Return a pointer to i-th external data object. More...

Data *const &	external_data_pt (const unsigned &i) const
	Return a pointer to i-th external data object (const version) More...

unsigned long	eqn_number (const unsigned &ieqn_local) const

int	local_eqn_number (const unsigned long &ieqn_global) const

unsigned	add_external_data (Data *const &data_pt, const bool &fd=true)

bool	external_data_fd (const unsigned &i) const

void	exclude_external_data_fd (const unsigned &i)

void	include_external_data_fd (const unsigned &i)

void	flush_external_data ()
	Flush all external data. More...

void	flush_external_data (Data *const &data_pt)
	Flush the object addressed by data_pt from the external data array. More...

unsigned	ninternal_data () const
	Return the number of internal data objects. More...

unsigned	nexternal_data () const
	Return the number of external data objects. More...

unsigned	ndof () const
	Return the number of equations/dofs in the element. More...

void	dof_vector (const unsigned &t, Vector< double > &dof)
	Return the vector of dof values at time level t. More...

void	dof_pt_vector (Vector< double * > &dof_pt)
	Return the vector of pointers to dof values. More...

void	set_internal_data_time_stepper (const unsigned &i, TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)

void	assign_internal_eqn_numbers (unsigned long &global_number, Vector< double * > &Dof_pt)

void	describe_dofs (std::ostream &out, const std::string &current_string) const

void	add_internal_value_pt_to_map (std::map< unsigned, double * > &map_of_value_pt)

virtual void	assign_local_eqn_numbers (const bool &store_local_dof_pt)

virtual void	complete_setup_of_dependencies ()

virtual void	get_residuals (Vector< double > &residuals)

virtual void	get_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)

virtual void	get_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)

virtual void	get_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)

virtual void	get_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)

virtual void	get_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)

virtual void	get_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)

virtual void	get_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)

virtual void	get_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)

virtual void	get_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)

virtual void	compute_norm (Vector< double > &norm)

virtual void	compute_norm (double &norm)

virtual unsigned	ndof_types () const

virtual void	get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const

Public Member Functions inherited from oomph::GeomObject
	GeomObject ()
	Default constructor. More...

	GeomObject (const unsigned &ndim)

	GeomObject (const unsigned &nlagrangian, const unsigned &ndim)

	GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)

	GeomObject (const GeomObject &dummy)=delete
	Broken copy constructor. More...

void	operator= (const GeomObject &)=delete
	Broken assignment operator. More...

virtual	~GeomObject ()
	(Empty) destructor More...

unsigned	nlagrangian () const
	Access function to # of Lagrangian coordinates. More...

unsigned	ndim () const
	Access function to # of Eulerian coordinates. More...

void	set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
	Set # of Lagrangian and Eulerian coordinates. More...

TimeStepper *&	time_stepper_pt ()

TimeStepper *	time_stepper_pt () const

virtual void	position (const double &t, const Vector< double > &zeta, Vector< double > &r) const

virtual void	dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const

virtual void	d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const

virtual void	d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const

Public Member Functions inherited from oomph::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...

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_spherical_nst (const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0

virtual double	dshape_and_dtest_eulerian_at_knot_spherical_nst (const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0

virtual void	pshape_spherical_nst (const Vector< double > &s, Shape &psi) const =0
	Compute the pressure shape functions at local coordinate s. More...

virtual void	pshape_spherical_nst (const Vector< double > &s, Shape &psi, Shape &test) const =0

virtual void	get_body_force_spherical_nst (const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)

virtual double	get_source_spherical_nst (double time, const unsigned &ipt, const Vector< double > &x)

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

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)

template<>
double	invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
	Zero-d specialisation of function to calculate inverse of jacobian mapping. More...

template<>
double	invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
	One-d specialisation of function to calculate inverse of jacobian mapping. More...

template<>
double	invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
	Two-d specialisation of function to calculate inverse of jacobian mapping. More...

template<>
double	invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const

template<>
void	dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const

template<>
void	dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const

template<>
void	dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const

template<>
void	dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const

template<>
void	d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const

template<>
void	d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const

template<>
void	d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const

template<>
void	d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const

template<>
void	transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const

template<>
void	transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const

template<>
void	transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const

template<>
void	transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const

Protected Member Functions inherited from oomph::GeneralisedElement
unsigned	add_internal_data (Data *const &data_pt, const bool &fd=true)

bool	internal_data_fd (const unsigned &i) const

void	exclude_internal_data_fd (const unsigned &i)

void	include_internal_data_fd (const unsigned &i)

void	clear_global_eqn_numbers ()

void	add_global_eqn_numbers (std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)

virtual void	assign_internal_and_external_local_eqn_numbers (const bool &store_local_dof_pt)

virtual void	assign_additional_local_eqn_numbers ()

int	internal_local_eqn (const unsigned &i, const unsigned &j) const

int	external_local_eqn (const unsigned &i, const unsigned &j)

void	fill_in_jacobian_from_internal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)

void	fill_in_jacobian_from_internal_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)

void	fill_in_jacobian_from_external_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)

void	fill_in_jacobian_from_external_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)

virtual void	update_before_internal_fd ()

virtual void	reset_after_internal_fd ()

virtual void	update_in_internal_fd (const unsigned &i)

virtual void	reset_in_internal_fd (const unsigned &i)

virtual void	update_before_external_fd ()

virtual void	reset_after_external_fd ()

virtual void	update_in_external_fd (const unsigned &i)

virtual void	reset_in_external_fd (const unsigned &i)

virtual void	fill_in_contribution_to_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)

virtual void	fill_in_contribution_to_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)

virtual void	fill_in_contribution_to_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)

virtual void	fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)

virtual void	fill_in_contribution_to_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)

virtual void	fill_in_contribution_to_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)

virtual void	fill_in_contribution_to_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)

Protected 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...

SphericalNavierStokesBodyForceFctPt	Body_force_fct_pt
	Pointer to body force function. More...

SphericalNavierStokesSourceFctPt	Source_fct_pt
	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 = 0.0
	Navier–Stokes equations static data. More...

static double	Default_Physical_Ratio_Value = 1.0
	Navier–Stokes equations static data. More...

static Vector< double >	Default_Gravity_vector
	Static default value for the gravity vector. More...

Additional Inherited Members
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 Navier–Stokes equations, in axisymmetric spherical polar coordinates. This contains the generic maths – any concrete implementation must be derived from this.

We also provide all functions required to use this element in FSI problems, by deriving it from the FSIFluidElement base class.

Member Typedef Documentation

◆ SphericalNavierStokesBodyForceFctPt

typedef void(* oomph::SphericalNavierStokesEquations::SphericalNavierStokesBodyForceFctPt) (const double &time, const Vector< double > &x, Vector< double > &body_force)

Function pointer to body force function fct(t,x,f(x)) x is a Vector!

◆ SphericalNavierStokesSourceFctPt

typedef double(* oomph::SphericalNavierStokesEquations::SphericalNavierStokesSourceFctPt) (const double &time, const Vector< double > &x)

Function pointer to source function fct(t,x) x is a Vector!

Constructor & Destructor Documentation

◆ SphericalNavierStokesEquations()

oomph::SphericalNavierStokesEquations::SphericalNavierStokesEquations ( )

inline

Constructor: NULL the body force and source function and make sure the ALE terms are included by default.

       : Body_force_fct_pt(0), Source_fct_pt(0), ALE_is_disabled(false)
     {
       // Set all the Physical parameter pointers to the default value zero
       Re_pt = &Default_Physical_Constant_Value;
       ReSt_pt = &Default_Physical_Constant_Value;
       ReInvFr_pt = &Default_Physical_Constant_Value;
       ReInvRo_pt = &Default_Physical_Constant_Value;
       G_pt = &Default_Gravity_vector;
       // Set the Physical ratios to the default value of 1
       Viscosity_Ratio_pt = &Default_Physical_Ratio_Value;
       Density_Ratio_pt = &Default_Physical_Ratio_Value;
     }

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

◆ actual()

Vector<double> oomph::SphericalNavierStokesEquations::actual ( const Vector< double > & x )

inline

     {
       const double Re = this->re();
  
       double r = x[0];
       double theta = x[1];
       Vector<double> ans(4, 0.0);
  
       ans[2] = r * sin(theta);
       ans[3] = 0.5 * Re * r * r * sin(theta) * sin(theta);
  
       return (ans);
     }

References UniformPSDSelfTest::r, GlobalPhysicalVariables::Re, re(), sin(), BiharmonicTestFunctions2::theta, and plotDoE::x.

◆ actual_dr()

Vector<double> oomph::SphericalNavierStokesEquations::actual_dr ( const Vector< double > & x )

inline

     {
       const double Re = this->re();
  
       double r = x[0];
       double theta = x[1];
       Vector<double> ans(4, 0.0);
  
       ans[2] = sin(theta);
       ans[3] = Re * r * sin(theta) * sin(theta);
  
       return (ans);
     }

References UniformPSDSelfTest::r, GlobalPhysicalVariables::Re, re(), sin(), BiharmonicTestFunctions2::theta, and plotDoE::x.

◆ actual_dth()

Vector<double> oomph::SphericalNavierStokesEquations::actual_dth ( const Vector< double > & x )

inline

     {
       const double Re = this->re();
  
       double r = x[0];
       double theta = x[1];
       Vector<double> ans(4, 0.0);
  
       ans[2] = r * cos(theta);
       ans[3] = Re * r * r * sin(theta) * cos(theta);
  
       return (ans);
     }

References cos(), UniformPSDSelfTest::r, GlobalPhysicalVariables::Re, re(), sin(), BiharmonicTestFunctions2::theta, and plotDoE::x.

◆ body_force_fct_pt() [1/2]

SphericalNavierStokesBodyForceFctPt& oomph::SphericalNavierStokesEquations::body_force_fct_pt ( )

inline

Access function for the body-force pointer.

     {
       return Body_force_fct_pt;
     }

References Body_force_fct_pt.

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

◆ body_force_fct_pt() [2/2]

SphericalNavierStokesBodyForceFctPt oomph::SphericalNavierStokesEquations::body_force_fct_pt ( ) const

inline

Access function for the body-force pointer. Const version.

     {
       return Body_force_fct_pt;
     }

References Body_force_fct_pt.

◆ compute_error() [1/2]

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

   {
     error = 0.0;
     norm = 0.0;
  
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Vector for coordintes
     Vector<double> x(2);
  
     // Set the value of n_intpt
     unsigned n_intpt = integral_pt()->nweight();
  
  
     outfile << "ZONE" << std::endl;
  
     // Exact solution Vector (u,v,w,p)
     Vector<double> exact_soln(4);
  
     // Loop over the integration points
     for (unsigned ipt = 0; ipt < n_intpt; ipt++)
     {
       // Assign values of s
       for (unsigned i = 0; i < 2; i++)
       {
         s[i] = integral_pt()->knot(ipt, i);
       }
  
       // Get the integral weight
       double w = integral_pt()->weight(ipt);
  
       // Get jacobian of mapping
       double J = J_eulerian(s);
  
       // Premultiply the weights and the Jacobian
       double W = w * J;
  
       // Get x position as Vector
       interpolated_x(s, x);
  
       // Get the exact solution at this point
       (*exact_soln_pt)(x, exact_soln);
  
       // Velocity error
       for (unsigned i = 0; i < 3; i++)
       {
         norm += exact_soln[i] * exact_soln[i] * W;
         error += (exact_soln[i] - interpolated_u_spherical_nst(s, i)) *
                  (exact_soln[i] - interpolated_u_spherical_nst(s, i)) * W;
       }
  
       // Output x,y,...,u_exact
       for (unsigned i = 0; i < 2; i++)
       {
         outfile << x[i] << " ";
       }
  
       // Output x,y,[z],u_error,v_error,[w_error]
       for (unsigned i = 0; i < 3; i++)
       {
         outfile << exact_soln[i] - interpolated_u_spherical_nst(s, i) << " ";
       }
       outfile << std::endl;
     }
   }

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

◆ compute_error() [2/2]

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

   {
     error = 0.0;
     norm = 0.0;
  
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Vector for coordintes
     Vector<double> x(2);
  
     // Set the value of n_intpt
     unsigned n_intpt = integral_pt()->nweight();
  
     outfile << "ZONE" << std::endl;
  
     // Exact solution Vector (u,v,[w],p)
     Vector<double> exact_soln(4);
  
     // Loop over the integration points
     for (unsigned ipt = 0; ipt < n_intpt; ipt++)
     {
       // Assign values of s
       for (unsigned i = 0; i < 2; i++)
       {
         s[i] = integral_pt()->knot(ipt, i);
       }
  
       // Get the integral weight
       double w = integral_pt()->weight(ipt);
  
       // Get jacobian of mapping
       double J = J_eulerian(s);
  
       // Premultiply the weights and the Jacobian
       double W = w * J;
  
       // Get x position as Vector
       interpolated_x(s, x);
  
       // Get exact solution at this point
       (*exact_soln_pt)(time, x, exact_soln);
  
       // Velocity error
       for (unsigned i = 0; i < 3; i++)
       {
         norm += exact_soln[i] * exact_soln[i] * W;
         error += (exact_soln[i] - interpolated_u_spherical_nst(s, i)) *
                  (exact_soln[i] - interpolated_u_spherical_nst(s, i)) * W;
       }
  
       // Output x,y,...,u_exact
       for (unsigned i = 0; i < 2; i++)
       {
         outfile << x[i] << " ";
       }
  
       // Output x,y,[z],u_error,v_error,[w_error]
       for (unsigned i = 0; i < 3; i++)
       {
         outfile << exact_soln[i] - interpolated_u_spherical_nst(s, i) << " ";
       }
  
       outfile << std::endl;
     }
   }

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

◆ compute_error_e()

void oomph::SphericalNavierStokesEquations::compute_error_e	(	std::ostream &	outfile,
		FiniteElement::SteadyExactSolutionFctPt	exact_soln_pt,
		FiniteElement::SteadyExactSolutionFctPt	exact_soln_dr_pt,
		FiniteElement::SteadyExactSolutionFctPt	exact_soln_dtheta_pt,
		double &	u_error,
		double &	u_norm,
		double &	p_error,
		double &	p_norm
	)

Validate against exact solution. Solution is provided direct from exact_soln function. Plot at a given number of plot points and compute the energy error and energy norm of the velocity solution over the element.

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

Exact solution Vectors (u,v,[w],p) -

need two extras to deal with the derivatives in the energy norm

   {
     u_error = 0.0;
     p_error = 0.0;
     u_norm = 0.0;
     p_norm = 0.0;
  
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Vector for coordinates
     Vector<double> x(2);
  
     // Set the value of n_intpt
     unsigned n_intpt = integral_pt()->nweight();
  
  
     Vector<double> exact_soln(4);
     Vector<double> exact_soln_dr(4);
     Vector<double> exact_soln_dth(4);
  
  
     // Loop over the integration points
     for (unsigned ipt = 0; ipt < n_intpt; ipt++)
     {
       // Assign values of s
       for (unsigned i = 0; i < 2; i++)
       {
         s[i] = integral_pt()->knot(ipt, i);
       }
  
       // Get the integral weight
       double w = integral_pt()->weight(ipt);
  
       // Get jacobian of mapping
       double J = J_eulerian(s);
  
       // Premultiply the weights and the Jacobian
       double W = w * J;
  
       // Get x position as Vector
       interpolated_x(s, x);
  
       // Get exact solution at the point x using the 'actual' functions
       (*exact_soln_pt)(x, exact_soln);
       (*exact_soln_dr_pt)(x, exact_soln_dr);
       (*exact_soln_dtheta_pt)(x, exact_soln_dth);
  
       // exact_soln = actual(x);
       // exact_soln_dr = actual_dr(x);
       // exact_soln_dth = actual_dth(x);
  
  
       double r = x[0];
       double theta = x[1];
  
       W *= r * r * sin(theta);
  
  
       // Velocity error in energy norm
       // Owing to the additional terms arising from the phi-components, we can't
       // loop over the velocities. Calculate each section separately instead.
  
       // Norm calculation involves the double contraction of two tensors, so we
       // take the summation of nine separate components for clarity. Fortunately
       // most of the components zero out in the axisymmetric case.
       //---------------------------------------------
       // Entry (1,1) in dyad matrix
       // No contribution
       // Entry (1,2) in dyad matrix
       // No contribution
       // Entry (1,3) in dyad matrix
       u_norm += (1 / (r * r)) * exact_soln[2] * exact_soln[2] * W;
       //----------------------------------------------
       // Entry (2,1) in dyad matrix
       // No contribution
       // Entry (2,2) in dyad matrix
       // No contribution
       // Entry (2,3) in dyad matrix
       u_norm += cot(theta) * cot(theta) * (1 / (r * r)) * exact_soln[2] *
                 exact_soln[2] * W;
       //-----------------------------------------------
       // Entry (3,1) in dyad matrix
       u_norm += exact_soln_dr[2] * exact_soln_dr[2] * W;
       // Entry (3,2) in dyad matrix
       u_norm += (1 / (r * r)) * exact_soln_dth[2] * exact_soln_dth[2] * W;
       // Entry (3,3) in dyad matrix
       // No contribution
       //-----------------------------------------------
       // End of norm calculation
  
  
       // Error calculation involves the summation of 9 squared components,
       // stated separately here for clarity
       //---------------------------------------------
       // Entry (1,1) in dyad matrix
       u_error += interpolated_dudx_spherical_nst(s, 0, 0) *
                  interpolated_dudx_spherical_nst(s, 0, 0) * W;
       // Entry (1,2) in dyad matrix
       u_error += (1 / (r * r)) *
                  (interpolated_u_spherical_nst(s, 1) -
                   interpolated_dudx_spherical_nst(s, 0, 1)) *
                  (interpolated_u_spherical_nst(s, 1) -
                   interpolated_dudx_spherical_nst(s, 0, 1)) *
                  W;
       // Entry (1,3) in dyad matrix
       u_error += (1 / (r * r)) *
                  (interpolated_u_spherical_nst(s, 2) - exact_soln[2]) *
                  (interpolated_u_spherical_nst(s, 2) - exact_soln[2]) * W;
       //---------------------------------------------
       // Entry (2,1) in dyad matrix
       u_error += interpolated_dudx_spherical_nst(s, 1, 0) *
                  interpolated_dudx_spherical_nst(s, 1, 0) * W;
       // Entry (2,2) in dyad matrix
       u_error += (1 / (r * r)) *
                  (interpolated_dudx_spherical_nst(s, 1, 1) +
                   interpolated_u_spherical_nst(s, 0)) *
                  (interpolated_dudx_spherical_nst(s, 1, 1) +
                   interpolated_u_spherical_nst(s, 0)) *
                  W;
       // Entry (2,3) in dyad matrix
       u_error += (1 / (r * r)) * cot(theta) * cot(theta) *
                  (interpolated_u_spherical_nst(s, 2) - exact_soln[2]) *
                  (interpolated_u_spherical_nst(s, 2) - exact_soln[2]) * W;
       //---------------------------------------------
       // Entry (3,1) in dyad matrix
       u_error += (exact_soln_dr[2] - interpolated_dudx_spherical_nst(s, 2, 0)) *
                  (exact_soln_dr[2] - interpolated_dudx_spherical_nst(s, 2, 0)) *
                  W;
       // Entry (3,2) in dyad matrix
       u_error +=
         (1 / (r * r)) *
         (exact_soln_dth[2] - interpolated_dudx_spherical_nst(s, 2, 1)) *
         (exact_soln_dth[2] - interpolated_dudx_spherical_nst(s, 2, 1)) * W;
       // Entry (3,3) in dyad matrix
       u_error += (1 / (r * r)) *
                  (interpolated_u_spherical_nst(s, 0) -
                   cot(theta) * interpolated_u_spherical_nst(s, 1)) *
                  (interpolated_u_spherical_nst(s, 0) -
                   cot(theta) * interpolated_u_spherical_nst(s, 1)) *
                  W;
       //---------------------------------------------
       // End of velocity error calculation
  
       // Pressure error in L2 norm - energy norm not required here as the
       // pressure only appears undifferentiated in the NS SPC weak form.
  
       p_norm += exact_soln[3] * exact_soln[3] * W;
  
       p_error += (exact_soln[3] - interpolated_p_spherical_nst(s)) *
                  (exact_soln[3] - interpolated_p_spherical_nst(s)) * W;
  
  
       // Output r and theta coordinates
       for (unsigned i = 0; i < 2; i++)
       {
         outfile << x[i] << " ";
       }
  
       // Output the velocity error details in the energy norm
       outfile << u_error << "  " << u_norm << "  ";
       // Output the pressure error details in the L2 norm
       outfile << p_error << "  " << p_norm << "  ";
       outfile << std::endl;
  
     } // End loop over integration points
  
   } // End of function statement

References cot(), ProblemParameters::exact_soln(), i, oomph::FiniteElement::integral_pt(), interpolated_dudx_spherical_nst(), interpolated_p_spherical_nst(), interpolated_u_spherical_nst(), oomph::FiniteElement::interpolated_x(), J, oomph::FiniteElement::J_eulerian(), oomph::Integral::knot(), oomph::Integral::nweight(), UniformPSDSelfTest::r, s, sin(), BiharmonicTestFunctions2::theta, w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and plotDoE::x.

◆ compute_shear_stress()

void oomph::SphericalNavierStokesEquations::compute_shear_stress ( std::ostream & outfile )

   {
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Vector for coordinates
     Vector<double> x(2);
  
     // Number of plot points
     unsigned Npts = 3;
  
     // Declaration of the shear stress variable
     double shear = 0.0;
  
  
     // Assign values of s
     for (unsigned i = 0; i < Npts; i++)
     {
       s[0] = 1.0;
       s[1] = -1 + 2.0 * i / (Npts - 1);
  
  
       // Get position as global coordinate
       interpolated_x(s, x);
  
       // Output r and theta coordinates
       outfile << x[0] << " ";
       outfile << x[1] << " ";
  
       // Output the shear stress at the current coordinate
       double r = x[0];
       shear = interpolated_dudx_spherical_nst(s, 2, 0) -
               (1 / r) * interpolated_u_spherical_nst(s, 2);
       outfile << shear << "  ";
  
       outfile << std::endl;
     }
   }

References i, interpolated_dudx_spherical_nst(), interpolated_u_spherical_nst(), oomph::FiniteElement::interpolated_x(), UniformPSDSelfTest::r, s, plotDoE::shear, and plotDoE::x.

◆ cot()

double oomph::SphericalNavierStokesEquations::cot ( const double & th ) const

inline

Include a cot function to simplify the NS momentum and jacobian expressions

     {
       return (1 / tan(th));
     }

References Eigen::bfloat16_impl::tan().

Referenced by compute_error_e(), and strain_rate().

◆ d_kin_energy_dt()

double oomph::SphericalNavierStokesEquations::d_kin_energy_dt ( ) const

Get integral of time derivative of kinetic energy over element.

Get integral of time derivative of kinetic energy over element:

   {
     throw OomphLibError("Not tested for spherical Navier-Stokes elements",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
  
  
     // Initialise
     double d_kin_en_dt = 0.0;
  
     // Set the value of n_intpt
     const unsigned n_intpt = integral_pt()->nweight();
  
     // Set the Vector to hold local coordinates
     Vector<double> s(2);
  
     // Get the number of nodes
     const unsigned n_node = this->nnode();
  
     // Storage for the shape function
     Shape psi(n_node);
     DShape dpsidx(n_node, 2);
  
     // Get the value at which the velocities are stored
     unsigned u_index[3];
     for (unsigned i = 0; i < 3; i++)
     {
       u_index[i] = this->u_index_spherical_nst(i);
     }
  
     // Loop over the integration points
     for (unsigned ipt = 0; ipt < n_intpt; ipt++)
     {
       // Get the jacobian of the mapping
       double J = dshape_eulerian_at_knot(ipt, psi, dpsidx);
  
       // Get the integral weight
       double w = integral_pt()->weight(ipt);
  
       // Now work out the velocity and the time derivative
       Vector<double> interpolated_u(3, 0.0), interpolated_dudt(3, 0.0);
  
       // Loop over the shape functions
       for (unsigned l = 0; l < n_node; l++)
       {
         // Loop over the velocity components
         for (unsigned i = 0; i < 3; i++)
         {
           interpolated_u[i] += nodal_value(l, u_index[i]) * psi(l);
           interpolated_dudt[i] += du_dt_spherical_nst(l, u_index[i]) * psi(l);
         }
       }
  
       // Get mesh velocity bit
       if (!ALE_is_disabled)
       {
         Vector<double> mesh_velocity(2, 0.0);
         DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
  
         // Loop over nodes again
         for (unsigned l = 0; l < n_node; l++)
         {
           // Mesh can only move in the first two directions
           for (unsigned i = 0; i < 2; i++)
           {
             mesh_velocity[i] += this->dnodal_position_dt(l, i) * psi(l);
           }
  
           // Now du/dx bit
           for (unsigned i = 0; i < 3; i++)
           {
             for (unsigned j = 0; j < 2; j++)
             {
               interpolated_dudx(i, j) +=
                 this->nodal_value(l, u_index[i]) * dpsidx(l, j);
             }
           }
         }
  
         // Subtract mesh velocity from du_dt
         for (unsigned i = 0; i < 3; i++)
         {
           for (unsigned k = 0; k < 2; k++)
           {
             interpolated_dudt[i] -= mesh_velocity[k] * interpolated_dudx(i, k);
           }
         }
       }
  
  
       // Loop over directions and add up u du/dt  terms
       double sum = 0.0;
       for (unsigned i = 0; i < 3; i++)
       {
         sum += interpolated_u[i] * interpolated_dudt[i];
       }
  
       d_kin_en_dt += sum * w * J;
     }
  
     return d_kin_en_dt;
   }

References ALE_is_disabled, oomph::FiniteElement::dnodal_position_dt(), oomph::FiniteElement::dshape_eulerian_at_knot(), du_dt_spherical_nst(), i, oomph::FiniteElement::integral_pt(), J, j, k, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, u_index_spherical_nst(), w, and oomph::Integral::weight().

◆ density_ratio()

const double& oomph::SphericalNavierStokesEquations::density_ratio ( ) const

inline

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

     {
       return *Density_Ratio_pt;
     }

References Density_Ratio_pt.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), and fill_in_generic_residual_contribution_spherical_nst().

◆ density_ratio_pt()

double*& oomph::SphericalNavierStokesEquations::density_ratio_pt ( )

inline

Pointer to Density ratio.

     {
       return Density_Ratio_pt;
     }

References Density_Ratio_pt.

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

◆ disable_ALE()

void oomph::SphericalNavierStokesEquations::disable_ALE ( )

inlinevirtual

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

Reimplemented from oomph::FiniteElement.

     {
       ALE_is_disabled = true;
     }

References ALE_is_disabled.

Referenced by oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::disable_ALE().

◆ dissipation() [1/2]

double oomph::SphericalNavierStokesEquations::dissipation ( ) const

Return integral of dissipation over element.

   {
     // Initialise
     double diss = 0.0;
  
     // Set the value of n_intpt
     unsigned n_intpt = integral_pt()->nweight();
  
     // Set the Vector to hold local coordinates
     Vector<double> s(2);
  
     // Loop over the integration points
     for (unsigned ipt = 0; ipt < n_intpt; ipt++)
     {
       // Assign values of s
       for (unsigned i = 0; i < 2; i++)
       {
         s[i] = integral_pt()->knot(ipt, i);
       }
  
       // Get the integral weight
       double w = integral_pt()->weight(ipt);
  
       // Get Jacobian of mapping
       double J = J_eulerian(s);
  
       // Get strain rate matrix
       DenseMatrix<double> strainrate(3, 3);
       strain_rate(s, strainrate);
  
       // Initialise
       double local_diss = 0.0;
       for (unsigned i = 0; i < 3; i++)
       {
         for (unsigned j = 0; j < 3; j++)
         {
           local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
         }
       }
  
       diss += local_diss * w * J;
     }
  
     return diss;
   }

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

Referenced by full_output(), and output().

◆ dissipation() [2/2]

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

Return dissipation at local coordinate s.

   {
     // Get strain rate matrix
     DenseMatrix<double> strainrate(3, 3);
     strain_rate(s, strainrate);
  
     // Initialise
     double local_diss = 0.0;
     for (unsigned i = 0; i < 3; i++)
     {
       for (unsigned j = 0; j < 3; j++)
       {
         local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
       }
     }
  
     return local_diss;
   }

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

◆ dshape_and_dtest_eulerian_at_knot_spherical_nst()

virtual double oomph::SphericalNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_spherical_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::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), and fill_in_generic_residual_contribution_spherical_nst().

◆ dshape_and_dtest_eulerian_spherical_nst()

virtual double oomph::SphericalNavierStokesEquations::dshape_and_dtest_eulerian_spherical_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::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

◆ du_dt_spherical_nst()

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

     {
       // Get the data's timestepper
       TimeStepper* time_stepper_pt = this->node_pt(n)->time_stepper_pt();
  
       // Initialise dudt
       double dudt = 0.0;
  
       // Loop over the timesteps, if there is a non Steady timestepper
       if (time_stepper_pt->type() != "Steady")
       {
         // Find the index at which the dof is stored
         const unsigned u_nodal_index = this->u_index_spherical_nst(i);
  
         // Number of timsteps (past & present)
         const unsigned n_time = time_stepper_pt->ntstorage();
         // Loop over the timesteps
         for (unsigned t = 0; t < n_time; t++)
         {
           dudt +=
             time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
         }
       }
  
       return dudt;
     }

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

Referenced by d_kin_energy_dt(), oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), fill_in_generic_residual_contribution_spherical_nst(), and full_output().

◆ enable_ALE()

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

     {
       ALE_is_disabled = false;
     }

References ALE_is_disabled.

Referenced by oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::enable_ALE().

◆ extract_velocity()

void oomph::SphericalNavierStokesEquations::extract_velocity ( std::ostream & outfile )

   {
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Vector for coordinates
     Vector<double> x(2);
  
     // Number of plot points
     unsigned Npts = 3;
  
  
     // Assign values of s
     for (unsigned i = 0; i < Npts; i++)
     {
       s[1] = 1.0;
       s[0] = -1 + 2.0 * i / (Npts - 1);
  
  
       // Get position as global coordinate
       interpolated_x(s, x);
       double r = x[0];
       double theta = x[1];
  
       // Output r and theta coordinates
       outfile << r << " ";
       outfile << theta << " ";
  
       // Output the velocity at the current coordinate
       outfile << interpolated_u_spherical_nst(s, 0) << "  ";
       outfile << interpolated_u_spherical_nst(s, 1) << "  ";
       outfile << interpolated_u_spherical_nst(s, 2) << "  ";
       outfile << interpolated_p_spherical_nst(s) << "  ";
  
       // Output shear stress along the given radius
       double shear = interpolated_dudx_spherical_nst(s, 2, 0) -
                      (1 / r) * interpolated_u_spherical_nst(s, 2);
       outfile << shear << "  ";
  
       // Coordinate setup for norm
       double J = J_eulerian(s);
       if (i != 1)
       {
         double w = integral_pt()->weight(i);
         double W = w * J;
         W *= r * r * sin(theta);
  
         // Output the velocity integral at current point
         double u_int = 0.0;
         for (unsigned j = 0; j < 3; j++)
         {
           u_int += interpolated_u_spherical_nst(s, j) * W;
         }
         outfile << u_int << " ";
  
         // Output the pressure integral at current point
         double p_int = interpolated_p_spherical_nst(s) * W;
         outfile << p_int;
         outfile << std::endl;
       }
  
       else
       {
         double w = integral_pt()->weight(5);
         double W = w * J;
         W *= r * r * sin(theta);
  
         // Output the velocity integral at current point
         double u_int = 0.0;
         for (unsigned j = 0; j < 3; j++)
         {
           u_int += interpolated_u_spherical_nst(s, j) * W;
         }
         outfile << u_int << " ";
  
         // Output the pressure integral at current point
         double p_int = interpolated_p_spherical_nst(s) * W;
         outfile << p_int;
         outfile << std::endl;
       }
     }
   }

References i, oomph::FiniteElement::integral_pt(), interpolated_dudx_spherical_nst(), interpolated_p_spherical_nst(), interpolated_u_spherical_nst(), oomph::FiniteElement::interpolated_x(), J, j, oomph::FiniteElement::J_eulerian(), UniformPSDSelfTest::r, s, plotDoE::shear, sin(), BiharmonicTestFunctions2::theta, w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and plotDoE::x.

◆ fill_in_contribution_to_jacobian()

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

inlinevirtual

Add the elemental contribution to the jacobian matrix. and the residuals vector. Note that this function will NOT initialise the residuals vector or the jacobian matrix. It must be called after the residuals vector and jacobian matrix have been initialised to zero. The default is to use finite differences to calculate the jacobian

Reimplemented from oomph::FiniteElement.

     {
       // Call the generic routine with the flag set to 1
       fill_in_generic_residual_contribution_spherical_nst(
         residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
     }

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

Referenced by oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::fill_in_contribution_to_jacobian().

◆ fill_in_contribution_to_jacobian_and_mass_matrix()

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

     {
       // Call the generic routine with the flag set to 2
       fill_in_generic_residual_contribution_spherical_nst(
         residuals, jacobian, mass_matrix, 2);
     }

References fill_in_generic_residual_contribution_spherical_nst().

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

◆ fill_in_contribution_to_residuals()

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

inlinevirtual

Compute the element's residual Vector.

Reimplemented from oomph::GeneralisedElement.

     {
       // Call the generic residuals function with flag set to 0
       // and using a dummy matrix argument
       fill_in_generic_residual_contribution_spherical_nst(
         residuals,
         GeneralisedElement::Dummy_matrix,
         GeneralisedElement::Dummy_matrix,
         0);
     }

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

Referenced by oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::fill_in_contribution_to_residuals().

◆ fill_in_generic_residual_contribution_spherical_nst()

void oomph::SphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst	(	Vector< double > &	residuals,
		DenseMatrix< double > &	jacobian,
		DenseMatrix< double > &	mass_matrix,
		unsigned	flag
	)

protectedvirtual

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

Reimplemented in oomph::RefineableSphericalNavierStokesEquations.

   {
     // Return immediately if there are no dofs
     if (ndof() == 0) return;
  
     // Find out how many nodes there are
     const unsigned n_node = nnode();
  
     // Get continuous time from timestepper of first node
     double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
  
     // Find out how many pressure dofs there are
     const unsigned n_pres = npres_spherical_nst();
  
     // Find the indices at which the local velocities are stored
     unsigned u_nodal_index[3];
     for (unsigned i = 0; i < 3; i++)
     {
       u_nodal_index[i] = u_index_spherical_nst(i);
     }
  
     // Set up memory for the shape and test functions
     Shape psif(n_node), testf(n_node);
     DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
  
     // Set up memory for pressure shape and test functions
     Shape psip(n_pres), testp(n_pres);
  
     // Number of integration points
     const unsigned n_intpt = integral_pt()->nweight();
  
     // Set the Vector to hold local coordinates
     Vector<double> s(2);
  
     // Get Physical Variables from Element
     // Reynolds number must be multiplied by the density ratio
     const double scaled_re = re() * density_ratio();
     const double scaled_re_st = re_st() * density_ratio();
     const double scaled_re_inv_ro = re_invro() * density_ratio();
     // double scaled_re_inv_fr = re_invfr()*density_ratio();
     // double visc_ratio = viscosity_ratio();
     Vector<double> G = g();
  
     // Integers to store the local equations and unknowns
     int local_eqn = 0, local_unknown = 0;
  
     // Loop over the integration points
     for (unsigned ipt = 0; ipt < n_intpt; ipt++)
     {
       // Assign values of s
       for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
       // Get the integral weight
       double w = integral_pt()->weight(ipt);
  
       // Call the derivatives of the shape and test functions
       double J = dshape_and_dtest_eulerian_at_knot_spherical_nst(
         ipt, psif, dpsifdx, testf, dtestfdx);
  
       // Call the pressure shape and test functions
       pshape_spherical_nst(s, psip, testp);
  
       // Premultiply the weights and the Jacobian
       const double W = w * J;
  
       // Calculate local values of the pressure and velocity components
       // Allocate
       double interpolated_p = 0.0;
       Vector<double> interpolated_u(3, 0.0);
       Vector<double> interpolated_x(2, 0.0);
       Vector<double> mesh_velocity(2, 0.0);
       Vector<double> dudt(3, 0.0);
       DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
  
       // Calculate pressure
       for (unsigned l = 0; l < n_pres; l++)
       {
         interpolated_p += p_spherical_nst(l) * psip(l);
       }
  
  
       // Calculate velocities and derivatives:
  
       // Loop over nodes
       for (unsigned l = 0; l < n_node; l++)
       {
         // cache the shape functions
         double psi_ = psif(l);
         // Loop over the positions separately
         for (unsigned i = 0; i < 2; i++)
         {
           interpolated_x[i] += raw_nodal_position(l, i) * psi_;
         }
  
         // Loop over velocity directions (three of these)
         for (unsigned i = 0; i < 3; i++)
         {
           // Get the nodal value
           const double u_value = raw_nodal_value(l, u_nodal_index[i]);
           interpolated_u[i] += u_value * psi_;
           dudt[i] += du_dt_spherical_nst(l, i) * psi_;
  
           // Loop over derivative directions
           for (unsigned j = 0; j < 2; j++)
           {
             interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
           }
         }
       }
  
       if (!ALE_is_disabled)
       {
         OomphLibWarning("ALE is not tested for spherical Navier Stokes",
                         "SphericalNS::fill_in...",
                         OOMPH_EXCEPTION_LOCATION);
         // Loop over nodes
         for (unsigned l = 0; l < n_node; l++)
         {
           // Loop over directions (only DIM (2) of these)
           for (unsigned i = 0; i < 2; i++)
           {
             mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif(l);
           }
         }
       }
  
       // Get the user-defined body force terms
       Vector<double> body_force(3);
       get_body_force_spherical_nst(time, ipt, s, interpolated_x, body_force);
  
       // Get the user-defined source function
       // double source = get_source_spherical_nst(time(),ipt,interpolated_x);
  
  
       // MOMENTUM EQUATIONS
       //------------------
  
       const double r = interpolated_x[0];
       const double r2 = r * r;
       const double sin_theta = sin(interpolated_x[1]);
       const double cos_theta = cos(interpolated_x[1]);
       const double cot_theta = cos_theta / sin_theta;
  
       const double u_r = interpolated_u[0];
       const double u_theta = interpolated_u[1];
       const double u_phi = interpolated_u[2];
  
       // Loop over the test functions
       for (unsigned l = 0; l < n_node; l++)
       {
         // Cannot loop over the velocity components,
         // so address each of them separately
  
         // FIRST: r-component momentum equation
  
         unsigned i = 0;
  
         // If it's not a boundary condition
         local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
         if (local_eqn >= 0)
         {
           // Convective r-terms
           double conv = u_r * interpolated_dudx(0, 0) * r;
  
           // Convective theta-terms
           conv += u_theta * interpolated_dudx(0, 1);
  
           // Remaining convective terms
           conv -= (u_theta * u_theta + u_phi * u_phi);
  
           // Add the time-derivative term and the convective terms
           // to the sum
           double sum = (scaled_re_st * r2 * dudt[0] + r * scaled_re * conv) *
                        sin_theta * testf[l];
  
           // Subtract the mesh velocity terms
           if (!ALE_is_disabled)
           {
             sum -= scaled_re_st *
                    (mesh_velocity[0] * interpolated_dudx(0, 0) * r +
                     mesh_velocity[1] * (interpolated_dudx(0, 1) - u_theta)) *
                    r * sin_theta * testf[l];
           }
  
           // Add the Coriolis term
           sum -= 2.0 * scaled_re_inv_ro * sin_theta * u_phi * r2 * sin_theta *
                  testf[l];
  
           // r-derivative test function component of stress tensor
           sum += (-interpolated_p + 2.0 * interpolated_dudx(0, 0)) * r2 *
                  sin_theta * dtestfdx(l, 0);
  
           // theta-derivative test function component of stress tensor
           sum +=
             (r * interpolated_dudx(1, 0) - u_theta + interpolated_dudx(0, 1)) *
             sin_theta * dtestfdx(l, 1);
  
           // Undifferentiated test function component of stress tensor
           sum += 2.0 *
                  ((-r * interpolated_p + interpolated_dudx(1, 1) + 2.0 * u_r) *
                     sin_theta +
                   u_theta * cos_theta) *
                  testf(l);
  
           // Make the residuals negative for consistency with the
           // other Navier-Stokes equations
           residuals[local_eqn] -= sum * W;
  
           // CALCULATE THE JACOBIAN
           if (flag)
           {
             // Loop over the velocity shape functions again
             for (unsigned l2 = 0; l2 < n_node; l2++)
             {
               // Cannot loop over the velocity components
               // --- need to compute these separately instead
  
               unsigned i2 = 0;
  
               // If at a non-zero degree of freedom add in the entry
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
               if (local_unknown >= 0)
               {
                 // Convective r-term contribution
                 double jac_conv = r * (psif(l2) * interpolated_dudx(0, 0) +
                                        u_r * dpsifdx(l2, 0));
  
                 // Convective theta-term contribution
                 jac_conv += u_theta * dpsifdx(l2, 1);
  
                 // Add the time-derivative contribution and the convective
                 // contribution to the sum
                 double jac_sum =
                   (scaled_re_st * node_pt(l2)->time_stepper_pt()->weight(1, 0) *
                      psif(l2) * r2 +
                    scaled_re * jac_conv * r) *
                   testf(l);
  
                 // Subtract the mesh-velocity terms
                 if (!ALE_is_disabled)
                 {
                   jac_sum -= scaled_re_st *
                              (mesh_velocity[0] * dpsifdx(l2, 0) * r +
                               mesh_velocity[1] * dpsifdx(l2, 1)) *
                              r * sin_theta * testf(l);
                 }
  
  
                 // Contribution from the r-derivative test function part
                 // of stress tensor
                 jac_sum += 2.0 * dpsifdx(l2, 0) * dtestfdx(l, 0) * r2;
  
                 // Contribution from the theta-derivative test function part
                 // of stress tensor
                 jac_sum += dpsifdx(l2, 1) * dtestfdx(l, 1);
  
  
                 // Contribution from the undifferentiated test function part
                 // of stress tensor
                 jac_sum += 4.0 * psif[l2] * testf(l);
  
                 // Add the total contribution to the jacobian multiplied
                 // by the jacobian weight
                 jacobian(local_eqn, local_unknown) -= jac_sum * sin_theta * W;
  
                 // Mass matrix term
                 if (flag == 2)
                 {
                   mass_matrix(local_eqn, local_unknown) +=
                     scaled_re_st * psif[l2] * testf[l] * r2 * sin_theta * W;
                 }
               }
               // End of i2=0 section
  
  
               i2 = 1;
  
               // If at a non-zero degree of freedom add in the entry
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
               if (local_unknown >= 0)
               {
                 // No time derivative contribution
  
                 // Convective contribution
                 double jac_sum = scaled_re *
                                  (interpolated_dudx(0, 1) - 2.0 * u_theta) *
                                  psif(l2) * r * sin_theta * testf(l);
  
                 // Mesh velocity terms
                 if (!ALE_is_disabled)
                 {
                   jac_sum += scaled_re_st * mesh_velocity[1] * psif(l2) * r *
                              sin_theta * testf(l);
                 }
  
                 // Contribution from the theta-derivative test function
                 // part of stress tensor
                 jac_sum +=
                   (r * dpsifdx(l2, 0) - psif(l2)) * dtestfdx(l, 1) * sin_theta;
  
                 // Contribution from the undifferentiated test function
                 // part of stress tensor
                 jac_sum += 2.0 *
                            (dpsifdx(l2, 1) * sin_theta + psif(l2) * cos_theta) *
                            testf(l);
  
                 // Add the full contribution to the jacobian matrix
                 jacobian(local_eqn, local_unknown) -= jac_sum * W;
  
               } // End of i2=1 section
  
  
               i2 = 2;
  
               // If at a non-zero degree of freedom add in the entry
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
               if (local_unknown >= 0)
               {
                 // Single convective-term contribution
                 jacobian(local_eqn, local_unknown) += 2.0 * scaled_re * u_phi *
                                                       psif(l2) * r * sin_theta *
                                                       testf[l] * W;
  
                 // Coriolis term
                 jacobian(local_eqn, local_unknown) +=
                   2.0 * scaled_re_inv_ro * sin_theta * psif(l2) * r2 *
                   sin_theta * testf[l] * W;
               }
               // End of i2=2 section
             }
             // End of the l2 loop over the test functions
  
             // Now loop over pressure shape functions
             // This is the contribution from pressure gradient
             for (unsigned l2 = 0; l2 < n_pres; l2++)
             {
               /*If we are at a non-zero degree of freedom in the entry*/
               local_unknown = p_local_eqn(l2);
               if (local_unknown >= 0)
               {
                 jacobian(local_eqn, local_unknown) +=
                   psip(l2) * (r2 * dtestfdx(l, 0) + 2.0 * r * testf(l)) *
                   sin_theta * W;
               }
             }
           } // End of Jacobian calculation
  
         } // End of if not boundary condition statement
           // End of r-component momentum equation
  
  
         // SECOND: theta-component momentum equation
  
         i = 1;
  
         // IF it's not a boundary condition
         local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
         if (local_eqn >= 0)
         {
           // All convective terms
           double conv = (u_r * interpolated_dudx(1, 0) * r +
                          u_theta * interpolated_dudx(1, 1) + u_r * u_theta) *
                           sin_theta -
                         u_phi * u_phi * cos_theta;
  
           // Add the time-derivative and convective terms to the
           // residuals sum
           double sum =
             (scaled_re_st * dudt[1] * r2 * sin_theta + scaled_re * r * conv) *
             testf(l);
  
           // Subtract the mesh velocity terms
           if (!ALE_is_disabled)
           {
             sum -= scaled_re_st *
                    (mesh_velocity[0] * interpolated_dudx(1, 0) * r +
                     mesh_velocity[1] * (interpolated_dudx(1, 1) + u_r)) *
                    r * sin_theta * testf(l);
           }
  
           // Add the Coriolis term
           sum -= 2.0 * scaled_re_inv_ro * cos_theta * u_phi * r2 * sin_theta *
                  testf[l];
  
           // r-derivative test function component of stress tensor
           sum +=
             (r * interpolated_dudx(1, 0) - u_theta + interpolated_dudx(0, 1)) *
             r * sin_theta * dtestfdx(l, 0);
  
           // theta-derivative test function component of stress tensor
           sum +=
             (-r * interpolated_p + 2.0 * interpolated_dudx(1, 1) + 2.0 * u_r) *
             dtestfdx(l, 1) * sin_theta;
  
           // Undifferentiated test function component of stress tensor
           sum -=
             ((r * interpolated_dudx(1, 0) - u_theta + interpolated_dudx(0, 1)) *
                sin_theta -
              (-r * interpolated_p + 2.0 * u_r + 2.0 * u_theta * cot_theta) *
                cos_theta) *
             testf(l);
  
           // Add the sum to the residuals,
           //(Negative for consistency)
           residuals[local_eqn] -= sum * W;
  
           // CALCULATE THE JACOBIAN
           if (flag)
           {
             // Loop over the velocity shape functions again
             for (unsigned l2 = 0; l2 < n_node; l2++)
             {
               // Cannot loop over the velocity components
               // --- need to compute these separately instead
  
               unsigned i2 = 0;
  
               // If at a non-zero degree of freedom add in the entry
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
               if (local_unknown >= 0)
               {
                 // No time derivative contribution
  
                 // Convective terms
                 double jac_sum = scaled_re *
                                  (r2 * interpolated_dudx(1, 0) + r * u_theta) *
                                  psif(l2) * sin_theta * testf(l);
  
                 // Mesh velocity terms
                 if (!ALE_is_disabled)
                 {
                   jac_sum -= scaled_re_st * mesh_velocity[1] * psif(l2) * r *
                              sin_theta * testf(l);
                 }
  
                 // Contribution from the r-derivative
                 // test function part of stress tensor
                 jac_sum += dpsifdx(l2, 1) * dtestfdx(l, 0) * sin_theta * r;
  
                 // Contribution from the theta-derivative test function
                 // part of stress tensor
                 jac_sum += 2.0 * psif(l2) * dtestfdx(l, 1) * sin_theta;
  
                 // Contribution from the undifferentiated test function
                 // part of stress tensor
                 jac_sum -=
                   (dpsifdx(l2, 1) * sin_theta - 2.0 * psif(l2) * cos_theta) *
                   testf(l);
  
                 // Add the sum to the jacobian
                 jacobian(local_eqn, local_unknown) -= jac_sum * W;
               }
               // End of i2=0 section
  
  
               i2 = 1;
  
               // If at a non-zero degree of freedom add in the entry
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
               if (local_unknown >= 0)
               {
                 // Contribution from the convective terms
                 double jac_conv = r * u_r * dpsifdx(l2, 0) +
                                   u_theta * dpsifdx(l2, 1) +
                                   (interpolated_dudx(1, 1) + u_r) * psif(l2);
  
                 // Add the time-derivative term and the
                 // convective terms
                 double jac_sum =
                   (scaled_re_st * node_pt(l2)->time_stepper_pt()->weight(1, 0) *
                      psif(l2) * r2 +
                    scaled_re * r * jac_conv) *
                   testf(l) * sin_theta;
  
  
                 // Mesh velocity terms
                 if (!ALE_is_disabled)
                 {
                   jac_sum -= scaled_re_st *
                              (mesh_velocity[0] * dpsifdx(l2, 0) * r +
                               mesh_velocity[1] * dpsifdx(l2, 1)) *
                              r * sin_theta * testf(l);
                 }
  
                 // Contribution from the r-derivative test function
                 // part of stress tensor
                 jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) * r *
                            dtestfdx(l, 0) * sin_theta;
  
                 // Contribution from the theta-derivative test function
                 // part of stress tensor
                 jac_sum += 2.0 * dpsifdx(l2, 1) * dtestfdx(l, 1) * sin_theta;
  
                 // Contribution from the undifferentiated test function
                 // part of stress tensor
                 jac_sum -= ((r * dpsifdx(l2, 0) - psif(l2)) * sin_theta -
                             2.0 * psif(l2) * cot_theta * cos_theta) *
                            testf(l);
  
                 // Add the contribution to the jacobian
                 jacobian(local_eqn, local_unknown) -= jac_sum * W;
  
                 // Mass matrix term
                 if (flag == 2)
                 {
                   mass_matrix(local_eqn, local_unknown) +=
                     scaled_re_st * psif[l2] * testf[l] * r2 * sin_theta * W;
                 }
               }
               // End of i2=1 section
  
  
               i2 = 2;
  
               // If at a non-zero degree of freedom add in the entry
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
               if (local_unknown >= 0)
               {
                 // Only a convective term contribution
                 jacobian(local_eqn, local_unknown) += scaled_re * 2.0 * r *
                                                       psif(l2) * u_phi *
                                                       cos_theta * testf(l) * W;
  
                 // Coriolis term
                 jacobian(local_eqn, local_unknown) +=
                   2.0 * scaled_re_inv_ro * cos_theta * psif(l2) * r2 *
                   sin_theta * testf[l] * W;
               }
               // End of i2=2 section
             }
             // End of the l2 loop over the test functions
  
             /*Now loop over pressure shape functions*/
             /*This is the contribution from pressure gradient*/
             for (unsigned l2 = 0; l2 < n_pres; l2++)
             {
               /*If we are at a non-zero degree of freedom in the entry*/
               local_unknown = p_local_eqn(l2);
               if (local_unknown >= 0)
               {
                 jacobian(local_eqn, local_unknown) +=
                   psip(l2) * r *
                   (dtestfdx(l, 1) * sin_theta + cos_theta * testf[l]) * W;
               }
             }
           } /*End of Jacobian calculation*/
  
         } // End of if not boundary condition statement
           // End of theta-component of momentum
  
  
         // THIRD: phi-component momentum equation
  
         i = 2;
  
         // IF it's not a boundary condition
         local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
         if (local_eqn >= 0)
         {
           // Convective r-terms
           double conv = u_r * interpolated_dudx(2, 0) * r * sin_theta;
  
           // Convective theta-terms
           conv += u_theta * interpolated_dudx(2, 1) * sin_theta;
  
           // Remaining convective terms
           conv += u_phi * (u_r * sin_theta + u_theta * cos_theta);
  
           // Add the time-derivative term and the convective terms to the sum
           double sum =
             (scaled_re_st * r2 * dudt[2] * sin_theta + scaled_re * conv * r) *
             testf(l);
  
           // Mesh velocity terms
           if (!ALE_is_disabled)
           {
             sum -= scaled_re_st *
                    (mesh_velocity[0] * interpolated_dudx(2, 0) * r +
                     mesh_velocity[1] * interpolated_dudx(2, 1)) *
                    r * sin_theta * testf(l);
           }
  
           // Add the Coriolis term
           sum += 2.0 * scaled_re_inv_ro *
                  (cos_theta * u_theta + sin_theta * u_r) * r2 * sin_theta *
                  testf[l];
  
           // r-derivative test function component of stress tensor
           sum += (r2 * interpolated_dudx(2, 0) - r * u_phi) * dtestfdx(l, 0) *
                  sin_theta;
  
           // theta-derivative test function component of stress tensor
           sum += (interpolated_dudx(2, 1) * sin_theta - u_phi * cos_theta) *
                  dtestfdx(l, 1);
  
           // Undifferentiated test function component of stress tensor
           sum -= ((r * interpolated_dudx(2, 0) - u_phi) * sin_theta +
                   (interpolated_dudx(2, 1) - u_phi * cot_theta) * cos_theta) *
                  testf(l);
  
           // Add the sum to the residuals
           residuals[local_eqn] -= sum * W;
  
  
           // CALCULATE THE JACOBIAN
           if (flag)
           {
             // Loop over the velocity shape functions again
             for (unsigned l2 = 0; l2 < n_node; l2++)
             {
               // Cannot loop over the velocity components
               // -- need to compute these separately instead
  
               unsigned i2 = 0;
  
               // If at a non-zero degree of freedom add in the entry
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
               if (local_unknown >= 0)
               {
                 // Contribution from convective terms
                 jacobian(local_eqn, local_unknown) -=
                   scaled_re * (r * interpolated_dudx(2, 0) + u_phi) * psif(l2) *
                   testf(l) * r * sin_theta * W;
  
                 // Coriolis term
                 jacobian(local_eqn, local_unknown) -=
                   2.0 * scaled_re_inv_ro * sin_theta * psif(l2) * r2 *
                   sin_theta * testf[l] * W;
               }
               // End of i2=0 section
  
               i2 = 1;
  
               // If at a non-zero degree of freedom add in the entry
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
               if (local_unknown >= 0)
               {
                 // Contribution from convective terms
                 jacobian(local_eqn, local_unknown) -=
                   scaled_re *
                   (interpolated_dudx(2, 1) * sin_theta + u_phi * cos_theta) *
                   r * psif(l2) * testf(l) * W;
  
  
                 // Coriolis term
                 jacobian(local_eqn, local_unknown) -=
                   2.0 * scaled_re_inv_ro * cos_theta * psif(l2) * r2 *
                   sin_theta * testf[l] * W;
               }
               // End of i2=1 section
  
  
               i2 = 2;
  
               // If at a non-zero degree of freedom add in the entry
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
               if (local_unknown >= 0)
               {
                 // Convective terms
                 double jac_conv = r * u_r * dpsifdx(l2, 0) * sin_theta;
  
                 // Convective theta-term contribution
                 jac_conv += u_theta * dpsifdx(l2, 1) * sin_theta;
  
                 // Contribution from the remaining convective terms
                 jac_conv += (u_r * sin_theta + u_theta * cos_theta) * psif(l2);
  
                 // Add the convective and time derivatives
                 double jac_sum =
                   (scaled_re_st * node_pt(l2)->time_stepper_pt()->weight(1, 0) *
                      psif(l2) * r2 * sin_theta +
                    scaled_re * r * jac_conv) *
                   testf(l);
  
  
                 // Mesh velocity terms
                 if (!ALE_is_disabled)
                 {
                   jac_sum -= scaled_re_st *
                              (mesh_velocity[0] * dpsifdx(l2, 0) * r +
                               mesh_velocity[1] * dpsifdx(l2, 1)) *
                              r * sin_theta * testf(l);
                 }
  
                 // Contribution from the r-derivative test function
                 // part of stress tensor
                 jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) * dtestfdx(l, 0) *
                            r * sin_theta;
  
                 // Contribution from the theta-derivative test function
                 // part of stress tensor
                 jac_sum += (dpsifdx(l2, 1) * sin_theta - psif(l2) * cos_theta) *
                            dtestfdx(l, 1);
  
                 // Contribution from the undifferentiated test function
                 // part of stress tensor
                 jac_sum -=
                   ((r * dpsifdx(l2, 0) - psif(l2)) * sin_theta +
                    (dpsifdx(l2, 1) - psif(l2) * cot_theta) * cos_theta) *
                   testf(l);
  
                 // Add to the jacobian
                 jacobian(local_eqn, local_unknown) -= jac_sum * W;
  
                 // Mass matrix term
                 if (flag == 2)
                 {
                   mass_matrix(local_eqn, local_unknown) +=
                     scaled_re_st * psif(l2) * testf[l] * r2 * sin_theta * W;
                 }
               }
               // End of i2=2 section
             }
             // End of the l2 loop over the test functions
  
             // We assume phi-derivatives to be zero
             // No phi-contribution to the pressure gradient to include here
  
           } // End of Jacobian calculation
  
         } // End of if not boundary condition statement
           // End of phi-component of momentum
  
       } // End of loop over shape functions
  
  
       // CONTINUITY EQUATION
       //-------------------
  
       // Loop over the shape functions
       for (unsigned l = 0; l < n_pres; l++)
       {
         local_eqn = p_local_eqn(l);
         // If not a boundary conditions
         // Can't loop over the velocity components, so put these in separately
         if (local_eqn >= 0)
         {
           residuals[local_eqn] += ((2.0 * u_r + r * interpolated_dudx(0, 0) +
                                     interpolated_dudx(1, 1)) *
                                      sin_theta +
                                    u_theta * cos_theta) *
                                   r * testp(l) * W;
  
           // CALCULATE THE JACOBIAN
           if (flag)
           {
             // Loop over the velocity shape functions
             for (unsigned l2 = 0; l2 < n_node; l2++)
             {
               // Can't loop over velocity components, so put these in separately
  
               // Contribution from the r-component
               unsigned i2 = 0;
               /*If we're at a non-zero degree of freedom add it in*/
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
  
               if (local_unknown >= 0)
               {
                 jacobian(local_eqn, local_unknown) +=
                   (2.0 * psif(l2) + r * dpsifdx(l2, 0)) * r * sin_theta *
                   testp(l) * W;
               }
               // End of contribution from r-component
  
  
               // Contribution from the theta-component
               i2 = 1;
               /*If we're at a non-zero degree of freedom add it in*/
               local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
  
               if (local_unknown >= 0)
               {
                 jacobian(local_eqn, local_unknown) +=
                   r * (dpsifdx(l2, 1) * sin_theta + psif(l2) * cos_theta) *
                   testp(l) * W;
               }
               // End of contribution from theta-component
  
             } // End of loop over l2
           } // End of Jacobian calculation
  
         } // End of if not boundary condition
       } // End of loop over l
     }
   }

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

◆ fix_pressure()

virtual void oomph::SphericalNavierStokesEquations::fix_pressure	(	const unsigned &	p_dof,
		const double &	p_value
	)

pure virtual

Pin p_dof-th pressure dof and set it to value specified by p_value.

Implemented in oomph::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

◆ full_output() [1/2]

void oomph::SphericalNavierStokesEquations::full_output ( std::ostream & outfile )

inline

Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format. Default number of plot points

     {
       unsigned nplot = 5;
       full_output(outfile, nplot);
     }

Referenced by oomph::QSphericalCrouzeixRaviartElement::full_output().

◆ full_output() [2/2]

void oomph::SphericalNavierStokesEquations::full_output	(	std::ostream &	outfile,
		const unsigned &	nplot
	)

Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format. nplot points in each coordinate direction

Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format. Specified number of plot points in each coordinate direction

   {
     throw OomphLibError(
       "Probably Broken", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
  
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Acceleration
     Vector<double> dudt(3);
  
     // Mesh elocity
     Vector<double> mesh_veloc(3, 0.0);
  
     // Tecplot header info
     outfile << tecplot_zone_string(nplot);
  
     // Find out how many nodes there are
     unsigned n_node = nnode();
  
     // Set up memory for the shape functions
     Shape psif(n_node);
     DShape dpsifdx(n_node, 2);
  
     // Loop over plot points
     unsigned num_plot_points = nplot_points(nplot);
     for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
     {
       // Get local coordinates of plot point
       get_s_plot(iplot, nplot, s);
  
       // Call the derivatives of the shape and test functions
       dshape_eulerian(s, psif, dpsifdx);
  
       // Allocate storage
       Vector<double> mesh_veloc(3);
       Vector<double> dudt(3);
       Vector<double> dudt_ALE(3);
       DenseMatrix<double> interpolated_dudx(3, 2);
  
       // Initialise everything to zero
       for (unsigned i = 0; i < 3; i++)
       {
         mesh_veloc[i] = 0.0;
         dudt[i] = 0.0;
         dudt_ALE[i] = 0.0;
         for (unsigned j = 0; j < 2; j++)
         {
           interpolated_dudx(i, j) = 0.0;
         }
       }
  
       // Calculate velocities and derivatives
  
  
       // Loop over directions
       for (unsigned i = 0; i < 3; i++)
       {
         // Get the index at which velocity i is stored
         unsigned u_nodal_index = u_index_spherical_nst(i);
         // Loop over nodes
         for (unsigned l = 0; l < n_node; l++)
         {
           dudt[i] += du_dt_spherical_nst(l, u_nodal_index) * psif(l);
           mesh_veloc[i] += dnodal_position_dt(l, i) * psif(l);
  
           // Loop over derivative directions for velocity gradients
           for (unsigned j = 0; j < 2; j++)
           {
             interpolated_dudx(i, j) +=
               nodal_value(l, u_nodal_index) * dpsifdx(l, j);
           }
         }
       }
  
  
       // Get dudt in ALE form (incl mesh veloc)
       for (unsigned i = 0; i < 3; i++)
       {
         dudt_ALE[i] = dudt[i];
         for (unsigned k = 0; k < 2; k++)
         {
           dudt_ALE[i] -= mesh_veloc[k] * interpolated_dudx(i, k);
         }
       }
  
  
       // Coordinates
       for (unsigned i = 0; i < 2; i++)
       {
         outfile << interpolated_x(s, i) << " ";
       }
  
       // Velocities
       for (unsigned i = 0; i < 3; i++)
       {
         outfile << interpolated_u_spherical_nst(s, i) << " ";
       }
  
       // Pressure
       outfile << interpolated_p_spherical_nst(s) << " ";
  
       // Accelerations
       for (unsigned i = 0; i < 3; i++)
       {
         outfile << dudt_ALE[i] << " ";
       }
  
       // Dissipation
       outfile << dissipation(s) << " ";
  
  
       outfile << std::endl;
     }
  
     // Write tecplot footer (e.g. FE connectivity lists)
     write_tecplot_zone_footer(outfile, nplot);
   }

◆ g()

const Vector<double>& oomph::SphericalNavierStokesEquations::g ( ) const

inline

Vector of gravitational components.

     {
       return *G_pt;
     }

References G_pt.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), fill_in_generic_residual_contribution_spherical_nst(), oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::fill_in_off_diagonal_jacobian_blocks_analytic(), and oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::get_body_force_spherical_nst().

◆ g_pt()

Vector<double>*& oomph::SphericalNavierStokesEquations::g_pt ( )

inline

Pointer to Vector of gravitational components.

     {
       return G_pt;
     }

References G_pt.

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

◆ get_body_force_spherical_nst()

virtual void oomph::SphericalNavierStokesEquations::get_body_force_spherical_nst	(	const double &	time,
		const unsigned &	ipt,
		const Vector< double > &	s,
		const Vector< double > &	x,
		Vector< double > &	result
	)

inlineprotectedvirtual

Calculate the body force at a given time and local and/or Eulerian position. This function is virtual so that it can be overloaded in multi-physics elements where the body force might depend on another variable.

Reimplemented in oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement.

     {
       // If the function pointer is zero return zero
       if (Body_force_fct_pt == 0)
       {
         // Loop over three spatial dimensions and set body forces to zero
         for (unsigned i = 0; i < 3; i++)
         {
           result[i] = 0.0;
         }
       }
       // Otherwise call the function
       else
       {
         (*Body_force_fct_pt)(time, x, result);
       }
     }

References Body_force_fct_pt, i, and plotDoE::x.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), and fill_in_generic_residual_contribution_spherical_nst().

◆ get_load()

void oomph::SphericalNavierStokesEquations::get_load	(	const Vector< double > &	s,
		const Vector< double > &	N,
		Vector< double > &	load
	)

inlinevirtual

This implements a pure virtual function defined in the FSIFluidElement class. The function computes the traction (on the viscous scale), at the element's local coordinate s, that the fluid element exerts onto an adjacent solid element. The number of arguments is imposed by the interface defined in the FSIFluidElement – only the unit normal N (pointing into the fluid!) is actually used in the computation.

Implements oomph::FSIFluidElement.

     {
       // Note: get_traction() computes the traction onto the fluid
       // if N is the outer unit normal onto the fluid; here we're
       // exepcting N to point into the fluid so we're getting the
       // traction onto the adjacent wall instead!
       get_traction(s, N, load);
     }

References get_traction(), load(), N, and s.

◆ get_pressure_and_velocity_mass_matrix_diagonal()

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

   {
 #ifdef PARANOID
     if ((which_one == 0) || (which_one == 1))
     {
       throw OomphLibError("Computation of diagonal of pressure mass matrix is "
                           "not impmented yet.\n",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Resize and initialise
     veloc_mass_diag.assign(ndof(), 0.0);
  
     // find out how many nodes there are
     const unsigned n_node = nnode();
  
     // find number of coordinates
     const unsigned n_value = 3;
  
     // find the indices at which the local velocities are stored
     Vector<unsigned> u_nodal_index(n_value);
     for (unsigned i = 0; i < n_value; i++)
     {
       u_nodal_index[i] = this->u_index_spherical_nst(i);
     }
  
     // Set up memory for test functions
     Shape test(n_node);
  
     // Number of integration points
     unsigned n_intpt = integral_pt()->nweight();
  
     // Integer to store the local equations no
     int local_eqn = 0;
  
     // Loop over the integration points
     for (unsigned ipt = 0; ipt < n_intpt; ipt++)
     {
       // Get the integral weight
       double w = integral_pt()->weight(ipt);
  
       // Get determinant of Jacobian of the mapping
       double J = J_eulerian_at_knot(ipt);
  
       // Premultiply weights and Jacobian
       double W = w * J;
  
       // Get the velocity test functions - there is no explicit
       // function to give the test function so use shape
       shape_at_knot(ipt, test);
  
       // Need to get the position to sort out the jacobian properly
       double r = 0.0;
       double theta = 0.0;
       for (unsigned l = 0; l < n_node; l++)
       {
         r += this->nodal_position(l, 0) * test(l);
         theta += this->nodal_position(l, 1) * test(l);
       }
  
       // Multiply by the geometric factor
       W *= r * r * sin(theta);
  
       // Loop over the veclocity test functions
       for (unsigned l = 0; l < n_node; l++)
       {
         // Loop over the velocity components
         for (unsigned i = 0; i < n_value; i++)
         {
           local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
  
           // If not a boundary condition
           if (local_eqn >= 0)
           {
             // Add the contribution
             veloc_mass_diag[local_eqn] += test[l] * test[l] * W;
           } // End of if not boundary condition statement
         } // End of loop over dimension
       } // End of loop over test functions
     }
   }

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(), sin(), Eigen::test, BiharmonicTestFunctions2::theta, u_index_spherical_nst(), w, oomph::QuadTreeNames::W, and oomph::Integral::weight().

◆ get_source_spherical_nst()

virtual double oomph::SphericalNavierStokesEquations::get_source_spherical_nst	(	double	time,
		const unsigned &	ipt,
		const Vector< double > &	x
	)

inlineprotectedvirtual

Calculate the source fct at given time and Eulerian position

     {
       // If the function pointer is zero return zero
       if (Source_fct_pt == 0)
       {
         return 0;
       }
       // Otherwise call the function
       else
       {
         return (*Source_fct_pt)(time, x);
       }
     }

References Source_fct_pt, and plotDoE::x.

◆ get_traction()

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

Compute traction (on the viscous scale) exerted onto the fluid at local coordinate s. N has to be outer unit normal to the fluid.

   {
     // Get velocity gradients
     DenseMatrix<double> strainrate(3, 3);
     strain_rate(s, strainrate);
  
     // Get pressure
     double press = interpolated_p_spherical_nst(s);
  
     // Loop over traction components
     for (unsigned i = 0; i < 3; i++)
     {
       // If we are in the r and theta direction
       // initialise the traction
       if (i < 2)
       {
         traction[i] = -press * N[i];
       }
       // Otherwise it's zero because the normal cannot have a component
       // in the phi direction
       else
       {
         traction[i] = 0.0;
       }
       // Loop over the possible traction directions
       for (unsigned j = 0; j < 2; j++)
       {
         traction[i] += 2.0 * strainrate(i, j) * N[j];
       }
     }
   }

References i, interpolated_p_spherical_nst(), j, N, s, and strain_rate().

Referenced by get_load().

◆ get_vorticity()

void oomph::SphericalNavierStokesEquations::get_vorticity	(	const Vector< double > &	s,
		Vector< double > &	vorticity
	)		const

Compute the vorticity vector at local coordinate s.

Compute 2D vorticity vector at local coordinate s (return in one and only component of vorticity vector

   {
     throw OomphLibError("Not tested for spherical Navier-Stokes elements",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
  
 #ifdef PARANOID
     if (vorticity.size() != 1)
     {
       std::ostringstream error_message;
       error_message << "Computation of vorticity in 2D requires a 1D vector\n"
                     << "which contains the only non-zero component of the\n"
                     << "vorticity vector. You've passed a vector of size "
                     << vorticity.size() << std::endl;
  
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Specify spatial dimension
     unsigned DIM = 2;
  
     // Velocity gradient matrix
     DenseMatrix<double> dudx(DIM);
  
     // Find out how many nodes there are in the element
     unsigned n_node = nnode();
  
     // Set up memory for the shape and test functions
     Shape psi(n_node);
     DShape dpsidx(n_node, DIM);
  
     // Call the derivatives of the shape functions
     dshape_eulerian(s, psi, dpsidx);
  
     // Initialise to zero
     for (unsigned i = 0; i < DIM; i++)
     {
       for (unsigned j = 0; j < DIM; j++)
       {
         dudx(i, j) = 0.0;
       }
     }
  
     // Loop over veclocity components
     for (unsigned i = 0; i < DIM; i++)
     {
       // Get the index at which the i-th velocity is stored
       unsigned u_nodal_index = u_index_spherical_nst(i);
       // Loop over derivative directions
       for (unsigned j = 0; j < DIM; j++)
       {
         // Loop over nodes
         for (unsigned l = 0; l < n_node; l++)
         {
           dudx(i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
         }
       }
     }
  
     // Z-component of vorticity
     vorticity[0] = dudx(1, 0) - dudx(0, 1);
   }

References DIM, oomph::FiniteElement::dshape_eulerian(), i, j, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, and u_index_spherical_nst().

Referenced by output_vorticity().

◆ interpolated_dudx_spherical_nst()

double oomph::SphericalNavierStokesEquations::interpolated_dudx_spherical_nst	(	const Vector< double > &	s,
		const unsigned &	i,
		const unsigned &	j
	)		const

inline

Return matrix entry dudx(i,j) of the FE interpolated velocity derivative at local coordinate s

     {
       // Find number of nodes
       unsigned n_node = nnode();
       // Local shape function
       Shape psi(n_node);
       DShape dpsidx(n_node, 2);
       // Find values of shape function
       (void)dshape_eulerian(s, psi, dpsidx);
  
       // Get nodal index at which i-th velocity is stored
       unsigned u_nodal_index = u_index_spherical_nst(i);
  
       // Initialise value of u
       double interpolated_dudx = 0.0;
       // Loop over the local nodes and sum
       for (unsigned l = 0; l < n_node; l++)
       {
         interpolated_dudx += nodal_value(l, u_nodal_index) * dpsidx(l, j);
       }
  
       return (interpolated_dudx);
     }

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

Referenced by compute_error_e(), compute_shear_stress(), and extract_velocity().

◆ interpolated_p_spherical_nst()

double oomph::SphericalNavierStokesEquations::interpolated_p_spherical_nst ( const Vector< double > & s ) const

inline

Return FE interpolated pressure at local coordinate s.

     {
       // Find number of nodes
       unsigned n_pres = npres_spherical_nst();
       // Local shape function
       Shape psi(n_pres);
       // Find values of shape function
       pshape_spherical_nst(s, psi);
  
       // Initialise value of p
       double interpolated_p = 0.0;
       // Loop over the local nodes and sum
       for (unsigned l = 0; l < n_pres; l++)
       {
         interpolated_p += p_spherical_nst(l) * psi[l];
       }
  
       return (interpolated_p);
     }

References npres_spherical_nst(), p_spherical_nst(), pshape_spherical_nst(), and s.

Referenced by compute_error_e(), extract_velocity(), full_output(), oomph::RefineableQSphericalCrouzeixRaviartElement::further_build(), oomph::RefineableQSphericalTaylorHoodElement::get_interpolated_values(), get_traction(), oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::output(), output(), and pressure_integral().

◆ interpolated_u_spherical_nst() [1/2]

double oomph::SphericalNavierStokesEquations::interpolated_u_spherical_nst	(	const Vector< double > &	s,
		const unsigned &	i
	)		const

inline

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

     {
       // Find number of nodes
       unsigned n_node = nnode();
       // Local shape function
       Shape psi(n_node);
       // Find values of shape function
       shape(s, psi);
  
       // Get nodal index at which i-th velocity is stored
       unsigned u_nodal_index = u_index_spherical_nst(i);
  
       // Initialise value of u
       double interpolated_u = 0.0;
       // Loop over the local nodes and sum
       for (unsigned l = 0; l < n_node; l++)
       {
         interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
       }
  
       return (interpolated_u);
     }

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

◆ interpolated_u_spherical_nst() [2/2]

void oomph::SphericalNavierStokesEquations::interpolated_u_spherical_nst	(	const Vector< double > &	s,
		Vector< double > &	veloc
	)		const

inline

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

     {
       // Find number of nodes
       unsigned n_node = nnode();
       // Local shape function
       Shape psi(n_node);
       // Find values of shape function
       shape(s, psi);
  
       for (unsigned i = 0; i < 3; i++)
       {
         // Index at which the nodal value is stored
         unsigned u_nodal_index = u_index_spherical_nst(i);
         // Initialise value of u
         veloc[i] = 0.0;
         // Loop over the local nodes and sum
         for (unsigned l = 0; l < n_node; l++)
         {
           veloc[i] += nodal_value(l, u_nodal_index) * psi[l];
         }
       }
     }

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

Referenced by compute_error(), compute_error_e(), compute_shear_stress(), extract_velocity(), full_output(), oomph::RefineableQSphericalTaylorHoodElement::get_interpolated_values(), oomph::RefineableQSphericalCrouzeixRaviartElement::get_interpolated_values(), oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::get_wind_spherical_adv_diff(), kin_energy(), oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::output(), and output().

◆ kin_energy()

double oomph::SphericalNavierStokesEquations::kin_energy ( ) const

Get integral of kinetic energy over element.

Get integral of kinetic energy over element:

   {
     // Initialise
     double kin_en = 0.0;
  
     // Set the value of n_intpt
     unsigned n_intpt = integral_pt()->nweight();
  
     // Set the Vector to hold local coordinates
     Vector<double> s(2);
  
     // Loop over the integration points
     for (unsigned ipt = 0; ipt < n_intpt; ipt++)
     {
       // Assign values of s
       for (unsigned i = 0; i < 2; i++)
       {
         s[i] = integral_pt()->knot(ipt, i);
       }
  
       // Get the integral weight
       double w = integral_pt()->weight(ipt);
  
       // Get Jacobian of mapping
       double J = J_eulerian(s);
  
       // Loop over directions
       double veloc_squared = 0.0;
       for (unsigned i = 0; i < 3; i++)
       {
         veloc_squared += interpolated_u_spherical_nst(s, i) *
                          interpolated_u_spherical_nst(s, i);
       }
  
       kin_en += 0.5 * veloc_squared * w * J;
     }
  
     return kin_en;
   }

References i, oomph::FiniteElement::integral_pt(), interpolated_u_spherical_nst(), J, oomph::FiniteElement::J_eulerian(), oomph::Integral::knot(), oomph::Integral::nweight(), s, w, and oomph::Integral::weight().

◆ npres_spherical_nst()

virtual unsigned oomph::SphericalNavierStokesEquations::npres_spherical_nst ( ) const

pure virtual

Function to return number of pressure degrees of freedom.

Implemented in oomph::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), fill_in_generic_residual_contribution_spherical_nst(), and interpolated_p_spherical_nst().

◆ output() [1/4]

void oomph::SphericalNavierStokesEquations::output ( FILE * file_pt )

inlinevirtual

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

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

     {
       unsigned nplot = 5;
       output(file_pt, nplot);
     }

References output().

◆ output() [2/4]

void oomph::SphericalNavierStokesEquations::output	(	FILE *	file_pt,
		const unsigned &	nplot
	)

virtual

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

C-style output function: x,y,[z],u,v,[w],p in tecplot format. Specified number of plot points in each coordinate direction.

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

   {
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Tecplot header info
     fprintf(file_pt, "%s", tecplot_zone_string(nplot).c_str());
  
     // Loop over plot points
     unsigned num_plot_points = nplot_points(nplot);
     for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
     {
       // Get local coordinates of plot point
       get_s_plot(iplot, nplot, s);
  
       // Coordinates
       for (unsigned i = 0; i < 2; i++)
       {
         fprintf(file_pt, "%g ", interpolated_x(s, i));
       }
  
       // Velocities
       for (unsigned i = 0; i < 3; i++)
       {
         fprintf(file_pt, "%g ", interpolated_u_spherical_nst(s, i));
       }
  
       // Pressure
       fprintf(file_pt, "%g \n", interpolated_p_spherical_nst(s));
     }
     fprintf(file_pt, "\n");
  
     // Write tecplot footer (e.g. FE connectivity lists)
     write_tecplot_zone_footer(file_pt, nplot);
   }

References oomph::FiniteElement::get_s_plot(), i, interpolated_p_spherical_nst(), interpolated_u_spherical_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::SphericalNavierStokesEquations::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::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

     {
       unsigned nplot = 5;
       output(outfile, nplot);
     }

Referenced by output(), oomph::QSphericalCrouzeixRaviartElement::output(), and oomph::QSphericalTaylorHoodElement::output().

◆ output() [4/4]

void oomph::SphericalNavierStokesEquations::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: x,y,[z],u,v,[w],p in tecplot format. Specified number of plot points in each coordinate direction.

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

   {
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Tecplot header info
     outfile << tecplot_zone_string(nplot);
  
     // Loop over plot points
     unsigned num_plot_points = nplot_points(nplot);
     for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
     {
       // Get local coordinates of plot point
       get_s_plot(iplot, nplot, s);
  
       // Coordinates in Cartesian Coordinates
       // for(unsigned i=0;i<DIM;i++)
       // {
       //  outfile << interpolated_x(s,i) << " ";
       // }
  
       // Output the coordinates in Cartesian coordintes
       double r = interpolated_x(s, 0);
       double theta = interpolated_x(s, 1);
       // The actual Cartesian values
       outfile << r * sin(theta) << " " << r * cos(theta) << " ";
  
       // Now the x and y velocities
       double u_r = interpolated_u_spherical_nst(s, 0);
       double u_theta = interpolated_u_spherical_nst(s, 1);
  
       outfile << u_r * sin(theta) + u_theta * cos(theta) << " "
               << u_r * cos(theta) - u_theta * sin(theta) << " ";
  
       // Now the polar coordinates
       // outfile << r << " " << theta << " ";
  
       // Velocities r, theta, phi
       for (unsigned i = 0; i < 3; i++)
       {
         outfile << interpolated_u_spherical_nst(s, i) << " ";
       }
  
       // Pressure
       outfile << interpolated_p_spherical_nst(s) << " ";
  
       // Dissipation
       outfile << dissipation(s) << " ";
  
       outfile << std::endl;
     }
  
     // Write tecplot footer (e.g. FE connectivity lists)
     write_tecplot_zone_footer(outfile, nplot);
   }

References cos(), dissipation(), oomph::FiniteElement::get_s_plot(), i, interpolated_p_spherical_nst(), interpolated_u_spherical_nst(), oomph::FiniteElement::interpolated_x(), oomph::FiniteElement::nplot_points(), UniformPSDSelfTest::r, s, sin(), oomph::FiniteElement::tecplot_zone_string(), BiharmonicTestFunctions2::theta, Global_Parameters::u_r(), Global_Parameters::u_theta(), and oomph::FiniteElement::write_tecplot_zone_footer().

◆ output_fct() [1/2]

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

   {
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Vector for coordinates
     Vector<double> x(2);
  
     // Tecplot header info
     outfile << tecplot_zone_string(nplot);
  
     // Exact solution Vector
     Vector<double> exact_soln;
  
     // Loop over plot points
     unsigned num_plot_points = nplot_points(nplot);
     for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
     {
       // Get local coordinates of plot point
       get_s_plot(iplot, nplot, s);
  
       // Get x position as Vector
       interpolated_x(s, x);
  
       // Get exact solution at this point
       (*exact_soln_pt)(time, x, exact_soln);
  
       // Output x,y,...
       for (unsigned i = 0; i < 3; i++)
       {
         outfile << x[i] << " ";
       }
  
       // Output "exact solution"
       for (unsigned i = 0; i < exact_soln.size(); i++)
       {
         outfile << exact_soln[i] << " ";
       }
  
       outfile << std::endl;
     }
  
     // Write tecplot footer (e.g. FE connectivity lists)
     write_tecplot_zone_footer(outfile, nplot);
   }

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

   {
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Vector for coordintes
     Vector<double> x(2);
  
     // Tecplot header info
     outfile << tecplot_zone_string(nplot);
  
     // Exact solution Vector
     Vector<double> exact_soln;
  
  
     // Loop over plot points
     unsigned num_plot_points = nplot_points(nplot);
     for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
     {
       // Get local coordinates of plot point
       get_s_plot(iplot, nplot, s);
  
       // Get x position as Vector
       interpolated_x(s, x);
  
       // Get exact solution at this point
       (*exact_soln_pt)(x, exact_soln);
  
       // Output x,y,...
       for (unsigned i = 0; i < 3; i++)
       {
         outfile << x[i] << " ";
       }
  
       // Output "exact solution"
       for (unsigned i = 0; i < exact_soln.size(); i++)
       {
         outfile << exact_soln[i] << " ";
       }
  
       outfile << std::endl;
     }
  
     // Write tecplot footer (e.g. FE connectivity lists)
     write_tecplot_zone_footer(outfile, nplot);
   }

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

   {
     // Find number of nodes
     unsigned n_node = nnode();
  
     // Local shape function
     Shape psi(n_node);
  
     // Vectors of local coordinates and coords and velocities
     Vector<double> s(2);
     Vector<double> interpolated_x(2);
     Vector<double> interpolated_u(3);
  
     // Tecplot header info
     outfile << tecplot_zone_string(nplot);
  
     // Loop over plot points
     unsigned num_plot_points = nplot_points(nplot);
     for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
     {
       // Get local coordinates of plot point
       get_s_plot(iplot, nplot, s);
  
       // Get shape functions
       shape(s, psi);
  
       // Loop over directions
       for (unsigned i = 0; i < 3; i++)
       {
         interpolated_x[i] = 0.0;
         interpolated_u[i] = 0.0;
         // Get the index at which velocity i is stored
         unsigned u_nodal_index = u_index_spherical_nst(i);
         // Loop over the local nodes and sum
         for (unsigned l = 0; l < n_node; l++)
         {
           interpolated_u[i] += nodal_value(t, l, u_nodal_index) * psi[l];
           interpolated_x[i] += nodal_position(t, l, i) * psi[l];
         }
       }
  
       // Coordinates
       for (unsigned i = 0; i < 2; i++)
       {
         outfile << interpolated_x[i] << " ";
       }
  
       // Velocities
       for (unsigned i = 0; i < 3; i++)
       {
         outfile << interpolated_u[i] << " ";
       }
  
       outfile << std::endl;
     }
  
     // Write tecplot footer (e.g. FE connectivity lists)
     write_tecplot_zone_footer(outfile, nplot);
   }

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_spherical_nst(), and oomph::FiniteElement::write_tecplot_zone_footer().

◆ output_vorticity()

void oomph::SphericalNavierStokesEquations::output_vorticity	(	std::ostream &	outfile,
		const unsigned &	nplot
	)

Output function: x,y,[z], [omega_x,omega_y,[and/or omega_z]] in tecplot format. nplot points in each coordinate direction

Output function for vorticity. x,y,[z],[omega_x,omega_y,[and/or omega_z]] in tecplot format. Specified number of plot points in each coordinate direction.

   {
     // Vector of local coordinates
     Vector<double> s(2);
  
     // Create vorticity vector of the required size
     Vector<double> vorticity(1);
  
     // Tecplot header info
     outfile << tecplot_zone_string(nplot);
  
     // Loop over plot points
     unsigned num_plot_points = nplot_points(nplot);
     for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
     {
       // Get local coordinates of plot point
       get_s_plot(iplot, nplot, s);
  
       // Coordinates
       for (unsigned i = 0; i < 2; i++)
       {
         outfile << interpolated_x(s, i) << " ";
       }
  
       // Get vorticity
       get_vorticity(s, vorticity);
  
       outfile << vorticity[0] << std::endl;
       ;
     }
  
     // Write tecplot footer (e.g. FE connectivity lists)
     write_tecplot_zone_footer(outfile, nplot);
   }

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

◆ p_local_eqn()

virtual int oomph::SphericalNavierStokesEquations::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::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), and fill_in_generic_residual_contribution_spherical_nst().

◆ p_nodal_index_spherical_nst()

virtual int oomph::SphericalNavierStokesEquations::p_nodal_index_spherical_nst ( ) const

inlinevirtual

Return the index at which the pressure is stored if it is stored at the nodes. If not stored at the nodes this will return a negative number.

Reimplemented in oomph::QSphericalTaylorHoodElement.

     {
       return Pressure_not_stored_at_node;
     }

References Pressure_not_stored_at_node.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst().

◆ p_spherical_nst()

virtual double oomph::SphericalNavierStokesEquations::p_spherical_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::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), fill_in_generic_residual_contribution_spherical_nst(), and interpolated_p_spherical_nst().

◆ pressure_integral()

double oomph::SphericalNavierStokesEquations::pressure_integral ( ) const

Integral of pressure over element.

Return pressure integrated over the element.

   {
     // Initialise
     double press_int = 0;
  
     // Set the value of n_intpt
     unsigned n_intpt = integral_pt()->nweight();
  
     // Set the Vector to hold local coordinates
     Vector<double> s(2);
  
     // Loop over the integration points
     for (unsigned ipt = 0; ipt < n_intpt; ipt++)
     {
       // Assign values of s
       for (unsigned i = 0; i < 2; i++)
       {
         s[i] = integral_pt()->knot(ipt, i);
       }
  
       // Get the integral weight
       double w = integral_pt()->weight(ipt);
  
       // Get Jacobian of mapping
       double J = J_eulerian(s);
  
       // Premultiply the weights and the Jacobian
       double W = w * J;
  
       // Get pressure
       double press = interpolated_p_spherical_nst(s);
  
       // Add
       press_int += press * W;
     }
  
     return press_int;
   }

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

◆ pshape_spherical_nst() [1/2]

virtual void oomph::SphericalNavierStokesEquations::pshape_spherical_nst	(	const Vector< double > &	s,
		Shape &	psi
	)		const

protectedpure virtual

Compute the pressure shape functions at local coordinate s.

Implemented in oomph::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), fill_in_generic_residual_contribution_spherical_nst(), and interpolated_p_spherical_nst().

◆ pshape_spherical_nst() [2/2]

virtual void oomph::SphericalNavierStokesEquations::pshape_spherical_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::QSphericalTaylorHoodElement, and oomph::QSphericalCrouzeixRaviartElement.

◆ re()

const double& oomph::SphericalNavierStokesEquations::re ( ) const

inline

Reynolds number.

     {
       return *Re_pt;
     }

References Re_pt.

Referenced by actual(), actual_dr(), actual_dth(), oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), and fill_in_generic_residual_contribution_spherical_nst().

◆ re_invfr()

const double& oomph::SphericalNavierStokesEquations::re_invfr ( ) const

inline

Global inverse Froude number.

     {
       return *ReInvFr_pt;
     }

References ReInvFr_pt.

◆ re_invfr_pt()

double*& oomph::SphericalNavierStokesEquations::re_invfr_pt ( )

inline

Pointer to global inverse Froude number.

     {
       return ReInvFr_pt;
     }

References ReInvFr_pt.

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

◆ re_invro()

const double& oomph::SphericalNavierStokesEquations::re_invro ( ) const

inline

Global Reynolds number multiplied by inverse Rossby number.

     {
       return *ReInvRo_pt;
     }

References ReInvRo_pt.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), and fill_in_generic_residual_contribution_spherical_nst().

◆ re_invro_pt()

double*& oomph::SphericalNavierStokesEquations::re_invro_pt ( )

inline

Pointer to global inverse Froude number.

     {
       return ReInvRo_pt;
     }

References ReInvRo_pt.

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

◆ re_pt()

double*& oomph::SphericalNavierStokesEquations::re_pt ( )

inline

Pointer to Reynolds number.

     {
       return Re_pt;
     }

References Re_pt.

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

◆ re_st()

const double& oomph::SphericalNavierStokesEquations::re_st ( ) const

inline

Product of Reynolds and Strouhal number (=Womersley number)

     {
       return *ReSt_pt;
     }

References ReSt_pt.

Referenced by oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), and fill_in_generic_residual_contribution_spherical_nst().

◆ re_st_pt()

double*& oomph::SphericalNavierStokesEquations::re_st_pt ( )

inline

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

     {
       return ReSt_pt;
     }

References ReSt_pt.

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

◆ source_fct_pt() [1/2]

SphericalNavierStokesSourceFctPt& oomph::SphericalNavierStokesEquations::source_fct_pt ( )

inline

Access function for the source-function pointer.

     {
       return Source_fct_pt;
     }

References Source_fct_pt.

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

◆ source_fct_pt() [2/2]

SphericalNavierStokesSourceFctPt oomph::SphericalNavierStokesEquations::source_fct_pt ( ) const

inline

Access function for the source-function pointer. Const version.

     {
       return Source_fct_pt;
     }

References Source_fct_pt.

◆ strain_rate()

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

Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)

Get strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)

   {
     // Velocity components
     Vector<double> interpolated_u(3, 0.0);
     // Coordinates
     double interpolated_r = 0.0;
     double interpolated_theta = 0.0;
     // Velocity gradient matrix
     DenseMatrix<double> dudx(3, 2, 0.0);
  
  
     // Find out how many nodes there are in the element
     unsigned n_node = nnode();
  
     // Set up memory for the shape and test functions
     Shape psi(n_node);
     DShape dpsidx(n_node, 2);
  
     // Call the derivatives of the shape functions
     dshape_eulerian(s, psi, dpsidx);
  
     // Get the values of the positions
     for (unsigned l = 0; l < n_node; l++)
     {
       interpolated_r += this->nodal_position(l, 0) * psi(l);
       interpolated_theta += this->nodal_position(l, 1) * psi(l);
     }
  
  
     // Loop over velocity components
     for (unsigned i = 0; i < 3; i++)
     {
       // Get the index at which the i-th velocity is stored
       unsigned u_nodal_index = u_index_spherical_nst(i);
       // Loop over nodes
       for (unsigned l = 0; l < n_node; l++)
       {
         // Get the velocity value
         const double nodal_u = this->nodal_value(l, u_nodal_index);
         // Add in the velocities
         interpolated_u[i] += nodal_u * psi(l);
         // Loop over the derivative directions
         for (unsigned j = 0; j < 2; j++)
         {
           dudx(i, j) += nodal_u * dpsidx(l, j);
         }
       }
     }
  
     // Assign those strain rates without any negative powers of the radius
     strainrate(0, 0) = dudx(0, 0);
     strainrate(0, 1) = 0.0;
     strainrate(0, 2) = 0.0;
     ;
     strainrate(1, 1) = 0.0;
     strainrate(1, 2) = 0.0;
     strainrate(2, 2) = 0.0;
  
  
     // Add the negative powers of the radius, unless we are
     // at the origin
     if (std::fabs(interpolated_r) > 1.0e-15)
     {
       const double pi = MathematicalConstants::Pi;
       double inverse_r = 1.0 / interpolated_r;
  
       // Should we include the cot terms (default no)
       bool include_cot_terms = false;
       double cot_theta = 0.0;
       // If we in the legal range then include cot
       if ((std::fabs(interpolated_theta) > 1.0e-15) &&
           (std::fabs(pi - interpolated_theta) > 1.0e-15))
       {
         include_cot_terms = true;
         cot_theta = this->cot(interpolated_theta);
       }
  
       strainrate(0, 1) =
         0.5 * (dudx(1, 0) + (dudx(0, 1) - interpolated_u[1]) * inverse_r);
       strainrate(0, 2) = 0.5 * (dudx(2, 0) - interpolated_u[2] * inverse_r);
       strainrate(1, 1) = (dudx(1, 1) + interpolated_u[0]) * inverse_r;
  
       if (include_cot_terms)
       {
         strainrate(1, 2) =
           0.5 * (dudx(2, 1) - interpolated_u[2] * cot_theta) * inverse_r;
         strainrate(2, 2) =
           (interpolated_u[0] + interpolated_u[1] * cot_theta) * inverse_r;
       }
     }
  
     // Sort out the symmetries
     strainrate(1, 0) = strainrate(0, 1);
     strainrate(2, 0) = strainrate(0, 2);
     strainrate(2, 1) = strainrate(1, 2);
   }

References cot(), oomph::FiniteElement::dshape_eulerian(), boost::multiprecision::fabs(), i, j, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), constants::pi, oomph::MathematicalConstants::Pi, s, and u_index_spherical_nst().

Referenced by dissipation(), get_traction(), and oomph::RefineableSphericalNavierStokesEquations::get_Z2_flux().

◆ u_index_spherical_nst()

virtual unsigned oomph::SphericalNavierStokesEquations::u_index_spherical_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.

     {
       return i;
     }

References i.

Referenced by d_kin_energy_dt(), du_dt_spherical_nst(), oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), fill_in_generic_residual_contribution_spherical_nst(), oomph::RefineableBuoyantQSphericalCrouzeixRaviartElement::fill_in_off_diagonal_jacobian_blocks_analytic(), full_output(), oomph::RefineableQSphericalTaylorHoodElement::get_interpolated_values(), oomph::RefineableQSphericalCrouzeixRaviartElement::get_interpolated_values(), get_pressure_and_velocity_mass_matrix_diagonal(), get_vorticity(), oomph::QSphericalCrouzeixRaviartElement::identify_load_data(), oomph::QSphericalTaylorHoodElement::identify_load_data(), interpolated_dudx_spherical_nst(), interpolated_u_spherical_nst(), output_veloc(), strain_rate(), and u_spherical_nst().

◆ u_spherical_nst() [1/2]

double oomph::SphericalNavierStokesEquations::u_spherical_nst	(	const unsigned &	n,
		const unsigned &	i
	)		const

inline

Velocity i at local node n. Uses suitably interpolated value for hanging nodes. The use of u_index_spherical_nst() permits the use of this element as the basis for multi-physics elements. The default is to assume that the i-th velocity component is stored at the i-th location of the node

     {
       return nodal_value(n, u_index_spherical_nst(i));
     }

References i, n, oomph::FiniteElement::nodal_value(), and u_index_spherical_nst().

◆ u_spherical_nst() [2/2]

double oomph::SphericalNavierStokesEquations::u_spherical_nst	(	const unsigned &	t,
		const unsigned &	n,
		const unsigned &	i
	)		const

inline

Velocity i at local node n at timestep t (t=0: present; t>0: previous). Uses suitably interpolated value for hanging nodes.

     {
       return nodal_value(t, n, u_index_spherical_nst(i));
     }

References i, n, oomph::FiniteElement::nodal_value(), plotPSD::t, and u_index_spherical_nst().

◆ viscosity_ratio()

const double& oomph::SphericalNavierStokesEquations::viscosity_ratio ( ) const

inline

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

     {
       return *Viscosity_Ratio_pt;
     }

References Viscosity_Ratio_pt.

◆ viscosity_ratio_pt()

double*& oomph::SphericalNavierStokesEquations::viscosity_ratio_pt ( )

inline

Pointer to Viscosity Ratio.

     {
       return Viscosity_Ratio_pt;
     }

References Viscosity_Ratio_pt.

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

Member Data Documentation

◆ ALE_is_disabled

bool oomph::SphericalNavierStokesEquations::ALE_is_disabled

protected

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

Referenced by d_kin_energy_dt(), disable_ALE(), enable_ALE(), oomph::RefineableSphericalNavierStokesEquations::fill_in_generic_residual_contribution_spherical_nst(), fill_in_generic_residual_contribution_spherical_nst(), and oomph::RefineableSphericalNavierStokesEquations::further_build().

◆ Body_force_fct_pt

SphericalNavierStokesBodyForceFctPt oomph::SphericalNavierStokesEquations::Body_force_fct_pt

protected

Pointer to body force function.

Referenced by body_force_fct_pt(), oomph::RefineableSphericalNavierStokesEquations::further_build(), and get_body_force_spherical_nst().

◆ Default_Gravity_vector

Vector< double > oomph::SphericalNavierStokesEquations::Default_Gravity_vector

staticprivate

Static default value for the gravity vector.

Navier-Stokes equations default gravity vector.

Referenced by SphericalNavierStokesEquations().

◆ Default_Physical_Constant_Value

double oomph::SphericalNavierStokesEquations::Default_Physical_Constant_Value = 0.0

staticprivate

Navier–Stokes equations static data.

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

Referenced by SphericalNavierStokesEquations().

◆ Default_Physical_Ratio_Value

double oomph::SphericalNavierStokesEquations::Default_Physical_Ratio_Value = 1.0

staticprivate

Navier–Stokes equations static data.

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

Referenced by SphericalNavierStokesEquations().

◆ Density_Ratio_pt

double* oomph::SphericalNavierStokesEquations::Density_Ratio_pt

protected

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

Referenced by density_ratio(), density_ratio_pt(), oomph::RefineableSphericalNavierStokesEquations::further_build(), and SphericalNavierStokesEquations().

◆ G_pt

Vector<double>* oomph::SphericalNavierStokesEquations::G_pt

protected

Pointer to global gravity Vector.

Referenced by oomph::RefineableSphericalNavierStokesEquations::further_build(), g(), g_pt(), and SphericalNavierStokesEquations().

◆ Gamma

Vector< double > oomph::SphericalNavierStokesEquations::Gamma

static

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

Navier–Stokes equations static data.

◆ Pressure_not_stored_at_node

int oomph::SphericalNavierStokesEquations::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_spherical_nst().

◆ Re_pt

double* oomph::SphericalNavierStokesEquations::Re_pt

protected

Pointer to global Reynolds number.

Referenced by oomph::RefineableSphericalNavierStokesEquations::further_build(), re(), re_pt(), and SphericalNavierStokesEquations().

◆ ReInvFr_pt

double* oomph::SphericalNavierStokesEquations::ReInvFr_pt

protected

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

Referenced by oomph::RefineableSphericalNavierStokesEquations::further_build(), re_invfr(), re_invfr_pt(), and SphericalNavierStokesEquations().

◆ ReInvRo_pt

double* oomph::SphericalNavierStokesEquations::ReInvRo_pt

protected

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

Referenced by oomph::RefineableSphericalNavierStokesEquations::further_build(), re_invro(), re_invro_pt(), and SphericalNavierStokesEquations().

◆ ReSt_pt

double* oomph::SphericalNavierStokesEquations::ReSt_pt

protected

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

Referenced by oomph::RefineableSphericalNavierStokesEquations::further_build(), re_st(), re_st_pt(), and SphericalNavierStokesEquations().

◆ Source_fct_pt

SphericalNavierStokesSourceFctPt oomph::SphericalNavierStokesEquations::Source_fct_pt

protected

Pointer to volumetric source function.

Referenced by oomph::RefineableSphericalNavierStokesEquations::further_build(), get_source_spherical_nst(), and source_fct_pt().

◆ Viscosity_Ratio_pt

double* oomph::SphericalNavierStokesEquations::Viscosity_Ratio_pt

protected

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

Referenced by oomph::RefineableSphericalNavierStokesEquations::further_build(), SphericalNavierStokesEquations(), viscosity_ratio(), and viscosity_ratio_pt().

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

Public Types

Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Static Private Attributes

Additional Inherited Members

Detailed Description

Member Typedef Documentation

◆ SphericalNavierStokesBodyForceFctPt

◆ SphericalNavierStokesSourceFctPt

Constructor & Destructor Documentation

◆ SphericalNavierStokesEquations()

Member Function Documentation

◆ actual()

◆ actual_dr()

◆ actual_dth()

◆ body_force_fct_pt() [1/2]

◆ body_force_fct_pt() [2/2]

◆ compute_error() [1/2]

◆ compute_error() [2/2]

◆ compute_error_e()

◆ compute_shear_stress()

◆ cot()

◆ d_kin_energy_dt()

◆ density_ratio()

◆ density_ratio_pt()

◆ disable_ALE()

◆ dissipation() [1/2]

◆ dissipation() [2/2]

◆ dshape_and_dtest_eulerian_at_knot_spherical_nst()

◆ dshape_and_dtest_eulerian_spherical_nst()

◆ du_dt_spherical_nst()

◆ enable_ALE()

◆ extract_velocity()

◆ fill_in_contribution_to_jacobian()

◆ fill_in_contribution_to_jacobian_and_mass_matrix()

◆ fill_in_contribution_to_residuals()

◆ fill_in_generic_residual_contribution_spherical_nst()

◆ fix_pressure()

◆ full_output() [1/2]

◆ full_output() [2/2]

◆ g()

◆ g_pt()

◆ get_body_force_spherical_nst()

◆ get_load()

◆ get_pressure_and_velocity_mass_matrix_diagonal()

◆ get_source_spherical_nst()

◆ get_traction()

◆ get_vorticity()

◆ interpolated_dudx_spherical_nst()

◆ interpolated_p_spherical_nst()

◆ interpolated_u_spherical_nst() [1/2]

◆ interpolated_u_spherical_nst() [2/2]

◆ kin_energy()

◆ npres_spherical_nst()

◆ output() [1/4]

◆ output() [2/4]

◆ output() [3/4]

◆ output() [4/4]

◆ output_fct() [1/2]

◆ output_fct() [2/2]

◆ output_veloc()

◆ output_vorticity()

◆ p_local_eqn()

◆ p_nodal_index_spherical_nst()

◆ p_spherical_nst()

◆ pressure_integral()

◆ pshape_spherical_nst() [1/2]

◆ pshape_spherical_nst() [2/2]

◆ re()

◆ re_invfr()

◆ re_invfr_pt()

◆ re_invro()

◆ re_invro_pt()

◆ re_pt()

◆ re_st()

◆ re_st_pt()

◆ source_fct_pt() [1/2]

◆ source_fct_pt() [2/2]