Public Member Functions | Private Member Functions | Private Attributes | List of all members
oomph::SolidICProblem Class Reference

#include <elastic_problems.h>

+ Inheritance diagram for oomph::SolidICProblem:

Public Member Functions

 SolidICProblem ()
 
 SolidICProblem (const SolidICProblem &)=delete
 Broken copy constructor. More...
 
void operator= (const SolidICProblem &)=delete
 Broken assignment operator. More...
 
void actions_after_newton_solve ()
 Update after solve (empty) More...
 
void actions_before_newton_solve ()
 Update the problem specs before solve. (empty) More...
 
void set_static_initial_condition (Problem *problem_pt, Mesh *mesh_pt, SolidInitialCondition *ic_pt, const double &time)
 
void set_static_initial_condition (Problem *problem_pt, Mesh *mesh_pt, SolidInitialCondition *ic_pt)
 
template<class TIMESTEPPER >
void set_newmark_initial_condition_directly (Problem *problem_pt, Mesh *mesh_pt, TIMESTEPPER *timestepper_pt, SolidInitialCondition *ic_pt, const double &dt)
 
template<class TIMESTEPPER >
void set_newmark_initial_condition_consistently (Problem *problem_pt, Mesh *mesh_pt, TIMESTEPPER *timestepper_pt, SolidInitialCondition *ic_pt, const double &dt, SolidFiniteElement::MultiplierFctPt multiplier_fct_pt=0)
 
double & max_residual_after_consistent_newton_ic ()
 
- Public Member Functions inherited from oomph::Problem
virtual void debug_hook_fct (const unsigned &i)
 
void set_analytic_dparameter (double *const &parameter_pt)
 
void unset_analytic_dparameter (double *const &parameter_pt)
 
bool is_dparameter_calculated_analytically (double *const &parameter_pt)
 
void set_analytic_hessian_products ()
 
void unset_analytic_hessian_products ()
 
bool are_hessian_products_calculated_analytically ()
 
void set_pinned_values_to_zero ()
 
bool distributed () const
 
virtual void actions_before_adapt ()
 
virtual void actions_after_adapt ()
 Actions that are to be performed after a mesh adaptation. More...
 
OomphCommunicator * communicator_pt ()
 access function to the oomph-lib communicator More...
 
const OomphCommunicator * communicator_pt () const
 access function to the oomph-lib communicator, const version More...
 
 Problem ()
 
 Problem (const Problem &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const Problem &)=delete
 Broken assignment operator. More...
 
virtual ~Problem ()
 Virtual destructor to clean up memory. More...
 
Mesh *& mesh_pt ()
 Return a pointer to the global mesh. More...
 
Mesh *const & mesh_pt () const
 Return a pointer to the global mesh (const version) More...
 
Mesh *& mesh_pt (const unsigned &imesh)
 
Mesh *const & mesh_pt (const unsigned &imesh) const
 Return a pointer to the i-th submesh (const version) More...
 
unsigned nsub_mesh () const
 Return number of submeshes. More...
 
unsigned add_sub_mesh (Mesh *const &mesh_pt)
 
void flush_sub_meshes ()
 
void build_global_mesh ()
 
void rebuild_global_mesh ()
 
LinearSolver *& linear_solver_pt ()
 Return a pointer to the linear solver object. More...
 
LinearSolver *const & linear_solver_pt () const
 Return a pointer to the linear solver object (const version) More...
 
LinearSolver *& mass_matrix_solver_for_explicit_timestepper_pt ()
 
LinearSolver * mass_matrix_solver_for_explicit_timestepper_pt () const
 
EigenSolver *& eigen_solver_pt ()
 Return a pointer to the eigen solver object. More...
 
EigenSolver *const & eigen_solver_pt () const
 Return a pointer to the eigen solver object (const version) More...
 
Time *& time_pt ()
 Return a pointer to the global time object. More...
 
Time * time_pt () const
 Return a pointer to the global time object (const version). More...
 
double & time ()
 Return the current value of continuous time. More...
 
double time () const
 Return the current value of continuous time (const version) More...
 
TimeStepper *& time_stepper_pt ()
 
const TimeStepper * time_stepper_pt () const
 
TimeStepper *& time_stepper_pt (const unsigned &i)
 Return a pointer to the i-th timestepper. More...
 
ExplicitTimeStepper *& explicit_time_stepper_pt ()
 Return a pointer to the explicit timestepper. More...
 
unsigned long set_timestepper_for_all_data (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data=false)
 
virtual void shift_time_values ()
 Shift all values along to prepare for next timestep. More...
 
AssemblyHandler *& assembly_handler_pt ()
 Return a pointer to the assembly handler object. More...
 
AssemblyHandler *const & assembly_handler_pt () const
 Return a pointer to the assembly handler object (const version) More...
 
double & minimum_dt ()
 Access function to min timestep in adaptive timestepping. More...
 
double & maximum_dt ()
 Access function to max timestep in adaptive timestepping. More...
 
unsigned & max_newton_iterations ()
 Access function to max Newton iterations before giving up. More...
 
void problem_is_nonlinear (const bool &prob_lin)
 Access function to Problem_is_nonlinear. More...
 
double & max_residuals ()
 
bool & time_adaptive_newton_crash_on_solve_fail ()
 Access function for Time_adaptive_newton_crash_on_solve_fail. More...
 
double & newton_solver_tolerance ()
 
void add_time_stepper_pt (TimeStepper *const &time_stepper_pt)
 
void set_explicit_time_stepper_pt (ExplicitTimeStepper *const &explicit_time_stepper_pt)
 
void initialise_dt (const double &dt)
 
void initialise_dt (const Vector< double > &dt)
 
Data *& global_data_pt (const unsigned &i)
 Return a pointer to the the i-th global data object. More...
 
void add_global_data (Data *const &global_data_pt)
 
void flush_global_data ()
 
LinearAlgebraDistribution *const & dof_distribution_pt () const
 Return the pointer to the dof distribution (read-only) More...
 
unsigned long ndof () const
 Return the number of dofs. More...
 
unsigned ntime_stepper () const
 Return the number of time steppers. More...
 
unsigned nglobal_data () const
 Return the number of global data values. More...
 
unsigned self_test ()
 Self-test: Check meshes and global data. Return 0 for OK. More...
 
void enable_store_local_dof_pt_in_elements ()
 
void disable_store_local_dof_pt_in_elements ()
 
unsigned long assign_eqn_numbers (const bool &assign_local_eqn_numbers=true)
 
void describe_dofs (std::ostream &out= *(oomph_info.stream_pt())) const
 
void enable_discontinuous_formulation ()
 
void disable_discontinuous_formulation ()
 
void get_dofs (DoubleVector &dofs) const
 
void get_dofs (const unsigned &t, DoubleVector &dofs) const
 Return vector of the t'th history value of all dofs. More...
 
void set_dofs (const DoubleVector &dofs)
 Set the values of the dofs. More...
 
void set_dofs (const unsigned &t, DoubleVector &dofs)
 Set the history values of the dofs. More...
 
void set_dofs (const unsigned &t, Vector< double * > &dof_pt)
 
void add_to_dofs (const double &lambda, const DoubleVector &increment_dofs)
 Add lambda x incremenet_dofs[l] to the l-th dof. More...
 
double * global_dof_pt (const unsigned &i)
 
double & dof (const unsigned &i)
 i-th dof in the problem More...
 
double dof (const unsigned &i) const
 i-th dof in the problem (const version) More...
 
double *& dof_pt (const unsigned &i)
 Pointer to i-th dof in the problem. More...
 
double * dof_pt (const unsigned &i) const
 Pointer to i-th dof in the problem (const version) More...
 
virtual void get_inverse_mass_matrix_times_residuals (DoubleVector &Mres)
 
virtual void get_dvaluesdt (DoubleVector &f)
 
virtual void get_residuals (DoubleVector &residuals)
 Get the total residuals Vector for the problem. More...
 
virtual void get_jacobian (DoubleVector &residuals, DenseDoubleMatrix &jacobian)
 
virtual void get_jacobian (DoubleVector &residuals, CRDoubleMatrix &jacobian)
 
virtual void get_jacobian (DoubleVector &residuals, CCDoubleMatrix &jacobian)
 
virtual void get_jacobian (DoubleVector &residuals, SumOfMatrices &jacobian)
 
void get_fd_jacobian (DoubleVector &residuals, DenseMatrix< double > &jacobian)
 Get the full Jacobian by finite differencing. More...
 
void get_derivative_wrt_global_parameter (double *const &parameter_pt, DoubleVector &result)
 
void get_hessian_vector_products (DoubleVectorWithHaloEntries const &Y, Vector< DoubleVectorWithHaloEntries > const &C, Vector< DoubleVectorWithHaloEntries > &product)
 
void solve_eigenproblem (const unsigned &n_eval, Vector< std::complex< double >> &eigenvalue, Vector< DoubleVector > &eigenvector, const bool &steady=true)
 Solve the eigenproblem. More...
 
