static_two_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 interface hydrostatics problem.
31 // The system consists of equal volumes of two fluids with
32 // the same physical properties, but a non-zero surface tension at
33 // their interface, in a domain of height 2 and half-width 0.5.
34 // The program solves for the interface position as the contact angle
35 // at the wall, alpha, decreases from pi/2. The resulting shapes should all be
36 // circular arcs and the pressure jump across the interface should be
37 // cos(alpha)/0.5 = 2 cos(alpha)/Ca.
38 
39 //OOMPH-LIB include files
40 #include "generic.h"
41 #include "navier_stokes.h"
42 #include "fluid_interface.h"
43 #include "constitutive.h"
44 #include "solid.h"
45 
46 // The mesh
47 #include "meshes/two_layer_spine_mesh.h"
48 #include "meshes/rectangular_quadmesh.h"
49 
50 //Use the std namespace
51 using namespace std;
52 
53 using namespace oomph;
54 
55 //The global physical variables
56 namespace Global_Physical_Variables
57 {
58  /// Pseudo-solid Poisson ratio
59  double Nu=0.1;
60 
61  ///Direction of the wall normal vector
62  Vector<double> Wall_normal;
63 
64  /// \short Function that specifies the wall unit normal
65  void wall_unit_normal_fct(const Vector<double> &x,
66  Vector<double> &normal)
67  {
68  normal=Wall_normal;
69  }
70 
71 }
72 
73 //============================================================================
74 ///A Problem class that solves the Navier--Stokes equations + free surface
75 ///in a 2D geometry using a spine-based node update
76 //============================================================================
77 template<class ELEMENT>
78 class CapProblem : public Problem
79 {
80 
81 public:
82 
83  //Constructor:
84  //Nx: Number of elements in the x (horizontal) direction
85  //Nh1: Number of elements in the y (vertical) direction in the lower fluid
86  //Nh2: Number of elements in the y (vertical) direction in the upper fluid
87  CapProblem(const unsigned &Nx, const unsigned &Nh1,
88  const unsigned &Nh2);
89 
90  /// Peform a parameter study: Solve problem for a range of contact angles
91  void parameter_study();
92 
93  /// Finish full specification of the elements
94  void finish_problem_setup();
95 
96  /// Pointer to the mesh
97  TwoLayerSpineMesh<SpineElement<ELEMENT> >* Bulk_mesh_pt;
98 
99  /// Pointer to the surface mesh
100  Mesh* Surface_mesh_pt;
101 
102  /// Pointer to the point mesh
103  Mesh* Point_mesh_pt;
104 
105  /// Update the spine mesh after every Newton step
106  void actions_before_newton_convergence_check() {Bulk_mesh_pt->node_update();}
107 
108  /// Create the volume constraint elements
109  void create_volume_constraint_elements();
110 
111  /// Doc the solution
112  void doc_solution(DocInfo& doc_info);
113 
114 private:
115 
116  /// The Capillary number
117  double Ca;
118 
119  /// The volume of the fluid
120  double Volume;
121 
122  /// The contact angle
123  double Angle;
124 
125  /// The volume constraint mesh
126  Mesh* Volume_constraint_mesh_pt;
127 
128  /// Trace file
129  ofstream Trace_file;
130 
131  /// Data that is traded for the volume constraint
132  Data* Traded_pressure_data_pt;
133 
134 };
135 
136 
137 
138 //======================================================================
139 /// Constructor
140 //======================================================================
141 template<class ELEMENT>
143 (const unsigned &Nx, const unsigned &Nh1, const unsigned &Nh2) :
144  Ca(2.1), //Initialise value of Ca to some random value
145  Volume(0.5), //Initialise the value of the volume
146  Angle(0.5*MathematicalConstants::Pi) //Initialise the contact angle
147 {
151 
152  //Construct our mesh
153  Bulk_mesh_pt = new TwoLayerSpineMesh<SpineElement<ELEMENT> >
154  (Nx,Nh1,Nh2,0.5,1.0,1.0);
155 
156  Surface_mesh_pt = new Mesh;
157 
158  //Loop over the horizontal elements
159  unsigned n_interface_lower = Bulk_mesh_pt->ninterface_lower();
160  for(unsigned i=0;i<n_interface_lower;i++)
161  {
162  //Construct a new 1D line element adjacent to the interface
163  FiniteElement *interface_element_pt = new
164  SpineLineFluidInterfaceElement<SpineElement<ELEMENT> >(
165  Bulk_mesh_pt->interface_lower_boundary_element_pt(i),
166  Bulk_mesh_pt->interface_lower_face_index_at_boundary(i));
167 
168  this->Surface_mesh_pt->add_element_pt(interface_element_pt);
169  }
170 
171 //Create the point mesh
172 Point_mesh_pt = new Mesh;
173  {
174  //Make the point (contact) element from the last surface element
175  FiniteElement* point_element_pt =
176  dynamic_cast<SpineLineFluidInterfaceElement<
177  SpineElement<ELEMENT> >*>(Surface_mesh_pt->element_pt(n_interface_lower-1))
178  ->make_bounding_element(1);
179 
180  //Add it to the mesh
181  this->Point_mesh_pt->add_element_pt(point_element_pt);
182 }
183 
184  //Hijack one of the pressure values in the upper fluid. Its value
185  //will affect the residual of that element but it will not
186  //be determined by it!
187  Traded_pressure_data_pt = dynamic_cast<ELEMENT*>(
188  Bulk_mesh_pt->upper_layer_element_pt(0))->hijack_internal_value(0,0);
189 
190  // Loop over the elements on the interface to pass pointer to Ca and
191  // the pointer to the Data item that contains the single
192  // (pressure) value that is "traded" for the volume constraint
193  unsigned n_interface = Surface_mesh_pt->nelement();
194  for(unsigned e=0;e<n_interface;e++)
195  {
196  //Cast to a 1D element
197  SpineLineFluidInterfaceElement<SpineElement<ELEMENT> >*el_pt=
198  dynamic_cast<SpineLineFluidInterfaceElement<
199  SpineElement<ELEMENT> >*>
200  (Surface_mesh_pt->element_pt(e));
201  //Set the Capillary number
202  el_pt->ca_pt() = &Ca;
203  }
204 
205 
206  //Set the boundary conditions
207 
208  //Pin the velocities on all external boundaries apart from the symmetry
209  //line (boundaries 4 and 5) where only the horizontal velocity is pinned
210  for (unsigned b=0;b<6;b++)
211  {
212  //Find the number of nodes on the boundary
213  unsigned n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
214  //Loop over the nodes on the boundary
215  for(unsigned n=0;n<n_boundary_node;n++)
216  {
217  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
218  if ((b!=4) && (b!=5))
219  {
220  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
221  }
222  }
223  }
224  // Set the contact angle boundary condition for the rightmost element
225  dynamic_cast<FluidInterfaceBoundingElement*>(
226  Point_mesh_pt->element_pt(0))->set_contact_angle(&Angle);
227 
228  dynamic_cast<FluidInterfaceBoundingElement*>(
229  Point_mesh_pt->element_pt(0))->ca_pt() = &Ca;
230 
231  dynamic_cast<FluidInterfaceBoundingElement*>(
232  Point_mesh_pt->element_pt(0))->wall_unit_normal_fct_pt() =
234 
235  // Pin a single pressure value: Set the pressure dof 0 in element 0
236  // of the lower layer to zero.
237  dynamic_cast<ELEMENT*>(Bulk_mesh_pt->lower_layer_element_pt(0))
238  ->fix_pressure(0,0.0);
239 
240  this->create_volume_constraint_elements();
241 
242  this->add_sub_mesh(Bulk_mesh_pt);
243  this->add_sub_mesh(Surface_mesh_pt);
244  this->add_sub_mesh(Point_mesh_pt);
245  this->add_sub_mesh(Volume_constraint_mesh_pt);
246 
247  this->build_global_mesh();
248 
249  //Setup all the equation numbering and look-up schemes
250  cout << "Number of unknowns: " << assign_eqn_numbers() << std::endl;
251 
252 }
253 
254 
255 
256 //======================================================================
257 /// Create the volume constraint elements
258 //========================================================================
259 template<class ELEMENT>
261 {
262  //The single volume constraint element
263  Volume_constraint_mesh_pt = new Mesh;
264  VolumeConstraintElement* vol_constraint_element =
265  new VolumeConstraintElement(&Volume,Traded_pressure_data_pt,0);
266  Volume_constraint_mesh_pt->add_element_pt(vol_constraint_element);
267 
268  //Loop over lower boundaries (or just 1, why?)
269  unsigned lower_boundaries[3]={0,1,5};
270  for(unsigned i=0;i<3;i++)
271  {
272  //Read out the actual boundaries
273  unsigned b = lower_boundaries[i];
274 
275  // How many bulk fluid elements are adjacent to boundary b?
276  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
277 
278  // Loop over the bulk fluid elements adjacent to boundary b?
279  for(unsigned e=0;e<n_element;e++)
280  {
281  // Get pointer to the bulk fluid element that is
282  // adjacent to boundary b
283  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
284  Bulk_mesh_pt->boundary_element_pt(b,e));
285 
286  //Find the index of the face of element e along boundary b
287  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
288 
289  // Create new element
290  SpineLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
291  new SpineLineVolumeConstraintBoundingElement<ELEMENT>(
292  bulk_elem_pt,face_index);
293 
294  //Set the "master" volume control element
295  el_pt->set_volume_constraint_element(vol_constraint_element);
296 
297  // Add it to the mesh
298  Volume_constraint_mesh_pt->add_element_pt(el_pt);
299  }
300  }
301 
302  //Loop over one side of the interface
303  {
304  // How many bulk fluid elements are adjacent to boundary b?
305  unsigned n_element = Bulk_mesh_pt->ninterface_lower();
306 
307  // Loop over the bulk fluid elements adjacent to boundary b?
308  for(unsigned e=0;e<n_element;e++)
309  {
310  // Get pointer to the bulk fluid element that is
311  // adjacent to boundary b
312  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
313  Bulk_mesh_pt->interface_lower_boundary_element_pt(e));
314 
315  //Find the index of the face of element e along boundary b
316  int face_index = Bulk_mesh_pt->interface_lower_face_index_at_boundary(e);
317 
318  // Create new element
319  SpineLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
320  new SpineLineVolumeConstraintBoundingElement<ELEMENT>(
321  bulk_elem_pt,face_index);
322 
323  //Set the "master" volume control element
324  el_pt->set_volume_constraint_element(vol_constraint_element);
325 
326  // Add it to the mesh
327  Volume_constraint_mesh_pt->add_element_pt(el_pt);
328  }
329  }
330 
331 }
332 
333 
334 
335 //======================================================================
336 /// Perform a parameter study
337 //======================================================================
338 template<class ELEMENT>
340 {
341 
342  // Create DocInfo object (allows checking if output directory exists)
343  DocInfo doc_info;
344  doc_info.set_directory("RESLT");
345  doc_info.number()=0;
346 
347 
348  // Open trace file
349  char filename[100];
350  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
351  Trace_file.open(filename);
352  Trace_file << "VARIABLES=\"<greek>a</greek><sub>prescribed</sub>\",";
353  Trace_file << "\"h<sub>left</sub>\",\"h<sub>right</sub>\",";
354  Trace_file << "\"<greek>a</greek><sub>left</sub>\",";
355  Trace_file << "\"<greek>a</greek><sub>right</sub>\",";
356  Trace_file << "\"p<sub>upper</sub>\",";
357  Trace_file << "\"p<sub>lower</sub>\",";
358  Trace_file << "\"p<sub>exact</sub>\"";
359  Trace_file << std::endl;
360 
361  // Doc initial state
362  doc_solution(doc_info);
363 
364  // Bump up counter
365  doc_info.number()++;
366 
367  //Gradually increase the contact angle
368  for(unsigned i=0;i<6;i++)
369  {
370  //Solve the problem
371  steady_newton_solve();
372 
373  //Output result
374  doc_solution(doc_info);
375 
376  // Bump up counter
377  doc_info.number()++;
378 
379  //Decrease the contact angle
380  Angle -= 5.0*MathematicalConstants::Pi/180.0;
381  }
382 
383 }
384 
385 
386 
387 
388 //========================================================================
389 /// Doc the solution
390 //========================================================================
391 template<class ELEMENT>
392 void CapProblem<ELEMENT>::doc_solution(DocInfo& doc_info)
393 {
394 
395  ofstream some_file;
396  char filename[100];
397 
398  // Number of plot points
399  unsigned npts=5;
400 
401 
402  //Output domain
403  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
404  doc_info.number());
405  some_file.open(filename);
406  Bulk_mesh_pt->output(some_file,npts);
407  Surface_mesh_pt->output(some_file,npts);
408  some_file.close();
409 
410 
411  // Number of interface elements
412 // unsigned ninterface=Bulk_mesh_pt->ninterface_element();
413 
414  // Number of spines
415  unsigned nspine=Bulk_mesh_pt->nspine();
416 
417  // Doc
418  Trace_file << Angle*180.0/MathematicalConstants::Pi;
419  Trace_file << " " << Bulk_mesh_pt->spine_pt(0)->height();
420  Trace_file << " " << Bulk_mesh_pt->spine_pt(nspine-1)->height();
421  Trace_file << " "
422  << dynamic_cast<ELEMENT*>(
423  Bulk_mesh_pt->upper_layer_element_pt(0))->p_nst(0);
424  Trace_file << " "
425  << dynamic_cast<ELEMENT*>(
426  Bulk_mesh_pt->lower_layer_element_pt(0))->p_nst(0);
427  Trace_file << " " << 2.0*cos(Angle)/Ca;
428  Trace_file << std::endl;
429 
430 }
431 
432 
433 
434 //==start_of_specific_mesh_class==========================================
435 /// Two layer mesh which employs a pseudo-solid node-update strategy.
436 /// This class is essentially a wrapper to an
437 /// ElasticRectangularQuadMesh, with an additional boundary
438 /// to represent the interface between the two fluid layers. In addition,
439 /// the mesh paritions the elements into those above and below
440 /// the interface and relabels boundaries so that we can impose
441 /// a volume constraint on the lower or upper fluid.
442 ///
443 /// 3
444 /// ---------------------------------------
445 /// | |
446 /// 4 | | 2
447 /// | 6 |
448 /// ---------------------------------------
449 /// | |
450 /// 5 | | 1
451 /// | |
452 /// ---------------------------------------
453 /// 0
454 //========================================================================
455 template <class ELEMENT>
457  public ElasticRectangularQuadMesh<ELEMENT>
458 {
459 
460 public:
461 
462  /// \short Constructor: Pass number of elements in x-direction, number of
463  /// elements in y-direction in bottom and top layer, respectively,
464  /// axial length and height of top and bottom layers, a boolean
465  /// flag to make the mesh periodic in the x-direction, and pointer
466  /// to timestepper (defaults to Steady timestepper)
467  ElasticTwoLayerMesh(const unsigned &nx,
468  const unsigned &ny1,
469  const unsigned &ny2,
470  const double &lx,
471  const double &h1,
472  const double &h2,
473  const bool& periodic_in_x=false,
474  TimeStepper* time_stepper_pt=
475  &Mesh::Default_TimeStepper)
476  :
477  RectangularQuadMesh<ELEMENT>(nx,ny1+ny2,lx,h1+h2,
478  periodic_in_x,time_stepper_pt),
479  ElasticRectangularQuadMesh<ELEMENT>(nx,ny1+ny2,lx,h1+h2,
480  periodic_in_x,time_stepper_pt)
481  {
482  //Set up the pointers to elements in the upper and lower fluid
483  //Calculate numbers of nodes in upper and lower regions
484  unsigned long n_lower = nx*ny1;
485  unsigned long n_upper = nx*ny2;
486  //Loop over lower elements and push back onto the array
487  Lower_layer_element_pt.resize(n_lower);
488  for(unsigned e=0;e<n_lower;e++)
489  {
490  Lower_layer_element_pt[e] = this->finite_element_pt(e);
491  }
492  //Loop over upper elements and push back onto the array
493  Upper_layer_element_pt.resize(n_upper);
494  for(unsigned e=0;e<n_upper;e++)
495  {
496  Upper_layer_element_pt[e] = this->finite_element_pt(n_lower + e);
497  }
498  //end of upper and lower fluid element assignment
499 
500  //Set the elements adjacent to the interface on both sides
501  Interface_lower_boundary_element_pt.resize(nx);
502  Interface_upper_boundary_element_pt.resize(nx);
503  {
504  unsigned count_lower=nx*(ny1-1);
505  unsigned count_upper= count_lower + nx;
506  for(unsigned e=0;e<nx;e++)
507  {
508  Interface_lower_boundary_element_pt[e] =
509  this->finite_element_pt(count_lower); ++count_lower;
510  Interface_upper_boundary_element_pt[e] =
511  this->finite_element_pt(count_upper); ++count_upper;
512  }
513  } //end of bulk elements next to interface setup
514 
515 
516  // Reset the number of boundaries
517  this->set_nboundary(7);
518 
519  //Relabel the boundary nodes
520  //Storage for the boundary coordinates that will be transferred directly
521  Vector<double> b_coord;
522  {
523  //Store the interface level
524  const double y_interface = h1 - this->Ymin;
525 
526  //Nodes on boundary 3 are now on boundaries 4 and 5
527  unsigned n_boundary_node = this->nboundary_node(3);
528  //Loop over the nodes remove them from boundary 3
529  //and add them to boundary 4 or 5 depending on their y coordinate
530  for(unsigned n=0;n<n_boundary_node;n++)
531  {
532  //Cache pointer to the node
533  Node* const nod_pt = this->boundary_node_pt(3,n);
534  //Get the boundary coordinates if set
535  if(this->boundary_coordinate_exists(3))
536  {
537  b_coord.resize(nod_pt->ncoordinates_on_boundary(3));
538  nod_pt->get_coordinates_on_boundary(3,b_coord);
539  //Indicate that the boundary coordinates are to be set on the
540  //new boundaries
541  this->Boundary_coordinate_exists[4]=true;
542  this->Boundary_coordinate_exists[5]=true;
543  }
544 
545  //Now remove the node from the old boundary
546  nod_pt->remove_from_boundary(3);
547 
548  //Find the height of the node
549  double y = nod_pt->x(1);
550  //If it's above or equal to the interface, it's on boundary 4
551  if(y >= y_interface)
552  {
553  this->add_boundary_node(4,nod_pt);
554  //Add the boundary coordinate if it has been set up
555  if(this->Boundary_coordinate_exists[4])
556  {
557  nod_pt->set_coordinates_on_boundary(4,b_coord);
558  }
559  }
560  //otherwise it's on boundary 5
561  if(y <= y_interface)
562  {
563  this->add_boundary_node(5,nod_pt);
564  //Add the boundary coordinate if it has been set up
565  if(this->Boundary_coordinate_exists[5])
566  {
567  nod_pt->set_coordinates_on_boundary(5,b_coord);
568  }
569  }
570  }
571 
572  //Now clear the boundary node information from boundary 3
573  this->Boundary_node_pt[3].clear();
574 
575  //Relabel the nodes on boundary 2 to be on boundary 3
576  n_boundary_node = this->nboundary_node(2);
577  //Loop over the nodes remove them from boundary 2
578  //and add them to boundary 3
579  for(unsigned n=0;n<n_boundary_node;n++)
580  {
581  //Cache pointer to the node
582  Node* const nod_pt = this->boundary_node_pt(2,n);
583  //Get the boundary coordinates if set
584  if(this->boundary_coordinate_exists(2))
585  {
586  b_coord.resize(nod_pt->ncoordinates_on_boundary(2));
587  nod_pt->get_coordinates_on_boundary(2,b_coord);
588  this->Boundary_coordinate_exists[3]=true;
589  }
590 
591  //Now remove the node from the boundary 2
592  nod_pt->remove_from_boundary(2);
593  //and add to boundary 3
594  this->add_boundary_node(3,nod_pt);
595  if(this->Boundary_coordinate_exists[3])
596  {
597  nod_pt->set_coordinates_on_boundary(3,b_coord);
598  }
599  }
600 
601  //Clear the information from boundary 2
602  this->Boundary_node_pt[2].clear();
603 
604  //Storage for nodes that should be removed from boundary 1
605  std::list<Node*> nodes_to_be_removed;
606 
607  //Nodes on boundary 1 are now on boundaries 1 and 2
608  n_boundary_node = this->nboundary_node(1);
609  //Loop over the nodes remove some of them from boundary 1
610  for(unsigned n=0;n<n_boundary_node;n++)
611  {
612  //Cache pointer to the node
613  Node* const nod_pt = this->boundary_node_pt(1,n);
614 
615  //Find the height of the node
616  double y = nod_pt->x(1);
617  //If it's above or equal to the interface it's on boundary 2
618  if(y >= y_interface)
619  {
620  //Get the boundary coordinates if set
621  if(this->boundary_coordinate_exists(1))
622  {
623  b_coord.resize(nod_pt->ncoordinates_on_boundary(1));
624  nod_pt->get_coordinates_on_boundary(1,b_coord);
625  this->Boundary_coordinate_exists[2]=true;
626  }
627 
628  //Now remove the node from the boundary 1 if above interace
629  if(y > y_interface)
630  {
631  nodes_to_be_removed.push_back(nod_pt);
632  }
633  //Always add to boundary 2
634  this->add_boundary_node(2,nod_pt);
635  //Add the boundary coordinate if it has been set up
636  if(this->Boundary_coordinate_exists[2])
637  {
638  nod_pt->set_coordinates_on_boundary(2,b_coord);
639  }
640  }
641  }
642 
643  //Loop over the nodes that are to be removed from boundary 1 and remove
644  //them
645  for(std::list<Node*>::iterator it=nodes_to_be_removed.begin();
646  it!=nodes_to_be_removed.end();++it)
647  {
648  this->remove_boundary_node(1,*it);
649  }
650  nodes_to_be_removed.clear();
651  } //end of boundary relabelling
652 
653  //Add the nodes to the interface
654 
655  //Read out number of linear points in the element
656  unsigned n_p = dynamic_cast<ELEMENT*>
657  (this->finite_element_pt(0))->nnode_1d();
658 
659  //Add the nodes on the interface to the boundary 6
660  //Storage for boundary coordinates (x-coordinate)
661  b_coord.resize(1);
662  this->Boundary_coordinate_exists[6];
663  //Starting index of the nodes
664  unsigned n_start=0;
665  for(unsigned e=0;e<nx;e++)
666  {
667  //If we are past the
668  if(e>0) {n_start=1;}
669  //Get the layer of elements just above the interface
670  FiniteElement* el_pt = this->finite_element_pt(nx*ny1+e);
671  //The first n_p nodes lie on the boundary
672  for(unsigned n=n_start;n<n_p;n++)
673  {
674  Node* nod_pt = el_pt->node_pt(n);
675  this->convert_to_boundary_node(nod_pt);
676  this->add_boundary_node(6,nod_pt);
677  b_coord[0] = nod_pt->x(0);
678  nod_pt->set_coordinates_on_boundary(6,b_coord);
679  }
680  }
681 
682  // Set up the boundary element information
683  this->setup_boundary_element_info();
684  }
685 
686  ///Access functions for pointers to elements in upper layer
687  FiniteElement* &upper_layer_element_pt(const unsigned long &i)
688  {return Upper_layer_element_pt[i];}
689 
690  ///Access functions for pointers to elements in bottom layer
691  FiniteElement* &lower_layer_element_pt(const unsigned long &i)
692  {return Lower_layer_element_pt[i];}
693 
694  ///Number of elements in upper layer
695  unsigned long nupper() const {return Upper_layer_element_pt.size();}
696 
697  ///Number of elements in top layer
698  unsigned long nlower() const {return Lower_layer_element_pt.size();}
699 
700  ///Access functions for pointers to elements in upper layer
701  FiniteElement* &interface_upper_boundary_element_pt(const unsigned long &i)
702  {return Interface_upper_boundary_element_pt[i];}
703 
704  ///Access functions for pointers to elements in bottom layer
705  FiniteElement* &interface_lower_boundary_element_pt(const unsigned long &i)
706  {return Interface_lower_boundary_element_pt[i];}
707 
708  ///Number of elements in upper layer
709  unsigned long ninterface_upper() const
710  {return Interface_upper_boundary_element_pt.size();}
711 
712  ///Number of elements in top layer
713  unsigned long ninterface_lower() const
714  {return Interface_lower_boundary_element_pt.size();}
715 
716  ///\short Index of the face of the elements next to the interface
717  ///in the upper region (always -2)
719  {return -2;}
720 
721  ///\short Index of the face of the elements next to the interface in
722  /// the lower region (always 2)
724  {return 2;}
725 
726 private:
727 
728  /// Vector of pointers to element in the upper layer
729  Vector <FiniteElement *> Lower_layer_element_pt;
730 
731  /// Vector of pointers to element in the lower layer
732  Vector <FiniteElement *> Upper_layer_element_pt;
733 
734  /// \short Vector of pointers to the elements adjacent to the interface
735  /// on the lower layer
736  Vector <FiniteElement*> Interface_lower_boundary_element_pt;
737 
738  /// \short Vector of pointers to the element adjacent to the interface
739  /// on the upper layer
740  Vector<FiniteElement *> Interface_upper_boundary_element_pt;
741 
742 
743 }; // End of specific mesh class
744 
745 
746 
747 //===========start_of_pseudo_elastic_class====================================
748 ///\short A class that solves the Navier--Stokes equations
749 ///to compute the shape of a static interface between two fluids in a
750 ///rectangular container with an imposed contact angle at the boundary.
751 //============================================================================
752 template<class ELEMENT>
753 class PseudoSolidCapProblem : public Problem
754 {
755 public:
756 
757  //Constructor: Arguements are the number of elements in the x
758  //direction and in each of the layers
760  const unsigned &Nx, const unsigned &Nh1, const unsigned &Nh2);
761 
762  /// Destructor: clean up memory allocated by the object
764 
765  /// Peform a parameter study: Solve problem for a range of contact angles
766  /// Pass name of output directory as a string
767  void parameter_study(const string& dir_name);
768 
769  /// Doc the solution
770  void doc_solution(DocInfo& doc_info);
771 
772 
773 private:
774 
775  /// Create the free surface elements
776  void create_free_surface_elements();
777 
778  /// Create the volume constraint elements
779  void create_volume_constraint_elements();
780 
781  /// Create the contact angle element
782  void create_contact_angle_element();
783 
784  /// The Capillary number
785  double Ca;
786 
787  /// The prescribed volume of the fluid
788  double Volume;
789 
790  /// The external pressure
791  double Pext;
792 
793  /// The contact angle
794  double Angle;
795 
796  /// Constitutive law used to determine the mesh deformation
797  ConstitutiveLaw *Constitutive_law_pt;
798 
799  // Pointer to the (single valued) Data item that
800  // will contain the pressure value that we're
801  // trading for the volume constraint
802  Data* Traded_pressure_data_pt;
803 
804  /// Trace file
805  ofstream Trace_file;
806 
807  ///Storage for the bulk mesh
809 
810  /// Storage for the free surface mesh
811  Mesh* Free_surface_mesh_pt;
812 
813  /// Storage for the element bounding the free surface
814  Mesh* Free_surface_bounding_mesh_pt;
815 
816  /// Storage for the elements that compute the enclosed volume
817  Mesh* Volume_computation_mesh_pt;
818 
819  /// Storage for the volume constraint
820  Mesh* Volume_constraint_mesh_pt;
821 
822 }; //end_of_pseudo_solid_problem_class
823 
824 
825 
826 //============start_of_constructor=====================================
827 /// Constructor: Pass boolean flag to indicate if the volume
828 /// constraint is applied by hijacking an internal pressure
829 /// or the external pressure
830 //======================================================================
831 template<class ELEMENT>
833  const unsigned &Nx, const unsigned &Nh1, const unsigned &Nh2) :
834  Ca(2.1), //Initialise value of Ca to some random value
835  Volume(0.5), //Initialise the value of the volume
836  Pext(1.23), //Initialise the external pressure to some random value
837  Angle(0.5*MathematicalConstants::Pi) //Initialise the contact angle
838 {
839  //Set the wall normal
843 
844  //Construct mesh
845  Bulk_mesh_pt = new ElasticTwoLayerMesh<ELEMENT>(Nx,Nh1,Nh2,0.5,1.0,1.0);
846 
847  //Hijack one of the pressure values in the fluid and use it
848  //as the pressure whose value is determined by the volume constraint.
849  //(Its value will affect the residual of that element but it will not
850  //be determined by it, i.e. it's hijacked).
851  Traded_pressure_data_pt = dynamic_cast<ELEMENT*>(
852  Bulk_mesh_pt->upper_layer_element_pt(0))->hijack_internal_value(0,0);
853 
854  //Set the constituive law
856  new GeneralisedHookean(&Global_Physical_Variables::Nu);
857 
858  //Loop over the elements to set the consitutive law
859  unsigned n_bulk = Bulk_mesh_pt->nelement();
860  for(unsigned e=0;e<n_bulk;e++)
861  {
862  ELEMENT* el_pt =
863  dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
864 
865  el_pt->constitutive_law_pt() = Constitutive_law_pt;
866  }
867 
868  //Set the boundary conditions
869 
870  //Fluid velocity conditions
871  //Pin the velocities on all external boundaries apart from the symmetry
872  //line boundaries 4 and 5) where only the horizontal velocity is pinned
873  for (unsigned b=0;b<6;b++)
874  {
875  //Find the number of nodes on the boundary
876  unsigned n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
877  //Loop over the nodes on the boundary
878  for(unsigned n=0;n<n_boundary_node;n++)
879  {
880  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
881  if((b!=4) && (b!=5))
882  {
883  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
884  }
885  }
886  } //end_of_fluid_boundary_conditions
887 
888  //PesudoSolid boundary conditions,
889  for(unsigned b=0;b<6;b++)
890  {
891  //Find the number of nodes on the boundary
892  unsigned n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
893  //Loop over the nodes on the boundary
894  for(unsigned n=0;n<n_boundary_node;n++)
895  {
896  //Pin vertical displacement on the bottom and top
897  if((b==0) || (b==3))
898  {
899  static_cast<SolidNode*>(Bulk_mesh_pt->boundary_node_pt(b,n))
900  ->pin_position(1);
901  }
902  //Pin horizontal displacement on the sides
903  if((b==1) || (b==2) || (b==4) || (b==5))
904  {
905  static_cast<SolidNode*>(Bulk_mesh_pt->boundary_node_pt(b,n))
906  ->pin_position(0);
907  }
908  }
909  } //end_of_solid_boundary_conditions
910 
911  //Constrain all nodes only to move vertically (not horizontally)
912  {
913  unsigned n_node = Bulk_mesh_pt->nnode();
914  for(unsigned n=0;n<n_node;n++)
915  {
916  static_cast<SolidNode*>(Bulk_mesh_pt->node_pt(n))->pin_position(0);
917  }
918  } //end_of_constraint
919 
920  // Pin a single pressure value in the lower fluid:
921  // Set the pressure dof 0 in element 0 to zero.
922  dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(0))->fix_pressure(0,0.0);
923 
924  //Create the free surface elements
926 
927  //Create the volume constraint elements
929 
930  //Need to make the bounding element
932 
933  //Now need to add all the meshes
934  this->add_sub_mesh(Bulk_mesh_pt);
935  this->add_sub_mesh(Free_surface_mesh_pt);
936  this->add_sub_mesh(Volume_computation_mesh_pt);
937  this->add_sub_mesh(Volume_constraint_mesh_pt);
938  this->add_sub_mesh(Free_surface_bounding_mesh_pt);
939 
940  //and build the global mesh
941  this->build_global_mesh();
942 
943  //Setup all the equation numbering and look-up schemes
944  cout << "Number of unknowns: " << assign_eqn_numbers() << std::endl;
945 
946 } //end_of_constructor
947 
948 
949 //==========================================================================
950 /// Destructor. Make sure to clean up all allocated memory, so that multiple
951 /// instances of the problem don't lead to excessive memory usage.
952 //==========================================================================
953 template<class ELEMENT>
955 {
956  //Delete the contact angle element
957  delete Free_surface_bounding_mesh_pt->element_pt(0);
958  Free_surface_bounding_mesh_pt->flush_element_and_node_storage();
959  delete Free_surface_bounding_mesh_pt;
960  //Delete the volume constraint mesh
961  delete Volume_constraint_mesh_pt;
962  //Delete the surface volume computation elements
963  unsigned n_element = Volume_computation_mesh_pt->nelement();
964  for(unsigned e=0;e<n_element;e++)
965  {delete Volume_computation_mesh_pt->element_pt(e);}
966  //Now flush the storage
967  Volume_computation_mesh_pt->flush_element_and_node_storage();
968  //Now delete the mesh
969  delete Volume_computation_mesh_pt;
970  //Delete the free surface elements
971  n_element = Free_surface_mesh_pt->nelement();
972  for(unsigned e=0;e<n_element;e++)
973  {delete Free_surface_mesh_pt->element_pt(e);}
974  //Now flush the storage
975  Free_surface_mesh_pt->flush_element_and_node_storage();
976  //Now delete the mesh
977  delete Free_surface_mesh_pt;
978 
979  //Delete the constitutive law
980  delete Constitutive_law_pt;
981 
982  //Then delete the bulk mesh
983  delete Bulk_mesh_pt;
984 }
985 
986 //============create_free_surface_element================================
987 /// Create the free surface elements
988 //========================================================================
989 template<class ELEMENT>
991 {
992  //Allocate storage for the free surface mesh
993  Free_surface_mesh_pt = new Mesh;
994 
995  // How many bulk fluid elements are adjacent to the interface
996  unsigned n_element = Bulk_mesh_pt->ninterface_lower();
997  //Loop over them
998  for(unsigned e=0;e<n_element;e++)
999  {
1000  // Create new element
1001  ElasticLineFluidInterfaceElement<ELEMENT>* el_pt =
1002  new ElasticLineFluidInterfaceElement<ELEMENT>(
1003  Bulk_mesh_pt->interface_lower_boundary_element_pt(e),
1004  Bulk_mesh_pt->interface_lower_face_index_at_boundary(e));
1005 
1006  // Add it to the mesh
1007  Free_surface_mesh_pt->add_element_pt(el_pt);
1008 
1009  //Add the capillary number
1010  el_pt->ca_pt() = &Ca;
1011  }
1012 
1013 }
1014 
1015 
1016 //============start_of_create_volume_constraint_elements==================
1017 /// Create the volume constraint elements
1018 //========================================================================
1019 template<class ELEMENT>
1021 {
1022  //Build the single volume constraint element
1023  Volume_constraint_mesh_pt = new Mesh;
1024  VolumeConstraintElement* vol_constraint_element =
1025  new VolumeConstraintElement(&Volume,Traded_pressure_data_pt,0);
1026  Volume_constraint_mesh_pt->add_element_pt(vol_constraint_element);
1027 
1028  //Now create the volume computation elements
1029  Volume_computation_mesh_pt = new Mesh;
1030 
1031  //Loop over lower boundaries (or just 1, why?)
1032  unsigned lower_boundaries[3]={0,1,5};
1033  //Loop over all boundaries
1034  for(unsigned i=0;i<3;i++)
1035  {
1036  //Read out the actual boundary
1037  unsigned b= lower_boundaries[i];
1038  // How many bulk fluid elements are adjacent to boundary b?
1039  unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
1040 
1041  // Loop over the bulk fluid elements adjacent to boundary b?
1042  for(unsigned e=0;e<n_element;e++)
1043  {
1044  // Get pointer to the bulk fluid element that is
1045  // adjacent to boundary b
1046  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
1047  Bulk_mesh_pt->boundary_element_pt(b,e));
1048 
1049  //Find the index of the face of element e along boundary b
1050  int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
1051 
1052  // Create new element
1053  ElasticLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
1054  new ElasticLineVolumeConstraintBoundingElement<ELEMENT>(
1055  bulk_elem_pt,face_index);
1056 
1057  //Set the "master" volume control element
1058  el_pt->set_volume_constraint_element(vol_constraint_element);
1059 
1060  // Add it to the mesh
1061  Volume_computation_mesh_pt->add_element_pt(el_pt);
1062  }
1063  }
1064 
1065 
1066  //Loop over one side of the interface
1067  {
1068  // How many bulk fluid elements are adjacent to boundary b?
1069  unsigned n_element = Bulk_mesh_pt->ninterface_lower();
1070 
1071  // Loop over the bulk fluid elements adjacent to boundary b?
1072  for(unsigned e=0;e<n_element;e++)
1073  {
1074  // Get pointer to the bulk fluid element that is
1075  // adjacent to boundary b
1076  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
1077  Bulk_mesh_pt->interface_lower_boundary_element_pt(e));
1078 
1079  //Find the index of the face of element e along boundary b
1080  int face_index = Bulk_mesh_pt->interface_lower_face_index_at_boundary(e);
1081 
1082  // Create new element
1083  ElasticLineVolumeConstraintBoundingElement<ELEMENT>* el_pt =
1084  new ElasticLineVolumeConstraintBoundingElement<ELEMENT>(
1085  bulk_elem_pt,face_index);
1086 
1087  //Set the "master" volume control element
1088  el_pt->set_volume_constraint_element(vol_constraint_element);
1089 
1090  // Add it to the mesh
1091  Volume_constraint_mesh_pt->add_element_pt(el_pt);
1092  }
1093  }
1094 
1095 
1096 } //end_of_create_volume_constraint_elements
1097 
1098 //==========start_of_create_contact_angle_elements========================
1099 /// Create the contact angle element
1100 //========================================================================
1101 template<class ELEMENT>
1103 {
1104  Free_surface_bounding_mesh_pt = new Mesh;
1105 
1106  //Find the element at the end of the free surface
1107  //The elements are assigned in order of increasing x coordinate
1108  unsigned n_free_surface = Free_surface_mesh_pt->nelement();
1109 
1110  //Make the bounding element for the contact angle constraint
1111  //which works because the order of elements in the mesh is known
1112  FluidInterfaceBoundingElement* el_pt =
1113  dynamic_cast<ElasticLineFluidInterfaceElement<ELEMENT>*>
1114  (Free_surface_mesh_pt->element_pt(n_free_surface-1))->
1115  make_bounding_element(1);
1116 
1117  //Set the contact angle (strong imposition)
1118  el_pt->set_contact_angle(&Angle);
1119 
1120  //Set the capillary number
1121  el_pt->ca_pt() = &Ca;
1122 
1123  //Set the wall normal of the external boundary
1124  el_pt->wall_unit_normal_fct_pt()
1126 
1127  //Add the element to the mesh
1128  Free_surface_bounding_mesh_pt->add_element_pt(el_pt);
1129 
1130 } //end_of_create_contact_angle_element
1131 
1132 
1133 
1134 //================start_of_parameter_study===========================
1135 /// Perform a parameter study. Pass name of output directory as
1136 /// a string
1137 //======================================================================
1138 template<class ELEMENT>
1139 void PseudoSolidCapProblem<ELEMENT>::parameter_study(const string& dir_name)
1140 {
1141  // Create DocInfo object (allows checking if output directory exists)
1142  DocInfo doc_info;
1143  doc_info.set_directory(dir_name);
1144  doc_info.number()=0;
1145 
1146  // Open trace file
1147  char filename[100];
1148  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
1149  Trace_file.open(filename);
1150  Trace_file << "VARIABLES=\"<greek>a</greek><sub>prescribed</sub>\",";
1151  Trace_file << "\"h<sub>left</sub>\",\"h<sub>right</sub>\",";
1152  Trace_file << "\"p<sub>fluid</sub>-p<sub>ext</sub>\",";
1153  Trace_file << "\"<greek>D</greek>p<sub>exact</sub>\"";
1154  Trace_file << std::endl;
1155 
1156  // Doc initial state
1157  doc_solution(doc_info);
1158 
1159  // Bump up counter
1160  doc_info.number()++;
1161 
1162  //Solve the problem for six different contact angles
1163  for(unsigned i=0;i<6;i++)
1164  {
1165  //Solve the problem
1166  steady_newton_solve();
1167 
1168  //Output result
1169  doc_solution(doc_info);
1170 
1171  // Bump up counter
1172  doc_info.number()++;
1173 
1174  //Decrease the contact angle
1175  Angle -= 5.0*MathematicalConstants::Pi/180.0;
1176  }
1177 
1178 }
1179 
1180 
1181 
1182 
1183 //==============start_of_doc_solution=====================================
1184 /// Doc the solution
1185 //========================================================================
1186 template<class ELEMENT>
1187 void PseudoSolidCapProblem<ELEMENT>::doc_solution(DocInfo& doc_info)
1188 {
1189  //Output stream
1190  ofstream some_file;
1191  char filename[100];
1192 
1193  // Number of plot points
1194  unsigned npts=5;
1195 
1196  //Output domain
1197  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
1198  doc_info.number());
1199  some_file.open(filename);
1200  Bulk_mesh_pt->output(some_file,npts);
1201  some_file.close();
1202 
1203 
1204  // Number of interface elements
1205  unsigned ninterface=Free_surface_mesh_pt->nelement();
1206  //Find number of nodes in the last interface element
1207  unsigned np = Free_surface_mesh_pt->finite_element_pt(ninterface-1)->nnode();
1208  // Document the contact angle (in degrees), the height of the interface at
1209  // the centre of the container, the height at the wall, the computed
1210  // pressure drop across the interface and
1211  // the analytic prediction of the pressure drop.
1212  Trace_file << Angle*180.0/MathematicalConstants::Pi;
1213  Trace_file << " " << Free_surface_mesh_pt->finite_element_pt(0)
1214  ->node_pt(0)->x(1)
1215  << " "
1216  << Free_surface_mesh_pt->finite_element_pt(ninterface-1)
1217  ->node_pt(np-1)->x(1);
1218  Trace_file << " "
1219  << dynamic_cast<ELEMENT*>(
1220  Bulk_mesh_pt->upper_layer_element_pt(0))->p_nst(0);
1221  Trace_file << " "
1222  << dynamic_cast<ELEMENT*>(
1223  Bulk_mesh_pt->lower_layer_element_pt(0))->p_nst(0);
1224  Trace_file << " " << 2.0*cos(Angle)/Ca;
1225  Trace_file << std::endl;
1226 
1227 } //end_of_doc_solution
1228 
1229 
1230 //======================================================================
1231 ///Main driver: Build problem and initiate parameter study
1232 //======================================================================
1233 int main()
1234 {
1235  {
1236  //Construct the problem with 4 x (4+4) elements
1238 
1239  //Solve the problem :)
1240  problem.parameter_study();
1241  }
1242 
1243  {
1244  //Construct the problem with 4 x (4+4) elements
1245  PseudoSolidCapProblem<Hijacked<
1246  PseudoSolidNodeUpdateElement<QCrouzeixRaviartElement<2>,
1247  QPVDElementWithPressure<2> > > >
1248  problem(4,4,4);
1249 
1250  //Solve the problem :)
1251  problem.parameter_study("RESLT_elastic");
1252  }
1253 
1254 }
void parameter_study(const string &dir_name)
unsigned long ninterface_lower() const
Number of elements in top layer.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
TwoLayerSpineMesh< SpineElement< ELEMENT > > * Bulk_mesh_pt
Pointer to the mesh.
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.
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.
ElasticTwoLayerMesh< ELEMENT > * Bulk_mesh_pt
Storage for the bulk mesh.
Vector< FiniteElement * > Lower_layer_element_pt
Vector of pointers to element in the upper layer.
unsigned long ninterface_upper() const
Number of elements in upper layer.
FiniteElement *& upper_layer_element_pt(const unsigned long &i)
Access functions for pointers to elements in upper layer.
void parameter_study(const string &dir_name)
Mesh * Free_surface_bounding_mesh_pt
Storage for the element bounding the free surface.
void create_volume_constraint_elements()
Create the volume constraint elements.
int main()
Main driver: Build problem and initiate parameter study.
Mesh * Free_surface_mesh_pt
Storage for the free surface mesh.
void create_volume_constraint_elements()
Create the volume constraint elements.
int interface_lower_face_index_at_boundary(const unsigned &e)
Index of the face of the elements next to the interface in the lower region (always 2) ...
unsigned long nlower() const
Number of elements in top layer.
Vector< FiniteElement * > Upper_layer_element_pt
Vector of pointers to element in the lower layer.
FiniteElement *& interface_upper_boundary_element_pt(const unsigned long &i)
Access functions for pointers to elements in upper layer.
FiniteElement *& interface_lower_boundary_element_pt(const unsigned long &i)
Access functions for pointers to elements in bottom layer.
void wall_unit_normal_fct(const Vector< double > &x, Vector< double > &normal)
Function that specifies the wall unit normal.
int interface_upper_face_index_at_boundary(const unsigned &e)
Index of the face of the elements next to the interface in the upper region (always -2) ...
Mesh * Volume_constraint_mesh_pt
Storage for the volume constraint.
CapProblem(const bool &hijack_internal)
~PseudoSolidCapProblem()
Destructor: clean up memory allocated by the object.
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...
Vector< FiniteElement * > Interface_upper_boundary_element_pt
Vector of pointers to the element adjacent to the interface on the upper layer.
void doc_solution(DocInfo &doc_info)
Doc the solution.
Vector< FiniteElement * > Interface_lower_boundary_element_pt
Vector of pointers to the elements adjacent to the interface on the lower layer.
Mesh * Volume_computation_mesh_pt
Storage for the elements that compute the enclosed volume.
ElasticTwoLayerMesh(const unsigned &nx, const unsigned &ny1, const unsigned &ny2, const double &lx, const double &h1, const double &h2, const bool &periodic_in_x=false, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction, number of elements in y-direction in bottom and ...
unsigned long nupper() const
Number of elements in upper layer.
double Nu
Pseudo-solid Poisson ratio.
Mesh * Bulk_mesh_pt
Storage for the bulk mesh.
FiniteElement *& lower_layer_element_pt(const unsigned long &i)
Access functions for pointers to elements in bottom layer.