refineable_t_junction.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 for refineable Young Laplace 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/rectangular_quadmesh.h"
36 
37 // Namespaces
38 using namespace std;
39 using namespace oomph;
40 
41 
42 
43 //======start_of_namespace========================================
44 /// Namespace for "global" problem parameters
45 //================================================================
46 namespace GlobalParameters
47 {
48 
49  /// Cos of contact angle
50  double Cos_gamma=cos(MathematicalConstants::Pi/6.0);
51 
52  /// Height control value for displacement control
53  double Controlled_height = 0.0;
54 
55  /// Length of domain
56  double L_x = 1.0;
57 
58  /// Width of domain
59  double L_y = 5.0;
60 
61  // Spine basis
62  //------------
63 
64  /// Spine basis: The position vector to the basis of the spine
65  /// as a function of the two coordinates x_1 and x_2, and its
66  /// derivatives w.r.t. to these coordinates.
67  /// dspine_B[i][j] = d spine_B[j] / dx_i
68  /// Spines start in the (x_1,x_2) plane at (x_1,x_2).
69  void spine_base_function(const Vector<double>& x,
70  Vector<double>& spine_B,
71  Vector< Vector<double> >& dspine_B)
72  {
73 
74  // Bspines and derivatives
75  spine_B[0] = x[0];
76  spine_B[1] = x[1];
77  spine_B[2] = 0.0 ;
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 ;
84 
85  } // End of bspine functions
86 
87 
88 
89  // Spines rotate in the y-direction
90  //---------------------------------
91 
92  /// Min. spine angle against horizontal plane
93  double Alpha_min = MathematicalConstants::Pi/2.0*1.5;
94 
95  /// Max. spine angle against horizontal plane
96  double Alpha_max = MathematicalConstants::Pi/2.0*0.5;
97 
98  /// Spine: The spine vector field as a function of the two
99  /// coordinates x_1 and x_2, and its derivatives w.r.t. to these coordinates:
100  /// dspine[i][j] = d spine[j] / dx_i
101  void spine_function(const Vector<double>& x,
102  Vector<double>& spine,
103  Vector< Vector<double> >& dspine)
104  {
105 
106  /// Spines (and derivatives) are independent of x[0] and rotate
107  /// in the x[1]-direction
108  spine[0]=0.0;
109  dspine[0][0]=0.0;
110  dspine[1][0]=0.0;
111 
112  spine[1]=cos(Alpha_min+(Alpha_max-Alpha_min)*x[1]/L_y);
113  dspine[0][1]=0.0;
114  dspine[1][1]=-sin(Alpha_min+(Alpha_max-Alpha_min)*x[1]/L_y)
115  *(Alpha_max-Alpha_min)/L_y;
116 
117  spine[2]=sin(Alpha_min+(Alpha_max-Alpha_min)*x[1]/L_y);
118  dspine[0][2]=0.0;
119  dspine[1][2]=cos(Alpha_min+(Alpha_max-Alpha_min)*x[1]/L_y)
120  *(Alpha_max-Alpha_min)/L_y;
121 
122  } // End spine function
123 
124 
125 } // end of namespace
126 
127 
128 
129 
130 
131 
132 //====== start_of_problem_class=======================================
133 /// 2D RefineableYoungLaplace problem on rectangular domain, discretised with
134 /// 2D QRefineableYoungLaplace elements. The specific type of element is
135 /// specified via the template parameter.
136 //====================================================================
137 template<class ELEMENT>
138 class RefineableYoungLaplaceProblem : public Problem
139 {
140 
141 public:
142 
143  /// Constructor:
145 
146  /// Destructor (empty)
148 
149  /// Update the problem specs before solve: Empty
151 
152  /// Update the problem after solve: Empty
154 
155  /// Actions before adapt: Wipe the mesh of contact angle elements
157  {
158  // Kill the contact angle elements and wipe contact angle mesh
159  if (Contact_angle_mesh_pt!=0) delete_contact_angle_elements();
160 
161  // Rebuild the Problem's global mesh from its various sub-meshes
162  rebuild_global_mesh();
163  }
164 
165 
166  /// Actions after adapt: Rebuild the mesh of contact angle elements
168  {
169  create_contact_angle_elements(1);
170  create_contact_angle_elements(3);
171 
172  // Set function pointers for contact-angle elements
173  unsigned nel=Contact_angle_mesh_pt->nelement();
174  for (unsigned e=0;e<nel;e++)
175  {
176  // Upcast from GeneralisedElement to YoungLaplace contact angle
177  // element
178  YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
179  dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*>(
180  Contact_angle_mesh_pt->element_pt(e));
181 
182  // Set the pointer to the prescribed contact angle
183  el_pt->prescribed_cos_gamma_pt() = &GlobalParameters::Cos_gamma;
184  }
185 
186  // Rebuild the Problem's global mesh from its various sub-meshes
187  rebuild_global_mesh();
188 
189  }
190 
191  /// Doc the solution. DocInfo object stores flags/labels for where the
192  /// output gets written to and the trace file
193  void doc_solution(DocInfo& doc_info, ofstream& trace_file);
194 
195 private:
196 
197  /// Create YoungLaplace contact angle elements on the
198  /// b-th boundary of the bulk mesh and add them to contact angle mesh
199  void create_contact_angle_elements(const unsigned& b);
200 
201  /// Delete contact angle elements
202  void delete_contact_angle_elements();
203 
204  /// Pointer to the "bulk" mesh
205  RefineableRectangularQuadMesh<ELEMENT>* Bulk_mesh_pt;
206 
207  /// Pointer to the contact angle mesh
209 
210  /// Pointer to mesh containing the height control element
212 
213  /// Node at which the height (displacement along spine) is controlled/doced
215 
216  /// Pointer to Data object that stores the prescribed curvature
217  Data* Kappa_pt;
218 
219 }; // end of problem class
220 
221 
222 //=====start_of_constructor===============================================
223 /// Constructor for RefineableYoungLaplace problem
224 //========================================================================
225 template<class ELEMENT>
227 {
228 
229  // Setup bulk mesh
230  //----------------
231 
232  // # of elements in x-direction
233  unsigned n_x=8;
234 
235  // # of elements in y-direction
236  unsigned n_y=8;
237 
238  // Domain length in x-direction
239  double l_x=GlobalParameters::L_x;
240 
241  // Domain length in y-direction
242  double l_y=GlobalParameters::L_y;
243 
244  // Build and assign mesh
245  Bulk_mesh_pt=new RefineableRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
246 
247  // Create/set error estimator
248  Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
249 
250  // Set targets for spatial adaptivity
251  Bulk_mesh_pt->max_permitted_error()=1.0e-4;
252  Bulk_mesh_pt->min_permitted_error()=1.0e-6;
253 
254  // Check that we've got an even number of elements otherwise
255  // out counting doesn't work...
256  if ((n_x%2!=0)||(n_y%2!=0))
257  {
258  cout << "n_x n_y should be even" << endl;
259  abort();
260  }
261 
262  // This is the element that contains the central node:
263  ELEMENT* prescribed_height_element_pt= dynamic_cast<ELEMENT*>(
264  Bulk_mesh_pt->element_pt(n_y*n_x/2+n_x/2));
265 
266  // The central node is node 0 in that element
267  Control_node_pt= static_cast<Node*>(prescribed_height_element_pt->node_pt(0));
268 
269  std::cout << "Controlling height at (x,y) : (" << Control_node_pt->x(0)
270  << "," << Control_node_pt->x(1) << ")" << "\n" << endl;
271 
272  // Create a height control element and store the
273  // pointer to the Kappa Data created by this object
274  HeightControlElement* height_control_element_pt=new HeightControlElement(
275  Control_node_pt,&GlobalParameters::Controlled_height);
276 
277  // Add to mesh
278  Height_control_mesh_pt = new Mesh;
279  Height_control_mesh_pt->add_element_pt(height_control_element_pt);
280 
281  // Store curvature data
282  Kappa_pt=height_control_element_pt->kappa_pt();
283 
284 
285  // Contact angle elements
286  //-----------------------
287 
288  // Create prescribed-contact-angle elements from all elements that are
289  // adjacent to boundary 1 and 3 and add them to their own mesh
290 
291  // set up new mesh
292  Contact_angle_mesh_pt=new Mesh;
293 
294  // creation of contact angle elements
295  create_contact_angle_elements(1);
296  create_contact_angle_elements(3);
297 
298 
299  // Add various meshes and build the global mesh
300  //----------------------------------------------
301  add_sub_mesh(Bulk_mesh_pt);
302  add_sub_mesh(Height_control_mesh_pt);
303  add_sub_mesh(Contact_angle_mesh_pt);
304  build_global_mesh();
305 
306 
307  // Boundary conditions
308  //--------------------
309 
310  // Set the boundary conditions for this problem: All nodes are
311  // free by default -- only need to pin the ones that have Dirichlet conditions
312  // here.
313  unsigned n_bound = Bulk_mesh_pt->nboundary();
314  for(unsigned b=0;b<n_bound;b++)
315  {
316  // Pin all boundaries for three cases and only boundaries
317  // 0 and 2 in all others:
318  if ((b==0)||(b==2))
319  {
320  unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
321  for (unsigned n=0;n<n_node;n++)
322  {
323  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
324  }
325  }
326  } // end bcs
327 
328  // Complete build of elements
329  //---------------------------
330 
331  // Complete the build of all elements so they are fully functional
332  unsigned n_bulk=Bulk_mesh_pt->nelement();
333  for(unsigned i=0;i<n_bulk;i++)
334  {
335  // Upcast from GeneralsedElement to the present element
336  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(i));
337 
338  //Set the spine function pointers
339  el_pt->spine_base_fct_pt() = GlobalParameters::spine_base_function;
340  el_pt->spine_fct_pt() = GlobalParameters::spine_function;
341 
342  // Set the curvature data for the element
343  el_pt->set_kappa(Kappa_pt);
344  }
345 
346  // Set function pointers for contact-angle elements
347  unsigned nel=Contact_angle_mesh_pt->nelement();
348  for (unsigned e=0;e<nel;e++)
349  {
350  // Upcast from GeneralisedElement to YoungLaplace contact angle
351  // element
352  YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
353  dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*>(
354  Contact_angle_mesh_pt->element_pt(e));
355 
356  // Set the pointer to the prescribed contact angle
357  el_pt->prescribed_cos_gamma_pt() = &GlobalParameters::Cos_gamma;
358  }
359 
360 
361  // Setup equation numbering scheme
362  cout <<"\nNumber of equations: " << assign_eqn_numbers() << endl;
363  cout << "\n********************************************\n" << endl;
364 
365 } // end of constructor
366 
367 
368 //============start_of_create_contact_angle_elements=====================
369 /// Create YoungLaplace contact angle elements on the b-th boundary of the
370 /// bulk mesh and add them to the contact angle mesh
371 //=======================================================================
372 template<class ELEMENT>
374  const unsigned &b)
375 {
376  // How many bulk elements are adjacent to boundary b?
377  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
378 
379  // Loop over the bulk elements adjacent to boundary b?
380  for(unsigned e=0;e<n_element;e++)
381  {
382  // Get pointer to the bulk element that is adjacent to boundary b
383  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
384  Bulk_mesh_pt->boundary_element_pt(b,e));
385 
386  // What is the index of the face of the bulk element at the boundary
387  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
388 
389  // Build the corresponding contact angle element
390  YoungLaplaceContactAngleElement<ELEMENT>* contact_angle_element_pt = new
391  YoungLaplaceContactAngleElement<ELEMENT>(bulk_elem_pt,face_index);
392 
393  //Add the contact angle element to the contact angle mesh
394  Contact_angle_mesh_pt->add_element_pt(contact_angle_element_pt);
395 
396  } //end of loop over bulk elements adjacent to boundary b
397 
398 } // end of create_contact_angle_elements
399 
400 
401 //============start_of_delete_contact_angle_elements=====================
402 /// Delete YoungLaplace contact angle elements
403 //=======================================================================
404 template<class ELEMENT>
406 {
407 
408  // How many contact angle elements are there?
409  unsigned n_element = Contact_angle_mesh_pt->nelement();
410 
411  // Loop over the surface elements
412  for(unsigned e=0;e<n_element;e++)
413  {
414  // Kill surface element
415  delete Contact_angle_mesh_pt->element_pt(e);
416  }
417 
418  // Wipe the mesh
419  Contact_angle_mesh_pt->flush_element_and_node_storage();
420 
421 
422 } // end of delete_contact_angle_elements
423 
424 
425 
426 //===============start_of_doc=============================================
427 /// Doc the solution: doc_info contains labels/output directory etc.
428 //========================================================================
429 template<class ELEMENT>
431  ofstream& trace_file)
432 {
433 
434  // Output kappa vs height
435  //-----------------------
436  trace_file << -1.0*Kappa_pt->value(0) << " ";
437  trace_file << Control_node_pt->value(0) ;
438  trace_file << endl;
439 
440  // Number of plot points: npts x npts
441  unsigned npts=5;
442 
443  // Output full solution
444  //---------------------
445  ofstream some_file;
446  char filename[100];
447  //YoungLaplaceEquations::Output_meniscus_and_spines=false;
448  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
449  doc_info.number());
450  some_file.open(filename);
451  Bulk_mesh_pt->output(some_file,npts);
452  some_file.close();
453 
454  // Output contact angle
455  //---------------------
456 
457  ofstream tangent_file;
458  sprintf(filename,"%s/tangent_to_contact_line%i.dat",
459  doc_info.directory().c_str(),
460  doc_info.number());
461  tangent_file.open(filename);
462 
463  ofstream normal_file;
464  sprintf(filename,"%s/normal_to_contact_line%i.dat",
465  doc_info.directory().c_str(),
466  doc_info.number());
467  normal_file.open(filename);
468 
469 
470  ofstream contact_angle_file;
471  sprintf(filename,"%s/contact_angle%i.dat",
472  doc_info.directory().c_str(),
473  doc_info.number());
474  contact_angle_file.open(filename);
475 
476  // Tangent and normal vectors to contact line
477  Vector<double> tangent(3);
478  Vector<double> normal(3);
479  Vector<double> r_contact(3);
480 
481  // How many contact angle elements are there?
482  unsigned n_element = Contact_angle_mesh_pt->nelement();
483 
484  // Loop over the surface elements
485  for(unsigned e=0;e<n_element;e++)
486  {
487 
488  tangent_file << "ZONE" << std::endl;
489  normal_file << "ZONE" << std::endl;
490  contact_angle_file << "ZONE" << std::endl;
491 
492  // Upcast from GeneralisedElement to YoungLaplace contact angle element
493  YoungLaplaceContactAngleElement<ELEMENT>* el_pt =
494  dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*>(
495  Contact_angle_mesh_pt->element_pt(e));
496 
497  // Loop over a few points in the contact angle element
498  Vector<double> s(1);
499  for (unsigned i=0;i<npts;i++)
500  {
501  s[0]=-1.0+2.0*double(i)/double(npts-1);
502 
503  dynamic_cast<ELEMENT*>(el_pt->bulk_element_pt())->
504  position(el_pt->local_coordinate_in_bulk(s),r_contact);
505 
506  el_pt->contact_line_vectors(s,tangent,normal);
507  tangent_file << r_contact[0] << " "
508  << r_contact[1] << " "
509  << r_contact[2] << " "
510  << tangent[0] << " "
511  << tangent[1] << " "
512  << tangent[2] << " " << std::endl;
513 
514  normal_file << r_contact[0] << " "
515  << r_contact[1] << " "
516  << r_contact[2] << " "
517  << normal[0] << " "
518  << normal[1] << " "
519  << normal[2] << " " << std::endl;
520 
521  contact_angle_file << r_contact[1] << " "
522  << el_pt->actual_cos_contact_angle(s)
523  << std::endl;
524  }
525 
526 
527  } // end of loop over both boundaries
528 
529  tangent_file.close();
530  normal_file.close();
531  contact_angle_file.close();
532 
533 
534 cout << "\n********************************************" << endl << endl;
535 
536 } // end of doc
537 
538 
539 
540 //===============start_of_main============================================
541 /// Drive code
542 //========================================================================
543 int main()
544 {
545 
546  // Create label for output
547  DocInfo doc_info;
548 
549  // Trace file
550  ofstream trace_file;
551 
552  // Set output directory
553  doc_info.set_directory("RESLT");
554 
555  // Open a trace file
556  char filename[100];
557  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
558  trace_file.open(filename);
559 
560  // Tecplot header for trace file: kappa and height value
561  trace_file << "VARIABLES=\"<GREEK>k</GREEK>\",\"h\"" << std::endl;
562  trace_file << "ZONE" << std::endl;
563 
564 
565  //Set up the problem
566  //------------------
567 
568  // Create the problem with 2D nine-node elements from the
569  // RefineableQYoungLaplaceElement family.
571 
572  // Perform one uniform refinement
573  problem.refine_uniformly();
574 
575  //Output the solution
576  problem.doc_solution(doc_info,trace_file);
577 
578  //Increment counter for solutions
579  doc_info.number()++;
580 
581  // Parameter incrementation
582  //-------------------------
583  double increment=0.1;
584 
585  // Loop over steps
586  unsigned nstep=2; // 10;
587  for (unsigned istep=0;istep<nstep;istep++)
588  {
590 
591  // Solve the problem
592  unsigned max_adapt=1;
593  problem.newton_solve(max_adapt);
594 
595  //Output the solution
596  problem.doc_solution(doc_info,trace_file);
597 
598  //Increment counter for solutions
599  doc_info.number()++;
600  }
601 
602  // Close output file
603  trace_file.close();
604 
605 } //end of main
606 
607 
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...
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.
Definition: barrel.cc:45
double Alpha_max
Max. spine angle against horizontal plane.
Definition: barrel.cc:99
double L_x
Length and width of the domain.
double Controlled_height
Height control value.
Definition: barrel.cc:51
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
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 Alpha_min
Min. spine angle against horizontal plane.
Definition: barrel.cc:96
double Cos_gamma
Cos of contact angle.
int main()
Drive code.