ThreeDimBethertonProblem< ELEMENT > Class Template Reference

Brethreton problem in square tube domain. More...

+ Inheritance diagram for ThreeDimBethertonProblem< ELEMENT >:

Public Member Functions

 ThreeDimBethertonProblem ()
 Constructor. More...
 
 ~ThreeDimBethertonProblem ()
 Destructor to clean up memory. More...
 
void fix_pressure (const unsigned int &e, const unsigned int &pdof, const double &pvalue)
 Fix pressure in element e at pressure dof pdof and set to pvalue. More...
 
void actions_before_newton_convergence_check ()
 Remesh before any convergence checke. More...
 
double get_outflow ()
 
double get_lambda ()
 
void multiply_aspect_ratio (const double factor, const bool &mult_width)
 
STSpineMesh< Hijacked< ELEMENT >, SpineSurfaceFluidInterfaceElement< ELEMENT > > * mesh_pt ()
 
void doc_solution (DocInfo &doc_info)
 Doc the solution. 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
 
virtual void actions_before_adapt ()
 
virtual void actions_after_adapt ()
 Actions that are to be performed after a mesh adaptation. More...
 
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 ()
 

Public Attributes

double Fixed_spine_height
 
- Public Attributes inherited from oomph::Problem
bool Shut_up_in_newton_solve
 

Private Attributes

Vector< FiniteElement * > Inlet_traction_element_pt
 
Vector< FiniteElement * > Outlet_traction_element_pt
 
Vector< FiniteElement * > Interface_line_element_pt
 
FiniteElementFix_spine_height_element_pt
 
DataBubble_pressure_data_pt
 
ofstream Trace_file
 Trace file. More...
 
SpineTrace_spine_r_h_top
 
SpineTrace_spine_r_h_middle
 
SpineTrace_spine_r_h_bottom
 
SpineTrace_spine_r_d_top
 
SpineTrace_spine_r_d_bottom
 

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...
 
- 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_solve ()
 
virtual void actions_after_newton_solve ()
 
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 ()
 
virtual double global_temporal_error_norm ()
 
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 ThreeDimBethertonProblem< ELEMENT >

Brethreton problem in square tube domain.

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

Constructor & Destructor Documentation

◆ ThreeDimBethertonProblem()

template<class ELEMENT >
ThreeDimBethertonProblem< ELEMENT >::ThreeDimBethertonProblem

Constructor.

Constructor for RectangularDrivenCavity problem.

