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
 
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_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...
 
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
 
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)
 
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 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: