TwoLayerInterfaceProblem< ELEMENT > Class Template Reference

Problem class to simulate viscous drop propagating along 2D channel. More...

+ Inheritance diagram for TwoLayerInterfaceProblem< ELEMENT >:

Public Member Functions

 TwoLayerInterfaceProblem ()
 Constructor. More...
 
 ~TwoLayerInterfaceProblem ()
 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 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...
 
void set_initial_condition ()
 Set the initial conditions. More...
 
void deform_interface (const double &epsilon, const unsigned &n_periods)
 Function to deform the interface. More...
 
void fix_pressure (const unsigned &e, const unsigned &pdof, const double &pvalue)
 Fix pressure in element e at pressure dof pdof and set to pvalue. 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  {
  Inflow_boundary_id =0 , Upper_wall_boundary_id =1 , Outflow_boundary_id =2 , Bottom_wall_boundary_id =3 ,
  Interface_boundary_id =4
}
 Enumeration of channel 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

MeshFree_surface_mesh_pt
 Pointers to mesh of free surface elements. More...
 
RefineableSolidTriangleMesh< ELEMENT > * Fluid_mesh_pt
 Pointer to Fluid_mesh. More...
 
Vector< TriangleMeshPolygon * > Drop_polygon_pt
 Vector storing pointer to the drop polygons. More...
 
TriangleMeshPolygonOuter_boundary_polyline_pt
 Triangle mesh polygon for outer boundary. More...
 
DataDrop_pressure_data_pt
 Pointer to a global drop pressure datum. More...
 
double Initial_value_for_drop_pressure
 Backed up drop pressure between adaptations. More...
 
ELEMENT * Hijacked_element_pt
 Pointer to hijacked element. More...
 
bool Use_volume_constraint
 Bool to indicate if volume constraint is applied (only for steady run) 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_convergence_check ()
 
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 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 TwoLayerInterfaceProblem< ELEMENT >

Problem class to simulate viscous drop propagating along 2D channel.

Member Enumeration Documentation

◆ anonymous enum

template<class ELEMENT >
anonymous enum
private

Enumeration of channel boundaries.

Enumerator
Inflow_boundary_id 
Upper_wall_boundary_id 
Outflow_boundary_id 
Bottom_wall_boundary_id 
Interface_boundary_id 
559  {
565  };
@ Interface_boundary_id
Definition: refineable_two_layer_interface.cc:564
@ Outflow_boundary_id
Definition: refineable_two_layer_interface.cc:562
@ Inflow_boundary_id
Definition: refineable_two_layer_interface.cc:560
@ Bottom_wall_boundary_id
Definition: refineable_two_layer_interface.cc:563
@ Upper_wall_boundary_id
Definition: refineable_two_layer_interface.cc:561

Constructor & Destructor Documentation

◆ TwoLayerInterfaceProblem()

template<class ELEMENT >
TwoLayerInterfaceProblem< ELEMENT >::TwoLayerInterfaceProblem

Constructor.