327 {
328 
331  double length_tip = Global_Physical_Variables::Length_tip;
332  double length_liq = Global_Physical_Variables::Length_liq;
333  double length_can = Global_Physical_Variables::Length_can;
335  double rat_press = Global_Physical_Variables::Rat_press;
336  this->Max_residuals = 10.0;
337  this->Max_newton_iterations = 15;
338 
339  unsigned ncan = Global_Physical_Variables::Ncan;
340  unsigned ntip = Global_Physical_Variables::Ntip;
341  unsigned nliq = Global_Physical_Variables::Nliq;
342 
343 // Build and assign mesh
344  Problem::mesh_pt() = new
346  (ncan,ntip,nliq, alpha,length_can,length_tip,length_liq,height,radius);
347 
348  cout<<"Nnodes boundary 5: "<<mesh_pt()->nboundary_node(5)<<endl;
349  cout<<"Nnodes spines : "<<mesh_pt()->nspine()<<endl;
350 
351  // Complete the build of all elements so they are fully functional
352 
353  // Set the fixed spine element
354  unsigned num_nodes = mesh_pt()->nnode();
355  SpineNode* spinenode_fixed=0;
356 
357  //Look fot the spine
358  for(unsigned inod = 0;inod< num_nodes;inod++)
359  {
360  SpineNode* spnd_pt = dynamic_cast<SpineNode*>(mesh_pt()->node_pt(inod));
362  double length_tip = Global_Physical_Variables::Length_tip;
363  double distance =
364  sqrt( (spnd_pt->x(0))* (spnd_pt->x(0)) +
365  (spnd_pt->x(1) + (length_tip * radius))*
366  (spnd_pt->x(1) + (length_tip*radius)) +
367  (spnd_pt->x(2) - height/2)* (spnd_pt->x(2) - height/2) );
368 
369  if(distance < 10E-8)
370  {
371  cout<<"Spine to be pinned found"<<endl;
372  spinenode_fixed = spnd_pt;
373  }
374  }
375 
376  //Create a fixed element using the central spine
377  FixSpineHeightElement* fix_spine_element_pt =
378  new FixSpineHeightElement(spinenode_fixed);
379 
380  //Set the fixed spine height
381  Fixed_spine_height = spinenode_fixed->spine_pt()->height();
382 
383  //Add the fixed height to it
384  fix_spine_element_pt->height_pt() = &Fixed_spine_height;
385 
386  //Store the element in the problem
387  Fix_spine_height_element_pt = fix_spine_element_pt;
388 
389  //Create a bubble presure data
390  Bubble_pressure_data_pt = new Data(1);
391 
392  //Set values
393  Bubble_pressure_data_pt->set_value(0,rat_press);
394 
395  //This will be global data
397 
398  //Set the pressure data
399  fix_spine_element_pt->set_traded_pressure_data(Bubble_pressure_data_pt);
400 
401  //Add the Fixed element to the mesh
402  mesh_pt()->add_element_pt(fix_spine_element_pt);
403 
404  cout<<"Bubble pressure "<<Bubble_pressure_data_pt->value(0)<<endl;
405 
406  //Find number of bulk elements in mesh
407  unsigned n_bulk = mesh_pt()->nbulk();
408 
409  // Loop over the elements to set up element-specific
410  // things that cannot be handled by constructor
411  for(unsigned e=0;e<n_bulk;e++)
412  {
413  // Upcast from GeneralisedElement to the present element
414  Hijacked<ELEMENT>* el_pt =
415  dynamic_cast<Hijacked<ELEMENT>* >(mesh_pt()->bulk_element_pt(e));
416 
417  //Set the Reynolds number
418  el_pt->re_pt() = &Global_Physical_Variables::Re;
419 
420  //Set the ReInvFr number
421  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
422 
423  //Set the gravity
424  el_pt->g_pt() = &Global_Physical_Variables::G;
425 
426  //Insist that shape derivatives are calculated by "full" finite
427  //differences (which is the silent default)
428  //el_pt->evaluate_shape_derivs_by_direct_fd();
429  }
430 
431 
432  //Find number of interface elements in mesh
433  unsigned n_interface = mesh_pt()->ninterface_element();
434  unsigned nhi = 0;
435 
436  for(unsigned e=0;e<n_interface;e++)
437  {
438  // Upcast from GeneralisedElement to the present element
441  (mesh_pt()->interface_element_pt(e));
443 
444  // Set the external pressure data
446 
447  // Quite inefficient way of Hijack the nodes
448  unsigned nnodes = el_pt->nnode();
449  for(unsigned i =0; i<nnodes;i++)
450  {
451  if( ( (el_pt->node_pt(i)->x(1)- length_can)*
452  (el_pt->node_pt(i)->x(1)- length_can) )< 10E-10 )
453  {
454  el_pt->hijack_nodal_value(i,1);
455  nhi++;
456  }
457  }
458 
459  //Insist that shape derivatives are calculated by "full" finite
460  //differences (which is the silent default)
461  //el_pt->evaluate_shape_derivs_by_direct_fd();
462  }
463 
464  cout<<nhi<<" hijacked nodes at the surface elements."<<endl;
465 
466 
467  //Create and set the inlet elements
468  unsigned long ninlet = mesh_pt()->nbulkinlet();
469  for(unsigned long i = 0;i<ninlet;i++)
470  {
472  inlet_element_pt =
474  mesh_pt()->bulk_inlet_element_pt(i),
475  mesh_pt()->face_index_inlet());
476 
477  //Add Elements to the local storage
478  Inlet_traction_element_pt.push_back(inlet_element_pt);
479 
480  //Set the traction function
481  inlet_element_pt->traction_fct_pt() =
483 
484  //Insist that shape derivatives are calculated by "full" finite
485  //differences (which is the silent default)
486  //inlet_element_pt->evaluate_shape_derivs_by_direct_fd();
487 
488 
489  //Add element to the mesh
490  mesh_pt()->add_element_pt(inlet_element_pt);
491  }
492 
493  //Create and set the outlet elements
494  unsigned long noutlet = mesh_pt()->nbulkoutlet();
495  for(unsigned long i = 0;i<noutlet;i++)
496  {
497  SpineGravityTractionElement<Hijacked<ELEMENT> >* outlet_element_pt =
499  mesh_pt()->bulk_outlet_element_pt(i),
500  mesh_pt()->face_index_outlet());
501 
502  //Set the ReInvFr number
503  outlet_element_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
504 
505  //Set the gravity
506  outlet_element_pt->g_pt() = &Global_Physical_Variables::G;
507 
508  //Add Elements to the local storage
509  Outlet_traction_element_pt.push_back(outlet_element_pt);
510 
511  //Insist that shape derivatives are calculated by "full" finite
512  //differences (which is the silent default)
513  //outlet_element_pt->evaluate_shape_derivs_by_direct_fd();
514 
515  //Add element to the mesh
516  mesh_pt()->add_element_pt(outlet_element_pt);
517  }
518 
519 
520  //Create and set the line elements elements
521  //(the coordinates s constants are the same as in the outlet elements
522  unsigned long nlineel = mesh_pt()->ninterfaceline();
523  for(unsigned long i = 0;i<nlineel;i++)
524  {
525  //Cast the return type to the edge element and also
526  //cast the interface element to the specific interface element
530  ( mesh_pt()->interface_line_element_pt(i))->make_bounding_element(
531  mesh_pt()->face_index_outlet()));
532 
533  //Set the capillary number
534  line_element_pt->ca_pt() = &Global_Physical_Variables::Ca;
535 
536  //Add Elements to the local storagea
537  Interface_line_element_pt.push_back(line_element_pt);
538 
539  //Add element to the mesh
540  mesh_pt()->add_element_pt(line_element_pt);
541 
542  //Hijack the nodes
543  for(unsigned e=0;e<line_element_pt->nnode();e++)
544  line_element_pt ->hijack_nodal_value(e,1);
545  }
546 
547  unsigned n_element = mesh_pt()->nelement();
548  cout<<"We have "<<n_bulk<<" bulk elemens, "<<n_interface<<" interface elements, "<<n_element<<" generic elements in this mesh"<<endl;
549 
550  mesh_pt()->node_update();
551 
552  //Now set the boundary conditions by pinning those variables for
553  //which we apply Dirichlet boundary conditions
554 
555 // Just pin the boundary Dirichlet conditions
556  unsigned ibound = 0;
557  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
558 
559  //Base of the tube (pin all values)
560  for (unsigned inod=0;inod<num_nod;inod++)
561  {
562  // Loop over values (ux uy and uz velocities)
563  for (unsigned i=0;i<3;i++)
564  {
565  mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
566  }
567  }
568 
569  //Inlet boundary (pin the transverse velocity components)
570  ibound = 1;
571  num_nod= mesh_pt()->nboundary_node(ibound);
572  for (unsigned inod=0;inod<num_nod;inod++)
573  {
574  // Loop over values (ux uy and uz velocities)
575  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
576  mesh_pt()->boundary_node_pt(ibound,inod)->pin(2);
577  }
578 
579  //Outer wall of the tube (pin all velocitites)
580  ibound = 2;
581  num_nod= mesh_pt()->nboundary_node(ibound);
582  for (unsigned inod=0;inod<num_nod;inod++)
583  {
584  // Loop over values (ux uy and uz velocities)
585  for (unsigned i=0;i<3;i++)
586  {
587  mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
588  }
589  }
590 
591  //Outlet (leave free)
592  //ibound = 3;
593  //num_nod= mesh_pt()->nboundary_node(ibound);
594  //for (unsigned inod=0;inod<num_nod;inod++)
595  // {
596  // }
597 
598  //Symmetry boundary (apply symmetry conditions in x)
599  ibound = 4;
600  num_nod= mesh_pt()->nboundary_node(ibound);
601  for (unsigned inod=0;inod<num_nod;inod++)
602  {
603  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
604  }
605 
606  //Free surface (leave free
607  //ibound = 5;
608  //num_nod= mesh_pt()->nboundary_node(ibound);
609  //for (unsigned inod=0;inod<num_nod;inod++)
610  // {
611  // }
612 
613  //Top of the tube (pin all velocities)
614  ibound = 6;
615  num_nod= mesh_pt()->nboundary_node(ibound);
616  for (unsigned inod=0;inod<num_nod;inod++)
617  {
618  // Loop over values (ux uy and uz velocities)
619  for (unsigned i=0;i<3;i++)
620  {
621  mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
622  }
623  }
624 
625  // Setup equation numbering scheme
626  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
627 
628 //Set initial guess that all velocities are uniform in the y direction
629  const unsigned n_node = mesh_pt()->nnode();
630  for(unsigned n=0;n<n_node;n++)
631  {
632  mesh_pt()->node_pt(n)->set_value(0,0.0);
633  mesh_pt()->node_pt(n)->set_value(1,1.0);
634  mesh_pt()->node_pt(n)->set_value(2,0.0);
635  }
636 
637  //Set boundary conditions
638  //Set the no-slip conditions on all tube boundaries
639  unsigned num_bound = mesh_pt()->nboundary();
640  for(unsigned ibound=0;ibound<num_bound;ibound++)
641  {
642  //If we are not on the symmetry boundary (4)
643  //Outlet (3), Inlet(1) or free surface (5), set the no-slip conditions
644  if((ibound != 4) && (ibound != 3) && (ibound != 1) && (ibound !=5) )
645  {
646  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
647  for (unsigned inod=0;inod<num_nod;inod++)
648  {
649  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,0.0);
650  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,1.0);
651  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(2,0.0);
652  }
653  }
654  }
655 
656  //Prepare output information
657 
658 //Open trace file
659  char filename[100];
660  sprintf(filename,"RESLT/trace.dat");
661  Trace_file.open(filename);
662  Trace_file << "VARIABLES=";
663  Trace_file << "\"Ca\",\t \"Bo\", \t \"Alpha\", \t \"m\", \t \"p<sub>bubble</sub>\", \t \"r<sub>h_top</sub>\", \t \"r<sub>d_top</sub>\", \t \"r<sub>h_middle</sub>\", \t \"r<sub>d_bottom</sub>\", \t \"r<sub>h_bottom</sub>\",\t \"lambda\" ";
664 
665  Trace_file<<endl;
666 
667 //Find the trace spines
668 
669 // Coordinates of the node which indicate the spine to trace
670  double x_trace;
671  double y_trace;
672  double z_trace;
673 
674  // All the nodes to trace are in the last slice
675  y_trace = length_can;
676 
677 //Find the spine of the top (horizontal spine)
678  x_trace = 0.0;
679  z_trace = height;
680 
681 
682  //Loop over all nodes
683  for(unsigned inod = 0;inod<n_node;inod++)
684  {
685  SpineNode* spnd_pt = dynamic_cast<SpineNode*>(mesh_pt()->node_pt(inod));
686 
687  double distance =sqrt(
688  (spnd_pt->x(0) - x_trace)* (spnd_pt->x(0) - x_trace) +
689  (spnd_pt->x(1) - y_trace)* (spnd_pt->x(1) - y_trace) +
690  (spnd_pt->x(2) - z_trace)* (spnd_pt->x(2) - z_trace) );
691  if(distance < 10E-8)
692  {
693  cout<<"Top horizontal Spine found"<<endl;
694  Trace_spine_r_h_top = spnd_pt->spine_pt();
695  }
696  }
697 
698 
699 //Find the spine of the top (diagonal spine)
700  x_trace = alpha/2;
701  z_trace = height;
702 
703  for(unsigned inod = 0;inod<n_node;inod++)
704  {
705  SpineNode* spnd_pt = dynamic_cast<SpineNode*>(mesh_pt()->node_pt(inod));
706 
707  double distance =sqrt(
708  (spnd_pt->x(0) - x_trace)* (spnd_pt->x(0) - x_trace) +
709  (spnd_pt->x(1) - y_trace)* (spnd_pt->x(1) - y_trace) +
710  (spnd_pt->x(2) - z_trace)* (spnd_pt->x(2) - z_trace) );
711  if(distance < 10E-8)
712  {
713  cout<<"Top diagonal Spine found"<<endl;
714  Trace_spine_r_d_top = spnd_pt->spine_pt();
715  }
716  }
717 
718 
719 //Find the spine of the middle (horizontal spine)
720  x_trace = alpha/2;
721  z_trace = height/2;
722 
723  for(unsigned inod = 0;inod<n_node;inod++)
724  {
725  SpineNode* spnd_pt = dynamic_cast<SpineNode*>(mesh_pt()->node_pt(inod));
726 
727  double distance =sqrt(
728  (spnd_pt->x(0) - x_trace)* (spnd_pt->x(0) - x_trace) +
729  (spnd_pt->x(1) - y_trace)* (spnd_pt->x(1) - y_trace) +
730  (spnd_pt->x(2) - z_trace)* (spnd_pt->x(2) - z_trace) );
731  if(distance < 10E-8)
732  {
733  cout<<"Middle horizontal Spine found"<<endl;
734  Trace_spine_r_h_middle = spnd_pt->spine_pt();
735  }
736  }
737 
738 
739 //Find the spine of the bottom (horizontal spine)
740  x_trace = 0.0;
741  z_trace = 0.0;
742 
743 
744 
745  for(unsigned inod = 0;inod<n_node;inod++)
746  {
747  SpineNode* spnd_pt = dynamic_cast<SpineNode*>(mesh_pt()->node_pt(inod));
748 
749  double distance =sqrt(
750  (spnd_pt->x(0) - x_trace)* (spnd_pt->x(0) - x_trace) +
751  (spnd_pt->x(1) - y_trace)* (spnd_pt->x(1) - y_trace) +
752  (spnd_pt->x(2) - z_trace)* (spnd_pt->x(2) - z_trace) );
753  if(distance < 10E-8)
754  {
755  cout<<"Bottom horizontal Spine found"<<endl;
756  Trace_spine_r_h_bottom = spnd_pt->spine_pt();
757  }
758  }
759 
760 
761 //Find the spine of the bottom (diagonal spine)
762  x_trace = alpha/2;
763  z_trace = 0.0;
764 
765 
766 
767  for(unsigned inod = 0;inod<n_node;inod++)
768  {
769  SpineNode* spnd_pt = dynamic_cast<SpineNode*>(mesh_pt()->node_pt(inod));
770 
771  double distance =sqrt(
772  (spnd_pt->x(0) - x_trace)* (spnd_pt->x(0) - x_trace) +
773  (spnd_pt->x(1) - y_trace)* (spnd_pt->x(1) - y_trace) +
774  (spnd_pt->x(2) - z_trace)* (spnd_pt->x(2) - z_trace) );
775  if(distance < 10E-8)
776  {
777  cout<<"Bottom diagonal Spine found"<<endl;
778  Trace_spine_r_d_bottom = spnd_pt->spine_pt();
779  }
780 
781  }
782 }
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
int i
Definition: BiCGSTAB_step_by_step.cpp:9
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Point element that is used to set the bubble pressure.
Definition: elastic_bretherton.cc:994
double *& height_pt()
Access function to the prescribed spine height.
Definition: elastic_bretherton.cc:1040
void set_traded_pressure_data(Data *traded_pressure_data_pt)
Definition: elastic_bretherton.cc:1065
Definition: st_mesh.h:54
Definition: elastic_bretherton.cc:344
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
Definition: elastic_bretherton.cc:522
double *& re_invfr_pt()
Pointer to Reynolds number divided by Froude number.
Definition: elastic_bretherton.cc:513
Spine * Trace_spine_r_d_top
Definition: three_d_breth.cc:315
Data * Bubble_pressure_data_pt
Definition: three_d_breth.cc:295
Spine * Trace_spine_r_h_top
Definition: three_d_breth.cc:309
Spine * Trace_spine_r_d_bottom
Definition: three_d_breth.cc:317
Spine * Trace_spine_r_h_bottom
Definition: three_d_breth.cc:313
ofstream Trace_file
Trace file.
Definition: three_d_breth.cc:305
Spine * Trace_spine_r_h_middle
Definition: three_d_breth.cc:311
Vector< FiniteElement * > Inlet_traction_element_pt
Definition: three_d_breth.cc:283
double Fixed_spine_height
Definition: three_d_breth.cc:300
FiniteElement * Fix_spine_height_element_pt
Definition: three_d_breth.cc:292
STSpineMesh< Hijacked< ELEMENT >, SpineSurfaceFluidInterfaceElement< ELEMENT > > * mesh_pt()
Definition: three_d_breth.cc:268
Vector< FiniteElement * > Interface_line_element_pt
Definition: three_d_breth.cc:289
Vector< FiniteElement * > Outlet_traction_element_pt
Definition: three_d_breth.cc:286
Definition: nodes.h:86
void set_value(const unsigned &i, const double &value_)
Definition: nodes.h:271
double value(const unsigned &i) const
Definition: nodes.h:293
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
double *& ca_pt()
Access function to the pointer specifying the capillary number.
Definition: interface_elements.h:172
double *& ca_pt()
Pointer to the Capillary number.
Definition: interface_elements.h:492
void set_external_pressure_data(Data *external_pressure_data_pt)
Definition: interface_elements.h:539
Definition: hijacked_elements.h:132
Data * hijack_nodal_value(const unsigned &n, const unsigned &i, const bool &return_data=true)
Definition: hijacked_elements.h:214
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
unsigned Max_newton_iterations
Maximum number of Newton iterations.
Definition: problem.h:599
void add_global_data(Data *const &global_data_pt)
Definition: problem.h:1654
double Max_residuals
Definition: problem.h:610
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Definition: problem.cc:1989
Definition: spines.h:477
Spine version of the LineFluidInterfaceBoundingElement.
Definition: specific_node_update_interface_elements.h:489
Definition: spines.h:328
Spine *& spine_pt()
Access function to spine.
Definition: spines.h:372
Definition: specific_node_update_interface_elements.h:628
double & height()
Access function to spine height.
Definition: spines.h:150
RealScalar alpha
Definition: level1_cplx_impl.h:151
double E
Elastic modulus.
Definition: TwenteMeshGluing.cpp:68
double Rat_press
Definition: three_d_breth.cc:89
unsigned Ntip
Definition: three_d_breth.cc:83
unsigned Ncan
Definition: three_d_breth.cc:81
double Radius
Radius of the pipe.
Definition: pipe.cc:55
double Length_liq
Definition: three_d_breth.cc:77
double Ca
Capillary number.
Definition: fibre.cc:61
double Length_can
Definition: three_d_breth.cc:79
double Alpha
Scaling factor for wall thickness (to be used in an exercise)
Definition: unsteady_ring.cc:73
double ReInvFr
Product of Reynolds number and inverse of Froude number.
Definition: fibre.cc:58
double height(const double &x)
Height of domain.
Definition: simple_spine_channel.cc:429
double Re
Reynolds number.
Definition: fibre.cc:55
double Length_tip
Definition: three_d_breth.cc:75
double Height
Definition: three_d_breth.cc:73
Vector< double > G(3)
Direction of gravity.
Definition: spherical_shell_convection.cc:62
unsigned Nliq
Definition: three_d_breth.cc:84
void hydrostatic_pressure(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Function that prescribes the hydrostatic pressure field at the outlet.
Definition: elastic_bretherton.cc:152
string filename
Definition: MergeRestartFiles.py:39
radius
Definition: UniformPSDSelfTest.py:15

References alpha, Global_Physical_Variables::Alpha, Global_Physical_Variables::Ca, oomph::FluidInterfaceBoundingElement::ca_pt(), oomph::FluidInterfaceElement::ca_pt(), Global_Physical_Variables::E, e(), MergeRestartFiles::filename, Global_Physical_Variables::G, SpineGravityTractionElement< ELEMENT >::g_pt(), oomph::Spine::height(), Global_Physical_Variables::height(), Global_Physical_Variables::Height, FixSpineHeightElement::height_pt(), oomph::Hijacked< ELEMENT >::hijack_nodal_value(), Global_Physical_Variables::hydrostatic_pressure(), i, Global_Physical_Variables::Length_can, Global_Physical_Variables::Length_liq, Global_Physical_Variables::Length_tip, oomph::Locate_zeta_helpers::Max_newton_iterations, n, Global_Physical_Variables::Ncan, Global_Physical_Variables::Nliq, oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), Global_Physical_Variables::Ntip, UniformPSDSelfTest::radius, Global_Physical_Variables::Radius, Global_Physical_Variables::Rat_press, Global_Physical_Variables::Re, SpineGravityTractionElement< ELEMENT >::re_invfr_pt(), Global_Physical_Variables::ReInvFr, oomph::FluidInterfaceElement::set_external_pressure_data(), FixSpineHeightElement::set_traded_pressure_data(), oomph::SpineNode::spine_pt(), sqrt(), oomph::Problem_Parameter::Trace_file, and oomph::Node::x().

