fsi_collapsible_channel.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 #include <iostream>
31 
32 // Generic oomph-lib includes
33 #include "generic.h"
34 #include "navier_stokes.h"
35 #include "beam.h"
36 
37 // The wall mesh
39 
40 //Include the fluid mesh
42 
43 using namespace std;
44 
45 using namespace oomph;
46 
47 
48 //==========start_of_BL_Squash =========================================
49 /// Namespace to define the mapping [0,1] -> [0,1] that re-distributes
50 /// nodal points across the channel width.
51 //======================================================================
52 namespace BL_Squash
53 {
54 
55  /// Boundary layer width
56  double Delta=0.1;
57 
58  /// Fraction of points in boundary layer
59  double Fract_in_BL=0.5;
60 
61  /// \short Mapping [0,1] -> [0,1] that re-distributes
62  /// nodal points across the channel width
63  double squash_fct(const double& s)
64  {
65  // Default return
66  double y=s;
67  if (s<0.5*Fract_in_BL)
68  {
69  y=Delta*2.0*s/Fract_in_BL;
70  }
71  else if (s>1.0-0.5*Fract_in_BL)
72  {
73  y=2.0*Delta/Fract_in_BL*s+1.0-2.0*Delta/Fract_in_BL;
74  }
75  else
76  {
77  y=(1.0-2.0*Delta)/(1.0-Fract_in_BL)*s+
78  (Delta-0.5*Fract_in_BL)/(1.0-Fract_in_BL);
79  }
80 
81  return y;
82  }
83 }// end of BL_Squash
84 
85 
86 
87 
88 
89 
90 //////////////////////////////////////////////////////////////////////////
91 //////////////////////////////////////////////////////////////////////////
92 //////////////////////////////////////////////////////////////////////////
93 
94 
95 //====start_of_underformed_wall============================================
96 /// Undeformed wall is a steady, straight 1D line in 2D space
97 /// \f[ x = X_0 + \zeta \f]
98 /// \f[ y = H \f]
99 //=========================================================================
100 class UndeformedWall : public GeomObject
101 {
102 
103 public:
104 
105  /// \short Constructor: arguments are the starting point and the height
106  /// above y=0.
107  UndeformedWall(const double& x0, const double& h): GeomObject(1,2)
108  {
109  X0=x0;
110  H=h;
111  }
112 
113 
114  /// \short Position vector at Lagrangian coordinate zeta
115  void position(const Vector<double>& zeta, Vector<double>& r) const
116  {
117  // Position Vector
118  r[0] = zeta[0]+X0;
119  r[1] = H;
120  }
121 
122 
123  /// \short Parametrised position on object: r(zeta). Evaluated at
124  /// previous timestep. t=0: current time; t>0: previous
125  /// timestep. Calls steady version.
126  void position(const unsigned& t, const Vector<double>& zeta,
127  Vector<double>& r) const
128  {
129  // Use the steady version
130  position(zeta,r);
131 
132  } // end of position
133 
134 
135  /// \short Posn vector and its 1st & 2nd derivatives
136  /// w.r.t. to coordinates:
137  /// \f$ \frac{dR_i}{d \zeta_\alpha}\f$ = drdzeta(alpha,i).
138  /// \f$ \frac{d^2R_i}{d \zeta_\alpha d \zeta_\beta}\f$ =
139  /// ddrdzeta(alpha,beta,i). Evaluated at current time.
140  void d2position(const Vector<double>& zeta,
141  Vector<double>& r,
142  DenseMatrix<double> &drdzeta,
143  RankThreeTensor<double> &ddrdzeta) const
144  {
145  // Position vector
146  r[0] = zeta[0]+X0;
147  r[1] = H;
148 
149  // Tangent vector
150  drdzeta(0,0)=1.0;
151  drdzeta(0,1)=0.0;
152 
153  // Derivative of tangent vector
154  ddrdzeta(0,0,0)=0.0;
155  ddrdzeta(0,0,1)=0.0;
156 
157  } // end of d2position
158 
159  private :
160 
161  /// x position of the undeformed beam's left end.
162  double X0;
163 
164  /// Height of the undeformed wall above y=0.
165  double H;
166 
167 }; //end_of_undeformed_wall
168 
169 
170 
171 
172 
173 
174 
175 //====start_of_physical_parameters=====================
176 /// Namespace for phyical parameters
177 //======================================================
178 namespace Global_Physical_Variables
179 {
180  /// Reynolds number
181  double Re=50.0;
182 
183  /// Womersley = Reynolds times Strouhal
184  double ReSt=50.0;
185 
186  /// Default pressure on the left boundary
187  double P_up=0.0;
188 
189  /// Traction applied on the fluid at the left (inflow) boundary
190  void prescribed_traction(const double& t,
191  const Vector<double>& x,
192  const Vector<double>& n,
193  Vector<double>& traction)
194  {
195  traction.resize(2);
196  traction[0]=P_up;
197  traction[1]=0.0;
198 
199  } //end traction
200 
201  /// Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.
202  double H=1.0e-2;
203 
204  /// 2nd Piola Kirchhoff pre-stress. As in Jensen & Heil (2003) paper.
205  double Sigma0=1.0e3;
206 
207  /// External pressure
208  double P_ext=0.0;
209 
210  /// \short Load function: Apply a constant external pressure to the wall.
211  /// Note: This is the load without the fluid contribution!
212  /// Fluid load gets added on by FSIWallElement.
213  void load(const Vector<double>& xi, const Vector<double>& x,
214  const Vector<double>& N, Vector<double>& load)
215  {
216  for(unsigned i=0;i<2;i++)
217  {
218  load[i] = -P_ext*N[i];
219  }
220  } //end of load
221 
222 
223  /// \short Fluid structure interaction parameter: Ratio of stresses used for
224  /// non-dimensionalisation of fluid to solid stresses.
225  double Q=1.0e-5;
226 
227 
228 } // end of namespace
229 
230 
231 
232 
233 //====start_of_problem_class==========================================
234 ///Problem class
235 //====================================================================
236 template <class ELEMENT>
237 class FSICollapsibleChannelProblem : public Problem
238 {
239 
240  public :
241 
242 /// \short Constructor: The arguments are the number of elements and
243 /// the lengths of the domain.
244  FSICollapsibleChannelProblem(const unsigned& nup,
245  const unsigned& ncollapsible,
246  const unsigned& ndown,
247  const unsigned& ny,
248  const double& lup,
249  const double& lcollapsible,
250  const double& ldown,
251  const double& ly);
252 
253  /// Destructor (empty)
255 
256 
257 #ifdef MACRO_ELEMENT_NODE_UPDATE
258 
259  /// Access function for the specific bulk (fluid) mesh
261  {
262  // Upcast from pointer to the Mesh base class to the specific
263  // element type that we're using here.
264  return dynamic_cast<
266  (Bulk_mesh_pt);
267  }
268 
269 #else
270 
271  /// Access function for the specific bulk (fluid) mesh
273  {
274  // Upcast from pointer to the Mesh base class to the specific
275  // element type that we're using here.
276  return dynamic_cast<
278  (Bulk_mesh_pt);
279  }
280 
281 #endif
282 
283 
284  /// Access function for the wall mesh
286  {
287  return Wall_mesh_pt;
288 
289  } // end of access to wall mesh
290 
291 
292  /// \short Update the problem specs before solve (empty)
294 
295  /// Update the problem after solve (empty)
297 
298  /// \short Update before checking Newton convergence: Update the
299  /// nodal positions in the fluid mesh in response to possible
300  /// changes in the wall shape
302  {
303  Bulk_mesh_pt->node_update();
304  }
305 
306  /// Doc the solution
307  void doc_solution(DocInfo& doc_info,ofstream& trace_file);
308 
309  /// Apply initial conditions
310  void set_initial_condition();
311 
312 private :
313 
314  /// Create the prescribed traction elements on boundary b
315  void create_traction_elements(const unsigned &b,
316  Mesh* const &bulk_mesh_pt,
317  Mesh* const &traction_mesh_pt);
318 
319  ///Number of elements in the x direction in the upstream part of the channel
320  unsigned Nup;
321 
322  /// \short Number of elements in the x direction in the collapsible part of
323  /// the channel
324  unsigned Ncollapsible;
325 
326  ///Number of elements in the x direction in the downstream part of the channel
327  unsigned Ndown;
328 
329  ///Number of elements across the channel
330  unsigned Ny;
331 
332  ///x-length in the upstream part of the channel
333  double Lup;
334 
335  ///x-length in the collapsible part of the channel
336  double Lcollapsible;
337 
338  ///x-length in the downstream part of the channel
339  double Ldown;
340 
341  ///Transverse length
342  double Ly;
343 
344 #ifdef MACRO_ELEMENT_NODE_UPDATE
345 
346  /// Pointer to the "bulk" mesh
348 
349 #else
350 
351  /// Pointer to the "bulk" mesh
353 
354 #endif
355 
356 
357  /// \short Pointer to the "surface" mesh that applies the traction at the
358  /// inflow
360 
361  /// Pointer to the "wall" mesh
363 
364  ///Pointer to the left control node
366 
367  ///Pointer to right control node
369 
370  /// Pointer to control node on the wall
372 
373 };//end of problem class
374 
375 
376 
377 
378 //=====start_of_constructor======================================
379 /// Constructor for the collapsible channel problem
380 //===============================================================
381 template <class ELEMENT>
383  const unsigned& nup,
384  const unsigned& ncollapsible,
385  const unsigned& ndown,
386  const unsigned& ny,
387  const double& lup,
388  const double& lcollapsible,
389  const double& ldown,
390  const double& ly)
391 {
392  // Store problem parameters
393  Nup=nup;
394  Ncollapsible=ncollapsible;
395  Ndown=ndown;
396  Ny=ny;
397  Lup=lup;
398  Lcollapsible=lcollapsible;
399  Ldown=ldown;
400  Ly=ly;
401 
402 
403  // Overwrite maximum allowed residual to accomodate bad initial guesses
404  Problem::Max_residuals=1000.0;
405 
406  // Allocate the timestepper -- this constructs the Problem's
407  // time object with a sufficient amount of storage to store the
408  // previous timsteps.
409  add_time_stepper_pt(new BDF<2>);
410 
411  // Geometric object that represents the undeformed wall:
412  // A straight line at height y=ly; starting at x=lup.
413  UndeformedWall* undeformed_wall_pt=new UndeformedWall(lup,ly);
414 
415  //Create the "wall" mesh with FSI Hermite elements
417  //(2*Ncollapsible+5,Lcollapsible,undeformed_wall_pt);
418  (Ncollapsible,Lcollapsible,undeformed_wall_pt);
419 
420 
421 
422  // Build a geometric object (one Lagrangian, two Eulerian coordinates)
423  // from the wall mesh
424  MeshAsGeomObject* wall_geom_object_pt=
425  new MeshAsGeomObject(Wall_mesh_pt);
426 
427 
428 #ifdef MACRO_ELEMENT_NODE_UPDATE
429 
430  //Build bulk (fluid) mesh
431  Bulk_mesh_pt =
433  (nup, ncollapsible, ndown, ny,
434  lup, lcollapsible, ldown, ly,
435  wall_geom_object_pt,
436  time_stepper_pt());
437 
438  // Set a non-trivial boundary-layer-squash function...
439  Bulk_mesh_pt->bl_squash_fct_pt() = &BL_Squash::squash_fct;
440 
441  // ... and update the nodal positions accordingly
442  Bulk_mesh_pt->node_update();
443 
444 #else
445 
446  //Build bulk (fluid) mesh
447  Bulk_mesh_pt =
449  (nup, ncollapsible, ndown, ny,
450  lup, lcollapsible, ldown, ly,
451  wall_geom_object_pt,
453  time_stepper_pt());
454 
455 #endif
456 
457 
458  // Create "surface mesh" that will contain only the prescribed-traction
459  // elements. The constructor just creates the mesh without
460  // giving it any elements, nodes, etc.
461  Applied_fluid_traction_mesh_pt = new Mesh;
462 
463  // Create prescribed-traction elements from all elements that are
464  // adjacent to boundary 5 (left boundary), but add them to a separate mesh.
465  create_traction_elements(5,Bulk_mesh_pt,Applied_fluid_traction_mesh_pt);
466 
467  // Add the sub meshes to the problem
468  add_sub_mesh(Bulk_mesh_pt);
469  add_sub_mesh(Applied_fluid_traction_mesh_pt);
470  add_sub_mesh(Wall_mesh_pt);
471 
472  // Combine all submeshes into a single Mesh
473  build_global_mesh();
474 
475 
476  // Complete build of fluid mesh
477  //-----------------------------
478 
479  // Loop over the elements to set up element-specific
480  // things that cannot be handled by constructor
481  unsigned n_element=Bulk_mesh_pt->nelement();
482  for(unsigned e=0;e<n_element;e++)
483  {
484  // Upcast from GeneralisedElement to the present element
485  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
486 
487  //Set the Reynolds number
488  el_pt->re_pt() = &Global_Physical_Variables::Re;
489 
490  // Set the Womersley number
491  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
492 
493  } // end loop over elements
494 
495 
496 
497  // Apply boundary conditions for fluid
498  //------------------------------------
499 
500  //Pin the velocity on the boundaries
501  //x and y-velocities pinned along boundary 0 (bottom boundary) :
502  unsigned ibound=0;
503  unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
504  for (unsigned inod=0;inod<num_nod;inod++)
505  {
506  for(unsigned i=0;i<2;i++)
507  {
508  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
509  }
510  }
511 
512  //x and y-velocities pinned along boundaries 2, 3, 4 (top boundaries) :
513  for(ibound=2;ibound<5;ibound++)
514  {
515  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
516  for (unsigned inod=0;inod<num_nod;inod++)
517  {
518  for(unsigned i=0;i<2;i++)
519  {
520  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
521  }
522  }
523  }
524 
525  //y-velocity pinned along boundary 1 (right boundary):
526  ibound=1;
527  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
528  for (unsigned inod=0;inod<num_nod;inod++)
529  {
530  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
531  }
532 
533 
534  //y-velocity pinned along boundary 5 (left boundary):
535  ibound=5;
536  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
537  for (unsigned inod=0;inod<num_nod;inod++)
538  {
539  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
540  }
541 //end of pin_velocity
542 
543  // Complete build of applied traction elements
544  //--------------------------------------------
545 
546  // Loop over the traction elements to pass pointer to prescribed
547  // traction function
548  unsigned n_el=Applied_fluid_traction_mesh_pt->nelement();
549  for(unsigned e=0;e<n_el;e++)
550  {
551  // Upcast from GeneralisedElement to NavierStokes traction element
552  NavierStokesTractionElement<ELEMENT> *el_pt =
553  dynamic_cast< NavierStokesTractionElement<ELEMENT>*>(
554  Applied_fluid_traction_mesh_pt->element_pt(e));
555 
556  // Set the pointer to the prescribed traction function
557  el_pt->traction_fct_pt() = &Global_Physical_Variables::prescribed_traction;
558  }
559 
560 
561 
562 
563  // Complete build of wall elements
564  //--------------------------------
565 
566  //Loop over the elements to set physical parameters etc.
567  n_element = wall_mesh_pt()->nelement();
568  for(unsigned e=0;e<n_element;e++)
569  {
570  // Upcast to the specific element type
571  FSIHermiteBeamElement *elem_pt =
572  dynamic_cast<FSIHermiteBeamElement*>(wall_mesh_pt()->element_pt(e));
573 
574  // Set physical parameters for each element:
575  elem_pt->sigma0_pt() = &Global_Physical_Variables::Sigma0;
576  elem_pt->h_pt() = &Global_Physical_Variables::H;
577 
578  // Set the load vector for each element
579  elem_pt->load_vector_fct_pt() = &Global_Physical_Variables::load;
580 
581  // Function that specifies the load ratios
582  elem_pt->q_pt() = &Global_Physical_Variables::Q;
583 
584  // Set the undeformed shape for each element
585  elem_pt->undeformed_beam_pt() = undeformed_wall_pt;
586 
587 
588  // The normal on the wall elements as computed by the FSIHermiteElements
589  // points away from the fluid rather than into the fluid (as assumed
590  // by default)
591  elem_pt->set_normal_pointing_out_of_fluid();
592 
593  } // end of loop over elements
594 
595 
596 
597  // Boundary conditions for wall mesh
598  //----------------------------------
599 
600  // Set the boundary conditions: Each end of the beam is fixed in space
601  // Loop over the boundaries (ends of the beam)
602  for(unsigned b=0;b<2;b++)
603  {
604  // Pin displacements in both x and y directions
605  wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(0);
606  wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(1);
607  }
608 
609 
610 
611 
612  //Choose control nodes
613  //---------------------
614 
615  // Left boundary
616  ibound=5;
617  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
618  unsigned control_nod=num_nod/2;
619  Left_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
620 
621  // Right boundary
622  ibound=1;
623  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
624  control_nod=num_nod/2;
625  Right_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
626 
627 
628  // Set the pointer to the control node on the wall
629  Wall_node_pt=wall_mesh_pt()->node_pt(Ncollapsible/2);
630 
631 
632 
633 
634  // Setup FSI
635  //----------
636 
637  // The velocity of the fluid nodes on the wall (fluid mesh boundary 3)
638  // is set by the wall motion -- hence the no-slip condition must be
639  // re-applied whenever a node update is performed for these nodes.
640  // Such tasks may be performed automatically by the auxiliary node update
641  // function specified by a function pointer:
642  ibound=3;
643  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
644  for (unsigned inod=0;inod<num_nod;inod++)
645  {
646  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->
647  set_auxiliary_node_update_fct_pt(
648  FSI_functions::apply_no_slip_on_moving_wall);
649  }
650 
651 
652  // Work out which fluid dofs affect the residuals of the wall elements:
653  // We pass the boundary between the fluid and solid meshes and
654  // pointers to the meshes. The interaction boundary is boundary 3 of the
655  // 2D fluid mesh.
656  FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
657  (this,3,Bulk_mesh_pt,Wall_mesh_pt);
658 
659  // Setup equation numbering scheme
660  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
661 
662 
663 }//end of constructor
664 
665 
666 
667 
668 //====start_of_doc_solution===================================================
669 /// Doc the solution
670 //============================================================================
671 template <class ELEMENT>
673  ofstream& trace_file)
674 {
675 
676 
677 #ifdef MACRO_ELEMENT_NODE_UPDATE
678 
679  // Doc fsi
680  if (CommandLineArgs::Argc>1)
681  {
682  FSI_functions::doc_fsi<MacroElementNodeUpdateNode>(Bulk_mesh_pt,Wall_mesh_pt,doc_info);
683  }
684 
685 #else
686 
687  // Doc fsi
688  if (CommandLineArgs::Argc>1)
689  {
690  FSI_functions::doc_fsi<AlgebraicNode>(Bulk_mesh_pt,Wall_mesh_pt,doc_info);
691  }
692 
693 #endif
694 
695  ofstream some_file;
696  char filename[100];
697 
698  // Number of plot points
699  unsigned npts;
700  npts=5;
701 
702  // Output fluid solution
703  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
704  doc_info.number());
705  some_file.open(filename);
706  bulk_mesh_pt()->output(some_file,npts);
707  some_file.close();
708 
709  // Document the wall shape
710  sprintf(filename,"%s/beam%i.dat",doc_info.directory().c_str(),
711  doc_info.number());
712  some_file.open(filename);
713  wall_mesh_pt()->output(some_file,npts);
714  some_file.close();
715 
716 
717  // Write trace file
718  trace_file << time_pt()->time() << " "
719  << Wall_node_pt->x(1) << " "
720  << Left_node_pt->value(0) << " "
721  << Right_node_pt->value(0) << " "
723  << std::endl;
724 
725 } // end_of_doc_solution
726 
727 
728 
729 
730 //=====start_of_create_traction_elements======================================
731 /// Create the traction elements
732 //============================================================================
733 template <class ELEMENT>
735  const unsigned &b, Mesh* const &bulk_mesh_pt, Mesh* const &traction_mesh_pt)
736 {
737 
738  // How many bulk elements are adjacent to boundary b?
739  unsigned n_element = bulk_mesh_pt->nboundary_element(b);
740 
741  // Loop over the bulk elements adjacent to boundary b?
742  for(unsigned e=0;e<n_element;e++)
743  {
744  // Get pointer to the bulk element that is adjacent to boundary b
745  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>
746  (bulk_mesh_pt->boundary_element_pt(b,e));
747 
748  //What is the index of the face of element e along boundary b
749  int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
750 
751  // Build the corresponding prescribed-traction element
752  NavierStokesTractionElement<ELEMENT>* flux_element_pt =
753  new NavierStokesTractionElement<ELEMENT>(bulk_elem_pt,face_index);
754 
755  //Add the prescribed-traction element to the surface mesh
756  traction_mesh_pt->add_element_pt(flux_element_pt);
757 
758  } //end of loop over bulk elements adjacent to boundary b
759 
760 } // end of create_traction_elements
761 
762 
763 
764 
765 //====start_of_apply_initial_condition========================================
766 /// Apply initial conditions
767 //============================================================================
768 template <class ELEMENT>
770 {
771  // Check that timestepper is from the BDF family
772  if (time_stepper_pt()->type()!="BDF")
773  {
774  std::ostringstream error_stream;
775  error_stream << "Timestepper has to be from the BDF family!\n"
776  << "You have specified a timestepper from the "
777  << time_stepper_pt()->type() << " family" << std::endl;
778 
779  throw OomphLibError(error_stream.str(),
780  OOMPH_CURRENT_FUNCTION,
781  OOMPH_EXCEPTION_LOCATION);
782  }
783 
784  // Update the mesh
785  bulk_mesh_pt()->node_update();
786 
787  // Loop over the nodes to set initial guess everywhere
788  unsigned num_nod = bulk_mesh_pt()->nnode();
789  for (unsigned n=0;n<num_nod;n++)
790  {
791  // Get nodal coordinates
792  Vector<double> x(2);
793  x[0]=bulk_mesh_pt()->node_pt(n)->x(0);
794  x[1]=bulk_mesh_pt()->node_pt(n)->x(1);
795 
796  // Assign initial condition: Steady Poiseuille flow
797  bulk_mesh_pt()->node_pt(n)->set_value(0,6.0*(x[1]/Ly)*(1.0-(x[1]/Ly)));
798  bulk_mesh_pt()->node_pt(n)->set_value(1,0.0);
799  }
800 
801  // Assign initial values for an impulsive start
802  bulk_mesh_pt()->assign_initial_values_impulsive();
803 
804 } // end of set_initial_condition
805 
806 
807 
808 
809 
810 //============start_of_main====================================================
811 /// Driver code for a collapsible channel problem with FSI.
812 /// Presence of command line arguments indicates validation run with
813 /// coarse resolution and small number of timesteps.
814 //=============================================================================
815 int main(int argc, char* argv[])
816 {
817 
818  // Store command line arguments
819  CommandLineArgs::setup(argc,argv);
820 
821  // Reduction in resolution for validation run?
822  unsigned coarsening_factor=1;
823  if (CommandLineArgs::Argc>1)
824  {
825  coarsening_factor=4;
826  }
827 
828  // Number of elements in the domain
829  unsigned nup=20/coarsening_factor;
830  unsigned ncollapsible=40/coarsening_factor;
831  unsigned ndown=40/coarsening_factor;
832  unsigned ny=16/coarsening_factor;
833 
834  // Length of the domain
835  double lup=5.0;
836  double lcollapsible=10.0;
837  double ldown=10.0;
838  double ly=1.0;
839 
840  // Set external pressure (on the wall stiffness scale).
842 
843  // Pressure on the left boundary: This is consistent with steady
844  // Poiseuille flow
845  Global_Physical_Variables::P_up=12.0*(lup+lcollapsible+ldown);
846 
847 #ifdef MACRO_ELEMENT_NODE_UPDATE
848 
849 #ifdef TAYLOR_HOOD
850 
851  // Build the problem with QTaylorHoodElements
853  <MacroElementNodeUpdateElement<QTaylorHoodElement<2> > >
854  problem(nup, ncollapsible, ndown, ny,
855  lup, lcollapsible, ldown, ly);
856 
857 #else
858 
859  // Build the problem with QCrouzeixRaviartElements
861  <MacroElementNodeUpdateElement<QCrouzeixRaviartElement<2> > >
862  problem(nup, ncollapsible, ndown, ny,
863  lup, lcollapsible, ldown, ly);
864 
865 #endif
866 
867 #else
868 
869 #ifdef TAYLOR_HOOD
870 
871  // Build the problem with QCrouzeixRaviartElements
873  <AlgebraicElement<QTaylorHoodElement<2> > >
874  problem(nup, ncollapsible, ndown, ny,
875  lup, lcollapsible, ldown, ly);
876 
877 #else
878 
879  // Build the problem with QCrouzeixRaviartElements
881  <AlgebraicElement<QCrouzeixRaviartElement<2> > >
882  problem(nup, ncollapsible, ndown, ny,
883  lup, lcollapsible, ldown, ly);
884 
885 #endif
886 
887 #endif
888  // Timestep. Note: Preliminary runs indicate that the period of
889  // the oscillation is about 1 so this gives us 40 steps per period.
890  double dt=1.0/40.0;
891 
892  // Initial time for the simulation
893  double t_min=0.0;
894 
895  // Maximum time for simulation
896  double t_max=3.5;
897 
898  // Initialise timestep
899  problem.time_pt()->time()=t_min;
900  problem.initialise_dt(dt);
901 
902  // Apply initial condition
903  problem.set_initial_condition();
904 
905  //Set output directory
906  DocInfo doc_info;
907  doc_info.set_directory("RESLT");
908 
909  // Open a trace file
910  ofstream trace_file;
911  char filename[100];
912  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
913  trace_file.open(filename);
914 
915  // Output the initial solution
916  problem.doc_solution(doc_info, trace_file);
917 
918  // Increment step number
919  doc_info.number()++;
920 
921  // Find number of timesteps (reduced for validation)
922  unsigned nstep = unsigned((t_max-t_min)/dt);
923  if (CommandLineArgs::Argc>1)
924  {
925  nstep=3;
926  }
927 
928  // Timestepping loop
929  for (unsigned istep=0;istep<nstep;istep++)
930  {
931  // Solve the problem
932  problem.unsteady_newton_solve(dt);
933 
934  // Outpt the solution
935  problem.doc_solution(doc_info, trace_file);
936 
937  // Step number
938  doc_info.number()++;
939  }
940 
941 
942  // Close trace file.
943  trace_file.close();
944 
945 }//end of main
946 
void actions_before_newton_convergence_check()
Update before checking Newton convergence: Update the nodal positions in the fluid mesh in response t...
OneDLagrangianMesh< FSIHermiteBeamElement > * wall_mesh_pt()
Access function for the wall mesh.
void actions_before_newton_solve()
Update the problem specs before solve (empty)
~FSICollapsibleChannelProblem()
Destructor (empty)
MacroElementNodeUpdateCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
double H
Height of the undeformed wall above y=0.
AlgebraicCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
double X0
x position of the undeformed beam's left end.
void set_initial_condition()
Apply initial conditions.
void load(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Load function: Apply a constant external pressure to the wall. Note: This is the load without the flu...
Node * Wall_node_pt
Pointer to control node on the wall.
AlgebraicCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.
double Sigma0
2nd Piola Kirchhoff pre-stress. As in Jensen & Heil (2003) paper.
double Q
Fluid structure interaction parameter: Ratio of stresses used for non-dimensionalisation of fluid to ...
Node * Right_node_pt
Pointer to right control node.
unsigned Nup
Number of elements in the x direction in the upstream part of the channel.
double P_up
Default pressure on the left boundary.
double Lcollapsible
x-length in the collapsible part of the channel
double ReSt
Womersley = Reynolds times Strouhal.
MacroElementNodeUpdateCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc the solution.
FSICollapsibleChannelProblem(const unsigned &nup, const unsigned &ncollapsible, const unsigned &ndown, const unsigned &ny, const double &lup, const double &lcollapsible, const double &ldown, const double &ly)
Constructor: The arguments are the number of elements and the lengths of the domain.
unsigned Ncollapsible
Number of elements in the x direction in the collapsible part of the channel.
double Delta
Boundary layer width.
void actions_after_newton_solve()
Update the problem after solve (empty)
void position(const Vector< double > &zeta, Vector< double > &r) const
Position vector at Lagrangian coordinate zeta.
Mesh * Applied_fluid_traction_mesh_pt
Pointer to the "surface" mesh that applies the traction at the inflow.
Node * Left_node_pt
Pointer to the left control node.
double Fract_in_BL
Fraction of points in boundary layer.
void create_traction_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &traction_mesh_pt)
Create the prescribed traction elements on boundary b.
UndeformedWall(const double &x0, const double &h)
Constructor: arguments are the starting point and the height above y=0.
void prescribed_traction(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Traction applied on the fluid at the left (inflow) boundary.
Collapsible channel mesh with algebraic node update.
int main(int argc, char *argv[])
double H
Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.
double Ldown
x-length in the downstream part of the channel
double Lup
x-length in the upstream part of the channel
virtual CollapsibleChannelDomain::BLSquashFctPt & bl_squash_fct_pt()
Function pointer for function that squashes the mesh near the walls. Default trivial mapping (the ide...
double squash_fct(const double &s)
Mapping [0,1] -> [0,1] that re-distributes nodal points across the channel width. ...
void position(const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
Parametrised position on object: r(zeta). Evaluated at previous timestep. t=0: current time; t>0: pre...
double P_ext
External pressure.
unsigned Ndown
Number of elements in the x direction in the downstream part of the channel.
void d2position(const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
Posn vector and its 1st & 2nd derivatives w.r.t. to coordinates: = drdzeta(alpha,i). = ddrdzeta(alpha,beta,i). Evaluated at current time.
unsigned Ny
Number of elements across the channel.
OneDLagrangianMesh< FSIHermiteBeamElement > * Wall_mesh_pt
Pointer to the "wall" mesh.