oomph::RefineableSphericalNavierStokesEquations Class Referenceabstract

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

#include <refineable_spherical_navier_stokes_elements.h>

+ Inheritance diagram for oomph::RefineableSphericalNavierStokesEquations:

Public Member Functions

 RefineableSphericalNavierStokesEquations ()
 Empty Constructor. More...
 
unsigned num_Z2_flux_terms ()
 Number of 'flux' terms for Z2 error estimation. More...
 
void get_Z2_flux (const Vector< double > &s, Vector< double > &flux)
 
double geometric_jacobian (const Vector< double > &x)
 Fill in the geometric Jacobian, which in this case is r*r*sin(theta) More...
 
void further_build ()
 Further build: pass the pointers down to the sons. More...
 
- Public Member Functions inherited from oomph::SphericalNavierStokesEquations
double cot (const double &th) const
 
 SphericalNavierStokesEquations ()
 
const doublere () const
 Reynolds number. More...
 
const doublere_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
const doubleviscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to Viscosity Ratio. More...
 
const doubledensity_ratio () const
 
double *& density_ratio_pt ()
 Pointer to Density ratio. More...
 
const doublere_invfr () const
 Global inverse Froude number. More...
 
double *& re_invfr_pt ()
 Pointer to global inverse Froude number. More...
 
const doublere_invro () const
 Global Reynolds number multiplied by inverse Rossby number. More...
 
double *& re_invro_pt ()
 Pointer to global inverse Froude number. More...
 
const Vector< double > & g () const
 Vector of gravitational components. More...
 
Vector< double > *& g_pt ()
 Pointer to Vector of gravitational components. More...
 
SphericalNavierStokesBodyForceFctPtbody_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...
 
SphericalNavierStokesSourceFctPtsource_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< doubleactual (const Vector< double > &x)
 
Vector< doubleactual_dr (const Vector< double > &x)
 
Vector< doubleactual_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)
 
MacroElementmacro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void describe_nodal_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void assign_all_generic_local_eqn_numbers (const bool &store_local_dof_pt)
 
Node *& node_pt (const unsigned &n)
 Return a pointer to the local node n. More...
 
Node *const & node_pt (const unsigned &n) const
 Return a pointer to the local node n (const version) More...
 
unsigned nnode () const
 Return the number of nodes. More...
 
virtual unsigned nnode_1d () const
 
double raw_nodal_position (const unsigned &n, const unsigned &i) const
 
double raw_nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position (const unsigned &n, const unsigned &i) const
 
double nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double dnodal_position_dt (const unsigned &n, const unsigned &i) const
 Return the i-th component of nodal velocity: dx/dt at local node n. More...
 
double dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual unsigned required_nvalue (const unsigned &n) const
 
unsigned nnodal_position_type () const
 
bool has_hanging_nodes () const
 
unsigned nodal_dimension () const
 Return the required Eulerian dimension of the nodes in this element. More...
 
virtual Nodeconstruct_node (const unsigned &n)
 
virtual Nodeconstruct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual Nodeconstruct_boundary_node (const unsigned &n)
 
virtual Nodeconstruct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
int get_node_number (Node *const &node_pt) const
 
virtual Nodeget_node_at_local_coordinate (const Vector< double > &s) const
 
double raw_nodal_value (const unsigned &n, const unsigned &i) const
 
double raw_nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
unsigned dim () const
 
virtual ElementGeometry::ElementGeometry element_geometry () const
 Return the geometry type of the element (either Q or T usually). More...
 
