inclined_plane.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 //A demo driver that solves the classic fluid flow problem of flow
31 //of a fluid film along an inclined plane. Stability analysis performed by
32 //Yih (1963), Benjamin (1957), ... and Blyth & Pozrikidis (2004).
33 
34 //This is an example of the subtleties involved in even a seemingly simple
35 //free surface problem.
36 
37 //Standard C++ library includes
38 #include <iostream>
39 #include <fstream>
40 #include <string>
41 #include <cmath>
42 
43 //Finite-Element library routines
44 #include "generic.h"
45 #include "navier_stokes.h"
46 #include "solid.h"
47 #include "fluid_interface.h"
48 #include "meshes/simple_rectangular_quadmesh.h"
49 
50 using namespace std;
51 
52 using namespace oomph;
53 
54 //The global physical variables
55 namespace Global_Physical_Variables
56 {
57  ///Reynolds number, based on the average velocity within the fluid film
58  double Re=0.0;
59 
60  ///The product of Reynolds number and inverse Froude number
61  ///is set to two in this problem, which gives the free surface velocity
62  ///to be sin(alpha). [Set to three in order to get the same scale as
63  ///used by Yih, Benjamin, etc]
64  double ReInvFr=2.0;
65 
66  ///Angle of incline of the slope (45 degrees)
67  double Alpha = 1.0*atan(1.0);
68 
69  ///\short The Vector direction of gravity, set in main()
70  Vector<double> G(2,0.0);
71 
72  ///The Capillary number
73  double Ca= 1.0;
74 
75  ///Set the wavenumber
76  double K = 0.1;
77 
78  ///Set the number of waves desired in the domain
79  double N_wave = 3;
80 
81  ///The length of the domain to fit the desired number of waves
82  double Length = 2*N_wave*4.0*atan(1.0)/K;
83 
84  ///Direction of the wall normal vector (at the inlet)
85  Vector<double> Wall_normal;
86 
87  /// \short Function that specifies the wall unit normal at the inlet
88  void wall_unit_normal_inlet_fct(const Vector<double> &x,
89  Vector<double> &normal)
90  {
91  normal=Wall_normal;
92  }
93 
94  /// \short Function that specified the wall unit normal at the outlet
95  void wall_unit_normal_outlet_fct(const Vector<double> &x,
96  Vector<double> &normal)
97  {
98  //Set the normal
99  normal = Wall_normal;
100  //and flip the sign
101  unsigned n_dim = normal.size();
102  for(unsigned i=0;i<n_dim;++i) {normal[i] *= -1.0;}
103  }
104 
105  ///The contact angle that is imposed at the inlet (pi)
106  double Inlet_Angle = 2.0*atan(1.0);
107 
108 
109  /// Function that prescribes the hydrostatic pressure field at the outlet
110  void hydrostatic_pressure_outlet(const double& time, const Vector<double> &x,
111  const Vector<double> &n,
112  Vector<double> &traction)
113  {
114  traction[0] = ReInvFr*G[1]*(1.0 - x[1]);
115  traction[1] = 0.0;
116  }
117 
118  /// Function that prescribes hydrostatic pressure field at the inlet
119  void hydrostatic_pressure_inlet(const double& time, const Vector<double> &x,
120  const Vector<double> &n,
121  Vector<double> &traction)
122  {
123  traction[0] = -ReInvFr*G[1]*(1.0 - x[1]);
124  traction[1] = 0.0;
125  }
126  //end of traction functions
127 
128  ///Constitutive law used to determine the mesh deformation
129  ConstitutiveLaw *Constitutive_law_pt;
130 
131  /// Pseudo-solid Poisson ratio
132  double Nu=0.1;
133 
134 }
135 
136 //=====================================================================
137 ///\short Generic problem class that will form the base class for both
138 ///spine and elastic mesh-updates of the problem.
139 ///Templated by the bulk element and interface element types
140 //====================================================================
141 template<class ELEMENT, class INTERFACE_ELEMENT>
142 class InclinedPlaneProblem : public Problem
143 {
144 
145 protected:
146 
147  ///Bulk fluid mesh
149 
150  ///Mesh for the traction elements that are added at inlet and outlet
152 
153  ///Mesh for the free surface elements
155 
156  ///Mesh for the point elements at each end of the free surface
158 
159  ///Prefix for output files
160  std::string Output_prefix;
161 
162 public:
163 
164  ///Generic Constructor (empty)
165  InclinedPlaneProblem(const unsigned &nx, const unsigned &ny,
166  const double &length) :
167  Output_prefix("Unset") { }
168 
169  ///Solve the steady problem
170  void solve_steady();
171 
172  ///Take n_tsteps timesteps of size dt
173  void timestep(const double &dt, const unsigned &n_tsteps);
174 
175  /// \short Actions before the timestep
176  /// (update the time-dependent boundary conditions)
178  {
179  //Read out the current time
180  double time = this->time_pt()->time();
181  //Now add a temporary sinusoidal suction and blowing to the base
182  //Amplitude of the perturbation
183  double epsilon = 0.01;
184  //Loop over the nodes on the base
185  unsigned n_node = this->Bulk_mesh_pt->nboundary_node(0);
186  for(unsigned n=0;n<n_node;n++)
187  {
188  Node* nod_pt = this->Bulk_mesh_pt->boundary_node_pt(0,n);
189  double arg = Global_Physical_Variables::K*nod_pt->x(0);
190  double value = sin(arg)*epsilon*time*exp(-time);
191  nod_pt->set_value(1,value);
192  }
193  } //end_of_actions_before_implicit_timestep
194 
195  ///\short Function to add the traction boundary elements to boundaries
196  /// 3(inlet) and 1(outlet) of the mesh
198  {
199  //Create a new (empty mesh)
200  Traction_mesh_pt = new Mesh;
201  //Inlet boundary conditions (boundary 3)
202  {
203  unsigned b = 3;
204  //Find the number of elements adjacent to mesh boundary
205  unsigned n_boundary_element = Bulk_mesh_pt->nboundary_element(b);
206  //Loop over these elements and create the traction elements
207  for(unsigned e=0;e<n_boundary_element;e++)
208  {
209  NavierStokesTractionElement<ELEMENT> *surface_element_pt =
210  new NavierStokesTractionElement<ELEMENT>
211  (Bulk_mesh_pt->boundary_element_pt(b,e),
212  Bulk_mesh_pt->face_index_at_boundary(b,e));
213  //Add the elements to the mesh
214  Traction_mesh_pt->add_element_pt(surface_element_pt);
215  //Set the traction function
216  surface_element_pt->traction_fct_pt() =
218  }
219  }
220 
221  //Outlet boundary conditions (boundary 1)
222  {
223  unsigned b=1;
224  //Find the number of elements adjacent to mesh boundary
225  unsigned n_boundary_element = Bulk_mesh_pt->nboundary_element(b);
226  //Loop over these elements and create the traction elements
227  for(unsigned e=0;e<n_boundary_element;e++)
228  {
229  NavierStokesTractionElement<ELEMENT> *surface_element_pt =
230  new NavierStokesTractionElement<ELEMENT>
231  (Bulk_mesh_pt->boundary_element_pt(b,e),
232  Bulk_mesh_pt->face_index_at_boundary(b,e));
233  //Add the elements to the mesh
234  Traction_mesh_pt->add_element_pt(surface_element_pt);
235  //Set the traction function
236  surface_element_pt->traction_fct_pt() =
238  }
239  }
240  } //end of make_traction_elements
241 
242  //Make the free surface elements on the top surface
244  {
245  //Create the (empty) meshes
246  Surface_mesh_pt = new Mesh;
247  Point_mesh_pt = new Mesh;
248 
249  //The free surface is on the boundary 2
250  unsigned b = 2;
251  unsigned n_boundary_element = Bulk_mesh_pt->nboundary_element(b);
252  //Loop over the elements and create the appropriate interface elements
253  for(unsigned e=0;e<n_boundary_element;e++)
254  {
255  INTERFACE_ELEMENT *surface_element_pt =
256  new INTERFACE_ELEMENT
257  (Bulk_mesh_pt->boundary_element_pt(b,e),
258  Bulk_mesh_pt->face_index_at_boundary(b,e));
259  //Add elements to the mesh
260  Surface_mesh_pt->add_element_pt(surface_element_pt);
261  //Assign the capillary number to the free surface
262  surface_element_pt->ca_pt() =
264 
265  //Make a point element from left-hand side of the
266  //first surface element (note that this relies on knowledge of
267  //the element order within the mesh)
268  if(e==0)
269  {
270  FluidInterfaceBoundingElement* point_element_pt =
271  surface_element_pt->make_bounding_element(-1);
272  //Add element to the point mesh
273  Point_mesh_pt->add_element_pt(point_element_pt);
274  //Set the capillary number
275  point_element_pt->ca_pt() = &Global_Physical_Variables::Ca;
276  //Set the wall normal
277  point_element_pt->wall_unit_normal_fct_pt() =
279  //Set the contact angle (using the strong version of the constraint)
280  point_element_pt->set_contact_angle(
282  }
283 
284  //Make another point element from the right-hand side of the
285  //last surface element (note that this relies on knowledge of
286  //the element order within the mesh)
287  if(e==n_boundary_element-1)
288  {
289  FluidInterfaceBoundingElement* point_element_pt =
290  surface_element_pt->make_bounding_element(1);
291  //Add element to the mesh
292  Point_mesh_pt->add_element_pt(point_element_pt);
293  //Set the capillary number
294  point_element_pt->ca_pt() = &Global_Physical_Variables::Ca;
295  // Set the function that specifies the wall normal
296  point_element_pt->wall_unit_normal_fct_pt() =
298  }
299  }
300  } //end of make_free_surface_elements
301 
302  ///\short Complete the build of the problem setting all standard
303  ///parameters and boundary conditions
305  {
306  using namespace Global_Physical_Variables;
307 
308  //Complete the build of the fluid elements by passing physical parameters
309  //Find the number of bulk elements
310  unsigned n_element = Bulk_mesh_pt->nelement();
311  //Loop over all the fluid elements
312  for(unsigned e=0;e<n_element;e++)
313  {
314  //Cast to a fluid element
315  ELEMENT *temp_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
316 
317  //Set the Reynolds number
318  temp_pt->re_pt() = &Re;
319  //The Strouhal number is 1, so ReSt = Re
320  temp_pt->re_st_pt() = &Re;
321  //Set the Reynolds number / Froude number
322  temp_pt->re_invfr_pt() = &ReInvFr;
323  //Set the direction of gravity
324  temp_pt->g_pt() = &G;
325  }
326 
327  //------------Set the boundary conditions for this problem----------
328 
329  {
330  //Determine whether we are solving an elastic problem or not
331  bool elastic = false;
332  if(dynamic_cast<SolidNode*>(Bulk_mesh_pt->node_pt(0))) {elastic=true;}
333 
334  //Loop over the bottom of the mesh (the wall of the channel)
335  unsigned n_node = Bulk_mesh_pt->nboundary_node(0);
336  for(unsigned j=0;j<n_node;j++)
337  {
338  //Pin the u- and v- velocities
339  Bulk_mesh_pt->boundary_node_pt(0,j)->pin(0);
340  Bulk_mesh_pt->boundary_node_pt(0,j)->pin(1);
341 
342  //If we are formulating the elastic problem pin both positions
343  //of nodes
344  if(elastic)
345  {
346  static_cast<SolidNode*>(Bulk_mesh_pt->boundary_node_pt(0,j))
347  ->pin_position(0);
348  static_cast<SolidNode*>(Bulk_mesh_pt->boundary_node_pt(0,j))
349  ->pin_position(1);
350  }
351  }
352 
353  //Loop over the inlet and set the Dirichlet condition
354  //of no vertical velocity
355  n_node = Bulk_mesh_pt->nboundary_node(3);
356  for(unsigned j=0;j<n_node;j++)
357  {
358  Bulk_mesh_pt->boundary_node_pt(3,j)->pin(1);
359 
360  //If elastic pin horizontal position of nodes
361  if(elastic)
362  {
363  static_cast<SolidNode*>(Bulk_mesh_pt->boundary_node_pt(3,j))
364  ->pin_position(0);
365  }
366  }
367 
368  //Loop over the outlet and set the Dirichlet condition
369  //of no vertical velocity
370  n_node = Bulk_mesh_pt->nboundary_node(1);
371  for(unsigned j=0;j<n_node;j++)
372  {
373  Bulk_mesh_pt->boundary_node_pt(1,j)->pin(1);
374 
375  //If elastic pin horizontal position
376  if(elastic)
377  {
378  static_cast<SolidNode*>(Bulk_mesh_pt->boundary_node_pt(1,j))
379  ->pin_position(0);
380  }
381  }
382  }
383 
384  //Attach the boundary conditions to the mesh
385  std::cout << assign_eqn_numbers() << " in the main problem" << std::endl;
386  } //end of complete_build
387 
388  ///Generic desructor to clean up the memory allocated in the problem
390  {
391  //Clear node storage before the mesh can be deleted,
392  //to avoid double deletion
393  this->Point_mesh_pt->flush_node_storage();
394  //Then delete the mesh
395  delete this->Point_mesh_pt;
396  //Clear node storage and then delete mesh
397  this->Surface_mesh_pt->flush_node_storage();
398  delete this->Surface_mesh_pt;
399  //Clear node storage and then delete mesh
400  this->Traction_mesh_pt->flush_node_storage();
401  delete this->Traction_mesh_pt;
402  //Delete the bulk mesh (no need to clear node storage)
403  delete this->Bulk_mesh_pt;
404  //Delete the time stepper
405  delete this->time_stepper_pt();
406  }
407 
408 };
409 
410 
411 //-------------------------------------------------------------------------
412 ///Solve the steady problem
413 //-------------------------------------------------------------------------
414 template<class ELEMENT,class INTERFACE_ELEMENT>
416 {
417  //Load the namespace
418  using namespace Global_Physical_Variables;
419 
420  //Initially set all nodes to the Nusselt flat-film solution
421  {
422  unsigned n_node = Bulk_mesh_pt->nnode();
423  for(unsigned n=0;n<n_node;n++)
424  {
425  double y = Bulk_mesh_pt->node_pt(n)->x(1);
426  //Top row
427  Bulk_mesh_pt->node_pt(n)->set_value(0,0.5*ReInvFr*sin(Alpha)*(2.0*y - y*y));
428  }
429  }
430 
431  //Do one steady solve
432  steady_newton_solve();
433 
434  //Output the full flow field
435  std::string filename = Output_prefix;;
436  filename.append("_output.dat");
437  ofstream file(filename.c_str());
438  Bulk_mesh_pt->output(file,5);
439  file.close();
440 } //end of solve_steady
441 
442 
443 //----------------------------------------------------------------------
444 /// Perform n_tsteps timesteps of size dt
445 //----------------------------------------------------------------------
446 template<class ELEMENT, class INTERFACE_ELEMENT>
448 timestep(const double &dt, const unsigned &n_tsteps)
449 {
450  //Need to use the Global variables here
451  using namespace Global_Physical_Variables;
452 
453  //Open an output file
454  std::string filename = Output_prefix;
455  filename.append("_time_trace.dat");
456  ofstream trace(filename.c_str());
457  //Counter that will be used to output the full flowfield
458  //at certain timesteps
459  int counter=0;
460 
461  //Initial output of the time and the value of the vertical position at the
462  //left and right-hand end of the free surface
463  trace << time_pt()->time() << " "
464  << Bulk_mesh_pt->boundary_node_pt(2,0)->value(1)
465  << " "
466  << Bulk_mesh_pt->
467  boundary_node_pt(2, Bulk_mesh_pt->nboundary_node(2)-1)->x(1)
468  << " "
469  << std::endl;
470 
471  //Loop over the desired number of timesteps
472  for(unsigned t=1;t<=n_tsteps;t++)
473  {
474  //Increase the counter
475  counter++;
476  cout << std::endl;
477  cout << "--------------TIMESTEP " << t<< " ------------------" << std::endl;
478 
479  //Take a timestep of size dt
480  unsteady_newton_solve(dt);
481 
482  //Uncomment to get full solution output
483  if(counter==2) //Change this number to get output every n steps
484  {
485  std::ofstream file;
486  std::ostringstream filename;
487  filename << Output_prefix << "_step" << Re << "_" << t << ".dat";
488  file.open(filename.str().c_str());
489  Bulk_mesh_pt->output(file,5);
490  file.close();
491 
492  counter=0;
493  }
494 
495  //Always output the interface
496  {
497  std::ofstream file;
498  std::ostringstream filename;
499  filename << Output_prefix << "_interface_" << Re << "_" << t << ".dat";
500  file.open(filename.str().c_str());
501  Surface_mesh_pt->output(file,5);
502  file.close();
503  }
504 
505  //Output the time and value of the vertical position of the free surface
506  //at the left- and right-hand ends
507  trace << time_pt()->time() << " "
508  << Bulk_mesh_pt->boundary_node_pt(2,0)->x(1) << " "
509  <<
510  Bulk_mesh_pt->
511  boundary_node_pt(2,Bulk_mesh_pt->nboundary_node(2)-1)->x(1) << " "
512  << std::endl;
513  }
514 } //end of timestep
515 
516 
517 //====================================================================
518 // Spine formulation of the problem
519 //===================================================================
520 
521 
522 //======================================================================
523 /// Create a spine mesh for the problem
524 //======================================================================
525 template <class ELEMENT>
527  public SimpleRectangularQuadMesh<ELEMENT>,
528  public SpineMesh
529 {
530 public:
531  SpineInclinedPlaneMesh(const unsigned &nx, const unsigned &ny,
532  const double &lx, const double &ly,
533  TimeStepper* time_stepper_pt) :
534  SimpleRectangularQuadMesh<ELEMENT>
535  (nx,ny,lx,ly,time_stepper_pt), SpineMesh()
536  {
537  //Find the number of linear points in the element
538  unsigned n_p = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
539  //Reserve storage for the number of spines
540  Spine_pt.reserve((n_p-1)*nx + 1);
541 
542  //Create single pointer to a spine
543  Spine* new_spine_pt=0;
544 
545  //Now loop over the elements horizontally
546  for(unsigned long j=0;j<nx;j++)
547  {
548  //In most elements, we don't assign a spine to the last column,
549  //beacuse that will be done by the next element
550  unsigned n_pmax = n_p-1;
551  //In the last element, however, we must assign the final spine
552  if(j==nx-1) {n_pmax = n_p;}
553 
554  //Loop over all nodes horizontally
555  for(unsigned l2=0;l2<n_pmax;l2++)
556  {
557  //Create a new spine with unit height and add to the mesh
558  new_spine_pt=new Spine(1.0);
559  Spine_pt.push_back(new_spine_pt);
560 
561  // Get the node
562  SpineNode* nod_pt=element_node_pt(j,l2);
563  //Set the pointer to spine
564  nod_pt->spine_pt() = new_spine_pt;
565  //Set the fraction
566  nod_pt->fraction() = 0.0;
567  // Pointer to the mesh that implements the update fct
568  nod_pt->spine_mesh_pt() = this;
569 
570  //Loop vertically along the spine
571  //Loop over the elements
572  for(unsigned long i=0;i<ny;i++)
573  {
574  //Loop over the vertical nodes, apart from the first
575  for(unsigned l1=1;l1<n_p;l1++)
576  {
577  // Get the node
578  SpineNode* nod_pt=element_node_pt(i*nx+j,l1*n_p+l2);
579  //Set the pointer to the spine
580  nod_pt->spine_pt() = new_spine_pt;
581  //Set the fraction
582  nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(ny);
583  // Pointer to the mesh that implements the update fct
584  nod_pt->spine_mesh_pt() = this;
585  }
586  }
587  }
588  } //End of horizontal loop over elements
589  } //end of constructor
590 
591  /// \short General node update function implements pure virtual function
592  /// defined in SpineMesh base class and performs specific node update
593  /// actions: along vertical spines
594  virtual void spine_node_update(SpineNode* spine_node_pt)
595  {
596  //Get fraction along the spine
597  double W = spine_node_pt->fraction();
598  //Get spine height
599  double H = spine_node_pt->h();
600  //Set the value of y
601  spine_node_pt->x(1) = W*H;
602  }
603 };
604 
605 
606 
607 //============================================================================
608 //Specific class for inclined plane problem using spines
609 //============================================================================
610 template<class ELEMENT, class TIMESTEPPER>
612  public InclinedPlaneProblem<ELEMENT,SpineLineFluidInterfaceElement<ELEMENT> >
613 {
614 public:
615 
616  //Constructor
617  SpineInclinedPlaneProblem(const unsigned &nx, const unsigned &ny,
618  const double &length):
619  InclinedPlaneProblem<ELEMENT,SpineLineFluidInterfaceElement<ELEMENT> >
620  (nx,ny,length)
621  {
622  //Set the name
623  this->Output_prefix = "spine";
624 
625  //Create our one and only timestepper, with adaptive timestepping
626  this->add_time_stepper_pt(new TIMESTEPPER);
627 
628  //Create the bulk mesh
629  this->Bulk_mesh_pt = new SpineInclinedPlaneMesh<ELEMENT>(
630  nx,ny,length,1.0,this->time_stepper_pt());
631 
632  //Create the traction elements
633  this->make_traction_elements();
634  //Create the free surface elements
635  this->make_free_surface_elements();
636 
637  //Add all sub meshes to the problem
638  this->add_sub_mesh(this->Bulk_mesh_pt);
639  this->add_sub_mesh(this->Traction_mesh_pt);
640  this->add_sub_mesh(this->Surface_mesh_pt);
641  this->add_sub_mesh(this->Point_mesh_pt);
642  //Create the global mesh
643  this->build_global_mesh();
644 
645  //Complete the build of the problem
646  this->complete_build();
647  }
648 
649  /// Spine heights/lengths are unknowns in the problem so their
650  /// values get corrected during each Newton step. However,
651  /// changing their value does not automatically change the
652  /// nodal positions, so we need to update all of them
654  {this->Bulk_mesh_pt->node_update();}
655 
656 };
657 
658 
659 //====================================================================
660 // Elastic formulation of the problem
661 //===================================================================
662 
663 
664 //======================================================================
665 /// Create an Elastic mesh for the problem
666 //======================================================================
667 template <class ELEMENT>
669  public SimpleRectangularQuadMesh<ELEMENT>,
670  public SolidMesh
671 {
672  //Public functions
673  public:
674  ElasticInclinedPlaneMesh(const unsigned &nx, const unsigned &ny,
675  const double &lx, const double &ly,
676  TimeStepper* time_stepper_pt) :
677  SimpleRectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,time_stepper_pt), SolidMesh()
678  {
679  //Make the current configuration the undeformed one
680  set_lagrangian_nodal_coordinates();
681  }
682 };
683 
684 
685 
686 
687 //============================================================================
688 //Specific class for inclined plane problem using pseudo-elastic formulation
689 //============================================================================
690 template<class ELEMENT, class TIMESTEPPER>
692  public InclinedPlaneProblem<ELEMENT,ElasticLineFluidInterfaceElement<ELEMENT> >
693 {
694 public:
695  //Constructor
696  ElasticInclinedPlaneProblem(const unsigned &nx, const unsigned &ny,
697  const double &length) :
698  InclinedPlaneProblem<ELEMENT,ElasticLineFluidInterfaceElement<ELEMENT> >
699  (nx,ny,length)
700  {
701  //Set the name
702  this->Output_prefix = "elastic";
703 
704  //Create our one and only timestepper, with adaptive timestepping
705  this->add_time_stepper_pt(new TIMESTEPPER);
706 
707  //Create the bulk mesh
708  this->Bulk_mesh_pt = new ElasticInclinedPlaneMesh<ELEMENT>(
709  nx,ny,length,1.0,this->time_stepper_pt());
710 
711  //Set the consititutive law for the elements
712  unsigned n_element = this->Bulk_mesh_pt->nelement();
713  //Loop over all the fluid elements
714  for(unsigned e=0;e<n_element;e++)
715  {
716  //Cast to a fluid element
717  ELEMENT *temp_pt = dynamic_cast<ELEMENT*>(
718  this->Bulk_mesh_pt->element_pt(e));
719  //Set the constitutive law
720  temp_pt->constitutive_law_pt() =
722  }
723 
724  //Create the traction elements
725  this->make_traction_elements();
726  //Create the free surface element
727  this->make_free_surface_elements();
728 
729  //Add all sub meshes to the problem
730  this->add_sub_mesh(this->Bulk_mesh_pt);
731  this->add_sub_mesh(this->Traction_mesh_pt);
732  this->add_sub_mesh(this->Surface_mesh_pt);
733  this->add_sub_mesh(this->Point_mesh_pt);
734  //Create the global mesh
735  this->build_global_mesh();
736 
737  //Complete the rest of the build
738  this->complete_build();
739  } //end of constructor
740 
741  ///\short Update Lagrangian positions after each timestep
742  ///(updated-lagrangian approach)
744  {
745  //Now loop over all the nodes and reset their Lagrangian coordinates
746  unsigned n_node = this->Bulk_mesh_pt->nnode();
747  for(unsigned n=0;n<n_node;n++)
748  {
749  //Cast node to an elastic node
750  SolidNode* temp_pt =
751  static_cast<SolidNode*>(this->Bulk_mesh_pt->node_pt(n));
752  for(unsigned j=0;j<2;j++) {temp_pt->xi(j) = temp_pt->x(j);}
753  }
754  } //end of actions_after_implicit_timestep
755 
756 };
757 
758 
759 //start of main
760 int main(int argc, char **argv)
761 {
762  using namespace Global_Physical_Variables;
763 
764  //Set the constitutive law for the mesh deformation
765  Constitutive_law_pt = new GeneralisedHookean(&Global_Physical_Variables::Nu);
766 
767 #ifdef CR_ELEMENT
768 #define FLUID_ELEMENT QCrouzeixRaviartElement<2>
769 #else
770 #define FLUID_ELEMENT QTaylorHoodElement<2>
771 #endif
772 
773  //Initialise physical parameters
774  //Scale Reynolds number to be independent of alpha.
775  Re = 4.0/sin(Alpha);
776 
777  //Set the direction of gravity
778  G[0] = sin(Alpha);
779  G[1] = -cos(Alpha);
780 
781  //The wall normal to the inlet is in the negative x direction
782  Wall_normal.resize(2);
783  Wall_normal[0] = -1.0;
784  Wall_normal[1] = 0.0;
785 
786  //Spine problem
787  {
788  //Create the problem
790  problem(30,4,Global_Physical_Variables::Length);
791 
792  //Solve the steady problem
793  problem.solve_steady();
794 
795  //Prepare the problem for timestepping
796  //(assume that it's been at the flat-film solution for all previous time)
797  double dt = 0.1;
798  problem.assign_initial_values_impulsive(dt);
799 
800  //Timestep it
801  problem.timestep(dt,2);
802  } //End of spine problem
803 
804 
805  //Elastic problem
806  {
807  //Create the problem
809  PseudoSolidNodeUpdateElement<FLUID_ELEMENT,QPVDElement<2,3> >, BDF<2> >
810  problem(30,4,Global_Physical_Variables::Length);
811 
812  //Solve the steady problem
813  problem.solve_steady();
814 
815  //Prepare the problem for timestepping
816  //(assume that it's been at the flat-film solution for all previous time)
817  double dt = 0.1;
818  problem.assign_initial_values_impulsive(dt);
819 
820  //Timestep it
821  problem.timestep(dt,2);
822  } //End of elastic problem
823 }
void hydrostatic_pressure_inlet(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Function that prescribes hydrostatic pressure field at the inlet.
double K
Set the wavenumber.
void timestep(const double &dt, const unsigned &n_tsteps)
Take n_tsteps timesteps of size dt.
double N_wave
Set the number of waves desired in the domain.
Mesh * Point_mesh_pt
Mesh for the point elements at each end of the free surface.
void actions_before_implicit_timestep()
Actions before the timestep (update the time-dependent boundary conditions)
Mesh * Bulk_mesh_pt
Bulk fluid mesh.
Mesh * Surface_mesh_pt
Mesh for the free surface elements.
Vector< double > Wall_normal
Direction of the wall normal vector (at the inlet)
Generic problem class that will form the base class for both spine and elastic mesh-updates of the pr...
std::string Output_prefix
Prefix for output files.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
virtual void spine_node_update(SpineNode *spine_node_pt)
General node update function implements pure virtual function defined in SpineMesh base class and per...
void actions_before_newton_convergence_check()
SpineInclinedPlaneProblem(const unsigned &nx, const unsigned &ny, const double &length)
void actions_after_implicit_timestep()
Update Lagrangian positions after each timestep (updated-lagrangian approach)
int main(int argc, char **argv)
ElasticInclinedPlaneProblem(const unsigned &nx, const unsigned &ny, const double &length)
void wall_unit_normal_outlet_fct(const Vector< double > &x, Vector< double > &normal)
Function that specified the wall unit normal at the outlet.
double Re
Reynolds number, based on the average velocity within the fluid film.
Vector< double > G(2, 0.0)
The Vector direction of gravity, set in main()
void wall_unit_normal_inlet_fct(const Vector< double > &x, Vector< double > &normal)
Function that specifies the wall unit normal at the inlet.
~InclinedPlaneProblem()
Generic desructor to clean up the memory allocated in the problem.
void complete_build()
Complete the build of the problem setting all standard parameters and boundary conditions.
Mesh * Traction_mesh_pt
Mesh for the traction elements that are added at inlet and outlet.
double Length
The length of the domain to fit the desired number of waves.
void make_traction_elements()
Function to add the traction boundary elements to boundaries 3(inlet) and 1(outlet) of the mesh...
ElasticInclinedPlaneMesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly, TimeStepper *time_stepper_pt)
double Inlet_Angle
The contact angle that is imposed at the inlet (pi)
Create a spine mesh for the problem.
void hydrostatic_pressure_outlet(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Function that prescribes the hydrostatic pressure field at the outlet.
SpineInclinedPlaneMesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly, TimeStepper *time_stepper_pt)
InclinedPlaneProblem(const unsigned &nx, const unsigned &ny, const double &length)
Generic Constructor (empty)
double Nu
Pseudo-solid Poisson ratio.
void solve_steady()
Solve the steady problem.
Create an Elastic mesh for the problem.
double Ca
The Capillary number.
double Alpha
Angle of incline of the slope (45 degrees)