32 #include "young_laplace.h"
35 #include "meshes/rectangular_quadmesh.h"
39 using namespace oomph;
49 template<
class ELEMENT>
71 if (Contact_angle_mesh_pt!=0) delete_contact_angle_elements();
74 rebuild_global_mesh();
85 create_contact_angle_elements(1);
86 create_contact_angle_elements(3);
89 unsigned nel=Contact_angle_mesh_pt->nelement();
90 for (
unsigned e=0;e<nel;e++)
94 YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
95 dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*
>(
96 Contact_angle_mesh_pt->element_pt(e));
104 rebuild_global_mesh();
109 void increment_parameters();
125 RefineableRectangularQuadMesh<ELEMENT>* Bulk_mesh_pt;
128 Mesh* Contact_angle_mesh_pt;
131 Mesh* Height_control_mesh_pt;
137 Node* Control_node_pt;
145 template<
class ELEMENT>
169 cout <<
"Lx = " << l_x <<
" and Ly = " << l_y << endl;
172 Bulk_mesh_pt=
new RefineableRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
175 Bulk_mesh_pt->spatial_error_estimator_pt()=
new Z2ErrorEstimator;
178 Bulk_mesh_pt->max_permitted_error()=1.0e-4;
179 Bulk_mesh_pt->min_permitted_error()=1.0e-6;
182 add_sub_mesh(Bulk_mesh_pt);
189 ELEMENT* prescribed_height_element_pt=
dynamic_cast<ELEMENT*
>(
195 Control_node_pt=
static_cast<Node*
>(
196 prescribed_height_element_pt->node_pt(0));
198 cout <<
"Controlling height at (x,y) : (" << Control_node_pt->x(0)
199 <<
"," << Control_node_pt->x(1) <<
")" << endl;
203 Height_control_element_pt=0;
204 Height_control_mesh_pt=0;
207 Height_control_element_pt=
new HeightControlElement(
211 Height_control_element_pt->kappa_pt()->
215 Height_control_mesh_pt =
new Mesh;
216 Height_control_mesh_pt->add_element_pt(Height_control_element_pt);
219 add_sub_mesh(Height_control_mesh_pt);
237 Contact_angle_mesh_pt=0;
242 Contact_angle_mesh_pt=
new Mesh;
245 create_contact_angle_elements(1);
246 create_contact_angle_elements(3);
249 add_sub_mesh(Contact_angle_mesh_pt);
264 unsigned n_bound = Bulk_mesh_pt->nboundary();
265 for(
unsigned b=0;b<n_bound;b++)
273 unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
274 for (
unsigned n=0;n<n_node;n++)
276 Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
285 unsigned n_bulk=Bulk_mesh_pt->nelement();
286 for(
unsigned i=0;i<n_bulk;i++)
289 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(Bulk_mesh_pt->element_pt(i));
308 unsigned nel=Contact_angle_mesh_pt->nelement();
309 for (
unsigned e=0;e<nel;e++)
313 YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
314 dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*
>(
315 Contact_angle_mesh_pt->element_pt(e));
323 cout <<
"\nNumber of equations: " << assign_eqn_numbers() << endl;
324 cout <<
"\n********************************************\n" << endl;
333 template<
class ELEMENT>
338 unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
341 for(
unsigned e=0;e<n_element;e++)
344 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
345 Bulk_mesh_pt->boundary_element_pt(b,e));
348 int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
351 YoungLaplaceContactAngleElement<ELEMENT>* contact_angle_element_pt =
new
352 YoungLaplaceContactAngleElement<ELEMENT>(bulk_elem_pt,face_index);
355 Contact_angle_mesh_pt->add_element_pt(contact_angle_element_pt);
365 template<
class ELEMENT>
370 unsigned n_element = Contact_angle_mesh_pt->nelement();
373 for(
unsigned e=0;e<n_element;e++)
376 delete Contact_angle_mesh_pt->element_pt(e);
380 Contact_angle_mesh_pt->flush_element_and_node_storage();
390 template<
class ELEMENT>
401 cout <<
"Solving for Prescribed KAPPA Value = " ;
409 cout <<
"Solving for Prescribed HEIGHT Value = " ;
419 template<
class ELEMENT>
421 ofstream& trace_file)
428 trace_file << Control_node_pt->value(0) ;
439 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
441 some_file.open(filename);
442 Bulk_mesh_pt->output(some_file,npts);
453 ofstream tangent_file;
454 sprintf(filename,
"%s/tangent_to_contact_line%i.dat",
455 doc_info.directory().c_str(),
457 tangent_file.open(filename);
459 ofstream normal_file;
460 sprintf(filename,
"%s/normal_to_contact_line%i.dat",
461 doc_info.directory().c_str(),
463 normal_file.open(filename);
466 ofstream contact_angle_file;
467 sprintf(filename,
"%s/contact_angle%i.dat",
468 doc_info.directory().c_str(),
470 contact_angle_file.open(filename);
473 Vector<double> tangent(3);
474 Vector<double> normal(3);
475 Vector<double> r_contact(3);
478 unsigned n_element = Contact_angle_mesh_pt->nelement();
481 for(
unsigned e=0;e<n_element;e++)
484 tangent_file <<
"ZONE" << std::endl;
485 normal_file <<
"ZONE" << std::endl;
486 contact_angle_file <<
"ZONE" << std::endl;
489 YoungLaplaceContactAngleElement<ELEMENT>* el_pt =
490 dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*
>(
491 Contact_angle_mesh_pt->element_pt(e));
495 for (
unsigned i=0;i<npts;i++)
497 s[0]=-1.0+2.0*double(i)/double(npts-1);
499 dynamic_cast<ELEMENT*
>(el_pt->bulk_element_pt())->
500 position(el_pt->local_coordinate_in_bulk(s),r_contact);
502 el_pt->contact_line_vectors(s,tangent,normal);
503 tangent_file << r_contact[0] <<
" "
504 << r_contact[1] <<
" "
505 << r_contact[2] <<
" "
508 << tangent[2] <<
" " << std::endl;
510 normal_file << r_contact[0] <<
" "
511 << r_contact[1] <<
" "
512 << r_contact[2] <<
" "
515 << normal[2] <<
" " << std::endl;
517 contact_angle_file << r_contact[1] <<
" "
518 << el_pt->actual_cos_contact_angle(s)
525 tangent_file.close();
527 contact_angle_file.close();
531 cout <<
"\n********************************************" << endl << endl;
541 void run_it(
const string& output_directory)
555 doc_info.set_directory(output_directory);
559 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
560 trace_file.open(filename);
564 "VARIABLES=\"<GREEK>k</GREEK>\",\"<GREEK>k</GREEK>_{ex}\",\"h\""
566 trace_file <<
"ZONE" << std::endl;
575 problem.refine_uniformly();
594 unsigned max_adapt=1;
595 problem.newton_solve(max_adapt);
616 int main(
int argc,
char* argv[])
620 CommandLineArgs::setup(argc,argv);
628 if (CommandLineArgs::Argc==1)
631 <<
"Running every case with limited number of steps for validation"
640 case_lo=atoi(argv[1]);
641 case_hi=atoi(argv[1]);
649 for (
unsigned my_case=case_lo;my_case<=case_hi;my_case++)
659 <<
"//////////////////////////////////////////////////////////\n"
660 <<
"All pinned solution \n"
661 <<
"//////////////////////////////////////////////////////////\n\n";
667 run_it(
"RESLT_adapt_all_pinned");
674 <<
"//////////////////////////////////////////////////////////\n"
675 <<
"Barrel-shaped solution with spine \n"
676 <<
"/////////////////////////////////////////////////////////\n\n";
681 run_it(
"RESLT_adapt_barrel_shape");
688 <<
"//////////////////////////////////////////////////////////\n"
689 <<
"Barrel-shaped solution without spines \n"
690 <<
"/////////////////////////////////////////////////////////\n\n";
695 run_it(
"RESLT_adapt_barrel_shape_without_spines");
702 <<
"//////////////////////////////////////////////////////////\n"
703 <<
"T-junction solution \n"
704 <<
"//////////////////////////////////////////////////////////\n\n";
714 run_it(
"RESLT_adapt_T_junction");
721 std::cout <<
"Wrong case! Options are:\n"
722 <<
"0: adaptive All pinned\n"
723 <<
"1: adaptive Barrel with spines\n"
724 <<
"2: adaptive Barrel without spines\n"
725 <<
"3: adaptive T_junction\n"
2D RefineableYoungLaplace problem on rectangular domain, discretised with 2D QRefineableYoungLaplace ...
void create_contact_angle_elements(const unsigned &b)
Create YoungLaplace contact angle elements on the b-th boundary of the bulk mesh and add them to cont...
void actions_after_newton_solve()
Update the problem after solve: Empty.
~RefineableYoungLaplaceProblem()
Destructor (empty)
void actions_before_newton_solve()
Update the problem specs before solve: Empty.
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to and the tra...
RefineableYoungLaplaceProblem()
Constructor:
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of contact angle elements.
void delete_contact_angle_elements()
Delete contact angle elements.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of contact angle elements.
HeightControlElement * Height_control_element_pt
Pointer to height control element.
void increment_parameters()
Increase the problem parameters before each solve.
double L_x
Length and width of the domain.
double Controlled_height
Height control value.
unsigned Control_element
Number of element in bulk mesh at which height control is applied. Initialise to 0 – will be overwrit...
double get_exact_kappa()
Exact kappa.
bool Use_height_control
Use height control (true) or not (false)?
double Controlled_height_increment
Increment for height control.
double Kappa_increment
Increment for prescribed curvature.
bool Use_spines
Use spines (true) or not (false)
void spine_function(const Vector< double > &x, Vector< double > &spine, Vector< Vector< double > > &dspine)
Spine: The spine vector field as a function of the two coordinates x_1 and x_2, and its derivatives w...
unsigned Nsteps
Number of steps.
unsigned N_x
Number of elements in the mesh.
Data * Kappa_pt
Pointer to Data object that stores the prescribed curvature.
void spine_base_function(const Vector< double > &x, Vector< double > &spine_B, Vector< Vector< double > > &dspine_B)
Spine basis: The position vector to the basis of the spine as a function of the two coordinates x_1 a...
double L_y
Width of domain.
double Gamma
Contact angle and its cos (dependent parameter – is reassigned)
@ T_junction_with_nonzero_contact_angle
double Kappa_initial
Initial value for kappa.
double Cos_gamma
Cos of contact angle.
int Case
What case are we considering: Choose one from the enumeration Cases.
void setup_dependent_parameters_and_sanity_check()
Setup dependent parameters and perform sanity check.
int main(int argc, char *argv[])
Driver code for 2D RefineableYoungLaplace problem. Input arguments: none (for validation) or case (0,...
void run_it(const string &output_directory)
Run code for current setting of parameter values – specify name of output directory.