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

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

#include <refineable_axisym_navier_stokes_elements.h>

+ Inheritance diagram for oomph::RefineableAxisymmetricNavierStokesEquations:

Public Member Functions

 RefineableAxisymmetricNavierStokesEquations ()
 Empty Constructor. More...
 
unsigned num_Z2_flux_terms ()
 Number of 'flux' terms for Z2 error estimation. More...
 
void get_Z2_flux (const Vector< double > &s, Vector< double > &flux)
 
double geometric_jacobian (const Vector< double > &x)
 Fill in the geometric Jacobian, which in this case is r. More...
 
void further_build ()
 Further build: pass the pointers down to the sons. More...
 
void dinterpolated_u_axi_nst_ddata (const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
 
virtual void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
- Public Member Functions inherited from oomph::AxisymmetricNavierStokesEquations
 AxisymmetricNavierStokesEquations ()
 Constructor: NULL the body force and source function. More...
 
const double & re () const
 Reynolds number. More...
 
const double & re_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
const double & re_invfr () const
 Global inverse Froude number. More...
 
double *& re_invfr_pt ()
 Pointer to global inverse Froude number. More...
 
const double & re_invro () const
 Global Reynolds number multiplied by inverse Rossby number. More...
 
double *& re_invro_pt ()
 Pointer to global inverse Froude number. More...
 
const Vector< double > & g () const
 Vector of gravitational components. More...
 
Vector< double > *& g_pt ()
 Pointer to Vector of gravitational components. More...
 
const double & density_ratio () const
 
double *& density_ratio_pt ()
 Pointer to Density ratio. More...
 
const double & viscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to Viscosity Ratio. More...
 
virtual unsigned npres_axi_nst () const =0
 Function to return number of pressure degrees of freedom. More...
 
virtual unsigned u_index_axi_nst (const unsigned &i) const
 
unsigned u_index_nst (const unsigned &i) const
 
unsigned n_u_nst () const
 
double du_dt_axi_nst (const unsigned &n, const unsigned &i) const
 
void disable_ALE ()
 
void enable_ALE ()
 
virtual double p_axi_nst (const unsigned &n_p) const =0
 
virtual int p_nodal_index_axi_nst () const
 Which nodal value represents the pressure? More...
 
double pressure_integral () const
 Integral of pressure over element. More...
 
double dissipation () const
 Return integral of dissipation over element. More...
 
double dissipation (const Vector< double > &s) const
 Return dissipation at local coordinate s. More...
 
double kin_energy () const
 Get integral of kinetic energy over element. More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 Strain-rate tensor: \( e_{ij} \) where \( i,j = r,z,\theta \) (in that order) More...
 
void traction (const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
 
void get_pressure_and_velocity_mass_matrix_diagonal (Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
 
unsigned nscalar_paraview () const
 
void scalar_value_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
 
std::string scalar_name_paraview (const unsigned &i) const
 
void point_output_data (const Vector< double > &s, Vector< double > &data)
 
void output (std::ostream &outfile)
 
void output (std::ostream &outfile, const unsigned &nplot)
 
void output (FILE *file_pt)
 
void output (FILE *file_pt, const unsigned &nplot)
 
void output_veloc (std::ostream &outfile, const unsigned &nplot, const unsigned &t)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 
void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
void fill_in_contribution_to_residuals (Vector< double > &residuals)
 Compute the element's residual Vector. More...
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void fill_in_contribution_to_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
void fill_in_contribution_to_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 Compute the element's residual Vector. More...
 
void fill_in_contribution_to_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
void interpolated_u_axi_nst (const Vector< double > &s, Vector< double > &veloc) const
 Compute vector of FE interpolated velocity u at local coordinate s. More...
 
double interpolated_u_axi_nst (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated velocity u[i] at local coordinate s. More...
 
double interpolated_u_axi_nst (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 Return FE interpolated velocity u[i] at local coordinate s. More...
 
double interpolated_p_axi_nst (const Vector< double > &s) const
 Return FE interpolated pressure at local coordinate s. More...
 
double interpolated_duds_axi_nst (const Vector< double > &s, const unsigned &i, const unsigned &j) const
 
double interpolated_dudx_axi_nst (const Vector< double > &s, const unsigned &i, const unsigned &j) const
 
double interpolated_dudt_axi_nst (const Vector< double > &s, const unsigned &i) const
 
double interpolated_d_dudx_dX_axi_nst (const Vector< double > &s, const unsigned &p, const unsigned &q, const unsigned &i, const unsigned &k) const
 
double compute_physical_size () const
 Compute the volume of the element. More...
 
- Public Member Functions inherited from oomph::FiniteElement
void set_dimension (const unsigned &dim)
 
void set_nodal_dimension (const unsigned &nodal_dim)
 
void set_nnodal_position_type (const unsigned &nposition_type)
 Set the number of types required to interpolate the coordinate. More...
 
void set_n_node (const unsigned &n)
 
int nodal_local_eqn (const unsigned &n, const unsigned &i) const
 
double dJ_eulerian_at_knot (const unsigned &ipt, Shape &psi, DenseMatrix< double > &djacobian_dX) const
 
 FiniteElement ()
 Constructor. More...
 
virtual ~FiniteElement ()
 
 FiniteElement (const FiniteElement &)=delete
 Broken copy constructor. More...
 
virtual bool local_coord_is_valid (const Vector< double > &s)
 Broken assignment operator. More...
 
virtual void move_local_coord_back_into_element (Vector< double > &s) const
 
void get_centre_of_gravity_and_max_radius_in_terms_of_zeta (Vector< double > &cog, double &max_radius) const
 
virtual void local_coordinate_of_node (const unsigned &j, Vector< double > &s) const
 
virtual void local_fraction_of_node (const unsigned &j, Vector< double > &s_fraction)
 
virtual double local_one_d_fraction_of_node (const unsigned &n1d, const unsigned &i)
 
virtual void set_macro_elem_pt (MacroElement *macro_elem_pt)
 
MacroElement * macro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void describe_nodal_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void assign_all_generic_local_eqn_numbers (const bool &store_local_dof_pt)
 
Node *& node_pt (const unsigned &n)
 Return a pointer to the local node n. More...
 
Node *const & node_pt (const unsigned &n) const
 Return a pointer to the local node n (const version) More...
 
unsigned nnode () const
 Return the number of nodes. More...
 
virtual unsigned nnode_1d () const
 
double raw_nodal_position (const unsigned &n, const unsigned &i) const
 
double raw_nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position (const unsigned &n, const unsigned &i) const
 
double nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double dnodal_position_dt (const unsigned &n, const unsigned &i) const
 Return the i-th component of nodal velocity: dx/dt at local node n. More...
 
double dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual unsigned required_nvalue (const unsigned &n) const
 
unsigned nnodal_position_type () const
 
bool has_hanging_nodes () const
 
unsigned nodal_dimension () const
 Return the required Eulerian dimension of the nodes in this element. More...
 
virtual 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_geometric_data (std::set< Data * > &geometric_data_pt)
 
virtual double s_min () const
 Min value of local coordinate. More...
 
virtual double s_max () const
 Max. value of local coordinate. More...
 
double size () const
 
void point_output (std::ostream &outfile, const Vector< double > &s)
 
virtual unsigned nplot_points_paraview (const unsigned &nplot) const
 
virtual unsigned nsub_elements_paraview (const unsigned &nplot) const
 
void output_paraview (std::ofstream &file_out, const unsigned &nplot) const
 
virtual void write_paraview_output_offset_information (std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
 
virtual void write_paraview_type (std::ofstream &file_out, const unsigned &nplot) const
 
virtual void write_paraview_offsets (std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
 
virtual void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
 
virtual void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
virtual void output (const unsigned &t, std::ostream &outfile, const unsigned &n_plot) const
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, const SolutionFunctorBase &exact_soln) const
 Output a time-dependent exact solution over the element. More...
 
virtual void get_s_plot (const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
 
virtual std::string tecplot_zone_string (const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (std::ostream &outfile, const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (FILE *file_pt, const unsigned &nplot) const
 
virtual unsigned nplot_points (const unsigned &nplot) const
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_abs_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error)
 
void integrate_fct (FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt, Vector< double > &integral)
 Integrate Vector-valued function over element. More...
 
void integrate_fct (FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt, const double &time, Vector< double > &integral)
 Integrate Vector-valued time-dep function over element. More...
 
virtual void build_face_element (const int &face_index, FaceElement *face_element_pt)
 
virtual unsigned self_test ()
 
virtual unsigned get_bulk_node_number (const int &face_index, const unsigned &i) const
 
virtual int face_outer_unit_normal_sign (const int &face_index) const
 Get the sign of the outer unit normal on the face given by face_index. More...
 
virtual unsigned nnode_on_face () const
 
void face_node_number_error_check (const unsigned &i) const
 Range check for face node numbers. More...
 
virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt (const int &face_index) const
 
virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt (const int &face_index) const
 
- Public Member Functions inherited from oomph::GeneralisedElement
 GeneralisedElement ()
 Constructor: Initialise all pointers and all values to zero. More...
 
virtual ~GeneralisedElement ()
 Virtual destructor to clean up any memory allocated by the object. More...
 
 GeneralisedElement (const GeneralisedElement &)=delete
 Broken copy constructor. More...
 
void operator= (const GeneralisedElement &)=delete
 Broken assignment operator. More...
 
Data *& internal_data_pt (const unsigned &i)
 Return a pointer to i-th internal data object. More...
 
Data *const & internal_data_pt (const unsigned &i) const
 Return a pointer to i-th internal data object (const version) More...
 
Data *& external_data_pt (const unsigned &i)
 Return a pointer to i-th external data object. More...
 
Data *const & external_data_pt (const unsigned &i) const
 Return a pointer to i-th external data object (const version) More...
 
unsigned long eqn_number (const unsigned &ieqn_local) const
 
int local_eqn_number (const unsigned long &ieqn_global) const
 
unsigned add_external_data (Data *const &data_pt, const bool &fd=true)
 
bool external_data_fd (const unsigned &i) const
 
void exclude_external_data_fd (const unsigned &i)
 
void include_external_data_fd (const unsigned &i)
 
void flush_external_data ()
 Flush all external data. More...
 
void flush_external_data (Data *const &data_pt)
 Flush the object addressed by data_pt from the external data array. More...
 
unsigned ninternal_data () const
 Return the number of internal data objects. More...
 
unsigned nexternal_data () const
 Return the number of external data objects. More...
 
unsigned ndof () const
 Return the number of equations/dofs in the element. More...
 
void dof_vector (const unsigned &t, Vector< double > &dof)
 Return the vector of dof values at time level t. More...
 
void dof_pt_vector (Vector< double * > &dof_pt)
 Return the vector of pointers to dof values. More...
 
void set_internal_data_time_stepper (const unsigned &i, TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void assign_internal_eqn_numbers (unsigned long &global_number, Vector< double * > &Dof_pt)
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void add_internal_value_pt_to_map (std::map< unsigned, double * > &map_of_value_pt)
 
virtual void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void complete_setup_of_dependencies ()
 
virtual void get_residuals (Vector< double > &residuals)
 
virtual void get_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
virtual void get_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void get_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void get_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void get_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void get_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void get_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void get_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void get_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_norm (double &norm)
 
virtual unsigned ndof_types () const
 
virtual void get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
 
- Public Member Functions inherited from oomph::GeomObject
 GeomObject ()
 Default constructor. More...
 
 GeomObject (const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)
 
 GeomObject (const GeomObject &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const GeomObject &)=delete
 Broken assignment operator. More...
 
virtual ~GeomObject ()
 (Empty) destructor More...
 
unsigned nlagrangian () const
 Access function to # of Lagrangian coordinates. More...
 
unsigned ndim () const
 Access function to # of Eulerian coordinates. More...
 
void set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
 Set # of Lagrangian and Eulerian coordinates. More...
 
TimeStepper *& time_stepper_pt ()
 
TimeStepper * time_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 
- Public Member Functions inherited from oomph::NavierStokesElementWithDiagonalMassMatrices
 NavierStokesElementWithDiagonalMassMatrices ()
 Empty constructor. More...
 
virtual ~NavierStokesElementWithDiagonalMassMatrices ()
 Virtual destructor. More...
 
 NavierStokesElementWithDiagonalMassMatrices (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken copy constructor. More...
 
void operator= (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken assignment operator. More...
 
- Public Member Functions inherited from oomph::RefineableElement
 RefineableElement ()
 
virtual ~RefineableElement ()
 Destructor, delete the allocated storage for the hanging equations. More...
 
 RefineableElement (const RefineableElement &)=delete
 Broken copy constructor. More...
 
void operator= (const RefineableElement &)=delete
 Broken assignment operator. More...
 
Tree * tree_pt ()
 Access function: Pointer to quadtree representation of this element. More...
 
void set_tree_pt (Tree *my_tree_pt)
 Set pointer to quadtree representation of this element. More...
 
virtual unsigned required_nsons () const
 
bool refinement_is_enabled ()
 Flag to indicate suppression of any refinement. More...
 
void disable_refinement ()
 Suppress of any refinement for this element. More...
 
void enable_refinement ()
 Emnable refinement for this element. More...
 
template<class ELEMENT >
void split (Vector< ELEMENT * > &son_pt) const
 
int local_hang_eqn (Node *const &node_pt, const unsigned &i)
 
virtual void build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)=0
 
void set_refinement_level (const int &refine_level)
 Set the refinement level. More...
 
unsigned refinement_level () const
 Return the Refinement level. More...
 
void select_for_refinement ()
 Select the element for refinement. More...
 
void deselect_for_refinement ()
 Deselect the element for refinement. More...
 
void select_sons_for_unrefinement ()
 Unrefinement will be performed by merging the four sons of this element. More...
 
void deselect_sons_for_unrefinement ()
 
bool to_be_refined ()
 Has the element been selected for refinement? More...
 
bool sons_to_be_unrefined ()
 Has the element been selected for unrefinement? More...
 
virtual void rebuild_from_sons (Mesh *&mesh_pt)=0
 
virtual void unbuild ()
 
virtual void deactivate_element ()
 
virtual bool nodes_built ()
 Return true if all the nodes have been built, false if not. More...
 
long number () const
 Element number (for debugging/plotting) More...
 
void set_number (const long &mynumber)
 Set element number (for debugging/plotting) More...
 
virtual unsigned ncont_interpolated_values () const =0
 
virtual void get_interpolated_values (const Vector< double > &s, Vector< double > &values)
 
virtual void get_interpolated_values (const unsigned &t, const Vector< double > &s, Vector< double > &values)=0
 
virtual Node * interpolating_node_pt (const unsigned &n, const int &value_id)
 
virtual double local_one_d_fraction_of_interpolating_node (const unsigned &n1d, const unsigned &i, const int &value_id)
 
virtual Node * get_interpolating_node_at_local_coordinate (const Vector< double > &s, const int &value_id)
 
virtual unsigned ninterpolating_node (const int &value_id)
 
virtual unsigned ninterpolating_node_1d (const int &value_id)
 
virtual void interpolating_basis (const Vector< double > &s, Shape &psi, const int &value_id) const
 
virtual void check_integrity (double &max_error)=0
 
void identify_field_data_for_interactions (std::set< std::pair< Data *, unsigned >> &paired_field_data)
 
void assign_nodal_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual RefineableElement * root_element_pt ()
 
virtual RefineableElement * father_element_pt () const
 Return a pointer to the father element. More...
 
void get_father_at_refinement_level (unsigned &refinement_level, RefineableElement *&father_at_reflevel_pt)
 
virtual void initial_setup (Tree *const &adopted_father_pt=0, const unsigned &initial_p_order=0)
 
virtual void pre_build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt)
 Pre-build the element. More...
 
virtual void setup_hanging_nodes (Vector< std::ofstream * > &output_stream)
 
virtual void further_setup_hanging_nodes ()
 
void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
unsigned nshape_controlling_nodes ()
 
std::map< Node *, unsigned > shape_controlling_node_lookup ()
 
- Public Member Functions inherited from oomph::ElementWithZ2ErrorEstimator
 ElementWithZ2ErrorEstimator ()
 Default empty constructor. More...
 
 ElementWithZ2ErrorEstimator (const ElementWithZ2ErrorEstimator &)=delete
 Broken copy constructor. More...
 
void operator= (const ElementWithZ2ErrorEstimator &)=delete
 Broken assignment operator. More...
 
virtual unsigned ncompound_fluxes ()
 
virtual void compute_exact_Z2_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
 
virtual void get_Z2_compound_flux_indices (Vector< unsigned > &flux_index)
 
virtual unsigned nvertex_node () const =0
 Number of vertex nodes in the element. More...
 
virtual Node * vertex_node_pt (const unsigned &j) const =0
 
virtual unsigned nrecovery_order ()=0
 Order of recovery shape functions. More...
 

Static Public Member Functions

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

Protected Member Functions

virtual Node * pressure_node_pt (const unsigned &n_p)
 
virtual void unpin_elemental_pressure_dofs ()=0
 Unpin all pressure dofs in the element. More...
 
virtual void pin_elemental_redundant_nodal_pressure_dofs ()
 
- Protected Member Functions inherited from oomph::AxisymmetricNavierStokesEquations
virtual int p_local_eqn (const unsigned &n) const =0
 
virtual double dshape_and_dtest_eulerian_axi_nst (const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_axi_nst (const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_axi_nst (const unsigned &ipt, Shape &psi, DShape &dpsidx, RankFourTensor< double > &d_dpsidx_dX, Shape &test, DShape &dtestdx, RankFourTensor< double > &d_dtestdx_dX, DenseMatrix< double > &djacobian_dX) const =0
 
virtual void pshape_axi_nst (const Vector< double > &s, Shape &psi) const =0
 Compute the pressure shape functions at local coordinate s. More...
 
virtual void pshape_axi_nst (const Vector< double > &s, Shape &psi, Shape &test) const =0
 
virtual void get_body_force_axi_nst (const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
 Calculate the body force fct at a given time and Eulerian position. More...
 
virtual void get_body_force_gradient_axi_nst (const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
 
double get_source_fct (const double &time, const unsigned &ipt, const Vector< double > &x)
 Calculate the source fct at given time and Eulerian position. More...
 
virtual void get_source_fct_gradient (const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
 
virtual void get_viscosity_ratio_axisym_nst (const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, double &visc_ratio)
 
- Protected Member Functions inherited from oomph::FiniteElement
template<unsigned DIM>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double invert_jacobian_mapping (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &inverse_jacobian) const
 
virtual void dJ_eulerian_dnodal_coordinates (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<unsigned DIM>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
virtual void d_dshape_eulerian_dnodal_coordinates (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<unsigned DIM>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
virtual void transform_derivatives (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
void transform_derivatives_diagonal (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
virtual void transform_second_derivatives (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
void fill_in_jacobian_from_nodal_by_fd (DenseMatrix< double > &jacobian)
 
virtual void update_before_nodal_fd ()
 
virtual void reset_after_nodal_fd ()
 
virtual void update_in_nodal_fd (const unsigned &i)
 
virtual void reset_in_nodal_fd (const unsigned &i)
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Zero-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 One-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Two-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
- Protected Member Functions inherited from oomph::GeneralisedElement
unsigned add_internal_data (Data *const &data_pt, const bool &fd=true)
 
bool internal_data_fd (const unsigned &i) const
 
void exclude_internal_data_fd (const unsigned &i)
 
void include_internal_data_fd (const unsigned &i)
 
void clear_global_eqn_numbers ()
 
void add_global_eqn_numbers (std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)
 
virtual void assign_internal_and_external_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void assign_additional_local_eqn_numbers ()
 
int internal_local_eqn (const unsigned &i, const unsigned &j) const
 
int external_local_eqn (const unsigned &i, const unsigned &j)
 
void fill_in_jacobian_from_internal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_internal_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
virtual void update_before_internal_fd ()
 
virtual void reset_after_internal_fd ()
 
virtual void update_in_internal_fd (const unsigned &i)
 
virtual void reset_in_internal_fd (const unsigned &i)
 
virtual void update_before_external_fd ()
 
virtual void reset_after_external_fd ()
 
virtual void update_in_external_fd (const unsigned &i)
 
virtual void reset_in_external_fd (const unsigned &i)
 
virtual void fill_in_contribution_to_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void fill_in_contribution_to_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void fill_in_contribution_to_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
- Protected Member Functions inherited from oomph::RefineableElement
void assemble_local_to_eulerian_jacobian (const DShape &dpsids, DenseMatrix< double > &jacobian) const
 
void assemble_local_to_eulerian_jacobian2 (const DShape &d2psids, DenseMatrix< double > &jacobian2) const
 
void assemble_eulerian_base_vectors (const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
 
double local_to_eulerian_mapping_diagonal (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
void assign_hanging_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign the local equation numbers for hanging node variables. More...
 
virtual void fill_in_jacobian_from_nodal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 

Private Member Functions

void fill_in_generic_residual_contribution_axi_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 
void fill_in_generic_dresidual_contribution_axi_nst (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
 
void fill_in_contribution_to_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 

Additional Inherited Members

- Public Types inherited from oomph::FiniteElement
typedef void(* SteadyExactSolutionFctPt) (const Vector< double > &, Vector< double > &)
 
typedef void(* UnsteadyExactSolutionFctPt) (const double &, const Vector< double > &, Vector< double > &)
 
- Public Attributes inherited from oomph::AxisymmetricNavierStokesEquations
void(*&)(const double &time, const Vector< double > &x, Vector< double > &f) axi_nst_body_force_fct_pt ()
 Access function for the body-force pointer. More...
 
double(*&)(const double &time, const Vector< double > &x) source_fct_pt ()
 Access function for the source-function pointer. More...
 
- Static Public Attributes inherited from oomph::AxisymmetricNavierStokesEquations
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
 
- Static Protected Member Functions inherited from oomph::RefineableElement
static void check_value_id (const int &n_continuously_interpolated_values, const int &value_id)
 
- Protected Attributes inherited from oomph::AxisymmetricNavierStokesEquations
double * Viscosity_Ratio_pt
 
double * Density_Ratio_pt
 
double * Re_pt
 Pointer to global Reynolds number. More...
 
double * ReSt_pt
 Pointer to global Reynolds number x Strouhal number (=Womersley) More...
 
double * ReInvFr_pt
 
double * ReInvRo_pt
 
Vector< double > * G_pt
 Pointer to global gravity Vector. More...
 
void(* Body_force_fct_pt )(const double &time, const Vector< double > &x, Vector< double > &result)
 Pointer to body force function. More...
 
double(* Source_fct_pt )(const double &time, const Vector< double > &x)
 Pointer to volumetric source function. More...
 
bool ALE_is_disabled
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElement * Macro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepper * Geom_object_time_stepper_pt
 
- Protected Attributes inherited from oomph::RefineableElement
Tree * Tree_pt
 A pointer to a general tree object. More...
 
unsigned Refine_level
 Refinement level. More...
 
bool To_be_refined
 Flag for refinement. More...
 
bool Refinement_is_enabled
 Flag to indicate suppression of any refinement. More...
 
bool Sons_to_be_unrefined
 Flag for unrefinement. More...
 
long Number
 Global element number – for plotting/validation purposes. More...
 
- Static Protected Attributes inherited from oomph::FiniteElement
static const unsigned Default_Initial_Nvalue = 0
 Default value for the number of values at a node. More...
 
static const double Node_location_tolerance = 1.0e-14
 
static const unsigned N2deriv [] = {0, 1, 3, 6}
 
- Static Protected Attributes inherited from oomph::GeneralisedElement
static DenseMatrix< double > Dummy_matrix
 
static std::deque< double * > Dof_pt_deque
 
- Static Protected Attributes inherited from oomph::RefineableElement
static double Max_integrity_tolerance = 1.0e-8
 Max. allowed discrepancy in element integrity check. More...
 

Detailed Description

Refineable version of the Axisymmetric Navier–Stokes equations.

Constructor & Destructor Documentation

◆ RefineableAxisymmetricNavierStokesEquations()

oomph::RefineableAxisymmetricNavierStokesEquations::RefineableAxisymmetricNavierStokesEquations ( )
inline

Empty Constructor.

68  : AxisymmetricNavierStokesEquations(),
69  RefineableElement(),
70  ElementWithZ2ErrorEstimator()
71  {
72  }
oomph::AxisymmetricNavierStokesEquations::AxisymmetricNavierStokesEquations
AxisymmetricNavierStokesEquations()
Constructor: NULL the body force and source function.
Definition: axisym_navier_stokes_elements.h:375
oomph::ElementWithZ2ErrorEstimator::ElementWithZ2ErrorEstimator
ElementWithZ2ErrorEstimator()
Default empty constructor.
Definition: error_estimator.h:82
oomph::RefineableElement::RefineableElement
RefineableElement()
Definition: refineable_elements.h:188

Member Function Documentation

◆ dinterpolated_u_axi_nst_ddata()

void oomph::RefineableAxisymmetricNavierStokesEquations::dinterpolated_u_axi_nst_ddata ( const Vector< double > &  s,
const unsigned &  i,
Vector< double > &  du_ddata,
Vector< unsigned > &  global_eqn_number 
)
inlinevirtual

Compute the derivatives of the i-th component of velocity at point s with respect to all data that can affect its value. In addition, return the global equation numbers corresponding to the data. Overload the non-refineable version to take account of hanging node information

Reimplemented from oomph::AxisymmetricNavierStokesEquations.

179  {
180  // Find number of nodes
181  unsigned n_node = this->nnode();
182  // Local shape function
183  Shape psi(n_node);
184  // Find values of shape function at the given local coordinate
185  this->shape(s, psi);
186 
187  // Find the index at which the velocity component is stored
188  const unsigned u_nodal_index = this->u_index_axi_nst(i);
189 
190  // Storage for hang info pointer
191  HangInfo* hang_info_pt = 0;
192  // Storage for global equation
193  int global_eqn = 0;
194 
195  // Find the number of dofs associated with interpolated u
196  unsigned n_u_dof = 0;
197  for (unsigned l = 0; l < n_node; l++)
198  {
199  unsigned n_master = 1;
200 
201  // Local bool (is the node hanging)
202  bool is_node_hanging = this->node_pt(l)->is_hanging();
203 
204  // If the node is hanging, get the number of master nodes
205  if (is_node_hanging)
206  {
207  hang_info_pt = this->node_pt(l)->hanging_pt();
208  n_master = hang_info_pt->nmaster();
209  }
210  // Otherwise there is just one master node, the node itself
211  else
212  {
213  n_master = 1;
214  }
215 
216  // Loop over the master nodes
217  for (unsigned m = 0; m < n_master; m++)
218  {
219  // Get the equation number
220  if (is_node_hanging)
221  {
222  // Get the equation number from the master node
223  global_eqn =
224  hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
225  }
226  else
227  {
228  // Global equation number
229  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
230  }
231 
232  // If it's positive add to the count
233  if (global_eqn >= 0)
234  {
235  ++n_u_dof;
236  }
237  }
238  }
239 
240  // Now resize the storage schemes
241  du_ddata.resize(n_u_dof, 0.0);
242  global_eqn_number.resize(n_u_dof, 0);
243 
244  // Loop over th nodes again and set the derivatives
245  unsigned count = 0;
246  // Loop over the local nodes and sum
247  for (unsigned l = 0; l < n_node; l++)
248  {
249  unsigned n_master = 1;
250  double hang_weight = 1.0;
251 
252  // Local bool (is the node hanging)
253  bool is_node_hanging = this->node_pt(l)->is_hanging();
254 
255  // If the node is hanging, get the number of master nodes
256  if (is_node_hanging)
257  {
258  hang_info_pt = this->node_pt(l)->hanging_pt();
259  n_master = hang_info_pt->nmaster();
260  }
261  // Otherwise there is just one master node, the node itself
262  else
263  {
264  n_master = 1;
265  }
266 
267  // Loop over the master nodes
268  for (unsigned m = 0; m < n_master; m++)
269  {
270  // If the node is hanging get weight from master node
271  if (is_node_hanging)
272  {
273  // Get the hang weight from the master node
274  hang_weight = hang_info_pt->master_weight(m);
275  }
276  else
277  {
278  // Node contributes with full weight
279  hang_weight = 1.0;
280  }
281 
282  // Get the equation number
283  if (is_node_hanging)
284  {
285  // Get the equation number from the master node
286  global_eqn =
287  hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
288  }
289  else
290  {
291  // Local equation number
292  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
293  }
294 
295  if (global_eqn >= 0)
296  {
297  // Set the global equation number
298  global_eqn_number[count] = global_eqn;
299  // Set the derivative with respect to the unknown
300  du_ddata[count] = psi[l] * hang_weight;
301  // Increase the counter
302  ++count;
303  }
304  }
305  }
306  }
oomph::AxisymmetricNavierStokesEquations::u_index_axi_nst
virtual unsigned u_index_axi_nst(const unsigned &i) const
Definition: axisym_navier_stokes_elements.h:506
oomph::Data::eqn_number
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::shape
virtual void shape(const Vector< double > &s, Shape &psi) const =0
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::Node::hanging_pt
HangInfo *const & hanging_pt() const
Definition: nodes.h:1228
oomph::Node::is_hanging
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
m
int * m
Definition: level2_cplx_impl.h:294

References oomph::Data::eqn_number(), oomph::Node::hanging_pt(), oomph::Node::is_hanging(), m, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::FiniteElement::shape(), and oomph::AxisymmetricNavierStokesEquations::u_index_axi_nst().

Referenced by RefineableQAxisymAdvectionDiffusionBoussinesqElement::get_dwind_axi_adv_diff_dexternal_element_data().

◆ fill_in_contribution_to_hessian_vector_products()

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

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

Reimplemented from oomph::AxisymmetricNavierStokesEquations.

391  {
392  throw OomphLibError("Not yet implemented\n",
393  OOMPH_CURRENT_FUNCTION,
394  OOMPH_EXCEPTION_LOCATION);
395  }
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 OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ fill_in_generic_dresidual_contribution_axi_nst()

void oomph::RefineableAxisymmetricNavierStokesEquations::fill_in_generic_dresidual_contribution_axi_nst ( double *const &  parameter_pt,
Vector< double > &  dres_dparam,
DenseMatrix< double > &  djac_dparam,
DenseMatrix< double > &  dmass_matrix_dparam,
unsigned  flag 
)
inlineprivatevirtual

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

Reimplemented from oomph::AxisymmetricNavierStokesEquations.

379  {
380  throw OomphLibError("Not yet implemented\n",
381  OOMPH_CURRENT_FUNCTION,
382  OOMPH_EXCEPTION_LOCATION);
383  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ fill_in_generic_residual_contribution_axi_nst()

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

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

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

Reimplemented from oomph::AxisymmetricNavierStokesEquations.

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

References oomph::AxisymmetricNavierStokesEquations::ALE_is_disabled, Global_Parameters::body_force(), oomph::AxisymmetricNavierStokesEquations::density_ratio(), DIM, oomph::FiniteElement::dnodal_position_dt(), oomph::AxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_axi_nst(), oomph::AxisymmetricNavierStokesEquations::du_dt_axi_nst(), G, oomph::AxisymmetricNavierStokesEquations::g(), oomph::AxisymmetricNavierStokesEquations::Gamma, oomph::AxisymmetricNavierStokesEquations::get_body_force_axi_nst(), oomph::AxisymmetricNavierStokesEquations::get_source_fct(), oomph::Node::hanging_pt(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), oomph::Node::is_hanging(), J, j, k, oomph::Integral::knot(), oomph::RefineableElement::local_hang_eqn(), m, m2(), oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), GlobalParameters::Nintpt, oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::FiniteElement::node_pt(), oomph::AxisymmetricNavierStokesEquations::npres_axi_nst(), oomph::Integral::nweight(), oomph::AxisymmetricNavierStokesEquations::p_axi_nst(), oomph::AxisymmetricNavierStokesEquations::p_local_eqn(), oomph::AxisymmetricNavierStokesEquations::p_nodal_index_axi_nst(), pressure_node_pt(), oomph::AxisymmetricNavierStokesEquations::pshape_axi_nst(), UniformPSDSelfTest::r, oomph::AxisymmetricNavierStokesEquations::re(), oomph::AxisymmetricNavierStokesEquations::re_invfr(), oomph::AxisymmetricNavierStokesEquations::re_invro(), oomph::AxisymmetricNavierStokesEquations::re_st(), s, TestProblem::source(), oomph::Time::time(), oomph::TimeStepper::time_pt(), oomph::Data::time_stepper_pt(), oomph::AxisymmetricNavierStokesEquations::u_index_axi_nst(), oomph::AxisymmetricNavierStokesEquations::viscosity_ratio(), w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and oomph::TimeStepper::weight().

◆ further_build()

void oomph::RefineableAxisymmetricNavierStokesEquations::further_build ( )
inlinevirtual

Further build: pass the pointers down to the sons.

Reimplemented from oomph::RefineableElement.

Reimplemented in oomph::RefineableAxisymmetricQCrouzeixRaviartElement.

137  {
138  // Find the father element
139  RefineableAxisymmetricNavierStokesEquations* cast_father_element_pt =
140  dynamic_cast<RefineableAxisymmetricNavierStokesEquations*>(
141  this->father_element_pt());
142 
143  // Set the viscosity ratio pointer
144  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
145  // Set the density ratio pointer
146  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
147  // Set pointer to global Reynolds number
148  this->Re_pt = cast_father_element_pt->re_pt();
149  // Set pointer to global Reynolds number x Strouhal number (=Womersley)
150  this->ReSt_pt = cast_father_element_pt->re_st_pt();
151  // Set pointer to global Reynolds number x inverse Froude number
152  this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
153  // Set pointer to the global Reynolds number x inverse Rossby number
154  this->ReInvRo_pt = cast_father_element_pt->re_invro_pt();
155  // Set pointer to global gravity Vector
156  this->G_pt = cast_father_element_pt->g_pt();
157 
158  // Set pointer to body force function
159  this->Body_force_fct_pt =
160  cast_father_element_pt->axi_nst_body_force_fct_pt();
161 
162  // Set pointer to volumetric source function
163  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
164 
165  // Set the ALE_is_disabled flag
166  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
167  }
oomph::AxisymmetricNavierStokesEquations::Source_fct_pt
double(* Source_fct_pt)(const double &time, const Vector< double > &x)
Pointer to volumetric source function.
Definition: axisym_navier_stokes_elements.h:168
oomph::AxisymmetricNavierStokesEquations::G_pt
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: axisym_navier_stokes_elements.h:160
oomph::AxisymmetricNavierStokesEquations::ReInvFr_pt
double * ReInvFr_pt
Definition: axisym_navier_stokes_elements.h:153
oomph::AxisymmetricNavierStokesEquations::Re_pt
double * Re_pt
Pointer to global Reynolds number.
Definition: axisym_navier_stokes_elements.h:146
oomph::AxisymmetricNavierStokesEquations::ReSt_pt
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: axisym_navier_stokes_elements.h:149
oomph::AxisymmetricNavierStokesEquations::Density_Ratio_pt
double * Density_Ratio_pt
Definition: axisym_navier_stokes_elements.h:141
oomph::AxisymmetricNavierStokesEquations::Body_force_fct_pt
void(* Body_force_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to body force function.
Definition: axisym_navier_stokes_elements.h:163
oomph::AxisymmetricNavierStokesEquations::ReInvRo_pt
double * ReInvRo_pt
Definition: axisym_navier_stokes_elements.h:157
oomph::AxisymmetricNavierStokesEquations::Viscosity_Ratio_pt
double * Viscosity_Ratio_pt
Definition: axisym_navier_stokes_elements.h:137
oomph::RefineableAxisymmetricNavierStokesEquations::RefineableAxisymmetricNavierStokesEquations
RefineableAxisymmetricNavierStokesEquations()
Empty Constructor.
Definition: refineable_axisym_navier_stokes_elements.h:67
oomph::RefineableElement::father_element_pt
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Definition: refineable_elements.h:539

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

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

◆ geometric_jacobian()

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

Fill in the geometric Jacobian, which in this case is r.

Reimplemented from oomph::ElementWithZ2ErrorEstimator.

130  {
131  return x[0];
132  }
plotDoE.x
list x
Definition: plotDoE.py:28

References plotDoE::x.

◆ get_dresidual_dnodal_coordinates()

void oomph::RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates ( RankThreeTensor< double > &  dresidual_dnodal_coordinates)
virtual

Compute derivatives of elemental residual vector with respect to nodal coordinates. This function computes these terms analytically and overwrites the default implementation in the FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}

Compute derivatives of elemental residual vector with respect to nodal coordinates. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij} Overloads the FD-based version in the FiniteElement base class.

Reimplemented from oomph::AxisymmetricNavierStokesEquations.

984  {
985  // Create an Oomph Lib warning
986  std::string warning_message = "This function has not been tested.\n";
987  std::string function = "RefineableAxisymmetricNavierStokesEquations::\n";
988  function += "get_dresidual_dnodal_coordinates(...)";
989  OomphLibWarning(warning_message, function, OOMPH_EXCEPTION_LOCATION);
990 
991  // Return immediately if there are no dofs
992  if (ndof() == 0)
993  {
994  return;
995  }
996 
997  // Determine number of nodes in element
998  const unsigned n_node = nnode();
999 
1000  // Get continuous time from timestepper of first node
1001  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1002 
1003  // Determine number of pressure dofs in element
1004  const unsigned n_pres = this->npres_axi_nst();
1005 
1006  // Find the indices at which the local velocities are stored
1007  unsigned u_nodal_index[3];
1008  for (unsigned i = 0; i < 3; i++)
1009  {
1010  u_nodal_index[i] = this->u_index_axi_nst(i);
1011  }
1012 
1013  // Which nodal value represents the pressure? (Negative if pressure
1014  // is not based on nodal interpolation).
1015  const int p_index = this->p_nodal_index_axi_nst();
1016 
1017  // Local array of booleans that are true if the l-th pressure value is
1018  // hanging (avoid repeated virtual function calls)
1019  bool pressure_dof_is_hanging[n_pres];
1020 
1021  // If the pressure is stored at a node
1022  if (p_index >= 0)
1023  {
1024  // Read out whether the pressure is hanging
1025  for (unsigned l = 0; l < n_pres; ++l)
1026  {
1027  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
1028  }
1029  }
1030  // Otherwise the pressure is not stored at a node and so cannot hang
1031  else
1032  {
1033  for (unsigned l = 0; l < n_pres; ++l)
1034  {
1035  pressure_dof_is_hanging[l] = false;
1036  }
1037  }
1038 
1039  // Set up memory for the shape and test functions
1040  // Note that there are only two dimensions, r and z, in this problem
1041  Shape psif(n_node), testf(n_node);
1042  DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
1043 
1044  // Set up memory for pressure shape and test functions
1045  Shape psip(n_pres), testp(n_pres);
1046 
1047  // Determine number of shape controlling nodes
1048  const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1049 
1050  // Deriatives of shape fct derivatives w.r.t. nodal coords
1051  RankFourTensor<double> d_dpsifdx_dX(2, n_shape_controlling_node, n_node, 2);
1052  RankFourTensor<double> d_dtestfdx_dX(
1053  2, n_shape_controlling_node, n_node, 2);
1054 
1055  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1056  DenseMatrix<double> dJ_dX(2, n_shape_controlling_node);
1057 
1058  // Derivatives of derivative of u w.r.t. nodal coords
1059  // Note that the entry d_dudx_dX(p,q,i,k) contains the derivative w.r.t.
1060  // nodal coordinate X_pq of du_i/dx_k. Since there are three velocity
1061  // components, i can take on the values 0, 1 and 2, while k and p only
1062  // take on the values 0 and 1 since there are only two spatial dimensions.
1063  RankFourTensor<double> d_dudx_dX(2, n_shape_controlling_node, 3, 2);
1064 
1065  // Derivatives of nodal velocities w.r.t. nodal coords:
1066  // Assumption: Interaction only local via no-slip so
1067  // X_pq only affects U_iq.
1068  // Note that the entry d_U_dX(p,q,i) contains the derivative w.r.t. nodal
1069  // coordinate X_pq of nodal value U_iq. In other words this entry is the
1070  // derivative of the i-th velocity component w.r.t. the p-th spatial
1071  // coordinate, taken at the q-th local node.
1072  RankThreeTensor<double> d_U_dX(2, n_shape_controlling_node, 3, 0.0);
1073 
1074  // Determine the number of integration points
1075  const unsigned n_intpt = integral_pt()->nweight();
1076 
1077  // Vector to hold local coordinates (two dimensions)
1078  Vector<double> s(2);
1079 
1080  // Get physical variables from element
1081  // (Reynolds number must be multiplied by the density ratio)
1082  const double scaled_re = this->re() * this->density_ratio();
1083  const double scaled_re_st = this->re_st() * this->density_ratio();
1084  const double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
1085  const double scaled_re_inv_ro = this->re_invro() * this->density_ratio();
1086  const double visc_ratio = this->viscosity_ratio();
1087  const Vector<double> G = this->g();
1088 
1089  // FD step
1090  double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
1091 
1092  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1093  // Assumption: Interaction only local via no-slip so
1094  // X_ij only affects U_ij.
1095  bool element_has_node_with_aux_node_update_fct = false;
1096 
1097  std::map<Node*, unsigned> local_shape_controlling_node_lookup =
1098  shape_controlling_node_lookup();
1099 
1100  // FD loop over shape-controlling nodes
1101  for (std::map<Node*, unsigned>::iterator it =
1102  local_shape_controlling_node_lookup.begin();
1103  it != local_shape_controlling_node_lookup.end();
1104  it++)
1105  {
1106  // Get pointer to q-th local shape-controlling node
1107  Node* nod_pt = it->first;
1108 
1109  // Get its number
1110  unsigned q = it->second;
1111 
1112  // Only compute if there's a node-update fct involved
1113  if (nod_pt->has_auxiliary_node_update_fct_pt())
1114  {
1115  element_has_node_with_aux_node_update_fct = true;
1116 
1117  // This functionality has not been tested so issue a warning
1118  // to this effect
1119  std::ostringstream warning_stream;
1120  warning_stream << "\nThe functionality to evaluate the additional"
1121  << "\ncontribution to the deriv of the residual eqn"
1122  << "\nw.r.t. the nodal coordinates which comes about"
1123  << "\nif a node's values are updated using an auxiliary"
1124  << "\nnode update function has NOT been tested for"
1125  << "\nrefineable axisymmetric Navier-Stokes elements."
1126  << "\nUse at your own risk" << std::endl;
1127  OomphLibWarning(warning_stream.str(),
1128  "RefineableAxisymmetricNavierStokesEquations::get_"
1129  "dresidual_dnodal_coordinates",
1130  OOMPH_EXCEPTION_LOCATION);
1131 
1132  // Current nodal velocity
1133  Vector<double> u_ref(3);
1134  for (unsigned i = 0; i < 3; i++)
1135  {
1136  u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
1137  }
1138 
1139  // FD
1140  for (unsigned p = 0; p < 2; p++)
1141  {
1142  // Make backup
1143  double backup = nod_pt->x(p);
1144 
1145  // Do FD step. No node update required as we're
1146  // attacking the coordinate directly...
1147  nod_pt->x(p) += eps_fd;
1148 
1149  // Do auxiliary node update (to apply no slip)
1150  nod_pt->perform_auxiliary_node_update_fct();
1151 
1152  // Loop over velocity components
1153  for (unsigned i = 0; i < 3; i++)
1154  {
1155  // Evaluate
1156  d_U_dX(p, q, i) =
1157  (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
1158  }
1159 
1160  // Reset
1161  nod_pt->x(p) = backup;
1162 
1163  // Do auxiliary node update (to apply no slip)
1164  nod_pt->perform_auxiliary_node_update_fct();
1165  }
1166  }
1167  }
1168 
1169  // Integer to store the local equation number
1170  int local_eqn = 0;
1171 
1172  // Pointers to hang info object
1173  HangInfo* hang_info_pt = 0;
1174 
1175  // Loop over the integration points
1176  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1177  {
1178  // Assign values of s
1179  for (unsigned i = 0; i < 2; i++)
1180  {
1181  s[i] = integral_pt()->knot(ipt, i);
1182  }
1183 
1184  // Get the integral weight
1185  const double w = integral_pt()->weight(ipt);
1186 
1187  // Call the derivatives of the shape and test functions
1188  const double J =
1189  this->dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,
1190  psif,
1191  dpsifdx,
1192  d_dpsifdx_dX,
1193  testf,
1194  dtestfdx,
1195  d_dtestfdx_dX,
1196  dJ_dX);
1197 
1198  // Call the pressure shape and test functions
1199  this->pshape_axi_nst(s, psip, testp);
1200 
1201  // Allocate storage for the position and the derivative of the
1202  // mesh positions w.r.t. time
1203  Vector<double> interpolated_x(2, 0.0);
1204  Vector<double> mesh_velocity(2, 0.0);
1205 
1206  // Allocate storage for the pressure, velocity components and their
1207  // time and spatial derivatives
1208  double interpolated_p = 0.0;
1209  Vector<double> interpolated_u(3, 0.0);
1210  Vector<double> dudt(3, 0.0);
1211  DenseMatrix<double> interpolated_dudx(3, 2, 0.0);
1212 
1213  // Calculate pressure at integration point
1214  for (unsigned l = 0; l < n_pres; l++)
1215  {
1216  interpolated_p += this->p_axi_nst(l) * psip[l];
1217  }
1218 
1219  // Calculate velocities and derivatives at integration point
1220  // ---------------------------------------------------------
1221 
1222  // Loop over nodes
1223  for (unsigned l = 0; l < n_node; l++)
1224  {
1225  // Cache the shape function
1226  const double psif_ = psif(l);
1227 
1228  // Loop over the two coordinate directions
1229  for (unsigned i = 0; i < 2; i++)
1230  {
1231  interpolated_x[i] += nodal_position(l, i) * psif_;
1232  }
1233 
1234  // Loop over the three velocity directions
1235  for (unsigned i = 0; i < 3; i++)
1236  {
1237  // Get the nodal value
1238  const double u_value = nodal_value(l, u_nodal_index[i]);
1239  interpolated_u[i] += u_value * psif_;
1240  dudt[i] += this->du_dt_axi_nst(l, i) * psif_;
1241 
1242  // Loop over derivative directions
1243  for (unsigned j = 0; j < 2; j++)
1244  {
1245  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1246  }
1247  }
1248  }
1249 
1250  // Get the mesh velocity if ALE is enabled
1251  if (!this->ALE_is_disabled)
1252  {
1253  // Loop over nodes
1254  for (unsigned l = 0; l < n_node; l++)
1255  {
1256  // Loop over the two coordinate directions
1257  for (unsigned i = 0; i < 2; i++)
1258  {
1259  mesh_velocity[i] += this->dnodal_position_dt(l, i) * psif[l];
1260  }
1261  }
1262  }
1263 
1264  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1265 
1266  // Loop over shape-controlling nodes
1267  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1268  {
1269  // Loop over the two coordinate directions
1270  for (unsigned p = 0; p < 2; p++)
1271  {
1272  // Loop over the three velocity components
1273  for (unsigned i = 0; i < 3; i++)
1274  {
1275  // Loop over the two coordinate directions
1276  for (unsigned k = 0; k < 2; k++)
1277  {
1278  double aux = 0.0;
1279 
1280  // Loop over nodes and add contribution
1281  for (unsigned j = 0; j < n_node; j++)
1282  {
1283  aux +=
1284  nodal_value(j, u_nodal_index[i]) * d_dpsifdx_dX(p, q, j, k);
1285  }
1286  d_dudx_dX(p, q, i, k) = aux;
1287  }
1288  }
1289  }
1290  }
1291 
1292  // Get weight of actual nodal position/value in computation of mesh
1293  // velocity from positional/value time stepper
1294  const double pos_time_weight =
1295  node_pt(0)->position_time_stepper_pt()->weight(1, 0);
1296  const double val_time_weight =
1297  node_pt(0)->time_stepper_pt()->weight(1, 0);
1298 
1299  // Get the user-defined body force terms
1300  Vector<double> body_force(3);
1301  this->get_body_force_axi_nst(time, ipt, s, interpolated_x, body_force);
1302 
1303  // Get the user-defined source function
1304  const double source = this->get_source_fct(time, ipt, interpolated_x);
1305 
1306  // Get gradient of body force function
1307  DenseMatrix<double> d_body_force_dx(3, 2, 0.0);
1308  this->get_body_force_gradient_axi_nst(
1309  time, ipt, s, interpolated_x, d_body_force_dx);
1310 
1311  // Get gradient of source function
1312  Vector<double> source_gradient(2, 0.0);
1313  this->get_source_fct_gradient(time, ipt, interpolated_x, source_gradient);
1314 
1315  // Cache r (the first position component)
1316  const double r = interpolated_x[0];
1317 
1318  // Assemble shape derivatives
1319  // --------------------------
1320 
1321  // ==================
1322  // MOMENTUM EQUATIONS
1323  // ==================
1324 
1325  // Number of master nodes and storage for the weight of the shape function
1326  unsigned n_master = 1;
1327  double hang_weight = 1.0;
1328 
1329  // Loop over the test functions
1330  for (unsigned l = 0; l < n_node; l++)
1331  {
1332  // Cache the test function
1333  const double testf_ = testf[l];
1334 
1335  // Local boolean to indicate whether the node is hanging
1336  bool is_node_hanging = node_pt(l)->is_hanging();
1337 
1338  // If the node is hanging
1339  if (is_node_hanging)
1340  {
1341  hang_info_pt = node_pt(l)->hanging_pt();
1342 
1343  // Read out number of master nodes from hanging data
1344  n_master = hang_info_pt->nmaster();
1345  }
1346  // Otherwise the node is its own master
1347  else
1348  {
1349  n_master = 1;
1350  }
1351 
1352  // Loop over the master nodes
1353  for (unsigned m = 0; m < n_master; m++)
1354  {
1355  // --------------------------------
1356  // FIRST (RADIAL) MOMENTUM EQUATION
1357  // --------------------------------
1358 
1359  // Get the equation number
1360  // If the node is hanging
1361  if (is_node_hanging)
1362  {
1363  // Get the equation number from the master node
1364  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1365  u_nodal_index[0]);
1366  // Get the hang weight from the master node
1367  hang_weight = hang_info_pt->master_weight(m);
1368  }
1369  // If the node is not hanging
1370  else
1371  {
1372  // Local equation number
1373  local_eqn = this->nodal_local_eqn(l, u_nodal_index[0]);
1374 
1375  // Node contributes with full weight
1376  hang_weight = 1.0;
1377  }
1378 
1379  // If it's not a boundary condition
1380  if (local_eqn >= 0)
1381  {
1382  // Loop over the two coordinate directions
1383  for (unsigned p = 0; p < 2; p++)
1384  {
1385  // Loop over shape controlling nodes
1386  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1387  {
1388  // Residual x deriv of Jacobian
1389  // ----------------------------
1390 
1391  // Add the user-defined body force terms
1392  double sum = r * body_force[0] * testf_;
1393 
1394  // Add the gravitational body force term
1395  sum += r * scaled_re_inv_fr * testf_ * G[0];
1396 
1397  // Add the pressure gradient term
1398  sum += interpolated_p * (testf_ + r * dtestfdx(l, 0));
1399 
1400  // Add the stress tensor terms
1401  // The viscosity ratio needs to go in here to ensure
1402  // continuity of normal stress is satisfied even in flows
1403  // with zero pressure gradient!
1404  sum -= visc_ratio * r * (1.0 + Gamma[0]) *
1405  interpolated_dudx(0, 0) * dtestfdx(l, 0);
1406 
1407  sum -= visc_ratio * r *
1408  (interpolated_dudx(0, 1) +
1409  Gamma[0] * interpolated_dudx(1, 0)) *
1410  dtestfdx(l, 1);
1411 
1412  sum -= visc_ratio * (1.0 + Gamma[0]) * interpolated_u[0] *
1413  testf_ / r;
1414 
1415  // Add the du/dt term
1416  sum -= scaled_re_st * r * dudt[0] * testf_;
1417 
1418  // Add the convective terms
1419  sum -= scaled_re *
1420  (r * interpolated_u[0] * interpolated_dudx(0, 0) -
1421  interpolated_u[2] * interpolated_u[2] +
1422  r * interpolated_u[1] * interpolated_dudx(0, 1)) *
1423  testf_;
1424 
1425  // Add the mesh velocity terms
1426  if (!ALE_is_disabled)
1427  {
1428  for (unsigned k = 0; k < 2; k++)
1429  {
1430  sum += scaled_re_st * r * mesh_velocity[k] *
1431  interpolated_dudx(0, k) * testf_;
1432  }
1433  }
1434 
1435  // Add the Coriolis term
1436  sum += 2.0 * r * scaled_re_inv_ro * interpolated_u[2] * testf_;
1437 
1438  // Multiply through by deriv of Jacobian and integration weight
1439  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1440  sum * dJ_dX(p, q) * w * hang_weight;
1441 
1442  // Derivative of residual x Jacobian
1443  // ---------------------------------
1444 
1445  // Body force
1446  sum = r * d_body_force_dx(0, p) * psif[q] * testf_;
1447  if (p == 0)
1448  {
1449  sum += body_force[0] * psif[q] * testf_;
1450  }
1451 
1452  // Gravitational body force
1453  if (p == 0)
1454  {
1455  sum += scaled_re_inv_fr * G[0] * psif[q] * testf_;
1456  }
1457 
1458  // Pressure gradient term
1459  sum += r * interpolated_p * d_dtestfdx_dX(p, q, l, 0);
1460  if (p == 0)
1461  {
1462  sum += interpolated_p * psif[q] * dtestfdx(l, 0);
1463  }
1464 
1465  // Viscous terms
1466  sum -=
1467  r * visc_ratio *
1468  ((1.0 + Gamma[0]) *
1469  (d_dudx_dX(p, q, 0, 0) * dtestfdx(l, 0) +
1470  interpolated_dudx(0, 0) * d_dtestfdx_dX(p, q, l, 0)) +
1471  (d_dudx_dX(p, q, 0, 1) + Gamma[0] * d_dudx_dX(p, q, 1, 0)) *
1472  dtestfdx(l, 1) +
1473  (interpolated_dudx(0, 1) +
1474  Gamma[0] * interpolated_dudx(1, 0)) *
1475  d_dtestfdx_dX(p, q, l, 1));
1476  if (p == 0)
1477  {
1478  sum -=
1479  visc_ratio *
1480  ((1.0 + Gamma[0]) *
1481  (interpolated_dudx(0, 0) * psif[q] * dtestfdx(l, 0) -
1482  interpolated_u[0] * psif[q] * testf_ / (r * r)) +
1483  (interpolated_dudx(0, 1) +
1484  Gamma[0] * interpolated_dudx(1, 0)) *
1485  psif[q] * dtestfdx(l, 1));
1486  }
1487 
1488  // Convective terms, including mesh velocity
1489  for (unsigned k = 0; k < 2; k++)
1490  {
1491  double tmp = scaled_re * interpolated_u[k];
1492  if (!ALE_is_disabled)
1493  {
1494  tmp -= scaled_re_st * mesh_velocity[k];
1495  }
1496  sum -= r * tmp * d_dudx_dX(p, q, 0, k) * testf_;
1497  if (p == 0)
1498  {
1499  sum -= tmp * interpolated_dudx(0, k) * psif[q] * testf_;
1500  }
1501  }
1502  if (!ALE_is_disabled)
1503  {
1504  sum += scaled_re_st * r * pos_time_weight *
1505  interpolated_dudx(0, p) * psif[q] * testf_;
1506  }
1507 
1508  // du/dt term
1509  if (p == 0)
1510  {
1511  sum -= scaled_re_st * dudt[0] * psif[q] * testf_;
1512  }
1513 
1514  // Coriolis term
1515  if (p == 0)
1516  {
1517  sum += 2.0 * scaled_re_inv_ro * interpolated_u[2] * psif[q] *
1518  testf_;
1519  }
1520 
1521  // Multiply through by Jacobian and integration weight
1522  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1523  sum * J * w * hang_weight;
1524 
1525  } // End of loop over shape controlling nodes q
1526  } // End of loop over coordinate directions p
1527 
1528  // Derivs w.r.t. to nodal velocities
1529  // ---------------------------------
1530  if (element_has_node_with_aux_node_update_fct)
1531  {
1532  // Loop over local nodes
1533  for (unsigned q_local = 0; q_local < n_node; q_local++)
1534  {
1535  // Number of master nodes and storage for the weight of
1536  // the shape function
1537  unsigned n_master2 = 1;
1538  double hang_weight2 = 1.0;
1539  HangInfo* hang_info2_pt = 0;
1540 
1541  // Local boolean to indicate whether the node is hanging
1542  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1543 
1544  Node* actual_shape_controlling_node_pt = node_pt(q_local);
1545 
1546  // If the node is hanging
1547  if (is_node_hanging2)
1548  {
1549  hang_info2_pt = node_pt(q_local)->hanging_pt();
1550 
1551  // Read out number of master nodes from hanging data
1552  n_master2 = hang_info2_pt->nmaster();
1553  }
1554  // Otherwise the node is its own master
1555  else
1556  {
1557  n_master2 = 1;
1558  }
1559 
1560  // Loop over the master nodes
1561  for (unsigned mm = 0; mm < n_master2; mm++)
1562  {
1563  if (is_node_hanging2)
1564  {
1565  actual_shape_controlling_node_pt =
1566  hang_info2_pt->master_node_pt(mm);
1567  hang_weight2 = hang_info2_pt->master_weight(mm);
1568  }
1569 
1570  // Find the corresponding number
1571  unsigned q = local_shape_controlling_node_lookup
1572  [actual_shape_controlling_node_pt];
1573 
1574  // Loop over the two coordinate directions
1575  for (unsigned p = 0; p < 2; p++)
1576  {
1577  // Contribution from deriv of first velocity component
1578  double tmp = scaled_re_st * r * val_time_weight *
1579  psif[q_local] * testf_;
1580  tmp += r * scaled_re * interpolated_dudx(0, 0) *
1581  psif[q_local] * testf_;
1582 
1583  for (unsigned k = 0; k < 2; k++)
1584  {
1585  double tmp2 = scaled_re * interpolated_u[k];
1586  if (!ALE_is_disabled)
1587  {
1588  tmp2 -= scaled_re_st * mesh_velocity[k];
1589  }
1590  tmp += r * tmp2 * dpsifdx(q_local, k) * testf_;
1591  }
1592 
1593  tmp += r * visc_ratio * (1.0 + Gamma[0]) *
1594  dpsifdx(q_local, 0) * dtestfdx(l, 0);
1595  tmp +=
1596  r * visc_ratio * dpsifdx(q_local, 1) * dtestfdx(l, 1);
1597  tmp += visc_ratio * (1.0 + Gamma[0]) * psif[q_local] *
1598  testf_ / r;
1599 
1600  // Multiply through by dU_0q/dX_pq
1601  double sum = -d_U_dX(p, q_local, 0) * tmp;
1602 
1603  // Contribution from deriv of second velocity component
1604  tmp = scaled_re * r * interpolated_dudx(0, 1) *
1605  psif[q_local] * testf_;
1606  tmp += r * visc_ratio * Gamma[0] * dpsifdx(q_local, 0) *
1607  dtestfdx(l, 1);
1608 
1609  // Multiply through by dU_1q/dX_pq
1610  sum -= d_U_dX(p, q, 1) * tmp;
1611 
1612  // Contribution from deriv of third velocity component
1613  tmp =
1614  2.0 *
1615  (r * scaled_re_inv_ro + scaled_re * interpolated_u[2]) *
1616  psif[q_local] * testf_;
1617 
1618  // Multiply through by dU_2q/dX_pq
1619  sum += d_U_dX(p, q, 2) * tmp;
1620 
1621  // Multiply through by Jacobian and integration weight
1622  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1623  sum * J * w * hang_weight * hang_weight2;
1624  }
1625  } // End of loop over master nodes
1626  } // End of loop over local nodes
1627  } // End of if(element_has_node_with_aux_node_update_fct)
1628  } // End of if it's not a boundary condition
1629 
1630  // --------------------------------
1631  // SECOND (AXIAL) MOMENTUM EQUATION
1632  // --------------------------------
1633 
1634  // Get the equation number
1635  // If the node is hanging
1636  if (is_node_hanging)
1637  {
1638  // Get the equation number from the master node
1639  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1640  u_nodal_index[1]);
1641  // Get the hang weight from the master node
1642  hang_weight = hang_info_pt->master_weight(m);
1643  }
1644  // If the node is not hanging
1645  else
1646  {
1647  // Local equation number
1648  local_eqn = this->nodal_local_eqn(l, u_nodal_index[1]);
1649 
1650  // Node contributes with full weight
1651  hang_weight = 1.0;
1652  }
1653 
1654  // If it's not a boundary condition
1655  if (local_eqn >= 0)
1656  {
1657  // Loop over the two coordinate directions
1658  for (unsigned p = 0; p < 2; p++)
1659  {
1660  // Loop over shape controlling nodes
1661  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1662  {
1663  // Residual x deriv of Jacobian
1664  // ----------------------------
1665 
1666  // Add the user-defined body force terms
1667  double sum = r * body_force[1] * testf_;
1668 
1669  // Add the gravitational body force term
1670  sum += r * scaled_re_inv_fr * testf_ * G[1];
1671 
1672  // Add the pressure gradient term
1673  sum += r * interpolated_p * dtestfdx(l, 1);
1674 
1675  // Add the stress tensor terms
1676  // The viscosity ratio needs to go in here to ensure
1677  // continuity of normal stress is satisfied even in flows
1678  // with zero pressure gradient!
1679  sum -= visc_ratio * r *
1680  (interpolated_dudx(1, 0) +
1681  Gamma[1] * interpolated_dudx(0, 1)) *
1682  dtestfdx(l, 0);
1683 
1684  sum -= visc_ratio * r * (1.0 + Gamma[1]) *
1685  interpolated_dudx(1, 1) * dtestfdx(l, 1);
1686 
1687  // Add the du/dt term
1688  sum -= scaled_re_st * r * dudt[1] * testf_;
1689 
1690  // Add the convective terms
1691  sum -= scaled_re *
1692  (r * interpolated_u[0] * interpolated_dudx(1, 0) +
1693  r * interpolated_u[1] * interpolated_dudx(1, 1)) *
1694  testf_;
1695 
1696  // Add the mesh velocity terms
1697  if (!ALE_is_disabled)
1698  {
1699  for (unsigned k = 0; k < 2; k++)
1700  {
1701  sum += scaled_re_st * r * mesh_velocity[k] *
1702  interpolated_dudx(1, k) * testf_;
1703  }
1704  }
1705 
1706  // Multiply through by deriv of Jacobian and integration weight
1707  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1708  sum * dJ_dX(p, q) * w * hang_weight;
1709 
1710  // Derivative of residual x Jacobian
1711  // ---------------------------------
1712 
1713  // Body force
1714  sum = r * d_body_force_dx(1, p) * psif[q] * testf_;
1715  if (p == 0)
1716  {
1717  sum += body_force[1] * psif[q] * testf_;
1718  }
1719 
1720  // Gravitational body force
1721  if (p == 0)
1722  {
1723  sum += scaled_re_inv_fr * G[1] * psif[q] * testf_;
1724  }
1725 
1726  // Pressure gradient term
1727  sum += r * interpolated_p * d_dtestfdx_dX(p, q, l, 1);
1728  if (p == 0)
1729  {
1730  sum += interpolated_p * psif[q] * dtestfdx(l, 1);
1731  }
1732 
1733  // Viscous terms
1734  sum -=
1735  r * visc_ratio *
1736  ((d_dudx_dX(p, q, 1, 0) + Gamma[1] * d_dudx_dX(p, q, 0, 1)) *
1737  dtestfdx(l, 0) +
1738  (interpolated_dudx(1, 0) +
1739  Gamma[1] * interpolated_dudx(0, 1)) *
1740  d_dtestfdx_dX(p, q, l, 0) +
1741  (1.0 + Gamma[1]) * d_dudx_dX(p, q, 1, 1) * dtestfdx(l, 1) +
1742  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
1743  d_dtestfdx_dX(p, q, l, 1));
1744  if (p == 0)
1745  {
1746  sum -=
1747  visc_ratio * ((interpolated_dudx(1, 0) +
1748  Gamma[1] * interpolated_dudx(0, 1)) *
1749  psif[q] * dtestfdx(l, 0) +
1750  (1.0 + Gamma[1]) * interpolated_dudx(1, 1) *
1751  psif[q] * dtestfdx(l, 1));
1752  }
1753 
1754  // Convective terms, including mesh velocity
1755  for (unsigned k = 0; k < 2; k++)
1756  {
1757  double tmp = scaled_re * interpolated_u[k];
1758  if (!ALE_is_disabled)
1759  {
1760  tmp -= scaled_re_st * mesh_velocity[k];
1761  }
1762  sum -= r * tmp * d_dudx_dX(p, q, 1, k) * testf_;
1763  if (p == 0)
1764  {
1765  sum -= tmp * interpolated_dudx(1, k) * psif[q] * testf_;
1766  }
1767  }
1768  if (!ALE_is_disabled)
1769  {
1770  sum += scaled_re_st * r * pos_time_weight *
1771  interpolated_dudx(1, p) * psif[q] * testf_;
1772  }
1773 
1774  // du/dt term
1775  if (p == 0)
1776  {
1777  sum -= scaled_re_st * dudt[1] * psif[q] * testf_;
1778  }
1779 
1780  // Multiply through by Jacobian and integration weight
1781  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1782  sum * J * w * hang_weight;
1783 
1784  } // End of loop over shape controlling nodes q
1785  } // End of loop over coordinate directions p
1786 
1787  // Derivs w.r.t. to nodal velocities
1788  // ---------------------------------
1789  if (element_has_node_with_aux_node_update_fct)
1790  {
1791  // Loop over local nodes
1792  for (unsigned q_local = 0; q_local < n_node; q_local++)
1793  {
1794  // Number of master nodes and storage for the weight of
1795  // the shape function
1796  unsigned n_master2 = 1;
1797  double hang_weight2 = 1.0;
1798  HangInfo* hang_info2_pt = 0;
1799 
1800  // Local boolean to indicate whether the node is hanging
1801  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1802 
1803  Node* actual_shape_controlling_node_pt = node_pt(q_local);
1804 
1805  // If the node is hanging
1806  if (is_node_hanging2)
1807  {
1808  hang_info2_pt = node_pt(q_local)->hanging_pt();
1809 
1810  // Read out number of master nodes from hanging data
1811  n_master2 = hang_info2_pt->nmaster();
1812  }
1813  // Otherwise the node is its own master
1814  else
1815  {
1816  n_master2 = 1;
1817  }
1818 
1819  // Loop over the master nodes
1820  for (unsigned mm = 0; mm < n_master2; mm++)
1821  {
1822  if (is_node_hanging2)
1823  {
1824  actual_shape_controlling_node_pt =
1825  hang_info2_pt->master_node_pt(mm);
1826  hang_weight2 = hang_info2_pt->master_weight(mm);
1827  }
1828 
1829  // Find the corresponding number
1830  unsigned q = local_shape_controlling_node_lookup
1831  [actual_shape_controlling_node_pt];
1832 
1833  // Loop over the two coordinate directions
1834  for (unsigned p = 0; p < 2; p++)
1835  {
1836  // Contribution from deriv of first velocity component
1837  double tmp = scaled_re * r * interpolated_dudx(1, 0) *
1838  psif[q_local] * testf_;
1839  tmp += r * visc_ratio * Gamma[1] * dpsifdx(q_local, 1) *
1840  dtestfdx(l, 0);
1841 
1842  // Multiply through by dU_0q/dX_pq
1843  double sum = -d_U_dX(p, q, 0) * tmp;
1844 
1845  // Contribution from deriv of second velocity component
1846  tmp = scaled_re_st * r * val_time_weight * psif[q_local] *
1847  testf_;
1848  tmp += r * scaled_re * interpolated_dudx(1, 1) *
1849  psif[q_local] * testf_;
1850 
1851  for (unsigned k = 0; k < 2; k++)
1852  {
1853  double tmp2 = scaled_re * interpolated_u[k];
1854  if (!ALE_is_disabled)
1855  {
1856  tmp2 -= scaled_re_st * mesh_velocity[k];
1857  }
1858  tmp += r * tmp2 * dpsifdx(q_local, k) * testf_;
1859  }
1860  tmp +=
1861  r * visc_ratio *
1862  (dpsifdx(q_local, 0) * dtestfdx(l, 0) +
1863  (1.0 + Gamma[1]) * dpsifdx(q_local, 1) * dtestfdx(l, 1));
1864 
1865  // Multiply through by dU_1q/dX_pq
1866  sum -= d_U_dX(p, q, 1) * tmp;
1867 
1868  // Multiply through by Jacobian and integration weight
1869  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1870  sum * J * w * hang_weight * hang_weight2;
1871  }
1872  } // End of loop over master nodes
1873  } // End of loop over local nodes
1874  } // End of if(element_has_node_with_aux_node_update_fct)
1875  } // End of if it's not a boundary condition
1876 
1877  // -----------------------------------
1878  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
1879  // -----------------------------------
1880 
1881  // Get the equation number
1882  // If the node is hanging
1883  if (is_node_hanging)
1884  {
1885  // Get the equation number from the master node
1886  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1887  u_nodal_index[2]);
1888  // Get the hang weight from the master node
1889  hang_weight = hang_info_pt->master_weight(m);
1890  }
1891  // If the node is not hanging
1892  else
1893  {
1894  // Local equation number
1895  local_eqn = this->nodal_local_eqn(l, u_nodal_index[2]);
1896 
1897  // Node contributes with full weight
1898  hang_weight = 1.0;
1899  }
1900 
1901  // If it's not a boundary condition
1902  if (local_eqn >= 0)
1903  {
1904  // Loop over the two coordinate directions
1905  for (unsigned p = 0; p < 2; p++)
1906  {
1907  // Loop over shape controlling nodes
1908  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1909  {
1910  // Residual x deriv of Jacobian
1911  // ----------------------------
1912 
1913  // Add the user-defined body force terms
1914  double sum = r * body_force[2] * testf_;
1915 
1916  // Add the gravitational body force term
1917  sum += r * scaled_re_inv_fr * testf_ * G[2];
1918 
1919  // There is NO pressure gradient term
1920 
1921  // Add in the stress tensor terms
1922  // The viscosity ratio needs to go in here to ensure
1923  // continuity of normal stress is satisfied even in flows
1924  // with zero pressure gradient!
1925  sum -=
1926  visc_ratio *
1927  (r * interpolated_dudx(2, 0) - Gamma[0] * interpolated_u[2]) *
1928  dtestfdx(l, 0);
1929 
1930  sum -=
1931  visc_ratio * r * interpolated_dudx(2, 1) * dtestfdx(l, 1);
1932 
1933  sum -= visc_ratio *
1934  ((interpolated_u[2] / r) -
1935  Gamma[0] * interpolated_dudx(2, 0)) *
1936  testf_;
1937 
1938  // Add the du/dt term
1939  sum -= scaled_re_st * r * dudt[2] * testf_;
1940 
1941  // Add the convective terms
1942  sum -= scaled_re *
1943  (r * interpolated_u[0] * interpolated_dudx(2, 0) +
1944  interpolated_u[0] * interpolated_u[2] +
1945  r * interpolated_u[1] * interpolated_dudx(2, 1)) *
1946  testf_;
1947 
1948  // Add the mesh velocity terms
1949  if (!ALE_is_disabled)
1950  {
1951  for (unsigned k = 0; k < 2; k++)
1952  {
1953  sum += scaled_re_st * r * mesh_velocity[k] *
1954  interpolated_dudx(2, k) * testf_;
1955  }
1956  }
1957 
1958  // Add the Coriolis term
1959  sum -= 2.0 * r * scaled_re_inv_ro * interpolated_u[0] * testf_;
1960 
1961  // Multiply through by deriv of Jacobian and integration weight
1962  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1963  sum * dJ_dX(p, q) * w * hang_weight;
1964 
1965  // Derivative of residual x Jacobian
1966  // ---------------------------------
1967 
1968  // Body force
1969  sum = r * d_body_force_dx(2, p) * psif[q] * testf_;
1970  if (p == 0)
1971  {
1972  sum += body_force[2] * psif[q] * testf_;
1973  }
1974 
1975  // Gravitational body force
1976  if (p == 0)
1977  {
1978  sum += scaled_re_inv_fr * G[2] * psif[q] * testf_;
1979  }
1980 
1981  // There is NO pressure gradient term
1982 
1983  // Viscous terms
1984  sum -= visc_ratio * r *
1985  (d_dudx_dX(p, q, 2, 0) * dtestfdx(l, 0) +
1986  interpolated_dudx(2, 0) * d_dtestfdx_dX(p, q, l, 0) +
1987  d_dudx_dX(p, q, 2, 1) * dtestfdx(l, 1) +
1988  interpolated_dudx(2, 1) * d_dtestfdx_dX(p, q, l, 1));
1989 
1990  sum += visc_ratio * Gamma[0] * d_dudx_dX(p, q, 2, 0) * testf_;
1991 
1992  if (p == 0)
1993  {
1994  sum -= visc_ratio *
1995  (interpolated_dudx(2, 0) * psif[q] * dtestfdx(l, 0) +
1996  interpolated_dudx(2, 1) * psif[q] * dtestfdx(l, 1) +
1997  interpolated_u[2] * psif[q] * testf_ / (r * r));
1998  }
1999 
2000  // Convective terms, including mesh velocity
2001  for (unsigned k = 0; k < 2; k++)
2002  {
2003  double tmp = scaled_re * interpolated_u[k];
2004  if (!ALE_is_disabled)
2005  {
2006  tmp -= scaled_re_st * mesh_velocity[k];
2007  }
2008  sum -= r * tmp * d_dudx_dX(p, q, 2, k) * testf_;
2009  if (p == 0)
2010  {
2011  sum -= tmp * interpolated_dudx(2, k) * psif[q] * testf_;
2012  }
2013  }
2014  if (!ALE_is_disabled)
2015  {
2016  sum += scaled_re_st * r * pos_time_weight *
2017  interpolated_dudx(2, p) * psif[q] * testf_;
2018  }
2019 
2020  // du/dt term
2021  if (p == 0)
2022  {
2023  sum -= scaled_re_st * dudt[2] * psif[q] * testf_;
2024  }
2025 
2026  // Coriolis term
2027  if (p == 0)
2028  {
2029  sum -= 2.0 * scaled_re_inv_ro * interpolated_u[0] * psif[q] *
2030  testf_;
2031  }
2032 
2033  // Multiply through by Jacobian and integration weight
2034  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2035  sum * J * w * hang_weight;
2036 
2037  } // End of loop over shape controlling nodes q
2038  } // End of loop over coordinate directions p
2039 
2040  // Derivs w.r.t. to nodal velocities
2041  // ---------------------------------
2042  if (element_has_node_with_aux_node_update_fct)
2043  {
2044  // Loop over local nodes
2045  for (unsigned q_local = 0; q_local < n_node; q_local++)
2046  {
2047  // Number of master nodes and storage for the weight of
2048  // the shape function
2049  unsigned n_master2 = 1;
2050  double hang_weight2 = 1.0;
2051  HangInfo* hang_info2_pt = 0;
2052 
2053  // Local boolean to indicate whether the node is hanging
2054  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2055 
2056  Node* actual_shape_controlling_node_pt = node_pt(q_local);
2057 
2058  // If the node is hanging
2059  if (is_node_hanging2)
2060  {
2061  hang_info2_pt = node_pt(q_local)->hanging_pt();
2062 
2063  // Read out number of master nodes from hanging data
2064  n_master2 = hang_info2_pt->nmaster();
2065  }
2066  // Otherwise the node is its own master
2067  else
2068  {
2069  n_master2 = 1;
2070  }
2071 
2072  // Loop over the master nodes
2073  for (unsigned mm = 0; mm < n_master2; mm++)
2074  {
2075  if (is_node_hanging2)
2076  {
2077  actual_shape_controlling_node_pt =
2078  hang_info2_pt->master_node_pt(mm);
2079  hang_weight2 = hang_info2_pt->master_weight(mm);
2080  }
2081 
2082  // Find the corresponding number
2083  unsigned q = local_shape_controlling_node_lookup
2084  [actual_shape_controlling_node_pt];
2085 
2086  // Loop over the two coordinate directions
2087  for (unsigned p = 0; p < 2; p++)
2088  {
2089  // Contribution from deriv of first velocity component
2090  double tmp = (2.0 * r * scaled_re_inv_ro +
2091  r * scaled_re * interpolated_dudx(2, 0) -
2092  scaled_re * interpolated_u[2]) *
2093  psif[q_local] * testf_;
2094 
2095  // Multiply through by dU_0q/dX_pq
2096  double sum = -d_U_dX(p, q, 0) * tmp;
2097 
2098  // Contribution from deriv of second velocity component
2099  // Multiply through by dU_1q/dX_pq
2100  sum -= d_U_dX(p, q, 1) * scaled_re * r *
2101  interpolated_dudx(2, 1) * psif[q_local] * testf_;
2102 
2103  // Contribution from deriv of third velocity component
2104  tmp = scaled_re_st * r * val_time_weight * psif[q_local] *
2105  testf_;
2106  tmp -=
2107  scaled_re * interpolated_u[0] * psif[q_local] * testf_;
2108 
2109  for (unsigned k = 0; k < 2; k++)
2110  {
2111  double tmp2 = scaled_re * interpolated_u[k];
2112  if (!ALE_is_disabled)
2113  {
2114  tmp2 -= scaled_re_st * mesh_velocity[k];
2115  }
2116  tmp += r * tmp2 * dpsifdx(q_local, k) * testf_;
2117  }
2118  tmp +=
2119  visc_ratio *
2120  (r * dpsifdx(q_local, 0) - Gamma[0] * psif[q_local]) *
2121  dtestfdx(l, 0);
2122  tmp +=
2123  r * visc_ratio * dpsifdx(q_local, 1) * dtestfdx(l, 1);
2124  tmp +=
2125  visc_ratio *
2126  ((psif[q_local] / r) - Gamma[0] * dpsifdx(q_local, 0)) *
2127  testf_;
2128 
2129  // Multiply through by dU_2q/dX_pq
2130  sum -= d_U_dX(p, q, 2) * tmp;
2131 
2132  // Multiply through by Jacobian and integration weight
2133  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2134  sum * J * w * hang_weight * hang_weight2;
2135  }
2136  } // End of loop over master nodes
2137  } // End of loop over local nodes
2138  } // End of if(element_has_node_with_aux_node_update_fct)
2139  } // End of if it's not a boundary condition
2140  } // End of loop over master nodes
2141  } // End of loop over test functions
2142 
2143 
2144  // ===================
2145  // CONTINUITY EQUATION
2146  // ===================
2147 
2148  // Loop over the shape functions
2149  for (unsigned l = 0; l < n_pres; l++)
2150  {
2151  // If the pressure dof is hanging
2152  if (pressure_dof_is_hanging[l])
2153  {
2154  // Pressure dof is hanging so it must be nodal-based
2155  // Get the hang info object
2156  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
2157 
2158  // Get the number of master nodes from the pressure node
2159  n_master = hang_info_pt->nmaster();
2160  }
2161  // Otherwise the node is its own master
2162  else
2163  {
2164  n_master = 1;
2165  }
2166 
2167  // Loop over the master nodes
2168  for (unsigned m = 0; m < n_master; m++)
2169  {
2170  // Get the number of the unknown
2171  // If the pressure dof is hanging
2172  if (pressure_dof_is_hanging[l])
2173  {
2174  // Get the local equation from the master node
2175  local_eqn =
2176  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
2177  // Get the weight for the node
2178  hang_weight = hang_info_pt->master_weight(m);
2179  }
2180  else
2181  {
2182  local_eqn = this->p_local_eqn(l);
2183  hang_weight = 1.0;
2184  }
2185 
2186  // Cache the test function
2187  const double testp_ = testp[l];
2188 
2189  // If not a boundary conditions
2190  if (local_eqn >= 0)
2191  {
2192  // Loop over the two coordinate directions
2193  for (unsigned p = 0; p < 2; p++)
2194  {
2195  // Loop over shape controlling nodes
2196  for (unsigned q = 0; q < n_shape_controlling_node; q++)
2197  {
2198  // Residual x deriv of Jacobian
2199  //-----------------------------
2200 
2201  // Source term
2202  double aux = -r * source;
2203 
2204  // Gradient terms
2205  aux += (interpolated_u[0] + r * interpolated_dudx(0, 0) +
2206  r * interpolated_dudx(1, 1));
2207 
2208  // Multiply through by deriv of Jacobian and integration weight
2209  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2210  aux * dJ_dX(p, q) * testp_ * w * hang_weight;
2211 
2212  // Derivative of residual x Jacobian
2213  // ---------------------------------
2214 
2215  // Gradient of source function
2216  aux = -r * source_gradient[p] * psif[q];
2217  if (p == 0)
2218  {
2219  aux -= source * psif[q];
2220  }
2221 
2222  // Gradient terms
2223  aux += r * (d_dudx_dX(p, q, 0, 0) + d_dudx_dX(p, q, 1, 1));
2224  if (p == 0)
2225  {
2226  aux += (interpolated_dudx(0, 0) + interpolated_dudx(1, 1)) *
2227  psif[q];
2228  }
2229 
2230  // Multiply through by Jacobian and integration weight
2231  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2232  aux * testp_ * J * w * hang_weight;
2233  }
2234  }
2235 
2236  // Derivs w.r.t. to nodal velocities
2237  // ---------------------------------
2238  if (element_has_node_with_aux_node_update_fct)
2239  {
2240  // Loop over local nodes
2241  for (unsigned q_local = 0; q_local < n_node; q_local++)
2242  {
2243  // Number of master nodes and storage for the weight of
2244  // the shape function
2245  unsigned n_master2 = 1;
2246  double hang_weight2 = 1.0;
2247  HangInfo* hang_info2_pt = 0;
2248 
2249  // Local boolean to indicate whether the node is hanging
2250  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2251 
2252  Node* actual_shape_controlling_node_pt = node_pt(q_local);
2253 
2254  // If the node is hanging
2255  if (is_node_hanging2)
2256  {
2257  hang_info2_pt = node_pt(q_local)->hanging_pt();
2258 
2259  // Read out number of master nodes from hanging data
2260  n_master2 = hang_info2_pt->nmaster();
2261  }
2262  // Otherwise the node is its own master
2263  else
2264  {
2265  n_master2 = 1;
2266  }
2267 
2268  // Loop over the master nodes
2269  for (unsigned mm = 0; mm < n_master2; mm++)
2270  {
2271  if (is_node_hanging2)
2272  {
2273  actual_shape_controlling_node_pt =
2274  hang_info2_pt->master_node_pt(mm);
2275  hang_weight2 = hang_info2_pt->master_weight(mm);
2276  }
2277 
2278  // Find the corresponding number
2279  unsigned q = local_shape_controlling_node_lookup
2280  [actual_shape_controlling_node_pt];
2281 
2282  // Loop over the two coordinate directions
2283  for (unsigned p = 0; p < 2; p++)
2284  {
2285  double aux = d_U_dX(p, q, 0) *
2286  (psif[q_local] + r * dpsifdx(q_local, 0)) +
2287  d_U_dX(p, q, 1) * r * dpsifdx(q_local, 1);
2288 
2289  // Multiply through by Jacobian, test function and
2290  // integration weight
2291  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2292  aux * testp_ * J * w * hang_weight * hang_weight2;
2293  }
2294  } // End of loop over (mm) master nodes
2295  } // End of loop over local nodes q_local
2296  } // End of if(element_has_node_with_aux_node_update_fct)
2297  } // End of if it's not a boundary condition
2298  } // End of loop over (m) master nodes
2299  } // End of loop over shape functions for continuity eqn
2300 
2301  } // End of loop over integration points
2302 
2303  } // End of get_dresidual_dnodal_coordinates(...)
p
float * p
Definition: Tutorial_Map_using.cpp:9
oomph::AxisymmetricNavierStokesEquations::get_source_fct_gradient
virtual void get_source_fct_gradient(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
Definition: axisym_navier_stokes_elements.h:305
oomph::AxisymmetricNavierStokesEquations::get_body_force_gradient_axi_nst
virtual void get_body_force_gradient_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
Definition: axisym_navier_stokes_elements.h:247
oomph::GeneralisedElement::Default_fd_jacobian_step
static double Default_fd_jacobian_step
Definition: elements.h:1198
oomph::GeneralisedElement::ndof
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
oomph::Node::position_time_stepper_pt
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
oomph::RefineableElement::nshape_controlling_nodes
unsigned nshape_controlling_nodes()
Definition: refineable_elements.h:627
oomph::RefineableElement::shape_controlling_node_lookup
std::map< Node *, unsigned > shape_controlling_node_lookup()
Definition: refineable_elements.h:636
tmp
Eigen::Matrix< Scalar, Dynamic, Dynamic, ColMajor > tmp
Definition: level3_impl.h:365
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286

References oomph::AxisymmetricNavierStokesEquations::ALE_is_disabled, Global_Parameters::body_force(), oomph::GeneralisedElement::Default_fd_jacobian_step, oomph::AxisymmetricNavierStokesEquations::density_ratio(), oomph::FiniteElement::dnodal_position_dt(), oomph::AxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_axi_nst(), oomph::AxisymmetricNavierStokesEquations::du_dt_axi_nst(), G, oomph::AxisymmetricNavierStokesEquations::g(), oomph::AxisymmetricNavierStokesEquations::Gamma, oomph::AxisymmetricNavierStokesEquations::get_body_force_axi_nst(), oomph::AxisymmetricNavierStokesEquations::get_body_force_gradient_axi_nst(), oomph::AxisymmetricNavierStokesEquations::get_source_fct(), oomph::AxisymmetricNavierStokesEquations::get_source_fct_gradient(), oomph::Node::hanging_pt(), oomph::Node::has_auxiliary_node_update_fct_pt(), i, oomph::FiniteElement::integral_pt(), oomph::FiniteElement::interpolated_x(), oomph::Node::is_hanging(), J, j, k, oomph::Integral::knot(), oomph::RefineableElement::local_hang_eqn(), m, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::GeneralisedElement::ndof(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_local_eqn(), oomph::FiniteElement::nodal_position(), oomph::FiniteElement::nodal_value(), oomph::FiniteElement::node_pt(), oomph::AxisymmetricNavierStokesEquations::npres_axi_nst(), oomph::RefineableElement::nshape_controlling_nodes(), oomph::Integral::nweight(), OOMPH_EXCEPTION_LOCATION, p, oomph::AxisymmetricNavierStokesEquations::p_axi_nst(), oomph::AxisymmetricNavierStokesEquations::p_local_eqn(), oomph::AxisymmetricNavierStokesEquations::p_nodal_index_axi_nst(), oomph::Node::perform_auxiliary_node_update_fct(), oomph::Node::position_time_stepper_pt(), pressure_node_pt(), oomph::AxisymmetricNavierStokesEquations::pshape_axi_nst(), Eigen::numext::q, UniformPSDSelfTest::r, oomph::AxisymmetricNavierStokesEquations::re(), oomph::AxisymmetricNavierStokesEquations::re_invfr(), oomph::AxisymmetricNavierStokesEquations::re_invro(), oomph::AxisymmetricNavierStokesEquations::re_st(), s, oomph::RefineableElement::shape_controlling_node_lookup(), TestProblem::source(), oomph::Global_string_for_annotation::string(), oomph::Time::time(), oomph::TimeStepper::time_pt(), oomph::Data::time_stepper_pt(), tmp, oomph::AxisymmetricNavierStokesEquations::u_index_axi_nst(), oomph::Data::value_pt(), oomph::AxisymmetricNavierStokesEquations::viscosity_ratio(), w, oomph::Integral::weight(), oomph::TimeStepper::weight(), and oomph::Node::x().

◆ get_Z2_flux()

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

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

Implements oomph::ElementWithZ2ErrorEstimator.

84  {
85  // Specify the number of velocity dimensions
86  unsigned DIM = 3;
87 
88 #ifdef PARANOID
89  unsigned num_entries = DIM + ((DIM * DIM) - DIM) / 2;
90  if (flux.size() < num_entries)
91  {
92  std::ostringstream error_message;
93  error_message << "The flux vector has the wrong number of entries, "
94  << flux.size() << ", whereas it should be " << num_entries
95  << std::endl;
96  throw OomphLibError(error_message.str(),
97  OOMPH_CURRENT_FUNCTION,
98  OOMPH_EXCEPTION_LOCATION);
99  }
100 #endif
101 
102 
103  // Get strain rate matrix
104  DenseMatrix<double> strainrate(DIM);
105  this->strain_rate(s, strainrate);
106 
107  // Pack into flux Vector
108  unsigned icount = 0;
109 
110  // Start with diagonal terms
111  for (unsigned i = 0; i < DIM; i++)
112  {
113  flux[icount] = strainrate(i, i);
114  icount++;
115  }
116 
117  // Off diagonals row by row
118  for (unsigned i = 0; i < DIM; i++)
119  {
120  for (unsigned j = i + 1; j < DIM; j++)
121  {
122  flux[icount] = strainrate(i, j);
123  icount++;
124  }
125  }
126  }
oomph::AxisymmetricNavierStokesEquations::strain_rate
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
Definition: axisym_navier_stokes_elements.cc:725
ProblemParameters::flux
void flux(const double &time, const Vector< double > &x, double &flux)
Get flux applied along boundary x=0.
Definition: pretend_melt.cc:59

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

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

◆ num_Z2_flux_terms()

unsigned oomph::RefineableAxisymmetricNavierStokesEquations::num_Z2_flux_terms ( )
inlinevirtual

Number of 'flux' terms for Z2 error estimation.

Implements oomph::ElementWithZ2ErrorEstimator.

76  {
77  // 3 diagonal strain rates, 3 off diagonal
78  return 6;
79  }

Referenced by oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::get_Z2_compound_flux_indices(), oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::get_Z2_flux(), and oomph::RefineableBuoyantQAxisymCrouzeixRaviartElement::num_Z2_flux_terms().

◆ pin_elemental_redundant_nodal_pressure_dofs()

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

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

Reimplemented in oomph::RefineableAxisymmetricQTaylorHoodElement.

63 {}

Referenced by pin_redundant_nodal_pressures().

◆ pin_redundant_nodal_pressures()

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

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

  • RefineableAxisymmetricNavierStokesEquations:: pin_all_nodal_pressure_dofs()

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

321  {
322  // Loop over all elements to brutally pin all nodal pressure degrees of
323  // freedom
324  unsigned n_element = element_pt.size();
325  for (unsigned e = 0; e < n_element; e++)
326  {
327  dynamic_cast<RefineableAxisymmetricNavierStokesEquations*>(
328  element_pt[e])
329  ->pin_elemental_redundant_nodal_pressure_dofs();
330  }
331  }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::RefineableAxisymmetricNavierStokesEquations::pin_elemental_redundant_nodal_pressure_dofs
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Definition: refineable_axisym_navier_stokes_elements.h:63

References e(), and pin_elemental_redundant_nodal_pressure_dofs().

◆ pressure_node_pt()

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

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

Reimplemented in oomph::RefineableAxisymmetricQTaylorHoodElement.

54  {
55  return NULL;
56  }

Referenced by fill_in_generic_residual_contribution_axi_nst(), and get_dresidual_dnodal_coordinates().

◆ unpin_all_pressure_dofs()

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

Unpin all pressure dofs in elements listed in vector.

336  {
337  // Loop over all elements to brutally unpin all nodal pressure degrees of
338  // freedom and internal pressure degrees of freedom
339  unsigned n_element = element_pt.size();
340  for (unsigned e = 0; e < n_element; e++)
341  {
342  dynamic_cast<RefineableAxisymmetricNavierStokesEquations*>(
343  element_pt[e])
344  ->unpin_elemental_pressure_dofs();
345  }
346  }
oomph::RefineableAxisymmetricNavierStokesEquations::unpin_elemental_pressure_dofs
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.

References e(), and unpin_elemental_pressure_dofs().

◆ unpin_elemental_pressure_dofs()

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

Unpin all pressure dofs in the element.

Implemented in oomph::RefineableAxisymmetricQCrouzeixRaviartElement, and oomph::RefineableAxisymmetricQTaylorHoodElement.

Referenced by unpin_all_pressure_dofs().

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