◆ ~ThreeDimBethertonProblem()

template<class ELEMENT >
ThreeDimBethertonProblem< ELEMENT >::~ThreeDimBethertonProblem

Destructor to clean up memory.

Destructor for RectangularDrivenCavity problem.

790 {
791 
792  // Timestepper gets killed in general problem destructor
793 
794  // Mesh gets killed in general problem destructor
795 
796 }

Member Function Documentation

◆ actions_before_newton_convergence_check()

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

Remesh before any convergence checke.

Reimplemented from oomph::Problem.

136  {
137  mesh_pt()->node_update();
138 
139  // This driver code cannot be allowed to use the analytical form of
140  // get_dresidual_dnodal_coordinates(...) that is implemented in the
141  // NavierStokesEquations class, since the elemental residuals have
142  // contributions from external data which is not taken into account
143  // by that routine. We therefore force the bulk elements to use the
144  // fully-finite differenced version.
145  const unsigned n_element = mesh_pt()->nelement();
146  for(unsigned e=0;e<n_element;e++)
147  {
148  if(dynamic_cast<ElementWithMovingNodes*>(mesh_pt()->element_pt(e)))
149  {
150  ElementWithMovingNodes* el_pt =
151  dynamic_cast<ElementWithMovingNodes*>(mesh_pt()->element_pt(e));
153  }
154  }
155  }
Definition: element_with_moving_nodes.h:65
void evaluate_shape_derivs_by_direct_fd()
Evaluate shape derivatives by direct finite differencing.
Definition: element_with_moving_nodes.h:210

