MeltContactProblem< ELEMENT > Class Template Reference

Problem class. More...

+ Inheritance diagram for MeltContactProblem< ELEMENT >:

Public Member Functions

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

Private Member Functions

void create_melt_elements ()
 Create melt elements. More...
 
void delete_melt_elements ()
 Delete melt elements. More...
 
void create_contact_elements ()
 Create contact elements. More...
 
void delete_contact_elements ()
 Delete contact elements. More...
 
void create_flux_elements ()
 Create flux elements. More...
 
void delete_flux_elements ()
 Delete flux elements. More...
 
void complete_problem_setup ()
 

Private Attributes

RefineableSolidTriangleMesh< ELEMENT > * Bulk_mesh_pt
 Pointer to bulk mesh. More...
 
Mesh * Surface_contact_mesh_pt
 Pointer to the "surface" mesh. More...
 
Mesh * Left_surface_melt_mesh_pt
 Pointer to the "surface" mesh. More...
 
Mesh * Surface_flux_mesh_pt
 Pointer to the "surface" mesh. More...
 
unsigned Left_melt_boundary_id
 ID of left melt boundary. More...
 
unsigned Right_melt_boundary_id
 ID of right melt boundary. More...
 
unsigned Lower_contact_boundary_id
 ID of lower contact boundary. More...
 
unsigned Upper_contact_boundary_id
 ID of upper 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...
 
ofstream Trace_file
 Trace file. More...
 
DocInfo Doc_info
 
BackupMeshForProjection< TElement< 1, 3 > > * Backed_up_left_surface_melt_mesh_pt
 

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_step ()
 
virtual void actions_after_newton_step ()
 
virtual void actions_after_implicit_timestep ()
 
virtual void actions_after_implicit_timestep_and_error_estimation ()
 
virtual void actions_before_explicit_timestep ()
 Actions that should be performed before each explicit time step. More...
 
virtual void actions_after_explicit_timestep ()
 Actions that should be performed after each explicit time step. More...
 
virtual void actions_before_read_unstructured_meshes ()
 
virtual void actions_after_read_unstructured_meshes ()
 
virtual void actions_after_change_in_global_parameter (double *const &parameter_pt)
 
virtual void actions_after_change_in_bifurcation_parameter ()
 
virtual void actions_after_parameter_increase (double *const &parameter_pt)
 
double & dof_derivative (const unsigned &i)
 
double & dof_current (const unsigned &i)
 
virtual void set_initial_condition ()
 
unsigned newton_solve_continuation (double *const &parameter_pt)
 
unsigned newton_solve_continuation (double *const &parameter_pt, DoubleVector &z)
 
void calculate_continuation_derivatives (double *const &parameter_pt)
 
void calculate_continuation_derivatives (const DoubleVector &z)
 
void calculate_continuation_derivatives_fd (double *const &parameter_pt)
 
bool does_pointer_correspond_to_problem_data (double *const &parameter_pt)
 
void set_consistent_pinned_values_for_continuation ()
 
- Protected Attributes inherited from oomph::Problem
Vector< Problem * > Copy_of_problem_pt
 
std::map< double *, bool > Calculate_dparameter_analytic
 
bool Calculate_hessian_products_analytic
 
LinearAlgebraDistribution * Dof_distribution_pt
 
Vector< double * > Dof_pt
 Vector of pointers to dofs. More...
 
DoubleVectorWithHaloEntries Element_count_per_dof
 
double Relaxation_factor
 
double Newton_solver_tolerance
 
unsigned Max_newton_iterations
 Maximum number of Newton iterations. More...
 
unsigned Nnewton_iter_taken
 
Vector< double > Max_res
 Maximum residuals at start and after each newton iteration. More...
 
double Max_residuals
 
bool Time_adaptive_newton_crash_on_solve_fail
 
bool Jacobian_reuse_is_enabled
 Is re-use of Jacobian in Newton iteration enabled? Default: false. More...
 
bool Jacobian_has_been_computed
 
bool Problem_is_nonlinear
 
bool Pause_at_end_of_sparse_assembly
 
bool Doc_time_in_distribute
 
unsigned Sparse_assembly_method
 
unsigned Sparse_assemble_with_arrays_initial_allocation
 
unsigned Sparse_assemble_with_arrays_allocation_increment
 
Vector< Vector< unsigned > > Sparse_assemble_with_arrays_previous_allocation
 
double Numerical_zero_for_sparse_assembly
 
double FD_step_used_in_get_hessian_vector_products
 
bool Mass_matrix_reuse_is_enabled
 
bool Mass_matrix_has_been_computed
 
bool Discontinuous_element_formulation
 
double Minimum_dt
 Minimum desired dt: if dt falls below this value, exit. More...
 
double Maximum_dt
 Maximum desired dt. More...
 
double DTSF_max_increase
 
double DTSF_min_decrease
 
double Minimum_dt_but_still_proceed
 
bool Scale_arc_length
 Boolean to control whether arc-length should be scaled. More...
 
double Desired_proportion_of_arc_length
 Proportion of the arc-length to taken by the parameter. More...
 
double Theta_squared
 
int Sign_of_jacobian
 Storage for the sign of the global Jacobian. More...
 
double Continuation_direction
 
double Parameter_derivative
 Storage for the derivative of the global parameter wrt arc-length. More...
 
double Parameter_current
 Storage for the present value of the global parameter. More...
 
bool Use_continuation_timestepper
 Boolean to control original or new storage of dof stuff. More...
 
unsigned Dof_derivative_offset
 
unsigned Dof_current_offset
 
Vector< double > Dof_derivative
 Storage for the derivative of the problem variables wrt arc-length. More...
 
Vector< double > Dof_current
 Storage for the present values of the variables. More...
 
double Ds_current
 Storage for the current step value. More...
 
unsigned Desired_newton_iterations_ds
 
double Minimum_ds
 Minimum desired value of arc-length. More...
 
