oomph::NavierStokesEquations< DIM > Class Template Referenceabstract

#include <navier_stokes_elements.h>

+ Inheritance diagram for oomph::NavierStokesEquations< DIM >:

Public Types

typedef void(* NavierStokesBodyForceFctPt) (const double &time, const Vector< double > &x, Vector< double > &body_force)
 
typedef double(* NavierStokesSourceFctPt) (const double &time, const Vector< double > &x)
 
typedef double(* NavierStokesPressureAdvDiffSourceFctPt) (const Vector< double > &x)
 
- Public Types inherited from oomph::FiniteElement
typedef void(* SteadyExactSolutionFctPt) (const Vector< double > &, Vector< double > &)
 
typedef void(* UnsteadyExactSolutionFctPt) (const double &, const Vector< double > &, Vector< double > &)
 

Public Member Functions

 NavierStokesEquations ()
 
const doublere () const
 Reynolds number. More...
 
const doublere_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
const doubleviscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to Viscosity Ratio. More...
 
const doubledensity_ratio () const
 
double *& density_ratio_pt ()
 Pointer to Density ratio. More...
 
const doublere_invfr () const
 Global inverse Froude number. More...
 
double *& re_invfr_pt ()
 Pointer to global inverse Froude number. More...
 
const Vector< double > & g () const
 Vector of gravitational components. More...
 
Vector< double > *& g_pt ()
 Pointer to Vector of gravitational components. More...
 
NavierStokesBodyForceFctPtbody_force_fct_pt ()
 Access function for the body-force pointer. More...
 
NavierStokesBodyForceFctPt body_force_fct_pt () const
 Access function for the body-force pointer. Const version. More...
 
NavierStokesSourceFctPtsource_fct_pt ()
 Access function for the source-function pointer. More...
 
NavierStokesSourceFctPt source_fct_pt () const
 Access function for the source-function pointer. Const version. More...
 
NavierStokesPressureAdvDiffSourceFctPtsource_fct_for_pressure_adv_diff ()
 
NavierStokesPressureAdvDiffSourceFctPt source_fct_for_pressure_adv_diff () const
 
intpinned_fp_pressure_eqn ()
 
virtual unsigned npres_nst () const =0
 Function to return number of pressure degrees of freedom. More...
 
virtual void pshape_nst (const Vector< double > &s, Shape &psi) const =0
 Compute the pressure shape functions at local coordinate s. More...
 
virtual void pshape_nst (const Vector< double > &s, Shape &psi, Shape &test) const =0
 
virtual double dpshape_and_dptest_eulerian_nst (const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const =0
 
double u_nst (const unsigned &n, const unsigned &i) const
 
double u_nst (const unsigned &t, const unsigned &n, const unsigned &i) const
 
virtual unsigned u_index_nst (const unsigned &i) const
 
unsigned n_u_nst () const
 
double du_dt_nst (const unsigned &n, const unsigned &i) const
 
void disable_ALE ()
 
void enable_ALE ()
 
virtual double p_nst (const unsigned &n_p) const =0
 
virtual double p_nst (const unsigned &t, const unsigned &n_p) const =0
 Pressure at local pressure "node" n_p at time level t. More...
 
virtual void fix_pressure (const unsigned &p_dof, const double &p_value)=0
 Pin p_dof-th pressure dof and set it to value specified by p_value. More...
 
virtual int p_nodal_index_nst () const
 
double pressure_integral () const
 Integral of pressure over element. More...
 
double dissipation () const
 Return integral of dissipation over element. More...
 
double dissipation (const Vector< double > &s) const
 Return dissipation at local coordinate s. More...
 
void get_vorticity (const Vector< double > &s, Vector< double > &vorticity) const
 Compute the vorticity vector at local coordinate s. More...
 
void get_vorticity (const Vector< double > &s, double &vorticity) const
 Compute the scalar vorticity at local coordinate s (2D) More...
 
double kin_energy () const
 Get integral of kinetic energy over element. More...
 
double d_kin_energy_dt () const
 Get integral of time derivative of kinetic energy over element. More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i) More...
 
