ContactProblem< ELEMENT > Class Template Reference

Problem class. More...

+ Inheritance diagram for ContactProblem< ELEMENT >:

Public Member Functions

 ContactProblem ()
 Constructor. More...
 
 ~ContactProblem ()
 Destructor (empty) More...
 
void actions_before_implicit_timestep ()
 Actions before timestep. More...
 
void actions_after_newton_solve ()
 Update the problem specs after solve (empty) More...
 
void actions_before_newton_solve ()
 Update the problem specs before solve (empty) More...
 
void actions_before_adapt ()
 Actions before adapt: wipe contact elements. More...
 
void actions_after_adapt ()
 
void switch_to_displ_control ()
 Switch to displ control. More...
 
void switch_to_force_control ()
 Switch to force control. More...
 
void doc_solution ()
 Doc the solution. More...
 
 ContactProblem ()
 Constructor. More...
 
 ~ContactProblem ()
 Destructor (empty) More...
 
void actions_after_newton_solve ()
 Update the problem specs after solve (empty) More...
 
void actions_before_newton_solve ()
 Update the problem specs before solve (empty) More...
 
void actions_before_adapt ()
 Actions before adapt: wipe contact elements. More...
 
void actions_after_adapt ()
 
void switch_to_displ_control ()
 Switch to displ control. More...
 
void switch_to_force_control ()
 Switch to force control. More...
 
void doc_solution ()
 Doc the solution. More...
 
 ContactProblem ()
 Constructor. More...
 
 ~ContactProblem ()
 Destructor (empty) More...
 
void actions_after_newton_solve ()
 Update the problem specs after solve (empty) More...
 
void actions_before_newton_solve ()
 Update the problem specs before solve (empty) More...
 
void actions_before_implicit_timestep ()
 Actions before next timestep. More...
 
void actions_before_adapt ()
 Actions before adapt: wipe contact elements. More...
 
void actions_after_adapt ()
 
void doc_solution ()
 Doc the solution. More...
 
double global_temporal_error_norm ()
 Dummy global error norm for adaptive time-stepping. More...
 
 ContactProblem ()
 Constructor. More...
 
 ~ContactProblem ()
 Destructor (empty) More...
 
void actions_after_newton_solve ()
 Update the problem specs after solve (empty) More...
 
void actions_before_newton_solve ()
 Update the problem specs before solve (empty) More...
 
void actions_before_adapt ()
 Actions before adapt: wipe contact elements. More...
 
void actions_after_adapt ()
 
void switch_to_displ_control ()
 Switch to displ control. More...
 
void switch_to_force_control ()
 Switch to force control. More...
 
void doc_solution ()
 Doc the solution. More...
 
- 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
 
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 Types

enum  {
  Bottom_boundary_id , Right_boundary_id , Contact_boundary_id , Left_boundary_id ,
  Right_top_boundary_id , Left_top_boundary_id , Boulder_top_boundary_id , Boulder_bottom_boundary_id
}
 
enum  {
  Bottom_boundary_id , Right_boundary_id , Contact_boundary_id , Left_boundary_id ,
  Right_top_boundary_id , Left_top_boundary_id
}
 

Private Member Functions

void create_displ_imposition_elements ()
 
void delete_displ_imposition_elements ()
 
void create_contact_heat_elements_on_boulder ()
 Create "contact" heat flux elements on surface of boulder. More...
 
void create_imposed_heat_flux_elements_on_boulder ()
 Create imposed heat flux elements on surface of boulder. More...
 
void create_contact_elements ()
 Create contact elements. More...
 
void delete_contact_elements ()
 Delete contact elements. More...
 
void complete_problem_setup ()
 
void create_displ_imposition_elements ()
 
void delete_displ_imposition_elements ()
 
void create_contact_elements ()
 Create contact elements. More...
 
void delete_contact_elements ()
 Delete contact elements. More...
 
void complete_problem_setup ()
 
void create_contact_elements ()
 Create contact elements. More...
 
void delete_contact_elements ()
 Delete contact elements. More...
 
void complete_problem_setup ()
 
void create_contact_elements ()
 Create contact elements. More...
 
void delete_contact_elements ()
 Delete contact elements. More...
 
void complete_problem_setup ()
 

Private Attributes

RefineableTriangleMesh< ProjectableUnsteadyHeatElement< TUnsteadyHeatElement< 2, 3 > > > * Boulder_mesh_pt
 Pointer to boulder mesh. More...
 
RefineableSolidTriangleMesh< ELEMENT > * Bulk_mesh_pt
 Pointer to bulk mesh. More...
 
Mesh * Surface_contact_mesh_pt
 Pointer to the "surface" contact mesh. More...
 
SolidMesh * Displ_imposition_mesh_pt
 
Mesh * Boulder_surface_heat_flux_mesh_pt
 Mesh of elements that imposed unit heat flux on top of boulder. More...
 
Mesh * Penetrator_mesh_pt
 Penetrator mesh. More...
 
Mesh * Boulder_surface_contact_mesh_pt
 
ofstream Trace_file
 Trace file. More...
 
DocInfo Doc_info
 
BackupMeshForProjection< TElement< 1, 3 > > * Backed_up_surface_contact_mesh_pt
 
SolidNode * Control_node_pt
 
double Xc_old
 x coordinate of old control node More...
 
unsigned Contact_id
 
double X_ll
 x coordinate of lower left corner More...
 
double X_ur
 x coordinate of upper right corner More...
 
double Y_ll
 y coordinate of lower left corner More...
 
double Y_ur
 y coordinate of upper right corner More...
 
TriangleMeshPolyLine * Contact_boundary_pt
 Contact boundary in its poly line representation. More...
 
double Maximum_element_length_on_contact_boundary
 Max. element length on contact boundary. More...
 
double Maximum_element_length_on_boulder_boundary
 
unsigned Contact_boundary_id
 ID of contact boundary. More...
 
unsigned Bottom_boundary_id
 ID of bottom boundary. More...
 
unsigned Left_boundary_id
 ID of left boundary. More...
 
unsigned Right_boundary_id
 ID of right boundary. More...
 
unsigned Left_top_boundary_id
 ID of top left boundary. More...
 
unsigned Right_top_boundary_id
 ID of top right boundary. More...
 

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_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 ()
 
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

template<class ELEMENT>
class ContactProblem< ELEMENT >

Problem class.

///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////

Member Enumeration Documentation

◆ anonymous enum

template<class ELEMENT >
anonymous enum
private
Enumerator
Bottom_boundary_id 

ID of bottom boundary.

Right_boundary_id 

ID of right boundary.

Contact_boundary_id 

ID of contact boundary.

Left_boundary_id 

ID of left boundary.

Right_top_boundary_id 
Left_top_boundary_id 
Boulder_top_boundary_id 
Boulder_bottom_boundary_id 
3947  {
3956  };
@ Boulder_bottom_boundary_id
Definition: heated_linear_solid_contact_with_gravity.cc:3955
@ Boulder_top_boundary_id
Definition: heated_linear_solid_contact_with_gravity.cc:3954
@ Right_boundary_id
ID of right boundary.
Definition: heated_linear_solid_contact_with_gravity.cc:3949
@ Bottom_boundary_id
ID of bottom boundary.
Definition: heated_linear_solid_contact_with_gravity.cc:3948
@ Left_top_boundary_id
Definition: heated_linear_solid_contact_with_gravity.cc:3953
@ Right_top_boundary_id
Definition: heated_linear_solid_contact_with_gravity.cc:3952
@ Left_boundary_id
ID of left boundary.
Definition: heated_linear_solid_contact_with_gravity.cc:3951
@ Contact_boundary_id
ID of contact boundary.
Definition: heated_linear_solid_contact_with_gravity.cc:3950

◆ anonymous enum

template<class ELEMENT >
anonymous enum
private
Enumerator
Bottom_boundary_id 

ID of bottom boundary.

Right_boundary_id 

ID of right boundary.

Contact_boundary_id 

ID of contact boundary.

Left_boundary_id 

ID of left boundary.

Right_top_boundary_id 
Left_top_boundary_id 

Constructor & Destructor Documentation

◆ ContactProblem() [1/4]

template<class ELEMENT >
ContactProblem< ELEMENT >::ContactProblem

Constructor.

Constructor for contact problem in square domain.

