Public Types | Public Member Functions | Protected Member Functions | Protected Attributes | Private Attributes | Static Private Attributes | List of all members
oomph::PoroelasticityEquations< DIM > Class Template Referenceabstract

#include <poroelasticity_elements.h>

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

Public Types

typedef void(* SourceFctPt) (const double &time, const Vector< double > &x, Vector< double > &f)
 Source function pointer typedef. More...
 
typedef void(* MassSourceFctPt) (const double &time, const Vector< double > &x, double &f)
 Mass source function pointer typedef. More...
 
- Public Types inherited from oomph::FiniteElement
typedef void(* SteadyExactSolutionFctPt) (const Vector< double > &, Vector< double > &)
 
typedef void(* UnsteadyExactSolutionFctPt) (const double &, const Vector< double > &, Vector< double > &)
 

Public Member Functions

 PoroelasticityEquations ()
 Constructor. More...
 
const double & lambda_sq () const
 Access function for timescale ratio (nondim density) More...
 
double *& lambda_sq_pt ()
 Access function for pointer to timescale ratio (nondim density) More...
 
const double & density_ratio () const
 Access function for the density ratio. More...
 
double *& density_ratio_pt ()
 Access function for pointer to the density ratio. More...
 
const double & k_inv () const
 Access function for the nondim inverse permeability. More...
 
double *& k_inv_pt ()
 Access function for pointer to the nondim inverse permeability. More...
 
const double & alpha () const
 Access function for alpha. More...
 
double *& alpha_pt ()
 Access function for pointer to alpha. More...
 
const double & porosity () const
 Access function for the porosity. More...
 
double *& porosity_pt ()
 Access function for pointer to the porosity. More...
 
SourceFctPt & force_solid_fct_pt ()
 Access function: Pointer to solid force function. More...
 
SourceFctPt force_solid_fct_pt () const
 Access function: Pointer to solid force function (const version) More...
 
SourceFctPt & force_fluid_fct_pt ()
 Access function: Pointer to fluid force function. More...
 
SourceFctPt force_fluid_fct_pt () const
 Access function: Pointer to fluid force function (const version) More...
 
MassSourceFctPt & mass_source_fct_pt ()
 Access function: Pointer to mass source function. More...
 
MassSourceFctPt mass_source_fct_pt () const
 Access function: Pointer to mass source function (const version) More...
 
void force_solid (const double &time, const Vector< double > &x, Vector< double > &b) const
 
void force_fluid (const double &time, const Vector< double > &x, Vector< double > &b) const
 
void mass_source (const double &time, const Vector< double > &x, double &b) const
 
ElasticityTensor *& elasticity_tensor_pt ()
 Return the pointer to the elasticity_tensor. More...
 
const double E (const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l) const
 Access function to the entries in the elasticity tensor. More...
 
void get_stress (const Vector< double > &s, DenseMatrix< double > &sigma) const
 
void get_strain (const Vector< double > &s, DenseMatrix< double > &strain) const
 Return the strain tensor. More...
 
virtual unsigned required_nvalue (const unsigned &n) const =0
 Number of values required at node n. More...
 
virtual unsigned u_index (const unsigned &n) const =0
 Return the nodal index of the n-th solid displacement unknown. More...
 
virtual int q_edge_local_eqn (const unsigned &n) const =0
 Return the equation number of the n-th edge (flux) degree of freedom. More...
 
virtual int q_internal_local_eqn (const unsigned &n) const =0
 
virtual unsigned q_edge_index (const unsigned &n) const =0
 Return the nodal index at which the nth edge unknown is stored. More...
 
virtual unsigned q_internal_index () const =0
 
virtual unsigned q_edge_node_number (const unsigned &n) const =0
 Return the number of the node where the nth edge unknown is stored. More...
 
virtual double q_edge (const unsigned &n) const =0
 Return the values of the edge (flux) degrees of freedom. More...
 
virtual double q_edge (const unsigned &t, const unsigned &n) const =0
 
virtual double q_internal (const unsigned &n) const =0
 Return the values of the internal (moment) degrees of freedom. More...
 
virtual double q_internal (const unsigned &t, const unsigned &n) const =0
 
virtual unsigned nq_basis () const =0
 Return the total number of computational basis functions for q. More...
 
virtual unsigned nq_basis_edge () const =0
 Return the number of edge basis functions for q. More...
 
virtual void get_q_basis_local (const Vector< double > &s, Shape &q_basis) const =0
 Returns the local form of the q basis at local coordinate s. More...
 
virtual void get_div_q_basis_local (const Vector< double > &s, Shape &div_q_basis_ds) const =0
 
void get_q_basis (const Vector< double > &s, Shape &q_basis) const
 Returns the transformed basis at local coordinate s. More...
 
virtual unsigned nedge_gauss_point () const =0
 Returns the number of gauss points along each edge of the element. More...
 
virtual double edge_gauss_point (const unsigned &edge, const unsigned &n) const =0
 Returns the nth gauss point along an edge. More...
 
virtual void edge_gauss_point_global (const unsigned &edge, const unsigned &n, Vector< double > &x) const =0
 Returns the global coordinates of the nth gauss point along an edge. More...
 
virtual void pin_q_internal_value (const unsigned &n)=0
 Pin the nth internal q value. More...
 
virtual int p_local_eqn (const unsigned &n) const =0
 Return the equation number of the n-th pressure degree of freedom. More...
 
virtual double p_value (unsigned &n) const =0
 Return the nth pressure value. More...
 
virtual unsigned np_basis () const =0
 Return the total number of pressure basis functions. More...
 
virtual void get_p_basis (const Vector< double > &s, Shape &p_basis) const =0
 Return the pressure basis. More...
 
virtual void pin_p_value (const unsigned &n, const double &p)=0
 Pin the nth pressure value. More...
 
virtual void scale_basis (Shape &basis) const =0
 Scale the edge basis to allow arbitrary edge mappings. More...
 
double transform_basis (const Vector< double > &s, const Shape &q_basis_local, Shape &psi, DShape &dpsi, Shape &q_basis) const
 
double transform_basis (const Vector< double > &s, const Shape &q_basis_local, Shape &psi, Shape &q_basis) const
 
void fill_in_contribution_to_residuals (Vector< double > &residuals)
 Fill in contribution to residuals for the Darcy equations. More...
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 Fill in the Jacobian matrix for the Newton method. More...
 
void interpolated_u (const Vector< double > &s, Vector< double > &disp) const
 Calculate the FE representation of u. More...
 
double interpolated_u (const Vector< double > &s, const unsigned &i) const
 Calculate the FE representation of the i-th component of u. More...
 
void interpolated_q (const Vector< double > &s, Vector< double > &u) const
 Calculate the FE representation of q. More...
 
double interpolated_q (const Vector< double > &s, const unsigned i) const
 Calculate the FE representation of the i-th component of q. More...
 
void interpolated_div_q (const Vector< double > &s, double &div_q) const
 Calculate the FE representation of div u. More...
 
double interpolated_div_q (const Vector< double > &s)
 Calculate the FE representation of div q and return it. More...
 
void interpolated_p (const Vector< double > &s, double &p) const
 Calculate the FE representation of p. More...
 
double interpolated_p (const Vector< double > &s) const
 Calculate the FE representation of p and return it. More...
 
double du_dt (const unsigned &n, const unsigned &i) const
 du/dt at local node n More...
 
double d2u_dt2 (const unsigned &n, const unsigned &i) const
 d^2u/dt^2 at local node n More...
 
double dq_edge_dt (const unsigned &n) const
 dq_edge/dt for the n-th edge degree of freedom More...
 
double dq_internal_dt (const unsigned &n) const
 dq_internal/dt for the n-th internal degree of freedom More...
 
void set_q_internal_timestepper (TimeStepper *const time_stepper_pt)
 Set the timestepper of the q internal data object. More...
 
unsigned self_test ()
 
void output (std::ostream &outfile)
 Output with default number of plot points. More...
 
void output (std::ostream &outfile, const unsigned &nplot)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
 
void output_fct (std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
 
void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
- Public Member Functions inherited from oomph::FiniteElement
void set_dimension (const unsigned &dim)
 
void set_nodal_dimension (const unsigned &nodal_dim)
 
void set_nnodal_position_type (const unsigned &nposition_type)
 Set the number of types required to interpolate the coordinate. More...
 
void set_n_node (const unsigned &n)
 
int nodal_local_eqn (const unsigned &n, const unsigned &i) const
 
double dJ_eulerian_at_knot (const unsigned &ipt, Shape &psi, DenseMatrix< double > &djacobian_dX) const
 
 FiniteElement ()
 Constructor. More...
 
virtual ~FiniteElement ()
 
 FiniteElement (const FiniteElement &)=delete
 Broken copy constructor. More...
 
virtual bool local_coord_is_valid (const Vector< double > &s)
 Broken assignment operator. More...
 
virtual void move_local_coord_back_into_element (Vector< double > &s) const
 
void get_centre_of_gravity_and_max_radius_in_terms_of_zeta (Vector< double > &cog, double &max_radius) const
 
virtual void local_coordinate_of_node (const unsigned &j, Vector< double > &s) const
 
virtual void local_fraction_of_node (const unsigned &j, Vector< double > &s_fraction)
 
virtual double local_one_d_fraction_of_node (const unsigned &n1d, const unsigned &i)
 
virtual void set_macro_elem_pt (MacroElement *macro_elem_pt)
 
MacroElement * macro_elem_pt ()
 Access function to pointer to macro element. More...
 
void get_x (const Vector< double > &s, Vector< double > &x) const
 
void get_x (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void get_x_from_macro_element (const Vector< double > &s, Vector< double > &x) const
 
virtual void get_x_from_macro_element (const unsigned &t, const Vector< double > &s, Vector< double > &x)
 
virtual void set_integration_scheme (Integral *const &integral_pt)
 Set the spatial integration scheme. More...
 
Integral *const & integral_pt () const
 Return the pointer to the integration scheme (const version) More...
 
virtual void shape (const Vector< double > &s, Shape &psi) const =0
 
virtual void shape_at_knot (const unsigned &ipt, Shape &psi) const
 
virtual void dshape_local (const Vector< double > &s, Shape &psi, DShape &dpsids) const
 
virtual void dshape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids) const
 
virtual void d2shape_local (const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual void d2shape_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
 
virtual double J_eulerian (const Vector< double > &s) const
 
virtual double J_eulerian_at_knot (const unsigned &ipt) const
 
void check_J_eulerian_at_knots (bool &passed) const
 
void check_jacobian (const double &jacobian) const
 
double dshape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx) const
 
virtual double dshape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, DenseMatrix< double > &djacobian_dX, RankFourTensor< double > &d_dpsidx_dX) const
 
double d2shape_eulerian (const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual double d2shape_eulerian_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
 
virtual void assign_nodal_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void describe_nodal_local_dofs (std::ostream &out, const std::string &current_string) const
 
virtual void assign_all_generic_local_eqn_numbers (const bool &store_local_dof_pt)
 
Node *& node_pt (const unsigned &n)
 Return a pointer to the local node n. More...
 
Node *const & node_pt (const unsigned &n) const
 Return a pointer to the local node n (const version) More...
 
unsigned nnode () const
 Return the number of nodes. More...
 
virtual unsigned nnode_1d () const
 
double raw_nodal_position (const unsigned &n, const unsigned &i) const
 
double raw_nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &i) const
 
double raw_dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double raw_dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position (const unsigned &n, const unsigned &i) const
 
double nodal_position (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double dnodal_position_dt (const unsigned &n, const unsigned &i) const
 Return the i-th component of nodal velocity: dx/dt at local node n. More...
 
double dnodal_position_dt (const unsigned &n, const unsigned &j, const unsigned &i) const
 
double nodal_position_gen (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double nodal_position_gen (const unsigned &t, const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &n, const unsigned &k, const unsigned &i) const
 
double dnodal_position_gen_dt (const unsigned &j, const unsigned &n, const unsigned &k, const unsigned &i) const
 
virtual void get_dresidual_dnodal_coordinates (RankThreeTensor< double > &dresidual_dnodal_coordinates)
 
virtual void disable_ALE ()
 
virtual void enable_ALE ()
 
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 Node * vertex_node_pt (const unsigned &j) const
 
virtual Node * construct_node (const unsigned &n)
 
virtual Node * construct_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
virtual Node * construct_boundary_node (const unsigned &n)
 
virtual Node * construct_boundary_node (const unsigned &n, TimeStepper *const &time_stepper_pt)
 
int get_node_number (Node *const &node_pt) const
 
virtual Node * get_node_at_local_coordinate (const Vector< double > &s) const
 
double raw_nodal_value (const unsigned &n, const unsigned &i) const
 
double raw_nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &n, const unsigned &i) const
 
double nodal_value (const unsigned &t, const unsigned &n, const unsigned &i) const
 
unsigned dim () const
 
virtual ElementGeometry::ElementGeometry element_geometry () const
 Return the geometry type of the element (either Q or T usually). More...
 
virtual double interpolated_x (const Vector< double > &s, const unsigned &i) const
 Return FE interpolated coordinate x[i] at local coordinate s. More...
 
virtual double interpolated_x (const unsigned &t, const Vector< double > &s, const unsigned &i) const
 
virtual void interpolated_x (const Vector< double > &s, Vector< double > &x) const
 Return FE interpolated position x[] at local coordinate s as Vector. More...
 
virtual void interpolated_x (const unsigned &t, const Vector< double > &s, Vector< double > &x) const
 
virtual double interpolated_dxdt (const Vector< double > &s, const unsigned &i, const unsigned &t)
 
virtual void interpolated_dxdt (const Vector< double > &s, const unsigned &t, Vector< double > &dxdt)
 
unsigned ngeom_data () const
 
Data * geom_data_pt (const unsigned &j)
 
void position (const Vector< double > &zeta, Vector< double > &r) const
 
void position (const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
 
void dposition_dt (const Vector< double > &zeta, const unsigned &t, Vector< double > &drdt)
 
virtual double zeta_nodal (const unsigned &n, const unsigned &k, const unsigned &i) const
 
void interpolated_zeta (const Vector< double > &s, Vector< double > &zeta) const
 
void locate_zeta (const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
 
virtual void node_update ()
 
virtual void identify_field_data_for_interactions (std::set< std::pair< Data *, unsigned >> &paired_field_data)
 
virtual void identify_geometric_data (std::set< Data * > &geometric_data_pt)
 
virtual double s_min () const
 Min value of local coordinate. More...
 
virtual double s_max () const
 Max. value of local coordinate. More...
 
double size () const
 
virtual double compute_physical_size () const
 
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 (FILE *file_pt)
 
virtual void output (FILE *file_pt, const unsigned &n_plot)
 
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_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 get_bulk_node_number (const int &face_index, const unsigned &i) const
 
virtual int face_outer_unit_normal_sign (const int &face_index) const
 Get the sign of the outer unit normal on the face given by face_index. More...
 
virtual unsigned nnode_on_face () const
 
void face_node_number_error_check (const unsigned &i) const
 Range check for face node numbers. More...
 
virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt (const int &face_index) const
 
virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt (const int &face_index) const
 
- Public Member Functions inherited from oomph::GeneralisedElement
 GeneralisedElement ()
 Constructor: Initialise all pointers and all values to zero. More...
 
virtual ~GeneralisedElement ()
 Virtual destructor to clean up any memory allocated by the object. More...
 
 GeneralisedElement (const GeneralisedElement &)=delete
 Broken copy constructor. More...
 
void operator= (const GeneralisedElement &)=delete
 Broken assignment operator. More...
 
Data *& internal_data_pt (const unsigned &i)
 Return a pointer to i-th internal data object. More...
 
Data *const & internal_data_pt (const unsigned &i) const
 Return a pointer to i-th internal data object (const version) More...
 
Data *& external_data_pt (const unsigned &i)
 Return a pointer to i-th external data object. More...
 
Data *const & external_data_pt (const unsigned &i) const
 Return a pointer to i-th external data object (const version) More...
 
unsigned long eqn_number (const unsigned &ieqn_local) const
 
int local_eqn_number (const unsigned long &ieqn_global) const
 
unsigned add_external_data (Data *const &data_pt, const bool &fd=true)
 
bool external_data_fd (const unsigned &i) const
 
void exclude_external_data_fd (const unsigned &i)
 
void include_external_data_fd (const unsigned &i)
 
void flush_external_data ()
 Flush all external data. More...
 
void flush_external_data (Data *const &data_pt)
 Flush the object addressed by data_pt from the external data array. More...
 
unsigned ninternal_data () const
 Return the number of internal data objects. More...
 
unsigned nexternal_data () const
 Return the number of external data objects. More...
 
unsigned ndof () const
 Return the number of equations/dofs in the element. More...
 
void dof_vector (const unsigned &t, Vector< double > &dof)
 Return the vector of dof values at time level t. More...
 
void dof_pt_vector (Vector< double * > &dof_pt)
 Return the vector of pointers to dof values. More...
 
void set_internal_data_time_stepper (const unsigned &i, TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void assign_internal_eqn_numbers (unsigned long &global_number, Vector< double * > &Dof_pt)
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void add_internal_value_pt_to_map (std::map< unsigned, double * > &map_of_value_pt)
 
virtual void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 
virtual void complete_setup_of_dependencies ()
 
virtual void get_residuals (Vector< double > &residuals)
 
virtual void get_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian)
 
virtual void get_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
 
virtual void get_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
 
virtual void get_dresiduals_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam)
 
virtual void get_djacobian_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
 
virtual void get_djacobian_and_dmass_matrix_dparameter (double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
 
virtual void get_hessian_vector_products (Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
 
virtual void get_inner_products (Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
 
virtual void get_inner_product_vectors (Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_norm (double &norm)
 
virtual unsigned ndof_types () const
 
virtual void get_dof_numbers_for_unknowns (std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
 
- Public Member Functions inherited from oomph::GeomObject
 GeomObject ()
 Default constructor. More...
 
 GeomObject (const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim)
 
 GeomObject (const unsigned &nlagrangian, const unsigned &ndim, TimeStepper *time_stepper_pt)
 
 GeomObject (const GeomObject &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const GeomObject &)=delete
 Broken assignment operator. More...
 
virtual ~GeomObject ()
 (Empty) destructor More...
 
unsigned nlagrangian () const
 Access function to # of Lagrangian coordinates. More...
 
unsigned ndim () const
 Access function to # of Eulerian coordinates. More...
 
void set_nlagrangian_and_ndim (const unsigned &n_lagrangian, const unsigned &n_dim)
 Set # of Lagrangian and Eulerian coordinates. More...
 
TimeStepper *& time_stepper_pt ()
 
TimeStepper * time_stepper_pt () const
 
virtual void position (const double &t, const Vector< double > &zeta, Vector< double > &r) const
 
virtual void dposition (const Vector< double > &zeta, DenseMatrix< double > &drdzeta) const
 
virtual void d2position (const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
 
virtual void d2position (const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
 

Protected Member Functions

virtual double shape_basis_test_local (const Vector< double > &s, Shape &psi, DShape &dpsi, Shape &u_basis, Shape &u_test, DShape &du_basis_dx, DShape &du_test_dx, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const =0
 
virtual double shape_basis_test_local_at_knot (const unsigned &ipt, Shape &psi, DShape &dpsi, Shape &u_basis, Shape &u_test, DShape &du_basis_dx, DShape &du_test_dx, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const =0
 
virtual void fill_in_generic_residual_contribution (Vector< double > &residuals, DenseMatrix< double > &jacobian, bool flag)
 Fill in residuals and, if flag==true, jacobian. More...
 
- 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_jacobian_and_mass_matrix (Vector< double > &residuals, DenseMatrix< double > &jacobian, 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

ElasticityTensor * Elasticity_tensor_pt
 Pointer to the elasticity tensor. More...
 
- Protected Attributes inherited from oomph::FiniteElement
MacroElement * Macro_elem_pt
 Pointer to the element's macro element (NULL by default) More...
 
- Protected Attributes inherited from oomph::GeomObject
unsigned NLagrangian
 Number of Lagrangian (intrinsic) coordinates. More...
 
unsigned Ndim
 Number of Eulerian coordinates. More...
 
TimeStepper * Geom_object_time_stepper_pt
 

Private Attributes

SourceFctPt Force_solid_fct_pt
 Pointer to solid source function. More...
 
SourceFctPt Force_fluid_fct_pt
 Pointer to fluid source function. More...
 
MassSourceFctPt Mass_source_fct_pt
 Pointer to the mass source function. More...
 
double * Lambda_sq_pt
 Timescale ratio (non-dim. density) More...
 
double * Density_ratio_pt
 Density ratio. More...
 
double * K_inv_pt
 1/k More...
 
double * Alpha_pt
 Alpha. More...
 
double * Porosity_pt
 Porosity. More...
 

Static Private Attributes

static double Default_lambda_sq_value = 1.0
 Static default value for timescale ratio (1.0 – for natural scaling) More...
 
static double Default_density_ratio_value = 1.0
 Static default value for the density ratio. More...
 
static double Default_k_inv_value = 1.0
 Static default value for 1/k. More...
 
static double Default_alpha_value = 1.0
 Static default value for alpha. More...
 
static double Default_porosity_value = 1.0
 Static default value for the porosity. More...
 

Additional Inherited Members

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

Detailed Description

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

Class implementing the generic maths of the poroelasticity equations: linear elasticity coupled with Darcy equations (using Raviart-Thomas elements with both edge and internal degrees of freedom)

Member Typedef Documentation

◆ MassSourceFctPt

template<unsigned DIM>
typedef void(* oomph::PoroelasticityEquations< DIM >::MassSourceFctPt) (const double &time, const Vector< double > &x, double &f)

Mass source function pointer typedef.

◆ SourceFctPt

template<unsigned DIM>
typedef void(* oomph::PoroelasticityEquations< DIM >::SourceFctPt) (const double &time, const Vector< double > &x, Vector< double > &f)

Source function pointer typedef.

Constructor & Destructor Documentation

◆ PoroelasticityEquations()

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

Constructor.

60  : Elasticity_tensor_pt(0),
61  Force_solid_fct_pt(0),
62  Force_fluid_fct_pt(0),
63  Mass_source_fct_pt(0),
64  Lambda_sq_pt(&Default_lambda_sq_value),
65  Density_ratio_pt(&Default_density_ratio_value),
66  K_inv_pt(&Default_k_inv_value),
67  Alpha_pt(&Default_alpha_value),
68  Porosity_pt(&Default_porosity_value)
69  {
70  }
oomph::PoroelasticityEquations::Force_fluid_fct_pt
SourceFctPt Force_fluid_fct_pt
Pointer to fluid source function.
Definition: poroelasticity_elements.h:790
oomph::PoroelasticityEquations::Mass_source_fct_pt
MassSourceFctPt Mass_source_fct_pt
Pointer to the mass source function.
Definition: poroelasticity_elements.h:793
oomph::PoroelasticityEquations::Porosity_pt
double * Porosity_pt
Porosity.
Definition: poroelasticity_elements.h:808
oomph::PoroelasticityEquations::Default_porosity_value
static double Default_porosity_value
Static default value for the porosity.
Definition: poroelasticity_elements.h:823
oomph::PoroelasticityEquations::Force_solid_fct_pt
SourceFctPt Force_solid_fct_pt
Pointer to solid source function.
Definition: poroelasticity_elements.h:787
oomph::PoroelasticityEquations::Default_density_ratio_value
static double Default_density_ratio_value
Static default value for the density ratio.
Definition: poroelasticity_elements.h:814
oomph::PoroelasticityEquations::Alpha_pt
double * Alpha_pt
Alpha.
Definition: poroelasticity_elements.h:805
oomph::PoroelasticityEquations::Default_lambda_sq_value
static double Default_lambda_sq_value
Static default value for timescale ratio (1.0 – for natural scaling)
Definition: poroelasticity_elements.h:811
oomph::PoroelasticityEquations::Default_alpha_value
static double Default_alpha_value
Static default value for alpha.
Definition: poroelasticity_elements.h:820
oomph::PoroelasticityEquations::Elasticity_tensor_pt
ElasticityTensor * Elasticity_tensor_pt
Pointer to the elasticity tensor.
Definition: poroelasticity_elements.h:783
oomph::PoroelasticityEquations::Density_ratio_pt
double * Density_ratio_pt
Density ratio.
Definition: poroelasticity_elements.h:799
oomph::PoroelasticityEquations::Default_k_inv_value
static double Default_k_inv_value
Static default value for 1/k.
Definition: poroelasticity_elements.h:817
oomph::PoroelasticityEquations::K_inv_pt
double * K_inv_pt
1/k
Definition: poroelasticity_elements.h:802
oomph::PoroelasticityEquations::Lambda_sq_pt
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
Definition: poroelasticity_elements.h:796

Member Function Documentation

◆ alpha()

template<unsigned DIM>
const double& oomph::PoroelasticityEquations< DIM >::alpha ( ) const
inline

Access function for alpha.

110  {
111  return *Alpha_pt;
112  }

References oomph::PoroelasticityEquations< DIM >::Alpha_pt.

◆ alpha_pt()

template<unsigned DIM>
double*& oomph::PoroelasticityEquations< DIM >::alpha_pt ( )
inline

Access function for pointer to alpha.

116  {
117  return Alpha_pt;
118  }

References oomph::PoroelasticityEquations< DIM >::Alpha_pt.

◆ compute_error() [1/2]

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

Compute the error between the FE solution and the exact solution using the H(div) norm for q and L^2 norm for p

Compute the error between the FE solution and the exact solution using the H(div) norm for q and L^2 norm for u and p

Reimplemented from oomph::FiniteElement.

467  {
468  for (unsigned i = 0; i < 3; i++)
469  {
470  error[i] = 0.0;
471  norm[i] = 0.0;
472  }
473 
474  // Vector of local coordinates
475  Vector<double> s(DIM);
476 
477  // Vector for coordinates
478  Vector<double> x(DIM);
479 
480  // Set the value of n_intpt
481  unsigned n_intpt = this->integral_pt()->nweight();
482 
483  outfile << "ZONE" << std::endl;
484 
485  // Exact solution Vector (u,v,[w])
486  Vector<double> exact_soln(2 * DIM + 2);
487 
488  // Loop over the integration points
489  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
490  {
491  // Assign values of s
492  for (unsigned i = 0; i < DIM; i++)
493  {
494  s[i] = this->integral_pt()->knot(ipt, i);
495  }
496 
497  // Get the integral weight
498  double w = this->integral_pt()->weight(ipt);
499 
500  // Get jacobian of mapping
501  double J = this->J_eulerian(s);
502 
503  // Premultiply the weights and the Jacobian
504  double W = w * J;
505 
506  // Get x position as Vector
507  this->interpolated_x(s, x);
508 
509  // Get exact solution at this point
510  (*exact_soln_pt)(x, exact_soln);
511 
512  // Displacement error
513  for (unsigned i = 0; i < DIM; i++)
514  {
515  norm[0] += exact_soln[i] * exact_soln[i] * W;
516  // Error due to q_i
517  error[0] += (exact_soln[i] - this->interpolated_u(s, i)) *
518  (exact_soln[i] - this->interpolated_u(s, i)) * W;
519  }
520 
521  // Flux error
522  for (unsigned i = 0; i < DIM; i++)
523  {
524  norm[1] += exact_soln[DIM + i] * exact_soln[DIM + i] * W;
525  // Error due to q_i
526  error[1] += (exact_soln[DIM + i] - this->interpolated_q(s, i)) *
527  (exact_soln[DIM + i] - this->interpolated_q(s, i)) * W;
528  }
529 
530  // Flux divergence error
531  norm[1] += exact_soln[2 * DIM] * exact_soln[2 * DIM] * W;
532  error[1] += (exact_soln[2 * DIM] - interpolated_div_q(s)) *
533  (exact_soln[2 * DIM] - interpolated_div_q(s)) * W;
534 
535  // Pressure error
536  norm[2] += exact_soln[2 * DIM + 1] * exact_soln[2 * DIM + 1] * W;
537  error[2] += (exact_soln[2 * DIM + 1] - this->interpolated_p(s)) *
538  (exact_soln[2 * DIM + 1] - this->interpolated_p(s)) * W;
539 
540  // Output x,y,[z]
541  for (unsigned i = 0; i < DIM; i++)
542  {
543  outfile << x[i] << " ";
544  }
545 
546  // Output u_1_error,u_2_error,...
547  for (unsigned i = 0; i < DIM; i++)
548  {
549  outfile << exact_soln[i] - this->interpolated_u(s, i) << " ";
550  }
551 
552  // Output q_1_error,q_2_error,...
553  for (unsigned i = 0; i < DIM; i++)
554  {
555  outfile << exact_soln[DIM + i] - this->interpolated_q(s, i) << " ";
556  }
557 
558  // Output p_error
559  outfile << exact_soln[2 * DIM + 1] - this->interpolated_p(s) << " ";
560 
561  outfile << std::endl;
562  }
563  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
J
JacobiRotation< float > J
Definition: Jacobi_makeJacobi.cpp:3
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
oomph::FiniteElement::J_eulerian
virtual double J_eulerian(const Vector< double > &s) const
Definition: elements.cc:4103
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::PoroelasticityEquations::interpolated_p
void interpolated_p(const Vector< double > &s, double &p) const
Calculate the FE representation of p.
Definition: poroelasticity_elements.h:540
oomph::PoroelasticityEquations::interpolated_q
void interpolated_q(const Vector< double > &s, Vector< double > &u) const
Calculate the FE representation of q.
Definition: poroelasticity_elements.h:437
oomph::PoroelasticityEquations::interpolated_div_q
void interpolated_div_q(const Vector< double > &s, double &div_q) const
Calculate the FE representation of div u.
Definition: poroelasticity_elements.h:482
oomph::PoroelasticityEquations::interpolated_u
void interpolated_u(const Vector< double > &s, Vector< double > &disp) const
Calculate the FE representation of u.
Definition: poroelasticity_elements.h:382
s
RealScalar s
Definition: level1_cplx_impl.h:130
DIM
#define DIM
Definition: linearised_navier_stokes_elements.h:44
ProblemParameters::exact_soln
void exact_soln(const double &time, const Vector< double > &x, Vector< double > &soln)
Definition: unstructured_two_d_curved.cc:301
calibrate.error
int error
Definition: calibrate.py:297
oomph::QuadTreeNames::W
@ W
Definition: quadtree.h:63
plotDoE.x
list x
Definition: plotDoE.py:28

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

◆ compute_error() [2/2]

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

Compute the error between the FE solution and the exact solution using the H(div) norm for q and L^2 norm for p. Unsteady version

Compute the error between the FE solution and the exact solution using the H(div) norm for u and L^2 norm for p. Unsteady version

Reimplemented from oomph::FiniteElement.

574  {
575  for (unsigned i = 0; i < 3; i++)
576  {
577  error[i] = 0.0;
578  norm[i] = 0.0;
579  }
580 
581  // Vector of local coordinates
582  Vector<double> s(DIM);
583 
584  // Vector for coordinates
585  Vector<double> x(DIM);
586 
587  // Set the value of n_intpt
588  unsigned n_intpt = this->integral_pt()->nweight();
589 
590  outfile << "ZONE" << std::endl;
591 
592  // Exact solution Vector (u,v,[w])
593  Vector<double> exact_soln(2 * DIM + 2);
594 
595  // Loop over the integration points
596  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
597  {
598  // Assign values of s
599  for (unsigned i = 0; i < DIM; i++)
600  {
601  s[i] = this->integral_pt()->knot(ipt, i);
602  }
603 
604  // Get the integral weight
605  double w = this->integral_pt()->weight(ipt);
606 
607  // Get jacobian of mapping
608  double J = this->J_eulerian(s);
609 
610  // Premultiply the weights and the Jacobian
611  double W = w * J;
612 
613  // Get x position as Vector
614  this->interpolated_x(s, x);
615 
616  // Get exact solution at this point
617  (*exact_soln_pt)(time, x, exact_soln);
618 
619  // Displacement error
620  for (unsigned i = 0; i < DIM; i++)
621  {
622  norm[0] += exact_soln[i] * exact_soln[i] * W;
623  // Error due to q_i
624  error[0] += (exact_soln[i] - this->interpolated_u(s, i)) *
625  (exact_soln[i] - this->interpolated_u(s, i)) * W;
626  }
627 
628  // Flux error
629  for (unsigned i = 0; i < DIM; i++)
630  {
631  norm[1] += exact_soln[DIM + i] * exact_soln[DIM + i] * W;
632  // Error due to q_i
633  error[1] += (exact_soln[DIM + i] - this->interpolated_q(s, i)) *
634  (exact_soln[DIM + i] - this->interpolated_q(s, i)) * W;
635  }
636 
637  // Flux divergence error
638  norm[1] += exact_soln[2 * DIM] * exact_soln[2 * DIM] * W;
639  error[1] += (exact_soln[2 * DIM] - interpolated_div_q(s)) *
640  (exact_soln[2 * DIM] - interpolated_div_q(s)) * W;
641 
642  // Pressure error
643  norm[2] += exact_soln[2 * DIM + 1] * exact_soln[2 * DIM + 1] * W;
644  error[2] += (exact_soln[2 * DIM + 1] - this->interpolated_p(s)) *
645  (exact_soln[2 * DIM + 1] - this->interpolated_p(s)) * W;
646 
647  // Output x,y,[z]
648  for (unsigned i = 0; i < DIM; i++)
649  {
650  outfile << x[i] << " ";
651  }
652 
653  // Output u_1_error,u_2_error,...
654  for (unsigned i = 0; i < DIM; i++)
655  {
656  outfile << exact_soln[i] - this->interpolated_u(s, i) << " ";
657  }
658 
659  // Output q_1_error,q_2_error,...
660  for (unsigned i = 0; i < DIM; i++)
661  {
662  outfile << exact_soln[DIM + i] - this->interpolated_q(s, i) << " ";
663  }
664 
665  // Output p_error
666  outfile << exact_soln[2 * DIM + 1] - this->interpolated_p(s) << " ";
667 
668  outfile << std::endl;
669  }
670  }

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

◆ d2u_dt2()

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::d2u_dt2 ( const unsigned &  n,
const unsigned &  i 
) const
inline

d^2u/dt^2 at local node n

606  {
607  // Get the timestepper
608  TimeStepper* time_stepper_pt = node_pt(n)->time_stepper_pt();
609 
610  // Storage for the derivative - initialise to 0
611  double d2u_dt2 = 0.0;
612 
613  // If we are doing an unsteady solve then calculate the derivative
614  if (!time_stepper_pt->is_steady())
615  {
616  // Get the nodal index
617  const unsigned u_nodal_index = u_index(i);
618 
619  // Get the number of values required to represent history
620  const unsigned n_time = time_stepper_pt->ntstorage();
621 
622  // Loop over history values
623  for (unsigned t = 0; t < n_time; t++)
624  {
625  // Add the contribution to the derivative
626  d2u_dt2 +=
627  time_stepper_pt->weight(2, t) * nodal_value(t, n, u_nodal_index);
628  }
629  }
630 
631  return d2u_dt2;
632  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::nodal_value
double nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2593
oomph::GeomObject::time_stepper_pt
TimeStepper *& time_stepper_pt()
Definition: geom_objects.h:192
oomph::PoroelasticityEquations::u_index
virtual unsigned u_index(const unsigned &n) const =0
Return the nodal index of the n-th solid displacement unknown.
oomph::PoroelasticityEquations::d2u_dt2
double d2u_dt2(const unsigned &n, const unsigned &i) const
d^2u/dt^2 at local node n
Definition: poroelasticity_elements.h:605
oomph::TimeStepper::ntstorage
unsigned ntstorage() const
Definition: timesteppers.h:601
oomph::TimeStepper::weight
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
oomph::TimeStepper::is_steady
bool is_steady() const
Definition: timesteppers.h:389
plotPSD.t
t
Definition: plotPSD.py:36

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

◆ density_ratio()

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

Access function for the density ratio.

86  {
87  return *Density_ratio_pt;
88  }

References oomph::PoroelasticityEquations< DIM >::Density_ratio_pt.

◆ density_ratio_pt()

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

Access function for pointer to the density ratio.

92  {
93  return Density_ratio_pt;
94  }

References oomph::PoroelasticityEquations< DIM >::Density_ratio_pt.

◆ dq_edge_dt()

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::dq_edge_dt ( const unsigned &  n) const
inline

dq_edge/dt for the n-th edge degree of freedom

636  {
637  unsigned node_num = q_edge_node_number(n);
638 
639  // get the timestepper
640  TimeStepper* time_stepper_pt = node_pt(node_num)->time_stepper_pt();
641 
642  // storage for the derivative - initialise to 0
643  double dq_dt = 0.0;
644 
645  // if we are doing an unsteady solve then calculate the derivative
646  if (!time_stepper_pt->is_steady())
647  {
648  // get the number of values required to represent history
649  const unsigned n_time = time_stepper_pt->ntstorage();
650 
651  // loop over history values
652  for (unsigned t = 0; t < n_time; t++)
653  {
654  // add the contribution to the derivative
655  dq_dt += time_stepper_pt->weight(1, t) * q_edge(t, n);
656  }
657  }
658 
659  return dq_dt;
660  }
oomph::PoroelasticityEquations::q_edge_node_number
virtual unsigned q_edge_node_number(const unsigned &n) const =0
Return the number of the node where the nth edge unknown is stored.
oomph::PoroelasticityEquations::q_edge
virtual double q_edge(const unsigned &n) const =0
Return the values of the edge (flux) degrees of freedom.

References oomph::TimeStepper::is_steady(), n, oomph::FiniteElement::node_pt(), oomph::TimeStepper::ntstorage(), oomph::PoroelasticityEquations< DIM >::q_edge(), oomph::PoroelasticityEquations< DIM >::q_edge_node_number(), plotPSD::t, oomph::GeomObject::time_stepper_pt(), oomph::Data::time_stepper_pt(), and oomph::TimeStepper::weight().

◆ dq_internal_dt()

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::dq_internal_dt ( const unsigned &  n) const
inline

dq_internal/dt for the n-th internal degree of freedom

664  {
665  // get the internal data index for q
666  unsigned internal_index = q_internal_index();
667 
668  // get the timestepper
669  TimeStepper* time_stepper_pt =
670  internal_data_pt(internal_index)->time_stepper_pt();
671 
672  // storage for the derivative - initialise to 0
673  double dq_dt = 0.0;
674 
675  // if we are doing an unsteady solve then calculate the derivative
676  if (!time_stepper_pt->is_steady())
677  {
678  // get the number of values required to represent history
679  const unsigned n_time = time_stepper_pt->ntstorage();
680 
681  // loop over history values
682  for (unsigned t = 0; t < n_time; t++)
683  {
684  // add the contribution to the derivative
685  dq_dt += time_stepper_pt->weight(1, t) * q_internal(t, n);
686  }
687  }
688 
689  return dq_dt;
690  }
oomph::GeneralisedElement::internal_data_pt
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
oomph::PoroelasticityEquations::q_internal_index
virtual unsigned q_internal_index() const =0
oomph::PoroelasticityEquations::q_internal
virtual double q_internal(const unsigned &n) const =0
Return the values of the internal (moment) degrees of freedom.

References oomph::GeneralisedElement::internal_data_pt(), oomph::TimeStepper::is_steady(), n, oomph::TimeStepper::ntstorage(), oomph::PoroelasticityEquations< DIM >::q_internal(), oomph::PoroelasticityEquations< DIM >::q_internal_index(), plotPSD::t, oomph::GeomObject::time_stepper_pt(), oomph::Data::time_stepper_pt(), and oomph::TimeStepper::weight().

◆ du_dt()

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::du_dt ( const unsigned &  n,
const unsigned &  i 
) const
inline

du/dt at local node n

576  {
577  // Get the timestepper
578  TimeStepper* time_stepper_pt = node_pt(n)->time_stepper_pt();
579 
580  // Storage for the derivative - initialise to 0
581  double du_dt = 0.0;
582 
583  // If we are doing an unsteady solve then calculate the derivative
584  if (!time_stepper_pt->is_steady())
585  {
586  // Get the nodal index
587  const unsigned u_nodal_index = u_index(i);
588 
589  // Get the number of values required to represent history
590  const unsigned n_time = time_stepper_pt->ntstorage();
591 
592  // Loop over history values
593  for (unsigned t = 0; t < n_time; t++)
594  {
595  // Add the contribution to the derivative
596  du_dt +=
597  time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
598  }
599  }
600 
601  return du_dt;
602  }
oomph::PoroelasticityEquations::du_dt
double du_dt(const unsigned &n, const unsigned &i) const
du/dt at local node n
Definition: poroelasticity_elements.h:575

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

◆ E()

template<unsigned DIM>
const double oomph::PoroelasticityEquations< DIM >::E ( const unsigned &  i,
const unsigned &  j,
const unsigned &  k,
const unsigned &  l 
) const
inline

Access function to the entries in the elasticity tensor.

240  {
241  return (*Elasticity_tensor_pt)(i, j, k, l);
242  }
k
char char char int int * k
Definition: level2_impl.h:374
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::PoroelasticityEquations< DIM >::Elasticity_tensor_pt, i, j, and k.

◆ edge_gauss_point()

template<unsigned DIM>
virtual double oomph::PoroelasticityEquations< DIM >::edge_gauss_point ( const unsigned &  edge,
const unsigned &  n 
) const
pure virtual

Returns the nth gauss point along an edge.

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ edge_gauss_point_global()

template<unsigned DIM>
virtual void oomph::PoroelasticityEquations< DIM >::edge_gauss_point_global ( const unsigned &  edge,
const unsigned &  n,
Vector< double > &  x 
) const
pure virtual

Returns the global coordinates of the nth gauss point along an edge.

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ elasticity_tensor_pt()

template<unsigned DIM>
ElasticityTensor*& oomph::PoroelasticityEquations< DIM >::elasticity_tensor_pt ( )
inline

Return the pointer to the elasticity_tensor.

231  {
232  return Elasticity_tensor_pt;
233  }

References oomph::PoroelasticityEquations< DIM >::Elasticity_tensor_pt.

◆ fill_in_contribution_to_jacobian()

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

Fill in the Jacobian matrix for the Newton method.

Reimplemented from oomph::FiniteElement.

377  {
378  this->fill_in_generic_residual_contribution(residuals, jacobian, 1);
379  }
oomph::PoroelasticityEquations::fill_in_generic_residual_contribution
virtual void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, bool flag)
Fill in residuals and, if flag==true, jacobian.
Definition: poroelasticity_elements.cc:674

References oomph::PoroelasticityEquations< DIM >::fill_in_generic_residual_contribution().

◆ fill_in_contribution_to_residuals()

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

Fill in contribution to residuals for the Darcy equations.

Reimplemented from oomph::GeneralisedElement.

369  {
370  this->fill_in_generic_residual_contribution(
371  residuals, GeneralisedElement::Dummy_matrix, 0);
372  }
oomph::GeneralisedElement::Dummy_matrix
static DenseMatrix< double > Dummy_matrix
Definition: elements.h:227

References oomph::GeneralisedElement::Dummy_matrix, and oomph::PoroelasticityEquations< DIM >::fill_in_generic_residual_contribution().

◆ fill_in_generic_residual_contribution()

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::fill_in_generic_residual_contribution ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
bool  flag 
)
protectedvirtual

Fill in residuals and, if flag==true, jacobian.

676  {
677  // Get the number of geometric nodes, total number of basis functions,
678  // and number of edges basis functions
679  const unsigned n_node = nnode();
680  const unsigned n_q_basis = nq_basis();
681  const unsigned n_q_basis_edge = nq_basis_edge();
682  const unsigned n_p_basis = np_basis();
683 
684  // Storage for the geometric and computational bases and the test functions
685  Shape psi(n_node), u_basis(n_node), u_test(n_node), q_basis(n_q_basis, DIM),
686  q_test(n_q_basis, DIM), p_basis(n_p_basis), p_test(n_p_basis),
687  div_q_basis_ds(n_q_basis), div_q_test_ds(n_q_basis);
688 
689  DShape dpsidx(n_node, DIM), du_basis_dx(n_node, DIM),
690  du_test_dx(n_node, DIM);
691 
692  // Get the number of integration points
693  unsigned n_intpt = integral_pt()->nweight();
694 
695  // Storage for the local coordinates
696  Vector<double> s(DIM);
697 
698  // Storage for the elasticity source function
699  Vector<double> f_solid(DIM);
700 
701  // Storage for the source function
702  Vector<double> f_fluid(DIM);
703 
704  // Storage for the mass source function
705  double mass_source_local;
706 
707  // Storage for Lambda_sq
708  double lambda_sq = this->lambda_sq();
709 
710  // Get the value of 1/k
711  double k_inv = this->k_inv();
712 
713  // Get the value of alpha
714  double alpha = this->alpha();
715 
716  // Get the value of the porosity
717  double porosity = this->porosity();
718 
719  // Get the density ratio
720  double density_ratio = this->density_ratio();
721 
722  // Precompute the ratio of fluid density to combined density
723  double rho_f_over_rho =
724  density_ratio / (1 + porosity * (density_ratio - 1));
725 
726  // Get continuous time from timestepper of first node
727  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
728 
729  // Local equation and unknown numbers
730  int local_eqn = 0, local_unknown = 0;
731 
732  // Loop over the integration points
733  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
734  {
735  // Find the local coordinates at the integration point
736  for (unsigned i = 0; i < DIM; i++)
737  {
738  s[i] = integral_pt()->knot(ipt, i);
739  }
740 
741  // Get the weight of the intetgration point
742  double w = integral_pt()->weight(ipt);
743 
744  // Call the basis functions and test functions and get the
745  // (geometric) jacobian of the current element
746  double J = shape_basis_test_local_at_knot(ipt,
747  psi,
748  dpsidx,
749  u_basis,
750  u_test,
751  du_basis_dx,
752  du_test_dx,
753  q_basis,
754  q_test,
755  p_basis,
756  p_test,
757  div_q_basis_ds,
758  div_q_test_ds);
759 
760  // Storage for interpolated values
761  Vector<double> interpolated_x(DIM, 0.0);
762  Vector<double> interpolated_u(DIM, 0.0);
763  DenseMatrix<double> interpolated_du_dx(DIM, DIM, 0.0);
764  double interpolated_div_du_dt_dx = 0.0;
765  Vector<double> interpolated_d2u_dt2(DIM, 0.0);
766  Vector<double> interpolated_q(DIM, 0.0);
767  double interpolated_div_q_ds = 0.0;
768  Vector<double> interpolated_dq_dt(DIM, 0.0);
769  double interpolated_p = 0.0;
770 
771  // loop over the vector basis functions to find interpolated x
772  for (unsigned l = 0; l < n_node; l++)
773  {
774  // Loop over the geometric basis functions
775  for (unsigned i = 0; i < DIM; i++)
776  {
777  interpolated_x[i] += nodal_position(l, i) * psi(l);
778 
779  interpolated_d2u_dt2[i] += this->d2u_dt2(l, i) * u_basis(l);
780 
781  // Get the nodal displacements
782  const double u_value = this->raw_nodal_value(l, u_index(i));
783  interpolated_u[i] += u_value * u_basis(l);
784 
785  // Loop over derivative directions
786  for (unsigned j = 0; j < DIM; j++)
787  {
788  interpolated_du_dx(i, j) += u_value * du_basis_dx(l, j);
789  }
790 
791  // divergence of the time derivative of the solid displacement
792  interpolated_div_du_dt_dx += this->du_dt(l, i) * du_basis_dx(l, i);
793  }
794  }
795 
796  // loop over the nodes and use the geometric basis functions to find the
797  // interpolated flux and its time derivative
798  for (unsigned i = 0; i < DIM; i++)
799  {
800  // Loop over the edge basis vectors
801  for (unsigned l = 0; l < n_q_basis_edge; l++)
802  {
803  interpolated_q[i] += q_edge(l) * q_basis(l, i);
804  interpolated_dq_dt[i] += dq_edge_dt(l) * q_basis(l, i);
805  }
806  // Loop over the internal basis vectors
807  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
808  {
809  interpolated_q[i] += q_internal(l - n_q_basis_edge) * q_basis(l, i);
810  interpolated_dq_dt[i] +=
811  dq_internal_dt(l - n_q_basis_edge) * q_basis(l, i);
812  }
813  }
814 
815  // loop over the pressure basis and find the interpolated pressure
816  for (unsigned l = 0; l < n_p_basis; l++)
817  {
818  interpolated_p += p_value(l) * p_basis(l);
819  }
820 
821  // loop over the u edge divergence basis and the u internal divergence
822  // basis to find interpolated div u
823  for (unsigned l = 0; l < n_q_basis_edge; l++)
824  {
825  interpolated_div_q_ds += q_edge(l) * div_q_basis_ds(l);
826  }
827  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
828  {
829  interpolated_div_q_ds +=
830  q_internal(l - n_q_basis_edge) * div_q_basis_ds(l);
831  }
832 
833  // Get the solid forcing
834  this->force_solid(time, interpolated_x, f_solid);
835 
836  // Get the fluid forcing
837  this->force_fluid(time, interpolated_x, f_fluid);
838 
839  // Get the mass source function
840  this->mass_source(time, interpolated_x, mass_source_local);
841 
842  // Premultiply the weights and the Jacobian
843  double W = w * J;
844 
845  // Linear elasticity:
846 
847  for (unsigned l = 0; l < n_node; l++)
848  {
849  for (unsigned a = 0; a < DIM; a++)
850  {
851  local_eqn = this->nodal_local_eqn(l, u_index(a));
852 
853  if (local_eqn >= 0)
854  {
855  residuals[local_eqn] +=
856  (lambda_sq * (interpolated_d2u_dt2[a] +
857  rho_f_over_rho * interpolated_dq_dt[a]) -
858  f_solid[a]) *
859  u_test(l) * W;
860 
861  // Stress term
862  for (unsigned b = 0; b < DIM; b++)
863  {
864  if (a == b)
865  {
866  residuals[local_eqn] -=
867  alpha * interpolated_p * du_test_dx(l, b) * W;
868  }
869 
870  for (unsigned c = 0; c < DIM; c++)
871  {
872  for (unsigned d = 0; d < DIM; d++)
873  {
874  // Add the stress terms to the residuals
875  residuals[local_eqn] += this->E(a, b, c, d) *
876  interpolated_du_dx(c, d) *
877  du_test_dx(l, b) * W;
878  }
879  }
880  }
881  // Jacobian entries
882  if (flag)
883  {
884  // d(u_eqn_l,a)/d(U_l2,c)
885  for (unsigned l2 = 0; l2 < n_node; l2++)
886  {
887  for (unsigned c = 0; c < DIM; c++)
888  {
889  local_unknown = this->nodal_local_eqn(l2, u_index(c));
890  if (local_unknown >= 0)
891  {
892  if (a == c)
893  {
894  jacobian(local_eqn, local_unknown) +=
895  lambda_sq *
896  this->node_pt(l2)->time_stepper_pt()->weight(2, 0) *
897  u_basis(l2) * u_test(l) * W;
898  }
899 
900  for (unsigned b = 0; b < DIM; b++)
901  {
902  for (unsigned d = 0; d < DIM; d++)
903  {
904  jacobian(local_eqn, local_unknown) +=
905  this->E(a, b, c, d) * du_basis_dx(l2, d) *
906  du_test_dx(l, b) * W;
907  }
908  }
909  }
910  }
911  }
912 
913  // d(u_eqn_l,a)/d(Q_l2)
914  for (unsigned l2 = 0; l2 < n_q_basis; l2++)
915  {
916  TimeStepper* timestepper_pt = 0;
917 
918  if (l2 < n_q_basis_edge)
919  {
920  local_unknown = q_edge_local_eqn(l2);
921  timestepper_pt =
922  this->node_pt(q_edge_node_number(l2))->time_stepper_pt();
923  }
924  else // n_q_basis_edge <= l < n_basis
925  {
926  local_unknown = q_internal_local_eqn(l2 - n_q_basis_edge);
927  timestepper_pt = this->internal_data_pt(q_internal_index())
928  ->time_stepper_pt();
929  }
930 
931  if (local_unknown >= 0)
932  {
933  jacobian(local_eqn, local_unknown) +=
934  lambda_sq * rho_f_over_rho * timestepper_pt->weight(1, 0) *
935  q_basis(l2, a) * u_test(l) * W;
936  }
937  }
938 
939  // d(u_eqn_l,a)/d(P_l2)
940  for (unsigned l2 = 0; l2 < n_p_basis; l2++)
941  {
942  local_unknown = p_local_eqn(l2);
943  if (local_unknown >= 0)
944  {
945  jacobian(local_eqn, local_unknown) -=
946  alpha * p_basis(l2) * du_test_dx(l, a) * W;
947  }
948  }
949  } // End of Jacobian entries
950  } // End of if not boundary condition
951  } // End of loop over dimensions
952  } // End of loop over u test functions
953 
954 
955  // Darcy:
956  // Loop over the test functions
957  for (unsigned l = 0; l < n_q_basis; l++)
958  {
959  if (l < n_q_basis_edge)
960  {
961  local_eqn = q_edge_local_eqn(l);
962  }
963  else // n_q_basis_edge <= l < n_basis
964  {
965  local_eqn = q_internal_local_eqn(l - n_q_basis_edge);
966  }
967 
968  // If it's not a boundary condition
969  if (local_eqn >= 0)
970  {
971  for (unsigned i = 0; i < DIM; i++)
972  {
973  residuals[local_eqn] +=
974  (k_inv * interpolated_q[i] - rho_f_over_rho * f_fluid[i] +
975  rho_f_over_rho * lambda_sq *
976  (interpolated_d2u_dt2[i] +
977  (1.0 / porosity) * interpolated_dq_dt[i])) *
978  q_test(l, i) * w * J;
979  }
980 
981  // deliberately no jacobian factor in this integral
982  residuals[local_eqn] -= (interpolated_p * div_q_test_ds(l)) * w;
983 
984  // Jacobian entries
985  if (flag)
986  {
987  // d(q_eqn_l)/d(U_l2,c)
988  for (unsigned l2 = 0; l2 < n_node; l2++)
989  {
990  for (unsigned c = 0; c < DIM; c++)
991  {
992  local_unknown = this->nodal_local_eqn(l2, u_index(c));
993  if (local_unknown >= 0)
994  {
995  jacobian(local_eqn, local_unknown) +=
996  rho_f_over_rho * lambda_sq *
997  this->node_pt(l2)->time_stepper_pt()->weight(2, 0) *
998  u_basis(l2) * q_test(l, c) * W;
999  }
1000  }
1001  }
1002 
1003  // d(q_eqn_l)/d(Q_l2)
1004  for (unsigned l2 = 0; l2 < n_q_basis; l2++)
1005  {
1006  TimeStepper* timestepper_pt = 0;
1007 
1008  if (l2 < n_q_basis_edge)
1009  {
1010  local_unknown = q_edge_local_eqn(l2);
1011  timestepper_pt =
1012  this->node_pt(q_edge_node_number(l2))->time_stepper_pt();
1013  }
1014  else // n_q_basis_edge <= l < n_basis
1015  {
1016  local_unknown = q_internal_local_eqn(l2 - n_q_basis_edge);
1017  timestepper_pt =
1018  this->internal_data_pt(q_internal_index())->time_stepper_pt();
1019  }
1020 
1021  if (local_unknown >= 0)
1022  {
1023  for (unsigned c = 0; c < DIM; c++)
1024  {
1025  jacobian(local_eqn, local_unknown) +=
1026  q_basis(l2, c) * q_test(l, c) *
1027  (k_inv + rho_f_over_rho * lambda_sq *
1028  timestepper_pt->weight(1, 0) / porosity) *
1029  W;
1030  }
1031  }
1032  }
1033 
1034  // d(q_eqn_l)/d(P_l2)
1035  for (unsigned l2 = 0; l2 < n_p_basis; l2++)
1036  {
1037  local_unknown = p_local_eqn(l2);
1038 
1039  // deliberately no jacobian factor in this term
1040  jacobian(local_eqn, local_unknown) -=
1041  p_basis(l2) * div_q_test_ds(l) * w;
1042  }
1043  } // End of Jacobian entries
1044  } // End of if not boundary condition
1045  } // End of loop over q test functions
1046 
1047  // loop over pressure test functions
1048  for (unsigned l = 0; l < n_p_basis; l++)
1049  {
1050  // get the local equation number
1051  local_eqn = p_local_eqn(l);
1052 
1053  // If it's not a boundary condition
1054  if (local_eqn >= 0)
1055  {
1056  // solid divergence term
1057  residuals[local_eqn] +=
1058  alpha * interpolated_div_du_dt_dx * p_test(l) * w * J;
1059  // fluid divergence term - deliberately no jacobian factor in this
1060  // term
1061  residuals[local_eqn] += interpolated_div_q_ds * p_test(l) * w;
1062  // mass source term
1063  residuals[local_eqn] -= mass_source_local * p_test(l) * w * J;
1064 
1065  // Jacobian entries
1066  if (flag)
1067  {
1068  // d(p_eqn_l)/d(U_l2,c)
1069  for (unsigned l2 = 0; l2 < n_node; l2++)
1070  {
1071  for (unsigned c = 0; c < DIM; c++)
1072  {
1073  local_unknown = this->nodal_local_eqn(l2, u_index(c));
1074 
1075  if (local_unknown >= 0)
1076  {
1077  jacobian(local_eqn, local_unknown) +=
1078  alpha * this->node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1079  du_basis_dx(l2, c) * p_test(l) * W;
1080  }
1081  }
1082  }
1083 
1084  // d(p_eqn_l)/d(Q_l2)
1085  for (unsigned l2 = 0; l2 < n_q_basis; l2++)
1086  {
1087  if (l2 < n_q_basis_edge)
1088  {
1089  local_unknown = q_edge_local_eqn(l2);
1090  }
1091  else // n_q_basis_edge <= l < n_basis
1092  {
1093  local_unknown = q_internal_local_eqn(l2 - n_q_basis_edge);
1094  }
1095 
1096  if (local_unknown >= 0)
1097  {
1098  jacobian(local_eqn, local_unknown) +=
1099  p_test(l) * div_q_basis_ds(l2) * w;
1100  }
1101  }
1102  } // End of Jacobian entries
1103  } // End of if not boundary condition
1104  } // End of loop over p test functions
1105  } // End of loop over integration points
1106  }
b
Scalar * b
Definition: benchVecAdd.cpp:17
oomph::FiniteElement::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Definition: elements.h:1432
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::FiniteElement::nodal_position
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317
oomph::FiniteElement::raw_nodal_value
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Definition: elements.h:2576
oomph::PoroelasticityEquations::dq_internal_dt
double dq_internal_dt(const unsigned &n) const
dq_internal/dt for the n-th internal degree of freedom
Definition: poroelasticity_elements.h:663
oomph::PoroelasticityEquations::porosity
const double & porosity() const
Access function for the porosity.
Definition: poroelasticity_elements.h:121
oomph::PoroelasticityEquations::q_edge_local_eqn
virtual int q_edge_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th edge (flux) degree of freedom.
oomph::PoroelasticityEquations::shape_basis_test_local_at_knot
virtual double shape_basis_test_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsi, Shape &u_basis, Shape &u_test, DShape &du_basis_dx, DShape &du_test_dx, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const =0
oomph::PoroelasticityEquations::k_inv
const double & k_inv() const
Access function for the nondim inverse permeability.
Definition: poroelasticity_elements.h:97
oomph::PoroelasticityEquations::force_solid
void force_solid(const double &time, const Vector< double > &x, Vector< double > &b) const
Definition: poroelasticity_elements.h:170
oomph::PoroelasticityEquations::mass_source
void mass_source(const double &time, const Vector< double > &x, double &b) const
Definition: poroelasticity_elements.h:214
oomph::PoroelasticityEquations::alpha
const double & alpha() const
Access function for alpha.
Definition: poroelasticity_elements.h:109
oomph::PoroelasticityEquations::q_internal_local_eqn
virtual int q_internal_local_eqn(const unsigned &n) const =0
oomph::PoroelasticityEquations::np_basis
virtual unsigned np_basis() const =0
Return the total number of pressure basis functions.
oomph::PoroelasticityEquations::dq_edge_dt
double dq_edge_dt(const unsigned &n) const
dq_edge/dt for the n-th edge degree of freedom
Definition: poroelasticity_elements.h:635
oomph::PoroelasticityEquations::density_ratio
const double & density_ratio() const
Access function for the density ratio.
Definition: poroelasticity_elements.h:85
oomph::PoroelasticityEquations::nq_basis_edge
virtual unsigned nq_basis_edge() const =0
Return the number of edge basis functions for q.
oomph::PoroelasticityEquations::p_local_eqn
virtual int p_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th pressure degree of freedom.
oomph::PoroelasticityEquations::force_fluid
void force_fluid(const double &time, const Vector< double > &x, Vector< double > &b) const
Definition: poroelasticity_elements.h:192
oomph::PoroelasticityEquations::p_value
virtual double p_value(unsigned &n) const =0
Return the nth pressure value.
oomph::PoroelasticityEquations::lambda_sq
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
Definition: poroelasticity_elements.h:73
oomph::PoroelasticityEquations::E
const double E(const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l) const
Access function to the entries in the elasticity tensor.
Definition: poroelasticity_elements.h:236
oomph::PoroelasticityEquations::nq_basis
virtual unsigned nq_basis() const =0
Return the total number of computational basis functions for q.
oomph::TimeStepper::time_pt
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
oomph::Time::time
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
a
const Scalar * a
Definition: level2_cplx_impl.h:32
calibrate.c
int c
Definition: calibrate.py:100

References DIM, and s.

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

◆ force_fluid()

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::force_fluid ( const double &  time,
const Vector< double > &  x,
Vector< double > &  b 
) const
inline

Indirect access to the fluid forcing function - returns 0 if no forcing function has been set

195  {
196  // If no function has been set, return zero vector
197  if (Force_fluid_fct_pt == 0)
198  {
199  // Get spatial dimension of element
200  unsigned n = dim();
201  for (unsigned i = 0; i < n; i++)
202  {
203  b[i] = 0.0;
204  }
205  }
206  else
207  {
208  (*Force_fluid_fct_pt)(time, x, b);
209  }
210  }
oomph::FiniteElement::dim
unsigned dim() const
Definition: elements.h:2611

References b, oomph::FiniteElement::dim(), oomph::PoroelasticityEquations< DIM >::Force_fluid_fct_pt, i, n, and plotDoE::x.

◆ force_fluid_fct_pt() [1/2]

template<unsigned DIM>
SourceFctPt& oomph::PoroelasticityEquations< DIM >::force_fluid_fct_pt ( )
inline

Access function: Pointer to fluid force function.

146  {
147  return Force_fluid_fct_pt;
148  }

References oomph::PoroelasticityEquations< DIM >::Force_fluid_fct_pt.

◆ force_fluid_fct_pt() [2/2]

template<unsigned DIM>
SourceFctPt oomph::PoroelasticityEquations< DIM >::force_fluid_fct_pt ( ) const
inline

Access function: Pointer to fluid force function (const version)

152  {
153  return Force_fluid_fct_pt;
154  }

References oomph::PoroelasticityEquations< DIM >::Force_fluid_fct_pt.

◆ force_solid()

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::force_solid ( const double &  time,
const Vector< double > &  x,
Vector< double > &  b 
) const
inline

Indirect access to the solid force function - returns 0 if no forcing function has been set

173  {
174  // If no function has been set, return zero vector
175  if (Force_solid_fct_pt == 0)
176  {
177  // Get spatial dimension of element
178  unsigned n = dim();
179  for (unsigned i = 0; i < n; i++)
180  {
181  b[i] = 0.0;
182  }
183  }
184  else
185  {
186  (*Force_solid_fct_pt)(time, x, b);
187  }
188  }

References b, oomph::FiniteElement::dim(), oomph::PoroelasticityEquations< DIM >::Force_solid_fct_pt, i, n, and plotDoE::x.

◆ force_solid_fct_pt() [1/2]

template<unsigned DIM>
SourceFctPt& oomph::PoroelasticityEquations< DIM >::force_solid_fct_pt ( )
inline

Access function: Pointer to solid force function.

134  {
135  return Force_solid_fct_pt;
136  }

References oomph::PoroelasticityEquations< DIM >::Force_solid_fct_pt.

◆ force_solid_fct_pt() [2/2]

template<unsigned DIM>
SourceFctPt oomph::PoroelasticityEquations< DIM >::force_solid_fct_pt ( ) const
inline

Access function: Pointer to solid force function (const version)

140  {
141  return Force_solid_fct_pt;
142  }

References oomph::PoroelasticityEquations< DIM >::Force_solid_fct_pt.

◆ get_div_q_basis_local()

template<unsigned DIM>
virtual void oomph::PoroelasticityEquations< DIM >::get_div_q_basis_local ( const Vector< double > &  s,
Shape &  div_q_basis_ds 
) const
pure virtual

Returns the local form of the q basis and dbasis/ds at local coordinate s

Implemented in oomph::TPoroelasticityElement< ORDER >, oomph::TPoroelasticityElement< ORDER >, and oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::interpolated_div_q().

◆ get_p_basis()

template<unsigned DIM>
virtual void oomph::PoroelasticityEquations< DIM >::get_p_basis ( const Vector< double > &  s,
Shape &  p_basis 
) const
pure virtual

Return the pressure basis.

Implemented in oomph::TPoroelasticityElement< ORDER >, oomph::TPoroelasticityElement< ORDER >, and oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::interpolated_p().

◆ get_q_basis()

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::get_q_basis ( const Vector< double > &  s,
Shape &  q_basis 
) const
inline

Returns the transformed basis at local coordinate s.

305  {
306  const unsigned n_node = this->nnode();
307  Shape psi(n_node, DIM);
308  const unsigned n_q_basis = this->nq_basis();
309  Shape q_basis_local(n_q_basis, DIM);
310  this->get_q_basis_local(s, q_basis_local);
311  (void)this->transform_basis(s, q_basis_local, psi, q_basis);
312  }
oomph::PoroelasticityEquations::transform_basis
double transform_basis(const Vector< double > &s, const Shape &q_basis_local, Shape &psi, DShape &dpsi, Shape &q_basis) const
Definition: poroelasticity_elements.cc:196
oomph::PoroelasticityEquations::get_q_basis_local
virtual void get_q_basis_local(const Vector< double > &s, Shape &q_basis) const =0
Returns the local form of the q basis at local coordinate s.

References DIM, oomph::PoroelasticityEquations< DIM >::get_q_basis_local(), oomph::FiniteElement::nnode(), oomph::PoroelasticityEquations< DIM >::nq_basis(), and oomph::PoroelasticityEquations< DIM >::transform_basis().

Referenced by oomph::PoroelasticityEquations< DIM >::interpolated_q().

◆ get_q_basis_local()

template<unsigned DIM>
virtual void oomph::PoroelasticityEquations< DIM >::get_q_basis_local ( const Vector< double > &  s,
Shape &  q_basis 
) const
pure virtual

Returns the local form of the q basis at local coordinate s.

Implemented in oomph::TPoroelasticityElement< ORDER >, oomph::TPoroelasticityElement< ORDER >, and oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::get_q_basis().

◆ get_strain()

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::get_strain ( const Vector< double > &  s,
DenseMatrix< double > &  strain 
) const

Return the strain tensor.

Compute the strain tensor at local coordinate s.

Now fill in the entries of the strain tensor

56  {
57 #ifdef PARANOID
58  if ((strain.ncol() != DIM) || (strain.nrow() != DIM))
59  {
60  std::ostringstream error_message;
61  error_message << "Strain matrix is " << strain.ncol() << " x "
62  << strain.nrow() << ", but dimension of the equations is "
63  << DIM << std::endl;
64  throw OomphLibError(
65  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
66  }
67 
68  /*
69  //Find out how many position types there are
70  unsigned n_position_type = this->nnodal_position_type();
71 
72  if(n_position_type != 1)
73  {
74  throw OomphLibError(
75  "PoroElasticity is not yet implemented for more than one position type",
76  OOMPH_CURRENT_FUNCTION,
77  OOMPH_EXCEPTION_LOCATION);
78  }
79  */
80 #endif
81 
82 
83  // Find out how many nodes there are in the element
84  unsigned n_node = nnode();
85 
86  // Find the indices at which the local velocities are stored
87  unsigned u_nodal_index[DIM];
88  for (unsigned i = 0; i < DIM; i++)
89  {
90  u_nodal_index[i] = u_index(i);
91  }
92 
93  // Set up memory for the shape and derivative functions
94  Shape psi(n_node);
95  DShape dpsidx(n_node, DIM);
96 
97  // Call the derivatives of the shape functions
98  (void)dshape_eulerian(s, psi, dpsidx);
99 
100  // Calculate interpolated values of the derivative of global position
101  DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
102 
103  // Loop over nodes
104  for (unsigned l = 0; l < n_node; l++)
105  {
106  // Loop over velocity components
107  for (unsigned i = 0; i < DIM; i++)
108  {
109  // Get the nodal value
110  const double u_value = this->nodal_value(l, u_nodal_index[i]);
111 
112  // Loop over derivative directions
113  for (unsigned j = 0; j < DIM; j++)
114  {
115  interpolated_dudx(i, j) += u_value * dpsidx(l, j);
116  }
117  }
118  }
119 
121  for (unsigned i = 0; i < DIM; i++)
122  {
123  // Do upper half of matrix
124  // Note that j must be signed here for the comparison test to work
125  // Also i must be cast to an int
126  for (int j = (DIM - 1); j >= static_cast<int>(i); j--)
127  {
128  // Off diagonal terms
129  if (static_cast<int>(i) != j)
130  {
131  strain(i, j) =
132  0.5 * (interpolated_dudx(i, j) + interpolated_dudx(j, i));
133  }
134  // Diagonal terms will including growth factor when it comes back in
135  else
136  {
137  strain(i, i) = interpolated_dudx(i, i);
138  }
139  }
140  // Matrix is symmetric so just copy lower half
141  for (int j = (i - 1); j >= 0; j--)
142  {
143  strain(i, j) = strain(j, i);
144  }
145  }
146  }
oomph::DenseMatrix::nrow
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:485
oomph::DenseMatrix::ncol
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:491
oomph::FiniteElement::dshape_eulerian
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Definition: elements.cc:3298
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References DIM, i, j, oomph::DenseMatrix< T >::ncol(), oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

◆ get_stress()

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::get_stress ( const Vector< double > &  s,
DenseMatrix< double > &  stress 
) const

Return the Cauchy stress tensor, as calculated from the elasticity tensor at specified local coordinate

Compute the Cauchy stress tensor at local coordinate s for displacement formulation.

156  {
157 #ifdef PARANOID
158  if ((stress.ncol() != DIM) || (stress.nrow() != DIM))
159  {
160  std::ostringstream error_message;
161  error_message << "Stress matrix is " << stress.ncol() << " x "
162  << stress.nrow() << ", but dimension of the equations is "
163  << DIM << std::endl;
164  throw OomphLibError(
165  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
166  }
167 #endif
168 
169  // Get strain
170  DenseMatrix<double> strain(DIM, DIM);
171  this->get_strain(s, strain);
172 
173  // Now fill in the entries of the stress tensor without exploiting
174  // symmetry -- sorry too lazy. This fct is only used for
175  // postprocessing anyway...
176  for (unsigned i = 0; i < DIM; i++)
177  {
178  for (unsigned j = 0; j < DIM; j++)
179  {
180  stress(i, j) = 0.0;
181  for (unsigned k = 0; k < DIM; k++)
182  {
183  for (unsigned l = 0; k < DIM; k++)
184  {
185  stress(i, j) += this->E(i, j, k, l) * strain(k, l);
186  }
187  }
188  }
189  }
190  }
oomph::PoroelasticityEquations::get_strain
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
Definition: poroelasticity_elements.cc:54

References DIM, Global_Physical_Variables::E, i, j, k, oomph::DenseMatrix< T >::ncol(), oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

◆ interpolated_div_q() [1/2]

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::interpolated_div_q ( const Vector< double > &  s)
inline

Calculate the FE representation of div q and return it.

528  {
529  // Temporary storage for div u
530  double div_q = 0;
531 
532  // Get the intepolated divergence
533  interpolated_div_q(s, div_q);
534 
535  // Return it
536  return div_q;
537  }

References oomph::PoroelasticityEquations< DIM >::interpolated_div_q(), and s.

◆ interpolated_div_q() [2/2]

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::interpolated_div_q ( const Vector< double > &  s,
double &  div_q 
) const
inline

Calculate the FE representation of div u.

483  {
484  // Zero the divergence
485  div_q = 0;
486 
487  // Get the number of nodes, q basis function, and q edge basis functions
488  unsigned n_node = nnode();
489  const unsigned n_q_basis = nq_basis();
490  const unsigned n_q_basis_edge = nq_basis_edge();
491 
492  // Storage for the divergence basis
493  Shape div_q_basis_ds(n_q_basis);
494 
495  // Storage for the geometric basis and it's derivatives
496  Shape psi(n_node);
497  DShape dpsi(n_node, DIM);
498 
499  // Call the geometric shape functions and their derivatives
500  this->dshape_local(s, psi, dpsi);
501 
502  // Storage for the inverse of the geometric jacobian (just so we can call
503  // the local to eulerian mapping)
504  DenseMatrix<double> inverse_jacobian(DIM);
505 
506  // Get the determinant of the geometric mapping
507  double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
508 
509  // Get the divergence basis (wrt local coords) at local coords s
510  get_div_q_basis_local(s, div_q_basis_ds);
511 
512  // Add the contribution to the divergence from the edge basis functions
513  for (unsigned l = 0; l < n_q_basis_edge; l++)
514  {
515  div_q += 1.0 / det * div_q_basis_ds(l) * q_edge(l);
516  }
517 
518  // Add the contribution to the divergence from the internal basis
519  // functions
520  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
521  {
522  div_q += 1.0 / det * div_q_basis_ds(l) * q_internal(l - n_q_basis_edge);
523  }
524  }
oomph::FiniteElement::local_to_eulerian_mapping
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Definition: elements.h:1508
oomph::FiniteElement::dshape_local
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Definition: elements.h:1981
oomph::PoroelasticityEquations::get_div_q_basis_local
virtual void get_div_q_basis_local(const Vector< double > &s, Shape &div_q_basis_ds) const =0

References DIM, oomph::FiniteElement::dshape_local(), oomph::PoroelasticityEquations< DIM >::get_div_q_basis_local(), oomph::FiniteElement::local_to_eulerian_mapping(), oomph::FiniteElement::nnode(), oomph::PoroelasticityEquations< DIM >::nq_basis(), oomph::PoroelasticityEquations< DIM >::nq_basis_edge(), oomph::PoroelasticityEquations< DIM >::q_edge(), oomph::PoroelasticityEquations< DIM >::q_internal(), and s.

Referenced by oomph::PoroelasticityEquations< DIM >::interpolated_div_q().

◆ interpolated_p() [1/2]

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::interpolated_p ( const Vector< double > &  s) const
inline

Calculate the FE representation of p and return it.

563  {
564  // Temporary storage for p
565  double p = 0;
566 
567  // Get the interpolated pressure
568  interpolated_p(s, p);
569 
570  // Return it
571  return p;
572  }
p
float * p
Definition: Tutorial_Map_using.cpp:9

References oomph::PoroelasticityEquations< DIM >::interpolated_p(), p, and s.

◆ interpolated_p() [2/2]

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::interpolated_p ( const Vector< double > &  s,
double &  p 
) const
inline

Calculate the FE representation of p.

541  {
542  // Get the number of p basis functions
543  unsigned n_p_basis = np_basis();
544 
545  // Storage for the p basis
546  Shape p_basis(n_p_basis);
547 
548  // Call the p basis
549  get_p_basis(s, p_basis);
550 
551  // Zero the pressure
552  p = 0;
553 
554  // Add the contribution to the pressure from each basis function
555  for (unsigned l = 0; l < n_p_basis; l++)
556  {
557  p += p_value(l) * p_basis(l);
558  }
559  }
oomph::PoroelasticityEquations::get_p_basis
virtual void get_p_basis(const Vector< double > &s, Shape &p_basis) const =0
Return the pressure basis.

References oomph::PoroelasticityEquations< DIM >::get_p_basis(), oomph::PoroelasticityEquations< DIM >::np_basis(), p, oomph::PoroelasticityEquations< DIM >::p_value(), and s.

Referenced by oomph::PoroelasticityEquations< DIM >::interpolated_p().

◆ interpolated_q() [1/2]

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

Calculate the FE representation of the i-th component of q.

461  {
462  unsigned n_q_basis = nq_basis();
463  unsigned n_q_basis_edge = nq_basis_edge();
464 
465  Shape q_basis(n_q_basis, DIM);
466 
467  get_q_basis(s, q_basis);
468  double q_i = 0.0;
469  for (unsigned l = 0; l < n_q_basis_edge; l++)
470  {
471  q_i += q_edge(l) * q_basis(l, i);
472  }
473  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
474  {
475  q_i += q_internal(l - n_q_basis_edge) * q_basis(l, i);
476  }
477 
478  return q_i;
479  }
oomph::PoroelasticityEquations::get_q_basis
void get_q_basis(const Vector< double > &s, Shape &q_basis) const
Returns the transformed basis at local coordinate s.
Definition: poroelasticity_elements.h:304

References DIM, oomph::PoroelasticityEquations< DIM >::get_q_basis(), i, oomph::PoroelasticityEquations< DIM >::nq_basis(), oomph::PoroelasticityEquations< DIM >::nq_basis_edge(), oomph::PoroelasticityEquations< DIM >::q_edge(), oomph::PoroelasticityEquations< DIM >::q_internal(), and s.

◆ interpolated_q() [2/2]

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::interpolated_q ( const Vector< double > &  s,
Vector< double > &  u 
) const
inline

Calculate the FE representation of q.

438  {
439  unsigned n_q_basis = nq_basis();
440  unsigned n_q_basis_edge = nq_basis_edge();
441 
442  Shape q_basis(n_q_basis, DIM);
443 
444  get_q_basis(s, q_basis);
445  for (unsigned i = 0; i < DIM; i++)
446  {
447  u[i] = 0.0;
448  for (unsigned l = 0; l < n_q_basis_edge; l++)
449  {
450  u[i] += q_edge(l) * q_basis(l, i);
451  }
452  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
453  {
454  u[i] += q_internal(l - n_q_basis_edge) * q_basis(l, i);
455  }
456  }
457  }

References DIM, oomph::PoroelasticityEquations< DIM >::get_q_basis(), i, oomph::PoroelasticityEquations< DIM >::nq_basis(), oomph::PoroelasticityEquations< DIM >::nq_basis_edge(), oomph::PoroelasticityEquations< DIM >::q_edge(), oomph::PoroelasticityEquations< DIM >::q_internal(), and s.

◆ interpolated_u() [1/2]

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

Calculate the FE representation of the i-th component of u.

411  {
412  // Find number of nodes
413  unsigned n_node = nnode();
414 
415  // Local shape function
416  Shape psi(n_node);
417 
418  // Find values of shape function
419  shape(s, psi);
420 
421  // Get nodal index at which i-th velocity is stored
422  unsigned u_nodal_index = u_index(i);
423 
424  // Initialise value of u
425  double interpolated_u = 0.0;
426 
427  // Loop over the local nodes and sum
428  for (unsigned l = 0; l < n_node; l++)
429  {
430  interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
431  }
432 
433  return (interpolated_u);
434  }
oomph::FiniteElement::shape
virtual void shape(const Vector< double > &s, Shape &psi) const =0

References i, oomph::PoroelasticityEquations< DIM >::interpolated_u(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nodal_value(), s, oomph::FiniteElement::shape(), and oomph::PoroelasticityEquations< DIM >::u_index().

◆ interpolated_u() [2/2]

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::interpolated_u ( const Vector< double > &  s,
Vector< double > &  disp 
) const
inline

Calculate the FE representation of u.

383  {
384  // Find number of nodes
385  unsigned n_node = nnode();
386 
387  // Local shape function
388  Shape psi(n_node);
389 
390  // Find values of shape function
391  shape(s, psi);
392 
393  for (unsigned i = 0; i < DIM; i++)
394  {
395  // Index at which the nodal value is stored
396  unsigned u_nodal_index = u_index(i);
397 
398  // Initialise value of u
399  disp[i] = 0.0;
400 
401  // Loop over the local nodes and sum
402  for (unsigned l = 0; l < n_node; l++)
403  {
404  disp[i] += nodal_value(l, u_nodal_index) * psi[l];
405  }
406  }
407  }

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

Referenced by oomph::PoroelasticityEquations< DIM >::interpolated_u().

◆ k_inv()

template<unsigned DIM>
const double& oomph::PoroelasticityEquations< DIM >::k_inv ( ) const
inline

Access function for the nondim inverse permeability.

98  {
99  return *K_inv_pt;
100  }

References oomph::PoroelasticityEquations< DIM >::K_inv_pt.

◆ k_inv_pt()

template<unsigned DIM>
double*& oomph::PoroelasticityEquations< DIM >::k_inv_pt ( )
inline

Access function for pointer to the nondim inverse permeability.

104  {
105  return K_inv_pt;
106  }

References oomph::PoroelasticityEquations< DIM >::K_inv_pt.

◆ lambda_sq()

template<unsigned DIM>
const double& oomph::PoroelasticityEquations< DIM >::lambda_sq ( ) const
inline

Access function for timescale ratio (nondim density)

74  {
75  return *Lambda_sq_pt;
76  }

References oomph::PoroelasticityEquations< DIM >::Lambda_sq_pt.

◆ lambda_sq_pt()

template<unsigned DIM>
double*& oomph::PoroelasticityEquations< DIM >::lambda_sq_pt ( )
inline

Access function for pointer to timescale ratio (nondim density)

80  {
81  return Lambda_sq_pt;
82  }

References oomph::PoroelasticityEquations< DIM >::Lambda_sq_pt.

◆ mass_source()

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::mass_source ( const double &  time,
const Vector< double > &  x,
double &  b 
) const
inline

Indirect access to the mass source function - returns 0 if no mass source function has been set

217  {
218  // If no function has been set, return zero vector
219  if (Mass_source_fct_pt == 0)
220  {
221  b = 0.0;
222  }
223  else
224  {
225  (*Mass_source_fct_pt)(time, x, b);
226  }
227  }

References b, oomph::PoroelasticityEquations< DIM >::Mass_source_fct_pt, and plotDoE::x.

◆ mass_source_fct_pt() [1/2]

template<unsigned DIM>
MassSourceFctPt& oomph::PoroelasticityEquations< DIM >::mass_source_fct_pt ( )
inline

Access function: Pointer to mass source function.

158  {
159  return Mass_source_fct_pt;
160  }

References oomph::PoroelasticityEquations< DIM >::Mass_source_fct_pt.

◆ mass_source_fct_pt() [2/2]

template<unsigned DIM>
MassSourceFctPt oomph::PoroelasticityEquations< DIM >::mass_source_fct_pt ( ) const
inline

Access function: Pointer to mass source function (const version)

164  {
165  return Mass_source_fct_pt;
166  }

References oomph::PoroelasticityEquations< DIM >::Mass_source_fct_pt.

◆ nedge_gauss_point()

template<unsigned DIM>
virtual unsigned oomph::PoroelasticityEquations< DIM >::nedge_gauss_point ( ) const
pure virtual

Returns the number of gauss points along each edge of the element.

Implemented in oomph::TPoroelasticityElement< ORDER >, oomph::TPoroelasticityElement< ORDER >, and oomph::TPoroelasticityElement< ORDER >.

◆ np_basis()

template<unsigned DIM>
virtual unsigned oomph::PoroelasticityEquations< DIM >::np_basis ( ) const
pure virtual

Return the total number of pressure basis functions.

Implemented in oomph::TPoroelasticityElement< ORDER >, oomph::TPoroelasticityElement< ORDER >, and oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::interpolated_p().

◆ nq_basis()

template<unsigned DIM>
virtual unsigned oomph::PoroelasticityEquations< DIM >::nq_basis ( ) const
pure virtual

Return the total number of computational basis functions for q.

Implemented in oomph::TPoroelasticityElement< ORDER >, oomph::TPoroelasticityElement< ORDER >, and oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::get_q_basis(), oomph::PoroelasticityEquations< DIM >::interpolated_div_q(), and oomph::PoroelasticityEquations< DIM >::interpolated_q().

◆ nq_basis_edge()

template<unsigned DIM>
virtual unsigned oomph::PoroelasticityEquations< DIM >::nq_basis_edge ( ) const
pure virtual

Return the number of edge basis functions for q.

Implemented in oomph::TPoroelasticityElement< ORDER >, oomph::TPoroelasticityElement< ORDER >, and oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::interpolated_div_q(), and oomph::PoroelasticityEquations< DIM >::interpolated_q().

◆ output() [1/2]

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

Output with default number of plot points.

Reimplemented from oomph::FiniteElement.

707  {
708  unsigned nplot = 5;
709  output(outfile, nplot);
710  }
oomph::PoroelasticityEquations::output
void output(std::ostream &outfile)
Output with default number of plot points.
Definition: poroelasticity_elements.h:706

Referenced by oomph::TPoroelasticityElement< ORDER >::output().

◆ output() [2/2]

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

Output FE representation of soln: x,y,u1,u2,div_q,p at Nplot^DIM plot points

Output FE representation of soln: x,y,u1,u2,q1,q2,div_q,p at Nplot^DIM plot points

Reimplemented from oomph::FiniteElement.

262  {
263  // Vector of local coordinates
264  Vector<double> s(DIM);
265 
266  // Tecplot header info
267  outfile << tecplot_zone_string(nplot);
268 
269  // Loop over plot points
270  unsigned num_plot_points = nplot_points(nplot);
271 
272  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
273  {
274  // Get local coordinates of plot point
275  get_s_plot(iplot, nplot, s);
276 
277  // Output the components of the position
278  for (unsigned i = 0; i < DIM; i++)
279  {
280  outfile << interpolated_x(s, i) << " ";
281  }
282 
283  // Output the components of the FE representation of u at s
284  for (unsigned i = 0; i < DIM; i++)
285  {
286  outfile << interpolated_u(s, i) << " ";
287  }
288 
289  // Output the components of the FE representation of q at s
290  for (unsigned i = 0; i < DIM; i++)
291  {
292  outfile << interpolated_q(s, i) << " ";
293  }
294 
295  // Output FE representation of div u at s
296  outfile << interpolated_div_q(s) << " ";
297 
298  // Output FE representation of p at s
299  outfile << interpolated_p(s) << " ";
300 
301  const unsigned n_node = this->nnode();
302 
303  Shape psi(n_node);
304  shape(s, psi);
305 
306  Vector<double> interpolated_du_dt(DIM, 0.0);
307 
308  for (unsigned i = 0; i < DIM; i++)
309  {
310  for (unsigned l = 0; l < n_node; l++)
311  {
312  interpolated_du_dt[i] += du_dt(l, i) * psi(l);
313  }
314  outfile << interpolated_du_dt[i] << " ";
315  }
316 
317  Vector<double> interpolated_d2u_dt2(DIM, 0.0);
318 
319  for (unsigned i = 0; i < DIM; i++)
320  {
321  for (unsigned l = 0; l < n_node; l++)
322  {
323  interpolated_d2u_dt2[i] += d2u_dt2(l, i) * psi(l);
324  }
325  outfile << interpolated_d2u_dt2[i] << " ";
326  }
327 
328  const unsigned n_q_basis = this->nq_basis();
329  const unsigned n_q_basis_edge = this->nq_basis_edge();
330 
331  Shape q_basis(n_q_basis, DIM), q_basis_local(n_q_basis, DIM);
332  this->get_q_basis_local(s, q_basis_local);
333  (void)this->transform_basis(s, q_basis_local, psi, q_basis);
334 
335  Vector<double> interpolated_dq_dt(DIM, 0.0);
336 
337  for (unsigned i = 0; i < DIM; i++)
338  {
339  for (unsigned l = 0; l < n_q_basis_edge; l++)
340  {
341  interpolated_dq_dt[i] += dq_edge_dt(l) * q_basis(l, i);
342  }
343  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
344  {
345  interpolated_dq_dt[i] +=
346  dq_internal_dt(l - n_q_basis_edge) * q_basis(l, i);
347  }
348  outfile << interpolated_dq_dt[i] << " ";
349  }
350 
351  outfile << std::endl;
352  }
353 
354  // Write tecplot footer (e.g. FE connectivity lists)
355  this->write_tecplot_zone_footer(outfile, nplot);
356  }
oomph::FiniteElement::tecplot_zone_string
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Definition: elements.h:3161
oomph::FiniteElement::get_s_plot
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
oomph::FiniteElement::nplot_points
virtual unsigned nplot_points(const unsigned &nplot) const
Definition: elements.h:3186
oomph::FiniteElement::write_tecplot_zone_footer
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Definition: elements.h:3174

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

◆ output_fct() [1/2]

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

Output FE representation of exact soln: x,y,u1,u2,div_q,p at Nplot^DIM plot points. Unsteady version

Reimplemented from oomph::FiniteElement.

416  {
417  // Vector of local coordinates
418  Vector<double> s(DIM);
419 
420  // Vector for coordintes
421  Vector<double> x(DIM);
422 
423  // Tecplot header info
424  outfile << this->tecplot_zone_string(nplot);
425 
426  // Exact solution Vector
427  Vector<double> exact_soln(2 * DIM + 2);
428 
429  // Loop over plot points
430  unsigned num_plot_points = this->nplot_points(nplot);
431 
432  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
433  {
434  // Get local coordinates of plot point
435  this->get_s_plot(iplot, nplot, s);
436 
437  // Get x position as Vector
438  this->interpolated_x(s, x);
439 
440  // Get exact solution at this point
441  (*exact_soln_pt)(time, x, exact_soln);
442 
443  // Output x,y,q_exact,p_exact,div_q_exact
444  for (unsigned i = 0; i < DIM; i++)
445  {
446  outfile << x[i] << " ";
447  }
448  for (unsigned i = 0; i < 2 * DIM + 2; i++)
449  {
450  outfile << exact_soln[i] << " ";
451  }
452  outfile << std::endl;
453  }
454 
455  // Write tecplot footer (e.g. FE connectivity lists)
456  this->write_tecplot_zone_footer(outfile, nplot);
457  }

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

◆ output_fct() [2/2]

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

Output FE representation of exact soln: x,y,u1,u2,div_q,p at Nplot^DIM plot points

Output FE representation of exact soln: x,y,u1,u2,q1,q2,div_q,p at Nplot^DIM plot points

Reimplemented from oomph::FiniteElement.

365  {
366  // Vector of local coordinates
367  Vector<double> s(DIM);
368 
369  // Vector for coordintes
370  Vector<double> x(DIM);
371 
372  // Tecplot header info
373  outfile << this->tecplot_zone_string(nplot);
374 
375  // Exact solution Vector
376  Vector<double> exact_soln(2 * DIM + 2);
377 
378  // Loop over plot points
379  unsigned num_plot_points = this->nplot_points(nplot);
380 
381  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
382  {
383  // Get local coordinates of plot point
384  this->get_s_plot(iplot, nplot, s);
385 
386  // Get x position as Vector
387  this->interpolated_x(s, x);
388 
389  // Get exact solution at this point
390  (*exact_soln_pt)(x, exact_soln);
391 
392  // Output x,y,q_exact,p_exact,div_q_exact
393  for (unsigned i = 0; i < DIM; i++)
394  {
395  outfile << x[i] << " ";
396  }
397  for (unsigned i = 0; i < 2 * DIM + 2; i++)
398  {
399  outfile << exact_soln[i] << " ";
400  }
401  outfile << std::endl;
402  }
403 
404  // Write tecplot footer (e.g. FE connectivity lists)
405  this->write_tecplot_zone_footer(outfile, nplot);
406  }

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

◆ p_local_eqn()

template<unsigned DIM>
virtual int oomph::PoroelasticityEquations< DIM >::p_local_eqn ( const unsigned &  n) const
pure virtual

Return the equation number of the n-th pressure degree of freedom.

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ p_value()

template<unsigned DIM>
virtual double oomph::PoroelasticityEquations< DIM >::p_value ( unsigned &  n) const
pure virtual

Return the nth pressure value.

Implemented in oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::interpolated_p().

◆ pin_p_value()

template<unsigned DIM>
virtual void oomph::PoroelasticityEquations< DIM >::pin_p_value ( const unsigned &  n,
const double &  p 
)
pure virtual

Pin the nth pressure value.

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ pin_q_internal_value()

template<unsigned DIM>
virtual void oomph::PoroelasticityEquations< DIM >::pin_q_internal_value ( const unsigned &  n)
pure virtual

Pin the nth internal q value.

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ porosity()

template<unsigned DIM>
const double& oomph::PoroelasticityEquations< DIM >::porosity ( ) const
inline

Access function for the porosity.

122  {
123  return *Porosity_pt;
124  }

References oomph::PoroelasticityEquations< DIM >::Porosity_pt.

◆ porosity_pt()

template<unsigned DIM>
double*& oomph::PoroelasticityEquations< DIM >::porosity_pt ( )
inline

Access function for pointer to the porosity.

128  {
129  return Porosity_pt;
130  }

References oomph::PoroelasticityEquations< DIM >::Porosity_pt.

◆ q_edge() [1/2]

template<unsigned DIM>
virtual double oomph::PoroelasticityEquations< DIM >::q_edge ( const unsigned &  n) const
pure virtual

Return the values of the edge (flux) degrees of freedom.

Implemented in oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::dq_edge_dt(), oomph::PoroelasticityEquations< DIM >::interpolated_div_q(), and oomph::PoroelasticityEquations< DIM >::interpolated_q().

◆ q_edge() [2/2]

template<unsigned DIM>
virtual double oomph::PoroelasticityEquations< DIM >::q_edge ( const unsigned &  t,
const unsigned &  n 
) const
pure virtual

Return the values of the edge (flux) degrees of freedom at time history level t

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ q_edge_index()

template<unsigned DIM>
virtual unsigned oomph::PoroelasticityEquations< DIM >::q_edge_index ( const unsigned &  n) const
pure virtual

Return the nodal index at which the nth edge unknown is stored.

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ q_edge_local_eqn()

template<unsigned DIM>
virtual int oomph::PoroelasticityEquations< DIM >::q_edge_local_eqn ( const unsigned &  n) const
pure virtual

Return the equation number of the n-th edge (flux) degree of freedom.

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ q_edge_node_number()

template<unsigned DIM>
virtual unsigned oomph::PoroelasticityEquations< DIM >::q_edge_node_number ( const unsigned &  n) const
pure virtual

Return the number of the node where the nth edge unknown is stored.

Implemented in oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::dq_edge_dt().

◆ q_internal() [1/2]

template<unsigned DIM>
virtual double oomph::PoroelasticityEquations< DIM >::q_internal ( const unsigned &  n) const
pure virtual

Return the values of the internal (moment) degrees of freedom.

Implemented in oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::dq_internal_dt(), oomph::PoroelasticityEquations< DIM >::interpolated_div_q(), and oomph::PoroelasticityEquations< DIM >::interpolated_q().

◆ q_internal() [2/2]

template<unsigned DIM>
virtual double oomph::PoroelasticityEquations< DIM >::q_internal ( const unsigned &  t,
const unsigned &  n 
) const
pure virtual

Return the values of the internal (moment) degrees of freedom at time history level t

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ q_internal_index()

template<unsigned DIM>
virtual unsigned oomph::PoroelasticityEquations< DIM >::q_internal_index ( ) const
pure virtual

Return the index of the internal data where the q_internal degrees of freedom are stored

Implemented in oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::dq_internal_dt(), and oomph::PoroelasticityEquations< DIM >::set_q_internal_timestepper().

◆ q_internal_local_eqn()

template<unsigned DIM>
virtual int oomph::PoroelasticityEquations< DIM >::q_internal_local_eqn ( const unsigned &  n) const
pure virtual

Return the equation number of the n-th internal (moment) degree of freedom

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ required_nvalue()

template<unsigned DIM>
virtual unsigned oomph::PoroelasticityEquations< DIM >::required_nvalue ( const unsigned &  n) const
pure virtual

Number of values required at node n.

Reimplemented from oomph::FiniteElement.

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ scale_basis()

template<unsigned DIM>
virtual void oomph::PoroelasticityEquations< DIM >::scale_basis ( Shape &  basis) const
pure virtual

Scale the edge basis to allow arbitrary edge mappings.

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ self_test()

template<unsigned DIM>
unsigned oomph::PoroelasticityEquations< DIM >::self_test ( )
inlinevirtual

Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.

Self-test: Have all internal values been classified as pinned/unpinned? Has pointer to spatial integration scheme been set? Return 0 if OK.

Reimplemented from oomph::FiniteElement.

701  {
702  return 0;
703  }

◆ set_q_internal_timestepper()

template<unsigned DIM>
void oomph::PoroelasticityEquations< DIM >::set_q_internal_timestepper ( TimeStepper *const  time_stepper_pt)
inline

Set the timestepper of the q internal data object.

694  {
695  unsigned q_index = q_internal_index();
696 
697  this->internal_data_pt(q_index)->set_time_stepper(time_stepper_pt, false);
698  }
oomph::Data::set_time_stepper
void set_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Definition: nodes.cc:406

References oomph::GeneralisedElement::internal_data_pt(), oomph::PoroelasticityEquations< DIM >::q_internal_index(), oomph::Data::set_time_stepper(), and oomph::GeomObject::time_stepper_pt().

◆ shape_basis_test_local()

template<unsigned DIM>
virtual double oomph::PoroelasticityEquations< DIM >::shape_basis_test_local ( const Vector< double > &  s,
Shape &  psi,
DShape &  dpsi,
Shape &  u_basis,
Shape &  u_test,
DShape &  du_basis_dx,
DShape &  du_test_dx,
Shape &  q_basis,
Shape &  q_test,
Shape &  p_basis,
Shape &  p_test,
Shape &  div_q_basis_ds,
Shape &  div_q_test_ds 
) const
protectedpure virtual

Returns the geometric basis, and the q, p and divergence basis functions and test functions at local coordinate s

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ shape_basis_test_local_at_knot()

template<unsigned DIM>
virtual double oomph::PoroelasticityEquations< DIM >::shape_basis_test_local_at_knot ( const unsigned &  ipt,
Shape &  psi,
DShape &  dpsi,
Shape &  u_basis,
Shape &  u_test,
DShape &  du_basis_dx,
DShape &  du_test_dx,
Shape &  q_basis,
Shape &  q_test,
Shape &  p_basis,
Shape &  p_test,
Shape &  div_q_basis_ds,
Shape &  div_q_test_ds 
) const
protectedpure virtual

Returns the geometric basis, and the q, p and divergence basis functions and test functions at integration point ipt

Implemented in oomph::TPoroelasticityElement< ORDER >.

◆ transform_basis() [1/2]

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::transform_basis ( const Vector< double > &  s,
const Shape &  q_basis_local,
Shape &  psi,
DShape &  dpsi,
Shape &  q_basis 
) const

Performs a div-conserving transformation of the vector basis functions from the reference element to the actual element

202  {
203  // Call the (geometric) shape functions and their derivatives
204  //(void)this->dshape_eulerian(s,psi,dpsi);
205  this->dshape_local(s, psi, dpsi);
206 
207  // Storage for the (geometric) jacobian and its inverse
208  DenseMatrix<double> jacobian(DIM), inverse_jacobian(DIM);
209 
210  // Get the jacobian of the geometric mapping and its determinant
211  double det = local_to_eulerian_mapping(dpsi, jacobian, inverse_jacobian);
212 
213  // Transform the derivative of the geometric basis so that it's w.r.t.
214  // global coordinates
215  this->transform_derivatives(inverse_jacobian, dpsi);
216 
217  // Get the number of computational basis vectors
218  const unsigned n_q_basis = this->nq_basis();
219 
220  // Loop over the basis vectors
221  for (unsigned l = 0; l < n_q_basis; l++)
222  {
223  // Loop over the spatial components
224  for (unsigned i = 0; i < DIM; i++)
225  {
226  // Zero the basis
227  q_basis(l, i) = 0.0;
228  }
229  }
230 
231  // Loop over the spatial components
232  for (unsigned i = 0; i < DIM; i++)
233  {
234  // And again
235  for (unsigned j = 0; j < DIM; j++)
236  {
237  // Get the element of the jacobian (must transpose it due to different
238  // conventions) and divide by the determinant
239  double jac_trans = jacobian(j, i) / det;
240 
241  // Loop over the computational basis vectors
242  for (unsigned l = 0; l < n_q_basis; l++)
243  {
244  // Apply the appropriate div-conserving mapping
245  q_basis(l, i) += jac_trans * q_basis_local(l, j);
246  }
247  }
248  }
249 
250  // Scale the basis by the ratio of the length of the edge to the length of
251  // the corresponding edge on the reference element
252  scale_basis(q_basis);
253 
254  return det;
255  }
oomph::FiniteElement::transform_derivatives
virtual void transform_derivatives(const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
Definition: elements.cc:2833
oomph::PoroelasticityEquations::scale_basis
virtual void scale_basis(Shape &basis) const =0
Scale the edge basis to allow arbitrary edge mappings.

References DIM, i, j, and s.

Referenced by oomph::PoroelasticityEquations< DIM >::get_q_basis(), and oomph::PoroelasticityEquations< DIM >::transform_basis().

◆ transform_basis() [2/2]

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::transform_basis ( const Vector< double > &  s,
const Shape &  q_basis_local,
Shape &  psi,
Shape &  q_basis 
) const
inline

Performs a div-conserving transformation of the vector basis functions from the reference element to the actual element

361  {
362  const unsigned n_node = this->nnode();
363  DShape dpsi(n_node, DIM);
364  return transform_basis(s, q_basis_local, psi, dpsi, q_basis);
365  }

References DIM, oomph::FiniteElement::nnode(), s, and oomph::PoroelasticityEquations< DIM >::transform_basis().

◆ u_index()

template<unsigned DIM>
virtual unsigned oomph::PoroelasticityEquations< DIM >::u_index ( const unsigned &  n) const
pure virtual

Return the nodal index of the n-th solid displacement unknown.

Implemented in oomph::TPoroelasticityElement< ORDER >.

Referenced by oomph::PoroelasticityEquations< DIM >::d2u_dt2(), oomph::PoroelasticityEquations< DIM >::du_dt(), and oomph::PoroelasticityEquations< DIM >::interpolated_u().

Member Data Documentation

◆ Alpha_pt

template<unsigned DIM>
double* oomph::PoroelasticityEquations< DIM >::Alpha_pt
private

Alpha.

Referenced by oomph::PoroelasticityEquations< DIM >::alpha(), and oomph::PoroelasticityEquations< DIM >::alpha_pt().

◆ Default_alpha_value

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::Default_alpha_value = 1.0
staticprivate

Static default value for alpha.

◆ Default_density_ratio_value

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::Default_density_ratio_value = 1.0
staticprivate

Static default value for the density ratio.

◆ Default_k_inv_value

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::Default_k_inv_value = 1.0
staticprivate

Static default value for 1/k.

◆ Default_lambda_sq_value

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::Default_lambda_sq_value = 1.0
staticprivate

Static default value for timescale ratio (1.0 – for natural scaling)

◆ Default_porosity_value

template<unsigned DIM>
double oomph::PoroelasticityEquations< DIM >::Default_porosity_value = 1.0
staticprivate

Static default value for the porosity.

◆ Density_ratio_pt

template<unsigned DIM>
double* oomph::PoroelasticityEquations< DIM >::Density_ratio_pt
private

Density ratio.

Referenced by oomph::PoroelasticityEquations< DIM >::density_ratio(), and oomph::PoroelasticityEquations< DIM >::density_ratio_pt().

◆ Elasticity_tensor_pt

template<unsigned DIM>
ElasticityTensor* oomph::PoroelasticityEquations< DIM >::Elasticity_tensor_pt
protected

Pointer to the elasticity tensor.

Referenced by oomph::PoroelasticityEquations< DIM >::E(), and oomph::PoroelasticityEquations< DIM >::elasticity_tensor_pt().

◆ Force_fluid_fct_pt

template<unsigned DIM>
SourceFctPt oomph::PoroelasticityEquations< DIM >::Force_fluid_fct_pt
private

Pointer to fluid source function.

Referenced by oomph::PoroelasticityEquations< DIM >::force_fluid(), and oomph::PoroelasticityEquations< DIM >::force_fluid_fct_pt().

◆ Force_solid_fct_pt

template<unsigned DIM>
SourceFctPt oomph::PoroelasticityEquations< DIM >::Force_solid_fct_pt
private

Pointer to solid source function.

Referenced by oomph::PoroelasticityEquations< DIM >::force_solid(), and oomph::PoroelasticityEquations< DIM >::force_solid_fct_pt().

◆ K_inv_pt

template<unsigned DIM>
double* oomph::PoroelasticityEquations< DIM >::K_inv_pt
private

1/k

Referenced by oomph::PoroelasticityEquations< DIM >::k_inv(), and oomph::PoroelasticityEquations< DIM >::k_inv_pt().

◆ Lambda_sq_pt

template<unsigned DIM>
double* oomph::PoroelasticityEquations< DIM >::Lambda_sq_pt
private

Timescale ratio (non-dim. density)

Referenced by oomph::PoroelasticityEquations< DIM >::lambda_sq(), and oomph::PoroelasticityEquations< DIM >::lambda_sq_pt().

◆ Mass_source_fct_pt

template<unsigned DIM>
MassSourceFctPt oomph::PoroelasticityEquations< DIM >::Mass_source_fct_pt
private

Pointer to the mass source function.

Referenced by oomph::PoroelasticityEquations< DIM >::mass_source(), and oomph::PoroelasticityEquations< DIM >::mass_source_fct_pt().

◆ Porosity_pt

template<unsigned DIM>
double* oomph::PoroelasticityEquations< DIM >::Porosity_pt
private

Porosity.

Referenced by oomph::PoroelasticityEquations< DIM >::porosity(), and oomph::PoroelasticityEquations< DIM >::porosity_pt().

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