References e(), and oomph::ElementWithMovingNodes::evaluate_shape_derivs_by_direct_fd().

◆ doc_solution()

template<class ELEMENT >
void ThreeDimBethertonProblem< ELEMENT >::doc_solution ( DocInfo doc_info)

Doc the solution.

806 {
807  ofstream some_file;
808  char filename[100];
809 
811  double width = Global_Physical_Variables::Alpha/2;
812  double diagonal = sqrt(height*height + width*width ) ;
813  double scale = 1/height; //scale for the radius
814 
815  double outflow = this->get_outflow();
816 
817  double capillary_number = Global_Physical_Variables::Ca;
818 
819  double wet_fraction = outflow/width/(height*2);
820 //The inflow is calculated in half channel
821 
822  double p_drop = Bubble_pressure_data_pt->value(0);
823 
824 // Get finger width at end of the computational domain
825  double lambda = get_lambda();
826 
827  Trace_file<<capillary_number;
828 
830 
832 
833  Trace_file<<" \t"<<wet_fraction;
834 
835  Trace_file<<" \t"<<p_drop;;
836 
837  Trace_file<<" \t"<<(height- Trace_spine_r_h_top->height() )*scale ;
838  Trace_file<<" \t"<<(diagonal- Trace_spine_r_d_top->height() )*scale ;
839  Trace_file<<" \t"<<(width - Trace_spine_r_h_middle->height() )*scale ;
840  Trace_file<<" \t"<<(diagonal - Trace_spine_r_d_bottom->height() )*scale ;
841  Trace_file<<" \t"<<(height - Trace_spine_r_h_bottom->height() )*scale ;
842 
843  Trace_file<<" \t"<<lambda ;
844 
845  Trace_file<<endl;
846 
847  // Number of plot points
848  unsigned npts;
849  npts=5;
850 
851  mesh_pt()->node_update();
852 
853  // Output solution
854  sprintf(filename,"%s/soln%i.dat",
855  doc_info.directory().c_str(),doc_info.number());
856  some_file.open(filename);
857  unsigned n_bulk = mesh_pt()->nbulk();
858  for(unsigned e=0;e<n_bulk;e++)
859  {
860  Hijacked<ELEMENT>* el_pt =
861  dynamic_cast<Hijacked<ELEMENT>* >(mesh_pt()->bulk_element_pt(e));
862  el_pt->output(some_file,npts);
863  }
864  some_file.close();
865 }
cout<< "The eigenvalues of A are:"<< endl<< ces.eigenvalues()<< endl;cout<< "The matrix of eigenvectors, V, is:"<< endl<< ces.eigenvectors()<< endl<< endl;complex< float > lambda
Definition: ComplexEigenSolver_compute.cpp:9
double get_lambda()
Definition: three_d_breth.cc:172
double get_outflow()
Definition: three_d_breth.cc:159
std::string directory() const
Output directory.
Definition: oomph_utilities.h:524
unsigned & number()
Number used (e.g.) for labeling output files.
Definition: oomph_utilities.h:554
void diagonal(const MatrixType &m)
Definition: diagonal.cpp:13
double Bo
Bond number.
Definition: elastic_bretherton.cc:97

