Public Member Functions | Static Public Member Functions | Protected Member Functions | List of all members
oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM > Class Template Referenceabstract

#include <generalised_newtonian_refineable_navier_stokes_elements.h>

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

Public Member Functions

 RefineableGeneralisedNewtonianNavierStokesEquations ()
 Constructor. More...
 
virtual Node * pressure_node_pt (const unsigned &n_p)
 
void get_pressure_and_velocity_mass_matrix_diagonal (Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
 
unsigned num_Z2_flux_terms ()
 Number of 'flux' terms for Z2 error estimation. More...
 
void get_Z2_flux (const Vector< double > &s, Vector< double > &flux)
 
void further_build ()
 Further build, pass the pointers down to the sons. More...
 
void dinterpolated_u_nst_ddata (const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
 
- Public Member Functions inherited from oomph::GeneralisedNewtonianNavierStokesEquations< DIM >
 GeneralisedNewtonianNavierStokesEquations ()
 
const double & re () const
 Reynolds number. More...
 
const double & re_st () const
 Product of Reynolds and Strouhal number (=Womersley number) More...
 
double *& re_pt ()
 Pointer to Reynolds number. More...
 
double *& re_st_pt ()
 Pointer to product of Reynolds and Strouhal number (=Womersley number) More...
 
const double & viscosity_ratio () const
 
double *& viscosity_ratio_pt ()
 Pointer to Viscosity Ratio. More...
 
const double & density_ratio () const
 
double *& density_ratio_pt ()
 Pointer to Density ratio. More...
 
const double & re_invfr () const
 Global inverse Froude number. More...
 
double *& re_invfr_pt ()
 Pointer to global inverse Froude number. More...
 
const Vector< double > & g () const
 Vector of gravitational components. More...
 
Vector< double > *& g_pt ()
 Pointer to Vector of gravitational components. More...
 
NavierStokesBodyForceFctPt & body_force_fct_pt ()
 Access function for the body-force pointer. More...
 
NavierStokesBodyForceFctPt body_force_fct_pt () const
 Access function for the body-force pointer. Const version. More...
 
NavierStokesSourceFctPt & source_fct_pt ()
 Access function for the source-function pointer. More...
 
NavierStokesSourceFctPt source_fct_pt () const
 Access function for the source-function pointer. Const version. More...
 
GeneralisedNewtonianConstitutiveEquation< DIM > *& constitutive_eqn_pt ()
 Access function for the constitutive equation pointer. More...
 
void use_current_strainrate_to_compute_second_invariant ()
 Switch to fully implict evaluation of non-Newtonian const eqn. More...
 
void use_extrapolated_strainrate_to_compute_second_invariant ()
 Use extrapolation for non-Newtonian const eqn. More...
 
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
 
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...
 
void max_and_min_invariant_and_viscosity (double &min_invariant, double &max_invariant, double &min_viscosity, double &max_viscosity) const
 
double dissipation () const
 Return integral of dissipation over element. More...
 
double dissipation (const Vector< double > &s) const
 Return dissipation at local coordinate s. More...
 
void get_vorticity (const Vector< double > &s, Vector< double > &vorticity) const
 Compute the vorticity vector at local coordinate s. More...
 
double kin_energy () const
 Get integral of kinetic energy over element. More...
 
double d_kin_energy_dt () const
 Get integral of time derivative of kinetic energy over element. More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i) More...
 
void strain_rate (const unsigned &t, const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 
virtual void extrapolated_strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 
virtual void extrapolated_strain_rate (const unsigned &ipt, DenseMatrix< double > &strain_rate) const
 
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)
 
unsigned nscalar_paraview () const
 
void scalar_value_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
 
std::string scalar_name_paraview (const unsigned &i) const
 
void output (std::ostream &outfile)
 
void output (std::ostream &outfile, const unsigned &nplot)
 
void output (FILE *file_pt)
 
void output (FILE *file_pt, const unsigned &nplot)
 
void full_output (std::ostream &outfile)
 
void full_output (std::ostream &outfile, const unsigned &nplot)
 
void output_veloc (std::ostream &outfile, const unsigned &nplot, const unsigned &t)
 
void output_vorticity (std::ostream &outfile, const unsigned &nplot)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 
void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
void 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 pin_all_non_pressure_dofs (std::map< Data *, std::vector< int >> &eqn_number_backup)
 Pin all non-pressure dofs and backup eqn numbers. More...
 
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 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...
 
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, 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)
 
MacroElement * macro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void describe_nodal_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void assign_all_generic_local_eqn_numbers (const bool &store_local_dof_pt)
 
Node *& node_pt (const unsigned &n)
 Return a pointer to the local node n. More...
 
Node *const & node_pt (const unsigned &n) const
 Return a pointer to the local node n (const version) More...
 
unsigned nnode () const
 Return the number of nodes. More...
 
virtual unsigned nnode_1d () const
 
double raw_nodal_position (const unsigned &n, const unsigned &i) const
 
double raw_nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position (const unsigned &n, const unsigned &i) const
 
double nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double dnodal_position_dt (const unsigned &n, const unsigned &i) const
 Return the i-th component of nodal velocity: dx/dt at local node n. More...
 
double dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual unsigned required_nvalue (const unsigned &n) const
 
unsigned nnodal_position_type () const
 
bool has_hanging_nodes () const
 
unsigned nodal_dimension () const
 Return the required Eulerian dimension of the nodes in this element. More...
 
virtual Node * construct_node (const unsigned &n)
 
virtual Node * construct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual Node * construct_boundary_node (const unsigned &n)
 
virtual Node * construct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
int get_node_number (Node *const &node_pt) const
 
virtual Node * get_node_at_local_coordinate (const Vector< double > &s) const
 
double raw_nodal_value (const unsigned &n, const unsigned &i) const
 
double raw_nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
unsigned dim () const
 
virtual ElementGeometry::ElementGeometry element_geometry () const
 Return the geometry type of the element (either Q or T usually). More...
 