bool Bifurcation_detection
 Boolean to control bifurcation detection via determinant of Jacobian. More...
 
bool Bisect_to_find_bifurcation
 Boolean to control wheter bisection is used to located bifurcation. More...
 
bool First_jacobian_sign_change
 Boolean to indicate whether a sign change has occured in the Jacobian. More...
 
bool Arc_length_step_taken
 Boolean to indicate whether an arc-length step has been taken. More...
 
bool Use_finite_differences_for_continuation_derivatives
 
OomphCommunicator * Communicator_pt
 The communicator for this problem. More...
 
bool Always_take_one_newton_step
 
double Timestep_reduction_factor_after_nonconvergence
 
bool Keep_temporal_error_below_tolerance
 
- Static Protected Attributes inherited from oomph::Problem
static ContinuationStorageScheme Continuation_time_stepper
 Storage for the single static continuation timestorage object. More...
 

Detailed Description

template<class ELEMENT>
class MeltContactProblem< ELEMENT >

Problem class.

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

Constructor & Destructor Documentation

◆ MeltContactProblem()

template<class ELEMENT >
MeltContactProblem< ELEMENT >::MeltContactProblem

Constructor.

Constructor for contact problem in square domain.

Identify penetrator as region 1

751 {
752 
753  // Output directory
754  Doc_info.set_directory("RESLT");
755 
756  // Output number
757  Doc_info.number()=0;
758 
759  // Open trace file
760  Trace_file.open("RESLT/trace.dat");
761 
762  // Allow for crap initial guess
763  Problem::Max_residuals=10000.0;
764 
765  // Allocate the timestepper -- this constructs the Problem's
766  // time object with a sufficient amount of storage to store the
767  // previous timsteps.
769 
770  // Pointer to the closed curve that defines the outer boundary
771  TriangleMeshClosedCurve* closed_curve_pt=0;
772 
773  // Build outer boundary as Polygon
774 
775  // The boundary, represented by polylines
776  Vector<TriangleMeshCurveSection*> boundary_polyline_pt;
777 
778  // Store coordinates of left and right point at which penetrator
779  // penetrates flat surface
780  Vector<double> right_penetrator(2);
781  Vector<double> left_penetrator(2);
782 
783  // Vertex coordinates on boundary
784  Vector<Vector<double> > bound_coords(2);
785 
786  // Left boundary
787  bound_coords[0].resize(2);
788  bound_coords[0][0]=0.0;
789  bound_coords[0][1]=1.0;
790 
791  bound_coords[1].resize(2);
792  bound_coords[1][0]=0.0;
793  bound_coords[1][1]=0.0;
794 
795  // Build the boundary polyline
797  boundary_polyline_pt.push_back(new TriangleMeshPolyLine(bound_coords,
799 
800  // Bottom boundary
801  bound_coords[0].resize(2);
802  bound_coords[0][0]=0.0;
803  bound_coords[0][1]=0.0;
804 
805  bound_coords[1].resize(2);
806  bound_coords[1][0]=1.0;
807  bound_coords[1][1]=0.0;
808 
809  // Build the boundary polyline
811  boundary_polyline_pt.push_back(new TriangleMeshPolyLine(bound_coords,
813 
814 
815  // Right boundary
816  bound_coords[0].resize(2);
817  bound_coords[0][0]=1.0;
818  bound_coords[0][1]=0.0;
819 
820  bound_coords[1].resize(2);
821  bound_coords[1][0]=1.0;
822  bound_coords[1][1]=1.0;
823 
824  // Build the boundary polyline
826  boundary_polyline_pt.push_back(new TriangleMeshPolyLine(bound_coords,
828 
829 
830  // Half opening angle of "boulder"
831  double half_opening_angle=acos((ProblemParameters::Centre[1]-1.0)/
833 
834 
835  // Points at which the penetrator meets the flat surface
836  right_penetrator[0]=ProblemParameters::Centre[0]+
837  ProblemParameters::Radius*sin(half_opening_angle);
838  right_penetrator[1]=1.0;
839 
840  left_penetrator[0]=ProblemParameters::Centre[0]-
841  ProblemParameters::Radius*sin(half_opening_angle);
842  left_penetrator[1]=1.0;
843 
844  // Right melt boundary
845  bound_coords[0].resize(2);
846  bound_coords[0][0]=1.0;
847  bound_coords[0][1]=1.0;
848 
849  bound_coords[1].resize(2);
850  bound_coords[1][0]=right_penetrator[0];
851  bound_coords[1][1]=right_penetrator[1];
852 
853 
854  // Build boundary poly line
856  boundary_polyline_pt.push_back
857  (new TriangleMeshPolyLine(bound_coords,Right_melt_boundary_id));
858 
859 
860 
861  // Vertex coordinates on lower contact boundary
862  unsigned nvertex_contact=20;
863  Vector<Vector<double> > circle_bound_coords(nvertex_contact);
864 
865  // First point
866  circle_bound_coords[0].resize(2);
867  circle_bound_coords[0][0]=right_penetrator[0];
868  circle_bound_coords[0][1]=right_penetrator[1];
869 
870  // Internal points
871  for (unsigned j=1;j<nvertex_contact-1;j++)
872  {
873  double phi=1.5*MathematicalConstants::Pi+half_opening_angle
874  -2.0*half_opening_angle*double(j)/double(nvertex_contact-1);
875 
876  circle_bound_coords[j].resize(2);
877  circle_bound_coords[j][0]=ProblemParameters::Centre[0]+
879 
880  circle_bound_coords[j][1]=ProblemParameters::Centre[1]+
882  }
883 
884  // Last point
885  circle_bound_coords[nvertex_contact-1].resize(2);
886  circle_bound_coords[nvertex_contact-1][0]=left_penetrator[0];
887  circle_bound_coords[nvertex_contact-1][1]=left_penetrator[1];
888 
889 
890  // Build boundary poly line
892  TriangleMeshPolyLine* lower_contact_boundary_pt=
893  new TriangleMeshPolyLine(circle_bound_coords,Lower_contact_boundary_id);
894 
895  // Store as internal poly line
896  Vector<TriangleMeshCurveSection*> internal_polyline_pt(1);
897  internal_polyline_pt[0]=lower_contact_boundary_pt;
898 
899  // Vertex coordinates on upper contact boundary
900  unsigned nvertex_contact_upper=unsigned(double(nvertex_contact)*
902  half_opening_angle);
903  circle_bound_coords.clear();
904  circle_bound_coords.resize(nvertex_contact_upper);
905 
906  // First point: recycle
907  circle_bound_coords[0].resize(2);
908  circle_bound_coords[0][0]=right_penetrator[0];
909  circle_bound_coords[0][1]=right_penetrator[1];
910  for (unsigned j=1;j<nvertex_contact_upper-1;j++)
911  {
912  double phi=-(0.5*MathematicalConstants::Pi-half_opening_angle)
913  +(2.0*MathematicalConstants::Pi-2.0*half_opening_angle)*
914  double(j)/double(nvertex_contact_upper-1);
915 
916  circle_bound_coords[j].resize(2);
917  circle_bound_coords[j][0]=ProblemParameters::Centre[0]+
919 
920  circle_bound_coords[j][1]=ProblemParameters::Centre[1]+
922  }
923 
924  // Last point
925  circle_bound_coords[nvertex_contact_upper-1].resize(2);
926  circle_bound_coords[nvertex_contact_upper-1][0]=left_penetrator[0];
927  circle_bound_coords[nvertex_contact_upper-1][1]=left_penetrator[1];
928 
929  // Build boundary poly line
931  TriangleMeshPolyLine* upper_contact_boundary_pt=
932  new TriangleMeshPolyLine(circle_bound_coords,Upper_contact_boundary_id);
933  boundary_polyline_pt.push_back(upper_contact_boundary_pt);
934 
935  // Left melt boundary
936  unsigned npt_left=3;
937  Vector<Vector<double> > left_melt_bound_coords(npt_left);
938  left_melt_bound_coords[0].resize(2);
939  left_melt_bound_coords[0][0]=left_penetrator[0];
940  left_melt_bound_coords[0][1]=left_penetrator[1];
941  for (unsigned j=1;j<npt_left-1;j++)
942  {
943  left_melt_bound_coords[j].resize(2);
944  left_melt_bound_coords[j][0]=left_penetrator[0]*double(j)/double(npt_left-1);
945  left_melt_bound_coords[j][1]=left_penetrator[1]-
946  (1.0-left_penetrator[1])*double(j)/double(npt_left-1);
947  }
948  left_melt_bound_coords[npt_left-1].resize(2);
949  left_melt_bound_coords[npt_left-1][0]=0.0;
950  left_melt_bound_coords[npt_left-1][1]=1.0;
951 
952 
953  // Build boundary poly line
955  TriangleMeshPolyLine* left_melt_boundary_pt=
956  new TriangleMeshPolyLine(left_melt_bound_coords,Left_melt_boundary_id);
957  boundary_polyline_pt.push_back(left_melt_boundary_pt);
958 
959  // Keep number of elements on boundary approx fixed.
960  double maximum_length=4.0*left_penetrator[0]/double(npt_left-1); //hierher why 4
961  left_melt_boundary_pt->set_maximum_length(maximum_length);
962 
963 
964  // Connect first vertex (vertex 0) of upper contact boundary
965  // with initial vertex of lower contact boundary
966  unsigned vertex_id_of_connection=0;
967  lower_contact_boundary_pt->connect_initial_vertex_to_polyline(
968  dynamic_cast<TriangleMeshPolyLine*>(upper_contact_boundary_pt),
969  vertex_id_of_connection);
970 
971 
972 //#define USE_BROKEN
973 #ifdef USE_BROKEN
974 
975  // hierher this misses out final node on lower contact boundary
976  // in boundary lookup scheme!
977 
978  // Connect last vertex of upper contact boundary
979  // with final vertex of lower contact boundary
980  vertex_id_of_connection=nvertex_contact_upper-1;
981  lower_contact_boundary_pt->connect_final_vertex_to_polyline(
982  dynamic_cast<TriangleMeshPolyLine*>(upper_contact_boundary_pt),
983  vertex_id_of_connection);
984 
985 #else
986 
987  // Connect first vertex left melt boundary
988  // with final vertex of lower contact boundary
989  vertex_id_of_connection=0;
990  lower_contact_boundary_pt->connect_final_vertex_to_polyline(
991  dynamic_cast<TriangleMeshPolyLine*>(left_melt_boundary_pt),
992  vertex_id_of_connection);
993 
994 #endif
995 
996  // Create open curve that defines boulder/ice interface
997  Vector<TriangleMeshOpenCurve*> inner_boundary_pt;
998  inner_boundary_pt.push_back(new TriangleMeshOpenCurve(internal_polyline_pt));
999 
1000 
1001  // Create the triangle mesh polygon for outer boundary
1002  //----------------------------------------------------
1003  TriangleMeshPolygon *outer_polygon =
1004  new TriangleMeshPolygon(boundary_polyline_pt);
1005 
1006  // Set the pointer
1007  closed_curve_pt = outer_polygon;
1008 
1009  // Now build the mesh
1010  //===================
1011 
1012  // Use the TriangleMeshParameters object for helping on the manage of the
1013  // TriangleMesh parameters
1014  TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
1015 
1016  // Specify the maximum area element
1017  double uniform_element_area=0.002;
1018  triangle_mesh_parameters.element_area() = uniform_element_area;
1019 
1020  // Specify the internal open boundary
1021  triangle_mesh_parameters.internal_open_curves_pt()=inner_boundary_pt;
1022 
1024  triangle_mesh_parameters.add_region_coordinates(1,ProblemParameters::Centre);
1025 
1026 
1027 #ifdef ADAPTIVE
1028 
1029  // Create the mesh
1030  Bulk_mesh_pt=new RefineableSolidTriangleMesh<ELEMENT>(triangle_mesh_parameters,
1031  time_stepper_pt());
1032 
1033  // Set error estimator for bulk mesh
1035  Bulk_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
1036 
1037  // // Set element size limits
1038  // Bulk_mesh_pt->max_element_size()=0.2;
1039  // Bulk_mesh_pt->min_element_size()=0.002;
1040 
1041 #else
1042 
1043  // Build mesh
1044  Bulk_mesh_pt=new SolidTriangleMesh<ELEMENT>(triangle_mesh_parameters,
1045  time_stepper_pt());
1046 
1047 #endif
1048 
1049  Bulk_mesh_pt->output("mesh.dat");
1050  Bulk_mesh_pt->output_boundaries("boundaries.dat");
1051 
1052  // Create the surface mesh as an empty mesh
1054 
1055  // Build 'em
1057 
1058  // Create the surface mesh as an empty mesh
1060 
1061  // Build 'em
1063 
1064  // Create the surface mesh as an empty mesh
1066 
1067  // Build 'em
1069 
1070  // Set boundary condition and complete the build of all elements
1072 
1073  // Add the sub meshes to the problem
1078 
1079  // Combine all submeshes into a single global Mesh
1081 
1082  // Set the initial conditions
1083  unsigned nnod = Bulk_mesh_pt->nnode();
1084  for(unsigned j=0;j<nnod;j++)
1085  {
1086  Bulk_mesh_pt->node_pt(j)->set_value(0,ExactSolution::U0);
1087  }
1088 
1089  // Do equation numbering
1090  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
1091 
1092 } // end of constructor
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
AnnoyingScalar acos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:138
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
unsigned Lower_contact_boundary_id
ID of lower contact boundary.
Definition: circular_boulder_melt.cc:718
unsigned Bottom_boundary_id
ID of bottom boundary.
Definition: circular_boulder_melt.cc:724
unsigned Left_boundary_id
ID of left boundary.
Definition: circular_boulder_melt.cc:727
Mesh * Left_surface_melt_mesh_pt
Pointer to the "surface" mesh.
Definition: circular_boulder_melt.cc:706
void complete_problem_setup()
Definition: circular_boulder_melt.cc:566
unsigned Upper_contact_boundary_id
ID of upper contact boundary.
Definition: circular_boulder_melt.cc:721
unsigned Right_boundary_id
ID of right boundary.
Definition: circular_boulder_melt.cc:730
void create_contact_elements()
Create contact elements.
Definition: circular_boulder_melt.cc:446
void create_flux_elements()
Create flux elements.
Definition: circular_boulder_melt.cc:507
void create_melt_elements()
Create melt elements.
Definition: circular_boulder_melt.cc:386
unsigned Left_melt_boundary_id
ID of left melt boundary.
Definition: circular_boulder_melt.cc:712
Mesh * Surface_flux_mesh_pt
Pointer to the "surface" mesh.
Definition: circular_boulder_melt.cc:709
RefineableSolidTriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to bulk mesh.
Definition: circular_boulder_melt.cc:693
ofstream Trace_file
Trace file.
Definition: circular_boulder_melt.cc:733
Mesh * Surface_contact_mesh_pt
Pointer to the "surface" mesh.
Definition: circular_boulder_melt.cc:703
unsigned Right_melt_boundary_id
ID of right melt boundary.
Definition: circular_boulder_melt.cc:715
DocInfo Doc_info
Definition: circular_boulder_melt.cc:736
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: mesh.h:67
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
Unstructured Triangle Mesh upgraded to solid mesh.
Definition: triangle_mesh.template.h:3896
Base class defining a closed curve for the Triangle mesh generation.
Definition: unstructured_two_d_mesh_geometry_base.h:1339
void connect_initial_vertex_to_polyline(TriangleMeshPolyLine *polyline_pt, const unsigned &vertex_number, const double &tolerance_for_connection=1.0e-14)
Definition: unstructured_two_d_mesh_geometry_base.cc:1120
void connect_final_vertex_to_polyline(TriangleMeshPolyLine *polyline_pt, const unsigned &vertex_number, const double &tolerance_for_connection=1.0e-14)
Definition: unstructured_two_d_mesh_geometry_base.cc:1189
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:1642
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
double Pi
Definition: two_d_biharmonic.cc:235
double U0
Constant/initial temperature.
Definition: circular_boulder_melt.cc:99
Z2ErrorEstimator * error_estimator_pt
Definition: MortaringCantileverCompareToNonMortaring.cpp:190
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
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References acos(), oomph::TriangleMeshParameters::add_region_coordinates(), ProblemParameters::Centre, oomph::TriangleMeshCurveSection::connect_final_vertex_to_polyline(), oomph::TriangleMeshCurveSection::connect_initial_vertex_to_polyline(), cos(), GlobalParameters::Doc_info, oomph::TriangleMeshParameters::element_area(), MeshRefinement::error_estimator_pt, oomph::TriangleMeshParameters::internal_open_curves_pt(), j, oomph::DocInfo::number(), BiharmonicTestFunctions2::Pi, ProblemParameters::Radius, oomph::DocInfo::set_directory(), oomph::TriangleMeshCurveSection::set_maximum_length(), sin(), oomph::Problem_Parameter::Trace_file, and ExactSolution::U0.

◆ ~MeltContactProblem()

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

Destructor (empty)

261 {}

Member Function Documentation

◆ actions_after_adapt()

template<class ELEMENT >
void MeltContactProblem< 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.

332  {
333  // Create contact elements
335 
336  // Create melt elements
338 
339  // Create flux elements
341 
342  // Rebuild the Problem's global mesh from its various sub-meshes
344 
345  // Rebuild elements
347 
348 // hierher
349  // // Now project from backup of original contact mesh to new one
350  // Backed_up_left_surface_melt_mesh_pt->project_onto_new_mesh(
351  // Left_surface_melt_mesh_pt);
352 
353  // // Kill backed up mesh
354  // delete Backed_up_left_surface_melt_mesh_pt;
355  // Backed_up_left_surface_melt_mesh_pt=0; // hierher add to constructor
356 
357  // Output flux with melt
358  ofstream some_file;
359  char filename[100];
360  sprintf(filename,"flux_with_melt_after.dat");
361  some_file.open(filename);
362  unsigned nel=Left_surface_melt_mesh_pt->nelement();
363  for (unsigned e=0;e<nel;e++)
364  {
365  dynamic_cast<SurfaceMeltElement<ELEMENT>* >(
366  Left_surface_melt_mesh_pt->element_pt(e))->output_melt(some_file);
367  }
368  some_file.close();
369 
370  oomph_info << "hierher doc_solution from after adapt\n";
371  doc_solution();
372 
373 
374  }
Array< double, 1, 3 > e(1./3., 0.5, 2.)
void doc_solution()
Doc the solution.
Definition: circular_boulder_melt.cc:1101
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
void rebuild_global_mesh()
Definition: problem.cc:1533
Definition: heat_transfer_and_melt_elements.h:2249
string filename
Definition: MergeRestartFiles.py:39
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

References e(), MergeRestartFiles::filename, and oomph::oomph_info.

◆ actions_after_newton_solve()

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

Update the problem specs after solve (empty)

Reimplemented from oomph::Problem.

264 {}

◆ actions_before_adapt()

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

Actions before adapt: wipe contact elements.

Reimplemented from oomph::Problem.

293  {
294 
295 // hierher
296  // // Make backup of surface mesh
297  // Backed_up_left_surface_melt_mesh_pt=
298  // new BackupMeshForProjection<TElement<1,3> >( //FaceGeometry<ELEMENT> >(
299  // Left_surface_melt_mesh_pt,Left_melt_boundary_id,
300  // Left_melt_boundary_id); //Note: we labeled the bc ID for face elements
301  // // by boundary id
302 
303 
304  // Output flux with melt
305  ofstream some_file;
306  char filename[100];
307  sprintf(filename,"flux_with_melt_before.dat");
308  some_file.open(filename);
309  unsigned nel=Left_surface_melt_mesh_pt->nelement();
310  for (unsigned e=0;e<nel;e++)
311  {
312  dynamic_cast<SurfaceMeltElement<ELEMENT>* >(
313  Left_surface_melt_mesh_pt->element_pt(e))->output_melt(some_file);
314  }
315  some_file.close();
316 
317 
318  // Kill the elements and wipe surface mesh
322 
323  // Rebuild the Problem's global mesh from its various sub-meshes
325  }
void delete_melt_elements()
Delete melt elements.
Definition: circular_boulder_melt.cc:427
void delete_flux_elements()
Delete flux elements.
Definition: circular_boulder_melt.cc:547
void delete_contact_elements()
Delete contact elements.
Definition: circular_boulder_melt.cc:490

References e(), and MergeRestartFiles::filename.

◆ actions_before_implicit_timestep()

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

Actions before next timestep.

Reimplemented from oomph::Problem.

272  {
273  // Amplitude of oscillation
274  double amplitude=0.2; //-0.2;
275 
277  (1.0+0.05*time_pt()->time())*
278  0.5*(1.0-cos(2.0*MathematicalConstants::Pi*time_pt()->time()));
279 
280  oomph_info << "Solving for y_c = "
281  << ProblemParameters::Centre[1] << std::endl;
282  }
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_initial
Initial y position of centre of penetrator.
Definition: circular_boulder_melt.cc:82

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

◆ actions_before_newton_convergence_check()

template<class ELEMENT >
void MeltContactProblem< ELEMENT >::actions_before_newton_convergence_check ( )
inlinevirtual

Newton convergence check – empty.

Reimplemented from oomph::Problem.

287  {
288  }

◆ actions_before_newton_solve()

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

Update the problem specs before solve (empty)

Reimplemented from oomph::Problem.

267 {}

◆ complete_problem_setup()

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

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

567  {
568 
569  // Set (pseudo-)solid mechanics properties for all elements
570  //---------------------------------------------------------
571  unsigned n_element = Bulk_mesh_pt->nelement();
572  for(unsigned e=0;e<n_element;e++)
573  {
574  //Cast to a solid element
575  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
576 
577  // Set the constitutive law
578  el_pt->constitutive_law_pt() =
580 
581  // Set density to zero
582  el_pt->lambda_sq_pt()=&ProblemParameters::Lambda_sq;
583 
584  // Disable inertia
585  el_pt->disable_inertia();
586 
587  // Set source function
588  el_pt->source_fct_pt() =
590  }
591 
592  // Apply boundary conditions for solid
593  //------------------------------------
594 
595  // Bottom: completely pinned
596  unsigned b=Bottom_boundary_id;
597  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
598  for (unsigned j=0;j<nnod;j++)
599  {
600  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
601  nod_pt->pin_position(0);
602  nod_pt->pin_position(1);
603  }
604 
605  // Sides: Symmetry bcs
607  nnod=Bulk_mesh_pt->nboundary_node(b);
608  for (unsigned j=0;j<nnod;j++)
609  {
610  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
611  nod_pt->pin_position(0);
612  }
614  nnod=Bulk_mesh_pt->nboundary_node(b);
615  for (unsigned j=0;j<nnod;j++)
616  {
617  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
618  nod_pt->pin_position(0);
619  }
620 
621  // Assign the Lagrangian coordinates -- sensible
622  // because we've completely rebuilt the mesh
623  // and haven't copied across any Lagrange multipliers
624  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
625  n_element=Left_surface_melt_mesh_pt->nelement();
626  for (unsigned e=0;e<n_element;e++)
627  {
629  dynamic_cast<SurfaceMeltElement<ELEMENT>*>(
632  }
633 
634  // Loop over the contact elements, pass pointer to penetrator and make sticky
635  //---------------------------------------------------------------------------
636  n_element=Surface_contact_mesh_pt->nelement();
637  for(unsigned e=0;e<n_element;e++)
638  {
639  // Upcast from GeneralisedElement
643 
644  // Set pointer to penetrator
646 
647  // Make it sticky to enforce permanent contact
648  el_pt->enable_stick();
649  }
650 
651 
652  // Loop over the flux elements to pass pointer to prescribed flux function
653  //------------------------------------------------------------------------
654  n_element=Left_surface_melt_mesh_pt->nelement();
655  for(unsigned e=0;e<n_element;e++)
656  {
657  // Upcast from GeneralisedElement to SurfaceMeltElement flux element
659  dynamic_cast<SurfaceMeltElement<ELEMENT>*>(
661 
662  // Pass validation flux
663  el_pt->flux_fct_pt()=
665 
666  // Set melt temperature
668  }
669 
670 
671 
672  // Loop over the flux elements to pass pointer to prescribed flux function
673  //------------------------------------------------------------------------
674  n_element=Surface_flux_mesh_pt->nelement();
675  for(unsigned e=0;e<n_element;e++)
676  {
677  // Upcast from GeneralisedElement
681 
682  // Pass validation flux
684  }
685 
686  }
Scalar * b
Definition: benchVecAdd.cpp:17
Definition: contact_elements.h:1086
Definition: nodes.h:1686
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
void enable_stick()
Definition: contact_elements.h:776
void set_lagrange_multiplier_pressure_to_zero()
Definition: heat_transfer_and_melt_elements.h:2433
double *& melt_temperature_pt()
Pointer to (non-default) melt temperature.
Definition: heat_transfer_and_melt_elements.h:2402
UnsteadyHeatPrescribedFluxFctPt & flux_fct_pt()
Access function for the prescribed-flux function pointer.
Definition: heat_transfer_and_melt_elements.h:147
Definition: heat_transfer_and_melt_elements.h:1522
void prescribed_flux_for_unsteady_heat_validation(const double &t, const Vector< double > &x, const Vector< double > &n, const double &u, double &flux)
Definition: circular_boulder_melt.cc:219
void flux_into_bulk(const double &t, const Vector< double > &x, const Vector< double > &n, const double &u, double &flux)
Definition: circular_boulder_melt.cc:185
void get_source_for_unsteady_heat_validation(const double &t, const Vector< double > &x, double &source)
Source function to make it an exact solution.
Definition: circular_boulder_melt.cc:133
Penetrator * Penetrator_pt
Penetrator.
Definition: circular_boulder_melt.cc:88
double Melt_temperature
Melt-temperature.
Definition: circular_boulder_melt.cc:67
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

References b, ProblemParameters::Constitutive_law_pt, e(), oomph::SurfaceContactElementBase< ELEMENT >::enable_stick(), oomph::TemplateFreeUnsteadyHeatBaseFaceElement::flux_fct_pt(), ExactSolution::flux_into_bulk(), ExactSolution::get_source_for_unsteady_heat_validation(), j, ProblemParameters::Lambda_sq, ProblemParameters::Melt_temperature, oomph::SurfaceMeltElement< ELEMENT >::melt_temperature_pt(), ProblemParameters::Penetrator_pt, oomph::SolidNode::pin_position(), ExactSolution::prescribed_flux_for_unsteady_heat_validation(), oomph::SurfaceMeltElement< ELEMENT >::set_lagrange_multiplier_pressure_to_zero(), and oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt().

◆ create_contact_elements()

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

Create contact elements.

447  {
448  oomph_info << "not doinig contact hierher\n";
449  return;
450 
451  // Which boundaries are in contact with penetrator?
452  Vector<unsigned> contact_boundaries;
453  contact_boundaries.push_back(Upper_contact_boundary_id);
454  contact_boundaries.push_back(Lower_contact_boundary_id);
455  unsigned n=contact_boundaries.size();
456  for (unsigned bb=0;bb<n;bb++)
457  {
458  // How many bulk elements are adjacent to boundary b?
459  unsigned b=contact_boundaries[bb];
460  unsigned n_element = Bulk_mesh_pt->nboundary_element_in_region(b,1);
461 
462  // Loop over the bulk elements adjacent to boundary b?
463  for(unsigned e=0;e<n_element;e++)
464  {
465  // Get pointer to the bulk element that is adjacent to boundary b
466  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
467  Bulk_mesh_pt->boundary_element_in_region_pt(b,1,e));
468 
469  //What is the face index of element e along boundary b
470  int face_index = Bulk_mesh_pt->face_index_at_boundary_in_region(b,1,e);
471 
472  // Id of additional dofs created by this element
473  // Use upper boundary ID for both
474  unsigned id=Upper_contact_boundary_id;
475 
476  // Build the corresponding contact element
477  NonlinearSurfaceContactElement<ELEMENT>* contact_element_pt = new
478  NonlinearSurfaceContactElement<ELEMENT>(bulk_elem_pt,face_index,id);
479 
480  //Add the contact element to the surface mesh
481  Surface_contact_mesh_pt->add_element_pt(contact_element_pt);
482 
483  } //end of loop over bulk elements adjacent to boundary b
484  }
485  }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:617

References b, e(), n, and oomph::oomph_info.

◆ create_flux_elements()

template<class ELEMENT >
void MeltContactProblem< ELEMENT >::create_flux_elements ( )
inlineprivate

Create flux elements.

508  {
509 
510  oomph_info << "not doinig flux hierher\n";
511  return;
512 
513  // Which boundaries are in contact with penetrator?
514  Vector<unsigned> flux_boundaries;
515  flux_boundaries.push_back(Upper_contact_boundary_id);
516  unsigned n=flux_boundaries.size();
517  for (unsigned bb=0;bb<n;bb++)
518  {
519  // How many bulk elements are adjacent to boundary b?
520  unsigned b=flux_boundaries[bb];
521  unsigned n_element = Bulk_mesh_pt->nboundary_element_in_region(b,1);
522 
523  // Loop over the bulk elements adjacent to boundary b?
524  for(unsigned e=0;e<n_element;e++)
525  {
526  // Get pointer to the bulk element that is adjacent to boundary b
527  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
528  Bulk_mesh_pt->boundary_element_in_region_pt(b,1,e));
529 
530  //What is the face index of element e along boundary b
531  int face_index = Bulk_mesh_pt->face_index_at_boundary_in_region(b,1,e);
532 
533  // Build the corresponding prescribed-flux element
534  UnsteadyHeatBaseFaceElement<ELEMENT>* flux_element_pt = new
535  UnsteadyHeatBaseFaceElement<ELEMENT>(bulk_elem_pt,face_index);
536 
537  // Add to mesh
538  Surface_flux_mesh_pt->add_element_pt(flux_element_pt);
539 
540  } //end of loop over bulk elements adjacent to boundary b
541  }
542  }

References b, e(), n, and oomph::oomph_info.

◆ create_melt_elements()

template<class ELEMENT >
void MeltContactProblem< ELEMENT >::create_melt_elements ( )
inlineprivate

Create melt elements.

387  {
388  Vector<unsigned> melt_boundary_id;
389  melt_boundary_id.push_back(Left_melt_boundary_id);
390  // hierher
391  melt_boundary_id.push_back(Right_melt_boundary_id);
392  //melt_boundary_id.push_back(Upper_contact_boundary_id);
393 
394  unsigned nb=melt_boundary_id.size();
395  for (unsigned bb=0;bb<nb;bb++)
396  {
397  unsigned b=melt_boundary_id[bb];
398  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
399 
400  // Loop over the bulk elements adjacent to boundary b?
401  for(unsigned e=0;e<n_element;e++)
402  {
403  // Get pointer to the bulk element that is adjacent to boundary b
404  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
405  Bulk_mesh_pt->boundary_element_pt(b,e));
406 
407  //What is the face index of element e along boundary b
408  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
409 
410  // ID of additional dofs created by this element use boundary ID
411  unsigned id=b;
412 
413  // Build the corresponding prescribed-flux element
414  SurfaceMeltElement<ELEMENT>* flux_element_pt = new
415  SurfaceMeltElement<ELEMENT>(bulk_elem_pt,face_index,id);
416 
417  //Add the prescribed-flux element to the surface mesh
418  Left_surface_melt_mesh_pt->add_element_pt(flux_element_pt);
419 
420  } //end of loop over bulk elements adjacent to boundary b
421  }
422  }
int nb
Definition: level2_impl.h:286

References b, e(), and nb.

◆ delete_contact_elements()

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

Delete contact elements.

491  {
492  // How many surface elements are in the surface mesh
493  unsigned n_element = Surface_contact_mesh_pt->nelement();
494 
495  // Loop over the surface elements
496  for(unsigned e=0;e<n_element;e++)
497  {
498  // Kill surface element
500  }
501 
502  // Wipe the mesh
504  }
void flush_element_and_node_storage()
Definition: mesh.h:407

References e().

◆ delete_flux_elements()

template<class ELEMENT >
void MeltContactProblem< ELEMENT >::delete_flux_elements ( )
inlineprivate

Delete flux elements.

548  {
549  // How many surface elements are in the surface mesh
550  unsigned n_element = Surface_flux_mesh_pt->nelement();
551 
552  // Loop over the surface elements
553  for(unsigned e=0;e<n_element;e++)
554  {
555  // Kill surface element
557  }
558 
559  // Wipe the mesh
561  }

References e().

◆ delete_melt_elements()

template<class ELEMENT >
void MeltContactProblem< ELEMENT >::delete_melt_elements ( )
inlineprivate

Delete melt elements.

428  {
429  // How many surface elements are in the surface mesh
430  unsigned n_element = Left_surface_melt_mesh_pt->nelement();
431 
432  // Loop over the surface elements
433  for(unsigned e=0;e<n_element;e++)
434  {
435  // Kill surface element
437  }
438 
439  // Wipe the mesh
441  }

References e().

◆ doc_solution()

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

Doc the solution.

1102 {
1103 
1104  oomph_info << "outputting for step: " << Doc_info.number() << std::endl;
1105 
1106 
1107  ofstream some_file;
1108  char filename[100];
1109 
1110  // Number of plot points
1111  unsigned npts;
1112  npts=5;
1113 
1114  // Output solution
1115  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
1116  Doc_info.number());
1117  some_file.open(filename);
1118  Bulk_mesh_pt->output(some_file,npts);
1119  some_file.close();
1120 
1121  // Output solution coarsely (only element vertices for easier
1122  // mesh visualisation)
1123  sprintf(filename,"%s/coarse_soln%i.dat",Doc_info.directory().c_str(),
1124  Doc_info.number());
1125  some_file.open(filename);
1126  Bulk_mesh_pt->output(some_file,2);
1127  some_file.close();
1128 
1129 
1130  // Output contact elements
1131  sprintf(filename,"%s/contact%i.dat",Doc_info.directory().c_str(),
1132  Doc_info.number());
1133  some_file.open(filename);
1134  unsigned nel=Surface_contact_mesh_pt->nelement();
1135  for (unsigned e=0;e<nel;e++)
1136  {
1138  Surface_contact_mesh_pt->element_pt(e))->output(some_file);
1139  }
1140  some_file.close();
1141 
1142  // Output flux with melt
1143  sprintf(filename,"%s/flux_with_melt%i.dat",Doc_info.directory().c_str(),
1144  Doc_info.number());
1145  some_file.open(filename);
1147  for (unsigned e=0;e<nel;e++)
1148  {
1149  dynamic_cast<SurfaceMeltElement<ELEMENT>* >(
1150  Left_surface_melt_mesh_pt->element_pt(e))->output_melt(some_file);
1151  }
1152  some_file.close();
1153 
1154 
1155  // Output exact solution
1156  sprintf(filename,"%s/exact_soln%i.dat",Doc_info.directory().c_str(),
1157  Doc_info.number());
1158  some_file.open(filename);
1159  Bulk_mesh_pt->output_fct(
1160  some_file,npts,time_pt()->time(),
1162  some_file.close();
1163 
1164  // Output exact position of melting line
1165  sprintf(filename,"%s/exact_height%i.dat",Doc_info.directory().c_str(),
1166  Doc_info.number());
1167  some_file.open(filename);
1168  unsigned nplot=100;
1169  for (unsigned j=0;j<nplot;j++)
1170  {
1171  double x=double(j)/double(nplot-1);
1172  some_file << x << " "
1174  << std::endl;
1175  }
1176  some_file.close();
1177 
1178  // Output exact solution
1179  sprintf(filename,"%s/exact_soln%i.dat",Doc_info.directory().c_str(),
1180  Doc_info.number());
1181  some_file.open(filename);
1182  Bulk_mesh_pt->output_fct(
1183  some_file,npts,time_pt()->time(),
1185  some_file.close();
1186 
1187  // Output exact position of melting line
1188  sprintf(filename,"%s/exact_height%i.dat",Doc_info.directory().c_str(),
1189  Doc_info.number());
1190  some_file.open(filename);
1191  nplot=100;
1192  for (unsigned j=0;j<nplot;j++)
1193  {
1194  double x=double(j)/double(nplot-1);
1195  some_file << x << " "
1197  << std::endl;
1198  }
1199  some_file.close();
1200 
1201 
1202  // Output penetrator
1203  sprintf(filename,"%s/penetrator%i.dat",Doc_info.directory().c_str(),
1204  Doc_info.number());
1205  some_file.open(filename);
1206  unsigned n=100;
1208  some_file.close();
1209 
1210  // Output Number of Newton iterations in form that can be visualised
1211  // as vector in paraview
1212  sprintf(filename,"%s/newton_iter%i.dat",Doc_info.directory().c_str(),
1213  Doc_info.number());
1214  some_file.open(filename);
1215  some_file << "0 0 0 " << Nnewton_iter_taken << std::endl;
1216  some_file.close();
1217 
1218  // Write norm of solution to trace file
1219  double norm=0.0;
1220  //Bulk_mesh_pt->compute_norm(norm);
1221  Trace_file << norm << std::endl;
1222 
1223  //Increment counter for solutions
1224  Doc_info.number()++;
1225 
1226 } // end of doc_solution
std::string directory() const
Output directory.
Definition: oomph_utilities.h:524
virtual void output(std::ostream &outfile, const unsigned &nplot) const =0
Output coordinates of penetrator at nplot plot points.
unsigned Nnewton_iter_taken
Definition: problem.h:603
double melting_surface_height(const double &t, const double &x)
Height of melting surface.
Definition: circular_boulder_melt.cc:156
void get_exact_u_for_unsteady_heat_validation(const double &t, const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
Definition: circular_boulder_melt.cc:112
list x
Definition: plotDoE.py:28
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.
Definition: overloaded_element_body.h:490

References oomph::DocInfo::directory(), GlobalParameters::Doc_info, e(), MergeRestartFiles::filename, ExactSolution::get_exact_u_for_unsteady_heat_validation(), j, ExactSolution::melting_surface_height(), n, oomph::DocInfo::number(), oomph::oomph_info, output(), oomph::Penetrator::output(), ProblemParameters::Penetrator_pt, oomph::Problem_Parameter::Trace_file, and plotDoE::x.

◆ global_temporal_error_norm()

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

Dummy global error norm for adaptive time-stepping.

Reimplemented from oomph::Problem.

380 {return 0.0;}

Member Data Documentation

◆ Backed_up_left_surface_melt_mesh_pt

template<class ELEMENT >
BackupMeshForProjection<TElement<1,3> >* MeltContactProblem< ELEMENT >::Backed_up_left_surface_melt_mesh_pt
private

Backup of Left_surface_melt_mesh_pt so the Lagrange multipliers and melt rate can be projected across

◆ Bottom_boundary_id

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

ID of bottom boundary.

◆ Bulk_mesh_pt

template<class ELEMENT >
RefineableSolidTriangleMesh<ELEMENT>* MeltContactProblem< ELEMENT >::Bulk_mesh_pt
private

Pointer to bulk mesh.

◆ Doc_info

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

◆ Left_boundary_id

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

ID of left boundary.

◆ Left_melt_boundary_id

template<class ELEMENT >
unsigned MeltContactProblem< ELEMENT >::Left_melt_boundary_id
private

ID of left melt boundary.

◆ Left_surface_melt_mesh_pt

template<class ELEMENT >
Mesh* MeltContactProblem< ELEMENT >::Left_surface_melt_mesh_pt
private

Pointer to the "surface" mesh.

◆ Lower_contact_boundary_id

template<class ELEMENT >
unsigned MeltContactProblem< ELEMENT >::Lower_contact_boundary_id
private

ID of lower contact boundary.

◆ Right_boundary_id

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

ID of right boundary.

◆ Right_melt_boundary_id

template<class ELEMENT >
unsigned MeltContactProblem< ELEMENT >::Right_melt_boundary_id
private

ID of right melt boundary.

◆ Surface_contact_mesh_pt

template<class ELEMENT >
Mesh* MeltContactProblem< ELEMENT >::Surface_contact_mesh_pt
private

Pointer to the "surface" mesh.

◆ Surface_flux_mesh_pt

template<class ELEMENT >
Mesh* MeltContactProblem< ELEMENT >::Surface_flux_mesh_pt
private

Pointer to the "surface" mesh.

◆ Trace_file

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

Trace file.

◆ Upper_contact_boundary_id

template<class ELEMENT >
unsigned MeltContactProblem< ELEMENT >::Upper_contact_boundary_id
private

ID of upper contact boundary.


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