void solve_eigenproblem (const unsigned &n_eval, Vector< std::complex< double >> &eigenvalue, const bool &steady=true)
 
virtual void get_eigenproblem_matrices (CRDoubleMatrix &mass_matrix, CRDoubleMatrix &main_matrix, const double &shift=0.0)
 
void assign_eigenvector_to_dofs (DoubleVector &eigenvector)
 Assign the eigenvector passed to the function to the dofs. More...
 
void add_eigenvector_to_dofs (const double &epsilon, const DoubleVector &eigenvector)
 
void store_current_dof_values ()
 Store the current values of the degrees of freedom. More...
 
void restore_dof_values ()
 Restore the stored values of the degrees of freedom. More...
 
void enable_jacobian_reuse ()
 
void disable_jacobian_reuse ()
 Disable recycling of Jacobian in Newton iteration. More...
 
bool jacobian_reuse_is_enabled ()
 Is recycling of Jacobian in Newton iteration enabled? More...
 
bool & use_predictor_values_as_initial_guess ()
 
void newton_solve ()
 Use Newton method to solve the problem. More...
 
void enable_globally_convergent_newton_method ()
 enable globally convergent Newton method More...
 
void disable_globally_convergent_newton_method ()
 disable globally convergent Newton method More...
 
void newton_solve (unsigned const &max_adapt)
 
void steady_newton_solve (unsigned const &max_adapt=0)
 
void copy (Problem *orig_problem_pt)
 
virtual Problem * make_copy ()
 
virtual void read (std::ifstream &restart_file, bool &unsteady_restart)
 
virtual void read (std::ifstream &restart_file)
 
virtual void dump (std::ofstream &dump_file) const
 
void dump (const std::string &dump_file_name) const
 
void delete_all_external_storage ()
 
virtual void symmetrise_eigenfunction_for_adaptive_pitchfork_tracking ()
 
double * bifurcation_parameter_pt () const
 
void get_bifurcation_eigenfunction (Vector< DoubleVector > &eigenfunction)
 
void activate_fold_tracking (double *const &parameter_pt, const bool &block_solve=true)
 
void activate_bifurcation_tracking (double *const &parameter_pt, const DoubleVector &eigenvector, const bool &block_solve=true)
 
void activate_bifurcation_tracking (double *const &parameter_pt, const DoubleVector &eigenvector, const DoubleVector &normalisation, const bool &block_solve=true)
 
void activate_pitchfork_tracking (double *const &parameter_pt, const DoubleVector &symmetry_vector, const bool &block_solve=true)
 
void activate_hopf_tracking (double *const &parameter_pt, const bool &block_solve=true)
 
void activate_hopf_tracking (double *const &parameter_pt, const double &omega, const DoubleVector &null_real, const DoubleVector &null_imag, const bool &block_solve=true)
 
void deactivate_bifurcation_tracking ()
 
void reset_assembly_handler_to_default ()
 Reset the system to the standard non-augemented state. More...
 
double arc_length_step_solve (double *const &parameter_pt, const double &ds, const unsigned &max_adapt=0)
 
double arc_length_step_solve (Data *const &data_pt, const unsigned &data_index, const double &ds, const unsigned &max_adapt=0)
 
void reset_arc_length_parameters ()
 
int & sign_of_jacobian ()
 
void explicit_timestep (const double &dt, const bool &shift_values=true)
 Take an explicit timestep of size dt. More...
 
void unsteady_newton_solve (const double &dt)
 
void unsteady_newton_solve (const double &dt, const bool &shift_values)
 
void unsteady_newton_solve (const double &dt, const unsigned &max_adapt, const bool &first, const bool &shift=true)
 
double doubly_adaptive_unsteady_newton_solve (const double &dt, const double &epsilon, const unsigned &max_adapt, const bool &first, const bool &shift=true)
 
double doubly_adaptive_unsteady_newton_solve (const double &dt, const double &epsilon, const unsigned &max_adapt, const unsigned &suppress_resolve_after_spatial_adapt_flag, const bool &first, const bool &shift=true)
 
double adaptive_unsteady_newton_solve (const double &dt_desired, const double &epsilon)
 
double adaptive_unsteady_newton_solve (const double &dt_desired, const double &epsilon, const bool &shift_values)
 
void assign_initial_values_impulsive ()
 
void assign_initial_values_impulsive (const double &dt)
 
void calculate_predictions ()
 Calculate predictions. More...
 
void enable_mass_matrix_reuse ()
 
void disable_mass_matrix_reuse ()
 
bool mass_matrix_reuse_is_enabled ()
 Return whether the mass matrix is being reused. More...
 
void refine_uniformly (const Vector< unsigned > &nrefine_for_mesh)
 
