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
 
OomphCommunicatorcommunicator_pt ()
 access function to the oomph-lib communicator More...
 
const OomphCommunicatorcommunicator_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 ()
 
LinearSolvermass_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...
 
Timetime_pt () const
 Return a pointer to the global time object (const version). More...
 
doubletime ()
 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 TimeSteppertime_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...
 
doubleminimum_dt ()
 Access function to min timestep in adaptive timestepping. More...
 
doublemaximum_dt ()
 Access function to max timestep in adaptive timestepping. More...
 
unsignedmax_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...
 
doublemax_residuals ()
 
booltime_adaptive_newton_crash_on_solve_fail ()
 Access function for Time_adaptive_newton_crash_on_solve_fail. More...
 
doublenewton_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...
 
doubleglobal_dof_pt (const unsigned &i)
 
doubledof (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...
 
doubledof_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...
 
booluse_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 Problemmake_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 ()
 
doublebifurcation_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 ()
 
intsign_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...
 
MeshSurface_contact_mesh_pt
 Pointer to the "surface" contact mesh. More...
 
SolidMeshDispl_imposition_mesh_pt
 
MeshBoulder_surface_heat_flux_mesh_pt
 Mesh of elements that imposed unit heat flux on top of boulder. More...
 
MeshPenetrator_mesh_pt
 Penetrator mesh. More...
 
MeshBoulder_surface_contact_mesh_pt
 
ofstream Trace_file
 Trace file. More...
 
DocInfo Doc_info
 
BackupMeshForProjection< TElement< 1, 3 > > * Backed_up_surface_contact_mesh_pt
 
SolidNodeControl_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...
 
TriangleMeshPolyLineContact_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)
 
doubledof_derivative (const unsigned &i)
 
doubledof_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 *, boolCalculate_dparameter_analytic
 
bool Calculate_hessian_products_analytic
 
LinearAlgebraDistributionDof_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< doubleMax_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< doubleDof_derivative
 Storage for the derivative of the problem variables wrt arc-length. More...
 
Vector< doubleDof_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
 
OomphCommunicatorCommunicator_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: