32 #include "young_laplace.h"
35 #include "meshes/rectangular_quadmesh.h"
39 using namespace oomph;
50 double Cos_gamma=cos(MathematicalConstants::Pi/6.0);
70 Vector<double>& spine_B,
71 Vector< Vector<double> >& dspine_B)
78 dspine_B[0][0] = 1.0 ;
79 dspine_B[1][0] = 0.0 ;
80 dspine_B[0][1] = 0.0 ;
81 dspine_B[1][1] = 1.0 ;
82 dspine_B[0][2] = 0.0 ;
83 dspine_B[1][2] = 0.0 ;
93 double Alpha_min = MathematicalConstants::Pi/2.0*1.5;
96 double Alpha_max = MathematicalConstants::Pi/2.0*0.5;
102 Vector<double>& spine,
103 Vector< Vector<double> >& dspine)
137 template<
class ELEMENT>
159 if (Contact_angle_mesh_pt!=0) delete_contact_angle_elements();
162 rebuild_global_mesh();
169 create_contact_angle_elements(1);
170 create_contact_angle_elements(3);
173 unsigned nel=Contact_angle_mesh_pt->nelement();
174 for (
unsigned e=0;e<nel;e++)
178 YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
179 dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*
>(
180 Contact_angle_mesh_pt->element_pt(e));
187 rebuild_global_mesh();
193 void doc_solution(DocInfo& doc_info, ofstream& trace_file);
199 void create_contact_angle_elements(
const unsigned& b);
202 void delete_contact_angle_elements();
225 template<
class ELEMENT>
245 Bulk_mesh_pt=
new RefineableRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
248 Bulk_mesh_pt->spatial_error_estimator_pt()=
new Z2ErrorEstimator;
251 Bulk_mesh_pt->max_permitted_error()=1.0e-4;
252 Bulk_mesh_pt->min_permitted_error()=1.0e-6;
256 if ((n_x%2!=0)||(n_y%2!=0))
258 cout <<
"n_x n_y should be even" << endl;
263 ELEMENT* prescribed_height_element_pt=
dynamic_cast<ELEMENT*
>(
264 Bulk_mesh_pt->element_pt(n_y*n_x/2+n_x/2));
267 Control_node_pt=
static_cast<Node*
>(prescribed_height_element_pt->node_pt(0));
269 std::cout <<
"Controlling height at (x,y) : (" << Control_node_pt->x(0)
270 <<
"," << Control_node_pt->x(1) <<
")" <<
"\n" << endl;
274 HeightControlElement* height_control_element_pt=
new HeightControlElement(
278 Height_control_mesh_pt =
new Mesh;
279 Height_control_mesh_pt->add_element_pt(height_control_element_pt);
282 Kappa_pt=height_control_element_pt->kappa_pt();
292 Contact_angle_mesh_pt=
new Mesh;
295 create_contact_angle_elements(1);
296 create_contact_angle_elements(3);
301 add_sub_mesh(Bulk_mesh_pt);
302 add_sub_mesh(Height_control_mesh_pt);
303 add_sub_mesh(Contact_angle_mesh_pt);
313 unsigned n_bound = Bulk_mesh_pt->nboundary();
314 for(
unsigned b=0;b<n_bound;b++)
320 unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
321 for (
unsigned n=0;n<n_node;n++)
323 Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
332 unsigned n_bulk=Bulk_mesh_pt->nelement();
333 for(
unsigned i=0;i<n_bulk;i++)
336 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(Bulk_mesh_pt->element_pt(i));
347 unsigned nel=Contact_angle_mesh_pt->nelement();
348 for (
unsigned e=0;e<nel;e++)
352 YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
353 dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*
>(
354 Contact_angle_mesh_pt->element_pt(e));
362 cout <<
"\nNumber of equations: " << assign_eqn_numbers() << endl;
363 cout <<
"\n********************************************\n" << endl;
372 template<
class ELEMENT>
377 unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
380 for(
unsigned e=0;e<n_element;e++)
383 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
384 Bulk_mesh_pt->boundary_element_pt(b,e));
387 int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
390 YoungLaplaceContactAngleElement<ELEMENT>* contact_angle_element_pt =
new
391 YoungLaplaceContactAngleElement<ELEMENT>(bulk_elem_pt,face_index);
394 Contact_angle_mesh_pt->add_element_pt(contact_angle_element_pt);
404 template<
class ELEMENT>
409 unsigned n_element = Contact_angle_mesh_pt->nelement();
412 for(
unsigned e=0;e<n_element;e++)
415 delete Contact_angle_mesh_pt->element_pt(e);
419 Contact_angle_mesh_pt->flush_element_and_node_storage();
429 template<
class ELEMENT>
431 ofstream& trace_file)
436 trace_file << -1.0*
Kappa_pt->value(0) <<
" ";
437 trace_file << Control_node_pt->value(0) ;
448 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
450 some_file.open(filename);
451 Bulk_mesh_pt->output(some_file,npts);
457 ofstream tangent_file;
458 sprintf(filename,
"%s/tangent_to_contact_line%i.dat",
459 doc_info.directory().c_str(),
461 tangent_file.open(filename);
463 ofstream normal_file;
464 sprintf(filename,
"%s/normal_to_contact_line%i.dat",
465 doc_info.directory().c_str(),
467 normal_file.open(filename);
470 ofstream contact_angle_file;
471 sprintf(filename,
"%s/contact_angle%i.dat",
472 doc_info.directory().c_str(),
474 contact_angle_file.open(filename);
477 Vector<double> tangent(3);
478 Vector<double> normal(3);
479 Vector<double> r_contact(3);
482 unsigned n_element = Contact_angle_mesh_pt->nelement();
485 for(
unsigned e=0;e<n_element;e++)
488 tangent_file <<
"ZONE" << std::endl;
489 normal_file <<
"ZONE" << std::endl;
490 contact_angle_file <<
"ZONE" << std::endl;
493 YoungLaplaceContactAngleElement<ELEMENT>* el_pt =
494 dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*
>(
495 Contact_angle_mesh_pt->element_pt(e));
499 for (
unsigned i=0;i<npts;i++)
501 s[0]=-1.0+2.0*double(i)/double(npts-1);
503 dynamic_cast<ELEMENT*
>(el_pt->bulk_element_pt())->
504 position(el_pt->local_coordinate_in_bulk(s),r_contact);
506 el_pt->contact_line_vectors(s,tangent,normal);
507 tangent_file << r_contact[0] <<
" "
508 << r_contact[1] <<
" "
509 << r_contact[2] <<
" "
512 << tangent[2] <<
" " << std::endl;
514 normal_file << r_contact[0] <<
" "
515 << r_contact[1] <<
" "
516 << r_contact[2] <<
" "
519 << normal[2] <<
" " << std::endl;
521 contact_angle_file << r_contact[1] <<
" "
522 << el_pt->actual_cos_contact_angle(s)
529 tangent_file.close();
531 contact_angle_file.close();
534 cout <<
"\n********************************************" << endl << endl;
553 doc_info.set_directory(
"RESLT");
557 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
558 trace_file.open(filename);
561 trace_file <<
"VARIABLES=\"<GREEK>k</GREEK>\",\"h\"" << std::endl;
562 trace_file <<
"ZONE" << std::endl;
573 problem.refine_uniformly();
583 double increment=0.1;
587 for (
unsigned istep=0;istep<nstep;istep++)
592 unsigned max_adapt=1;
593 problem.newton_solve(max_adapt);
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.
Data * Kappa_pt
Pointer to Data object that stores the prescribed curvature.
~RefineableYoungLaplaceProblem()
Destructor (empty)
RefineableRectangularQuadMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
void actions_before_newton_solve()
Update the problem specs before solve: Empty.
Mesh * Contact_angle_mesh_pt
Pointer to the contact angle mesh.
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.
Mesh * Height_control_mesh_pt
Pointer to mesh containing the height control element.
Node * Control_node_pt
Node at which the height (displacement along spine) is controlled/doced.
Namespace for "global" problem parameters.
double Alpha_max
Max. spine angle against horizontal plane.
double L_x
Length and width of the domain.
double Controlled_height
Height control value.
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...
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 Alpha_min
Min. spine angle against horizontal plane.
double Cos_gamma
Cos of contact angle.