References Global_Physical_Variables::Alpha, Global_Physical_Variables::Bo, Global_Physical_Variables::Ca, diagonal(), oomph::DocInfo::directory(), e(), MergeRestartFiles::filename, Global_Physical_Variables::height(), Global_Physical_Variables::Height, lambda, oomph::DocInfo::number(), sqrt(), and oomph::Problem_Parameter::Trace_file.

◆ fix_pressure()

template<class ELEMENT >
void ThreeDimBethertonProblem< ELEMENT >::fix_pressure ( const unsigned int e,
const unsigned int pdof,
const double pvalue 
)
inline

Fix pressure in element e at pressure dof pdof and set to pvalue.

128  {
129  //Cast to full element type and fix the pressure at that element
130  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
131  fix_pressure(pdof,pvalue);
132  }
void fix_pressure(const unsigned int &e, const unsigned int &pdof, const double &pvalue)
Fix pressure in element e at pressure dof pdof and set to pvalue.
Definition: three_d_breth.cc:126

References e().

◆ get_lambda()

template<class ELEMENT >
double ThreeDimBethertonProblem< ELEMENT >::get_lambda ( )
inline
173  {
174  //Initial max
175  double max_x = 0.0;
176  //Loop over all in interface elements
177  unsigned n_line_interface = Interface_line_element_pt.size();
178  for(unsigned e=0;e<n_line_interface;e++)
179  {
180  //Cache the element
181  FiniteElement* const element_pt = Interface_line_element_pt[e];
182  //Loop over all nodes
183  const unsigned n_node = element_pt->nnode();
184  for(unsigned n=0;n<n_node;n++)
185  {
186  if(element_pt->node_pt(n)->x(0) > max_x)
187  {max_x = element_pt->node_pt(n)->x(0);}
188  }
189  }
190  //Now scale by the half-width
191  return max_x/(Global_Physical_Variables::Alpha/2);
192  }
Definition: elements.h:1313

