young_laplace.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 //LIC//
8 //LIC// This library is free software; you can redistribute it and/or
9 //LIC// modify it under the terms of the GNU Lesser General Public
10 //LIC// License as published by the Free Software Foundation; either
11 //LIC// version 2.1 of the License, or (at your option) any later version.
12 //LIC//
13 //LIC// This library is distributed in the hope that it will be useful,
14 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 //LIC// Lesser General Public License for more details.
17 //LIC//
18 //LIC// You should have received a copy of the GNU Lesser General Public
19 //LIC// License along with this library; if not, write to the Free Software
20 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 //LIC// 02110-1301 USA.
22 //LIC//
23 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 //LIC//
25 //LIC//====================================================================
26 // Driver code for a 2D YoungLaplace problem
27 
28 // Generic routines
29 #include "generic.h"
30 
31 // The YoungLaplace equations
32 #include "young_laplace.h"
33 
34 // The mesh
35 #include "meshes/simple_rectangular_quadmesh.h"
36 
37 
38 using namespace std;
39 using namespace oomph;
40 
41 
42 // Namespace (shared with refineable version)
44 
45 
46 
47 //====== start_of_problem_class=======================================
48 /// 2D YoungLaplace problem on rectangular domain, discretised with
49 /// 2D QYoungLaplace elements. The specific type of element is
50 /// specified via the template parameter.
51 //====================================================================
52 template<class ELEMENT>
53 class YoungLaplaceProblem : public Problem
54 {
55 
56 public:
57 
58  /// Constructor:
60 
61  /// Destructor (empty)
63 
64  /// Update the problem specs before solve
66 
67  /// Update the problem after solve: Empty
69 
70  /// Doc the solution. DocInfo object stores flags/labels for where the
71  /// output gets written to and the trace file
72  void doc_solution(DocInfo& doc_info, ofstream& trace_file);
73 
74 private:
75 
76  /// Create YoungLaplace contact angle elements on the
77  /// b-th boundary of the problem's mesh and add them to mesh
78  void create_contact_angle_elements(const unsigned& b);
79 
80  /// Number of YoungLaplace "bulk" elements (We're attaching the
81  /// contact angle elements to the bulk mesh --> only the first
82  /// N_bulk_elements elements in the mesh are bulk elements!)
83  unsigned N_bulk_elements;
84 
85  /// Number of last FaceElement on boundary 1
87 
88  /// Number of last FaceElement on boundary 3
90 
91  /// Node at which the height (displacement along spine) is controlled/doced
92  Node* Control_node_pt;
93 
94 }; // end of problem class
95 
96 
97 //=====start_of_constructor===============================================
98 /// Constructor for YoungLaplace problem
99 //========================================================================
100 template<class ELEMENT>
102 {
103 
104  // Setup dependent parameters in namespace
106 
107  // Setup mesh
108  //-----------
109 
110  // # of elements in x-direction
111  unsigned n_x=GlobalParameters::N_x;
112 
113  // # of elements in y-direction
114  unsigned n_y=GlobalParameters::N_y;
115 
116  // Domain length in x-direction
117  double l_x=GlobalParameters::L_x;
118 
119  // Domain length in y-direction
120  double l_y=GlobalParameters::L_y;
121 
122  // Print Size of the mesh
123  cout << "Lx = " << l_x << " and Ly = " << l_y << endl;
124 
125  // Build and assign mesh
126  Problem::mesh_pt()=new SimpleRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
127 
128 
129  // Choose the prescribed height element
130  ELEMENT* prescribed_height_element_pt= dynamic_cast<ELEMENT*>(
131  mesh_pt()->element_pt(GlobalParameters::Control_element));
132 
133  // ...and the associated control node (node 0 in that element)
134  // (we're storing this node even if there's no height-control, for
135  // output purposes...)
136  Control_node_pt= static_cast<Node*>(
137  prescribed_height_element_pt->node_pt(0));
138 
139  // If needed, create a height control element and store the
140  // pointer to the Kappa Data created by this object
141  HeightControlElement* height_control_element_pt=0;
143  {
144  cout << "Controlling height at (x,y) : (" << Control_node_pt->x(0)
145  << "," << Control_node_pt->x(1) << ")" << "\n" << endl;
146 
147  height_control_element_pt=new HeightControlElement(
148  Control_node_pt,&GlobalParameters::Controlled_height);
149 
150  GlobalParameters::Kappa_pt=height_control_element_pt->kappa_pt();
151 
152  // Assign initial value for kappa value
153  height_control_element_pt->kappa_pt()
154  ->set_value(0,GlobalParameters::Kappa_initial);
155  }
156  //...otherwise create a kappa data item from scratch and pin its
157  // single unknown value
158  else
159  {
160  GlobalParameters::Kappa_pt=new Data(1);
163  }
164 
165  // Number of elements in the bulk mesh
166  N_bulk_elements = mesh_pt()->nelement();
167 
168  // Create contact angle elements from all elements that are
169  // adjacent to boundary 1 and 3 and add them to the (single) global mesh
172  {
173  create_contact_angle_elements(1);
174  Last_element_on_boundary1=mesh_pt()->nelement();
175 
176  create_contact_angle_elements(3);
177  Last_element_on_boundary3=mesh_pt()->nelement();
178  }
179  else
180  {
181  Last_element_on_boundary1=0;
182  Last_element_on_boundary3=0;
183  }
184 
185  // Boundary conditions
186  //--------------------
187 
188  // Set the boundary conditions for this problem: All nodes are
189  // free by default -- only need to pin the ones that have Dirichlet conditions
190  // here.
191  unsigned n_bound = mesh_pt()->nboundary();
192  for(unsigned b=0;b<n_bound;b++)
193  {
194  // Pin all boundaries for two cases and only boundaries
195  // 0 and 2 in all others:
197  (b==0)||
198  (b==2))
199  {
200  unsigned n_node = mesh_pt()->nboundary_node(b);
201  for (unsigned n=0;n<n_node;n++)
202  {
203  mesh_pt()->boundary_node_pt(b,n)->pin(0);
204  }
205  }
206  }
207 
208  // Complete build of elements
209  //---------------------------
210 
211  // Complete the build of all elements so they are fully functional
212  for(unsigned i=0;i<N_bulk_elements;i++)
213  {
214  // Upcast from GeneralsedElement to the present element
215  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
216 
218  {
219  //Set the spine function pointers
220  el_pt->spine_base_fct_pt() = GlobalParameters::spine_base_function;
221  el_pt->spine_fct_pt() = GlobalParameters::spine_function;
222  }
223 
224  // Set the curvature data for the element
225  el_pt->set_kappa(GlobalParameters::Kappa_pt);
226 
227  }
228 
229  // Set function pointers for contact-angle elements
232  {
233  // Loop over the contact-angle elements (located at the "end" of the
234  // mesh) to pass function pointer to contact angle function.
235  for (unsigned e=N_bulk_elements;e<Last_element_on_boundary1;e++)
236  {
237  // Upcast from GeneralisedElement to YoungLaplace contact angle element
238  YoungLaplaceContactAngleElement<ELEMENT>* el_pt =
239  dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*>(
240  mesh_pt()->element_pt(e));
241 
242  // Set the pointer to the prescribed contact angle
243  el_pt->prescribed_cos_gamma_pt() = &GlobalParameters::Cos_gamma;
244  }
245  for (unsigned e=Last_element_on_boundary1;e<Last_element_on_boundary3;e++)
246  {
247  // Upcast from GeneralisedElement to YoungLaplace contact angle element
248  YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
249  dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*>(
250  mesh_pt()->element_pt(e));
251 
252  // Set the pointer to the prescribed contact angle
253  el_pt->prescribed_cos_gamma_pt() = &GlobalParameters::Cos_gamma;
254  }
255  }
256 
257  // Height Element setup
258  //---------------------
259 
260  /// Add height control element to mesh at the very end
262  {
263  mesh_pt()->add_element_pt(height_control_element_pt);
264  }
265 
266 
267  // Setup equation numbering scheme
268  cout <<"\nNumber of equations: " << assign_eqn_numbers() << endl;
269  cout << "\n********************************************\n" << endl;
270 
271 } // end of constructor
272 
273 
274 
275 //============start_of_create_contact_angle_elements=====================
276 /// Create YoungLaplace contact angle elements on the b-th boundary of the Mesh
277 //=======================================================================
278 template<class ELEMENT>
280  const unsigned &b)
281 {
282  // How many bulk elements are adjacent to boundary b?
283  unsigned n_element = mesh_pt()->nboundary_element(b);
284 
285  // Loop over the bulk elements adjacent to boundary b?
286  for(unsigned e=0;e<n_element;e++)
287  {
288  // Get pointer to the bulk element that is adjacent to boundary b
289  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
290  mesh_pt()->boundary_element_pt(b,e));
291 
292  // What is the index of the face of the bulk element at the boundary
293  int face_index = mesh_pt()->face_index_at_boundary(b,e);
294 
295  // Build the corresponding prescribed contact angle element
296  YoungLaplaceContactAngleElement<ELEMENT>* contact_angle_element_pt = new
297  YoungLaplaceContactAngleElement<ELEMENT>(bulk_elem_pt,face_index);
298 
299  //Add the prescribed contact angle element to the mesh
300  mesh_pt()->add_element_pt(contact_angle_element_pt);
301 
302  } //end of loop over bulk elements adjacent to boundary b
303 
304 } // end of create_contact_angle_elements
305 
306 
307 
308 
309 //========================================start_of_actions_before_solve===
310 /// Update the problem specs before solve: (Re-)set boundary conditions
311 /// to the values from the exact solution.
312 //========================================================================
313 template<class ELEMENT>
315 {
316 
317  // Increment kappa or height value
319  {
320  double kappa=GlobalParameters::Kappa_pt->value(0);
322  GlobalParameters::Kappa_pt->set_value(0,kappa);
323 
324  cout << "Solving for Prescribed KAPPA Value = " ;
325  cout << GlobalParameters::Kappa_pt->value(0) << "\n" << endl;
326  }
327  else
328  {
331 
332  cout << "Solving for Prescribed HEIGHT Value = " ;
333  cout << GlobalParameters::Controlled_height << "\n" << endl;
334  }
335 
336 } // end of actions before solve
337 
338 
339 
340 
341 //===============start_of_doc=============================================
342 /// Doc the solution: doc_info contains labels/output directory etc.
343 //========================================================================
344 template<class ELEMENT>
345 void YoungLaplaceProblem<ELEMENT>::doc_solution(DocInfo& doc_info,
346  ofstream& trace_file)
347 {
348 
349  // Output kappa vs height of the apex
350  //------------------------------------
351  trace_file << -1.0*GlobalParameters::Kappa_pt->value(0) << " ";
352  trace_file << GlobalParameters::get_exact_kappa() << " ";
353  trace_file << Control_node_pt->value(0) ;
354  trace_file << endl;
355 
356  // Number of plot points: npts x npts
357  unsigned npts=5;
358 
359  // Output full solution
360  //---------------------
361  ofstream some_file;
362  char filename[100];
363  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
364  doc_info.number());
365  some_file.open(filename);
366  for(unsigned i=0;i<N_bulk_elements;i++)
367  {
368  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
369  el_pt->output(some_file,npts);
370  }
371  some_file.close();
372 
373  // Output contact angle
374  //---------------------
375 
376  //Doc contact angle stuff
379  {
380 
381  ofstream tangent_file;
382  sprintf(filename,"%s/tangent_to_contact_line%i.dat",
383  doc_info.directory().c_str(),
384  doc_info.number());
385  tangent_file.open(filename);
386 
387  ofstream normal_file;
388  sprintf(filename,"%s/normal_to_contact_line%i.dat",
389  doc_info.directory().c_str(),
390  doc_info.number());
391  normal_file.open(filename);
392 
393 
394  ofstream contact_angle_file;
395  sprintf(filename,"%s/contact_angle%i.dat",
396  doc_info.directory().c_str(),
397  doc_info.number());
398  contact_angle_file.open(filename);
399 
400  // Tangent and normal vectors to contact line
401  Vector<double> tangent(3);
402  Vector<double> normal(3);
403  Vector<double> r_contact(3);
404 
405 
406  // Loop over both boundaries
407  for (unsigned bnd=0;bnd<2;bnd++)
408  {
409  unsigned el_lo=N_bulk_elements;
410  unsigned el_hi=Last_element_on_boundary1;
411  if (1==bnd)
412  {
413  el_lo=Last_element_on_boundary1;
414  el_hi=Last_element_on_boundary3;
415  }
416 
417  // Loop over the contact-angle elements (located at the "end" of the
418  // mesh)
419  for (unsigned e=el_lo;e<el_hi;e++)
420  {
421 
422  tangent_file << "ZONE" << std::endl;
423  normal_file << "ZONE" << std::endl;
424  contact_angle_file << "ZONE" << std::endl;
425 
426  // Upcast from GeneralisedElement to YoungLaplace contact angle element
427  YoungLaplaceContactAngleElement<ELEMENT>* el_pt =
428  dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*>(
429  mesh_pt()->element_pt(e));
430 
431  // Loop over a few points in the contact angle element
432  Vector<double> s(1);
433  for (unsigned i=0;i<npts;i++)
434  {
435  s[0]=-1.0+2.0*double(i)/double(npts-1);
436 
437  dynamic_cast<ELEMENT*>(el_pt->bulk_element_pt())->
438  position(el_pt->local_coordinate_in_bulk(s),r_contact);
439 
440  el_pt->contact_line_vectors(s,tangent,normal);
441  tangent_file << r_contact[0] << " "
442  << r_contact[1] << " "
443  << r_contact[2] << " "
444  << tangent[0] << " "
445  << tangent[1] << " "
446  << tangent[2] << " " << std::endl;
447 
448  normal_file << r_contact[0] << " "
449  << r_contact[1] << " "
450  << r_contact[2] << " "
451  << normal[0] << " "
452  << normal[1] << " "
453  << normal[2] << " " << std::endl;
454 
455 
456  contact_angle_file << r_contact[1] << " "
457  << el_pt->actual_cos_contact_angle(s)
458  << std::endl;
459  }
460  }
461 
462  } // end of loop over both boundaries
463 
464  tangent_file.close();
465  normal_file.close();
466  contact_angle_file.close();
467 
468  }
469 
470 cout << "\n********************************************" << endl << endl;
471 
472 } // end of doc
473 
474 //========================================================================
475 /// Run code for current setting of parameter values -- specify name
476 /// of output directory
477 //========================================================================
478 void run_it(const string& output_directory)
479 {
480 
481  // Create label for output
482  //------------------------
483  DocInfo doc_info;
484 
485  // Set outputs
486  //------------
487 
488  // Trace file
489  ofstream trace_file;
490 
491  // Set output directory
492  doc_info.set_directory(output_directory);
493 
494  //Open a trace file
495  char filename[100];
496  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
497  trace_file.open(filename);
498 
499  // Write kappa, exact kappa if exists and height values
500  trace_file
501  << "VARIABLES=\"<GREEK>k</GREEK>\",\"<GREEK>k</GREEK>_{ex}\",\"h\""
502  << std::endl;
503  trace_file << "ZONE" << std::endl;
504 
505  //Set up the problem
506  //------------------
507 
508  // Create the problem with 2D nine-node elements from the
509  // QYoungLaplaceElement family.
511 
512  //Output the solution
513  problem.doc_solution(doc_info,trace_file);
514 
515  //Increment counter for solutions
516  doc_info.number()++;
517 
518 
519  // Parameter incrementation
520  //-------------------------
521 
522  // Loop over steps
523  for (unsigned istep=0;istep<GlobalParameters::Nsteps;istep++)
524  {
525 
526  // Solve the problem
527  problem.newton_solve();
528 
529  //Output the solution
530  problem.doc_solution(doc_info,trace_file);
531 
532  //Increment counter for solutions
533  doc_info.number()++;
534 
535  }
536 
537  // Close output file
538  trace_file.close();
539 
540 } //end of run_it()
541 
542 
543 
544 //===== start_of_main=====================================================
545 /// Driver code for 2D YoungLaplace problem. Input arguments: none
546 /// (for validation) or case (0,1,2 for all pinned, barrel and T junction)
547 /// and number of steps.
548 //========================================================================
549 int main(int argc, char* argv[])
550 {
551 
552  // Store command line arguments
553  CommandLineArgs::setup(argc,argv);
554 
555  // Cases to run (By default (validation) run all
556  unsigned case_lo=0;
557  unsigned case_hi=2;
558 
559  // No command line args: Running every case with
560  // limited number of steps
561  if (CommandLineArgs::Argc==1)
562  {
563  std::cout
564  << "Running every case with limited number of steps for validation"
565  << std::endl;
566 
567  // Number of steps
569  }
570  else
571  {
572  // Which case to run?
573  case_lo=atoi(argv[1]);
574  case_hi=atoi(argv[1]);
575 
576  // Number of steps
577  GlobalParameters::Nsteps=atoi(argv[2]);
578  }
579 
580 
581  // Loop over chosen case(s)
582  //-------------------------
583  for (unsigned my_case=case_lo;my_case<=case_hi;my_case++)
584  {
585 
586  // Choose
587  switch (my_case)
588  {
589 
590  case 0:
591 
592  cout << endl << endl
593  << "//////////////////////////////////////////////////////////\n"
594  << "All pinned solution \n"
595  << "//////////////////////////////////////////////////////////\n\n";
596 
598 
599  // Run with spines
601  run_it("RESLT_all_pinned");
602 
603  break;
604 
605  case 1:
606 
607  cout << endl << endl
608  << "//////////////////////////////////////////////////////////\n"
609  << "Barrel-shaped solution \n"
610  << "/////////////////////////////////////////////////////////\n\n";
611 
613 
614  // Run with spines
616  run_it("RESLT_barrel_shape");
617 
618  break;
619 
620  case 2:
621 
622  cout << endl << endl
623  << "//////////////////////////////////////////////////////////\n"
624  << "T-junction solution \n"
625  << "//////////////////////////////////////////////////////////\n\n";
626 
629 
630  // Adjust spine orientation
632  GlobalParameters::Gamma=MathematicalConstants::Pi/6.0;
633 
634  // Run with spines
636  run_it("RESLT_T_junction");
637 
638  break;
639 
640  default:
641 
642  std::cout << "Wrong case! Options are:\n"
643  << "0: All pinned\n"
644  << "1: Barrel \n"
645  << "2: T_junction\n"
646  << std::endl;
647  assert(false);
648 
649  }
650 
651  }
652 
653 
654 } //end of main
655 
656 
657 
658 
659 
660 
2D YoungLaplace problem on rectangular domain, discretised with 2D QYoungLaplace elements....
Definition: barrel.cc:140
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.
Definition: barrel.cc:51
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.
Definition: barrel.cc:54
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...
Definition: barrel.cc:104
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...
Definition: barrel.cc:72
double L_y
Width of domain.
double Gamma
Contact angle and its cos (dependent parameter – is reassigned)
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.