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
 
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 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...
 
MeshSurface_contact_mesh_pt
 Pointer to the "surface" mesh. More...
 
MeshLeft_surface_melt_mesh_pt
 Pointer to the "surface" mesh. More...
 
MeshSurface_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)
 
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 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: