AxisymmetricVibratingShellProblem< ELEMENT > Class Template Reference

Problem class to simulate the settling of a viscous drop. More...

+ Inheritance diagram for AxisymmetricVibratingShellProblem< ELEMENT >:

Public Member Functions

 AxisymmetricVibratingShellProblem ()
 Constructor. More...
 
 ~AxisymmetricVibratingShellProblem ()
 Destructor. More...
 
void actions_before_adapt ()
 Actions before adapt: Wipe the mesh of free surface elements. More...
 
void actions_after_adapt ()
 Actions after adapt: Rebuild the mesh of free surface elements. More...
 
void actions_after_newton_solve ()
 Update the after solve (empty) More...
 
void actions_before_newton_solve ()
 Update the problem specs before solve. More...
 
void actions_before_newton_convergence_check ()
 Allow inverted elements during Newton iteration. More...
 
void actions_after_newton_convergence_check ()
 Don't allow inverted elements in general. More...
 
void complete_problem_setup ()
 Set boundary conditions and complete the build of all elements. More...
 
void doc_solution (const std::string &comment="")
 Doc the solution. More...
 
void compute_error_estimate (double &max_err, double &min_err)
 Compute the error estimates and assign to elements for plotting. More...
 
double global_temporal_error_norm ()
 Error norm to determine the next time step. More...
 
double & next_dt ()
 Access the next suggested timestep. More...
 
double height_central_node ()
 Get the height at the centre node. More...
 
void update_latest_fixed_point_iteration_guess_for_strain_rate_for_all_elements ()
 Update latest guess for strain rate. More...
 
void enable_fixed_point_iteration_for_strain_rate_for_all_elements ()
 
void disable_fixed_point_iteration_for_strain_rate_for_all_elements ()
 Disable use of fixed point iteration. More...
 
void enable_aitken_extrapolation_for_all_elements ()
 Enable use of Aitken extrapolation for all elements. More...
 
void disable_aitken_extrapolation_for_all_elements ()
 Disable use of Aitken extrapolation. More...
 
void set_nprev_for_extrapolation_of_strain_rate_for_all_elements (const unsigned &nprev)
 
double calculate_error_of_fixed_point_iteration ()
 get the error of the fixed point iteration More...
 
double strainrate_norm ()
 
void dump_it (ofstream &dump_file)
 Dump problem data to allow for later restart. More...
 
void restart ()
 Restart. 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 Types

enum  { Internal_boundary1_id =0 , Free_surface_boundary_id =1 , Symmetry_line_id =2 , Shell_wall_boundary_id =3 }
 Enumeration of boundaries. More...
 

Private Member Functions

void create_free_surface_elements ()
 Create free surface elements. More...
 
void delete_free_surface_elements ()
 Delete free surface elements. More...
 

Private Attributes

Node * Central_node_on_free_surface_pt
 Pointer to free surface node on symmetry line. More...
 
double Next_dt
 
Mesh * Free_surface_mesh_pt
 Pointers to mesh of free surface elements. More...
 
RefineableSolidTriangleMesh< ELEMENT > * Fluid_mesh_pt
 Pointer to Fluid_mesh. More...
 
Data * External_pressure_data_pt
 Pointer to a global external pressure datum. More...
 
TriangleMeshClosedCurve * Outer_boundary_polyline_pt
 Triangle mesh polygon for outer boundary. More...
 
Vector< TriangleMeshOpenCurve * > Internal_open_boundary_pt
 
double Pext
 External pressure. 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_before_implicit_timestep ()
 
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 AxisymmetricVibratingShellProblem< ELEMENT >

Problem class to simulate the settling of a viscous drop.

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

Member Enumeration Documentation

◆ anonymous enum

template<class ELEMENT >
anonymous enum
private

Enumeration of boundaries.

Enumerator
Internal_boundary1_id 
Free_surface_boundary_id 
Symmetry_line_id 
Shell_wall_boundary_id 
704  {
709  };
@ Symmetry_line_id
Definition: axisym_vibrating_shell_non_newtonian.cc:707
@ Shell_wall_boundary_id
Definition: axisym_vibrating_shell_non_newtonian.cc:708
@ Free_surface_boundary_id
Definition: axisym_vibrating_shell_non_newtonian.cc:706
@ Internal_boundary1_id
Definition: axisym_vibrating_shell_non_newtonian.cc:705

Constructor & Destructor Documentation

◆ AxisymmetricVibratingShellProblem()

Constructor.

Get pointer to the free surface node

Get the node's boundary coordinate

Get the node's coordinate in the spline representation

Get pointer to the free surface node

Get the node's r-coordinate

Set external pressure to zero

721 {
722 
725  Max_residuals=10.0;
727 
728  //enable_globally_convergent_newton_method();
729 
730  // Allocate the timestepper -- this constructs the Problem's
731  // time object with a sufficient amount of storage to store the
732  // previous timsteps.
734  {
735  case 1:
736  this->add_time_stepper_pt(new BDF<1>(true));
737  oomph_info << "Using BDF1\n";
738  break;
739 
740  case 2:
741  this->add_time_stepper_pt(new BDF<2>(true));
742  oomph_info << "Using BDF2\n";
743  break;
744 
745  case 4:
746  this->add_time_stepper_pt(new BDF<4>(true));
747  oomph_info << "Using BDF4\n";
748  break;
749 
750  default:
751  oomph_info << "Wrong BDF type: " << Problem_Parameter::BDF_type
752  << std::endl;
753  break;
754  }
755 
756  // Build the boundary segments for outer boundary, consisting of
757  //--------------------------------------------------------------
758  // 2 separate polylines
759  //------------------------
760 
761  Vector<TriangleMeshCurveSection*> boundary_curve_section_pt(3);
762  Vector<TriangleMeshPolyLine*> internal_polyline_pt(2);
763 
764  // Each polyline only has two vertices -- provide storage for their
765  // coordinates
766  Vector<Vector<double> > vertex_coord;
767 
768  Ellipse* Shell_wall_pt = new Ellipse(Problem_Parameter::Shell_radius,
770 
771 
772  double zeta_start=0.0;
773  double zeta_end=MathematicalConstants::Pi/2.0;
774 
775  // number of points along the shell wall
776  unsigned npoints_wall=80;
777 
778  // get the increment in boundary coordinate
779  double unit_zeta = (zeta_end-zeta_start)/double(npoints_wall-1);
780 
781  // resize the vector storing the vertices
782  vertex_coord.resize(npoints_wall);
783 
784  // boundary coordinate
785  Vector<double> zeta(1,0.0);
786 
787  // coordinates of point on the boundary
788  Vector<double> coord(2,0.0);
789 
790  // Create points on boundary
791  for(unsigned ipoint=0; ipoint<npoints_wall;ipoint++)
792  {
793  // Get the coordinates
794  zeta[0]=zeta_start+unit_zeta*double(ipoint);
795  Shell_wall_pt->position(zeta,coord);
796 
797  // resize vertex coordinate
798  vertex_coord[ipoint].resize(2);
799 
800  // Shift
801  vertex_coord[ipoint][0]=coord[0];
802  //std::cout<<vertex_coord[ipoint][0]<< " ";
803  vertex_coord[ipoint][1]=coord[1];
804  //std::cout<<vertex_coord[ipoint][1]<< std::endl;
805  }
806 
807  boundary_curve_section_pt[0]=
808  new TriangleMeshPolyLine(vertex_coord, Shell_wall_boundary_id);
809 
810  ofstream file;
811  file.open((Problem_Parameter::Directory+"/boundary_section0.dat").c_str());
812  boundary_curve_section_pt[0]->output(file);
813  file.close();
814 
815  // Resize the vertex coordinates vector
816  vertex_coord.resize(2);
817 
818  for(unsigned i=0;i<2;i++)
819  {
820  vertex_coord[i].resize(2);
821  }
822 
823  // Build the symmetry polyline
824  vertex_coord[0][0]=0.0;
825  vertex_coord[0][1]=Problem_Parameter::Shell_radius;
826  vertex_coord[1][0]=0.0;
827  vertex_coord[1][1]=Problem_Parameter::free_surface_profile(0.0); // obacht
828  boundary_curve_section_pt[1]= new TriangleMeshPolyLine(vertex_coord,
830 
831  file.open((Problem_Parameter::Directory+"/boundary_section1.dat").c_str());
832  boundary_curve_section_pt[1]->output(file);
833  file.close();
834 
835  // number of points along the shell wall
836  unsigned npoints_fs=80;
837 
838  // get the increment in boundary coordinate
839  double increment = 1.0/double(npoints_fs-1);
840 
841  // resize the vector storing the vertices
842  vertex_coord.resize(npoints_fs);
843 
844  // Create points on boundary
845  for(unsigned ipoint=0; ipoint<npoints_fs;ipoint++)
846  {
847  // resize vertex coordinate
848  vertex_coord[ipoint].resize(2);
849 
850  if(ipoint==0)
851  {
852  vertex_coord[ipoint][0]=0.0;
853  //std::cout<<vertex_coord[ipoint][0]<< " ";
854  vertex_coord[ipoint][1]=
856  //std::cout<<vertex_coord[ipoint][1]<< std::endl;
857  continue;
858  }
859  else if(ipoint==npoints_fs-1)
860  {
861  vertex_coord[ipoint][0]=Problem_Parameter::Shell_radius;
862  //std::cout<<vertex_coord[ipoint][0]<< " ";
863  vertex_coord[ipoint][1]=0.0;
864  //std::cout<<vertex_coord[ipoint][1]<< std::endl;
865  continue;
866  }
867 
868  // Get the coordinates
869  coord[0]=double(ipoint)*increment;
870  coord[1]=Problem_Parameter::free_surface_profile(coord[0]);
871  //-1.0*(-0.05*tanh(10.0*coord[0]-7.5)+0.05); // obacht
872 
873  // Shift
874  vertex_coord[ipoint][0]=coord[0];
875  //std::cout<<vertex_coord[ipoint][0]<< " ";
876  vertex_coord[ipoint][1]=coord[1];
877  //std::cout<<vertex_coord[ipoint][1]<< std::endl;
878  }
879 
880  // Build the boundary polyline
881  boundary_curve_section_pt[2]=
883 
884  file.open((Problem_Parameter::Directory+"/boundary_section2.dat").c_str());
885  boundary_curve_section_pt[2]->output(file);
886  file.close();
887 
888  // Create the triangle mesh polygon for outer boundary
890  new TriangleMeshClosedCurve(boundary_curve_section_pt);
891 
892  // Use the TriangleMeshParameters object for gathering all
893  // the necessary arguments for the TriangleMesh object
894  TriangleMeshParameters triangle_mesh_parameters(
896 
897  // Define the maximum element areas
898  triangle_mesh_parameters.element_area() =
900 
901  // enable the boundary refinement
902  triangle_mesh_parameters.enable_boundary_refinement();
903 
904  // Create the mesh
905  Fluid_mesh_pt =
906  new RefineableSolidTriangleMesh<ELEMENT>(
907  triangle_mesh_parameters,this->time_stepper_pt());
908 
909 
910  // Increase number of bins to reduce bias/drift in area targets
911  Fluid_mesh_pt->nbin_x_for_area_transfer()=500;
912  Fluid_mesh_pt->nbin_y_for_area_transfer()=500;
913 
914 // Fluid_mesh_pt->output((Problem_Parameter::Directory+"/Mesh_before_snapping.dat").c_str());
915 
916  //----------------------------------------------------------------
917  // Snap nodes manually onto the curved boundary
918  //----------------------------------------------------------------
919 
920  // loop over the nodes on the wall boundary
921  unsigned n_boundary_node = Fluid_mesh_pt->
922  nboundary_node(Shell_wall_boundary_id);
923 
924  for(unsigned n=0;n<n_boundary_node;n++)
925  {
927  Node* nod_pt = Fluid_mesh_pt->boundary_node_pt(Shell_wall_boundary_id,n);
928 
930  Vector<double> zeta(1);
932 
934  Vector<double> new_x(2,0.0);
935  // updating zeta
936  zeta[0]=zeta_start+zeta[0]*(zeta_end-zeta_start);
937  Shell_wall_pt->position(zeta,new_x);
938  nod_pt->x(0) = new_x[0];
939  nod_pt->x(1) = new_x[1];
940 
941  }
942 
943  n_boundary_node=Fluid_mesh_pt->nboundary_node(Free_surface_boundary_id);
944 
945  for(unsigned n=0;n<n_boundary_node;n++)
946  {
948  Node* nod_pt = Fluid_mesh_pt->boundary_node_pt(Free_surface_boundary_id,n);
949 
952  {
953  continue;
954  }
955 
957  double r=nod_pt->x(0);
958 
959  // calculate node's z-coordinate
960  double new_z=Problem_Parameter::free_surface_profile(r); // obacht
961 
962  // move node
963  nod_pt->x(1) = new_z;
964 
965  }
966 
967  //-----------------------------------------------------------------
968  // End of snapping
969  //-----------------------------------------------------------------
970 
971  //Fluid_mesh_pt->output((Problem_Parameter::Directory+"/Mesh_after_snapping.dat").c_str());
972  Fluid_mesh_pt->output_boundaries((Problem_Parameter::Directory+"/boundaries.dat").c_str());
973 
974  // Set error estimator for bulk mesh
976  Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
977 
978  // Set targets for spatial adaptivity
979  Fluid_mesh_pt->max_permitted_error()=5.0e-3; // was 0.005; 0.0005
980  Fluid_mesh_pt->min_permitted_error()=1.0e-3; // was 0.001; 0.0001
982  Fluid_mesh_pt->min_element_size()=1.0e-5; // was 0.001 or 0.00001
983 
984 
986  Pext=0.0;
987 
988  //Create a Data object whose single value stores the
989  //external pressure
991 
992  // Set external pressure
994 
995  // The external pressure is pinned -- the external pressure
996  // sets the pressure throughout the domain -- we do not have
997  // the liberty to fix another pressure value!
999 
1000  // Construct the mesh of free surface elements
1003 
1004  // Set boundary condition and complete the build of all bulk elements
1006 
1007  // Combine meshes
1008  //---------------
1009 
1010  // Add Fluid_mesh_pt sub meshes
1011  this->add_sub_mesh(Fluid_mesh_pt);
1012 
1013  // Add Free_surface sub meshes
1014  this->add_sub_mesh(this->Free_surface_mesh_pt);
1015 
1016  // Build global mesh
1017  this->build_global_mesh();
1018 
1019  // Set lagrangian coordinates for pseudo-solid
1020  Fluid_mesh_pt->set_lagrangian_nodal_coordinates();
1021 
1022  // Use mumps
1023  //linear_solver_pt()=new MumpsSolver;
1024 
1025  // Setup equation numbering scheme
1026  cout <<"Number of equations: " << this->assign_eqn_numbers() << std::endl;
1027 
1028 
1029 } // end_of_constructor
int i
Definition: BiCGSTAB_step_by_step.cpp:9
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
Data * External_pressure_data_pt
Pointer to a global external pressure datum.
Definition: axisym_vibrating_shell_non_newtonian.cc:695
void create_free_surface_elements()
Create free surface elements.
Definition: axisym_vibrating_shell_non_newtonian.cc:1282
TriangleMeshClosedCurve * Outer_boundary_polyline_pt
Triangle mesh polygon for outer boundary.
Definition: axisym_vibrating_shell_non_newtonian.cc:698
RefineableSolidTriangleMesh< ELEMENT > * Fluid_mesh_pt
Pointer to Fluid_mesh.
Definition: axisym_vibrating_shell_non_newtonian.cc:692
void complete_problem_setup()
Set boundary conditions and complete the build of all elements.
Definition: axisym_vibrating_shell_non_newtonian.cc:1076
double Pext
External pressure.
Definition: axisym_vibrating_shell_non_newtonian.cc:712
Mesh * Free_surface_mesh_pt
Pointers to mesh of free surface elements.
Definition: axisym_vibrating_shell_non_newtonian.cc:689
Definition: timesteppers.h:1165
Definition: nodes.h:86
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
void set_value(const unsigned &i, const double &value_)
Definition: nodes.h:271
Definition: geom_objects.h:644
void position(const Vector< double > &zeta, Vector< double > &r) const
Position Vector at Lagrangian coordinate zeta.
Definition: geom_objects.h:745
Definition: mesh.h:67
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
virtual bool is_on_boundary() const
Definition: nodes.h:1373
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Definition: nodes.cc:2379
bool Always_take_one_newton_step
Definition: problem.h:2197
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
unsigned Max_newton_iterations
Maximum number of Newton iterations.
Definition: problem.h:599
void build_global_mesh()
Definition: problem.cc:1493
double Minimum_dt_but_still_proceed
Definition: problem.h:725
double Max_residuals
Definition: problem.h:610
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Definition: problem.cc:1989
TimeStepper *& time_stepper_pt()
Definition: problem.h:1524
Base class defining a closed curve for the Triangle mesh generation.
Definition: unstructured_two_d_mesh_geometry_base.h:1339
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
Definition: oomph-lib/src/generic/Vector.h:58
Definition: error_estimator.h:266
double Pi
Definition: two_d_biharmonic.cc:235
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
Z2ErrorEstimator * error_estimator_pt
Definition: MortaringCantileverCompareToNonMortaring.cpp:190
r
Definition: UniformPSDSelfTest.py:20
std::string Directory
Output directory.
Definition: axisym_vibrating_shell_non_newtonian.cc:64
double Uniform_element_area
Uniform initial element area.
Definition: axisym_vibrating_shell_non_newtonian.cc:61
double Shell_radius
Radius of the shell.
Definition: axisym_vibrating_shell_non_newtonian.cc:143
unsigned BDF_type
Which bdf timestepper do we use?
Definition: axisym_vibrating_shell_non_newtonian.cc:54
double free_surface_profile(const double r)
Definition: axisym_vibrating_shell_non_newtonian.cc:128
double Dt_min
Minimum allowable timestep.
Definition: axisym_vibrating_shell_non_newtonian.cc:103
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

References oomph::Problem_Parameter::BDF_type, oomph::Problem_Parameter::Directory, oomph::Problem_Parameter::Dt_min, oomph::TriangleMeshParameters::element_area(), oomph::TriangleMeshParameters::enable_boundary_refinement(), MeshRefinement::error_estimator_pt, oomph::Problem_Parameter::free_surface_profile(), oomph::Node::get_coordinates_on_boundary(), i, oomph::Node::is_on_boundary(), oomph::Locate_zeta_helpers::Max_newton_iterations, n, oomph::oomph_info, Global_Physical_Variables::Pext, BiharmonicTestFunctions2::Pi, oomph::Ellipse::position(), UniformPSDSelfTest::r, oomph::Problem_Parameter::Shell_radius, oomph::Problem_Parameter::Uniform_element_area, oomph::Node::x(), and Eigen::zeta().

◆ ~AxisymmetricVibratingShellProblem()

Destructor.

1037 {
1038  // Fluid timestepper
1039  delete this->time_stepper_pt(0);
1040 
1041  // Kill data associated with outer boundary
1043  for (unsigned j=0;j<n;j++)
1044  {
1046  }
1048 
1049 
1050 
1051  // Flush element of free surface elements
1053  delete Free_surface_mesh_pt;
1054 
1055  // Delete error estimator
1056  delete Fluid_mesh_pt->spatial_error_estimator_pt();
1057 
1058  // Delete fluid mesh
1059  delete Fluid_mesh_pt;
1060 
1061  // Delete pressure data
1063 
1064  // Kill const eqn
1066 
1068 
1069 } // end_of_destructor
void delete_free_surface_elements()
Delete free surface elements.
Definition: axisym_vibrating_shell_non_newtonian.cc:671
virtual TriangleMeshCurveSection * curve_section_pt(const unsigned &i) const
Pointer to i-th constituent curve section.
Definition: unstructured_two_d_mesh_geometry_base.h:1308
virtual unsigned ncurve_section() const
Number of constituent curves.
Definition: unstructured_two_d_mesh_geometry_base.h:1172
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
Definition: jeffery_orbit.cc:82
GeneralisedNewtonianConstitutiveEquation< 3 > * Const_eqn_pt
Fluid constitutive equation.
Definition: axisym_vibrating_shell_non_newtonian.cc:213
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::Problem_Parameter::Const_eqn_pt, Problem_Parameter::Constitutive_law_pt, j, and n.

Member Function Documentation

◆ actions_after_adapt()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::actions_after_adapt ( )
inlinevirtual

Actions after adapt: Rebuild the mesh of free surface elements.

Reimplemented from oomph::Problem.

306  {
307 
308  // Create the elements that impose the displacement constraint
310 
311  // Rebuild the Problem's global mesh from its various sub-meshes
312  this->rebuild_global_mesh();
313 
314  // Setup the problem again -- remember that fluid mesh has been
315  // completely rebuilt and its element's don't have any
316  // pointers to Re etc. yet
318 
319  }// end of actions_after_adapt
void rebuild_global_mesh()
Definition: problem.cc:1533

◆ actions_after_newton_convergence_check()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::actions_after_newton_convergence_check ( )
inline

Don't allow inverted elements in general.

352  {
353  unsigned n_element = Fluid_mesh_pt->nelement();
354  for(unsigned e=0;e<n_element;e++)
355  {
356  // Upcast from GeneralisedElement to the present element
357  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
358 
359  el_pt->Accept_negative_jacobian=false;
360  }
361 
362  }
Array< double, 1, 3 > e(1./3., 0.5, 2.)

References e().

◆ actions_after_newton_solve()

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

Update the after solve (empty)

Reimplemented from oomph::Problem.

323 {}

◆ actions_before_adapt()

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

Actions before adapt: Wipe the mesh of free surface elements.

Reimplemented from oomph::Problem.

294  {
295  // Kill the elements and wipe surface mesh
297 
298  // Rebuild the Problem's global mesh from its various sub-meshes
299  this->rebuild_global_mesh();
300 
301  }// end of actions_before_adapt

◆ actions_before_newton_convergence_check()

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

Allow inverted elements during Newton iteration.

Reimplemented from oomph::Problem.

337  {
338  unsigned n_element = Fluid_mesh_pt->nelement();
339  for(unsigned e=0;e<n_element;e++)
340  {
341  // Upcast from GeneralisedElement to the present element
342  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
343 
344  // accept negative jacobian in solution process
345  el_pt->Accept_negative_jacobian=true;
346  }
347 
348  }

References e().

◆ actions_before_newton_solve()

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

Update the problem specs before solve.

Reimplemented from oomph::Problem.

327  {
328  //Reset the Lagrangian coordinates of the nodes to be the current
329  //Eulerian coordinates -- this makes the current configuration
330  //stress free
331  Fluid_mesh_pt->set_lagrangian_nodal_coordinates();
332  }

◆ calculate_error_of_fixed_point_iteration()

template<class ELEMENT >
double AxisymmetricVibratingShellProblem< ELEMENT >::calculate_error_of_fixed_point_iteration ( )
inline

get the error of the fixed point iteration

474  {
475  // Get elemental max/min/norms
476  double norm_squared=0.0;
477  double latest_guess_norm_squared=0.0;
478  double error_norm_squared=0.0;
479 
480  unsigned nel=Fluid_mesh_pt->nelement();
481  for (unsigned e=0;e<nel;e++)
482  {
483  // Get element
484  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(
485  Fluid_mesh_pt->element_pt(e));
486 
487  // Get norms of invariant
488  double el_norm_squared=0.0;
489  double el_latest_guess_norm_squared=0.0;
490  double el_error_norm_squared=0.0;
491  el_pt->square_of_norm_of_fixed_point(
492  el_norm_squared,
493  el_latest_guess_norm_squared,
494  el_error_norm_squared);
495 
496  // Add it...
497  norm_squared+=el_norm_squared;
498  latest_guess_norm_squared+=el_latest_guess_norm_squared;
499  error_norm_squared+=el_error_norm_squared;
500  }
501 
502  oomph_info << "Norm of current strain rate invariant: "
503  << sqrt(norm_squared) << std::endl;
504  oomph_info << "Norm of latest fixed point iteration guess for "
505  << "strain rate invariant: "
506  << sqrt(latest_guess_norm_squared) << std::endl;
507  oomph_info << "Norm of error in fixed point iteration "
508  << "strain rate invariant: "
509  << sqrt(error_norm_squared) << " equivalent to " ;
510  if (sqrt(norm_squared)!=0.0)
511  {
512  oomph_info << sqrt(error_norm_squared)/sqrt(norm_squared)*100.0 << " %";
513  }
514  oomph_info << std::endl;
515 
516  if (sqrt(norm_squared)!=0.0)
517  {
518  return sqrt(error_norm_squared)/sqrt(norm_squared)*100.0;
519  }
520  else
521  {
522  return 0.0;
523  }
524 
525  }
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134

References e(), oomph::oomph_info, and sqrt().

◆ complete_problem_setup()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::complete_problem_setup

Set boundary conditions and complete the build of all elements.

Set up the problem: apply BC and make bulk elements fully functional.

Get pointer to the free surface node

1077 {
1078  // Re-set the boundary conditions for fluid problem: All nodes are
1079  // free by default -- just pin the ones that have Dirichlet conditions
1080  // here.
1081  unsigned nbound=Fluid_mesh_pt->nboundary();
1082  for(unsigned ibound=0;ibound<nbound;ibound++)
1083  {
1084  unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
1085  for (unsigned inod=0;inod<num_nod;inod++)
1086  {
1087  // Get node
1088  Node* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
1089 
1090  //Pin both velocities at the shell wall boundary
1091  if(ibound==Shell_wall_boundary_id)
1092  {
1093  nod_pt->pin(0);
1094  nod_pt->pin(1);
1095  }
1096 
1097  // pin horizontal velocity at symmetry boundary
1098  if(ibound==Symmetry_line_id)
1099  {
1100  nod_pt->pin(0);
1101  }
1102 
1103  // pin Lagrange multiplier at the intersection of the shell wall
1104  // boundary and the free surface
1105  if( (nod_pt->is_on_boundary(Shell_wall_boundary_id)) &&
1107  {
1108  // Get the number of values at this node
1109  unsigned n_value=nod_pt->nvalue();
1110 
1111  // check that it is the corner node with 5 values
1112  // (u,v,w,p,delta)
1113  if(n_value != 5)
1114  {
1115  oomph_info <<" Here!\n";
1116  abort();
1117  }
1118 
1119  nod_pt->pin(n_value-1);
1120  }
1121 
1122  // Pin pseudo-solid positions apart from free surface boundary which
1123  // we allow to move
1124  SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
1125  if(ibound==Shell_wall_boundary_id)
1126  {
1127  solid_node_pt->pin_position(0);
1128  solid_node_pt->pin_position(1);
1129  }
1130  else if(ibound==Symmetry_line_id)
1131  {
1132  solid_node_pt->pin_position(0);
1133  }
1134 
1135  }
1136 
1137  } // end loop over boundaries
1138 
1139  // Complete the build of all elements so they are fully functional
1140  // Remember that adaptation for triangle meshes involves a complete
1141  // regneration of the mesh (rather than splitting as in tree-based
1142  // meshes where such parameters can be passed down from the father
1143  // element!)
1144  unsigned n_element = Fluid_mesh_pt->nelement();
1145  for(unsigned e=0;e<n_element;e++)
1146  {
1147  // Upcast from GeneralisedElement to the present element
1148  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
1149 
1150  // Modify the tolerance for when the mapping is considered singular
1151  el_pt->Tolerance_for_singular_jacobian=1.0e-26;
1152 
1153  // Set the Reynolds number
1154  el_pt->re_pt() = &Problem_Parameter::Re;
1155 
1156  // Set the Womersley number
1157  el_pt->re_st_pt() = &Problem_Parameter::Re_St;
1158 
1159  // Set the body force
1160  el_pt->axi_nst_body_force_fct_pt() = &Problem_Parameter::oscillation;
1161 
1162  // Set the product of Reynolds number and inverse Froude number
1163  el_pt->re_invfr_pt() = &Problem_Parameter::ReInvFr;
1164 
1165  // Set the constitutive law for pseudo-elastic mesh deformation
1166  el_pt->constitutive_law_pt()=Problem_Parameter::Constitutive_law_pt;
1167 
1168  // Assign Constitutive equation
1169  el_pt->constitutive_eqn_pt() = Problem_Parameter::Const_eqn_pt;
1170 
1171  // Use extrapolated strain rate when determining viscosity
1172  el_pt->use_extrapolated_strainrate_to_compute_second_invariant();
1173 
1174  // for all nodes, pin azimuthal velocity (fewer dofs)
1175  unsigned n_node=el_pt->nnode();
1176  for(unsigned inod=0;inod<n_node;inod++)
1177  {
1178  el_pt->node_pt(inod)->pin(2);
1179  el_pt->node_pt(inod)->set_value(2,0.0);
1180  }
1181 
1182  } // end of functional elements
1183 
1184 
1185  // Setup extrapolation
1188 
1189 
1190  // Re-apply boundary values on Dirichlet boundary conditions
1191  // (Boundary conditions are ignored when the solution is transferred
1192  // from the old to the new mesh by projection; this leads to a slight
1193  // change in the boundary values (which are, of course, never changed,
1194  // unlike the actual unknowns for which the projected values only
1195  // serve as an initial guess)
1196 
1197  // Set velocity and history values of velocity on walls
1198  nbound=this->Fluid_mesh_pt->nboundary();
1199  for(unsigned ibound=0;ibound<nbound;++ibound)
1200  {
1201  if( (ibound==Shell_wall_boundary_id) ||
1202  (ibound==Symmetry_line_id) )
1203  {
1204  // Loop over nodes on this boundary
1205  unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(ibound);
1206  for (unsigned inod=0;inod<num_nod;inod++)
1207  {
1208  // Get node
1209  Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(ibound,inod);
1210 
1211  // Get number of previous (history) values
1212  unsigned n_prev=nod_pt->time_stepper_pt()->nprev_values();
1213 
1214  // Velocity is and was zero at all previous times
1215  for (unsigned t=0;t<=n_prev;t++)
1216  {
1217  if(ibound==Shell_wall_boundary_id)
1218  {
1219  nod_pt->set_value(t,1,0.0);
1220  nod_pt->set_value(t,0,0.0);
1221  }
1222  else
1223  {
1224  nod_pt->set_value(t,0,0.0);
1225  }
1226 
1227  // Move nodes on symmetry line exactly onto r=0 (for all times)
1228  if(ibound==Symmetry_line_id)
1229  {
1230  nod_pt->x(t,0)=0.0;
1231  }
1232  }
1233  }
1234  }
1235  }
1236 
1237 
1238  // Update pointer to central node
1239  oomph_info << "Updating central node" << std::endl;
1241  const unsigned n_boundary_node = Fluid_mesh_pt->
1242  nboundary_node(Free_surface_boundary_id);
1243  for(unsigned n=0;n<n_boundary_node;n++)
1244  {
1246  Node* nod_pt = Fluid_mesh_pt->boundary_node_pt(Free_surface_boundary_id,n);
1247 
1248  // Is this the one?
1249  if (nod_pt->is_on_boundary(Symmetry_line_id))
1250  {
1252  {
1253  oomph_info << "Odd -- more than one free surface node on sym line?\n";
1254  abort();
1255  }
1257  }
1258  }
1259 
1261  {
1262  oomph_info << "Odd -- not found the free surface node on sym line...\n";
1263  abort();
1264  }
1265  else
1266  {
1267  oomph_info << "Updated central node" << std::endl;
1268  oomph_info << "Central node now at: "
1269  << Central_node_on_free_surface_pt->x(0) << " "
1270  << Central_node_on_free_surface_pt->x(1) << " "
1271  << std::endl;
1272  }
1273 
1274 } // end of complete_problem_setup
void set_nprev_for_extrapolation_of_strain_rate_for_all_elements(const unsigned &nprev)
Definition: axisym_vibrating_shell_non_newtonian.cc:447
Node * Central_node_on_free_surface_pt
Pointer to free surface node on symmetry line.
Definition: axisym_vibrating_shell_non_newtonian.cc:661
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
Definition: nodes.h:1686
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1816
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
double ReInvFr
Product of Reynolds number and inverse of Froude number.
Definition: refineable_two_layer_interface.cc:303
double Re
Reynolds number.
Definition: jeffery_orbit.cc:59
void oscillation(const double &time, const Vector< double > &x, Vector< double > &force)
Function describing the oscillation.
Definition: axisym_vibrating_shell_non_newtonian.cc:106
unsigned Nprev_for_extrapolation_of_strain_rate
Definition: axisym_vibrating_shell_non_newtonian.cc:58
double Re_St
Reynolds multiplied Strouhal.
Definition: axisym_vibrating_shell_non_newtonian.cc:91
t
Definition: plotPSD.py:36

References oomph::Problem_Parameter::Const_eqn_pt, Problem_Parameter::Constitutive_law_pt, e(), oomph::Node::is_on_boundary(), n, oomph::Problem_Parameter::Nprev_for_extrapolation_of_strain_rate, oomph::TimeStepper::nprev_values(), oomph::Data::nvalue(), oomph::oomph_info, oomph::Problem_Parameter::oscillation(), oomph::Data::pin(), oomph::SolidNode::pin_position(), Problem_Parameter::Re, oomph::Problem_Parameter::Re_St, Problem_Parameter::ReInvFr, oomph::Data::set_value(), plotPSD::t, oomph::Data::time_stepper_pt(), and oomph::Node::x().

◆ compute_error_estimate()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::compute_error_estimate ( double &  max_err,
double &  min_err 
)

Compute the error estimates and assign to elements for plotting.

Compute error estimates and assign to elements for plotting.

1657 {
1658  // Get error estimator
1659  ErrorEstimator* err_est_pt=Fluid_mesh_pt->spatial_error_estimator_pt();
1660 
1661  // Get/output error estimates
1662  unsigned nel=Fluid_mesh_pt->nelement();
1663  Vector<double> elemental_error(nel);
1664 
1665  // We need a dynamic cast, get_element_errors from the Fluid_mesh_pt
1666  // Dynamic cast is used because get_element_errors require a Mesh* ans
1667  // not a SolidMesh*
1668  Mesh* fluid_mesh_pt=dynamic_cast<Mesh*>(Fluid_mesh_pt);
1669  err_est_pt->get_element_errors(fluid_mesh_pt,
1670  elemental_error);
1671 
1672  // Set errors for post-processing and find extrema
1673  max_err=0.0;
1674  min_err=DBL_MAX;
1675  for (unsigned e=0;e<nel;e++)
1676  {
1677  dynamic_cast<MyTaylorHoodElement*>(Fluid_mesh_pt->element_pt(e))->
1678  set_error(elemental_error[e]);
1679 
1680  max_err=std::max(max_err,elemental_error[e]);
1681  min_err=std::min(min_err,elemental_error[e]);
1682  }
1683 
1684 }
Base class for spatial error estimators.
Definition: error_estimator.h:40
void get_element_errors(Mesh *&mesh_pt, Vector< double > &elemental_error)
Definition: error_estimator.h:56
Overload TaylorHood element to modify output.
Definition: pressure_driven_torus.cc:98
#define min(a, b)
Definition: datatypes.h:22
#define max(a, b)
Definition: datatypes.h:23
void set_error(const double &error)
Set error value for post-processing.
Definition: overloaded_element_body.h:432

References e(), oomph::ErrorEstimator::get_element_errors(), max, min, and set_error().

◆ create_free_surface_elements()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::create_free_surface_elements
private

Create free surface elements.

Create elements that impose the kinematic and dynamic bcs for the pseudo-solid fluid mesh

1283 {
1284  // How many bulk fluid elements are adjacent to boundary b?
1285  unsigned n_element =
1286  Fluid_mesh_pt->nboundary_element(Free_surface_boundary_id);
1287 
1288  // Loop over the bulk fluid elements adjacent to boundary b?
1289  for(unsigned e=0;e<n_element;e++)
1290  {
1291  // Get pointer to the bulk fluid element that is
1292  // adjacent to boundary b
1293  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
1294  Fluid_mesh_pt->boundary_element_pt(Free_surface_boundary_id,e));
1295 
1296  //Find the index of the face of element e along boundary b
1297  int face_index =
1298  Fluid_mesh_pt->face_index_at_boundary(Free_surface_boundary_id,e);
1299 
1300  // Create new element
1303  bulk_elem_pt,face_index);
1304 
1305  // Add it to the mesh
1307 
1308  //Add the appropriate boundary number
1310 
1311  //Specify the capillary number
1312  el_pt->ca_pt() = &Problem_Parameter::Ca;
1313 
1314  //Specify the Strouhal number
1315  el_pt->st_pt() = &Problem_Parameter::St;
1316 
1317  // Specify the bubble pressure (pointer to Data object and
1318  // index of value within that Data object that corresponds
1319  // to the traded pressure
1321 
1322  }
1323 
1324 } // end of create_free_surface_elements
Specialise the Elastic update case to axisymmetric equations.
Definition: specific_node_update_interface_elements.h:1257
void set_boundary_number_in_bulk_mesh(const unsigned &b)
Set function for the boundary number in bulk mesh.
Definition: elements.h:4482
double *& ca_pt()
Pointer to the Capillary number.
Definition: interface_elements.h:492
double *& st_pt()
The pointer to the Strouhal number.
Definition: interface_elements.h:504
void set_external_pressure_data(Data *external_pressure_data_pt)
Definition: interface_elements.h:539
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:617
double St
Strouhal number.
Definition: jeffery_orbit.cc:62
double Ca
Capillary number.
Definition: refineable_two_layer_interface.cc:314

References Problem_Parameter::Ca, oomph::FluidInterfaceElement::ca_pt(), e(), oomph::FaceElement::set_boundary_number_in_bulk_mesh(), oomph::FluidInterfaceElement::set_external_pressure_data(), Problem_Parameter::St, and oomph::FluidInterfaceElement::st_pt().

◆ delete_free_surface_elements()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::delete_free_surface_elements ( )
inlineprivate

Delete free surface elements.

672  {
673  // How many surface elements are in the surface mesh
674  unsigned n_element = Free_surface_mesh_pt->nelement();
675 
676  // Loop over the surface elements
677  for(unsigned e=0;e<n_element;e++)
678  {
679  // Kill surface element
681  }
682 
683  // Wipe the mesh
685 
686  } // end of delete_free_surface_elements
void flush_element_and_node_storage()
Definition: mesh.h:407
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

References e().

◆ disable_aitken_extrapolation_for_all_elements()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::disable_aitken_extrapolation_for_all_elements ( )
inline

Disable use of Aitken extrapolation.

436  {
437  unsigned n_element = Fluid_mesh_pt->nelement();
438  for(unsigned e=0;e<n_element;e++)
439  {
440  // Upcast from GeneralisedElement to the present element
441  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
442  el_pt->disable_aitken_extrapolation();
443  }
444  }

References e().

◆ disable_fixed_point_iteration_for_strain_rate_for_all_elements()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::disable_fixed_point_iteration_for_strain_rate_for_all_elements ( )
inline

Disable use of fixed point iteration.

412  {
413  unsigned n_element = Fluid_mesh_pt->nelement();
414  for(unsigned e=0;e<n_element;e++)
415  {
416  // Upcast from GeneralisedElement to the present element
417  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
418  el_pt->disable_fixed_point_iteration_for_strain_rate();
419  }
420  }

References e().

◆ doc_solution()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::doc_solution ( const std::string &  comment = "")

Doc the solution.

1390 {
1391 
1392  ofstream some_file;
1393  char filename[1000];
1394 
1395  oomph_info << "Docing trace step: "
1396  << Problem_Parameter::Doc_info_trace.number() << std::endl;
1397 
1398  // Compute errors and assign to each element for plotting
1399  double max_err=0.0;
1400  double min_err=0.0;
1401  compute_error_estimate(max_err,min_err);
1402 
1403  // Write restart file
1405  ("--suppress_restart_files") &&
1407  ("--validation") )
1408  {
1409  sprintf(filename,"%s/restart%i.dat",
1410  Problem_Parameter::Doc_info_soln.directory().c_str(),
1412  ofstream dump_file;
1413  dump_file.open(filename);
1414  dump_file.precision(20);
1415  dump_it(dump_file);
1416  dump_file.close();
1417  }
1418 
1419  // Get body force
1420  Vector<double> body_force(3,0.0);
1421  Vector<double> x(3,0.0);
1423 
1424  // only output the actual solution when it's the right time
1425  // hierher if(this->time_pt()->time() >= double(Count_doc)*Dt_doc)
1426  {
1427  oomph_info << "Docing soln step: "
1428  << Problem_Parameter::Doc_info_soln.number() << std::endl;
1429 
1430  // Number of plot points
1431  unsigned npts=5;
1432 
1433 
1434  // Actual solution
1435  sprintf(filename,"%s/soln%i.dat",
1436  Problem_Parameter::Doc_info_soln.directory().c_str(),
1438  some_file.open(filename);
1439  some_file << dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(0))
1440  ->variable_identifier();
1441  this->Fluid_mesh_pt->output(some_file,npts);
1442  some_file.close();
1443 
1444 
1445  // Actual solution on the free surface
1446  sprintf(filename,"%s/free_surface_soln%i.dat",
1447  Problem_Parameter::Doc_info_soln.directory().c_str(),
1449  some_file.open(filename);
1450  unsigned nel=Free_surface_mesh_pt->nelement();
1451  for (unsigned e=0;e<nel;e++)
1452  {
1454  Free_surface_mesh_pt->element_pt(e))->output(some_file,npts);
1455  }
1456  some_file.close();
1457 
1458 
1459  // Output free surface
1460  sprintf(filename,"%s/free_surface%i.dat",
1461  Problem_Parameter::Doc_info_soln.directory().c_str(),
1463  some_file.open(filename);
1464  unsigned num_nod=Fluid_mesh_pt->nboundary_node(Free_surface_boundary_id);
1465  if (num_nod>0)
1466  {
1467  for (unsigned inod=0;inod<num_nod;inod++)
1468  {
1469  Fluid_mesh_pt->boundary_node_pt(Free_surface_boundary_id, inod)
1470  ->output(some_file);
1471  }
1472  }
1473  some_file.close();
1474 
1475 
1476  // Coarse solution (mesh)
1477  unsigned npts_coarse=2;
1478  sprintf(filename,"%s/coarse_soln%i.dat",
1479  Problem_Parameter::Doc_info_soln.directory().c_str(),
1481  some_file.open(filename);
1482  some_file << dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(0))
1483  ->variable_identifier();
1484  this->Fluid_mesh_pt->output(some_file,npts_coarse);
1485  some_file.close();
1486 
1487  // Doc body force
1488  sprintf(filename,"%s/body_force%i.dat",
1489  Problem_Parameter::Doc_info_soln.directory().c_str(),
1491  some_file.open(filename);
1492  some_file << "-0.5 0.0 0.0 "
1493  << body_force[0] << " "
1494  << body_force[1] << " "
1495  << body_force[2] << " "
1496  << std::endl;
1497  some_file.close();
1498 
1499 
1500  // Increment the doc_info number
1502 
1503  }
1504 
1505  // Assemble square of L2 norm
1506  double square_of_l2_norm=0.0;
1507  unsigned nel=Fluid_mesh_pt->nelement();
1508  for (unsigned e=0;e<nel;e++)
1509  {
1511  dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(e))->
1513  }
1514  Problem_Parameter::Norm_file << sqrt(square_of_l2_norm) << std::endl;
1515 
1516 
1517  // Output boundaries
1518  sprintf(filename,"%s/boundaries%i.dat",
1519  Problem_Parameter::Doc_info_soln.directory().c_str(),
1521  some_file.open(filename);
1522  this->Fluid_mesh_pt->output_boundaries(some_file);
1523  some_file.close();
1524 
1525  // Get max/min area
1526  double max_area=0.0;
1527  double min_area=0.0;
1528  Fluid_mesh_pt->max_and_min_element_size(max_area, min_area);
1529 
1530 
1531  // Compute current volume and error in extrapolation from fluid elements
1532  double current_vol=0.0;
1533  double norm_squared=0.0;
1534  double visc_norm_squared=0.0;
1535  double extrapolated_norm_squared=0.0;
1536  double error_norm_squared=0.0;
1537  double min_invariant=DBL_MAX;
1538  double max_invariant=-DBL_MAX;
1539  double min_viscosity=DBL_MAX;
1540  double max_viscosity=-DBL_MAX;
1541  for (unsigned e=0;e<nel;e++)
1542  {
1543  // Get element
1544  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(
1545  Fluid_mesh_pt->element_pt(e));
1546 
1547  // Add to physical size (actual volume)
1548  current_vol+=el_pt->compute_physical_size();
1549 
1550  // Check size (in computational coordinates
1551  double size=el_pt->size();
1552 
1553  // Get norms of invariant
1554  double el_norm_squared=0.0;
1555  double el_extrapolated_norm_squared=0.0;
1556  double el_error_norm_squared=0.0;
1557  double test_size=0.0;
1558  test_size=el_pt->square_of_norm_of_strain_invariant(
1559  el_norm_squared,
1560  el_extrapolated_norm_squared,
1561  el_error_norm_squared);
1562 
1563  if (std::fabs(test_size-size)>1.0e-10)
1564  {
1565  oomph_info << "Trouble: "
1566  << test_size << " "
1567  << size << std::endl;
1568  }
1569 
1570  // Get norms of viscosity
1571  double el_visc_norm_squared=0.0;
1572  double el_visc_extrapolated_norm_squared=0.0;
1573  double el_visc_error_norm_squared=0.0;
1574  el_pt->square_of_norm_of_viscosity(
1575  el_visc_norm_squared,
1576  el_visc_extrapolated_norm_squared,
1577  el_visc_error_norm_squared);
1578 
1579  // Add it...
1580  norm_squared+=el_norm_squared;
1581  visc_norm_squared+=el_visc_norm_squared;
1582  extrapolated_norm_squared+=el_extrapolated_norm_squared;
1583  error_norm_squared+=el_error_norm_squared;
1584 
1585  // Get viscosity extrema
1586  double el_min_invariant=0.0;
1587  double el_max_invariant=0.0;
1588  double el_min_viscosity=0.0;
1589  double el_max_viscosity=0.0;
1590  el_pt->max_and_min_invariant_and_viscosity(el_min_invariant,
1591  el_max_invariant,
1592  el_min_viscosity,
1593  el_max_viscosity);
1594 
1595  // Update overall extrema
1596  min_invariant=std::min(min_invariant,el_min_invariant);
1597  max_invariant=std::max(max_invariant,el_max_invariant);
1598  min_viscosity=std::min(min_viscosity,el_min_viscosity);
1599  max_viscosity=std::max(max_viscosity,el_max_viscosity);
1600  }
1601 
1602  oomph_info << "Norm of strain rate invariant: "
1603  << sqrt(norm_squared) << std::endl;
1604  oomph_info << "Norm of extrapolated strain rate invariant: "
1605  << sqrt(extrapolated_norm_squared) << std::endl;
1606  oomph_info << "Norm of error in extrapolated strain rate invariant: "
1607  << sqrt(error_norm_squared) << " equivalent to " ;
1608  if (sqrt(norm_squared)!=0.0)
1609  {
1610  oomph_info << sqrt(error_norm_squared)/sqrt(norm_squared)*100.0 << " %";
1611  }
1612  oomph_info << std::endl;
1613 
1614 
1615  oomph_info << "min_invariant = " << min_invariant << "\n"
1616  << "max_invariant = " << max_invariant << "\n"
1617  << "min_viscosity = " << min_viscosity << "\n"
1618  << "max_viscosity = " << max_viscosity
1619  << std::endl;
1620 
1621  // Write trace file
1623  << this->time_pt()->time() << " " // 1
1624  << Central_node_on_free_surface_pt->x(1) << " " // 2
1625  << body_force[1] << " " // 3
1626  << current_vol << " " // 4
1627  << Fluid_mesh_pt->nelement() << " " // 5
1628  << max_area << " " // 6
1629  << min_area << " " // 7
1630  << max_err << " " // spatial error 8
1631  << min_err << " " // spatial error 9
1632  << sqrt(norm_squared) << " " // strain invariant 10
1633  << sqrt(extrapolated_norm_squared) << " " // strain invariant 11
1634  << sqrt(error_norm_squared) << " " // strain invariant 12
1635  << min_invariant << " " // 13
1636  << max_invariant << " " // 14
1637  << min_viscosity << " " // 15
1638  << max_viscosity << " " // 16
1639  << max_viscosity << " " // 17
1640  << this->time_pt()->dt() << " " // 18
1641  << global_temporal_error_norm() << " " // temporal error measure 19
1642  << Problem_Parameter::Doc_info_trace.number() << " " // 20
1643  << std::endl;
1644 
1645  // Increment the doc_info number
1647 
1648 } // end_of_doc_full_solution
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
double global_temporal_error_norm()
Error norm to determine the next time step.
Definition: axisym_vibrating_shell_non_newtonian.cc:1331
void dump_it(ofstream &dump_file)
Dump problem data to allow for later restart.
Definition: axisym_vibrating_shell_non_newtonian.cc:561
void compute_error_estimate(double &max_err, double &min_err)
Compute the error estimates and assign to elements for plotting.
Definition: axisym_vibrating_shell_non_newtonian.cc:1656
unsigned & number()
Number used (e.g.) for labeling output files.
Definition: oomph_utilities.h:554
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 & dt(const unsigned &t=0)
Definition: timesteppers.h:136
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Definition: axisym_linear_elasticity/cylinder/cylinder.cc:96
string filename
Definition: MergeRestartFiles.py:39
ofstream Norm_file
Definition: refineable_two_layer_interface.cc:341
ofstream Trace_file
Trace file.
Definition: refineable_two_layer_interface.cc:335
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
bool command_line_flag_has_been_set(const std::string &flag)
Definition: oomph_utilities.cc:501
DocInfo Doc_info_soln
Doc info solutions.
Definition: axisym_vibrating_shell_non_newtonian.cc:73
DocInfo Doc_info_trace
Doc info trace file.
Definition: axisym_vibrating_shell_non_newtonian.cc:70
list x
Definition: plotDoE.py:28
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.
Definition: overloaded_element_body.h:490
double square_of_l2_norm()
Get square of L2 norm of velocity.
Definition: overloaded_element_body.h:1031

References Global_Parameters::body_force(), oomph::CommandLineArgs::command_line_flag_has_been_set(), oomph::Problem_Parameter::Doc_info_soln, oomph::Problem_Parameter::Doc_info_trace, e(), boost::multiprecision::fabs(), MergeRestartFiles::filename, max, min, Problem_Parameter::Norm_file, oomph::DocInfo::number(), oomph::oomph_info, oomph::Problem_Parameter::oscillation(), output(), size, sqrt(), square_of_l2_norm(), Problem_Parameter::Trace_file, and plotDoE::x.

◆ dump_it()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::dump_it ( ofstream &  dump_file)
inline

Dump problem data to allow for later restart.

562  {
563 
564  // Write doc numbers
566  << " # current doc number for trace" << std::endl;
568  << " # current doc number for soln" << std::endl;
569 
570  // Next timestep (required for restart from temporally adaptive
571  // run
572  dump_file << Next_dt << " # next timestep " << std::endl;
573 
574  // Dump the refinement pattern and the generic problem data
575  Problem::dump(dump_file);
576  }
double Next_dt
Definition: axisym_vibrating_shell_non_newtonian.cc:665

References oomph::Problem_Parameter::Doc_info_soln, oomph::Problem_Parameter::Doc_info_trace, and oomph::DocInfo::number().

◆ enable_aitken_extrapolation_for_all_elements()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::enable_aitken_extrapolation_for_all_elements ( )
inline

Enable use of Aitken extrapolation for all elements.

424  {
425  unsigned n_element = Fluid_mesh_pt->nelement();
426  for(unsigned e=0;e<n_element;e++)
427  {
428  // Upcast from GeneralisedElement to the present element
429  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
430  el_pt->enable_aitken_extrapolation();
431  }
432  }

References e().

◆ enable_fixed_point_iteration_for_strain_rate_for_all_elements()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::enable_fixed_point_iteration_for_strain_rate_for_all_elements ( )
inline

Enable use of fixed point iteration for all elements

400  {
401  unsigned n_element = Fluid_mesh_pt->nelement();
402  for(unsigned e=0;e<n_element;e++)
403  {
404  // Upcast from GeneralisedElement to the present element
405  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
406  el_pt->enable_fixed_point_iteration_for_strain_rate();
407  }
408  }

References e().

◆ global_temporal_error_norm()

template<class ELEMENT >
double AxisymmetricVibratingShellProblem< ELEMENT >::global_temporal_error_norm
virtual

Error norm to determine the next time step.

Calculate the global temporal error norm.

Reimplemented from oomph::Problem.

1332 {
1333 
1334  //oomph_info << "hierher bypassed global temporal error norm\n";
1335  //return 0.0;
1336 
1337  // Initialise
1338  double global_error = 0.0;
1339 
1340  //Find out how many nodes there are in the problem
1341  const unsigned n_node = Fluid_mesh_pt->nnode();
1342 
1343 
1344  oomph_info << "in here with node = " << n_node << std::endl;
1345 
1346  //Loop over the nodes and calculate the errors in the positions
1347  for(unsigned n=0;n<n_node;n++)
1348  {
1349  //Find number of dimensions of the node
1350  const unsigned n_dim = Fluid_mesh_pt->node_pt(n)->ndim();
1351  //Set the position error to zero
1352  double node_position_error = 0.0;
1353  //Loop over the dimensions
1354  for(unsigned i=0;i<n_dim;i++)
1355  {
1356  //Get position error
1357  double error =
1358  Fluid_mesh_pt->node_pt(n)->position_time_stepper_pt()->
1359  temporal_error_in_position(Fluid_mesh_pt->node_pt(n),i);
1360 
1361  //Add the square of the individual error to the position error
1362  node_position_error += error*error;
1363  }
1364 
1365  //Divide the position error by the number of dimensions
1366  node_position_error /= n_dim;
1367 
1368  //Now add to the global error
1369  global_error += node_position_error;
1370  }
1371 
1372  //Now the global error must be divided by the number of nodes
1373  global_error /= n_node;
1374 
1375  oomph_info << "done global error = " << global_error << std::endl;
1376 
1377  //Return the square root of the errr
1378  return sqrt(global_error);
1379 
1380 
1381 }
int error
Definition: calibrate.py:297

References calibrate::error, i, n, oomph::oomph_info, and sqrt().

◆ height_central_node()

template<class ELEMENT >
double AxisymmetricVibratingShellProblem< ELEMENT >::height_central_node ( )
inline

Get the height at the centre node.

381  {
383  }

◆ next_dt()

template<class ELEMENT >
double& AxisymmetricVibratingShellProblem< ELEMENT >::next_dt ( )
inline

Access the next suggested timestep.

377 {return Next_dt;}

◆ restart()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::restart ( )
inline

Restart.

581  {
582  // Pointer to restart file
583  ifstream* restart_file_pt=0;
584 
585  // Open restart file from stem
586  restart_file_pt=new ifstream(Problem_Parameter::Restart_file.c_str(),
587  ios_base::in);
588  if (restart_file_pt!=0)
589  {
590  oomph_info << "Have opened "
592  << " for restart. " << std::endl;
593  }
594  else
595  {
596  std::ostringstream error_stream;
597  error_stream
598  << "ERROR while trying to open "
600  << " for restart." << std::endl;
601 
602  throw OomphLibError(
603  error_stream.str(),
604  "restart()",
606  }
607 
608 
609  // Read restart data:
610  //-------------------
611  if (restart_file_pt!=0)
612  {
613 
614  // Doc number for trace
615  //---------------------
616  // Read line up to termination sign
617  string input_string;
618  getline(*restart_file_pt,input_string,'#');
619 
620  // Ignore rest of line
621  restart_file_pt->ignore(80,'\n');
622 
623  // Doc number
624  Problem_Parameter::Doc_info_trace.number()=unsigned(atoi(input_string.c_str()));
625 
626 
627  // Doc number for solution
628  //---------------------
629  // Read line up to termination sign
630  getline(*restart_file_pt,input_string,'#');
631 
632  // Ignore rest of line
633  restart_file_pt->ignore(80,'\n');
634 
635  // Doc number
636  Problem_Parameter::Doc_info_soln.number()=unsigned(atoi(input_string.c_str()));
637 
638 
639  // Next timestep (required for restart from temporally adaptive run
640  //-----------------------------------------------------------------
641 
642  // Read line up to termination sign
643  getline(*restart_file_pt,input_string,'#');
644 
645  // Ignore rest of line
646  restart_file_pt->ignore(80,'\n');
647 
648  // Next timestep
649  Next_dt=double(atof(input_string.c_str()));
650 
651 
652  // Refine the mesh and read in the generic problem data
653  Problem::read(*restart_file_pt);
654  }
655  }
Definition: oomph_definitions.h:222
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::enable_if_t< PacketLoad, PacketType > read(const TensorMapper &tensorMapper, const StorageIndex &NCIndex, const StorageIndex &CIndex, const StorageIndex &ld)
read, a template function used for loading the data from global memory. This function is used to guar...
Definition: TensorContractionSycl.h:162
std::string Restart_file
Name of restart file.
Definition: axisym_vibrating_shell_non_newtonian.cc:67
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61

References oomph::Problem_Parameter::Doc_info_soln, oomph::Problem_Parameter::Doc_info_trace, oomph::DocInfo::number(), OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Eigen::TensorSycl::internal::read(), and oomph::Problem_Parameter::Restart_file.

◆ set_nprev_for_extrapolation_of_strain_rate_for_all_elements()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::set_nprev_for_extrapolation_of_strain_rate_for_all_elements ( const unsigned &  nprev)
inline
449  {
450  unsigned n_element = Fluid_mesh_pt->nelement();
451  for(unsigned e=0;e<n_element;e++)
452  {
453  // Upcast from GeneralisedElement to the present element
454  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
455 
456  // Use extrapolated strain rate when determining viscosity?
458  ("--use_current_strainrate_for_viscosity"))
459  {
460  el_pt->use_extrapolated_strainrate_to_compute_second_invariant();
461 
462  // Set extrapolation order
463  el_pt->nprev_for_extrapolation_of_strain_rate()=nprev;
464  }
465  else
466  {
467  el_pt->use_current_strainrate_to_compute_second_invariant();
468  }
469  }
470  }

References oomph::CommandLineArgs::command_line_flag_has_been_set(), and e().

◆ strainrate_norm()

template<class ELEMENT >
double AxisymmetricVibratingShellProblem< ELEMENT >::strainrate_norm ( )
inline
529  {
530  // Compute current volume and error in extrapolation from fluid elements
531  double norm_squared=0.0;
532  double extrapolated_norm_squared=0.0;
533  double error_norm_squared=0.0;
534  unsigned nel=Fluid_mesh_pt->nelement();
535  for (unsigned e=0;e<nel;e++)
536  {
537  // Get element
538  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(
539  Fluid_mesh_pt->element_pt(e));
540 
541  // Get norms of invariant
542  double el_norm_squared=0.0;
543  double el_extrapolated_norm_squared=0.0;
544  double el_error_norm_squared=0.0;
545  double test_size=0.0;
546  test_size=el_pt->square_of_norm_of_strain_invariant(
547  el_norm_squared,
548  el_extrapolated_norm_squared,
549  el_error_norm_squared);
550 
551  // Add it...
552  norm_squared+=el_norm_squared;
553  extrapolated_norm_squared+=el_extrapolated_norm_squared;
554  error_norm_squared+=el_error_norm_squared;
555  }
556 
557  return sqrt(norm_squared);
558  }

References e(), and sqrt().

◆ update_latest_fixed_point_iteration_guess_for_strain_rate_for_all_elements()

template<class ELEMENT >
void AxisymmetricVibratingShellProblem< ELEMENT >::update_latest_fixed_point_iteration_guess_for_strain_rate_for_all_elements ( )
inline

Update latest guess for strain rate.

387  {
388  unsigned n_element = Fluid_mesh_pt->nelement();
389  for(unsigned e=0;e<n_element;e++)
390  {
391  // Upcast from GeneralisedElement to the present element
392  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
393  el_pt->update_latest_fixed_point_iteration_guess_for_strain_rate();
394  }
395  }

References e().

Member Data Documentation

◆ Central_node_on_free_surface_pt

template<class ELEMENT >
Node* AxisymmetricVibratingShellProblem< ELEMENT >::Central_node_on_free_surface_pt
private

Pointer to free surface node on symmetry line.

◆ External_pressure_data_pt

template<class ELEMENT >
Data* AxisymmetricVibratingShellProblem< ELEMENT >::External_pressure_data_pt
private

Pointer to a global external pressure datum.

◆ Fluid_mesh_pt

template<class ELEMENT >
RefineableSolidTriangleMesh<ELEMENT>* AxisymmetricVibratingShellProblem< ELEMENT >::Fluid_mesh_pt
private

Pointer to Fluid_mesh.

◆ Free_surface_mesh_pt

template<class ELEMENT >
Mesh* AxisymmetricVibratingShellProblem< ELEMENT >::Free_surface_mesh_pt
private

Pointers to mesh of free surface elements.

◆ Internal_open_boundary_pt

template<class ELEMENT >
Vector<TriangleMeshOpenCurve* > AxisymmetricVibratingShellProblem< ELEMENT >::Internal_open_boundary_pt
private

◆ Next_dt

template<class ELEMENT >
double AxisymmetricVibratingShellProblem< ELEMENT >::Next_dt
private

Suggestion for the next timestep (allows it to be written to or read from a restart file)

◆ Outer_boundary_polyline_pt

template<class ELEMENT >
TriangleMeshClosedCurve* AxisymmetricVibratingShellProblem< ELEMENT >::Outer_boundary_polyline_pt
private

Triangle mesh polygon for outer boundary.

◆ Pext

template<class ELEMENT >
double AxisymmetricVibratingShellProblem< ELEMENT >::Pext
private

External pressure.


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