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

#include <polar_navier_stokes_elements.h>

+ Inheritance diagram for oomph::PolarNavierStokesEquations:

Public Types

typedef void(* NavierStokesBodyForceFctPt) (const double &time, const Vector< double > &x, Vector< double > &body_force)
 
typedef double(* NavierStokesSourceFctPt) (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

 PolarNavierStokesEquations ()
 Constructor: NULL the body force and source function. More...
 
const double & re () const
 Reynolds number. More...
 
const double & alpha () const
 Alpha. More...
 
const double & re_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& alpha_pt ()
 Pointer to Alpha. 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 Vector< double > & g () const
 Vector of gravitational components. More...
 
Vector< double > *& g_pt ()
 Pointer to Vector of gravitational components. More...
 
NavierStokesBodyForceFctPt & body_force_fct_pt ()
 Access function for the body-force pointer. More...
 
NavierStokesBodyForceFctPt body_force_fct_pt () const
 Access function for the body-force pointer. Const version. More...
 
NavierStokesSourceFctPt & source_fct_pt ()
 Access function for the source-function pointer. More...
 
NavierStokesSourceFctPt source_fct_pt () const
 Access function for the source-function pointer. Const version. More...
 
virtual unsigned npres_pnst () const =0
 Function to return number of pressure degrees of freedom. More...
 
virtual double u_pnst (const unsigned &n, const unsigned &i) const =0
 
virtual double u_pnst (const unsigned &t, const unsigned &n, const unsigned &i) const =0
 
virtual unsigned u_index_pnst (const unsigned &i) const
 
double du_dt_pnst (const unsigned &n, const unsigned &i) const
 
virtual double p_pnst (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_pnst ()
 
virtual Node * pressure_node_pt (const unsigned &n_p)
 
double pressure_integral () const
 Integral of pressure over element. More...
 
double dissipation () const
 Return integral of dissipation over element. More...
 
double dissipation (const Vector< double > &s) const
 Return dissipation at local coordinate s. More...
 
double kin_energy () const
 Get integral of kinetic energy over element. More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 Strain-rate tensor: Now returns polar strain. More...
 
void strain_rate_by_r (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 Function to return polar strain multiplied by r. 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 > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
 
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_fct (std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 
void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
void fill_in_contribution_to_residuals (Vector< double > &residuals)
 Compute the element's residual Vector. More...
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void fill_in_contribution_to_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void fill_in_generic_residual_contribution (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 
void interpolated_u_pnst (const Vector< double > &s, Vector< double > &veloc) const
 Compute vector of FE interpolated velocity u at local coordinate s. More...
 
double interpolated_u_pnst (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated velocity u[i] at local coordinate s. More...
 
double interpolated_p_pnst (const Vector< double > &s) const
 Return FE interpolated pressure at local coordinate s. More...
 
double interpolated_dudx_pnst (const Vector< double > &s, const unsigned &i, const unsigned &j) const
 
- 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 void disable_ALE ()
 
virtual void enable_ALE ()
 
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
 

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)=0
 
virtual double dshape_and_dtest_eulerian_pnst (const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_pnst (const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual void pshape_pnst (const Vector< double > &s, Shape &psi) const =0
 Compute the pressure shape functions at local coordinate s. More...
 
virtual void pshape_pnst (const Vector< double > &s, Shape &psi, Shape &test) const =0
 
void get_body_force (double time, const Vector< double > &x, Vector< double > &result)
 Calculate the body force at a given time and Eulerian position. More...
 
double get_source_fct (double time, const Vector< double > &x)
 
- 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 * Alpha_pt
 Pointer to the angle alpha. More...
 
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
 
Vector< double > * G_pt
 Pointer to global gravity Vector. More...
 
NavierStokesBodyForceFctPt Body_force_fct_pt
 Pointer to body force function. More...
 
NavierStokesSourceFctPt Source_fct_pt
 Pointer to volumetric source function. More...
 
- 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 polar Navier–Stokes equations, This contains the generic maths – any concrete implementation must be derived from this.

Member Typedef Documentation

◆ NavierStokesBodyForceFctPt

typedef void(* oomph::PolarNavierStokesEquations::NavierStokesBodyForceFctPt) (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!

◆ NavierStokesSourceFctPt

typedef double(* oomph::PolarNavierStokesEquations::NavierStokesSourceFctPt) (const double &time, const Vector< double > &x)

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

Constructor & Destructor Documentation

◆ PolarNavierStokesEquations()

oomph::PolarNavierStokesEquations::PolarNavierStokesEquations ( )
inline

Constructor: NULL the body force and source function.

245  : Body_force_fct_pt(0), Source_fct_pt(0)
246  {
247  // Set all the Physical parameter pointers to the default value zero
248  Re_pt = &Default_Physical_Constant_Value;
249  ReSt_pt = &Default_Physical_Constant_Value;
250  ReInvFr_pt = &Default_Physical_Constant_Value;
251  G_pt = &Default_Gravity_vector;
252  // Set the Physical ratios to the default value of 1
253  Viscosity_Ratio_pt = &Default_Physical_Ratio_Value;
254  Density_Ratio_pt = &Default_Physical_Ratio_Value;
255  }
oomph::PolarNavierStokesEquations::Default_Physical_Constant_Value
static double Default_Physical_Constant_Value
Navier–Stokes equations static data.
Definition: polar_navier_stokes_elements.h:128
oomph::PolarNavierStokesEquations::ReSt_pt
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: polar_navier_stokes_elements.h:157
oomph::PolarNavierStokesEquations::Viscosity_Ratio_pt
double * Viscosity_Ratio_pt
Definition: polar_navier_stokes_elements.h:142
oomph::PolarNavierStokesEquations::G_pt
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: polar_navier_stokes_elements.h:164
oomph::PolarNavierStokesEquations::Body_force_fct_pt
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
Definition: polar_navier_stokes_elements.h:167
oomph::PolarNavierStokesEquations::Default_Physical_Ratio_Value
static double Default_Physical_Ratio_Value
Navier–Stokes equations static data.
Definition: polar_navier_stokes_elements.h:132
oomph::PolarNavierStokesEquations::Re_pt
double * Re_pt
Pointer to global Reynolds number.
Definition: polar_navier_stokes_elements.h:154
oomph::PolarNavierStokesEquations::Default_Gravity_vector
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
Definition: polar_navier_stokes_elements.h:135
oomph::PolarNavierStokesEquations::Density_Ratio_pt
double * Density_Ratio_pt
Definition: polar_navier_stokes_elements.h:146
oomph::PolarNavierStokesEquations::Source_fct_pt
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
Definition: polar_navier_stokes_elements.h:170
oomph::PolarNavierStokesEquations::ReInvFr_pt
double * ReInvFr_pt
Definition: polar_navier_stokes_elements.h:161

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

Member Function Documentation

◆ alpha()

const double& oomph::PolarNavierStokesEquations::alpha ( ) const
inline

Alpha.

273  {
274  return *Alpha_pt;
275  }
oomph::PolarNavierStokesEquations::Alpha_pt
double * Alpha_pt
Pointer to the angle alpha.
Definition: polar_navier_stokes_elements.h:151

References Alpha_pt.

Referenced by fill_in_generic_residual_contribution(), oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution(), output(), strain_rate(), and strain_rate_by_r().

◆ alpha_pt()

double*& oomph::PolarNavierStokesEquations::alpha_pt ( )
inline

Pointer to Alpha.

291  {
292  return Alpha_pt;
293  }

References Alpha_pt.

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

◆ body_force_fct_pt() [1/2]

NavierStokesBodyForceFctPt& oomph::PolarNavierStokesEquations::body_force_fct_pt ( )
inline

Access function for the body-force pointer.

353  {
354  return Body_force_fct_pt;
355  }

References Body_force_fct_pt.

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

◆ body_force_fct_pt() [2/2]

NavierStokesBodyForceFctPt oomph::PolarNavierStokesEquations::body_force_fct_pt ( ) const
inline

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

359  {
360  return Body_force_fct_pt;
361  }

References Body_force_fct_pt.

◆ compute_error() [1/2]

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

144  {
145  error = 0.0;
146  norm = 0.0;
147 
148  // Vector of local coordinates
149  Vector<double> s(2);
150 
151  // Vector for coordintes
152  Vector<double> x(2);
153 
154  // Set the value of n_intpt
155  unsigned n_intpt = integral_pt()->nweight();
156 
157  outfile << "ZONE" << std::endl;
158 
159  // Exact solution Vector (u,v,[w],p)
160  Vector<double> exact_soln(2 + 1);
161 
162  // Loop over the integration points
163  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
164  {
165  // Assign values of s
166  for (unsigned i = 0; i < 2; i++)
167  {
168  s[i] = integral_pt()->knot(ipt, i);
169  }
170 
171  // Get the integral weight
172  double w = integral_pt()->weight(ipt);
173 
174  // Get jacobian of mapping
175  double J = J_eulerian(s);
176 
177  // Premultiply the weights and the Jacobian
178  double W = w * J;
179 
180  // Get x position as Vector
181  interpolated_x(s, x);
182 
183  // Get exact solution at this point
184  (*exact_soln_pt)(x, exact_soln);
185 
186  // Velocity error
187  for (unsigned i = 0; i < 2; i++)
188  {
189  norm += exact_soln[i] * exact_soln[i] * W;
190  error += (exact_soln[i] - interpolated_u_pnst(s, i)) *
191  (exact_soln[i] - interpolated_u_pnst(s, i)) * W;
192  }
193 
194  // Output x,y,...,u_exact
195  for (unsigned i = 0; i < 2; i++)
196  {
197  outfile << x[i] << " ";
198  }
199 
200  // Output x,y,[z],u_error,v_error,[w_error]
201  for (unsigned i = 0; i < 2; i++)
202  {
203  outfile << exact_soln[i] - interpolated_u_pnst(s, i) << " ";
204  }
205 
206  outfile << std::endl;
207  }
208  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
J
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
oomph::FiniteElement::J_eulerian
virtual double J_eulerian(const Vector< double > &s) const
Definition: elements.cc:4103
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::PolarNavierStokesEquations::interpolated_u_pnst
void interpolated_u_pnst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
Definition: polar_navier_stokes_elements.h:662
s
RealScalar s
Definition: level1_cplx_impl.h:130
ProblemParameters::exact_soln
void exact_soln(const double &time, const Vector< double > &x, Vector< double > &soln)
Definition: unstructured_two_d_curved.cc:301
calibrate.error
int error
Definition: calibrate.py:297
oomph::QuadTreeNames::W
@ W
Definition: quadtree.h:63
plotDoE.x
list x
Definition: plotDoE.py:28

References calibrate::error, ProblemParameters::exact_soln(), i, oomph::FiniteElement::integral_pt(), interpolated_u_pnst(), 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::PolarNavierStokesEquations::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.

67  {
68  error = 0.0;
69  norm = 0.0;
70 
71  // Vector of local coordinates
72  Vector<double> s(2);
73 
74  // Vector for coordintes
75  Vector<double> x(2);
76 
77  // Set the value of n_intpt
78  unsigned n_intpt = integral_pt()->nweight();
79 
80  outfile << "ZONE" << std::endl;
81 
82  // Exact solution Vector (u,v,[w],p)
83  Vector<double> exact_soln(2 + 1);
84 
85  // Loop over the integration points
86  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
87  {
88  // Assign values of s
89  for (unsigned i = 0; i < 2; i++)
90  {
91  s[i] = integral_pt()->knot(ipt, i);
92  }
93 
94  // Get the integral weight
95  double w = integral_pt()->weight(ipt);
96 
97  // Get jacobian of mapping
98  double J = J_eulerian(s);
99 
100  // Premultiply the weights and the Jacobian
101  double W = w * J;
102 
103  // Get x position as Vector
104  interpolated_x(s, x);
105 
106  // Get exact solution at this point
107  (*exact_soln_pt)(time, x, exact_soln);
108 
109  // Velocity error
110  for (unsigned i = 0; i < 2; i++)
111  {
112  norm += exact_soln[i] * exact_soln[i] * W;
113  error += (exact_soln[i] - interpolated_u_pnst(s, i)) *
114  (exact_soln[i] - interpolated_u_pnst(s, i)) * W;
115  }
116 
117  // Output x,y,...,u_exact
118  for (unsigned i = 0; i < 2; i++)
119  {
120  outfile << x[i] << " ";
121  }
122 
123  // Output x,y,[z],u_error,v_error,[w_error]
124  for (unsigned i = 0; i < 2; i++)
125  {
126  outfile << exact_soln[i] - interpolated_u_pnst(s, i) << " ";
127  }
128 
129  outfile << std::endl;
130  }
131  }

References calibrate::error, ProblemParameters::exact_soln(), i, oomph::FiniteElement::integral_pt(), interpolated_u_pnst(), 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.

◆ density_ratio()

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

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

317  {
318  return *Density_Ratio_pt;
319  }

References Density_Ratio_pt.

◆ density_ratio_pt()

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

Pointer to Density ratio.

323  {
324  return Density_Ratio_pt;
325  }

References Density_Ratio_pt.

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

◆ dissipation() [1/2]

double oomph::PolarNavierStokesEquations::dissipation ( ) const

Return integral of dissipation over element.

722  {
723  // Initialise
724  double diss = 0.0;
725 
726  // Set the value of n_intpt
727  unsigned n_intpt = integral_pt()->nweight();
728 
729  // Set the Vector to hold local coordinates
730  Vector<double> s(2);
731 
732  // Loop over the integration points
733  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
734  {
735  // Assign values of s
736  for (unsigned i = 0; i < 2; i++)
737  {
738  s[i] = integral_pt()->knot(ipt, i);
739  }
740 
741  // Get the integral weight
742  double w = integral_pt()->weight(ipt);
743 
744  // Get Jacobian of mapping
745  double J = J_eulerian(s);
746 
747  // Get strain rate matrix
748  DenseMatrix<double> strainrate(2, 2);
749  // DenseDoubleMatrix strainrate(2,2);
750  strain_rate(s, strainrate);
751 
752  // Initialise
753  double local_diss = 0.0;
754  for (unsigned i = 0; i < 2; i++)
755  {
756  for (unsigned j = 0; j < 2; j++)
757  {
758  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
759  }
760  }
761 
762  diss += local_diss * w * J;
763  }
764 
765  return diss;
766  }
oomph::PolarNavierStokesEquations::strain_rate
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: Now returns polar strain.
Definition: polar_navier_stokes_elements.cc:822
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

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

Referenced by full_output().

◆ dissipation() [2/2]

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

Return dissipation at local coordinate s.

799  {
800  // Get strain rate matrix
801  DenseMatrix<double> strainrate(2, 2);
802  // DenseDoubleMatrix strainrate(2,2);
803  strain_rate(s, strainrate);
804 
805  // Initialise
806  double local_diss = 0.0;
807  for (unsigned i = 0; i < 2; i++)
808  {
809  for (unsigned j = 0; j < 2; j++)
810  {
811  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
812  }
813  }
814 
815  return local_diss;
816  }

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

◆ dshape_and_dtest_eulerian_at_knot_pnst()

virtual double oomph::PolarNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_pnst ( 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::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

Referenced by fill_in_generic_residual_contribution(), and oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution().

◆ dshape_and_dtest_eulerian_pnst()

virtual double oomph::PolarNavierStokesEquations::dshape_and_dtest_eulerian_pnst ( 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::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

◆ du_dt_pnst()

double oomph::PolarNavierStokesEquations::du_dt_pnst ( 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.

402  {
403  // Get the data's timestepper
404  TimeStepper* time_stepper_pt = node_pt(n)->time_stepper_pt();
405 
406  // Number of timsteps (past & present)
407  unsigned n_time = time_stepper_pt->ntstorage();
408 
409  // Initialise dudt
410  double dudt = 0.0;
411  // Loop over the timesteps, if there is a non Steady timestepper
412  if (time_stepper_pt->type() != "Steady")
413  {
414  for (unsigned t = 0; t < n_time; t++)
415  {
416  dudt += time_stepper_pt->weight(1, t) * u_pnst(t, n, i);
417  }
418  }
419 
420  return dudt;
421  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::GeomObject::time_stepper_pt
TimeStepper *& time_stepper_pt()
Definition: geom_objects.h:192
oomph::PolarNavierStokesEquations::u_pnst
virtual double u_pnst(const unsigned &n, const unsigned &i) const =0
oomph::TimeStepper::ntstorage
unsigned ntstorage() const
Definition: timesteppers.h:601
oomph::TimeStepper::weight
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
oomph::TimeStepper::type
std::string type() const
Definition: timesteppers.h:490
plotPSD.t
t
Definition: plotPSD.py:36

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

Referenced by full_output().

◆ fill_in_contribution_to_jacobian()

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

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

Reimplemented from oomph::FiniteElement.

577  {
578  // Call the generic routine with the flag set to 1
579  fill_in_generic_residual_contribution(
580  residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
581 
582  /*
583  // Only needed for test purposes
584  //Create a new matrix
585  unsigned n_dof = ndof();
586  DenseMatrix<double> jacobian_fd(n_dof,n_dof,0.0);
587  fill_in_jacobian_from_internal_by_fd(residuals,jacobian_fd);
588  fill_in_jacobian_from_nodal_by_fd(residuals,jacobian_fd);
589 
590  if(n_dof==21)
591  {
592  for(unsigned i=0;i<n_dof;i++)
593  {
594  for(unsigned j=0;j<n_dof;j++)
595  {
596  if(std::fabs(jacobian(i,j) - jacobian_fd(i,j)) > 1.0e-3)
597  {
598  std::cout << i << " " << j << " " << std::fabs(jacobian(i,j) -
599  jacobian_fd(i,j)) << std::endl;
600  }
601  }
602  }
603  }
604  // Only needed for test purposes
605 
606  // Only needed for test purposes
607  //Create a new matrix
608  unsigned n_dof = ndof();
609  DenseMatrix<double>
610  jacobian_new(n_dof,n_dof,0.0),jacobian_old(n_dof,n_dof,0.0);
611  Vector<double> residuals_new(n_dof,0.0),residuals_old(n_dof,0.0);
612 
613  //Call the generic routine with the flag set to 1
614  PolarNavierStokesEquations::fill_in_generic_residual_contribution(residuals_old,jacobian_old,
615  GeneralisedElement::Dummy_matrix,1);
616  //Call the generic routine with the flag set to 1
617  fill_in_generic_residual_contribution(residuals_new,jacobian_new,
618  GeneralisedElement::Dummy_matrix,1);
619 
620  if(n_dof==21)
621  {
622  for(unsigned i=0;i<n_dof;i++)
623  {
624  if(std::fabs(residuals_new[i] - residuals_old[i])>1.0e-8) std::cout << i
625  << " " << std::fabs(residuals_new[i] - residuals_old[i]) << std::endl;
626  //std::cout << residuals_new[i] << std::endl;
627  for(unsigned j=0;j<n_dof;j++)
628  {
629  if(std::fabs(jacobian_new(i,j) - jacobian_old(i,j)) > 1.0e-8)
630  {
631  std::cout << i << " " << j << " " << std::fabs(jacobian_new(i,j) -
632  jacobian_old(i,j)) << std::endl;
633  }
634  }
635  }
636  }
637  // Only needed for test purposes
638  */
639  } // End of fill_in_contribution_to_jacobian
oomph::GeneralisedElement::Dummy_matrix
static DenseMatrix< double > Dummy_matrix
Definition: elements.h:227
oomph::PolarNavierStokesEquations::fill_in_generic_residual_contribution
virtual void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Definition: polar_navier_stokes_elements.cc:1100

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

◆ fill_in_contribution_to_jacobian_and_mass_matrix()

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

Compute the element's residual Vector and the jacobian matrix Plus the mass matrix especially for eigenvalue problems

Reimplemented from oomph::GeneralisedElement.

647  {
648  // Call the generic routine with the flag set to 2
649  fill_in_generic_residual_contribution(
650  residuals, jacobian, mass_matrix, 2);
651  }

References fill_in_generic_residual_contribution().

◆ fill_in_contribution_to_residuals()

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

Compute the element's residual Vector.

Reimplemented from oomph::GeneralisedElement.

565  {
566  // Call the generic residuals function with flag set to 0
567  fill_in_generic_residual_contribution(residuals,
568  GeneralisedElement::Dummy_matrix,
569  GeneralisedElement::Dummy_matrix,
570  0);
571  }

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

◆ fill_in_generic_residual_contribution()

void oomph::PolarNavierStokesEquations::fill_in_generic_residual_contribution ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
DenseMatrix< double > &  mass_matrix,
unsigned  flag 
)
virtual

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

/////////////////////////////////////////////////////////////////// / The finished version of the new equations ///// /////////////////////////////////////////////////////////////////// Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobian as well. flag=2 for Residuals, Jacobian and mass_matrix

This is now my new version with Jacobian and dimensionless phi

Reimplemented in oomph::RefineablePolarNavierStokesEquations.

1105  {
1106  // Find out how many nodes there are
1107  unsigned n_node = nnode();
1108 
1109  // Find out how many pressure dofs there are
1110  unsigned n_pres = npres_pnst();
1111 
1112  // Find the indices at which the local velocities are stored
1113  unsigned u_nodal_index[2];
1114  for (unsigned i = 0; i < 2; i++)
1115  {
1116  u_nodal_index[i] = u_index_pnst(i);
1117  }
1118 
1119  // Set up memory for the shape and test functions
1120  Shape psif(n_node), testf(n_node);
1121  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
1122 
1123  // Set up memory for pressure shape and test functions
1124  Shape psip(n_pres), testp(n_pres);
1125 
1126  // Number of integration points
1127  unsigned n_intpt = integral_pt()->nweight();
1128 
1129  // Set the Vector to hold local coordinates
1130  Vector<double> s(2);
1131 
1132  // Get the reynolds number and Alpha
1133  const double Re = re();
1134  const double Alpha = alpha();
1135  const double Re_St = re_st();
1136 
1137  // Integers to store the local equations and unknowns
1138  int local_eqn = 0, local_unknown = 0;
1139 
1140  // Loop over the integration points
1141  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1142  {
1143  // Assign values of s
1144  for (unsigned i = 0; i < 2; i++) s[i] = integral_pt()->knot(ipt, i);
1145  // Get the integral weight
1146  double w = integral_pt()->weight(ipt);
1147 
1148  // Call the derivatives of the shape and test functions
1149  double J = dshape_and_dtest_eulerian_at_knot_pnst(
1150  ipt, psif, dpsifdx, testf, dtestfdx);
1151 
1152  // Call the pressure shape and test functions
1153  pshape_pnst(s, psip, testp);
1154 
1155  // Premultiply the weights and the Jacobian
1156  double W = w * J;
1157 
1158  // Calculate local values of the pressure and velocity components
1159  // Allocate
1160  double interpolated_p = 0.0;
1161  Vector<double> interpolated_u(2);
1162  Vector<double> interpolated_x(2);
1163  // Vector<double> mesh_velocity(2);
1164  // Vector<double> dudt(2);
1165  DenseMatrix<double> interpolated_dudx(2, 2);
1166 
1167  // Initialise to zero
1168  for (unsigned i = 0; i < 2; i++)
1169  {
1170  // dudt[i] = 0.0;
1171  // mesh_velocity[i] = 0.0;
1172  interpolated_u[i] = 0.0;
1173  interpolated_x[i] = 0.0;
1174  for (unsigned j = 0; j < 2; j++)
1175  {
1176  interpolated_dudx(i, j) = 0.0;
1177  }
1178  }
1179 
1180  // Calculate pressure
1181  for (unsigned l = 0; l < n_pres; l++)
1182  interpolated_p += p_pnst(l) * psip[l];
1183 
1184  // Calculate velocities and derivatives:
1185 
1186  // Loop over nodes
1187  for (unsigned l = 0; l < n_node; l++)
1188  {
1189  // Loop over directions
1190  for (unsigned i = 0; i < 2; i++)
1191  {
1192  // Get the nodal value
1193  double u_value = this->nodal_value(l, u_nodal_index[i]);
1194  interpolated_u[i] += u_value * psif[l];
1195  interpolated_x[i] += this->nodal_position(l, i) * psif[l];
1196  // dudt[i]+=du_dt_pnst(l,i)*psif[l];
1197  // mesh_velocity[i] +=dx_dt(l,i)*psif[l];
1198 
1199  // Loop over derivative directions
1200  for (unsigned j = 0; j < 2; j++)
1201  {
1202  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1203  }
1204  }
1205  }
1206 
1207  // MOMENTUM EQUATIONS
1208  //------------------
1209 
1210  // Loop over the test functions
1211  for (unsigned l = 0; l < n_node; l++)
1212  {
1213  // Can't loop over velocity components as don't have identical
1214  // contributions Do seperately for i = {0,1} instead
1215  unsigned i = 0;
1216  {
1217  /*IF it's not a boundary condition*/
1218  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
1219  if (local_eqn >= 0)
1220  {
1221  // Add the testf[l] term of the stress tensor
1222  residuals[local_eqn] +=
1223  ((interpolated_p / interpolated_x[0]) -
1224  ((1. + Gamma[i]) / pow(interpolated_x[0], 2.)) *
1225  ((1. / Alpha) * interpolated_dudx(1, 1) + interpolated_u[0])) *
1226  testf[l] * interpolated_x[0] * Alpha * W;
1227 
1228  // Add the dtestfdx(l,0) term of the stress tensor
1229  residuals[local_eqn] +=
1230  (interpolated_p - (1. + Gamma[i]) * interpolated_dudx(0, 0)) *
1231  dtestfdx(l, 0) * interpolated_x[0] * Alpha * W;
1232 
1233  // Add the dtestfdx(l,1) term of the stress tensor
1234  residuals[local_eqn] -=
1235  ((1. / (interpolated_x[0] * Alpha)) * interpolated_dudx(0, 1) -
1236  (interpolated_u[1] / interpolated_x[0]) +
1237  Gamma[i] * interpolated_dudx(1, 0)) *
1238  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
1239  interpolated_x[0] * Alpha * W;
1240 
1241  // Convective terms
1242  residuals[local_eqn] -=
1243  Re *
1244  (interpolated_u[0] * interpolated_dudx(0, 0) +
1245  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
1246  interpolated_dudx(0, 1) -
1247  (pow(interpolated_u[1], 2.) / interpolated_x[0])) *
1248  testf[l] * interpolated_x[0] * Alpha * W;
1249 
1250 
1251  // CALCULATE THE JACOBIAN
1252  if (flag)
1253  {
1254  // Loop over the velocity shape functions again
1255  for (unsigned l2 = 0; l2 < n_node; l2++)
1256  {
1257  // Again can't loop over velocity components due to loss of
1258  // symmetry
1259  unsigned i2 = 0;
1260  {
1261  // If at a non-zero degree of freedom add in the entry
1262  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
1263  if (local_unknown >= 0)
1264  {
1265  // Add contribution to Elemental Matrix
1266  jacobian(local_eqn, local_unknown) -=
1267  (1. + Gamma[i]) *
1268  (psif[l2] / pow(interpolated_x[0], 2.)) * testf[l] *
1269  interpolated_x[0] * Alpha * W;
1270 
1271  jacobian(local_eqn, local_unknown) -=
1272  (1. + Gamma[i]) * dpsifdx(l2, 0) * dtestfdx(l, 0) *
1273  interpolated_x[0] * Alpha * W;
1274 
1275  jacobian(local_eqn, local_unknown) -=
1276  (1. / (interpolated_x[0] * Alpha)) * dpsifdx(l2, 1) *
1277  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
1278  interpolated_x[0] * Alpha * W;
1279 
1280  // Now add in the inertial terms
1281  jacobian(local_eqn, local_unknown) -=
1282  Re *
1283  (psif[l2] * interpolated_dudx(0, 0) +
1284  interpolated_u[0] * dpsifdx(l2, 0) +
1285  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
1286  dpsifdx(l2, 1)) *
1287  testf[l] * interpolated_x[0] * Alpha * W;
1288 
1289  // extra bit for mass matrix
1290  if (flag == 2)
1291  {
1292  mass_matrix(local_eqn, local_unknown) +=
1293  Re_St * psif[l2] * testf[l] * interpolated_x[0] *
1294  Alpha * W;
1295  }
1296 
1297  } // End of (Jacobian's) if not boundary condition statement
1298  } // End of i2=0 section
1299 
1300  i2 = 1;
1301  {
1302  // If at a non-zero degree of freedom add in the entry
1303  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
1304  if (local_unknown >= 0)
1305  {
1306  // Add contribution to Elemental Matrix
1307  jacobian(local_eqn, local_unknown) -=
1308  ((1. + Gamma[i]) / (pow(interpolated_x[0], 2.) * Alpha)) *
1309  dpsifdx(l2, 1) * testf[l] * interpolated_x[0] * Alpha * W;
1310 
1311  jacobian(local_eqn, local_unknown) -=
1312  (-(psif[l2] / interpolated_x[0]) +
1313  Gamma[i] * dpsifdx(l2, 0)) *
1314  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
1315  interpolated_x[0] * Alpha * W;
1316 
1317  // Now add in the inertial terms
1318  jacobian(local_eqn, local_unknown) -=
1319  Re *
1320  ((psif[l2] / (interpolated_x[0] * Alpha)) *
1321  interpolated_dudx(0, 1) -
1322  2 * interpolated_u[1] * (psif[l2] / interpolated_x[0])) *
1323  testf[l] * interpolated_x[0] * Alpha * W;
1324 
1325  } // End of (Jacobian's) if not boundary condition statement
1326  } // End of i2=1 section
1327 
1328  } // End of l2 loop
1329 
1330  /*Now loop over pressure shape functions*/
1331  /*This is the contribution from pressure gradient*/
1332  for (unsigned l2 = 0; l2 < n_pres; l2++)
1333  {
1334  /*If we are at a non-zero degree of freedom in the entry*/
1335  local_unknown = p_local_eqn(l2);
1336  if (local_unknown >= 0)
1337  {
1338  jacobian(local_eqn, local_unknown) +=
1339  (psip[l2] / interpolated_x[0]) * testf[l] *
1340  interpolated_x[0] * Alpha * W;
1341 
1342  jacobian(local_eqn, local_unknown) +=
1343  psip[l2] * dtestfdx(l, 0) * interpolated_x[0] * Alpha * W;
1344  }
1345  }
1346  } /*End of Jacobian calculation*/
1347 
1348  } // End of if not boundary condition statement
1349  } // End of i=0 section
1350 
1351  i = 1;
1352  {
1353  /*IF it's not a boundary condition*/
1354  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
1355  if (local_eqn >= 0)
1356  {
1357  // Add the testf[l] term of the stress tensor
1358  residuals[local_eqn] +=
1359  ((1. / (pow(interpolated_x[0], 2.) * Alpha)) *
1360  interpolated_dudx(0, 1) -
1361  (interpolated_u[1] / pow(interpolated_x[0], 2.)) +
1362  Gamma[i] * (1. / interpolated_x[0]) * interpolated_dudx(1, 0)) *
1363  testf[l] * interpolated_x[0] * Alpha * W;
1364 
1365  // Add the dtestfdx(l,0) term of the stress tensor
1366  residuals[local_eqn] -=
1367  (interpolated_dudx(1, 0) +
1368  Gamma[i] *
1369  ((1. / (interpolated_x[0] * Alpha)) * interpolated_dudx(0, 1) -
1370  (interpolated_u[1] / interpolated_x[0]))) *
1371  dtestfdx(l, 0) * interpolated_x[0] * Alpha * W;
1372 
1373  // Add the dtestfdx(l,1) term of the stress tensor
1374  residuals[local_eqn] +=
1375  (interpolated_p -
1376  (1. + Gamma[i]) *
1377  ((1. / (interpolated_x[0] * Alpha)) * interpolated_dudx(1, 1) +
1378  (interpolated_u[0] / interpolated_x[0]))) *
1379  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
1380  interpolated_x[0] * Alpha * W;
1381 
1382  // Convective terms
1383  residuals[local_eqn] -=
1384  Re *
1385  (interpolated_u[0] * interpolated_dudx(1, 0) +
1386  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
1387  interpolated_dudx(1, 1) +
1388  ((interpolated_u[0] * interpolated_u[1]) / interpolated_x[0])) *
1389  testf[l] * interpolated_x[0] * Alpha * W;
1390 
1391  // CALCULATE THE JACOBIAN
1392  if (flag)
1393  {
1394  // Loop over the velocity shape functions again
1395  for (unsigned l2 = 0; l2 < n_node; l2++)
1396  {
1397  // Again can't loop over velocity components due to loss of
1398  // symmetry
1399  unsigned i2 = 0;
1400  {
1401  // If at a non-zero degree of freedom add in the entry
1402  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
1403  if (local_unknown >= 0)
1404  {
1405  // Add contribution to Elemental Matrix
1406  jacobian(local_eqn, local_unknown) +=
1407  (1. / (pow(interpolated_x[0], 2.) * Alpha)) *
1408  dpsifdx(l2, 1) * testf[l] * interpolated_x[0] * Alpha * W;
1409 
1410  jacobian(local_eqn, local_unknown) -=
1411  Gamma[i] * (1. / (interpolated_x[0] * Alpha)) *
1412  dpsifdx(l2, 1) * dtestfdx(l, 0) * interpolated_x[0] *
1413  Alpha * W;
1414 
1415  jacobian(local_eqn, local_unknown) -=
1416  (1 + Gamma[i]) * (psif[l2] / interpolated_x[0]) *
1417  (1. / (interpolated_x[0] * Alpha)) * dtestfdx(l, 1) *
1418  interpolated_x[0] * Alpha * W;
1419 
1420  // Now add in the inertial terms
1421  jacobian(local_eqn, local_unknown) -=
1422  Re *
1423  (psif[l2] * interpolated_dudx(1, 0) +
1424  (psif[l2] * interpolated_u[1] / interpolated_x[0])) *
1425  testf[l] * interpolated_x[0] * Alpha * W;
1426 
1427  } // End of (Jacobian's) if not boundary condition statement
1428  } // End of i2=0 section
1429 
1430  i2 = 1;
1431  {
1432  // If at a non-zero degree of freedom add in the entry
1433  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
1434  if (local_unknown >= 0)
1435  {
1436  // Add contribution to Elemental Matrix
1437  jacobian(local_eqn, local_unknown) +=
1438  (-(psif[l2] / pow(interpolated_x[0], 2.)) +
1439  Gamma[i] * (1. / interpolated_x[0]) * dpsifdx(l2, 0)) *
1440  testf[l] * interpolated_x[0] * Alpha * W;
1441 
1442  jacobian(local_eqn, local_unknown) -=
1443  (dpsifdx(l2, 0) -
1444  Gamma[i] * (psif[l2] / interpolated_x[0])) *
1445  dtestfdx(l, 0) * interpolated_x[0] * Alpha * W;
1446 
1447  jacobian(local_eqn, local_unknown) -=
1448  (1. + Gamma[i]) * (1. / (interpolated_x[0] * Alpha)) *
1449  dpsifdx(l2, 1) * (1. / (interpolated_x[0] * Alpha)) *
1450  dtestfdx(l, 1) * interpolated_x[0] * Alpha * W;
1451 
1452  // Now add in the inertial terms
1453  jacobian(local_eqn, local_unknown) -=
1454  Re *
1455  (interpolated_u[0] * dpsifdx(l2, 0) +
1456  (psif[l2] / (interpolated_x[0] * Alpha)) *
1457  interpolated_dudx(1, 1) +
1458  (interpolated_u[1] / (interpolated_x[0] * Alpha)) *
1459  dpsifdx(l2, 1) +
1460  (interpolated_u[0] * psif[l2] / interpolated_x[0])) *
1461  testf[l] * interpolated_x[0] * Alpha * W;
1462 
1463  // extra bit for mass matrix
1464  if (flag == 2)
1465  {
1466  mass_matrix(local_eqn, local_unknown) +=
1467  Re_St * psif[l2] * testf[l] * interpolated_x[0] *
1468  Alpha * W;
1469  }
1470 
1471  } // End of (Jacobian's) if not boundary condition statement
1472  } // End of i2=1 section
1473 
1474  } // End of l2 loop
1475 
1476  /*Now loop over pressure shape functions*/
1477  /*This is the contribution from pressure gradient*/
1478  for (unsigned l2 = 0; l2 < n_pres; l2++)
1479  {
1480  /*If we are at a non-zero degree of freedom in the entry*/
1481  local_unknown = p_local_eqn(l2);
1482  if (local_unknown >= 0)
1483  {
1484  jacobian(local_eqn, local_unknown) +=
1485  (psip[l2] / interpolated_x[0]) * (1. / Alpha) *
1486  dtestfdx(l, 1) * interpolated_x[0] * Alpha * W;
1487  }
1488  }
1489  } /*End of Jacobian calculation*/
1490 
1491  } // End of if not boundary condition statement
1492 
1493  } // End of i=1 section
1494 
1495  } // End of loop over shape functions
1496 
1497  // CONTINUITY EQUATION
1498  //-------------------
1499 
1500  // Loop over the shape functions
1501  for (unsigned l = 0; l < n_pres; l++)
1502  {
1503  local_eqn = p_local_eqn(l);
1504  // If not a boundary conditions
1505  if (local_eqn >= 0)
1506  {
1507  residuals[local_eqn] +=
1508  (interpolated_dudx(0, 0) + (interpolated_u[0] / interpolated_x[0]) +
1509  (1. / (interpolated_x[0] * Alpha)) * interpolated_dudx(1, 1)) *
1510  testp[l] * interpolated_x[0] * Alpha * W;
1511 
1512 
1513  /*CALCULATE THE JACOBIAN*/
1514  if (flag)
1515  {
1516  /*Loop over the velocity shape functions*/
1517  for (unsigned l2 = 0; l2 < n_node; l2++)
1518  {
1519  unsigned i2 = 0;
1520  {
1521  /*If we're at a non-zero degree of freedom add it in*/
1522  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
1523  if (local_unknown >= 0)
1524  {
1525  jacobian(local_eqn, local_unknown) +=
1526  (dpsifdx(l2, 0) + (psif[l2] / interpolated_x[0])) *
1527  testp[l] * interpolated_x[0] * Alpha * W;
1528  }
1529  } // End of i2=0 section
1530 
1531  i2 = 1;
1532  {
1533  /*If we're at a non-zero degree of freedom add it in*/
1534  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
1535  if (local_unknown >= 0)
1536  {
1537  jacobian(local_eqn, local_unknown) +=
1538  (1. / (interpolated_x[0] * Alpha)) * dpsifdx(l2, 1) *
1539  testp[l] * interpolated_x[0] * Alpha * W;
1540  }
1541  } // End of i2=1 section
1542 
1543  } /*End of loop over l2*/
1544  } /*End of Jacobian calculation*/
1545 
1546  } // End of if not boundary condition
1547  } // End of loop over l
1548 
1549 
1550  } // End of loop over integration points
1551 
1552  } // End of add_generic_residual_contribution
oomph::FiniteElement::nodal_value
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
oomph::FiniteElement::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::FiniteElement::nodal_position
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317
oomph::PolarNavierStokesEquations::pshape_pnst
virtual void pshape_pnst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
oomph::PolarNavierStokesEquations::Gamma
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition: polar_navier_stokes_elements.h:261
oomph::PolarNavierStokesEquations::p_local_eqn
virtual int p_local_eqn(const unsigned &n)=0
oomph::PolarNavierStokesEquations::alpha
const double & alpha() const
Alpha.
Definition: polar_navier_stokes_elements.h:272
oomph::PolarNavierStokesEquations::p_pnst
virtual double p_pnst(const unsigned &n_p) const =0
oomph::PolarNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_pnst
virtual double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
oomph::PolarNavierStokesEquations::u_index_pnst
virtual unsigned u_index_pnst(const unsigned &i) const
Definition: polar_navier_stokes_elements.h:394
oomph::PolarNavierStokesEquations::re
const double & re() const
Reynolds number.
Definition: polar_navier_stokes_elements.h:266
oomph::PolarNavierStokesEquations::re_st
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: polar_navier_stokes_elements.h:278
oomph::PolarNavierStokesEquations::npres_pnst
virtual unsigned npres_pnst() const =0
Function to return number of pressure degrees of freedom.
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
GlobalPhysicalVariables::Re
double Re
Reynolds number.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:61
TanhSolnForAdvectionDiffusion::Alpha
double Alpha
Parameter for steepness of step.
Definition: two_d_adv_diff_adapt.cc:53
oomph::Problem_Parameter::Re_St
double Re_St
Reynolds multiplied Strouhal.
Definition: axisym_vibrating_shell_non_newtonian.cc:91

References TanhSolnForAdvectionDiffusion::Alpha, alpha(), dshape_and_dtest_eulerian_at_knot_pnst(), Gamma, i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), J, j, oomph::Integral::knot(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), npres_pnst(), oomph::Integral::nweight(), p_local_eqn(), p_pnst(), Eigen::bfloat16_impl::pow(), pshape_pnst(), GlobalPhysicalVariables::Re, re(), oomph::Problem_Parameter::Re_St, re_st(), s, u_index_pnst(), w, oomph::QuadTreeNames::W, and oomph::Integral::weight().

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::PolarNavierStokesEquations::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::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

◆ full_output() [1/2]

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

511  {
512  unsigned nplot = 5;
513  full_output(outfile, nplot);
514  }
oomph::PolarNavierStokesEquations::full_output
void full_output(std::ostream &outfile)
Definition: polar_navier_stokes_elements.h:510

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

◆ full_output() [2/2]

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

609  {
610  // Vector of local coordinates
611  Vector<double> s(2);
612 
613  // Acceleration
614  Vector<double> dudt(2);
615 
616  // Mesh elocity
617  Vector<double> mesh_veloc(2);
618 
619  // Tecplot header info
620  outfile << tecplot_zone_string(nplot);
621 
622  // Find out how many nodes there are
623  unsigned n_node = nnode();
624 
625  // Set up memory for the shape functions
626  Shape psif(n_node);
627  DShape dpsifdx(n_node, 2);
628 
629  // Loop over plot points
630  unsigned num_plot_points = nplot_points(nplot);
631  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
632  {
633  // Get local coordinates of plot point
634  get_s_plot(iplot, nplot, s);
635 
636  // Call the derivatives of the shape and test functions
637  dshape_eulerian(s, psif, dpsifdx);
638 
639  // Allocate storage
640  Vector<double> mesh_veloc(2);
641  Vector<double> dudt(2);
642  Vector<double> dudt_ALE(2);
643  DenseMatrix<double> interpolated_dudx(2, 2);
644  // DenseDoubleMatrix interpolated_dudx(2,2);
645 
646  // Initialise everything to zero
647  for (unsigned i = 0; i < 2; i++)
648  {
649  mesh_veloc[i] = 0.0;
650  dudt[i] = 0.0;
651  dudt_ALE[i] = 0.0;
652  for (unsigned j = 0; j < 2; j++)
653  {
654  interpolated_dudx(i, j) = 0.0;
655  }
656  }
657 
658  // Calculate velocities and derivatives
659 
660  // Loop over nodes
661  for (unsigned l = 0; l < n_node; l++)
662  {
663  // Loop over directions
664  for (unsigned i = 0; i < 2; i++)
665  {
666  dudt[i] += du_dt_pnst(l, i) * psif[l];
667  mesh_veloc[i] += dnodal_position_dt(l, i) * psif[l];
668 
669  // Loop over derivative directions for velocity gradients
670  for (unsigned j = 0; j < 2; j++)
671  {
672  interpolated_dudx(i, j) += u_pnst(l, i) * dpsifdx(l, j);
673  }
674  }
675  }
676 
677 
678  // Get dudt in ALE form (incl mesh veloc)
679  for (unsigned i = 0; i < 2; i++)
680  {
681  dudt_ALE[i] = dudt[i];
682  for (unsigned k = 0; k < 2; k++)
683  {
684  dudt_ALE[i] -= mesh_veloc[k] * interpolated_dudx(i, k);
685  }
686  }
687 
688 
689  // Coordinates
690  for (unsigned i = 0; i < 2; i++)
691  {
692  outfile << interpolated_x(s, i) << " ";
693  }
694 
695  // Velocities
696  for (unsigned i = 0; i < 2; i++)
697  {
698  outfile << interpolated_u_pnst(s, i) << " ";
699  }
700 
701  // Pressure
702  outfile << interpolated_p_pnst(s) << " ";
703 
704  // Accelerations
705  for (unsigned i = 0; i < 2; i++)
706  {
707  outfile << dudt_ALE[i] << " ";
708  }
709 
710  // Dissipation
711  outfile << dissipation(s) << " ";
712 
713 
714  outfile << std::endl;
715  }
716  }
oomph::FiniteElement::tecplot_zone_string
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Definition: elements.h:3161
oomph::FiniteElement::get_s_plot
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Definition: elements.h:3148
oomph::FiniteElement::nplot_points
virtual unsigned nplot_points(const unsigned &nplot) const
Definition: elements.h:3186
oomph::FiniteElement::dshape_eulerian
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Definition: elements.cc:3298
oomph::FiniteElement::dnodal_position_dt
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2333
oomph::PolarNavierStokesEquations::dissipation
double dissipation() const
Return integral of dissipation over element.
Definition: polar_navier_stokes_elements.cc:721
oomph::PolarNavierStokesEquations::du_dt_pnst
double du_dt_pnst(const unsigned &n, const unsigned &i) const
Definition: polar_navier_stokes_elements.h:401
oomph::PolarNavierStokesEquations::interpolated_p_pnst
double interpolated_p_pnst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
Definition: polar_navier_stokes_elements.h:706
k
char char char int int * k
Definition: level2_impl.h:374

References dissipation(), oomph::FiniteElement::dnodal_position_dt(), oomph::FiniteElement::dshape_eulerian(), du_dt_pnst(), oomph::FiniteElement::get_s_plot(), i, interpolated_p_pnst(), interpolated_u_pnst(), oomph::FiniteElement::interpolated_x(), j, k, oomph::FiniteElement::nnode(), oomph::FiniteElement::nplot_points(), s, oomph::FiniteElement::tecplot_zone_string(), and u_pnst().

◆ g()

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

Vector of gravitational components.

341  {
342  return *G_pt;
343  }

References G_pt.

◆ g_pt()

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

Pointer to Vector of gravitational components.

347  {
348  return G_pt;
349  }

References G_pt.

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

◆ get_body_force()

void oomph::PolarNavierStokesEquations::get_body_force ( double  time,
const Vector< double > &  x,
Vector< double > &  result 
)
inlineprotected

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

210  {
211  // If the function pointer is zero return zero
212  if (Body_force_fct_pt == 0)
213  {
214  // Loop over dimensions and set body forces to zero
215  for (unsigned i = 0; i < 2; i++)
216  {
217  result[i] = 0.0;
218  }
219  }
220  // Otherwise call the function
221  else
222  {
223  (*Body_force_fct_pt)(time, x, result);
224  }
225  }

References Body_force_fct_pt, i, and plotDoE::x.

◆ get_load()

void oomph::PolarNavierStokesEquations::get_load ( const Vector< double > &  s,
const Vector< double > &  xi,
const Vector< double > &  x,
const Vector< double > &  N,
Vector< double > &  load 
)
inline

The potential loading on an external SolidElement is always provided by the traction function

477  {
478  get_traction(s, N, load);
479  }
load
void load(Archive &ar, ParticleHandler &handl)
Definition: Particles.h:21
oomph::PolarNavierStokesEquations::get_traction
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
Definition: polar_navier_stokes_elements.cc:772
N
@ N
Definition: constructor.cpp:22

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

◆ get_source_fct()

double oomph::PolarNavierStokesEquations::get_source_fct ( double  time,
const Vector< double > &  x 
)
inlineprotected

Calculate the source fct at given time and Eulerian position

230  {
231  // If the function pointer is zero return zero
232  if (Source_fct_pt == 0)
233  {
234  return 0;
235  }
236  // Otherwise call the function
237  else
238  {
239  return (*Source_fct_pt)(time, x);
240  }
241  }

References Source_fct_pt, and plotDoE::x.

◆ get_traction()

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

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

775  {
776  // Get velocity gradients
777  DenseMatrix<double> strainrate(2, 2);
778  // DenseDoubleMatrix strainrate(2,2);
779  strain_rate(s, strainrate);
780 
781  // Get pressure
782  double press = interpolated_p_pnst(s);
783 
784  // Loop over traction components
785  for (unsigned i = 0; i < 2; i++)
786  {
787  traction[i] = -press * N[i];
788  for (unsigned j = 0; j < 2; j++)
789  {
790  traction[i] += 2.0 * strainrate(i, j) * N[j];
791  }
792  }
793  }

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

Referenced by get_load().

◆ interpolated_dudx_pnst()

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

//////////////////////////////////////////////////////////////// My own work: /// Return FE interpolated velocity derivative du[i]/dx[j] /// at local coordinate s /// ////////////////////////////////////////////////////////////////

734  {
735  // Find number of nodes
736  unsigned n_node = nnode();
737 
738  // Set up memory for the shape and test functions
739  Shape psif(n_node);
740  DShape dpsifdx(n_node, 2);
741 
742  // double J =
743  dshape_eulerian(s, psif, dpsifdx);
744 
745  // Initialise to zero
746  double interpolated_dudx = 0.0;
747 
748  // Calculate velocity derivative:
749 
750  // Loop over nodes
751  for (unsigned l = 0; l < n_node; l++)
752  {
753  interpolated_dudx += u_pnst(l, i) * dpsifdx(l, j);
754  }
755 
756  return (interpolated_dudx);
757  }

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

◆ interpolated_p_pnst()

double oomph::PolarNavierStokesEquations::interpolated_p_pnst ( const Vector< double > &  s) const
inline

Return FE interpolated pressure at local coordinate s.

707  {
708  // Find number of nodes
709  unsigned n_pres = npres_pnst();
710  // Local shape function
711  Shape psi(n_pres);
712  // Find values of shape function
713  pshape_pnst(s, psi);
714 
715  // Initialise value of p
716  double interpolated_p = 0.0;
717  // Loop over the local nodes and sum
718  for (unsigned l = 0; l < n_pres; l++)
719  {
720  interpolated_p += p_pnst(l) * psi[l];
721  }
722 
723  return (interpolated_p);
724  }

References npres_pnst(), p_pnst(), pshape_pnst(), and s.

Referenced by full_output(), oomph::RefineablePolarCrouzeixRaviartElement::further_build(), oomph::RefineablePolarTaylorHoodElement::get_interpolated_values(), get_traction(), output(), and pressure_integral().

◆ interpolated_u_pnst() [1/2]

double oomph::PolarNavierStokesEquations::interpolated_u_pnst ( const Vector< double > &  s,
const unsigned &  i 
) const
inline

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

686  {
687  // Find number of nodes
688  unsigned n_node = nnode();
689  // Local shape function
690  Shape psi(n_node);
691  // Find values of shape function
692  shape(s, psi);
693 
694  // Initialise value of u
695  double interpolated_u = 0.0;
696  // Loop over the local nodes and sum
697  for (unsigned l = 0; l < n_node; l++)
698  {
699  interpolated_u += u_pnst(l, i) * psi[l];
700  }
701 
702  return (interpolated_u);
703  }
oomph::FiniteElement::shape
virtual void shape(const Vector< double > &s, Shape &psi) const =0

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

◆ interpolated_u_pnst() [2/2]

void oomph::PolarNavierStokesEquations::interpolated_u_pnst ( const Vector< double > &  s,
Vector< double > &  veloc 
) const
inline

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

664  {
665  // Find number of nodes
666  unsigned n_node = nnode();
667  // Local shape function
668  Shape psi(n_node);
669  // Find values of shape function
670  shape(s, psi);
671 
672  for (unsigned i = 0; i < 2; i++)
673  {
674  // Initialise value of u
675  veloc[i] = 0.0;
676  // Loop over the local nodes and sum
677  for (unsigned l = 0; l < n_node; l++)
678  {
679  veloc[i] += u_pnst(l, i) * psi[l];
680  }
681  }
682  }

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

Referenced by compute_error(), full_output(), oomph::RefineablePolarTaylorHoodElement::get_interpolated_values(), oomph::RefineablePolarCrouzeixRaviartElement::get_interpolated_values(), kin_energy(), and output().

◆ kin_energy()

double oomph::PolarNavierStokesEquations::kin_energy ( ) const

Get integral of kinetic energy over element.

Get integral of kinetic energy over element:

1007  {
1008  // Initialise
1009  double kin_en = 0.0;
1010 
1011  // Set the value of n_intpt
1012  unsigned n_intpt = integral_pt()->nweight();
1013 
1014  // Set the Vector to hold local coordinates
1015  Vector<double> s(2);
1016 
1017  // Loop over the integration points
1018  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1019  {
1020  // Assign values of s
1021  for (unsigned i = 0; i < 2; i++)
1022  {
1023  s[i] = integral_pt()->knot(ipt, i);
1024  }
1025 
1026  // Get the integral weight
1027  double w = integral_pt()->weight(ipt);
1028 
1029  // Get Jacobian of mapping
1030  double J = J_eulerian(s);
1031 
1032  // Loop over directions
1033  double veloc_squared = 0.0;
1034  for (unsigned i = 0; i < 2; i++)
1035  {
1036  veloc_squared += interpolated_u_pnst(s, i) * interpolated_u_pnst(s, i);
1037  }
1038 
1039  kin_en += 0.5 * veloc_squared * w * J;
1040  }
1041 
1042  return kin_en;
1043  }

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

◆ npres_pnst()

virtual unsigned oomph::PolarNavierStokesEquations::npres_pnst ( ) const
pure virtual

Function to return number of pressure degrees of freedom.

Implemented in oomph::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

Referenced by fill_in_generic_residual_contribution(), oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution(), and interpolated_p_pnst().

◆ output() [1/4]

void oomph::PolarNavierStokesEquations::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::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

498  {
499  unsigned nplot = 5;
500  output(file_pt, nplot);
501  }
oomph::PolarNavierStokesEquations::output
void output(std::ostream &outfile)
Definition: polar_navier_stokes_elements.h:485

References output().

◆ output() [2/4]

void oomph::PolarNavierStokesEquations::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::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

564  {
565  // Vector of local coordinates
566  Vector<double> s(2);
567 
568  // Tecplot header info
569  // outfile << tecplot_zone_string(nplot);
570  fprintf(file_pt, "%s", tecplot_zone_string(nplot).c_str());
571 
572  // Loop over plot points
573  unsigned num_plot_points = nplot_points(nplot);
574  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
575  {
576  // Get local coordinates of plot point
577  get_s_plot(iplot, nplot, s);
578 
579  // Coordinates
580  for (unsigned i = 0; i < 2; i++)
581  {
582  // outfile << interpolated_x(s,i) << " ";
583  fprintf(file_pt, "%g ", interpolated_x(s, i));
584  }
585 
586  // Velocities
587  for (unsigned i = 0; i < 2; i++)
588  {
589  // outfile << interpolated_u_pnst(s,i) << " ";
590  fprintf(file_pt, "%g ", interpolated_u_pnst(s, i));
591  }
592 
593  // Pressure
594  // outfile << interpolated_p_pnst(s) << " ";
595  fprintf(file_pt, "%g \n", interpolated_p_pnst(s));
596  }
597  fprintf(file_pt, "\n");
598  }

References oomph::FiniteElement::get_s_plot(), i, interpolated_p_pnst(), interpolated_u_pnst(), oomph::FiniteElement::interpolated_x(), oomph::FiniteElement::nplot_points(), s, and oomph::FiniteElement::tecplot_zone_string().

◆ output() [3/4]

void oomph::PolarNavierStokesEquations::output ( std::ostream &  outfile)
inlinevirtual

Output functionget_vels(const Vector<double>& x_to_get, Vector<double>& vels): x,y,[z],u,v,[w],p in tecplot format. Default number of plot points

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

486  {
487  unsigned nplot = 5;
488  output(outfile, nplot);
489  }

Referenced by output(), oomph::PolarCrouzeixRaviartElement::output(), and oomph::PolarTaylorHoodElement::output().

◆ output() [4/4]

void oomph::PolarNavierStokesEquations::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::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

396  {
397  // Vector of local coordinates
398  Vector<double> s(2);
399 
400  // Tecplot header info
401  outfile << tecplot_zone_string(nplot);
402 
403  // Loop over plot points
404  unsigned num_plot_points = nplot_points(nplot);
405  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
406  {
407  // Get local coordinates of plot point
408  get_s_plot(iplot, nplot, s);
409 
410  // Work out global physical coordinate
411  const double Alpha = alpha();
412  double r = interpolated_x(s, 0);
413  double phi = interpolated_x(s, 1);
414  double theta = Alpha * phi;
415 
416  // Coordinates
417  outfile << r * cos(theta) << " " << r * sin(theta) << " ";
418 
419  // Velocities
420  outfile << interpolated_u_pnst(s, 0) * cos(theta) -
421  interpolated_u_pnst(s, 1) * sin(theta)
422  << " ";
423  outfile << interpolated_u_pnst(s, 0) * sin(theta) +
424  interpolated_u_pnst(s, 1) * cos(theta)
425  << " ";
426 
427  // Pressure
428  outfile << interpolated_p_pnst(s) << " ";
429 
430  // Radial and Azimuthal velocities
431  outfile << interpolated_u_pnst(s, 0) << " " << interpolated_u_pnst(s, 1)
432  << " ";
433  // comment start here
434  /*
435  double similarity_solution,dsimilarity_solution;
436  get_similarity_solution(theta,Alpha,similarity_solution,dsimilarity_solution);
437  // Similarity solution:
438  outfile << similarity_solution/(Alpha*r) << " ";
439  // Error from similarity solution:
440  outfile << interpolated_u_pnst(s,0) - similarity_solution/(Alpha*r) << "
441  ";
442  */
443  /*
444  //Work out the Stokes flow similarity solution
445  double mult = (Alpha/(sin(2.*Alpha)-2.*Alpha*cos(2.*Alpha)));
446  double similarity_solution = mult*(cos(2.*Alpha*phi)-cos(2.*Alpha));
447  // Similarity solution:
448  outfile << similarity_solution/(Alpha*r) << " ";
449  // Error from similarity solution:
450  outfile << interpolated_u_pnst(s,0) - similarity_solution/(Alpha*r) << "
451  ";
452  //comment end here
453  */
454  outfile << 0 << " " << 1 << " ";
455 
456  // r and theta for better plotting
457  outfile << r << " " << phi << " ";
458 
459  outfile << std::endl;
460  }
461  outfile << std::endl;
462  }
cos
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
sin
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
BiharmonicTestFunctions2::theta
double theta
Definition: two_d_biharmonic.cc:236
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20

References TanhSolnForAdvectionDiffusion::Alpha, alpha(), cos(), oomph::FiniteElement::get_s_plot(), interpolated_p_pnst(), interpolated_u_pnst(), oomph::FiniteElement::interpolated_x(), oomph::FiniteElement::nplot_points(), UniformPSDSelfTest::r, s, sin(), oomph::FiniteElement::tecplot_zone_string(), and BiharmonicTestFunctions2::theta.

◆ output_fct() [1/2]

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

276  {
277  // Vector of local coordinates
278  Vector<double> s(2);
279 
280  // Vector for coordintes
281  Vector<double> x(2);
282 
283  // Tecplot header info
284  outfile << tecplot_zone_string(nplot);
285 
286  // Exact solution Vector
287  Vector<double> exact_soln;
288 
289  // Loop over plot points
290  unsigned num_plot_points = nplot_points(nplot);
291  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
292  {
293  // Get local coordinates of plot point
294  get_s_plot(iplot, nplot, s);
295 
296  // Get x position as Vector
297  interpolated_x(s, x);
298 
299  // Get exact solution at this point
300  (*exact_soln_pt)(time, x, exact_soln);
301 
302  // Output x,y,...
303  for (unsigned i = 0; i < 2; i++)
304  {
305  outfile << x[i] << " ";
306  }
307 
308  // Output "exact solution"
309  for (unsigned i = 0; i < exact_soln.size(); i++)
310  {
311  outfile << exact_soln[i] << " ";
312  }
313 
314  outfile << std::endl;
315  }
316 
317  // Write tecplot footer (e.g. FE connectivity lists)
318  write_tecplot_zone_footer(outfile, nplot);
319  }
oomph::FiniteElement::write_tecplot_zone_footer
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Definition: elements.h:3174

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

220  {
221  // Vector of local coordinates
222  Vector<double> s(2);
223 
224  // Vector for coordintes
225  Vector<double> x(2);
226 
227  // Tecplot header info
228  outfile << tecplot_zone_string(nplot);
229 
230  // Exact solution Vector
231  Vector<double> exact_soln;
232 
233  // Loop over plot points
234  unsigned num_plot_points = nplot_points(nplot);
235  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
236  {
237  // Get local coordinates of plot point
238  get_s_plot(iplot, nplot, s);
239 
240  // Get x position as Vector
241  interpolated_x(s, x);
242 
243  // Get exact solution at this point
244  (*exact_soln_pt)(x, exact_soln);
245 
246  // Output x,y,...
247  for (unsigned i = 0; i < 2; i++)
248  {
249  outfile << x[i] << " ";
250  }
251 
252  // Output "exact solution"
253  for (unsigned i = 0; i < exact_soln.size(); i++)
254  {
255  outfile << exact_soln[i] << " ";
256  }
257 
258  outfile << std::endl;
259  }
260 
261  // Write tecplot footer (e.g. FE connectivity lists)
262  write_tecplot_zone_footer(outfile, nplot);
263  }

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

331  {
332  // Find number of nodes
333  unsigned n_node = nnode();
334 
335  // Local shape function
336  Shape psi(n_node);
337 
338  // Vectors of local coordinates and coords and velocities
339  Vector<double> s(2);
340  Vector<double> interpolated_x(2);
341  Vector<double> interpolated_u(2);
342 
343  // Tecplot header info
344  outfile << tecplot_zone_string(nplot);
345 
346  // Loop over plot points
347  unsigned num_plot_points = nplot_points(nplot);
348  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
349  {
350  // Get local coordinates of plot point
351  get_s_plot(iplot, nplot, s);
352 
353  // Get shape functions
354  shape(s, psi);
355 
356  // Loop over directions
357  for (unsigned i = 0; i < 2; i++)
358  {
359  interpolated_x[i] = 0.0;
360  interpolated_u[i] = 0.0;
361  // Loop over the local nodes and sum
362  for (unsigned l = 0; l < n_node; l++)
363  {
364  interpolated_u[i] += u_pnst(t, l, i) * psi[l];
365  interpolated_x[i] += nodal_position(t, l, i) * psi[l];
366  }
367  }
368 
369  // Coordinates
370  for (unsigned i = 0; i < 2; i++)
371  {
372  outfile << interpolated_x[i] << " ";
373  }
374 
375  // Velocities
376  for (unsigned i = 0; i < 2; i++)
377  {
378  outfile << interpolated_u[i] << " ";
379  }
380 
381  outfile << std::endl;
382  }
383  // Write tecplot footer (e.g. FE connectivity lists)
384  write_tecplot_zone_footer(outfile, nplot);
385  }

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

◆ p_local_eqn()

virtual int oomph::PolarNavierStokesEquations::p_local_eqn ( const unsigned &  n)
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::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

Referenced by fill_in_generic_residual_contribution(), and oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution().

◆ p_nodal_index_pnst()

virtual int oomph::PolarNavierStokesEquations::p_nodal_index_pnst ( )
inlinevirtual

Which nodal value represents the pressure? (Default: negative, indicating that pressure is not based on nodal interpolation).

Reimplemented in oomph::PolarTaylorHoodElement.

433  {
434  return Pressure_not_stored_at_node;
435  }
oomph::PolarNavierStokesEquations::Pressure_not_stored_at_node
static int Pressure_not_stored_at_node
Definition: polar_navier_stokes_elements.h:124

References Pressure_not_stored_at_node.

Referenced by oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution().

◆ p_pnst()

virtual double oomph::PolarNavierStokesEquations::p_pnst ( const unsigned &  n_p) const
pure virtual

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

Implemented in oomph::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

Referenced by fill_in_generic_residual_contribution(), oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution(), and interpolated_p_pnst().

◆ pressure_integral()

double oomph::PolarNavierStokesEquations::pressure_integral ( ) const

Integral of pressure over element.

Return pressure integrated over the element.

1049  {
1050  // Initialise
1051  double press_int = 0;
1052 
1053  // Set the value of n_intpt
1054  unsigned n_intpt = integral_pt()->nweight();
1055 
1056  // Set the Vector to hold local coordinates
1057  Vector<double> s(2);
1058 
1059  // Loop over the integration points
1060  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1061  {
1062  // Assign values of s
1063  for (unsigned i = 0; i < 2; i++)
1064  {
1065  s[i] = integral_pt()->knot(ipt, i);
1066  }
1067 
1068  // Get the integral weight
1069  double w = integral_pt()->weight(ipt);
1070 
1071  // Get Jacobian of mapping
1072  double J = J_eulerian(s);
1073 
1074  // Premultiply the weights and the Jacobian
1075  double W = w * J;
1076 
1077  // Get pressure
1078  double press = interpolated_p_pnst(s);
1079 
1080  // Add
1081  press_int += press * W;
1082  }
1083 
1084  return press_int;
1085  }

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

◆ pressure_node_pt()

virtual Node* oomph::PolarNavierStokesEquations::pressure_node_pt ( const unsigned &  n_p)
inlinevirtual

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

Reimplemented in oomph::RefineablePolarTaylorHoodElement, oomph::RefineablePolarNavierStokesEquations, and oomph::PolarTaylorHoodElement.

440  {
441  return NULL;
442  }

◆ pshape_pnst() [1/2]

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

Compute the pressure shape functions at local coordinate s.

Implemented in oomph::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

Referenced by fill_in_generic_residual_contribution(), oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution(), and interpolated_p_pnst().

◆ pshape_pnst() [2/2]

virtual void oomph::PolarNavierStokesEquations::pshape_pnst ( 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::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

◆ re()

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

Reynolds number.

267  {
268  return *Re_pt;
269  }

References Re_pt.

Referenced by fill_in_generic_residual_contribution(), and oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution().

◆ re_invfr()

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

Global inverse Froude number.

329  {
330  return *ReInvFr_pt;
331  }

References ReInvFr_pt.

◆ re_invfr_pt()

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

Pointer to global inverse Froude number.

335  {
336  return ReInvFr_pt;
337  }

References ReInvFr_pt.

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

◆ re_pt()

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

Pointer to Reynolds number.

285  {
286  return Re_pt;
287  }

References Re_pt.

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

◆ re_st()

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

Product of Reynolds and Strouhal number (=Womersley number)

279  {
280  return *ReSt_pt;
281  }

References ReSt_pt.

Referenced by fill_in_generic_residual_contribution(), and oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution().

◆ re_st_pt()

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

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

297  {
298  return ReSt_pt;
299  }

References ReSt_pt.

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

◆ source_fct_pt() [1/2]

NavierStokesSourceFctPt& oomph::PolarNavierStokesEquations::source_fct_pt ( )
inline

Access function for the source-function pointer.

365  {
366  return Source_fct_pt;
367  }

References Source_fct_pt.

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

◆ source_fct_pt() [2/2]

NavierStokesSourceFctPt oomph::PolarNavierStokesEquations::source_fct_pt ( ) const
inline

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

371  {
372  return Source_fct_pt;
373  }

References Source_fct_pt.

◆ strain_rate()

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

Strain-rate tensor: Now returns polar strain.

Get strain-rate tensor: Slightly more complicated in polar coordinates (see eg. Aris)

824  {
825 #ifdef PARANOID
826  if ((strainrate.ncol() != 2) || (strainrate.nrow() != 2))
827  {
828  std::ostringstream error_stream;
829  error_stream << "Wrong size " << strainrate.ncol() << " "
830  << strainrate.nrow() << std::endl;
831  throw OomphLibError(
832  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
833  }
834 #endif
835 
836  // Velocity gradient matrix
837  DenseMatrix<double> dudx(2);
838 
839  // Find out how many nodes there are in the element
840  unsigned n_node = nnode();
841 
842  // Set up memory for the shape and test functions
843  Shape psif(n_node);
844  DShape dpsifdx(n_node, 2);
845 
846  // Call the derivatives of the shape functions
847  dshape_eulerian(s, psif, dpsifdx);
848 
849  // Find the indices at which the local velocities are stored
850  unsigned u_nodal_index[2];
851  for (unsigned i = 0; i < 2; i++)
852  {
853  u_nodal_index[i] = u_index_pnst(i);
854  }
855 
856  // Calculate local values of the velocity components
857  // Allocate
858  Vector<double> interpolated_u(2);
859  Vector<double> interpolated_x(2);
860  DenseMatrix<double> interpolated_dudx(2, 2);
861 
862  // Initialise to zero
863  for (unsigned i = 0; i < 2; i++)
864  {
865  interpolated_u[i] = 0.0;
866  interpolated_x[i] = 0.0;
867  for (unsigned j = 0; j < 2; j++)
868  {
869  interpolated_dudx(i, j) = 0.0;
870  }
871  }
872 
873  // Calculate velocities and derivatives:
874  // Loop over nodes
875  for (unsigned l = 0; l < n_node; l++)
876  {
877  // Loop over directions
878  for (unsigned i = 0; i < 2; i++)
879  {
880  // Get the nodal value
881  double u_value = this->nodal_value(l, u_nodal_index[i]);
882  interpolated_u[i] += u_value * psif[l];
883  interpolated_x[i] += this->nodal_position(l, i) * psif[l];
884 
885  // Loop over derivative directions
886  for (unsigned j = 0; j < 2; j++)
887  {
888  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
889  }
890  }
891  }
892 
893  const double Alpha = alpha();
894  // Can no longer loop over strain components
895  strainrate(0, 0) = interpolated_dudx(0, 0);
896  strainrate(1, 1) =
897  (1. / (Alpha * interpolated_x[0])) * interpolated_dudx(1, 1) +
898  (interpolated_u[0] / interpolated_x[0]);
899  strainrate(1, 0) =
900  0.5 * (interpolated_dudx(1, 0) - (interpolated_u[1] / interpolated_x[0]) +
901  (1. / (Alpha * interpolated_x[0])) * interpolated_dudx(0, 1));
902  strainrate(0, 1) = strainrate(1, 0);
903  }
oomph::DenseMatrix::nrow
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:485
oomph::DenseMatrix::ncol
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:491
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References TanhSolnForAdvectionDiffusion::Alpha, alpha(), oomph::FiniteElement::dshape_eulerian(), i, oomph::FiniteElement::interpolated_x(), j, oomph::DenseMatrix< T >::ncol(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, and u_index_pnst().

Referenced by dissipation(), and get_traction().

◆ strain_rate_by_r()

void oomph::PolarNavierStokesEquations::strain_rate_by_r ( const Vector< double > &  s,
DenseMatrix< double > &  strainrate 
) const

Function to return polar strain multiplied by r.

Return polar strain-rate tensor multiplied by r Slightly more complicated in polar coordinates (see eg. Aris)

911  {
912 #ifdef PARANOID
913  if ((strainrate.ncol() != 2) || (strainrate.nrow() != 2))
914  {
915  std::ostringstream error_stream;
916  error_stream << "Wrong size " << strainrate.ncol() << " "
917  << strainrate.nrow() << std::endl;
918  throw OomphLibError(
919  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
920  }
921 #endif
922 
923  // Velocity gradient matrix
924  DenseMatrix<double> dudx(2);
925 
926  // Find out how many nodes there are in the element
927  unsigned n_node = nnode();
928 
929  // Set up memory for the shape and test functions
930  Shape psif(n_node);
931  DShape dpsifdx(n_node, 2);
932 
933  // Call the derivatives of the shape functions
934  dshape_eulerian(s, psif, dpsifdx);
935 
936  // Find the indices at which the local velocities are stored
937  unsigned u_nodal_index[2];
938  for (unsigned i = 0; i < 2; i++)
939  {
940  u_nodal_index[i] = u_index_pnst(i);
941  }
942 
943  // Calculate local values of the velocity components
944  // Allocate
945  Vector<double> interpolated_u(2);
946  Vector<double> interpolated_x(2);
947  DenseMatrix<double> interpolated_dudx(2, 2);
948 
949  // Initialise to zero
950  for (unsigned i = 0; i < 2; i++)
951  {
952  interpolated_u[i] = 0.0;
953  interpolated_x[i] = 0.0;
954  for (unsigned j = 0; j < 2; j++)
955  {
956  interpolated_dudx(i, j) = 0.0;
957  }
958  }
959 
960  // Calculate velocities and derivatives:
961  // Loop over nodes
962  for (unsigned l = 0; l < n_node; l++)
963  {
964  // Loop over directions
965  for (unsigned i = 0; i < 2; i++)
966  {
967  // Get the nodal value
968  double u_value = this->nodal_value(l, u_nodal_index[i]);
969  interpolated_u[i] += u_value * psif[l];
970  interpolated_x[i] += this->nodal_position(l, i) * psif[l];
971 
972  // Loop over derivative directions
973  for (unsigned j = 0; j < 2; j++)
974  {
975  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
976  }
977  }
978  }
979 
980  const double Alpha = alpha();
981  // Can no longer loop over strain components
982  strainrate(0, 0) = interpolated_dudx(0, 0);
983  strainrate(1, 1) =
984  (1. / (Alpha * interpolated_x[0])) * interpolated_dudx(1, 1) +
985  (interpolated_u[0] / interpolated_x[0]);
986  strainrate(1, 0) =
987  0.5 * (interpolated_dudx(1, 0) - (interpolated_u[1] / interpolated_x[0]) +
988  (1. / (Alpha * interpolated_x[0])) * interpolated_dudx(0, 1));
989  strainrate(0, 1) = strainrate(1, 0);
990 
991  strainrate(0, 0) *= interpolated_x[0];
992  strainrate(1, 1) *= interpolated_x[0];
993  strainrate(1, 0) *= interpolated_x[0];
994  strainrate(0, 1) *= interpolated_x[0];
995  /*
996  strainrate(0,0)*=interpolated_x[0];
997  strainrate(1,1)*=interpolated_x[0];
998  strainrate(1,0)*=interpolated_x[0];
999  strainrate(0,1)*=interpolated_x[0];
1000  */
1001  }

References TanhSolnForAdvectionDiffusion::Alpha, alpha(), oomph::FiniteElement::dshape_eulerian(), i, oomph::FiniteElement::interpolated_x(), j, oomph::DenseMatrix< T >::ncol(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, and u_index_pnst().

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

◆ u_index_pnst()

virtual unsigned oomph::PolarNavierStokesEquations::u_index_pnst ( 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.

395  {
396  return i;
397  }

References i.

Referenced by fill_in_generic_residual_contribution(), oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution(), oomph::RefineablePolarTaylorHoodElement::get_interpolated_values(), oomph::RefineablePolarCrouzeixRaviartElement::get_interpolated_values(), oomph::RefineablePolarTaylorHoodElement::insert_load_data(), oomph::RefineablePolarCrouzeixRaviartElement::insert_load_data(), strain_rate(), and strain_rate_by_r().

◆ u_pnst() [1/2]

virtual double oomph::PolarNavierStokesEquations::u_pnst ( const unsigned &  n,
const unsigned &  i 
) const
pure virtual

Velocity i at local node n. Uses suitably interpolated value for hanging nodes.

Implemented in oomph::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

Referenced by du_dt_pnst(), full_output(), interpolated_dudx_pnst(), interpolated_u_pnst(), and output_veloc().

◆ u_pnst() [2/2]

virtual double oomph::PolarNavierStokesEquations::u_pnst ( const unsigned &  t,
const unsigned &  n,
const unsigned &  i 
) const
pure virtual

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

Implemented in oomph::PolarTaylorHoodElement, and oomph::PolarCrouzeixRaviartElement.

◆ viscosity_ratio()

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

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

304  {
305  return *Viscosity_Ratio_pt;
306  }

References Viscosity_Ratio_pt.

◆ viscosity_ratio_pt()

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

Pointer to Viscosity Ratio.

310  {
311  return Viscosity_Ratio_pt;
312  }

References Viscosity_Ratio_pt.

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

Member Data Documentation

◆ Alpha_pt

double* oomph::PolarNavierStokesEquations::Alpha_pt
protected

Pointer to the angle alpha.

Referenced by alpha(), alpha_pt(), and oomph::RefineablePolarNavierStokesEquations::further_build().

◆ Body_force_fct_pt

NavierStokesBodyForceFctPt oomph::PolarNavierStokesEquations::Body_force_fct_pt
protected

Pointer to body force function.

Referenced by body_force_fct_pt(), oomph::RefineablePolarNavierStokesEquations::further_build(), and get_body_force().

◆ Default_Gravity_vector

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

Static default value for the gravity vector.

Navier-Stokes equations default gravity vector.

Referenced by PolarNavierStokesEquations().

◆ Default_Physical_Constant_Value

double oomph::PolarNavierStokesEquations::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 PolarNavierStokesEquations().

◆ Default_Physical_Ratio_Value

double oomph::PolarNavierStokesEquations::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 PolarNavierStokesEquations().

◆ Density_Ratio_pt

double* oomph::PolarNavierStokesEquations::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::RefineablePolarNavierStokesEquations::further_build(), and PolarNavierStokesEquations().

◆ G_pt

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

Pointer to global gravity Vector.

Referenced by oomph::RefineablePolarNavierStokesEquations::further_build(), g(), g_pt(), and PolarNavierStokesEquations().

◆ Gamma

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

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

Start of what would've been navier_stokes_elements.cc //.

Navier–Stokes equations static data

Referenced by fill_in_generic_residual_contribution(), and oomph::RefineablePolarNavierStokesEquations::fill_in_generic_residual_contribution().

◆ Pressure_not_stored_at_node

int oomph::PolarNavierStokesEquations::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_pnst().

◆ Re_pt

double* oomph::PolarNavierStokesEquations::Re_pt
protected

Pointer to global Reynolds number.

Referenced by oomph::RefineablePolarNavierStokesEquations::further_build(), PolarNavierStokesEquations(), re(), and re_pt().

◆ ReInvFr_pt

double* oomph::PolarNavierStokesEquations::ReInvFr_pt
protected

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

Referenced by oomph::RefineablePolarNavierStokesEquations::further_build(), PolarNavierStokesEquations(), re_invfr(), and re_invfr_pt().

◆ ReSt_pt

double* oomph::PolarNavierStokesEquations::ReSt_pt
protected

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

Referenced by oomph::RefineablePolarNavierStokesEquations::further_build(), PolarNavierStokesEquations(), re_st(), and re_st_pt().

◆ Source_fct_pt

NavierStokesSourceFctPt oomph::PolarNavierStokesEquations::Source_fct_pt
protected

Pointer to volumetric source function.

Referenced by oomph::RefineablePolarNavierStokesEquations::further_build(), get_source_fct(), and source_fct_pt().

◆ Viscosity_Ratio_pt

double* oomph::PolarNavierStokesEquations::Viscosity_Ratio_pt
protected

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

Referenced by oomph::RefineablePolarNavierStokesEquations::further_build(), PolarNavierStokesEquations(), viscosity_ratio(), and viscosity_ratio_pt().

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