static_single_layer.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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Driver code for a 2D free-surface hydrostatics problem.
31 // The system consists of a layer of fluid
32 // in a domain of height 1 and half-width 0.5.
33 // The program solves for the interface position as the contact angle
34 // at the wall, alpha, decreases from pi/2. The resulting shapes should all be
35 // circular arcs and the pressure jump across the interface should be
36 // cos(alpha)/0.5 = 2 cos(alpha)/Ca.
37 
38 //OOMPH-LIB include files
39 #include "generic.h"
40 #include "navier_stokes.h"
41 #include "fluid_interface.h"
42 #include "constitutive.h"
43 #include "solid.h"
44 
45 // The mesh
46 #include "meshes/single_layer_spine_mesh.h"
47 
48 //Use the std namespace
49 using namespace std;
50 
51 using namespace oomph;
52 
53 //====start_of_namespace=================================
54 /// Namespace for phyical parameters
55 //=======================================================
56 namespace Global_Physical_Variables
57 {
58 
59  /// Pseudo-solid Poisson ratio
60  double Nu=0.1;
61 
62  ///Direction of the wall normal vector
63  Vector<double> Wall_normal;
64 
65  /// \short Function that specifies the wall unit normal
66  void wall_unit_normal_fct(const Vector<double> &x,
67  Vector<double> &normal)
68  {
69  normal=Wall_normal;
70  }
71 
72 } //end_of_namespace
73 
74 
75 
76 //============================================================================
77 ///A Problem class that solves the Navier--Stokes equations + free surface
78 ///in a 2D geometry using a spine-based node update
79 //============================================================================
80 template<class ELEMENT>
81 class CapProblem : public Problem
82 {
83 
84 public:
85 
86  //Constructor: Boolean flag indicates if volume constraint is
87  //applied by hijacking internal or external pressure
88  CapProblem(const bool& hijack_internal);
89 
90  //Destructor: clean up all allocated memory
91  ~CapProblem();
92 
93  /// Perform a parameter study: Solve problem for a range of contact angles
94  /// Pass name of output directory as a string
95  void parameter_study(const string& dir_name);
96 
97  /// Update the spine mesh after every Newton step
98  void actions_before_newton_convergence_check() {Bulk_mesh_pt->node_update();}
99 
100  /// Create the volume constraint elements
101  void create_volume_constraint_elements();
102 
103  /// Doc the solution
104  void doc_solution(DocInfo& doc_info);
105 
106 private:
107 
108  /// The Capillary number
109  double Ca;
110 
111  /// The volume of the fluid
112  double Volume;
113 
114  /// The external pressure
115  double Pext;
116 
117  /// The contact angle
118  double Angle;
119 
120  /// The bulk mesh of fluid elements
121  SingleLayerSpineMesh<SpineElement<ELEMENT> >* Bulk_mesh_pt;
122 
123  /// The mesh for the interface elements
125 
126  /// The mesh for the element at the contact point
128 
129  /// The volume constraint mesh
131 
132  /// Trace file
133  ofstream Trace_file;
134 
135  /// Data object whose single value stores the external pressure
137 
138  /// Data that is traded for the volume constraint
140 
141 };
142 
143 
144 
145 //======================================================================
146 /// Constructor: Pass boolean flag to indicate if the volume
147 /// constraint is applied by hijacking an internal pressure
148 /// or the external pressure
149 //======================================================================
150 template<class ELEMENT>
151 CapProblem<ELEMENT>::CapProblem(const bool& hijack_internal) :
152  Ca(2.1), //Initialise value of Ca to some random value
153  Volume(0.5), //Initialise the value of the volume
154  Pext(1.23), //Initialise the external pressure to some random value
155  Angle(0.5*MathematicalConstants::Pi) //Initialise the contact angle
156 {
157  //Set the wall normal
161 
162  // Number of elements in the horizontal direction
163  unsigned nx=4;
164 
165  // Number of elements in the vertical direction
166  unsigned nh=4;
167 
168  // Halfwidth of domain
169  double half_width=0.5;
170 
171  //Construct bulk mesh
172  Bulk_mesh_pt =
173  new SingleLayerSpineMesh<SpineElement<ELEMENT> >(nx,nh,half_width,1.0);
174 
175  //Create the surface mesh that will contain the interface elements
176  //First create storage, but with no elements or nodes
177  Surface_mesh_pt = new Mesh;
178  //Also create storage for a point mesh that contain the single element
179  //responsible for enforcing the contact angle condition
180  Point_mesh_pt = new Mesh;
181 
182  //Loop over the horizontal elements adjacent to the upper surface
183  //(boundary 2) and create the surface elements
184  unsigned n_boundary_element = Bulk_mesh_pt->nboundary_element(2);
185  for(unsigned e=0;e<n_boundary_element;e++)
186  {
187  //Construct a new 1D line element adjacent to boundary 2
188  FiniteElement *interface_element_pt =
189  new SpineLineFluidInterfaceElement<SpineElement<ELEMENT> >(
190  Bulk_mesh_pt->boundary_element_pt(2,e),
191  Bulk_mesh_pt->face_index_at_boundary(2,e));
192 
193  //Push it back onto the stack
194  this->Surface_mesh_pt->add_element_pt(interface_element_pt);
195 
196  //Find the (single) node that is on the solid boundary (boundary 1)
197  unsigned n_node = interface_element_pt->nnode();
198  //We only need to check the right-hand nodes (because I know the
199  //ordering of the nodes within the element)
200  if(interface_element_pt->node_pt(n_node-1)->is_on_boundary(1))
201  {
202  //Make the point (contact) element from right-hand edge of the element
203  FiniteElement* point_element_pt =
204  dynamic_cast<SpineLineFluidInterfaceElement<
205  SpineElement<ELEMENT> >*>(interface_element_pt)
206  ->make_bounding_element(1);
207 
208  //Add it to the mesh
209  this->Point_mesh_pt->add_element_pt(point_element_pt);
210  }
211  }
212 
213  //Create a Data object whose single value stores the
214  //external pressure
215  External_pressure_data_pt = new Data(1);
216 
217  // Set external pressure
218  External_pressure_data_pt->set_value(0,Pext);
219 
220  // Which pressure are we trading for the volume constraint: We
221  // can either hijack an internal pressure or use the external pressure.
222  if (hijack_internal)
223  {
224  // The external pressure is pinned -- the external pressure
225  // sets the pressure throughout the domain -- we do not have
226  // the liberty to fix another pressure value!
228 
229  //Hijack one of the pressure values in the fluid and use it
230  //as the pressure whose value is determined by the volume constraint.
231  //(Its value will affect the residual of that element but it will not
232  //be determined by it, i.e. it's hijacked).
233  Traded_pressure_data_pt = dynamic_cast<ELEMENT*>(
234  Bulk_mesh_pt->element_pt(0))->hijack_internal_value(0,0);
235  }
236  else
237  {
238  // Regard the external pressure as an unknown and add
239  // it to the problem's global data so it gets included
240  // in the equation numbering. Note that, at the moment,
241  // there's no equation that determines its value!
242  add_global_data(External_pressure_data_pt);
243 
244  // Declare the external pressure to be the pressure determined
245  // by the volume constraint, i.e. the pressure that's "traded":
247 
248  // Since the external pressure is "traded" for the volume constraint,
249  // it no longer sets the overall pressure, and we
250  // can add an arbitrary constant to all pressures. To make
251  // the solution unique, we pin a single pressure value in the bulk:
252  // We arbitrarily set the pressure dof 0 in element 0 to zero.
253  dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(0))->fix_pressure(0,0.0);
254  }
255 
256 
257 
258  // Loop over the elements on the interface to pass pointer to Ca and
259  // the pointers to the Data items that contains the single
260  // (pressure) value that is "traded" for the volume constraint
261  // and the external pressure
262  unsigned n_interface = Surface_mesh_pt->nelement();
263  for(unsigned e=0;e<n_interface;e++)
264  {
265  //Cast to a 1D element
266  SpineLineFluidInterfaceElement<SpineElement<ELEMENT> >*el_pt=
267  dynamic_cast<SpineLineFluidInterfaceElement<SpineElement<ELEMENT> >*>
268  (Surface_mesh_pt->element_pt(e));
269 
270  //Set the Capillary number
271  el_pt->ca_pt() = &Ca;
272 
273  //Pass the Data item that contains the single external pressure value
274  el_pt->set_external_pressure_data(External_pressure_data_pt);
275  }
276 
277 
278  //Set the boundary conditions
279 
280  //Pin the velocities on all boundaries apart from the free surface
281  //(boundary 2) where all velocities are free, and apart from the symmetry
282  //line (boundary 3) where only the horizontal velocity is pinned
283  unsigned n_bound=Bulk_mesh_pt->nboundary();
284  for (unsigned b=0;b<n_bound;b++)
285  {
286  if (b!=2)
287  {
288  //Find the number of nodes on the boundary
289  unsigned n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
290  //Loop over the nodes on the boundary
291  for(unsigned n=0;n<n_boundary_node;n++)
292  {
293  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
294  if (b!=3)
295  {
296  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
297  }
298  }
299  }
300  }
301 
302  // Set the contact angle boundary condition for the rightmost element
303  // (pass pointer to double that specifies the contact angle)
304  FluidInterfaceBoundingElement *contact_angle_element_pt
305  = dynamic_cast<FluidInterfaceBoundingElement*>(
306  Point_mesh_pt->element_pt(0));
307 
308  contact_angle_element_pt->set_contact_angle(&Angle);
309  contact_angle_element_pt->ca_pt() = &Ca;
310  contact_angle_element_pt->wall_unit_normal_fct_pt()
312 
313  //Now add the volume constraint
315 
316  this->add_sub_mesh(Bulk_mesh_pt);
317  this->add_sub_mesh(Surface_mesh_pt);
318  this->add_sub_mesh(Point_mesh_pt);
319  this->add_sub_mesh(Volume_constraint_mesh_pt);
320 
321  this->build_global_mesh();
322 
323  //Setup all the equation numbering and look-up schemes
324  cout << "Number of unknowns: " << assign_eqn_numbers() << std::endl;
325 
326 }
327 
328 
329 //==========================================================================
330 /// Destructor. Make sure to clean up all allocated memory, so that multiple
331 /// instances of the problem don't lead to excessive memory usage.
332 //==========================================================================
333 template<class ELEMENT>
335 {
336  //Loop over the volume mesh and delete the elements
337  unsigned n_element = Volume_constraint_mesh_pt->nelement();
338  for(unsigned e=0;e<n_element;e++)
339  {delete Volume_constraint_mesh_pt->element_pt(e);}
340  //Now flush the storage
341  Volume_constraint_mesh_pt->flush_element_and_node_storage();
342  //Now delete the mesh
343  delete Volume_constraint_mesh_pt;
344 
345  //Delete the traded pressure if not the same as the external pressure
346  if(Traded_pressure_data_pt!=External_pressure_data_pt)
347  {delete Traded_pressure_data_pt;}
348  //Next delete the external data
349  delete External_pressure_data_pt;
350 
351  //Loop over the point mesh and delete the elements
352  n_element = Point_mesh_pt->nelement();
353  for(unsigned e=0;e<n_element;e++)
354  {delete Point_mesh_pt->element_pt(e);}
355  //Now flush the storage
356  Point_mesh_pt->flush_element_and_node_storage();
357  //Then delete the mesh
358  delete Point_mesh_pt;
359 
360  //Loop over the surface mesh and delete the elements
361  n_element = Surface_mesh_pt->nelement();
362  for(unsigned e=0;e<n_element;e++)
363  {delete Surface_mesh_pt->element_pt(e);}
364  //Now flush the storage
365  Surface_mesh_pt->flush_element_and_node_storage();
366  //Then delete the mesh
367  delete Surface_mesh_pt;
368 
369  //Then delete the bulk mesh
370  delete Bulk_mesh_pt;
371 }
372 
373 
374 //=======================================================================
375 /// Create the volume constraint elements
376 //========================================================================
377 template<class ELEMENT>
379 {
380  //The single volume constraint element
381  Volume_constraint_mesh_pt = new Mesh;
382  VolumeConstraintElement* vol_constraint_element =
383  new VolumeConstraintElement(&Volume,Traded_pressure_data_pt,0);
384  Volume_constraint_mesh_pt->add_element_pt(vol_constraint_element);
385 
386  //Loop over all boundaries (or just 1 and 2 why?)
387  for(unsigned b=0;b<4;b++)
388  {
389  // How many bulk fluid elements are adjacent to boundary b?
390  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
391 
392  // Loop over the bulk fluid elements adjacent to boundary b?
393  for(unsigned e=0;e<n_element;e++)
394  {
395  // Get pointer to the bulk fluid element that is
396  // adjacent to boundary b
397  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
398  Bulk_mesh_pt->boundary_element_pt(b,e));
399 
400  //Find the index of the face of element e along boundary b
401  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
402 
403  // Create new element
404  SpineLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
405  new SpineLineVolumeConstraintBoundingElement<ELEMENT>(
406  bulk_elem_pt,face_index);
407 
408  //Set the "master" volume control element
409  el_pt->set_volume_constraint_element(vol_constraint_element);
410 
411  // Add it to the mesh
412  Volume_constraint_mesh_pt->add_element_pt(el_pt);
413  }
414  }
415 }
416 
417 
418 //======================================================================
419 /// Perform a parameter study. Pass name of output directory as
420 /// a string
421 //======================================================================
422 template<class ELEMENT>
423 void CapProblem<ELEMENT>::parameter_study(const string& dir_name)
424 {
425 
426  // Create DocInfo object (allows checking if output directory exists)
427  DocInfo doc_info;
428  doc_info.set_directory(dir_name);
429  doc_info.number()=0;
430 
431 
432  // Open trace file
433  char filename[100];
434  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
435  Trace_file.open(filename);
436  Trace_file << "VARIABLES=\"<greek>a</greek><sub>prescribed</sub>\",";
437  Trace_file << "\"h<sub>left</sub>\",\"h<sub>right</sub>\",";
438  Trace_file << "\"p<sub>fluid</sub>-p<sub>ext</sub>\",";
439  Trace_file << "\"<greek>D</greek>p<sub>exact</sub>\"";
440  Trace_file << std::endl;
441 
442  //Gradually increase the contact angle
443  for(unsigned i=0;i<6;i++)
444  {
445  //Solve the problem
446  steady_newton_solve();
447 
448  //Output result
449  doc_solution(doc_info);
450 
451  // Bump up counter
452  doc_info.number()++;
453 
454  //Decrease the contact angle
455  Angle -= 5.0*MathematicalConstants::Pi/180.0;
456  }
457 
458 }
459 
460 
461 
462 
463 //========================================================================
464 /// Doc the solution
465 //========================================================================
466 template<class ELEMENT>
467 void CapProblem<ELEMENT>::doc_solution(DocInfo& doc_info)
468 {
469 
470  ofstream some_file;
471  char filename[100];
472 
473  // Number of plot points
474  unsigned npts=5;
475 
476 
477  //Output domain
478  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
479  doc_info.number());
480  some_file.open(filename);
481  Bulk_mesh_pt->output(some_file,npts);
482  Surface_mesh_pt->output(some_file,npts);
483  Point_mesh_pt->output(some_file,npts);
484  some_file.close();
485 
486  // Number of spines
487  unsigned nspine=Bulk_mesh_pt->nspine();
488 
489  // Doc
490  Trace_file << Angle*180.0/MathematicalConstants::Pi;
491  Trace_file << " " << Bulk_mesh_pt->spine_pt(0)->height();
492  Trace_file << " " << Bulk_mesh_pt->spine_pt(nspine-1)->height();
493  Trace_file << " "
494  << dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(0))
495  ->p_nst(0)- External_pressure_data_pt->value(0);
496  Trace_file << " " << -2.0*cos(Angle)/Ca;
497  Trace_file << std::endl;
498 
499 }
500 
501 
502 //===========start_of_pseudo_elastic_class====================================
503 ///\short A class that solves the Navier--Stokes equations
504 ///to compute the shape of a static interface in a rectangular container
505 ///with imposed contact angle at the boundary.
506 //============================================================================
507 template<class ELEMENT>
508 class PseudoSolidCapProblem : public Problem
509 {
510 public:
511 
512  ///Constructor: Boolean flag indicates if volume constraint is
513  ///applied by hijacking internal or external pressure
514  PseudoSolidCapProblem(const bool& hijack_internal);
515 
516  /// Destructor: clean up memory allocated by the object
518 
519  /// Peform a parameter study: Solve problem for a range of contact angles
520  /// Pass name of output directory as a string
521  void parameter_study(const string& dir_name);
522 
523  /// Doc the solution
524  void doc_solution(DocInfo& doc_info);
525 
526 private:
527 
528  /// Create the free surface elements
530 
531  /// Create the volume constraint elements
533 
534  /// Create the contact angle element
536 
537  /// The Capillary number
538  double Ca;
539 
540  /// The prescribed volume of the fluid
541  double Volume;
542 
543  /// The external pressure
544  double Pext;
545 
546  /// The contact angle
547  double Angle;
548 
549  /// Constitutive law used to determine the mesh deformation
550  ConstitutiveLaw *Constitutive_law_pt;
551 
552  /// Data object whose single value stores the external pressure
554 
555  // Pointer to the (single valued) Data item that
556  // will contain the pressure value that we're
557  // trading for the volume constraint
559 
560  /// Trace file
561  ofstream Trace_file;
562 
563  ///Storage for the bulk mesh
565 
566  /// Storage for the free surface mesh
568 
569  /// Storage for the element bounding the free surface
571 
572  /// Storage for the elements that compute the enclosed volume
574 
575  /// Storage for the volume constraint
577 
578 }; //end_of_pseudo_solid_problem_class
579 
580 
581 
582 //============start_of_constructor=====================================
583 /// Constructor: Pass boolean flag to indicate if the volume
584 /// constraint is applied by hijacking an internal pressure
585 /// or the external pressure
586 //======================================================================
587 template<class ELEMENT>
589  Ca(2.1), //Initialise value of Ca to some random value
590  Volume(0.5), //Initialise the value of the volume
591  Pext(1.23), //Initialise the external pressure to some random value
592  Angle(0.5*MathematicalConstants::Pi) //Initialise the contact angle
593 {
594  //Set the wall normal
598 
599  // Number of elements in the horizontal direction
600  unsigned nx=4;
601 
602  // Number of elements in the vertical direction
603  unsigned nh=4;
604 
605  // Halfwidth of domain
606  double half_width=0.5;
607 
608  //Construct mesh
609  Bulk_mesh_pt = new ElasticRectangularQuadMesh<ELEMENT>(nx,nh,half_width,1.0);
610 
611  //Create a Data object whose single value stores the
612  //external pressure
613  External_pressure_data_pt = new Data(1);
614 
615  // Set external pressure
616  External_pressure_data_pt->set_value(0,Pext);
617 
618  // Which pressure are we trading for the volume constraint: We
619  // can either hijack an internal pressure or use the external pressure.
620  if (hijack_internal)
621  {
622  // The external pressure is pinned -- the external pressure
623  // sets the pressure throughout the domain -- we do not have
624  // the liberty to fix another pressure value!
626 
627  //Hijack one of the pressure values in the fluid and use it
628  //as the pressure whose value is determined by the volume constraint.
629  //(Its value will affect the residual of that element but it will not
630  //be determined by it, i.e. it's hijacked).
631  Traded_pressure_data_pt = dynamic_cast<ELEMENT*>(
632  Bulk_mesh_pt->element_pt(0))->hijack_internal_value(0,0);
633  }
634  else
635  {
636  // Regard the external pressure is an unknown and add
637  // it to the problem's global data so it gets included
638  // in the equation numbering. Note that, at the moment,
639  // there's no equation that determines its value!
640  add_global_data(External_pressure_data_pt);
641 
642  // Declare the external pressure to be the pressure determined
643  // by the volume constraint, i.e. the pressure that's "traded":
645 
646  // Since the external pressure is "traded" for the volume constraint,
647  // it no longer sets the overall pressure, and we
648  // can add an arbitrary constant to all pressures. To make
649  // the solution unique, we pin a single pressure value in the bulk:
650  // We arbitrarily set the pressure dof 0 in element 0 to zero.
651  dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(0))->fix_pressure(0,0.0);
652  }
653 
654  //Set the constituive law
656  new GeneralisedHookean(&Global_Physical_Variables::Nu);
657 
658  //Loop over the elements to set the consitutive law and jacobian
659  unsigned n_bulk = Bulk_mesh_pt->nelement();
660  for(unsigned e=0;e<n_bulk;e++)
661  {
662  ELEMENT* el_pt =
663  dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
664 
665  el_pt->constitutive_law_pt() = Constitutive_law_pt;
666  }
667 
668  //Set the boundary conditions
669 
670  //Fluid velocity conditions
671  //Pin the velocities on all boundaries apart from the free surface
672  //(boundary 2) where all velocities are free, and apart from the symmetry
673  //line (boundary 3) where only the horizontal velocity is pinned
674  unsigned n_bound=Bulk_mesh_pt->nboundary();
675  for (unsigned b=0;b<n_bound;b++)
676  {
677  if (b!=2)
678  {
679  //Find the number of nodes on the boundary
680  unsigned n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
681  //Loop over the nodes on the boundary
682  for(unsigned n=0;n<n_boundary_node;n++)
683  {
684  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
685  if (b!=3)
686  {
687  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
688  }
689  }
690  }
691  } //end_of_fluid_boundary_conditions
692 
693  //PesudoSolid boundary conditions
694  for (unsigned b=0;b<n_bound;b++)
695  {
696  if (b!=2)
697  {
698  //Find the number of nodes on the boundary
699  unsigned n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
700  //Loop over the nodes on the boundary
701  for(unsigned n=0;n<n_boundary_node;n++)
702  {
703  //Pin vertical displacement on the bottom
704  if(b==0)
705  {
706  static_cast<SolidNode*>(Bulk_mesh_pt->boundary_node_pt(b,n))
707  ->pin_position(1);
708  }
709  if((b==1) || (b==3))
710  {
711  //Pin horizontal displacement on the sizes
712  static_cast<SolidNode*>(Bulk_mesh_pt->boundary_node_pt(b,n))
713  ->pin_position(0);
714  }
715  }
716  } //end_of_solid_boundary_conditions
717  }
718 
719  //Constrain all nodes only to move vertically (not horizontally)
720  {
721  unsigned n_node = Bulk_mesh_pt->nnode();
722  for(unsigned n=0;n<n_node;n++)
723  {
724  static_cast<SolidNode*>(Bulk_mesh_pt->node_pt(n))->pin_position(0);
725  }
726  } //end_of_constraint
727 
728  //Create the free surface elements
730 
731  //Create the volume constraint elements
733 
734  //Need to make the bounding element
736 
737  //Now need to add all the meshes
738  this->add_sub_mesh(Bulk_mesh_pt);
739  this->add_sub_mesh(Free_surface_mesh_pt);
740  this->add_sub_mesh(Volume_computation_mesh_pt);
741  this->add_sub_mesh(Volume_constraint_mesh_pt);
742  this->add_sub_mesh(Free_surface_bounding_mesh_pt);
743 
744  //and build the global mesh
745  this->build_global_mesh();
746 
747  //Setup all the equation numbering and look-up schemes
748  cout << "Number of unknowns: " << assign_eqn_numbers() << std::endl;
749 
750 } //end_of_constructor
751 
752 
753 //==========================================================================
754 /// Destructor. Make sure to clean up all allocated memory, so that multiple
755 /// instances of the problem don't lead to excessive memory usage.
756 //==========================================================================
757 template<class ELEMENT>
759 {
760  //Delete the contact angle element
761  delete Free_surface_bounding_mesh_pt->element_pt(0);
762  Free_surface_bounding_mesh_pt->flush_element_and_node_storage();
763  delete Free_surface_bounding_mesh_pt;
764  //Delete the volume constraint mesh
765  delete Volume_constraint_mesh_pt;
766  //Delete the surface volume computation elements
767  unsigned n_element = Volume_computation_mesh_pt->nelement();
768  for(unsigned e=0;e<n_element;e++)
769  {delete Volume_computation_mesh_pt->element_pt(e);}
770  //Now flush the storage
771  Volume_computation_mesh_pt->flush_element_and_node_storage();
772  //Now delete the mesh
773  delete Volume_computation_mesh_pt;
774  //Delete the free surface elements
775  n_element = Free_surface_mesh_pt->nelement();
776  for(unsigned e=0;e<n_element;e++)
777  {delete Free_surface_mesh_pt->element_pt(e);}
778  //Now flush the storage
779  Free_surface_mesh_pt->flush_element_and_node_storage();
780  //Now delete the mesh
781  delete Free_surface_mesh_pt;
782 
783  //Delete the constitutive law
784  delete Constitutive_law_pt;
785 
786  //If not the same as the external pressure, delete the traded pressure
787  if(Traded_pressure_data_pt!=External_pressure_data_pt)
788  {delete Traded_pressure_data_pt;}
789  //Next delete the external data
790  delete External_pressure_data_pt;
791  //Then delete the bulk mesh
792  delete Bulk_mesh_pt;
793 }
794 
795 //============create_free_surface_element================================
796 /// Create the free surface elements
797 //========================================================================
798 template<class ELEMENT>
800 {
801  //Allocate storage for the free surface mesh
802  Free_surface_mesh_pt = new Mesh;
803 
804  //Loop over boundary 2 and create the interface elements
805  unsigned b=2;
806 
807  // How many bulk fluid elements are adjacent to boundary b?
808  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
809 
810  // Loop over the bulk fluid elements adjacent to boundary b?
811  for(unsigned e=0;e<n_element;e++)
812  {
813  // Get pointer to the bulk fluid element that is
814  // adjacent to boundary b
815  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
816  Bulk_mesh_pt->boundary_element_pt(b,e));
817 
818  //Find the index of the face of element e along boundary b
819  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
820 
821  // Create new element
822  ElasticLineFluidInterfaceElement<ELEMENT>* el_pt =
823  new ElasticLineFluidInterfaceElement<ELEMENT>(
824  bulk_elem_pt,face_index);
825 
826  // Add it to the mesh
827  Free_surface_mesh_pt->add_element_pt(el_pt);
828 
829  //Add the appropriate boundary number
830  el_pt->set_boundary_number_in_bulk_mesh(b);
831 
832  //Add the capillary number
833  el_pt->ca_pt() = &Ca;
834 
835  //Add the external pressure data
836  el_pt->set_external_pressure_data(External_pressure_data_pt);
837  }
838 
839 }
840 
841 
842 //============start_of_create_volume_constraint_elements==================
843 /// Create the volume constraint elements
844 //========================================================================
845 template<class ELEMENT>
847 {
848  //Build the single volume constraint element
849  Volume_constraint_mesh_pt = new Mesh;
850  VolumeConstraintElement* vol_constraint_element =
851  new VolumeConstraintElement(&Volume,Traded_pressure_data_pt,0);
852  Volume_constraint_mesh_pt->add_element_pt(vol_constraint_element);
853 
854  //Now create the volume computation elements
855  Volume_computation_mesh_pt = new Mesh;
856 
857  //Loop over all boundaries
858  for(unsigned b=0;b<4;b++)
859  {
860  // How many bulk fluid elements are adjacent to boundary b?
861  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
862 
863  // Loop over the bulk fluid elements adjacent to boundary b?
864  for(unsigned e=0;e<n_element;e++)
865  {
866  // Get pointer to the bulk fluid element that is
867  // adjacent to boundary b
868  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
869  Bulk_mesh_pt->boundary_element_pt(b,e));
870 
871  //Find the index of the face of element e along boundary b
872  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
873 
874  // Create new element
875  ElasticLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
876  new ElasticLineVolumeConstraintBoundingElement<ELEMENT>(
877  bulk_elem_pt,face_index);
878 
879  //Set the "master" volume control element
880  el_pt->set_volume_constraint_element(vol_constraint_element);
881 
882  // Add it to the mesh
883  Volume_computation_mesh_pt->add_element_pt(el_pt);
884  }
885  }
886 } //end_of_create_volume_constraint_elements
887 
888 //==========start_of_create_contact_angle_elements========================
889 /// Create the contact angle element
890 //========================================================================
891 template<class ELEMENT>
893 {
894  Free_surface_bounding_mesh_pt = new Mesh;
895 
896  //Find the element at the end of the free surface
897  //The elements are assigned in order of increasing x coordinate
898  unsigned n_free_surface = Free_surface_mesh_pt->nelement();
899 
900  //Make the bounding element for the contact angle constraint
901  //which works because the order of elements in the mesh is known
902  FluidInterfaceBoundingElement* el_pt =
903  dynamic_cast<ElasticLineFluidInterfaceElement<ELEMENT>*>
904  (Free_surface_mesh_pt->element_pt(n_free_surface-1))->
905  make_bounding_element(1);
906 
907  //Set the contact angle (strong imposition)
908  el_pt->set_contact_angle(&Angle);
909 
910  //Set the capillary number
911  el_pt->ca_pt() = &Ca;
912 
913  //Set the wall normal of the external boundary
914  el_pt->wall_unit_normal_fct_pt()
916 
917  //Add the element to the mesh
918  Free_surface_bounding_mesh_pt->add_element_pt(el_pt);
919 
920 } //end_of_create_contact_angle_element
921 
922 
923 
924 //================start_of_parameter_study===========================
925 /// Perform a parameter study. Pass name of output directory as
926 /// a string
927 //======================================================================
928 template<class ELEMENT>
930 {
931  // Create DocInfo object (allows checking if output directory exists)
932  DocInfo doc_info;
933  doc_info.set_directory(dir_name);
934  doc_info.number()=0;
935 
936  // Open trace file
937  char filename[100];
938  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
939  Trace_file.open(filename);
940  Trace_file << "VARIABLES=\"<greek>a</greek><sub>prescribed</sub>\",";
941  Trace_file << "\"h<sub>left</sub>\",\"h<sub>right</sub>\",";
942  Trace_file << "\"p<sub>fluid</sub>-p<sub>ext</sub>\",";
943  Trace_file << "\"<greek>D</greek>p<sub>exact</sub>\"";
944  Trace_file << std::endl;
945 
946  //Solve the problem for six different contact angles
947  for(unsigned i=0;i<6;i++)
948  {
949  //Solve the problem
950  steady_newton_solve();
951 
952  //Output result
953  doc_solution(doc_info);
954 
955  // Bump up counter
956  doc_info.number()++;
957 
958  //Decrease the contact angle
959  Angle -= 5.0*MathematicalConstants::Pi/180.0;
960  }
961 
962 }
963 
964 
965 
966 
967 //==============start_of_doc_solution=====================================
968 /// Doc the solution
969 //========================================================================
970 template<class ELEMENT>
972 {
973  //Output stream
974  ofstream some_file;
975  char filename[100];
976 
977  // Number of plot points
978  unsigned npts=5;
979 
980  //Output domain
981  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
982  doc_info.number());
983  some_file.open(filename);
984  Bulk_mesh_pt->output(some_file,npts);
985  some_file.close();
986 
987 
988  // Number of interface elements
989  unsigned ninterface=Free_surface_mesh_pt->nelement();
990  //Find number of nodes in the last interface element
991  unsigned np = Free_surface_mesh_pt->finite_element_pt(ninterface-1)->nnode();
992  // Document the contact angle (in degrees), the height of the interface at
993  // the centre of the container, the height at the wall, the computed
994  // pressure drop across the interface and
995  // the analytic prediction of the pressure drop.
996  Trace_file << Angle*180.0/MathematicalConstants::Pi;
997  Trace_file << " " << Free_surface_mesh_pt->finite_element_pt(0)
998  ->node_pt(0)->x(1)
999  << " "
1000  << Free_surface_mesh_pt->finite_element_pt(ninterface-1)
1001  ->node_pt(np-1)->x(1);
1002  Trace_file << " "
1003  << dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(0))->p_nst(0)-
1004  External_pressure_data_pt->value(0);
1005  Trace_file << " " << -2.0*cos(Angle)/Ca;
1006  Trace_file << std::endl;
1007 
1008 } //end_of_doc_solution
1009 
1010 
1011 
1012 
1013 //===start_of_main=======================================================
1014 ///Main driver: Build problem and initiate parameter study
1015 //======================================================================
1016 int main()
1017 {
1018 
1019  // Solve the problem twice, once hijacking an internal, once
1020  // hijacking the external pressure
1021  for (unsigned i=0;i<2;i++)
1022  {
1023 
1024  bool hijack_internal=false;
1025  if (i==1) hijack_internal=true;
1026  //Construct the problem
1027  CapProblem<Hijacked<QCrouzeixRaviartElement<2> > > problem(hijack_internal);
1028 
1029  string dir_name="RESLT_hijacked_external";
1030  if (i==1) dir_name="RESLT_hijacked_internal";
1031 
1032  //Do parameter study
1033  problem.parameter_study(dir_name);
1034 
1035  }
1036 
1037 
1038  // Solve the pseudosolid problem twice, once hijacking an internal, once
1039  // hijacking the external pressure
1040  for (unsigned i=0;i<2;i++)
1041  {
1042  bool hijack_internal=false;
1043  if (i==1) hijack_internal=true;
1044  //Construct the problem
1045  PseudoSolidCapProblem<Hijacked<
1046  PseudoSolidNodeUpdateElement<QCrouzeixRaviartElement<2>,
1047  QPVDElementWithPressure<2> > > > problem(hijack_internal);
1048 
1049  string dir_name="RESLT_elastic_hijacked_external";
1050  if (i==1) dir_name="RESLT_elastic_hijacked_internal";
1051 
1052  //Do parameter study
1053  problem.parameter_study(dir_name);
1054  }
1055 
1056 } //end_of_main
void parameter_study(const string &dir_name)
double Pext
The external pressure.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
void create_contact_angle_element()
Create the contact angle element.
Vector< double > Wall_normal
Direction of the wall normal vector.
void actions_before_newton_convergence_check()
Update the spine mesh after every Newton step.
double Pext
The external pressure.
void parameter_study()
Peform a parameter study: Solve problem for a range of contact angles.
void create_free_surface_elements()
Create the free surface elements.
void parameter_study(const string &dir_name)
double Angle
The contact angle.
Mesh * Free_surface_bounding_mesh_pt
Storage for the element bounding the free surface.
void create_volume_constraint_elements()
Create the volume constraint elements.
Mesh * Surface_mesh_pt
The mesh for the interface elements.
Mesh * Free_surface_mesh_pt
Storage for the free surface mesh.
ofstream Trace_file
Trace file.
void create_volume_constraint_elements()
Create the volume constraint elements.
double Volume
The prescribed volume of the fluid.
Data * External_pressure_data_pt
Data object whose single value stores the external pressure.
void wall_unit_normal_fct(const Vector< double > &x, Vector< double > &normal)
Function that specifies the wall unit normal.
SingleLayerSpineMesh< SpineElement< ELEMENT > > * Bulk_mesh_pt
The bulk mesh of fluid elements.
Mesh * Volume_constraint_mesh_pt
Storage for the volume constraint.
CapProblem(const bool &hijack_internal)
~PseudoSolidCapProblem()
Destructor: clean up memory allocated by the object.
ofstream Trace_file
Trace file.
PseudoSolidCapProblem(const bool &hijack_internal)
void doc_solution(DocInfo &doc_info)
Doc the solution.
A class that solves the Navier–Stokes equations to compute the shape of a static interface in a recta...
Mesh * Volume_constraint_mesh_pt
The volume constraint mesh.
void doc_solution(DocInfo &doc_info)
Doc the solution.
int main()
Main driver: Build problem and initiate parameter study.
double Ca
The Capillary number.
double Volume
The volume of the fluid.
double Ca
The Capillary number.
Mesh * Volume_computation_mesh_pt
Storage for the elements that compute the enclosed volume.
Data * External_pressure_data_pt
Data object whose single value stores the external pressure.
Data * Traded_pressure_data_pt
Data that is traded for the volume constraint.
double Nu
Pseudo-solid Poisson ratio.
Mesh * Bulk_mesh_pt
Storage for the bulk mesh.
double Angle
The contact angle.
Mesh * Point_mesh_pt
The mesh for the element at the contact point.