oomph::PVDEquationsWithPressure< DIM > Class Template Referenceabstract

#include <solid_elements.h>

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

Public Member Functions

 PVDEquationsWithPressure ()
 Constructor, by default the element is NOT incompressible. More...
 
void get_stress (const Vector< double > &s, DenseMatrix< double > &sigma)
 
bool is_incompressible () const
 Return whether the material is incompressible. More...
 
void set_incompressible ()
 Set the material to be incompressible. More...
 
void set_compressible ()
 Set the material to be compressible. More...
 
virtual double solid_p (const unsigned &l)=0
 Return the lth solid pressure. More...
 
virtual void set_solid_p (const unsigned &l, const double &p_value)=0
 Set the lth solid pressure to p_value. More...
 
void fill_in_contribution_to_residuals (Vector< double > &residuals)
 Fill in the residuals. 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)
 
double interpolated_solid_p (const Vector< double > &s)
 Return the interpolated_solid_pressure. More...
 
void output (std::ostream &outfile)
 Output: x,y,[z],xi0,xi1,[xi2],p,gamma. More...
 
void output (std::ostream &outfile, const unsigned &n_plot)
 Output: x,y,[z],xi0,xi1,[xi2],p,gamma. More...
 
void output (FILE *file_pt)
 C-style output: x,y,[z],xi0,xi1,[xi2],p,gamma. More...
 
void output (FILE *file_pt, const unsigned &n_plot)
 C-style output: x,y,[z],xi0,xi1,[xi2],p,gamma. More...
 
void extended_output (std::ostream &outfile, const unsigned &n_plot)
 Output: x,y,[z],xi0,xi1,[xi2],gamma strain and stress components. More...
 
void get_mass_matrix_diagonal (Vector< double > &mass_diag)
 
unsigned ndof_types () const
 
void get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
 
- Public Member Functions inherited from oomph::PVDEquationsBase< DIM >
 PVDEquationsBase ()
 
ConstitutiveLaw *& constitutive_law_pt ()
 Return the constitutive law pointer. More...
 
const doublelambda_sq () const
 Access function for timescale ratio (nondim density) More...
 
double *& lambda_sq_pt ()
 Access function for pointer to timescale ratio (nondim density) More...
 
IsotropicGrowthFctPtisotropic_growth_fct_pt ()
 Access function: Pointer to isotropic growth function. More...
 
PrestressFctPtprestress_fct_pt ()
 Access function: Pointer to pre-stress function. More...
 
IsotropicGrowthFctPt isotropic_growth_fct_pt () const
 Access function: Pointer to isotropic growth function (const version) More...
 
BodyForceFctPtbody_force_fct_pt ()
 Access function: Pointer to body force function. More...
 
BodyForceFctPt body_force_fct_pt () const
 Access function: Pointer to body force function (const version) More...
 
void enable_inertia ()
 Switch on solid inertia. More...
 
void disable_inertia ()
 Switch off solid inertia. More...
 
bool is_inertia_enabled () const
 Access function to flag that switches inertia on/off (const version) More...
 
virtual unsigned npres_solid () const
 
virtual int solid_p_local_eqn (const unsigned &i) const
 
virtual int solid_p_nodal_index () const
 
virtual void unpin_elemental_solid_pressure_dofs ()=0
 Unpin all solid pressure dofs in the element. More...
 
virtual void pin_elemental_redundant_nodal_solid_pressures ()
 Pin the element's redundant solid pressures (needed for refinement) More...
 
void get_strain (const Vector< double > &s, DenseMatrix< double > &strain) const
 Return the strain tensor. More...
 
void get_energy (double &pot_en, double &kin_en)
 Get potential (strain) and kinetic energy. More...
 
void get_deformed_covariant_basis_vectors (const Vector< double > &s, DenseMatrix< double > &def_covariant_basis)
 
void get_principal_stress (const Vector< double > &s, DenseMatrix< double > &principal_stress_vector, Vector< double > &principal_stress)
 
virtual void get_isotropic_growth (const unsigned &ipt, const Vector< double > &s, const Vector< double > &xi, double &gamma) const
 
void body_force (const Vector< double > &xi, Vector< double > &b) const
 
void enable_evaluate_jacobian_by_fd ()
 Set Jacobian to be evaluated by FD? Else: Analytically. More...
 
void disable_evaluate_jacobian_by_fd ()
 Set Jacobian to be evaluated analytically Else: by FD. More...
 
bool is_jacobian_evaluated_by_fd () const
 Return the flag indicating whether the jacobian is evaluated by fd. More...
 
double prestress (const unsigned &i, const unsigned &j, const Vector< double > xi)
 
- Public Member Functions inherited from oomph::SolidFiniteElement
void set_lagrangian_dimension (const unsigned &lagrangian_dimension)
 
virtual bool has_internal_solid_data ()
 
 SolidFiniteElement ()
 Constructor: Set defaults. More...
 
virtual ~SolidFiniteElement ()
 Destructor to clean up any allocated memory. More...
 
 SolidFiniteElement (const SolidFiniteElement &)=delete
 Broken copy constructor. More...
 
unsigned ngeom_data () const
 Broken assignment operator. More...
 
Datageom_data_pt (const unsigned &j)
 
void identify_geometric_data (std::set< Data * > &geometric_data_pt)
 
double zeta_nodal (const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual void get_x_and_xi (const Vector< double > &s, Vector< double > &x_fe, Vector< double > &x, Vector< double > &xi_fe, Vector< double > &xi) const
 
virtual void set_macro_elem_pt (MacroElement *macro_elem_pt)
 
virtual void set_macro_elem_pt (MacroElement *macro_elem_pt, MacroElement *undeformed_macro_elem_pt)
 
void set_undeformed_macro_elem_pt (MacroElement *undeformed_macro_elem_pt)
 
MacroElementundeformed_macro_elem_pt ()
 Access function to pointer to "undeformed" macro element. More...
 
double dshape_lagrangian (const Vector< double > &s, Shape &psi, DShape &dpsidxi) const
 
virtual double dshape_lagrangian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidxi) const
 
double d2shape_lagrangian (const Vector< double > &s, Shape &psi, DShape &dpsidxi, DShape &d2psidxi) const
 
virtual double d2shape_lagrangian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidxi, DShape &d2psidxi) const
 
unsigned lagrangian_dimension () const
 
unsigned nnodal_lagrangian_type () const
 
Nodeconstruct_node (const unsigned &n)
 Construct the local node n and return a pointer to it. More...
 
Nodeconstruct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
Nodeconstruct_boundary_node (const unsigned &n)
 
Nodeconstruct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual void assign_all_generic_local_eqn_numbers (const bool &store_local_dof_pt)
 
void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
double raw_lagrangian_position (const unsigned &n, const unsigned &i) const
 
double raw_lagrangian_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double lagrangian_position (const unsigned &n, const unsigned &i) const
 Return i-th Lagrangian coordinate at local node n. More...
 
