SolarRadiationProblem< ELEMENT > Class Template Reference

Problem class. More...

+ Inheritance diagram for SolarRadiationProblem< ELEMENT >:

Public Member Functions

 SolarRadiationProblem ()
 Constructor. More...
 
 ~SolarRadiationProblem ()
 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. More...
 
void actions_before_adapt ()
 Actions before adapt: wipe contact elements. More...
 
void actions_after_adapt ()
 
void update_limiting_angles ()
 Update limiting angles for diffuse radiation. More...
 
void setup_shielding_nodes ()
 
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_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 * 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
 
Vector< Node * > Ordered_shielding_node_pt
 Storage for ordered shielding nodes. More...
 
Vector< FiniteElement * > Ordered_shielding_face_element_pt
 Face elements on potentially sun-exposed boundaries. More...
 

Additional Inherited Members

- Public Types inherited from oomph::Problem
typedef void(* SpatialErrorEstimatorFctPt) (Mesh *&mesh_pt, Vector< double > &elemental_error)
 Function pointer for spatial error estimator. More...
 
typedef void(* SpatialErrorEstimatorWithDocFctPt) (Mesh *&mesh_pt, Vector< double > &elemental_error, DocInfo &doc_info)
 Function pointer for spatial error estimator with doc. More...
 
- Public Attributes inherited from oomph::Problem
bool Shut_up_in_newton_solve
 
- Static Public Attributes inherited from oomph::Problem
static bool Suppress_warning_about_actions_before_read_unstructured_meshes
 
- Protected Types inherited from oomph::Problem
enum  Assembly_method {
  Perform_assembly_using_vectors_of_pairs , Perform_assembly_using_two_vectors , Perform_assembly_using_maps , Perform_assembly_using_lists ,
  Perform_assembly_using_two_arrays
}
 Enumerated flags to determine which sparse assembly method is used. More...
 
- Protected Member Functions inherited from oomph::Problem
unsigned setup_element_count_per_dof ()
 
virtual void sparse_assemble_row_or_column_compressed (Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
 
virtual void actions_before_newton_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 SolarRadiationProblem< ELEMENT >

Problem class.

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

Constructor & Destructor Documentation

◆ SolarRadiationProblem()

template<class ELEMENT >
SolarRadiationProblem< ELEMENT >::SolarRadiationProblem

Constructor.

Constructor for contact problem in square domain.

Identify penetrator as region 1

779 {
780 
781  // Output directory
782  Doc_info.set_directory("RESLT");
783 
784  // Output number
785  Doc_info.number()=0;
786 
787  // Open trace file
788  Trace_file.open("RESLT/trace.dat");
789 
790  // Allow for crap initial guess
791  Problem::Max_residuals=10000.0;
792 
793  // Allocate the timestepper -- this constructs the Problem's
794  // time object with a sufficient amount of storage to store the
795  // previous timsteps.
797 
798  // Pointer to the closed curve that defines the outer boundary
799  TriangleMeshClosedCurve* closed_curve_pt=0;
800 
801  // Build outer boundary as Polygon
802 
803  // The boundary, represented by polylines
804  Vector<TriangleMeshCurveSection*> boundary_polyline_pt;
805 
806  // Store coordinates of left and right point at which penetrator
807  // penetrates flat surface
808  Vector<double> right_penetrator(2);
809  Vector<double> left_penetrator(2);
810 
811  // Vertex coordinates on boundary
812  Vector<Vector<double> > bound_coords(2);
813 
814  // Left boundary
815  bound_coords[0].resize(2);
816  bound_coords[0][0]=0.0;
817  bound_coords[0][1]=1.0;
818 
819  bound_coords[1].resize(2);
820  bound_coords[1][0]=0.0;
821  bound_coords[1][1]=0.0;
822 
823  // Build the boundary polyline
825  boundary_polyline_pt.push_back(new TriangleMeshPolyLine(bound_coords,
827 
828  // Bottom boundary
829  bound_coords[0].resize(2);
830  bound_coords[0][0]=0.0;
831  bound_coords[0][1]=0.0;
832 
833  bound_coords[1].resize(2);
834  bound_coords[1][0]=1.0;
835  bound_coords[1][1]=0.0;
836 
837  // Build the boundary polyline
839  boundary_polyline_pt.push_back(new TriangleMeshPolyLine(bound_coords,
841 
842 
843  // Right boundary
844  bound_coords[0].resize(2);
845  bound_coords[0][0]=1.0;
846  bound_coords[0][1]=0.0;
847 
848  bound_coords[1].resize(2);
849  bound_coords[1][0]=1.0;
850  bound_coords[1][1]=1.0;
851 
852  // Build the boundary polyline
854  boundary_polyline_pt.push_back(new TriangleMeshPolyLine(bound_coords,
856 
857 
858  // Half opening angle of "boulder"
859  double half_opening_angle=acos((ProblemParameters::Centre[1]-1.0)/
861 
862 
863  // Points at which the penetrator meets the flat surface
864  right_penetrator[0]=ProblemParameters::Centre[0]+
865  ProblemParameters::Radius*sin(half_opening_angle);
866  right_penetrator[1]=1.0;
867 
868  left_penetrator[0]=ProblemParameters::Centre[0]-
869  ProblemParameters::Radius*sin(half_opening_angle);
870  left_penetrator[1]=1.0;
871 
872  // Right melt boundary
873  bound_coords[0].resize(2);
874  bound_coords[0][0]=1.0;
875  bound_coords[0][1]=1.0;
876 
877  bound_coords[1].resize(2);
878  bound_coords[1][0]=right_penetrator[0];
879  bound_coords[1][1]=right_penetrator[1];
880 
881 
882  // Build boundary poly line
884  boundary_polyline_pt.push_back
885  (new TriangleMeshPolyLine(bound_coords,Right_melt_boundary_id));
886 
887 
888 
889  // Vertex coordinates on lower contact boundary
890  unsigned nvertex_contact=20;
891  Vector<Vector<double> > circle_bound_coords(nvertex_contact);
892 
893  // First point
894  circle_bound_coords[0].resize(2);
895  circle_bound_coords[0][0]=right_penetrator[0];
896  circle_bound_coords[0][1]=right_penetrator[1];
897 
898  // Internal points
899  for (unsigned j=1;j<nvertex_contact-1;j++)
900  {
901  double phi=1.5*MathematicalConstants::Pi+half_opening_angle
902  -2.0*half_opening_angle*double(j)/double(nvertex_contact-1);
903 
904  circle_bound_coords[j].resize(2);
905  circle_bound_coords[j][0]=ProblemParameters::Centre[0]+
907 
908  circle_bound_coords[j][1]=ProblemParameters::Centre[1]+
910  }
911 
912  // Last point
913  circle_bound_coords[nvertex_contact-1].resize(2);
914  circle_bound_coords[nvertex_contact-1][0]=left_penetrator[0];
915  circle_bound_coords[nvertex_contact-1][1]=left_penetrator[1];
916 
917 
918  // Build boundary poly line
920  TriangleMeshPolyLine* lower_contact_boundary_pt=
921  new TriangleMeshPolyLine(circle_bound_coords,Lower_contact_boundary_id);
922 
923  // Store as internal poly line
924  Vector<TriangleMeshCurveSection*> internal_polyline_pt(1);
925  internal_polyline_pt[0]=lower_contact_boundary_pt;
926 
927  // Vertex coordinates on upper contact boundary
928  unsigned nvertex_contact_upper=unsigned(double(nvertex_contact)*
930  half_opening_angle);
931  circle_bound_coords.clear();
932  circle_bound_coords.resize(nvertex_contact_upper);
933 
934  // First point: recycle
935  circle_bound_coords[0].resize(2);
936  circle_bound_coords[0][0]=right_penetrator[0];
937  circle_bound_coords[0][1]=right_penetrator[1];
938  for (unsigned j=1;j<nvertex_contact_upper-1;j++)
939  {
940  double phi=-(0.5*MathematicalConstants::Pi-half_opening_angle)
941  +(2.0*MathematicalConstants::Pi-2.0*half_opening_angle)*
942  double(j)/double(nvertex_contact_upper-1);
943 
944  circle_bound_coords[j].resize(2);
945  circle_bound_coords[j][0]=ProblemParameters::Centre[0]+
947 
948  circle_bound_coords[j][1]=ProblemParameters::Centre[1]+
950  }
951 
952  // Last point
953  circle_bound_coords[nvertex_contact_upper-1].resize(2);
954  circle_bound_coords[nvertex_contact_upper-1][0]=left_penetrator[0];
955  circle_bound_coords[nvertex_contact_upper-1][1]=left_penetrator[1];
956 
957  // Build boundary poly line
959  TriangleMeshPolyLine* upper_contact_boundary_pt=
960  new TriangleMeshPolyLine(circle_bound_coords,Upper_contact_boundary_id);
961  boundary_polyline_pt.push_back(upper_contact_boundary_pt);
962 
963  // Left melt boundary
964  Vector<Vector<double> > left_melt_bound_coords;
965  Vector<double> tmp_vector(2);
966  tmp_vector[0]=left_penetrator[0];
967  tmp_vector[1]=left_penetrator[1];
968  left_melt_bound_coords.push_back(tmp_vector);
969 
970 
971  // Funky bit
972  {
973  tmp_vector[0]=0.2;
974  tmp_vector[1]=1.0;
975  left_melt_bound_coords.push_back(tmp_vector);
976 
977  tmp_vector[0]=0.2;
978  tmp_vector[1]=2.2;
979  left_melt_bound_coords.push_back(tmp_vector);
980 
981  tmp_vector[0]=0.5;
982  tmp_vector[1]=2.2;
983  left_melt_bound_coords.push_back(tmp_vector);
984 
985  tmp_vector[0]=0.5;
986  tmp_vector[1]=1.7;
987  left_melt_bound_coords.push_back(tmp_vector);
988 
989  tmp_vector[0]=0.4;
990  tmp_vector[1]=1.7;
991  left_melt_bound_coords.push_back(tmp_vector);
992 
993  tmp_vector[0]=0.4;
994  tmp_vector[1]=2.0;
995  left_melt_bound_coords.push_back(tmp_vector);
996 
997  tmp_vector[0]=0.35;
998  tmp_vector[1]=2.0;
999  left_melt_bound_coords.push_back(tmp_vector);
1000 
1001  tmp_vector[0]=0.35;
1002  tmp_vector[1]=1.65;
1003  left_melt_bound_coords.push_back(tmp_vector);
1004 
1005  tmp_vector[0]=0.55;
1006  tmp_vector[1]=1.65;
1007  left_melt_bound_coords.push_back(tmp_vector);
1008 
1009 
1010  tmp_vector[0]=0.55;
1011  tmp_vector[1]=2.25;
1012  left_melt_bound_coords.push_back(tmp_vector);
1013 
1014  tmp_vector[0]=0.15;
1015  tmp_vector[1]=2.25;
1016  left_melt_bound_coords.push_back(tmp_vector);
1017 
1018  tmp_vector[0]=0.15;
1019  tmp_vector[1]=1.0;
1020  left_melt_bound_coords.push_back(tmp_vector);
1021  }
1022 
1023 
1024  tmp_vector[0]=0.0;
1025  tmp_vector[1]=1.0;
1026  left_melt_bound_coords.push_back(tmp_vector);
1027 
1028  // Build boundary poly line
1030  TriangleMeshPolyLine* left_melt_boundary_pt=
1031  new TriangleMeshPolyLine(left_melt_bound_coords,Left_melt_boundary_id);
1032  boundary_polyline_pt.push_back(left_melt_boundary_pt);
1033 
1034  // Connect first vertex (vertex 0) of upper contact boundary
1035  // with initial vertex of lower contact boundary
1036  unsigned vertex_id_of_connection=0;
1037  lower_contact_boundary_pt->connect_initial_vertex_to_polyline(
1038  dynamic_cast<TriangleMeshPolyLine*>(upper_contact_boundary_pt),
1039  vertex_id_of_connection);
1040 
1041 
1042 //#define USE_BROKEN
1043 #ifdef USE_BROKEN
1044 
1045  // hierher this misses out final node on lower contact boundary
1046  // in boundary lookup scheme!
1047 
1048  // Connect last vertex of upper contact boundary
1049  // with final vertex of lower contact boundary
1050  vertex_id_of_connection=nvertex_contact_upper-1;
1051  lower_contact_boundary_pt->connect_final_vertex_to_polyline(
1052  dynamic_cast<TriangleMeshPolyLine*>(upper_contact_boundary_pt),
1053  vertex_id_of_connection);
1054 
1055 #else
1056 
1057  // Connect first vertex left melt boundary
1058  // with final vertex of lower contact boundary
1059  vertex_id_of_connection=0;
1060  lower_contact_boundary_pt->connect_final_vertex_to_polyline(
1061  dynamic_cast<TriangleMeshPolyLine*>(left_melt_boundary_pt),
1062  vertex_id_of_connection);
1063 
1064 #endif
1065 
1066  // Create open curve that defines boulder/ice interface
1067  Vector<TriangleMeshOpenCurve*> inner_boundary_pt;
1068  inner_boundary_pt.push_back(new TriangleMeshOpenCurve(internal_polyline_pt));
1069 
1070 
1071  // Create the triangle mesh polygon for outer boundary
1072  //----------------------------------------------------
1073  TriangleMeshPolygon *outer_polygon =
1074  new TriangleMeshPolygon(boundary_polyline_pt);
1075 
1076  // Set the pointer
1077  closed_curve_pt = outer_polygon;
1078 
1079  // Now build the mesh
1080  //===================
1081 
1082  // Use the TriangleMeshParameters object for helping on the manage of the
1083  // TriangleMesh parameters
1084  TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
1085 
1086  // Specify the maximum area element
1087  double uniform_element_area=0.002;
1088  triangle_mesh_parameters.element_area() = uniform_element_area;
1089 
1090  // Specify the internal open boundary
1091  triangle_mesh_parameters.internal_open_curves_pt()=inner_boundary_pt;
1092 
1094  triangle_mesh_parameters.add_region_coordinates(1,ProblemParameters::Centre);
1095 
1096 
1097 #ifdef ADAPTIVE
1098 
1099  // Create the mesh
1100  Bulk_mesh_pt=new RefineableSolidTriangleMesh<ELEMENT>(triangle_mesh_parameters,
1101  time_stepper_pt());
1102 
1103  // Set error estimator for bulk mesh
1105  Bulk_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
1106 
1107  // // Set element size limits
1108  // Bulk_mesh_pt->max_element_size()=0.2;
1109  // Bulk_mesh_pt->min_element_size()=0.002;
1110 
1111 #else
1112 
1113  // Build mesh
1114  Bulk_mesh_pt=new SolidTriangleMesh<ELEMENT>(triangle_mesh_parameters,
1115  time_stepper_pt());
1116 
1117 #endif
1118 
1119  Bulk_mesh_pt->output("mesh.dat");
1120  Bulk_mesh_pt->output_boundaries("boundaries.dat");
1121 
1122  // Create the surface mesh as an empty mesh
1124 
1125  // Build 'em
1127 
1128  // Create the surface mesh as an empty mesh
1130 
1131  // Build 'em
1133 
1134  // Set boundary condition and complete the build of all elements
1136 
1137  // Add the sub meshes to the problem
1141 
1142  // Combine all submeshes into a single global Mesh
1144 
1145  // Set the initial conditions
1146  unsigned nnod = Bulk_mesh_pt->nnode();
1147  for(unsigned j=0;j<nnod;j++)
1148  {
1149  Bulk_mesh_pt->node_pt(j)->set_value(0,ExactSolution::U0);
1150  }
1151 
1152  // Do equation numbering
1153  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
1154 
1155 } // 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
void create_contact_elements()
Create contact elements.
Definition: circular_boulder_solar_radiation.cc:483
Mesh * Surface_flux_mesh_pt
Pointer to the "surface" mesh.
Definition: circular_boulder_solar_radiation.cc:736
unsigned Right_melt_boundary_id
ID of right melt boundary.
Definition: circular_boulder_solar_radiation.cc:742
unsigned Lower_contact_boundary_id
ID of lower contact boundary.
Definition: circular_boulder_solar_radiation.cc:745
unsigned Right_boundary_id
ID of right boundary.
Definition: circular_boulder_solar_radiation.cc:757
unsigned Left_boundary_id
ID of left boundary.
Definition: circular_boulder_solar_radiation.cc:754
DocInfo Doc_info
Definition: circular_boulder_solar_radiation.cc:763
void complete_problem_setup()
Definition: circular_boulder_solar_radiation.cc:598
Mesh * Surface_contact_mesh_pt
Pointer to the "surface" mesh.
Definition: circular_boulder_solar_radiation.cc:733
ofstream Trace_file
Trace file.
Definition: circular_boulder_solar_radiation.cc:760
void create_flux_elements()
Create flux elements.
Definition: circular_boulder_solar_radiation.cc:540
unsigned Bottom_boundary_id
ID of bottom boundary.
Definition: circular_boulder_solar_radiation.cc:751
RefineableSolidTriangleMesh< ELEMENT > * Bulk_mesh_pt
Pointer to bulk mesh.
Definition: circular_boulder_solar_radiation.cc:723
unsigned Upper_contact_boundary_id
ID of upper contact boundary.
Definition: circular_boulder_solar_radiation.cc:748
unsigned Left_melt_boundary_id
ID of left melt boundary.
Definition: circular_boulder_solar_radiation.cc:739
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
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(), sin(), oomph::Problem_Parameter::Trace_file, and ExactSolution::U0.

◆ ~SolarRadiationProblem()

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

Destructor (empty)

263 {}

Member Function Documentation

◆ actions_after_adapt()

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

328  {
329  // Create elements
332 
333  // Rebuild the Problem's global mesh from its various sub-meshes
335 
336  // Rebuild elements
338  }
void rebuild_global_mesh()
Definition: problem.cc:1533

◆ actions_after_newton_solve()

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

Update the problem specs after solve (empty)

Reimplemented from oomph::Problem.

266 {}

◆ actions_before_adapt()

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

Actions before adapt: wipe contact elements.

Reimplemented from oomph::Problem.

314  {
315  // Kill the elements and wipe surface mesh
318 
319  // Rebuild the Problem's global mesh from its various sub-meshes
321  }
void delete_contact_elements()
Delete contact elements.
Definition: circular_boulder_solar_radiation.cc:523
void delete_flux_elements()
Delete flux elements.
Definition: circular_boulder_solar_radiation.cc:579

◆ actions_before_implicit_timestep()

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

Actions before next timestep.

Reimplemented from oomph::Problem.

273  {
274  // Amplitude of oscillation
275  double amplitude=0.2; //-0.2;
276 
278  (1.0+0.05*time_pt()->time())*
279  0.5*(1.0-cos(2.0*MathematicalConstants::Pi*time_pt()->time()));
280 
281  oomph_info << "Solving for y_c = "
282  << ProblemParameters::Centre[1] << std::endl;
283 
284  // If update of limiting angles is skipped during Newton iteration
285  // update it at least before every timestep
287  "--skip_update_limiting_angles_during_newton_iteration"))
288  {
290  }
291 
292  }
void update_limiting_angles()
Update limiting angles for diffuse radiation.
Definition: circular_boulder_solar_radiation.cc:342
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
bool command_line_flag_has_been_set(const std::string &flag)
Definition: oomph_utilities.cc:501
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

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

◆ actions_before_newton_convergence_check()

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

Newton convergence check.

Reimplemented from oomph::Problem.

297  {
298  // Sometimes the Newton iteration
299  // stagnates, probably because of the shadow jumping
300  // over integration points. Not updating limiting angles
301  // is a bit naughty but the error incurred should be
302  // small-ish if timesteps (i.e. change in surface shape remain
303  // small between solves).
305  "--skip_update_limiting_angles_during_newton_iteration"))
306  {
308  }
309  }

References oomph::CommandLineArgs::command_line_flag_has_been_set().

◆ actions_before_newton_solve()

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

Update the problem specs before solve (empty)

Reimplemented from oomph::Problem.

269 {}

◆ complete_problem_setup()

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

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

599  {
600 
601 
602  // Set material properties for elements in rock region
603  //---------------------------------------------------
604  unsigned r=1; // hierher enumerate globally
605  unsigned nel=Bulk_mesh_pt->nregion_element(r);
606  for (unsigned e=0;e<nel;e++)
607  {
608  ELEMENT *el_pt=dynamic_cast<ELEMENT*>(
609  Bulk_mesh_pt->region_element_pt(r,e));
610 
611  // Non-dim thermal inertia
612  el_pt->alpha_pt()=&ProblemParameters::Alpha_rock;
613 
614  // Non-dim thermal inertia
615  el_pt->beta_pt()=&ProblemParameters::Beta_rock;
616  }
617 
618 
619  // Set (pseudo-)solid mechanics properties for all elements
620  //---------------------------------------------------------
621  unsigned n_element = Bulk_mesh_pt->nelement();
622  for(unsigned e=0;e<n_element;e++)
623  {
624  //Cast to a solid element
625  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
626 
627  // Set the constitutive law
628  el_pt->constitutive_law_pt() =
630 
631  // Set density to zero
632  el_pt->lambda_sq_pt()=&ProblemParameters::Lambda_sq;
633 
634  // Disable inertia
635  el_pt->disable_inertia();
636 
637  // Set source function
638  //el_pt->source_fct_pt() =
639  // &ExactSolution::get_source_for_unsteady_heat_validation;
640  }
641 
642  // Apply boundary conditions for solid
643  //------------------------------------
644 
645  // Bottom: completely pinned
646  unsigned b=Bottom_boundary_id;
647  unsigned nnod=Bulk_mesh_pt->nboundary_node(b);
648  for (unsigned j=0;j<nnod;j++)
649  {
650  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
651  nod_pt->pin_position(0);
652  nod_pt->pin_position(1);
653  }
654 
655  // Sides: Symmetry bcs
657  nnod=Bulk_mesh_pt->nboundary_node(b);
658  for (unsigned j=0;j<nnod;j++)
659  {
660  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
661  nod_pt->pin_position(0);
662  }
664  nnod=Bulk_mesh_pt->nboundary_node(b);
665  for (unsigned j=0;j<nnod;j++)
666  {
667  SolidNode* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
668  nod_pt->pin_position(0);
669  }
670 
671  // Assign the Lagrangian coordinates -- sensible
672  // because we've completely rebuilt the mesh
673  // and haven't copied across any Lagrange multipliers
674  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
675 
676  // Loop over the contact elements, pass pointer to penetrator and make sticky
677  //---------------------------------------------------------------------------
678  n_element=Surface_contact_mesh_pt->nelement();
679  for(unsigned e=0;e<n_element;e++)
680  {
681  // Upcast from GeneralisedElement
685 
686  // Set pointer to penetrator
688 
689  // Make it sticky to enforce permanent contact
690  el_pt->enable_stick();
691  }
692 
693 
694  // Loop over the flux elements to pass pointer to prescribed flux function
695  //------------------------------------------------------------------------
696  n_element=Surface_flux_mesh_pt->nelement();
697  for(unsigned e=0;e<n_element;e++)
698  {
699  // Upcast from GeneralisedElement
703 
704  // Set atmospheric radiation fct
707  }
708 
709 
710  // Setup shielding nodes
712 
713  // hierher
715 
716  }
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Scalar * b
Definition: benchVecAdd.cpp:17
void setup_shielding_nodes()
Definition: circular_boulder_solar_radiation.cc:360
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
Definition: contact_elements.h:1086
void(*&)(const double &time, double &solar_flux_magnitude, Vector< double > &solar_flux_unit_vector, double &total_diffuse_radiation) atmospheric_radiation_fct_pt()
Reference to the atmospheric radiation function pointer.
Definition: temporary_stefan_boltzmann_elements.h:209
Definition: nodes.h:1686
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1816
Definition: temporary_stefan_boltzmann_elements.h:1572
void set_penetrator_pt(Penetrator *penetrator_pt)
Set pointer to penetrator.
Definition: contact_elements.h:844
void enable_stick()
Definition: contact_elements.h:776
Penetrator * Penetrator_pt
Penetrator.
Definition: circular_boulder_melt.cc:88
void atmospheric_radiation(const double &time, double &solar_flux_magnitude, Vector< double > &solar_flux_unit_vector, double &total_diffuse_radiation)
Definition: circular_boulder_solar_radiation.cc:101
double Beta_rock
Nondim thermal conductivity for rock.
Definition: circular_boulder_solar_radiation.cc:123
double Alpha_rock
Nondim thermal inertia for rock.
Definition: circular_boulder_solar_radiation.cc:120
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
r
Definition: UniformPSDSelfTest.py:20

References ProblemParameters::Alpha_rock, ProblemParameters::atmospheric_radiation(), oomph::SolarRadiationBase::atmospheric_radiation_fct_pt, b, ProblemParameters::Beta_rock, ProblemParameters::Constitutive_law_pt, e(), oomph::SurfaceContactElementBase< ELEMENT >::enable_stick(), j, ProblemParameters::Lambda_sq, ProblemParameters::Penetrator_pt, oomph::SolidNode::pin_position(), UniformPSDSelfTest::r, and oomph::SurfaceContactElementBase< ELEMENT >::set_penetrator_pt().

◆ create_contact_elements()

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

Create contact elements.

484  {
485  // Which boundaries are in contact with penetrator?
486  Vector<unsigned> contact_boundaries;
487  contact_boundaries.push_back(Upper_contact_boundary_id);
488  contact_boundaries.push_back(Lower_contact_boundary_id);
489  unsigned n=contact_boundaries.size();
490  for (unsigned bb=0;bb<n;bb++)
491  {
492  // How many bulk elements are adjacent to boundary b?
493  unsigned b=contact_boundaries[bb];
494  unsigned n_element = Bulk_mesh_pt->nboundary_element_in_region(b,1);
495 
496  // Loop over the bulk elements adjacent to boundary b?
497  for(unsigned e=0;e<n_element;e++)
498  {
499  // Get pointer to the bulk element that is adjacent to boundary b
500  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
501  Bulk_mesh_pt->boundary_element_in_region_pt(b,1,e));
502 
503  //What is the face index of element e along boundary b
504  int face_index = Bulk_mesh_pt->face_index_at_boundary_in_region(b,1,e);
505 
506  // Id of additional dofs created by this element
507  unsigned id=0; // hierher make member
508 
509  // Build the corresponding contact element
510  NonlinearSurfaceContactElement<ELEMENT>* contact_element_pt = new
511  NonlinearSurfaceContactElement<ELEMENT>(bulk_elem_pt,face_index,id);
512 
513  //Add the contact element to the surface mesh
514  Surface_contact_mesh_pt->add_element_pt(contact_element_pt);
515 
516  } //end of loop over bulk elements adjacent to boundary b
517  }
518  }
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(), and n.

◆ create_flux_elements()

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

Create flux elements.

541  {
542  // Which boundaries are subject to incoming radiation
543  Vector<unsigned> flux_boundaries;
544  flux_boundaries.push_back(Upper_contact_boundary_id);
545  flux_boundaries.push_back(Left_melt_boundary_id);
546  flux_boundaries.push_back(Right_melt_boundary_id);
547  unsigned n=flux_boundaries.size();
548  for (unsigned bb=0;bb<n;bb++)
549  {
550  // How many bulk elements are adjacent to boundary b?
551  unsigned b=flux_boundaries[bb];
552  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
553 
554  // Loop over the bulk elements adjacent to boundary b?
555  for(unsigned e=0;e<n_element;e++)
556  {
557  // Get pointer to the bulk element that is adjacent to boundary b
558  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
559  Bulk_mesh_pt->boundary_element_pt(b,e));
560 
561  //What is the face index of element e along boundary b
562  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
563 
564  // Create flux element
567  face_index);
568 
569  // Add to mesh
571 
572  } //end of loop over bulk elements adjacent to boundary b
573  }
574  }

References b, e(), and n.

◆ delete_contact_elements()

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

Delete contact elements.

524  {
525  // How many surface elements are in the surface mesh
526  unsigned n_element = Surface_contact_mesh_pt->nelement();
527 
528  // Loop over the surface elements
529  for(unsigned e=0;e<n_element;e++)
530  {
531  // Kill surface element
533  }
534 
535  // Wipe the mesh
537  }
void flush_element_and_node_storage()
Definition: mesh.h:407

References e().

◆ delete_flux_elements()

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

Delete flux elements.

580  {
581  // How many surface elements are in the surface mesh
582  unsigned n_element = Surface_flux_mesh_pt->nelement();
583 
584  // Loop over the surface elements
585  for(unsigned e=0;e<n_element;e++)
586  {
587  // Kill surface element
589  }
590 
591  // Wipe the mesh
593  }

References e().

◆ doc_solution()

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

Doc the solution.

1165 {
1166 
1167  oomph_info << "outputting step: " << Doc_info.number() << std::endl;
1168 
1169  ofstream some_file;
1170  char filename[100];
1171 
1172  // Number of plot points
1173  unsigned npts;
1174  npts=5;
1175 
1176  // Output solution
1177  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
1178  Doc_info.number());
1179  some_file.open(filename);
1180  Bulk_mesh_pt->output(some_file,npts);
1181  some_file.close();
1182 
1183  // Output solution coarsely (only element vertices for easier
1184  // mesh visualisation)
1185  sprintf(filename,"%s/coarse_soln%i.dat",Doc_info.directory().c_str(),
1186  Doc_info.number());
1187  some_file.open(filename);
1188  Bulk_mesh_pt->output(some_file,2);
1189  some_file.close();
1190 
1191  // Output contact elements
1192  sprintf(filename,"%s/contact%i.dat",Doc_info.directory().c_str(),
1193  Doc_info.number());
1194  some_file.open(filename);
1195  unsigned nel=Surface_contact_mesh_pt->nelement();
1196  for (unsigned e=0;e<nel;e++)
1197  {
1199  Surface_contact_mesh_pt->element_pt(e))->output(some_file);
1200  }
1201  some_file.close();
1202 
1203  // Output exact solution
1204  sprintf(filename,"%s/exact_soln%i.dat",Doc_info.directory().c_str(),
1205  Doc_info.number());
1206  some_file.open(filename);
1207  Bulk_mesh_pt->output_fct(
1208  some_file,npts,time_pt()->time(),
1210  some_file.close();
1211 
1212  // Output Solar radiation
1213  double solar_flux_magnitude=0.0;
1214  Vector<double> solar_flux_unit_vector(2);
1215  double total_diffuse_radiation=0.0;
1217  solar_flux_magnitude,
1218  solar_flux_unit_vector,
1219  total_diffuse_radiation);
1220  sprintf(filename,"%s/solar_radiation%i.dat",Doc_info.directory().c_str(),
1221  Doc_info.number());
1222  some_file.open(filename);
1223  some_file << "0.0 3.5 "
1224  << solar_flux_magnitude << " "
1225  << solar_flux_unit_vector[0] << " "
1226  << solar_flux_unit_vector[1] << " "
1227  << "\n";
1228  some_file.close();
1229 
1230 
1231  // Output atmospheric radiation along all exposed surfaces
1232  sprintf(filename,"%s/atmospheric_radiation%i.dat",
1233  Doc_info.directory().c_str(),
1234  Doc_info.number());
1235  some_file.open(filename);
1237  for (unsigned e=0;e<nel;e++)
1238  {
1241  element_pt(e))->output_atmospheric_radiation(some_file);
1242  }
1243  some_file.close();
1244 
1245 
1246  // Output illumination angles for all integration points
1247  sprintf(filename,"%s/illumination_angles%i.dat",
1248  Doc_info.directory().c_str(),
1249  Doc_info.number());
1250  some_file.open(filename);
1252  for (unsigned e=0;e<nel;e++)
1253  {
1256  element_pt(e))->output_limiting_angles(some_file);
1257  }
1258  some_file.close();
1259 
1260  bool plot_individual=true;
1261  if (plot_individual)
1262  {
1263 
1264  // Output cone of diffuse radiation for all integration points
1265  // (suitable for paraview animation -- each element is in a separate
1266  // file)
1267  double radius=2.0;
1269  unsigned count=0;
1270  for (unsigned e=0;e<nel;e++)
1271  {
1272  sprintf(filename,"%s/diffuse_radiation_cone%i_%i.dat",
1273  Doc_info.directory().c_str(),
1274  Doc_info.number(),count++);
1275  some_file.open(filename);
1278  ->output_diffuse_radiation_cone(some_file,radius);
1279  some_file.close();
1280  }
1281  count=0;
1282  for (unsigned e=0;e<nel;e++)
1283  {
1284  sprintf(filename,"%s/diffuse_radiation_cone_max_angle%i_%i.dat",
1285  Doc_info.directory().c_str(),
1286  Doc_info.number(),count++);
1287  some_file.open(filename);
1290  ->output_diffuse_radiation_cone_max_angle(some_file,radius);
1291  some_file.close();
1292  }
1293  count=0;
1294  for (unsigned e=0;e<nel;e++)
1295  {
1296  sprintf(filename,"%s/diffuse_radiation_cone_min_angle%i_%i.dat",
1297  Doc_info.directory().c_str(),
1298  Doc_info.number(),count++);
1299  some_file.open(filename);
1302  ->output_diffuse_radiation_cone_min_angle(some_file,radius);
1303  some_file.close();
1304  }
1305  }
1306 
1307 
1308  // Output penetrator
1309  sprintf(filename,"%s/penetrator%i.dat",Doc_info.directory().c_str(),
1310  Doc_info.number());
1311  some_file.open(filename);
1312  unsigned n=100;
1314  some_file.close();
1315 
1316  // Output Number of Newton iterations in form that can be visualised
1317  // as vector in paraview
1318  sprintf(filename,"%s/newton_iter%i.dat",Doc_info.directory().c_str(),
1319  Doc_info.number());
1320  some_file.open(filename);
1321  some_file << "0 0 0 " << Nnewton_iter_taken << std::endl;
1322  some_file.close();
1323 
1324  // Write norm of solution to trace file
1325  double norm=0.0;
1326  //Bulk_mesh_pt->compute_norm(norm);
1327  Trace_file << norm << std::endl;
1328 
1329  //Increment counter for solutions
1330  Doc_info.number()++;
1331 
1332 } // end of doc_solution
Vector< FiniteElement * > Ordered_shielding_face_element_pt
Face elements on potentially sun-exposed boundaries.
Definition: circular_boulder_solar_radiation.cc:769
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
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
string filename
Definition: MergeRestartFiles.py:39
radius
Definition: UniformPSDSelfTest.py:15
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.
Definition: overloaded_element_body.h:490

