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

#include <refineable_navier_stokes_elements.h>

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

Public Member Functions

 RefineableNavierStokesEquations ()
 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::NavierStokesEquations< DIM >
 NavierStokesEquations ()
 
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...
 
NavierStokesPressureAdvDiffSourceFctPt & source_fct_for_pressure_adv_diff ()
 
NavierStokesPressureAdvDiffSourceFctPt source_fct_for_pressure_adv_diff () const
 
int & pinned_fp_pressure_eqn ()
 
virtual unsigned npres_nst () const =0
 Function to return number of pressure degrees of freedom. More...
 
virtual void pshape_nst (const Vector< double > &s, Shape &psi) const =0
 Compute the pressure shape functions at local coordinate s. More...
 
virtual void pshape_nst (const Vector< double > &s, Shape &psi, Shape &test) const =0
 
virtual double dpshape_and_dptest_eulerian_nst (const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const =0
 
double u_nst (const unsigned &n, const unsigned &i) const
 
double u_nst (const unsigned &t, const unsigned &n, const unsigned &i) const
 
virtual unsigned u_index_nst (const unsigned &i) const
 
unsigned n_u_nst () const
 
double du_dt_nst (const unsigned &n, const unsigned &i) const
 
void disable_ALE ()
 
void enable_ALE ()
 
virtual double p_nst (const unsigned &n_p) const =0
 
virtual double p_nst (const unsigned &t, const unsigned &n_p) const =0
 Pressure at local pressure "node" n_p at time level t. More...
 
virtual void fix_pressure (const unsigned &p_dof, const double &p_value)=0
 Pin p_dof-th pressure dof and set it to value specified by p_value. More...
 
virtual int p_nodal_index_nst () const
 
double pressure_integral () const
 Integral of pressure over element. More...
 
double dissipation () const
 Return integral of dissipation over element. More...
 
double dissipation (const Vector< double > &s) const
 Return dissipation at local coordinate s. More...
 
void get_vorticity (const Vector< double > &s, Vector< double > &vorticity) const
 Compute the vorticity vector at local coordinate s. More...
 
void get_vorticity (const Vector< double > &s, double &vorticity) const
 Compute the scalar vorticity at local coordinate s (2D) More...
 
double kin_energy () const
 Get integral of kinetic energy over element. More...
 
double d_kin_energy_dt () const
 Get integral of time derivative of kinetic energy over element. More...
 
void strain_rate (const Vector< double > &s, DenseMatrix< double > &strain_rate) const
 Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i) More...
 
void get_traction (const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
 
void get_traction (const Vector< double > &s, const Vector< double > &N, Vector< double > &traction_p, Vector< double > &traction_visc_n, Vector< double > &traction_visc_t)
 
void get_load (const Vector< double > &s, const Vector< double > &N, Vector< double > &load)
 
unsigned nscalar_paraview () const
 
void scalar_value_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
 
void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
std::string scalar_name_paraview (const unsigned &i) const
 
void output (std::ostream &outfile)
 
void output (std::ostream &outfile, const unsigned &nplot)
 
void output (FILE *file_pt)
 
void output (FILE *file_pt, const unsigned &nplot)
 
void full_output (std::ostream &outfile)
 
void full_output (std::ostream &outfile, const unsigned &nplot)
 
void output_veloc (std::ostream &outfile, const unsigned &nplot, const unsigned &t)
 
void output_vorticity (std::ostream &outfile, const unsigned &nplot)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 
void compute_norm (double &norm)
 Compute norm of solution: square of the L2 norm of the velocities. More...
 
void compute_norm (Vector< double > &norm)
 Compute the vector norm of the FEM solution. More...
 
void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
void fill_in_contribution_to_residuals (Vector< double > &residuals)
 Compute the element's residual Vector. More...
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void fill_in_contribution_to_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
void fill_in_contribution_to_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 Compute the element's residual Vector. More...
 
void fill_in_contribution_to_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
void fill_in_pressure_advection_diffusion_residuals (Vector< double > &residuals)
 
void fill_in_pressure_advection_diffusion_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void pin_all_non_pressure_dofs (std::map< Data *, std::vector< int >> &eqn_number_backup)
 Pin all non-pressure dofs and backup eqn numbers. More...
 
virtual void build_fp_press_adv_diff_robin_bc_element (const unsigned &face_index)=0
 
void output_pressure_advection_diffusion_robin_elements (std::ostream &outfile)
 
void delete_pressure_advection_diffusion_robin_elements ()
 
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...
 
double interpolated_dudx_nst (const Vector< double > &s, const unsigned &i, const unsigned &j) const
 
void point_output_data (const Vector< double > &s, Vector< double > &data)
 
void get_vorticity (const Vector< double > &s, Vector< double > &vorticity) const
 
void get_vorticity (const Vector< double > &s, double &vorticity) const
 
void get_vorticity (const Vector< double > &s, Vector< double > &vorticity) const
 Compute 3D vorticity vector at local coordinate s. More...
 
- Public Member Functions inherited from oomph::FSIFluidElement
 FSIFluidElement ()
 Constructor. More...
 
 FSIFluidElement (const FSIFluidElement &)=delete
 Broken copy constructor. More...
 
void operator= (const FSIFluidElement &)=delete
 Broken assignment operator. More...
 
virtual void identify_load_data (std::set< std::pair< Data *, unsigned >> &paired_load_data)=0
 
virtual void identify_pressure_data (std::set< std::pair< Data *, unsigned >> &paired_pressure_data)=0
 
- Public Member Functions inherited from oomph::FiniteElement
void set_dimension (const unsigned &dim)
 
void set_nodal_dimension (const unsigned &nodal_dim)
 
void set_nnodal_position_type (const unsigned &nposition_type)
 Set the number of types required to interpolate the coordinate. More...
 
void set_n_node (const unsigned &n)
 
int nodal_local_eqn (const unsigned &n, const unsigned &i) const
 
double dJ_eulerian_at_knot (const unsigned &ipt, Shape &psi, DenseMatrix< double > &djacobian_dX) const
 
 FiniteElement ()
 Constructor. More...
 
virtual ~FiniteElement ()
 
 FiniteElement (const FiniteElement &)=delete
 Broken copy constructor. More...
 
virtual bool local_coord_is_valid (const Vector< double > &s)
 Broken assignment operator. More...
 
virtual void move_local_coord_back_into_element (Vector< double > &s) const
 
void get_centre_of_gravity_and_max_radius_in_terms_of_zeta (Vector< double > &cog, double &max_radius) const
 
virtual void local_coordinate_of_node (const unsigned &j, Vector< double > &s) const
 
virtual void local_fraction_of_node (const unsigned &j, Vector< double > &s_fraction)
 
virtual double local_one_d_fraction_of_node (const unsigned &n1d, const unsigned &i)
 
virtual void set_macro_elem_pt (MacroElement *macro_elem_pt)
 
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 output (const unsigned &t, std::ostream &outfile, const unsigned &n_plot) const
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, const SolutionFunctorBase &exact_soln) const
 Output a time-dependent exact solution over the element. More...
 