void refine_uniformly (const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
 
void refine_uniformly_and_prune (const Vector< unsigned > &nrefine_for_mesh)
 
void refine_uniformly_and_prune (const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
 
void refine_uniformly (DocInfo &doc_info)
 
void refine_uniformly_and_prune (DocInfo &doc_info)
 
void refine_uniformly ()
 
void refine_uniformly (const unsigned &i_mesh, DocInfo &doc_info)
 Do uniform refinement for submesh i_mesh with documentation. More...
 
void refine_uniformly (const unsigned &i_mesh)
 Do uniform refinement for submesh i_mesh without documentation. More...
 
void p_refine_uniformly (const Vector< unsigned > &nrefine_for_mesh)
 
void p_refine_uniformly (const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
 
void p_refine_uniformly_and_prune (const Vector< unsigned > &nrefine_for_mesh)
 
void p_refine_uniformly_and_prune (const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
 
void p_refine_uniformly (DocInfo &doc_info)
 
void p_refine_uniformly_and_prune (DocInfo &doc_info)
 
void p_refine_uniformly ()
 
void p_refine_uniformly (const unsigned &i_mesh, DocInfo &doc_info)
 Do uniform p-refinement for submesh i_mesh with documentation. More...
 
void p_refine_uniformly (const unsigned &i_mesh)
 Do uniform p-refinement for submesh i_mesh without documentation. More...
 
void refine_selected_elements (const Vector< unsigned > &elements_to_be_refined)
 
void refine_selected_elements (const Vector< RefineableElement * > &elements_to_be_refined_pt)
 
void refine_selected_elements (const unsigned &i_mesh, const Vector< unsigned > &elements_to_be_refined)
 
void refine_selected_elements (const unsigned &i_mesh, const Vector< RefineableElement * > &elements_to_be_refined_pt)
 
void refine_selected_elements (const Vector< Vector< unsigned >> &elements_to_be_refined)
 
void refine_selected_elements (const Vector< Vector< RefineableElement * >> &elements_to_be_refined_pt)
 
void p_refine_selected_elements (const Vector< unsigned > &elements_to_be_refined)
 
void p_refine_selected_elements (const Vector< PRefineableElement * > &elements_to_be_refined_pt)
 
void p_refine_selected_elements (const unsigned &i_mesh, const Vector< unsigned > &elements_to_be_refined)
 
void p_refine_selected_elements (const unsigned &i_mesh, const Vector< PRefineableElement * > &elements_to_be_refined_pt)
 
void p_refine_selected_elements (const Vector< Vector< unsigned >> &elements_to_be_refined)
 
void p_refine_selected_elements (const Vector< Vector< PRefineableElement * >> &elements_to_be_refined_pt)
 
unsigned unrefine_uniformly ()
 
unsigned unrefine_uniformly (const unsigned &i_mesh)
 
void p_unrefine_uniformly (DocInfo &doc_info)
 
void p_unrefine_uniformly (const unsigned &i_mesh, DocInfo &doc_info)
 Do uniform p-unrefinement for submesh i_mesh without documentation. More...
 
void adapt (unsigned &n_refined, unsigned &n_unrefined)
 
void adapt ()
 
void p_adapt (unsigned &n_refined, unsigned &n_unrefined)
 
void p_adapt ()
 
void adapt_based_on_error_estimates (unsigned &n_refined, unsigned &n_unrefined, Vector< Vector< double >> &elemental_error)
 
void adapt_based_on_error_estimates (Vector< Vector< double >> &elemental_error)
 
void get_all_error_estimates (Vector< Vector< double >> &elemental_error)
 
void doc_errors (DocInfo &doc_info)
 Get max and min error for all elements in submeshes. More...
 
void doc_errors ()
 Get max and min error for all elements in submeshes. More...
 
void enable_info_in_newton_solve ()
 
void disable_info_in_newton_solve ()
 Disable the output of information when in the newton solver. More...
 
- Public Member Functions inherited from oomph::ExplicitTimeSteppableObject
 ExplicitTimeSteppableObject ()
 Empty constructor. More...
 
 ExplicitTimeSteppableObject (const ExplicitTimeSteppableObject &)=delete
 Broken copy constructor. More...
 
void operator= (const ExplicitTimeSteppableObject &)=delete
 Broken assignment operator. More...
 
virtual ~ExplicitTimeSteppableObject ()
 Empty destructor. More...
 
virtual void actions_before_explicit_stage ()
 
virtual void actions_after_explicit_stage ()
 

Private Member Functions

void backup_original_state ()
 Backup original state of all data associated with mesh. More...
 
void reset_original_state ()
 Reset original state of all data associated with mesh. More...
 
void setup_problem ()
 

Private Attributes

SolidInitialCondition * IC_pt
 Pointer to initial condition object. More...
 
Vector< int > Backup_pinned
 Vector to store pinned status of all data. More...
 
Vector< Vector< Data * > > Backup_ext_data
 Vector of Vectors to store pointers to exernal data in the elements. More...
 
double Max_residual_after_consistent_newton_ic
 

Additional Inherited Members

- Public Types inherited from oomph::Problem
typedef void(* SpatialErrorEstimatorFctPt) (Mesh *&mesh_pt, Vector< double > &elemental_error)
 Function pointer for spatial error estimator. More...
 
typedef void(* SpatialErrorEstimatorWithDocFctPt) (Mesh *&mesh_pt, Vector< double > &elemental_error, DocInfo &doc_info)
 Function pointer for spatial error estimator with doc. More...
 
- Public Attributes inherited from oomph::Problem
bool Shut_up_in_newton_solve
 
- Static Public Attributes inherited from oomph::Problem
static bool Suppress_warning_about_actions_before_read_unstructured_meshes
 
- Protected Types inherited from oomph::Problem
enum  Assembly_method {
  Perform_assembly_using_vectors_of_pairs , Perform_assembly_using_two_vectors , Perform_assembly_using_maps , Perform_assembly_using_lists ,
  Perform_assembly_using_two_arrays
}
 Enumerated flags to determine which sparse assembly method is used. More...
 
- Protected Member Functions inherited from oomph::Problem
unsigned setup_element_count_per_dof ()
 
virtual void sparse_assemble_row_or_column_compressed (Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
 
virtual void actions_before_newton_convergence_check ()
 
virtual void actions_before_newton_step ()
 
virtual void actions_after_newton_step ()
 
virtual void actions_before_implicit_timestep ()
 
virtual void actions_after_implicit_timestep ()
 
virtual void actions_after_implicit_timestep_and_error_estimation ()
 
virtual void actions_before_explicit_timestep ()
 Actions that should be performed before each explicit time step. More...
 
virtual void actions_after_explicit_timestep ()
 Actions that should be performed after each explicit time step. More...
 
virtual void actions_before_read_unstructured_meshes ()
 
virtual void actions_after_read_unstructured_meshes ()
 
virtual void actions_after_change_in_global_parameter (double *const &parameter_pt)
 
virtual void actions_after_change_in_bifurcation_parameter ()
 
virtual void actions_after_parameter_increase (double *const &parameter_pt)
 
double & dof_derivative (const unsigned &i)
 
double & dof_current (const unsigned &i)
 
virtual void set_initial_condition ()
 
virtual double global_temporal_error_norm ()
 
unsigned newton_solve_continuation (double *const &parameter_pt)
 
unsigned newton_solve_continuation (double *const &parameter_pt, DoubleVector &z)
 
void calculate_continuation_derivatives (double *const &parameter_pt)
 
void calculate_continuation_derivatives (const DoubleVector &z)
 
void calculate_continuation_derivatives_fd (double *const &parameter_pt)
 
bool does_pointer_correspond_to_problem_data (double *const &parameter_pt)
 
void set_consistent_pinned_values_for_continuation ()
 
- Protected Attributes inherited from oomph::Problem
Vector< Problem * > Copy_of_problem_pt
 
std::map< double *, bool > Calculate_dparameter_analytic
 
bool Calculate_hessian_products_analytic
 
LinearAlgebraDistribution * Dof_distribution_pt
 
Vector< double * > Dof_pt
 Vector of pointers to dofs. More...
 
DoubleVectorWithHaloEntries Element_count_per_dof
 
double Relaxation_factor
 
double Newton_solver_tolerance
 
unsigned Max_newton_iterations
 Maximum number of Newton iterations. More...
 
unsigned Nnewton_iter_taken
 
Vector< double > Max_res
 Maximum residuals at start and after each newton iteration. More...
 
double Max_residuals
 
bool Time_adaptive_newton_crash_on_solve_fail
 
bool Jacobian_reuse_is_enabled
 Is re-use of Jacobian in Newton iteration enabled? Default: false. More...
 
bool Jacobian_has_been_computed
 
bool Problem_is_nonlinear
 
bool Pause_at_end_of_sparse_assembly
 
bool Doc_time_in_distribute
 
unsigned Sparse_assembly_method
 
unsigned Sparse_assemble_with_arrays_initial_allocation
 
unsigned Sparse_assemble_with_arrays_allocation_increment
 
Vector< Vector< unsigned > > Sparse_assemble_with_arrays_previous_allocation
 
double Numerical_zero_for_sparse_assembly
 
double FD_step_used_in_get_hessian_vector_products
 
bool Mass_matrix_reuse_is_enabled
 
bool Mass_matrix_has_been_computed
 
bool Discontinuous_element_formulation
 
double Minimum_dt
 Minimum desired dt: if dt falls below this value, exit. More...
 
double Maximum_dt
 Maximum desired dt. More...
 
double DTSF_max_increase
 
double DTSF_min_decrease
 
double Minimum_dt_but_still_proceed
 
bool Scale_arc_length
 Boolean to control whether arc-length should be scaled. More...
 
double Desired_proportion_of_arc_length
 Proportion of the arc-length to taken by the parameter. More...
 
double Theta_squared
 
int Sign_of_jacobian
 Storage for the sign of the global Jacobian. More...
 
double Continuation_direction
 
double Parameter_derivative
 Storage for the derivative of the global parameter wrt arc-length. More...
 
double Parameter_current
 Storage for the present value of the global parameter. More...
 
bool Use_continuation_timestepper
 Boolean to control original or new storage of dof stuff. More...
 
unsigned Dof_derivative_offset
 
unsigned Dof_current_offset
 
Vector< double > Dof_derivative
 Storage for the derivative of the problem variables wrt arc-length. More...
 
Vector< double > Dof_current
 Storage for the present values of the variables. More...
 
double Ds_current
 Storage for the current step value. More...
 
unsigned Desired_newton_iterations_ds
 
double Minimum_ds
 Minimum desired value of arc-length. More...
 
bool Bifurcation_detection
 Boolean to control bifurcation detection via determinant of Jacobian. More...
 
bool Bisect_to_find_bifurcation
 Boolean to control wheter bisection is used to located bifurcation. More...
 
bool First_jacobian_sign_change
 Boolean to indicate whether a sign change has occured in the Jacobian. More...
 
bool Arc_length_step_taken
 Boolean to indicate whether an arc-length step has been taken. More...
 
bool Use_finite_differences_for_continuation_derivatives
 
OomphCommunicator * Communicator_pt
 The communicator for this problem. More...
 
bool Always_take_one_newton_step
 
double Timestep_reduction_factor_after_nonconvergence
 
bool Keep_temporal_error_below_tolerance
 
- Static Protected Attributes inherited from oomph::Problem
static ContinuationStorageScheme Continuation_time_stepper
 Storage for the single static continuation timestorage object. More...
 

Detailed Description

///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// IC problem for an elastic body discretised on a given (sub)-mesh. We switch the elements' residuals and Jacobians to the system of equations that forces the wall shape to become that of a specified "initial condition object".

Boundary conditions for all data items associated with the mesh's elements' are temporarily overwritten so that all positional dofs (and only those!) become free while all other data (nodal data and the element's internal and external data) are either pinned (for nodal and internal data) or completely removed (for pointers to external data). The complete removal of the (pointers to the) external data is necessary because in FSI problems, positional data of certain elements can feature as external data of others. Hence pinning them in their incarnation as external data would also pin them in their incarnation as positional data.

All data and its pinned-status are restored at the end of the call to setup_ic().

Constructor & Destructor Documentation

◆ SolidICProblem() [1/2]

oomph::SolidICProblem::SolidICProblem ( )
inline

Constructor. Initialise pointer to IC object to NULL. Create a dummy mesh that can be deleted when static problem finally goes out of scope at end of program execution.

Create dummy mesh

94  : IC_pt(0)
95  {
97  mesh_pt() = new DummyMesh;
98 
99  // Default value for checking of consistent assignment of Newmark IC
100  Max_residual_after_consistent_newton_ic = 1.0e-12;
101  }
oomph::Problem::mesh_pt
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
oomph::SolidICProblem::Max_residual_after_consistent_newton_ic
double Max_residual_after_consistent_newton_ic
Definition: elastic_problems.h:201
oomph::SolidICProblem::IC_pt
SolidInitialCondition * IC_pt
Pointer to initial condition object.
Definition: elastic_problems.h:190

References Max_residual_after_consistent_newton_ic, and oomph::Problem::mesh_pt().

◆ SolidICProblem() [2/2]

oomph::SolidICProblem::SolidICProblem ( const SolidICProblem &  )
delete

Broken copy constructor.

Member Function Documentation

◆ actions_after_newton_solve()

void oomph::SolidICProblem::actions_after_newton_solve ( )
inlinevirtual

Update after solve (empty)

Reimplemented from oomph::Problem.

111 {}

◆ actions_before_newton_solve()

void oomph::SolidICProblem::actions_before_newton_solve ( )
inlinevirtual

Update the problem specs before solve. (empty)

Reimplemented from oomph::Problem.

114 {}

◆ backup_original_state()

void oomph::SolidICProblem::backup_original_state ( )
private

Backup original state of all data associated with mesh.

Backup pinned status of all data associated with the mesh. Also backup the (pointers to the) elements' external data.

192  {
193  // Find out how many nodes there are
194  unsigned long n_node = mesh_pt()->nnode();
195 
196  // Flush vector which holds backup of pinned status
197  Backup_pinned.clear();
198 
199  // Flush vector which holds vectors with backup of (pointers to) external
200  // data
201  Backup_ext_data.clear();
202 
203  // Loop over all the nodes
204  for (unsigned n = 0; n < n_node; n++)
205  {
206  // Cast to an elastic node
207  SolidNode* node_pt = dynamic_cast<SolidNode*>(mesh_pt()->node_pt(n));
208 
209 #ifdef PARANOID
210  if (node_pt == 0)
211  {
212  throw OomphLibError("Wasn't able to cast to SolidNode\n",
213  OOMPH_CURRENT_FUNCTION,
214  OOMPH_EXCEPTION_LOCATION);
215  }
216 #endif
217 
218  // Get spatial dimension of node
219  unsigned ndim = node_pt->ndim();
220 
221  // Find out how many positional dofs there are
222  SolidFiniteElement* elem_pt =
223  dynamic_cast<SolidFiniteElement*>(mesh_pt()->element_pt(0));
224 
225 #ifdef PARANOID
226  if (elem_pt == 0)
227  {
228  throw OomphLibError("Wasn't able to cast to SolidFiniteElement\n",
229  OOMPH_CURRENT_FUNCTION,
230  OOMPH_EXCEPTION_LOCATION);
231  }
232 #endif
233 
234  unsigned ntype = elem_pt->nnodal_position_type();
235 
236 
237  // Loop over coordinate directions
238  for (unsigned i = 0; i < ndim; i++)
239  {
240  // Loop over type of dof
241  for (unsigned k = 0; k < ntype; k++)
242  {
243  // Backup pinned status
244  Backup_pinned.push_back(node_pt->position_is_pinned(k, i));
245  }
246  }
247 
248  // Loop over nodal values
249  unsigned nval = node_pt->nvalue();
250  for (unsigned ival = 0; ival < nval; ival++)
251  {
252  // Backup pinned status
253  Backup_pinned.push_back(node_pt->is_pinned(ival));
254  }
255  }
256 
257  // Loop over the elements
258  unsigned Nelement = mesh_pt()->nelement();
259  Backup_ext_data.resize(Nelement);
260  for (unsigned i = 0; i < Nelement; i++)
261  {
262  // Cast to proper element type
263  SolidFiniteElement* elem_pt =
264  dynamic_cast<SolidFiniteElement*>(mesh_pt()->element_pt(i));
265 
266 #ifdef PARANOID
267  if (elem_pt == 0)
268  {
269  throw OomphLibError("Wasn't able to cast to SolidFiniteElement\n",
270  OOMPH_CURRENT_FUNCTION,
271  OOMPH_EXCEPTION_LOCATION);
272  }
273 #endif
274 
275  // Find out number of external data
276  unsigned next = elem_pt->nexternal_data();
277  Backup_ext_data[i].resize(next);
278 
279  // Loop over external data
280  for (unsigned j = 0; j < next; j++)
281  {
282  Data* data_pt = elem_pt->external_data_pt(j);
283 
284  // Backup the pointer to external data
285  Backup_ext_data[i][j] = data_pt;
286  }
287 
288  // Find out number of internal data
289  unsigned nint = elem_pt->ninternal_data();
290 
291  // Loop over internal data
292  for (unsigned j = 0; j < nint; j++)
293  {
294  Data* data_pt = elem_pt->internal_data_pt(j);
295 
296  // Loop over internal values
297  unsigned nval = data_pt->nvalue();
298  for (unsigned ival = 0; ival < nval; ival++)
299  {
300  // Backup pinned status
301  Backup_pinned.push_back(data_pt->is_pinned(ival));
302  }
303  }
304 
305 
306 #ifdef PARANOID
307  // If there is internal solid data, complain
308  if (elem_pt->has_internal_solid_data())
309  {
310  std::string error_message =
311  "Automatic assignment of initial conditions doesn't work yet\n";
312  error_message +=
313  "for elasticity elements with internal solid dofs (pressures)\n";
314 
315  throw OomphLibError(
316  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
317  }
318 #endif
319  // for(unsigned i=0;i<nint_solid;i++)
320  // {
321  // Data* data_pt=elem_pt->internal_solid_data_pt(i);
322 
323  // // Loop over values
324  // unsigned nval=data_pt->nvalue();
325  // for (unsigned ival=0;ival<nval;ival++)
326  // {
327  // // Backup pinned status
328  // Backup_pinned.push_back(data_pt->is_pinned(ival));
329  // }
330  // }
331  }
332 
333  // Record number of dofs whose status was pinned
334  // oomph_info << "Number of backed up values " << Backup_pinned.size() <<
335  // std::endl;
336  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
oomph::Mesh::element_pt
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
oomph::Mesh::nnode
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
oomph::Mesh::node_pt
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
oomph::Mesh::nelement
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
oomph::SolidICProblem::Backup_ext_data
Vector< Vector< Data * > > Backup_ext_data
Vector of Vectors to store pointers to exernal data in the elements.
Definition: elastic_problems.h:196
oomph::SolidICProblem::Backup_pinned
Vector< int > Backup_pinned
Vector to store pinned status of all data.
Definition: elastic_problems.h:193
k
char char char int int * k
Definition: level2_impl.h:374
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References Backup_ext_data, Backup_pinned, oomph::Mesh::element_pt(), oomph::GeneralisedElement::external_data_pt(), oomph::SolidFiniteElement::has_internal_solid_data(), i, oomph::GeneralisedElement::internal_data_pt(), oomph::Data::is_pinned(), j, k, oomph::Problem::mesh_pt(), n, oomph::Node::ndim(), oomph::Mesh::nelement(), oomph::GeneralisedElement::nexternal_data(), oomph::GeneralisedElement::ninternal_data(), oomph::FiniteElement::nnodal_position_type(), oomph::Mesh::nnode(), oomph::Mesh::node_pt(), oomph::Data::nvalue(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::SolidNode::position_is_pinned(), and oomph::Global_string_for_annotation::string().

Referenced by set_newmark_initial_condition_consistently(), set_newmark_initial_condition_directly(), and set_static_initial_condition().

◆ max_residual_after_consistent_newton_ic()

double& oomph::SolidICProblem::max_residual_after_consistent_newton_ic ( )
inline

Max. tolerated residual after application of consistent Newmark IC. Used to check if we have specified the correct timescale ratio (non-dim density).

173  {
174  return Max_residual_after_consistent_newton_ic;
175  }

References Max_residual_after_consistent_newton_ic.

◆ operator=()

void oomph::SolidICProblem::operator= ( const SolidICProblem &  )
delete

Broken assignment operator.

◆ reset_original_state()

void oomph::SolidICProblem::reset_original_state ( )
private

Reset original state of all data associated with mesh.

Reset pinned status of all data and re-instate the pointers to the elements' external data.

344  {
345  // Find out how many nodes there are
346  unsigned long n_node = mesh_pt()->nnode();
347 
348  // Initialise counter for backed up dofs
349  unsigned count = 0;
350 
351  // Loop over all the nodes
352  for (unsigned n = 0; n < n_node; n++)
353  {
354  // Cast to an elastic node
355  SolidNode* node_pt = dynamic_cast<SolidNode*>(mesh_pt()->node_pt(n));
356 
357 #ifdef PARANOID
358  if (node_pt == 0)
359  {
360  throw OomphLibError("Wasn't able to cast to SolidNode\n",
361  OOMPH_CURRENT_FUNCTION,
362  OOMPH_EXCEPTION_LOCATION);
363  }
364 #endif
365 
366  // Get spatial dimension of node
367  unsigned ndim = node_pt->ndim();
368 
369  // Find out how many positional dofs there are
370  SolidFiniteElement* elem_pt =
371  dynamic_cast<SolidFiniteElement*>(mesh_pt()->element_pt(0));
372 
373 #ifdef PARANOID
374  if (elem_pt == 0)
375  {
376  throw OomphLibError("Wasn't able to cast to SolidFiniteElement\n",
377  OOMPH_CURRENT_FUNCTION,
378  OOMPH_EXCEPTION_LOCATION);
379  }
380 #endif
381 
382  unsigned ntype = elem_pt->nnodal_position_type();
383 
384  // Loop over coordinate directions
385  for (unsigned i = 0; i < ndim; i++)
386  {
387  // Loop over type of dof
388  for (unsigned k = 0; k < ntype; k++)
389  {
390  // Reset pinned status (positional dofs were all unpinned)
391  if (Backup_pinned[count])
392  {
393  node_pt->pin_position(k, i);
394  }
395  count++;
396  }
397  }
398 
399  // Loop over nodal values
400  unsigned nval = node_pt->nvalue();
401  for (unsigned ival = 0; ival < nval; ival++)
402  {
403  // Reset pinned status (nodal values were all pinned)
404  if (Backup_pinned[count])
405  {
406  node_pt->unpin(ival);
407  }
408  }
409  }
410 
411 
412  // Loop over the elements
413  unsigned Nelement = mesh_pt()->nelement();
414  for (unsigned i = 0; i < Nelement; i++)
415  {
416  // Cast to proper element type
417  SolidFiniteElement* elem_pt =
418  dynamic_cast<SolidFiniteElement*>(mesh_pt()->element_pt(i));
419 
420 #ifdef PARANOID
421  if (elem_pt == 0)
422  {
423  throw OomphLibError("Wasn't able to cast to SolidFiniteElement\n",
424  OOMPH_CURRENT_FUNCTION,
425  OOMPH_EXCEPTION_LOCATION);
426  }
427 #endif
428 
429  // Switch back to normal Jacobian
430  dynamic_cast<SolidFiniteElement*>(elem_pt)
431  ->disable_solve_for_consistent_newmark_accel();
432 
433  // Switch off flag for setting initial condition
434  dynamic_cast<SolidFiniteElement*>(elem_pt)->solid_ic_pt() = 0;
435 
436  // IF it's an FSI wall element then turn on external stuff again
437  if (FSIWallElement* fsi_elem_pt = dynamic_cast<FSIWallElement*>(elem_pt))
438  {
439  fsi_elem_pt->include_external_load_data();
440  }
441 
442  // Find out number of external data
443  unsigned next = Backup_ext_data[i].size();
444 
445  // Loop over external data
446  for (unsigned j = 0; j < next; j++)
447  {
448  // Backed up external data
449  Data* data_pt = Backup_ext_data[i][j];
450 
451  // Add external data
452  elem_pt->add_external_data(data_pt);
453  }
454 
455  // Find out number of internal data
456  unsigned nint = elem_pt->ninternal_data();
457 
458  // Loop over internal data
459  for (unsigned j = 0; j < nint; j++)
460  {
461  Data* data_pt = elem_pt->internal_data_pt(j);
462 
463  // Loop over internal values
464  unsigned nval = data_pt->nvalue();
465  for (unsigned ival = 0; ival < nval; ival++)
466  {
467  // Restore pinned status (values were all pinned)
468  if (!Backup_pinned[count])
469  {
470  data_pt->unpin(ival);
471  }
472  }
473  }
474 
475 
476 #ifdef PARANOID
477  // If there is internal solid data, complain
478  if (elem_pt->has_internal_solid_data())
479  {
480  std::string error_message =
481  "Automatic assignment of initial conditions doesn't work yet\n";
482  error_message +=
483  "for elasticity elements with internal solid dofs (pressures)\n";
484 
485  throw OomphLibError(
486  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
487  }
488 #endif
489  // for(unsigned i=0;i<nint_solid;i++)
490  // {
491  // Data* data_pt=elem_pt->internal_solid_data_pt(i);
492 
493  // // Loop over values
494  // unsigned nval=data_pt->nvalue();
495  // for (unsigned ival=0;ival<nval;ival++)
496  // {
497  // // Restore pinned status (values were all pinned)
498  // if (!Backup_pinned[count])
499  // {
500  // data_pt->unpin(ival);
501  // }
502  // }
503  // }
504  }
505 
506  // Check number of recovered pinned values
507  // oomph_info << "Recovered pin values " << count << std::endl;
508 
509  // Flush vector which holds backup of pinned status
510  Backup_pinned.clear();
511 
512  // Flush vector which holds vectors with backup of (pointers to) external
513  // data
514  Backup_ext_data.clear();
515  }

References oomph::GeneralisedElement::add_external_data(), Backup_ext_data, Backup_pinned, oomph::Mesh::element_pt(), oomph::SolidFiniteElement::has_internal_solid_data(), i, oomph::GeneralisedElement::internal_data_pt(), j, k, oomph::Problem::mesh_pt(), n, oomph::Node::ndim(), oomph::Mesh::nelement(), oomph::GeneralisedElement::ninternal_data(), oomph::FiniteElement::nnodal_position_type(), oomph::Mesh::nnode(), oomph::Mesh::node_pt(), oomph::Data::nvalue(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::SolidNode::pin_position(), oomph::Global_string_for_annotation::string(), and oomph::Data::unpin().

Referenced by set_newmark_initial_condition_consistently(), set_newmark_initial_condition_directly(), and set_static_initial_condition().

◆ set_newmark_initial_condition_consistently()

template<class TIMESTEPPER >
void oomph::SolidICProblem::set_newmark_initial_condition_consistently ( Problem *  problem_pt,
Mesh *  wall_mesh_pt,
TIMESTEPPER *  timestepper_pt,
SolidInitialCondition *  ic_pt,
const double &  dt,
SolidFiniteElement::MultiplierFctPt  multiplier_fct_pt = 0 
)

Setup initial condition for time-integration with Newmark's method. Past displacements and velocities are assigned directly (consistent with the idea that a second order ODE can take ICs up to 1st order, while the history value for the previous acceleration is determined by the condition that the weak equation is satisfied at the initial time.) The multiplier function needs to specify the factor that multiplies the inertia terms – typically this is a constant, given by the ratio \( \Lambda^2 \) of the problem's intrinsic timescale to the time used to non-dimensionalise the equations. If the function (pointer) is not specified the multiplier is assumed to be equal to 1.0 – appropriate for a non-dimensionalisation based on the problem's intrinsic timescale.

Setup initial condition for time-integration with Newmark's method. Past displacements and velocities are assigned directly (consistent with the idea that a second order ODE can take ICs up to 1st order, while the history value for the previous acceleration is determined by the condition that that the weak equation is satisfied at the initial time. The multiplier function needs to specify the factor that multiplies the inertia terms – typically this is a constant, given by the ratio \( \Lambda^2 \) of the problem's intrinsic timescale to the time used to non-dimensionalise the equations. If the function (pointer) is not specified, the multiplier is assumed to be equal to 1.0 – appropriate for a non-dimensionalisation based on the problem's intrinsic timescale.

Pointer to member type solver

363  {
364 #ifdef PARANOID
365  if (timestepper_pt->type() != "Newmark")
366  {
367  std::ostringstream error_message;
368  error_message
369  << "SolidICProblem::set_newmark_initial_condition_consistently()\n"
370  << "can only be called for Newmark type timestepper whereas\n "
371  << "you've called it for " << timestepper_pt->type() << std::endl;
372 
373  throw OomphLibError(
374  error_message.str(),
375  "SolidICProblem::set_newmark_initial_condition_consistently()",
376  OOMPH_EXCEPTION_LOCATION);
377  }
378 #endif
379 
380  // Set value of dt
381  timestepper_pt->time_pt()->dt() = dt;
382 
383  // Set the weights
384  timestepper_pt->set_weights();
385 
386  // Delete dummy mesh
387  delete mesh_pt();
388 
389  // Set pointer to mesh
390  mesh_pt() = wall_mesh_pt;
391 
392  // Set pointer to initial condition object
393  IC_pt = ic_pt;
394 
395  // Backup the pinned status of all dofs and remove external data
396  // of all elements
397  backup_original_state();
398 
399  // Now alter the pinned status so that the IC problem for the
400  // positional variables can be solved; setup equation numbering
401  // scheme
402  setup_problem();
403 
404  // Number of history values
405  unsigned ntstorage =
406  IC_pt->geom_object_pt()->time_stepper_pt()->ntstorage();
407 
408  // Set values at previous time
409  //----------------------------
410 
411  // Loop over number of previous times stored
412  unsigned nprevtime =
413  IC_pt->geom_object_pt()->time_stepper_pt()->nprev_values();
414 
415  // Backup previous times:
416  Vector<double> prev_time(nprevtime + 1);
417  for (unsigned i = 0; i <= nprevtime; i++)
418  {
419  prev_time[i] =
420  IC_pt->geom_object_pt()->time_stepper_pt()->time_pt()->time(i);
421  }
422 
423  // Loop over previous times & set values themselves
424  //-------------------------------------------------
425  for (unsigned i = 1; i <= nprevtime; i++)
426  {
427  // Set time for geometric object that specifies initial condition
428  // [Note: this acts on time everywhere!]
429  IC_pt->geom_object_pt()->time_stepper_pt()->time_pt()->time() =
430  prev_time[i];
431 
432  // Set flag to ensure that the t_deriv-th time derivative
433  // of the prescribed solution gets stored in displacements
434  IC_pt->ic_time_deriv() = 0;
435 
436  // Solve the problem for initial shape: After this solve
437  // The nodes's current positions represent the position at
438  // previous time level i.
439  newton_solve();
440 
441  // Loop over all the nodes
442  unsigned n_node = mesh_pt()->nnode();
443  for (unsigned n = 0; n < n_node; n++)
444  {
445  // Get the variable position data
446  Data* position_data_pt = dynamic_cast<SolidNode*>(mesh_pt()->node_pt(n))
447  ->variable_position_pt();
448  // Get number of values
449  unsigned nval = position_data_pt->nvalue();
450 
451  // Assign values at previous times into their corresponding
452  // slots
453  for (unsigned ival = 0; ival < nval; ival++)
454  {
455  position_data_pt->set_value(
456  i, ival, position_data_pt->value(0, ival));
457  }
458  }
459  }
460 
461  // Set veloc (1st time deriv) at previous time and store in appropriate
462  //---------------------------------------------------------------------
463  // history value. Also assign zero value for 2nd time deriv. at
464  //-------------------------------------------------------------
465  // previous time
466  //--------------
467 
468  // Set time for geometric object that specifies initial condition
469  // to previous time [Note: this acts on time everywhere!]
470  IC_pt->geom_object_pt()->time_stepper_pt()->time_pt()->time() =
471  prev_time[1];
472 
473  // Set flag to ensure that the t_deriv-th time derivative
474  // of the prescribed solution gets stored in displacements
475  IC_pt->ic_time_deriv() = 1;
476 
477  // Solve the problem for initial shape: After this solve
478  // The nodes's current positions represent the time derivatives of at
479  // the positons at previous time level.
480  newton_solve();
481 
482  // Loop over all the nodes
483  unsigned n_node = mesh_pt()->nnode();
484  for (unsigned n = 0; n < n_node; n++)
485  {
486  // Get the position data
487  Data* position_data_pt =
488  dynamic_cast<SolidNode*>(mesh_pt()->node_pt(n))->variable_position_pt();
489 
490  // Get number of values
491  unsigned nval = position_data_pt->nvalue();
492 
493  // Assign values at previous times into their corresponding
494  // slots (last but one history value of Newmark scheme)
495  for (unsigned ival = 0; ival < nval; ival++)
496  {
497  position_data_pt->set_value(
498  ntstorage - 2, ival, position_data_pt->value(0, ival));
499  position_data_pt->set_value(ntstorage - 1, ival, 0.0);
500  }
501  }
502 
503 
504  // Set values at current time
505  //---------------------------
506 
507  // Reset time to current value
508  // [Note: this acts on time everywhere!]
509  IC_pt->geom_object_pt()->time_stepper_pt()->time_pt()->time() =
510  prev_time[0];
511 
512  // Set flag to ensure that the t_deriv-th time derivative
513  // of the prescribed solution gets stored in displacements
514  IC_pt->ic_time_deriv() = 0;
515 
516  // Solve the problem for initial shape
517  newton_solve();
518 
519 
520  // Now solve for the correction to the Newmark accelerations
521  //----------------------------------------------------------
522  // at previous time:
523  //------------------
524 
525  // Loop over the elements
526  unsigned Nelement = mesh_pt()->nelement();
527  for (unsigned i = 0; i < Nelement; i++)
528  {
529  // Cast to proper element type
530  SolidFiniteElement* elem_pt =
531  dynamic_cast<SolidFiniteElement*>(mesh_pt()->element_pt(i));
532 
533  // Switch system to the one that determines the Newmark accelerations
534  // by setting the Jacobian to the mass matrix
535  elem_pt->enable_solve_for_consistent_newmark_accel();
536 
537  // Set pointer to multiplier function
538  elem_pt->multiplier_fct_pt() = multiplier_fct_pt;
539 
540  // Switch off pointer to initial condition object
541  elem_pt->solid_ic_pt() = 0;
542  }
543 
544  // Correction vector
545  DoubleVector correction;
546 
548  typedef void (LinearSolver::*SolverMemPtr)(Problem* const& problem,
549  DoubleVector& result);
550  SolverMemPtr solver_mem_pt = &LinearSolver::solve;
551 
552  // Now do the linear solve
553  (linear_solver_pt()->*solver_mem_pt)(this, correction);
554 
555  // Update discrete 2nd deriv at previous time so that it's consistent
556  // with PDE at current time
557 
558  // Loop over all the nodes
559  for (unsigned n = 0; n < n_node; n++)
560  {
561  // Get the pointer to the position data
562  Data* position_data_pt =
563  dynamic_cast<SolidNode*>(mesh_pt()->node_pt(n))->variable_position_pt();
564 
565  // Get number of values
566  unsigned nval = position_data_pt->nvalue();
567 
568  // Assign values for the history value that corresponds to the
569  // previous accel in Newmark scheme so that the PDE is satsified
570  // at current time
571  for (unsigned ival = 0; ival < nval; ival++)
572  {
573  // Global equation number
574  int ieqn = position_data_pt->eqn_number(ival);
575 
576 #ifdef PARANOID
577  if (ieqn < 0)
578  {
579  throw OomphLibError(
580  "No positional dofs should be pinned at this stage!",
581  OOMPH_CURRENT_FUNCTION,
582  OOMPH_EXCEPTION_LOCATION);
583  }
584 #endif
585 
586  // Update the value
587  *(position_data_pt->value_pt(ntstorage - 1, ival)) -= correction[ieqn];
588  }
589  }
590 
591 
592 #ifdef PARANOID
593  // Check the residual
594  DoubleVector residuals;
595  get_residuals(residuals);
596  double res_max = residuals.max();
597  oomph_info
598  << "Max. residual after assigning consistent initial conditions: "
599  << res_max << std::endl;
600  if (res_max > Max_residual_after_consistent_newton_ic)
601  {
602  std::ostringstream error_message;
603  error_message << "Residual is bigger than allowed! [Current tolerance: "
604  << Max_residual_after_consistent_newton_ic << "]\n\n";
605  error_message << "This is probably because you've not specified the "
606  << "correct multiplier \n(the product of growth factor "
607  << "and timescale ratio [the non-dim density]). \nPlease "
608  << "check the Solid Mechanics Theory Tutorial for "
609  << "details. \n\n"
610  << "If you're sure that the residual is OK, overwrite "
611  << "the default tolerance using\n";
612  error_message
613  << "SolidICProblem::max_residual_after_consistent_newton_ic()"
614  << std::endl
615  << "or recompile without the PARANOID flag." << std::endl;
616 
617  throw OomphLibError(
618  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
619  }
620 #endif
621 
622  // Reset the pinned status and re-attach the external data to the elements
623  reset_original_state();
624 
625  // Set pointer to dummy mesh so there's something that can be deleted
626  // when static problem finally goes out of scope.
627  mesh_pt() = new DummyMesh;
628 
629  // We have temporarily over-written equation numbers -- need
630  // to reset them now
631  oomph_info << "Number of equations in big problem: "
632  << problem_pt->assign_eqn_numbers() << std::endl;
633  }
DoubleVector
std::vector< double > DoubleVector
loads clump configuration
Definition: ClumpInput.h:26
oomph::GeomObject::time_stepper_pt
TimeStepper *& time_stepper_pt()
Definition: geom_objects.h:192
oomph::LinearSolver::solve
virtual void solve(Problem *const &problem_pt, DoubleVector &result)=0
oomph::Problem::get_residuals
virtual void get_residuals(DoubleVector &residuals)
Get the total residuals Vector for the problem.
Definition: problem.cc:3714
oomph::Problem::linear_solver_pt
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1466
oomph::Problem::newton_solve
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8783
oomph::Problem::Problem
Problem()
Definition: problem.cc:69
oomph::SolidICProblem::backup_original_state
void backup_original_state()
Backup original state of all data associated with mesh.
Definition: elastic_problems.cc:191
oomph::SolidICProblem::setup_problem
void setup_problem()
Definition: elastic_problems.cc:48
oomph::SolidICProblem::reset_original_state
void reset_original_state()
Reset original state of all data associated with mesh.
Definition: elastic_problems.cc:343
oomph::SolidInitialCondition::geom_object_pt
GeomObject *& geom_object_pt()
Definition: elements.h:3511
oomph::SolidInitialCondition::ic_time_deriv
unsigned & ic_time_deriv()
Which time derivative are we currently assigning?
Definition: elements.h:3517
oomph::TimeStepper::nprev_values
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
oomph::TimeStepper::ntstorage
unsigned ntstorage() const
Definition: timesteppers.h:601
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
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319
problem
Constructor for SteadyAxisymAdvectionDiffusion problem
Definition: steady_axisym_advection_diffusion.cc:213

References oomph::Problem::assign_eqn_numbers(), backup_original_state(), oomph::Mesh::element_pt(), oomph::SolidFiniteElement::enable_solve_for_consistent_newmark_accel(), oomph::Data::eqn_number(), oomph::SolidInitialCondition::geom_object_pt(), oomph::Problem::get_residuals(), i, IC_pt, oomph::SolidInitialCondition::ic_time_deriv(), oomph::Problem::linear_solver_pt(), oomph::DoubleVector::max(), Max_residual_after_consistent_newton_ic, oomph::Problem::mesh_pt(), oomph::SolidFiniteElement::multiplier_fct_pt(), n, oomph::Mesh::nelement(), oomph::Problem::newton_solve(), oomph::Mesh::nnode(), oomph::Mesh::node_pt(), oomph::TimeStepper::nprev_values(), oomph::TimeStepper::ntstorage(), oomph::Data::nvalue(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, problem, reset_original_state(), oomph::Data::set_value(), setup_problem(), oomph::SolidFiniteElement::solid_ic_pt(), oomph::LinearSolver::solve(), oomph::Time::time(), oomph::TimeStepper::time_pt(), oomph::GeomObject::time_stepper_pt(), oomph::Data::value(), and oomph::Data::value_pt().

◆ set_newmark_initial_condition_directly()

template<class TIMESTEPPER >
void oomph::SolidICProblem::set_newmark_initial_condition_directly ( Problem *  problem_pt,
Mesh *  wall_mesh_pt,
TIMESTEPPER *  timestepper_pt,
SolidInitialCondition *  ic_pt,
const double &  dt 
)

Setup initial condition for time-integration with Newmark's method. History values are assigned to that the velocity and accelerations determined by the Newmark scheme are exact at the initial time.

218  {
219 #ifdef PARANOID
220  if (timestepper_pt->type() != "Newmark")
221  {
222  std::ostringstream error_message;
223  error_message
224  << "SolidICProblem::set_newmark_initial_condition_directly()\n"
225  << "can only be called for Newmark type timestepper whereas\n "
226  << "you've called it for " << timestepper_pt->type() << std::endl;
227 
228  throw OomphLibError(
229  error_message.str(),
230  "SolidICProblem::set_newmark_initial_condition_directly()",
231  OOMPH_EXCEPTION_LOCATION);
232  }
233 #endif
234 
235  // Set value of dt
236  timestepper_pt->time_pt()->dt() = dt;
237 
238  // Set the weights
239  timestepper_pt->set_weights();
240 
241  // Delete dummy mesh
242  delete mesh_pt();
243 
244  // Set pointer to mesh
245  mesh_pt() = wall_mesh_pt;
246 
247  // Set pointer to initial condition object
248  IC_pt = ic_pt;
249 
250  // Backup the pinned status of all dofs and remove external data
251  // of all elements
252  backup_original_state();
253 
254  // Now alter the pinned status so that the IC problem for the
255  // positional variables can be solved; setup equation numbering
256  // scheme
257  setup_problem();
258 
259 
260  // Store times at which we need to assign ic:
261  double current_time = timestepper_pt->time_pt()->time();
262  double previous_time = timestepper_pt->time_pt()->time(1);
263 
264  // Stage 1: Set values and time derivs at current time
265  //----------------------------------------------------
266 
267  // [Note: this acts on time everywhere!]
268  IC_pt->geom_object_pt()->time_stepper_pt()->time_pt()->time() =
269  current_time;
270 
271  // Loop over time-derivatives
272  for (unsigned t_deriv = 0; t_deriv <= 2; t_deriv++)
273  {
274  // Set flag to ensure that the t_deriv-th time derivative
275  // of the prescribed solution gets stored in displacements
276  IC_pt->ic_time_deriv() = t_deriv;
277 
278  // Solve the problem for initial shape
279  newton_solve();
280 
281  // Loop over all the nodes
282  unsigned n_node = mesh_pt()->nnode();
283  for (unsigned n = 0; n < n_node; n++)
284  {
285  // Assign current time derivative in its temporary storage
286  // position
287  timestepper_pt->assign_initial_data_values_stage1(
288  t_deriv,
289  dynamic_cast<SolidNode*>(mesh_pt()->node_pt(n))
290  ->variable_position_pt());
291  }
292  }
293 
294 
295  // Stage 2: Now get position at previous time and adjust previous
296  //---------------------------------------------------------------
297  // values of veloc and accel in Newmark scheme so that current
298  //---------------------------------------------------------------
299  // veloc and accel (specified in step 1) are represented exactly.
300  //---------------------------------------------------------------
301 
302  // [Note: this acts on time everywhere and needs to be reset!]
303  IC_pt->geom_object_pt()->time_stepper_pt()->time_pt()->time() =
304  previous_time;
305 
306  // Set flag to ensure that the t_deriv-th time derivative
307  // of the prescribed solution gets stored in displacements
308  IC_pt->ic_time_deriv() = 0;
309 
310  // Solve the problem for initial shape
311  newton_solve();
312 
313  // Loop over all the nodes and make the final adjustments
314  unsigned n_node = mesh_pt()->nnode();
315  for (unsigned n = 0; n < n_node; n++)
316  {
317  timestepper_pt->assign_initial_data_values_stage2(
318  dynamic_cast<SolidNode*>(mesh_pt()->node_pt(n))
319  ->variable_position_pt());
320  }
321 
322  // Reset time
323  IC_pt->geom_object_pt()->time_stepper_pt()->time_pt()->time() =
324  current_time;
325 
326  // Reset the pinned status and re-attach the external data to the elements
327  reset_original_state();
328 
329  // Set pointer to dummy mesh so there's something that can be deleted
330  // when static problem finally goes out of scope.
331  mesh_pt() = new DummyMesh;
332 
333  // We have temporarily over-written equation numbers -- need
334  // to reset them now
335  oomph_info << "Number of equations in big problem: "
336  << problem_pt->assign_eqn_numbers() << std::endl;
337  }

References oomph::Problem::assign_eqn_numbers(), backup_original_state(), oomph::SolidInitialCondition::geom_object_pt(), IC_pt, oomph::SolidInitialCondition::ic_time_deriv(), oomph::Problem::mesh_pt(), n, oomph::Problem::newton_solve(), oomph::Mesh::nnode(), OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, reset_original_state(), setup_problem(), oomph::Time::time(), oomph::TimeStepper::time_pt(), and oomph::GeomObject::time_stepper_pt().

◆ set_static_initial_condition() [1/2]

void oomph::SolidICProblem::set_static_initial_condition ( Problem *  problem_pt,
Mesh *  mesh_pt,
SolidInitialCondition *  ic_pt 
)
inline

Force the elastic structure that is discretised on the specified mesh to deform in the shape of the initial condition object (wrapper)

129  {
130  double time;
131  set_static_initial_condition(problem_pt, mesh_pt, ic_pt, time);
132  }
oomph::Problem::time
double & time()
Return the current value of continuous time.
Definition: problem.cc:11531
oomph::SolidICProblem::set_static_initial_condition
void set_static_initial_condition(Problem *problem_pt, Mesh *mesh_pt, SolidInitialCondition *ic_pt, const double &time)
Definition: elastic_problems.cc:522

References oomph::Problem::mesh_pt(), set_static_initial_condition(), and oomph::Problem::time().

◆ set_static_initial_condition() [2/2]

void oomph::SolidICProblem::set_static_initial_condition ( Problem *  problem_pt,
Mesh *  wall_mesh_pt,
SolidInitialCondition *  ic_pt,
const double &  time 
)

Force the elastic structure that is discretised on the specified mesh to deform in the shape of the initial condition object (evaluated at the time specified)

IC problem for wall: Deform wall into the static initial shape described by the IC object at given time.

527  {
528  // Tell this sub-problem it is distributed if the main problem is
529  // distributed
530 #ifdef OOMPH_HAS_MPI
531  if (problem_pt->problem_has_been_distributed())
532  {
533  Problem_has_been_distributed = true;
534  }
535  // This (sub-)problem needs to know about the oomph communicator
536  delete Communicator_pt;
537  Communicator_pt = new OomphCommunicator(problem_pt->communicator_pt());
538 #endif
539 
540 
541  // Backup value of time
542  double backup_time = 0.0;
543 
544  // Set value of time for IC object (needs to be backed up and
545  // restored since it might be pointed to by other objects)
546  TimeStepper* timestepper_pt = ic_pt->geom_object_pt()->time_stepper_pt();
547  if (timestepper_pt != 0)
548  {
549  backup_time = timestepper_pt->time_pt()->time();
550  timestepper_pt->time_pt()->time() = time;
551  }
552 
553  // Delete dummy mesh
554  delete mesh_pt();
555 
556  // Set pointer to mesh
557  mesh_pt() = wall_mesh_pt;
558 
559  // Set pointer to initial condition object
560  IC_pt = ic_pt;
561 
562  // Backup the pinned status of all dofs and remove external data
563  // of all elements
564  backup_original_state();
565 
566  // Now alter the pinned status so that the IC problem for the
567  // positional variables can be solved; setup equation numbering
568  // scheme
569  setup_problem();
570 
571  // Assign displacements
572  IC_pt->ic_time_deriv() = 0;
573 
574  // Solve the problem for initial shape
575  newton_solve();
576 
577  // Impulsive start
578  assign_initial_values_impulsive();
579 
580  // Reset the pinned status and re-attach the external data to the elements
581  reset_original_state();
582 
583  // Reset time
584  if (timestepper_pt != 0)
585  {
586  timestepper_pt->time_pt()->time() = backup_time;
587  }
588 
589  // Set pointer to dummy mesh so there's something that can be deleted
590  // when static problem finally goes out of scope.
591  mesh_pt() = new DummyMesh;
592 
593  // We have temporarily over-written equation numbers -- need
594  // to reset them now
595  oomph_info << "Number of equations in big problem: "
596  << problem_pt->assign_eqn_numbers() << std::endl;
597  }
oomph::Problem::assign_initial_values_impulsive
void assign_initial_values_impulsive()
Definition: problem.cc:11499
oomph::Problem::Communicator_pt
OomphCommunicator * Communicator_pt
The communicator for this problem.
Definition: problem.h:1242

References oomph::Problem::assign_eqn_numbers(), oomph::Problem::assign_initial_values_impulsive(), backup_original_state(), oomph::Problem::Communicator_pt, oomph::Problem::communicator_pt(), oomph::SolidInitialCondition::geom_object_pt(), IC_pt, oomph::SolidInitialCondition::ic_time_deriv(), oomph::Problem::mesh_pt(), oomph::Problem::newton_solve(), oomph::oomph_info, reset_original_state(), setup_problem(), oomph::Problem::time(), oomph::Time::time(), oomph::TimeStepper::time_pt(), and oomph::GeomObject::time_stepper_pt().

Referenced by set_static_initial_condition().

◆ setup_problem()

void oomph::SolidICProblem::setup_problem ( )
private

Change pinned status of all data associated with mesh so that the IC problem can be solved.

/////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////// Setup IC problem by:

  • Pinning all nodal values in the mesh
  • Pinning the internal data of all elements.
  • Freeing/unpinnning all positional data.
  • Flushing the pointers to the elements' external data.
  • Setting the pointer to the IC object for all elements to ensure that the correct residuals/Jacobians are computed.
49  {
50  // Find out how many nodes there are
51  unsigned long n_node = mesh_pt()->nnode();
52 
53  // Loop over all the nodes
54  for (unsigned n = 0; n < n_node; n++)
55  {
56  // Cast to an elastic node
57  SolidNode* node_pt = dynamic_cast<SolidNode*>(mesh_pt()->node_pt(n));
58 
59 #ifdef PARANOID
60  if (node_pt == 0)
61  {
62  throw OomphLibError("Wasn't able to cast to SolidNode\n",
63  OOMPH_CURRENT_FUNCTION,
64  OOMPH_EXCEPTION_LOCATION);
65  }
66 #endif
67 
68  // Get spatial dimension of node
69  unsigned ndim = node_pt->ndim();
70 
71  // Find out how many positional dofs there are
72  SolidFiniteElement* elem_pt =
73  dynamic_cast<SolidFiniteElement*>(mesh_pt()->element_pt(0));
74 
75 #ifdef PARANOID
76  if (elem_pt == 0)
77  {
78  throw OomphLibError("Wasn't able to cast to SolidFiniteElement\n",
79  OOMPH_CURRENT_FUNCTION,
80  OOMPH_EXCEPTION_LOCATION);
81  }
82 #endif
83 
84  unsigned ntype = elem_pt->nnodal_position_type();
85 
86  // Loop over coordinate directions
87  for (unsigned i = 0; i < ndim; i++)
88  {
89  // Loop over type of dof
90  for (unsigned k = 0; k < ntype; k++)
91  {
92  // Unpin them
93  node_pt->unpin_position(k, i);
94  }
95  }
96 
97  // Loop over nodal values
98  unsigned nval = node_pt->nvalue();
99  for (unsigned ival = 0; ival < nval; ival++)
100  {
101  // Pin them
102  node_pt->pin(ival);
103  }
104  }
105 
106 
107  // Loop over the elements
108  unsigned Nelement = mesh_pt()->nelement();
109  for (unsigned i = 0; i < Nelement; i++)
110  {
111  // Cast to proper element type
112  SolidFiniteElement* elem_pt =
113  dynamic_cast<SolidFiniteElement*>(mesh_pt()->element_pt(i));
114 
115 #ifdef PARANOID
116  if (elem_pt == 0)
117  {
118  throw OomphLibError("Wasn't able to cast to SolidFiniteElement\n",
119  OOMPH_CURRENT_FUNCTION,
120  OOMPH_EXCEPTION_LOCATION);
121  }
122 #endif
123 
124 
125  // Set flag for setting initial condition
126  elem_pt->solid_ic_pt() = IC_pt;
127 
128  // We've backed up the element's external data: Flush it
129  elem_pt->flush_external_data();
130 
131  // IF it's an FSI wall element then kill external stuff
132  if (FSIWallElement* fsi_elem_pt = dynamic_cast<FSIWallElement*>(elem_pt))
133  {
134  fsi_elem_pt->exclude_external_load_data();
135  }
136 
137  // Find out number of internal data
138  unsigned nint = elem_pt->ninternal_data();
139 
140  // Loop over internal data
141  for (unsigned j = 0; j < nint; j++)
142  {
143  Data* data_pt = elem_pt->internal_data_pt(j);
144 
145  // Loop over internal values
146  unsigned nval = data_pt->nvalue();
147  for (unsigned ival = 0; ival < nval; ival++)
148  {
149  // Pin internal values
150  data_pt->pin(ival);
151  }
152  }
153 
154 #ifdef PARANOID
155  // Is there internal solid data
156  if (elem_pt->has_internal_solid_data())
157  {
158  std::string error_message =
159  "Automatic assignment of initial conditions doesn't work yet\n";
160  error_message +=
161  "for elasticity elements with internal solid dofs (pressures)\n";
162 
163  throw OomphLibError(
164  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
165  }
166 #endif
167  // for(unsigned i=0;i<nint_solid;i++)
168  // {
169  // Data* data_pt=elem_pt->internal_solid_data_pt(i);
170 
171  // // Loop over values
172  // unsigned nval=data_pt->nvalue();
173  // for (unsigned ival=0;ival<nval;ival++)
174  // {
175  // // Pin internal values
176  // data_pt->pin(ival);
177  // }
178  // }
179  }
180 
181  // Setup equation numbers for IC problem
182  oomph_info << "# of dofs for wall initial guess" << assign_eqn_numbers()
183  << std::endl;
184  }
oomph::Problem::assign_eqn_numbers
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Definition: problem.cc:1989

References oomph::Problem::assign_eqn_numbers(), oomph::Mesh::element_pt(), oomph::GeneralisedElement::flush_external_data(), oomph::SolidFiniteElement::has_internal_solid_data(), i, IC_pt, oomph::GeneralisedElement::internal_data_pt(), j, k, oomph::Problem::mesh_pt(), n, oomph::Node::ndim(), oomph::Mesh::nelement(), oomph::GeneralisedElement::ninternal_data(), oomph::FiniteElement::nnodal_position_type(), oomph::Mesh::nnode(), oomph::Mesh::node_pt(), oomph::Data::nvalue(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::Data::pin(), oomph::SolidFiniteElement::solid_ic_pt(), oomph::Global_string_for_annotation::string(), and oomph::SolidNode::unpin_position().

Referenced by set_newmark_initial_condition_consistently(), set_newmark_initial_condition_directly(), and set_static_initial_condition().

Member Data Documentation

◆ Backup_ext_data

Vector<Vector<Data*> > oomph::SolidICProblem::Backup_ext_data
private

Vector of Vectors to store pointers to exernal data in the elements.

Referenced by backup_original_state(), and reset_original_state().

◆ Backup_pinned

Vector<int> oomph::SolidICProblem::Backup_pinned
private

Vector to store pinned status of all data.

Referenced by backup_original_state(), and reset_original_state().

◆ IC_pt

SolidInitialCondition* oomph::SolidICProblem::IC_pt
private

Pointer to initial condition object.

Referenced by set_newmark_initial_condition_consistently(), set_newmark_initial_condition_directly(), set_static_initial_condition(), and setup_problem().

◆ Max_residual_after_consistent_newton_ic

double oomph::SolidICProblem::Max_residual_after_consistent_newton_ic
private

Max. tolerated residual after application of consistent Newmark IC. Used to check if we have specified the correct timescale ratio (non-dim density).

Referenced by max_residual_after_consistent_newton_ic(), set_newmark_initial_condition_consistently(), and SolidICProblem().

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