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