References ProblemParameters::atmospheric_radiation(), oomph::DocInfo::directory(), GlobalParameters::Doc_info, e(), MergeRestartFiles::filename, ExactSolution::get_exact_u_for_unsteady_heat_validation(), n, oomph::DocInfo::number(), oomph::oomph_info, output(), oomph::Penetrator::output(), ProblemParameters::Penetrator_pt, UniformPSDSelfTest::radius, and oomph::Problem_Parameter::Trace_file.

◆ global_temporal_error_norm()

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

Dummy global error norm for adaptive time-stepping.

Reimplemented from oomph::Problem.

477 {return 0.0;}

◆ setup_shielding_nodes()

template<class ELEMENT >
void SolarRadiationProblem< ELEMENT >::setup_shielding_nodes ( )
inline
361  {
362  // Wipe
365 
366  // Left horizon
367  // hierher Ordered_shielding_node_pt.push_back(Left_horizon_node_pt);
368 
369  // Associate elements with vertex nodes
370  std::map<Node*,std::set<FiniteElement*> > adj_el_pt;
371  unsigned nel=Surface_flux_mesh_pt->nelement();
372 
373  oomph_info << "Number of elements: " << nel << std::endl;
374 
375  for (unsigned e=0;e<nel;e++)
376  {
378  unsigned nnod=el_pt->nnode();
379 
380  // Associate element with first vertex node
381  Node* nod_pt=el_pt->node_pt(0);
382  adj_el_pt[nod_pt].insert(el_pt);
383 
384  // Associate element with last vertex node
385  nod_pt=el_pt->node_pt(nnod-1);
386  adj_el_pt[nod_pt].insert(el_pt);
387  }
388 
389  // Find node that is only associated with single element
390  unsigned nfound=0;
391  Node* first_node_pt=0;
392  FiniteElement* first_element_pt=0;
393  for (std::map<Node*,std::set<FiniteElement*> >::iterator it=
394  adj_el_pt.begin();
395  it!=adj_el_pt.end();it++)
396  {
397  if ((*it).second.size()==1)
398  {
399  if (nfound==0)
400  {
401  first_node_pt=(*it).first;
402  first_element_pt=*(((*it).second).begin());
403  nfound++;
404  }
405  else
406  {
407  Node* other_node_pt=(*it).first;
408  if (other_node_pt->x(0)<=first_node_pt->x(0))
409  {
410  first_node_pt=other_node_pt;
411  first_element_pt=*(((*it).second).begin());
412  }
413  break;
414  }
415  }
416  }
417 
418 
419  // Current element/node
420  FiniteElement* current_element_pt=first_element_pt;
421  Node* current_node_pt=first_node_pt;
422 
423  oomph_info << "FIRST NODE: "
424  << first_node_pt->x(0) << " "
425  << first_node_pt->x(1) << " " << std::endl;
426 
427  // Add it
428  Ordered_shielding_node_pt.push_back(current_node_pt);
429  Ordered_shielding_face_element_pt.push_back(current_element_pt);
430 
431  // Keep going until we reach the end
432  unsigned n_associated_els=2;
433  while (n_associated_els==2)
434  {
435  // Next node
436  Node* next_node_pt=current_element_pt->node_pt(0);
437  if (next_node_pt==current_node_pt)
438  {
439  unsigned nnod=current_element_pt->nnode();
440  next_node_pt=current_element_pt->node_pt(nnod-1);
441  }
442 
443  // Get next element
444  std::set<FiniteElement*>::iterator it=adj_el_pt[next_node_pt].begin();
445  FiniteElement* next_element_pt=(*it);
446  if (next_element_pt==current_element_pt)
447  {
448  it++;
449  next_element_pt=(*it);
450  }
451 
452  // Add it
453  Ordered_shielding_node_pt.push_back(next_node_pt);
454  Ordered_shielding_face_element_pt.push_back(next_element_pt);
455 
456  // bump up
457  current_element_pt=next_element_pt;
458  current_node_pt=next_node_pt;
459 
460  n_associated_els=adj_el_pt[current_node_pt].size();
461  }
462 
463  // Kill last one
465 
466  oomph_info << "Done setup shielding nodes: "
467  << Ordered_shielding_node_pt.size() << std::endl;
468 
469  // Right horizon
470  // hierher Ordered_shielding_node_pt.push_back(Right_horizon_node_pt);
471  }
Vector< Node * > Ordered_shielding_node_pt
Storage for ordered shielding nodes.
Definition: circular_boulder_solar_radiation.cc:766
Definition: elements.h:1313
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060

References e(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::oomph_info, and oomph::Node::x().

◆ update_limiting_angles()

template<class ELEMENT >
void SolarRadiationProblem< ELEMENT >::update_limiting_angles ( )
inline

Update limiting angles for diffuse radiation.

343  {
344  // Update limiting angles for diffuse radiation, given the
345  // Vector of pointers to face elements that make up the "upper boundary"
346  // that can potentially shield the integration points from diffuse
347  // radiation
348  unsigned nel=Surface_flux_mesh_pt->nelement();
349  for (unsigned e=0;e<nel;e++)
350  {
354  }
355  }

References e().

Member Data Documentation

◆ Bottom_boundary_id

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

ID of bottom boundary.

◆ Bulk_mesh_pt

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

Pointer to bulk mesh.

◆ Doc_info

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

◆ Left_boundary_id

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

ID of left boundary.

◆ Left_melt_boundary_id

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

ID of left melt boundary.

◆ Lower_contact_boundary_id

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

ID of lower contact boundary.

◆ Ordered_shielding_face_element_pt

template<class ELEMENT >
Vector<FiniteElement*> SolarRadiationProblem< ELEMENT >::Ordered_shielding_face_element_pt
private

Face elements on potentially sun-exposed boundaries.

◆ Ordered_shielding_node_pt

template<class ELEMENT >
Vector<Node*> SolarRadiationProblem< ELEMENT >::Ordered_shielding_node_pt
private

Storage for ordered shielding nodes.

◆ Right_boundary_id

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

ID of right boundary.

◆ Right_melt_boundary_id

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

ID of right melt boundary.

◆ Surface_contact_mesh_pt

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

Pointer to the "surface" mesh.

◆ Surface_flux_mesh_pt

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

Pointer to the "surface" mesh.

◆ Trace_file

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

Trace file.

◆ Upper_contact_boundary_id

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

ID of upper contact boundary.


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