32 #include "young_laplace.h"
35 #include "meshes/simple_rectangular_quadmesh.h"
39 using namespace oomph;
52 template<
class ELEMENT>
78 void create_contact_angle_elements(
const unsigned& b);
92 Node* Control_node_pt;
100 template<
class ELEMENT>
123 cout <<
"Lx = " << l_x <<
" and Ly = " << l_y << endl;
126 Problem::mesh_pt()=
new SimpleRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
130 ELEMENT* prescribed_height_element_pt=
dynamic_cast<ELEMENT*
>(
136 Control_node_pt=
static_cast<Node*
>(
137 prescribed_height_element_pt->node_pt(0));
141 HeightControlElement* height_control_element_pt=0;
144 cout <<
"Controlling height at (x,y) : (" << Control_node_pt->x(0)
145 <<
"," << Control_node_pt->x(1) <<
")" <<
"\n" << endl;
147 height_control_element_pt=
new HeightControlElement(
153 height_control_element_pt->kappa_pt()
166 N_bulk_elements = mesh_pt()->nelement();
173 create_contact_angle_elements(1);
174 Last_element_on_boundary1=mesh_pt()->nelement();
176 create_contact_angle_elements(3);
177 Last_element_on_boundary3=mesh_pt()->nelement();
181 Last_element_on_boundary1=0;
182 Last_element_on_boundary3=0;
191 unsigned n_bound = mesh_pt()->nboundary();
192 for(
unsigned b=0;b<n_bound;b++)
200 unsigned n_node = mesh_pt()->nboundary_node(b);
201 for (
unsigned n=0;n<n_node;n++)
203 mesh_pt()->boundary_node_pt(b,n)->pin(0);
212 for(
unsigned i=0;i<N_bulk_elements;i++)
215 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(mesh_pt()->element_pt(i));
235 for (
unsigned e=N_bulk_elements;e<Last_element_on_boundary1;e++)
238 YoungLaplaceContactAngleElement<ELEMENT>* el_pt =
239 dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*
>(
240 mesh_pt()->element_pt(e));
245 for (
unsigned e=Last_element_on_boundary1;e<Last_element_on_boundary3;e++)
248 YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
249 dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*
>(
250 mesh_pt()->element_pt(e));
263 mesh_pt()->add_element_pt(height_control_element_pt);
268 cout <<
"\nNumber of equations: " << assign_eqn_numbers() << endl;
269 cout <<
"\n********************************************\n" << endl;
278 template<
class ELEMENT>
283 unsigned n_element = mesh_pt()->nboundary_element(b);
286 for(
unsigned e=0;e<n_element;e++)
289 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
290 mesh_pt()->boundary_element_pt(b,e));
293 int face_index = mesh_pt()->face_index_at_boundary(b,e);
296 YoungLaplaceContactAngleElement<ELEMENT>* contact_angle_element_pt =
new
297 YoungLaplaceContactAngleElement<ELEMENT>(bulk_elem_pt,face_index);
300 mesh_pt()->add_element_pt(contact_angle_element_pt);
313 template<
class ELEMENT>
324 cout <<
"Solving for Prescribed KAPPA Value = " ;
332 cout <<
"Solving for Prescribed HEIGHT Value = " ;
344 template<
class ELEMENT>
346 ofstream& trace_file)
353 trace_file << Control_node_pt->value(0) ;
363 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
365 some_file.open(filename);
366 for(
unsigned i=0;i<N_bulk_elements;i++)
368 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(mesh_pt()->element_pt(i));
369 el_pt->output(some_file,npts);
381 ofstream tangent_file;
382 sprintf(filename,
"%s/tangent_to_contact_line%i.dat",
383 doc_info.directory().c_str(),
385 tangent_file.open(filename);
387 ofstream normal_file;
388 sprintf(filename,
"%s/normal_to_contact_line%i.dat",
389 doc_info.directory().c_str(),
391 normal_file.open(filename);
394 ofstream contact_angle_file;
395 sprintf(filename,
"%s/contact_angle%i.dat",
396 doc_info.directory().c_str(),
398 contact_angle_file.open(filename);
401 Vector<double> tangent(3);
402 Vector<double> normal(3);
403 Vector<double> r_contact(3);
407 for (
unsigned bnd=0;bnd<2;bnd++)
409 unsigned el_lo=N_bulk_elements;
410 unsigned el_hi=Last_element_on_boundary1;
413 el_lo=Last_element_on_boundary1;
414 el_hi=Last_element_on_boundary3;
419 for (
unsigned e=el_lo;e<el_hi;e++)
422 tangent_file <<
"ZONE" << std::endl;
423 normal_file <<
"ZONE" << std::endl;
424 contact_angle_file <<
"ZONE" << std::endl;
427 YoungLaplaceContactAngleElement<ELEMENT>* el_pt =
428 dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*
>(
429 mesh_pt()->element_pt(e));
433 for (
unsigned i=0;i<npts;i++)
435 s[0]=-1.0+2.0*double(i)/double(npts-1);
437 dynamic_cast<ELEMENT*
>(el_pt->bulk_element_pt())->
438 position(el_pt->local_coordinate_in_bulk(s),r_contact);
440 el_pt->contact_line_vectors(s,tangent,normal);
441 tangent_file << r_contact[0] <<
" "
442 << r_contact[1] <<
" "
443 << r_contact[2] <<
" "
446 << tangent[2] <<
" " << std::endl;
448 normal_file << r_contact[0] <<
" "
449 << r_contact[1] <<
" "
450 << r_contact[2] <<
" "
453 << normal[2] <<
" " << std::endl;
456 contact_angle_file << r_contact[1] <<
" "
457 << el_pt->actual_cos_contact_angle(s)
464 tangent_file.close();
466 contact_angle_file.close();
470 cout <<
"\n********************************************" << endl << endl;
478 void run_it(
const string& output_directory)
492 doc_info.set_directory(output_directory);
496 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
497 trace_file.open(filename);
501 <<
"VARIABLES=\"<GREEK>k</GREEK>\",\"<GREEK>k</GREEK>_{ex}\",\"h\""
503 trace_file <<
"ZONE" << std::endl;
527 problem.newton_solve();
549 int main(
int argc,
char* argv[])
553 CommandLineArgs::setup(argc,argv);
561 if (CommandLineArgs::Argc==1)
564 <<
"Running every case with limited number of steps for validation"
573 case_lo=atoi(argv[1]);
574 case_hi=atoi(argv[1]);
583 for (
unsigned my_case=case_lo;my_case<=case_hi;my_case++)
593 <<
"//////////////////////////////////////////////////////////\n"
594 <<
"All pinned solution \n"
595 <<
"//////////////////////////////////////////////////////////\n\n";
601 run_it(
"RESLT_all_pinned");
608 <<
"//////////////////////////////////////////////////////////\n"
609 <<
"Barrel-shaped solution \n"
610 <<
"/////////////////////////////////////////////////////////\n\n";
616 run_it(
"RESLT_barrel_shape");
623 <<
"//////////////////////////////////////////////////////////\n"
624 <<
"T-junction solution \n"
625 <<
"//////////////////////////////////////////////////////////\n\n";
636 run_it(
"RESLT_T_junction");
642 std::cout <<
"Wrong case! Options are:\n"
2D YoungLaplace problem on rectangular domain, discretised with 2D QYoungLaplace elements....
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...
unsigned Last_element_on_boundary1
Number of last FaceElement on boundary 1.
void create_contact_angle_elements(const unsigned &b)
Create YoungLaplace contact angle elements on the b-th boundary of the problem's mesh and add them to...
unsigned N_bulk_elements
Number of YoungLaplace "bulk" elements (We're attaching the contact angle elements to the bulk mesh...
YoungLaplaceProblem()
Constructor:
unsigned Last_element_on_boundary3
Number of last FaceElement on boundary 3.
void actions_after_newton_solve()
Update the problem after solve: Empty.
void actions_before_newton_solve()
Update the problem specs before solve.
~YoungLaplaceProblem()
Destructor (empty)
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 YoungLaplace 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.