40 #include "navier_stokes.h" 41 #include "fluid_interface.h" 42 #include "constitutive.h" 46 #include "meshes/single_layer_spine_mesh.h" 51 using namespace oomph;
67 Vector<double> &normal)
80 template<
class ELEMENT>
95 void parameter_study(
const string& dir_name);
101 void create_volume_constraint_elements();
104 void doc_solution(DocInfo& doc_info);
150 template<
class ELEMENT>
155 Angle(0.5*MathematicalConstants::Pi)
169 double half_width=0.5;
173 new SingleLayerSpineMesh<SpineElement<ELEMENT> >(nx,nh,half_width,1.0);
184 unsigned n_boundary_element =
Bulk_mesh_pt->nboundary_element(2);
185 for(
unsigned e=0;e<n_boundary_element;e++)
188 FiniteElement *interface_element_pt =
189 new SpineLineFluidInterfaceElement<SpineElement<ELEMENT> >(
197 unsigned n_node = interface_element_pt->nnode();
200 if(interface_element_pt->node_pt(n_node-1)->is_on_boundary(1))
203 FiniteElement* point_element_pt =
204 dynamic_cast<SpineLineFluidInterfaceElement<
205 SpineElement<ELEMENT>
>*>(interface_element_pt)
206 ->make_bounding_element(1);
234 Bulk_mesh_pt->element_pt(0))->hijack_internal_value(0,0);
253 dynamic_cast<ELEMENT*
>(
Bulk_mesh_pt->element_pt(0))->fix_pressure(0,0.0);
263 for(
unsigned e=0;e<n_interface;e++)
266 SpineLineFluidInterfaceElement<SpineElement<ELEMENT> >*el_pt=
267 dynamic_cast<SpineLineFluidInterfaceElement<SpineElement<ELEMENT>
>*>
271 el_pt->ca_pt() = &
Ca;
284 for (
unsigned b=0;b<n_bound;b++)
289 unsigned n_boundary_node =
Bulk_mesh_pt->nboundary_node(b);
291 for(
unsigned n=0;n<n_boundary_node;n++)
304 FluidInterfaceBoundingElement *contact_angle_element_pt
305 =
dynamic_cast<FluidInterfaceBoundingElement*
>(
308 contact_angle_element_pt->set_contact_angle(&
Angle);
309 contact_angle_element_pt->ca_pt() = &
Ca;
310 contact_angle_element_pt->wall_unit_normal_fct_pt()
321 this->build_global_mesh();
324 cout <<
"Number of unknowns: " << assign_eqn_numbers() << std::endl;
333 template<
class ELEMENT>
338 for(
unsigned e=0;e<n_element;e++)
353 for(
unsigned e=0;e<n_element;e++)
362 for(
unsigned e=0;e<n_element;e++)
377 template<
class ELEMENT>
382 VolumeConstraintElement* vol_constraint_element =
387 for(
unsigned b=0;b<4;b++)
390 unsigned n_element =
Bulk_mesh_pt->nboundary_element(b);
393 for(
unsigned e=0;e<n_element;e++)
397 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
401 int face_index =
Bulk_mesh_pt->face_index_at_boundary(b,e);
404 SpineLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
405 new SpineLineVolumeConstraintBoundingElement<ELEMENT>(
406 bulk_elem_pt,face_index);
409 el_pt->set_volume_constraint_element(vol_constraint_element);
422 template<
class ELEMENT>
428 doc_info.set_directory(dir_name);
434 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
436 Trace_file <<
"VARIABLES=\"<greek>a</greek><sub>prescribed</sub>\",";
437 Trace_file <<
"\"h<sub>left</sub>\",\"h<sub>right</sub>\",";
438 Trace_file <<
"\"p<sub>fluid</sub>-p<sub>ext</sub>\",";
439 Trace_file <<
"\"<greek>D</greek>p<sub>exact</sub>\"";
443 for(
unsigned i=0;i<6;i++)
446 steady_newton_solve();
455 Angle -= 5.0*MathematicalConstants::Pi/180.0;
466 template<
class ELEMENT>
478 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
480 some_file.open(filename);
507 template<
class ELEMENT>
529 void create_free_surface_elements();
535 void create_contact_angle_element();
587 template<
class ELEMENT>
592 Angle(0.5*MathematicalConstants::Pi)
606 double half_width=0.5;
609 Bulk_mesh_pt =
new ElasticRectangularQuadMesh<ELEMENT>(nx,nh,half_width,1.0);
632 Bulk_mesh_pt->element_pt(0))->hijack_internal_value(0,0);
651 dynamic_cast<ELEMENT*
>(
Bulk_mesh_pt->element_pt(0))->fix_pressure(0,0.0);
660 for(
unsigned e=0;e<n_bulk;e++)
675 for (
unsigned b=0;b<n_bound;b++)
680 unsigned n_boundary_node =
Bulk_mesh_pt->nboundary_node(b);
682 for(
unsigned n=0;n<n_boundary_node;n++)
694 for (
unsigned b=0;b<n_bound;b++)
699 unsigned n_boundary_node =
Bulk_mesh_pt->nboundary_node(b);
701 for(
unsigned n=0;n<n_boundary_node;n++)
706 static_cast<SolidNode*
>(
Bulk_mesh_pt->boundary_node_pt(b,n))
712 static_cast<SolidNode*
>(
Bulk_mesh_pt->boundary_node_pt(b,n))
722 for(
unsigned n=0;n<n_node;n++)
724 static_cast<SolidNode*
>(
Bulk_mesh_pt->node_pt(n))->pin_position(0);
745 this->build_global_mesh();
748 cout <<
"Number of unknowns: " << assign_eqn_numbers() << std::endl;
757 template<
class ELEMENT>
768 for(
unsigned e=0;e<n_element;e++)
776 for(
unsigned e=0;e<n_element;e++)
798 template<
class ELEMENT>
808 unsigned n_element =
Bulk_mesh_pt->nboundary_element(b);
811 for(
unsigned e=0;e<n_element;e++)
815 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
819 int face_index =
Bulk_mesh_pt->face_index_at_boundary(b,e);
822 ElasticLineFluidInterfaceElement<ELEMENT>* el_pt =
823 new ElasticLineFluidInterfaceElement<ELEMENT>(
824 bulk_elem_pt,face_index);
830 el_pt->set_boundary_number_in_bulk_mesh(b);
833 el_pt->ca_pt() = &
Ca;
845 template<
class ELEMENT>
850 VolumeConstraintElement* vol_constraint_element =
858 for(
unsigned b=0;b<4;b++)
861 unsigned n_element =
Bulk_mesh_pt->nboundary_element(b);
864 for(
unsigned e=0;e<n_element;e++)
868 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
872 int face_index =
Bulk_mesh_pt->face_index_at_boundary(b,e);
875 ElasticLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
876 new ElasticLineVolumeConstraintBoundingElement<ELEMENT>(
877 bulk_elem_pt,face_index);
880 el_pt->set_volume_constraint_element(vol_constraint_element);
891 template<
class ELEMENT>
902 FluidInterfaceBoundingElement* el_pt =
903 dynamic_cast<ElasticLineFluidInterfaceElement<ELEMENT>*
> 905 make_bounding_element(1);
908 el_pt->set_contact_angle(&
Angle);
911 el_pt->ca_pt() = &
Ca;
914 el_pt->wall_unit_normal_fct_pt()
928 template<
class ELEMENT>
933 doc_info.set_directory(dir_name);
938 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
940 Trace_file <<
"VARIABLES=\"<greek>a</greek><sub>prescribed</sub>\",";
941 Trace_file <<
"\"h<sub>left</sub>\",\"h<sub>right</sub>\",";
942 Trace_file <<
"\"p<sub>fluid</sub>-p<sub>ext</sub>\",";
943 Trace_file <<
"\"<greek>D</greek>p<sub>exact</sub>\"";
947 for(
unsigned i=0;i<6;i++)
950 steady_newton_solve();
959 Angle -= 5.0*MathematicalConstants::Pi/180.0;
970 template<
class ELEMENT>
981 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
983 some_file.open(filename);
1001 ->node_pt(np-1)->x(1);
1003 <<
dynamic_cast<ELEMENT*
>(
Bulk_mesh_pt->element_pt(0))->p_nst(0)-
1021 for (
unsigned i=0;i<2;i++)
1024 bool hijack_internal=
false;
1025 if (i==1) hijack_internal=
true;
1029 string dir_name=
"RESLT_hijacked_external";
1030 if (i==1) dir_name=
"RESLT_hijacked_internal";
1040 for (
unsigned i=0;i<2;i++)
1042 bool hijack_internal=
false;
1043 if (i==1) hijack_internal=
true;
1046 PseudoSolidNodeUpdateElement<QCrouzeixRaviartElement<2>,
1047 QPVDElementWithPressure<2> > > > problem(hijack_internal);
1049 string dir_name=
"RESLT_elastic_hijacked_external";
1050 if (i==1) dir_name=
"RESLT_elastic_hijacked_internal";
1053 problem.parameter_study(dir_name);
void parameter_study(const string &dir_name)
double Pext
The external pressure.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
void create_contact_angle_element()
Create the contact angle element.
Vector< double > Wall_normal
Direction of the wall normal vector.
void actions_before_newton_convergence_check()
Update the spine mesh after every Newton step.
Data * Traded_pressure_data_pt
double Pext
The external pressure.
void parameter_study()
Peform a parameter study: Solve problem for a range of contact angles.
void create_free_surface_elements()
Create the free surface elements.
void parameter_study(const string &dir_name)
double Angle
The contact angle.
Mesh * Free_surface_bounding_mesh_pt
Storage for the element bounding the free surface.
void create_volume_constraint_elements()
Create the volume constraint elements.
Mesh * Surface_mesh_pt
The mesh for the interface elements.
Mesh * Free_surface_mesh_pt
Storage for the free surface mesh.
ofstream Trace_file
Trace file.
void create_volume_constraint_elements()
Create the volume constraint elements.
double Volume
The prescribed volume of the fluid.
Data * External_pressure_data_pt
Data object whose single value stores the external pressure.
void wall_unit_normal_fct(const Vector< double > &x, Vector< double > &normal)
Function that specifies the wall unit normal.
SingleLayerSpineMesh< SpineElement< ELEMENT > > * Bulk_mesh_pt
The bulk mesh of fluid elements.
Mesh * Volume_constraint_mesh_pt
Storage for the volume constraint.
CapProblem(const bool &hijack_internal)
~PseudoSolidCapProblem()
Destructor: clean up memory allocated by the object.
ofstream Trace_file
Trace file.
PseudoSolidCapProblem(const bool &hijack_internal)
void doc_solution(DocInfo &doc_info)
Doc the solution.
Namespace for phyical parameters.
A class that solves the Navier–Stokes equations to compute the shape of a static interface in a rect...
Mesh * Volume_constraint_mesh_pt
The volume constraint mesh.
void doc_solution(DocInfo &doc_info)
Doc the solution.
int main()
Main driver: Build problem and initiate parameter study.
double Ca
The Capillary number.
double Volume
The volume of the fluid.
double Ca
The Capillary number.
Mesh * Volume_computation_mesh_pt
Storage for the elements that compute the enclosed volume.
Data * External_pressure_data_pt
Data object whose single value stores the external pressure.
Data * Traded_pressure_data_pt
Data that is traded for the volume constraint.
double Nu
Pseudo-solid Poisson ratio.
Mesh * Bulk_mesh_pt
Storage for the bulk mesh.
double Angle
The contact angle.
Mesh * Point_mesh_pt
The mesh for the element at the contact point.