4016 {
4017 
4018  // Allow for crap initial guess or convergence problems...
4020 
4021  // Initialise
4022  Control_node_pt=0;
4023 
4024  // Initialise x coordinate of old control node (will be overwritten
4025  // when needed)
4026  Xc_old=0.0;
4027 
4028  // ID of additional nodal values created by contact elements to store
4029  // contact pressure/Lagr. mult.
4030  Contact_id=1;
4031 
4032  // Initialise
4034 
4035  // Output directory
4036  Doc_info.set_directory("RESLT");
4037 
4038  // Output number
4039  Doc_info.number()=0;
4040 
4041  // Open trace file
4042  Trace_file.open("RESLT/trace.dat");
4043 
4044  // Set boundaries of domain
4045  X_ll=0.0;;
4046  X_ur=1.0;
4047  Y_ll=0.0;
4048  Y_ur=1.0;
4049 
4050  // Allocate the timestepper -- this constructs the Problem's
4051  // time object with a sufficient amount of storage to store the
4052  // previous timsteps.
4054 
4055  // Create boulder mesh
4056  //--------------------
4057 
4058  // Create circle representing outer boundary
4059  Circle* boulder_boundary_circle_pt=new Circle(ProblemParameters::Centre[0],
4062 
4063  // Provide storage for pointers to the two parts of the curvilinear boundary
4064  Vector<TriangleMeshCurveSection*> boulder_outer_curvilinear_boundary_pt(2);
4065 
4066  // First bit
4067  double zeta_start=0.0;
4068  double zeta_end=MathematicalConstants::Pi;
4069  unsigned n_seg_boulder=1000;
4070  boulder_outer_curvilinear_boundary_pt[0]=new TriangleMeshCurviLine(
4071  boulder_boundary_circle_pt,zeta_start,zeta_end,n_seg_boulder,
4073 
4074  // Second bit
4075  zeta_start=MathematicalConstants::Pi;
4076  zeta_end=2.0*MathematicalConstants::Pi;
4077  boulder_outer_curvilinear_boundary_pt[1]=new TriangleMeshCurviLine(
4078  boulder_boundary_circle_pt,zeta_start,zeta_end,n_seg_boulder,
4080 
4081  // Combine to curvilinear boundary and define the
4082  // outer boundary
4083  TriangleMeshClosedCurve* boulder_outer_boundary_pt=
4084  new TriangleMeshClosedCurve(boulder_outer_curvilinear_boundary_pt);
4085 
4086  // Use the TriangleMeshParameters object for helping on the manage
4087  // of the TriangleMesh parameters. The only parameter that needs to take
4088  // is the outer boundary.
4089  TriangleMeshParameters boulder_triangle_mesh_parameters(
4090  boulder_outer_boundary_pt);
4091 
4092  // Target element size in boulder mesh
4093  boulder_triangle_mesh_parameters.element_area() =
4095 
4096  // Build boulder mesh
4100  boulder_triangle_mesh_parameters,
4101  time_stepper_pt());
4102 
4103  // Set error estimator for bulk mesh
4104  Z2ErrorEstimator* boulder_error_estimator_pt=new Z2ErrorEstimator;
4105  Boulder_mesh_pt->spatial_error_estimator_pt()=boulder_error_estimator_pt;
4106 
4107  // Set element size limits (faire: just uncommented both of these)
4110 
4111 
4112 #ifdef STRUCTURED_MESH
4113 
4114  // Create the mesh
4115 
4116  // # of elements in x-direction
4117  unsigned n_x=11;
4118 
4119  // # of elements in y-direction
4120  unsigned n_y=10;
4121 
4122  // Domain length in x-direction
4123  double l_x=1.0;
4124 
4125  // Domain length in y-direction
4126  double l_y=1.0;
4127 
4128  //Now create the mesh
4130  (n_x,n_y,l_x,l_y,time_stepper_pt());
4131 
4132 #else
4133 
4134  // Pointer to the closed curve that defines the outer boundary
4135  TriangleMeshClosedCurve* closed_curve_pt=0;
4136 
4137  // Build outer boundary as Polygon
4138 
4139  // The boundary is bounded by five distinct boundaries, each
4140  // represented by its own polyline
4141  Vector<TriangleMeshCurveSection*> boundary_polyline_pt(6);
4142 
4143  // Vertex coordinates on boundary
4144  Vector<Vector<double> > bound_coords(2);
4145 
4146  // Left boundary
4147  bound_coords[0].resize(2);
4148  bound_coords[0][0]=X_ll;
4149  bound_coords[0][1]=Y_ur;
4150 
4151  bound_coords[1].resize(2);
4152  bound_coords[1][0]=X_ll;
4153  bound_coords[1][1]=Y_ll;
4154 
4155  // Build the boundary polyline
4156  boundary_polyline_pt[0]=new TriangleMeshPolyLine(bound_coords,
4158 
4159  // Bottom boundary
4160  bound_coords[0].resize(2);
4161  bound_coords[0][0]=X_ll;
4162  bound_coords[0][1]=Y_ll;
4163 
4164  bound_coords[1].resize(2);
4165  bound_coords[1][0]=X_ur;
4166  bound_coords[1][1]=Y_ll;
4167 
4168  // Build the boundary polyline
4169  boundary_polyline_pt[1]=new TriangleMeshPolyLine(bound_coords,
4171 
4172  // Right boundary
4173  bound_coords[0].resize(2);
4174  bound_coords[0][0]=X_ur;
4175  bound_coords[0][1]=Y_ll;
4176 
4177  bound_coords[1].resize(2);
4178  bound_coords[1][0]=X_ur;
4179  bound_coords[1][1]=Y_ur;
4180 
4181  // Build the boundary polyline
4182  boundary_polyline_pt[2]=new TriangleMeshPolyLine(bound_coords,
4184 
4185 
4186  // Right top boundary
4187  unsigned npt_right=10;
4188  Vector<Vector<double> > right_top_bound_coords(npt_right);
4189  right_top_bound_coords[0].resize(2);
4190  right_top_bound_coords[0][0]=X_ur;
4191  right_top_bound_coords[0][1]=Y_ur;
4192  for (unsigned j=1;j<npt_right;j++)
4193  {
4194  right_top_bound_coords[j].resize(2);
4196  double(npt_right-1);
4197  double y=Y_ur;
4198  right_top_bound_coords[j][0]=x;
4199  right_top_bound_coords[j][1]=y;
4200  }
4201 
4202 
4203  // Build boundary poly line
4204  TriangleMeshPolyLine* right_top_boundary_pt=
4205  new TriangleMeshPolyLine(right_top_bound_coords,
4207  boundary_polyline_pt[3]=right_top_boundary_pt;
4208 
4209 
4210  // Contact boundary
4211  unsigned npt_contact=200; // hierher 20;
4212  Vector<Vector<double> > contact_bound_coords(npt_contact);
4213  contact_bound_coords[0].resize(2);
4214  contact_bound_coords[0][0]=right_top_bound_coords[npt_right-1][0];
4215  contact_bound_coords[0][1]=right_top_bound_coords[npt_right-1][1];
4216  for (unsigned j=1;j<npt_contact;j++)
4217  {
4218  contact_bound_coords[j].resize(2);
4221  ProblemParameters::X_contact_end_left)*double(j)/double(npt_contact-1);
4222  double y=Y_ur;
4223  contact_bound_coords[j][0]=x;
4224  contact_bound_coords[j][1]=y;
4225  }
4226 
4227  // Build boundary poly line
4229  new TriangleMeshPolyLine(contact_bound_coords,
4231  boundary_polyline_pt[4]=Contact_boundary_pt;
4232 
4233  // Keep elements near boundary nice and small
4237 
4238  // Left top boundary
4239  unsigned npt_left=15;
4240  Vector<Vector<double> > top_left_bound_coords(npt_left);
4241  top_left_bound_coords[0].resize(2);
4242  top_left_bound_coords[0][0]=contact_bound_coords[npt_contact-1][0];
4243  top_left_bound_coords[0][1]=contact_bound_coords[npt_contact-1][1];
4244  for (unsigned j=1;j<npt_left-1;j++)
4245  {
4246  top_left_bound_coords[j].resize(2);
4249  double(j)/double(npt_left-1);
4250  double y=Y_ur;
4251  top_left_bound_coords[j][0]=x;
4252  top_left_bound_coords[j][1]=y;
4253  }
4254  top_left_bound_coords[npt_left-1].resize(2);
4255  top_left_bound_coords[npt_left-1][0]=X_ll;
4256  top_left_bound_coords[npt_left-1][1]=Y_ur;
4257 
4258  // Build boundary poly line
4259  TriangleMeshPolyLine* top_left_boundary_pt=
4260  new TriangleMeshPolyLine(top_left_bound_coords,
4262  boundary_polyline_pt[5]=top_left_boundary_pt;
4263 
4264  // Create the triangle mesh polygon for outer boundary
4265  //----------------------------------------------------
4266  TriangleMeshPolygon *outer_polygon =
4267  new TriangleMeshPolygon(boundary_polyline_pt);
4268 
4269  // Set the pointer
4270  closed_curve_pt = outer_polygon;
4271 
4272  // Now build the mesh
4273  //===================
4274 
4275  // Use the TriangleMeshParameters object for helping on the manage of the
4276  // TriangleMesh parameters
4277  TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
4278 
4279  // Specify the maximum area element
4280  double uniform_element_area=ProblemParameters::El_area;
4281  triangle_mesh_parameters.element_area() = uniform_element_area;
4282 
4283  // Create the mesh
4284  Bulk_mesh_pt=
4285  new RefineableSolidTriangleMesh<ELEMENT>(triangle_mesh_parameters,
4286  time_stepper_pt());
4287 
4288 #endif
4289 
4290 
4291 #ifdef STRUCTURED_MESH
4292 
4293  Vector<double> lower_left(2);
4294  lower_left[0]=0.3;
4295  lower_left[1]=0.8;
4296  Vector<double> upper_right(2);
4297  upper_right[0]=0.7;
4298  upper_right[1]=1.0;
4299  unsigned central_node_number=4;
4300  Bulk_mesh_pt->spatial_error_estimator_pt()=
4302  lower_left,
4303  upper_right,
4304  central_node_number);
4305 
4306 #else
4307 
4308  // Set error estimator for bulk mesh
4310  Bulk_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
4311 
4312  // Set element size limits
4313  Bulk_mesh_pt->max_element_size()=ProblemParameters::El_area;
4314  Bulk_mesh_pt->min_element_size()=0.1*ProblemParameters::El_area;
4315 
4316 #endif
4317 
4318  // Create "surface" mesh on boulder that applies contact heat flux boundary
4319  // condition to boulder
4322 
4323  // Create "surface" mesh on boulder that applies flux boundary
4324  // condition to boulder
4327 
4328  // Create contact elements
4331 
4332  // Create elements that enforce prescribed boundary motion
4333  // by Lagrange multipliers
4336 
4337  // Create mesh for penetrator element
4339 
4340  // Set boundary condition and complete the build of all elements
4342 
4343  // Initial assigment
4345 
4346  // Add the sub meshes to the problem
4354 
4355  // Combine all submeshes into a single global Mesh
4357 
4358  // hierher
4359  //Problem::Newton_solver_tolerance=1.0e-6;
4360 
4361  // hierher
4362  //linear_solver_pt()=new FD_LU;
4363 
4364  // Do equation numbering
4365  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
4366 
4367 } // end of constructor
double X_ll
x coordinate of lower left corner
Definition: heated_linear_solid_contact_with_gravity.cc:3990
unsigned Contact_id
Definition: heated_linear_solid_contact_with_gravity.cc:3987
void create_imposed_heat_flux_elements_on_boulder()
Create imposed heat flux elements on surface of boulder.
Definition: heated_linear_solid_contact_with_gravity.cc:3402
SolidNode * Control_node_pt
Definition: heated_linear_solid_contact_with_gravity.cc:3980
DocInfo Doc_info
Definition: heated_linear_solid_contact_with_gravity.cc:3962
RefineableTriangleMesh< ProjectableUnsteadyHeatElement< TUnsteadyHeatElement< 2, 3 > > > * Boulder_mesh_pt
Pointer to boulder mesh.
Definition: heated_linear_solid_contact_with_gravity.cc:3915
Mesh * Surface_contact_mesh_pt
Pointer to the "surface" contact mesh.
Definition: heated_linear_solid_contact_with_gravity.cc:3930
SolidMesh * Displ_imposition_mesh_pt
Definition: heated_linear_solid_contact_with_gravity.cc:3934
RefineableSolidTriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to bulk mesh.
Definition: heated_linear_solid_contact_with_gravity.cc:3925
double Y_ll
y coordinate of lower left corner
Definition: heated_linear_solid_contact_with_gravity.cc:3996
TriangleMeshPolyLine * Contact_boundary_pt
Contact boundary in its poly line representation.
Definition: heated_linear_solid_contact_with_gravity.cc:4002
void create_contact_elements()
Create contact elements.
Definition: heated_linear_solid_contact_with_gravity.cc:3437
void create_contact_heat_elements_on_boulder()
Create "contact" heat flux elements on surface of boulder.
Definition: heated_linear_solid_contact_with_gravity.cc:3367
double Maximum_element_length_on_contact_boundary
Max. element length on contact boundary.
Definition: heated_linear_solid_contact_with_gravity.cc:4005
Mesh * Penetrator_mesh_pt
Penetrator mesh.
Definition: heated_linear_solid_contact_with_gravity.cc:3940
double X_ur
x coordinate of upper right corner
Definition: heated_linear_solid_contact_with_gravity.cc:3993
double Xc_old
x coordinate of old control node
Definition: heated_linear_solid_contact_with_gravity.cc:3983
Mesh * Boulder_surface_contact_mesh_pt
Definition: heated_linear_solid_contact_with_gravity.cc:3944
ofstream Trace_file
Trace file.
Definition: heated_linear_solid_contact_with_gravity.cc:3959
void complete_problem_setup()
Definition: heated_linear_solid_contact_with_gravity.cc:3517
double Y_ur
y coordinate of upper right corner
Definition: heated_linear_solid_contact_with_gravity.cc:3999
BackupMeshForProjection< TElement< 1, 3 > > * Backed_up_surface_contact_mesh_pt
Definition: heated_linear_solid_contact_with_gravity.cc:3974
Mesh * Boulder_surface_heat_flux_mesh_pt
Mesh of elements that imposed unit heat flux on top of boulder.
Definition: heated_linear_solid_contact_with_gravity.cc:3937
void create_displ_imposition_elements()
Definition: heated_linear_solid_contact_with_gravity.cc:3125
Definition: geom_objects.h:873
void set_directory(const std::string &directory)
Definition: oomph_utilities.cc:298
unsigned & number()
Number used (e.g.) for labeling output files.
Definition: oomph_utilities.h:554
Definition: error_estimator.h:435
Definition: rectangular_quadmesh.template.h:573
Definition: mesh.h:67
Definition: timesteppers.h:1074
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
void add_time_stepper_pt(TimeStepper *const &time_stepper_pt)
Definition: problem.cc:1545
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Definition: problem.h:1330
void build_global_mesh()
Definition: problem.cc:1493
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Definition: problem.cc:1989
TimeStepper *& time_stepper_pt()
Definition: problem.h:1524
UnsteadyHeat upgraded to become projectable.
Definition: unsteady_heat_elements.h:699
ErrorEstimator *& spatial_error_estimator_pt()
Access to spatial error estimator.
Definition: refineable_mesh.h:143
Unstructured refineable Triangle Mesh.
Definition: triangle_mesh.template.h:2249
double & max_element_size()
Max element size allowed during adaptation.
Definition: triangle_mesh.template.h:2477
double & min_element_size()
Min element size allowed during adaptation.
Definition: triangle_mesh.template.h:2483
Definition: mesh.h:2562
Definition: Tunsteady_heat_elements.h:62
Base class defining a closed curve for the Triangle mesh generation.
Definition: unstructured_two_d_mesh_geometry_base.h:1339
void set_maximum_length(const double &maximum_length)
Definition: unstructured_two_d_mesh_geometry_base.h:284
Definition: unstructured_two_d_mesh_geometry_base.h:662
Definition: triangle_mesh.template.h:94
Class defining a polyline for use in Triangle Mesh generation.
Definition: unstructured_two_d_mesh_geometry_base.h:868
Class defining a closed polygon for the Triangle mesh generation.
Definition: unstructured_two_d_mesh_geometry_base.h:1451
Definition: oomph-lib/src/generic/Vector.h:58
Definition: error_estimator.h:266
Scalar * y
Definition: level1_cplx_impl.h:128
double Pi
Definition: two_d_biharmonic.cc:235
Z2ErrorEstimator * error_estimator_pt
Definition: MortaringCantileverCompareToNonMortaring.cpp:190
double El_area
Initial/max element area.
Definition: heated_linear_solid_contact_with_gravity.cc:2857
Vector< double > Centre
Position of centre of penetrator.
Definition: circular_boulder_melt.cc:85
double Radius
Radius of penetrator.
Definition: circular_boulder_melt.cc:79
double X_contact_end_right
Right end of contact region (for unstructured mesh only)
Definition: heated_linear_solid_contact_with_gravity.cc:2796
double X_contact_end_left
Left end of contact region (for unstructured mesh only)
Definition: heated_linear_solid_contact_with_gravity.cc:2793
double Y_c
... OR THESE...
Definition: heated_linear_solid_contact_with_gravity.cc:2846
unsigned Max_newton_iterations
Maximum number of newton iterations.
Definition: elements.cc:1654
list x
Definition: plotDoE.py:28
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References ProblemParameters::Centre, GlobalParameters::Doc_info, ProblemParameters::El_area, oomph::TriangleMeshParameters::element_area(), MeshRefinement::error_estimator_pt, j, oomph::Locate_zeta_helpers::Max_newton_iterations, oomph::DocInfo::number(), BiharmonicTestFunctions2::Pi, ProblemParameters::Radius, oomph::DocInfo::set_directory(), oomph::Problem_Parameter::Trace_file, plotDoE::x, ProblemParameters::X_contact_end_left, ProblemParameters::X_contact_end_right, y, and ProblemParameters::Y_c.

◆ ~ContactProblem() [1/4]

template<class ELEMENT >
ContactProblem< ELEMENT >::~ContactProblem ( )
inline

Destructor (empty)

2900 {}

◆ ContactProblem() [2/4]

template<class ELEMENT >
ContactProblem< ELEMENT >::ContactProblem ( )

Constructor.

◆ ~ContactProblem() [2/4]

template<class ELEMENT >
ContactProblem< ELEMENT >::~ContactProblem ( )
inline

Destructor (empty)

923 {}

◆ ContactProblem() [3/4]

template<class ELEMENT >
ContactProblem< ELEMENT >::ContactProblem ( )

Constructor.

◆ ~ContactProblem() [3/4]

template<class ELEMENT >
ContactProblem< ELEMENT >::~ContactProblem ( )
inline

Destructor (empty)

112 {}

◆ ContactProblem() [4/4]

template<class ELEMENT >
ContactProblem< ELEMENT >::ContactProblem ( )

Constructor.

◆ ~ContactProblem() [4/4]

template<class ELEMENT >
ContactProblem< ELEMENT >::~ContactProblem ( )
inline

Destructor (empty)

705 {}

Member Function Documentation

◆ actions_after_adapt() [1/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_after_adapt ( )
inlinevirtual

Actions after adapt: Setup the problem again – remember that the mesh has been completely rebuilt and its element's don't have any pointers to source fcts etc. yet

Reimplemented from oomph::Problem.

3023  {
3024  // Create "contact" heat flux elements on surface of boulder
3026 
3027  // Create normal heat flux elements on surface of boulder
3029 
3030  // Create contact elements
3032 
3033  // Create elements that impose displacement of melt line
3035 
3036  // Now project from backup of original contact mesh to new one
3039 
3040  // For maximum stability: Reset the current nodal positions to be
3041  // the "stress-free" ones -- this assignment means that the
3042  // parameter study no longer corresponds to a physical experiment
3043  // but is what we'd do if we wanted to use the solid solve
3044  // to update a fluid mesh in an FSI problem, say.
3045  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
3046 
3047  // Rebuild elements
3049 
3050  // Rebuild the Problem's global mesh from its various sub-meshes
3052 
3053  // Kill backed up mesh
3056 
3057  // Output impose displ elements
3058  ofstream some_file;
3059  char filename[100];
3060  sprintf(filename,"impose_after.dat");
3061  some_file.open(filename);
3062  Displ_imposition_mesh_pt->output(some_file);
3063  some_file.close();
3064 
3065  }
void project_onto_new_mesh(Mesh *new_mesh_pt)
Definition: face_mesh_project.h:394
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
void rebuild_global_mesh()
Definition: problem.cc:1533
string filename
Definition: MergeRestartFiles.py:39

References ContactProblem< ELEMENT >::Backed_up_surface_contact_mesh_pt, ContactProblem< ELEMENT >::Bulk_mesh_pt, ContactProblem< ELEMENT >::complete_problem_setup(), ContactProblem< ELEMENT >::create_contact_elements(), ContactProblem< ELEMENT >::create_contact_heat_elements_on_boulder(), ContactProblem< ELEMENT >::create_displ_imposition_elements(), ContactProblem< ELEMENT >::create_imposed_heat_flux_elements_on_boulder(), ContactProblem< ELEMENT >::Displ_imposition_mesh_pt, MergeRestartFiles::filename, oomph::Mesh::output(), oomph::BackupMeshForProjection< GEOMETRIC_ELEMENT >::project_onto_new_mesh(), oomph::Problem::rebuild_global_mesh(), and ContactProblem< ELEMENT >::Surface_contact_mesh_pt.

◆ actions_after_adapt() [2/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_after_adapt ( )
inlinevirtual

Actions after adapt: Setup the problem again – remember that the mesh has been completely rebuilt and its element's don't have any pointers to source fcts etc. yet

Reimplemented from oomph::Problem.

1003  {
1004  // Create contact elements
1006 
1007  // Create elements that impose displacement of melt line
1009 
1010  // Now project from backup of original contact mesh to new one
1013 
1014  // For maximum stability: Reset the current nodal positions to be
1015  // the "stress-free" ones -- this assignment means that the
1016  // parameter study no longer corresponds to a physical experiment
1017  // but is what we'd do if we wanted to use the solid solve
1018  // to update a fluid mesh in an FSI problem, say.
1019  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
1020 
1021  // Rebuild elements
1023 
1024  // Rebuild the Problem's global mesh from its various sub-meshes
1026 
1027  // Kill backed up mesh
1030 
1031  // Output impose displ elements
1032  ofstream some_file;
1033  char filename[100];
1034  sprintf(filename,"impose_after.dat");
1035  some_file.open(filename);
1036  Displ_imposition_mesh_pt->output(some_file);
1037  some_file.close();
1038 
1039  }

References MergeRestartFiles::filename.

◆ actions_after_adapt() [3/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_after_adapt ( )
inlinevirtual

Actions after adapt: Setup the problem again – remember that the mesh has been completely rebuilt and its element's don't have any pointers to source fcts etc. yet

Reimplemented from oomph::Problem.

171  {
172  // Create contact elements
174 
175  // Rebuild the Problem's global mesh from its various sub-meshes
177 
178  // Rebuild elements
180 
181  // Now project from backup of original contact mesh to new one
182  oomph_info << "Projecting contact pressure.\n";
185 
186  // Kill backed up mesh
189 
190  // // Output contact elements
191  // ofstream some_file;
192  // char filename[100];
193  // sprintf(filename,"contact_after.dat");
194  // some_file.open(filename);
195  // unsigned nel=Surface_contact_mesh_pt->nelement();
196  // for (unsigned e=0;e<nel;e++)
197  // {
198  // dynamic_cast<NonlinearSurfaceContactElement<ELEMENT>* >(
199  // Surface_contact_mesh_pt->element_pt(e))->output(some_file);
200  // }
201  // some_file.close();
202  // //pause("done");
203 
204  }
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

References oomph::oomph_info.

◆ actions_after_adapt() [4/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_after_adapt ( )
inlinevirtual

Actions after adapt: Setup the problem again – remember that the mesh has been completely rebuilt and its element's don't have any pointers to source fcts etc. yet

Reimplemented from oomph::Problem.

766  {
767  // Create contact elements
769 
770  // Now project from backup of original contact mesh to new one
773 
774  // Rebuild elements
776 
777  // Rebuild the Problem's global mesh from its various sub-meshes
779 
780  // Kill backed up mesh
783 
784  // // Output contact elements
785  // ofstream some_file;
786  // char filename[100];
787  // sprintf(filename,"contact_after.dat");
788  // some_file.open(filename);
789  // unsigned nel=Surface_contact_mesh_pt->nelement();
790  // for (unsigned e=0;e<nel;e++)
791  // {
792  // dynamic_cast<NonlinearSurfaceContactElement<ELEMENT>* >(
793  // Surface_contact_mesh_pt->element_pt(e))->output(some_file);
794  // }
795  // some_file.close();
796 
797  }

◆ actions_after_newton_solve() [1/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_after_newton_solve ( )
inlinevirtual

Update the problem specs after solve (empty)

Reimplemented from oomph::Problem.

2923 {}

◆ actions_after_newton_solve() [2/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_after_newton_solve ( )
inlinevirtual

Update the problem specs after solve (empty)

Reimplemented from oomph::Problem.

926 {}

◆ actions_after_newton_solve() [3/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_after_newton_solve ( )
inlinevirtual

Update the problem specs after solve (empty)

Reimplemented from oomph::Problem.

115 {}

◆ actions_after_newton_solve() [4/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_after_newton_solve ( )
inlinevirtual

Update the problem specs after solve (empty)

Reimplemented from oomph::Problem.

708 {}

◆ actions_before_adapt() [1/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_adapt ( )
inlinevirtual

Actions before adapt: wipe contact elements.

Reimplemented from oomph::Problem.

2952  {
2953  // Backup x coordinate of old control node
2954  Xc_old=Control_node_pt->x(0);
2955 
2956 #ifdef STRUCTURED_MESH
2957 
2958  // Make backup of surface mesh
2962 #else
2963 
2964  //double max_el_length=Maximum_element_length_on_contact_boundary*
2965  //ProblemParameters::Element_length_factor;
2966 
2967  //Impose max_el_length from el_area, ie derive average element length from El_area then
2968  //apply multiplier
2969  double max_el_length = ProblemParameters::Element_length_factor*
2971 
2972  Bulk_mesh_pt->boundary_polyline_pt(Contact_boundary_id)
2973  ->set_maximum_length(max_el_length);
2974 
2976  ->set_maximum_length(max_el_length);
2977 
2978  // Make backup of surface mesh
2982 #endif
2983 
2984 
2985  // // Output contact elements
2986  // ofstream some_file;
2987  // char filename[100];
2988  // sprintf(filename,"contact_before.dat");
2989  // some_file.open(filename);
2990  // unsigned nel=Surface_contact_mesh_pt->nelement();
2991  // for (unsigned e=0;e<nel;e++)
2992  // {
2993  // dynamic_cast<HeatedLinearSurfaceContactElement<ELEMENT>* >(
2994  // Surface_contact_mesh_pt->element_pt(e))->output(some_file);
2995  // }
2996  // some_file.close();
2997 
2998 
2999  // Output impose displ elements
3000  ofstream some_file;
3001  char filename[100];
3002  sprintf(filename,"impose_before.dat");
3003  some_file.open(filename);
3004  Displ_imposition_mesh_pt->output(some_file);
3005  some_file.close();
3006 
3007  // // Kill the elements and wipe surface mesh
3010 
3011  // Wipe the mesh
3013 
3014  // Rebuild the Problem's global mesh from its various sub-meshes
3016  }
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
void delete_contact_elements()
Delete contact elements.
Definition: heated_linear_solid_contact_with_gravity.cc:3470
void delete_displ_imposition_elements()
Definition: heated_linear_solid_contact_with_gravity.cc:3349
Definition: face_mesh_project.h:223
void flush_element_and_node_storage()
Definition: mesh.h:407
TriangleMeshPolyLine * boundary_polyline_pt(const unsigned &b)
Definition: unstructured_two_d_mesh_geometry_base.h:1937
double Element_length_factor
Factor for element length on contact boundary.
Definition: heated_linear_solid_contact_with_gravity.cc:2861

References ContactProblem< ELEMENT >::Backed_up_surface_contact_mesh_pt, ContactProblem< ELEMENT >::Boulder_bottom_boundary_id, ContactProblem< ELEMENT >::Boulder_mesh_pt, oomph::UnstructuredTwoDMeshGeometryBase::boundary_polyline_pt(), ContactProblem< ELEMENT >::Bulk_mesh_pt, ContactProblem< ELEMENT >::Contact_boundary_id, ContactProblem< ELEMENT >::Control_node_pt, ContactProblem< ELEMENT >::delete_contact_elements(), ContactProblem< ELEMENT >::delete_displ_imposition_elements(), ContactProblem< ELEMENT >::Displ_imposition_mesh_pt, ProblemParameters::El_area, ProblemParameters::Element_length_factor, MergeRestartFiles::filename, oomph::Mesh::flush_element_and_node_storage(), oomph::Mesh::output(), ContactProblem< ELEMENT >::Penetrator_mesh_pt, oomph::Problem::rebuild_global_mesh(), oomph::TriangleMeshCurveSection::set_maximum_length(), sqrt(), ContactProblem< ELEMENT >::Surface_contact_mesh_pt, oomph::Node::x(), and ContactProblem< ELEMENT >::Xc_old.

◆ actions_before_adapt() [2/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_adapt ( )
inlinevirtual

Actions before adapt: wipe contact elements.

Reimplemented from oomph::Problem.

942  {
943  // Backup x coordinate of old control node
945 
946 #ifdef STRUCTURED_MESH
947 
948  // Make backup of surface mesh
952 #else
953 
954  Bulk_mesh_pt->boundary_polyline_pt(Contact_boundary_id)->set_maximum_length(
957 
958  // Make backup of surface mesh
962 #endif
963 
964 
965  // // Output contact elements
966  // ofstream some_file;
967  // char filename[100];
968  // sprintf(filename,"contact_before.dat");
969  // some_file.open(filename);
970  // unsigned nel=Surface_contact_mesh_pt->nelement();
971  // for (unsigned e=0;e<nel;e++)
972  // {
973  // dynamic_cast<LinearSurfaceContactElement<ELEMENT>* >(
974  // Surface_contact_mesh_pt->element_pt(e))->output(some_file);
975  // }
976  // some_file.close();
977 
978 
979  // Output impose displ elements
980  ofstream some_file;
981  char filename[100];
982  sprintf(filename,"impose_before.dat");
983  some_file.open(filename);
984  Displ_imposition_mesh_pt->output(some_file);
985  some_file.close();
986 
987  // // Kill the elements and wipe surface mesh
990 
991  // Wipe the mesh
993 
994  // Rebuild the Problem's global mesh from its various sub-meshes
996  }

References ProblemParameters::Element_length_factor, MergeRestartFiles::filename, and oomph::Mesh::output().

◆ actions_before_adapt() [3/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_adapt ( )
inlinevirtual

Actions before adapt: wipe contact elements.

Reimplemented from oomph::Problem.

140  {
141  // Make backup of surface mesh
145 
146  // // Output contact elements
147  // ofstream some_file;
148  // char filename[100];
149  // sprintf(filename,"contact_before.dat");
150  // some_file.open(filename);
151  // unsigned nel=Surface_contact_mesh_pt->nelement();
152  // for (unsigned e=0;e<nel;e++)
153  // {
154  // dynamic_cast<NonlinearSurfaceContactElement<ELEMENT>* >(
155  // Surface_contact_mesh_pt->element_pt(e))->output(some_file);
156  // }
157  // some_file.close();
158 
159  // // Kill the elements and wipe surface mesh
161 
162  // Rebuild the Problem's global mesh from its various sub-meshes
164  }

◆ actions_before_adapt() [4/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_adapt ( )
inlinevirtual

Actions before adapt: wipe contact elements.

Reimplemented from oomph::Problem.

716  {
717  // Backup x coordinate of old control node
719 
720 #ifdef STRUCTURED_MESH
721 
722  // Make backup of surface mesh
726 #else
727 
728  Bulk_mesh_pt->boundary_polyline_pt(Contact_boundary_id)->set_maximum_length(
731 
732  // Make backup of surface mesh
736 #endif
737 
738  // // Output contact elements
739  // ofstream some_file;
740  // char filename[100];
741  // sprintf(filename,"contact_before.dat");
742  // some_file.open(filename);
743  // unsigned nel=Surface_contact_mesh_pt->nelement();
744  // for (unsigned e=0;e<nel;e++)
745  // {
746  // dynamic_cast<NonlinearSurfaceContactElement<ELEMENT>* >(
747  // Surface_contact_mesh_pt->element_pt(e))->output(some_file);
748  // }
749  // some_file.close();
750 
751  // // Kill the elements and wipe surface mesh
753 
754  // Wipe the mesh
756 
757  // Rebuild the Problem's global mesh from its various sub-meshes
759  }

References ProblemParameters::Element_length_factor.

◆ actions_before_implicit_timestep() [1/2]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_implicit_timestep ( )
inlinevirtual

Actions before timestep.

Reimplemented from oomph::Problem.

2904  {
2905  double dyc=0.00024;
2906  ProblemParameters::Centre[1]-=dyc;
2907  oomph_info << "Re-solving imposed circle pos for yc="
2909  << std::endl;
2910  unsigned b=Contact_boundary_id;
2911  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
2912  for (unsigned j=0;j<nnod;j++)
2913  {
2914  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
2915  // hierher index read out from element, though
2916  // this shouldn't actually be set at all. Kill
2918  }
2919  }
Scalar * b
Definition: benchVecAdd.cpp:17
void set_value(const unsigned &i, const double &value_)
Definition: nodes.h:271
Definition: nodes.h:1686
double T_contact
hierher
Definition: heated_linear_solid_contact_with_gravity.cc:2782

References b, ContactProblem< ELEMENT >::Bulk_mesh_pt, ProblemParameters::Centre, ContactProblem< ELEMENT >::Contact_boundary_id, j, oomph::oomph_info, oomph::Data::set_value(), and ProblemParameters::T_contact.

◆ actions_before_implicit_timestep() [2/2]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_implicit_timestep ( )
inlinevirtual

Actions before next timestep.

Reimplemented from oomph::Problem.

123  {
124  // Amplitude of oscillation
125  double amplitude=0.1;
126 
127  // Update position of centre -- amplitude of oscillation increases
128  // by 5% per period...
129  double time=time_pt()->time();
131  (1.0+0.05*time)*0.5*(1.0-cos(2.0*MathematicalConstants::Pi*time));
132 
133  oomph_info << "Solving for y_c = "
135  << " for time: " << time << std::endl;
136  }
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
double & time()
Return the current value of continuous time.
Definition: problem.cc:11531
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
double Y_c_max
Initial/max y-position.
Definition: solid_contact.cc:80

References ProblemParameters::Centre, cos(), oomph::oomph_info, BiharmonicTestFunctions2::Pi, and ProblemParameters::Y_c_max.

◆ actions_before_newton_solve() [1/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_newton_solve ( )
inlinevirtual

Update the problem specs before solve (empty)

Reimplemented from oomph::Problem.

2927  {
2928  // For maximum stability: Reset the current nodal positions to be
2929  // the "stress-free" ones -- this assignment means that the
2930  // parameter study no longer corresponds to a physical experiment
2931  // but is what we'd do if we wanted to use the solid solve
2932  // to update a fluid mesh in an FSI problem, say.
2933  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
2934 
2935  // DoubleVector r;
2936  // CRDoubleMatrix jac;
2937 
2938  // oomph_info << "SETTING UP JAC FOR OUTPUT OF MOST RECENT JAC\n";
2939  // get_jacobian(r,jac);
2940  // oomph_info << "DONE SETTING UP JAC FOR OUTPUT OF MOST RECENT JAC\n";
2941  // jac.sparse_indexed_output("most_recent_jacobian.dat");
2942 
2943  // ofstream descr_file;
2944  // descr_file.open("most_recent_description.dat");
2945  // describe_dofs(descr_file);
2946  // descr_file.close();
2947  }

References ContactProblem< ELEMENT >::Bulk_mesh_pt.

◆ actions_before_newton_solve() [2/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_newton_solve ( )
inlinevirtual

Update the problem specs before solve (empty)

Reimplemented from oomph::Problem.

930  {
931  // For maximum stability: Reset the current nodal positions to be
932  // the "stress-free" ones -- this assignment means that the
933  // parameter study no longer corresponds to a physical experiment
934  // but is what we'd do if we wanted to use the solid solve
935  // to update a fluid mesh in an FSI problem, say.
936  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
937  }

◆ actions_before_newton_solve() [3/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_newton_solve ( )
inlinevirtual

Update the problem specs before solve (empty)

Reimplemented from oomph::Problem.

118 {}

◆ actions_before_newton_solve() [4/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::actions_before_newton_solve ( )
inlinevirtual

Update the problem specs before solve (empty)

Reimplemented from oomph::Problem.

711 {}

◆ complete_problem_setup() [1/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::complete_problem_setup ( )
inlineprivate

Helper function to (re-)set boundary condition and complete the build of all elements

3518  {
3519 
3520  // Set (pseudo-)solid mechanics properties for all elements
3521  //---------------------------------------------------------
3522  unsigned n_element = Bulk_mesh_pt->nelement();
3523  for(unsigned e=0;e<n_element;e++)
3524  {
3525  //Cast to a solid element
3526  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
3527 
3528  // Set the constitutive law
3529  el_pt->constitutive_law_pt() =
3531 
3532  // Set density to zero
3533  el_pt->LinearElasticityEquations<2>::lambda_sq_pt()=
3535 
3536  // Set density to zero
3537  el_pt->PVDEquationsBase<2>::lambda_sq_pt()=
3539 
3540  // Set the elasticity tensor
3541  el_pt->elasticity_tensor_pt() = &ProblemParameters::E;
3542 
3543  // Disable inertia
3544  el_pt->LinearElasticityEquations<2>::disable_inertia();
3545 
3546  // Disable inertia
3547  el_pt->PVDEquationsBase<2>::disable_inertia();
3548  }
3549 
3550  // Apply boundary conditions for solid
3551  //------------------------------------
3552 
3553  // Bottom: completely pinned
3554  unsigned b=Bottom_boundary_id;
3555  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
3556  for (unsigned j=0;j<nnod;j++)
3557  {
3558  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
3559  nod_pt->pin_position(0);
3560  nod_pt->pin_position(1);
3561  nod_pt->pin(0);
3562  nod_pt->pin(1);
3563  }
3564 
3565  // Sides: Symmetry bcs
3567  nnod=Bulk_mesh_pt->nboundary_node(b);
3568  for (unsigned j=0;j<nnod;j++)
3569  {
3570  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
3571  nod_pt->pin_position(0);
3572  nod_pt->pin(0);
3573  }
3575  nnod=Bulk_mesh_pt->nboundary_node(b);
3576  for (unsigned j=0;j<nnod;j++)
3577  {
3578  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
3579  nod_pt->pin_position(0);
3580  nod_pt->pin(0);
3581  }
3582 
3583 
3584  // Enforce contact at most central or most heavily loaded node
3585  //------------------------------------------------------------
3586  {
3587 
3588  // Update angle?
3589  bool update_angle=false;
3590  double phi_old=0.0;
3592  dynamic_cast<HeatedCircularPenetratorElement*>(
3594  if (pen_el_pt!=0)
3595  {
3596  update_angle=true;
3597  phi_old=pen_el_pt->angle();
3598  }
3599 
3600 
3601  // Find closest/most loaded node
3602  double x_c=0.5;
3603  Control_node_pt=0;
3604  SolidNode* most_central_node_pt=0;
3605  SolidNode* most_loaded_node_pt=0;
3606  double dist_min=DBL_MAX;
3607  double load_max=0.0;
3608  unsigned b=Contact_boundary_id;
3609  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
3610  for (unsigned j=0;j<nnod;j++)
3611  {
3612  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
3613 
3614  // Find closest node
3615  double dist=std::fabs(nod_pt->x(0)-x_c);
3616  if (dist<dist_min)
3617  {
3618  dist_min=dist;
3619  most_central_node_pt=nod_pt;
3620  }
3621 
3622  // Find most loaded node
3623  BoundaryNodeBase *bnod_pt=dynamic_cast<BoundaryNodeBase*>(nod_pt);
3624  unsigned index_of_contact_pressure=
3626  if (nod_pt->value(index_of_contact_pressure)>load_max)
3627  {
3628  load_max=nod_pt->value(index_of_contact_pressure);
3629  most_loaded_node_pt=nod_pt;
3630  }
3631  }
3632 
3633  // No load
3634  if (most_loaded_node_pt==0)
3635  {
3636  oomph_info << "Choosing most central node as control node\n";
3637  Control_node_pt=most_central_node_pt;
3638  }
3639  else
3640  {
3641  oomph_info << "Choosing most loaded node as control node\n";
3642  Control_node_pt=most_loaded_node_pt;
3643  }
3644  oomph_info << "Control node located at: "
3645  << Control_node_pt->x(0) << " "
3646  << Control_node_pt->x(1) << " "
3647  << std::endl;
3648 
3649 
3650  // Update angle
3651  if (update_angle)
3652  {
3653  double phi_new=asin((Xc_old-Control_node_pt->x(0))/
3654  ProblemParameters::Radius+sin(phi_old));
3655  dynamic_cast<HeatedCircularPenetratorElement*>(
3656  ProblemParameters::Penetrator_pt)->set_angle(phi_new);
3657  oomph_info << "Old/new angle: " << phi_old << " "
3658  << phi_new << std::endl;
3659  }
3660 
3661 
3662  //...............................................................
3663 
3664 
3665 
3666  // Set target vertical position of control node to its current
3667  // position. This is fine as initial assignment; it's then over-written
3668  // before the next solve (if the displacement is controlled).
3669  // If/when this function is called during adaptation it's fine too
3670  // because it doesn't impose another displacement increment -- it simply
3671  // maintains the position that was obtained on the previous solve.
3672  unsigned index_of_vertical_displacement=1;
3674  Control_node_pt->value(index_of_vertical_displacement);
3675 
3676  // Index of nodal value at control node that stores the traded
3677  // contact pressure
3678  unsigned index_of_traded_contact_pressure=UINT_MAX;
3679 
3680  // Impose position of centre when penetrator position directly imposed
3681  //--------------------------------------------------------------------
3683  {
3684  // Delete old one
3686 
3687  // Make new one
3691  }
3692  // Compute penetrator position as part of the solution, either by
3693  // --------------------------------------------------------------
3694  // prescribing nodal position or weight
3695  // ------------------------------------
3696  else //--
3697  {
3698  // Loop over face elements to identify the ones that contain
3699  // the control node
3700  unsigned nel=Surface_contact_mesh_pt->nelement();
3701  bool found=false;
3702  for (unsigned e=0;e<nel;e++)
3703  {
3707  unsigned nnod=el_pt->nnode();
3708  for (unsigned j=0;j<nnod;j++)
3709  {
3710  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(el_pt->node_pt(j));
3711  if (solid_nod_pt==Control_node_pt)
3712  {
3713  // Got it!
3714  found=true;
3715 
3716  // Find index at which contact pressure/Lagrange multiplier
3717  // is stored
3718  BoundaryNodeBase *bnod_pt =
3719  dynamic_cast<BoundaryNodeBase*>(solid_nod_pt);
3720 
3721  // Get the index of the first nodal value associated with
3722  // this FaceElement
3723  unsigned new_index_of_traded_contact_pressure=
3725 
3726  // Copy across
3727 #ifdef PARANOID
3728  if (index_of_traded_contact_pressure!=UINT_MAX)
3729  {
3730  if (new_index_of_traded_contact_pressure!=
3731  index_of_traded_contact_pressure)
3732  {
3733  std::stringstream junk;
3734  junk << "Inconsistency in identification of index of traded"
3735  << "contact pressure: "
3736  << new_index_of_traded_contact_pressure
3737  << " != " << index_of_traded_contact_pressure;
3738  throw OomphLibError(
3739  junk.str(),
3742  }
3743  }
3744 #endif
3745  index_of_traded_contact_pressure=
3746  new_index_of_traded_contact_pressure;
3747  }
3748  }
3749  }
3750  if (!found)
3751  {
3752  std::stringstream junk;
3753  junk << "Control node not found!";
3754  throw OomphLibError(
3755  junk.str(),
3758  }
3759 
3760  // Back up old penetrator (if it existed)
3761  HeatedCircularPenetratorElement* old_penetrator_pt=
3762  dynamic_cast<HeatedCircularPenetratorElement*>(
3764 
3765  // Make new one
3766  unsigned index_of_horizontal_displacement=0;
3767  unsigned index_of_vertical_displacement=1;
3770  index_of_traded_contact_pressure,
3771  index_of_horizontal_displacement,
3772  index_of_vertical_displacement,
3774 
3775  // Displacement or weight imposed?
3776  bool impose_displ=true;
3777  if (old_penetrator_pt!=0)
3778  {
3779  if (!old_penetrator_pt->yc_is_imposed())
3780  {
3781  impose_displ=false;
3782  }
3783  }
3784  if (impose_displ)
3785  {
3786  dynamic_cast<HeatedCircularPenetratorElement*>(
3788  }
3789  else
3790  {
3791  dynamic_cast<HeatedCircularPenetratorElement*>(
3792  ProblemParameters::Penetrator_pt)->impose_weight(
3794  }
3795 
3796 
3797  // Angle or horizontal force imposed?
3798  bool impose_angle=true;
3799  if (old_penetrator_pt!=0)
3800  {
3801  if (!old_penetrator_pt->rotation_angle_is_imposed())
3802  {
3803  impose_angle=false;
3804  }
3805  }
3806  if (impose_angle)
3807  {
3808  dynamic_cast<HeatedCircularPenetratorElement*>(
3809  ProblemParameters::Penetrator_pt)->impose_rotation_angle(
3811  }
3812  else
3813  {
3814  dynamic_cast<HeatedCircularPenetratorElement*>(
3815  ProblemParameters::Penetrator_pt)->impose_horizontal_force(
3817  dynamic_cast<HeatedCircularPenetratorElement*>(
3819  old_penetrator_pt->angle());
3820  }
3821 
3822  // Now kill the old one!
3823  delete old_penetrator_pt;
3824 
3825  // Add to mesh
3827  dynamic_cast<HeatedCircularPenetratorElement*>(
3829  }
3830 
3831  // Pass contact elements to penetrator element and declare their
3832  // positions and Lagrange multiplier (contact pressure) values
3833  // to be external data.
3835  dynamic_cast<HeatedCircularPenetratorElement*>(
3837  if (el_pt!=0)
3838  {
3840  }
3841 
3842  } // end of penetrator position computed as part of the solution
3843  // (directly or indirectly)
3844 
3845  // Loop over the contact elements to pass pointer to penetrator
3846  //-------------------------------------------------------------
3847  n_element=Surface_contact_mesh_pt->nelement();
3848  for(unsigned e=0;e<n_element;e++)
3849  {
3850  // Upcast from GeneralisedElement
3854 
3855  // Set pointer to penetrator
3857  }
3858 
3859 
3860 
3861  // Find external elements (heated contact elements) that provide
3862  // the heat flux onto the boulder to maintain continuity of
3863  // temperature
3864  Multi_domain_functions::
3865  Accept_failed_locate_zeta_in_setup_multi_domain_interaction=true;
3866 
3870 
3871  Multi_domain_functions::
3872  Accept_failed_locate_zeta_in_setup_multi_domain_interaction=false;
3873 
3874  // Find external elements (penetrator flux elements) that provide
3875  // the target temperature
3876  Multi_domain_functions::
3877  Accept_failed_locate_zeta_in_setup_multi_domain_interaction=true;
3878 
3883 
3884  Multi_domain_functions::
3885  Accept_failed_locate_zeta_in_setup_multi_domain_interaction=false;
3886 
3887 
3888  // hierher setup interaction
3889 
3890  // // Loop over the heat contact elements to pass pointer to flux
3891  // //------------------------------------------------------------
3892  // n_element=Boulder_surface_contact_mesh_pt->nelement();
3893  // for(unsigned e=0;e<n_element;e++)
3894  // {
3895  // // Upcast from GeneralisedElement
3896  // HeatedPenetratorFluxElement
3897  // <ProjectableUnsteadyHeatElement<TUnsteadyHeatElement<2,3> > >* el_pt =
3898  // dynamic_cast<HeatedPenetratorFluxElement
3899  // <ProjectableUnsteadyHeatElement<TUnsteadyHeatElement<2,3> > >*>(
3900  // Boulder_surface_contact_mesh_pt->element_pt(e));
3901 
3902 
3903  // // Set the pointer to the prescribed flux function
3904  // // hierher setup interaction el_pt->flux_fct_pt() = &ProblemParameters::unit_flux;
3905 
3906  // }
3907 
3908 
3909 
3910 
3911  }
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Definition: heated_linear_solid_contact_with_gravity.cc:2167
double angle() const
Angle of rotation around contact point.
Definition: heated_linear_solid_contact_with_gravity.cc:2301
bool rotation_angle_is_imposed()
Definition: heated_linear_solid_contact_with_gravity.cc:2487
bool yc_is_imposed()
Is vertical positon of control node imposed? If false then weight imposed.
Definition: heated_linear_solid_contact_with_gravity.cc:2463
void set_contact_element_mesh_pt(Mesh *contact_element_mesh_pt)
Definition: heated_linear_solid_contact_with_gravity.cc:2326
Definition: nodes.h:1996
unsigned index_of_first_value_assigned_by_face_element(const unsigned &face_id=0) const
Definition: nodes.h:2061
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Heated circular penetrator.
Definition: contact_elements.h:392
Definition: heated_linear_solid_contact_with_gravity.cc:601
Definition: heated_linear_solid_contact_with_gravity.cc:1434
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:617
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
double value(const unsigned &i) const
Definition: nodes.cc:2408
Definition: oomph_definitions.h:222
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1816
void set_penetrator_pt(Penetrator *penetrator_pt)
Set pointer to penetrator.
Definition: contact_elements.h:844
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 asin(const bfloat16 &a)
Definition: BFloat16.h:634
bool found
Definition: MergeRestartFiles.py:24
Penetrator * Penetrator_pt
Penetrator.
Definition: circular_boulder_melt.cc:88
double Horizontal_force
Horizontal force of penetrator.
Definition: heated_linear_solid_contact_with_gravity.cc:2841
double Rotation_angle
Target rotation angle about control node.
Definition: heated_linear_solid_contact_with_gravity.cc:2849
double Weight
NOTE: WE IMPOSE EITHER THESE ...
Definition: heated_linear_solid_contact_with_gravity.cc:2838
IsotropicElasticityTensor E(Nu)
The elasticity tensor.
bool Impose_position_of_centre
Definition: heated_linear_solid_contact_with_gravity.cc:2815
double Lambda_sq
Timescale ratio (non-dim density)
Definition: unstructured_two_d_curved.cc:54
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
Definition: circular_boulder_melt.cc:76
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
void setup_multi_domain_interaction(Problem *problem_pt, Mesh *const &mesh_pt, Mesh *const &external_mesh_pt, const unsigned &interaction_index=0)
Definition: multi_domain.template.cc:280
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References oomph::Mesh::add_element_pt(), HeatedCircularPenetratorElement::angle(), Eigen::bfloat16_impl::asin(), b, ContactProblem< ELEMENT >::Bottom_boundary_id, ContactProblem< ELEMENT >::Boulder_surface_contact_mesh_pt, ContactProblem< ELEMENT >::Bulk_mesh_pt, ProblemParameters::Centre, ProblemParameters::Constitutive_law_pt, ContactProblem< ELEMENT >::Contact_boundary_id, ContactProblem< ELEMENT >::Contact_id, ContactProblem< ELEMENT >::Control_node_pt, e(), ProblemParameters::E, oomph::Mesh::element_pt(), boost::multiprecision::fabs(), MergeRestartFiles::found, ProblemParameters::Horizontal_force, ProblemParameters::Impose_position_of_centre, oomph::BoundaryNodeBase::index_of_first_value_assigned_by_face_element(), j, ProblemParameters::Lambda_sq, ContactProblem< ELEMENT >::Left_boundary_id, oomph::Mesh::nelement(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, ContactProblem< ELEMENT >::Penetrator_mesh_pt, ProblemParameters::Penetrator_pt, oomph::Data::pin(), oomph::SolidNode::pin_position(), ProblemParameters::Radius, ContactProblem< ELEMENT >::Right_boundary_id, ProblemParameters::Rotation_angle, HeatedCircularPenetratorElement::rotation_angle_is_imposed(), HeatedCircularPenetratorElement::set_contact_element_mesh_pt(), oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt(), oomph::Multi_domain_functions::setup_multi_domain_interaction(), sin(), ContactProblem< ELEMENT >::Surface_contact_mesh_pt, oomph::Node::value(), ProblemParameters::Weight, oomph::Node::x(), ContactProblem< ELEMENT >::Xc_old, ProblemParameters::Y_c, and HeatedCircularPenetratorElement::yc_is_imposed().

Referenced by ContactProblem< ELEMENT >::actions_after_adapt().

◆ complete_problem_setup() [2/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::complete_problem_setup ( )
inlineprivate

Helper function to (re-)set boundary condition and complete the build of all elements

1393  {
1394 
1395  // Set (pseudo-)solid mechanics properties for all elements
1396  //---------------------------------------------------------
1397  unsigned n_element = Bulk_mesh_pt->nelement();
1398  for(unsigned e=0;e<n_element;e++)
1399  {
1400  //Cast to a solid element
1401  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
1402 
1403  // Set the constitutive law
1404  el_pt->constitutive_law_pt() =
1406 
1407  // Set density to zero
1408  el_pt->LinearElasticityEquations<2>::lambda_sq_pt()=
1410 
1411  // Set density to zero
1412  el_pt->PVDEquationsBase<2>::lambda_sq_pt()=
1414 
1415  // Set the elasticity tensor
1416  el_pt->elasticity_tensor_pt() = &ProblemParameters::E;
1417 
1418  // Disable inertia
1419  el_pt->LinearElasticityEquations<2>::disable_inertia();
1420 
1421  // Disable inertia
1422  el_pt->PVDEquationsBase<2>::disable_inertia();
1423 
1424  // Set body force for real solid
1425  el_pt->LinearElasticityEquations<2>::body_force_fct_pt() =
1427 
1428  }
1429 
1430  // Apply boundary conditions for solid
1431  //------------------------------------
1432 
1433  // Bottom: completely pinned
1434  unsigned b=Bottom_boundary_id;
1435  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
1436  for (unsigned j=0;j<nnod;j++)
1437  {
1438  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
1439  nod_pt->pin_position(0);
1440  nod_pt->pin_position(1);
1441  nod_pt->pin(0);
1442  nod_pt->pin(1);
1443  }
1444 
1445  // Sides: Symmetry bcs
1447  nnod=Bulk_mesh_pt->nboundary_node(b);
1448  for (unsigned j=0;j<nnod;j++)
1449  {
1450  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
1451  nod_pt->pin_position(0);
1452  nod_pt->pin(0);
1453  }
1455  nnod=Bulk_mesh_pt->nboundary_node(b);
1456  for (unsigned j=0;j<nnod;j++)
1457  {
1458  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
1459  nod_pt->pin_position(0);
1460  nod_pt->pin(0);
1461  }
1462 
1463  // Enforce contact at most central or most heavily loaded node
1464  //------------------------------------------------------------
1465  {
1466 
1467  // Update angle?
1468  bool update_angle=false;
1469  double phi_old=0.0;
1470  CircularPenetratorElement* pen_el_pt=
1471  dynamic_cast<CircularPenetratorElement*>(
1473  if (pen_el_pt!=0)
1474  {
1475  update_angle=true;
1476  phi_old=pen_el_pt->angle();
1477  }
1478 
1479 
1480  // Find closest/most loaded node
1481  double x_c=0.5;
1482  Control_node_pt=0;
1483  SolidNode* most_central_node_pt=0;
1484  SolidNode* most_loaded_node_pt=0;
1485  double dist_min=DBL_MAX;
1486  double load_max=0.0;
1487  unsigned b=Contact_boundary_id;
1488  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
1489  for (unsigned j=0;j<nnod;j++)
1490  {
1491  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
1492 
1493  // Find closest node
1494  double dist=std::fabs(nod_pt->x(0)-x_c);
1495  if (dist<dist_min)
1496  {
1497  dist_min=dist;
1498  most_central_node_pt=nod_pt;
1499  }
1500 
1501  // Find most loaded node
1502  BoundaryNodeBase *bnod_pt=dynamic_cast<BoundaryNodeBase*>(nod_pt);
1503  unsigned index_of_contact_pressure=
1505  if (nod_pt->value(index_of_contact_pressure)>load_max)
1506  {
1507  load_max=nod_pt->value(index_of_contact_pressure);
1508  most_loaded_node_pt=nod_pt;
1509  }
1510  }
1511 
1512  // No load
1513  if (most_loaded_node_pt==0)
1514  {
1515  oomph_info << "Choosing most central node as control node\n";
1516  Control_node_pt=most_central_node_pt;
1517  }
1518  else
1519  {
1520  oomph_info << "Choosing most loaded node as control node\n";
1521  Control_node_pt=most_loaded_node_pt;
1522  }
1523  oomph_info << "Control node located at: "
1524  << Control_node_pt->x(0) << " "
1525  << Control_node_pt->x(1) << " "
1526  << std::endl;
1527 
1528 
1529  // Update angle
1530  if (update_angle)
1531  {
1532  double phi_new=asin((Xc_old-Control_node_pt->x(0))/
1533  ProblemParameters::Radius+sin(phi_old));
1534  dynamic_cast<CircularPenetratorElement*>(
1535  ProblemParameters::Penetrator_pt)->set_angle(phi_new);
1536  oomph_info << "Old/new angle: " << phi_old << " "
1537  << phi_new << std::endl;
1538  }
1539 
1540 
1541  //...............................................................
1542 
1543 
1544 
1545  // Set target vertical position of control node to its current
1546  // position. This is fine as initial assignment; it's then over-written
1547  // before the next solve (if the displacement is controlled).
1548  // If/when this function is called during adaptation it's fine too
1549  // because it doesn't impose another displacement increment -- it simply
1550  // maintains the position that was obtained on the previous solve.
1551  unsigned index_of_vertical_displacement=1;
1553  Control_node_pt->value(index_of_vertical_displacement);
1554 
1555  // Index of nodal value at control node that stores the traded
1556  // contact pressure
1557  unsigned index_of_traded_contact_pressure=UINT_MAX;
1558 
1559  // Impose position of centre when penetrator position directly imposed
1560  //--------------------------------------------------------------------
1562  {
1563  // Delete old one
1565 
1566  // Make new one
1570  }
1571  // Compute penetrator position as part of the solution, either by
1572  // --------------------------------------------------------------
1573  // prescribing nodal position or weight
1574  // ------------------------------------
1575  else
1576  {
1577  // Loop over face elements to identify the ones that contain
1578  // the control node
1579  unsigned nel=Surface_contact_mesh_pt->nelement();
1580  bool found=false;
1581  for (unsigned e=0;e<nel;e++)
1582  {
1584  dynamic_cast<LinearSurfaceContactElement<ELEMENT>*>(
1586  unsigned nnod=el_pt->nnode();
1587  for (unsigned j=0;j<nnod;j++)
1588  {
1589  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(el_pt->node_pt(j));
1590  if (solid_nod_pt==Control_node_pt)
1591  {
1592  // Got it!
1593  found=true;
1594 
1595  // Find index at which contact pressure/Lagrange multiplier
1596  // is stored
1597  BoundaryNodeBase *bnod_pt =
1598  dynamic_cast<BoundaryNodeBase*>(solid_nod_pt);
1599 
1600  // Get the index of the first nodal value associated with
1601  // this FaceElement
1602  unsigned new_index_of_traded_contact_pressure=
1604 
1605  // Copy across
1606 #ifdef PARANOID
1607  if (index_of_traded_contact_pressure!=UINT_MAX)
1608  {
1609  if (new_index_of_traded_contact_pressure!=
1610  index_of_traded_contact_pressure)
1611  {
1612  std::stringstream junk;
1613  junk << "Inconsistency in identification of index of traded"
1614  << "contact pressure: "
1615  << new_index_of_traded_contact_pressure
1616  << " != " << index_of_traded_contact_pressure;
1617  throw OomphLibError(
1618  junk.str(),
1621  }
1622  }
1623 #endif
1624  index_of_traded_contact_pressure=
1625  new_index_of_traded_contact_pressure;
1626  }
1627  }
1628  }
1629  if (!found)
1630  {
1631  std::stringstream junk;
1632  junk << "Control node not found!";
1633  throw OomphLibError(
1634  junk.str(),
1637  }
1638 
1639  // Back up old penetrator (if it existed)
1640  CircularPenetratorElement* old_penetrator_pt=
1641  dynamic_cast<CircularPenetratorElement*>(
1643 
1644  // Make new one
1645  unsigned index_of_horizontal_displacement=0;
1646  unsigned index_of_vertical_displacement=1;
1649  index_of_traded_contact_pressure,
1650  index_of_horizontal_displacement,
1651  index_of_vertical_displacement,
1653 
1654  // Displacement or weight imposed?
1655  bool impose_displ=true;
1656  if (old_penetrator_pt!=0)
1657  {
1658  if (!old_penetrator_pt->yc_is_imposed())
1659  {
1660  impose_displ=false;
1661  }
1662  }
1663  if (impose_displ)
1664  {
1665  dynamic_cast<CircularPenetratorElement*>(
1667  }
1668  else
1669  {
1670  dynamic_cast<CircularPenetratorElement*>(
1671  ProblemParameters::Penetrator_pt)->impose_weight(
1673  }
1674 
1675 
1676  // Angle or horizontal force imposed?
1677  bool impose_angle=true;
1678  if (old_penetrator_pt!=0)
1679  {
1680  if (!old_penetrator_pt->rotation_angle_is_imposed())
1681  {
1682  impose_angle=false;
1683  }
1684  }
1685  if (impose_angle)
1686  {
1687  dynamic_cast<CircularPenetratorElement*>(
1688  ProblemParameters::Penetrator_pt)->impose_rotation_angle(
1690  }
1691  else
1692  {
1693  dynamic_cast<CircularPenetratorElement*>(
1694  ProblemParameters::Penetrator_pt)->impose_horizontal_force(
1696  dynamic_cast<CircularPenetratorElement*>(
1698  old_penetrator_pt->angle());
1699  }
1700 
1701  // Now kill the old one!
1702  delete old_penetrator_pt;
1703 
1704  // Add to mesh
1706  dynamic_cast<CircularPenetratorElement*>(
1708  }
1709 
1710  // Pass contact elements to penetrator element and declare their
1711  // positions and Lagrange multiplier (contact pressure) values
1712  // to be external data.
1715  if (el_pt!=0)
1716  {
1718  }
1719 
1720  } // end of penetrator position computed as part of the solution
1721  // (directly or indirectly)
1722 
1723  // Loop over the contact elements to pass pointer to penetrator
1724  //-------------------------------------------------------------
1725  n_element=Surface_contact_mesh_pt->nelement();
1726  for(unsigned e=0;e<n_element;e++)
1727  {
1728  // Upcast from GeneralisedElement
1730  dynamic_cast<LinearSurfaceContactElement<ELEMENT>*>(
1732 
1733  // Set pointer to penetrator
1735  }
1736 
1737  }
Definition: linear_solid_contact_with_gravity.cc:256
bool rotation_angle_is_imposed()
Definition: linear_solid_contact_with_gravity.cc:541
bool yc_is_imposed()
Is vertical positon of control node imposed? If false then weight imposed.
Definition: linear_solid_contact_with_gravity.cc:517
double angle() const
Angle of rotation around contact point.
Definition: linear_solid_contact_with_gravity.cc:354
void set_contact_element_mesh_pt(Mesh *contact_element_mesh_pt)
Definition: linear_solid_contact_with_gravity.cc:379
Penetrator – here implemented as a circle.
Definition: contact_elements.h:249
Definition: contact_elements.h:1773
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
The body force function.
Definition: linear_solid_contact_with_gravity.cc:894

References CircularPenetratorElement::angle(), Eigen::bfloat16_impl::asin(), b, ProblemParameters::body_force(), ProblemParameters::Centre, ProblemParameters::Constitutive_law_pt, e(), ProblemParameters::E, boost::multiprecision::fabs(), MergeRestartFiles::found, ProblemParameters::Horizontal_force, ProblemParameters::Impose_position_of_centre, oomph::BoundaryNodeBase::index_of_first_value_assigned_by_face_element(), j, ProblemParameters::Lambda_sq, oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, ProblemParameters::Penetrator_pt, oomph::Data::pin(), oomph::SolidNode::pin_position(), ProblemParameters::Radius, ProblemParameters::Rotation_angle, CircularPenetratorElement::rotation_angle_is_imposed(), CircularPenetratorElement::set_contact_element_mesh_pt(), oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt(), sin(), oomph::Node::value(), ProblemParameters::Weight, oomph::Node::x(), ProblemParameters::Y_c, and CircularPenetratorElement::yc_is_imposed().

◆ complete_problem_setup() [3/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::complete_problem_setup ( )
inlineprivate

Helper function to (re-)set boundary condition and complete the build of all elements

265  {
266 
267  // Set (pseudo-)solid mechanics properties for all elements
268  //---------------------------------------------------------
269  unsigned n_element = Bulk_mesh_pt->nelement();
270  for(unsigned e=0;e<n_element;e++)
271  {
272  //Cast to a solid element
273  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
274 
275  // Set the constitutive law
276  el_pt->constitutive_law_pt() =
278 
279  // Set density to zero
280  el_pt->lambda_sq_pt()=&ProblemParameters::Lambda_sq;
281 
282  // Disable inertia
283  el_pt->disable_inertia();
284  }
285 
286  // Apply boundary conditions for solid
287  //------------------------------------
288 
289  // Bottom: completely pinned
290  unsigned b=Bottom_boundary_id;
291  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
292  for (unsigned j=0;j<nnod;j++)
293  {
294  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
295  nod_pt->pin_position(0);
296  nod_pt->pin_position(1);
297  }
298 
299  // Sides: Symmetry bcs
301  nnod=Bulk_mesh_pt->nboundary_node(b);
302  for (unsigned j=0;j<nnod;j++)
303  {
304  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
305  nod_pt->pin_position(0);
306  }
308  nnod=Bulk_mesh_pt->nboundary_node(b);
309  for (unsigned j=0;j<nnod;j++)
310  {
311  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
312  nod_pt->pin_position(0);
313  }
314 
315 
316  // hierher
317  // if (!CommandLineArgs::command_line_flag_has_been_set("--proper_elasticity"))
318  // {
319  // // Assign the Lagrangian coordinates -- sensible
320  // // because we've completely rebuilt the mesh
321  // // and haven't copied across any Lagrange multipliers
322  // Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
323  // }
324 
325 
326  // Loop over the contact elements to pass pointer to penetrator
327  //-------------------------------------------------------------
328  n_element=Surface_contact_mesh_pt->nelement();
329  for(unsigned e=0;e<n_element;e++)
330  {
331  // Upcast from GeneralisedElement
335 
336  // Set pointer to penetrator
338  }
339  }
Definition: contact_elements.h:1086

References b, ProblemParameters::Constitutive_law_pt, e(), j, ProblemParameters::Lambda_sq, ProblemParameters::Penetrator_pt, oomph::SolidNode::pin_position(), and oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt().

◆ complete_problem_setup() [4/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::complete_problem_setup ( )
inlineprivate

Helper function to (re-)set boundary condition and complete the build of all elements

905  {
906 
907  // Set (pseudo-)solid mechanics properties for all elements
908  //---------------------------------------------------------
909  unsigned n_element = Bulk_mesh_pt->nelement();
910  for(unsigned e=0;e<n_element;e++)
911  {
912  //Cast to a solid element
913  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
914 
915  // Set the constitutive law
916  el_pt->constitutive_law_pt() =
918 
919  // Set density to zero
920  el_pt->lambda_sq_pt()=&ProblemParameters::Lambda_sq;
921 
922  // Disable inertia
923  el_pt->disable_inertia();
924 
925  // Set body force
926  el_pt->body_force_fct_pt() = &ProblemParameters::body_force;
927 
928  }
929 
930  // Apply boundary conditions for solid
931  //------------------------------------
932 
933  // Bottom: completely pinned
934  unsigned b=Bottom_boundary_id;
935  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
936  for (unsigned j=0;j<nnod;j++)
937  {
938  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
939  nod_pt->pin_position(0);
940  nod_pt->pin_position(1);
941  }
942 
943  // Sides: Symmetry bcs
945  nnod=Bulk_mesh_pt->nboundary_node(b);
946  for (unsigned j=0;j<nnod;j++)
947  {
948  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
949  nod_pt->pin_position(0);
950  }
952  nnod=Bulk_mesh_pt->nboundary_node(b);
953  for (unsigned j=0;j<nnod;j++)
954  {
955  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
956  nod_pt->pin_position(0);
957  }
958 
959 
960  // Enforce contact at most central or most heavily loaded node
961  //------------------------------------------------------------
962  {
963 
964  // Update angle?
965  bool update_angle=false;
966  double phi_old=0.0;
967  CircularPenetratorElement* pen_el_pt=
968  dynamic_cast<CircularPenetratorElement*>(
970  if (pen_el_pt!=0)
971  {
972  update_angle=true;
973  phi_old=pen_el_pt->angle();
974  }
975 
976  double x_c=0.5;
977  Control_node_pt=0;
978  SolidNode* most_central_node_pt=0;
979  SolidNode* most_loaded_node_pt=0;
980  double dist_min=DBL_MAX;
981  double load_max=0.0;
982  unsigned b=Contact_boundary_id;
983  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
984  for (unsigned j=0;j<nnod;j++)
985  {
986  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
987 
988  // Find closest node
989  double dist=std::fabs(nod_pt->x(0)-x_c);
990  if (dist<dist_min)
991  {
992  dist_min=dist;
993  most_central_node_pt=nod_pt;
994  }
995 
996  // Find most loaded node
997  BoundaryNodeBase *bnod_pt=dynamic_cast<BoundaryNodeBase*>(nod_pt);
998  unsigned index_of_contact_pressure=
1000  if (nod_pt->value(index_of_contact_pressure)>load_max)
1001  {
1002  load_max=nod_pt->value(index_of_contact_pressure);
1003  most_loaded_node_pt=nod_pt;
1004  }
1005  }
1006 
1007  // No load
1008  if (most_loaded_node_pt==0)
1009  {
1010  oomph_info << "Choosing most central node as control node\n";
1011  Control_node_pt=most_central_node_pt;
1012  }
1013  else
1014  {
1015  oomph_info << "Choosing most loaded node as control node\n";
1016  Control_node_pt=most_loaded_node_pt;
1017  }
1018  oomph_info << "Control node located at: "
1019  << Control_node_pt->x(0) << " "
1020  << Control_node_pt->x(1) << " "
1021  << std::endl;
1022 
1023 
1024  // Update angle
1025  if (update_angle)
1026  {
1027  double phi_new=asin((Xc_old-Control_node_pt->x(0))/
1028  ProblemParameters::Radius+sin(phi_old));
1029  dynamic_cast<CircularPenetratorElement*>(
1030  ProblemParameters::Penetrator_pt)->set_angle(phi_new);
1031  oomph_info << "Old/new angle: " << phi_old << " "
1032  << phi_new << std::endl;
1033  }
1034 
1035 
1036  // Set target vertical position of control node to its current
1037  // position. This is fine as initial assignment; it's then over-written
1038  // before the next solve (if the displacement is controlled).
1039  // If/when this function is called during adaptation it's fine too
1040  // because it doesn't impose another displacement increment -- it simply
1041  // maintains the position that was obtained on the previous solve.
1043 
1044  // Index of nodal value at control node that stores the traded
1045  // contact pressure
1046  unsigned index_of_traded_contact_pressure=UINT_MAX;
1047 
1048  // Impose position of centre when penetrator position directly imposed
1049  //--------------------------------------------------------------------
1051  {
1052  // Delete old one
1054 
1055  // Make new one
1059  }
1060  // Compute penetrator position as part of the solution, either by
1061  // --------------------------------------------------------------
1062  // prescribing nodal position or weight
1063  // ------------------------------------
1064  else
1065  {
1066  // Loop over face elements to identify the ones that contain
1067  // the control node
1068  unsigned nel=Surface_contact_mesh_pt->nelement();
1069  bool found=false;
1070  for (unsigned e=0;e<nel;e++)
1071  {
1075  unsigned nnod=el_pt->nnode();
1076  for (unsigned j=0;j<nnod;j++)
1077  {
1078  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(el_pt->node_pt(j));
1079  if (solid_nod_pt==Control_node_pt)
1080  {
1081  // Got it!
1082  found=true;
1083 
1084  // Find index at which contact pressure/Lagrange multiplier
1085  // is stored
1086  BoundaryNodeBase *bnod_pt =
1087  dynamic_cast<BoundaryNodeBase*>(solid_nod_pt);
1088 
1089  // Get the index of the first nodal value associated with
1090  // this FaceElement
1091  unsigned new_index_of_traded_contact_pressure=
1093 
1094  // Copy across
1095 #ifdef PARANOID
1096  if (index_of_traded_contact_pressure!=UINT_MAX)
1097  {
1098  if (new_index_of_traded_contact_pressure!=
1099  index_of_traded_contact_pressure)
1100  {
1101  std::stringstream junk;
1102  junk << "Inconsistency in identification of index of traded"
1103  << "contact pressure: "
1104  << new_index_of_traded_contact_pressure
1105  << " != " << index_of_traded_contact_pressure;
1106  throw OomphLibError(
1107  junk.str(),
1110  }
1111  }
1112 #endif
1113  index_of_traded_contact_pressure=
1114  new_index_of_traded_contact_pressure;
1115  }
1116  }
1117  }
1118  if (!found)
1119  {
1120  std::stringstream junk;
1121  junk << "Control node not found!";
1122  throw OomphLibError(
1123  junk.str(),
1126  }
1127 
1128  // Back up old penetrator (if it existed)
1129  CircularPenetratorElement* old_penetrator_pt=
1130  dynamic_cast<CircularPenetratorElement*>(
1132 
1133  // Make new one
1136  index_of_traded_contact_pressure,
1138 
1139  // Displacement or weight imposed?
1140  bool impose_displ=true;
1141  if (old_penetrator_pt!=0)
1142  {
1143  if (!old_penetrator_pt->yc_is_imposed())
1144  {
1145  impose_displ=false;
1146  }
1147  }
1148  if (impose_displ)
1149  {
1150  dynamic_cast<CircularPenetratorElement*>(
1152  }
1153  else
1154  {
1155  dynamic_cast<CircularPenetratorElement*>(
1156  ProblemParameters::Penetrator_pt)->impose_weight(
1158  }
1159 
1160 
1161  // Angle or horizontal force imposed?
1162  bool impose_angle=true;
1163  if (old_penetrator_pt!=0)
1164  {
1165  if (!old_penetrator_pt->rotation_angle_is_imposed())
1166  {
1167  impose_angle=false;
1168  }
1169  }
1170  if (impose_angle)
1171  {
1172  dynamic_cast<CircularPenetratorElement*>(
1173  ProblemParameters::Penetrator_pt)->impose_rotation_angle(
1175  }
1176  else
1177  {
1178  dynamic_cast<CircularPenetratorElement*>(
1179  ProblemParameters::Penetrator_pt)->impose_horizontal_force(
1181  dynamic_cast<CircularPenetratorElement*>(
1183  old_penetrator_pt->angle());
1184  }
1185 
1186  // Now kill the old one!
1187  delete old_penetrator_pt;
1188 
1189  // Add to mesh
1191  dynamic_cast<CircularPenetratorElement*>(
1193  }
1194 
1195  // Pass contact elements to penetrator element and declare their
1196  // positions and Lagrange multiplier (contact pressure) values
1197  // to be external data.
1200  if (el_pt!=0)
1201  {
1203  }
1204 
1205  } // end of penetrator position computed as part of the solution
1206  // (directly or indirectly)
1207 
1208 
1209  // Loop over the contact elements to pass pointer to penetrator
1210  //-------------------------------------------------------------
1211  n_element=Surface_contact_mesh_pt->nelement();
1212  for(unsigned e=0;e<n_element;e++)
1213  {
1214  // Upcast from GeneralisedElement
1218 
1219  // Set pointer to penetrator
1221  }
1222 
1223  }

References CircularPenetratorElement::angle(), Eigen::bfloat16_impl::asin(), b, ProblemParameters::body_force(), ProblemParameters::Centre, ProblemParameters::Constitutive_law_pt, e(), boost::multiprecision::fabs(), MergeRestartFiles::found, ProblemParameters::Horizontal_force, ProblemParameters::Impose_position_of_centre, oomph::BoundaryNodeBase::index_of_first_value_assigned_by_face_element(), j, ProblemParameters::Lambda_sq, oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, ProblemParameters::Penetrator_pt, oomph::SolidNode::pin_position(), ProblemParameters::Radius, ProblemParameters::Rotation_angle, CircularPenetratorElement::rotation_angle_is_imposed(), CircularPenetratorElement::set_contact_element_mesh_pt(), oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt(), sin(), oomph::Node::value(), ProblemParameters::Weight, oomph::Node::x(), ProblemParameters::Y_c, and CircularPenetratorElement::yc_is_imposed().

◆ create_contact_elements() [1/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::create_contact_elements ( )
inlineprivate

Create contact elements.

3438  {
3439  // How many bulk elements are adjacent to boundary b?
3440  unsigned b=Contact_boundary_id;
3441  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
3442 
3443  // Loop over the bulk elements adjacent to boundary b?
3444  for(unsigned e=0;e<n_element;e++)
3445  {
3446  // Get pointer to the bulk element that is adjacent to boundary b
3447  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
3448  Bulk_mesh_pt->boundary_element_pt(b,e));
3449 
3450  //What is the face index of element e along boundary b
3451  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
3452 
3453  // Build the corresponding contact element
3454  HeatedLinearSurfaceContactElement<ELEMENT>* contact_element_pt = new
3455  HeatedLinearSurfaceContactElement<ELEMENT>(bulk_elem_pt,face_index,
3456  Contact_id);
3457 
3458  // Pass pointer to penetrator
3460 
3461  //Add the contact element to the surface mesh
3462  Surface_contact_mesh_pt->add_element_pt(contact_element_pt);
3463 
3464  } //end of loop over bulk elements adjacent to boundary b
3465  }

References oomph::Mesh::add_element_pt(), b, ContactProblem< ELEMENT >::Bulk_mesh_pt, ContactProblem< ELEMENT >::Contact_boundary_id, ContactProblem< ELEMENT >::Contact_id, e(), ProblemParameters::Penetrator_pt, oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt(), and ContactProblem< ELEMENT >::Surface_contact_mesh_pt.

Referenced by ContactProblem< ELEMENT >::actions_after_adapt().

◆ create_contact_elements() [2/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::create_contact_elements ( )
inlineprivate

Create contact elements.

1344  {
1345  // How many bulk elements are adjacent to boundary b?
1346  unsigned b=Contact_boundary_id;
1347  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
1348 
1349  // Loop over the bulk elements adjacent to boundary b?
1350  for(unsigned e=0;e<n_element;e++)
1351  {
1352  // Get pointer to the bulk element that is adjacent to boundary b
1353  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
1354  Bulk_mesh_pt->boundary_element_pt(b,e));
1355 
1356  //What is the face index of element e along boundary b
1357  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
1358 
1359  // Build the corresponding contact element
1360  LinearSurfaceContactElement<ELEMENT>* contact_element_pt = new
1361  LinearSurfaceContactElement<ELEMENT>(bulk_elem_pt,face_index,
1362  Contact_id);
1363 
1364  //Add the contact element to the surface mesh
1365  Surface_contact_mesh_pt->add_element_pt(contact_element_pt);
1366 
1367  } //end of loop over bulk elements adjacent to boundary b
1368  }

References b, and e().

◆ create_contact_elements() [3/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::create_contact_elements ( )
inlineprivate

Create contact elements.

217  {
218  // How many bulk elements are adjacent to boundary b?
219  unsigned b=Contact_boundary_id;
220  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
221 
222  // Loop over the bulk elements adjacent to boundary b?
223  for(unsigned e=0;e<n_element;e++)
224  {
225  // Get pointer to the bulk element that is adjacent to boundary b
226  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
227  Bulk_mesh_pt->boundary_element_pt(b,e));
228 
229  //What is the face index of element e along boundary b
230  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
231 
232  // Build the corresponding contact element
233  NonlinearSurfaceContactElement<ELEMENT>* contact_element_pt = new
234  NonlinearSurfaceContactElement<ELEMENT>(bulk_elem_pt,face_index);
235 
236  //Add the contact element to the surface mesh
237  Surface_contact_mesh_pt->add_element_pt(contact_element_pt);
238 
239  } //end of loop over bulk elements adjacent to boundary b
240  }

References b, and e().

◆ create_contact_elements() [4/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::create_contact_elements ( )
inlineprivate

Create contact elements.

856  {
857  // How many bulk elements are adjacent to boundary b?
858  unsigned b=Contact_boundary_id;
859  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
860 
861  // Loop over the bulk elements adjacent to boundary b?
862  for(unsigned e=0;e<n_element;e++)
863  {
864  // Get pointer to the bulk element that is adjacent to boundary b
865  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
866  Bulk_mesh_pt->boundary_element_pt(b,e));
867 
868  //What is the face index of element e along boundary b
869  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
870 
871  // Build the corresponding contact element
872  NonlinearSurfaceContactElement<ELEMENT>* contact_element_pt = new
873  NonlinearSurfaceContactElement<ELEMENT>(bulk_elem_pt,face_index,
874  Contact_id);
875 
876  //Add the contact element to the surface mesh
877  Surface_contact_mesh_pt->add_element_pt(contact_element_pt);
878 
879  } //end of loop over bulk elements adjacent to boundary b
880  }

References b, and e().

◆ create_contact_heat_elements_on_boulder()

template<class ELEMENT >
void ContactProblem< ELEMENT >::create_contact_heat_elements_on_boulder ( )
inlineprivate

Create "contact" heat flux elements on surface of boulder.

3368  {
3369  // How many bulk elements are adjacent to boundary b?
3370  unsigned b=Boulder_bottom_boundary_id;
3371  unsigned n_element = Boulder_mesh_pt->nboundary_element(b);
3372 
3373  // Loop over the bulk elements adjacent to boundary b?
3374  for(unsigned e=0;e<n_element;e++)
3375  {
3376  // Get pointer to the bulk element that is adjacent to boundary b
3380 
3381  //What is the face index of element e along boundary b
3382  int face_index = Boulder_mesh_pt->face_index_at_boundary(b,e);
3383 
3384  // Build the corresponding contact element
3387  contact_element_pt = new
3390  bulk_elem_pt,face_index);
3391 
3392  // Pass pointer to penetrator
3393  contact_element_pt->set_penetrator_pt(ProblemParameters::Penetrator_pt);
3394 
3395  //Add the contact element to the surface mesh
3396  Boulder_surface_contact_mesh_pt->add_element_pt(contact_element_pt);
3397 
3398  } //end of loop over bulk elements adjacent to boundary b
3399  }
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
Definition: mesh.h:896
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840

References oomph::Mesh::add_element_pt(), b, ContactProblem< ELEMENT >::Boulder_bottom_boundary_id, ContactProblem< ELEMENT >::Boulder_mesh_pt, ContactProblem< ELEMENT >::Boulder_surface_contact_mesh_pt, oomph::Mesh::boundary_element_pt(), e(), oomph::Mesh::face_index_at_boundary(), oomph::Mesh::nboundary_element(), and ProblemParameters::Penetrator_pt.

Referenced by ContactProblem< ELEMENT >::actions_after_adapt().

◆ create_displ_imposition_elements() [1/2]

template<class ELEMENT >
void ContactProblem< ELEMENT >::create_displ_imposition_elements ( )
inlineprivate

Create elements that enforce prescribed boundary motion by Lagrange multipliers

3126  {
3127 
3128  Vector<unsigned> boundary_id;
3129  boundary_id.push_back(Contact_boundary_id);
3130 #ifndef STRUCTURED_MESH
3131  Node* left_contact_node_pt=0;
3132  Node* right_contact_node_pt=0;
3133  Node* left_left_top_node_pt=0;
3134  Node* right_left_top_node_pt=0;
3135  Node* left_right_top_node_pt=0;
3136  Node* right_right_top_node_pt=0;
3137  boundary_id.push_back(Left_top_boundary_id);
3138  boundary_id.push_back(Right_top_boundary_id);
3139 #endif
3140  unsigned nb=boundary_id.size();
3141  for (unsigned bb=0;bb<nb;bb++)
3142  {
3143  unsigned b=boundary_id[bb];
3144 
3145  // How many bulk elements are adjacent to boundary b?
3146  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
3147 
3148  // Loop over the bulk elements adjacent to boundary b?
3149  for(unsigned e=0;e<n_element;e++)
3150  {
3151  // Get pointer to the bulk element that is adjacent to boundary b
3152  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
3153  Bulk_mesh_pt->boundary_element_pt(b,e));
3154 
3155  //Find the index of the face of element e along boundary b
3156  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
3157 
3158  // Create new element
3161  bulk_elem_pt,face_index);
3162 
3163  // Add to mesh
3165 
3166 #ifdef STRUCTURED_MESH
3167 
3168  // Set the GeomObject that defines the boundary shape and
3169  // specify which bulk boundary we are attached to (needed to extract
3170  // the boundary coordinate from the bulk nodes)
3173 
3174 #else
3175 
3176  switch(b)
3177  {
3178  case Contact_boundary_id:
3179 
3180  // Set the GeomObject that defines the boundary shape and
3181  // specify which bulk boundary we are attached to (needed to extract
3182  // the boundary coordinate from the bulk nodes)
3185 
3186  {
3187  unsigned nnod=el_pt->nnode();
3188  for (unsigned j=0;j<nnod;j++)
3189  {
3190  Node* nod_pt=el_pt->node_pt(j);
3191  if (nod_pt->is_on_boundary(Left_top_boundary_id))
3192  {
3193  left_contact_node_pt=nod_pt;
3194  }
3195  else if (nod_pt->is_on_boundary(Right_top_boundary_id))
3196  {
3197  right_contact_node_pt=nod_pt;
3198  }
3199 
3200  }
3201  }
3202 
3203  break;
3204 
3205  case Left_top_boundary_id:
3206 
3207  // Set the GeomObject that defines the boundary shape and
3208  // specify which bulk boundary we are attached to (needed to extract
3209  // the boundary coordinate from the bulk nodes)
3212 
3213  {
3214  unsigned nnod=el_pt->nnode();
3215  for (unsigned j=0;j<nnod;j++)
3216  {
3217  Node* nod_pt=el_pt->node_pt(j);
3218  if (nod_pt->is_on_boundary(Left_boundary_id))
3219  {
3220  left_left_top_node_pt=nod_pt;
3221  }
3222  else if (nod_pt->is_on_boundary(Contact_boundary_id))
3223  {
3224  right_left_top_node_pt=nod_pt;
3225  }
3226  }
3227  }
3228 
3229  break;
3230 
3231  case Right_top_boundary_id:
3232 
3233  // Set the GeomObject that defines the boundary shape and
3234  // specify which bulk boundary we are attached to (needed to extract
3235  // the boundary coordinate from the bulk nodes)
3238 
3239  {
3240  unsigned nnod=el_pt->nnode();
3241  for (unsigned j=0;j<nnod;j++)
3242  {
3243  Node* nod_pt=el_pt->node_pt(j);
3244  if (nod_pt->is_on_boundary(Right_boundary_id))
3245  {
3246  right_right_top_node_pt=nod_pt;
3247  }
3248  else if (nod_pt->is_on_boundary(Contact_boundary_id))
3249  {
3250  left_right_top_node_pt=nod_pt;
3251  }
3252  }
3253  }
3254  break;
3255 
3256  default:
3257 
3258  // Never get here...
3259  oomph_info << "Never get here! b = " << b << "\n";
3260  abort();
3261  }
3262 
3263 
3264 #endif
3265 
3266  // Loop over the nodes
3267  unsigned nnod=el_pt->nnode();
3268  for (unsigned j=0;j<nnod;j++)
3269  {
3270  Node* nod_pt = el_pt->node_pt(j);
3271 
3272  // Is the node also on side boundaries?
3273  if ((nod_pt->is_on_boundary(Left_boundary_id))||
3274  (nod_pt->is_on_boundary(Right_boundary_id)))
3275  {
3276  // How many nodal values were used by the "bulk" element
3277  // that originally created this node?
3278  unsigned n_bulk_value=el_pt->nbulk_value(j);
3279 
3280  // hierher: this statement is unlikely to be true. Fix!
3281  // The remaining ones are Lagrange multipliers and we pin them.
3282  unsigned nval=nod_pt->nvalue();
3283  for (unsigned j=n_bulk_value;j<nval;j++)
3284  {
3285  nod_pt->pin(j);
3286  }
3287  }
3288  }
3289  }
3290  }
3291 
3292 #ifndef STRUCTURED_MESH
3293 
3294  Vector<double> zeta_left(1);
3295  Vector<double> zeta_right(1);
3296  left_left_top_node_pt->get_coordinates_on_boundary(Left_top_boundary_id,
3297  zeta_left);
3298  right_left_top_node_pt->get_coordinates_on_boundary(Left_top_boundary_id,
3299  zeta_right);
3300  if (zeta_left[0]>zeta_right[0])
3301  {
3303  oomph_info << "left top is reversed\n";
3304  }
3305  else
3306  {
3308  oomph_info << "left top is not reversed\n";
3309  }
3310 
3311  left_contact_node_pt->get_coordinates_on_boundary(Contact_boundary_id,
3312  zeta_left);
3313  right_contact_node_pt->get_coordinates_on_boundary(Contact_boundary_id,
3314  zeta_right);
3315  if (zeta_left[0]>zeta_right[0])
3316  {
3318  oomph_info << "contact is reversed\n";
3319  }
3320  else
3321  {
3323  oomph_info << "contact is not reversed\n";
3324  }
3325 
3326  left_right_top_node_pt->get_coordinates_on_boundary(Right_top_boundary_id,
3327  zeta_left);
3328  right_right_top_node_pt->get_coordinates_on_boundary(Right_top_boundary_id,
3329  zeta_right);
3330  if (zeta_left[0]>zeta_right[0])
3331  {
3333  oomph_info << "right top is reversed\n";
3334  }
3335  else
3336  {
3338  oomph_info << "right top is not reversed\n";
3339  }
3340 
3341 #endif
3342 
3343 
3344  } // end of create_displ_imposition_elements
void set_reversed()
Local coordinates are reversed.
Definition: heated_linear_solid_contact_with_gravity.cc:2076
void set_non_reversed()
Local coordinates are not reversed.
Definition: heated_linear_solid_contact_with_gravity.cc:2082
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
unsigned & nbulk_value(const unsigned &n)
Definition: elements.h:4845
Definition: solid_traction_elements.h:1107
void set_boundary_shape_geom_object_pt(GeomObject *boundary_shape_geom_object_pt, const unsigned &boundary_number_in_bulk_mesh)
Definition: solid_traction_elements.h:1209
Definition: nodes.h:906
virtual bool is_on_boundary() const
Definition: nodes.h:1373
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Definition: nodes.cc:2379
int nb
Definition: level2_impl.h:286
WarpedLine Boundary_geom_object(0.0)
GeomObject specifying the shape of the boundary: Initially it's flat.
WarpedLine Boundary_geom_object_contact(1.0e-10, X_contact_end_left, X_contact_end_right)
WarpedLine Boundary_geom_object_right(1.0e-10, X_contact_end_right, 1.0)
WarpedLine Boundary_geom_object_left(1.0e-10, 0.0, X_contact_end_left)

References oomph::Mesh::add_element_pt(), b, Global_Physical_Variables::Boundary_geom_object, ProblemParameters::Boundary_geom_object_contact, ProblemParameters::Boundary_geom_object_left, ProblemParameters::Boundary_geom_object_right, ContactProblem< ELEMENT >::Bulk_mesh_pt, ContactProblem< ELEMENT >::Contact_boundary_id, ContactProblem< ELEMENT >::Displ_imposition_mesh_pt, e(), oomph::Node::get_coordinates_on_boundary(), oomph::Node::is_on_boundary(), j, ContactProblem< ELEMENT >::Left_boundary_id, ContactProblem< ELEMENT >::Left_top_boundary_id, nb, oomph::FaceElement::nbulk_value(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::Data::nvalue(), oomph::oomph_info, oomph::Data::pin(), ContactProblem< ELEMENT >::Right_boundary_id, ContactProblem< ELEMENT >::Right_top_boundary_id, oomph::ImposeDisplacementByLagrangeMultiplierElement< ELEMENT >::set_boundary_shape_geom_object_pt(), WarpedLine::set_non_reversed(), and WarpedLine::set_reversed().

Referenced by ContactProblem< ELEMENT >::actions_after_adapt().

◆ create_displ_imposition_elements() [2/2]

template<class ELEMENT >
void ContactProblem< ELEMENT >::create_displ_imposition_elements ( )
inlineprivate

Create elements that enforce prescribed boundary motion by Lagrange multipliers

1103  {
1104 
1105  Vector<unsigned> boundary_id;
1106  boundary_id.push_back(Contact_boundary_id);
1107 #ifndef STRUCTURED_MESH
1108  Node* left_contact_node_pt=0;
1109  Node* right_contact_node_pt=0;
1110  Node* left_left_top_node_pt=0;
1111  Node* right_left_top_node_pt=0;
1112  Node* left_right_top_node_pt=0;
1113  Node* right_right_top_node_pt=0;
1114  boundary_id.push_back(Left_top_boundary_id);
1115  boundary_id.push_back(Right_top_boundary_id);
1116 #endif
1117  unsigned nb=boundary_id.size();
1118  for (unsigned bb=0;bb<nb;bb++)
1119  {
1120  unsigned b=boundary_id[bb];
1121 
1122  // How many bulk elements are adjacent to boundary b?
1123  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
1124 
1125  // Loop over the bulk elements adjacent to boundary b?
1126  for(unsigned e=0;e<n_element;e++)
1127  {
1128  // Get pointer to the bulk element that is adjacent to boundary b
1129  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
1130  Bulk_mesh_pt->boundary_element_pt(b,e));
1131 
1132  //Find the index of the face of element e along boundary b
1133  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
1134 
1135  // Create new element
1138  bulk_elem_pt,face_index);
1139 
1140  // Add to mesh
1142 
1143 #ifdef STRUCTURED_MESH
1144 
1145  // Set the GeomObject that defines the boundary shape and
1146  // specify which bulk boundary we are attached to (needed to extract
1147  // the boundary coordinate from the bulk nodes)
1150 
1151 #else
1152 
1153  switch(b)
1154  {
1155  case Contact_boundary_id:
1156 
1157  // Set the GeomObject that defines the boundary shape and
1158  // specify which bulk boundary we are attached to (needed to extract
1159  // the boundary coordinate from the bulk nodes)
1162 
1163  {
1164  unsigned nnod=el_pt->nnode();
1165  for (unsigned j=0;j<nnod;j++)
1166  {
1167  Node* nod_pt=el_pt->node_pt(j);
1168  if (nod_pt->is_on_boundary(Left_top_boundary_id))
1169  {
1170  left_contact_node_pt=nod_pt;
1171  }
1172  else if (nod_pt->is_on_boundary(Right_top_boundary_id))
1173  {
1174  right_contact_node_pt=nod_pt;
1175  }
1176  }
1177  }
1178 
1179  break;
1180 
1181  case Left_top_boundary_id:
1182 
1183  // Set the GeomObject that defines the boundary shape and
1184  // specify which bulk boundary we are attached to (needed to extract
1185  // the boundary coordinate from the bulk nodes)
1188 
1189  {
1190  unsigned nnod=el_pt->nnode();
1191  for (unsigned j=0;j<nnod;j++)
1192  {
1193  Node* nod_pt=el_pt->node_pt(j);
1194  if (nod_pt->is_on_boundary(Left_boundary_id))
1195  {
1196  left_left_top_node_pt=nod_pt;
1197  }
1198  else if (nod_pt->is_on_boundary(Contact_boundary_id))
1199  {
1200  right_left_top_node_pt=nod_pt;
1201  }
1202  }
1203  }
1204 
1205  break;
1206 
1207  case Right_top_boundary_id:
1208 
1209  // Set the GeomObject that defines the boundary shape and
1210  // specify which bulk boundary we are attached to (needed to extract
1211  // the boundary coordinate from the bulk nodes)
1214 
1215  {
1216  unsigned nnod=el_pt->nnode();
1217  for (unsigned j=0;j<nnod;j++)
1218  {
1219  Node* nod_pt=el_pt->node_pt(j);
1220  if (nod_pt->is_on_boundary(Right_boundary_id))
1221  {
1222  right_right_top_node_pt=nod_pt;
1223  }
1224  else if (nod_pt->is_on_boundary(Contact_boundary_id))
1225  {
1226  left_right_top_node_pt=nod_pt;
1227  }
1228  }
1229  }
1230  break;
1231 
1232  default:
1233 
1234  // Never get here...
1235  oomph_info << "Never get here! b = " << b << "\n";
1236  abort();
1237  }
1238 
1239 
1240 #endif
1241 
1242  // Loop over the nodes
1243  unsigned nnod=el_pt->nnode();
1244  for (unsigned j=0;j<nnod;j++)
1245  {
1246  Node* nod_pt = el_pt->node_pt(j);
1247 
1248  // Is the node also on side boundaries?
1249  if ((nod_pt->is_on_boundary(Left_boundary_id))||
1250  (nod_pt->is_on_boundary(Right_boundary_id)))
1251  {
1252  // How many nodal values were used by the "bulk" element
1253  // that originally created this node?
1254  unsigned n_bulk_value=el_pt->nbulk_value(j);
1255 
1256  // The remaining ones are Lagrange multipliers and we pin them.
1257  unsigned nval=nod_pt->nvalue();
1258  for (unsigned j=n_bulk_value;j<nval;j++)
1259  {
1260  nod_pt->pin(j);
1261  }
1262  }
1263  }
1264  }
1265  }
1266 
1267 #ifndef STRUCTURED_MESH
1268 
1269  Vector<double> zeta_left(1);
1270  Vector<double> zeta_right(1);
1271  left_left_top_node_pt->get_coordinates_on_boundary(Left_top_boundary_id,
1272  zeta_left);
1273  right_left_top_node_pt->get_coordinates_on_boundary(Left_top_boundary_id,
1274  zeta_right);
1275  if (zeta_left[0]>zeta_right[0])
1276  {
1278  oomph_info << "left top is reversed\n";
1279  }
1280  else
1281  {
1283  oomph_info << "left top is not reversed\n";
1284  }
1285 
1286  left_contact_node_pt->get_coordinates_on_boundary(Contact_boundary_id,
1287  zeta_left);
1288  right_contact_node_pt->get_coordinates_on_boundary(Contact_boundary_id,
1289  zeta_right);
1290  if (zeta_left[0]>zeta_right[0])
1291  {
1293  oomph_info << "contact is reversed\n";
1294  }
1295  else
1296  {
1298  oomph_info << "contact is not reversed\n";
1299  }
1300 
1301  left_right_top_node_pt->get_coordinates_on_boundary(Right_top_boundary_id,
1302  zeta_left);
1303  right_right_top_node_pt->get_coordinates_on_boundary(Right_top_boundary_id,
1304  zeta_right);
1305  if (zeta_left[0]>zeta_right[0])
1306  {
1308  oomph_info << "right top is reversed\n";
1309  }
1310  else
1311  {
1313  oomph_info << "right top is not reversed\n";
1314  }
1315 
1316 #endif
1317 
1318 
1319  } // end of create_displ_imposition_elements

References b, Global_Physical_Variables::Boundary_geom_object, ProblemParameters::Boundary_geom_object_contact, ProblemParameters::Boundary_geom_object_left, ProblemParameters::Boundary_geom_object_right, e(), oomph::Node::get_coordinates_on_boundary(), oomph::Node::is_on_boundary(), j, nb, oomph::FaceElement::nbulk_value(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::Data::nvalue(), oomph::oomph_info, oomph::Data::pin(), oomph::ImposeDisplacementByLagrangeMultiplierElement< ELEMENT >::set_boundary_shape_geom_object_pt(), WarpedLine::set_non_reversed(), and WarpedLine::set_reversed().

◆ create_imposed_heat_flux_elements_on_boulder()

template<class ELEMENT >
void ContactProblem< ELEMENT >::create_imposed_heat_flux_elements_on_boulder ( )
inlineprivate

Create imposed heat flux elements on surface of boulder.

3403  {
3404  // How many bulk elements are adjacent to boundary b?
3405  unsigned b=Boulder_top_boundary_id;
3406  unsigned n_element = Boulder_mesh_pt->nboundary_element(b);
3407 
3408  // Loop over the bulk elements adjacent to boundary b?
3409  for(unsigned e=0;e<n_element;e++)
3410  {
3411  // Get pointer to the bulk element that is adjacent to boundary b
3415 
3416  //What is the face index of element e along boundary b
3417  int face_index = Boulder_mesh_pt->face_index_at_boundary(b,e);
3418 
3419  // Build the corresponding flux element
3422  flux_element_pt = new
3425  bulk_elem_pt,face_index);
3426 
3427  // Pass pointer to heat flux
3428  flux_element_pt->flux_fct_pt() = &ProblemParameters::unit_flux;
3429 
3430  //Add the contact element to the surface mesh
3432 
3433  } //end of loop over bulk elements adjacent to boundary b
3434  }
Definition: unsteady_heat_flux_elements.h:56
void unit_flux(const double &time, const Vector< double > &x, double &flux)
hierher temp flux.
Definition: heated_linear_solid_contact_with_gravity.cc:2776

References oomph::Mesh::add_element_pt(), b, ContactProblem< ELEMENT >::Boulder_mesh_pt, ContactProblem< ELEMENT >::Boulder_surface_heat_flux_mesh_pt, ContactProblem< ELEMENT >::Boulder_top_boundary_id, oomph::Mesh::boundary_element_pt(), e(), oomph::Mesh::face_index_at_boundary(), oomph::Mesh::nboundary_element(), and ProblemParameters::unit_flux().

Referenced by ContactProblem< ELEMENT >::actions_after_adapt().

◆ delete_contact_elements() [1/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::delete_contact_elements ( )
inlineprivate

Delete contact elements.

3471  {
3472  // How many surface elements are in the surface mesh
3473  unsigned n_element = Surface_contact_mesh_pt->nelement();
3474 
3475  // Loop over the surface elements
3476  for(unsigned e=0;e<n_element;e++)
3477  {
3478  // Kill surface element
3480  }
3481 
3482  // Wipe the mesh
3484 
3485 
3486  // How many surface elements are in the surface mesh
3488 
3489  // Loop over the surface elements
3490  for(unsigned e=0;e<n_element;e++)
3491  {
3492  // Kill surface element
3494  }
3495 
3496  // Wipe the mesh
3498 
3499 
3500  // How many surface elements are in the surface mesh
3502 
3503  // Loop over the surface elements
3504  for(unsigned e=0;e<n_element;e++)
3505  {
3506  // Kill surface element
3508  }
3509 
3510  // Wipe the mesh
3512 
3513  }

References ContactProblem< ELEMENT >::Boulder_surface_contact_mesh_pt, ContactProblem< ELEMENT >::Boulder_surface_heat_flux_mesh_pt, e(), oomph::Mesh::element_pt(), oomph::Mesh::flush_element_and_node_storage(), oomph::Mesh::nelement(), and ContactProblem< ELEMENT >::Surface_contact_mesh_pt.

Referenced by ContactProblem< ELEMENT >::actions_before_adapt().

◆ delete_contact_elements() [2/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::delete_contact_elements ( )
inlineprivate

Delete contact elements.

1374  {
1375  // How many surface elements are in the surface mesh
1376  unsigned n_element = Surface_contact_mesh_pt->nelement();
1377 
1378  // Loop over the surface elements
1379  for(unsigned e=0;e<n_element;e++)
1380  {
1381  // Kill surface element
1383  }
1384 
1385  // Wipe the mesh
1387  }

References e().

◆ delete_contact_elements() [3/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::delete_contact_elements ( )
inlineprivate

Delete contact elements.

246  {
247  // How many surface elements are in the surface mesh
248  unsigned n_element = Surface_contact_mesh_pt->nelement();
249 
250  // Loop over the surface elements
251  for(unsigned e=0;e<n_element;e++)
252  {
253  // Kill surface element
255  }
256 
257  // Wipe the mesh
259  }

References e().

◆ delete_contact_elements() [4/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::delete_contact_elements ( )
inlineprivate

Delete contact elements.

886  {
887  // How many surface elements are in the surface mesh
888  unsigned n_element = Surface_contact_mesh_pt->nelement();
889 
890  // Loop over the surface elements
891  for(unsigned e=0;e<n_element;e++)
892  {
893  // Kill surface element
895  }
896 
897  // Wipe the mesh
899  }

References e().

◆ delete_displ_imposition_elements() [1/2]

template<class ELEMENT >
void ContactProblem< ELEMENT >::delete_displ_imposition_elements ( )
inlineprivate

Delete elements that enforce prescribed boundary motion by Lagrange multiplliers

3350  {
3351  // How many surface elements are in the surface mesh
3352  unsigned n_element = Displ_imposition_mesh_pt->nelement();
3353 
3354  // Loop over the surface elements
3355  for(unsigned e=0;e<n_element;e++)
3356  {
3357  // Kill surface element
3359  }
3360 
3361  // Wipe the mesh
3363  }

References ContactProblem< ELEMENT >::Displ_imposition_mesh_pt, e(), oomph::Mesh::element_pt(), oomph::Mesh::flush_element_and_node_storage(), and oomph::Mesh::nelement().

Referenced by ContactProblem< ELEMENT >::actions_before_adapt().

◆ delete_displ_imposition_elements() [2/2]

template<class ELEMENT >
void ContactProblem< ELEMENT >::delete_displ_imposition_elements ( )
inlineprivate

Delete elements that enforce prescribed boundary motion by Lagrange multiplliers

1325  {
1326  // How many surface elements are in the surface mesh
1327  unsigned n_element = Displ_imposition_mesh_pt->nelement();
1328 
1329  // Loop over the surface elements
1330  for(unsigned e=0;e<n_element;e++)
1331  {
1332  // Kill surface element
1334  }
1335 
1336  // Wipe the mesh
1338  }

References e().

◆ doc_solution() [1/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::doc_solution

Doc the solution.

4376 {
4377 
4378  oomph_info << "Outputting for step: " << Doc_info.number() << std::endl;
4379  ofstream some_file;
4380  char filename[100];
4381 
4382 
4383  // Write restart file
4384 /* sprintf(filename,"%s/restart%i.dat",doc_info.directory().c_str(),
4385  doc_info.number());
4386  some_file.open(filename);
4387  dump(some_file);
4388  some_file.close();*/
4389 
4390  // Number of plot points
4391  unsigned npts;
4392  npts=3;
4393 
4394  // Output solution
4395  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
4396  Doc_info.number());
4397  some_file.open(filename);
4398  Bulk_mesh_pt->output(some_file,npts);
4399  some_file.close();
4400 
4401  // Output solution coarsely (only element vertices for easier
4402  // mesh visualisation)
4403  sprintf(filename,"%s/coarse_soln%i.dat",Doc_info.directory().c_str(),
4404  Doc_info.number());
4405  some_file.open(filename);
4406  Bulk_mesh_pt->output(some_file,2);
4407  some_file.close();
4408 
4409  // Output solution
4410  sprintf(filename,"%s/boulder_soln%i.dat",Doc_info.directory().c_str(),
4411  Doc_info.number());
4412  some_file.open(filename);
4413  Boulder_mesh_pt->output(some_file,npts);
4414  some_file.close();
4415 
4416  // Output contact region on boulder
4417  sprintf(filename,"%s/boulder_contact%i.dat",Doc_info.directory().c_str(),
4418  Doc_info.number());
4419  some_file.open(filename);
4420  Boulder_surface_contact_mesh_pt->output(some_file,npts);
4421  some_file.close();
4422 
4423  // Output contact region on boulder
4424  sprintf(filename,"%s/boulder_heat_flux%i.dat",Doc_info.directory().c_str(),
4425  Doc_info.number());
4426  some_file.open(filename);
4427  Boulder_surface_heat_flux_mesh_pt->output(some_file,npts);
4428  some_file.close();
4429 
4430  // Output solution coarsely (only element vertices for easier
4431  // mesh visualisation)
4432  sprintf(filename,"%s/boulder_coarse_soln%i.dat",Doc_info.directory().c_str(),
4433  Doc_info.number());
4434  some_file.open(filename);
4435  Boulder_mesh_pt->output(some_file,2);
4436  some_file.close();
4437 
4438  // Output contact elements
4439  sprintf(filename,"%s/imposed_displ%i.dat",Doc_info.directory().c_str(),
4440  Doc_info.number());
4441  some_file.open(filename);
4442  Displ_imposition_mesh_pt->output(some_file);
4443  some_file.close();
4444 
4445  // Output contact elements
4446  sprintf(filename,"%s/contact%i.dat",Doc_info.directory().c_str(),
4447  Doc_info.number());
4448  some_file.open(filename);
4449  unsigned nel=Surface_contact_mesh_pt->nelement();
4450  for (unsigned e=0;e<nel;e++)
4451  {
4453  Surface_contact_mesh_pt->element_pt(e))->output(some_file,3);
4454  }
4455  some_file.close();
4456 
4457 
4458 
4459  // Output integration points of contact elements
4460  sprintf(filename,"%s/contact_integration_points%i.dat",
4461  Doc_info.directory().c_str(),
4462  Doc_info.number());
4463  some_file.open(filename);
4464  Vector<double> s(1);
4465  Vector<double> x(2);
4467  for (unsigned e=0;e<nel;e++)
4468  {
4472  unsigned nint=el_pt->integral_pt()->nweight();
4473  for (unsigned j=0;j<nint;j++)
4474  {
4475  s[0]=el_pt->integral_pt()->knot(j,0);
4476  el_pt->interpolated_x(s,x);
4477  some_file << x[0] << " " << x[1] << std::endl;
4478  }
4479  }
4480  some_file.close();
4481 
4482 
4483 
4484  // Output penetrator
4485  sprintf(filename,"%s/penetrator%i.dat",Doc_info.directory().c_str(),
4486  Doc_info.number());
4487  some_file.open(filename);
4488  unsigned n=500;
4490  some_file.close();
4491 
4492  // Output contact elements and assemble total resulting force
4493  Vector<double> total_contact_force(2,0.0);
4494  Vector<double> contact_force(2,0.0);
4495  sprintf(filename,"%s/contact%i.dat",Doc_info.directory().c_str(),
4496  Doc_info.number());
4497  some_file.open(filename);
4499  for (unsigned e=0;e<nel;e++)
4500  {
4504  el_pt->output(some_file,3);
4505  el_pt->resulting_contact_force(contact_force);
4506  total_contact_force[0]+=contact_force[0];
4507  total_contact_force[1]+=contact_force[1];
4508  }
4509  some_file.close();
4510 
4511  double radius_of_elastic_body=0.0;
4512 #ifdef STRUCTURED_MESH
4513  radius_of_elastic_body=
4515 #else
4516  radius_of_elastic_body=
4518 #endif
4519 
4520  // Get half-width of Hertz contact region
4521  double b_hertz=0.0;
4522  if (std::fabs(total_contact_force[1])>1.0e-16)
4523  {
4526  (-1.0/radius_of_elastic_body+
4528  (-total_contact_force[1]));
4529  }
4530  oomph_info << "b_hertz " << b_hertz << std::endl;
4531 
4532  double p_max_hertz = 0.0;
4533  if(b_hertz!=0.0)
4534  {
4535  p_max_hertz=2.0*total_contact_force[1]/
4536  (MathematicalConstants::Pi*b_hertz);
4537  }
4538 
4539  // Output Hertzian pressure contact distribution
4540  sprintf(filename,"%s/hertz%i.dat",Doc_info.directory().c_str(),
4541  Doc_info.number());
4542  some_file.open(filename);
4543  n=500;
4545  dynamic_cast<HeatedCircularPenetratorElement*>(
4547  Vector<double> centre(2);
4548  if (pen_el_pt!=0)
4549  {
4550  Vector<double> my_centre(pen_el_pt->centre());
4551  centre[0]=my_centre[0];
4552  centre[1]=my_centre[1];
4553  }
4554  else
4555  {
4556  centre[0]=ProblemParameters::Centre[0];
4557  centre[1]=ProblemParameters::Centre[1];
4558  }
4559  double x_c=centre[0];
4560  double width=2.0*b_hertz;
4561  for (unsigned j=0;j<n;j++)
4562  {
4563  double x=x_c-0.5*width+width*double(j)/double(n-1);
4564  double p=0.0;
4565  if (abs((x-x_c))<b_hertz)
4566  {
4567  p=p_max_hertz*sqrt(1.0-pow((x-x_c)/b_hertz,2));
4568  }
4569  some_file << x << " 0.0 " << p << std::endl;
4570  }
4571  some_file.close();
4572 
4573 
4574  double target_weight=0.0;
4575  if (pen_el_pt!=0)
4576  {
4577  target_weight=pen_el_pt->target_weight();
4578  }
4579  // double angle=0.0;
4580  // if (pen_el_pt!=0)
4581  // {
4582  // angle=pen_el_pt->angle();
4583  // }
4584 
4585 
4586  // Write trace file
4587  Trace_file << total_contact_force[0] << " "
4588  << total_contact_force[1] << " "
4589  << centre[0] << " "
4590  << centre[1] << " "
4591  << target_weight << " "
4592  //ALH/hierher Suppressed so that validation tests pass
4593  //For small changes in loading a different contact point can be selected
4594  //which means that different angles are computed
4595  // << angle << " "
4596  << std::endl;
4597 
4598 
4599  //Increment counter for solutions
4600  Doc_info.number()++;
4601 
4602 } // end of doc_solution
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
float * p
Definition: Tutorial_Map_using.cpp:9
double target_weight()
Target weight (returns zero if not imposed)
Definition: heated_linear_solid_contact_with_gravity.cc:2411
Vector< double > centre() const
Get centre of penetrator.
Definition: heated_linear_solid_contact_with_gravity.cc:2592
double radius() const
Return radius.
Definition: heated_linear_solid_contact_with_gravity.cc:2052
std::string directory() const
Output directory.
Definition: oomph_utilities.h:524
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Definition: elements.h:4528
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
Definition: heated_linear_solid_contact_with_gravity.cc:736
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
void resulting_contact_force(Vector< double > &contact_force)
Resulting contact force.
Definition: contact_elements.h:2371
virtual void output(std::ostream &outfile, const unsigned &nplot) const =0
Output coordinates of penetrator at nplot plot points.
RealScalar s
Definition: level1_cplx_impl.h:130
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
double Nu
Poisson's ratio.
Definition: unstructured_two_d_curved.cc:65
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.
Definition: overloaded_element_body.h:490

References abs(), Global_Physical_Variables::Boundary_geom_object, ProblemParameters::Boundary_geom_object_contact, ProblemParameters::Centre, HeatedCircularPenetratorElement::centre(), oomph::DocInfo::directory(), GlobalParameters::Doc_info, e(), boost::multiprecision::fabs(), MergeRestartFiles::filename, oomph::FiniteElement::integral_pt(), oomph::FaceElement::interpolated_x(), j, oomph::Integral::knot(), n, ProblemParameters::Nu, oomph::DocInfo::number(), oomph::Integral::nweight(), oomph::oomph_info, oomph::HeatedLinearSurfaceContactElement< ELEMENT >::output(), output(), oomph::Penetrator::output(), p, ProblemParameters::Penetrator_pt, BiharmonicTestFunctions2::Pi, Eigen::ArrayBase< Derived >::pow(), ProblemParameters::Radius, WarpedLine::radius(), oomph::LinearSurfaceContactElement< ELEMENT >::resulting_contact_force(), s, sqrt(), HeatedCircularPenetratorElement::target_weight(), oomph::Problem_Parameter::Trace_file, and plotDoE::x.

◆ doc_solution() [2/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::doc_solution ( )

Doc the solution.

◆ doc_solution() [3/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::doc_solution ( )

Doc the solution.

◆ doc_solution() [4/4]

template<class ELEMENT >
void ContactProblem< ELEMENT >::doc_solution ( )

Doc the solution.

◆ global_temporal_error_norm()

template<class ELEMENT >
double ContactProblem< ELEMENT >::global_temporal_error_norm ( )
inlinevirtual

Dummy global error norm for adaptive time-stepping.

Reimplemented from oomph::Problem.

210 {return 0.0;}

◆ switch_to_displ_control() [1/3]

template<class ELEMENT >
void ContactProblem< ELEMENT >::switch_to_displ_control ( )
inline

◆ switch_to_displ_control() [2/3]

template<class ELEMENT >
void ContactProblem< ELEMENT >::switch_to_displ_control ( )
inline

◆ switch_to_displ_control() [3/3]

template<class ELEMENT >
void ContactProblem< ELEMENT >::switch_to_displ_control ( )
inline

◆ switch_to_force_control() [1/3]

template<class ELEMENT >
void ContactProblem< ELEMENT >::switch_to_force_control ( )
inline

Switch to force control.

3080  {
3081  dynamic_cast<HeatedCircularPenetratorElement*>(
3082  ProblemParameters::Penetrator_pt)->impose_weight(
3084  dynamic_cast<HeatedCircularPenetratorElement*>(
3085  ProblemParameters::Penetrator_pt)->impose_horizontal_force(
3087  dynamic_cast<HeatedCircularPenetratorElement*>(
3088  ProblemParameters::Penetrator_pt)->set_equilibrium_target_forces();
3089 
3090  // Re-set contact mesh -- we now need to treat the positions
3091  // and penalty pressures of all nodes as external data of the
3092  // penetrator element!
3093  dynamic_cast<HeatedCircularPenetratorElement*>(
3094  ProblemParameters::Penetrator_pt)->set_contact_element_mesh_pt(
3096 
3097  // Reset penetrator because its equilibrium data (which is external
3098  // data for contact elements) has changed
3099  unsigned n_element=Surface_contact_mesh_pt->nelement();
3100  for(unsigned e=0;e<n_element;e++)
3101  {
3102  // Upcast from GeneralisedElement
3106 
3107  // Set pointer to penetrator
3109  }
3110 
3111  cout <<"New number of equations: " << assign_eqn_numbers() << std::endl;
3112 
3113  }

References oomph::Problem::assign_eqn_numbers(), e(), oomph::Mesh::element_pt(), ProblemParameters::Horizontal_force, oomph::Mesh::nelement(), ProblemParameters::Penetrator_pt, oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt(), ContactProblem< ELEMENT >::Surface_contact_mesh_pt, and ProblemParameters::Weight.

◆ switch_to_force_control() [2/3]

template<class ELEMENT >
void ContactProblem< ELEMENT >::switch_to_force_control ( )
inline

Switch to force control.

1057  {
1058  dynamic_cast<CircularPenetratorElement*>(
1059  ProblemParameters::Penetrator_pt)->impose_weight(
1061  dynamic_cast<CircularPenetratorElement*>(
1062  ProblemParameters::Penetrator_pt)->impose_horizontal_force(
1064  dynamic_cast<CircularPenetratorElement*>(
1065  ProblemParameters::Penetrator_pt)->set_equilibrium_target_forces();
1066 
1067  // Re-set contact mesh -- we now need to treat the positions
1068  // and penalty pressures of all nodes as external data of the
1069  // penetrator element!
1070  dynamic_cast<CircularPenetratorElement*>(
1071  ProblemParameters::Penetrator_pt)->set_contact_element_mesh_pt(
1073 
1074  // Reset penetrator because its equilibrium data (which is external
1075  // data for contact elements) has changed
1076  unsigned n_element=Surface_contact_mesh_pt->nelement();
1077  for(unsigned e=0;e<n_element;e++)
1078  {
1079  // Upcast from GeneralisedElement
1081  dynamic_cast<LinearSurfaceContactElement<ELEMENT>*>(
1083 
1084  // Set pointer to penetrator
1086  }
1087 
1088  cout <<"New number of equations: " << assign_eqn_numbers() << std::endl;
1089 
1090  }

References e(), ProblemParameters::Horizontal_force, ProblemParameters::Penetrator_pt, oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt(), and ProblemParameters::Weight.

◆ switch_to_force_control() [3/3]

template<class ELEMENT >
void ContactProblem< ELEMENT >::switch_to_force_control ( )
inline

Switch to force control.

813  {
814  dynamic_cast<CircularPenetratorElement*>(
815  ProblemParameters::Penetrator_pt)->impose_weight(
817  dynamic_cast<CircularPenetratorElement*>(
818  ProblemParameters::Penetrator_pt)->impose_horizontal_force(
820  dynamic_cast<CircularPenetratorElement*>(
821  ProblemParameters::Penetrator_pt)->set_equilibrium_target_forces();
822 
823  // Re-set contact mesh -- we now need to treat the positions
824  // and penalty pressures of all nodes as external data of the
825  // penetrator element!
826  dynamic_cast<CircularPenetratorElement*>(
827  ProblemParameters::Penetrator_pt)->set_contact_element_mesh_pt(
829 
830  // Reset penetrator because its equilibrium data (which is external
831  // data for contact elements) has changed
832  unsigned n_element=Surface_contact_mesh_pt->nelement();
833  for(unsigned e=0;e<n_element;e++)
834  {
835  // Upcast from GeneralisedElement
839 
840  // Set pointer to penetrator
842  }
843 
844  cout <<"New number of equations: " << assign_eqn_numbers() << std::endl;
845 
846  }

References e(), ProblemParameters::Horizontal_force, ProblemParameters::Penetrator_pt, oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt(), and ProblemParameters::Weight.

Member Data Documentation

◆ Backed_up_surface_contact_mesh_pt

template<class ELEMENT >
BackupMeshForProjection< TElement< 1, 3 > > * ContactProblem< ELEMENT >::Backed_up_surface_contact_mesh_pt
private

Backup of Surface_contact_mesh_pt so the Lagrange multipliers can be projected across

Referenced by ContactProblem< ELEMENT >::actions_after_adapt(), and ContactProblem< ELEMENT >::actions_before_adapt().

◆ Bottom_boundary_id

template<class ELEMENT >
unsigned ContactProblem< ELEMENT >::Bottom_boundary_id
private

ID of bottom boundary.

◆ Boulder_mesh_pt

◆ Boulder_surface_contact_mesh_pt

template<class ELEMENT >
Mesh* ContactProblem< ELEMENT >::Boulder_surface_contact_mesh_pt
private

Pointer to the "surface" mesh on boulder that applies flux boundary condition to boulder

Referenced by ContactProblem< ELEMENT >::complete_problem_setup(), ContactProblem< ELEMENT >::create_contact_heat_elements_on_boulder(), and ContactProblem< ELEMENT >::delete_contact_elements().

◆ Boulder_surface_heat_flux_mesh_pt

template<class ELEMENT >
Mesh* ContactProblem< ELEMENT >::Boulder_surface_heat_flux_mesh_pt
private

Mesh of elements that imposed unit heat flux on top of boulder.

Referenced by ContactProblem< ELEMENT >::create_imposed_heat_flux_elements_on_boulder(), and ContactProblem< ELEMENT >::delete_contact_elements().

◆ Bulk_mesh_pt

◆ Contact_boundary_id

template<class ELEMENT >
unsigned ContactProblem< ELEMENT >::Contact_boundary_id
private

ID of contact boundary.

◆ Contact_boundary_pt

template<class ELEMENT >
TriangleMeshPolyLine * ContactProblem< ELEMENT >::Contact_boundary_pt
private

Contact boundary in its poly line representation.

◆ Contact_id

template<class ELEMENT >
unsigned ContactProblem< ELEMENT >::Contact_id
private

ID of additional nodal values created by contact elements to store contact pressure/Lagr. mult.

Referenced by ContactProblem< ELEMENT >::complete_problem_setup(), and ContactProblem< ELEMENT >::create_contact_elements().

◆ Control_node_pt

template<class ELEMENT >
SolidNode * ContactProblem< ELEMENT >::Control_node_pt
private

Pointer to control node where Lagrange multiplier (contact pressure) is "pseudo-hijacked" to impose either displacement or weight constraint.

Referenced by ContactProblem< ELEMENT >::actions_before_adapt(), ContactProblem< ELEMENT >::complete_problem_setup(), and ContactProblem< ELEMENT >::switch_to_displ_control().

◆ Displ_imposition_mesh_pt

template<class ELEMENT >
SolidMesh * ContactProblem< ELEMENT >::Displ_imposition_mesh_pt
private

◆ Doc_info

template<class ELEMENT >
DocInfo ContactProblem< ELEMENT >::Doc_info
private

◆ Left_boundary_id

template<class ELEMENT >
unsigned ContactProblem< ELEMENT >::Left_boundary_id
private

ID of left boundary.

◆ Left_top_boundary_id

template<class ELEMENT >
unsigned ContactProblem< ELEMENT >::Left_top_boundary_id
private

ID of top left boundary.

◆ Maximum_element_length_on_boulder_boundary

template<class ELEMENT >
double ContactProblem< ELEMENT >::Maximum_element_length_on_boulder_boundary
private

◆ Maximum_element_length_on_contact_boundary

template<class ELEMENT >
double ContactProblem< ELEMENT >::Maximum_element_length_on_contact_boundary
private

Max. element length on contact boundary.

◆ Penetrator_mesh_pt

template<class ELEMENT >
Mesh * ContactProblem< ELEMENT >::Penetrator_mesh_pt
private

◆ Right_boundary_id

template<class ELEMENT >
unsigned ContactProblem< ELEMENT >::Right_boundary_id
private

ID of right boundary.

◆ Right_top_boundary_id

template<class ELEMENT >
unsigned ContactProblem< ELEMENT >::Right_top_boundary_id
private

ID of top right boundary.

◆ Surface_contact_mesh_pt

◆ Trace_file

template<class ELEMENT >
ofstream ContactProblem< ELEMENT >::Trace_file
private

Trace file.

◆ X_ll

template<class ELEMENT >
double ContactProblem< ELEMENT >::X_ll
private

x coordinate of lower left corner

◆ X_ur

template<class ELEMENT >
double ContactProblem< ELEMENT >::X_ur
private

x coordinate of upper right corner

◆ Xc_old

template<class ELEMENT >
double ContactProblem< ELEMENT >::Xc_old
private

◆ Y_ll

template<class ELEMENT >
double ContactProblem< ELEMENT >::Y_ll
private

y coordinate of lower left corner

◆ Y_ur

template<class ELEMENT >
double ContactProblem< ELEMENT >::Y_ur
private

y coordinate of upper right corner


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