virtual double interpolated_x (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated coordinate x[i] at local coordinate s. More...
 
virtual double interpolated_x (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_x (const Vector< double > &s, Vector< double > &x) const
 Return FE interpolated position x[] at local coordinate s as Vector. More...
 
virtual void interpolated_x (const unsigned &t, const Vector< double > &s, Vector< double > &x) const
 
virtual double interpolated_dxdt (const Vector< double > &s, const unsigned &i, const unsigned &t)
 
virtual void interpolated_dxdt (const Vector< double > &s, const unsigned &t, Vector< double > &dxdt)
 
unsigned ngeom_data () const
 
Data * geom_data_pt (const unsigned &j)
 
void position (const Vector< double > &zeta, Vector< double > &r) const
 
void position (const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
 
void dposition_dt (const Vector< double > &zeta, const unsigned &t, Vector< double > &drdt)
 
virtual double zeta_nodal (const unsigned &n, const unsigned &k, const unsigned &i) const
 
void interpolated_zeta (const Vector< double > &s, Vector< double > &zeta) const
 
void locate_zeta (const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
 
virtual void node_update ()
 
virtual void identify_geometric_data (std::set< Data * > &geometric_data_pt)
 
virtual double s_min () const
 Min value of local coordinate. More...
 
virtual double s_max () const
 Max. value of local coordinate. More...
 
double size () const
 
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 scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
virtual void output (const unsigned &t, std::ostream &outfile, const unsigned &n_plot) const
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, const SolutionFunctorBase &exact_soln) const
 Output a time-dependent exact solution over the element. More...
 
virtual void get_s_plot (const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
 
virtual std::string tecplot_zone_string (const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (std::ostream &outfile, const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (FILE *file_pt, const unsigned &nplot) const
 
virtual unsigned nplot_points (const unsigned &nplot) const
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_abs_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error)
 
void integrate_fct (FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt, Vector< double > &integral)
 Integrate Vector-valued function over element. More...
 
void integrate_fct (FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt, const double &time, Vector< double > &integral)
 Integrate Vector-valued time-dep function over element. More...
 
virtual void build_face_element (const int &face_index, FaceElement *face_element_pt)
 
virtual unsigned self_test ()
 
virtual unsigned get_bulk_node_number (const int &face_index, const unsigned &i) const
 
virtual int face_outer_unit_normal_sign (const int &face_index) const
 Get the sign of the outer unit normal on the face given by face_index. More...
 
virtual unsigned nnode_on_face () const
 
void face_node_number_error_check (const unsigned &i) const
 Range check for face node numbers. More...
 
virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt (const int &face_index) const
 
virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt (const int &face_index) const
 
- Public Member Functions inherited from oomph::GeneralisedElement
 GeneralisedElement ()
 Constructor: Initialise all pointers and all values to zero. More...
 
virtual ~GeneralisedElement ()
 Virtual destructor to clean up any memory allocated by the object. More...
 
 GeneralisedElement (const GeneralisedElement &)=delete
 Broken copy constructor. More...
 
void operator= (const GeneralisedElement &)=delete
 Broken assignment operator. More...
 
Data *& internal_data_pt (const unsigned &i)
 Return a pointer to i-th internal data object. More...
 
Data *const & internal_data_pt (const unsigned &i) const
 Return a pointer to i-th internal data object (const version) More...
 
Data *& external_data_pt (const unsigned &i)
 Return a pointer to i-th external data object. More...
 
Data *const & external_data_pt (const unsigned &i) const
 Return a pointer to i-th external data object (const version) More...
 
unsigned long eqn_number (const unsigned &ieqn_local) const
 
int local_eqn_number (const unsigned long &ieqn_global) const
 
unsigned add_external_data (Data *const &data_pt, const bool &fd=true)
 
bool external_data_fd (const unsigned &i) const
 
void exclude_external_data_fd (const unsigned &i)
 
void include_external_data_fd (const unsigned &i)
 
void flush_external_data ()
 Flush all external data. More...
 
void flush_external_data (Data *const &data_pt)
 Flush the object addressed by data_pt from the external data array. More...
 
unsigned ninternal_data () const
 Return the number of internal data objects. More...
 
unsigned nexternal_data () const
 Return the number of external data objects. More...
 
unsigned ndof () const
 Return the number of equations/dofs in the element. More...
 
void dof_vector (const unsigned &t, Vector< double > &dof)
 Return the vector of dof values at time level t. More...
 
void dof_pt_vector (Vector< double * > &dof_pt)
 Return the vector of pointers to dof values. More...
 
void set_internal_data_time_stepper (const unsigned &i, TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void assign_internal_eqn_numbers (unsigned long &global_number, Vector< double * > &Dof_pt)
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void add_internal_value_pt_to_map (std::map< unsigned, double * > &map_of_value_pt)
 
virtual void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void complete_setup_of_dependencies ()
 
virtual void get_residuals (Vector< double > &residuals)
 
virtual void get_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
virtual void get_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void get_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void get_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void get_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void get_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void get_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void get_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void get_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_norm (double &norm)
 
virtual unsigned ndof_types () const
 
virtual void get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
 
- Public Member Functions inherited from oomph::GeomObject
 GeomObject ()
 Default constructor. More...
 
 GeomObject (const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)
 
 GeomObject (const GeomObject &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const GeomObject &)=delete
 Broken assignment operator. More...
 
virtual ~GeomObject ()
 (Empty) destructor More...
 
unsigned nlagrangian () const
 Access function to # of Lagrangian coordinates. More...
 
unsigned ndim () const
 Access function to # of Eulerian coordinates. More...
 
void set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
 Set # of Lagrangian and Eulerian coordinates. More...
 
TimeStepper *& time_stepper_pt ()
 
TimeStepper * time_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 
- Public Member Functions inherited from oomph::GeneralisedNewtonianTemplateFreeNavierStokesEquationsBase
 GeneralisedNewtonianTemplateFreeNavierStokesEquationsBase ()
 Constructor (empty) More...
 
virtual ~GeneralisedNewtonianTemplateFreeNavierStokesEquationsBase ()
 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...
 
- Public Member Functions inherited from oomph::RefineableElement
 RefineableElement ()
 
virtual ~RefineableElement ()
 Destructor, delete the allocated storage for the hanging equations. More...
 
 RefineableElement (const RefineableElement &)=delete
 Broken copy constructor. More...
 
void operator= (const RefineableElement &)=delete
 Broken assignment operator. More...
 
Tree * tree_pt ()
 Access function: Pointer to quadtree representation of this element. More...
 
void set_tree_pt (Tree *my_tree_pt)
 Set pointer to quadtree representation of this element. More...
 
virtual unsigned required_nsons () const
 
bool refinement_is_enabled ()
 Flag to indicate suppression of any refinement. More...
 
void disable_refinement ()
 Suppress of any refinement for this element. More...
 
void enable_refinement ()
 Emnable refinement for this element. More...
 
template<class ELEMENT >
void split (Vector< ELEMENT * > &son_pt) const
 
int local_hang_eqn (Node *const &node_pt, const unsigned &i)
 
virtual void build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)=0
 
void set_refinement_level (const int &refine_level)
 Set the refinement level. More...
 
unsigned refinement_level () const
 Return the Refinement level. More...
 
void select_for_refinement ()
 Select the element for refinement. More...
 
void deselect_for_refinement ()
 Deselect the element for refinement. More...
 
void select_sons_for_unrefinement ()
 Unrefinement will be performed by merging the four sons of this element. More...
 
void deselect_sons_for_unrefinement ()
 
bool to_be_refined ()
 Has the element been selected for refinement? More...
 
bool sons_to_be_unrefined ()
 Has the element been selected for unrefinement? More...
 
virtual void rebuild_from_sons (Mesh *&mesh_pt)=0
 
virtual void unbuild ()
 
virtual void deactivate_element ()
 
virtual bool nodes_built ()
 Return true if all the nodes have been built, false if not. More...
 
long number () const
 Element number (for debugging/plotting) More...
 
void set_number (const long &mynumber)
 Set element number (for debugging/plotting) More...
 
virtual unsigned ncont_interpolated_values () const =0
 
virtual void get_interpolated_values (const Vector< double > &s, Vector< double > &values)
 
virtual void get_interpolated_values (const unsigned &t, const Vector< double > &s, Vector< double > &values)=0
 
virtual Node * interpolating_node_pt (const unsigned &n, const int &value_id)
 
virtual double local_one_d_fraction_of_interpolating_node (const unsigned &n1d, const unsigned &i, const int &value_id)
 
virtual Node * get_interpolating_node_at_local_coordinate (const Vector< double > &s, const int &value_id)
 
virtual unsigned ninterpolating_node (const int &value_id)
 
virtual unsigned ninterpolating_node_1d (const int &value_id)
 
virtual void interpolating_basis (const Vector< double > &s, Shape &psi, const int &value_id) const
 
virtual void check_integrity (double &max_error)=0
 
void identify_field_data_for_interactions (std::set< std::pair< Data *, unsigned >> &paired_field_data)
 
void assign_nodal_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual RefineableElement * root_element_pt ()
 
virtual RefineableElement * father_element_pt () const
 Return a pointer to the father element. More...
 
void get_father_at_refinement_level (unsigned &refinement_level, RefineableElement *&father_at_reflevel_pt)
 
virtual void initial_setup (Tree *const &adopted_father_pt=0, const unsigned &initial_p_order=0)
 
virtual void pre_build (Mesh *&mesh_pt, Vector< Node * > &new_node_pt)
 Pre-build the element. More...
 
virtual void setup_hanging_nodes (Vector< std::ofstream * > &output_stream)
 
virtual void further_setup_hanging_nodes ()
 
void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
unsigned nshape_controlling_nodes ()
 
std::map< Node *, unsigned > shape_controlling_node_lookup ()
 
- Public Member Functions inherited from oomph::ElementWithZ2ErrorEstimator
 ElementWithZ2ErrorEstimator ()
 Default empty constructor. More...
 
 ElementWithZ2ErrorEstimator (const ElementWithZ2ErrorEstimator &)=delete
 Broken copy constructor. More...
 
void operator= (const ElementWithZ2ErrorEstimator &)=delete
 Broken assignment operator. More...
 
virtual unsigned ncompound_fluxes ()
 
virtual void compute_exact_Z2_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
 
virtual void get_Z2_compound_flux_indices (Vector< unsigned > &flux_index)
 
virtual unsigned nvertex_node () const =0
 Number of vertex nodes in the element. More...
 
virtual Node * vertex_node_pt (const unsigned &j) const =0
 
virtual unsigned nrecovery_order ()=0
 Order of recovery shape functions. More...
 
virtual double geometric_jacobian (const Vector< double > &x)
 

Static Public Member Functions

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

Protected Member Functions

virtual void unpin_elemental_pressure_dofs ()=0
 Unpin all pressure dofs in the element. More...
 
virtual void pin_elemental_redundant_nodal_pressure_dofs ()
 
void fill_in_generic_residual_contribution_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
 
virtual void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
- Protected Member Functions inherited from oomph::GeneralisedNewtonianNavierStokesEquations< DIM >
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 double dpshape_and_dptest_eulerian_nst (const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) 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_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
template<unsigned DIM>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double invert_jacobian_mapping (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &inverse_jacobian) const
 
virtual void dJ_eulerian_dnodal_coordinates (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<unsigned DIM>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
virtual void d_dshape_eulerian_dnodal_coordinates (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<unsigned DIM>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
virtual void transform_derivatives (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
void transform_derivatives_diagonal (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
virtual void transform_second_derivatives (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
void fill_in_jacobian_from_nodal_by_fd (DenseMatrix< double > &jacobian)
 
virtual void update_before_nodal_fd ()
 
virtual void reset_after_nodal_fd ()
 
virtual void update_in_nodal_fd (const unsigned &i)
 
virtual void reset_in_nodal_fd (const unsigned &i)
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Zero-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 One-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Two-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
- Protected Member Functions inherited from oomph::GeneralisedElement
unsigned add_internal_data (Data *const &data_pt, const bool &fd=true)
 
bool internal_data_fd (const unsigned &i) const
 
void exclude_internal_data_fd (const unsigned &i)
 
void include_internal_data_fd (const unsigned &i)
 
void clear_global_eqn_numbers ()
 
void add_global_eqn_numbers (std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)
 
virtual void assign_internal_and_external_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void assign_additional_local_eqn_numbers ()
 
int internal_local_eqn (const unsigned &i, const unsigned &j) const
 
int external_local_eqn (const unsigned &i, const unsigned &j)
 
void fill_in_jacobian_from_internal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_internal_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
virtual void update_before_internal_fd ()
 
virtual void reset_after_internal_fd ()
 
virtual void update_in_internal_fd (const unsigned &i)
 
virtual void reset_in_internal_fd (const unsigned &i)
 
virtual void update_before_external_fd ()
 
virtual void reset_after_external_fd ()
 
virtual void update_in_external_fd (const unsigned &i)
 
virtual void reset_in_external_fd (const unsigned &i)
 
virtual void fill_in_contribution_to_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void fill_in_contribution_to_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void fill_in_contribution_to_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
- Protected Member Functions inherited from oomph::RefineableElement
void assemble_local_to_eulerian_jacobian (const DShape &dpsids, DenseMatrix< double > &jacobian) const
 
void assemble_local_to_eulerian_jacobian2 (const DShape &d2psids, DenseMatrix< double > &jacobian2) const
 
void assemble_eulerian_base_vectors (const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
 
double local_to_eulerian_mapping_diagonal (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
void assign_hanging_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign the local equation numbers for hanging node variables. More...
 
virtual void fill_in_jacobian_from_nodal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 

Additional Inherited Members

- Public Types inherited from oomph::GeneralisedNewtonianNavierStokesEquations< DIM >
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 > &)
 
- Static Public Attributes inherited from oomph::GeneralisedNewtonianNavierStokesEquations< DIM >
static Vector< double > Gamma
 Vector to decide whether the stress-divergence form is used or not. More...
 
- Static Public Attributes inherited from oomph::FiniteElement
static double Tolerance_for_singular_jacobian = 1.0e-16
 Tolerance below which the jacobian is considered singular. More...
 
static bool Accept_negative_jacobian = false
 
static bool Suppress_output_while_checking_for_inverted_elements
 
- Static Public Attributes inherited from oomph::GeneralisedElement
static bool Suppress_warning_about_repeated_internal_data
 
static bool Suppress_warning_about_repeated_external_data = true
 
static double Default_fd_jacobian_step = 1.0e-8
 
- Static Protected Member Functions inherited from oomph::RefineableElement
static void check_value_id (const int &n_continuously_interpolated_values, const int &value_id)
 
- Protected Attributes inherited from oomph::GeneralisedNewtonianNavierStokesEquations< DIM >
double * Viscosity_Ratio_pt
 
double * Density_Ratio_pt
 
double * Re_pt
 Pointer to global Reynolds number. More...
 
double * ReSt_pt
 Pointer to global Reynolds number x Strouhal number (=Womersley) More...
 
double * ReInvFr_pt
 
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...
 
GeneralisedNewtonianConstitutiveEquation< DIM > * Constitutive_eqn_pt
 Pointer to the generalised Newtonian constitutive equation. More...
 
bool Use_extrapolated_strainrate_to_compute_second_invariant
 
bool ALE_is_disabled
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElement * Macro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepper * Geom_object_time_stepper_pt
 
- Protected Attributes inherited from oomph::RefineableElement
Tree * Tree_pt
 A pointer to a general tree object. More...
 
unsigned Refine_level
 Refinement level. More...
 
bool To_be_refined
 Flag for refinement. More...
 
bool Refinement_is_enabled
 Flag to indicate suppression of any refinement. More...
 
bool Sons_to_be_unrefined
 Flag for unrefinement. More...
 
long Number
 Global element number – for plotting/validation purposes. More...
 
- Static Protected Attributes inherited from oomph::GeneralisedNewtonianNavierStokesEquations< DIM >
static bool Pre_multiply_by_viscosity_ratio = false
 
- Static Protected Attributes inherited from oomph::FiniteElement
static const unsigned Default_Initial_Nvalue = 0
 Default value for the number of values at a node. More...
 
static const double Node_location_tolerance = 1.0e-14
 
static const unsigned N2deriv [] = {0, 1, 3, 6}
 
- Static Protected Attributes inherited from oomph::GeneralisedElement
static DenseMatrix< double > Dummy_matrix
 
static std::deque< double * > Dof_pt_deque
 
- Static Protected Attributes inherited from oomph::RefineableElement
static double Max_integrity_tolerance = 1.0e-8
 Max. allowed discrepancy in element integrity check. More...
 

Detailed Description

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

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// Refineable version of the Navier–Stokes equations

Constructor & Destructor Documentation

◆ RefineableGeneralisedNewtonianNavierStokesEquations()

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

Constructor.

72  : GeneralisedNewtonianNavierStokesEquations<DIM>(),
73  RefineableElement(),
74  ElementWithZ2ErrorEstimator()
75  {
76  }
oomph::ElementWithZ2ErrorEstimator::ElementWithZ2ErrorEstimator
ElementWithZ2ErrorEstimator()
Default empty constructor.
Definition: error_estimator.h:82
oomph::RefineableElement::RefineableElement
RefineableElement()
Definition: refineable_elements.h:188

Member Function Documentation

◆ dinterpolated_u_nst_ddata()

template<unsigned DIM>
void oomph::RefineableGeneralisedNewtonianNavierStokesEquations< 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. Overload the non-refineable version to take account of hanging node information

Reimplemented from oomph::GeneralisedNewtonianNavierStokesEquations< DIM >.

227  {
228  // Find number of nodes
229  unsigned n_node = this->nnode();
230  // Local shape function
231  Shape psi(n_node);
232  // Find values of shape function at the given local coordinate
233  this->shape(s, psi);
234 
235  // Find the index at which the velocity component is stored
236  const unsigned u_nodal_index = this->u_index_nst(i);
237 
238  // Storage for hang info pointer
239  HangInfo* hang_info_pt = 0;
240  // Storage for global equation
241  int global_eqn = 0;
242 
243  // Find the number of dofs associated with interpolated u
244  unsigned n_u_dof = 0;
245  for (unsigned l = 0; l < n_node; l++)
246  {
247  unsigned n_master = 1;
248 
249  // Local bool (is the node hanging)
250  bool is_node_hanging = this->node_pt(l)->is_hanging();
251 
252  // If the node is hanging, get the number of master nodes
253  if (is_node_hanging)
254  {
255  hang_info_pt = this->node_pt(l)->hanging_pt();
256  n_master = hang_info_pt->nmaster();
257  }
258  // Otherwise there is just one master node, the node itself
259  else
260  {
261  n_master = 1;
262  }
263 
264  // Loop over the master nodes
265  for (unsigned m = 0; m < n_master; m++)
266  {
267  // Get the equation number
268  if (is_node_hanging)
269  {
270  // Get the equation number from the master node
271  global_eqn =
272  hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
273  }
274  else
275  {
276  // Global equation number
277  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
278  }
279 
280  // If it's positive add to the count
281  if (global_eqn >= 0)
282  {
283  ++n_u_dof;
284  }
285  }
286  }
287 
288  // Now resize the storage schemes
289  du_ddata.resize(n_u_dof, 0.0);
290  global_eqn_number.resize(n_u_dof, 0);
291 
292  // Loop over th nodes again and set the derivatives
293  unsigned count = 0;
294  // Loop over the local nodes and sum
295  for (unsigned l = 0; l < n_node; l++)
296  {
297  unsigned n_master = 1;
298  double hang_weight = 1.0;
299 
300  // Local bool (is the node hanging)
301  bool is_node_hanging = this->node_pt(l)->is_hanging();
302 
303  // If the node is hanging, get the number of master nodes
304  if (is_node_hanging)
305  {
306  hang_info_pt = this->node_pt(l)->hanging_pt();
307  n_master = hang_info_pt->nmaster();
308  }
309  // Otherwise there is just one master node, the node itself
310  else
311  {
312  n_master = 1;
313  }
314 
315  // Loop over the master nodes
316  for (unsigned m = 0; m < n_master; m++)
317  {
318  // If the node is hanging get weight from master node
319  if (is_node_hanging)
320  {
321  // Get the hang weight from the master node
322  hang_weight = hang_info_pt->master_weight(m);
323  }
324  else
325  {
326  // Node contributes with full weight
327  hang_weight = 1.0;
328  }
329 
330  // Get the equation number
331  if (is_node_hanging)
332  {
333  // Get the equation number from the master node
334  global_eqn =
335  hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
336  }
337  else
338  {
339  // Local equation number
340  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
341  }
342 
343  if (global_eqn >= 0)
344  {
345  // Set the global equation number
346  global_eqn_number[count] = global_eqn;
347  // Set the derivative with respect to the unknown
348  du_ddata[count] = psi[l] * hang_weight;
349  // Increase the counter
350  ++count;
351  }
352  }
353  }
354  }
oomph::Data::eqn_number
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::shape
virtual void shape(const Vector< double > &s, Shape &psi) const =0
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::GeneralisedNewtonianNavierStokesEquations::u_index_nst
virtual unsigned u_index_nst(const unsigned &i) const
Definition: generalised_newtonian_navier_stokes_elements.h:580
oomph::Node::hanging_pt
HangInfo *const & hanging_pt() const
Definition: nodes.h:1228
oomph::Node::is_hanging
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
m
int * m
Definition: level2_cplx_impl.h:294

References oomph::Data::eqn_number(), oomph::Node::hanging_pt(), oomph::Node::is_hanging(), m, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::FiniteElement::shape(), and oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::u_index_nst().

◆ fill_in_generic_residual_contribution_nst()

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

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

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

Reimplemented from oomph::GeneralisedNewtonianNavierStokesEquations< DIM >.

296  {
297  // Find out how many nodes there are
298  unsigned n_node = nnode();
299 
300  // Get continuous time from timestepper of first node
301  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
302 
303  // Find out how many pressure dofs there are
304  unsigned n_pres = this->npres_nst();
305 
306  // Get the indices at which the velocity components are stored
307  unsigned u_nodal_index[DIM];
308  for (unsigned i = 0; i < DIM; i++)
309  {
310  u_nodal_index[i] = this->u_index_nst(i);
311  }
312 
313  // Which nodal value represents the pressure? (Negative if pressure
314  // is not based on nodal interpolation).
315  int p_index = this->p_nodal_index_nst();
316 
317  // Local array of booleans that are true if the l-th pressure value is
318  // hanging (avoid repeated virtual function calls)
319  bool pressure_dof_is_hanging[n_pres];
320  // If the pressure is stored at a node
321  if (p_index >= 0)
322  {
323  // Read out whether the pressure is hanging
324  for (unsigned l = 0; l < n_pres; ++l)
325  {
326  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
327  }
328  }
329  // Otherwise the pressure is not stored at a node and so cannot hang
330  else
331  {
332  for (unsigned l = 0; l < n_pres; ++l)
333  {
334  pressure_dof_is_hanging[l] = false;
335  }
336  }
337 
338  // Set up memory for the shape and test functions
339  Shape psif(n_node), testf(n_node);
340  DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
341 
342 
343  // Set up memory for pressure shape and test functions
344  Shape psip(n_pres), testp(n_pres);
345 
346  // Set the value of n_intpt
347  unsigned n_intpt = integral_pt()->nweight();
348 
349  // Set the Vector to hold local coordinates
350  Vector<double> s(DIM);
351 
352  // Get Physical Variables from Element
353  // Reynolds number must be multiplied by the density ratio
354  double scaled_re = this->re() * this->density_ratio();
355  double scaled_re_st = this->re_st() * this->density_ratio();
356  double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
357  double visc_ratio = this->viscosity_ratio();
358  Vector<double> G = this->g();
359 
360  // Integers that store the local equations and unknowns
361  int local_eqn = 0, local_unknown = 0;
362 
363  // Pointers to hang info objects
364  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
365 
366  // Local boolean for ALE (or not)
367  bool ALE_is_disabled_flag = this->ALE_is_disabled;
368 
369  // Loop over the integration points
370  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
371  {
372  // Assign values of s
373  for (unsigned i = 0; i < DIM; i++)
374  {
375  s[i] = integral_pt()->knot(ipt, i);
376  }
377 
378  // Get the integral weight
379  double w = integral_pt()->weight(ipt);
380 
381  // Call the derivatives of the shape and test functions
382  double J = this->dshape_and_dtest_eulerian_at_knot_nst(
383  ipt, psif, dpsifdx, testf, dtestfdx);
384 
385  // Call the pressure shape and test functions
386  this->pshape_nst(s, psip, testp);
387 
388  // Premultiply the weights and the Jacobian
389  double W = w * J;
390 
391  // Calculate local values of the pressure and velocity components
392  //--------------------------------------------------------------
393  double interpolated_p = 0.0;
394  Vector<double> interpolated_u(DIM, 0.0);
395  Vector<double> interpolated_x(DIM, 0.0);
396  Vector<double> mesh_veloc(DIM, 0.0);
397  Vector<double> dudt(DIM, 0.0);
398  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
399 
400  // Calculate pressure
401  for (unsigned l = 0; l < n_pres; l++)
402  {
403  interpolated_p += this->p_nst(l) * psip[l];
404  }
405 
406 
407  // Calculate velocities and derivatives
408 
409  // Loop over nodes
410  for (unsigned l = 0; l < n_node; l++)
411  {
412  // Loop over directions
413  for (unsigned i = 0; i < DIM; i++)
414  {
415  // Get the nodal value
416  double u_value = this->nodal_value(l, u_nodal_index[i]);
417  interpolated_u[i] += u_value * psif[l];
418  interpolated_x[i] += this->nodal_position(l, i) * psif[l];
419  dudt[i] += this->du_dt_nst(l, i) * psif[l];
420 
421  // Loop over derivative directions for velocity gradients
422  for (unsigned j = 0; j < DIM; j++)
423  {
424  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
425  }
426  }
427  }
428 
429  if (!ALE_is_disabled_flag)
430  {
431  // Loop over nodes
432  for (unsigned l = 0; l < n_node; l++)
433  {
434  // Loop over directions
435  for (unsigned i = 0; i < DIM; i++)
436  {
437  mesh_veloc[i] += this->dnodal_position_dt(l, i) * psif[l];
438  }
439  }
440  }
441 
442  // The strainrate used to compute the second invariant
443  DenseMatrix<double> strainrate_to_compute_second_invariant(DIM, DIM, 0.0);
444 
445  // the strainrate used to calculate the second invariant
446  // can be either the current one or the one extrapolated from
447  // previous velocity values
448  if (!this->Use_extrapolated_strainrate_to_compute_second_invariant)
449  {
450  this->strain_rate(s, strainrate_to_compute_second_invariant);
451  }
452  else
453  {
454  this->extrapolated_strain_rate(ipt,
455  strainrate_to_compute_second_invariant);
456  }
457 
458  // Calculate the second invariant
459  double second_invariant = SecondInvariantHelper::second_invariant(
460  strainrate_to_compute_second_invariant);
461 
462  // Get the viscosity according to the constitutive equation
463  double viscosity = this->Constitutive_eqn_pt->viscosity(second_invariant);
464 
465  // Get the user-defined body force terms
466  Vector<double> body_force(DIM);
467  this->get_body_force_nst(time, ipt, s, interpolated_x, body_force);
468 
469  // Get the user-defined source function
470  double source = this->get_source_nst(time, ipt, interpolated_x);
471 
472  // The generalised Newtonian viscosity differentiated w.r.t.
473  // the unknown velocities
474  DenseMatrix<double> dviscosity_dunknown(DIM, n_node, 0.0);
475 
476  if (!this->Use_extrapolated_strainrate_to_compute_second_invariant)
477  {
478  // Calculate the derivate of the viscosity w.r.t. the second invariant
479  double dviscosity_dsecond_invariant =
480  this->Constitutive_eqn_pt->dviscosity_dinvariant(second_invariant);
481 
482  // calculate the strainrate
483  DenseMatrix<double> strainrate(DIM, DIM, 0.0);
484  this->strain_rate(s, strainrate);
485 
486  // pre-compute the derivative of the second invariant w.r.t. the
487  // entries in the rate of strain tensor
488  DenseMatrix<double> dinvariant_dstrainrate(DIM, DIM, 0.0);
489 
490  // setting up a Kronecker Delta matrix, which has ones at the diagonal
491  // and zeros off-diagonally
492  DenseMatrix<double> kroneckerdelta(DIM, DIM, 0.0);
493 
494  for (unsigned i = 0; i < DIM; i++)
495  {
496  for (unsigned j = 0; j < DIM; j++)
497  {
498  if (i == j)
499  {
500  // Set the diagonal entries of the Kronecker delta
501  kroneckerdelta(i, i) = 1.0;
502 
503  double tmp = 0.0;
504  for (unsigned k = 0; k < DIM; k++)
505  {
506  if (k != i)
507  {
508  tmp += strainrate(k, k);
509  }
510  }
511  dinvariant_dstrainrate(i, i) = tmp;
512  }
513  else
514  {
515  dinvariant_dstrainrate(i, j) = -strainrate(j, i);
516  }
517  }
518  }
519 
520  // a rank four tensor storing the derivative of the strainrate
521  // w.r.t. the unknowns
522  RankFourTensor<double> dstrainrate_dunknown(DIM, DIM, DIM, n_node);
523 
524  // loop over the nodes
525  for (unsigned l = 0; l < n_node; l++)
526  {
527  // loop over the velocity components
528  for (unsigned k = 0; k < DIM; k++)
529  {
530  // loop over the entries of the rate of strain tensor
531  for (unsigned i = 0; i < DIM; i++)
532  {
533  for (unsigned j = 0; j < DIM; j++)
534  {
535  // initialise the entry to zero
536  dstrainrate_dunknown(i, j, k, l) = 0.0;
537 
538  // loop over velocities and directions
539  for (unsigned m = 0; m < DIM; m++)
540  {
541  for (unsigned n = 0; n < DIM; n++)
542  {
543  dstrainrate_dunknown(i, j, k, l) +=
544  0.5 *
545  (kroneckerdelta(i, m) * kroneckerdelta(j, n) +
546  kroneckerdelta(i, n) * kroneckerdelta(j, m)) *
547  kroneckerdelta(m, k) * dpsifdx(l, n);
548  }
549  }
550  }
551  }
552  }
553  }
554 
555  // Now calculate the derivative of the viscosity w.r.t. the unknowns
556  // loop over the nodes
557  for (unsigned l = 0; l < n_node; l++)
558  {
559  // loop over the velocities
560  for (unsigned k = 0; k < DIM; k++)
561  {
562  // loop over the entries in the rate of strain tensor
563  for (unsigned i = 0; i < DIM; i++)
564  {
565  for (unsigned j = 0; j < DIM; j++)
566  {
567  dviscosity_dunknown(k, l) += dviscosity_dsecond_invariant *
568  dinvariant_dstrainrate(i, j) *
569  dstrainrate_dunknown(i, j, k, l);
570  }
571  }
572  }
573  }
574  }
575 
576 
577  // MOMENTUM EQUATIONS
578  //==================
579 
580  // Number of master nodes and storage for the weight of the shape function
581  unsigned n_master = 1;
582  double hang_weight = 1.0;
583 
584  // Loop over the nodes for the test functions/equations
585  //----------------------------------------------------
586  for (unsigned l = 0; l < n_node; l++)
587  {
588  // Local boolean to indicate whether the node is hanging
589  bool is_node_hanging = node_pt(l)->is_hanging();
590 
591  // If the node is hanging
592  if (is_node_hanging)
593  {
594  hang_info_pt = node_pt(l)->hanging_pt();
595 
596  // Read out number of master nodes from hanging data
597  n_master = hang_info_pt->nmaster();
598  }
599  // Otherwise the node is its own master
600  else
601  {
602  n_master = 1;
603  }
604 
605  // Loop over the master nodes
606  for (unsigned m = 0; m < n_master; m++)
607  {
608  // Loop over velocity components for equations
609  for (unsigned i = 0; i < DIM; i++)
610  {
611  // Get the equation number
612  // If the node is hanging
613  if (is_node_hanging)
614  {
615  // Get the equation number from the master node
616  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
617  u_nodal_index[i]);
618  // Get the hang weight from the master node
619  hang_weight = hang_info_pt->master_weight(m);
620  }
621  // If the node is not hanging
622  else
623  {
624  // Local equation number
625  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
626 
627  // Node contributes with full weight
628  hang_weight = 1.0;
629  }
630 
631  // If it's not a boundary condition...
632  if (local_eqn >= 0)
633  {
634  // Temporary variable to hold the residuals
635  double sum = 0.0;
636 
637  // Add the user-defined body force terms
638  sum += body_force[i];
639 
640  // Add the gravitational body force term
641  sum += scaled_re_inv_fr * G[i];
642 
643  // Add in the inertial term
644  sum -= scaled_re_st * dudt[i];
645 
646  // Convective terms, including mesh velocity
647  for (unsigned k = 0; k < DIM; k++)
648  {
649  double tmp = scaled_re * interpolated_u[k];
650  if (!ALE_is_disabled_flag)
651  {
652  tmp -= scaled_re_st * mesh_veloc[k];
653  }
654  sum -= tmp * interpolated_dudx(i, k);
655  }
656 
657  // Add the pressure gradient term
658  // Potentially pre-multiply by viscosity ratio
659  if (this->Pre_multiply_by_viscosity_ratio)
660  {
661  sum = visc_ratio * viscosity *
662  (sum * testf[l] + interpolated_p * dtestfdx(l, i)) * W *
663  hang_weight;
664  }
665  else
666  {
667  sum = (sum * testf[l] + interpolated_p * dtestfdx(l, i)) * W *
668  hang_weight;
669  }
670 
671  // Add in the stress tensor terms
672  // The viscosity ratio needs to go in here to ensure
673  // continuity of normal stress is satisfied even in flows
674  // with zero pressure gradient!
675  for (unsigned k = 0; k < DIM; k++)
676  {
677  sum -= visc_ratio * viscosity *
678  (interpolated_dudx(i, k) +
679  this->Gamma[i] * interpolated_dudx(k, i)) *
680  dtestfdx(l, k) * W * hang_weight;
681  }
682 
683  // Add contribution
684  residuals[local_eqn] += sum;
685 
686  // CALCULATE THE JACOBIAN
687  if (flag)
688  {
689  // Number of master nodes and weights
690  unsigned n_master2 = 1;
691  double hang_weight2 = 1.0;
692  // Loop over the velocity nodes for columns
693  for (unsigned l2 = 0; l2 < n_node; l2++)
694  {
695  // Local boolean to indicate whether the node is hanging
696  bool is_node2_hanging = node_pt(l2)->is_hanging();
697 
698  // If the node is hanging
699  if (is_node2_hanging)
700  {
701  hang_info2_pt = node_pt(l2)->hanging_pt();
702  // Read out number of master nodes from hanging data
703  n_master2 = hang_info2_pt->nmaster();
704  }
705  // Otherwise the node is its own master
706  else
707  {
708  n_master2 = 1;
709  }
710 
711  // Loop over the master nodes
712  for (unsigned m2 = 0; m2 < n_master2; m2++)
713  {
714  // Loop over the velocity components
715  for (unsigned i2 = 0; i2 < DIM; i2++)
716  {
717  // Get the number of the unknown
718  // If the node is hanging
719  if (is_node2_hanging)
720  {
721  // Get the equation number from the master node
722  local_unknown = this->local_hang_eqn(
723  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
724  // Get the hang weights
725  hang_weight2 = hang_info2_pt->master_weight(m2);
726  }
727  else
728  {
729  local_unknown =
730  this->nodal_local_eqn(l2, u_nodal_index[i2]);
731  hang_weight2 = 1.0;
732  }
733 
734  // If the unknown is non-zero
735  if (local_unknown >= 0)
736  {
737  // Add contribution to Elemental Matrix
738  jacobian(local_eqn, local_unknown) -=
739  visc_ratio * viscosity * this->Gamma[i] *
740  dpsifdx(l2, i) * dtestfdx(l, i2) * W * hang_weight *
741  hang_weight2;
742 
743  // Now add in the inertial terms
744  jacobian(local_eqn, local_unknown) -=
745  scaled_re * psif[l2] * interpolated_dudx(i, i2) *
746  testf[l] * W * hang_weight * hang_weight2;
747 
748  if (
749  !this
750  ->Use_extrapolated_strainrate_to_compute_second_invariant)
751  {
752  for (unsigned k = 0; k < DIM; k++)
753  {
754  jacobian(local_eqn, local_unknown) -=
755  visc_ratio * dviscosity_dunknown(i2, l2) *
756  (interpolated_dudx(i, k) +
757  this->Gamma[i] * interpolated_dudx(k, i)) *
758  dtestfdx(l, k) * W * hang_weight * hang_weight2;
759  }
760  }
761 
762  // Extra diagonal components if i2=i
763  if (i2 == i)
764  {
765  // Mass matrix entries
766  // Again note the positive sign because the mass
767  // matrix is taken on the other side of the equation
768  if (flag == 2)
769  {
770  mass_matrix(local_eqn, local_unknown) +=
771  scaled_re_st * psif[l2] * testf[l] * W *
772  hang_weight * hang_weight2;
773  }
774 
775  // du/dt term
776  jacobian(local_eqn, local_unknown) -=
777  scaled_re_st *
778  node_pt(l2)->time_stepper_pt()->weight(1, 0) *
779  psif[l2] * testf[l] * W * hang_weight *
780  hang_weight2;
781 
782  // Extra advective terms
783  for (unsigned k = 0; k < DIM; k++)
784  {
785  double tmp = scaled_re * interpolated_u[k];
786  if (!ALE_is_disabled_flag)
787  {
788  tmp -= scaled_re_st * mesh_veloc[k];
789  }
790 
791  jacobian(local_eqn, local_unknown) -=
792  tmp * dpsifdx(l2, k) * testf[l] * W *
793  hang_weight * hang_weight2;
794  }
795 
796  // Extra viscous terms
797  for (unsigned k = 0; k < DIM; k++)
798  {
799  jacobian(local_eqn, local_unknown) -=
800  visc_ratio * viscosity * dpsifdx(l2, k) *
801  dtestfdx(l, k) * W * hang_weight * hang_weight2;
802  }
803  }
804  }
805  }
806  }
807  }
808 
809  // Loop over the pressure shape functions
810  for (unsigned l2 = 0; l2 < n_pres; l2++)
811  {
812  // If the pressure dof is hanging
813  if (pressure_dof_is_hanging[l2])
814  {
815  hang_info2_pt =
816  this->pressure_node_pt(l2)->hanging_pt(p_index);
817  // Pressure dof is hanging so it must be nodal-based
818  // Get the number of master nodes from the pressure node
819  n_master2 = hang_info2_pt->nmaster();
820  }
821  // Otherwise the node is its own master
822  else
823  {
824  n_master2 = 1;
825  }
826 
827  // Loop over the master nodes
828  for (unsigned m2 = 0; m2 < n_master2; m2++)
829  {
830  // Get the number of the unknown
831  // If the pressure dof is hanging
832  if (pressure_dof_is_hanging[l2])
833  {
834  // Get the unknown from the master node
835  local_unknown = this->local_hang_eqn(
836  hang_info2_pt->master_node_pt(m2), p_index);
837  // Get the weight from the hanging object
838  hang_weight2 = hang_info2_pt->master_weight(m2);
839  }
840  else
841  {
842  local_unknown = this->p_local_eqn(l2);
843  hang_weight2 = 1.0;
844  }
845 
846  // If the unknown is not pinned
847  if (local_unknown >= 0)
848  {
849  if (this->Pre_multiply_by_viscosity_ratio)
850  {
851  jacobian(local_eqn, local_unknown) +=
852  visc_ratio * viscosity * psip[l2] * dtestfdx(l, i) *
853  W * hang_weight * hang_weight2;
854  }
855  else
856  {
857  jacobian(local_eqn, local_unknown) +=
858  psip[l2] * dtestfdx(l, i) * W * hang_weight *
859  hang_weight2;
860  }
861  }
862  }
863  }
864 
865  } // End of Jacobian calculation
866 
867  } // End of if not boundary condition statement
868 
869  } // End of loop over components of non-hanging node
870 
871  } // End of loop over master nodes
872 
873  } // End of loop over nodes for equations
874 
875 
876  // CONTINUITY EQUATION
877  //===================
878 
879  // Loop over the pressure shape functions
880  for (unsigned l = 0; l < n_pres; l++)
881  {
882  // If the pressure dof is hanging
883  if (pressure_dof_is_hanging[l])
884  {
885  // Pressure dof is hanging so it must be nodal-based
886  // Get the hang info object
887  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
888 
889  // Get the number of master nodes from the pressure node
890  n_master = hang_info_pt->nmaster();
891  }
892  // Otherwise the node is its own master
893  else
894  {
895  n_master = 1;
896  }
897 
898  // Loop over the master nodes
899  for (unsigned m = 0; m < n_master; m++)
900  {
901  // Get the number of the unknown
902  // If the pressure dof is hanging
903  if (pressure_dof_is_hanging[l])
904  {
905  // Get the local equation from the master node
906  local_eqn =
907  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
908  // Get the weight for the node
909  hang_weight = hang_info_pt->master_weight(m);
910  }
911  else
912  {
913  local_eqn = this->p_local_eqn(l);
914  hang_weight = 1.0;
915  }
916 
917  // If the equation is not pinned
918  if (local_eqn >= 0)
919  {
920  // Source term
921  residuals[local_eqn] -= source * testp[l] * W * hang_weight;
922 
923  // Loop over velocity components
924  for (unsigned k = 0; k < DIM; k++)
925  {
926  // Potentially pre-multiply by viscosity ratio
927  if (this->Pre_multiply_by_viscosity_ratio)
928  {
929  residuals[local_eqn] += visc_ratio * viscosity *
930  interpolated_dudx(k, k) * testp[l] * W *
931  hang_weight;
932  }
933  else
934  {
935  residuals[local_eqn] +=
936  interpolated_dudx(k, k) * testp[l] * W * hang_weight;
937  }
938  }
939 
940  // CALCULATE THE JACOBIAN
941  if (flag)
942  {
943  unsigned n_master2 = 1;
944  double hang_weight2 = 1.0;
945  // Loop over the velocity nodes for columns
946  for (unsigned l2 = 0; l2 < n_node; l2++)
947  {
948  // Local boolean to indicate whether the node is hanging
949  bool is_node2_hanging = node_pt(l2)->is_hanging();
950 
951  // If the node is hanging
952  if (is_node2_hanging)
953  {
954  hang_info2_pt = node_pt(l2)->hanging_pt();
955  // Read out number of master nodes from hanging data
956  n_master2 = hang_info2_pt->nmaster();
957  }
958  // Otherwise the node is its own master
959  else
960  {
961  n_master2 = 1;
962  }
963 
964  // Loop over the master nodes
965  for (unsigned m2 = 0; m2 < n_master2; m2++)
966  {
967  // Loop over the velocity components
968  for (unsigned i2 = 0; i2 < DIM; i2++)
969  {
970  // Get the number of the unknown
971  // If the node is hanging
972  if (is_node2_hanging)
973  {
974  // Get the equation number from the master node
975  local_unknown = this->local_hang_eqn(
976  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
977  hang_weight2 = hang_info2_pt->master_weight(m2);
978  }
979  else
980  {
981  local_unknown =
982  this->nodal_local_eqn(l2, u_nodal_index[i2]);
983  hang_weight2 = 1.0;
984  }
985 
986  // If the unknown is not pinned
987  if (local_unknown >= 0)
988  {
989  if (this->Pre_multiply_by_viscosity_ratio)
990  {
991  jacobian(local_eqn, local_unknown) +=
992  visc_ratio * viscosity * dpsifdx(l2, i2) * testp[l] *
993  W * hang_weight * hang_weight2;
994  }
995  else
996  {
997  jacobian(local_eqn, local_unknown) +=
998  dpsifdx(l2, i2) * testp[l] * W * hang_weight *
999  hang_weight2;
1000  }
1001  }
1002  }
1003  }
1004  }
1005 
1006  // NO PRESSURE CONTRIBUTION TO THE JACOBIAN
1007 
1008  } // End of jacobian calculation
1009  }
1010  }
1011  } // End of loop over pressure variables
1012 
1013  } // End of loop over integration points
1014  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
G
JacobiRotation< float > G
Definition: Jacobi_makeGivens.cpp:2
J
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
m2
MatrixType m2(n_dims)
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::FiniteElement::nodal_value
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
oomph::FiniteElement::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
oomph::FiniteElement::nodal_position
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317
oomph::FiniteElement::dnodal_position_dt
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2333
oomph::GeneralisedNewtonianNavierStokesEquations::get_source_nst
virtual double get_source_nst(const double &time, const unsigned &ipt, const Vector< double > &x)
Definition: generalised_newtonian_navier_stokes_elements.h:318
oomph::GeneralisedNewtonianNavierStokesEquations::Pre_multiply_by_viscosity_ratio
static bool Pre_multiply_by_viscosity_ratio
Definition: generalised_newtonian_navier_stokes_elements.h:159
oomph::GeneralisedNewtonianNavierStokesEquations::Use_extrapolated_strainrate_to_compute_second_invariant
bool Use_extrapolated_strainrate_to_compute_second_invariant
Definition: generalised_newtonian_navier_stokes_elements.h:197
oomph::GeneralisedNewtonianNavierStokesEquations::ALE_is_disabled
bool ALE_is_disabled
Definition: generalised_newtonian_navier_stokes_elements.h:202
oomph::GeneralisedNewtonianNavierStokesEquations::get_body_force_nst
virtual void get_body_force_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Definition: generalised_newtonian_navier_stokes_elements.h:250
oomph::GeneralisedNewtonianNavierStokesEquations::Constitutive_eqn_pt
GeneralisedNewtonianConstitutiveEquation< DIM > * Constitutive_eqn_pt
Pointer to the generalised Newtonian constitutive equation.
Definition: generalised_newtonian_navier_stokes_elements.h:193
oomph::GeneralisedNewtonianNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_nst
virtual double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
oomph::GeneralisedNewtonianNavierStokesEquations::strain_rate
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: generalised_newtonian_navier_stokes_elements.cc:1012
oomph::GeneralisedNewtonianNavierStokesEquations::viscosity_ratio
const double & viscosity_ratio() const
Definition: generalised_newtonian_navier_stokes_elements.h:455
oomph::GeneralisedNewtonianNavierStokesEquations::g
const Vector< double > & g() const
Vector of gravitational components.
Definition: generalised_newtonian_navier_stokes_elements.h:492
oomph::GeneralisedNewtonianNavierStokesEquations::extrapolated_strain_rate
virtual void extrapolated_strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Definition: generalised_newtonian_navier_stokes_elements.cc:1175
oomph::GeneralisedNewtonianNavierStokesEquations::du_dt_nst
double du_dt_nst(const unsigned &n, const unsigned &i) const
Definition: generalised_newtonian_navier_stokes_elements.h:594
oomph::GeneralisedNewtonianNavierStokesEquations::re_invfr
const double & re_invfr() const
Global inverse Froude number.
Definition: generalised_newtonian_navier_stokes_elements.h:480
oomph::GeneralisedNewtonianNavierStokesEquations::pshape_nst
virtual void pshape_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
oomph::GeneralisedNewtonianNavierStokesEquations::p_nodal_index_nst
virtual int p_nodal_index_nst() const
Definition: generalised_newtonian_navier_stokes_elements.h:650
oomph::GeneralisedNewtonianNavierStokesEquations::Gamma
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition: generalised_newtonian_navier_stokes_elements.h:425
oomph::GeneralisedNewtonianNavierStokesEquations::npres_nst
virtual unsigned npres_nst() const =0
Function to return number of pressure degrees of freedom.
oomph::GeneralisedNewtonianNavierStokesEquations::density_ratio
const double & density_ratio() const
Definition: generalised_newtonian_navier_stokes_elements.h:468
oomph::GeneralisedNewtonianNavierStokesEquations::re_st
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: generalised_newtonian_navier_stokes_elements.h:436
oomph::GeneralisedNewtonianNavierStokesEquations::p_nst
virtual double p_nst(const unsigned &n_p) const =0
oomph::GeneralisedNewtonianNavierStokesEquations::re
const double & re() const
Reynolds number.
Definition: generalised_newtonian_navier_stokes_elements.h:430
oomph::HangInfo::nmaster
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::RefineableElement::local_hang_eqn
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Definition: refineable_elements.h:278
oomph::RefineableGeneralisedNewtonianNavierStokesEquations::pressure_node_pt
virtual Node * pressure_node_pt(const unsigned &n_p)
Definition: generalised_newtonian_refineable_navier_stokes_elements.h:119
oomph::TimeStepper::weight
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
oomph::TimeStepper::time_pt
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
oomph::Time::time
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
s
RealScalar s
Definition: level1_cplx_impl.h:130
k
char char char int int * k
Definition: level2_impl.h:374
tmp
Eigen::Matrix< Scalar, Dynamic, Dynamic, ColMajor > tmp
Definition: level3_impl.h:365
DIM
#define DIM
Definition: linearised_navier_stokes_elements.h:44
Global_Parameters::body_force
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Definition: axisym_linear_elasticity/cylinder/cylinder.cc:96
TestProblem::source
void source(const Vector< double > &x, Vector< double > &f)
Source function.
Definition: unstructured_two_d_circle.cc:46
oomph::QuadTreeNames::W
@ W
Definition: quadtree.h:63
oomph::SecondInvariantHelper::second_invariant
double second_invariant(const DenseMatrix< double > &tensor)
Compute second invariant of tensor.
Definition: oomph_utilities.cc:163
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::ALE_is_disabled, Global_Parameters::body_force(), DIM, extrapolated_strain_rate(), G, GlobalPhysicalVariables::Gamma, i, J, j, k, m, m2(), oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), n, oomph::HangInfo::nmaster(), s, oomph::SecondInvariantHelper::second_invariant(), TestProblem::source(), tmp, w, and oomph::QuadTreeNames::W.

◆ further_build()

template<unsigned DIM>
void oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::further_build ( )
inlinevirtual

Further build, pass the pointers down to the sons.

Reimplemented from oomph::RefineableElement.

Reimplemented in oomph::RefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >, oomph::RefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >, oomph::PRefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >, oomph::RefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >, oomph::PRefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >, and oomph::PRefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >.

186  {
187  // Find the father element
188  RefineableGeneralisedNewtonianNavierStokesEquations<
189  DIM>* cast_father_element_pt =
190  dynamic_cast<RefineableGeneralisedNewtonianNavierStokesEquations<DIM>*>(
191  this->father_element_pt());
192 
193  // Set the viscosity ratio pointer
194  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
195  // Set the density ratio pointer
196  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
197  // Set pointer to global Reynolds number
198  this->Re_pt = cast_father_element_pt->re_pt();
199  // Set pointer to global Reynolds number x Strouhal number (=Womersley)
200  this->ReSt_pt = cast_father_element_pt->re_st_pt();
201  // Set pointer to global Reynolds number x inverse Froude number
202  this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
203  // Set pointer to global gravity Vector
204  this->G_pt = cast_father_element_pt->g_pt();
205 
206  // Set pointer to body force function
207  this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
208 
209  // Set pointer to volumetric source function
210  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
211 
212  // Set the ALE flag
213  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
214  }
oomph::GeneralisedNewtonianNavierStokesEquations::ReInvFr_pt
double * ReInvFr_pt
Definition: generalised_newtonian_navier_stokes_elements.h:181
oomph::GeneralisedNewtonianNavierStokesEquations::Viscosity_Ratio_pt
double * Viscosity_Ratio_pt
Definition: generalised_newtonian_navier_stokes_elements.h:165
oomph::GeneralisedNewtonianNavierStokesEquations::Density_Ratio_pt
double * Density_Ratio_pt
Definition: generalised_newtonian_navier_stokes_elements.h:169
oomph::GeneralisedNewtonianNavierStokesEquations::Body_force_fct_pt
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
Definition: generalised_newtonian_navier_stokes_elements.h:187
oomph::GeneralisedNewtonianNavierStokesEquations::Source_fct_pt
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
Definition: generalised_newtonian_navier_stokes_elements.h:190
oomph::GeneralisedNewtonianNavierStokesEquations::Re_pt
double * Re_pt
Pointer to global Reynolds number.
Definition: generalised_newtonian_navier_stokes_elements.h:174
oomph::GeneralisedNewtonianNavierStokesEquations::ReSt_pt
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: generalised_newtonian_navier_stokes_elements.h:177
oomph::GeneralisedNewtonianNavierStokesEquations::G_pt
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: generalised_newtonian_navier_stokes_elements.h:184
oomph::RefineableElement::father_element_pt
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Definition: refineable_elements.h:539
oomph::RefineableGeneralisedNewtonianNavierStokesEquations::RefineableGeneralisedNewtonianNavierStokesEquations
RefineableGeneralisedNewtonianNavierStokesEquations()
Constructor.
Definition: generalised_newtonian_refineable_navier_stokes_elements.h:71

References oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::ALE_is_disabled, oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::Body_force_fct_pt, oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::Density_Ratio_pt, DIM, oomph::RefineableElement::father_element_pt(), oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::G_pt, oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::Re_pt, oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::ReInvFr_pt, oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::ReSt_pt, oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::Source_fct_pt, and oomph::GeneralisedNewtonianNavierStokesEquations< DIM >::Viscosity_Ratio_pt.

Referenced by oomph::PRefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >::further_build(), and oomph::RefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >::further_build().

◆ get_dresidual_dnodal_coordinates()

template<unsigned DIM>
void oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::get_dresidual_dnodal_coordinates ( RankThreeTensor< double > &  dresidual_dnodal_coordinates)
protectedvirtual

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::GeneralisedNewtonianNavierStokesEquations< DIM >.

1027  {
1028  // Return immediately if there are no dofs
1029  if (ndof() == 0)
1030  {
1031  return;
1032  }
1033 
1034  // Determine number of nodes in element
1035  const unsigned n_node = nnode();
1036 
1037  // Get continuous time from timestepper of first node
1038  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1039 
1040  // Determine number of pressure dofs in element
1041  const unsigned n_pres = this->npres_nst();
1042 
1043  // Find the indices at which the local velocities are stored
1044  unsigned u_nodal_index[DIM];
1045  for (unsigned i = 0; i < DIM; i++)
1046  {
1047  u_nodal_index[i] = this->u_index_nst(i);
1048  }
1049 
1050  // Which nodal value represents the pressure? (Negative if pressure
1051  // is not based on nodal interpolation).
1052  const int p_index = this->p_nodal_index_nst();
1053 
1054  // Local array of booleans that are true if the l-th pressure value is
1055  // hanging (avoid repeated virtual function calls)
1056  bool pressure_dof_is_hanging[n_pres];
1057 
1058  // If the pressure is stored at a node
1059  if (p_index >= 0)
1060  {
1061  // Read out whether the pressure is hanging
1062  for (unsigned l = 0; l < n_pres; ++l)
1063  {
1064  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
1065  }
1066  }
1067  // Otherwise the pressure is not stored at a node and so cannot hang
1068  else
1069  {
1070  for (unsigned l = 0; l < n_pres; ++l)
1071  {
1072  pressure_dof_is_hanging[l] = false;
1073  }
1074  }
1075 
1076  // Set up memory for the shape and test functions
1077  Shape psif(n_node), testf(n_node);
1078  DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
1079 
1080  // Set up memory for pressure shape and test functions
1081  Shape psip(n_pres), testp(n_pres);
1082 
1083  // Determine number of shape controlling nodes
1084  const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1085 
1086  // Deriatives of shape fct derivatives w.r.t. nodal coords
1087  RankFourTensor<double> d_dpsifdx_dX(
1088  DIM, n_shape_controlling_node, n_node, DIM);
1089  RankFourTensor<double> d_dtestfdx_dX(
1090  DIM, n_shape_controlling_node, n_node, DIM);
1091 
1092  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1093  DenseMatrix<double> dJ_dX(DIM, n_shape_controlling_node);
1094 
1095  // Derivatives of derivative of u w.r.t. nodal coords
1096  RankFourTensor<double> d_dudx_dX(DIM, n_shape_controlling_node, DIM, DIM);
1097 
1098  // Derivatives of nodal velocities w.r.t. nodal coords:
1099  // Assumption: Interaction only local via no-slip so
1100  // X_ij only affects U_ij.
1101  DenseMatrix<double> d_U_dX(DIM, n_shape_controlling_node, 0.0);
1102 
1103  // Determine the number of integration points
1104  const unsigned n_intpt = integral_pt()->nweight();
1105 
1106  // Vector to hold local coordinates
1107  Vector<double> s(DIM);
1108 
1109  // Get physical variables from element
1110  // (Reynolds number must be multiplied by the density ratio)
1111  double scaled_re = this->re() * this->density_ratio();
1112  double scaled_re_st = this->re_st() * this->density_ratio();
1113  double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
1114  double visc_ratio = this->viscosity_ratio();
1115  Vector<double> G = this->g();
1116 
1117  // FD step
1118  double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
1119 
1120  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1121  // Assumption: Interaction only local via no-slip so
1122  // X_ij only affects U_ij.
1123  bool element_has_node_with_aux_node_update_fct = false;
1124 
1125  std::map<Node*, unsigned> local_shape_controlling_node_lookup =
1126  shape_controlling_node_lookup();
1127 
1128  // FD loop over shape-controlling nodes
1129  for (std::map<Node*, unsigned>::iterator it =
1130  local_shape_controlling_node_lookup.begin();
1131  it != local_shape_controlling_node_lookup.end();
1132  it++)
1133  {
1134  // Get node
1135  Node* nod_pt = it->first;
1136 
1137  // Get its number
1138  unsigned q = it->second;
1139 
1140  // Only compute if there's a node-update fct involved
1141  if (nod_pt->has_auxiliary_node_update_fct_pt())
1142  {
1143  element_has_node_with_aux_node_update_fct = true;
1144 
1145  // Current nodal velocity
1146  Vector<double> u_ref(DIM);
1147  for (unsigned i = 0; i < DIM; i++)
1148  {
1149  u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
1150  }
1151 
1152  // FD
1153  for (unsigned p = 0; p < DIM; p++)
1154  {
1155  // Make backup
1156  double backup = nod_pt->x(p);
1157 
1158  // Do FD step. No node update required as we're
1159  // attacking the coordinate directly...
1160  nod_pt->x(p) += eps_fd;
1161 
1162  // Do auxiliary node update (to apply no slip)
1163  nod_pt->perform_auxiliary_node_update_fct();
1164 
1165  // Evaluate
1166  d_U_dX(p, q) =
1167  (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
1168 
1169  // Reset
1170  nod_pt->x(p) = backup;
1171 
1172  // Do auxiliary node update (to apply no slip)
1173  nod_pt->perform_auxiliary_node_update_fct();
1174  }
1175  }
1176  }
1177 
1178  // Integer to store the local equation number
1179  int local_eqn = 0;
1180 
1181  // Pointers to hang info object
1182  HangInfo* hang_info_pt = 0;
1183 
1184  // Loop over the integration points
1185  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1186  {
1187  // Assign values of s
1188  for (unsigned i = 0; i < DIM; i++)
1189  {
1190  s[i] = integral_pt()->knot(ipt, i);
1191  }
1192 
1193  // Get the integral weight
1194  const double w = integral_pt()->weight(ipt);
1195 
1196  // Call the derivatives of the shape and test functions
1197  const double J =
1198  this->dshape_and_dtest_eulerian_at_knot_nst(ipt,
1199  psif,
1200  dpsifdx,
1201  d_dpsifdx_dX,
1202  testf,
1203  dtestfdx,
1204  d_dtestfdx_dX,
1205  dJ_dX);
1206 
1207  // Call the pressure shape and test functions
1208  this->pshape_nst(s, psip, testp);
1209 
1210  // Calculate local values of the pressure and velocity components
1211  // Allocate
1212  double interpolated_p = 0.0;
1213  Vector<double> interpolated_u(DIM, 0.0);
1214  Vector<double> interpolated_x(DIM, 0.0);
1215  Vector<double> mesh_velocity(DIM, 0.0);
1216  Vector<double> dudt(DIM, 0.0);
1217  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1218 
1219  // Calculate pressure
1220  for (unsigned l = 0; l < n_pres; l++)
1221  {
1222  interpolated_p += this->p_nst(l) * psip[l];
1223  }
1224 
1225  // Calculate velocities and derivatives:
1226 
1227  // Loop over nodes
1228  for (unsigned l = 0; l < n_node; l++)
1229  {
1230  // Loop over directions
1231  for (unsigned i = 0; i < DIM; i++)
1232  {
1233  // Get the nodal value
1234  const double u_value = nodal_value(l, u_nodal_index[i]);
1235  interpolated_u[i] += u_value * psif[l];
1236  interpolated_x[i] += nodal_position(l, i) * psif[l];
1237  dudt[i] += this->du_dt_nst(l, i) * psif[l];
1238 
1239  // Loop over derivative directions
1240  for (unsigned j = 0; j < DIM; j++)
1241  {
1242  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1243  }
1244  }
1245  }
1246 
1247  if (!this->ALE_is_disabled)
1248  {
1249  // Loop over nodes
1250  for (unsigned l = 0; l < n_node; l++)
1251  {
1252  // Loop over directions
1253  for (unsigned i = 0; i < DIM; i++)
1254  {
1255  mesh_velocity[i] += this->dnodal_position_dt(l, i) * psif[l];
1256  }
1257  }
1258  }
1259 
1260  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1261 
1262  // Loop over shape-controlling nodes
1263  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1264  {
1265  // Loop over coordinate directions
1266  for (unsigned p = 0; p < DIM; p++)
1267  {
1268  for (unsigned i = 0; i < DIM; i++)
1269  {
1270  for (unsigned k = 0; k < DIM; k++)
1271  {
1272  double aux = 0.0;
1273  for (unsigned j = 0; j < n_node; j++)
1274  {
1275  aux +=
1276  nodal_value(j, u_nodal_index[i]) * d_dpsifdx_dX(p, q, j, k);
1277  }
1278  d_dudx_dX(p, q, i, k) = aux;
1279  }
1280  }
1281  }
1282  }
1283 
1284  // Get weight of actual nodal position/value in computation of mesh
1285  // velocity from positional/value time stepper
1286  const double pos_time_weight =
1287  node_pt(0)->position_time_stepper_pt()->weight(1, 0);
1288  const double val_time_weight =
1289  node_pt(0)->time_stepper_pt()->weight(1, 0);
1290 
1291  // Get the user-defined body force terms
1292  Vector<double> body_force(DIM);
1293  this->get_body_force_nst(time, ipt, s, interpolated_x, body_force);
1294 
1295  // Get the user-defined source function
1296  const double source = this->get_source_nst(time, ipt, interpolated_x);
1297 
1298  // Get gradient of body force function
1299  DenseMatrix<double> d_body_force_dx(DIM, DIM, 0.0);
1300  this->get_body_force_gradient_nst(
1301  time, ipt, s, interpolated_x, d_body_force_dx);
1302 
1303  // Get gradient of source function
1304  Vector<double> source_gradient(DIM, 0.0);
1305  this->get_source_gradient_nst(time, ipt, interpolated_x, source_gradient);
1306 
1307 
1308  // Assemble shape derivatives
1309  //---------------------------
1310 
1311  // MOMENTUM EQUATIONS
1312  // ------------------
1313 
1314  // Number of master nodes and storage for the weight of the shape function
1315  unsigned n_master = 1;
1316  double hang_weight = 1.0;
1317 
1318  // Loop over the test functions
1319  for (unsigned l = 0; l < n_node; l++)
1320  {
1321  // Local boolean to indicate whether the node is hanging
1322  bool is_node_hanging = node_pt(l)->is_hanging();
1323 
1324  // If the node is hanging
1325  if (is_node_hanging)
1326  {
1327  hang_info_pt = node_pt(l)->hanging_pt();
1328 
1329  // Read out number of master nodes from hanging data
1330  n_master = hang_info_pt->nmaster();
1331  }
1332  // Otherwise the node is its own master
1333  else
1334  {
1335  n_master = 1;
1336  }
1337 
1338  // Loop over the master nodes
1339  for (unsigned m = 0; m < n_master; m++)
1340  {
1341  // Loop over coordinate directions
1342  for (unsigned i = 0; i < DIM; i++)
1343  {
1344  // Get the equation number
1345  // If the node is hanging
1346  if (is_node_hanging)
1347  {
1348  // Get the equation number from the master node
1349  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1350  u_nodal_index[i]);
1351  // Get the hang weight from the master node
1352  hang_weight = hang_info_pt->master_weight(m);
1353  }
1354  // If the node is not hanging
1355  else
1356  {
1357  // Local equation number
1358  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
1359 
1360  // Node contributes with full weight
1361  hang_weight = 1.0;
1362  }
1363 
1364  // IF it's not a boundary condition
1365  if (local_eqn >= 0)
1366  {
1367  // Loop over coordinate directions
1368  for (unsigned p = 0; p < DIM; p++)
1369  {
1370  // Loop over shape controlling nodes
1371  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1372  {
1373  // Residual x deriv of Jacobian
1374  // ----------------------------
1375 
1376  // Add the user-defined body force terms
1377  double sum = body_force[i] * testf[l];
1378 
1379  // Add the gravitational body force term
1380  sum += scaled_re_inv_fr * testf[l] * G[i];
1381 
1382  // Add the pressure gradient term
1383  sum += interpolated_p * dtestfdx(l, i);
1384 
1385  // Add in the stress tensor terms
1386  // The viscosity ratio needs to go in here to ensure
1387  // continuity of normal stress is satisfied even in flows
1388  // with zero pressure gradient!
1389  for (unsigned k = 0; k < DIM; k++)
1390  {
1391  sum -= visc_ratio *
1392  (interpolated_dudx(i, k) +
1393  this->Gamma[i] * interpolated_dudx(k, i)) *
1394  dtestfdx(l, k);
1395  }
1396 
1397  // Add in the inertial terms
1398 
1399  // du/dt term
1400  sum -= scaled_re_st * dudt[i] * testf[l];
1401 
1402  // Convective terms, including mesh velocity
1403  for (unsigned k = 0; k < DIM; k++)
1404  {
1405  double tmp = scaled_re * interpolated_u[k];
1406  if (!this->ALE_is_disabled)
1407  {
1408  tmp -= scaled_re_st * mesh_velocity[k];
1409  }
1410  sum -= tmp * interpolated_dudx(i, k) * testf[l];
1411  }
1412 
1413  // Multiply throsugh by deriv of Jacobian and integration
1414  // weight
1415  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1416  sum * dJ_dX(p, q) * w * hang_weight;
1417 
1418  // Derivative of residual x Jacobian
1419  // ---------------------------------
1420 
1421  // Body force
1422  sum = d_body_force_dx(i, p) * psif(q) * testf(l);
1423 
1424  // Pressure gradient term
1425  sum += interpolated_p * d_dtestfdx_dX(p, q, l, i);
1426 
1427  // Viscous term
1428  for (unsigned k = 0; k < DIM; k++)
1429  {
1430  sum -=
1431  visc_ratio * ((interpolated_dudx(i, k) +
1432  this->Gamma[i] * interpolated_dudx(k, i)) *
1433  d_dtestfdx_dX(p, q, l, k) +
1434  (d_dudx_dX(p, q, i, k) +
1435  this->Gamma[i] * d_dudx_dX(p, q, k, i)) *
1436  dtestfdx(l, k));
1437  }
1438 
1439  // Convective terms, including mesh velocity
1440  for (unsigned k = 0; k < DIM; k++)
1441  {
1442  double tmp = scaled_re * interpolated_u[k];
1443  if (!this->ALE_is_disabled)
1444  {
1445  tmp -= scaled_re_st * mesh_velocity[k];
1446  }
1447  sum -= tmp * d_dudx_dX(p, q, i, k) * testf(l);
1448  }
1449  if (!this->ALE_is_disabled)
1450  {
1451  sum += scaled_re_st * pos_time_weight * psif(q) *
1452  interpolated_dudx(i, p) * testf(l);
1453  }
1454 
1455  // Multiply through by Jacobian and integration weight
1456  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1457  sum * J * w * hang_weight;
1458 
1459  } // End of loop over shape controlling nodes q
1460  } // End of loop over coordinate directions p
1461 
1462 
1463  // Derivs w.r.t. to nodal velocities
1464  // ---------------------------------
1465  if (element_has_node_with_aux_node_update_fct)
1466  {
1467  // Loop over local nodes
1468  for (unsigned q_local = 0; q_local < n_node; q_local++)
1469  {
1470  // Number of master nodes and storage for the weight of
1471  // the shape function
1472  unsigned n_master2 = 1;
1473  double hang_weight2 = 1.0;
1474  HangInfo* hang_info2_pt = 0;
1475 
1476  // Local boolean to indicate whether the node is hanging
1477  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1478 
1479  Node* actual_shape_controlling_node_pt = node_pt(q_local);
1480 
1481  // If the node is hanging
1482  if (is_node_hanging2)
1483  {
1484  hang_info2_pt = node_pt(q_local)->hanging_pt();
1485 
1486  // Read out number of master nodes from hanging data
1487  n_master2 = hang_info2_pt->nmaster();
1488  }
1489  // Otherwise the node is its own master
1490  else
1491  {
1492  n_master2 = 1;
1493  }
1494 
1495  // Loop over the master nodes
1496  for (unsigned mm = 0; mm < n_master2; mm++)
1497  {
1498  if (is_node_hanging2)
1499  {
1500  actual_shape_controlling_node_pt =
1501  hang_info2_pt->master_node_pt(mm);
1502  hang_weight2 = hang_info2_pt->master_weight(mm);
1503  }
1504 
1505  // Find the corresponding number
1506  unsigned q = local_shape_controlling_node_lookup
1507  [actual_shape_controlling_node_pt];
1508 
1509  // Loop over coordinate directions
1510  for (unsigned p = 0; p < DIM; p++)
1511  {
1512  double sum = -visc_ratio * this->Gamma[i] *
1513  dpsifdx(q_local, i) * dtestfdx(l, p) -
1514  scaled_re * psif(q_local) *
1515  interpolated_dudx(i, p) * testf(l);
1516  if (i == p)
1517  {
1518  sum -= scaled_re_st * val_time_weight * psif(q_local) *
1519  testf(l);
1520  for (unsigned k = 0; k < DIM; k++)
1521  {
1522  sum -=
1523  visc_ratio * dpsifdx(q_local, k) * dtestfdx(l, k);
1524  double tmp = scaled_re * interpolated_u[k];
1525  if (!this->ALE_is_disabled)
1526  {
1527  tmp -= scaled_re_st * mesh_velocity[k];
1528  }
1529  sum -= tmp * dpsifdx(q_local, k) * testf(l);
1530  }
1531  }
1532 
1533  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1534  sum * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1535  }
1536  } // End of loop over master nodes
1537  } // End of loop over local nodes
1538  } // End of if(element_has_node_with_aux_node_update_fct)
1539 
1540 
1541  } // local_eqn>=0
1542  }
1543  }
1544  } // End of loop over test functions
1545 
1546 
1547  // CONTINUITY EQUATION
1548  // -------------------
1549 
1550  // Loop over the shape functions
1551  for (unsigned l = 0; l < n_pres; l++)
1552  {
1553  // If the pressure dof is hanging
1554  if (pressure_dof_is_hanging[l])
1555  {
1556  // Pressure dof is hanging so it must be nodal-based
1557  // Get the hang info object
1558  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1559 
1560  // Get the number of master nodes from the pressure node
1561  n_master = hang_info_pt->nmaster();
1562  }
1563  // Otherwise the node is its own master
1564  else
1565  {
1566  n_master = 1;
1567  }
1568 
1569  // Loop over the master nodes
1570  for (unsigned m = 0; m < n_master; m++)
1571  {
1572  // Get the number of the unknown
1573  // If the pressure dof is hanging
1574  if (pressure_dof_is_hanging[l])
1575  {
1576  // Get the local equation from the master node
1577  local_eqn =
1578  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1579  // Get the weight for the node
1580  hang_weight = hang_info_pt->master_weight(m);
1581  }
1582  else
1583  {
1584  local_eqn = this->p_local_eqn(l);
1585  hang_weight = 1.0;
1586  }
1587 
1588  // If not a boundary conditions
1589  if (local_eqn >= 0)
1590  {
1591  // Loop over coordinate directions
1592  for (unsigned p = 0; p < DIM; p++)
1593  {
1594  // Loop over nodes
1595  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1596  {
1597  // Residual x deriv of Jacobian
1598  //-----------------------------
1599 
1600  // Source term
1601  double aux = -source;
1602 
1603  // Loop over velocity components
1604  for (unsigned k = 0; k < DIM; k++)
1605  {
1606  aux += interpolated_dudx(k, k);
1607  }
1608 
1609  // Multiply through by deriv of Jacobian and integration weight
1610  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1611  aux * dJ_dX(p, q) * testp[l] * w * hang_weight;
1612 
1613 
1614  // Derivative of residual x Jacobian
1615  // ---------------------------------
1616 
1617  // Loop over velocity components
1618  aux = -source_gradient[p] * psif(q);
1619  for (unsigned k = 0; k < DIM; k++)
1620  {
1621  aux += d_dudx_dX(p, q, k, k);
1622  }
1623  // Multiply through by Jacobian and integration weight
1624  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1625  aux * testp[l] * J * w * hang_weight;
1626  }
1627  }
1628 
1629 
1630  // Derivs w.r.t. to nodal velocities
1631  // ---------------------------------
1632  if (element_has_node_with_aux_node_update_fct)
1633  {
1634  // Loop over local nodes
1635  for (unsigned q_local = 0; q_local < n_node; q_local++)
1636  {
1637  // Number of master nodes and storage for the weight of
1638  // the shape function
1639  unsigned n_master2 = 1;
1640  double hang_weight2 = 1.0;
1641  HangInfo* hang_info2_pt = 0;
1642 
1643  // Local boolean to indicate whether the node is hanging
1644  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1645 
1646  Node* actual_shape_controlling_node_pt = node_pt(q_local);
1647 
1648  // If the node is hanging
1649  if (is_node_hanging2)
1650  {
1651  hang_info2_pt = node_pt(q_local)->hanging_pt();
1652 
1653  // Read out number of master nodes from hanging data
1654  n_master2 = hang_info2_pt->nmaster();
1655  }
1656  // Otherwise the node is its own master
1657  else
1658  {
1659  n_master2 = 1;
1660  }
1661 
1662  // Loop over the master nodes
1663  for (unsigned mm = 0; mm < n_master2; mm++)
1664  {
1665  if (is_node_hanging2)
1666  {
1667  actual_shape_controlling_node_pt =
1668  hang_info2_pt->master_node_pt(mm);
1669  hang_weight2 = hang_info2_pt->master_weight(mm);
1670  }
1671 
1672  // Find the corresponding number
1673  unsigned q = local_shape_controlling_node_lookup
1674  [actual_shape_controlling_node_pt];
1675 
1676  // Loop over coordinate directions
1677  for (unsigned p = 0; p < DIM; p++)
1678  {
1679  double aux = dpsifdx(q_local, p) * testp(l);
1680  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1681  aux * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1682  }
1683  } // End of loop over (mm) master nodes
1684  } // End of loop over local nodes q_local
1685  } // End of if(element_has_node_with_aux_node_update_fct)
1686  } // End of if it's not a boundary condition
1687  } // End of loop over (m) master nodes
1688  } // End of loop over shape functions for continuity eqn
1689 
1690  } // End of loop over integration points
1691  }
p
float * p
Definition: Tutorial_Map_using.cpp:9
oomph::GeneralisedElement::Default_fd_jacobian_step
static double Default_fd_jacobian_step
Definition: elements.h:1198
oomph::GeneralisedElement::ndof
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
oomph::GeneralisedNewtonianNavierStokesEquations::get_body_force_gradient_nst
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: generalised_newtonian_navier_stokes_elements.h:277
oomph::GeneralisedNewtonianNavierStokesEquations::get_source_gradient_nst
virtual void get_source_gradient_nst(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
Definition: generalised_newtonian_navier_stokes_elements.h:340
oomph::Node::position_time_stepper_pt
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
oomph::RefineableElement::nshape_controlling_nodes
unsigned nshape_controlling_nodes()
Definition: refineable_elements.h:627
oomph::RefineableElement::shape_controlling_node_lookup
std::map< Node *, unsigned > shape_controlling_node_lookup()
Definition: refineable_elements.h:636
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019

References oomph::ALE_is_disabled, Global_Parameters::body_force(), oomph::GeneralisedElement::Default_fd_jacobian_step, DIM, G, GlobalPhysicalVariables::Gamma, oomph::Node::has_auxiliary_node_update_fct_pt(), i, J, j, k, m, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::HangInfo::nmaster(), 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_pressure_and_velocity_mass_matrix_diagonal()

template<unsigned DIM>
void oomph::RefineableGeneralisedNewtonianNavierStokesEquations< 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)

Reimplemented from oomph::GeneralisedNewtonianNavierStokesEquations< DIM >.

43  {
44  // Resize and initialise
45  unsigned n_dof = ndof();
46 
47  if ((which_one == 0) || (which_one == 1))
48  press_mass_diag.assign(n_dof, 0.0);
49  if ((which_one == 0) || (which_one == 2))
50  veloc_mass_diag.assign(n_dof, 0.0);
51 
52  // Pointers to hang info object
53  HangInfo* hang_info_pt = 0;
54 
55  // Number of master nodes and weight for shape fcts
56  unsigned n_master = 1;
57  double hang_weight = 1.0;
58 
59  // find out how many nodes there are
60  unsigned n_node = nnode();
61 
62  // Set up memory for veloc shape functions
63  Shape psi(n_node);
64 
65  // Find number of pressure dofs
66  unsigned n_pres = this->npres_nst();
67 
68  // Pressure shape function
69  Shape psi_p(n_pres);
70 
71  // Local coordinates
72  Vector<double> s(DIM);
73 
74  // find the indices at which the local velocities are stored
75  Vector<unsigned> u_nodal_index(DIM);
76  for (unsigned i = 0; i < DIM; i++)
77  {
78  u_nodal_index[i] = this->u_index_nst(i);
79  }
80 
81  // Which nodal value represents the pressure? (Negative if pressure
82  // is not based on nodal interpolation).
83  int p_index = this->p_nodal_index_nst();
84 
85  // Local array of booleans that are true if the l-th pressure value is
86  // hanging (avoid repeated virtual function calls)
87  bool pressure_dof_is_hanging[n_pres];
88 
89  // If the pressure is stored at a node
90  if (p_index >= 0)
91  {
92  // Read out whether the pressure is hanging
93  for (unsigned l = 0; l < n_pres; ++l)
94  {
95  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
96  }
97  }
98  // Otherwise the pressure is not stored at a node and so cannot hang
99  else
100  {
101  for (unsigned l = 0; l < n_pres; ++l)
102  {
103  pressure_dof_is_hanging[l] = false;
104  }
105  }
106 
107 
108  // Number of integration points
109  unsigned n_intpt = integral_pt()->nweight();
110 
111  // Integer to store the local equations no
112  int local_eqn = 0;
113 
114  // Loop over the integration points
115  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
116  {
117  // Get the integral weight
118  double w = integral_pt()->weight(ipt);
119 
120  // Get determinant of Jacobian of the mapping
121  double J = J_eulerian_at_knot(ipt);
122 
123  // Assign values of s
124  for (unsigned i = 0; i < DIM; i++)
125  {
126  s[i] = integral_pt()->knot(ipt, i);
127  }
128 
129  // Premultiply weights and Jacobian
130  double W = w * J;
131 
132 
133  // Do we want the velocity one?
134  if ((which_one == 0) || (which_one == 2))
135  {
136  // Get the velocity shape functions
137  shape_at_knot(ipt, psi);
138 
139 
140  // Number of master nodes and storage for the weight of the shape
141  // function
142  unsigned n_master = 1;
143  double hang_weight = 1.0;
144 
145  // Loop over the nodes for the test functions/equations
146  //----------------------------------------------------
147  for (unsigned l = 0; l < n_node; l++)
148  {
149  // Local boolean to indicate whether the node is hanging
150  bool is_node_hanging = node_pt(l)->is_hanging();
151 
152  // If the node is hanging
153  if (is_node_hanging)
154  {
155  hang_info_pt = node_pt(l)->hanging_pt();
156 
157  // Read out number of master nodes from hanging data
158  n_master = hang_info_pt->nmaster();
159  }
160  // Otherwise the node is its own master
161  else
162  {
163  n_master = 1;
164  }
165 
166  // Loop over the master nodes
167  for (unsigned m = 0; m < n_master; m++)
168  {
169  // Loop over velocity components for equations
170  for (unsigned i = 0; i < DIM; i++)
171  {
172  // Get the equation number
173  // If the node is hanging
174  if (is_node_hanging)
175  {
176  // Get the equation number from the master node
177  local_eqn = this->local_hang_eqn(
178  hang_info_pt->master_node_pt(m), u_nodal_index[i]);
179  // Get the hang weight from the master node
180  hang_weight = hang_info_pt->master_weight(m);
181  }
182  // If the node is not hanging
183  else
184  {
185  // Local equation number
186  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
187 
188  // Node contributes with full weight
189  hang_weight = 1.0;
190  }
191 
192  // If it's not a boundary condition...
193  if (local_eqn >= 0)
194  {
195  // //Loop over the veclocity shape functions
196  // for(unsigned l=0; l<n_node; l++)
197  // {
198  // //Loop over the velocity components
199  // for(unsigned i=0; i<DIM; i++)
200  // {
201  // local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
202 
203  // //If not a boundary condition
204  // if(local_eqn >= 0)
205  // {
206 
207 
208  // Add the contribution
209  veloc_mass_diag[local_eqn] += pow(psi[l] * hang_weight, 2) * W;
210  }
211  }
212  }
213  }
214  }
215 
216  // Do we want the pressure one?
217  if ((which_one == 0) || (which_one == 1))
218  {
219  // Get the pressure shape functions
220  this->pshape_nst(s, psi_p);
221 
222  // Loop over the pressure shape functions
223  for (unsigned l = 0; l < n_pres; l++)
224  {
225  // If the pressure dof is hanging
226  if (pressure_dof_is_hanging[l])
227  {
228  // Pressure dof is hanging so it must be nodal-based
229  // Get the hang info object
230  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
231 
232  // Get the number of master nodes from the pressure node
233  n_master = hang_info_pt->nmaster();
234  }
235  // Otherwise the node is its own master
236  else
237  {
238  n_master = 1;
239  }
240 
241  // Loop over the master nodes
242  for (unsigned m = 0; m < n_master; m++)
243  {
244  // Get the number of the unknown
245  // If the pressure dof is hanging
246  if (pressure_dof_is_hanging[l])
247  {
248  // Get the local equation from the master node
249  local_eqn =
250  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
251  // Get the weight for the node
252  hang_weight = hang_info_pt->master_weight(m);
253  }
254  else
255  {
256  local_eqn = this->p_local_eqn(l);
257  hang_weight = 1.0;
258  }
259 
260  // If the equation is not pinned
261  if (local_eqn >= 0)
262  {
263  // //Loop over the veclocity shape functions
264  // for(unsigned l=0; l<n_pres; l++)
265  // {
266  // // Get equation number
267  // local_eqn = p_local_eqn(l);
268 
269  // //If not a boundary condition
270  // if(local_eqn >= 0)
271  // {
272 
273 
274  // Add the contribution
275  press_mass_diag[local_eqn] += pow(psi_p[l] * hang_weight, 2) * W;
276  }
277  }
278  }
279  }
280  }
281  }
oomph::FiniteElement::J_eulerian_at_knot
virtual double J_eulerian_at_knot(const unsigned &ipt) const
Definition: elements.cc:4168
oomph::FiniteElement::shape_at_knot
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Definition: elements.cc:3220
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625

References DIM, i, J, m, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::HangInfo::nmaster(), Eigen::bfloat16_impl::pow(), s, w, and oomph::QuadTreeNames::W.

◆ get_Z2_flux()

template<unsigned DIM>
void oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::get_Z2_flux ( const Vector< double > &  s,
Vector< double > &  flux 
)
inlinevirtual

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

Implements oomph::ElementWithZ2ErrorEstimator.

144  {
145 #ifdef PARANOID
146  unsigned num_entries = DIM + (DIM * (DIM - 1)) / 2;
147  if (flux.size() < num_entries)
148  {
149  std::ostringstream error_message;
150  error_message << "The flux vector has the wrong number of entries, "
151  << flux.size() << ", whereas it should be at least "
152  << num_entries << std::endl;
153  throw OomphLibError(error_message.str(),
154  OOMPH_CURRENT_FUNCTION,
155  OOMPH_EXCEPTION_LOCATION);
156  }
157 #endif
158 
159  // Get strain rate matrix
160  DenseMatrix<double> strainrate(DIM);
161  this->strain_rate(s, strainrate);
162 
163  // Pack into flux Vector
164  unsigned icount = 0;
165 
166  // Start with diagonal terms
167  for (unsigned i = 0; i < DIM; i++)
168  {
169  flux[icount] = strainrate(i, i);
170  icount++;
171  }
172 
173  // Off diagonals row by row
174  for (unsigned i = 0; i < DIM; i++)
175  {
176  for (unsigned j = i + 1; j < DIM; j++)
177  {
178  flux[icount] = strainrate(i, j);
179  icount++;
180  }
181  }
182  }
ProblemParameters::flux
void flux(const double &time, const Vector< double > &x, double &flux)
Get flux applied along boundary x=0.
Definition: pretend_melt.cc:59
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

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

◆ num_Z2_flux_terms()

template<unsigned DIM>
unsigned oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::num_Z2_flux_terms ( )
inlinevirtual

Number of 'flux' terms for Z2 error estimation.

Implements oomph::ElementWithZ2ErrorEstimator.

136  {
137  // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
138  return DIM + (DIM * (DIM - 1)) / 2;
139  }

References DIM.

◆ pin_elemental_redundant_nodal_pressure_dofs()

template<unsigned DIM>
virtual void oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::pin_elemental_redundant_nodal_pressure_dofs ( )
inlineprotectedvirtual

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

Reimplemented in oomph::RefineableGeneralisedNewtonianQTaylorHoodElement< DIM >.

67 {}

Referenced by oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::pin_redundant_nodal_pressures().

◆ pin_redundant_nodal_pressures()

template<unsigned DIM>
static void oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::pin_redundant_nodal_pressures ( const Vector< GeneralisedElement * > &  element_pt)
inlinestatic

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

  • RefineableNavierStokesEquations:: pin_elemental_redundant_nodal_pressure_dofs()

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

91  {
92  // Loop over all elements and call the function that pins their
93  // unused nodal pressure data
94  unsigned n_element = element_pt.size();
95  for (unsigned e = 0; e < n_element; e++)
96  {
97  dynamic_cast<RefineableGeneralisedNewtonianNavierStokesEquations<DIM>*>(
98  element_pt[e])
99  ->pin_elemental_redundant_nodal_pressure_dofs();
100  }
101  }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::RefineableGeneralisedNewtonianNavierStokesEquations::pin_elemental_redundant_nodal_pressure_dofs
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Definition: generalised_newtonian_refineable_navier_stokes_elements.h:67

References e(), and oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::pin_elemental_redundant_nodal_pressure_dofs().

◆ pressure_node_pt()

template<unsigned DIM>
virtual Node* oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::pressure_node_pt ( const unsigned &  n_p)
inlinevirtual

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

Reimplemented in oomph::RefineableGeneralisedNewtonianQTaylorHoodElement< DIM >.

120  {
121  return NULL;
122  }

◆ unpin_all_pressure_dofs()

template<unsigned DIM>
static void oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::unpin_all_pressure_dofs ( const Vector< GeneralisedElement * > &  element_pt)
inlinestatic

Unpin all pressure dofs in elements listed in vector.

106  {
107  // Loop over all elements
108  unsigned n_element = element_pt.size();
109  for (unsigned e = 0; e < n_element; e++)
110  {
111  dynamic_cast<RefineableGeneralisedNewtonianNavierStokesEquations<DIM>*>(
112  element_pt[e])
113  ->unpin_elemental_pressure_dofs();
114  }
115  }
oomph::RefineableGeneralisedNewtonianNavierStokesEquations::unpin_elemental_pressure_dofs
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.

References e(), and oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::unpin_elemental_pressure_dofs().

◆ unpin_elemental_pressure_dofs()

template<unsigned DIM>
virtual void oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::unpin_elemental_pressure_dofs ( )
protectedpure virtual

Unpin all pressure dofs in the element.

Implemented in oomph::PRefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >, oomph::RefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM >, and oomph::RefineableGeneralisedNewtonianQTaylorHoodElement< DIM >.

Referenced by oomph::RefineableGeneralisedNewtonianNavierStokesEquations< DIM >::unpin_all_pressure_dofs().

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