virtual double interpolated_x (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated coordinate x[i] at local coordinate s. More...
 
virtual double interpolated_x (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_x (const Vector< double > &s, Vector< double > &x) const
 Return FE interpolated position x[] at local coordinate s as Vector. More...
 
virtual void interpolated_x (const unsigned &t, const Vector< double > &s, Vector< double > &x) const
 
virtual double interpolated_dxdt (const Vector< double > &s, const unsigned &i, const unsigned &t)
 
virtual void interpolated_dxdt (const Vector< double > &s, const unsigned &t, Vector< double > &dxdt)
 
unsigned ngeom_data () const
 
Datageom_data_pt (const unsigned &j)
 
void position (const Vector< double > &zeta, Vector< double > &r) const
 
void position (const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
 
void dposition_dt (const Vector< double > &zeta, const unsigned &t, Vector< double > &drdt)
 
virtual double zeta_nodal (const unsigned &n, const unsigned &k, const unsigned &i) const
 
void interpolated_zeta (const Vector< double > &s, Vector< double > &zeta) const
 
void locate_zeta (const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
 
virtual void node_update ()
 
virtual void identify_geometric_data (std::set< Data * > &geometric_data_pt)
 
virtual double s_min () const
 Min value of local coordinate. More...
 
virtual double s_max () const
 Max. value of local coordinate. More...
 
double size () const
 
virtual double compute_physical_size () const
 
virtual void point_output_data (const Vector< double > &s, Vector< double > &data)
 
void point_output (std::ostream &outfile, const Vector< double > &s)
 
virtual unsigned nplot_points_paraview (const unsigned &nplot) const
 
virtual unsigned nsub_elements_paraview (const unsigned &nplot) const
 
void output_paraview (std::ofstream &file_out, const unsigned &nplot) const
 
virtual void write_paraview_output_offset_information (std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
 
virtual void write_paraview_type (std::ofstream &file_out, const unsigned &nplot) const
 
virtual void write_paraview_offsets (std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
 
virtual unsigned nscalar_paraview () const
 
virtual void scalar_value_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
 
virtual void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
 
virtual void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
virtual std::string scalar_name_paraview (const unsigned &i) const
 
virtual void output (const unsigned &t, std::ostream &outfile, const unsigned &n_plot) const
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, const SolutionFunctorBase &exact_soln) const
 Output a time-dependent exact solution over the element. More...
 
virtual void get_s_plot (const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
 
virtual std::string tecplot_zone_string (const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (std::ostream &outfile, const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (FILE *file_pt, const unsigned &nplot) const
 
virtual unsigned nplot_points (const unsigned &nplot) const
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_abs_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error)
 
void integrate_fct (FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt, Vector< double > &integral)
 Integrate Vector-valued function over element. More...
 
void integrate_fct (FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt, const double &time, Vector< double > &integral)
 Integrate Vector-valued time-dep function over element. More...
 
virtual void build_face_element (const int &face_index, FaceElement *face_element_pt)
 
virtual unsigned self_test ()
 
virtual unsigned get_bulk_node_number (const int &face_index, const unsigned &i) const
 
virtual int face_outer_unit_normal_sign (const int &face_index) const
 Get the sign of the outer unit normal on the face given by face_index. More...
 
virtual unsigned nnode_on_face () const
 
void face_node_number_error_check (const unsigned &i) const
 Range check for face node numbers. More...
 
virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt (const int &face_index) const
 
virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt (const int &face_index) const
 
- Public Member Functions inherited from oomph::GeneralisedElement
 GeneralisedElement ()
 Constructor: Initialise all pointers and all values to zero. More...
 
virtual ~GeneralisedElement ()
 Virtual destructor to clean up any memory allocated by the object. More...
 
 GeneralisedElement (const GeneralisedElement &)=delete
 Broken copy constructor. More...
 
void operator= (const GeneralisedElement &)=delete
 Broken assignment operator. More...
 
Data *& internal_data_pt (const unsigned &i)
 Return a pointer to i-th internal data object. More...
 
Data *const & internal_data_pt (const unsigned &i) const
 Return a pointer to i-th internal data object (const version) More...
 
Data *& external_data_pt (const unsigned &i)
 Return a pointer to i-th external data object. More...
 
Data *const & external_data_pt (const unsigned &i) const
 Return a pointer to i-th external data object (const version) More...
 
unsigned long eqn_number (const unsigned &ieqn_local) const
 
int local_eqn_number (const unsigned long &ieqn_global) const
 
unsigned add_external_data (Data *const &data_pt, const bool &fd=true)
 
bool external_data_fd (const unsigned &i) const
 
void exclude_external_data_fd (const unsigned &i)
 
void include_external_data_fd (const unsigned &i)
 
void flush_external_data ()
 Flush all external data. More...
 
void flush_external_data (Data *const &data_pt)
 Flush the object addressed by data_pt from the external data array. More...
 
unsigned ninternal_data () const
 Return the number of internal data objects. More...
 
unsigned nexternal_data () const
 Return the number of external data objects. More...
 
unsigned ndof () const
 Return the number of equations/dofs in the element. More...
 
void dof_vector (const unsigned &t, Vector< double > &dof)
 Return the vector of dof values at time level t. More...
 
void dof_pt_vector (Vector< double * > &dof_pt)
 Return the vector of pointers to dof values. More...
 
void set_internal_data_time_stepper (const unsigned &i, TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void assign_internal_eqn_numbers (unsigned long &global_number, Vector< double * > &Dof_pt)
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void add_internal_value_pt_to_map (std::map< unsigned, double * > &map_of_value_pt)
 
virtual void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void complete_setup_of_dependencies ()
 
virtual void get_residuals (Vector< double > &residuals)
 
virtual void get_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
virtual void get_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void get_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void get_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void get_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void get_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void get_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void get_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void get_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_norm (double &norm)
 
virtual unsigned ndof_types () const
 
virtual void get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
 
- Public Member Functions inherited from oomph::GeomObject
 GeomObject ()
 Default constructor. More...
 
 GeomObject (const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)
 
 GeomObject (const GeomObject &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const GeomObject &)=delete
 Broken assignment operator. More...
 
virtual ~GeomObject ()
 (Empty) destructor More...
 
unsigned nlagrangian () const
 Access function to # of Lagrangian coordinates. More...
 
unsigned ndim () const
 Access function to # of Eulerian coordinates. More...
 
void set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
 Set # of Lagrangian and Eulerian coordinates. More...
 
TimeStepper *& time_stepper_pt ()
 
TimeSteppertime_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 
- Public Member Functions inherited from oomph::NavierStokesElementWithDiagonalMassMatrices
 NavierStokesElementWithDiagonalMassMatrices ()
 Empty constructor. More...
 
virtual ~NavierStokesElementWithDiagonalMassMatrices ()
 Virtual destructor. More...
 
 NavierStokesElementWithDiagonalMassMatrices (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken copy constructor. More...
 
void operator= (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken assignment operator. More...
 
- Public Member Functions inherited from oomph::RefineableElement
 RefineableElement ()
 
virtual ~RefineableElement ()
 Destructor, delete the allocated storage for the hanging equations. More...
 
 RefineableElement (const RefineableElement &)=delete
 Broken copy constructor. More...
 
void operator= (const RefineableElement &)=delete
 Broken assignment operator. More...
 
Treetree_pt ()
 Access function: Pointer to quadtree representation of this element. More...
 
void set_tree_pt (Tree *my_tree_pt)
 Set pointer to quadtree representation of this element. More...
 
virtual unsigned required_nsons () const
 
bool refinement_is_enabled ()
 Flag to indicate suppression of any refinement. More...
 
void disable_refinement ()
 Suppress of any refinement for this element. More...
 
void enable_refinement ()
 Emnable refinement for this element. More...
 
template<class ELEMENT >
void split (Vector< ELEMENT * > &son_pt) const
 
int local_hang_eqn (Node *const &node_pt, const unsigned &i)
 
virtual void build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)=0
 
void set_refinement_level (const int &refine_level)
 Set the refinement level. More...
 
unsigned refinement_level () const
 Return the Refinement level. More...
 
void select_for_refinement ()
 Select the element for refinement. More...
 
void deselect_for_refinement ()
 Deselect the element for refinement. More...
 
void select_sons_for_unrefinement ()
 Unrefinement will be performed by merging the four sons of this element. More...
 
void deselect_sons_for_unrefinement ()
 
bool to_be_refined ()
 Has the element been selected for refinement? More...
 
bool sons_to_be_unrefined ()
 Has the element been selected for unrefinement? More...
 
virtual void rebuild_from_sons (Mesh *&mesh_pt)=0
 
virtual void unbuild ()
 
virtual void deactivate_element ()
 
virtual bool nodes_built ()
 Return true if all the nodes have been built, false if not. More...
 
long number () const
 Element number (for debugging/plotting) More...
 
void set_number (const long &mynumber)
 Set element number (for debugging/plotting) More...
 
virtual unsigned ncont_interpolated_values () const =0
 
virtual void get_interpolated_values (const Vector< double > &s, Vector< double > &values)
 
virtual void get_interpolated_values (const unsigned &t, const Vector< double > &s, Vector< double > &values)=0
 
virtual Nodeinterpolating_node_pt (const unsigned &n, const int &value_id)
 
virtual double local_one_d_fraction_of_interpolating_node (const unsigned &n1d, const unsigned &i, const int &value_id)
 
virtual Nodeget_interpolating_node_at_local_coordinate (const Vector< double > &s, const int &value_id)
 
virtual unsigned ninterpolating_node (const int &value_id)
 
virtual unsigned ninterpolating_node_1d (const int &value_id)
 
virtual void interpolating_basis (const Vector< double > &s, Shape &psi, const int &value_id) const
 
virtual void check_integrity (double &max_error)=0
 
void identify_field_data_for_interactions (std::set< std::pair< Data *, unsigned >> &paired_field_data)
 
void assign_nodal_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual RefineableElementroot_element_pt ()
 
virtual RefineableElementfather_element_pt () const
 Return a pointer to the father element. More...
 
void get_father_at_refinement_level (unsigned &refinement_level, RefineableElement *&father_at_reflevel_pt)
 
virtual void initial_setup (Tree *const &adopted_father_pt=0, const unsigned &initial_p_order=0)
 
virtual void pre_build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt)
 Pre-build the element. More...
 
virtual void setup_hanging_nodes (Vector< std::ofstream * > &output_stream)
 
virtual void further_setup_hanging_nodes ()
 
void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
unsigned nshape_controlling_nodes ()
 
std::map< Node *, unsignedshape_controlling_node_lookup ()
 
- Public Member Functions inherited from oomph::ElementWithZ2ErrorEstimator
 ElementWithZ2ErrorEstimator ()
 Default empty constructor. More...
 
 ElementWithZ2ErrorEstimator (const ElementWithZ2ErrorEstimator &)=delete
 Broken copy constructor. More...
 
void operator= (const ElementWithZ2ErrorEstimator &)=delete
 Broken assignment operator. More...
 
virtual unsigned ncompound_fluxes ()
 
virtual void compute_exact_Z2_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
 
virtual void get_Z2_compound_flux_indices (Vector< unsigned > &flux_index)
 
virtual unsigned nvertex_node () const =0
 Number of vertex nodes in the element. More...
 
virtual Nodevertex_node_pt (const unsigned &j) const =0
 
virtual unsigned nrecovery_order ()=0
 Order of recovery shape functions. More...
 

Static Public Member Functions

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

Protected Member Functions

virtual Nodepressure_node_pt (const unsigned &n_p)
 
virtual void unpin_elemental_pressure_dofs ()=0
 Unpin all pressure dofs in the element. More...
 
virtual void pin_elemental_redundant_nodal_pressure_dofs ()
 
- Protected Member Functions inherited from oomph::SphericalNavierStokesEquations
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)
 
- Protected Member Functions inherited from oomph::FiniteElement
template<unsigned DIM>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double invert_jacobian_mapping (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &inverse_jacobian) const
 
virtual void dJ_eulerian_dnodal_coordinates (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<unsigned DIM>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
virtual void d_dshape_eulerian_dnodal_coordinates (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<unsigned DIM>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
virtual void transform_derivatives (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
void transform_derivatives_diagonal (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
virtual void transform_second_derivatives (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
void fill_in_jacobian_from_nodal_by_fd (DenseMatrix< double > &jacobian)
 
virtual void update_before_nodal_fd ()
 
virtual void reset_after_nodal_fd ()
 
virtual void update_in_nodal_fd (const unsigned &i)
 
virtual void reset_in_nodal_fd (const unsigned &i)
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Zero-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 One-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Two-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
- Protected Member Functions inherited from oomph::GeneralisedElement
unsigned add_internal_data (Data *const &data_pt, const bool &fd=true)
 
bool internal_data_fd (const unsigned &i) const
 
void exclude_internal_data_fd (const unsigned &i)
 
void include_internal_data_fd (const unsigned &i)
 
void clear_global_eqn_numbers ()
 
void add_global_eqn_numbers (std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)
 
virtual void assign_internal_and_external_local_eqn_numbers (const bool &store_local_dof_pt)
 
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 Member Functions inherited from oomph::RefineableElement
void assemble_local_to_eulerian_jacobian (const DShape &dpsids, DenseMatrix< double > &jacobian) const
 
void assemble_local_to_eulerian_jacobian2 (const DShape &d2psids, DenseMatrix< double > &jacobian2) const
 
void assemble_eulerian_base_vectors (const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
 
double local_to_eulerian_mapping_diagonal (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
void assign_hanging_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign the local equation numbers for hanging node variables. More...
 
virtual void fill_in_jacobian_from_nodal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 

Private Member Functions

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

Additional Inherited Members

- Public Types inherited from oomph::SphericalNavierStokesEquations
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 > &)
 
- Static Public Attributes inherited from oomph::SphericalNavierStokesEquations
static Vector< doubleGamma
 Vector to decide whether the stress-divergence form is used or not. More...
 
- Static Public Attributes inherited from oomph::FiniteElement
static double Tolerance_for_singular_jacobian = 1.0e-16
 Tolerance below which the jacobian is considered singular. More...
 
static bool Accept_negative_jacobian = false
 
static bool Suppress_output_while_checking_for_inverted_elements
 
- Static Public Attributes inherited from oomph::GeneralisedElement
static bool Suppress_warning_about_repeated_internal_data
 
static bool Suppress_warning_about_repeated_external_data = true
 
static double Default_fd_jacobian_step = 1.0e-8
 
- Static Protected Member Functions inherited from oomph::RefineableElement
static void check_value_id (const int &n_continuously_interpolated_values, const int &value_id)
 
- Protected Attributes inherited from oomph::SphericalNavierStokesEquations
doubleViscosity_Ratio_pt
 
doubleDensity_Ratio_pt
 
doubleRe_pt
 Pointer to global Reynolds number. More...
 
doubleReSt_pt
 Pointer to global Reynolds number x Strouhal number (=Womersley) More...
 
doubleReInvFr_pt
 
doubleReInvRo_pt
 
Vector< double > * G_pt
 Pointer to global gravity Vector. More...
 
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
MacroElementMacro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepperGeom_object_time_stepper_pt
 
- Protected Attributes inherited from oomph::RefineableElement
TreeTree_pt
 A pointer to a general tree object. More...
 
unsigned Refine_level
 Refinement level. More...
 
bool To_be_refined
 Flag for refinement. More...
 
bool Refinement_is_enabled
 Flag to indicate suppression of any refinement. More...
 
bool Sons_to_be_unrefined
 Flag for unrefinement. More...
 
long Number
 Global element number – for plotting/validation purposes. More...
 
- Static Protected Attributes inherited from oomph::FiniteElement
static const unsigned Default_Initial_Nvalue = 0
 Default value for the number of values at a node. More...
 
static const double Node_location_tolerance = 1.0e-14
 
static const unsigned N2deriv [] = {0, 1, 3, 6}
 
- Static Protected Attributes inherited from oomph::GeneralisedElement
static DenseMatrix< doubleDummy_matrix
 
static std::deque< double * > Dof_pt_deque
 
- Static Protected Attributes inherited from oomph::RefineableElement
static double Max_integrity_tolerance = 1.0e-8
 Max. allowed discrepancy in element integrity check. More...
 

Detailed Description

Refineable version of the Spherical Navier–Stokes equations.

Constructor & Destructor Documentation

◆ RefineableSphericalNavierStokesEquations()

oomph::RefineableSphericalNavierStokesEquations::RefineableSphericalNavierStokesEquations ( )
inline

Empty Constructor.

71  {
72  }
ElementWithZ2ErrorEstimator()
Default empty constructor.
Definition: error_estimator.h:82
RefineableElement()
Definition: refineable_elements.h:188
SphericalNavierStokesEquations()
Definition: spherical_navier_stokes_elements.h:226

Member Function Documentation

◆ fill_in_generic_residual_contribution_spherical_nst()

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

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

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

Reimplemented from oomph::SphericalNavierStokesEquations.

43  {
44  // Find out how many nodes there are
45  const unsigned n_node = nnode();
46 
47  // Get continuous time from timestepper of first node
48  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
49 
50  // Find out how many pressure dofs there are
51  const unsigned n_pres = npres_spherical_nst();
52 
53  // Get the local indices of the nodal coordinates
54  unsigned u_nodal_index[3];
55  for (unsigned i = 0; i < 3; ++i)
56  {
57  u_nodal_index[i] = this->u_index_spherical_nst(i);
58  }
59 
60  // Which nodal value represents the pressure? (Negative if pressure
61  // is not based on nodal interpolation).
62  const int p_index = this->p_nodal_index_spherical_nst();
63 
64  // Local array of booleans that are true if the l-th pressure value is
65  // hanging (avoid repeated virtual function calls)
66  bool pressure_dof_is_hanging[n_pres];
67  // If the pressure is stored at a node
68  if (p_index >= 0)
69  {
70  // Read out whether the pressure is hanging
71  for (unsigned l = 0; l < n_pres; ++l)
72  {
73  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
74  }
75  }
76  // Otherwise the pressure is not stored at a node and so cannot hang
77  else
78  {
79  for (unsigned l = 0; l < n_pres; ++l)
80  {
81  pressure_dof_is_hanging[l] = false;
82  }
83  }
84 
85 
86  // Set up memory for the shape and test functions
87  Shape psif(n_node), testf(n_node);
88  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
89 
90  // Set up memory for pressure shape and test functions
91  Shape psip(n_pres), testp(n_pres);
92 
93  // Set the value of n_intpt
94  const unsigned n_intpt = integral_pt()->nweight();
95 
96  // Set the Vector to hold local coordinates
97  Vector<double> s(2);
98 
99  // Get Physical Variables from Element
100  // Reynolds number must be multiplied by the density ratio
101  const double dens_ratio = density_ratio();
102  const double scaled_re = re() * dens_ratio;
103  const double scaled_re_st = re_st() * dens_ratio;
104  const double scaled_re_inv_ro = re_invro() * dens_ratio;
105  // const double scaled_re_inv_fr = re_invfr()*dens_ratio;
106  // const double visc_ratio = viscosity_ratio();
107  Vector<double> G = g();
108 
109  // Integers to store the local equation and unknown numbers
110  int local_eqn = 0, local_unknown = 0;
111 
112  // Local storage for pointers to hang info objects
113  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
114 
115  // Loop over the integration points
116  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
117  {
118  // Assign values of s
119  for (unsigned i = 0; i < 2; i++)
120  {
121  s[i] = integral_pt()->knot(ipt, i);
122  }
123 
124  // Get the integral weight
125  double w = integral_pt()->weight(ipt);
126 
127  // Call the derivatives of the shape and test functions
129  ipt, psif, dpsifdx, testf, dtestfdx);
130 
131  // Call the pressure shape and test functions
132  pshape_spherical_nst(s, psip, testp);
133 
134  // Premultiply the weights and the Jacobian
135  double W = w * J;
136 
137  // Calculate local values of the pressure and velocity components
138  //--------------------------------------------------------------
139  double interpolated_p = 0.0;
140  Vector<double> interpolated_u(3, 0.0);
141  Vector<double> interpolated_x(2, 0.0);
142  Vector<double> mesh_velocity(2, 0.0);
143  Vector<double> dudt(3, 0.0);
144  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
145 
146  // Calculate pressure
147  for (unsigned l = 0; l < n_pres; l++)
148  {
149  interpolated_p += this->p_spherical_nst(l) * psip(l);
150  }
151 
152  // Calculate velocities and derivatives
153 
154  // Loop over nodes
155  for (unsigned l = 0; l < n_node; l++)
156  {
157  // Cache the shape function
158  double psi_ = psif(l);
159  // Loop over positions separately
160  for (unsigned i = 0; i < 2; i++)
161  {
162  interpolated_x[i] += nodal_position(l, i) * psi_;
163  }
164 
165  // Loop over velocity directions (three of these)
166  for (unsigned i = 0; i < 3; i++)
167  {
168  const double u_value = this->nodal_value(l, u_nodal_index[i]);
169  interpolated_u[i] += u_value * psi_;
170  dudt[i] += du_dt_spherical_nst(l, i) * psi_;
171 
172  // Loop over derivative directions for gradients
173  for (unsigned j = 0; j < 2; j++)
174  {
175  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
176  }
177  }
178 
179  // Only bother to calculate the mesh velocity if we are using an ALE
180  // method
181  if (!ALE_is_disabled)
182  {
183  throw OomphLibError(
184  "ALE is not properly implemented for Refineable Spherical NS yet",
187 
188  // Loop over directions (only DIM (2) of these)
189  for (unsigned i = 0; i < 2; i++)
190  {
191  mesh_velocity[i] += this->dnodal_position_dt(l, i) * psi_;
192  }
193  }
194  } // End of loop over the nodes
195 
196  // Get the user-defined body force terms
197  Vector<double> body_force(3);
199  time, ipt, s, interpolated_x, body_force);
200 
201  // Get the user-defined source function
202  // double
203  // source=this->get_source_spherical_nst(time(),ipt,interpolated_x);
204 
205  // r is the first postition component
206  const double r = interpolated_x[0];
207  const double r2 = r * r;
208  // const double theta = interpolated_x[1];
209  const double sin_theta = sin(interpolated_x[1]);
210  const double cos_theta = cos(interpolated_x[1]);
211  const double cot_theta = cos_theta / sin_theta;
212 
213  const double u_r = interpolated_u[0];
214  const double u_theta = interpolated_u[1];
215  const double u_phi = interpolated_u[2];
216 
217  // Pre-calculate the scaling factor of the jacobian
218  // double W_pure = W;
219 
220  // W *= r*r*sin(theta);
221 
222  // MOMENTUM EQUATIONS
223  //==================
224  // Number of master nodes and storage for the weight of the shape function
225  unsigned n_master = 1;
226  double hang_weight = 1.0;
227 
228  // Loop over the nodes for the test functions/equations
229  //----------------------------------------------------
230  for (unsigned l = 0; l < n_node; l++)
231  {
232  // Local boolean that indicates the hanging status of the node
233  const bool is_node_hanging = node_pt(l)->is_hanging();
234 
235  // If the node is hanging
236  if (is_node_hanging)
237  {
238  // Get the hanging pointer
239  hang_info_pt = node_pt(l)->hanging_pt();
240  // Read out number of master nodes from hanging data
241  n_master = hang_info_pt->nmaster();
242  }
243  // Otherwise the node is its own master
244  else
245  {
246  n_master = 1;
247  }
248 
249  // Loop over the master nodes
250  for (unsigned m = 0; m < n_master; m++)
251  {
252  // Loop over velocity components for equations
253  for (unsigned i = 0; i < 2 + 1; i++)
254  {
255  // Get the equation number
256  // If the node is hanging
257  if (is_node_hanging)
258  {
259  // Get the equation number from the master node
260  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
261  u_nodal_index[i]);
262  // Get the hang weight from the master node
263  hang_weight = hang_info_pt->master_weight(m);
264  }
265  // If the node is not hanging
266  else
267  {
268  // Local equation number
269  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
270  // Node contributes with full weight
271  hang_weight = 1.0;
272  }
273 
274  // If it's not a boundary condition...
275  if (local_eqn >= 0)
276  {
277  // initialise for residual calculation
278  double sum = 0.0;
279 
280  switch (i)
281  {
282  // RADIAL MOMENTUM EQUATION
283  case 0:
284  {
285  // Convective r-terms
286  double conv = r * u_r * interpolated_dudx(0, 0);
287 
288  // Convective theta-terms
289  conv += u_theta * interpolated_dudx(0, 1);
290 
291  // Remaining convective terms
292  conv -= (u_theta * u_theta + u_phi * u_phi);
293 
294  // Add the time-derivative and convective terms
295  sum += (scaled_re_st * dudt[0] * r2 + scaled_re * r * conv) *
296  sin_theta * testf(l);
297 
298  // Add the user-defined body force term
299  sum -= r2 * sin_theta * body_force[0] * testf(l);
300 
301  // Add the Coriolis term
302  sum -= 2.0 * scaled_re_inv_ro * sin_theta * u_phi * r2 *
303  sin_theta * testf(l);
304 
305  // r-derivative test function component of stress tensor
306  sum += (-interpolated_p + 2 * interpolated_dudx(0, 0)) * r2 *
307  sin_theta * dtestfdx(l, 0);
308 
309  // theta-derivative test function component of stress tensor
310  sum += (r * interpolated_dudx(1, 0) - u_theta +
311  interpolated_dudx(0, 1)) *
312  sin_theta * dtestfdx(l, 1);
313 
314  // Undifferentiated test function component of stress tensor
315  sum += 2.0 *
316  ((-r * interpolated_p + +interpolated_dudx(1, 1) +
317  2.0 * u_r) *
318  sin_theta +
319  u_theta * cos_theta) *
320  testf(l);
321  }
322  break;
323 
324  // THETA-COMPONENT MOMENTUM EQUATION
325  case 1:
326  {
327  // All convective terms
328  double conv =
329  (u_r * interpolated_dudx(1, 0) * r +
330  u_theta * interpolated_dudx(1, 1) + u_r * u_theta) *
331  sin_theta -
332  u_phi * u_phi * cos_theta;
333 
334  // Add the time-derivative and convective terms to the
335  // residual
336  sum += (scaled_re_st * r2 * sin_theta * dudt[1] +
337  scaled_re * r * conv) *
338  testf(l);
339 
340  // Add the user-defined body force term
341  sum -= r2 * sin_theta * body_force[1] * testf(l);
342 
343  // Add the Coriolis term
344  sum -= 2.0 * scaled_re_inv_ro * cos_theta * u_phi * r2 *
345  sin_theta * testf(l);
346 
347  // r-derivative test function component of stress tensor
348  sum += (r * interpolated_dudx(1, 0) - u_theta +
349  interpolated_dudx(0, 1)) *
350  r * sin_theta * dtestfdx(l, 0);
351 
352  // theta-derivative test function component of stress tensor
353  sum += (-r * interpolated_p + 2.0 * interpolated_dudx(1, 1) +
354  2 * u_r) *
355  sin_theta * dtestfdx(l, 1);
356 
357  // Undifferentiated test function component of stress tensor
358  sum -= ((r * interpolated_dudx(1, 0) - u_theta +
359  interpolated_dudx(0, 1)) *
360  sin_theta -
361  (-r * interpolated_p + 2.0 * u_r +
362  2.0 * u_theta * cot_theta) *
363  cos_theta) *
364  testf(l);
365  }
366  break;
367 
368  // PHI-COMPONENT MOMENTUM EQUATION
369  case 2:
370 
371  {
372  // Convective r-terms
373  double conv = u_r * interpolated_dudx(2, 0) * r * sin_theta;
374 
375  // Convective theta-terms
376  conv += u_theta * interpolated_dudx(2, 1) * sin_theta;
377 
378  // Remaining convective terms
379  conv += u_phi * (u_r * sin_theta + u_theta * cos_theta);
380 
381  // Add the time-derivative and convective terms
382  sum += (scaled_re_st * r2 * dudt[2] * sin_theta +
383  scaled_re * conv * r) *
384  testf(l);
385 
386  sum -= r2 * sin_theta * body_force[2] * testf(l);
387 
388  // Add the Coriolis term
389  sum += 2.0 * scaled_re_inv_ro *
390  (cos_theta * u_theta + sin_theta * u_r) * r2 *
391  sin_theta * testf(l);
392 
393 
394  // r-derivative test function component of stress tensor
395  sum += (r2 * interpolated_dudx(2, 0) - r * u_phi) *
396  dtestfdx(l, 0) * sin_theta;
397 
398  // theta-derivative test function component of stress tensor
399  sum +=
400  (interpolated_dudx(2, 1) * sin_theta - u_phi * cos_theta) *
401  dtestfdx(l, 1);
402 
403  // Undifferentiated test function component of stress tensor
404  sum -= ((r * interpolated_dudx(2, 0) - u_phi) * sin_theta +
405  (interpolated_dudx(2, 1) - u_phi * cot_theta) *
406  cos_theta) *
407  testf(l);
408  }
409  break;
410  }
411 
412  // Add contribution
413  // Sign changed to be consistent with other NS implementations
414  residuals[local_eqn] -= sum * W * hang_weight;
415 
416  // CALCULATE THE JACOBIAN
417  if (flag)
418  {
419  // Number of master nodes and weights
420  unsigned n_master2 = 1;
421  double hang_weight2 = 1.0;
422  // Loop over the velocity nodes for columns
423  for (unsigned l2 = 0; l2 < n_node; l2++)
424  {
425  // Local boolean for hanging status
426  const bool is_node2_hanging = node_pt(l2)->is_hanging();
427 
428  // If the node is hanging
429  if (is_node2_hanging)
430  {
431  hang_info2_pt = node_pt(l2)->hanging_pt();
432  // Read out number of master nodes from hanging data
433  n_master2 = hang_info2_pt->nmaster();
434  }
435  // Otherwise the node is its own master
436  else
437  {
438  n_master2 = 1;
439  }
440 
441  // Loop over the master nodes
442  for (unsigned m2 = 0; m2 < n_master2; m2++)
443  {
444  // Loop over the velocity components
445  for (unsigned i2 = 0; i2 < 2 + 1; i2++)
446  {
447  // Get the number of the unknown
448  // If the node is hanging
449  if (is_node2_hanging)
450  {
451  // Get the equation number from the master node
452  local_unknown = this->local_hang_eqn(
453  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
454  hang_weight2 = hang_info2_pt->master_weight(m2);
455  }
456  else
457  {
458  local_unknown =
459  this->nodal_local_eqn(l2, u_nodal_index[i2]);
460  hang_weight2 = 1.0;
461  }
462 
463  // If the unknown is non-zero
464  if (local_unknown >= 0)
465  {
466  // Different results for i and i2
467  switch (i)
468  {
469  // RADIAL MOMENTUM EQUATION
470  case 0:
471  switch (i2)
472  {
473  // radial component
474  case 0:
475 
476  // Add the mass matrix entries
477  if (flag == 2)
478  {
479  mass_matrix(local_eqn, local_unknown) +=
480  scaled_re_st * psif(l2) * testf(l) * r2 *
481  sin_theta * W * hang_weight * hang_weight2;
482  }
483 
484  {
485  // Convective r-term contribution
486  double jac_conv =
487  r * (psif(l2) * interpolated_dudx(0, 0) +
488  u_r * dpsifdx(l2, 0));
489 
490  // Convective theta-term contribution
491  jac_conv += u_theta * dpsifdx(l2, 1);
492 
493  // Add the time-derivative contribution and
494  // the convective
495  // contribution to the sum
496  double jac_sum =
497  (scaled_re_st *
499  1, 0) *
500  psif(l2) * r2 +
501  scaled_re * jac_conv * r) *
502  testf(l);
503 
504 
505  // Contribution from the r-derivative test
506  // function
507  // part of stress tensor
508  jac_sum +=
509  2.0 * dpsifdx(l2, 0) * dtestfdx(l, 0) * r2;
510 
511  // Contribution from the theta-derivative
512  // test function part of stress tensor
513  jac_sum += dpsifdx(l2, 1) * dtestfdx(l, 1);
514 
515 
516  // Contribution from the undifferentiated
517  // test function part
518  // of stress tensor
519  jac_sum += 4.0 * psif[l2] * testf(l);
520 
521  // Add the total contribution to the
522  // jacobian multiplied
523  // by the jacobian weight
524  jacobian(local_eqn, local_unknown) -=
525  jac_sum * sin_theta * W * hang_weight *
526  hang_weight2;
527  }
528 
529  break;
530 
531  // axial component
532  case 1:
533  {
534  // No time derivative contribution
535 
536  // Convective contribution
537  double jac_sum =
538  scaled_re *
539  (interpolated_dudx(0, 1) - 2.0 * u_theta) *
540  psif(l2) * r * sin_theta * testf(l);
541 
542  // Contribution from the theta-derivative
543  // test function
544  // part of stress tensor
545  jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) *
546  dtestfdx(l, 1) * sin_theta;
547 
548  // Contribution from the undifferentiated
549  // test function
550  // part of stress tensor
551  jac_sum += 2.0 *
552  (dpsifdx(l2, 1) * sin_theta +
553  psif(l2) * cos_theta) *
554  testf(l);
555 
556  // Add the full contribution to the jacobian
557  // matrix
558  jacobian(local_eqn, local_unknown) -=
559  jac_sum * W * hang_weight * hang_weight2;
560 
561  } // End of i2 = 1 section
562 
563  break;
564 
565  // azimuthal component
566  case 2:
567  {
568  // Single convective-term contribution
569  jacobian(local_eqn, local_unknown) +=
570  2.0 * scaled_re * u_phi * psif[l2] * r *
571  sin_theta * testf[l] * W * hang_weight *
572  hang_weight2;
573 
574  // Add the Coriolis term
575  jacobian(local_eqn, local_unknown) +=
576  2.0 * scaled_re_inv_ro * sin_theta *
577  psif(l2) * r2 * sin_theta * testf[l] * W *
578  hang_weight * hang_weight2;
579  }
580 
581  break;
582  } // End of contribution radial momentum eqn
583  break;
584 
585  // AXIAL MOMENTUM EQUATION
586  case 1:
587  switch (i2)
588  {
589  // radial component
590  case 0:
591  {
592  // Convective terms
593  double jac_sum =
594  scaled_re *
595  (r2 * interpolated_dudx(1, 0) + r * u_theta) *
596  psif(l2) * sin_theta * testf(l);
597 
598  // Contribution from the r-derivative
599  // test function part of stress tensor
600  jac_sum += dpsifdx(l2, 1) * dtestfdx(l, 0) *
601  sin_theta * r;
602 
603  // Contribution from the theta-derivative
604  // test function
605  // part of stress tensor
606  jac_sum +=
607  2.0 * psif(l2) * dtestfdx(l, 1) * sin_theta;
608 
609  // Contribution from the undifferentiated
610  // test function
611  // part of stress tensor
612  jac_sum -= (dpsifdx(l2, 1) * sin_theta -
613  2.0 * psif(l2) * cos_theta) *
614  testf(l);
615 
616  // Add the sum to the jacobian
617  jacobian(local_eqn, local_unknown) -=
618  jac_sum * W * hang_weight * hang_weight2;
619  }
620 
621  break;
622 
623  // axial component
624  case 1:
625 
626  // Add the mass matrix terms
627  if (flag == 2)
628  {
629  mass_matrix(local_eqn, local_unknown) +=
630  scaled_re_st * psif[l2] * testf[l] * W *
631  r2 * sin_theta * hang_weight * hang_weight2;
632  }
633 
634 
635  {
636  // Contribution from the convective terms
637  double jac_conv =
638  r * u_r * dpsifdx(l2, 0) +
639  u_theta * dpsifdx(l2, 1) +
640  (interpolated_dudx(1, 1) + u_r) * psif(l2);
641 
642  // Add the time-derivative term and the
643  // convective terms
644  double jac_sum =
645  (scaled_re_st *
647  1, 0) *
648  psif(l2) * r2 +
649  scaled_re * r * jac_conv) *
650  testf(l) * sin_theta;
651 
652 
653  // Contribution from the r-derivative test
654  // function
655  // part of stress tensor
656  jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) *
657  r * dtestfdx(l, 0) * sin_theta;
658 
659  // Contribution from the theta-derivative
660  // test function
661  // part of stress tensor
662  jac_sum += 2.0 * dpsifdx(l2, 1) *
663  dtestfdx(l, 1) * sin_theta;
664 
665  // Contribution from the undifferentiated
666  // test function
667  // part of stress tensor
668  jac_sum -=
669  ((r * dpsifdx(l2, 0) - psif(l2)) *
670  sin_theta -
671  2.0 * psif(l2) * cot_theta * cos_theta) *
672  testf(l);
673 
674  // Add the contribution to the jacobian
675  jacobian(local_eqn, local_unknown) -=
676  jac_sum * W * hang_weight * hang_weight2;
677  }
678 
679  break;
680 
681  // azimuthal component
682  case 2:
683  {
684  // Only a convective term contribution
685  jacobian(local_eqn, local_unknown) +=
686  scaled_re * 2.0 * r * psif(l2) * u_phi *
687  cos_theta * testf(l) * W * hang_weight *
688  hang_weight2;
689 
690  // Add the Coriolis term
691  jacobian(local_eqn, local_unknown) +=
692  2.0 * scaled_re_inv_ro * cos_theta *
693  psif(l2) * r2 * sin_theta * testf[l] * W *
694  hang_weight * hang_weight2;
695  }
696 
697  break;
698  }
699  break;
700 
701  // AZIMUTHAL MOMENTUM EQUATION
702  case 2:
703  switch (i2)
704  {
705  // radial component
706  case 0:
707  {
708  // Contribution from convective terms
709  jacobian(local_eqn, local_unknown) -=
710  scaled_re *
711  (r * interpolated_dudx(2, 0) + u_phi) *
712  psif(l2) * testf(l) * r * sin_theta * W *
713  hang_weight * hang_weight2;
714 
715  // Coriolis term
716  jacobian(local_eqn, local_unknown) -=
717  2.0 * scaled_re_inv_ro * sin_theta *
718  psif(l2) * r2 * sin_theta * testf[l] * W *
719  hang_weight * hang_weight2;
720  }
721  break;
722 
723  // axial component
724  case 1:
725  {
726  // Contribution from convective terms
727  jacobian(local_eqn, local_unknown) -=
728  scaled_re *
729  (interpolated_dudx(2, 1) * sin_theta +
730  u_phi * cos_theta) *
731  r * psif(l2) * testf(l) * W * hang_weight *
732  hang_weight2;
733 
734  // Coriolis term
735  jacobian(local_eqn, local_unknown) -=
736  2.0 * scaled_re_inv_ro * cos_theta *
737  psif(l2) * r2 * sin_theta * testf[l] * W *
738  hang_weight * hang_weight2;
739  }
740 
741  break;
742 
743  // azimuthal component
744  case 2:
745 
746  // Add the mass matrix terms
747  if (flag == 2)
748  {
749  mass_matrix(local_eqn, local_unknown) +=
750  scaled_re_st * psif[l2] * testf[l] * W *
751  r2 * sin_theta * hang_weight * hang_weight2;
752  }
753 
754  {
755  // Convective terms
756  double jac_conv =
757  r * u_r * dpsifdx(l2, 0) * sin_theta;
758 
759  // Convective theta-term contribution
760  jac_conv +=
761  u_theta * dpsifdx(l2, 1) * sin_theta;
762 
763  // Contribution from the remaining convective
764  // terms
765  jac_conv +=
766  (u_r * sin_theta + u_theta * cos_theta) *
767  psif(l2);
768 
769  // Add the convective and time derivatives
770  double jac_sum =
771  (scaled_re_st *
773  1, 0) *
774  psif(l2) * r2 * sin_theta +
775  scaled_re * r * jac_conv) *
776  testf(l);
777 
778 
779  // Contribution from the r-derivative test
780  // function
781  // part of stress tensor
782  jac_sum += (r * dpsifdx(l2, 0) - psif(l2)) *
783  dtestfdx(l, 0) * r * sin_theta;
784 
785  // Contribution from the theta-derivative
786  // test function
787  // part of stress tensor
788  jac_sum += (dpsifdx(l2, 1) * sin_theta -
789  psif(l2) * cos_theta) *
790  dtestfdx(l, 1);
791 
792  // Contribution from the undifferentiated
793  // test function
794  // part of stress tensor
795  jac_sum -=
796  ((r * dpsifdx(l2, 0) - psif(l2)) *
797  sin_theta +
798  (dpsifdx(l2, 1) - psif(l2) * cot_theta) *
799  cos_theta) *
800  testf(l);
801 
802  // Add to the jacobian
803  jacobian(local_eqn, local_unknown) -=
804  jac_sum * W * hang_weight * hang_weight2;
805  }
806 
807  break;
808  }
809  break;
810  }
811  }
812  }
813  }
814  } // End of loop over the nodes
815 
816 
817  // Loop over the pressure shape functions
818  for (unsigned l2 = 0; l2 < n_pres; l2++)
819  {
820  // If the pressure dof is hanging
821  if (pressure_dof_is_hanging[l2])
822  {
823  // Pressure dof is hanging so it must be nodal-based
824  hang_info2_pt = pressure_node_pt(l2)->hanging_pt(p_index);
825 
826  // Get the number of master nodes from the pressure node
827  n_master2 = hang_info2_pt->nmaster();
828  }
829  // Otherwise the node is its own master
830  else
831  {
832  n_master2 = 1;
833  }
834 
835  // Loop over the master nodes
836  for (unsigned m2 = 0; m2 < n_master2; m2++)
837  {
838  // Get the number of the unknown
839  // If the pressure dof is hanging
840  if (pressure_dof_is_hanging[l2])
841  {
842  // Get the unknown from the node
843  local_unknown = local_hang_eqn(
844  hang_info2_pt->master_node_pt(m2), p_index);
845  // Get the weight from the hanging object
846  hang_weight2 = hang_info2_pt->master_weight(m2);
847  }
848  else
849  {
850  local_unknown = p_local_eqn(l2);
851  hang_weight2 = 1.0;
852  }
853 
854  // If the unknown is not pinned
855  if (local_unknown >= 0)
856  {
857  // Add in contributions to different equations
858  switch (i)
859  {
860  // RADIAL MOMENTUM EQUATION
861  case 0:
862  jacobian(local_eqn, local_unknown) +=
863  psip(l2) *
864  (r2 * dtestfdx(l, 0) + 2.0 * r * testf[l]) * W *
865  sin_theta * hang_weight * hang_weight2;
866 
867 
868  break;
869 
870  // AXIAL MOMENTUM EQUATION
871  case 1:
872  jacobian(local_eqn, local_unknown) +=
873  psip(l2) * r *
874  (dtestfdx(l, 1) * sin_theta +
875  cos_theta * testf(l)) *
876  W * hang_weight * hang_weight2;
877 
878  break;
879 
880  // AZIMUTHAL MOMENTUM EQUATION
881  case 2:
882  break;
883  }
884  }
885  }
886  } // End of loop over pressure dofs
887  } // End of Jacobian calculation
888  } // End of if not boundary condition statement
889  } // End of loop over velocity components
890  } // End of loop over master nodes
891  } // End of loop over nodes
892 
893 
894  // CONTINUITY EQUATION
895  //===================
896 
897  // Loop over the pressure shape functions
898  for (unsigned l = 0; l < n_pres; l++)
899  {
900  // If the pressure dof is hanging
901  if (pressure_dof_is_hanging[l])
902  {
903  // Pressure dof is hanging so it must be nodal-based
904  hang_info_pt = pressure_node_pt(l)->hanging_pt(p_index);
905  // Get the number of master nodes from the pressure node
906  n_master = hang_info_pt->nmaster();
907  }
908  // Otherwise the node is its own master
909  else
910  {
911  n_master = 1;
912  }
913 
914  // Loop over the master nodes
915  for (unsigned m = 0; m < n_master; m++)
916  {
917  // Get the number of the unknown
918  // If the pressure dof is hanging
919  if (pressure_dof_is_hanging[l])
920  {
921  local_eqn =
922  local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
923  hang_weight = hang_info_pt->master_weight(m);
924  }
925  else
926  {
927  local_eqn = p_local_eqn(l);
928  hang_weight = 1.0;
929  }
930 
931  // If the equation is not pinned
932  if (local_eqn >= 0)
933  {
934  // The entire continuity equation
935  residuals[local_eqn] += ((2.0 * u_r + r * interpolated_dudx(0, 0) +
936  interpolated_dudx(1, 1)) *
937  sin_theta +
938  u_theta * cos_theta) *
939  r * testp(l) * W * hang_weight;
940 
941  // CALCULATE THE JACOBIAN
942  //======================
943  if (flag)
944  {
945  unsigned n_master2 = 1;
946  double hang_weight2 = 1.0;
947  // Loop over the velocity nodes for columns
948  for (unsigned l2 = 0; l2 < n_node; l2++)
949  {
950  // Local boolean to indicate whether the node is hanging
951  bool is_node2_hanging = node_pt(l2)->is_hanging();
952 
953  // If the node is hanging
954  if (is_node2_hanging)
955  {
956  hang_info2_pt = node_pt(l2)->hanging_pt();
957  // Read out number of master nodes from hanging data
958  n_master2 = hang_info2_pt->nmaster();
959  }
960  // Otherwise the node is its own master
961  else
962  {
963  n_master2 = 1;
964  }
965 
966  // Loop over the master nodes
967  for (unsigned m2 = 0; m2 < n_master2; m2++)
968  {
969  // Loop over the velocity components
970  for (unsigned i2 = 0; i2 < 2 + 1; i2++)
971  {
972  // Get the number of the unknown
973  // If the node is hanging
974  if (is_node2_hanging)
975  {
976  // Get the equation number from the master node
977  local_unknown = local_hang_eqn(
978  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
979  hang_weight2 = hang_info2_pt->master_weight(m2);
980  }
981  else
982  {
983  local_unknown =
984  this->nodal_local_eqn(l2, u_nodal_index[i2]);
985  hang_weight2 = 1.0;
986  }
987 
988  // If the unknown is not pinned
989  if (local_unknown >= 0)
990  {
991  switch (i2)
992  {
993  // radial component
994  case 0:
995  jacobian(local_eqn, local_unknown) +=
996  (2.0 * psif(l2) + r * dpsifdx(l2, 0)) * r *
997  sin_theta * testp(l) * W * hang_weight *
998  hang_weight2;
999 
1000 
1001  break;
1002 
1003  // axial component
1004  case 1:
1005  jacobian(local_eqn, local_unknown) +=
1006  r *
1007  (dpsifdx(l2, 1) * sin_theta +
1008  psif(l2) * cos_theta) *
1009  testp(l) * W * hang_weight * hang_weight2;
1010 
1011 
1012  break;
1013 
1014  // azimuthal component
1015  case 2:
1016  break;
1017  }
1018  }
1019  }
1020  }
1021  } // End of loop over nodes
1022 
1023  // NO PRESSURE CONTRIBUTIONS TO CONTINUITY EQUATION
1024  } // End of Jacobian calculation
1025  }
1026  }
1027  } // End of loop over pressure nodes
1028 
1029  } // end of loop over integration points
1030  }
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
int i
Definition: BiCGSTAB_step_by_step.cpp:9
JacobiRotation< float > G
Definition: Jacobi_makeGivens.cpp:2
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
RowVector3d w
Definition: Matrix_resize_int.cpp:3
MatrixType m2(n_dims)
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2333
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
HangInfo *const & hanging_pt() const
Definition: nodes.h:1228
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Definition: refineable_elements.h:278
virtual Node * pressure_node_pt(const unsigned &n_p)
Definition: refineable_spherical_navier_stokes_elements.h:53
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.
virtual unsigned u_index_spherical_nst(const unsigned &i) const
Definition: spherical_navier_stokes_elements.h:387
const double & re() const
Reynolds number.
Definition: spherical_navier_stokes_elements.h:249
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
Definition: spherical_navier_stokes_elements.h:311
virtual double p_spherical_nst(const unsigned &n_p) const =0
const double & density_ratio() const
Definition: spherical_navier_stokes_elements.h:287
virtual void get_body_force_spherical_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Definition: spherical_navier_stokes_elements.h:166
double du_dt_spherical_nst(const unsigned &n, const unsigned &i) const
Definition: spherical_navier_stokes_elements.h:395
virtual int p_nodal_index_spherical_nst() const
Definition: spherical_navier_stokes_elements.h:448
const Vector< double > & g() const
Vector of gravitational components.
Definition: spherical_navier_stokes_elements.h:323
virtual int p_local_eqn(const unsigned &n) const =0
virtual unsigned npres_spherical_nst() const =0
Function to return number of pressure degrees of freedom.
bool ALE_is_disabled
Definition: spherical_navier_stokes_elements.h:124
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: spherical_navier_stokes_elements.h:255
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
RealScalar s
Definition: level1_cplx_impl.h:130
int * m
Definition: level2_cplx_impl.h:294
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Definition: axisym_linear_elasticity/cylinder/cylinder.cc:96
double u_r(const double &time, const Vector< double > &x)
Calculate the time dependent form of the r-component of displacement.
Definition: axisym_linear_elasticity/cylinder/cylinder.cc:130
double u_theta(const double &time, const Vector< double > &x)
Calculate the time dependent form of the theta-component of displacement.
Definition: axisym_linear_elasticity/cylinder/cylinder.cc:146
r
Definition: UniformPSDSelfTest.py:20
@ W
Definition: quadtree.h:63
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::SphericalNavierStokesEquations::ALE_is_disabled, Global_Parameters::body_force(), cos(), oomph::SphericalNavierStokesEquations::density_ratio(), oomph::FiniteElement::dnodal_position_dt(), oomph::SphericalNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_spherical_nst(), oomph::SphericalNavierStokesEquations::du_dt_spherical_nst(), G, oomph::SphericalNavierStokesEquations::g(), oomph::SphericalNavierStokesEquations::get_body_force_spherical_nst(), oomph::Node::hanging_pt(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), oomph::Node::is_hanging(), J, j, oomph::Integral::knot(), oomph::RefineableElement::local_hang_eqn(), m, m2(), oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::FiniteElement::node_pt(), oomph::SphericalNavierStokesEquations::npres_spherical_nst(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::SphericalNavierStokesEquations::p_local_eqn(), oomph::SphericalNavierStokesEquations::p_nodal_index_spherical_nst(), oomph::SphericalNavierStokesEquations::p_spherical_nst(), pressure_node_pt(), oomph::SphericalNavierStokesEquations::pshape_spherical_nst(), UniformPSDSelfTest::r, oomph::SphericalNavierStokesEquations::re(), oomph::SphericalNavierStokesEquations::re_invro(), oomph::SphericalNavierStokesEquations::re_st(), s, sin(), oomph::Time::time(), oomph::TimeStepper::time_pt(), oomph::Data::time_stepper_pt(), oomph::SphericalNavierStokesEquations::u_index_spherical_nst(), Global_Parameters::u_r(), Global_Parameters::u_theta(), w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and oomph::TimeStepper::weight().

◆ further_build()

void oomph::RefineableSphericalNavierStokesEquations::further_build ( )
inlinevirtual

Further build: pass the pointers down to the sons.

Reimplemented from oomph::RefineableElement.

Reimplemented in oomph::RefineableQSphericalCrouzeixRaviartElement.

133  {
134  // Find the father element
135  RefineableSphericalNavierStokesEquations* cast_father_element_pt =
137  this->father_element_pt());
138 
139  // Set the viscosity ratio pointer
140  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
141  // Set the density ratio pointer
142  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
143  // Set pointer to global Reynolds number
144  this->Re_pt = cast_father_element_pt->re_pt();
145  // Set pointer to global Reynolds number x Strouhal number (=Womersley)
146  this->ReSt_pt = cast_father_element_pt->re_st_pt();
147  // Set pointer to the global Reynolds number x inverse Rossby number
148  this->ReInvRo_pt = cast_father_element_pt->re_invro_pt();
149  // Set pointer to global Reynolds number x inverse Froude number
150  this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
151  // Set pointer to global gravity Vector
152  this->G_pt = cast_father_element_pt->g_pt();
153 
154  // Set pointer to body force function
155  this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
156 
157  // Set pointer to volumetric source function
158  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
159 
160  // Set the ALE flag
161  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
162  }
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Definition: refineable_elements.h:539
RefineableSphericalNavierStokesEquations()
Empty Constructor.
Definition: refineable_spherical_navier_stokes_elements.h:67
SphericalNavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
Definition: spherical_navier_stokes_elements.h:119
double * Density_Ratio_pt
Definition: spherical_navier_stokes_elements.h:94
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: spherical_navier_stokes_elements.h:102
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: spherical_navier_stokes_elements.h:113
double * ReInvRo_pt
Definition: spherical_navier_stokes_elements.h:110
double * Re_pt
Pointer to global Reynolds number.
Definition: spherical_navier_stokes_elements.h:99
double * Viscosity_Ratio_pt
Definition: spherical_navier_stokes_elements.h:90
double * ReInvFr_pt
Definition: spherical_navier_stokes_elements.h:106
SphericalNavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
Definition: spherical_navier_stokes_elements.h:116

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

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

◆ geometric_jacobian()

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

Fill in the geometric Jacobian, which in this case is r*r*sin(theta)

Reimplemented from oomph::ElementWithZ2ErrorEstimator.

126  {
127  return x[0] * x[0] * sin(x[1]);
128  }
list x
Definition: plotDoE.py:28

References sin(), and plotDoE::x.

◆ get_Z2_flux()

void oomph::RefineableSphericalNavierStokesEquations::get_Z2_flux ( const Vector< double > &  s,
Vector< double > &  flux 
)
inlinevirtual

Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.

Implements oomph::ElementWithZ2ErrorEstimator.

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

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

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

◆ num_Z2_flux_terms()

unsigned oomph::RefineableSphericalNavierStokesEquations::num_Z2_flux_terms ( )
inlinevirtual

◆ pin_elemental_redundant_nodal_pressure_dofs()

virtual void oomph::RefineableSphericalNavierStokesEquations::pin_elemental_redundant_nodal_pressure_dofs ( )
inlineprotectedvirtual

Pin unused nodal pressure dofs (empty by default, because by default pressure dofs are not associated with nodes)

Reimplemented in oomph::RefineableQSphericalTaylorHoodElement.

63 {}

Referenced by pin_redundant_nodal_pressures().

◆ pin_redundant_nodal_pressures()

static void oomph::RefineableSphericalNavierStokesEquations::pin_redundant_nodal_pressures ( const Vector< GeneralisedElement * > &  element_pt)
inlinestatic

Loop over all elements in Vector (which typically contains all the elements in a fluid mesh) and pin the nodal pressure degrees of freedom that are not being used. Function uses the member function

  • RefineableSphericalNavierStokesEquations:: pin_all_nodal_pressure_dofs()

which is empty by default and should be implemented for elements with nodal pressure degrees of freedom (e.g. for refineable Taylor-Hood.)

176  {
177  // Loop over all elements to brutally pin all nodal pressure degrees of
178  // freedom
179  unsigned n_element = element_pt.size();
180  for (unsigned e = 0; e < n_element; e++)
181  {
182  dynamic_cast<RefineableSphericalNavierStokesEquations*>(element_pt[e])
184  }
185  }
Array< double, 1, 3 > e(1./3., 0.5, 2.)
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Definition: refineable_spherical_navier_stokes_elements.h:63

References e(), and pin_elemental_redundant_nodal_pressure_dofs().

◆ pressure_node_pt()

virtual Node* oomph::RefineableSphericalNavierStokesEquations::pressure_node_pt ( const unsigned n_p)
inlineprotectedvirtual

Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interpolation).

Reimplemented in oomph::RefineableQSphericalTaylorHoodElement.

54  {
55  return NULL;
56  }

Referenced by fill_in_generic_residual_contribution_spherical_nst().

◆ unpin_all_pressure_dofs()

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

Unpin all pressure dofs in elements listed in vector.

190  {
191  // Loop over all elements to brutally unpin all nodal pressure degrees of
192  // freedom and internal pressure degrees of freedom
193  unsigned n_element = element_pt.size();
194  for (unsigned e = 0; e < n_element; e++)
195  {
196  dynamic_cast<RefineableSphericalNavierStokesEquations*>(element_pt[e])
198  }
199  }
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.

References e(), and unpin_elemental_pressure_dofs().

◆ unpin_elemental_pressure_dofs()

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

Unpin all pressure dofs in the element.

Implemented in oomph::RefineableQSphericalCrouzeixRaviartElement, and oomph::RefineableQSphericalTaylorHoodElement.

Referenced by unpin_all_pressure_dofs().


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