void get_traction (const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
 
void get_traction (const Vector< double > &s, const Vector< double > &N, Vector< double > &traction_p, Vector< double > &traction_visc_n, Vector< double > &traction_visc_t)
 
void get_load (const Vector< double > &s, const Vector< double > &N, Vector< double > &load)
 
void get_pressure_and_velocity_mass_matrix_diagonal (Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
 
unsigned nscalar_paraview () const
 
void scalar_value_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
 
void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
std::string scalar_name_paraview (const unsigned &i) const
 
void output (std::ostream &outfile)
 
void output (std::ostream &outfile, const unsigned &nplot)
 
void output (FILE *file_pt)
 
void output (FILE *file_pt, const unsigned &nplot)
 
void full_output (std::ostream &outfile)
 
void full_output (std::ostream &outfile, const unsigned &nplot)
 
void output_veloc (std::ostream &outfile, const unsigned &nplot, const unsigned &t)
 
void output_vorticity (std::ostream &outfile, const unsigned &nplot)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 
void compute_norm (double &norm)
 Compute norm of solution: square of the L2 norm of the velocities. More...
 
void compute_norm (Vector< double > &norm)
 Compute the vector norm of the FEM solution. More...
 
void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
void compute_error (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 fill_in_pressure_advection_diffusion_residuals (Vector< double > &residuals)
 
void fill_in_pressure_advection_diffusion_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void pin_all_non_pressure_dofs (std::map< Data *, std::vector< int >> &eqn_number_backup)
 Pin all non-pressure dofs and backup eqn numbers. More...
 
virtual void build_fp_press_adv_diff_robin_bc_element (const unsigned &face_index)=0
 
void output_pressure_advection_diffusion_robin_elements (std::ostream &outfile)
 
void delete_pressure_advection_diffusion_robin_elements ()
 
virtual void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
void interpolated_u_nst (const Vector< double > &s, Vector< double > &veloc) const
 Compute vector of FE interpolated velocity u at local coordinate s. More...
 
double interpolated_u_nst (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated velocity u[i] at local coordinate s. More...
 
double interpolated_u_nst (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void dinterpolated_u_nst_ddata (const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
 
virtual double interpolated_p_nst (const Vector< double > &s) const
 Return FE interpolated pressure at local coordinate s. More...
 
double interpolated_p_nst (const unsigned &t, const Vector< double > &s) const
 Return FE interpolated pressure at local coordinate s at time level t. More...
 
double interpolated_dudx_nst (const Vector< double > &s, const unsigned &i, const unsigned &j) const
 
void point_output_data (const Vector< double > &s, Vector< double > &data)
 
void get_vorticity (const Vector< double > &s, Vector< double > &vorticity) const
 
void get_vorticity (const Vector< double > &s, double &vorticity) const
 
void get_vorticity (const Vector< double > &s, Vector< double > &vorticity) const
 Compute 3D vorticity vector at local coordinate s. More...
 
- Public Member Functions inherited from oomph::FSIFluidElement
 FSIFluidElement ()
 Constructor. More...
 
 FSIFluidElement (const FSIFluidElement &)=delete
 Broken copy constructor. More...
 
void operator= (const FSIFluidElement &)=delete
 Broken assignment operator. More...
 
virtual void identify_load_data (std::set< std::pair< Data *, unsigned >> &paired_load_data)=0
 
virtual void identify_pressure_data (std::set< std::pair< Data *, unsigned >> &paired_pressure_data)=0
 
- Public Member Functions inherited from oomph::FiniteElement
void set_dimension (const unsigned &dim)
 
void set_nodal_dimension (const unsigned &nodal_dim)
 
void set_nnodal_position_type (const unsigned &nposition_type)
 Set the number of types required to interpolate the coordinate. More...
 
void set_n_node (const unsigned &n)
 
int nodal_local_eqn (const unsigned &n, const unsigned &i) const
 
double dJ_eulerian_at_knot (const unsigned &ipt, Shape &psi, DenseMatrix< double > &djacobian_dX) const
 
 FiniteElement ()
 Constructor. More...
 
virtual ~FiniteElement ()
 
 FiniteElement (const FiniteElement &)=delete
 Broken copy constructor. More...
 
virtual bool local_coord_is_valid (const Vector< double > &s)
 Broken assignment operator. More...
 
virtual void move_local_coord_back_into_element (Vector< double > &s) const
 
void get_centre_of_gravity_and_max_radius_in_terms_of_zeta (Vector< double > &cog, double &max_radius) const
 
virtual void local_coordinate_of_node (const unsigned &j, Vector< double > &s) const
 
virtual void local_fraction_of_node (const unsigned &j, Vector< double > &s_fraction)
 
virtual double local_one_d_fraction_of_node (const unsigned &n1d, const unsigned &i)
 
virtual void set_macro_elem_pt (MacroElement *macro_elem_pt)
 
MacroElementmacro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void assign_nodal_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void describe_nodal_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void assign_all_generic_local_eqn_numbers (const bool &store_local_dof_pt)
 
Node *& node_pt (const unsigned &n)
 Return a pointer to the local node n. More...
 
Node *const & node_pt (const unsigned &n) const
 Return a pointer to the local node n (const version) More...
 
unsigned nnode () const
 Return the number of nodes. More...
 
virtual unsigned nnode_1d () const
 
double raw_nodal_position (const unsigned &n, const unsigned &i) const
 
double raw_nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position (const unsigned &n, const unsigned &i) const
 
double nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double dnodal_position_dt (const unsigned &n, const unsigned &i) const
 Return the i-th component of nodal velocity: dx/dt at local node n. More...
 
double dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual unsigned required_nvalue (const unsigned &n) const
 
unsigned nnodal_position_type () const
 
bool has_hanging_nodes () const
 
unsigned nodal_dimension () const
 Return the required Eulerian dimension of the nodes in this element. More...
 
virtual unsigned nvertex_node () const
 
virtual Nodevertex_node_pt (const unsigned &j) const
 
virtual Nodeconstruct_node (const unsigned &n)
 
virtual Nodeconstruct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual Nodeconstruct_boundary_node (const unsigned &n)
 
virtual Nodeconstruct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
int get_node_number (Node *const &node_pt) const
 
virtual Nodeget_node_at_local_coordinate (const Vector< double > &s) const
 
double raw_nodal_value (const unsigned &n, const unsigned &i) const
 
double raw_nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
unsigned dim () const
 
virtual ElementGeometry::ElementGeometry element_geometry () const
 Return the geometry type of the element (either Q or T usually). More...
 
virtual double interpolated_x (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated coordinate x[i] at local coordinate s. More...
 
virtual double interpolated_x (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_x (const Vector< double > &s, Vector< double > &x) const
 Return FE interpolated position x[] at local coordinate s as Vector. More...
 
virtual void interpolated_x (const unsigned &t, const Vector< double > &s, Vector< double > &x) const
 
virtual double interpolated_dxdt (const Vector< double > &s, const unsigned &i, const unsigned &t)
 
virtual void interpolated_dxdt (const Vector< double > &s, const unsigned &t, Vector< double > &dxdt)
 
unsigned ngeom_data () const
 
Datageom_data_pt (const unsigned &j)
 
void position (const Vector< double > &zeta, Vector< double > &r) const
 
void position (const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
 
void dposition_dt (const Vector< double > &zeta, const unsigned &t, Vector< double > &drdt)
 
virtual double zeta_nodal (const unsigned &n, const unsigned &k, const unsigned &i) const
 
void interpolated_zeta (const Vector< double > &s, Vector< double > &zeta) const
 
void locate_zeta (const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
 
virtual void node_update ()
 
virtual void identify_field_data_for_interactions (std::set< std::pair< Data *, unsigned >> &paired_field_data)
 
virtual void identify_geometric_data (std::set< Data * > &geometric_data_pt)
 
virtual double s_min () const
 Min value of local coordinate. More...
 
virtual double s_max () const
 Max. value of local coordinate. More...
 
double size () const
 
virtual double compute_physical_size () const
 
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 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, 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 unsigned ndof_types () const
 
virtual void get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
 
- Public Member Functions inherited from oomph::GeomObject
 GeomObject ()
 Default constructor. More...
 
 GeomObject (const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)
 
 GeomObject (const GeomObject &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const GeomObject &)=delete
 Broken assignment operator. More...
 
virtual ~GeomObject ()
 (Empty) destructor More...
 
unsigned nlagrangian () const
 Access function to # of Lagrangian coordinates. More...
 
unsigned ndim () const
 Access function to # of Eulerian coordinates. More...
 
void set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
 Set # of Lagrangian and Eulerian coordinates. More...
 
TimeStepper *& time_stepper_pt ()
 
TimeSteppertime_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 
- Public Member Functions inherited from oomph::TemplateFreeNavierStokesEquationsBase
 TemplateFreeNavierStokesEquationsBase ()
 Constructor (empty) More...
 
virtual ~TemplateFreeNavierStokesEquationsBase ()
 Virtual destructor (empty) More...
 
virtual int p_local_eqn (const unsigned &n) const =0
 
- Public Member Functions inherited from oomph::NavierStokesElementWithDiagonalMassMatrices
 NavierStokesElementWithDiagonalMassMatrices ()
 Empty constructor. More...
 
virtual ~NavierStokesElementWithDiagonalMassMatrices ()
 Virtual destructor. More...
 
 NavierStokesElementWithDiagonalMassMatrices (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken copy constructor. More...
 
void operator= (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken assignment operator. More...
 

Static Public Attributes

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

Protected Member Functions

virtual double dshape_and_dtest_eulerian_nst (const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_nst (const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
 
virtual double dshape_and_dtest_eulerian_at_knot_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 get_body_force_nst (const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
 
virtual void get_body_force_gradient_nst (const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
 
virtual double get_source_nst (const double &time, const unsigned &ipt, const Vector< double > &x)
 
virtual void get_source_gradient_nst (const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
 
virtual void fill_in_generic_residual_contribution_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 
virtual void fill_in_generic_pressure_advection_diffusion_contribution_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
 
virtual void fill_in_generic_dresidual_contribution_nst (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
 
void fill_in_contribution_to_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
- Protected Member Functions inherited from oomph::FiniteElement
virtual void assemble_local_to_eulerian_jacobian (const DShape &dpsids, DenseMatrix< double > &jacobian) const
 
virtual void assemble_local_to_eulerian_jacobian2 (const DShape &d2psids, DenseMatrix< double > &jacobian2) const
 
virtual void assemble_eulerian_base_vectors (const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
 
template<unsigned DIM>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double invert_jacobian_mapping (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_eulerian_mapping_diagonal (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual void dJ_eulerian_dnodal_coordinates (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<unsigned DIM>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
virtual void d_dshape_eulerian_dnodal_coordinates (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<unsigned DIM>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
virtual void transform_derivatives (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
void transform_derivatives_diagonal (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
virtual void transform_second_derivatives (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
virtual void fill_in_jacobian_from_nodal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void fill_in_jacobian_from_nodal_by_fd (DenseMatrix< double > &jacobian)
 
virtual void update_before_nodal_fd ()
 
virtual void reset_after_nodal_fd ()
 
virtual void update_in_nodal_fd (const unsigned &i)
 
virtual void reset_in_nodal_fd (const unsigned &i)
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Zero-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 One-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Two-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
- Protected Member Functions inherited from oomph::GeneralisedElement
unsigned add_internal_data (Data *const &data_pt, const bool &fd=true)
 
bool internal_data_fd (const unsigned &i) const
 
void exclude_internal_data_fd (const unsigned &i)
 
void include_internal_data_fd (const unsigned &i)
 
void clear_global_eqn_numbers ()
 
void add_global_eqn_numbers (std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)
 
virtual void assign_internal_and_external_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void assign_additional_local_eqn_numbers ()
 
int internal_local_eqn (const unsigned &i, const unsigned &j) const
 
int external_local_eqn (const unsigned &i, const unsigned &j)
 
void fill_in_jacobian_from_internal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_internal_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
virtual void update_before_internal_fd ()
 
virtual void reset_after_internal_fd ()
 
virtual void update_in_internal_fd (const unsigned &i)
 
virtual void reset_in_internal_fd (const unsigned &i)
 
virtual void update_before_external_fd ()
 
virtual void reset_after_external_fd ()
 
virtual void update_in_external_fd (const unsigned &i)
 
virtual void reset_in_external_fd (const unsigned &i)
 
virtual void fill_in_contribution_to_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void fill_in_contribution_to_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void fill_in_contribution_to_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 

Protected Attributes

doubleViscosity_Ratio_pt
 
doubleDensity_Ratio_pt
 
doubleRe_pt
 Pointer to global Reynolds number. More...
 
doubleReSt_pt
 Pointer to global Reynolds number x Strouhal number (=Womersley) More...
 
doubleReInvFr_pt
 
Vector< double > * G_pt
 Pointer to global gravity Vector. More...
 
NavierStokesBodyForceFctPt Body_force_fct_pt
 Pointer to body force function. More...
 
NavierStokesSourceFctPt Source_fct_pt
 Pointer to volumetric source function. More...
 
NavierStokesPressureAdvDiffSourceFctPt Press_adv_diff_source_fct_pt
 
bool ALE_is_disabled
 
Vector< FpPressureAdvDiffRobinBCElementBase * > Pressure_advection_diffusion_robin_element_pt
 
int Pinned_fp_pressure_eqn
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElementMacro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepperGeom_object_time_stepper_pt
 

Static Private Attributes

static int Pressure_not_stored_at_node = -100
 
static double Default_Physical_Constant_Value = 0.0
 Navier–Stokes equations static data. More...
 
static double Default_Physical_Ratio_Value = 1.0
 Navier–Stokes equations static data. More...
 
static Vector< doubleDefault_Gravity_vector
 Static default value for the gravity vector. More...
 

Additional Inherited Members

- Static Protected Attributes inherited from oomph::FiniteElement
static const unsigned Default_Initial_Nvalue = 0
 Default value for the number of values at a node. More...
 
static const double Node_location_tolerance = 1.0e-14
 
static const unsigned N2deriv [] = {0, 1, 3, 6}
 
- Static Protected Attributes inherited from oomph::GeneralisedElement
static DenseMatrix< doubleDummy_matrix
 
static std::deque< double * > Dof_pt_deque
 

Detailed Description

template<unsigned DIM>
class oomph::NavierStokesEquations< DIM >

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// A class for elements that solve the cartesian Navier–Stokes equations, templated by the dimension DIM. This contains the generic maths – any concrete implementation must be derived from this.

We're solving:

\( { Re \left( St \frac{\partial u_i}{\partial t} + (u_j - u_j^{M}) \frac{\partial u_i}{\partial x_j} \right) = - \frac{\partial p}{\partial x_i} - R_\rho B_i(x_j) - \frac{Re}{Fr} G_i + \frac{\partial }{\partial x_j} \left[ R_\mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \right] } \)

and

\( { \frac{\partial u_i}{\partial x_i} = Q } \)

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

Member Typedef Documentation

◆ NavierStokesBodyForceFctPt

template<unsigned DIM>
typedef void(* oomph::NavierStokesEquations< DIM >::NavierStokesBodyForceFctPt) (const double &time, const Vector< double > &x, Vector< double > &body_force)

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

◆ NavierStokesPressureAdvDiffSourceFctPt

template<unsigned DIM>
typedef double(* oomph::NavierStokesEquations< DIM >::NavierStokesPressureAdvDiffSourceFctPt) (const Vector< double > &x)

Function pointer to source function fct(x) for the pressure advection diffusion equation (only used during validation!). x is a Vector!

◆ NavierStokesSourceFctPt

template<unsigned DIM>
typedef double(* oomph::NavierStokesEquations< DIM >::NavierStokesSourceFctPt) (const double &time, const Vector< double > &x)

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

Constructor & Destructor Documentation

◆ NavierStokesEquations()

template<unsigned DIM>
oomph::NavierStokesEquations< DIM >::NavierStokesEquations ( )
inline

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

678  : Body_force_fct_pt(0),
679  Source_fct_pt(0),
681  ALE_is_disabled(false),
683  {
684  // Set all the Physical parameter pointers to the default value zero
689  // Set the Physical ratios to the default value of 1
692  }
double * Viscosity_Ratio_pt
Definition: navier_stokes_elements.h:436
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
Definition: navier_stokes_elements.h:461
double * Re_pt
Pointer to global Reynolds number.
Definition: navier_stokes_elements.h:445
NavierStokesPressureAdvDiffSourceFctPt Press_adv_diff_source_fct_pt
Definition: navier_stokes_elements.h:465
int Pinned_fp_pressure_eqn
Definition: navier_stokes_elements.h:479
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: navier_stokes_elements.h:455
static double Default_Physical_Constant_Value
Navier–Stokes equations static data.
Definition: navier_stokes_elements.h:422
double * ReInvFr_pt
Definition: navier_stokes_elements.h:452
static double Default_Physical_Ratio_Value
Navier–Stokes equations static data.
Definition: navier_stokes_elements.h:426
double * Density_Ratio_pt
Definition: navier_stokes_elements.h:440
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
Definition: navier_stokes_elements.h:429
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
Definition: navier_stokes_elements.h:458
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: navier_stokes_elements.h:448
bool ALE_is_disabled
Definition: navier_stokes_elements.h:470

References oomph::NavierStokesEquations< DIM >::Default_Gravity_vector, oomph::NavierStokesEquations< DIM >::Default_Physical_Constant_Value, oomph::NavierStokesEquations< DIM >::Default_Physical_Ratio_Value, oomph::NavierStokesEquations< DIM >::Density_Ratio_pt, oomph::NavierStokesEquations< DIM >::G_pt, oomph::NavierStokesEquations< DIM >::Re_pt, oomph::NavierStokesEquations< DIM >::ReInvFr_pt, oomph::NavierStokesEquations< DIM >::ReSt_pt, and oomph::NavierStokesEquations< DIM >::Viscosity_Ratio_pt.

Member Function Documentation

◆ body_force_fct_pt() [1/2]

template<unsigned DIM>
NavierStokesBodyForceFctPt& oomph::NavierStokesEquations< DIM >::body_force_fct_pt ( )
inline

Access function for the body-force pointer.

778  {
779  return Body_force_fct_pt;
780  }

References oomph::NavierStokesEquations< DIM >::Body_force_fct_pt.

Referenced by oomph::RefineableNavierStokesEquations< DIM >::further_build().

◆ body_force_fct_pt() [2/2]

template<unsigned DIM>
NavierStokesBodyForceFctPt oomph::NavierStokesEquations< DIM >::body_force_fct_pt ( ) const
inline

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

784  {
785  return Body_force_fct_pt;
786  }

References oomph::NavierStokesEquations< DIM >::Body_force_fct_pt.

◆ build_fp_press_adv_diff_robin_bc_element()

template<unsigned DIM>
virtual void oomph::NavierStokesEquations< DIM >::build_fp_press_adv_diff_robin_bc_element ( const unsigned face_index)
pure virtual

◆ compute_error() [1/4]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::compute_error ( FiniteElement::SteadyExactSolutionFctPt  exact_soln_pt,
double error,
double norm 
)
virtual

Validate against exact solution. Solution is provided via function pointer. Compute L2 error and L2 norm of velocity solution over element.

Validate against exact velocity solution Solution is provided via a function pointer. Compute L2 error and L2 norm of velocity solution over element.

Reimplemented from oomph::FiniteElement.

512  {
513  // Initialise error norm value
514  error = 0.0;
515 
516  // Initialise exact solution norm value
517  norm = 0.0;
518 
519  // Vector of local coordinates
520  Vector<double> s(DIM);
521 
522  // Vector for coordintes
523  Vector<double> x(DIM);
524 
525  // Set the value of n_intpt
526  unsigned n_intpt = integral_pt()->nweight();
527 
528  // Exact solution Vector (u,v,[w],p)
529  Vector<double> exact_soln(DIM + 1);
530 
531  // Loop over the integration points
532  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
533  {
534  // Assign values of s
535  for (unsigned i = 0; i < DIM; i++)
536  {
537  s[i] = integral_pt()->knot(ipt, i);
538  }
539 
540  // Get the integral weight
541  double w = integral_pt()->weight(ipt);
542 
543  // Get jacobian of mapping
544  double J = J_eulerian(s);
545 
546  // Premultiply the weights and the Jacobian
547  double W = w * J;
548 
549  // Get x position as Vector
550  interpolated_x(s, x);
551 
552  // Get exact solution at this point
553  (*exact_soln_pt)(x, exact_soln);
554 
555  // Velocity error
556  for (unsigned i = 0; i < DIM; i++)
557  {
558  norm += exact_soln[i] * exact_soln[i] * W;
559  error += (exact_soln[i] - interpolated_u_nst(s, i)) *
560  (exact_soln[i] - interpolated_u_nst(s, i)) * W;
561  }
562  } // for(unsigned ipt=0;ipt<n_intpt;ipt++)
563  } // End of compute_error
int i
Definition: BiCGSTAB_step_by_step.cpp:9
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
RowVector3d w
Definition: Matrix_resize_int.cpp:3
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual double J_eulerian(const Vector< double > &s) const
Definition: elements.cc:4103
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
void interpolated_u_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
Definition: navier_stokes_elements.h:1505
RealScalar s
Definition: level1_cplx_impl.h:130
#define DIM
Definition: linearised_navier_stokes_elements.h:44
void exact_soln(const double &time, const Vector< double > &x, Vector< double > &soln)
Definition: unstructured_two_d_curved.cc:301
int error
Definition: calibrate.py:297
@ W
Definition: quadtree.h:63
list x
Definition: plotDoE.py:28

References DIM, calibrate::error, ProblemParameters::exact_soln(), i, J, s, w, oomph::QuadTreeNames::W, and plotDoE::x.

◆ compute_error() [2/4]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::compute_error ( FiniteElement::UnsteadyExactSolutionFctPt  exact_soln_pt,
const double time,
double error,
double norm 
)
virtual

Validate against exact solution. Solution is provided via function pointer. Compute L2 error and L2 norm of velocity solution over element.

Validate against exact velocity solution at given time. Solution is provided via function pointer. Compute L2 error and L2 norm of velocity solution over element.

Reimplemented from oomph::FiniteElement.

451  {
452  error = 0.0;
453  norm = 0.0;
454 
455  // Vector of local coordinates
456  Vector<double> s(DIM);
457 
458  // Vector for coordintes
459  Vector<double> x(DIM);
460 
461  // Set the value of n_intpt
462  unsigned n_intpt = integral_pt()->nweight();
463 
464  // Exact solution Vector (u,v,[w],p)
465  Vector<double> exact_soln(DIM + 1);
466 
467  // Loop over the integration points
468  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
469  {
470  // Assign values of s
471  for (unsigned i = 0; i < DIM; i++)
472  {
473  s[i] = integral_pt()->knot(ipt, i);
474  }
475 
476  // Get the integral weight
477  double w = integral_pt()->weight(ipt);
478 
479  // Get jacobian of mapping
480  double J = J_eulerian(s);
481 
482  // Premultiply the weights and the Jacobian
483  double W = w * J;
484 
485  // Get x position as Vector
486  interpolated_x(s, x);
487 
488  // Get exact solution at this point
489  (*exact_soln_pt)(time, x, exact_soln);
490 
491  // Velocity error
492  for (unsigned i = 0; i < DIM; i++)
493  {
494  norm += exact_soln[i] * exact_soln[i] * W;
495  error += (exact_soln[i] - interpolated_u_nst(s, i)) *
496  (exact_soln[i] - interpolated_u_nst(s, i)) * W;
497  }
498  }
499  } // End of compute_error

References DIM, calibrate::error, ProblemParameters::exact_soln(), i, J, s, w, oomph::QuadTreeNames::W, and plotDoE::x.

◆ compute_error() [3/4]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::compute_error ( std::ostream &  outfile,
FiniteElement::SteadyExactSolutionFctPt  exact_soln_pt,
double error,
double norm 
)
virtual

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

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

Reimplemented from oomph::FiniteElement.

373  {
374  error = 0.0;
375  norm = 0.0;
376 
377  // Vector of local coordinates
378  Vector<double> s(DIM);
379 
380  // Vector for coordintes
381  Vector<double> x(DIM);
382 
383  // Set the value of n_intpt
384  unsigned n_intpt = integral_pt()->nweight();
385 
386 
387  outfile << "ZONE" << std::endl;
388 
389  // Exact solution Vector (u,v,[w],p)
390  Vector<double> exact_soln(DIM + 1);
391 
392  // Loop over the integration points
393  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
394  {
395  // Assign values of s
396  for (unsigned i = 0; i < DIM; i++)
397  {
398  s[i] = integral_pt()->knot(ipt, i);
399  }
400 
401  // Get the integral weight
402  double w = integral_pt()->weight(ipt);
403 
404  // Get jacobian of mapping
405  double J = J_eulerian(s);
406 
407  // Premultiply the weights and the Jacobian
408  double W = w * J;
409 
410  // Get x position as Vector
411  interpolated_x(s, x);
412 
413  // Get exact solution at this point
414  (*exact_soln_pt)(x, exact_soln);
415 
416  // Velocity error
417  for (unsigned i = 0; i < DIM; i++)
418  {
419  norm += exact_soln[i] * exact_soln[i] * W;
420  error += (exact_soln[i] - interpolated_u_nst(s, i)) *
421  (exact_soln[i] - interpolated_u_nst(s, i)) * W;
422  }
423 
424  // Output x,y,...,u_exact
425  for (unsigned i = 0; i < DIM; i++)
426  {
427  outfile << x[i] << " ";
428  }
429 
430  // Output x,y,[z],u_error,v_error,[w_error]
431  for (unsigned i = 0; i < DIM; i++)
432  {
433  outfile << exact_soln[i] - interpolated_u_nst(s, i) << " ";
434  }
435  outfile << std::endl;
436  }
437  }

References DIM, calibrate::error, ProblemParameters::exact_soln(), i, J, s, w, oomph::QuadTreeNames::W, and plotDoE::x.

◆ compute_error() [4/4]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::compute_error ( std::ostream &  outfile,
FiniteElement::UnsteadyExactSolutionFctPt  exact_soln_pt,
const double time,
double error,
double norm 
)
virtual

Validate against exact solution at given time Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element

Validate against exact velocity solution at given time. Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element.

Reimplemented from oomph::FiniteElement.

295  {
296  error = 0.0;
297  norm = 0.0;
298 
299  // Vector of local coordinates
300  Vector<double> s(DIM);
301 
302  // Vector for coordintes
303  Vector<double> x(DIM);
304 
305  // Set the value of n_intpt
306  unsigned n_intpt = integral_pt()->nweight();
307 
308  outfile << "ZONE" << std::endl;
309 
310  // Exact solution Vector (u,v,[w],p)
311  Vector<double> exact_soln(DIM + 1);
312 
313  // Loop over the integration points
314  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
315  {
316  // Assign values of s
317  for (unsigned i = 0; i < DIM; i++)
318  {
319  s[i] = integral_pt()->knot(ipt, i);
320  }
321 
322  // Get the integral weight
323  double w = integral_pt()->weight(ipt);
324 
325  // Get jacobian of mapping
326  double J = J_eulerian(s);
327 
328  // Premultiply the weights and the Jacobian
329  double W = w * J;
330 
331  // Get x position as Vector
332  interpolated_x(s, x);
333 
334  // Get exact solution at this point
335  (*exact_soln_pt)(time, x, exact_soln);
336 
337  // Velocity error
338  for (unsigned i = 0; i < DIM; i++)
339  {
340  norm += exact_soln[i] * exact_soln[i] * W;
341  error += (exact_soln[i] - interpolated_u_nst(s, i)) *
342  (exact_soln[i] - interpolated_u_nst(s, i)) * W;
343  }
344 
345  // Output x,y,...,u_exact
346  for (unsigned i = 0; i < DIM; i++)
347  {
348  outfile << x[i] << " ";
349  }
350 
351  // Output x,y,[z],u_error,v_error,[w_error]
352  for (unsigned i = 0; i < DIM; i++)
353  {
354  outfile << exact_soln[i] - interpolated_u_nst(s, i) << " ";
355  }
356 
357  outfile << std::endl;
358  }
359  }

References DIM, calibrate::error, ProblemParameters::exact_soln(), i, J, s, w, oomph::QuadTreeNames::W, and plotDoE::x.

◆ compute_norm() [1/2]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::compute_norm ( double norm)
virtual

Compute norm of solution: square of the L2 norm of the velocities.

Compute norm of the solution.

Reimplemented from oomph::GeneralisedElement.

Reimplemented in oomph::PseudoSolidNodeUpdateElement< TCrouzeixRaviartElement< 2 >, TPVDBubbleEnrichedElement< 2, 3 > >, and oomph::PseudoSolidNodeUpdateElement< TTaylorHoodElement< 2 >, TPVDElement< 2, 3 > >.

187  {
188  // Initialise
189  norm = 0.0;
190 
191  // Vector of local coordinates
192  Vector<double> s(DIM);
193 
194  // Solution
195  double u = 0.0;
196 
197  // Find out how many nodes there are in the element
198  unsigned n_node = this->nnode();
199 
200  Shape psi(n_node);
201 
202  // Set the value of n_intpt
203  unsigned n_intpt = this->integral_pt()->nweight();
204 
205  // Loop over the integration points
206  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
207  {
208  // Assign values of s
209  for (unsigned i = 0; i < DIM; i++)
210  {
211  s[i] = this->integral_pt()->knot(ipt, i);
212  }
213 
214  // Get the integral weight
215  double w = this->integral_pt()->weight(ipt);
216 
217  // Get jacobian of mapping
218  double J = this->J_eulerian(s);
219 
220  // Premultiply the weights and the Jacobian
221  double W = w * J;
222 
223  for (unsigned i = 0; i < DIM; i++)
224  {
225  // Get FE function value
226  u = this->interpolated_u_nst(s, i);
227 
228  // Add to norm
229  norm += u * u * W;
230  }
231  }
232  }
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210

References DIM, i, J, s, w, and oomph::QuadTreeNames::W.

Referenced by oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >::compute_norm(), and oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::compute_norm().

◆ compute_norm() [2/2]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::compute_norm ( Vector< double > &  norm)
virtual

Compute the vector norm of the FEM solution.

Compute the vector norm of FEM solution.

Reimplemented from oomph::GeneralisedElement.

240  {
241  // Resize the solution norm vector
242  norm.resize(DIM + 1, 0.0);
243 
244  // Vector of local coordinates
245  Vector<double> s(DIM, 0.0);
246 
247  // Set the value of n_intpt
248  unsigned n_intpt = integral_pt()->nweight();
249 
250  // Loop over the integration points
251  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
252  {
253  // Assign values of s
254  for (unsigned i = 0; i < DIM; i++)
255  {
256  // Compute the i-th local coordinate of this knot point
257  s[i] = integral_pt()->knot(ipt, i);
258  }
259 
260  // Get the integral weight
261  double w = integral_pt()->weight(ipt);
262 
263  // Get Jacobian of mapping
264  double J = J_eulerian(s);
265 
266  // Premultiply the weights and the Jacobian
267  double W = w * J;
268 
269  // Compute the velocity norm
270  for (unsigned i = 0; i < DIM; i++)
271  {
272  // Update the solution norm value
273  norm[i] += std::pow(interpolated_u_nst(s, i), 2) * W;
274  }
275 
276  // Update the pressure norm value
277  norm[DIM] += std::pow(interpolated_p_nst(s), 2) * W;
278  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
279  } // End of compute_norm
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
Definition: navier_stokes_elements.h:1639
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625

References DIM, i, J, Eigen::bfloat16_impl::pow(), s, w, and oomph::QuadTreeNames::W.

◆ d_kin_energy_dt()

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::d_kin_energy_dt

Get integral of time derivative of kinetic energy over element.

Get integral of time derivative of kinetic energy over element:

1461  {
1462  // Initialise
1463  double d_kin_en_dt = 0.0;
1464 
1465  // Set the value of n_intpt
1466  const unsigned n_intpt = integral_pt()->nweight();
1467 
1468  // Set the Vector to hold local coordinates
1469  Vector<double> s(DIM);
1470 
1471  // Get the number of nodes
1472  const unsigned n_node = this->nnode();
1473 
1474  // Storage for the shape function
1475  Shape psi(n_node);
1476  DShape dpsidx(n_node, DIM);
1477 
1478  // Get the value at which the velocities are stored
1479  unsigned u_index[DIM];
1480  for (unsigned i = 0; i < DIM; i++)
1481  {
1482  u_index[i] = this->u_index_nst(i);
1483  }
1484 
1485  // Loop over the integration points
1486  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1487  {
1488  // Get the jacobian of the mapping
1489  double J = dshape_eulerian_at_knot(ipt, psi, dpsidx);
1490 
1491  // Get the integral weight
1492  double w = integral_pt()->weight(ipt);
1493 
1494  // Now work out the velocity and the time derivative
1495  Vector<double> interpolated_u(DIM, 0.0), interpolated_dudt(DIM, 0.0);
1496 
1497  // Loop over the shape functions
1498  for (unsigned l = 0; l < n_node; l++)
1499  {
1500  // Loop over the dimensions
1501  for (unsigned i = 0; i < DIM; i++)
1502  {
1503  interpolated_u[i] += nodal_value(l, u_index[i]) * psi(l);
1504  interpolated_dudt[i] += du_dt_nst(l, u_index[i]) * psi(l);
1505  }
1506  }
1507 
1508  // Get mesh velocity bit
1509  if (!ALE_is_disabled)
1510  {
1511  Vector<double> mesh_velocity(DIM, 0.0);
1512  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1513 
1514  // Loop over nodes again
1515  for (unsigned l = 0; l < n_node; l++)
1516  {
1517  for (unsigned i = 0; i < DIM; i++)
1518  {
1519  mesh_velocity[i] += this->dnodal_position_dt(l, i) * psi(l);
1520 
1521  for (unsigned j = 0; j < DIM; j++)
1522  {
1523  interpolated_dudx(i, j) +=
1524  this->nodal_value(l, u_index[i]) * dpsidx(l, j);
1525  }
1526  }
1527  }
1528 
1529  // Subtract mesh velocity from du_dt
1530  for (unsigned i = 0; i < DIM; i++)
1531  {
1532  for (unsigned k = 0; k < DIM; k++)
1533  {
1534  interpolated_dudt[i] -= mesh_velocity[k] * interpolated_dudx(i, k);
1535  }
1536  }
1537  }
1538 
1539 
1540  // Loop over directions and add up u du/dt terms
1541  double sum = 0.0;
1542  for (unsigned i = 0; i < DIM; i++)
1543  {
1544  sum += interpolated_u[i] * interpolated_dudt[i];
1545  }
1546 
1547  d_kin_en_dt += sum * w * J;
1548  }
1549 
1550  return d_kin_en_dt;
1551  }
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Definition: elements.cc:3325
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
virtual unsigned u_index_nst(const unsigned &i) const
Definition: navier_stokes_elements.h:866
double du_dt_nst(const unsigned &n, const unsigned &i) const
Definition: navier_stokes_elements.h:880
char char char int int * k
Definition: level2_impl.h:374
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References DIM, i, J, s, and w.

◆ delete_pressure_advection_diffusion_robin_elements()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::delete_pressure_advection_diffusion_robin_elements ( )
inlinevirtual

Delete the FaceElements that apply the Robin boundary condition to the pressure advection diffusion problem required by Fp preconditioner

Implements oomph::TemplateFreeNavierStokesEquationsBase.

1487  {
1488  unsigned nel = Pressure_advection_diffusion_robin_element_pt.size();
1489  for (unsigned e = 0; e < nel; e++)
1490  {
1492  }
1494  }
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Vector< FpPressureAdvDiffRobinBCElementBase * > Pressure_advection_diffusion_robin_element_pt
Definition: navier_stokes_elements.h:475

References e(), and oomph::NavierStokesEquations< DIM >::Pressure_advection_diffusion_robin_element_pt.

◆ density_ratio()

template<unsigned DIM>
const double& oomph::NavierStokesEquations< DIM >::density_ratio ( ) const
inline

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

742  {
743  return *Density_Ratio_pt;
744  }

References oomph::NavierStokesEquations< DIM >::Density_Ratio_pt.

◆ density_ratio_pt()

template<unsigned DIM>
double*& oomph::NavierStokesEquations< DIM >::density_ratio_pt ( )
inline

◆ dinterpolated_u_nst_ddata()

template<unsigned DIM>
virtual void oomph::NavierStokesEquations< DIM >::dinterpolated_u_nst_ddata ( const Vector< double > &  s,
const unsigned i,
Vector< double > &  du_ddata,
Vector< unsigned > &  global_eqn_number 
)
inlinevirtual

Compute the derivatives of the i-th component of velocity at point s with respect to all data that can affect its value. In addition, return the global equation numbers corresponding to the data. The function is virtual so that it can be overloaded in the refineable version

Reimplemented in oomph::RefineableNavierStokesEquations< DIM >.

1591  {
1592  // Find number of nodes
1593  unsigned n_node = nnode();
1594  // Local shape function
1595  Shape psi(n_node);
1596  // Find values of shape function
1597  shape(s, psi);
1598 
1599  // Find the index at which the velocity component is stored
1600  const unsigned u_nodal_index = u_index_nst(i);
1601 
1602  // Find the number of dofs associated with interpolated u
1603  unsigned n_u_dof = 0;
1604  for (unsigned l = 0; l < n_node; l++)
1605  {
1606  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
1607  // If it's positive add to the count
1608  if (global_eqn >= 0)
1609  {
1610  ++n_u_dof;
1611  }
1612  }
1613 
1614  // Now resize the storage schemes
1615  du_ddata.resize(n_u_dof, 0.0);
1616  global_eqn_number.resize(n_u_dof, 0);
1617 
1618  // Loop over th nodes again and set the derivatives
1619  unsigned count = 0;
1620  // Loop over the local nodes and sum
1621  for (unsigned l = 0; l < n_node; l++)
1622  {
1623  // Get the global equation number
1624  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
1625  if (global_eqn >= 0)
1626  {
1627  // Set the global equation number
1628  global_eqn_number[count] = global_eqn;
1629  // Set the derivative with respect to the unknown
1630  du_ddata[count] = psi[l];
1631  // Increase the counter
1632  ++count;
1633  }
1634  }
1635  }
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void shape(const Vector< double > &s, Shape &psi) const =0

References oomph::Data::eqn_number(), i, oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), s, oomph::FiniteElement::shape(), and oomph::NavierStokesEquations< DIM >::u_index_nst().

Referenced by QAdvectionDiffusionElementWithExternalElement< DIM >::get_dwind_adv_diff_dexternal_element_data().

◆ disable_ALE()

◆ dissipation() [1/2]

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::dissipation

Return integral of dissipation over element.

1047  {
1048  // Initialise
1049  double diss = 0.0;
1050 
1051  // Set the value of n_intpt
1052  unsigned n_intpt = integral_pt()->nweight();
1053 
1054  // Set the Vector to hold local coordinates
1055  Vector<double> s(DIM);
1056 
1057  // Loop over the integration points
1058  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1059  {
1060  // Assign values of s
1061  for (unsigned i = 0; i < DIM; i++)
1062  {
1063  s[i] = integral_pt()->knot(ipt, i);
1064  }
1065 
1066  // Get the integral weight
1067  double w = integral_pt()->weight(ipt);
1068 
1069  // Get Jacobian of mapping
1070  double J = J_eulerian(s);
1071 
1072  // Get strain rate matrix
1073  DenseMatrix<double> strainrate(DIM, DIM);
1074  strain_rate(s, strainrate);
1075 
1076  // Initialise
1077  double local_diss = 0.0;
1078  for (unsigned i = 0; i < DIM; i++)
1079  {
1080  for (unsigned j = 0; j < DIM; j++)
1081  {
1082  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
1083  }
1084  }
1085 
1086  diss += local_diss * w * J;
1087  }
1088 
1089  return diss;
1090  }
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
Definition: navier_stokes_elements.cc:1185

References DIM, i, J, j, s, and w.

◆ dissipation() [2/2]

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::dissipation ( const Vector< double > &  s) const

Return dissipation at local coordinate s.

1163  {
1164  // Get strain rate matrix
1165  DenseMatrix<double> strainrate(DIM, DIM);
1166  strain_rate(s, strainrate);
1167 
1168  // Initialise
1169  double local_diss = 0.0;
1170  for (unsigned i = 0; i < DIM; i++)
1171  {
1172  for (unsigned j = 0; j < DIM; j++)
1173  {
1174  local_diss += 2.0 * strainrate(i, j) * strainrate(i, j);
1175  }
1176  }
1177 
1178  return local_diss;
1179  }

◆ dpshape_and_dptest_eulerian_nst()

template<unsigned DIM>
virtual double oomph::NavierStokesEquations< DIM >::dpshape_and_dptest_eulerian_nst ( const Vector< double > &  s,
Shape ppsi,
DShape dppsidx,
Shape ptest,
DShape dptestdx 
) const
pure virtual

◆ dshape_and_dtest_eulerian_at_knot_nst() [1/2]

template<unsigned DIM>
virtual double oomph::NavierStokesEquations< DIM >::dshape_and_dtest_eulerian_at_knot_nst ( const unsigned ipt,
Shape psi,
DShape dpsidx,
RankFourTensor< double > &  d_dpsidx_dX,
Shape test,
DShape dtestdx,
RankFourTensor< double > &  d_dtestdx_dX,
DenseMatrix< double > &  djacobian_dX 
) const
protectedpure virtual

◆ dshape_and_dtest_eulerian_at_knot_nst() [2/2]

template<unsigned DIM>
virtual double oomph::NavierStokesEquations< DIM >::dshape_and_dtest_eulerian_at_knot_nst ( const unsigned ipt,
Shape psi,
DShape dpsidx,
Shape test,
DShape dtestdx 
) const
protectedpure virtual

◆ dshape_and_dtest_eulerian_nst()

template<unsigned DIM>
virtual double oomph::NavierStokesEquations< DIM >::dshape_and_dtest_eulerian_nst ( const Vector< double > &  s,
Shape psi,
DShape dpsidx,
Shape test,
DShape dtestdx 
) const
protectedpure virtual

◆ du_dt_nst()

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::du_dt_nst ( const unsigned n,
const unsigned i 
) const
inline

i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes.

881  {
882  // Get the data's timestepper
883  TimeStepper* time_stepper_pt = this->node_pt(n)->time_stepper_pt();
884 
885  // Initialise dudt
886  double dudt = 0.0;
887 
888  // Loop over the timesteps, if there is a non Steady timestepper
889  if (!time_stepper_pt->is_steady())
890  {
891  // Find the index at which the dof is stored
892  const unsigned u_nodal_index = this->u_index_nst(i);
893 
894  // Number of timsteps (past & present)
895  const unsigned n_time = time_stepper_pt->ntstorage();
896  // Loop over the timesteps
897  for (unsigned t = 0; t < n_time; t++)
898  {
899  dudt +=
900  time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
901  }
902  }
903 
904  return dudt;
905  }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
TimeStepper *& time_stepper_pt()
Definition: geom_objects.h:192
unsigned ntstorage() const
Definition: timesteppers.h:601
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
bool is_steady() const
Definition: timesteppers.h:389
t
Definition: plotPSD.py:36

References oomph::TimeStepper::is_steady(), n, oomph::FiniteElement::nodal_value(), oomph::FiniteElement::node_pt(), oomph::TimeStepper::ntstorage(), plotPSD::t, oomph::GeomObject::time_stepper_pt(), oomph::Data::time_stepper_pt(), oomph::NavierStokesEquations< DIM >::u_index_nst(), and oomph::TimeStepper::weight().

Referenced by oomph::MyTaylorHoodElement< DIM >::output().

◆ enable_ALE()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::enable_ALE ( )
inlinevirtual

(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative. Note: By default, ALE is enabled, at the expense of possibly creating unnecessary work in problems where the mesh is, in fact, stationary.

Reimplemented from oomph::FiniteElement.

919  {
920  ALE_is_disabled = false;
921  }

References oomph::NavierStokesEquations< DIM >::ALE_is_disabled.

Referenced by RefineableBuoyantQCrouzeixRaviartElement< DIM >::enable_ALE(), oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >::enable_ALE(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::enable_ALE(), oomph::BuoyantQCrouzeixRaviartElement< DIM >::enable_ALE(), and oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::enable_ALE().

◆ fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter ( double *const &  parameter_pt,
Vector< double > &  dres_dparam,
DenseMatrix< double > &  djac_dparam,
DenseMatrix< double > &  dmass_matrix_dparam 
)
inlinevirtual

Add the element's contribution to its residuals vector, jacobian matrix and mass matrix

Reimplemented from oomph::GeneralisedElement.

1330  {
1331  // Call the generic routine with the flag set to 2
1333  parameter_pt, dres_dparam, djac_dparam, dmass_matrix_dparam, 2);
1334  }
virtual void fill_in_generic_dresidual_contribution_nst(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
Definition: navier_stokes_elements.cc:2087

References oomph::NavierStokesEquations< DIM >::fill_in_generic_dresidual_contribution_nst().

◆ fill_in_contribution_to_djacobian_dparameter()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_djacobian_dparameter ( double *const &  parameter_pt,
Vector< double > &  dres_dparam,
DenseMatrix< double > &  djac_dparam 
)
inlinevirtual

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

Reimplemented from oomph::GeneralisedElement.

1313  {
1314  // Call the generic routine with the flag set to 1
1316  parameter_pt,
1317  dres_dparam,
1318  djac_dparam,
1320  1);
1321  }
static DenseMatrix< double > Dummy_matrix
Definition: elements.h:227

References oomph::GeneralisedElement::Dummy_matrix, and oomph::NavierStokesEquations< DIM >::fill_in_generic_dresidual_contribution_nst().

◆ fill_in_contribution_to_dresiduals_dparameter()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_dresiduals_dparameter ( double *const &  parameter_pt,
Vector< double > &  dres_dparam 
)
inlinevirtual

Compute the element's residual Vector.

Reimplemented from oomph::GeneralisedElement.

1296  {
1297  // Call the generic residuals function with flag set to 0
1298  // and using a dummy matrix argument
1300  parameter_pt,
1301  dres_dparam,
1304  0);
1305  }

References oomph::GeneralisedElement::Dummy_matrix, and oomph::NavierStokesEquations< DIM >::fill_in_generic_dresidual_contribution_nst().

◆ fill_in_contribution_to_hessian_vector_products()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_hessian_vector_products ( Vector< double > const &  Y,
DenseMatrix< double > const &  C,
DenseMatrix< double > &  product 
)
protectedvirtual

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

Reimplemented from oomph::GeneralisedElement.

2109  {
2110  throw OomphLibError("Not yet implemented\n",
2113  }
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ fill_in_contribution_to_jacobian()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_jacobian ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian 
)
inlinevirtual

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

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::PseudoSolidNodeUpdateElement< TCrouzeixRaviartElement< 2 >, TPVDBubbleEnrichedElement< 2, 3 > >, and oomph::PseudoSolidNodeUpdateElement< TTaylorHoodElement< 2 >, TPVDElement< 2, 3 > >.

1275  {
1276  // Call the generic routine with the flag set to 1
1278  residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
1279  }
virtual void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Definition: navier_stokes_elements.cc:1774

References oomph::GeneralisedElement::Dummy_matrix, and oomph::NavierStokesEquations< DIM >::fill_in_generic_residual_contribution_nst().

Referenced by RefineableBuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_jacobian(), RefineableQCrouzeixRaviartElementWithExternalElement< DIM >::fill_in_contribution_to_jacobian(), oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_jacobian(), RefineableQCrouzeixRaviartElementWithTwoExternalElement< DIM >::fill_in_contribution_to_jacobian(), QCrouzeixRaviartElementWithTwoExternalElement< DIM >::fill_in_contribution_to_jacobian(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_jacobian(), oomph::BuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_jacobian(), and oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_jacobian().

◆ fill_in_contribution_to_jacobian_and_mass_matrix()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_jacobian_and_mass_matrix ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
DenseMatrix< double > &  mass_matrix 
)
inlinevirtual

◆ fill_in_contribution_to_residuals()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_residuals ( Vector< double > &  residuals)
inlinevirtual

Compute the element's residual Vector.

Reimplemented from oomph::GeneralisedElement.

Reimplemented in oomph::PseudoSolidNodeUpdateElement< TCrouzeixRaviartElement< 2 >, TPVDBubbleEnrichedElement< 2, 3 > >, and oomph::PseudoSolidNodeUpdateElement< TTaylorHoodElement< 2 >, TPVDElement< 2, 3 > >.

1261  {
1262  // Call the generic residuals function with flag set to 0
1263  // and using a dummy matrix argument
1265  residuals,
1268  0);
1269  }

References oomph::GeneralisedElement::Dummy_matrix, and oomph::NavierStokesEquations< DIM >::fill_in_generic_residual_contribution_nst().

Referenced by RefineableBuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_residuals(), RefineableQCrouzeixRaviartElementWithExternalElement< DIM >::fill_in_contribution_to_residuals(), oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_residuals(), RefineableQCrouzeixRaviartElementWithTwoExternalElement< DIM >::fill_in_contribution_to_residuals(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_residuals(), oomph::BuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_residuals(), oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::fill_in_contribution_to_residuals(), oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >::fill_in_off_diagonal_jacobian_blocks_by_fd(), and oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::fill_in_off_diagonal_jacobian_blocks_by_fd().

◆ fill_in_generic_dresidual_contribution_nst()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_generic_dresidual_contribution_nst ( double *const &  parameter_pt,
Vector< double > &  dres_dparam,
DenseMatrix< double > &  djac_dparam,
DenseMatrix< double > &  dmass_matrix_dparam,
unsigned  flag 
)
protectedvirtual

Compute the derivatives of the residuals for the Navier–Stokes equations with respect to a parameter Flag=1 (or 0): do (or don't) compute the Jacobian as well. Flag=2: Fill in mass matrix too.

Compute the derivatives of the residuals for the Navier–Stokes equations with respect to a parameter; flag=2 or 1 or 0: do (or don't) compute the Jacobian and mass matrix as well

2093  {
2094  throw OomphLibError("Not yet implemented\n",
2097  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

Referenced by oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter(), oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_djacobian_dparameter(), and oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_dresiduals_dparameter().

◆ fill_in_generic_pressure_advection_diffusion_contribution_nst()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_generic_pressure_advection_diffusion_contribution_nst ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
unsigned  flag 
)
protectedvirtual

Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp preconditioner. flag=1(or 0): do (or don't) compute the Jacobian as well.

Reimplemented in oomph::RefineableNavierStokesEquations< DIM >.

1607  {
1608  // Return immediately if there are no dofs
1609  if (ndof() == 0) return;
1610 
1611  // Find out how many nodes there are
1612  unsigned n_node = nnode();
1613 
1614  // Find out how many pressure dofs there are
1615  unsigned n_pres = npres_nst();
1616 
1617  // Find the indices at which the local velocities are stored
1618  unsigned u_nodal_index[DIM];
1619  for (unsigned i = 0; i < DIM; i++)
1620  {
1621  u_nodal_index[i] = u_index_nst(i);
1622  }
1623 
1624  // Set up memory for the velocity shape fcts
1625  Shape psif(n_node);
1626  DShape dpsidx(n_node, DIM);
1627 
1628  // Set up memory for pressure shape and test functions
1629  Shape psip(n_pres), testp(n_pres);
1630  DShape dpsip(n_pres, DIM);
1631  DShape dtestp(n_pres, DIM);
1632 
1633  // Number of integration points
1634  unsigned n_intpt = integral_pt()->nweight();
1635 
1636  // Set the Vector to hold local coordinates
1637  Vector<double> s(DIM);
1638 
1639  // Get Physical Variables from Element
1640  // Reynolds number must be multiplied by the density ratio
1641  double scaled_re = re() * density_ratio();
1642 
1643  // Integers to store the local equations and unknowns
1644  int local_eqn = 0, local_unknown = 0;
1645 
1646  // Loop over the integration points
1647  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1648  {
1649  // Assign values of s
1650  for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
1651 
1652  // Get the integral weight
1653  double w = integral_pt()->weight(ipt);
1654 
1655  // Call the derivatives of the veloc shape functions
1656  // (Derivs not needed but they are free)
1657  double J = this->dshape_eulerian_at_knot(ipt, psif, dpsidx);
1658 
1659  // Call the pressure shape and test functions
1660  this->dpshape_and_dptest_eulerian_nst(s, psip, dpsip, testp, dtestp);
1661 
1662  // Premultiply the weights and the Jacobian
1663  double W = w * J;
1664 
1665  // Calculate local values of the pressure and velocity components
1666  // Allocate
1667  Vector<double> interpolated_u(DIM, 0.0);
1668  Vector<double> interpolated_x(DIM, 0.0);
1669  Vector<double> interpolated_dpdx(DIM, 0.0);
1670 
1671  // Calculate pressure gradient
1672  for (unsigned l = 0; l < n_pres; l++)
1673  {
1674  for (unsigned i = 0; i < DIM; i++)
1675  {
1676  interpolated_dpdx[i] += p_nst(l) * dpsip(l, i);
1677  }
1678  }
1679 
1680  // Calculate velocities
1681 
1682  // Loop over nodes
1683  for (unsigned l = 0; l < n_node; l++)
1684  {
1685  // Loop over directions
1686  for (unsigned i = 0; i < DIM; i++)
1687  {
1688  // Get the nodal value
1689  double u_value = raw_nodal_value(l, u_nodal_index[i]);
1690  interpolated_u[i] += u_value * psif[l];
1691  interpolated_x[i] += raw_nodal_position(l, i) * psif[l];
1692  }
1693  }
1694 
1695  // Source function (for validaton only)
1696  double source = 0.0;
1698  {
1700  }
1701 
1702  // Loop over the shape functions
1703  for (unsigned l = 0; l < n_pres; l++)
1704  {
1705  local_eqn = p_local_eqn(l);
1706 
1707  // If not a boundary conditions
1708  if (local_eqn >= 0)
1709  {
1710  residuals[local_eqn] -= source * testp[l] * W;
1711  for (unsigned k = 0; k < DIM; k++)
1712  {
1713  residuals[local_eqn] +=
1714  interpolated_dpdx[k] *
1715  (scaled_re * interpolated_u[k] * testp[l] + dtestp(l, k)) * W;
1716  }
1717 
1718  // Jacobian too?
1719  if (flag)
1720  {
1721  // Loop over the shape functions
1722  for (unsigned l2 = 0; l2 < n_pres; l2++)
1723  {
1724  local_unknown = p_local_eqn(l2);
1725 
1726  // If not a boundary conditions
1727  if (local_unknown >= 0)
1728  {
1729  if ((int(eqn_number(local_eqn)) != Pinned_fp_pressure_eqn) &&
1730  (int(eqn_number(local_unknown)) != Pinned_fp_pressure_eqn))
1731  {
1732  for (unsigned k = 0; k < DIM; k++)
1733  {
1734  jacobian(local_eqn, local_unknown) +=
1735  dtestp(l2, k) *
1736  (scaled_re * interpolated_u[k] * testp[l] +
1737  dtestp(l, k)) *
1738  W;
1739  }
1740  }
1741  else
1742  {
1743  if ((int(eqn_number(local_eqn)) == Pinned_fp_pressure_eqn) &&
1744  (int(eqn_number(local_unknown)) ==
1746  {
1747  jacobian(local_eqn, local_unknown) = 1.0;
1748  }
1749  }
1750  }
1751  }
1752  } /*End of Jacobian calculation*/
1753  } // End of if not boundary condition
1754  } // End of loop over l
1755  }
1756 
1757  // Now add boundary contributions from Robin BCs
1758  unsigned nrobin = Pressure_advection_diffusion_robin_element_pt.size();
1759  for (unsigned e = 0; e < nrobin; e++)
1760  {
1762  ->fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
1763  residuals, jacobian, flag);
1764  }
1765  }
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2576
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.cc:1686
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
unsigned long eqn_number(const unsigned &ieqn_local) const
Definition: elements.h:704
virtual double dpshape_and_dptest_eulerian_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const =0
virtual unsigned npres_nst() const =0
Function to return number of pressure degrees of freedom.
const double & re() const
Reynolds number.
Definition: navier_stokes_elements.h:703
virtual double p_nst(const unsigned &n_p) const =0
const double & density_ratio() const
Definition: navier_stokes_elements.h:741
virtual int p_local_eqn(const unsigned &n) const =0
return int(ret)+1
void source(const Vector< double > &x, Vector< double > &f)
Source function.
Definition: unstructured_two_d_circle.cc:46

References DIM, e(), i, int(), J, k, s, TestProblem::source(), w, and oomph::QuadTreeNames::W.

Referenced by oomph::NavierStokesEquations< DIM >::fill_in_pressure_advection_diffusion_jacobian(), and oomph::NavierStokesEquations< DIM >::fill_in_pressure_advection_diffusion_residuals().

◆ fill_in_generic_residual_contribution_nst()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_generic_residual_contribution_nst ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
DenseMatrix< double > &  mass_matrix,
unsigned  flag 
)
protectedvirtual

Compute the residuals for the Navier–Stokes equations. Flag=1 (or 0): do (or don't) compute the Jacobian as well. Flag=2: Fill in mass matrix too.

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

Reimplemented in oomph::RefineableNavierStokesEquations< DIM >.

1779  {
1780  // Return immediately if there are no dofs
1781  if (ndof() == 0) return;
1782 
1783  // Find out how many nodes there are
1784  unsigned n_node = nnode();
1785 
1786  // Get continuous time from timestepper of first node
1787  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1788 
1789  // Find out how many pressure dofs there are
1790  unsigned n_pres = npres_nst();
1791 
1792  // Find the indices at which the local velocities are stored
1793  unsigned u_nodal_index[DIM];
1794  for (unsigned i = 0; i < DIM; i++)
1795  {
1796  u_nodal_index[i] = u_index_nst(i);
1797  }
1798 
1799  // Set up memory for the shape and test functions
1800  Shape psif(n_node), testf(n_node);
1801  DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
1802 
1803  // Set up memory for pressure shape and test functions
1804  Shape psip(n_pres), testp(n_pres);
1805 
1806  // Number of integration points
1807  unsigned n_intpt = integral_pt()->nweight();
1808 
1809  // Set the Vector to hold local coordinates
1810  Vector<double> s(DIM);
1811 
1812  // Get Physical Variables from Element
1813  // Reynolds number must be multiplied by the density ratio
1814  double scaled_re = re() * density_ratio();
1815  double scaled_re_st = re_st() * density_ratio();
1816  double scaled_re_inv_fr = re_invfr() * density_ratio();
1817  double visc_ratio = viscosity_ratio();
1818  Vector<double> G = g();
1819 
1820  // Integers to store the local equations and unknowns
1821  int local_eqn = 0, local_unknown = 0;
1822 
1823  // Loop over the integration points
1824  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1825  {
1826  // Assign values of s
1827  for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
1828  // Get the integral weight
1829  double w = integral_pt()->weight(ipt);
1830 
1831  // Call the derivatives of the shape and test functions
1833  ipt, psif, dpsifdx, testf, dtestfdx);
1834 
1835  // Call the pressure shape and test functions
1836  pshape_nst(s, psip, testp);
1837 
1838  // Premultiply the weights and the Jacobian
1839  double W = w * J;
1840 
1841  // Calculate local values of the pressure and velocity components
1842  // Allocate
1843  double interpolated_p = 0.0;
1844  Vector<double> interpolated_u(DIM, 0.0);
1845  Vector<double> interpolated_x(DIM, 0.0);
1846  Vector<double> mesh_velocity(DIM, 0.0);
1847  Vector<double> dudt(DIM, 0.0);
1848  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1849 
1850  // Calculate pressure
1851  for (unsigned l = 0; l < n_pres; l++)
1852  interpolated_p += p_nst(l) * psip[l];
1853 
1854  // Calculate velocities and derivatives:
1855 
1856  // Loop over nodes
1857  for (unsigned l = 0; l < n_node; l++)
1858  {
1859  // Loop over directions
1860  for (unsigned i = 0; i < DIM; i++)
1861  {
1862  // Get the nodal value
1863  double u_value = raw_nodal_value(l, u_nodal_index[i]);
1864  interpolated_u[i] += u_value * psif[l];
1865  interpolated_x[i] += raw_nodal_position(l, i) * psif[l];
1866  dudt[i] += du_dt_nst(l, i) * psif[l];
1867 
1868  // Loop over derivative directions
1869  for (unsigned j = 0; j < DIM; j++)
1870  {
1871  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1872  }
1873  }
1874  }
1875 
1876  if (!ALE_is_disabled)
1877  {
1878  // Loop over nodes
1879  for (unsigned l = 0; l < n_node; l++)
1880  {
1881  // Loop over directions
1882  for (unsigned i = 0; i < DIM; i++)
1883  {
1884  mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif[l];
1885  }
1886  }
1887  }
1888 
1889  // Get the user-defined body force terms
1890  Vector<double> body_force(DIM);
1892 
1893  // Get the user-defined source function
1894  double source = get_source_nst(time, ipt, interpolated_x);
1895 
1896 
1897  // MOMENTUM EQUATIONS
1898  //------------------
1899 
1900  // Loop over the test functions
1901  for (unsigned l = 0; l < n_node; l++)
1902  {
1903  // Loop over the velocity components
1904  for (unsigned i = 0; i < DIM; i++)
1905  {
1906  /*IF it's not a boundary condition*/
1907  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
1908  if (local_eqn >= 0)
1909  {
1910  // Add the user-defined body force terms
1911  residuals[local_eqn] += body_force[i] * testf[l] * W;
1912 
1913  // Add the gravitational body force term
1914  residuals[local_eqn] += scaled_re_inv_fr * testf[l] * G[i] * W;
1915 
1916  // Add the pressure gradient term
1917  residuals[local_eqn] += interpolated_p * dtestfdx(l, i) * W;
1918 
1919  // Add in the stress tensor terms
1920  // The viscosity ratio needs to go in here to ensure
1921  // continuity of normal stress is satisfied even in flows
1922  // with zero pressure gradient!
1923  for (unsigned k = 0; k < DIM; k++)
1924  {
1925  residuals[local_eqn] -=
1926  visc_ratio *
1927  (interpolated_dudx(i, k) + Gamma[i] * interpolated_dudx(k, i)) *
1928  dtestfdx(l, k) * W;
1929  }
1930 
1931  // Add in the inertial terms
1932  // du/dt term
1933  residuals[local_eqn] -= scaled_re_st * dudt[i] * testf[l] * W;
1934 
1935 
1936  // Convective terms, including mesh velocity
1937  for (unsigned k = 0; k < DIM; k++)
1938  {
1939  double tmp = scaled_re * interpolated_u[k];
1940  if (!ALE_is_disabled) tmp -= scaled_re_st * mesh_velocity[k];
1941  residuals[local_eqn] -=
1942  tmp * interpolated_dudx(i, k) * testf[l] * W;
1943  }
1944 
1945  // CALCULATE THE JACOBIAN
1946  if (flag)
1947  {
1948  // Loop over the velocity shape functions again
1949  for (unsigned l2 = 0; l2 < n_node; l2++)
1950  {
1951  // Loop over the velocity components again
1952  for (unsigned i2 = 0; i2 < DIM; i2++)
1953  {
1954  // If at a non-zero degree of freedom add in the entry
1955  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
1956  if (local_unknown >= 0)
1957  {
1958  // Add contribution to Elemental Matrix
1959  jacobian(local_eqn, local_unknown) -=
1960  visc_ratio * Gamma[i] * dpsifdx(l2, i) * dtestfdx(l, i2) *
1961  W;
1962 
1963  // Extra component if i2 = i
1964  if (i2 == i)
1965  {
1966  /*Loop over velocity components*/
1967  for (unsigned k = 0; k < DIM; k++)
1968  {
1969  jacobian(local_eqn, local_unknown) -=
1970  visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) * W;
1971  }
1972  }
1973 
1974  // Now add in the inertial terms
1975  jacobian(local_eqn, local_unknown) -=
1976  scaled_re * psif[l2] * interpolated_dudx(i, i2) *
1977  testf[l] * W;
1978 
1979  // Extra component if i2=i
1980  if (i2 == i)
1981  {
1982  // Add the mass matrix term (only diagonal entries)
1983  // Note that this is positive because the mass matrix
1984  // is taken to the other side of the equation when
1985  // formulating the generalised eigenproblem.
1986  if (flag == 2)
1987  {
1988  mass_matrix(local_eqn, local_unknown) +=
1989  scaled_re_st * psif[l2] * testf[l] * W;
1990  }
1991 
1992  // du/dt term
1993  jacobian(local_eqn, local_unknown) -=
1994  scaled_re_st *
1995  node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1996  psif[l2] * testf[l] * W;
1997 
1998  // Loop over the velocity components
1999  for (unsigned k = 0; k < DIM; k++)
2000  {
2001  double tmp = scaled_re * interpolated_u[k];
2002  if (!ALE_is_disabled)
2003  tmp -= scaled_re_st * mesh_velocity[k];
2004  jacobian(local_eqn, local_unknown) -=
2005  tmp * dpsifdx(l2, k) * testf[l] * W;
2006  }
2007  }
2008  }
2009  }
2010  }
2011 
2012  /*Now loop over pressure shape functions*/
2013  /*This is the contribution from pressure gradient*/
2014  for (unsigned l2 = 0; l2 < n_pres; l2++)
2015  {
2016  /*If we are at a non-zero degree of freedom in the entry*/
2017  local_unknown = p_local_eqn(l2);
2018  if (local_unknown >= 0)
2019  {
2020  jacobian(local_eqn, local_unknown) +=
2021  psip[l2] * dtestfdx(l, i) * W;
2022  }
2023  }
2024  } /*End of Jacobian calculation*/
2025 
2026  } // End of if not boundary condition statement
2027 
2028  } // End of loop over dimension
2029  } // End of loop over shape functions
2030 
2031 
2032  // CONTINUITY EQUATION
2033  //-------------------
2034 
2035  // Loop over the shape functions
2036  for (unsigned l = 0; l < n_pres; l++)
2037  {
2038  local_eqn = p_local_eqn(l);
2039  // If not a boundary conditions
2040  if (local_eqn >= 0)
2041  {
2042  // Source term
2043  // residuals[local_eqn] -=source*testp[l]*W;
2044  double aux = -source;
2045 
2046  // Loop over velocity components
2047  for (unsigned k = 0; k < DIM; k++)
2048  {
2049  // residuals[local_eqn] += interpolated_dudx(k,k)*testp[l]*W;
2050  aux += interpolated_dudx(k, k);
2051  }
2052 
2053  residuals[local_eqn] += aux * testp[l] * W;
2054 
2055  /*CALCULATE THE JACOBIAN*/
2056  if (flag)
2057  {
2058  /*Loop over the velocity shape functions*/
2059  for (unsigned l2 = 0; l2 < n_node; l2++)
2060  {
2061  /*Loop over velocity components*/
2062  for (unsigned i2 = 0; i2 < DIM; i2++)
2063  {
2064  /*If we're at a non-zero degree of freedom add it in*/
2065  local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2066  if (local_unknown >= 0)
2067  {
2068  jacobian(local_eqn, local_unknown) +=
2069  dpsifdx(l2, i2) * testp[l] * W;
2070  }
2071  } /*End of loop over i2*/
2072  } /*End of loop over l2*/
2073  } /*End of Jacobian calculation*/
2074 
2075  } // End of if not boundary condition
2076  } // End of loop over l
2077  }
2078  }
JacobiRotation< float > G
Definition: Jacobi_makeGivens.cpp:2
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
double raw_dnodal_position_dt(const unsigned &n, const unsigned &i) const
Definition: elements.h:2256
virtual double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
virtual void pshape_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: navier_stokes_elements.h:709
virtual void get_body_force_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Definition: navier_stokes_elements.h:517
const double & re_invfr() const
Global inverse Froude number.
Definition: navier_stokes_elements.h:753
const double & viscosity_ratio() const
Definition: navier_stokes_elements.h:728
virtual double get_source_nst(const double &time, const unsigned &ipt, const Vector< double > &x)
Definition: navier_stokes_elements.h:585
const Vector< double > & g() const
Vector of gravitational components.
Definition: navier_stokes_elements.h:765
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition: navier_stokes_elements.h:698
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
Eigen::Matrix< Scalar, Dynamic, Dynamic, ColMajor > tmp
Definition: level3_impl.h:365
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Definition: axisym_linear_elasticity/cylinder/cylinder.cc:96

References oomph::ALE_is_disabled, Global_Parameters::body_force(), DIM, G, GlobalPhysicalVariables::Gamma, i, J, j, k, s, TestProblem::source(), tmp, w, and oomph::QuadTreeNames::W.

Referenced by oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_jacobian(), oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_jacobian_and_mass_matrix(), and oomph::NavierStokesEquations< DIM >::fill_in_contribution_to_residuals().

◆ fill_in_pressure_advection_diffusion_jacobian()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_pressure_advection_diffusion_jacobian ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian 
)
inlinevirtual

Compute the residuals and Jacobian for the associated pressure advection diffusion problem. Used by the Fp preconditioner.

Implements oomph::TemplateFreeNavierStokesEquationsBase.

1350  {
1352  residuals, jacobian, 1);
1353  }
virtual void fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Definition: navier_stokes_elements.cc:1605

References oomph::NavierStokesEquations< DIM >::fill_in_generic_pressure_advection_diffusion_contribution_nst().

◆ fill_in_pressure_advection_diffusion_residuals()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::fill_in_pressure_advection_diffusion_residuals ( Vector< double > &  residuals)
inlinevirtual

Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp preconditioner.

Implements oomph::TemplateFreeNavierStokesEquationsBase.

References oomph::GeneralisedElement::Dummy_matrix, and oomph::NavierStokesEquations< DIM >::fill_in_generic_pressure_advection_diffusion_contribution_nst().

◆ fix_pressure()

template<unsigned DIM>
virtual void oomph::NavierStokesEquations< DIM >::fix_pressure ( const unsigned p_dof,
const double p_value 
)
pure virtual

◆ full_output() [1/2]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::full_output ( std::ostream &  outfile)
inline

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

1181  {
1182  unsigned nplot = 5;
1183  full_output(outfile, nplot);
1184  }
void full_output(std::ostream &outfile)
Definition: navier_stokes_elements.h:1180

Referenced by oomph::QCrouzeixRaviartElement< DIM >::full_output(), oomph::TCrouzeixRaviartElement< DIM >::full_output(), and STSpineMesh< ELEMENT, INTERFACE_ELEMENT >::full_output().

◆ full_output() [2/2]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::full_output ( std::ostream &  outfile,
const unsigned nplot 
)

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

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

849  {
850  // Vector of local coordinates
851  Vector<double> s(DIM);
852 
853  // Acceleration
854  Vector<double> dudt(DIM);
855 
856  // Mesh elocity
857  Vector<double> mesh_veloc(DIM, 0.0);
858 
859  // Tecplot header info
860  outfile << tecplot_zone_string(nplot);
861 
862  // Find out how many nodes there are
863  unsigned n_node = nnode();
864 
865  // Set up memory for the shape functions
866  Shape psif(n_node);
867  DShape dpsifdx(n_node, DIM);
868 
869  // Loop over plot points
870  unsigned num_plot_points = nplot_points(nplot);
871  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
872  {
873  // Get local coordinates of plot point
874  get_s_plot(iplot, nplot, s);
875 
876  // Call the derivatives of the shape and test functions
877  dshape_eulerian(s, psif, dpsifdx);
878 
879  // Allocate storage
880  Vector<double> mesh_veloc(DIM);
881  Vector<double> dudt(DIM);
882  Vector<double> dudt_ALE(DIM);
883  DenseMatrix<double> interpolated_dudx(DIM, DIM);
884 
885  // Initialise everything to zero
886  for (unsigned i = 0; i < DIM; i++)
887  {
888  mesh_veloc[i] = 0.0;
889  dudt[i] = 0.0;
890  dudt_ALE[i] = 0.0;
891  for (unsigned j = 0; j < DIM; j++)
892  {
893  interpolated_dudx(i, j) = 0.0;
894  }
895  }
896 
897  // Calculate velocities and derivatives
898 
899 
900  // Loop over directions
901  for (unsigned i = 0; i < DIM; i++)
902  {
903  // Get the index at which velocity i is stored
904  unsigned u_nodal_index = u_index_nst(i);
905  // Loop over nodes
906  for (unsigned l = 0; l < n_node; l++)
907  {
908  dudt[i] += du_dt_nst(l, u_nodal_index) * psif[l];
909  mesh_veloc[i] += dnodal_position_dt(l, i) * psif[l];
910 
911  // Loop over derivative directions for velocity gradients
912  for (unsigned j = 0; j < DIM; j++)
913  {
914  interpolated_dudx(i, j) +=
915  nodal_value(l, u_nodal_index) * dpsifdx(l, j);
916  }
917  }
918  }
919 
920 
921  // Get dudt in ALE form (incl mesh veloc)
922  for (unsigned i = 0; i < DIM; i++)
923  {
924  dudt_ALE[i] = dudt[i];
925  for (unsigned k = 0; k < DIM; k++)
926  {
927  dudt_ALE[i] -= mesh_veloc[k] * interpolated_dudx(i, k);
928  }
929  }
930 
931 
932  // Coordinates
933  for (unsigned i = 0; i < DIM; i++)
934  {
935  outfile << interpolated_x(s, i) << " ";
936  }
937 
938  // Velocities
939  for (unsigned i = 0; i < DIM; i++)
940  {
941  outfile << interpolated_u_nst(s, i) << " ";
942  }
943 
944  // Pressure
945  outfile << interpolated_p_nst(s) << " ";
946 
947  // Accelerations
948  for (unsigned i = 0; i < DIM; i++)
949  {
950  outfile << dudt_ALE[i] << " ";
951  }
952 
953  // // Mesh velocity
954  // for(unsigned i=0;i<DIM;i++)
955  // {
956  // outfile << mesh_veloc[i] << " ";
957  // }
958 
959  // Dissipation
960  outfile << dissipation(s) << " ";
961 
962 
963  outfile << std::endl;
964  }
965 
966  // Write tecplot footer (e.g. FE connectivity lists)
967  write_tecplot_zone_footer(outfile, nplot);
968  }
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Definition: elements.h:3161
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Definition: elements.h:3186
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Definition: elements.h:3174
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Definition: elements.cc:3298
double dissipation() const
Return integral of dissipation over element.
Definition: navier_stokes_elements.cc:1046

References DIM, i, j, k, and s.

◆ g()

◆ g_pt()

template<unsigned DIM>
Vector<double>*& oomph::NavierStokesEquations< DIM >::g_pt ( )
inline

Pointer to Vector of gravitational components.

772  {
773  return G_pt;
774  }

References oomph::NavierStokesEquations< DIM >::G_pt.

Referenced by oomph::RefineableNavierStokesEquations< DIM >::further_build().

◆ get_body_force_gradient_nst()

template<unsigned DIM>
virtual void oomph::NavierStokesEquations< DIM >::get_body_force_gradient_nst ( const double time,
const unsigned ipt,
const Vector< double > &  s,
const Vector< double > &  x,
DenseMatrix< double > &  d_body_force_dx 
)
inlineprotectedvirtual

Get gradient of body force term at (Eulerian) position x. This function is virtual to allow overloading in multi-physics problems where the strength of the source function might be determined by another system of equations. Computed via function pointer (if set) or by finite differencing (default)

550  {
551  // hierher: Implement function pointer version
552  /* //If no gradient function has been set, FD it */
553  /* if(Body_force_fct_gradient_pt==0) */
554  {
555  // Reference value
556  Vector<double> body_force(DIM, 0.0);
557  get_body_force_nst(time, ipt, s, x, body_force);
558 
559  // FD it
561  Vector<double> body_force_pls(DIM, 0.0);
562  Vector<double> x_pls(x);
563  for (unsigned i = 0; i < DIM; i++)
564  {
565  x_pls[i] += eps_fd;
566  get_body_force_nst(time, ipt, s, x_pls, body_force_pls);
567  for (unsigned j = 0; j < DIM; j++)
568  {
569  d_body_force_dx(j, i) =
570  (body_force_pls[j] - body_force[j]) / eps_fd;
571  }
572  x_pls[i] = x[i];
573  }
574  }
575  /* else */
576  /* { */
577  /* // Get gradient */
578  /* (*Source_fct_gradient_pt)(time,x,gradient); */
579  /* } */
580  }
static double Default_fd_jacobian_step
Definition: elements.h:1198

References Global_Parameters::body_force(), oomph::GeneralisedElement::Default_fd_jacobian_step, DIM, oomph::NavierStokesEquations< DIM >::get_body_force_nst(), i, j, s, and plotDoE::x.

◆ get_body_force_nst()

template<unsigned DIM>
virtual void oomph::NavierStokesEquations< DIM >::get_body_force_nst ( const double time,
const unsigned ipt,
const Vector< double > &  s,
const Vector< double > &  x,
Vector< double > &  result 
)
inlineprotectedvirtual

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

Reimplemented in oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >, oomph::BuoyantQCrouzeixRaviartElement< DIM >, oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >, QCrouzeixRaviartElementWithTwoExternalElement< DIM >, RefineableQCrouzeixRaviartElementWithTwoExternalElement< DIM >, oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >, RefineableQCrouzeixRaviartElementWithExternalElement< DIM >, and RefineableBuoyantQCrouzeixRaviartElement< DIM >.

522  {
523  // If the function pointer is zero return zero
524  if (Body_force_fct_pt == 0)
525  {
526  // Loop over dimensions and set body forces to zero
527  for (unsigned i = 0; i < DIM; i++)
528  {
529  result[i] = 0.0;
530  }
531  }
532  // Otherwise call the function
533  else
534  {
535  (*Body_force_fct_pt)(time, x, result);
536  }
537  }

References oomph::NavierStokesEquations< DIM >::Body_force_fct_pt, DIM, i, and plotDoE::x.

Referenced by oomph::NavierStokesEquations< DIM >::get_body_force_gradient_nst().

◆ get_dresidual_dnodal_coordinates()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::get_dresidual_dnodal_coordinates ( RankThreeTensor< double > &  dresidual_dnodal_coordinates)
virtual

Compute derivatives of elemental residual vector with respect to nodal coordinates. Overwrites default implementation in 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 FE base class.

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::RefineableNavierStokesEquations< DIM >.

2125  {
2126  // Return immediately if there are no dofs
2127  if (ndof() == 0)
2128  {
2129  return;
2130  }
2131 
2132  // Determine number of nodes in element
2133  const unsigned n_node = nnode();
2134 
2135  // Get continuous time from timestepper of first node
2136  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
2137 
2138  // Determine number of pressure dofs in element
2139  const unsigned n_pres = npres_nst();
2140 
2141  // Find the indices at which the local velocities are stored
2142  unsigned u_nodal_index[DIM];
2143  for (unsigned i = 0; i < DIM; i++)
2144  {
2145  u_nodal_index[i] = u_index_nst(i);
2146  }
2147 
2148  // Set up memory for the shape and test functions
2149  Shape psif(n_node), testf(n_node);
2150  DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
2151 
2152  // Set up memory for pressure shape and test functions
2153  Shape psip(n_pres), testp(n_pres);
2154 
2155  // Deriatives of shape fct derivatives w.r.t. nodal coords
2156  RankFourTensor<double> d_dpsifdx_dX(DIM, n_node, n_node, DIM);
2157  RankFourTensor<double> d_dtestfdx_dX(DIM, n_node, n_node, DIM);
2158 
2159  // Derivative of Jacobian of mapping w.r.t. to nodal coords
2160  DenseMatrix<double> dJ_dX(DIM, n_node);
2161 
2162  // Derivatives of derivative of u w.r.t. nodal coords
2163  RankFourTensor<double> d_dudx_dX(DIM, n_node, DIM, DIM);
2164 
2165  // Derivatives of nodal velocities w.r.t. nodal coords:
2166  // Assumption: Interaction only local via no-slip so
2167  // X_ij only affects U_ij.
2168  DenseMatrix<double> d_U_dX(DIM, n_node, 0.0);
2169 
2170  // Determine the number of integration points
2171  const unsigned n_intpt = integral_pt()->nweight();
2172 
2173  // Vector to hold local coordinates
2174  Vector<double> s(DIM);
2175 
2176  // Get physical variables from element
2177  // (Reynolds number must be multiplied by the density ratio)
2178  double scaled_re = re() * density_ratio();
2179  double scaled_re_st = re_st() * density_ratio();
2180  double scaled_re_inv_fr = re_invfr() * density_ratio();
2181  double visc_ratio = viscosity_ratio();
2182  Vector<double> G = g();
2183 
2184  // FD step
2186 
2187  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
2188  // Assumption: Interaction only local via no-slip so
2189  // X_ij only affects U_ij.
2190  bool element_has_node_with_aux_node_update_fct = false;
2191  for (unsigned q = 0; q < n_node; q++)
2192  {
2193  Node* nod_pt = node_pt(q);
2194 
2195  // Only compute if there's a node-update fct involved
2196  if (nod_pt->has_auxiliary_node_update_fct_pt())
2197  {
2198  element_has_node_with_aux_node_update_fct = true;
2199 
2200  // Current nodal velocity
2201  Vector<double> u_ref(DIM);
2202  for (unsigned i = 0; i < DIM; i++)
2203  {
2204  u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
2205  }
2206 
2207  // FD
2208  for (unsigned p = 0; p < DIM; p++)
2209  {
2210  // Make backup
2211  double backup = nod_pt->x(p);
2212 
2213  // Do FD step. No node update required as we're
2214  // attacking the coordinate directly...
2215  nod_pt->x(p) += eps_fd;
2216 
2217  // Do auxiliary node update (to apply no slip)
2218  nod_pt->perform_auxiliary_node_update_fct();
2219 
2220  // Evaluate
2221  d_U_dX(p, q) =
2222  (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
2223 
2224  // Reset
2225  nod_pt->x(p) = backup;
2226 
2227  // Do auxiliary node update (to apply no slip)
2228  nod_pt->perform_auxiliary_node_update_fct();
2229  }
2230  }
2231  }
2232 
2233  // Integer to store the local equation number
2234  int local_eqn = 0;
2235 
2236  // Loop over the integration points
2237  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2238  {
2239  // Assign values of s
2240  for (unsigned i = 0; i < DIM; i++)
2241  {
2242  s[i] = integral_pt()->knot(ipt, i);
2243  }
2244 
2245  // Get the integral weight
2246  const double w = integral_pt()->weight(ipt);
2247 
2248  // Call the derivatives of the shape and test functions
2249  const double J = dshape_and_dtest_eulerian_at_knot_nst(ipt,
2250  psif,
2251  dpsifdx,
2252  d_dpsifdx_dX,
2253  testf,
2254  dtestfdx,
2255  d_dtestfdx_dX,
2256  dJ_dX);
2257 
2258  // Call the pressure shape and test functions
2259  pshape_nst(s, psip, testp);
2260 
2261  // Calculate local values of the pressure and velocity components
2262  // Allocate
2263  double interpolated_p = 0.0;
2264  Vector<double> interpolated_u(DIM, 0.0);
2265  Vector<double> interpolated_x(DIM, 0.0);
2266  Vector<double> mesh_velocity(DIM, 0.0);
2267  Vector<double> dudt(DIM, 0.0);
2268  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
2269 
2270  // Calculate pressure
2271  for (unsigned l = 0; l < n_pres; l++)
2272  {
2273  interpolated_p += p_nst(l) * psip[l];
2274  }
2275 
2276  // Calculate velocities and derivatives:
2277 
2278  // Loop over nodes
2279  for (unsigned l = 0; l < n_node; l++)
2280  {
2281  // Loop over directions
2282  for (unsigned i = 0; i < DIM; i++)
2283  {
2284  // Get the nodal value
2285  double u_value = raw_nodal_value(l, u_nodal_index[i]);
2286  interpolated_u[i] += u_value * psif[l];
2287  interpolated_x[i] += raw_nodal_position(l, i) * psif[l];
2288  dudt[i] += du_dt_nst(l, i) * psif[l];
2289 
2290  // Loop over derivative directions
2291  for (unsigned j = 0; j < DIM; j++)
2292  {
2293  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
2294  }
2295  }
2296  }
2297 
2298  if (!ALE_is_disabled)
2299  {
2300  // Loop over nodes
2301  for (unsigned l = 0; l < n_node; l++)
2302  {
2303  // Loop over directions
2304  for (unsigned i = 0; i < DIM; i++)
2305  {
2306  mesh_velocity[i] += this->raw_dnodal_position_dt(l, i) * psif[l];
2307  }
2308  }
2309  }
2310 
2311  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
2312  for (unsigned q = 0; q < n_node; q++)
2313  {
2314  // Loop over coordinate directions
2315  for (unsigned p = 0; p < DIM; p++)
2316  {
2317  for (unsigned i = 0; i < DIM; i++)
2318  {
2319  for (unsigned k = 0; k < DIM; k++)
2320  {
2321  double aux = 0.0;
2322  for (unsigned j = 0; j < n_node; j++)
2323  {
2324  aux += raw_nodal_value(j, u_nodal_index[i]) *
2325  d_dpsifdx_dX(p, q, j, k);
2326  }
2327  d_dudx_dX(p, q, i, k) = aux;
2328  }
2329  }
2330  }
2331  }
2332 
2333  // Get weight of actual nodal position/value in computation of mesh
2334  // velocity from positional/value time stepper
2335  const double pos_time_weight =
2337  const double val_time_weight =
2338  node_pt(0)->time_stepper_pt()->weight(1, 0);
2339 
2340  // Get the user-defined body force terms
2341  Vector<double> body_force(DIM);
2343 
2344  // Get the user-defined source function
2345  const double source = get_source_nst(time, ipt, interpolated_x);
2346 
2347  // Note: Can use raw values and avoid bypassing hanging information
2348  // because this is the non-refineable version!
2349 
2350  // Get gradient of body force function
2351  DenseMatrix<double> d_body_force_dx(DIM, DIM, 0.0);
2353  time, ipt, s, interpolated_x, d_body_force_dx);
2354 
2355  // Get gradient of source function
2356  Vector<double> source_gradient(DIM, 0.0);
2357  get_source_gradient_nst(time, ipt, interpolated_x, source_gradient);
2358 
2359 
2360  // Assemble shape derivatives
2361  // --------------------------
2362 
2363  // MOMENTUM EQUATIONS
2364  // ------------------
2365 
2366  // Loop over the test functions
2367  for (unsigned l = 0; l < n_node; l++)
2368  {
2369  // Loop over coordinate directions
2370  for (unsigned i = 0; i < DIM; i++)
2371  {
2372  // Get the local equation
2373  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
2374  ;
2375 
2376  // IF it's not a boundary condition
2377  if (local_eqn >= 0)
2378  {
2379  // Loop over coordinate directions
2380  for (unsigned p = 0; p < DIM; p++)
2381  {
2382  // Loop over nodes
2383  for (unsigned q = 0; q < n_node; q++)
2384  {
2385  // Residual x deriv of Jacobian
2386  //-----------------------------
2387 
2388  // Add the user-defined body force terms
2389  double sum = body_force[i] * testf[l];
2390 
2391  // Add the gravitational body force term
2392  sum += scaled_re_inv_fr * testf[l] * G[i];
2393 
2394  // Add the pressure gradient term
2395  sum += interpolated_p * dtestfdx(l, i);
2396 
2397  // Add in the stress tensor terms
2398  // The viscosity ratio needs to go in here to ensure
2399  // continuity of normal stress is satisfied even in flows
2400  // with zero pressure gradient!
2401  for (unsigned k = 0; k < DIM; k++)
2402  {
2403  sum -= visc_ratio *
2404  (interpolated_dudx(i, k) +
2405  Gamma[i] * interpolated_dudx(k, i)) *
2406  dtestfdx(l, k);
2407  }
2408 
2409  // Add in the inertial terms
2410 
2411  // du/dt term
2412  sum -= scaled_re_st * dudt[i] * testf[l];
2413 
2414  // Convective terms, including mesh velocity
2415  for (unsigned k = 0; k < DIM; k++)
2416  {
2417  double tmp = scaled_re * interpolated_u[k];
2418  if (!ALE_is_disabled)
2419  {
2420  tmp -= scaled_re_st * mesh_velocity[k];
2421  }
2422  sum -= tmp * interpolated_dudx(i, k) * testf[l];
2423  }
2424 
2425  // Multiply through by deriv of Jacobian and integration weight
2426  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2427  sum * dJ_dX(p, q) * w;
2428 
2429  // Derivative of residual x Jacobian
2430  //----------------------------------
2431 
2432  // Body force
2433  sum = d_body_force_dx(i, p) * psif(q) * testf(l);
2434 
2435  // Pressure gradient term
2436  sum += interpolated_p * d_dtestfdx_dX(p, q, l, i);
2437 
2438  // Viscous term
2439  for (unsigned k = 0; k < DIM; k++)
2440  {
2441  sum -= visc_ratio * ((interpolated_dudx(i, k) +
2442  Gamma[i] * interpolated_dudx(k, i)) *
2443  d_dtestfdx_dX(p, q, l, k) +
2444  (d_dudx_dX(p, q, i, k) +
2445  Gamma[i] * d_dudx_dX(p, q, k, i)) *
2446  dtestfdx(l, k));
2447  }
2448 
2449  // Convective terms, including mesh velocity
2450  for (unsigned k = 0; k < DIM; k++)
2451  {
2452  double tmp = scaled_re * interpolated_u[k];
2453  if (!ALE_is_disabled)
2454  {
2455  tmp -= scaled_re_st * mesh_velocity[k];
2456  }
2457  sum -= tmp * d_dudx_dX(p, q, i, k) * testf(l);
2458  }
2459  if (!ALE_is_disabled)
2460  {
2461  sum += scaled_re_st * pos_time_weight * psif(q) *
2462  interpolated_dudx(i, p) * testf(l);
2463  }
2464 
2465  // Multiply through by Jacobian and integration weight
2466  dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2467 
2468  // Derivs w.r.t. to nodal velocities
2469  // ---------------------------------
2470  if (element_has_node_with_aux_node_update_fct)
2471  {
2472  sum =
2473  -visc_ratio * Gamma[i] * dpsifdx(q, i) * dtestfdx(l, p) -
2474  scaled_re * psif(q) * interpolated_dudx(i, p) * testf(l);
2475  if (i == p)
2476  {
2477  sum -= scaled_re_st * val_time_weight * psif(q) * testf(l);
2478  for (unsigned k = 0; k < DIM; k++)
2479  {
2480  sum -= visc_ratio * dpsifdx(q, k) * dtestfdx(l, k);
2481  double tmp = scaled_re * interpolated_u[k];
2482  if (!ALE_is_disabled)
2483  tmp -= scaled_re_st * mesh_velocity[k];
2484  sum -= tmp * dpsifdx(q, k) * testf(l);
2485  }
2486  }
2487  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2488  sum * d_U_dX(p, q) * J * w;
2489  }
2490  }
2491  }
2492  }
2493  }
2494  } // End of loop over test functions
2495 
2496 
2497  // CONTINUITY EQUATION
2498  // -------------------
2499 
2500  // Loop over the shape functions
2501  for (unsigned l = 0; l < n_pres; l++)
2502  {
2503  local_eqn = p_local_eqn(l);
2504 
2505  // If not a boundary conditions
2506  if (local_eqn >= 0)
2507  {
2508  // Loop over coordinate directions
2509  for (unsigned p = 0; p < DIM; p++)
2510  {
2511  // Loop over nodes
2512  for (unsigned q = 0; q < n_node; q++)
2513  {
2514  // Residual x deriv of Jacobian
2515  //-----------------------------
2516 
2517  // Source term
2518  double aux = -source;
2519 
2520  // Loop over velocity components
2521  for (unsigned k = 0; k < DIM; k++)
2522  {
2523  aux += interpolated_dudx(k, k);
2524  }
2525 
2526  // Multiply through by deriv of Jacobian and integration weight
2527  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2528  aux * dJ_dX(p, q) * testp[l] * w;
2529 
2530  // Derivative of residual x Jacobian
2531  //----------------------------------
2532 
2533  // Loop over velocity components
2534  aux = -source_gradient[p] * psif(q);
2535  for (unsigned k = 0; k < DIM; k++)
2536  {
2537  aux += d_dudx_dX(p, q, k, k);
2538  }
2539  // Multiply through by Jacobian and integration weight
2540  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2541  aux * testp[l] * J * w;
2542 
2543  // Derivs w.r.t. to nodal velocities
2544  //---------------------------------
2545  if (element_has_node_with_aux_node_update_fct)
2546  {
2547  aux = dpsifdx(q, p) * testp(l);
2548  dresidual_dnodal_coordinates(local_eqn, p, q) +=
2549  aux * d_U_dX(p, q) * J * w;
2550  }
2551  }
2552  }
2553  }
2554  }
2555  } // End of loop over integration points
2556  }
float * p
Definition: Tutorial_Map_using.cpp:9
virtual void get_source_gradient_nst(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
Definition: navier_stokes_elements.h:607
virtual void get_body_force_gradient_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
Definition: navier_stokes_elements.h:544
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019

References oomph::ALE_is_disabled, Global_Parameters::body_force(), DIM, G, GlobalPhysicalVariables::Gamma, oomph::Node::has_auxiliary_node_update_fct_pt(), i, J, j, k, p, oomph::Node::perform_auxiliary_node_update_fct(), Eigen::numext::q, s, TestProblem::source(), tmp, oomph::Data::value_pt(), w, and oomph::Node::x().

◆ get_load()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::get_load ( const Vector< double > &  s,
const Vector< double > &  N,
Vector< double > &  load 
)
inlinevirtual

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

Implements oomph::FSIFluidElement.

995  {
996  // Note: get_traction() computes the traction onto the fluid
997  // if N is the outer unit normal onto the fluid; here we're
998  // exepcting N to point into the fluid so we're getting the
999  // traction onto the adjacent wall instead!
1000  get_traction(s, N, load);
1001  }
void load(Archive &ar, ParticleHandler &handl)
Definition: Particles.h:21
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
Definition: navier_stokes_elements.cc:1098
@ N
Definition: constructor.cpp:22

References oomph::NavierStokesEquations< DIM >::get_traction(), load(), N, and s.

◆ get_pressure_and_velocity_mass_matrix_diagonal()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::get_pressure_and_velocity_mass_matrix_diagonal ( Vector< double > &  press_mass_diag,
Vector< double > &  veloc_mass_diag,
const unsigned which_one = 0 
)
virtual

Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed, otherwise only the pressure (which_one=1) or the velocity mass matrix (which_one=2 – the LSC version of the preconditioner only needs that one)

Implements oomph::TemplateFreeNavierStokesEquationsBase.

Reimplemented in oomph::RefineableNavierStokesEquations< DIM >.

75  {
76  // Resize and initialise
77  unsigned n_dof = ndof();
78 
79  if ((which_one == 0) || (which_one == 1))
80  press_mass_diag.assign(n_dof, 0.0);
81  if ((which_one == 0) || (which_one == 2))
82  veloc_mass_diag.assign(n_dof, 0.0);
83 
84  // find out how many nodes there are
85  unsigned n_node = nnode();
86 
87  // find number of spatial dimensions
88  unsigned n_dim = this->dim();
89 
90  // Local coordinates
91  Vector<double> s(n_dim);
92 
93  // find the indices at which the local velocities are stored
94  Vector<unsigned> u_nodal_index(n_dim);
95  for (unsigned i = 0; i < n_dim; i++)
96  {
97  u_nodal_index[i] = this->u_index_nst(i);
98  }
99 
100  // Set up memory for veloc shape functions
101  Shape psi(n_node);
102 
103  // Find number of pressure dofs
104  unsigned n_pres = npres_nst();
105 
106  // Pressure shape function
107  Shape psi_p(n_pres);
108 
109  // Number of integration points
110  unsigned n_intpt = integral_pt()->nweight();
111 
112  // Integer to store the local equations no
113  int local_eqn = 0;
114 
115  // Loop over the integration points
116  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
117  {
118  // Get the integral weight
119  double w = integral_pt()->weight(ipt);
120 
121  // Get determinant of Jacobian of the mapping
122  double J = J_eulerian_at_knot(ipt);
123 
124  // Assign values of s
125  for (unsigned i = 0; i < n_dim; i++)
126  {
127  s[i] = integral_pt()->knot(ipt, i);
128  }
129 
130  // Premultiply weights and Jacobian
131  double W = w * J;
132 
133 
134  // Do we want the velocity one?
135  if ((which_one == 0) || (which_one == 2))
136  {
137  // Get the velocity shape functions
138  shape_at_knot(ipt, psi);
139 
140  // Loop over the veclocity shape functions
141  for (unsigned l = 0; l < n_node; l++)
142  {
143  // Loop over the velocity components
144  for (unsigned i = 0; i < n_dim; i++)
145  {
146  local_eqn = nodal_local_eqn(l, u_nodal_index[i]);
147 
148  // If not a boundary condition
149  if (local_eqn >= 0)
150  {
151  // Add the contribution
152  veloc_mass_diag[local_eqn] += pow(psi[l], 2) * W;
153  }
154  }
155  }
156  }
157 
158  // Do we want the pressure one?
159  if ((which_one == 0) || (which_one == 1))
160  {
161  // Get the pressure shape functions
162  pshape_nst(s, psi_p);
163 
164  // Loop over the veclocity shape functions
165  for (unsigned l = 0; l < n_pres; l++)
166  {
167  // Get equation number
168  local_eqn = p_local_eqn(l);
169 
170  // If not a boundary condition
171  if (local_eqn >= 0)
172  {
173  // Add the contribution
174  press_mass_diag[local_eqn] += pow(psi_p[l], 2) * W;
175  }
176  }
177  }
178  }
179  }
unsigned dim() const
Definition: elements.h:2611
virtual double J_eulerian_at_knot(const unsigned &ipt) const
Definition: elements.cc:4168
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Definition: elements.cc:3220

References i, J, Eigen::bfloat16_impl::pow(), s, w, and oomph::QuadTreeNames::W.

◆ get_source_gradient_nst()

template<unsigned DIM>
virtual void oomph::NavierStokesEquations< DIM >::get_source_gradient_nst ( const double time,
const unsigned ipt,
const Vector< double > &  x,
Vector< double > &  gradient 
)
inlineprotectedvirtual

Get gradient of source term at (Eulerian) position x. This function is virtual to allow overloading in multi-physics problems where the strength of the source function might be determined by another system of equations. Computed via function pointer (if set) or by finite differencing (default)

611  {
612  // hierher: Implement function pointer version
613  /* //If no gradient function has been set, FD it */
614  /* if(Source_fct_gradient_pt==0) */
615  {
616  // Reference value
617  double source = get_source_nst(time, ipt, x);
618 
619  // FD it
621  double source_pls = 0.0;
622  Vector<double> x_pls(x);
623  for (unsigned i = 0; i < DIM; i++)
624  {
625  x_pls[i] += eps_fd;
626  source_pls = get_source_nst(time, ipt, x_pls);
627  gradient[i] = (source_pls - source) / eps_fd;
628  x_pls[i] = x[i];
629  }
630  }
631  /* else */
632  /* { */
633  /* // Get gradient */
634  /* (*Source_fct_gradient_pt)(time,x,gradient); */
635  /* } */
636  }

References oomph::GeneralisedElement::Default_fd_jacobian_step, DIM, oomph::NavierStokesEquations< DIM >::get_source_nst(), i, TestProblem::source(), and plotDoE::x.

◆ get_source_nst()

template<unsigned DIM>
virtual double oomph::NavierStokesEquations< DIM >::get_source_nst ( const double time,
const unsigned ipt,
const Vector< double > &  x 
)
inlineprotectedvirtual

Calculate the source fct at given time and Eulerian position

588  {
589  // If the function pointer is zero return zero
590  if (Source_fct_pt == 0)
591  {
592  return 0;
593  }
594  // Otherwise call the function
595  else
596  {
597  return (*Source_fct_pt)(time, x);
598  }
599  }

References oomph::NavierStokesEquations< DIM >::Source_fct_pt, and plotDoE::x.

Referenced by oomph::NavierStokesEquations< DIM >::get_source_gradient_nst().

◆ get_traction() [1/2]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::get_traction ( const Vector< double > &  s,
const Vector< double > &  N,
Vector< double > &  traction 
)

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

1101  {
1102  // Get velocity gradients
1103  DenseMatrix<double> strainrate(DIM, DIM);
1104  strain_rate(s, strainrate);
1105 
1106  // Get pressure
1107  double press = interpolated_p_nst(s);
1108 
1109  // Loop over traction components
1110  for (unsigned i = 0; i < DIM; i++)
1111  {
1112  traction[i] = -press * N[i];
1113  for (unsigned j = 0; j < DIM; j++)
1114  {
1115  traction[i] += 2.0 * strainrate(i, j) * N[j];
1116  }
1117  }
1118  }

References DIM, i, j, N, and s.

Referenced by oomph::NavierStokesEquations< DIM >::get_load().

◆ get_traction() [2/2]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::get_traction ( const Vector< double > &  s,
const Vector< double > &  N,
Vector< double > &  traction_p,
Vector< double > &  traction_visc_n,
Vector< double > &  traction_visc_t 
)

Compute traction (on the viscous scale) exerted onto the fluid at local coordinate s, decomposed into pressure and normal and tangential viscous components. N has to be outer unit normal to the fluid.

1132  {
1133  Vector<double> traction_visc(DIM);
1134 
1135  // Get velocity gradients
1136  DenseMatrix<double> strainrate(DIM, DIM);
1137  strain_rate(s, strainrate);
1138 
1139  // Get pressure
1140  double press = interpolated_p_nst(s);
1141 
1142  // Loop over traction components
1143  for (unsigned i = 0; i < DIM; i++)
1144  {
1145  // pressure
1146  traction_p[i] = -press * N[i];
1147  for (unsigned j = 0; j < DIM; j++)
1148  {
1149  // viscous stress
1150  traction_visc[i] += 2.0 * strainrate(i, j) * N[j];
1151  }
1152  // Decompose viscous stress into normal and tangential components
1153  traction_visc_n[i] = traction_visc[i] * N[i];
1154  traction_visc_t[i] = traction_visc[i] - traction_visc_n[i];
1155  }
1156  }

References DIM, i, j, N, and s.

◆ get_vorticity() [1/5]

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

Compute vorticity vector at local coordinate s and return the one and only component of vorticity vector (only makes sense when solving the 2D N.St. equations)

1327  {
1328  // Create a vector to store the vorticity
1329  Vector<double> vort(1, 0.0);
1330 
1331  // Calculate the vorticity
1332  get_vorticity(s, vort);
1333 
1334  // Assign the vorticity
1335  vorticity = vort[0];
1336  }
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.

References s.

◆ get_vorticity() [2/5]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::get_vorticity ( const Vector< double > &  s,
double vorticity 
) const

Compute the scalar vorticity at local coordinate s (2D)

◆ get_vorticity() [3/5]

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

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

1258  {
1259 #ifdef PARANOID
1260  if (vorticity.size() != 1)
1261  {
1262  std::ostringstream error_message;
1263  error_message << "Computation of vorticity in 2D requires a 1D vector\n"
1264  << "which contains the only non-zero component of the\n"
1265  << "vorticity vector. You've passed a vector of size "
1266  << vorticity.size() << std::endl;
1267 
1268  throw OomphLibError(
1269  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1270  }
1271 #endif
1272 
1273  // Specify spatial dimension
1274  unsigned DIM = 2;
1275 
1276  // Velocity gradient matrix
1277  DenseMatrix<double> dudx(DIM);
1278 
1279  // Find out how many nodes there are in the element
1280  unsigned n_node = nnode();
1281 
1282  // Set up memory for the shape and test functions
1283  Shape psi(n_node);
1284  DShape dpsidx(n_node, DIM);
1285 
1286  // Call the derivatives of the shape functions
1287  dshape_eulerian(s, psi, dpsidx);
1288 
1289  // Initialise to zero
1290  for (unsigned i = 0; i < DIM; i++)
1291  {
1292  for (unsigned j = 0; j < DIM; j++)
1293  {
1294  dudx(i, j) = 0.0;
1295  }
1296  }
1297 
1298  // Loop over veclocity components
1299  for (unsigned i = 0; i < DIM; i++)
1300  {
1301  // Get the index at which the i-th velocity is stored
1302  unsigned u_nodal_index = u_index_nst(i);
1303  // Loop over derivative directions
1304  for (unsigned j = 0; j < DIM; j++)
1305  {
1306  // Loop over nodes
1307  for (unsigned l = 0; l < n_node; l++)
1308  {
1309  dudx(i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1310  }
1311  }
1312  }
1313 
1314  // Z-component of vorticity
1315  vorticity[0] = dudx(1, 0) - dudx(0, 1);
1316  }

References DIM, i, j, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

◆ get_vorticity() [4/5]

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

Compute 3D vorticity vector at local coordinate s.

1345  {
1346 #ifdef PARANOID
1347  if (vorticity.size() != 3)
1348  {
1349  std::ostringstream error_message;
1350  error_message << "Computation of vorticity in 3D requires a 3D vector\n"
1351  << "which contains the only non-zero component of the\n"
1352  << "vorticity vector. You've passed a vector of size "
1353  << vorticity.size() << std::endl;
1354 
1355  throw OomphLibError(
1356  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1357  }
1358 #endif
1359 
1360  // Specify spatial dimension
1361  unsigned DIM = 3;
1362 
1363  // Velocity gradient matrix
1364  DenseMatrix<double> dudx(DIM);
1365 
1366  // Find out how many nodes there are in the element
1367  unsigned n_node = nnode();
1368 
1369  // Set up memory for the shape and test functions
1370  Shape psi(n_node);
1371  DShape dpsidx(n_node, DIM);
1372 
1373  // Call the derivatives of the shape functions
1374  dshape_eulerian(s, psi, dpsidx);
1375 
1376  // Initialise to zero
1377  for (unsigned i = 0; i < DIM; i++)
1378  {
1379  for (unsigned j = 0; j < DIM; j++)
1380  {
1381  dudx(i, j) = 0.0;
1382  }
1383  }
1384 
1385  // Loop over veclocity components
1386  for (unsigned i = 0; i < DIM; i++)
1387  {
1388  // Get the index at which the i-th velocity component is stored
1389  unsigned u_nodal_index = u_index_nst(i);
1390  // Loop over derivative directions
1391  for (unsigned j = 0; j < DIM; j++)
1392  {
1393  // Loop over nodes
1394  for (unsigned l = 0; l < n_node; l++)
1395  {
1396  dudx(i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1397  }
1398  }
1399  }
1400 
1401  vorticity[0] = dudx(2, 1) - dudx(1, 2);
1402  vorticity[1] = dudx(0, 2) - dudx(2, 0);
1403  vorticity[2] = dudx(1, 0) - dudx(0, 1);
1404  }

References DIM, i, j, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

◆ get_vorticity() [5/5]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::get_vorticity ( const Vector< double > &  s,
Vector< double > &  vorticity 
) const

Compute the vorticity vector at local coordinate s.

◆ interpolated_dudx_nst()

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::interpolated_dudx_nst ( const Vector< double > &  s,
const unsigned i,
const unsigned j 
) const
inline

Return FE interpolated derivatives of velocity component u[i] w.r.t spatial global coordinate direction x[j] at local coordinate s

1687  {
1688  // Determine number of nodes
1689  const unsigned n_node = nnode();
1690 
1691  // Allocate storage for local shape function and its derivatives
1692  // with respect to space
1693  Shape psif(n_node);
1694  DShape dpsifdx(n_node, DIM);
1695 
1696  // Find values of shape function (ignore jacobian)
1697  (void)this->dshape_eulerian(s, psif, dpsifdx);
1698 
1699  // Get the index at which the velocity is stored
1700  const unsigned u_nodal_index = u_index_nst(i);
1701 
1702  // Initialise value of dudx
1703  double interpolated_dudx = 0.0;
1704 
1705  // Loop over the local nodes and sum
1706  for (unsigned l = 0; l < n_node; l++)
1707  {
1708  interpolated_dudx += nodal_value(l, u_nodal_index) * dpsifdx(l, j);
1709  }
1710 
1711  return (interpolated_dudx);
1712  }

References DIM, oomph::FiniteElement::dshape_eulerian(), i, j, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), and oomph::NavierStokesEquations< DIM >::u_index_nst().

Referenced by LinearisedQTaylorHoodMultiDomainElement::get_base_flow_dudx(), LinearisedQCrouzeixRaviartMultiDomainElement::get_base_flow_dudx(), RefineableLinearisedQTaylorHoodMultiDomainElement::get_base_flow_dudx(), and RefineableLinearisedQCrouzeixRaviartMultiDomainElement::get_base_flow_dudx().

◆ interpolated_p_nst() [1/2]

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::interpolated_p_nst ( const unsigned t,
const Vector< double > &  s 
) const
inline

Return FE interpolated pressure at local coordinate s at time level t.

1662  {
1663  // Find number of nodes
1664  unsigned n_pres = npres_nst();
1665  // Local shape function
1666  Shape psi(n_pres);
1667  // Find values of shape function
1668  pshape_nst(s, psi);
1669 
1670  // Initialise value of p
1671  double interpolated_p = 0.0;
1672  // Loop over the local nodes and sum
1673  for (unsigned l = 0; l < n_pres; l++)
1674  {
1675  interpolated_p += p_nst(t, l) * psi[l];
1676  }
1677 
1678  return (interpolated_p);
1679  }

References oomph::NavierStokesEquations< DIM >::npres_nst(), oomph::NavierStokesEquations< DIM >::p_nst(), oomph::NavierStokesEquations< DIM >::pshape_nst(), s, and plotPSD::t.

◆ interpolated_p_nst() [2/2]

template<unsigned DIM>
virtual double oomph::NavierStokesEquations< DIM >::interpolated_p_nst ( const Vector< double > &  s) const
inlinevirtual

Return FE interpolated pressure at local coordinate s.

1640  {
1641  // Find number of nodes
1642  unsigned n_pres = npres_nst();
1643  // Local shape function
1644  Shape psi(n_pres);
1645  // Find values of shape function
1646  pshape_nst(s, psi);
1647 
1648  // Initialise value of p
1649  double interpolated_p = 0.0;
1650  // Loop over the local nodes and sum
1651  for (unsigned l = 0; l < n_pres; l++)
1652  {
1653  interpolated_p += p_nst(l) * psi[l];
1654  }
1655 
1656  return (interpolated_p);
1657  }

References oomph::NavierStokesEquations< DIM >::npres_nst(), oomph::NavierStokesEquations< DIM >::p_nst(), oomph::NavierStokesEquations< DIM >::pshape_nst(), and s.

Referenced by oomph::PRefineableQCrouzeixRaviartElement< DIM >::further_build(), oomph::RefineableQCrouzeixRaviartElement< DIM >::further_build(), oomph::RefineableQTaylorHoodElement< DIM >::get_interpolated_values(), oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >::output(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::output(), oomph::MyTaylorHoodElement< DIM >::output(), oomph::BuoyantQCrouzeixRaviartElement< DIM >::output(), oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::output(), oomph::NavierStokesEquations< DIM >::point_output_data(), and oomph::NavierStokesEquations< DIM >::scalar_value_paraview().

◆ interpolated_u_nst() [1/3]

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::interpolated_u_nst ( const unsigned t,
const Vector< double > &  s,
const unsigned i 
) const
inline

Return FE interpolated velocity u[i] at local coordinate s at time level t (t=0: present; t>0: history)

1558  {
1559  // Find number of nodes
1560  unsigned n_node = nnode();
1561 
1562  // Local shape function
1563  Shape psi(n_node);
1564 
1565  // Find values of shape function
1566  shape(s, psi);
1567 
1568  // Get nodal index at which i-th velocity is stored
1569  unsigned u_nodal_index = u_index_nst(i);
1570 
1571  // Initialise value of u
1572  double interpolated_u = 0.0;
1573  // Loop over the local nodes and sum
1574  for (unsigned l = 0; l < n_node; l++)
1575  {
1576  interpolated_u += nodal_value(t, l, u_nodal_index) * psi[l];
1577  }
1578 
1579  return (interpolated_u);
1580  }

References i, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), s, oomph::FiniteElement::shape(), plotPSD::t, and oomph::NavierStokesEquations< DIM >::u_index_nst().

◆ interpolated_u_nst() [2/3]

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::interpolated_u_nst ( const Vector< double > &  s,
const unsigned i 
) const
inline

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

1531  {
1532  // Find number of nodes
1533  unsigned n_node = nnode();
1534  // Local shape function
1535  Shape psi(n_node);
1536  // Find values of shape function
1537  shape(s, psi);
1538 
1539  // Get nodal index at which i-th velocity is stored
1540  unsigned u_nodal_index = u_index_nst(i);
1541 
1542  // Initialise value of u
1543  double interpolated_u = 0.0;
1544  // Loop over the local nodes and sum
1545  for (unsigned l = 0; l < n_node; l++)
1546  {
1547  interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
1548  }
1549 
1550  return (interpolated_u);
1551  }

References i, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), s, oomph::FiniteElement::shape(), and oomph::NavierStokesEquations< DIM >::u_index_nst().

◆ interpolated_u_nst() [3/3]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::interpolated_u_nst ( const Vector< double > &  s,
Vector< double > &  veloc 
) const
inline

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

1507  {
1508  // Find number of nodes
1509  unsigned n_node = nnode();
1510  // Local shape function
1511  Shape psi(n_node);
1512  // Find values of shape function
1513  shape(s, psi);
1514 
1515  for (unsigned i = 0; i < DIM; i++)
1516  {
1517  // Index at which the nodal value is stored
1518  unsigned u_nodal_index = u_index_nst(i);
1519  // Initialise value of u
1520  veloc[i] = 0.0;
1521  // Loop over the local nodes and sum
1522  for (unsigned l = 0; l < n_node; l++)
1523  {
1524  veloc[i] += nodal_value(l, u_nodal_index) * psi[l];
1525  }
1526  }
1527  }

References DIM, i, oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), s, oomph::FiniteElement::shape(), and oomph::NavierStokesEquations< DIM >::u_index_nst().

Referenced by LinearisedQTaylorHoodMultiDomainElement::get_base_flow_u(), LinearisedQCrouzeixRaviartMultiDomainElement::get_base_flow_u(), RefineableLinearisedQTaylorHoodMultiDomainElement::get_base_flow_u(), RefineableLinearisedQCrouzeixRaviartMultiDomainElement::get_base_flow_u(), oomph::RefineableQTaylorHoodElement< DIM >::get_interpolated_values(), oomph::RefineableQCrouzeixRaviartElement< DIM >::get_interpolated_values(), oomph::PRefineableQCrouzeixRaviartElement< DIM >::get_interpolated_values(), RefineableQAdvectionDiffusionElementWithExternalElement< DIM >::get_wind_adv_diff(), QAdvectionDiffusionElementWithExternalElement< DIM >::get_wind_adv_diff(), oomph::BuoyantQCrouzeixRaviartElement< DIM >::get_wind_adv_diff(), oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::get_wind_adv_diff(), oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >::get_wind_adv_diff_react(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::get_wind_adv_diff_react(), oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >::output(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::output(), oomph::MyTaylorHoodElement< DIM >::output(), oomph::BuoyantQCrouzeixRaviartElement< DIM >::output(), oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::output(), oomph::NavierStokesEquations< DIM >::point_output_data(), and oomph::NavierStokesEquations< DIM >::scalar_value_paraview().

◆ kin_energy()

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::kin_energy

Get integral of kinetic energy over element.

Get integral of kinetic energy over element: Note that this is the "raw" kinetic energy in the sense that the density ratio has not been included. In problems with two or more fluids the user will have to remember to premultiply certain elements by the appropriate density ratio.

1417  {
1418  // Initialise
1419  double kin_en = 0.0;
1420 
1421  // Set the value of n_intpt
1422  unsigned n_intpt = integral_pt()->nweight();
1423 
1424  // Set the Vector to hold local coordinates
1425  Vector<double> s(DIM);
1426 
1427  // Loop over the integration points
1428  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1429  {
1430  // Assign values of s
1431  for (unsigned i = 0; i < DIM; i++)
1432  {
1433  s[i] = integral_pt()->knot(ipt, i);
1434  }
1435 
1436  // Get the integral weight
1437  double w = integral_pt()->weight(ipt);
1438 
1439  // Get Jacobian of mapping
1440  double J = J_eulerian(s);
1441 
1442  // Loop over directions
1443  double veloc_squared = 0.0;
1444  for (unsigned i = 0; i < DIM; i++)
1445  {
1446  veloc_squared += interpolated_u_nst(s, i) * interpolated_u_nst(s, i);
1447  }
1448 
1449  kin_en += 0.5 * veloc_squared * w * J;
1450  }
1451 
1452  return kin_en;
1453  }

References DIM, i, J, s, and w.

◆ n_u_nst()

template<unsigned DIM>
unsigned oomph::NavierStokesEquations< DIM >::n_u_nst ( ) const
inline

Return the number of velocity components Used in FluidInterfaceElements

874  {
875  return DIM;
876  }

References DIM.

◆ npres_nst()

◆ nscalar_paraview()

template<unsigned DIM>
unsigned oomph::NavierStokesEquations< DIM >::nscalar_paraview ( ) const
inlinevirtual

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

Reimplemented from oomph::FiniteElement.

1015  {
1016  return DIM + 1;
1017  }

References DIM.

◆ output() [1/4]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::output ( FILE *  file_pt)
inlinevirtual

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

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::TTaylorHoodElement< DIM >, oomph::TCrouzeixRaviartElement< DIM >, oomph::QTaylorHoodElement< DIM >, oomph::QTaylorHoodElement< 2 >, oomph::QCrouzeixRaviartElement< DIM >, oomph::PseudoSolidNodeUpdateElement< TCrouzeixRaviartElement< 2 >, TPVDBubbleEnrichedElement< 2, 3 > >, and oomph::PseudoSolidNodeUpdateElement< TTaylorHoodElement< 2 >, TPVDElement< 2, 3 > >.

1168  {
1169  unsigned nplot = 5;
1170  output(file_pt, nplot);
1171  }
void output(std::ostream &outfile)
Definition: navier_stokes_elements.h:1155

References oomph::NavierStokesEquations< DIM >::output().

◆ output() [2/4]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::output ( FILE *  file_pt,
const unsigned nplot 
)
virtual

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

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

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::QTaylorHoodElement< DIM >, oomph::QTaylorHoodElement< 2 >, oomph::QCrouzeixRaviartElement< DIM >, oomph::TTaylorHoodElement< DIM >, oomph::TCrouzeixRaviartElement< DIM >, oomph::PseudoSolidNodeUpdateElement< TCrouzeixRaviartElement< 2 >, TPVDBubbleEnrichedElement< 2, 3 > >, and oomph::PseudoSolidNodeUpdateElement< TTaylorHoodElement< 2 >, TPVDElement< 2, 3 > >.

804  {
805  // Vector of local coordinates
806  Vector<double> s(DIM);
807 
808  // Tecplot header info
809  fprintf(file_pt, "%s", tecplot_zone_string(nplot).c_str());
810 
811  // Loop over plot points
812  unsigned num_plot_points = nplot_points(nplot);
813  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
814  {
815  // Get local coordinates of plot point
816  get_s_plot(iplot, nplot, s);
817 
818  // Coordinates
819  for (unsigned i = 0; i < DIM; i++)
820  {
821  fprintf(file_pt, "%g ", interpolated_x(s, i));
822  }
823 
824  // Velocities
825  for (unsigned i = 0; i < DIM; i++)
826  {
827  fprintf(file_pt, "%g ", interpolated_u_nst(s, i));
828  }
829 
830  // Pressure
831  fprintf(file_pt, "%g \n", interpolated_p_nst(s));
832  }
833  fprintf(file_pt, "\n");
834 
835  // Write tecplot footer (e.g. FE connectivity lists)
836  write_tecplot_zone_footer(file_pt, nplot);
837  }

References DIM, i, and s.

◆ output() [3/4]

◆ output() [4/4]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::output ( std::ostream &  outfile,
const unsigned nplot 
)
virtual

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

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

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::TTaylorHoodElement< DIM >, oomph::TCrouzeixRaviartElement< DIM >, oomph::QTaylorHoodElement< DIM >, oomph::QTaylorHoodElement< 2 >, oomph::QCrouzeixRaviartElement< DIM >, oomph::PseudoSolidNodeUpdateElement< TCrouzeixRaviartElement< 2 >, TPVDBubbleEnrichedElement< 2, 3 > >, and oomph::PseudoSolidNodeUpdateElement< TTaylorHoodElement< 2 >, TPVDElement< 2, 3 > >.

758  {
759  // Vector of local coordinates
760  Vector<double> s(DIM);
761 
762  // Tecplot header info
763  outfile << tecplot_zone_string(nplot);
764 
765  // Loop over plot points
766  unsigned num_plot_points = nplot_points(nplot);
767  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
768  {
769  // Get local coordinates of plot point
770  get_s_plot(iplot, nplot, s);
771 
772  // Coordinates
773  for (unsigned i = 0; i < DIM; i++)
774  {
775  outfile << interpolated_x(s, i) << " ";
776  }
777 
778  // Velocities
779  for (unsigned i = 0; i < DIM; i++)
780  {
781  outfile << interpolated_u_nst(s, i) << " ";
782  }
783 
784  // Pressure
785  outfile << interpolated_p_nst(s) << " ";
786 
787  outfile << std::endl;
788  }
789  outfile << std::endl;
790 
791  // Write tecplot footer (e.g. FE connectivity lists)
792  write_tecplot_zone_footer(outfile, nplot);
793  }

References DIM, i, and s.

◆ output_fct() [1/2]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::output_fct ( std::ostream &  outfile,
const unsigned nplot,
const double time,
FiniteElement::UnsteadyExactSolutionFctPt  exact_soln_pt 
)
virtual

Output exact solution specified via function pointer at a given time and at a given number of plot points. Function prints as many components as are returned in solution Vector.

Output "exact" solution at a given time Solution is provided via function pointer. Plot at a given number of plot points. Function prints as many components as are returned in solution Vector.

Reimplemented from oomph::FiniteElement.

634  {
635  // Vector of local coordinates
636  Vector<double> s(DIM);
637 
638  // Vector for coordintes
639  Vector<double> x(DIM);
640 
641  // Tecplot header info
642  outfile << tecplot_zone_string(nplot);
643 
644  // Exact solution Vector
645  Vector<double> exact_soln;
646 
647  // Loop over plot points
648  unsigned num_plot_points = nplot_points(nplot);
649  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
650  {
651  // Get local coordinates of plot point
652  get_s_plot(iplot, nplot, s);
653 
654  // Get x position as Vector
655  interpolated_x(s, x);
656 
657  // Get exact solution at this point
658  (*exact_soln_pt)(time, x, exact_soln);
659 
660  // Output x,y,...
661  for (unsigned i = 0; i < DIM; i++)
662  {
663  outfile << x[i] << " ";
664  }
665 
666  // Output "exact solution"
667  for (unsigned i = 0; i < exact_soln.size(); i++)
668  {
669  outfile << exact_soln[i] << " ";
670  }
671 
672  outfile << std::endl;
673  }
674 
675  // Write tecplot footer (e.g. FE connectivity lists)
676  write_tecplot_zone_footer(outfile, nplot);
677  }

References DIM, s, and plotDoE::x.

◆ output_fct() [2/2]

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::output_fct ( std::ostream &  outfile,
const unsigned nplot,
FiniteElement::SteadyExactSolutionFctPt  exact_soln_pt 
)
virtual

Output exact solution specified via function pointer at a given number of plot points. Function prints as many components as are returned in solution Vector

Output "exact" solution Solution is provided via function pointer. Plot at a given number of plot points. Function prints as many components as are returned in solution Vector.

Reimplemented from oomph::FiniteElement.

577  {
578  // Vector of local coordinates
579  Vector<double> s(DIM);
580 
581  // Vector for coordintes
582  Vector<double> x(DIM);
583 
584  // Tecplot header info
585  outfile << tecplot_zone_string(nplot);
586 
587  // Exact solution Vector
588  Vector<double> exact_soln;
589 
590  // Loop over plot points
591  unsigned num_plot_points = nplot_points(nplot);
592  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
593  {
594  // Get local coordinates of plot point
595  get_s_plot(iplot, nplot, s);
596 
597  // Get x position as Vector
598  interpolated_x(s, x);
599 
600  // Get exact solution at this point
601  (*exact_soln_pt)(x, exact_soln);
602 
603  // Output x,y,...
604  for (unsigned i = 0; i < DIM; i++)
605  {
606  outfile << x[i] << " ";
607  }
608 
609  // Output "exact solution"
610  for (unsigned i = 0; i < exact_soln.size(); i++)
611  {
612  outfile << exact_soln[i] << " ";
613  }
614 
615  outfile << std::endl;
616  }
617 
618  // Write tecplot footer (e.g. FE connectivity lists)
619  write_tecplot_zone_footer(outfile, nplot);
620  }

References DIM, ProblemParameters::exact_soln(), i, s, and plotDoE::x.

◆ output_pressure_advection_diffusion_robin_elements()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::output_pressure_advection_diffusion_robin_elements ( std::ostream &  outfile)
inline

Output the FaceElements that apply the Robin boundary condition to the pressure advection diffusion problem required by Fp preconditioner

1451  {
1452  unsigned nel = Pressure_advection_diffusion_robin_element_pt.size();
1453  for (unsigned e = 0; e < nel; e++)
1454  {
1455  FaceElement* face_el_pt =
1457  outfile << "ZONE" << std::endl;
1458  Vector<double> unit_normal(DIM);
1459  Vector<double> x(DIM);
1460  Vector<double> s(DIM - 1);
1461  unsigned n = face_el_pt->integral_pt()->nweight();
1462  for (unsigned ipt = 0; ipt < n; ipt++)
1463  {
1464  for (unsigned i = 0; i < DIM - 1; i++)
1465  {
1466  s[i] = face_el_pt->integral_pt()->knot(ipt, i);
1467  }
1468  face_el_pt->interpolated_x(s, x);
1469  face_el_pt->outer_unit_normal(ipt, unit_normal);
1470  for (unsigned i = 0; i < DIM; i++)
1471  {
1472  outfile << x[i] << " ";
1473  }
1474  for (unsigned i = 0; i < DIM; i++)
1475  {
1476  outfile << unit_normal[i] << " ";
1477  }
1478  outfile << std::endl;
1479  }
1480  }
1481  }

References DIM, e(), i, oomph::FiniteElement::integral_pt(), oomph::FaceElement::interpolated_x(), oomph::Integral::knot(), n, oomph::Integral::nweight(), oomph::FaceElement::outer_unit_normal(), oomph::NavierStokesEquations< DIM >::Pressure_advection_diffusion_robin_element_pt, s, and plotDoE::x.

◆ output_veloc()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::output_veloc ( std::ostream &  outfile,
const unsigned nplot,
const unsigned t 
)

Output function: x,y,[z],u,v,[w] in tecplot format. nplot points in each coordinate direction at timestep t (t=0: present; t>0: previous timestep)

Output function: Velocities only x,y,[z],u,v,[w] in tecplot format at specified previous timestep (t=0: present; t>0: previous timestep). Specified number of plot points in each coordinate direction.

690  {
691  // Find number of nodes
692  unsigned n_node = nnode();
693 
694  // Local shape function
695  Shape psi(n_node);
696 
697  // Vectors of local coordinates and coords and velocities
698  Vector<double> s(DIM);
699  Vector<double> interpolated_x(DIM);
700  Vector<double> interpolated_u(DIM);
701 
702  // Tecplot header info
703  outfile << tecplot_zone_string(nplot);
704 
705  // Loop over plot points
706  unsigned num_plot_points = nplot_points(nplot);
707  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
708  {
709  // Get local coordinates of plot point
710  get_s_plot(iplot, nplot, s);
711 
712  // Get shape functions
713  shape(s, psi);
714 
715  // Loop over directions
716  for (unsigned i = 0; i < DIM; i++)
717  {
718  interpolated_x[i] = 0.0;
719  interpolated_u[i] = 0.0;
720  // Get the index at which velocity i is stored
721  unsigned u_nodal_index = u_index_nst(i);
722  // Loop over the local nodes and sum
723  for (unsigned l = 0; l < n_node; l++)
724  {
725  interpolated_u[i] += nodal_value(t, l, u_nodal_index) * psi[l];
726  interpolated_x[i] += nodal_position(t, l, i) * psi[l];
727  }
728  }
729 
730  // Coordinates
731  for (unsigned i = 0; i < DIM; i++)
732  {
733  outfile << interpolated_x[i] << " ";
734  }
735 
736  // Velocities
737  for (unsigned i = 0; i < DIM; i++)
738  {
739  outfile << interpolated_u[i] << " ";
740  }
741 
742  outfile << std::endl;
743  }
744 
745  // Write tecplot footer (e.g. FE connectivity lists)
746  write_tecplot_zone_footer(outfile, nplot);
747  }
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317

References DIM, i, s, oomph::OneDimLagrange::shape(), and plotPSD::t.

◆ output_vorticity()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::output_vorticity ( std::ostream &  outfile,
const unsigned nplot 
)

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

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

980  {
981  // Vector of local coordinates
982  Vector<double> s(DIM);
983 
984  // Create vorticity vector of the required size
985  Vector<double> vorticity;
986  if (DIM == 2)
987  {
988  vorticity.resize(1);
989  }
990  else if (DIM == 3)
991  {
992  vorticity.resize(3);
993  }
994  else
995  {
996  std::string error_message =
997  "Can't output vorticity in 1D - in fact, what do you mean\n";
998  error_message += "by the 1D Navier-Stokes equations?\n";
999 
1000  throw OomphLibError(
1002  }
1003 
1004  // Tecplot header info
1005  outfile << tecplot_zone_string(nplot);
1006 
1007  // Loop over plot points
1008  unsigned num_plot_points = nplot_points(nplot);
1009  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1010  {
1011  // Get local coordinates of plot point
1012  get_s_plot(iplot, nplot, s);
1013 
1014  // Coordinates
1015  for (unsigned i = 0; i < DIM; i++)
1016  {
1017  outfile << interpolated_x(s, i) << " ";
1018  }
1019 
1020  // Get vorticity
1021  get_vorticity(s, vorticity);
1022 
1023  if (DIM == 2)
1024  {
1025  outfile << vorticity[0];
1026  }
1027  else
1028  {
1029  outfile << vorticity[0] << " " << vorticity[1] << " " << vorticity[2]
1030  << " ";
1031  }
1032 
1033  outfile << std::endl;
1034  }
1035  outfile << std::endl;
1036 
1037  // Write tecplot footer (e.g. FE connectivity lists)
1038  write_tecplot_zone_footer(outfile, nplot);
1039  }
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286

◆ p_nodal_index_nst()

template<unsigned DIM>
virtual int oomph::NavierStokesEquations< DIM >::p_nodal_index_nst ( ) const
inlinevirtual

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

Implements oomph::TemplateFreeNavierStokesEquationsBase.

Reimplemented in oomph::TTaylorHoodElement< DIM >, oomph::QTaylorHoodElement< DIM >, and oomph::QTaylorHoodElement< 2 >.

937  {
939  }
static int Pressure_not_stored_at_node
Definition: navier_stokes_elements.h:418

References oomph::NavierStokesEquations< DIM >::Pressure_not_stored_at_node.

Referenced by oomph::NavierStokesEquations< DIM >::pin_all_non_pressure_dofs().

◆ p_nst() [1/2]

template<unsigned DIM>
virtual double oomph::NavierStokesEquations< DIM >::p_nst ( const unsigned n_p) const
pure virtual

◆ p_nst() [2/2]

template<unsigned DIM>
virtual double oomph::NavierStokesEquations< DIM >::p_nst ( const unsigned t,
const unsigned n_p 
) const
pure virtual

◆ pin_all_non_pressure_dofs()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::pin_all_non_pressure_dofs ( std::map< Data *, std::vector< int >> &  eqn_number_backup)
inlinevirtual

Pin all non-pressure dofs and backup eqn numbers.

Implements oomph::TemplateFreeNavierStokesEquationsBase.

1359  {
1360  // Loop over internal data and pin the values (having established that
1361  // pressure dofs aren't amongst those)
1362  unsigned nint = this->ninternal_data();
1363  for (unsigned j = 0; j < nint; j++)
1364  {
1365  Data* data_pt = this->internal_data_pt(j);
1366  if (eqn_number_backup[data_pt].size() == 0)
1367  {
1368  unsigned nvalue = data_pt->nvalue();
1369  eqn_number_backup[data_pt].resize(nvalue);
1370  for (unsigned i = 0; i < nvalue; i++)
1371  {
1372  // Backup
1373  eqn_number_backup[data_pt][i] = data_pt->eqn_number(i);
1374 
1375  // Pin everything
1376  data_pt->pin(i);
1377  }
1378  }
1379  }
1380 
1381  // Now deal with nodal values
1382  unsigned nnod = this->nnode();
1383  for (unsigned j = 0; j < nnod; j++)
1384  {
1385  Node* nod_pt = this->node_pt(j);
1386  if (eqn_number_backup[nod_pt].size() == 0)
1387  {
1388  unsigned nvalue = nod_pt->nvalue();
1389  eqn_number_backup[nod_pt].resize(nvalue);
1390  for (unsigned i = 0; i < nvalue; i++)
1391  {
1392  // Pin everything apart from the nodal pressure
1393  // value
1394  if (int(i) != this->p_nodal_index_nst())
1395  {
1396  // Backup
1397  eqn_number_backup[nod_pt][i] = nod_pt->eqn_number(i);
1398 
1399  // Pin
1400  nod_pt->pin(i);
1401  }
1402  // Else it's a pressure value
1403  else
1404  {
1405  // Exclude non-nodal pressure based elements
1406  if (this->p_nodal_index_nst() >= 0)
1407  {
1408  // Backup
1409  eqn_number_backup[nod_pt][i] = nod_pt->eqn_number(i);
1410  }
1411  }
1412  }
1413 
1414 
1415  // If it's a solid node deal with its positional data too
1416  SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
1417  if (solid_nod_pt != 0)
1418  {
1419  Data* solid_posn_data_pt = solid_nod_pt->variable_position_pt();
1420  if (eqn_number_backup[solid_posn_data_pt].size() == 0)
1421  {
1422  unsigned nvalue = solid_posn_data_pt->nvalue();
1423  eqn_number_backup[solid_posn_data_pt].resize(nvalue);
1424  for (unsigned i = 0; i < nvalue; i++)
1425  {
1426  // Backup
1427  eqn_number_backup[solid_posn_data_pt][i] =
1428  solid_posn_data_pt->eqn_number(i);
1429 
1430  // Pin
1431  solid_posn_data_pt->pin(i);
1432  }
1433  }
1434  }
1435  }
1436  }
1437  }
double size() const
Definition: elements.cc:4290
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:823
virtual int p_nodal_index_nst() const
Definition: navier_stokes_elements.h:936

References oomph::Data::eqn_number(), i, oomph::GeneralisedElement::internal_data_pt(), j, oomph::GeneralisedElement::ninternal_data(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::Data::nvalue(), oomph::NavierStokesEquations< DIM >::p_nodal_index_nst(), oomph::Data::pin(), oomph::FiniteElement::size(), and oomph::SolidNode::variable_position_pt().

◆ pinned_fp_pressure_eqn()

template<unsigned DIM>
int& oomph::NavierStokesEquations< DIM >::pinned_fp_pressure_eqn ( )
inlinevirtual

Global eqn number of pressure dof that's pinned in pressure adv diff problem

Implements oomph::TemplateFreeNavierStokesEquationsBase.

818  {
819  return Pinned_fp_pressure_eqn;
820  }

References oomph::NavierStokesEquations< DIM >::Pinned_fp_pressure_eqn.

◆ point_output_data()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::point_output_data ( const Vector< double > &  s,
Vector< double > &  data 
)
inlinevirtual

Output solution in data vector at local cordinates s: x,y [,z], u,v,[w], p

Reimplemented from oomph::FiniteElement.

1718  {
1719  // Dimension
1720  unsigned dim = s.size();
1721 
1722  // Resize data for values
1723  data.resize(2 * dim + 1);
1724 
1725  // Write values in the vector
1726  for (unsigned i = 0; i < dim; i++)
1727  {
1728  data[i] = interpolated_x(s, i);
1729  data[i + dim] = this->interpolated_u_nst(s, i);
1730  }
1731  data[2 * dim] = this->interpolated_p_nst(s);
1732  }
int data[]
Definition: Map_placement_new.cpp:1

References data, oomph::FiniteElement::dim(), i, oomph::NavierStokesEquations< DIM >::interpolated_p_nst(), oomph::NavierStokesEquations< DIM >::interpolated_u_nst(), oomph::FiniteElement::interpolated_x(), and s.

◆ pressure_integral()

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::pressure_integral

Integral of pressure over element.

Return pressure integrated over the element.

1559  {
1560  // Initialise
1561  double press_int = 0;
1562 
1563  // Set the value of n_intpt
1564  unsigned n_intpt = integral_pt()->nweight();
1565 
1566  // Set the Vector to hold local coordinates
1567  Vector<double> s(DIM);
1568 
1569  // Loop over the integration points
1570  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1571  {
1572  // Assign values of s
1573  for (unsigned i = 0; i < DIM; i++)
1574  {
1575  s[i] = integral_pt()->knot(ipt, i);
1576  }
1577 
1578  // Get the integral weight
1579  double w = integral_pt()->weight(ipt);
1580 
1581  // Get Jacobian of mapping
1582  double J = J_eulerian(s);
1583 
1584  // Premultiply the weights and the Jacobian
1585  double W = w * J;
1586 
1587  // Get pressure
1588  double press = interpolated_p_nst(s);
1589 
1590  // Add
1591  press_int += press * W;
1592  }
1593 
1594  return press_int;
1595  }

References DIM, i, J, s, w, and oomph::QuadTreeNames::W.

◆ pshape_nst() [1/2]

◆ pshape_nst() [2/2]

◆ re()

template<unsigned DIM>
const double& oomph::NavierStokesEquations< DIM >::re ( ) const
inline

Reynolds number.

704  {
705  return *Re_pt;
706  }

References oomph::NavierStokesEquations< DIM >::Re_pt.

◆ re_invfr()

template<unsigned DIM>
const double& oomph::NavierStokesEquations< DIM >::re_invfr ( ) const
inline

Global inverse Froude number.

754  {
755  return *ReInvFr_pt;
756  }

References oomph::NavierStokesEquations< DIM >::ReInvFr_pt.

◆ re_invfr_pt()

template<unsigned DIM>
double*& oomph::NavierStokesEquations< DIM >::re_invfr_pt ( )
inline

Pointer to global inverse Froude number.

760  {
761  return ReInvFr_pt;
762  }

References oomph::NavierStokesEquations< DIM >::ReInvFr_pt.

Referenced by oomph::RefineableNavierStokesEquations< DIM >::further_build().

◆ re_pt()

◆ re_st()

template<unsigned DIM>
const double& oomph::NavierStokesEquations< DIM >::re_st ( ) const
inline

Product of Reynolds and Strouhal number (=Womersley number)

710  {
711  return *ReSt_pt;
712  }

References oomph::NavierStokesEquations< DIM >::ReSt_pt.

◆ re_st_pt()

template<unsigned DIM>
double*& oomph::NavierStokesEquations< DIM >::re_st_pt ( )
inline

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

722  {
723  return ReSt_pt;
724  }

References oomph::NavierStokesEquations< DIM >::ReSt_pt.

Referenced by oomph::RefineableNavierStokesEquations< DIM >::further_build().

◆ scalar_name_paraview()

template<unsigned DIM>
std::string oomph::NavierStokesEquations< DIM >::scalar_name_paraview ( const unsigned i) const
inlinevirtual

Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.

Reimplemented from oomph::FiniteElement.

1128  {
1129  // Velocities
1130  if (i < DIM)
1131  {
1132  return "Velocity " + StringConversion::to_string(i);
1133  }
1134  // Preussre
1135  else if (i == DIM)
1136  {
1137  return "Pressure";
1138  }
1139  // Never get here
1140  else
1141  {
1142  std::stringstream error_stream;
1143  error_stream << "These Navier Stokes elements only store " << DIM + 1
1144  << " fields,\n"
1145  << "but i is currently " << i << std::endl;
1146  throw OomphLibError(
1147  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1148  // Dummy return
1149  return " ";
1150  }
1151  }
std::string to_string(T object, unsigned float_precision=8)
Definition: oomph_utilities.h:189

References DIM, i, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::StringConversion::to_string().

◆ scalar_value_fct_paraview()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::scalar_value_fct_paraview ( std::ofstream &  file_out,
const unsigned i,
const unsigned nplot,
const double time,
FiniteElement::UnsteadyExactSolutionFctPt  exact_soln_pt 
) const
inlinevirtual

Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specific element type.

Reimplemented from oomph::FiniteElement.

1073  {
1074 #ifdef PARANOID
1075  if (i > DIM)
1076  {
1077  // Create an output stream
1078  std::stringstream error_stream;
1079 
1080  // Create the error message
1081  error_stream << "These Navier Stokes elements only store " << DIM + 1
1082  << " fields, but i is currently " << i << std::endl;
1083 
1084  // Throw the error message
1085  throw OomphLibError(
1086  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1087  }
1088 #endif
1089 
1090  // Vector of local coordinates
1091  Vector<double> s(DIM + 1, 0.0);
1092 
1093  // Storage for the spatial coordinates
1094  Vector<double> spatial_coordinates(DIM, 0.0);
1095 
1096  // How many plot points do we have in total?
1097  unsigned num_plot_points = nplot_points_paraview(nplot);
1098 
1099  // Loop over plot points
1100  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1101  {
1102  // Get the local coordinates of the iplot-th plot point
1103  get_s_plot(iplot, nplot, s);
1104 
1105  // Loop over the spatial coordinates
1106  for (unsigned j = 0; j < DIM; j++)
1107  {
1108  // Assign the i-th spatial coordinate
1109  spatial_coordinates[j] = interpolated_x(s, j);
1110  }
1111 
1112  // Exact solution vector (here it's simply a scalar)
1113  Vector<double> exact_soln(DIM + 1, 0.0);
1114 
1115  // Get the exact solution at this point
1116  (*exact_soln_pt)(time, spatial_coordinates, exact_soln);
1117 
1118  // Output the interpolated solution value
1119  file_out << exact_soln[i] << std::endl;
1120  } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
1121  } // End of scalar_value_fct_paraview
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Definition: elements.h:2862

References DIM, ProblemParameters::exact_soln(), oomph::FiniteElement::get_s_plot(), i, oomph::FiniteElement::interpolated_x(), j, oomph::FiniteElement::nplot_points_paraview(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

◆ scalar_value_paraview()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::scalar_value_paraview ( std::ofstream &  file_out,
const unsigned i,
const unsigned nplot 
) const
inlinevirtual

Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specific element type.

Reimplemented from oomph::FiniteElement.

1024  {
1025  // Vector of local coordinates
1026  Vector<double> s(DIM);
1027 
1028 
1029  // Loop over plot points
1030  unsigned num_plot_points = nplot_points_paraview(nplot);
1031  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1032  {
1033  // Get local coordinates of plot point
1034  get_s_plot(iplot, nplot, s);
1035 
1036  // Velocities
1037  if (i < DIM)
1038  {
1039  file_out << interpolated_u_nst(s, i) << std::endl;
1040  }
1041 
1042  // Pressure
1043  else if (i == DIM)
1044  {
1045  file_out << interpolated_p_nst(s) << std::endl;
1046  }
1047 
1048  // Never get here
1049  else
1050  {
1051 #ifdef PARANOID
1052  std::stringstream error_stream;
1053  error_stream << "These Navier Stokes elements only store " << DIM + 1
1054  << " fields, "
1055  << "but i is currently " << i << std::endl;
1056  throw OomphLibError(error_stream.str(),
1059 #endif
1060  }
1061  }
1062  }

References DIM, oomph::FiniteElement::get_s_plot(), i, oomph::NavierStokesEquations< DIM >::interpolated_p_nst(), oomph::NavierStokesEquations< DIM >::interpolated_u_nst(), oomph::FiniteElement::nplot_points_paraview(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

◆ source_fct_for_pressure_adv_diff() [1/2]

template<unsigned DIM>
NavierStokesPressureAdvDiffSourceFctPt& oomph::NavierStokesEquations< DIM >::source_fct_for_pressure_adv_diff ( )
inline

Access function for the source-function pointer for pressure advection diffusion (used for validation only).

803  {
805  }

References oomph::NavierStokesEquations< DIM >::Press_adv_diff_source_fct_pt.

◆ source_fct_for_pressure_adv_diff() [2/2]

template<unsigned DIM>
NavierStokesPressureAdvDiffSourceFctPt oomph::NavierStokesEquations< DIM >::source_fct_for_pressure_adv_diff ( ) const
inline

Access function for the source-function pointer for pressure advection diffusion (used for validation only). Const version.

811  {
813  }

References oomph::NavierStokesEquations< DIM >::Press_adv_diff_source_fct_pt.

◆ source_fct_pt() [1/2]

template<unsigned DIM>
NavierStokesSourceFctPt& oomph::NavierStokesEquations< DIM >::source_fct_pt ( )
inline

Access function for the source-function pointer.

790  {
791  return Source_fct_pt;
792  }

References oomph::NavierStokesEquations< DIM >::Source_fct_pt.

Referenced by oomph::RefineableNavierStokesEquations< DIM >::further_build().

◆ source_fct_pt() [2/2]

template<unsigned DIM>
NavierStokesSourceFctPt oomph::NavierStokesEquations< DIM >::source_fct_pt ( ) const
inline

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

796  {
797  return Source_fct_pt;
798  }

References oomph::NavierStokesEquations< DIM >::Source_fct_pt.

◆ strain_rate()

template<unsigned DIM>
void oomph::NavierStokesEquations< DIM >::strain_rate ( const Vector< double > &  s,
DenseMatrix< double > &  strain_rate 
) const

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

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

1187  {
1188 #ifdef PARANOID
1189  if ((strainrate.ncol() != DIM) || (strainrate.nrow() != DIM))
1190  {
1191  std::ostringstream error_message;
1192  error_message << "The strain rate has incorrect dimensions "
1193  << strainrate.ncol() << " x " << strainrate.nrow()
1194  << " Not " << DIM << std::endl;
1195 
1196  throw OomphLibError(
1197  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1198  }
1199 #endif
1200 
1201  // Velocity gradient matrix
1202  DenseMatrix<double> dudx(DIM);
1203 
1204  // Find out how many nodes there are in the element
1205  unsigned n_node = nnode();
1206 
1207  // Set up memory for the shape and test functions
1208  Shape psi(n_node);
1209  DShape dpsidx(n_node, DIM);
1210 
1211  // Call the derivatives of the shape functions
1212  dshape_eulerian(s, psi, dpsidx);
1213 
1214  // Initialise to zero
1215  for (unsigned i = 0; i < DIM; i++)
1216  {
1217  for (unsigned j = 0; j < DIM; j++)
1218  {
1219  dudx(i, j) = 0.0;
1220  }
1221  }
1222 
1223  // Loop over veclocity components
1224  for (unsigned i = 0; i < DIM; i++)
1225  {
1226  // Get the index at which the i-th velocity is stored
1227  unsigned u_nodal_index = u_index_nst(i);
1228  // Loop over derivative directions
1229  for (unsigned j = 0; j < DIM; j++)
1230  {
1231  // Loop over nodes
1232  for (unsigned l = 0; l < n_node; l++)
1233  {
1234  dudx(i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1235  }
1236  }
1237  }
1238 
1239  // Loop over veclocity components
1240  for (unsigned i = 0; i < DIM; i++)
1241  {
1242  // Loop over derivative directions
1243  for (unsigned j = 0; j < DIM; j++)
1244  {
1245  strainrate(i, j) = 0.5 * (dudx(i, j) + dudx(j, i));
1246  }
1247  }
1248  }

Referenced by oomph::RefineableNavierStokesEquations< DIM >::get_Z2_flux(), oomph::TCrouzeixRaviartElement< DIM >::get_Z2_flux(), and oomph::TTaylorHoodElement< DIM >::get_Z2_flux().

◆ u_index_nst()

template<unsigned DIM>
virtual unsigned oomph::NavierStokesEquations< DIM >::u_index_nst ( const unsigned i) const
inlinevirtual

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

867  {
868  return i;
869  }

References i.

Referenced by oomph::NavierStokesEquations< DIM >::dinterpolated_u_nst_ddata(), oomph::RefineableNavierStokesEquations< DIM >::dinterpolated_u_nst_ddata(), oomph::NavierStokesEquations< DIM >::du_dt_nst(), oomph::BuoyantQCrouzeixRaviartElement< DIM >::fill_in_off_diagonal_jacobian_blocks_analytic(), oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::fill_in_off_diagonal_jacobian_blocks_analytic(), oomph::DoubleBuoyantQCrouzeixRaviartElement< DIM >::fill_in_off_diagonal_jacobian_blocks_by_fd(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::fill_in_off_diagonal_jacobian_blocks_by_fd(), oomph::RefineableQTaylorHoodElement< DIM >::get_interpolated_values(), oomph::RefineableQCrouzeixRaviartElement< DIM >::get_interpolated_values(), oomph::PRefineableQCrouzeixRaviartElement< DIM >::get_interpolated_values(), oomph::RefineableQTaylorHoodElement< DIM >::identify_load_data(), oomph::RefineableQCrouzeixRaviartElement< DIM >::identify_load_data(), oomph::NavierStokesEquations< DIM >::interpolated_dudx_nst(), oomph::NavierStokesEquations< DIM >::interpolated_u_nst(), oomph::MyTaylorHoodElement< DIM >::output(), oomph::MyTaylorHoodElement< DIM >::square_of_l2_norm(), and oomph::NavierStokesEquations< DIM >::u_nst().

◆ u_nst() [1/2]

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::u_nst ( const unsigned n,
const unsigned i 
) const
inline

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

849  {
850  return nodal_value(n, u_index_nst(i));
851  }

References i, n, oomph::FiniteElement::nodal_value(), and oomph::NavierStokesEquations< DIM >::u_index_nst().

◆ u_nst() [2/2]

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::u_nst ( const unsigned t,
const unsigned n,
const unsigned i 
) const
inline

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

856  {
857  return nodal_value(t, n, u_index_nst(i));
858  }

References i, n, oomph::FiniteElement::nodal_value(), plotPSD::t, and oomph::NavierStokesEquations< DIM >::u_index_nst().

◆ viscosity_ratio()

template<unsigned DIM>
const double& oomph::NavierStokesEquations< DIM >::viscosity_ratio ( ) const
inline

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

729  {
730  return *Viscosity_Ratio_pt;
731  }

References oomph::NavierStokesEquations< DIM >::Viscosity_Ratio_pt.

◆ viscosity_ratio_pt()

template<unsigned DIM>
double*& oomph::NavierStokesEquations< DIM >::viscosity_ratio_pt ( )
inline

Member Data Documentation

◆ ALE_is_disabled

template<unsigned DIM>
bool oomph::NavierStokesEquations< DIM >::ALE_is_disabled
protected

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

Referenced by oomph::NavierStokesEquations< DIM >::disable_ALE(), oomph::NavierStokesEquations< DIM >::enable_ALE(), and oomph::RefineableNavierStokesEquations< DIM >::further_build().

◆ Body_force_fct_pt

◆ Default_Gravity_vector

template<unsigned DIM>
Vector< double > oomph::NavierStokesEquations< DIM >::Default_Gravity_vector
staticprivate

Static default value for the gravity vector.

Navier-Stokes equations default gravity vector.

Referenced by oomph::NavierStokesEquations< DIM >::NavierStokesEquations().

◆ Default_Physical_Constant_Value

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::Default_Physical_Constant_Value = 0.0
staticprivate

Navier–Stokes equations static data.

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

Referenced by oomph::NavierStokesEquations< DIM >::NavierStokesEquations().

◆ Default_Physical_Ratio_Value

template<unsigned DIM>
double oomph::NavierStokesEquations< DIM >::Default_Physical_Ratio_Value = 1.0
staticprivate

Navier–Stokes equations static data.

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

Referenced by oomph::NavierStokesEquations< DIM >::NavierStokesEquations().

◆ Density_Ratio_pt

template<unsigned DIM>
double* oomph::NavierStokesEquations< DIM >::Density_Ratio_pt
protected

◆ G_pt

◆ Gamma

template<unsigned DIM>
Vector< double > oomph::NavierStokesEquations< DIM >::Gamma
static

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

Navier–Stokes equations static data.

////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////

◆ Pinned_fp_pressure_eqn

template<unsigned DIM>
int oomph::NavierStokesEquations< DIM >::Pinned_fp_pressure_eqn
protected

Global eqn number of pressure dof that's pinned in pressure advection diffusion problem (defaults to -1)

Referenced by oomph::NavierStokesEquations< DIM >::pinned_fp_pressure_eqn().

◆ Press_adv_diff_source_fct_pt

template<unsigned DIM>
NavierStokesPressureAdvDiffSourceFctPt oomph::NavierStokesEquations< DIM >::Press_adv_diff_source_fct_pt
protected

Pointer to source function pressure advection diffusion equation (only to be used during validation)

Referenced by oomph::NavierStokesEquations< DIM >::source_fct_for_pressure_adv_diff().

◆ Pressure_advection_diffusion_robin_element_pt

◆ Pressure_not_stored_at_node

template<unsigned DIM>
int oomph::NavierStokesEquations< DIM >::Pressure_not_stored_at_node = -100
staticprivate

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

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

Referenced by oomph::NavierStokesEquations< DIM >::p_nodal_index_nst().

◆ Re_pt

◆ ReInvFr_pt

template<unsigned DIM>
double* oomph::NavierStokesEquations< DIM >::ReInvFr_pt
protected

◆ ReSt_pt

◆ Source_fct_pt

◆ Viscosity_Ratio_pt

template<unsigned DIM>
double* oomph::NavierStokesEquations< DIM >::Viscosity_Ratio_pt
protected

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