576 {
577  // Output directory
579 
580  // Allocate the timestepper -- this constructs the Problem's
581  // time object with a sufficient amount of storage to store the
582  // previous timsteps.
583  this->add_time_stepper_pt(new BDF<2>);
584 
585  // Build the boundary segments for outer boundary, consisting of
586  //--------------------------------------------------------------
587  // four separate polylines
588  //------------------------
589  Vector<TriangleMeshCurveSection*> boundary_polyline_pt(4);
590 
591  // Each polyline only has three vertices -- provide storage for their
592  // coordinates
593  Vector<Vector<double> > vertex_coord(3);
594  for(unsigned i=0;i<3;i++)
595  {
596  vertex_coord[i].resize(2);
597  }
598 
599  const double H = Problem_Parameter::Height;
600  const double h = Problem_Parameter::Interface_Height*H;
601 
602  // First polyline: Inflow
603  vertex_coord[0][0]=0.0;
604  vertex_coord[0][1]=0.0;
605  vertex_coord[1][0]=0.0;
606  vertex_coord[1][1]= h;
607  vertex_coord[2][0]=0.0;
608  vertex_coord[2][1]= H;
609 
610  // Build the 1st boundary polyline
611  boundary_polyline_pt[0] = new TriangleMeshPolyLine(vertex_coord,
613 
614  // Second boundary polyline: Upper wall
615  vertex_coord[0][0]=0.0;
616  vertex_coord[0][1]=H;
617  vertex_coord[1][0]=0.5*Problem_Parameter::Length;
618  vertex_coord[1][1]=H;
619  vertex_coord[2][0]=Problem_Parameter::Length;
620  vertex_coord[2][1]=H;
621 
622  // Build the 2nd boundary polyline
623  boundary_polyline_pt[1] = new TriangleMeshPolyLine(vertex_coord,
625 
626  // Third boundary polyline: Outflow
627  vertex_coord[0][0]=Problem_Parameter::Length;
628  vertex_coord[0][1]=H;
629  vertex_coord[1][0]=Problem_Parameter::Length;
630  vertex_coord[1][1]=h;
631  vertex_coord[2][0]=Problem_Parameter::Length;
632  vertex_coord[2][1]=0.0;
633 
634  // Build the 3rd boundary polyline
635  boundary_polyline_pt[2] = new TriangleMeshPolyLine(vertex_coord,
637 
638  // Fourth boundary polyline: Bottom wall
639  vertex_coord[0][0]=Problem_Parameter::Length;
640  vertex_coord[0][1]=0.0;
641  vertex_coord[1][0]=0.5*Problem_Parameter::Length;
642  vertex_coord[1][1]=0.0;
643  vertex_coord[2][0]=0.0;
644  vertex_coord[2][1]=0.0;
645 
646  // Build the 4th boundary polyline
647  boundary_polyline_pt[3] = new TriangleMeshPolyLine(vertex_coord,
649 
650  // Create the triangle mesh polygon for outer boundary
651  Outer_boundary_polyline_pt = new TriangleMeshPolygon(boundary_polyline_pt);
652 
653  //Here we need to put the dividing internal line in
654  Vector<TriangleMeshOpenCurve *> interface_pt(1);
655  //Set the vertex coordinates
656  vertex_coord[0][0]=0.0;
657  vertex_coord[0][1]=h;
658  vertex_coord[1][0]=0.5*Problem_Parameter::Length;
659  vertex_coord[1][1]=h;
660  vertex_coord[2][0]=Problem_Parameter::Length;
661  vertex_coord[2][1]=h;
662 
663 //Create the internal line
664  TriangleMeshPolyLine* interface_polyline_pt =
665  new TriangleMeshPolyLine(vertex_coord,
667 
668  // Do the connection with the destination boundary, in this case
669  // the connection is done with the inflow boundary
670  interface_polyline_pt->connect_initial_vertex_to_polyline(
671  dynamic_cast<TriangleMeshPolyLine*>(boundary_polyline_pt[0]),1);
672 
673  // Do the connection with the destination boundary, in this case
674  // the connection is done with the outflow boundary
675  interface_polyline_pt->connect_final_vertex_to_polyline(
676  dynamic_cast<TriangleMeshPolyLine*>(boundary_polyline_pt[2]),1);
677 
678  Vector<TriangleMeshCurveSection*> interface_curve_pt(1);
679  interface_curve_pt[0] = interface_polyline_pt;
680 
681  interface_pt[0] = new TriangleMeshOpenCurve(interface_curve_pt);
682 
683  // Now build the mesh, based on the boundaries specified by
684  //---------------------------------------------------------
685  // polygons just created
686  //----------------------
687 
688  // Convert to "closed curve" objects
690 
691  unsigned n_internal_closed_boundaries = 0;
693  inner_boundaries_pt(n_internal_closed_boundaries);
694 
695  // Target area for initial mesh
696  double uniform_element_area=0.01;
697 
698  // Use the TriangleMeshParameter object for gathering all
699  // the necessary arguments for the TriangleMesh object
700  TriangleMeshParameters triangle_mesh_parameters(
701  outer_closed_curve_pt);
702 
703  //Define the inner boundaries
704  triangle_mesh_parameters.internal_closed_curve_pt() = inner_boundaries_pt;
705 
706  // Define the holes on the boundary
707  triangle_mesh_parameters.internal_open_curves_pt() = interface_pt;
708 
709  // Define the maximum element area
710  triangle_mesh_parameters.element_area() =
711  uniform_element_area;
712 
713  Vector<double> lower_region(2);
714  lower_region[0] = 0.5*Problem_Parameter::Length;
715  lower_region[1] = 0.5*Problem_Parameter::Interface_Height;
716 
717  // Define the region
718  triangle_mesh_parameters.add_region_coordinates(1, lower_region);
719 
720  // Create the mesh
721  Fluid_mesh_pt =
722  new RefineableSolidTriangleMesh<ELEMENT>(
723  triangle_mesh_parameters, this->time_stepper_pt());
724 
725  // Set error estimator for bulk mesh
727  Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
728 
729  // Set targets for spatial adaptivity
730  Fluid_mesh_pt->max_permitted_error()=0.005;
731  Fluid_mesh_pt->min_permitted_error()=0.001;
732  Fluid_mesh_pt->max_element_size()=0.2;
733  Fluid_mesh_pt->min_element_size()=0.001;
734 
735  // Use coarser mesh during validation
737  {
738  Fluid_mesh_pt->min_element_size()=0.01;
739  }
740 
741  // Output boundary and mesh initial mesh for information
742  this->Fluid_mesh_pt->output_boundaries("boundaries.dat");
743  this->Fluid_mesh_pt->output("mesh.dat");
744 
745  // Set boundary condition and complete the build of all elements
747 
748  // Construct the mesh of free surface elements
751 
752  // Combine meshes
753  //---------------
754 
755  // Add Fluid_mesh_pt sub meshes
757 
758  // Add Free_surface sub meshes
759  this->add_sub_mesh(this->Free_surface_mesh_pt);
760 
761  // Build global mesh
762  this->build_global_mesh();
763 
764  // Setup equation numbering scheme
765  cout <<"Number of equations: " << this->assign_eqn_numbers() << std::endl;
766 
767 } // end_of_constructor
int i
Definition: BiCGSTAB_step_by_step.cpp:9
MatrixXf H
Definition: HessenbergDecomposition_matrixH.cpp:4
Mesh * Free_surface_mesh_pt
Pointers to mesh of free surface elements.
Definition: refineable_two_layer_interface.cc:533
void create_free_surface_elements()
Create free surface elements.
Definition: refineable_two_layer_interface.cc:775
void complete_problem_setup()
Set boundary conditions and complete the build of all elements.
Definition: refineable_two_layer_interface.cc:825
RefineableSolidTriangleMesh< ELEMENT > * Fluid_mesh_pt
Pointer to Fluid_mesh.
Definition: refineable_two_layer_interface.cc:537
TriangleMeshPolygon * Outer_boundary_polyline_pt
Triangle mesh polygon for outer boundary.
Definition: refineable_two_layer_interface.cc:543
void set_directory(const std::string &directory)
Definition: oomph_utilities.cc:298
Definition: mesh.h:67
void add_time_stepper_pt(TimeStepper *const &time_stepper_pt)
Definition: problem.cc:1545
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Definition: problem.h:1330
void build_global_mesh()
Definition: problem.cc:1493
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Definition: problem.cc:1989
TimeStepper *& time_stepper_pt()
Definition: problem.h:1524
Base class defining a closed curve for the Triangle mesh generation.
Definition: unstructured_two_d_mesh_geometry_base.h:1339
void connect_initial_vertex_to_polyline(TriangleMeshPolyLine *polyline_pt, const unsigned &vertex_number, const double &tolerance_for_connection=1.0e-14)
Definition: unstructured_two_d_mesh_geometry_base.cc:1120
void connect_final_vertex_to_polyline(TriangleMeshPolyLine *polyline_pt, const unsigned &vertex_number, const double &tolerance_for_connection=1.0e-14)
Definition: unstructured_two_d_mesh_geometry_base.cc:1189
Definition: unstructured_two_d_mesh_geometry_base.h:1642
Definition: triangle_mesh.template.h:94
Class defining a polyline for use in Triangle Mesh generation.
Definition: unstructured_two_d_mesh_geometry_base.h:868
Class defining a closed polygon for the Triangle mesh generation.
Definition: unstructured_two_d_mesh_geometry_base.h:1451
Definition: oomph-lib/src/generic/Vector.h:58
Definition: error_estimator.h:266
Z2ErrorEstimator * error_estimator_pt
Definition: MortaringCantileverCompareToNonMortaring.cpp:190
double Interface_Height
Definition: refineable_two_layer_interface.cc:329
DocInfo Doc_info
Doc info.
Definition: refineable_two_layer_interface.cc:291
double Length
Length of the channel.
Definition: refineable_two_layer_interface.cc:323
double Height
Definition: refineable_two_layer_interface.cc:326
bool command_line_flag_has_been_set(const std::string &flag)
Definition: oomph_utilities.cc:501

References oomph::TriangleMeshParameters::add_region_coordinates(), oomph::CommandLineArgs::command_line_flag_has_been_set(), oomph::TriangleMeshCurveSection::connect_final_vertex_to_polyline(), oomph::TriangleMeshCurveSection::connect_initial_vertex_to_polyline(), Problem_Parameter::Doc_info, oomph::TriangleMeshParameters::element_area(), MeshRefinement::error_estimator_pt, H, Problem_Parameter::Height, i, Problem_Parameter::Interface_Height, oomph::TriangleMeshParameters::internal_closed_curve_pt(), oomph::TriangleMeshParameters::internal_open_curves_pt(), Problem_Parameter::Length, and oomph::DocInfo::set_directory().

◆ ~TwoLayerInterfaceProblem()

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

Destructor.

361  {
362  // Fluid timestepper
363  delete this->time_stepper_pt(0);
364 
365  // Kill data associated with outer boundary
367  for (unsigned j=0;j<n;j++)
368  {
370  }
372 
373  // Flush element of free surface elements
375  delete Free_surface_mesh_pt;
376 
377  // Delete error estimator
378  delete Fluid_mesh_pt->spatial_error_estimator_pt();
379 
380  // Delete fluid mesh
381  delete Fluid_mesh_pt;
382 
383  // Delete the global pressure drop data
384  delete Drop_pressure_data_pt;
385 
386  // Kill const eqn
388 
389  }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
void delete_free_surface_elements()
Delete free surface elements.
Definition: refineable_two_layer_interface.cc:512
Data * Drop_pressure_data_pt
Pointer to a global drop pressure datum.
Definition: refineable_two_layer_interface.cc:546
unsigned npolyline() const
Number of constituent polylines.
Definition: unstructured_two_d_mesh_geometry_base.h:1482
TriangleMeshPolyLine * polyline_pt(const unsigned &i) const
Pointer to i-th constituent polyline.
Definition: unstructured_two_d_mesh_geometry_base.h:1488
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
Definition: jeffery_orbit.cc:82
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References Problem_Parameter::Constitutive_law_pt, j, and n.

Member Function Documentation

◆ actions_after_adapt()

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

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

Reimplemented from oomph::Problem.

406  {
407  // Create the elements that impose the displacement constraint
409 
410  // Rebuild the Problem's global mesh from its various sub-meshes
411  this->rebuild_global_mesh();
412 
413  // Setup the problem again -- remember that fluid mesh has been
414  // completely rebuilt and its element's don't have any
415  // pointers to Re etc. yet
417 
418  }// end of actions_after_adapt
void rebuild_global_mesh()
Definition: problem.cc:1533

◆ actions_after_newton_solve()

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

Update the after solve (empty)

Reimplemented from oomph::Problem.

422 {}

◆ actions_before_adapt()

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

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

Reimplemented from oomph::Problem.

394  {
395  // Kill the elements and wipe surface mesh
397 
398  // Rebuild the Problem's global mesh from its various sub-meshes
399  this->rebuild_global_mesh();
400 
401  }// end of actions_before_adapt

◆ actions_before_newton_solve()

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

Update the problem specs before solve.

Reimplemented from oomph::Problem.

426  {
427  //Reset the Lagrangian coordinates of the nodes to be the current
428  //Eulerian coordinates -- this makes the current configuration
429  //stress free
430  Fluid_mesh_pt->set_lagrangian_nodal_coordinates();
431  }

◆ complete_problem_setup()

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

Set boundary conditions and complete the build of all elements.

826  {
827  // Re-set the boundary conditions for fluid problem: All nodes are
828  // free by default -- just pin the ones that have Dirichlet conditions
829  // here.
830  unsigned nbound=Fluid_mesh_pt->nboundary();
831  for(unsigned ibound=0;ibound<nbound;ibound++)
832  {
833  unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
834  for (unsigned inod=0;inod<num_nod;inod++)
835  {
836  // Get node
837  Node* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
838 
839  //Pin both velocities on side boundaries (1 and 3)
840  if((ibound==1) || (ibound==3))
841  {
842  nod_pt->pin(0);
843  nod_pt->pin(1);
844  }
845 
846  //If it's the inflow or outflow pin only the horizontal velocity
847  if((ibound==0) || (ibound==2)) {nod_pt->pin(0);}
848 
849  // Pin horizontal pseudo-solid positions
850  SolidNode* solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
851  solid_node_pt->pin_position(0);
852 
853  //Pin the vertical position on the upper and lower walls
854  if((ibound==1) || (ibound==3))
855  {
856  solid_node_pt->pin_position(1);
857  }
858  else
859  {
860  solid_node_pt->unpin_position(1);
861  }
862  } //End of loop over nodes
863  } // end loop over boundaries
864 
865  // Complete the build of all elements so they are fully functional
866  // Remember that adaptation for triangle meshes involves a complete
867  // regneration of the mesh (rather than splitting as in tree-based
868  // meshes where such parameters can be passed down from the father
869  // element!)
870  unsigned n_element = Fluid_mesh_pt->nelement();
871  for(unsigned e=0;e<n_element;e++)
872  {
873  // Upcast from GeneralisedElement to the present element
874  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
875 
876  // Set the Reynolds number
877  el_pt->re_pt() = &Problem_Parameter::Re;
878 
879  // Set the Womersley number (same as Re since St=1)
880  el_pt->re_st_pt() = &Problem_Parameter::ReSt;
881 
882  // Set the product of the Reynolds number and the inverse of the
883  // Froude number
884  el_pt->re_invfr_pt() = &Problem_Parameter::ReInvFr;
885 
886  // Set the direction of gravity
887  el_pt->g_pt() = &Problem_Parameter::G;
888 
889  // Set the constitutive law for pseudo-elastic mesh deformation
890  el_pt->constitutive_law_pt()=Problem_Parameter::Constitutive_law_pt;
891  }
892 
893  //For the elements in the upper region (region 0),
894  //set the viscosity ratio
895  n_element = Fluid_mesh_pt->nregion_element(1);
896  for(unsigned e=0;e<n_element;e++)
897  {
898  // Upcast from GeneralisedElement to the present element
899  ELEMENT* el_pt =
900  dynamic_cast<ELEMENT*>(Fluid_mesh_pt->region_element_pt(1,e));
901 
902  el_pt->viscosity_ratio_pt() = &Problem_Parameter::Viscosity_Ratio;
903 
904  el_pt->density_ratio_pt() = &Problem_Parameter::Density_Ratio;
905  }
906 
907 
908  // Re-apply boundary values on Dirichlet boundary conditions
909  // (Boundary conditions are ignored when the solution is transferred
910  // from the old to the new mesh by projection; this leads to a slight
911  // change in the boundary values (which are, of course, never changed,
912  // unlike the actual unknowns for which the projected values only
913  // serve as an initial guess)
914 
915  // Set velocity and history values of velocity on walls
916  nbound=this->Fluid_mesh_pt->nboundary();
917  for(unsigned ibound=0;ibound<nbound;++ibound)
918  {
919  // Loop over nodes on this boundary
920  unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(ibound);
921  for (unsigned inod=0;inod<num_nod;inod++)
922  {
923  // Get node
924  Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(ibound,inod);
925 
926  // Get number of previous (history) values
927  unsigned n_prev=nod_pt->time_stepper_pt()->nprev_values();
928 
929  //If we are on the upper or lower walls
930  // Velocity is and was zero at all previous times
931  if((ibound==1) || (ibound==3))
932  {
933  //Loop over time history values
934  for (unsigned t=0;t<=n_prev;t++)
935  {
936  nod_pt->set_value(t,0,0.0);
937  nod_pt->set_value(t,1,0.0);
938 
939  // Nodes have always been there...
940  nod_pt->x(t,0)=nod_pt->x(0,0);
941  nod_pt->x(t,1)=nod_pt->x(0,1);
942  }
943  }
944  //If we are on the side wals there is no horizontal velocity or
945  //change in horizontal position
946  if((ibound==0) || (ibound==2))
947  {
948  for (unsigned t=0;t<=n_prev;t++)
949  {
950  nod_pt->set_value(t,0,0.0);
951  // Nodes have always been there...
952  nod_pt->x(t,0)=nod_pt->x(0,0);
953  }
954  }
955  //But there is a change in vertical position!
956  }
957  } //End of loop over boundaries
958 
959  fix_pressure(0,0,0.0);
960 
961  }
Array< double, 1, 3 > e(1./3., 0.5, 2.)
void fix_pressure(const unsigned &e, const unsigned &pdof, const double &pvalue)
Fix pressure in element e at pressure dof pdof and set to pvalue.
Definition: refineable_two_layer_interface.cc:495
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
void set_value(const unsigned &i, const double &value_)
Definition: nodes.h:271
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
Definition: nodes.h:1686
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1816
void unpin_position(const unsigned &i)
Unpin the nodal position.
Definition: nodes.h:1829
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
double ReSt
Womersley number (Reynolds x Strouhal)
Definition: refineable_two_layer_interface.cc:300
Vector< double > G(2)
Direction of gravity.
double Density_Ratio
Definition: refineable_two_layer_interface.cc:311
double ReInvFr
Product of Reynolds number and inverse of Froude number.
Definition: refineable_two_layer_interface.cc:303
double Viscosity_Ratio
Definition: refineable_two_layer_interface.cc:307
double Re
Reynolds number.
Definition: jeffery_orbit.cc:59
t
Definition: plotPSD.py:36

References Problem_Parameter::Constitutive_law_pt, Problem_Parameter::Density_Ratio, e(), Problem_Parameter::G, oomph::TimeStepper::nprev_values(), oomph::Data::pin(), oomph::SolidNode::pin_position(), Problem_Parameter::Re, Problem_Parameter::ReInvFr, Problem_Parameter::ReSt, oomph::Data::set_value(), plotPSD::t, oomph::Data::time_stepper_pt(), oomph::SolidNode::unpin_position(), Problem_Parameter::Viscosity_Ratio, and oomph::Node::x().

◆ compute_error_estimate()

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

1035 {
1036  // Get error estimator
1037  ErrorEstimator* err_est_pt=Fluid_mesh_pt->spatial_error_estimator_pt();
1038 
1039  // Get/output error estimates
1040  unsigned nel=Fluid_mesh_pt->nelement();
1041  Vector<double> elemental_error(nel);
1042 
1043  // We need a dynamic cast, get_element_errors from the Fluid_mesh_pt
1044  // Dynamic cast is used because get_element_errors require a Mesh* ans
1045  // not a SolidMesh*
1046  Mesh* fluid_mesh_pt=dynamic_cast<Mesh*>(Fluid_mesh_pt);
1047  err_est_pt->get_element_errors(fluid_mesh_pt,
1048  elemental_error);
1049 
1050  // Set errors for post-processing and find extrema
1051  max_err=0.0;
1052  min_err=DBL_MAX;
1053  for (unsigned e=0;e<nel;e++)
1054  {
1055  dynamic_cast<MyCrouzeixRaviartElement*>(Fluid_mesh_pt->element_pt(e))->
1056  set_error(elemental_error[e]);
1057 
1058  max_err=std::max(max_err,elemental_error[e]);
1059  min_err=std::min(min_err,elemental_error[e]);
1060  }
1061 
1062 }
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 CrouzeixRaviart element to modify output.
Definition: refineable_two_layer_interface.cc:54
#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 TwoLayerInterfaceProblem< 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

776 {
777 
778  //Loop over the free surface boundaries
779  unsigned nb=Fluid_mesh_pt->nboundary();
780  for(unsigned b=4;b<nb;b++)
781  {
782  // Note: region is important
783  // How many bulk fluid elements are adjacent to boundary b in region 0?
784  unsigned n_element = Fluid_mesh_pt->nboundary_element_in_region(b,0);
785 
786  // Loop over the bulk fluid elements adjacent to boundary b?
787  for(unsigned e=0;e<n_element;e++)
788  {
789  // Get pointer to the bulk fluid element that is
790  // adjacent to boundary b in region 0
791  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
792  Fluid_mesh_pt->boundary_element_in_region_pt(b,0,e));
793 
794  //Find the index of the face of element e along boundary b in region 0
795  int face_index = Fluid_mesh_pt->face_index_at_boundary_in_region(b,0,e);
796 
797  // Create new element
800  bulk_elem_pt,face_index);
801 
802  // Add it to the mesh
804 
805  //Add the appropriate boundary number
807 
808  //Specify the Strouhal number
809  el_pt->st_pt() = &Problem_Parameter::St;
810 
811  //Specify the capillary number
812  el_pt->ca_pt() = &Problem_Parameter::Ca;
813  }
814  }
815 }
Scalar * b
Definition: benchVecAdd.cpp:17
Specialise the elastic update template class to concrete 1D case.
Definition: specific_node_update_interface_elements.h:1220
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 add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:617
int nb
Definition: level2_impl.h:286
double St
Strouhal number.
Definition: jeffery_orbit.cc:62
double Ca
Capillary number.
Definition: refineable_two_layer_interface.cc:314

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

◆ deform_interface()

template<class ELEMENT >
void TwoLayerInterfaceProblem< ELEMENT >::deform_interface ( const double epsilon,
const unsigned n_periods 
)
inline

Function to deform the interface.

471  {
472  // Determine number of nodes in the "bulk" mesh
473  const unsigned n_node = Fluid_mesh_pt->nnode();
474 
475  // Loop over all nodes in mesh
476  for(unsigned n=0;n<n_node;n++)
477  {
478  // Determine eulerian position of node
479  const double current_x_pos = Fluid_mesh_pt->node_pt(n)->x(0);
480  const double current_y_pos = Fluid_mesh_pt->node_pt(n)->x(1);
481 
482  // Determine new vertical position of node
483  const double new_y_pos = current_y_pos
484  + (1.0-fabs(1.0-current_y_pos))*epsilon
485  *(cos(2.0*n_periods*MathematicalConstants::Pi*current_x_pos/Problem_Parameter::Length));
486 
487  // Set new position
488  Fluid_mesh_pt->node_pt(n)->x(1) = new_y_pos;
489  }
490 } // End of deform_interface
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
double Pi
Definition: two_d_biharmonic.cc:235
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
double epsilon
Definition: osc_ring_sarah_asymptotics.h:43

References cos(), oomph::SarahBL::epsilon, boost::multiprecision::fabs(), Problem_Parameter::Length, n, and BiharmonicTestFunctions2::Pi.

◆ delete_free_surface_elements()

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

Delete free surface elements.

513  {
514  // How many surface elements are in the surface mesh
515  unsigned n_element = Free_surface_mesh_pt->nelement();
516 
517  // Loop over the surface elements
518  for(unsigned e=0;e<n_element;e++)
519  {
520  // Kill surface element
522  }
523 
524 
525  // Wipe the mesh
527 
528 
529  } // 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().

◆ doc_solution()

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

Doc the solution.

970 {
971 
972  // Output the time
973  cout << "Time is now " << time_pt()->time() << std::endl;
974 
975  double min_x_coordinate = 2.0;
976  unsigned min_boundary_node = 0;
977  //Find the left-most node on the boundary
978  unsigned n_bound = Fluid_mesh_pt->nboundary_node(4);
979  for(unsigned n=0;n<n_bound;++n)
980  {
981  Node* nod_pt = Fluid_mesh_pt->boundary_node_pt(4,n);
982  if(nod_pt->x(0) < min_x_coordinate)
983  {
984  min_x_coordinate = nod_pt->x(0);
985  min_boundary_node = n;
986  }
987  }
988 
989  // Document time and vertical position of left hand side of interface
990  // in trace file
992  << time_pt()->time() << " "
993  << Fluid_mesh_pt->boundary_node_pt(4,min_boundary_node)->x(1) << std::endl;
994 
995  ofstream some_file;
996  char filename[100];
997 
998  // Set number of plot points (in each coordinate direction)
999  const unsigned npts = 5;
1000 
1001  // Open solution output file
1002  sprintf(filename,"%s/soln%i.dat",
1003  Problem_Parameter::Doc_info.directory().c_str(),
1004  Problem_Parameter::Doc_info.number());
1005  some_file.open(filename);
1006 
1007  // Output solution to file
1008  Fluid_mesh_pt->output(some_file,npts);
1009 
1010  // Close solution output file
1011  some_file.close();
1012 
1013  // Open interface solution output file
1014  sprintf(filename,"%s/interface_soln%i.dat",
1015  Problem_Parameter::Doc_info.directory().c_str(),
1016  Problem_Parameter::Doc_info.number());
1017  some_file.open(filename);
1018 
1019  // Output solution to file
1020  Free_surface_mesh_pt->output(some_file,npts);
1021 
1022  // Close solution output file
1023  some_file.close();
1024 
1025  // Increment the doc_info number
1027 }
unsigned & number()
Number used (e.g.) for labeling output files.
Definition: oomph_utilities.h:554
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
string filename
Definition: MergeRestartFiles.py:39
ofstream Trace_file
Trace file.
Definition: refineable_two_layer_interface.cc:335

References Problem_Parameter::Doc_info, MergeRestartFiles::filename, n, oomph::DocInfo::number(), Problem_Parameter::Trace_file, and oomph::Node::x().

◆ fix_pressure()

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

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

498  {
499  // Fix the pressure at that element
500  dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e))->
501  fix_pressure(pdof,pvalue);
502  }

References e().

◆ set_initial_condition()

template<class ELEMENT >
void TwoLayerInterfaceProblem< ELEMENT >::set_initial_condition ( )
inlinevirtual

Set the initial conditions.

Reimplemented from oomph::Problem.

445 {
446  // Determine number of nodes in mesh
447  const unsigned n_node = Fluid_mesh_pt->nnode();
448 
449  // Loop over all nodes in mesh
450  for(unsigned n=0;n<n_node;n++)
451  {
452  // Loop over the two velocity components
453  for(unsigned i=0;i<2;i++)
454  {
455  // Set velocity component i of node n to zero
456  Fluid_mesh_pt->node_pt(n)->set_value(i,0.0);
457  }
458  }
459 
460  // Initialise the previous velocity values for timestepping
461  // corresponding to an impulsive start
463 
464 } // End of set_initial_condition
void assign_initial_values_impulsive()
Definition: problem.cc:11499

References i, and n.

Member Data Documentation

◆ Drop_polygon_pt

template<class ELEMENT >
Vector<TriangleMeshPolygon*> TwoLayerInterfaceProblem< ELEMENT >::Drop_polygon_pt
private

Vector storing pointer to the drop polygons.

◆ Drop_pressure_data_pt

template<class ELEMENT >
Data* TwoLayerInterfaceProblem< ELEMENT >::Drop_pressure_data_pt
private

Pointer to a global drop pressure datum.

◆ Fluid_mesh_pt

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

Pointer to Fluid_mesh.

◆ Free_surface_mesh_pt

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

Pointers to mesh of free surface elements.

◆ Hijacked_element_pt

template<class ELEMENT >
ELEMENT* TwoLayerInterfaceProblem< ELEMENT >::Hijacked_element_pt
private

Pointer to hijacked element.

◆ Initial_value_for_drop_pressure

template<class ELEMENT >
double TwoLayerInterfaceProblem< ELEMENT >::Initial_value_for_drop_pressure
private

Backed up drop pressure between adaptations.

◆ Outer_boundary_polyline_pt

template<class ELEMENT >
TriangleMeshPolygon* TwoLayerInterfaceProblem< ELEMENT >::Outer_boundary_polyline_pt
private

Triangle mesh polygon for outer boundary.

◆ Use_volume_constraint

template<class ELEMENT >
bool TwoLayerInterfaceProblem< ELEMENT >::Use_volume_constraint
private

Bool to indicate if volume constraint is applied (only for steady run)


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