elastic_two_layer_interface.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 for an adaptive two-dimensional two fluid interface problem,
31 // where the mesh is deformed using a pseudo-solid node-update strategy
32 
33 // Generic oomph-lib header
34 #include "generic.h"
35 
36 // Navier-Stokes headers
37 #include "navier_stokes.h"
38 
39 // Interface headers
40 #include "fluid_interface.h"
41 
42 // Constitutive law headers
43 #include "constitutive.h"
44 
45 // Solid headers
46 #include "solid.h"
47 
48 // The mesh
49 #include "meshes/rectangular_quadmesh.h"
50 
51 using namespace std;
52 
53 using namespace oomph;
54 
55 
56 //==start_of_namespace====================================================
57 /// Namespace for physical parameters
58 //========================================================================
59 namespace Global_Physical_Variables
60 {
61 
62  /// Reynolds number
63  double Re = 5.0;
64 
65  /// Strouhal number
66  double St = 1.0;
67 
68  /// Womersley number (Reynolds x Strouhal)
69  double ReSt = 5.0;
70 
71  /// Product of Reynolds number and inverse of Froude number
72  double ReInvFr = 5.0;
73 
74  /// \short Ratio of viscosity in upper fluid to viscosity in lower
75  /// fluid. Reynolds number etc. is based on viscosity in lower fluid.
76  double Viscosity_Ratio = 0.1;
77 
78  /// \short Ratio of density in upper fluid to density in lower
79  /// fluid. Reynolds number etc. is based on density in lower fluid.
80  double Density_Ratio = 0.5;
81 
82  /// Capillary number
83  double Ca = 0.01;
84 
85  /// Direction of gravity
86  Vector<double> G(2);
87 
88  /// Pseudo-solid Poisson ratio
89  double Nu = 0.1;
90 
91 } // End of namespace
92 
93 
94 //////////////////////////////////////////////////////////////////////////
95 //////////////////////////////////////////////////////////////////////////
96 //////////////////////////////////////////////////////////////////////////
97 
98 
99 //==start_of_specific_mesh_class==========================================
100 /// Two layer mesh which employs a pseudo-solid node-update strategy.
101 /// This class is essentially a wrapper to an
102 /// ElasticRefineableRectangularQuadMesh, with an additional boundary
103 /// to represent the interface between the two fluid layers.
104 //========================================================================
105 template <class ELEMENT>
107  public virtual ElasticRefineableRectangularQuadMesh<ELEMENT>
108 {
109 
110 public:
111 
112  /// \short Constructor: Pass number of elements in x-direction, number of
113  /// elements in y-direction in bottom and top layer, respectively,
114  /// axial length and height of top and bottom layers, a boolean
115  /// flag to make the mesh periodic in the x-direction, and pointer
116  /// to timestepper (defaults to Steady timestepper)
117  ElasticRefineableTwoLayerMesh(const unsigned &nx,
118  const unsigned &ny1,
119  const unsigned &ny2,
120  const double &lx,
121  const double &h1,
122  const double &h2,
123  const bool& periodic_in_x,
124  TimeStepper* time_stepper_pt=
125  &Mesh::Default_TimeStepper)
126  : RectangularQuadMesh<ELEMENT>(nx,ny1+ny2,lx,h1+h2,
127  periodic_in_x,time_stepper_pt),
128  ElasticRectangularQuadMesh<ELEMENT>(nx,ny1+ny2,lx,h1+h2,
129  periodic_in_x,time_stepper_pt),
130  ElasticRefineableRectangularQuadMesh<ELEMENT>(nx,ny1+ny2,lx,h1+h2,
131  periodic_in_x,
132  time_stepper_pt)
133  {
134  // ----------------------------------------------------
135  // Convert all nodes on the interface to boundary nodes
136  // ----------------------------------------------------
137 
138  // Set the number of boundaries to 5
139  this->set_nboundary(5);
140 
141  // Loop over horizontal elements
142  for(unsigned e=0;e<nx;e++)
143  {
144  // Get pointer to element in lower fluid adjacent to interface
145  FiniteElement* el_pt = this->finite_element_pt(nx*(ny1-1)+e);
146 
147  // Determine number of nodes in this element
148  const unsigned n_node = el_pt->nnode();
149 
150  // The last three nodes in this element are those on the interface.
151  // Loop over these nodes and convert them to boundary nodes.
152  for(unsigned n=0;n<3;n++)
153  {
154  Node* nod_pt = el_pt->node_pt(n_node-3+n);
155  this->convert_to_boundary_node(nod_pt);
156  this->add_boundary_node(4,nod_pt);
157  }
158  } // End of loop over horizontal elements
159 
160  // Set up the boundary element information
161  this->setup_boundary_element_info();
162  }
163 
164 }; // End of specific mesh class
165 
166 
167 //////////////////////////////////////////////////////////////////////////
168 //////////////////////////////////////////////////////////////////////////
169 //////////////////////////////////////////////////////////////////////////
170 
171 
172 //==start_of_problem_class================================================
173 /// Two fluid interface problem in a rectangular domain which is
174 /// periodic in the x direction
175 //========================================================================
176 template <class ELEMENT, class TIMESTEPPER>
177 class InterfaceProblem : public Problem
178 {
179 
180 public:
181 
182  /// Constructor
184 
185  /// Destructor (empty)
187 
188  /// Set initial conditions
189  void set_initial_condition();
190 
191  /// Set boundary conditions
192  void set_boundary_conditions();
193 
194  /// Document the solution
195  void doc_solution(DocInfo &doc_info);
196 
197  /// Do unsteady run up to maximum time t_max with given timestep dt
198  void unsteady_run(const double &t_max, const double &dt);
199 
200 private:
201 
202  /// No actions required before solve step
204 
205  /// No actions required after solve step
207 
208  /// \short Actions before the timestep: For maximum stability, reset
209  /// the current nodal positions to be the "stress-free" ones.
211  {
212  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
213  }
214 
215  /// Strip off the interface elements before adapting the bulk mesh
216  void actions_before_adapt();
217 
218  /// Rebuild the mesh of interface elements after adapting the bulk mesh
219  void actions_after_adapt();
220 
221  /// Create the 1d interface elements
222  void create_interface_elements();
223 
224  /// Delete the 1d interface elements
225  void delete_interface_elements();
226 
227  /// Deform the mesh/free surface to a prescribed function
228  void deform_free_surface(const double &epsilon, const unsigned &n_periods);
229 
230  /// Fix pressure in element e at pressure dof pdof and set to pvalue
231  void fix_pressure(const unsigned &e,
232  const unsigned &pdof,
233  const double &pvalue)
234  {
235  // Fix the pressure at that element
236  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
237  fix_pressure(pdof,pvalue);
238  }
239 
240  /// Pointer to the (specific) "bulk" mesh
242 
243  /// Pointer to the "surface" mesh
245 
246  // Pointer to the constitutive law used to determine the mesh deformation
247  ConstitutiveLaw* Constitutive_law_pt;
248 
249  /// Width of domain
250  double Lx;
251 
252  /// Trace file
253  ofstream Trace_file;
254 
255 }; // End of problem class
256 
257 
258 
259 //==start_of_constructor==================================================
260 /// Constructor for two fluid interface problem
261 //========================================================================
262 template <class ELEMENT, class TIMESTEPPER>
265 {
266  // Allocate the timestepper (this constructs the time object as well)
267  add_time_stepper_pt(new TIMESTEPPER);
268 
269  // Define number of elements in x direction
270  const unsigned n_x = 3;
271 
272  // Define number of elements in y direction in lower fluid (fluid 1)
273  const unsigned n_y1 = 3;
274 
275  // Define number of elements in y direction in upper fluid (fluid 2)
276  const unsigned n_y2 = 3;
277 
278  // Define width of domain and store as class member data
279  const double l_x = 1.0;
280  this->Lx = l_x;
281 
282  // Define height of lower fluid layer
283  const double h1 = 1.0;
284 
285  // Define height of upper fluid layer
286  const double h2 = 1.0;
287 
288  // Build and assign the "bulk" mesh (the "true" boolean flag tells
289  // the mesh constructor that the domain is periodic in x)
290  Bulk_mesh_pt = new ElasticRefineableTwoLayerMesh<ELEMENT>
291  (n_x,n_y1,n_y2,l_x,h1,h2,true,time_stepper_pt());
292 
293  // Create and set the error estimator for spatial adaptivity
294  Bulk_mesh_pt->spatial_error_estimator_pt() = new Z2ErrorEstimator;
295 
296  // Set the maximum refinement level for the mesh to 4
297  Bulk_mesh_pt->max_refinement_level() = 4;
298 
299  // Create the "surface" mesh that will contain only the interface
300  // elements. The constructor just creates the mesh without giving
301  // it any elements, nodes, etc.
302  Surface_mesh_pt = new Mesh;
303 
304  // Create interface elements at the boundary between the two fluids,
305  // and add them to the surface mesh
306  create_interface_elements();
307 
308  // Add the two sub meshes to the problem
309  add_sub_mesh(Bulk_mesh_pt);
310  add_sub_mesh(Surface_mesh_pt);
311 
312  // Combine all sub-meshes into a single mesh
313  build_global_mesh();
314 
315  // --------------------------------------------
316  // Set the boundary conditions for this problem
317  // --------------------------------------------
318 
319  // All values are free by default -- just pin the ones that have
320  // Dirichlet conditions here
321 
322  // Determine number of mesh boundaries
323  const unsigned n_boundary = Bulk_mesh_pt->nboundary();
324 
325  // Loop over mesh boundaries
326  for(unsigned b=0;b<n_boundary;b++)
327  {
328  // Determine number of nodes on boundary b
329  const unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
330 
331  // Loop over nodes on boundary b
332  for(unsigned n=0;n<n_node;n++)
333  {
334  // Fluid boundary conditions:
335  // --------------------------
336 
337  // Pin x-component of velocity (no slip/penetration)
338  // on all boundaries other than the interface (b=4)
339  if(b!=4) { Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0); }
340 
341  // Pin y-component of velocity on both solid boundaries (no penetration)
342  if(b==0 || b==2) { Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1); }
343 
344  // Solid boundary conditions:
345  // --------------------------
346 
347  // Pin vertical mesh position on both solid boundaries (no penetration)
348  if(b==0 || b==2) { Bulk_mesh_pt->boundary_node_pt(b,n)->pin_position(1); }
349 
350  } // End of loop over nodes on boundary b
351  } // End of loop over mesh boundaries
352 
353  // Pin horizontal position of all nodes
354  const unsigned n_node = Bulk_mesh_pt->nnode();
355  for(unsigned n=0;n<n_node;n++) { Bulk_mesh_pt->node_pt(n)->pin_position(0); }
356 
357  // Define a constitutive law for the solid equations: generalised Hookean
358  Constitutive_law_pt = new GeneralisedHookean(&Global_Physical_Variables::Nu);
359 
360  // ----------------------------------------------------------------
361  // Complete the problem setup to make the elements fully functional
362  // ----------------------------------------------------------------
363 
364  // Compute number of bulk elements in lower/upper fluids
365  const unsigned n_lower = n_x*n_y1;
366  const unsigned n_upper = n_x*n_y2;
367 
368  // Loop over bulk elements in lower fluid
369  for(unsigned e=0;e<n_lower;e++)
370  {
371  // Upcast from GeneralisedElement to the present element
372  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
373 
374  // Set the Reynolds number
375  el_pt->re_pt() = &Global_Physical_Variables::Re;
376 
377  // Set the Womersley number
378  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
379 
380  // Set the product of the Reynolds number and the inverse of the
381  // Froude number
382  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
383 
384  // Set the direction of gravity
385  el_pt->g_pt() = &Global_Physical_Variables::G;
386 
387  // Set the constitutive law
388  el_pt->constitutive_law_pt() = Constitutive_law_pt;
389 
390  } // End of loop over bulk elements in lower fluid
391 
392  // Loop over bulk elements in upper fluid
393  for(unsigned e=n_lower;e<(n_lower+n_upper);e++)
394  {
395  // Upcast from GeneralisedElement to the present element
396  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
397 
398  // Set the Reynolds number
399  el_pt->re_pt() = &Global_Physical_Variables::Re;
400 
401  // Set the Womersley number
402  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
403 
404  // Set the product of the Reynolds number and the inverse of the
405  // Froude number
406  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
407 
408  // Set the viscosity ratio
409  el_pt->viscosity_ratio_pt() = &Global_Physical_Variables::Viscosity_Ratio;
410 
411  // Set the density ratio
412  el_pt->density_ratio_pt() = &Global_Physical_Variables::Density_Ratio;
413 
414  // Set the direction of gravity
415  el_pt->g_pt() = &Global_Physical_Variables::G;
416 
417  // Set the constitutive law
418  el_pt->constitutive_law_pt() = Constitutive_law_pt;
419 
420  } // End of loop over bulk elements in upper fluid
421 
422  // Set the pressure in the first element at 'node' 0 to 0.0
423  fix_pressure(0,0,0.0);
424 
425  // Apply the boundary conditions
426  set_boundary_conditions();
427 
428  // Set up equation numbering scheme
429  cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
430 
431 } // End of constructor
432 
433 
434 
435 //==start_of_set_initial_condition========================================
436 /// \short Set initial conditions: Set all nodal velocities to zero and
437 /// initialise the previous velocities and nodal positions to correspond
438 /// to an impulsive start
439 //========================================================================
440 template <class ELEMENT, class TIMESTEPPER>
442 {
443  // Determine number of nodes in mesh
444  const unsigned n_node = mesh_pt()->nnode();
445 
446  // Loop over all nodes in mesh
447  for(unsigned n=0;n<n_node;n++)
448  {
449  // Loop over the two velocity components
450  for(unsigned i=0;i<2;i++)
451  {
452  // Set velocity component i of node n to zero
453  mesh_pt()->node_pt(n)->set_value(i,0.0);
454  }
455  }
456 
457  // Initialise the previous velocity values and nodal positions
458  // for timestepping corresponding to an impulsive start
459  assign_initial_values_impulsive();
460 
461 } // End of set_initial_condition
462 
463 
464 
465 //==start_of_set_boundary_conditions======================================
466 /// \short Set boundary conditions: Set both velocity components
467 /// to zero on the top and bottom (solid) walls and the horizontal
468 /// component only to zero on the side (periodic) boundaries
469 //========================================================================
470 template <class ELEMENT, class TIMESTEPPER>
472 {
473  // Determine number of mesh boundaries
474  const unsigned n_boundary = Bulk_mesh_pt->nboundary();
475 
476  // Loop over mesh boundaries
477  for(unsigned b=0;b<n_boundary;b++)
478  {
479  // Determine number of nodes on boundary b
480  const unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
481 
482  // Loop over nodes on boundary b
483  for(unsigned n=0;n<n_node;n++)
484  {
485  // Set x-component of the velocity to zero on all boundaries
486  // other than the interface (b=4)
487  if(b!=4) { Bulk_mesh_pt->boundary_node_pt(b,n)->set_value(0,0.0); }
488 
489  // Set y-component of the velocity to zero on both solid
490  // boundaries (top and bottom)
491  if(b==0 || b==2)
492  {
493  Bulk_mesh_pt->boundary_node_pt(b,n)->set_value(1,0.0);
494  }
495  } // End of loop over nodes on boundary b
496  } // End of loop over mesh boundaries
497 
498 } // End of set_boundary_conditions
499 
500 
501 
502 //==start_of_actions_before_adapt=========================================
503 /// Strip off the interface elements before adapting the bulk mesh
504 //========================================================================
505 template<class ELEMENT, class TIMESTEPPER>
507 {
508  // Delete the interface elements and wipe the surface mesh
509  delete_interface_elements();
510 
511  // Rebuild the Problem's global mesh from its various sub-meshes
512  rebuild_global_mesh();
513 
514 } // End of actions_before_adapt
515 
516 
517 
518 //==start_of_actions_after_adapt==========================================
519 /// Rebuild the mesh of interface elements after adapting the bulk mesh
520 //========================================================================
521 template<class ELEMENT, class TIMESTEPPER>
523 {
524  // Create the interface elements
525  this->create_interface_elements();
526 
527  // Rebuild the Problem's global mesh from its various sub-meshes
528  rebuild_global_mesh();
529 
530  // Pin horizontal displacement of all nodes
531  const unsigned n_node = Bulk_mesh_pt->nnode();
532  for(unsigned n=0;n<n_node;n++) { Bulk_mesh_pt->node_pt(n)->pin_position(0); }
533 
534  // Unpin all fluid pressure dofs
535  RefineableNavierStokesEquations<2>::
536  unpin_all_pressure_dofs(Bulk_mesh_pt->element_pt());
537 
538  // Pin redudant fluid pressure dofs
539  RefineableNavierStokesEquations<2>::
540  pin_redundant_nodal_pressures(Bulk_mesh_pt->element_pt());
541 
542  // Now set the pressure in the first element at 'node' 0 to 0.0
543  fix_pressure(0,0,0.0);
544 
545  // Pin the redundant solid pressures (if any)
546  PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
547  Bulk_mesh_pt->element_pt());
548 
549  // Reset the boundary conditions
550  set_boundary_conditions();
551 
552 } // End of actions_after_adapt
553 
554 
555 
556 //==start_of_create_interface_elements====================================
557 /// \short Create interface elements between the two fluids in the mesh
558 /// pointed to by Bulk_mesh_pt and add the elements to the Mesh object
559 /// pointed to by Surface_mesh_pt.
560 //========================================================================
561 template <class ELEMENT, class TIMESTEPPER>
563 {
564  // Determine number of bulk elements adjacent to interface (boundary 4)
565  const unsigned n_element = this->Bulk_mesh_pt->nboundary_element(4);
566 
567  // Loop over those elements adjacent to the interface
568  for(unsigned e=0;e<n_element;e++)
569  {
570  // Get pointer to the bulk element that is adjacent to the interface
571  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
572  this->Bulk_mesh_pt->boundary_element_pt(4,e));
573 
574  // We only want to attach interface elements to the bulk elements
575  // which are BELOW the interface, and so we filter out those above by
576  // referring to the viscosity_ratio_pt
577  if(bulk_elem_pt->viscosity_ratio_pt()
579  {
580  // Find index of the face of element e that corresponds to the interface
581  const int face_index = this->Bulk_mesh_pt->face_index_at_boundary(4,e);
582 
583  // Create the interface element
584  FiniteElement* interface_element_element_pt =
585  new ElasticLineFluidInterfaceElement<ELEMENT>(bulk_elem_pt,face_index);
586 
587  // Add the interface element to the surface mesh
588  this->Surface_mesh_pt->add_element_pt(interface_element_element_pt);
589  }
590  }
591 
592  // --------------------------------------------------------
593  // Complete the setup to make the elements fully functional
594  // --------------------------------------------------------
595 
596  // Determine number of 1D interface elements in mesh
597  const unsigned n_interface_element = this->Surface_mesh_pt->nelement();
598 
599  // Loop over the interface elements
600  for(unsigned e=0;e<n_interface_element;e++)
601  {
602  // Upcast from GeneralisedElement to the present element
603  ElasticLineFluidInterfaceElement<ELEMENT>* el_pt =
604  dynamic_cast<ElasticLineFluidInterfaceElement<ELEMENT>*>
605  (Surface_mesh_pt->element_pt(e));
606 
607  // Set the Strouhal number
608  el_pt->st_pt() = &Global_Physical_Variables::St;
609 
610  // Set the Capillary number
611  el_pt->ca_pt() = &Global_Physical_Variables::Ca;
612 
613  } // End of loop over interface elements
614 
615 } // End of create_interface_elements
616 
617 
618 
619 //==start_of_delete_interface_elements====================================
620 /// Delete the interface elements and wipe the surface mesh
621 //========================================================================
622 template <class ELEMENT, class TIMESTEPPER>
624 {
625  // Determine number of interface elements
626  const unsigned n_interface_element = Surface_mesh_pt->nelement();
627 
628  // Loop over interface elements and delete
629  for(unsigned e=0;e<n_interface_element;e++)
630  {
631  delete Surface_mesh_pt->element_pt(e);
632  }
633 
634  // Wipe the mesh
635  Surface_mesh_pt->flush_element_and_node_storage();
636 
637 } // End of delete_interface_elements
638 
639 
640 
641 //==start_of_deform_free_surface==========================================
642 /// Deform the mesh/free surface to a prescribed function
643 //========================================================================
644 template <class ELEMENT, class TIMESTEPPER>
646 deform_free_surface(const double &epsilon,const unsigned &n_periods)
647 {
648  // Determine number of nodes in the "bulk" mesh
649  const unsigned n_node = Bulk_mesh_pt->nnode();
650 
651  // Loop over all nodes in mesh
652  for(unsigned n=0;n<n_node;n++)
653  {
654  // Determine eulerian position of node
655  const double current_x_pos = Bulk_mesh_pt->node_pt(n)->x(0);
656  const double current_y_pos = Bulk_mesh_pt->node_pt(n)->x(1);
657 
658  // Determine new vertical position of node
659  const double new_y_pos = current_y_pos
660  + (1.0-fabs(1.0-current_y_pos))*epsilon
661  *(cos(2.0*n_periods*MathematicalConstants::Pi*current_x_pos/Lx));
662 
663  // Set new position
664  Bulk_mesh_pt->node_pt(n)->x(1) = new_y_pos;
665  }
666 } // End of deform_free_surface
667 
668 
669 
670 //==start_of_doc_solution=================================================
671 /// Doc the solution
672 //========================================================================
673 template <class ELEMENT, class TIMESTEPPER>
675 {
676 
677  // Output the time
678  cout << "Time is now " << time_pt()->time() << std::endl;
679 
680  // Upcast from GeneralisedElement to the present element
681  ElasticLineFluidInterfaceElement<ELEMENT>* el_pt =
682  dynamic_cast<ElasticLineFluidInterfaceElement<ELEMENT>*>
683  (Surface_mesh_pt->element_pt(0));
684 
685  // Document time and vertical position of left hand side of interface
686  // in trace file
687  Trace_file << time_pt()->time() << " "
688  << el_pt->node_pt(0)->x(1) << std::endl;
689 
690  ofstream some_file;
691  char filename[100];
692 
693  // Set number of plot points (in each coordinate direction)
694  const unsigned npts = 5;
695 
696  // Open solution output file
697  sprintf(filename,"%s/soln%i.dat",
698  doc_info.directory().c_str(),doc_info.number());
699  some_file.open(filename);
700 
701  // Output solution to file
702  Bulk_mesh_pt->output(some_file,npts);
703 
704  // Close solution output file
705  some_file.close();
706 
707  // Open interface solution output file
708  sprintf(filename,"%s/interface_soln%i.dat",
709  doc_info.directory().c_str(),doc_info.number());
710  some_file.open(filename);
711 
712  // Output solution to file
713  Surface_mesh_pt->output(some_file,npts);
714 
715  // Close solution output file
716  some_file.close();
717 
718 } // End of doc_solution
719 
720 
721 
722 //==start_of_unsteady_run=================================================
723 /// Perform run up to specified time t_max with given timestep dt
724 //========================================================================
725 template <class ELEMENT, class TIMESTEPPER>
727 unsteady_run(const double &t_max, const double &dt)
728 {
729 
730  // Set value of epsilon
731  const double epsilon = 0.1;
732 
733  // Set number of periods for cosine term
734  const unsigned n_periods = 1;
735 
736  // Deform the mesh/free surface
737  deform_free_surface(epsilon,n_periods);
738 
739  // Initialise DocInfo object
740  DocInfo doc_info;
741 
742  // Set output directory
743  doc_info.set_directory("RESLT");
744 
745  // Initialise counter for solutions
746  doc_info.number()=0;
747 
748  // Open trace file
749  char filename[100];
750  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
751  Trace_file.open(filename);
752 
753  // Initialise trace file
754  Trace_file << "time, interface height" << std::endl;
755 
756  // Initialise timestep
757  initialise_dt(dt);
758 
759  // Set initial condition
760  set_initial_condition();
761 
762  // Maximum number of spatial adaptations per timestep
763  unsigned max_adapt = 2;
764 
765  // Call refine_uniformly twice
766  for(unsigned i=0;i<2;i++) { refine_uniformly(); }
767 
768  // Doc initial solution
769  doc_solution(doc_info);
770 
771  // Increment counter for solutions
772  doc_info.number()++;
773 
774  // Determine number of timesteps
775  const unsigned n_timestep = unsigned(t_max/dt);
776 
777  // Are we on the first timestep? At this point, yes!
778  bool first_timestep = true;
779 
780  // Timestepping loop
781  for(unsigned t=1;t<=n_timestep;t++)
782  {
783  // Output current timestep to screen
784  cout << "\nTimestep " << t << " of " << n_timestep << std::endl;
785 
786  // Take one fixed timestep with spatial adaptivity
787  unsteady_newton_solve(dt,max_adapt,first_timestep);
788 
789  // No longer on first timestep, so set first_timestep flag to false
790  first_timestep = false;
791 
792  // Reset maximum number of adaptations for all future timesteps
793  max_adapt = 1;
794 
795  // Doc solution
796  doc_solution(doc_info);
797 
798  // Increment counter for solutions
799  doc_info.number()++;
800 
801  } // End of timestepping loop
802 
803 } // End of unsteady_run
804 
805 
806 //////////////////////////////////////////////////////////////////////////
807 //////////////////////////////////////////////////////////////////////////
808 //////////////////////////////////////////////////////////////////////////
809 
810 
811 //==start_of_main=========================================================
812 /// Driver code for two-dimensional two fluid interface problem
813 //========================================================================
814 int main(int argc, char* argv[])
815 {
816  // Store command line arguments
817  CommandLineArgs::setup(argc,argv);
818 
819  // -----------------------------------------------------------------
820  // Define possible command line arguments and parse the ones that
821  // were actually specified
822  // -----------------------------------------------------------------
823 
824  // Are we performing a validation run?
825  CommandLineArgs::specify_command_line_flag("--validation");
826 
827  // Parse command line
828  CommandLineArgs::parse_and_assign();
829 
830  // Doc what has actually been specified on the command line
831  CommandLineArgs::doc_specified_flags();
832 
833  // Check that definition of Womersley number is consistent with those
834  // of the Reynolds and Strouhal numbers
837  {
838  std::ostringstream error_stream;
839  error_stream << "Definition of Global_Physical_Variables::ReSt is "
840  << "inconsistant with those\n"
841  << "of Global_Physical_Variables::Re and "
842  << "Global_Physical_Variables::St." << std::endl;
843  throw OomphLibError(error_stream.str(),
844  OOMPH_CURRENT_FUNCTION,OOMPH_EXCEPTION_LOCATION);
845  }
846 
847  /// Maximum time
848  double t_max = 0.6;
849 
850  /// Duration of timestep
851  const double dt = 0.0025;
852 
853  // If we are doing validation run, use smaller number of timesteps
854  if(CommandLineArgs::command_line_flag_has_been_set("--validation"))
855  {
856  t_max = 0.005;
857  }
858 
859  // Set direction of gravity (vertically downwards)
862 
863  // Set up the elastic test problem with QCrouzeixRaviartElements,
864  // using the BDF<2> timestepper
865  InterfaceProblem<RefineablePseudoSolidNodeUpdateElement<
866  RefineableQCrouzeixRaviartElement<2>,RefineableQPVDElement<2,3> >, BDF<2> >
867  problem;
868 
869  // Run the unsteady simulation
870  problem.unsteady_run(t_max,dt);
871 
872 } // End of main
double Lx
Width of domain.
void actions_before_adapt()
Strip off the interface elements before adapting the bulk mesh.
Vector< double > G(2)
Direction of gravity.
Mesh * Surface_mesh_pt
Pointer to the "surface" mesh.
ElasticRefineableTwoLayerMesh(const unsigned &nx, const unsigned &ny1, const unsigned &ny2, const double &lx, const double &h1, const double &h2, const bool &periodic_in_x, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction, number of elements in y-direction in bottom and ...
void delete_interface_elements()
Delete the 1d interface elements.
ofstream Trace_file
Trace file.
void actions_before_newton_solve()
No actions required before solve step.
double ReSt
Womersley number (Reynolds x Strouhal)
void actions_after_adapt()
Rebuild the mesh of interface elements after adapting the bulk mesh.
void fix_pressure(const unsigned &e, const unsigned &pdof, const double &pvalue)
Fix pressure in element e at pressure dof pdof and set to pvalue.
~InterfaceProblem()
Destructor (empty)
double Density_Ratio
Ratio of density in upper fluid to density in lower fluid. Reynolds number etc. is based on density i...
void actions_after_newton_solve()
No actions required after solve step.
void set_boundary_conditions()
Set boundary conditions.
void set_initial_condition()
Set initial conditions.
void unsteady_run(const double &t_max, const double &dt)
Do unsteady run up to maximum time t_max with given timestep dt.
double Viscosity_Ratio
Ratio of viscosity in upper fluid to viscosity in lower fluid. Reynolds number etc. is based on viscosity in lower fluid.
void create_interface_elements()
Create the 1d interface elements.
void deform_free_surface(const double &epsilon, const unsigned &n_periods)
Deform the mesh/free surface to a prescribed function.
double ReInvFr
Product of Reynolds number and inverse of Froude number.
ElasticRefineableTwoLayerMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the (specific) "bulk" mesh.
void actions_before_implicit_timestep()
Actions before the timestep: For maximum stability, reset the current nodal positions to be the "stre...
int main(int argc, char *argv[])
Driver code for two-dimensional two fluid interface problem.
void doc_solution(DocInfo &doc_info)
Document the solution.
ConstitutiveLaw * Constitutive_law_pt
double Nu
Pseudo-solid Poisson ratio.