double lagrangian_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual double interpolated_xi (const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_xi (const Vector< double > &s, Vector< double > &xi) const
 
virtual void interpolated_dxids (const Vector< double > &s, DenseMatrix< double > &dxids) const
 
virtual void J_lagrangian (const Vector< double > &s) const
 
virtual double J_lagrangian_at_knot (const unsigned &ipt) const
 
SolidInitialCondition *& solid_ic_pt ()
 Pointer to object that describes the initial condition. More...
 
void enable_solve_for_consistent_newmark_accel ()
 
void disable_solve_for_consistent_newmark_accel ()
 Set to reset the problem being solved to be the standard problem. More...
 
MultiplierFctPtmultiplier_fct_pt ()
 
MultiplierFctPt multiplier_fct_pt () const
 
virtual void get_residuals_for_solid_ic (Vector< double > &residuals)
 
void fill_in_residuals_for_solid_ic (Vector< double > &residuals)
 
void fill_in_jacobian_for_solid_ic (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void fill_in_jacobian_for_newmark_accel (DenseMatrix< double > &jacobian)
 
void compute_norm (double &el_norm)
 
int position_local_eqn (const unsigned &n, const unsigned &k, const unsigned &j) const
 
- Public Member Functions inherited from oomph::FiniteElement
void set_dimension (const unsigned &dim)
 
void set_nodal_dimension (const unsigned &nodal_dim)
 
void set_nnodal_position_type (const unsigned &nposition_type)
 Set the number of types required to interpolate the coordinate. More...
 
void set_n_node (const unsigned &n)
 
int nodal_local_eqn (const unsigned &n, const unsigned &i) const
 
double dJ_eulerian_at_knot (const unsigned &ipt, Shape &psi, DenseMatrix< double > &djacobian_dX) const
 
 FiniteElement ()
 Constructor. More...
 
virtual ~FiniteElement ()
 
 FiniteElement (const FiniteElement &)=delete
 Broken copy constructor. More...
 
virtual bool local_coord_is_valid (const Vector< double > &s)
 Broken assignment operator. More...
 
virtual void move_local_coord_back_into_element (Vector< double > &s) const
 
void get_centre_of_gravity_and_max_radius_in_terms_of_zeta (Vector< double > &cog, double &max_radius) const
 
virtual void local_coordinate_of_node (const unsigned &j, Vector< double > &s) const
 
virtual void local_fraction_of_node (const unsigned &j, Vector< double > &s_fraction)
 
virtual double local_one_d_fraction_of_node (const unsigned &n1d, const unsigned &i)
 
MacroElementmacro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void assign_nodal_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void describe_nodal_local_dofs (std::ostream &out, const std::string &current_string) const
 
Node *& node_pt (const unsigned &n)
 Return a pointer to the local node n. More...
 
Node *const & node_pt (const unsigned &n) const
 Return a pointer to the local node n (const version) More...
 
unsigned nnode () const
 Return the number of nodes. More...
 
virtual unsigned nnode_1d () const
 
double raw_nodal_position (const unsigned &n, const unsigned &i) const
 
double raw_nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position (const unsigned &n, const unsigned &i) const
 
double nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double dnodal_position_dt (const unsigned &n, const unsigned &i) const
 Return the i-th component of nodal velocity: dx/dt at local node n. More...
 
double dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
virtual void disable_ALE ()
 
virtual void enable_ALE ()
 
virtual unsigned required_nvalue (const unsigned &n) const
 
unsigned nnodal_position_type () const
 
bool has_hanging_nodes () const
 
unsigned nodal_dimension () const
 Return the required Eulerian dimension of the nodes in this element. More...
 
virtual unsigned nvertex_node () const
 
virtual Nodevertex_node_pt (const unsigned &j) const
 
int get_node_number (Node *const &node_pt) const
 
virtual Nodeget_node_at_local_coordinate (const Vector< double > &s) const
 
double raw_nodal_value (const unsigned &n, const unsigned &i) const
 
double raw_nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
unsigned dim () const
 
virtual ElementGeometry::ElementGeometry element_geometry () const
 Return the geometry type of the element (either Q or T usually). More...
 
virtual double interpolated_x (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated coordinate x[i] at local coordinate s. More...
 
virtual double interpolated_x (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_x (const Vector< double > &s, Vector< double > &x) const
 Return FE interpolated position x[] at local coordinate s as Vector. More...
 
virtual void interpolated_x (const unsigned &t, const Vector< double > &s, Vector< double > &x) const
 
virtual double interpolated_dxdt (const Vector< double > &s, const unsigned &i, const unsigned &t)
 
virtual void interpolated_dxdt (const Vector< double > &s, const unsigned &t, Vector< double > &dxdt)
 
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)
 
void interpolated_zeta (const Vector< double > &s, Vector< double > &zeta) const
 
void locate_zeta (const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
 
virtual void node_update ()
 
virtual void identify_field_data_for_interactions (std::set< std::pair< Data *, unsigned >> &paired_field_data)
 
virtual double s_min () const
 Min value of local coordinate. More...
 
virtual double s_max () const
 Max. value of local coordinate. More...
 
double size () const
 
virtual double compute_physical_size () const
 
virtual void point_output_data (const Vector< double > &s, Vector< double > &data)
 
void point_output (std::ostream &outfile, const Vector< double > &s)
 
virtual unsigned nplot_points_paraview (const unsigned &nplot) const
 
virtual unsigned nsub_elements_paraview (const unsigned &nplot) const
 
void output_paraview (std::ofstream &file_out, const unsigned &nplot) const
 
virtual void write_paraview_output_offset_information (std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
 
virtual void write_paraview_type (std::ofstream &file_out, const unsigned &nplot) const
 
virtual void write_paraview_offsets (std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
 
virtual unsigned nscalar_paraview () const
 
virtual void scalar_value_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
 
virtual void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
 
virtual void scalar_value_fct_paraview (std::ofstream &file_out, const unsigned &i, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
virtual std::string scalar_name_paraview (const unsigned &i) const
 
virtual void output (const unsigned &t, std::ostream &outfile, const unsigned &n_plot) const
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 Output an exact solution over the element. More...
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 Output a time-dependent exact solution over the element. More...
 
virtual void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, const SolutionFunctorBase &exact_soln) const
 Output a time-dependent exact solution over the element. More...
 
virtual void get_s_plot (const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
 
virtual std::string tecplot_zone_string (const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (std::ostream &outfile, const unsigned &nplot) const
 
virtual void write_tecplot_zone_footer (FILE *file_pt, const unsigned &nplot) const
 
virtual unsigned nplot_points (const unsigned &nplot) const
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Calculate the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_abs_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error)
 
void integrate_fct (FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt, Vector< double > &integral)
 Integrate Vector-valued function over element. More...
 
void integrate_fct (FiniteElement::UnsteadyExactSolutionFctPt integrand_fct_pt, const double &time, Vector< double > &integral)
 Integrate Vector-valued time-dep function over element. More...
 
virtual void build_face_element (const int &face_index, FaceElement *face_element_pt)
 
virtual unsigned self_test ()
 
virtual unsigned get_bulk_node_number (const int &face_index, const unsigned &i) const
 
virtual int face_outer_unit_normal_sign (const int &face_index) const
 Get the sign of the outer unit normal on the face given by face_index. More...
 
virtual unsigned nnode_on_face () const
 
void face_node_number_error_check (const unsigned &i) const
 Range check for face node numbers. More...
 
virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt (const int &face_index) const
 
virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt (const int &face_index) const
 
- Public Member Functions inherited from oomph::GeneralisedElement
 GeneralisedElement ()
 Constructor: Initialise all pointers and all values to zero. More...
 
virtual ~GeneralisedElement ()
 Virtual destructor to clean up any memory allocated by the object. More...
 
 GeneralisedElement (const GeneralisedElement &)=delete
 Broken copy constructor. More...
 
void operator= (const GeneralisedElement &)=delete
 Broken assignment operator. More...
 
Data *& internal_data_pt (const unsigned &i)
 Return a pointer to i-th internal data object. More...
 
Data *const & internal_data_pt (const unsigned &i) const
 Return a pointer to i-th internal data object (const version) More...
 
Data *& external_data_pt (const unsigned &i)
 Return a pointer to i-th external data object. More...
 
Data *const & external_data_pt (const unsigned &i) const
 Return a pointer to i-th external data object (const version) More...
 
unsigned long eqn_number (const unsigned &ieqn_local) const
 
int local_eqn_number (const unsigned long &ieqn_global) const
 
unsigned add_external_data (Data *const &data_pt, const bool &fd=true)
 
bool external_data_fd (const unsigned &i) const
 
void exclude_external_data_fd (const unsigned &i)
 
void include_external_data_fd (const unsigned &i)
 
void flush_external_data ()
 Flush all external data. More...
 
void flush_external_data (Data *const &data_pt)
 Flush the object addressed by data_pt from the external data array. More...
 
unsigned ninternal_data () const
 Return the number of internal data objects. More...
 
unsigned nexternal_data () const
 Return the number of external data objects. More...
 
unsigned ndof () const
 Return the number of equations/dofs in the element. More...
 
void dof_vector (const unsigned &t, Vector< double > &dof)
 Return the vector of dof values at time level t. More...
 
void dof_pt_vector (Vector< double * > &dof_pt)
 Return the vector of pointers to dof values. More...
 
void set_internal_data_time_stepper (const unsigned &i, TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void assign_internal_eqn_numbers (unsigned long &global_number, Vector< double * > &Dof_pt)
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void add_internal_value_pt_to_map (std::map< unsigned, double * > &map_of_value_pt)
 
virtual void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void complete_setup_of_dependencies ()
 
virtual void get_residuals (Vector< double > &residuals)
 
virtual void get_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
virtual void get_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void get_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void get_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void get_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void get_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void get_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void get_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void get_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
virtual void compute_norm (Vector< double > &norm)
 
- Public Member Functions inherited from oomph::GeomObject
 GeomObject ()
 Default constructor. More...
 
 GeomObject (const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)
 
 GeomObject (const GeomObject &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const GeomObject &)=delete
 Broken assignment operator. More...
 
virtual ~GeomObject ()
 (Empty) destructor More...
 
unsigned nlagrangian () const
 Access function to # of Lagrangian coordinates. More...
 
unsigned ndim () const
 Access function to # of Eulerian coordinates. More...
 
void set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
 Set # of Lagrangian and Eulerian coordinates. More...
 
TimeStepper *& time_stepper_pt ()
 
TimeSteppertime_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 
- Public Member Functions inherited from oomph::SolidElementWithDiagonalMassMatrix
 SolidElementWithDiagonalMassMatrix ()
 Empty constructor. More...
 
virtual ~SolidElementWithDiagonalMassMatrix ()
 Virtual destructor. More...
 
 SolidElementWithDiagonalMassMatrix (const SolidElementWithDiagonalMassMatrix &)=delete
 Broken copy constructor. More...
 
void operator= (const SolidElementWithDiagonalMassMatrix &)=delete
 Broken assignment operator. More...
 

Protected Member Functions

void get_stress (const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma_dev, DenseMatrix< double > &Gcontra, double &gen_dil, double &inv_kappa)
 
void get_d_stress_dG_upper (const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, const double &gen_dil, const double &inv_kappa, const double &interpolated_solid_p, RankFourTensor< double > &d_sigma_dG, DenseMatrix< double > &d_gen_dil_dG)
 
virtual void solid_pshape (const Vector< double > &s, Shape &psi) const =0
 Return the solid pressure shape functions. More...
 
void solid_pshape_at_knot (const unsigned &ipt, Shape &psi) const
 Return the stored solid shape functions at the knots. More...
 
virtual void fill_in_generic_residual_contribution_pvd_with_pressure (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &flag)
 
void get_stress (const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma_dev, DenseMatrix< double > &Gcontra, double &detG)
 
void get_d_stress_dG_upper (const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, const double &detG, const double &interpolated_solid_p, RankFourTensor< double > &d_sigma_dG, DenseMatrix< double > &d_detG_dG)
 
- Protected Member Functions inherited from oomph::SolidFiniteElement
void fill_in_generic_jacobian_for_solid_ic (Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
 
void set_nnodal_lagrangian_type (const unsigned &nlagrangian_type)
 
virtual double local_to_lagrangian_mapping (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
double local_to_lagrangian_mapping (const DShape &dpsids, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_lagrangian_mapping_diagonal (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual void assign_solid_local_eqn_numbers (const bool &store_local_dof)
 Assign local equation numbers for the solid equations in the element. More...
 
void describe_solid_local_dofs (std::ostream &out, const std::string &current_string) const
 Classifies dofs locally for solid specific aspects. More...
 
virtual void fill_in_jacobian_from_solid_position_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void fill_in_jacobian_from_solid_position_by_fd (DenseMatrix< double > &jacobian)
 
virtual void update_before_solid_position_fd ()
 
virtual void reset_after_solid_position_fd ()
 
virtual void update_in_solid_position_fd (const unsigned &i)
 
virtual void reset_in_solid_position_fd (const unsigned &i)
 
- Protected Member Functions inherited from oomph::FiniteElement
virtual void assemble_local_to_eulerian_jacobian (const DShape &dpsids, DenseMatrix< double > &jacobian) const
 
virtual void assemble_local_to_eulerian_jacobian2 (const DShape &d2psids, DenseMatrix< double > &jacobian2) const
 
virtual void assemble_eulerian_base_vectors (const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
 
template<unsigned DIM>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double invert_jacobian_mapping (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
double local_to_eulerian_mapping (const DShape &dpsids, DenseMatrix< double > &inverse_jacobian) const
 
virtual double local_to_eulerian_mapping_diagonal (const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
virtual void dJ_eulerian_dnodal_coordinates (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<unsigned DIM>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
virtual void d_dshape_eulerian_dnodal_coordinates (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<unsigned DIM>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
virtual void transform_derivatives (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
void transform_derivatives_diagonal (const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
 
virtual void transform_second_derivatives (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<unsigned DIM>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
virtual void fill_in_jacobian_from_nodal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
void fill_in_jacobian_from_nodal_by_fd (DenseMatrix< double > &jacobian)
 
virtual void update_before_nodal_fd ()
 
virtual void reset_after_nodal_fd ()
 
virtual void update_in_nodal_fd (const unsigned &i)
 
virtual void reset_in_nodal_fd (const unsigned &i)
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Zero-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 One-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 Two-d specialisation of function to calculate inverse of jacobian mapping. More...
 
template<>
double invert_jacobian (const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void dJ_eulerian_dnodal_coordinates_templated_helper (const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void d_dshape_eulerian_dnodal_coordinates_templated_helper (const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_template (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
template<>
void transform_second_derivatives_diagonal (const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
 
- Protected Member Functions inherited from oomph::GeneralisedElement
unsigned add_internal_data (Data *const &data_pt, const bool &fd=true)
 
bool internal_data_fd (const unsigned &i) const
 
void exclude_internal_data_fd (const unsigned &i)
 
void include_internal_data_fd (const unsigned &i)
 
void clear_global_eqn_numbers ()
 
void add_global_eqn_numbers (std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)
 
virtual void assign_internal_and_external_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void assign_additional_local_eqn_numbers ()
 
int internal_local_eqn (const unsigned &i, const unsigned &j) const
 
int external_local_eqn (const unsigned &i, const unsigned &j)
 
void fill_in_jacobian_from_internal_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_internal_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
void fill_in_jacobian_from_external_by_fd (DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
 
virtual void update_before_internal_fd ()
 
virtual void reset_after_internal_fd ()
 
virtual void update_in_internal_fd (const unsigned &i)
 
virtual void reset_in_internal_fd (const unsigned &i)
 
virtual void update_before_external_fd ()
 
virtual void reset_after_external_fd ()
 
virtual void update_in_external_fd (const unsigned &i)
 
virtual void reset_in_external_fd (const unsigned &i)
 
virtual void fill_in_contribution_to_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void fill_in_contribution_to_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void fill_in_contribution_to_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void fill_in_contribution_to_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void fill_in_contribution_to_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void fill_in_contribution_to_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 

Protected Attributes

bool Incompressible
 Boolean to determine whether the solid is incompressible or not. More...
 
- Protected Attributes inherited from oomph::PVDEquationsBase< DIM >
IsotropicGrowthFctPt Isotropic_growth_fct_pt
 Pointer to isotropic growth function. More...
 
PrestressFctPt Prestress_fct_pt
 Pointer to prestress function. More...
 
ConstitutiveLawConstitutive_law_pt
 Pointer to the constitutive law. More...
 
doubleLambda_sq_pt
 Timescale ratio (non-dim. density) More...
 
bool Unsteady
 Flag that switches inertia on/off. More...
 
BodyForceFctPt Body_force_fct_pt
 Pointer to body force function. More...
 
bool Evaluate_jacobian_by_fd
 Use FD to evaluate Jacobian. More...
 
- Protected Attributes inherited from oomph::SolidFiniteElement
MacroElementUndeformed_macro_elem_pt
 Pointer to the element's "undeformed" macro element (NULL by default) More...
 
SolidInitialConditionSolid_ic_pt
 Pointer to object that specifies the initial condition. More...
 
bool Solve_for_consistent_newmark_accel_flag
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElementMacro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepperGeom_object_time_stepper_pt
 

Additional Inherited Members

- Public Types inherited from oomph::PVDEquationsBase< DIM >
typedef void(* IsotropicGrowthFctPt) (const Vector< double > &xi, double &gamma)
 
typedef double(* PrestressFctPt) (const unsigned &i, const unsigned &j, const Vector< double > &xi)
 
typedef void(* BodyForceFctPt) (const double &t, const Vector< double > &xi, Vector< double > &b)
 
- Public Types inherited from oomph::SolidFiniteElement
typedef double(* MultiplierFctPt) (const Vector< double > &xi)
 
- 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 Member Functions inherited from oomph::PVDEquationsBase< DIM >
static void pin_redundant_nodal_solid_pressures (const Vector< GeneralisedElement * > &element_pt)
 
static void unpin_all_solid_pressure_dofs (const Vector< GeneralisedElement * > &element_pt)
 Unpin all pressure dofs in elements listed in vector. 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 Attributes inherited from oomph::PVDEquationsBase< DIM >
static double Default_lambda_sq_value = 1.0
 Static default value for timescale ratio (1.0 – for natural scaling) 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< doubleDummy_matrix
 
static std::deque< double * > Dof_pt_deque
 

Detailed Description

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

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// A class for elements that solve the equations of solid mechanics, based on the principle of virtual displacements, with a contitutive equation that involves a pressure. This formulation is required in the case of incompressible materials, in which the additional constraint that volume must be conserved is applied. In this case, the Incompressible flag must be set to true. If the Incompressible flag is not set to true, we use the nearly-incompressible formulation of the constitutive equations.

Constructor & Destructor Documentation

◆ PVDEquationsWithPressure()

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

Constructor, by default the element is NOT incompressible.

866  : PVDEquationsBase<DIM>(), Incompressible(false)
867  {
868  }
bool Incompressible
Boolean to determine whether the solid is incompressible or not.
Definition: solid_elements.h:1237

Member Function Documentation

◆ extended_output()

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::extended_output ( std::ostream &  outfile,
const unsigned n_plot 
)

Output: x,y,[z],xi0,xi1,[xi2],gamma strain and stress components.

Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress components

2118  {
2119  Vector<double> x(DIM);
2120  Vector<double> xi(DIM);
2121  Vector<double> s(DIM);
2122  DenseMatrix<double> stress_or_strain(DIM, DIM);
2123 
2124  // Tecplot header info
2125  outfile << this->tecplot_zone_string(n_plot);
2126 
2127  // Loop over plot points
2128  unsigned num_plot_points = this->nplot_points(n_plot);
2129  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
2130  {
2131  // Get local coordinates of plot point
2132  this->get_s_plot(iplot, n_plot, s);
2133 
2134  // Get Eulerian and Lagrangian coordinates
2135  this->interpolated_x(s, x);
2136  this->interpolated_xi(s, xi);
2137 
2138  // Get isotropic growth
2139  double gamma;
2140  // Dummy integration point
2141  unsigned ipt = 0;
2142  this->get_isotropic_growth(ipt, s, xi, gamma);
2143 
2144  // Output the x,y,..
2145  for (unsigned i = 0; i < DIM; i++)
2146  {
2147  outfile << x[i] << " ";
2148  }
2149 
2150  // Output xi0,xi1,..
2151  for (unsigned i = 0; i < DIM; i++)
2152  {
2153  outfile << xi[i] << " ";
2154  }
2155 
2156  // Output growth
2157  outfile << gamma << " ";
2158 
2159  // Output pressure
2160  outfile << interpolated_solid_p(s) << " ";
2161 
2162  // get the strain
2163  this->get_strain(s, stress_or_strain);
2164  for (unsigned i = 0; i < DIM; i++)
2165  {
2166  for (unsigned j = 0; j <= i; j++)
2167  {
2168  outfile << stress_or_strain(j, i) << " ";
2169  }
2170  }
2171 
2172  // get the stress
2173  this->get_stress(s, stress_or_strain);
2174  for (unsigned i = 0; i < DIM; i++)
2175  {
2176  for (unsigned j = 0; j <= i; j++)
2177  {
2178  outfile << stress_or_strain(j, i) << " ";
2179  }
2180  }
2181 
2182 
2183  outfile << std::endl;
2184  }
2185 
2186 
2187  // Write tecplot footer (e.g. FE connectivity lists)
2188  this->write_tecplot_zone_footer(outfile, n_plot);
2189  outfile << std::endl;
2190  }
int i
Definition: BiCGSTAB_step_by_step.cpp:9
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Definition: elements.h:3161
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
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Definition: elements.h:3186
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Definition: elements.h:3174
virtual void get_isotropic_growth(const unsigned &ipt, const Vector< double > &s, const Vector< double > &xi, double &gamma) const
Definition: solid_elements.h:267
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
Definition: solid_elements.cc:47
void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma)
Definition: solid_elements.cc:2267
double interpolated_solid_p(const Vector< double > &s)
Return the interpolated_solid_pressure.
Definition: solid_elements.h:1011
virtual double interpolated_xi(const Vector< double > &s, const unsigned &i) const
Definition: elements.cc:7104
RealScalar s
Definition: level1_cplx_impl.h:130
#define DIM
Definition: linearised_navier_stokes_elements.h:44
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:116
list x
Definition: plotDoE.py:28
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References DIM, mathsFunc::gamma(), i, j, s, and plotDoE::x.

◆ fill_in_contribution_to_jacobian()

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

Fill in contribution to Jacobian (either by FD or analytically, for the positional variables; control this via evaluate_jacobian_by_fd(). Note: Jacobian entries arising from derivatives w.r.t. pressure terms are always computed analytically.

Reimplemented from oomph::SolidFiniteElement.

916  {
917  // Solve for the consistent acceleration in the Newmark scheme
918  // Note that this replaces solid entries only
920  (this->Solid_ic_pt != 0))
921  {
922  std::string error_message = "Can't assign consistent Newmark history\n";
923  error_message += " values for solid element with pressure dofs\n";
924 
925  throw OomphLibError(
927  }
928 
929  // FD
930  if (this->Evaluate_jacobian_by_fd)
931  {
932  // Call the generic routine with the flag set to 2: Computes residuals
933  // and derivatives w.r.t. to pressure variables
935  residuals, jacobian, GeneralisedElement::Dummy_matrix, 2);
936 
937  // Call the finite difference routine for the deriatives w.r.t.
938  // the positional variables
940  }
941  // Do it fully analytically
942  else
943  {
944  // Call the generic routine with the flag set to 1: Get residual
945  // and fully analytical Jacobian
947  residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
948  }
949  }
static DenseMatrix< double > Dummy_matrix
Definition: elements.h:227
bool Evaluate_jacobian_by_fd
Use FD to evaluate Jacobian.
Definition: solid_elements.h:430
virtual void fill_in_generic_residual_contribution_pvd_with_pressure(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &flag)
Definition: solid_elements.cc:1271
bool Solve_for_consistent_newmark_accel_flag
Definition: elements.h:4302
SolidInitialCondition * Solid_ic_pt
Pointer to object that specifies the initial condition.
Definition: elements.h:4131
virtual void fill_in_jacobian_from_solid_position_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Definition: elements.cc:6985
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References oomph::GeneralisedElement::Dummy_matrix, oomph::PVDEquationsBase< DIM >::Evaluate_jacobian_by_fd, oomph::PVDEquationsWithPressure< DIM >::fill_in_generic_residual_contribution_pvd_with_pressure(), oomph::SolidFiniteElement::fill_in_jacobian_from_solid_position_by_fd(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::SolidFiniteElement::Solid_ic_pt, oomph::SolidFiniteElement::Solve_for_consistent_newmark_accel_flag, and oomph::Global_string_for_annotation::string().

◆ fill_in_contribution_to_jacobian_and_mass_matrix()

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

Fill in contribution to Mass matrix and Jacobian (either by FD or analytically, for the positional variables; control this via evaluate_jacobian_by_fd(). Note: Jacobian entries arising from derivatives w.r.t. pressure terms are always computed analytically. Note that the Jacobian is multiplied by minus one to ensure that the mass matrix is positive semi-definite.

Reimplemented from oomph::GeneralisedElement.

962  {
963  // Solve for the consistent acceleration in the Newmark scheme
964  // Note that this replaces solid entries only
966  (this->Solid_ic_pt != 0))
967  {
968  std::string error_message = "Can't assign consistent Newmark history\n";
969  error_message += " values for solid element with pressure dofs\n";
970 
971  throw OomphLibError(
973  }
974 
975  // FD
976  if (this->Evaluate_jacobian_by_fd)
977  {
978  // Call the generic routine with the flag set to 4: Computes residuals
979  // and derivatives w.r.t. to pressure variables
980  // and the mass matrix
982  residuals, jacobian, mass_matrix, 4);
983 
984  // Call the finite difference routine for the deriatives w.r.t.
985  // the positional variables
987  }
988  // Do it fully analytically
989  else
990  {
991  // Call the generic routine with the flag set to 3: Get residual
992  // and fully analytical Jacobian
994  residuals, jacobian, mass_matrix, 3);
995  }
996 
997  // Multiply the residuals and jacobian by minus one
998  const unsigned n_dof = this->ndof();
999  for (unsigned i = 0; i < n_dof; i++)
1000  {
1001  residuals[i] *= -1.0;
1002  for (unsigned j = 0; j < n_dof; j++)
1003  {
1004  jacobian(i, j) *= -1.0;
1005  }
1006  }
1007  }
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835

References oomph::PVDEquationsBase< DIM >::Evaluate_jacobian_by_fd, oomph::PVDEquationsWithPressure< DIM >::fill_in_generic_residual_contribution_pvd_with_pressure(), oomph::SolidFiniteElement::fill_in_jacobian_from_solid_position_by_fd(), i, j, oomph::GeneralisedElement::ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::SolidFiniteElement::Solid_ic_pt, oomph::SolidFiniteElement::Solve_for_consistent_newmark_accel_flag, and oomph::Global_string_for_annotation::string().

◆ fill_in_contribution_to_residuals()

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

Fill in the residuals.

Reimplemented from oomph::GeneralisedElement.

900  {
901  // Call the generic residuals function with flag set to 0
902  // using a dummy matrix argument
904  residuals,
907  0);
908  }

References oomph::GeneralisedElement::Dummy_matrix, and oomph::PVDEquationsWithPressure< DIM >::fill_in_generic_residual_contribution_pvd_with_pressure().

◆ fill_in_generic_residual_contribution_pvd_with_pressure()

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

Returns the residuals for the discretised principle of virtual displacements, formulated in the incompressible/ near-incompressible case.

  • If flag==0, compute only the residual vector.
  • If flag==1, compute residual vector and fully analytical Jacobian
  • If flag==2, also compute the pressure-related entries in the Jacobian (all others need to be done by finite differencing.

Fill in element's contribution to the elemental residual vector and/or Jacobian matrix. flag=0: compute only residual vector flag=1: compute both, fully analytically flag=2: compute both, using FD for the derivatives w.r.t. to the discrete displacment dofs. flag=3: compute residuals, jacobian (full analytic) and mass matrix flag=4: compute residuals, jacobian (FD for derivatives w.r.t. displacements) and mass matrix

Reimplemented in oomph::RefineablePVDEquationsWithPressure< DIM >.

1276  {
1277 #ifdef PARANOID
1278  // Check if the constitutive equation requires the explicit imposition of an
1279  // incompressibility constraint
1281  (!Incompressible))
1282  {
1283  throw OomphLibError("The constitutive law requires the use of the "
1284  "incompressible formulation by setting the element's "
1285  "member function set_incompressible()",
1288  }
1289 #endif
1290 
1291 
1292  // Simply set up initial condition?
1293  if (this->Solid_ic_pt != 0)
1294  {
1295  this->get_residuals_for_solid_ic(residuals);
1296  return;
1297  }
1298 
1299  // Find out how many nodes there are
1300  const unsigned n_node = this->nnode();
1301 
1302  // Find out how many position types of dof there are
1303  const unsigned n_position_type = this->nnodal_position_type();
1304 
1305  // Find out how many pressure dofs there are
1306  const unsigned n_solid_pres = this->npres_solid();
1307 
1308  // Set up memory for the shape functions
1309  Shape psi(n_node, n_position_type);
1310  DShape dpsidxi(n_node, n_position_type, DIM);
1311 
1312  // Set up memory for the pressure shape functions
1313  Shape psisp(n_solid_pres);
1314 
1315  // Set the value of n_intpt
1316  const unsigned n_intpt = this->integral_pt()->nweight();
1317 
1318  // Set the vector to hold the local coordinates in the element
1319  Vector<double> s(DIM);
1320 
1321  // Timescale ratio (non-dim density)
1322  double lambda_sq = this->lambda_sq();
1323 
1324  // Time factor
1325  double time_factor = 0.0;
1326  if (lambda_sq > 0)
1327  {
1328  time_factor = this->node_pt(0)->position_time_stepper_pt()->weight(2, 0);
1329  }
1330 
1331  // Integers to hold the local equation and unknown numbers
1332  int local_eqn = 0, local_unknown = 0;
1333 
1334  // Loop over the integration points
1335  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1336  {
1337  // Assign the values of s
1338  for (unsigned i = 0; i < DIM; ++i)
1339  {
1340  s[i] = this->integral_pt()->knot(ipt, i);
1341  }
1342 
1343  // Get the integral weight
1344  double w = this->integral_pt()->weight(ipt);
1345 
1346  // Call the derivatives of the shape functions
1347  double J = this->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
1348 
1349  // Call the pressure shape functions
1350  solid_pshape_at_knot(ipt, psisp);
1351 
1352  // Storage for Lagrangian coordinates (initialised to zero)
1353  Vector<double> interpolated_xi(DIM, 0.0);
1354 
1355  // Deformed tangent vectors
1356  DenseMatrix<double> interpolated_G(DIM);
1357 
1358  // Setup memory for accelerations
1359  Vector<double> accel(DIM);
1360 
1361  // Initialise to zero
1362  for (unsigned i = 0; i < DIM; i++)
1363  {
1364  // Initialise acclerations
1365  accel[i] = 0.0;
1366  for (unsigned j = 0; j < DIM; j++)
1367  {
1368  interpolated_G(i, j) = 0.0;
1369  }
1370  }
1371 
1372  // Calculate displacements and derivatives and lagrangian coordinates
1373  for (unsigned l = 0; l < n_node; l++)
1374  {
1375  // Loop over positional dofs
1376  for (unsigned k = 0; k < n_position_type; k++)
1377  {
1378  // Loop over displacement components (deformed position)
1379  for (unsigned i = 0; i < DIM; i++)
1380  {
1381  // Calculate the lagrangian coordinates and the accelerations
1382  interpolated_xi[i] +=
1383  this->lagrangian_position_gen(l, k, i) * psi(l, k);
1384 
1385  // Only compute accelerations if inertia is switched on
1386  // otherwise the timestepper might not be able to
1387  // work out dx_gen_dt(2,...)
1388  if ((lambda_sq > 0.0) && (this->Unsteady))
1389  {
1390  accel[i] += this->dnodal_position_gen_dt(2, l, k, i) * psi(l, k);
1391  }
1392 
1393  // Loop over derivative directions
1394  for (unsigned j = 0; j < DIM; j++)
1395  {
1396  interpolated_G(j, i) +=
1397  this->nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
1398  }
1399  }
1400  }
1401  }
1402 
1403  // Get isotropic growth factor
1404  double gamma = 1.0;
1406 
1407  // Get body force at current time
1408  Vector<double> b(DIM);
1409  this->body_force(interpolated_xi, b);
1410 
1411  // We use Cartesian coordinates as the reference coordinate
1412  // system. In this case the undeformed metric tensor is always
1413  // the identity matrix -- stretched by the isotropic growth
1414  double diag_entry = pow(gamma, 2.0 / double(DIM));
1416  for (unsigned i = 0; i < DIM; i++)
1417  {
1418  for (unsigned j = 0; j < DIM; j++)
1419  {
1420  if (i == j)
1421  {
1422  g(i, j) = diag_entry;
1423  }
1424  else
1425  {
1426  g(i, j) = 0.0;
1427  }
1428  }
1429  }
1430 
1431  // Premultiply the undeformed volume ratio (from the isotropic
1432  // growth), the weights and the Jacobian
1433  double W = gamma * w * J;
1434 
1435  // Calculate the interpolated solid pressure
1436  double interpolated_solid_p = 0.0;
1437  for (unsigned l = 0; l < n_solid_pres; l++)
1438  {
1439  interpolated_solid_p += solid_p(l) * psisp[l];
1440  }
1441 
1442 
1443  // Declare and calculate the deformed metric tensor
1445 
1446  // Assign values of G
1447  for (unsigned i = 0; i < DIM; i++)
1448  {
1449  // Do upper half of matrix
1450  for (unsigned j = i; j < DIM; j++)
1451  {
1452  // Initialise G(i,j) to zero
1453  G(i, j) = 0.0;
1454  // Now calculate the dot product
1455  for (unsigned k = 0; k < DIM; k++)
1456  {
1457  G(i, j) += interpolated_G(i, k) * interpolated_G(j, k);
1458  }
1459  }
1460  // Matrix is symmetric so just copy lower half
1461  for (unsigned j = 0; j < i; j++)
1462  {
1463  G(i, j) = G(j, i);
1464  }
1465  }
1466 
1467  // Now calculate the deviatoric stress and all pressure-related
1468  // quantitites
1469  DenseMatrix<double> sigma(DIM, DIM), sigma_dev(DIM, DIM), Gup(DIM, DIM);
1470  double detG = 0.0;
1471  double gen_dil = 0.0;
1472  double inv_kappa = 0.0;
1473 
1474  // Get stress derivative by FD only needed for Jacobian
1475 
1476  // Stress etc derivatives
1477  RankFourTensor<double> d_stress_dG(DIM, DIM, DIM, DIM, 0.0);
1478  DenseMatrix<double> d_detG_dG(DIM, DIM, 0.0);
1479  DenseMatrix<double> d_gen_dil_dG(DIM, DIM, 0.0);
1480 
1481  // Derivative of metric tensor w.r.t. to nodal coords
1482  RankFiveTensor<double> d_G_dX(
1483  n_node, n_position_type, DIM, DIM, DIM, 0.0);
1484 
1485  // Get Jacobian too?
1486  if ((flag == 1) || (flag == 3))
1487  {
1488  // Derivative of metric tensor w.r.t. to discrete positional dofs
1489  // NOTE: Since G is symmetric we only compute the upper triangle
1490  // and DO NOT copy the entries across. Subsequent computations
1491  // must (and, in fact, do) therefore only operate with upper
1492  // triangular entries
1493  for (unsigned ll = 0; ll < n_node; ll++)
1494  {
1495  for (unsigned kk = 0; kk < n_position_type; kk++)
1496  {
1497  for (unsigned ii = 0; ii < DIM; ii++)
1498  {
1499  for (unsigned aa = 0; aa < DIM; aa++)
1500  {
1501  for (unsigned bb = aa; bb < DIM; bb++)
1502  {
1503  d_G_dX(ll, kk, ii, aa, bb) =
1504  interpolated_G(aa, ii) * dpsidxi(ll, kk, bb) +
1505  interpolated_G(bb, ii) * dpsidxi(ll, kk, aa);
1506  }
1507  }
1508  }
1509  }
1510  }
1511  }
1512 
1513 
1514  // Incompressible: Compute the deviatoric part of the stress tensor, the
1515  // contravariant deformed metric tensor and the determinant
1516  // of the deformed covariant metric tensor.
1517  if (Incompressible)
1518  {
1519  get_stress(g, G, sigma_dev, Gup, detG);
1520 
1521  // Get full stress
1522  for (unsigned a = 0; a < DIM; a++)
1523  {
1524  for (unsigned b = 0; b < DIM; b++)
1525  {
1526  sigma(a, b) = sigma_dev(a, b) - interpolated_solid_p * Gup(a, b);
1527  }
1528  }
1529 
1530  // Get Jacobian too?
1531  if ((flag == 1) || (flag == 3))
1532  {
1533  // Get the "upper triangular" entries of the derivatives of the stress
1534  // tensor with respect to G
1535  this->get_d_stress_dG_upper(
1536  g, G, sigma, detG, interpolated_solid_p, d_stress_dG, d_detG_dG);
1537  }
1538  }
1539  // Nearly incompressible: Compute the deviatoric part of the
1540  // stress tensor, the contravariant deformed metric tensor,
1541  // the generalised dilatation and the inverse bulk modulus.
1542  else
1543  {
1544  get_stress(g, G, sigma_dev, Gup, gen_dil, inv_kappa);
1545 
1546  // Get full stress
1547  for (unsigned a = 0; a < DIM; a++)
1548  {
1549  for (unsigned b = 0; b < DIM; b++)
1550  {
1551  sigma(a, b) = sigma_dev(a, b) - interpolated_solid_p * Gup(a, b);
1552  }
1553  }
1554 
1555  // Get Jacobian too?
1556  if ((flag == 1) || (flag == 3))
1557  {
1558  // Get the "upper triangular" entries of the derivatives of the stress
1559  // tensor with respect to G
1560  this->get_d_stress_dG_upper(g,
1561  G,
1562  sigma,
1563  gen_dil,
1564  inv_kappa,
1566  d_stress_dG,
1567  d_gen_dil_dG);
1568  }
1569  }
1570 
1571  // Add pre-stress
1572  for (unsigned i = 0; i < DIM; i++)
1573  {
1574  for (unsigned j = 0; j < DIM; j++)
1575  {
1576  sigma(i, j) += this->prestress(i, j, interpolated_xi);
1577  }
1578  }
1579 
1580  //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL
1581  // DISPLACEMENTS========
1582 
1583  // Loop over the test functions, nodes of the element
1584  for (unsigned l = 0; l < n_node; l++)
1585  {
1586  // Loop over the types of dof
1587  for (unsigned k = 0; k < n_position_type; k++)
1588  {
1589  // Offset for faster access
1590  const unsigned offset5 = dpsidxi.offset(l, k);
1591 
1592  // Loop over the displacement components
1593  for (unsigned i = 0; i < DIM; i++)
1594  {
1595  // Get the equation number
1596  local_eqn = this->position_local_eqn(l, k, i);
1597 
1598  /*IF it's not a boundary condition*/
1599  if (local_eqn >= 0)
1600  {
1601  // Initialise the contribution
1602  double sum = 0.0;
1603 
1604  // Acceleration and body force
1605  sum += (lambda_sq * accel[i] - b[i]) * psi(l, k);
1606 
1607  // Stress term
1608  for (unsigned a = 0; a < DIM; a++)
1609  {
1610  unsigned count = offset5;
1611  for (unsigned b = 0; b < DIM; b++)
1612  {
1613  // Add the stress terms to the residuals
1614  sum += sigma(a, b) * interpolated_G(a, i) *
1615  dpsidxi.raw_direct_access(count);
1616  ++count;
1617  }
1618  }
1619  residuals[local_eqn] += W * sum;
1620 
1621  // Add in the mass matrix terms
1622  if (flag > 2)
1623  {
1624  // Loop over the nodes of the element again
1625  for (unsigned ll = 0; ll < n_node; ll++)
1626  {
1627  // Loop of types of dofs again
1628  for (unsigned kk = 0; kk < n_position_type; kk++)
1629  {
1630  // Get the number of the unknown
1631  int local_unknown = this->position_local_eqn(ll, kk, i);
1632 
1633  /*IF it's not a boundary condition*/
1634  if (local_unknown >= 0)
1635  {
1636  mass_matrix(local_eqn, local_unknown) +=
1637  lambda_sq * psi(l, k) * psi(ll, kk) * W;
1638  }
1639  }
1640  }
1641  }
1642 
1643  // Add in the jacobian terms
1644  if ((flag == 1) || (flag == 3))
1645  {
1646  // Offset for faster access in general stress loop
1647  const unsigned offset1 = d_G_dX.offset(l, k, i);
1648 
1649  // Loop over the nodes of the element again
1650  for (unsigned ll = 0; ll < n_node; ll++)
1651  {
1652  // Loop of types of dofs again
1653  for (unsigned kk = 0; kk < n_position_type; kk++)
1654  {
1655  // Loop over the displacement components again
1656  for (unsigned ii = 0; ii < DIM; ii++)
1657  {
1658  // Get the number of the unknown
1659  int local_unknown = this->position_local_eqn(ll, kk, ii);
1660 
1661  /*IF it's not a boundary condition*/
1662  if (local_unknown >= 0)
1663  {
1664  // Offset for faster access in general stress loop
1665  const unsigned offset2 = d_G_dX.offset(ll, kk, ii);
1666  const unsigned offset4 = dpsidxi.offset(ll, kk);
1667 
1668  // General stress term
1669  //--------------------
1670  double sum = 0.0;
1671  unsigned count1 = offset1;
1672  for (unsigned a = 0; a < DIM; a++)
1673  {
1674  // Bump up direct access because we're only
1675  // accessing upper triangle
1676  count1 += a;
1677  for (unsigned b = a; b < DIM; b++)
1678  {
1679  double factor = d_G_dX.raw_direct_access(count1);
1680  if (a == b) factor *= 0.5;
1681 
1682  // Offset for faster access
1683  unsigned offset3 = d_stress_dG.offset(a, b);
1684  unsigned count2 = offset2;
1685  unsigned count3 = offset3;
1686 
1687  for (unsigned aa = 0; aa < DIM; aa++)
1688  {
1689  // Bump up direct access because we're only
1690  // accessing upper triangle
1691  count2 += aa;
1692  count3 += aa;
1693 
1694  // Only upper half of derivatives w.r.t. symm
1695  // tensor
1696  for (unsigned bb = aa; bb < DIM; bb++)
1697  {
1698  sum += factor *
1699  d_stress_dG.raw_direct_access(count3) *
1700  d_G_dX.raw_direct_access(count2);
1701  ++count2;
1702  ++count3;
1703  }
1704  }
1705  ++count1;
1706  }
1707  }
1708 
1709  // Multiply by weight and add contribution
1710  // (Add directly because this bit is nonsymmetric)
1711  jacobian(local_eqn, local_unknown) += sum * W;
1712 
1713  // Only upper triangle (no separate test for bc as
1714  // local_eqn is already nonnegative)
1715  if ((i == ii) && (local_unknown >= local_eqn))
1716  {
1717  // Initialise the contribution
1718  double sum = 0.0;
1719 
1720  // Inertia term
1721  sum +=
1722  lambda_sq * time_factor * psi(ll, kk) * psi(l, k);
1723 
1724  // Stress term
1725  unsigned count4 = offset4;
1726  for (unsigned a = 0; a < DIM; a++)
1727  {
1728  // Cache term
1729  const double factor =
1730  dpsidxi.raw_direct_access(count4); // ll ,kk
1731  ++count4;
1732 
1733  unsigned count5 = offset5;
1734  for (unsigned b = 0; b < DIM; b++)
1735  {
1736  sum += sigma(a, b) * factor *
1737  dpsidxi.raw_direct_access(count5); // l ,k
1738  ++count5;
1739  }
1740  }
1741 
1742  // Add to jacobian
1743  jacobian(local_eqn, local_unknown) += sum * W;
1744  // Add to lower triangular parts
1745  if (local_eqn != local_unknown)
1746  {
1747  jacobian(local_unknown, local_eqn) += sum * W;
1748  }
1749  }
1750  } // End of if not boundary condition
1751  }
1752  }
1753  }
1754  }
1755 
1756  // Derivatives w.r.t. pressure dofs
1757  if (flag > 0)
1758  {
1759  // Loop over the pressure dofs for unknowns
1760  for (unsigned l2 = 0; l2 < n_solid_pres; l2++)
1761  {
1762  local_unknown = this->solid_p_local_eqn(l2);
1763 
1764  // If it's not a boundary condition
1765  if (local_unknown >= 0)
1766  {
1767  // Add the pressure terms to the jacobian
1768  for (unsigned a = 0; a < DIM; a++)
1769  {
1770  for (unsigned b = 0; b < DIM; b++)
1771  {
1772  jacobian(local_eqn, local_unknown) -=
1773  psisp[l2] * Gup(a, b) * interpolated_G(a, i) *
1774  dpsidxi(l, k, b) * W;
1775  }
1776  }
1777  }
1778  }
1779  } // End of Jacobian terms
1780 
1781  } // End of if not boundary condition
1782  } // End of loop over coordinate directions
1783  } // End of loop over types of dof
1784  } // End of loop over shape functions
1785 
1786  //==============CONSTRAINT EQUATIONS FOR PRESSURE=====================
1787 
1788  // Now loop over the pressure degrees of freedom
1789  for (unsigned l = 0; l < n_solid_pres; l++)
1790  {
1791  local_eqn = this->solid_p_local_eqn(l);
1792 
1793  // Pinned (unlikely, actually) or real dof?
1794  if (local_eqn >= 0)
1795  {
1796  // For true incompressibility we need to conserve volume
1797  // so the determinant of the deformed metric tensor
1798  // needs to be equal to that of the undeformed one, which
1799  // is equal to the volumetric growth factor
1800  if (Incompressible)
1801  {
1802  residuals[local_eqn] += (detG - gamma) * psisp[l] * W;
1803 
1804 
1805  // Get Jacobian too?
1806  if ((flag == 1) || (flag == 3))
1807  {
1808  // Loop over the nodes of the element again
1809  for (unsigned ll = 0; ll < n_node; ll++)
1810  {
1811  // Loop of types of dofs again
1812  for (unsigned kk = 0; kk < n_position_type; kk++)
1813  {
1814  // Loop over the displacement components again
1815  for (unsigned ii = 0; ii < DIM; ii++)
1816  {
1817  // Get the number of the unknown
1818  int local_unknown = this->position_local_eqn(ll, kk, ii);
1819 
1820  /*IF it's not a boundary condition*/
1821  if (local_unknown >= 0)
1822  {
1823  // Offset for faster access
1824  const unsigned offset = d_G_dX.offset(ll, kk, ii);
1825 
1826  // General stress term
1827  double sum = 0.0;
1828  unsigned count = offset;
1829  for (unsigned aa = 0; aa < DIM; aa++)
1830  {
1831  // Bump up direct access because we're only
1832  // accessing upper triangle
1833  count += aa;
1834 
1835  // Only upper half
1836  for (unsigned bb = aa; bb < DIM; bb++)
1837  {
1838  sum += d_detG_dG(aa, bb) *
1839  d_G_dX.raw_direct_access(count) * psisp(l);
1840  ++count;
1841  }
1842  }
1843  jacobian(local_eqn, local_unknown) += sum * W;
1844  }
1845  }
1846  }
1847  }
1848 
1849  // No Jacobian terms due to pressure since it does not feature
1850  // in the incompressibility constraint
1851  }
1852  }
1853  // Nearly incompressible: (Neg.) pressure given by product of
1854  // bulk modulus and generalised dilatation
1855  else
1856  {
1857  residuals[local_eqn] +=
1858  (inv_kappa * interpolated_solid_p + gen_dil) * psisp[l] * W;
1859 
1860  // Add in the jacobian terms
1861  if ((flag == 1) || (flag == 3))
1862  {
1863  // Loop over the nodes of the element again
1864  for (unsigned ll = 0; ll < n_node; ll++)
1865  {
1866  // Loop of types of dofs again
1867  for (unsigned kk = 0; kk < n_position_type; kk++)
1868  {
1869  // Loop over the displacement components again
1870  for (unsigned ii = 0; ii < DIM; ii++)
1871  {
1872  // Get the number of the unknown
1873  int local_unknown = this->position_local_eqn(ll, kk, ii);
1874 
1875  /*IF it's not a boundary condition*/
1876  if (local_unknown >= 0)
1877  {
1878  // Offset for faster access
1879  const unsigned offset = d_G_dX.offset(ll, kk, ii);
1880 
1881  // General stress term
1882  double sum = 0.0;
1883  unsigned count = offset;
1884  for (unsigned aa = 0; aa < DIM; aa++)
1885  {
1886  // Bump up direct access because we're only
1887  // accessing upper triangle
1888  count += aa;
1889 
1890  // Only upper half
1891  for (unsigned bb = aa; bb < DIM; bb++)
1892  {
1893  sum += d_gen_dil_dG(aa, bb) *
1894  d_G_dX.raw_direct_access(count) * psisp(l);
1895  ++count;
1896  }
1897  }
1898  jacobian(local_eqn, local_unknown) += sum * W;
1899  }
1900  }
1901  }
1902  }
1903  }
1904 
1905  // Derivatives w.r.t. pressure dofs
1906  if (flag > 0)
1907  {
1908  // Loop over the pressure nodes again
1909  for (unsigned l2 = 0; l2 < n_solid_pres; l2++)
1910  {
1911  local_unknown = this->solid_p_local_eqn(l2);
1912  // If not pinnned
1913  if (local_unknown >= 0)
1914  {
1915  jacobian(local_eqn, local_unknown) +=
1916  inv_kappa * psisp[l2] * psisp[l] * W;
1917  }
1918  }
1919  } // End of jacobian terms
1920 
1921  } // End of else
1922 
1923  } // End of if not boundary condition
1924  }
1925 
1926  } // End of loop over integration points
1927  }
JacobiRotation< float > G
Definition: Jacobi_makeGivens.cpp:2
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
RowVector3d w
Definition: Matrix_resize_int.cpp:3
Scalar * b
Definition: benchVecAdd.cpp:17
virtual bool requires_incompressibility_constraint()=0
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nnodal_position_type() const
Definition: elements.h:2463
double dnodal_position_gen_dt(const unsigned &n, const unsigned &k, const unsigned &i) const
Definition: elements.h:2369
double nodal_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Definition: elements.h:2349
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
bool Unsteady
Flag that switches inertia on/off.
Definition: solid_elements.h:421
virtual int solid_p_local_eqn(const unsigned &i) const
Definition: solid_elements.h:178
ConstitutiveLaw * Constitutive_law_pt
Pointer to the constitutive law.
Definition: solid_elements.h:415
double prestress(const unsigned &i, const unsigned &j, const Vector< double > xi)
Definition: solid_elements.h:393
virtual unsigned npres_solid() const
Definition: solid_elements.h:171
void body_force(const Vector< double > &xi, Vector< double > &b) const
Definition: solid_elements.h:287
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
Definition: solid_elements.h:108
virtual double solid_p(const unsigned &l)=0
Return the lth solid pressure.
void get_d_stress_dG_upper(const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, const double &gen_dil, const double &inv_kappa, const double &interpolated_solid_p, RankFourTensor< double > &d_sigma_dG, DenseMatrix< double > &d_gen_dil_dG)
Definition: solid_elements.h:1177
void solid_pshape_at_knot(const unsigned &ipt, Shape &psi) const
Return the stored solid shape functions at the knots.
Definition: solid_elements.h:1220
double lagrangian_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Definition: elements.h:3912
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Definition: elements.h:4137
virtual void get_residuals_for_solid_ic(Vector< double > &residuals)
Definition: elements.h:4003
virtual double dshape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi) const
Definition: elements.cc:6737
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
const Scalar * a
Definition: level2_cplx_impl.h:32
char char char int int * k
Definition: level2_impl.h:374
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
int sigma
Definition: calibrate.py:179
@ W
Definition: quadtree.h:63

References a, b, Global_Parameters::body_force(), Constitutive::Constitutive_law_pt, DIM, G, mathsFunc::gamma(), i, J, j, k, oomph::RankFourTensor< T >::offset(), oomph::DShape::offset(), oomph::RankFiveTensor< T >::offset(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Eigen::bfloat16_impl::pow(), oomph::RankFourTensor< T >::raw_direct_access(), oomph::RankFiveTensor< T >::raw_direct_access(), oomph::DShape::raw_direct_access(), oomph::ConstitutiveLaw::requires_incompressibility_constraint(), s, calibrate::sigma, w, and oomph::QuadTreeNames::W.

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

◆ get_d_stress_dG_upper() [1/2]

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::get_d_stress_dG_upper ( const DenseMatrix< double > &  g,
const DenseMatrix< double > &  G,
const DenseMatrix< double > &  sigma,
const double detG,
const double interpolated_solid_p,
RankFourTensor< double > &  d_sigma_dG,
DenseMatrix< double > &  d_detG_dG 
)
inlineprotected

Return the derivative of the 2nd Piola Kirchhoff stress tensor, as calculated from the constitutive law in the incompresible formulation. Also return derivative of the determinant of the deformed covariant metric tensor (likely to be needed in the incompressibility constraint)

1295  {
1296 #ifdef PARANOID
1297  // If the pointer to the constitutive law hasn't been set, issue an error
1298  if (this->Constitutive_law_pt == 0)
1299  {
1300  // Write an error message
1301  std::string error_message =
1302  "Elements derived from PVDEquationsWithPressure \n";
1303  error_message += "must have a constitutive law:\n";
1304  error_message +=
1305  "set one using the constitutive_law_pt() member function";
1306  // Throw the error
1307  throw OomphLibError(
1309  }
1310 #endif
1311  // Only bother with the symmetric part by passing false as last entry
1313  g, G, sigma, detG, interpolated_solid_p, d_sigma_dG, d_detG_dG, false);
1314  }
virtual void calculate_d_second_piola_kirchhoff_stress_dG(const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, RankFourTensor< double > &d_sigma_dG, const bool &symmetrize_tensor=true)
Definition: constitutive_laws.cc:352

References oomph::ConstitutiveLaw::calculate_d_second_piola_kirchhoff_stress_dG(), oomph::PVDEquationsBase< DIM >::Constitutive_law_pt, G, oomph::PVDEquationsWithPressure< DIM >::interpolated_solid_p(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, calibrate::sigma, and oomph::Global_string_for_annotation::string().

◆ get_d_stress_dG_upper() [2/2]

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::get_d_stress_dG_upper ( const DenseMatrix< double > &  g,
const DenseMatrix< double > &  G,
const DenseMatrix< double > &  sigma,
const double gen_dil,
const double inv_kappa,
const double interpolated_solid_p,
RankFourTensor< double > &  d_sigma_dG,
DenseMatrix< double > &  d_gen_dil_dG 
)
inlineprotected

Return the derivative of the deviatoric part of the 2nd Piola Kirchhoff stress tensor, as calculated from the constitutive law in the nearly incompresible formulation. Also return the derivative of the generalised dilatation.

1186  {
1187 #ifdef PARANOID
1188  // If the pointer to the constitutive law hasn't been set, issue an error
1189  if (this->Constitutive_law_pt == 0)
1190  {
1191  // Write an error message
1192  std::string error_message =
1193  "Elements derived from PVDEquationsWithPressure \n";
1194  error_message += "must have a constitutive law:\n";
1195  error_message +=
1196  "set one using the constitutive_law_pt() member function";
1197  // Throw the error
1198  throw OomphLibError(
1200  }
1201 #endif
1202  // Only bother with the symmetric part by passing false as last entry
1204  g,
1205  G,
1206  sigma,
1207  gen_dil,
1208  inv_kappa,
1210  d_sigma_dG,
1211  d_gen_dil_dG,
1212  false);
1213  }

References oomph::ConstitutiveLaw::calculate_d_second_piola_kirchhoff_stress_dG(), oomph::PVDEquationsBase< DIM >::Constitutive_law_pt, G, oomph::PVDEquationsWithPressure< DIM >::interpolated_solid_p(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, calibrate::sigma, and oomph::Global_string_for_annotation::string().

◆ get_dof_numbers_for_unknowns()

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::get_dof_numbers_for_unknowns ( std::list< std::pair< unsigned long, unsigned >> &  dof_lookup_list) const
inlinevirtual

Create a list of pairs for all unknowns in this element, so that the first entry in each pair contains the global equation number of the unknown, while the second one contains the number of the "DOF" that this unknown is associated with. (Function can obviously only be called if the equation numbering scheme has been set up.) There are DIM+1 types of DOF: displacement compnents and pressure

Reimplemented from oomph::PVDEquationsBase< DIM >.

1079  {
1080  // temporary pair (used to store dof lookup prior to being added to list
1081  std::pair<unsigned, unsigned> dof_lookup;
1082 
1083  // number of nodes
1084  const unsigned n_node = this->nnode();
1085 
1086  // Get the number of position dofs and dimensions at the node
1087  const unsigned n_position_type = this->nnodal_position_type();
1088  const unsigned nodal_dim = this->nodal_dimension();
1089 
1090  // Integer storage for local unknown
1091  int local_unknown = 0;
1092 
1093  // Loop over the nodes
1094  for (unsigned n = 0; n < n_node; n++)
1095  {
1096  // Loop over position dofs
1097  for (unsigned k = 0; k < n_position_type; k++)
1098  {
1099  // Loop over dimension
1100  for (unsigned i = 0; i < nodal_dim; i++)
1101  {
1102  // If the variable is free
1103  local_unknown = this->position_local_eqn(n, k, i);
1104 
1105  // ignore pinned values
1106  if (local_unknown >= 0)
1107  {
1108  // store dof lookup in temporary pair: First entry in pair
1109  // is global equation number; second entry is dof type
1110  dof_lookup.first = this->eqn_number(local_unknown);
1111  dof_lookup.second = i;
1112 
1113  // add to list
1114  dof_lookup_list.push_front(dof_lookup);
1115  }
1116  }
1117  }
1118  }
1119 
1120  // Do solid pressure degrees of freedom
1121  unsigned np = this->npres_solid();
1122  for (unsigned j = 0; j < np; j++)
1123  {
1124  int local_unknown = this->solid_p_local_eqn(j);
1125  // ignore pinned values
1126  if (local_unknown >= 0)
1127  {
1128  // store dof lookup in temporary pair: First entry in pair
1129  // is global equation number; second entry is dof type
1130  dof_lookup.first = this->eqn_number(local_unknown);
1131  dof_lookup.second = DIM;
1132 
1133  // add to list
1134  dof_lookup_list.push_front(dof_lookup);
1135  }
1136  }
1137  }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
unsigned long eqn_number(const unsigned &ieqn_local) const
Definition: elements.h:704

References DIM, oomph::GeneralisedElement::eqn_number(), i, j, k, n, oomph::FiniteElement::nnodal_position_type(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_dimension(), oomph::PVDEquationsBase< DIM >::npres_solid(), oomph::SolidFiniteElement::position_local_eqn(), and oomph::PVDEquationsBase< DIM >::solid_p_local_eqn().

◆ get_mass_matrix_diagonal()

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::get_mass_matrix_diagonal ( Vector< double > &  mass_diag)
virtual

Compute the diagonal of the displacement mass matrix for LSC preconditioner

Compute the diagonal of the velocity mass matrix for LSC preconditioner.

Implements oomph::SolidElementWithDiagonalMassMatrix.

2200  {
2201  // Resize and initialise
2202  mass_diag.assign(this->ndof(), 0.0);
2203 
2204  // find out how many nodes there are
2205  unsigned n_node = this->nnode();
2206 
2207  // Find out how many position types of dof there are
2208  const unsigned n_position_type = this->nnodal_position_type();
2209 
2210  // Set up memory for the shape functions
2211  Shape psi(n_node, n_position_type);
2212  DShape dpsidxi(n_node, n_position_type, DIM);
2213 
2214  // Number of integration points
2215  unsigned n_intpt = this->integral_pt()->nweight();
2216 
2217  // Integer to store the local equations no
2218  int local_eqn = 0;
2219 
2220  // Loop over the integration points
2221  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2222  {
2223  // Get the integral weight
2224  double w = this->integral_pt()->weight(ipt);
2225 
2226  // Call the derivatives of the shape functions
2227  double J = this->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
2228 
2229  // Premultiply weights and Jacobian
2230  double W = w * J;
2231 
2232  // Loop over the nodes
2233  for (unsigned l = 0; l < n_node; l++)
2234  {
2235  // Loop over the types of dof
2236  for (unsigned k = 0; k < n_position_type; k++)
2237  {
2238  // Loop over the directions
2239  for (unsigned i = 0; i < DIM; i++)
2240  {
2241  // Get the equation number
2242  local_eqn = this->position_local_eqn(l, k, i);
2243 
2244  // If not a boundary condition
2245  if (local_eqn >= 0)
2246  {
2247  // Add the contribution
2248  mass_diag[local_eqn] += pow(psi(l, k), 2) * W;
2249  } // End of if not boundary condition statement
2250  } // End of loop over dimension
2251  } // End of dof type
2252  } // End of loop over basis functions
2253  }
2254  }

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

◆ get_stress() [1/3]

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::get_stress ( const DenseMatrix< double > &  g,
const DenseMatrix< double > &  G,
DenseMatrix< double > &  sigma_dev,
DenseMatrix< double > &  Gcontra,
double detG 
)
inlineprotected

Return the deviatoric part of the 2nd Piola Kirchhoff stress tensor, as calculated from the constitutive law in the incompresible formulation. Also return the contravariant deformed metric tensor, and the determinant of the deformed covariant metric tensor (likely to be needed in the incompressibility constraint)

1263  {
1264 #ifdef PARANOID
1265  // If the pointer to the constitutive law hasn't been set, issue an error
1266  if (this->Constitutive_law_pt == 0)
1267  {
1268  // Write an error message
1269  std::string error_message =
1270  "Elements derived from PVDEquationsWithPressure \n";
1271  error_message += "must have a constitutive law:\n";
1272  error_message +=
1273  "set one using the constitutive_law_pt() member function";
1274  // Throw the error
1275  throw OomphLibError(
1277  }
1278 #endif
1280  g, G, sigma_dev, Gcontra, detG);
1281  }
virtual void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)=0

References oomph::ConstitutiveLaw::calculate_second_piola_kirchhoff_stress(), oomph::PVDEquationsBase< DIM >::Constitutive_law_pt, G, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::Global_string_for_annotation::string().

◆ get_stress() [2/3]

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::get_stress ( const DenseMatrix< double > &  g,
const DenseMatrix< double > &  G,
DenseMatrix< double > &  sigma_dev,
DenseMatrix< double > &  Gcontra,
double gen_dil,
double inv_kappa 
)
inlineprotected

Return the deviatoric part of the 2nd Piola Kirchhoff stress tensor, as calculated from the constitutive law in the nearly incompresible formulation. Also return the contravariant deformed metric tensor, the generalised dilatation, and the inverse of the bulk modulus.

1151  {
1152 #ifdef PARANOID
1153  // If the pointer to the constitutive law hasn't been set, issue an error
1154  if (this->Constitutive_law_pt == 0)
1155  {
1156  // Write an error message
1157  std::string error_message =
1158  "Elements derived from PVDEquationsWithPressure \n";
1159  error_message += "must have a constitutive law:\n";
1160  error_message +=
1161  "set one using the constitutive_law_pt() member function";
1162  // Throw the error
1163  throw OomphLibError(
1165  }
1166 #endif
1168  g, G, sigma_dev, Gcontra, gen_dil, inv_kappa);
1169  }

References oomph::ConstitutiveLaw::calculate_second_piola_kirchhoff_stress(), oomph::PVDEquationsBase< DIM >::Constitutive_law_pt, G, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::Global_string_for_annotation::string().

◆ get_stress() [3/3]

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::get_stress ( const Vector< double > &  s,
DenseMatrix< double > &  sigma 
)
virtual

Return the 2nd Piola Kirchoff stress tensor, as calculated from the constitutive law at specified local coordinate

Compute the contravariant second Piola Kirchoff stress at a given local coordinate. Note: this replicates a lot of code that is already coontained in get_residuals() but without sacrificing efficiency (re-computing the shape functions several times) or creating helper functions with horrendous interfaces (to pass all the functions which shouldn't be recomputed) about this is unavoidable.

Implements oomph::PVDEquationsBase< DIM >.

2269  {
2270  // Find out how many nodes there are
2271  unsigned n_node = this->nnode();
2272 
2273  // Find out how many positional dofs there are
2274  unsigned n_position_type = this->nnodal_position_type();
2275 
2276  // Find out how many pressure dofs there are
2277  unsigned n_solid_pres = this->npres_solid();
2278 
2279  // Set up memory for the shape functions
2280  Shape psi(n_node, n_position_type);
2281  DShape dpsidxi(n_node, n_position_type, DIM);
2282 
2283  // Set up memory for the pressure shape functions
2284  Shape psisp(n_solid_pres);
2285 
2286  // Find values of shape function
2287  solid_pshape(s, psisp);
2288 
2289  // Call the derivatives of the shape functions (ignore Jacobian)
2290  (void)this->dshape_lagrangian(s, psi, dpsidxi);
2291 
2292  // Lagrangian coordinates
2293  Vector<double> xi(DIM);
2294  this->interpolated_xi(s, xi);
2295 
2296  // Get isotropic growth factor
2297  double gamma;
2298  // Dummy integration point
2299  unsigned ipt = 0;
2300  this->get_isotropic_growth(ipt, s, xi, gamma);
2301 
2302  // We use Cartesian coordinates as the reference coordinate
2303  // system. In this case the undeformed metric tensor is always
2304  // the identity matrix -- stretched by the isotropic growth
2305  double diag_entry = pow(gamma, 2.0 / double(DIM));
2307  for (unsigned i = 0; i < DIM; i++)
2308  {
2309  for (unsigned j = 0; j < DIM; j++)
2310  {
2311  if (i == j)
2312  {
2313  g(i, j) = diag_entry;
2314  }
2315  else
2316  {
2317  g(i, j) = 0.0;
2318  }
2319  }
2320  }
2321 
2322 
2323  // Calculate interpolated values of the derivative of global position
2324  // wrt lagrangian coordinates
2325  DenseMatrix<double> interpolated_G(DIM);
2326 
2327  // Initialise to zero
2328  for (unsigned i = 0; i < DIM; i++)
2329  {
2330  for (unsigned j = 0; j < DIM; j++)
2331  {
2332  interpolated_G(i, j) = 0.0;
2333  }
2334  }
2335 
2336  // Calculate displacements and derivatives
2337  for (unsigned l = 0; l < n_node; l++)
2338  {
2339  // Loop over positional dofs
2340  for (unsigned k = 0; k < n_position_type; k++)
2341  {
2342  // Loop over displacement components (deformed position)
2343  for (unsigned i = 0; i < DIM; i++)
2344  {
2345  // Loop over derivative directions
2346  for (unsigned j = 0; j < DIM; j++)
2347  {
2348  interpolated_G(j, i) +=
2349  this->nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
2350  }
2351  }
2352  }
2353  }
2354 
2355  // Declare and calculate the deformed metric tensor
2357 
2358  // Assign values of G
2359  for (unsigned i = 0; i < DIM; i++)
2360  {
2361  // Do upper half of matrix
2362  // Note that j must be signed here for the comparison test to work
2363  // Also i must be cast to an int
2364  for (int j = (DIM - 1); j >= static_cast<int>(i); j--)
2365  {
2366  // Initialise G(i,j) to zero
2367  G(i, j) = 0.0;
2368  // Now calculate the dot product
2369  for (unsigned k = 0; k < DIM; k++)
2370  {
2371  G(i, j) += interpolated_G(i, k) * interpolated_G(j, k);
2372  }
2373  }
2374  // Matrix is symmetric so just copy lower half
2375  for (int j = (i - 1); j >= 0; j--)
2376  {
2377  G(i, j) = G(j, i);
2378  }
2379  }
2380 
2381 
2382  // Calculate the interpolated solid pressure
2383  double interpolated_solid_p = 0.0;
2384  for (unsigned l = 0; l < n_solid_pres; l++)
2385  {
2386  interpolated_solid_p += solid_p(l) * psisp[l];
2387  }
2388 
2389  // Now calculate the deviatoric stress and all pressure-related
2390  // quantitites
2391  DenseMatrix<double> sigma_dev(DIM), Gup(DIM);
2392  double detG = 0.0;
2393  double gen_dil = 0.0;
2394  double inv_kappa = 0.0;
2395 
2396  // Incompressible: Compute the deviatoric part of the stress tensor, the
2397  // contravariant deformed metric tensor and the determinant
2398  // of the deformed covariant metric tensor.
2399 
2400  if (Incompressible)
2401  {
2402  get_stress(g, G, sigma_dev, Gup, detG);
2403  }
2404  // Nearly incompressible: Compute the deviatoric part of the
2405  // stress tensor, the contravariant deformed metric tensor,
2406  // the generalised dilatation and the inverse bulk modulus.
2407  else
2408  {
2409  get_stress(g, G, sigma_dev, Gup, gen_dil, inv_kappa);
2410  }
2411 
2412  // Get complete stress
2413  for (unsigned i = 0; i < DIM; i++)
2414  {
2415  for (unsigned j = 0; j < DIM; j++)
2416  {
2417  sigma(i, j) = -interpolated_solid_p * Gup(i, j) + sigma_dev(i, j);
2418  }
2419  }
2420  }
virtual void solid_pshape(const Vector< double > &s, Shape &psi) const =0
Return the solid pressure shape functions.
double dshape_lagrangian(const Vector< double > &s, Shape &psi, DShape &dpsidxi) const
Definition: elements.cc:6710

References DIM, G, mathsFunc::gamma(), i, j, k, Eigen::bfloat16_impl::pow(), s, and calibrate::sigma.

◆ interpolated_solid_p()

template<unsigned DIM>
double oomph::PVDEquationsWithPressure< DIM >::interpolated_solid_p ( const Vector< double > &  s)
inline

Return the interpolated_solid_pressure.

1012  {
1013  // Find number of nodes
1014  unsigned n_solid_pres = this->npres_solid();
1015  // Local shape function
1016  Shape psisp(n_solid_pres);
1017  // Find values of shape function
1018  solid_pshape(s, psisp);
1019 
1020  // Initialise value of solid_p
1021  double interpolated_solid_p = 0.0;
1022  // Loop over the local nodes and sum
1023  for (unsigned l = 0; l < n_solid_pres; l++)
1024  {
1025  interpolated_solid_p += solid_p(l) * psisp[l];
1026  }
1027 
1028  return (interpolated_solid_p);
1029  }

References oomph::PVDEquationsBase< DIM >::npres_solid(), s, oomph::PVDEquationsWithPressure< DIM >::solid_p(), and oomph::PVDEquationsWithPressure< DIM >::solid_pshape().

Referenced by oomph::RefineableQPVDElementWithPressure< DIM >::further_build(), oomph::PVDEquationsWithPressure< DIM >::get_d_stress_dG_upper(), and oomph::RefineableQPVDElementWithContinuousPressure< DIM >::get_interpolated_values().

◆ is_incompressible()

template<unsigned DIM>
bool oomph::PVDEquationsWithPressure< DIM >::is_incompressible ( ) const
inline

Return whether the material is incompressible.

876  {
877  return Incompressible;
878  }

References oomph::PVDEquationsWithPressure< DIM >::Incompressible.

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

◆ ndof_types()

template<unsigned DIM>
unsigned oomph::PVDEquationsWithPressure< DIM >::ndof_types ( ) const
inlinevirtual

returns the number of DOF types associated with this element: displacement components and pressure

Reimplemented from oomph::PVDEquationsBase< DIM >.

1065  {
1066  return DIM + 1;
1067  }

References DIM.

◆ output() [1/4]

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

C-style output: x,y,[z],xi0,xi1,[xi2],p,gamma.

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::TPVDElementWithContinuousPressure< DIM >, oomph::QPVDElementWithContinuousPressure< DIM >, and oomph::QPVDElementWithPressure< DIM >.

1045  {
1046  unsigned n_plot = 5;
1047  output(file_pt, n_plot);
1048  }
void output(std::ostream &outfile)
Output: x,y,[z],xi0,xi1,[xi2],p,gamma.
Definition: solid_elements.h:1033

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

◆ output() [2/4]

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

C-style output: x,y,[z],xi0,xi1,[xi2],p,gamma.

C-stsyle output: x,y,[z],xi0,xi1,[xi2],p,gamma.

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::TPVDElementWithContinuousPressure< DIM >, oomph::QPVDElementWithContinuousPressure< DIM >, and oomph::QPVDElementWithPressure< DIM >.

1992  {
1993  // Set output Vector
1994  Vector<double> s(DIM);
1995  Vector<double> x(DIM);
1996  Vector<double> xi(DIM);
1997 
1998  switch (DIM)
1999  {
2000  case 2:
2001  // Tecplot header info
2002  // outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
2003  fprintf(file_pt, "ZONE I=%i, J=%i\n", n_plot, n_plot);
2004 
2005  // Loop over element nodes
2006  for (unsigned l2 = 0; l2 < n_plot; l2++)
2007  {
2008  s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
2009  for (unsigned l1 = 0; l1 < n_plot; l1++)
2010  {
2011  s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
2012 
2013  // Get Eulerian and Lagrangian coordinates
2014  this->interpolated_x(s, x);
2015  this->interpolated_xi(s, xi);
2016 
2017  // Get isotropic growth
2018  double gamma;
2019  // Dummy integration point
2020  unsigned ipt = 0;
2021  this->get_isotropic_growth(ipt, s, xi, gamma);
2022 
2023  // Output the x,y,..
2024  for (unsigned i = 0; i < DIM; i++)
2025  {
2026  // outfile << x[i] << " ";
2027  fprintf(file_pt, "%g ", x[i]);
2028  }
2029  // Output xi0,xi1,..
2030  for (unsigned i = 0; i < DIM; i++)
2031  {
2032  // outfile << xi[i] << " ";
2033  fprintf(file_pt, "%g ", xi[i]);
2034  }
2035  // Output growth
2036  // outfile << gamma << " ";
2037  fprintf(file_pt, "%g ", gamma);
2038 
2039  // Output pressure
2040  // outfile << interpolated_solid_p(s) << " ";
2041  // outfile << std::endl;
2042  fprintf(file_pt, "%g \n", interpolated_solid_p(s));
2043  }
2044  }
2045 
2046  break;
2047 
2048  case 3:
2049  // Tecplot header info
2050  // outfile << "ZONE I=" << n_plot
2051  // << ", J=" << n_plot
2052  // << ", K=" << n_plot << std::endl;
2053  fprintf(file_pt, "ZONE I=%i, J=%i, K=%i \n", n_plot, n_plot, n_plot);
2054 
2055  // Loop over element nodes
2056  for (unsigned l3 = 0; l3 < n_plot; l3++)
2057  {
2058  s[2] = -1.0 + l3 * 2.0 / (n_plot - 1);
2059  for (unsigned l2 = 0; l2 < n_plot; l2++)
2060  {
2061  s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
2062  for (unsigned l1 = 0; l1 < n_plot; l1++)
2063  {
2064  s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
2065 
2066  // Get Eulerian and Lagrangian coordinates
2067  this->interpolated_x(s, x);
2068  this->interpolated_xi(s, xi);
2069 
2070  // Get isotropic growth
2071  double gamma;
2072  // Dummy integration point
2073  unsigned ipt = 0;
2074  this->get_isotropic_growth(ipt, s, xi, gamma);
2075 
2076  // Output the x,y,..
2077  for (unsigned i = 0; i < DIM; i++)
2078  {
2079  // outfile << x[i] << " ";
2080  fprintf(file_pt, "%g ", x[i]);
2081  }
2082  // Output xi0,xi1,..
2083  for (unsigned i = 0; i < DIM; i++)
2084  {
2085  // outfile << xi[i] << " ";
2086  fprintf(file_pt, "%g ", xi[i]);
2087  }
2088  // Output growth
2089  // outfile << gamma << " ";
2090  fprintf(file_pt, "%g ", gamma);
2091 
2092  // Output pressure
2093  // outfile << interpolated_solid_p(s) << " ";
2094  // outfile << std::endl;
2095  fprintf(file_pt, "%g \n", interpolated_solid_p(s));
2096  }
2097  }
2098  }
2099  break;
2100 
2101  default:
2102  std::ostringstream error_message;
2103  error_message << "No output routine for PVDEquationsWithPressure<"
2104  << DIM << "> elements. Write it yourself!" << std::endl;
2105  throw OomphLibError(error_message.str(),
2108  }
2109  }

References DIM, mathsFunc::gamma(), i, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, and plotDoE::x.

◆ output() [3/4]

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::output ( std::ostream &  outfile)
inlinevirtual

◆ output() [4/4]

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

Output: x,y,[z],xi0,xi1,[xi2],p,gamma.

Reimplemented from oomph::FiniteElement.

Reimplemented in oomph::TPVDElementWithContinuousPressure< DIM >, oomph::QPVDElementWithContinuousPressure< DIM >, and oomph::QPVDElementWithPressure< DIM >.

1936  {
1937  // Set output Vector
1938  Vector<double> s(DIM);
1939  Vector<double> x(DIM);
1940  Vector<double> xi(DIM);
1941 
1942  // Tecplot header info
1943  outfile << this->tecplot_zone_string(n_plot);
1944 
1945  // Loop over plot points
1946  unsigned num_plot_points = this->nplot_points(n_plot);
1947  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1948  {
1949  // Get local coordinates of plot point
1950  this->get_s_plot(iplot, n_plot, s);
1951 
1952  // Get Eulerian and Lagrangian coordinates
1953  this->interpolated_x(s, x);
1954  this->interpolated_xi(s, xi);
1955 
1956  // Get isotropic growth
1957  double gamma;
1958  // Dummy integration point
1959  unsigned ipt = 0;
1960  this->get_isotropic_growth(ipt, s, xi, gamma);
1961 
1962  // Output the x,y,..
1963  for (unsigned i = 0; i < DIM; i++)
1964  {
1965  outfile << x[i] << " ";
1966  }
1967 
1968  // Output xi0,xi1,..
1969  for (unsigned i = 0; i < DIM; i++)
1970  {
1971  outfile << xi[i] << " ";
1972  }
1973 
1974  // Output growth
1975  outfile << gamma << " ";
1976  // Output pressure
1977  outfile << interpolated_solid_p(s) << " ";
1978  outfile << "\n";
1979  }
1980 
1981  // Write tecplot footer (e.g. FE connectivity lists)
1982  this->write_tecplot_zone_footer(outfile, n_plot);
1983  }

References DIM, mathsFunc::gamma(), i, s, and plotDoE::x.

◆ set_compressible()

◆ set_incompressible()

◆ set_solid_p()

template<unsigned DIM>
virtual void oomph::PVDEquationsWithPressure< DIM >::set_solid_p ( const unsigned l,
const double p_value 
)
pure virtual

◆ solid_p()

◆ solid_pshape()

◆ solid_pshape_at_knot()

template<unsigned DIM>
void oomph::PVDEquationsWithPressure< DIM >::solid_pshape_at_knot ( const unsigned ipt,
Shape psi 
) const
inlineprotected

Return the stored solid shape functions at the knots.

1221  {
1222  // Find the dimension of the element
1223  unsigned Dim = this->dim();
1224  // Storage for local coordinates of the integration point
1225  Vector<double> s(Dim);
1226  // Set the local coordinates
1227  for (unsigned i = 0; i < Dim; i++)
1228  {
1229  s[i] = this->integral_pt()->knot(ipt, i);
1230  }
1231  // Get the shape function
1232  solid_pshape(s, psi);
1233  }
unsigned dim() const
Definition: elements.h:2611
static const unsigned Dim
Problem dimension.
Definition: two_d_tilted_square.cc:62

References Global_Variables::Dim, oomph::FiniteElement::dim(), i, oomph::FiniteElement::integral_pt(), oomph::Integral::knot(), s, and oomph::PVDEquationsWithPressure< DIM >::solid_pshape().

Member Data Documentation

◆ Incompressible


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