virtual void get_s_plot (const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
 
virtual std::string tecplot_zone_string (const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (std::ostream &outfile, const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (FILE *file_pt, const unsigned &nplot) const
 
virtual unsigned nplot_points (const unsigned &nplot) const
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_abs_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error)
 
void integrate_fct (FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt, Vector< double > &integral)
 Integrate Vector-valued function over element. More...
 
void integrate_fct (FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt, const double &time, Vector< double > &integral)
 Integrate Vector-valued time-dep function over element. More...
 
virtual void build_face_element (const int &face_index, FaceElement *face_element_pt)
 
virtual unsigned self_test ()
 
virtual unsigned get_bulk_node_number (const int &face_index, const unsigned &i) const
 
virtual int face_outer_unit_normal_sign (const int &face_index) const
 Get the sign of the outer unit normal on the face given by face_index. More...
 
virtual unsigned nnode_on_face () const
 
void face_node_number_error_check (const unsigned &i) const
 Range check for face node numbers. More...
 
virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt (const int &face_index) const
 
virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt (const int &face_index) const
 
- Public Member Functions inherited from oomph::GeneralisedElement
 GeneralisedElement ()
 Constructor: Initialise all pointers and all values to zero. More...
 
virtual ~GeneralisedElement ()
 Virtual destructor to clean up any memory allocated by the object. More...
 
 GeneralisedElement (const GeneralisedElement &)=delete
 Broken copy constructor. More...
 
void operator= (const GeneralisedElement &)=delete
 Broken assignment operator. More...
 
Data *& internal_data_pt (const unsigned &i)
 Return a pointer to i-th internal data object. More...
 
Data *const & internal_data_pt (const unsigned &i) const
 Return a pointer to i-th internal data object (const version) More...
 
Data *& external_data_pt (const unsigned &i)
 Return a pointer to i-th external data object. More...
 
Data *const & external_data_pt (const unsigned &i) const
 Return a pointer to i-th external data object (const version) More...
 
unsigned long eqn_number (const unsigned &ieqn_local) const
 
int local_eqn_number (const unsigned long &ieqn_global) const
 
unsigned add_external_data (Data *const &data_pt, const bool &fd=true)
 
bool external_data_fd (const unsigned &i) const
 
void exclude_external_data_fd (const unsigned &i)
 
void include_external_data_fd (const unsigned &i)
 
void flush_external_data ()
 Flush all external data. More...
 
void flush_external_data (Data *const &data_pt)
 Flush the object addressed by data_pt from the external data array. More...
 
unsigned ninternal_data () const
 Return the number of internal data objects. More...
 
unsigned nexternal_data () const
 Return the number of external data objects. More...
 
unsigned ndof () const
 Return the number of equations/dofs in the element. More...
 
void dof_vector (const unsigned &t, Vector< double > &dof)
 Return the vector of dof values at time level t. More...
 
void dof_pt_vector (Vector< double * > &dof_pt)
 Return the vector of pointers to dof values. More...
 
void set_internal_data_time_stepper (const unsigned &i, TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void assign_internal_eqn_numbers (unsigned long &global_number, Vector< double * > &Dof_pt)
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void add_internal_value_pt_to_map (std::map< unsigned, double * > &map_of_value_pt)
 
virtual void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void complete_setup_of_dependencies ()
 
virtual void get_residuals (Vector< double > &residuals)
 
virtual void get_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
virtual void get_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void get_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void get_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void get_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void get_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void get_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void get_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void get_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
virtual unsigned ndof_types () const
 
virtual void get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
 
- Public Member Functions inherited from oomph::GeomObject
 GeomObject ()
 Default constructor. More...
 
 GeomObject (const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)
 
 GeomObject (const GeomObject &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const GeomObject &)=delete
 Broken assignment operator. More...
 
virtual ~GeomObject ()
 (Empty) destructor More...
 
unsigned nlagrangian () const
 Access function to # of Lagrangian coordinates. More...
 
unsigned ndim () const
 Access function to # of Eulerian coordinates. More...
 
void set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
 Set # of Lagrangian and Eulerian coordinates. More...
 
TimeStepper *& time_stepper_pt ()
 
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::TemplateFreeNavierStokesEquationsBase
 TemplateFreeNavierStokesEquationsBase ()
 Constructor (empty) More...
 
virtual ~TemplateFreeNavierStokesEquationsBase ()
 Virtual destructor (empty) More...
 
virtual int p_local_eqn (const unsigned &n) const =0
 
- Public Member Functions inherited from oomph::NavierStokesElementWithDiagonalMassMatrices
 NavierStokesElementWithDiagonalMassMatrices ()
 Empty constructor. More...
 
virtual ~NavierStokesElementWithDiagonalMassMatrices ()
 Virtual destructor. More...
 
 NavierStokesElementWithDiagonalMassMatrices (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken copy constructor. More...
 
void operator= (const NavierStokesElementWithDiagonalMassMatrices &)=delete
 Broken assignment operator. More...
 
- 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 ()
 
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)
 
void fill_in_generic_pressure_advection_diffusion_contribution_nst (Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
 
virtual void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
- Protected Member Functions inherited from oomph::NavierStokesEquations< 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 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)
 
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::NavierStokesEquations< 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::NavierStokesEquations< 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::NavierStokesEquations< 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...
 
NavierStokesPressureAdvDiffSourceFctPt Press_adv_diff_source_fct_pt
 
bool ALE_is_disabled
 
Vector< FpPressureAdvDiffRobinBCElementBase * > Pressure_advection_diffusion_robin_element_pt
 
int Pinned_fp_pressure_eqn
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElement * Macro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepper * Geom_object_time_stepper_pt
 
- Protected Attributes inherited from oomph::RefineableElement
Tree * Tree_pt
 A pointer to a general tree object. More...
 
unsigned Refine_level
 Refinement level. More...
 
bool To_be_refined
 Flag for refinement. More...
 
bool Refinement_is_enabled
 Flag to indicate suppression of any refinement. More...
 
bool Sons_to_be_unrefined
 Flag for unrefinement. More...
 
long Number
 Global element number – for plotting/validation purposes. More...
 
- Static Protected Attributes inherited from oomph::FiniteElement
static const unsigned Default_Initial_Nvalue = 0
 Default value for the number of values at a node. More...
 
static const double Node_location_tolerance = 1.0e-14
 
static const unsigned N2deriv [] = {0, 1, 3, 6}
 
- Static Protected Attributes inherited from oomph::GeneralisedElement
static DenseMatrix< double > Dummy_matrix
 
static std::deque< double * > Dof_pt_deque
 
- Static Protected Attributes inherited from oomph::RefineableElement
static double Max_integrity_tolerance = 1.0e-8
 Max. allowed discrepancy in element integrity check. More...
 

Detailed Description

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

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

Constructor & Destructor Documentation

◆ RefineableNavierStokesEquations()

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

Constructor.

334  : NavierStokesEquations<DIM>(),
335  RefineableElement(),
336  ElementWithZ2ErrorEstimator()
337  {
338  }
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::RefineableNavierStokesEquations< 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::NavierStokesEquations< DIM >.

486  {
487  // Find number of nodes
488  unsigned n_node = this->nnode();
489  // Local shape function
490  Shape psi(n_node);
491  // Find values of shape function at the given local coordinate
492  this->shape(s, psi);
493 
494  // Find the index at which the velocity component is stored
495  const unsigned u_nodal_index = this->u_index_nst(i);
496 
497  // Storage for hang info pointer
498  HangInfo* hang_info_pt = 0;
499  // Storage for global equation
500  int global_eqn = 0;
501 
502  // Find the number of dofs associated with interpolated u
503  unsigned n_u_dof = 0;
504  for (unsigned l = 0; l < n_node; l++)
505  {
506  unsigned n_master = 1;
507 
508  // Local bool (is the node hanging)
509  bool is_node_hanging = this->node_pt(l)->is_hanging();
510 
511  // If the node is hanging, get the number of master nodes
512  if (is_node_hanging)
513  {
514  hang_info_pt = this->node_pt(l)->hanging_pt();
515  n_master = hang_info_pt->nmaster();
516  }
517  // Otherwise there is just one master node, the node itself
518  else
519  {
520  n_master = 1;
521  }
522 
523  // Loop over the master nodes
524  for (unsigned m = 0; m < n_master; m++)
525  {
526  // Get the equation number
527  if (is_node_hanging)
528  {
529  // Get the equation number from the master node
530  global_eqn =
531  hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
532  }
533  else
534  {
535  // Global equation number
536  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
537  }
538 
539  // If it's positive add to the count
540  if (global_eqn >= 0)
541  {
542  ++n_u_dof;
543  }
544  }
545  }
546 
547  // Now resize the storage schemes
548  du_ddata.resize(n_u_dof, 0.0);
549  global_eqn_number.resize(n_u_dof, 0);
550 
551  // Loop over th nodes again and set the derivatives
552  unsigned count = 0;
553  // Loop over the local nodes and sum
554  for (unsigned l = 0; l < n_node; l++)
555  {
556  unsigned n_master = 1;
557  double hang_weight = 1.0;
558 
559  // Local bool (is the node hanging)
560  bool is_node_hanging = this->node_pt(l)->is_hanging();
561 
562  // If the node is hanging, get the number of master nodes
563  if (is_node_hanging)
564  {
565  hang_info_pt = this->node_pt(l)->hanging_pt();
566  n_master = hang_info_pt->nmaster();
567  }
568  // Otherwise there is just one master node, the node itself
569  else
570  {
571  n_master = 1;
572  }
573 
574  // Loop over the master nodes
575  for (unsigned m = 0; m < n_master; m++)
576  {
577  // If the node is hanging get weight from master node
578  if (is_node_hanging)
579  {
580  // Get the hang weight from the master node
581  hang_weight = hang_info_pt->master_weight(m);
582  }
583  else
584  {
585  // Node contributes with full weight
586  hang_weight = 1.0;
587  }
588 
589  // Get the equation number
590  if (is_node_hanging)
591  {
592  // Get the equation number from the master node
593  global_eqn =
594  hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
595  }
596  else
597  {
598  // Local equation number
599  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
600  }
601 
602  if (global_eqn >= 0)
603  {
604  // Set the global equation number
605  global_eqn_number[count] = global_eqn;
606  // Set the derivative with respect to the unknown
607  du_ddata[count] = psi[l] * hang_weight;
608  // Increase the counter
609  ++count;
610  }
611  }
612  }
613  }
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::NavierStokesEquations::u_index_nst
virtual unsigned u_index_nst(const unsigned &i) const
Definition: navier_stokes_elements.h:866
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::NavierStokesEquations< DIM >::u_index_nst().

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

◆ fill_in_generic_pressure_advection_diffusion_contribution_nst()

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

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

Reimplemented from oomph::NavierStokesEquations< DIM >.

293  {
294  // Return immediately if there are no dofs
295  if (ndof() == 0) return;
296 
297  // Find out how many nodes there are
298  unsigned n_node = nnode();
299 
300  // Find out how many pressure dofs there are
301  unsigned n_pres = this->npres_nst();
302 
303  // Find the indices at which the local velocities are stored
304  unsigned u_nodal_index[DIM];
305  for (unsigned i = 0; i < DIM; i++)
306  {
307  u_nodal_index[i] = this->u_index_nst(i);
308  }
309 
310 
311  // Which nodal value represents the pressure? (Negative if pressure
312  // is not based on nodal interpolation).
313  int p_index = this->p_nodal_index_nst();
314 
315  // Local array of booleans that are true if the l-th pressure value is
316  // hanging (avoid repeated virtual function calls)
317  bool pressure_dof_is_hanging[n_pres];
318  // If the pressure is stored at a node
319  if (p_index >= 0)
320  {
321  // Read out whether the pressure is hanging
322  for (unsigned l = 0; l < n_pres; ++l)
323  {
324  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
325  }
326  }
327  // Otherwise the pressure is not stored at a node and so cannot hang
328  else
329  {
330  // pressure advection diffusion doesn't work for this one!
331  throw OomphLibError(
332  "Pressure advection diffusion does not work in this case\n",
333  OOMPH_CURRENT_FUNCTION,
334  OOMPH_EXCEPTION_LOCATION);
335 
336  for (unsigned l = 0; l < n_pres; ++l)
337  {
338  pressure_dof_is_hanging[l] = false;
339  }
340  }
341 
342  // Set up memory for the velocity shape fcts
343  Shape psif(n_node);
344  DShape dpsidx(n_node, DIM);
345 
346  // Set up memory for pressure shape and test functions
347  Shape psip(n_pres), testp(n_pres);
348  DShape dpsip(n_pres, DIM);
349  DShape dtestp(n_pres, DIM);
350 
351  // Number of integration points
352  unsigned n_intpt = integral_pt()->nweight();
353 
354  // Set the Vector to hold local coordinates
355  Vector<double> s(DIM);
356 
357  // Get Physical Variables from Element
358  // Reynolds number must be multiplied by the density ratio
359  double scaled_re = this->re() * this->density_ratio();
360 
361  // Integers to store the local equations and unknowns
362  int local_eqn = 0, local_unknown = 0;
363 
364  // Pointers to hang info objects
365  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
366 
367  // Loop over the integration points
368  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
369  {
370  // Assign values of s
371  for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
372 
373  // Get the integral weight
374  double w = integral_pt()->weight(ipt);
375 
376  // Call the derivatives of the veloc shape functions
377  // (Derivs not needed but they are free)
378  double J = this->dshape_eulerian_at_knot(ipt, psif, dpsidx);
379 
380  // Call the pressure shape and test functions
381  this->dpshape_and_dptest_eulerian_nst(s, psip, dpsip, testp, dtestp);
382 
383  // Premultiply the weights and the Jacobian
384  double W = w * J;
385 
386  // Calculate local values of the pressure and velocity components
387  // Allocate
388  Vector<double> interpolated_u(DIM, 0.0);
389  Vector<double> interpolated_x(DIM, 0.0);
390  Vector<double> interpolated_dpdx(DIM, 0.0);
391 
392  // Calculate pressure gradient
393  for (unsigned l = 0; l < n_pres; l++)
394  {
395  for (unsigned i = 0; i < DIM; i++)
396  {
397  interpolated_dpdx[i] += this->p_nst(l) * dpsip(l, i);
398  }
399  }
400 
401  // Calculate velocities
402 
403  // Loop over nodes
404  for (unsigned l = 0; l < n_node; l++)
405  {
406  // Loop over directions
407  for (unsigned i = 0; i < DIM; i++)
408  {
409  // Get the nodal value
410  double u_value = nodal_value(l, u_nodal_index[i]);
411  interpolated_u[i] += u_value * psif[l];
412  interpolated_x[i] += nodal_position(l, i) * psif[l];
413  }
414  }
415 
416  // Source function (for validaton only)
417  double source = 0.0;
418  if (this->Press_adv_diff_source_fct_pt != 0)
419  {
420  source = this->Press_adv_diff_source_fct_pt(interpolated_x);
421  }
422 
423 
424  // Number of master nodes and storage for the weight of the shape function
425  unsigned n_master = 1;
426  double hang_weight = 1.0;
427 
428 
429  // Loop over the pressure shape functions
430  for (unsigned l = 0; l < n_pres; l++)
431  {
432  // If the pressure dof is hanging
433  if (pressure_dof_is_hanging[l])
434  {
435  // Pressure dof is hanging so it must be nodal-based
436  // Get the hang info object
437  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
438 
439  // Get the number of master nodes from the pressure node
440  n_master = hang_info_pt->nmaster();
441  }
442  // Otherwise the node is its own master
443  else
444  {
445  n_master = 1;
446  }
447 
448  // Loop over the master nodes
449  for (unsigned m = 0; m < n_master; m++)
450  {
451  // Get the number of the unknown
452  // If the pressure dof is hanging
453  if (pressure_dof_is_hanging[l])
454  {
455  // Get the local equation from the master node
456  local_eqn =
457  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
458  // Get the weight for the node
459  hang_weight = hang_info_pt->master_weight(m);
460  }
461  else
462  {
463  local_eqn = this->p_local_eqn(l);
464  hang_weight = 1.0;
465  }
466 
467  // If the equation is not pinned
468  if (local_eqn >= 0)
469  {
470  residuals[local_eqn] -= source * testp[l] * W * hang_weight;
471  for (unsigned k = 0; k < DIM; k++)
472  {
473  residuals[local_eqn] +=
474  interpolated_dpdx[k] *
475  (scaled_re * interpolated_u[k] * testp[l] + dtestp(l, k)) * W *
476  hang_weight;
477  }
478 
479  // Jacobian too?
480  if (flag)
481  {
482  // Number of master nodes and weights
483  unsigned n_master2 = 1;
484  double hang_weight2 = 1.0;
485 
486  // Loop over the pressure shape functions
487  for (unsigned l2 = 0; l2 < n_pres; l2++)
488  {
489  // If the pressure dof is hanging
490  if (pressure_dof_is_hanging[l2])
491  {
492  hang_info2_pt =
493  this->pressure_node_pt(l2)->hanging_pt(p_index);
494  // Pressure dof is hanging so it must be nodal-based
495  // Get the number of master nodes from the pressure node
496  n_master2 = hang_info2_pt->nmaster();
497  }
498  // Otherwise the node is its own master
499  else
500  {
501  n_master2 = 1;
502  }
503 
504  // Loop over the master nodes
505  for (unsigned m2 = 0; m2 < n_master2; m2++)
506  {
507  // Get the number of the unknown
508  // If the pressure dof is hanging
509  if (pressure_dof_is_hanging[l2])
510  {
511  // Get the unknown from the master node
512  local_unknown = this->local_hang_eqn(
513  hang_info2_pt->master_node_pt(m2), p_index);
514  // Get the weight from the hanging object
515  hang_weight2 = hang_info2_pt->master_weight(m2);
516  }
517  else
518  {
519  local_unknown = this->p_local_eqn(l2);
520  hang_weight2 = 1.0;
521  }
522 
523  // If the unknown is not pinned
524  if (local_unknown >= 0)
525  {
526  if ((int(eqn_number(local_eqn)) !=
527  this->Pinned_fp_pressure_eqn) &&
528  (int(eqn_number(local_unknown)) !=
529  this->Pinned_fp_pressure_eqn))
530  {
531  for (unsigned k = 0; k < DIM; k++)
532  {
533  jacobian(local_eqn, local_unknown) +=
534  dtestp(l2, k) *
535  (scaled_re * interpolated_u[k] * testp[l] +
536  dtestp(l, k)) *
537  W * hang_weight * hang_weight2;
538  }
539  }
540  else
541  {
542  if ((int(eqn_number(local_eqn)) ==
543  this->Pinned_fp_pressure_eqn) &&
544  (int(eqn_number(local_unknown)) ==
545  this->Pinned_fp_pressure_eqn))
546  {
547  jacobian(local_eqn, local_unknown) = 1.0;
548  }
549  }
550  }
551  }
552  }
553  } /*End of Jacobian calculation*/
554  } // End of if not boundary condition
555  } // End of loop over master nodes
556  } // End of loop over l
557  } // end of integration loop
558 
559  // Now add boundary contributions from Robin BCs
560  unsigned nrobin =
561  this->Pressure_advection_diffusion_robin_element_pt.size();
562  for (unsigned e = 0; e < nrobin; e++)
563  {
564  this->Pressure_advection_diffusion_robin_element_pt[e]
565  ->fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
566  residuals, jacobian, flag);
567  }
568  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
e
Array< double, 1, 3 > e(1./3., 0.5, 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::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::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::dshape_eulerian_at_knot
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Definition: elements.cc:3325
oomph::GeneralisedElement::ndof
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
oomph::GeneralisedElement::eqn_number
unsigned long eqn_number(const unsigned &ieqn_local) const
Definition: elements.h:704
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::NavierStokesEquations::dpshape_and_dptest_eulerian_nst
virtual double dpshape_and_dptest_eulerian_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const =0
oomph::NavierStokesEquations::p_nodal_index_nst
virtual int p_nodal_index_nst() const
Definition: navier_stokes_elements.h:936
oomph::NavierStokesEquations::npres_nst
virtual unsigned npres_nst() const =0
Function to return number of pressure degrees of freedom.
oomph::NavierStokesEquations::Press_adv_diff_source_fct_pt
NavierStokesPressureAdvDiffSourceFctPt Press_adv_diff_source_fct_pt
Definition: navier_stokes_elements.h:465
oomph::NavierStokesEquations::Pinned_fp_pressure_eqn
int Pinned_fp_pressure_eqn
Definition: navier_stokes_elements.h:479
oomph::NavierStokesEquations::re
const double & re() const
Reynolds number.
Definition: navier_stokes_elements.h:703
oomph::NavierStokesEquations::p_nst
virtual double p_nst(const unsigned &n_p) const =0
oomph::NavierStokesEquations::density_ratio
const double & density_ratio() const
Definition: navier_stokes_elements.h:741
oomph::NavierStokesEquations::Pressure_advection_diffusion_robin_element_pt
Vector< FpPressureAdvDiffRobinBCElementBase * > Pressure_advection_diffusion_robin_element_pt
Definition: navier_stokes_elements.h:475
oomph::RefineableElement::local_hang_eqn
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Definition: refineable_elements.h:278
oomph::RefineableNavierStokesEquations::pressure_node_pt
virtual Node * pressure_node_pt(const unsigned &n_p)
Definition: refineable_navier_stokes_elements.h:379
oomph::TemplateFreeNavierStokesEquationsBase::p_local_eqn
virtual int p_local_eqn(const unsigned &n) const =0
s
RealScalar s
Definition: level1_cplx_impl.h:130
int
return int(ret)+1
k
char char char int int * k
Definition: level2_impl.h:374
DIM
#define DIM
Definition: linearised_navier_stokes_elements.h:44
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_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, e(), i, int(), J, k, m, m2(), oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::HangInfo::nmaster(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, TestProblem::source(), w, and oomph::QuadTreeNames::W.

◆ fill_in_generic_residual_contribution_nst()

template<unsigned DIM>
void oomph::RefineableNavierStokesEquations< 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::NavierStokesEquations< DIM >.

583  {
584  // Find out how many nodes there are
585  unsigned n_node = nnode();
586 
587  // Get continuous time from timestepper of first node
588  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
589 
590  // Find out how many pressure dofs there are
591  unsigned n_pres = this->npres_nst();
592 
593  // Get the indices at which the velocity components are stored
594  unsigned u_nodal_index[DIM];
595  for (unsigned i = 0; i < DIM; i++)
596  {
597  u_nodal_index[i] = this->u_index_nst(i);
598  }
599 
600  // Which nodal value represents the pressure? (Negative if pressure
601  // is not based on nodal interpolation).
602  int p_index = this->p_nodal_index_nst();
603 
604  // Local array of booleans that are true if the l-th pressure value is
605  // hanging (avoid repeated virtual function calls)
606  bool pressure_dof_is_hanging[n_pres];
607  // If the pressure is stored at a node
608  if (p_index >= 0)
609  {
610  // Read out whether the pressure is hanging
611  for (unsigned l = 0; l < n_pres; ++l)
612  {
613  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
614  }
615  }
616  // Otherwise the pressure is not stored at a node and so cannot hang
617  else
618  {
619  for (unsigned l = 0; l < n_pres; ++l)
620  {
621  pressure_dof_is_hanging[l] = false;
622  }
623  }
624 
625  // Set up memory for the shape and test functions
626  Shape psif(n_node), testf(n_node);
627  DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
628 
629 
630  // Set up memory for pressure shape and test functions
631  Shape psip(n_pres), testp(n_pres);
632 
633  // Set the value of n_intpt
634  unsigned n_intpt = integral_pt()->nweight();
635 
636  // Set the Vector to hold local coordinates
637  Vector<double> s(DIM);
638 
639  // Get Physical Variables from Element
640  // Reynolds number must be multiplied by the density ratio
641  double scaled_re = this->re() * this->density_ratio();
642  double scaled_re_st = this->re_st() * this->density_ratio();
643  double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
644  double visc_ratio = this->viscosity_ratio();
645  Vector<double> G = this->g();
646 
647  // Integers that store the local equations and unknowns
648  int local_eqn = 0, local_unknown = 0;
649 
650  // Pointers to hang info objects
651  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
652 
653  // Local boolean for ALE (or not)
654  bool ALE_is_disabled_flag = this->ALE_is_disabled;
655 
656  // Loop over the integration points
657  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
658  {
659  // Assign values of s
660  for (unsigned i = 0; i < DIM; i++)
661  {
662  s[i] = integral_pt()->knot(ipt, i);
663  }
664 
665  // Get the integral weight
666  double w = integral_pt()->weight(ipt);
667 
668  // Call the derivatives of the shape and test functions
669  double J = this->dshape_and_dtest_eulerian_at_knot_nst(
670  ipt, psif, dpsifdx, testf, dtestfdx);
671 
672  // Call the pressure shape and test functions
673  this->pshape_nst(s, psip, testp);
674 
675  // Premultiply the weights and the Jacobian
676  double W = w * J;
677 
678  // Calculate local values of the pressure and velocity components
679  //--------------------------------------------------------------
680  double interpolated_p = 0.0;
681  Vector<double> interpolated_u(DIM, 0.0);
682  Vector<double> interpolated_x(DIM, 0.0);
683  Vector<double> mesh_veloc(DIM, 0.0);
684  Vector<double> dudt(DIM, 0.0);
685  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
686 
687  // Calculate pressure
688  for (unsigned l = 0; l < n_pres; l++)
689  {
690  interpolated_p += this->p_nst(l) * psip[l];
691  }
692 
693 
694  // Calculate velocities and derivatives
695 
696  // Loop over nodes
697  for (unsigned l = 0; l < n_node; l++)
698  {
699  // Loop over directions
700  for (unsigned i = 0; i < DIM; i++)
701  {
702  // Get the nodal value
703  double u_value = this->nodal_value(l, u_nodal_index[i]);
704  interpolated_u[i] += u_value * psif[l];
705  interpolated_x[i] += this->nodal_position(l, i) * psif[l];
706  dudt[i] += this->du_dt_nst(l, i) * psif[l];
707 
708  // Loop over derivative directions for velocity gradients
709  for (unsigned j = 0; j < DIM; j++)
710  {
711  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
712  }
713  }
714  }
715 
716  if (!ALE_is_disabled_flag)
717  {
718  // Loop over nodes
719  for (unsigned l = 0; l < n_node; l++)
720  {
721  // Loop over directions
722  for (unsigned i = 0; i < DIM; i++)
723  {
724  mesh_veloc[i] += this->dnodal_position_dt(l, i) * psif[l];
725  }
726  }
727  }
728 
729  // Get the user-defined body force terms
730  Vector<double> body_force(DIM);
731  this->get_body_force_nst(time, ipt, s, interpolated_x, body_force);
732 
733  // Get the user-defined source function
734  double source = this->get_source_nst(time, ipt, interpolated_x);
735 
736  // MOMENTUM EQUATIONS
737  //==================
738 
739  // Number of master nodes and storage for the weight of the shape function
740  unsigned n_master = 1;
741  double hang_weight = 1.0;
742 
743  // Loop over the nodes for the test functions/equations
744  //----------------------------------------------------
745  for (unsigned l = 0; l < n_node; l++)
746  {
747  // Local boolean to indicate whether the node is hanging
748  bool is_node_hanging = node_pt(l)->is_hanging();
749 
750  // If the node is hanging
751  if (is_node_hanging)
752  {
753  hang_info_pt = node_pt(l)->hanging_pt();
754 
755  // Read out number of master nodes from hanging data
756  n_master = hang_info_pt->nmaster();
757  }
758  // Otherwise the node is its own master
759  else
760  {
761  n_master = 1;
762  }
763 
764  // Loop over the master nodes
765  for (unsigned m = 0; m < n_master; m++)
766  {
767  // Loop over velocity components for equations
768  for (unsigned i = 0; i < DIM; i++)
769  {
770  // Get the equation number
771  // If the node is hanging
772  if (is_node_hanging)
773  {
774  // Get the equation number from the master node
775  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
776  u_nodal_index[i]);
777  // Get the hang weight from the master node
778  hang_weight = hang_info_pt->master_weight(m);
779  }
780  // If the node is not hanging
781  else
782  {
783  // Local equation number
784  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
785 
786  // Node contributes with full weight
787  hang_weight = 1.0;
788  }
789 
790  // If it's not a boundary condition...
791  if (local_eqn >= 0)
792  {
793  // Temporary variable to hold the residuals
794  double sum = 0.0;
795 
796  // Add the user-defined body force terms
797  sum += body_force[i];
798 
799  // Add the gravitational body force term
800  sum += scaled_re_inv_fr * G[i];
801 
802  // Add in the inertial term
803  sum -= scaled_re_st * dudt[i];
804 
805  // Convective terms, including mesh velocity
806  for (unsigned k = 0; k < DIM; k++)
807  {
808  double tmp = scaled_re * interpolated_u[k];
809  if (!ALE_is_disabled_flag)
810  {
811  tmp -= scaled_re_st * mesh_veloc[k];
812  }
813  sum -= tmp * interpolated_dudx(i, k);
814  }
815 
816  // Add the pressure gradient term
817  sum = (sum * testf[l] + interpolated_p * dtestfdx(l, i)) * W *
818  hang_weight;
819 
820  // Add in the stress tensor terms
821  // The viscosity ratio needs to go in here to ensure
822  // continuity of normal stress is satisfied even in flows
823  // with zero pressure gradient!
824  for (unsigned k = 0; k < DIM; k++)
825  {
826  sum -= visc_ratio *
827  (interpolated_dudx(i, k) +
828  this->Gamma[i] * interpolated_dudx(k, i)) *
829  dtestfdx(l, k) * W * hang_weight;
830  }
831 
832  // Add contribution
833  residuals[local_eqn] += sum;
834 
835  // CALCULATE THE JACOBIAN
836  if (flag)
837  {
838  // Number of master nodes and weights
839  unsigned n_master2 = 1;
840  double hang_weight2 = 1.0;
841  // Loop over the velocity nodes for columns
842  for (unsigned l2 = 0; l2 < n_node; l2++)
843  {
844  // Local boolean to indicate whether the node is hanging
845  bool is_node2_hanging = node_pt(l2)->is_hanging();
846 
847  // If the node is hanging
848  if (is_node2_hanging)
849  {
850  hang_info2_pt = node_pt(l2)->hanging_pt();
851  // Read out number of master nodes from hanging data
852  n_master2 = hang_info2_pt->nmaster();
853  }
854  // Otherwise the node is its own master
855  else
856  {
857  n_master2 = 1;
858  }
859 
860  // Loop over the master nodes
861  for (unsigned m2 = 0; m2 < n_master2; m2++)
862  {
863  // Loop over the velocity components
864  for (unsigned i2 = 0; i2 < DIM; i2++)
865  {
866  // Get the number of the unknown
867  // If the node is hanging
868  if (is_node2_hanging)
869  {
870  // Get the equation number from the master node
871  local_unknown = this->local_hang_eqn(
872  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
873  // Get the hang weights
874  hang_weight2 = hang_info2_pt->master_weight(m2);
875  }
876  else
877  {
878  local_unknown =
879  this->nodal_local_eqn(l2, u_nodal_index[i2]);
880  hang_weight2 = 1.0;
881  }
882 
883  // If the unknown is non-zero
884  if (local_unknown >= 0)
885  {
886  // Add contribution to Elemental Matrix
887  jacobian(local_eqn, local_unknown) -=
888  visc_ratio * this->Gamma[i] * dpsifdx(l2, i) *
889  dtestfdx(l, i2) * W * hang_weight * hang_weight2;
890 
891  // Now add in the inertial terms
892  jacobian(local_eqn, local_unknown) -=
893  scaled_re * psif[l2] * interpolated_dudx(i, i2) *
894  testf[l] * W * hang_weight * hang_weight2;
895 
896  // Extra diagonal components if i2=i
897  if (i2 == i)
898  {
899  // Mass matrix entries
900  // Again note the positive sign because the mass
901  // matrix is taken on the other side of the equation
902  if (flag == 2)
903  {
904  mass_matrix(local_eqn, local_unknown) +=
905  scaled_re_st * psif[l2] * testf[l] * W *
906  hang_weight * hang_weight2;
907  }
908 
909  // du/dt term
910  jacobian(local_eqn, local_unknown) -=
911  scaled_re_st *
912  node_pt(l2)->time_stepper_pt()->weight(1, 0) *
913  psif[l2] * testf[l] * W * hang_weight *
914  hang_weight2;
915 
916  // Extra advective terms
917  for (unsigned k = 0; k < DIM; k++)
918  {
919  double tmp = scaled_re * interpolated_u[k];
920  if (!ALE_is_disabled_flag)
921  {
922  tmp -= scaled_re_st * mesh_veloc[k];
923  }
924 
925  jacobian(local_eqn, local_unknown) -=
926  tmp * dpsifdx(l2, k) * testf[l] * W *
927  hang_weight * hang_weight2;
928  }
929 
930  // Extra viscous terms
931  for (unsigned k = 0; k < DIM; k++)
932  {
933  jacobian(local_eqn, local_unknown) -=
934  visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) * W *
935  hang_weight * hang_weight2;
936  }
937  }
938  }
939  }
940  }
941  }
942 
943  // Loop over the pressure shape functions
944  for (unsigned l2 = 0; l2 < n_pres; l2++)
945  {
946  // If the pressure dof is hanging
947  if (pressure_dof_is_hanging[l2])
948  {
949  hang_info2_pt =
950  this->pressure_node_pt(l2)->hanging_pt(p_index);
951  // Pressure dof is hanging so it must be nodal-based
952  // Get the number of master nodes from the pressure node
953  n_master2 = hang_info2_pt->nmaster();
954  }
955  // Otherwise the node is its own master
956  else
957  {
958  n_master2 = 1;
959  }
960 
961  // Loop over the master nodes
962  for (unsigned m2 = 0; m2 < n_master2; m2++)
963  {
964  // Get the number of the unknown
965  // If the pressure dof is hanging
966  if (pressure_dof_is_hanging[l2])
967  {
968  // Get the unknown from the master node
969  local_unknown = this->local_hang_eqn(
970  hang_info2_pt->master_node_pt(m2), p_index);
971  // Get the weight from the hanging object
972  hang_weight2 = hang_info2_pt->master_weight(m2);
973  }
974  else
975  {
976  local_unknown = this->p_local_eqn(l2);
977  hang_weight2 = 1.0;
978  }
979 
980  // If the unknown is not pinned
981  if (local_unknown >= 0)
982  {
983  jacobian(local_eqn, local_unknown) +=
984  psip[l2] * dtestfdx(l, i) * W * hang_weight *
985  hang_weight2;
986  }
987  }
988  }
989 
990  } // End of Jacobian calculation
991 
992  } // End of if not boundary condition statement
993 
994  } // End of loop over components of non-hanging node
995 
996  } // End of loop over master nodes
997 
998  } // End of loop over nodes for equations
999 
1000 
1001  // CONTINUITY EQUATION
1002  //===================
1003 
1004  // Loop over the pressure shape functions
1005  for (unsigned l = 0; l < n_pres; l++)
1006  {
1007  // If the pressure dof is hanging
1008  if (pressure_dof_is_hanging[l])
1009  {
1010  // Pressure dof is hanging so it must be nodal-based
1011  // Get the hang info object
1012  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1013 
1014  // Get the number of master nodes from the pressure node
1015  n_master = hang_info_pt->nmaster();
1016  }
1017  // Otherwise the node is its own master
1018  else
1019  {
1020  n_master = 1;
1021  }
1022 
1023  // Loop over the master nodes
1024  for (unsigned m = 0; m < n_master; m++)
1025  {
1026  // Get the number of the unknown
1027  // If the pressure dof is hanging
1028  if (pressure_dof_is_hanging[l])
1029  {
1030  // Get the local equation from the master node
1031  local_eqn =
1032  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1033  // Get the weight for the node
1034  hang_weight = hang_info_pt->master_weight(m);
1035  }
1036  else
1037  {
1038  local_eqn = this->p_local_eqn(l);
1039  hang_weight = 1.0;
1040  }
1041 
1042  // If the equation is not pinned
1043  if (local_eqn >= 0)
1044  {
1045  // Source term
1046  residuals[local_eqn] -= source * testp[l] * W * hang_weight;
1047 
1048  // Loop over velocity components
1049  for (unsigned k = 0; k < DIM; k++)
1050  {
1051  residuals[local_eqn] +=
1052  interpolated_dudx(k, k) * testp[l] * W * hang_weight;
1053  }
1054 
1055  // CALCULATE THE JACOBIAN
1056  if (flag)
1057  {
1058  unsigned n_master2 = 1;
1059  double hang_weight2 = 1.0;
1060  // Loop over the velocity nodes for columns
1061  for (unsigned l2 = 0; l2 < n_node; l2++)
1062  {
1063  // Local boolean to indicate whether the node is hanging
1064  bool is_node2_hanging = node_pt(l2)->is_hanging();
1065 
1066  // If the node is hanging
1067  if (is_node2_hanging)
1068  {
1069  hang_info2_pt = node_pt(l2)->hanging_pt();
1070  // Read out number of master nodes from hanging data
1071  n_master2 = hang_info2_pt->nmaster();
1072  }
1073  // Otherwise the node is its own master
1074  else
1075  {
1076  n_master2 = 1;
1077  }
1078 
1079  // Loop over the master nodes
1080  for (unsigned m2 = 0; m2 < n_master2; m2++)
1081  {
1082  // Loop over the velocity components
1083  for (unsigned i2 = 0; i2 < DIM; i2++)
1084  {
1085  // Get the number of the unknown
1086  // If the node is hanging
1087  if (is_node2_hanging)
1088  {
1089  // Get the equation number from the master node
1090  local_unknown = this->local_hang_eqn(
1091  hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
1092  hang_weight2 = hang_info2_pt->master_weight(m2);
1093  }
1094  else
1095  {
1096  local_unknown =
1097  this->nodal_local_eqn(l2, u_nodal_index[i2]);
1098  hang_weight2 = 1.0;
1099  }
1100 
1101  // If the unknown is not pinned
1102  if (local_unknown >= 0)
1103  {
1104  jacobian(local_eqn, local_unknown) +=
1105  dpsifdx(l2, i2) * testp[l] * W * hang_weight *
1106  hang_weight2;
1107  }
1108  }
1109  }
1110  }
1111 
1112  // NO PRESSURE CONTRIBUTION TO THE JACOBIAN
1113 
1114  } // End of jacobian calculation
1115  }
1116  }
1117  } // End of loop over pressure variables
1118 
1119  } // End of loop over integration points
1120  }
G
JacobiRotation< float > G
Definition: Jacobi_makeGivens.cpp:2
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::FiniteElement::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
oomph::FiniteElement::dnodal_position_dt
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2333
oomph::HangInfo::nmaster
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
oomph::NavierStokesEquations::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::NavierStokesEquations::pshape_nst
virtual void pshape_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
oomph::NavierStokesEquations::re_st
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition: navier_stokes_elements.h:709
oomph::NavierStokesEquations::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: navier_stokes_elements.h:517
oomph::NavierStokesEquations::re_invfr
const double & re_invfr() const
Global inverse Froude number.
Definition: navier_stokes_elements.h:753
oomph::NavierStokesEquations::viscosity_ratio
const double & viscosity_ratio() const
Definition: navier_stokes_elements.h:728
oomph::NavierStokesEquations::get_source_nst
virtual double get_source_nst(const double &time, const unsigned &ipt, const Vector< double > &x)
Definition: navier_stokes_elements.h:585
oomph::NavierStokesEquations::du_dt_nst
double du_dt_nst(const unsigned &n, const unsigned &i) const
Definition: navier_stokes_elements.h:880
oomph::NavierStokesEquations::ALE_is_disabled
bool ALE_is_disabled
Definition: navier_stokes_elements.h:470
oomph::NavierStokesEquations::g
const Vector< double > & g() const
Vector of gravitational components.
Definition: navier_stokes_elements.h:765
oomph::NavierStokesEquations::Gamma
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition: navier_stokes_elements.h:698
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
tmp
Eigen::Matrix< Scalar, Dynamic, Dynamic, ColMajor > tmp
Definition: level3_impl.h:365
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
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

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

◆ further_build()

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

Further build, pass the pointers down to the sons.

Reimplemented from oomph::RefineableElement.

Reimplemented in oomph::RefineableQCrouzeixRaviartElement< DIM >, oomph::RefineableQCrouzeixRaviartElement< 2 >, oomph::RefineableQCrouzeixRaviartElement< DIM >, oomph::RefineableQCrouzeixRaviartElement< 2 >, oomph::PRefineableQCrouzeixRaviartElement< DIM >, oomph::RefineableQCrouzeixRaviartElement< DIM >, oomph::RefineableQCrouzeixRaviartElement< 2 >, oomph::PRefineableQCrouzeixRaviartElement< DIM >, and oomph::PRefineableQCrouzeixRaviartElement< DIM >.

446  {
447  // Find the father element
448  RefineableNavierStokesEquations<DIM>* cast_father_element_pt =
449  dynamic_cast<RefineableNavierStokesEquations<DIM>*>(
450  this->father_element_pt());
451 
452  // Set the viscosity ratio pointer
453  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
454  // Set the density ratio pointer
455  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
456  // Set pointer to global Reynolds number
457  this->Re_pt = cast_father_element_pt->re_pt();
458  // Set pointer to global Reynolds number x Strouhal number (=Womersley)
459  this->ReSt_pt = cast_father_element_pt->re_st_pt();
460  // Set pointer to global Reynolds number x inverse Froude number
461  this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
462  // Set pointer to global gravity Vector
463  this->G_pt = cast_father_element_pt->g_pt();
464 
465  // Set pointer to body force function
466  this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
467 
468  // Set pointer to volumetric source function
469  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
470 
471  // Set the ALE flag
472  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
473  }
oomph::NavierStokesEquations::Viscosity_Ratio_pt
double * Viscosity_Ratio_pt
Definition: navier_stokes_elements.h:436
oomph::NavierStokesEquations::Source_fct_pt
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
Definition: navier_stokes_elements.h:461
oomph::NavierStokesEquations::Re_pt
double * Re_pt
Pointer to global Reynolds number.
Definition: navier_stokes_elements.h:445
oomph::NavierStokesEquations::G_pt
Vector< double > * G_pt
Pointer to global gravity Vector.
Definition: navier_stokes_elements.h:455
oomph::NavierStokesEquations::ReInvFr_pt
double * ReInvFr_pt
Definition: navier_stokes_elements.h:452
oomph::NavierStokesEquations::Density_Ratio_pt
double * Density_Ratio_pt
Definition: navier_stokes_elements.h:440
oomph::NavierStokesEquations::Body_force_fct_pt
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
Definition: navier_stokes_elements.h:458
oomph::NavierStokesEquations::ReSt_pt
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Definition: navier_stokes_elements.h:448
oomph::RefineableElement::father_element_pt
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Definition: refineable_elements.h:539

References oomph::NavierStokesEquations< DIM >::ALE_is_disabled, oomph::NavierStokesEquations< DIM >::Body_force_fct_pt, oomph::NavierStokesEquations< DIM >::body_force_fct_pt(), oomph::NavierStokesEquations< DIM >::Density_Ratio_pt, oomph::NavierStokesEquations< DIM >::density_ratio_pt(), oomph::RefineableElement::father_element_pt(), oomph::NavierStokesEquations< DIM >::G_pt, oomph::NavierStokesEquations< DIM >::g_pt(), oomph::NavierStokesEquations< DIM >::re_invfr_pt(), oomph::NavierStokesEquations< DIM >::Re_pt, oomph::NavierStokesEquations< DIM >::re_pt(), oomph::NavierStokesEquations< DIM >::re_st_pt(), oomph::NavierStokesEquations< DIM >::ReInvFr_pt, oomph::NavierStokesEquations< DIM >::ReSt_pt, oomph::NavierStokesEquations< DIM >::Source_fct_pt, oomph::NavierStokesEquations< DIM >::source_fct_pt(), oomph::NavierStokesEquations< DIM >::Viscosity_Ratio_pt, and oomph::NavierStokesEquations< DIM >::viscosity_ratio_pt().

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

◆ get_dresidual_dnodal_coordinates()

template<unsigned DIM>
void oomph::RefineableNavierStokesEquations< 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::RefineableElement.

1132  {
1133  // Return immediately if there are no dofs
1134  if (ndof() == 0)
1135  {
1136  return;
1137  }
1138 
1139  // Determine number of nodes in element
1140  const unsigned n_node = nnode();
1141 
1142  // Get continuous time from timestepper of first node
1143  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1144 
1145  // Determine number of pressure dofs in element
1146  const unsigned n_pres = this->npres_nst();
1147 
1148  // Find the indices at which the local velocities are stored
1149  unsigned u_nodal_index[DIM];
1150  for (unsigned i = 0; i < DIM; i++)
1151  {
1152  u_nodal_index[i] = this->u_index_nst(i);
1153  }
1154 
1155  // Which nodal value represents the pressure? (Negative if pressure
1156  // is not based on nodal interpolation).
1157  const int p_index = this->p_nodal_index_nst();
1158 
1159  // Local array of booleans that are true if the l-th pressure value is
1160  // hanging (avoid repeated virtual function calls)
1161  bool pressure_dof_is_hanging[n_pres];
1162 
1163  // If the pressure is stored at a node
1164  if (p_index >= 0)
1165  {
1166  // Read out whether the pressure is hanging
1167  for (unsigned l = 0; l < n_pres; ++l)
1168  {
1169  pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
1170  }
1171  }
1172  // Otherwise the pressure is not stored at a node and so cannot hang
1173  else
1174  {
1175  for (unsigned l = 0; l < n_pres; ++l)
1176  {
1177  pressure_dof_is_hanging[l] = false;
1178  }
1179  }
1180 
1181  // Set up memory for the shape and test functions
1182  Shape psif(n_node), testf(n_node);
1183  DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
1184 
1185  // Set up memory for pressure shape and test functions
1186  Shape psip(n_pres), testp(n_pres);
1187 
1188  // Determine number of shape controlling nodes
1189  const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1190 
1191  // Deriatives of shape fct derivatives w.r.t. nodal coords
1192  RankFourTensor<double> d_dpsifdx_dX(
1193  DIM, n_shape_controlling_node, n_node, DIM);
1194  RankFourTensor<double> d_dtestfdx_dX(
1195  DIM, n_shape_controlling_node, n_node, DIM);
1196 
1197  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1198  DenseMatrix<double> dJ_dX(DIM, n_shape_controlling_node);
1199 
1200  // Derivatives of derivative of u w.r.t. nodal coords
1201  RankFourTensor<double> d_dudx_dX(DIM, n_shape_controlling_node, DIM, DIM);
1202 
1203  // Derivatives of nodal velocities w.r.t. nodal coords:
1204  // Assumption: Interaction only local via no-slip so
1205  // X_ij only affects U_ij.
1206  DenseMatrix<double> d_U_dX(DIM, n_shape_controlling_node, 0.0);
1207 
1208  // Determine the number of integration points
1209  const unsigned n_intpt = integral_pt()->nweight();
1210 
1211  // Vector to hold local coordinates
1212  Vector<double> s(DIM);
1213 
1214  // Get physical variables from element
1215  // (Reynolds number must be multiplied by the density ratio)
1216  double scaled_re = this->re() * this->density_ratio();
1217  double scaled_re_st = this->re_st() * this->density_ratio();
1218  double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
1219  double visc_ratio = this->viscosity_ratio();
1220  Vector<double> G = this->g();
1221 
1222  // FD step
1223  double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
1224 
1225  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1226  // Assumption: Interaction only local via no-slip so
1227  // X_ij only affects U_ij.
1228  bool element_has_node_with_aux_node_update_fct = false;
1229 
1230  std::map<Node*, unsigned> local_shape_controlling_node_lookup =
1231  shape_controlling_node_lookup();
1232 
1233  // FD loop over shape-controlling nodes
1234  for (std::map<Node*, unsigned>::iterator it =
1235  local_shape_controlling_node_lookup.begin();
1236  it != local_shape_controlling_node_lookup.end();
1237  it++)
1238  {
1239  // Get node
1240  Node* nod_pt = it->first;
1241 
1242  // Get its number
1243  unsigned q = it->second;
1244 
1245  // Only compute if there's a node-update fct involved
1246  if (nod_pt->has_auxiliary_node_update_fct_pt())
1247  {
1248  element_has_node_with_aux_node_update_fct = true;
1249 
1250  // Current nodal velocity
1251  Vector<double> u_ref(DIM);
1252  for (unsigned i = 0; i < DIM; i++)
1253  {
1254  u_ref[i] = *(nod_pt->value_pt(u_nodal_index[i]));
1255  }
1256 
1257  // FD
1258  for (unsigned p = 0; p < DIM; p++)
1259  {
1260  // Make backup
1261  double backup = nod_pt->x(p);
1262 
1263  // Do FD step. No node update required as we're
1264  // attacking the coordinate directly...
1265  nod_pt->x(p) += eps_fd;
1266 
1267  // Do auxiliary node update (to apply no slip)
1268  nod_pt->perform_auxiliary_node_update_fct();
1269 
1270  // Evaluate
1271  d_U_dX(p, q) =
1272  (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
1273 
1274  // Reset
1275  nod_pt->x(p) = backup;
1276 
1277  // Do auxiliary node update (to apply no slip)
1278  nod_pt->perform_auxiliary_node_update_fct();
1279  }
1280  }
1281  }
1282 
1283  // Integer to store the local equation number
1284  int local_eqn = 0;
1285 
1286  // Pointers to hang info object
1287  HangInfo* hang_info_pt = 0;
1288 
1289  // Loop over the integration points
1290  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1291  {
1292  // Assign values of s
1293  for (unsigned i = 0; i < DIM; i++)
1294  {
1295  s[i] = integral_pt()->knot(ipt, i);
1296  }
1297 
1298  // Get the integral weight
1299  const double w = integral_pt()->weight(ipt);
1300 
1301  // Call the derivatives of the shape and test functions
1302  const double J =
1303  this->dshape_and_dtest_eulerian_at_knot_nst(ipt,
1304  psif,
1305  dpsifdx,
1306  d_dpsifdx_dX,
1307  testf,
1308  dtestfdx,
1309  d_dtestfdx_dX,
1310  dJ_dX);
1311 
1312  // Call the pressure shape and test functions
1313  this->pshape_nst(s, psip, testp);
1314 
1315  // Calculate local values of the pressure and velocity components
1316  // Allocate
1317  double interpolated_p = 0.0;
1318  Vector<double> interpolated_u(DIM, 0.0);
1319  Vector<double> interpolated_x(DIM, 0.0);
1320  Vector<double> mesh_velocity(DIM, 0.0);
1321  Vector<double> dudt(DIM, 0.0);
1322  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1323 
1324  // Calculate pressure
1325  for (unsigned l = 0; l < n_pres; l++)
1326  {
1327  interpolated_p += this->p_nst(l) * psip[l];
1328  }
1329 
1330  // Calculate velocities and derivatives:
1331 
1332  // Loop over nodes
1333  for (unsigned l = 0; l < n_node; l++)
1334  {
1335  // Loop over directions
1336  for (unsigned i = 0; i < DIM; i++)
1337  {
1338  // Get the nodal value
1339  const double u_value = nodal_value(l, u_nodal_index[i]);
1340  interpolated_u[i] += u_value * psif[l];
1341  interpolated_x[i] += nodal_position(l, i) * psif[l];
1342  dudt[i] += this->du_dt_nst(l, i) * psif[l];
1343 
1344  // Loop over derivative directions
1345  for (unsigned j = 0; j < DIM; j++)
1346  {
1347  interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
1348  }
1349  }
1350  }
1351 
1352  if (!this->ALE_is_disabled)
1353  {
1354  // Loop over nodes
1355  for (unsigned l = 0; l < n_node; l++)
1356  {
1357  // Loop over directions
1358  for (unsigned i = 0; i < DIM; i++)
1359  {
1360  mesh_velocity[i] += this->dnodal_position_dt(l, i) * psif[l];
1361  }
1362  }
1363  }
1364 
1365  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1366 
1367  // Loop over shape-controlling nodes
1368  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1369  {
1370  // Loop over coordinate directions
1371  for (unsigned p = 0; p < DIM; p++)
1372  {
1373  for (unsigned i = 0; i < DIM; i++)
1374  {
1375  for (unsigned k = 0; k < DIM; k++)
1376  {
1377  double aux = 0.0;
1378  for (unsigned j = 0; j < n_node; j++)
1379  {
1380  aux +=
1381  nodal_value(j, u_nodal_index[i]) * d_dpsifdx_dX(p, q, j, k);
1382  }
1383  d_dudx_dX(p, q, i, k) = aux;
1384  }
1385  }
1386  }
1387  }
1388 
1389  // Get weight of actual nodal position/value in computation of mesh
1390  // velocity from positional/value time stepper
1391  const double pos_time_weight =
1392  node_pt(0)->position_time_stepper_pt()->weight(1, 0);
1393  const double val_time_weight =
1394  node_pt(0)->time_stepper_pt()->weight(1, 0);
1395 
1396  // Get the user-defined body force terms
1397  Vector<double> body_force(DIM);
1398  this->get_body_force_nst(time, ipt, s, interpolated_x, body_force);
1399 
1400  // Get the user-defined source function
1401  const double source = this->get_source_nst(time, ipt, interpolated_x);
1402 
1403  // Get gradient of body force function
1404  DenseMatrix<double> d_body_force_dx(DIM, DIM, 0.0);
1405  this->get_body_force_gradient_nst(
1406  time, ipt, s, interpolated_x, d_body_force_dx);
1407 
1408  // Get gradient of source function
1409  Vector<double> source_gradient(DIM, 0.0);
1410  this->get_source_gradient_nst(time, ipt, interpolated_x, source_gradient);
1411 
1412 
1413  // Assemble shape derivatives
1414  //---------------------------
1415 
1416  // MOMENTUM EQUATIONS
1417  // ------------------
1418 
1419  // Number of master nodes and storage for the weight of the shape function
1420  unsigned n_master = 1;
1421  double hang_weight = 1.0;
1422 
1423  // Loop over the test functions
1424  for (unsigned l = 0; l < n_node; l++)
1425  {
1426  // Local boolean to indicate whether the node is hanging
1427  bool is_node_hanging = node_pt(l)->is_hanging();
1428 
1429  // If the node is hanging
1430  if (is_node_hanging)
1431  {
1432  hang_info_pt = node_pt(l)->hanging_pt();
1433 
1434  // Read out number of master nodes from hanging data
1435  n_master = hang_info_pt->nmaster();
1436  }
1437  // Otherwise the node is its own master
1438  else
1439  {
1440  n_master = 1;
1441  }
1442 
1443  // Loop over the master nodes
1444  for (unsigned m = 0; m < n_master; m++)
1445  {
1446  // Loop over coordinate directions
1447  for (unsigned i = 0; i < DIM; i++)
1448  {
1449  // Get the equation number
1450  // If the node is hanging
1451  if (is_node_hanging)
1452  {
1453  // Get the equation number from the master node
1454  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1455  u_nodal_index[i]);
1456  // Get the hang weight from the master node
1457  hang_weight = hang_info_pt->master_weight(m);
1458  }
1459  // If the node is not hanging
1460  else
1461  {
1462  // Local equation number
1463  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
1464 
1465  // Node contributes with full weight
1466  hang_weight = 1.0;
1467  }
1468 
1469  // IF it's not a boundary condition
1470  if (local_eqn >= 0)
1471  {
1472  // Loop over coordinate directions
1473  for (unsigned p = 0; p < DIM; p++)
1474  {
1475  // Loop over shape controlling nodes
1476  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1477  {
1478  // Residual x deriv of Jacobian
1479  // ----------------------------
1480 
1481  // Add the user-defined body force terms
1482  double sum = body_force[i] * testf[l];
1483 
1484  // Add the gravitational body force term
1485  sum += scaled_re_inv_fr * testf[l] * G[i];
1486 
1487  // Add the pressure gradient term
1488  sum += interpolated_p * dtestfdx(l, i);
1489 
1490  // Add in the stress tensor terms
1491  // The viscosity ratio needs to go in here to ensure
1492  // continuity of normal stress is satisfied even in flows
1493  // with zero pressure gradient!
1494  for (unsigned k = 0; k < DIM; k++)
1495  {
1496  sum -= visc_ratio *
1497  (interpolated_dudx(i, k) +
1498  this->Gamma[i] * interpolated_dudx(k, i)) *
1499  dtestfdx(l, k);
1500  }
1501 
1502  // Add in the inertial terms
1503 
1504  // du/dt term
1505  sum -= scaled_re_st * dudt[i] * testf[l];
1506 
1507  // Convective terms, including mesh velocity
1508  for (unsigned k = 0; k < DIM; k++)
1509  {
1510  double tmp = scaled_re * interpolated_u[k];
1511  if (!this->ALE_is_disabled)
1512  {
1513  tmp -= scaled_re_st * mesh_velocity[k];
1514  }
1515  sum -= tmp * interpolated_dudx(i, k) * testf[l];
1516  }
1517 
1518  // Multiply throsugh by deriv of Jacobian and integration
1519  // weight
1520  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1521  sum * dJ_dX(p, q) * w * hang_weight;
1522 
1523  // Derivative of residual x Jacobian
1524  // ---------------------------------
1525 
1526  // Body force
1527  sum = d_body_force_dx(i, p) * psif(q) * testf(l);
1528 
1529  // Pressure gradient term
1530  sum += interpolated_p * d_dtestfdx_dX(p, q, l, i);
1531 
1532  // Viscous term
1533  for (unsigned k = 0; k < DIM; k++)
1534  {
1535  sum -=
1536  visc_ratio * ((interpolated_dudx(i, k) +
1537  this->Gamma[i] * interpolated_dudx(k, i)) *
1538  d_dtestfdx_dX(p, q, l, k) +
1539  (d_dudx_dX(p, q, i, k) +
1540  this->Gamma[i] * d_dudx_dX(p, q, k, i)) *
1541  dtestfdx(l, k));
1542  }
1543 
1544  // Convective terms, including mesh velocity
1545  for (unsigned k = 0; k < DIM; k++)
1546  {
1547  double tmp = scaled_re * interpolated_u[k];
1548  if (!this->ALE_is_disabled)
1549  {
1550  tmp -= scaled_re_st * mesh_velocity[k];
1551  }
1552  sum -= tmp * d_dudx_dX(p, q, i, k) * testf(l);
1553  }
1554  if (!this->ALE_is_disabled)
1555  {
1556  sum += scaled_re_st * pos_time_weight * psif(q) *
1557  interpolated_dudx(i, p) * testf(l);
1558  }
1559 
1560  // Multiply through by Jacobian and integration weight
1561  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1562  sum * J * w * hang_weight;
1563 
1564  } // End of loop over shape controlling nodes q
1565  } // End of loop over coordinate directions p
1566 
1567 
1568  // Derivs w.r.t. to nodal velocities
1569  // ---------------------------------
1570  if (element_has_node_with_aux_node_update_fct)
1571  {
1572  // Loop over local nodes
1573  for (unsigned q_local = 0; q_local < n_node; q_local++)
1574  {
1575  // Number of master nodes and storage for the weight of
1576  // the shape function
1577  unsigned n_master2 = 1;
1578  double hang_weight2 = 1.0;
1579  HangInfo* hang_info2_pt = 0;
1580 
1581  // Local boolean to indicate whether the node is hanging
1582  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1583 
1584  Node* actual_shape_controlling_node_pt = node_pt(q_local);
1585 
1586  // If the node is hanging
1587  if (is_node_hanging2)
1588  {
1589  hang_info2_pt = node_pt(q_local)->hanging_pt();
1590 
1591  // Read out number of master nodes from hanging data
1592  n_master2 = hang_info2_pt->nmaster();
1593  }
1594  // Otherwise the node is its own master
1595  else
1596  {
1597  n_master2 = 1;
1598  }
1599 
1600  // Loop over the master nodes
1601  for (unsigned mm = 0; mm < n_master2; mm++)
1602  {
1603  if (is_node_hanging2)
1604  {
1605  actual_shape_controlling_node_pt =
1606  hang_info2_pt->master_node_pt(mm);
1607  hang_weight2 = hang_info2_pt->master_weight(mm);
1608  }
1609 
1610  // Find the corresponding number
1611  unsigned q = local_shape_controlling_node_lookup
1612  [actual_shape_controlling_node_pt];
1613 
1614  // Loop over coordinate directions
1615  for (unsigned p = 0; p < DIM; p++)
1616  {
1617  double sum = -visc_ratio * this->Gamma[i] *
1618  dpsifdx(q_local, i) * dtestfdx(l, p) -
1619  scaled_re * psif(q_local) *
1620  interpolated_dudx(i, p) * testf(l);
1621  if (i == p)
1622  {
1623  sum -= scaled_re_st * val_time_weight * psif(q_local) *
1624  testf(l);
1625  for (unsigned k = 0; k < DIM; k++)
1626  {
1627  sum -=
1628  visc_ratio * dpsifdx(q_local, k) * dtestfdx(l, k);
1629  double tmp = scaled_re * interpolated_u[k];
1630  if (!this->ALE_is_disabled)
1631  {
1632  tmp -= scaled_re_st * mesh_velocity[k];
1633  }
1634  sum -= tmp * dpsifdx(q_local, k) * testf(l);
1635  }
1636  }
1637 
1638  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1639  sum * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1640  }
1641  } // End of loop over master nodes
1642  } // End of loop over local nodes
1643  } // End of if(element_has_node_with_aux_node_update_fct)
1644 
1645 
1646  } // local_eqn>=0
1647  }
1648  }
1649  } // End of loop over test functions
1650 
1651 
1652  // CONTINUITY EQUATION
1653  // -------------------
1654 
1655  // Loop over the shape functions
1656  for (unsigned l = 0; l < n_pres; l++)
1657  {
1658  // If the pressure dof is hanging
1659  if (pressure_dof_is_hanging[l])
1660  {
1661  // Pressure dof is hanging so it must be nodal-based
1662  // Get the hang info object
1663  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1664 
1665  // Get the number of master nodes from the pressure node
1666  n_master = hang_info_pt->nmaster();
1667  }
1668  // Otherwise the node is its own master
1669  else
1670  {
1671  n_master = 1;
1672  }
1673 
1674  // Loop over the master nodes
1675  for (unsigned m = 0; m < n_master; m++)
1676  {
1677  // Get the number of the unknown
1678  // If the pressure dof is hanging
1679  if (pressure_dof_is_hanging[l])
1680  {
1681  // Get the local equation from the master node
1682  local_eqn =
1683  this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1684  // Get the weight for the node
1685  hang_weight = hang_info_pt->master_weight(m);
1686  }
1687  else
1688  {
1689  local_eqn = this->p_local_eqn(l);
1690  hang_weight = 1.0;
1691  }
1692 
1693  // If not a boundary conditions
1694  if (local_eqn >= 0)
1695  {
1696  // Loop over coordinate directions
1697  for (unsigned p = 0; p < DIM; p++)
1698  {
1699  // Loop over nodes
1700  for (unsigned q = 0; q < n_shape_controlling_node; q++)
1701  {
1702  // Residual x deriv of Jacobian
1703  //-----------------------------
1704 
1705  // Source term
1706  double aux = -source;
1707 
1708  // Loop over velocity components
1709  for (unsigned k = 0; k < DIM; k++)
1710  {
1711  aux += interpolated_dudx(k, k);
1712  }
1713 
1714  // Multiply through by deriv of Jacobian and integration weight
1715  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1716  aux * dJ_dX(p, q) * testp[l] * w * hang_weight;
1717 
1718 
1719  // Derivative of residual x Jacobian
1720  // ---------------------------------
1721 
1722  // Loop over velocity components
1723  aux = -source_gradient[p] * psif(q);
1724  for (unsigned k = 0; k < DIM; k++)
1725  {
1726  aux += d_dudx_dX(p, q, k, k);
1727  }
1728  // Multiply through by Jacobian and integration weight
1729  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1730  aux * testp[l] * J * w * hang_weight;
1731  }
1732  }
1733 
1734 
1735  // Derivs w.r.t. to nodal velocities
1736  // ---------------------------------
1737  if (element_has_node_with_aux_node_update_fct)
1738  {
1739  // Loop over local nodes
1740  for (unsigned q_local = 0; q_local < n_node; q_local++)
1741  {
1742  // Number of master nodes and storage for the weight of
1743  // the shape function
1744  unsigned n_master2 = 1;
1745  double hang_weight2 = 1.0;
1746  HangInfo* hang_info2_pt = 0;
1747 
1748  // Local boolean to indicate whether the node is hanging
1749  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1750 
1751  Node* actual_shape_controlling_node_pt = node_pt(q_local);
1752 
1753  // If the node is hanging
1754  if (is_node_hanging2)
1755  {
1756  hang_info2_pt = node_pt(q_local)->hanging_pt();
1757 
1758  // Read out number of master nodes from hanging data
1759  n_master2 = hang_info2_pt->nmaster();
1760  }
1761  // Otherwise the node is its own master
1762  else
1763  {
1764  n_master2 = 1;
1765  }
1766 
1767  // Loop over the master nodes
1768  for (unsigned mm = 0; mm < n_master2; mm++)
1769  {
1770  if (is_node_hanging2)
1771  {
1772  actual_shape_controlling_node_pt =
1773  hang_info2_pt->master_node_pt(mm);
1774  hang_weight2 = hang_info2_pt->master_weight(mm);
1775  }
1776 
1777  // Find the corresponding number
1778  unsigned q = local_shape_controlling_node_lookup
1779  [actual_shape_controlling_node_pt];
1780 
1781  // Loop over coordinate directions
1782  for (unsigned p = 0; p < DIM; p++)
1783  {
1784  double aux = dpsifdx(q_local, p) * testp(l);
1785  dresidual_dnodal_coordinates(local_eqn, p, q) +=
1786  aux * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1787  }
1788  } // End of loop over (mm) master nodes
1789  } // End of loop over local nodes q_local
1790  } // End of if(element_has_node_with_aux_node_update_fct)
1791  } // End of if it's not a boundary condition
1792  } // End of loop over (m) master nodes
1793  } // End of loop over shape functions for continuity eqn
1794 
1795  } // End of loop over integration points
1796  }
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::NavierStokesEquations::get_source_gradient_nst
virtual void get_source_gradient_nst(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
Definition: navier_stokes_elements.h:607
oomph::NavierStokesEquations::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: navier_stokes_elements.h:544
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::RefineableNavierStokesEquations< 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::NavierStokesEquations< 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::RefineableNavierStokesEquations< 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.

404  {
405 #ifdef PARANOID
406  unsigned num_entries = DIM + (DIM * (DIM - 1)) / 2;
407  if (flux.size() < num_entries)
408  {
409  std::ostringstream error_message;
410  error_message << "The flux vector has the wrong number of entries, "
411  << flux.size() << ", whereas it should be at least "
412  << num_entries << std::endl;
413  throw OomphLibError(error_message.str(),
414  OOMPH_CURRENT_FUNCTION,
415  OOMPH_EXCEPTION_LOCATION);
416  }
417 #endif
418 
419  // Get strain rate matrix
420  DenseMatrix<double> strainrate(DIM);
421  this->strain_rate(s, strainrate);
422 
423  // Pack into flux Vector
424  unsigned icount = 0;
425 
426  // Start with diagonal terms
427  for (unsigned i = 0; i < DIM; i++)
428  {
429  flux[icount] = strainrate(i, i);
430  icount++;
431  }
432 
433  // Off diagonals row by row
434  for (unsigned i = 0; i < DIM; i++)
435  {
436  for (unsigned j = i + 1; j < DIM; j++)
437  {
438  flux[icount] = strainrate(i, j);
439  icount++;
440  }
441  }
442  }
oomph::NavierStokesEquations::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: navier_stokes_elements.cc:1185
ProblemParameters::flux
void flux(const double &time, const Vector< double > &x, double &flux)
Get flux applied along boundary x=0.
Definition: pretend_melt.cc:59

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

Referenced by RefineableBuoyantQCrouzeixRaviartElement< DIM >::get_Z2_flux(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::get_Z2_flux(), and oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::get_Z2_flux().

◆ num_Z2_flux_terms()

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

Number of 'flux' terms for Z2 error estimation.

Implements oomph::ElementWithZ2ErrorEstimator.

396  {
397  // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
398  return DIM + (DIM * (DIM - 1)) / 2;
399  }

References DIM.

Referenced by RefineableBuoyantQCrouzeixRaviartElement< DIM >::get_Z2_compound_flux_indices(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::get_Z2_compound_flux_indices(), oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::get_Z2_compound_flux_indices(), RefineableBuoyantQCrouzeixRaviartElement< DIM >::get_Z2_flux(), oomph::RefineableDoubleBuoyantQCrouzeixRaviartElement< DIM >::get_Z2_flux(), and oomph::RefineableBuoyantQCrouzeixRaviartElement< DIM >::get_Z2_flux().

◆ pin_elemental_redundant_nodal_pressure_dofs()

template<unsigned DIM>
virtual void oomph::RefineableNavierStokesEquations< 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::RefineableQTaylorHoodElement< DIM >, and oomph::RefineableQTaylorHoodElement< 2 >.

329 {}

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

◆ pin_redundant_nodal_pressures()

template<unsigned DIM>
static void oomph::RefineableNavierStokesEquations< 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.)

353  {
354  // Loop over all elements and call the function that pins their
355  // unused nodal pressure data
356  unsigned n_element = element_pt.size();
357  for (unsigned e = 0; e < n_element; e++)
358  {
359  dynamic_cast<RefineableNavierStokesEquations<DIM>*>(element_pt[e])
360  ->pin_elemental_redundant_nodal_pressure_dofs();
361  }
362  }
oomph::RefineableNavierStokesEquations::pin_elemental_redundant_nodal_pressure_dofs
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Definition: refineable_navier_stokes_elements.h:329

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

Referenced by FlowAroundCylinderProblem< ELEMENT >::actions_after_adapt(), RefineablePorousChannelProblem< ELEMENT >::actions_after_adapt(), MovingBlockProblem< ELEMENT >::actions_after_adapt(), and RefineablePorousChannelProblem< ELEMENT >::RefineablePorousChannelProblem().

◆ pressure_node_pt()

template<unsigned DIM>
virtual Node* oomph::RefineableNavierStokesEquations< 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::RefineableQTaylorHoodElement< DIM >, and oomph::RefineableQTaylorHoodElement< 2 >.

380  {
381  return NULL;
382  }

◆ unpin_all_pressure_dofs()

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

Unpin all pressure dofs in elements listed in vector.

367  {
368  // Loop over all elements
369  unsigned n_element = element_pt.size();
370  for (unsigned e = 0; e < n_element; e++)
371  {
372  dynamic_cast<RefineableNavierStokesEquations<DIM>*>(element_pt[e])
373  ->unpin_elemental_pressure_dofs();
374  }
375  }
oomph::RefineableNavierStokesEquations::unpin_elemental_pressure_dofs
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.

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

Referenced by FlowAroundCylinderProblem< ELEMENT >::actions_after_adapt(), RefineablePorousChannelProblem< ELEMENT >::actions_after_adapt(), MovingBlockProblem< ELEMENT >::actions_after_adapt(), RefineableConvectionProblem< NST_ELEMENT, AD_ELEMENT >::actions_after_adapt(), and RefineableDDConvectionProblem< NST_ELEMENT, AD_ELEMENT >::actions_after_adapt().

◆ unpin_elemental_pressure_dofs()

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

Unpin all pressure dofs in the element.

Implemented in oomph::PRefineableQCrouzeixRaviartElement< DIM >, oomph::RefineableQCrouzeixRaviartElement< DIM >, oomph::RefineableQCrouzeixRaviartElement< 2 >, oomph::RefineableQTaylorHoodElement< DIM >, and oomph::RefineableQTaylorHoodElement< 2 >.

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

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