References Global_Physical_Variables::Alpha, e(), n, oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), and oomph::Node::x().

◆ get_outflow()

template<class ELEMENT >
double ThreeDimBethertonProblem< ELEMENT >::get_outflow ( )
inline
160  {
161  double flow = 0;
162  unsigned int nel = Outlet_traction_element_pt.size();
163  for(unsigned int i =0;i<nel;i++)
164  {
165  flow += dynamic_cast<SpineGravityTractionElement<Hijacked<ELEMENT> >* >
166  (Outlet_traction_element_pt[i])->flow();
167  }
168  return flow;
169  }

References i.

◆ mesh_pt()

template<class ELEMENT >
STSpineMesh<Hijacked<ELEMENT>, SpineSurfaceFluidInterfaceElement<ELEMENT> >* ThreeDimBethertonProblem< ELEMENT >::mesh_pt ( )
inline
269  {
270  // Upcast from pointer to the Mesh base class to the specific
271  // element type that we're using here.
272  return dynamic_cast< STSpineMesh<Hijacked<ELEMENT>,
273  SpineSurfaceFluidInterfaceElement<ELEMENT> >* >(Problem::mesh_pt());
274  }

◆ multiply_aspect_ratio()

template<class ELEMENT >
void ThreeDimBethertonProblem< ELEMENT >::multiply_aspect_ratio ( const double  factor,
const bool mult_width 
)
inline
197  {
198  // until now we have just worked with a square channel
199  // lets change the aspect ratio
200  // we will change also the y aspect ratio
201  // so that the tip is also streched in the y direction can be changed
202 
203  //Loop over all nodes in the mesh
204  const unsigned n_node = mesh_pt()->nnode();
205  for(unsigned n=0;n<n_node;n++)
206  {
207  //Multiply the x-coordinate and y-coordinate by the specified factor
208  Node* const node_pt = mesh_pt()->node_pt(n);
209  if(mult_width){node_pt->x(0) = node_pt->x(0) * factor ;}
210  node_pt->x(1) = node_pt->x(1) * factor;
211  }
212 
213  // Update also the origin of the spines to the new position
214  // of the sidewalls
215  const unsigned n_spine = mesh_pt()->nspine();
216  for(unsigned s = 0;s<n_spine;s++)
217  {
218  Spine* spine_pt = mesh_pt()->spine_pt(s);
219  const double x_spine = spine_pt->geom_parameter(0);
220  const double y_spine = spine_pt->geom_parameter(1);
221  if(mult_width){spine_pt->geom_parameter(0) = x_spine * factor;}
222  spine_pt->geom_parameter(1) = y_spine * factor;
223  }
224 
225 
226  // Update the spines heights according to the new geometry
227  // This is done by looping over all nodes on the free surface
228  // and calculating the distance from the new side wall
229  const unsigned n_boundary_node = mesh_pt()->nboundary_node(5);
230 
231  for(unsigned n=0;n<n_boundary_node;n++)
232  {
233  //Get the root of the spine
234  SpineNode* spine_node_pt =
235  dynamic_cast<SpineNode*>(mesh_pt()->boundary_node_pt(5,n));
236  Spine* spine_pt = spine_node_pt->spine_pt();
237  double x_spine = spine_pt->geom_parameter(0);
238  double y_spine = spine_pt->geom_parameter(1);
239  double z_spine = spine_pt->geom_parameter(2);
240 
241  //Find the distance of the node on the free surface from the root
242  //of the spine
243  double distance = sqrt(
244  (spine_node_pt->x(0) - x_spine )*(spine_node_pt->x(0) - x_spine )+
245  (spine_node_pt->x(1) - y_spine)*(spine_node_pt->x(1) - y_spine)+
246  (spine_node_pt->x(2) - z_spine)*(spine_node_pt->x(2) - z_spine) );
247  //Set the spine length to be the new distance
248  spine_node_pt->h() = distance;
249  }
250 
251  //Change also the distance of the fixed spine by the factor
253  // At the end we change the value of Alpha
254  if(mult_width){Global_Physical_Variables::Alpha =
256  if(!mult_width){ Global_Physical_Variables::Radius =
264  }
Definition: nodes.h:906
double & h()
Access function to spine height.
Definition: spines.h:397
Definition: spines.h:64
double & geom_parameter(const unsigned &i)
Definition: spines.h:287
RealScalar s
Definition: level1_cplx_impl.h:130

References Global_Physical_Variables::Alpha, oomph::Spine::geom_parameter(), oomph::SpineNode::h(), Global_Physical_Variables::Length_can, Global_Physical_Variables::Length_liq, Global_Physical_Variables::Length_tip, n, Global_Physical_Variables::Radius, s, oomph::SpineNode::spine_pt(), sqrt(), and oomph::Node::x().

Member Data Documentation

◆ Bubble_pressure_data_pt

template<class ELEMENT >
Data* ThreeDimBethertonProblem< ELEMENT >::Bubble_pressure_data_pt
private

◆ Fix_spine_height_element_pt

template<class ELEMENT >
FiniteElement* ThreeDimBethertonProblem< ELEMENT >::Fix_spine_height_element_pt
private

◆ Fixed_spine_height

template<class ELEMENT >
double ThreeDimBethertonProblem< ELEMENT >::Fixed_spine_height

◆ Inlet_traction_element_pt

template<class ELEMENT >
Vector<FiniteElement*> ThreeDimBethertonProblem< ELEMENT >::Inlet_traction_element_pt
private

◆ Interface_line_element_pt

template<class ELEMENT >
Vector<FiniteElement*> ThreeDimBethertonProblem< ELEMENT >::Interface_line_element_pt
private

◆ Outlet_traction_element_pt

template<class ELEMENT >
Vector<FiniteElement*> ThreeDimBethertonProblem< ELEMENT >::Outlet_traction_element_pt
private

◆ Trace_file

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

Trace file.

◆ Trace_spine_r_d_bottom

template<class ELEMENT >
Spine* ThreeDimBethertonProblem< ELEMENT >::Trace_spine_r_d_bottom
private

◆ Trace_spine_r_d_top

template<class ELEMENT >
Spine* ThreeDimBethertonProblem< ELEMENT >::Trace_spine_r_d_top
private

◆ Trace_spine_r_h_bottom

template<class ELEMENT >
Spine* ThreeDimBethertonProblem< ELEMENT >::Trace_spine_r_h_bottom
private

◆ Trace_spine_r_h_middle

template<class ELEMENT >
Spine* ThreeDimBethertonProblem< ELEMENT >::Trace_spine_r_h_middle
private

◆ Trace_spine_r_h_top

template<class ELEMENT >
Spine* ThreeDimBethertonProblem< ELEMENT >::Trace_spine_r_h_top
private

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