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...
 
doublenext_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
 
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 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

NodeCentral_node_on_free_surface_pt
 Pointer to free surface node on symmetry line. More...
 
double Next_dt
 
MeshFree_surface_mesh_pt
 Pointers to mesh of free surface elements. More...
 
RefineableSolidTriangleMesh< ELEMENT > * Fluid_mesh_pt
 Pointer to Fluid_mesh. More...
 
DataExternal_pressure_data_pt
 Pointer to a global external pressure datum. More...
 
TriangleMeshClosedCurveOuter_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)
 
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 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: