fsi_collapsible_channel_adapt.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 
49 
50 //==========start_of_BL_Squash =========================================
51 /// Namespace to define the mapping [0,1] -> [0,1] that re-distributes
52 /// nodal points across the channel width.
53 //======================================================================
54 namespace BL_Squash
55 {
56 
57  /// Boundary layer width
58  double Delta=0.1;
59 
60  /// Fraction of points in boundary layer
61  double Fract_in_BL=0.5;
62 
63  /// \short Mapping [0,1] -> [0,1] that re-distributes
64  /// nodal points across the channel width
65  double squash_fct(const double& s)
66  {
67  // Default return
68  double y=s;
69  if (s<0.5*Fract_in_BL)
70  {
71  y=Delta*2.0*s/Fract_in_BL;
72  }
73  else if (s>1.0-0.5*Fract_in_BL)
74  {
75  y=2.0*Delta/Fract_in_BL*s+1.0-2.0*Delta/Fract_in_BL;
76  }
77  else
78  {
79  y=(1.0-2.0*Delta)/(1.0-Fract_in_BL)*s+
80  (Delta-0.5*Fract_in_BL)/(1.0-Fract_in_BL);
81  }
82 
83  return y;
84  }
85 }// end of BL_Squash
86 
87 
88 
89 //////////////////////////////////////////////////////////////////////////
90 //////////////////////////////////////////////////////////////////////////
91 //////////////////////////////////////////////////////////////////////////
92 
93 
94 //====start_of_underformed_wall============================================
95 /// Undeformed wall is a steady, straight 1D line in 2D space
96 /// \f[ x = X_0 + \zeta \f]
97 /// \f[ y = H \f]
98 //=========================================================================
99 class UndeformedWall : public GeomObject
100 {
101 
102 public:
103 
104  /// \short Constructor: arguments are the starting point and the height
105  /// above y=0.
106  UndeformedWall(const double& x0, const double& h): GeomObject(1,2)
107  {
108  X0=x0;
109  H=h;
110  }
111 
112 
113  /// \short Position vector at Lagrangian coordinate zeta
114  void position(const Vector<double>& zeta, Vector<double>& r) const
115  {
116  // Position Vector
117  r[0] = zeta[0]+X0;
118  r[1] = H;
119  }
120 
121 
122  /// \short Parametrised position on object: r(zeta). Evaluated at
123  /// previous timestep. t=0: current time; t>0: previous
124  /// timestep. Calls steady version.
125  void position(const unsigned& t, const Vector<double>& zeta,
126  Vector<double>& r) const
127  {
128  // Use the steady version
129  position(zeta,r);
130 
131  } // end of position
132 
133 
134  /// \short Posn vector and its 1st & 2nd derivatives
135  /// w.r.t. to coordinates:
136  /// \f$ \frac{dR_i}{d \zeta_\alpha}\f$ = drdzeta(alpha,i).
137  /// \f$ \frac{d^2R_i}{d \zeta_\alpha d \zeta_\beta}\f$ =
138  /// ddrdzeta(alpha,beta,i). Evaluated at current time.
139  virtual void d2position(const Vector<double>& zeta,
140  Vector<double>& r,
141  DenseMatrix<double> &drdzeta,
142  RankThreeTensor<double> &ddrdzeta) const
143  {
144  // Position vector
145  r[0] = zeta[0]+X0;
146  r[1] = H;
147 
148  // Tangent vector
149  drdzeta(0,0)=1.0;
150  drdzeta(0,1)=0.0;
151 
152  // Derivative of tangent vector
153  ddrdzeta(0,0,0)=0.0;
154  ddrdzeta(0,0,1)=0.0;
155 
156  } // end of d2position
157 
158  private :
159 
160  /// x position of the undeformed beam's left end.
161  double X0;
162 
163  /// Height of the undeformed wall above y=0.
164  double H;
165 
166 }; //end_of_undeformed_wall
167 
168 
169 
170 
171 
172 
173 
174 //====start_of_physical_parameters=====================
175 /// Namespace for phyical parameters
176 //======================================================
177 namespace Global_Physical_Variables
178 {
179  /// Reynolds number
180  double Re=50.0;
181 
182  /// Womersley = Reynolds times Strouhal
183  double ReSt=50.0;
184 
185  /// Default pressure on the left boundary
186  double P_up=0.0;
187 
188  /// Traction applied on the fluid at the left (inflow) boundary
189  void prescribed_traction(const double& t,
190  const Vector<double>& x,
191  const Vector<double>& n,
192  Vector<double>& traction)
193  {
194  traction.resize(2);
195  traction[0]=P_up;
196  traction[1]=0.0;
197 
198  } //end traction
199 
200  /// Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.
201  double H=1.0e-2;
202 
203  /// 2nd Piola Kirchhoff pre-stress. As in Jensen & Heil (2003) paper.
204  double Sigma0=1.0e3;
205 
206  /// External pressure
207  double P_ext=0.0;
208 
209  /// \short Load function: Apply a constant external pressure to the wall.
210  /// Note: This is the load without the fluid contribution!
211  /// Fluid load gets added on by FSIWallElement.
212  void load(const Vector<double>& xi, const Vector<double>& x,
213  const Vector<double>& N, Vector<double>& load)
214  {
215  for(unsigned i=0;i<2;i++)
216  {
217  load[i] = -P_ext*N[i];
218  }
219  } //end of load
220 
221 
222  /// \short Fluid structure interaction parameter: Ratio of stresses used for
223  /// non-dimensionalisation of fluid to solid stresses.
224  double Q=1.0e-5;
225 
226 
227 } // end of namespace
228 
229 
230 
231 
232 //====start_of_problem_class==========================================
233 ///Problem class
234 //====================================================================
235 template <class ELEMENT>
236 class FSICollapsibleChannelProblem : public Problem
237 {
238 
239  public :
240 
241 /// \short Constructor: The arguments are the number of elements and
242 /// the lengths of the domain.
243  FSICollapsibleChannelProblem(const unsigned& nup,
244  const unsigned& ncollapsible,
245  const unsigned& ndown,
246  const unsigned& ny,
247  const double& lup,
248  const double& lcollapsible,
249  const double& ldown,
250  const double& ly);
251 
252  /// Destructor (empty)
254 
255 #ifdef MACRO_ELEMENT_NODE_UPDATE
256 
257  /// Access function for the specific bulk (fluid) mesh
259  {
260  // Upcast from pointer to the Mesh base class to the specific
261  // element type that we're using here.
262  return dynamic_cast<
264  (Bulk_mesh_pt);
265  }
266 
267 #else
268 
269  /// Access function for the specific bulk (fluid) mesh
271  {
272  // Upcast from pointer to the Mesh base class to the specific
273  // element type that we're using here.
274  return dynamic_cast<
276  (Bulk_mesh_pt);
277  }
278 
279 #endif
280 
281 
282  /// Access function for the wall mesh
284  {
285  return Wall_mesh_pt;
286 
287  } // end of access to wall mesh
288 
289 
290  /// Actions before adapt: Wipe the mesh of prescribed traction elements
291  void actions_before_adapt();
292 
293  /// Actions after adapt: Rebuild the mesh of prescribed traction elements
294  /// and reset FSI
295  void actions_after_adapt();
296 
297  /// \short Update the problem specs before solve (empty)
299 
300  /// Update the problem after solve (empty)
302 
303  /// \short Update before checking Newton convergence: Update the
304  /// nodal positions in the fluid mesh in response to possible
305  /// changes in the wall shape
307  {
308  Bulk_mesh_pt->node_update();
309  }
310 
311  /// Doc the solution
312  void doc_solution(DocInfo& doc_info,ofstream& trace_file);
313 
314  /// Apply initial conditions
315  void set_initial_condition();
316 
317 private :
318 
319  /// Create the prescribed traction elements on boundary b
320  void create_traction_elements(const unsigned &b,
321  Mesh* const &bulk_mesh_pt,
322  Mesh* const &traction_mesh_pt);
323 
324  /// Delete prescribed traction elements from the surface mesh
325  void delete_traction_elements(Mesh* const &traction_mesh_pt);
326 
327  ///Number of elements in the x direction in the upstream part of the channel
328  unsigned Nup;
329 
330  /// \short Number of elements in the x direction in the collapsible part of
331  /// the channel
332  unsigned Ncollapsible;
333 
334  ///Number of elements in the x direction in the downstream part of the channel
335  unsigned Ndown;
336 
337  ///Number of elements across the channel
338  unsigned Ny;
339 
340  ///x-length in the upstream part of the channel
341  double Lup;
342 
343  ///x-length in the collapsible part of the channel
344  double Lcollapsible;
345 
346  ///x-length in the downstream part of the channel
347  double Ldown;
348 
349  ///Transverse length
350  double Ly;
351 
352 #ifdef MACRO_ELEMENT_NODE_UPDATE
353 
354  /// Pointer to the "bulk" mesh
356 
357 #else
358 
359  /// Pointer to the "bulk" mesh
361 
362 #endif
363 
364  /// \short Pointer to the "surface" mesh that applies the traction at the
365  /// inflow
366  Mesh* Applied_fluid_traction_mesh_pt;
367 
368  /// Pointer to the "wall" mesh
370 
371  ///Pointer to the left control node
372  Node* Left_node_pt;
373 
374  ///Pointer to right control node
375  Node* Right_node_pt;
376 
377  /// Pointer to control node on the wall
378  Node* Wall_node_pt;
379 
380 };//end of problem class
381 
382 
383 
384 
385 //=====start_of_constructor======================================
386 /// Constructor for the collapsible channel problem
387 //===============================================================
388 template <class ELEMENT>
390  const unsigned& nup,
391  const unsigned& ncollapsible,
392  const unsigned& ndown,
393  const unsigned& ny,
394  const double& lup,
395  const double& lcollapsible,
396  const double& ldown,
397  const double& ly)
398 {
399  // Store problem parameters
400  Nup=nup;
401  Ncollapsible=ncollapsible;
402  Ndown=ndown;
403  Ny=ny;
404  Lup=lup;
405  Lcollapsible=lcollapsible;
406  Ldown=ldown;
407  Ly=ly;
408 
409 
410  // Overwrite maximum allowed residual to accomodate bad initial guesses
411  Problem::Max_residuals=1000.0;
412 
413  // Allocate the timestepper for the Navier-Stokes equations
414  BDF<2>* fluid_time_stepper_pt=new BDF<2>;
415 
416  // Add the fluid timestepper to the Problem's collection of timesteppers.
417  add_time_stepper_pt(fluid_time_stepper_pt);
418 
419  // Create a dummy Steady timestepper that stores two history values
420  Steady<2>* wall_time_stepper_pt = new Steady<2>;
421 
422  // Add the wall timestepper to the Problem's collection of timesteppers.
423  add_time_stepper_pt(wall_time_stepper_pt);
424 
425  // Geometric object that represents the undeformed wall:
426  // A straight line at height y=ly; starting at x=lup.
427  UndeformedWall* undeformed_wall_pt=new UndeformedWall(lup,ly);
428 
429  //Create the "wall" mesh with FSI Hermite beam elements, passing the
430  //dummy wall timestepper to the constructor
432  (Ncollapsible,Lcollapsible,undeformed_wall_pt,wall_time_stepper_pt);
433 
434  // Build a geometric object (one Lagrangian, two Eulerian coordinates)
435  // from the wall mesh
436  MeshAsGeomObject* wall_geom_object_pt=
437  new MeshAsGeomObject(Wall_mesh_pt);
438 
439 #ifdef MACRO_ELEMENT_NODE_UPDATE
440 
441  //Build bulk (fluid) mesh
442  Bulk_mesh_pt =
444  (nup, ncollapsible, ndown, ny,
445  lup, lcollapsible, ldown, ly,
446  wall_geom_object_pt,
447  fluid_time_stepper_pt);
448 
449  // Set a non-trivial boundary-layer-squash function...
450  Bulk_mesh_pt->bl_squash_fct_pt() = &BL_Squash::squash_fct;
451 
452  // ... and update the nodal positions accordingly
453  Bulk_mesh_pt->node_update();
454 
455 #else
456 
457  //Build bulk (fluid) mesh
458  Bulk_mesh_pt =
460  (nup, ncollapsible, ndown, ny,
461  lup, lcollapsible, ldown, ly,
462  wall_geom_object_pt,
464  fluid_time_stepper_pt);
465 
466 #endif
467 
468  if (CommandLineArgs::Argc>1)
469  {
470  // One round of uniform refinement
471  Bulk_mesh_pt->refine_uniformly();
472  }
473 
474  // Create "surface mesh" that will contain only the prescribed-traction
475  // elements. The constructor just creates the mesh without
476  // giving it any elements, nodes, etc.
477  Applied_fluid_traction_mesh_pt = new Mesh;
478 
479  // Create prescribed-traction elements from all elements that are
480  // adjacent to boundary 5 (left boundary), but add them to a separate mesh.
481  create_traction_elements(5,Bulk_mesh_pt,Applied_fluid_traction_mesh_pt);
482 
483  // Add the sub meshes to the problem
484  add_sub_mesh(Bulk_mesh_pt);
485  add_sub_mesh(Applied_fluid_traction_mesh_pt);
486  add_sub_mesh(Wall_mesh_pt);
487 
488  // Combine all submeshes into a single Mesh
489  build_global_mesh();
490 
491  //Set errror estimator
492  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
493  bulk_mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
494 
495 
496  // Complete build of fluid mesh
497  //-----------------------------
498 
499  // Loop over the elements to set up element-specific
500  // things that cannot be handled by constructor
501  unsigned n_element=Bulk_mesh_pt->nelement();
502  for(unsigned e=0;e<n_element;e++)
503  {
504  // Upcast from GeneralisedElement to the present element
505  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
506 
507  //Set the Reynolds number
508  el_pt->re_pt() = &Global_Physical_Variables::Re;
509 
510  // Set the Womersley number
511  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
512 
513  } // end loop over elements
514 
515 
516  // Pin redudant pressure dofs
517  RefineableNavierStokesEquations<2>::
518  pin_redundant_nodal_pressures(Bulk_mesh_pt->element_pt());
519 
520 
521  // Apply boundary conditions for fluid
522  //------------------------------------
523 
524 
525  //Pin the velocity on the boundaries
526  //x and y-velocities pinned along boundary 0 (bottom boundary) :
527  unsigned ibound=0;
528  unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
529  for (unsigned inod=0;inod<num_nod;inod++)
530  {
531  for(unsigned i=0;i<2;i++)
532  {
533  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
534  }
535  }
536 
537  //x and y-velocities pinned along boundaries 2, 3, 4 (top boundaries) :
538  for(ibound=2;ibound<5;ibound++)
539  {
540  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
541  for (unsigned inod=0;inod<num_nod;inod++)
542  {
543  for(unsigned i=0;i<2;i++)
544  {
545  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
546  }
547  }
548  }
549 
550  //y-velocity pinned along boundary 1 (right boundary):
551  ibound=1;
552  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
553  for (unsigned inod=0;inod<num_nod;inod++)
554  {
555  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
556  }
557 
558 
559  //y-velocity pinned along boundary 5 (left boundary):
560  ibound=5;
561  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
562  for (unsigned inod=0;inod<num_nod;inod++)
563  {
564 
565  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
566 
567  }//end of pin_velocity
568 
569 
570 
571  // Complete build of applied traction elements
572  //--------------------------------------------
573 
574  // Loop over the traction elements to pass pointer to prescribed
575  // traction function
576  unsigned n_el=Applied_fluid_traction_mesh_pt->nelement();
577  for(unsigned e=0;e<n_el;e++)
578  {
579  // Upcast from GeneralisedElement to NavierStokes traction element
580  NavierStokesTractionElement<ELEMENT> *el_pt =
581  dynamic_cast< NavierStokesTractionElement<ELEMENT>*>(
582  Applied_fluid_traction_mesh_pt->element_pt(e));
583 
584  // Set the pointer to the prescribed traction function
585  el_pt->traction_fct_pt() = &Global_Physical_Variables::prescribed_traction;
586 
587  }
588 
589 
590 
591 
592  // Complete build of wall elements
593  //--------------------------------
594 
595  //Loop over the elements to set physical parameters etc.
596  n_element = wall_mesh_pt()->nelement();
597  for(unsigned e=0;e<n_element;e++)
598  {
599  // Upcast to the specific element type
600  FSIHermiteBeamElement *elem_pt =
601  dynamic_cast<FSIHermiteBeamElement*>(wall_mesh_pt()->element_pt(e));
602 
603  // Set physical parameters for each element:
604  elem_pt->sigma0_pt() = &Global_Physical_Variables::Sigma0;
605  elem_pt->h_pt() = &Global_Physical_Variables::H;
606 
607  // Set the load vector for each element
608  elem_pt->load_vector_fct_pt() = &Global_Physical_Variables::load;
609 
610  // Function that specifies the load ratios
611  elem_pt->q_pt() = &Global_Physical_Variables::Q;
612 
613  // Set the undeformed shape for each element
614  elem_pt->undeformed_beam_pt() = undeformed_wall_pt;
615 
616  // The normal on the wall elements as computed by the FSIHermiteElements
617  // points away from the fluid rather than into the fluid (as assumed
618  // by default)
619  elem_pt->set_normal_pointing_out_of_fluid();
620 
621  } // end of loop over elements
622 
623 
624 
625  // Boundary conditions for wall mesh
626  //----------------------------------
627 
628  // Set the boundary conditions: Each end of the beam is fixed in space
629  // Loop over the boundaries (ends of the beam)
630  for(unsigned b=0;b<2;b++)
631  {
632  // Pin displacements in both x and y directions
633  wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(0);
634  wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(1);
635  }
636 
637 
638 
639  //Choose control nodes
640  //---------------------
641 
642  // Left boundary
643  ibound=5;
644  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
645  unsigned control_nod=num_nod/2;
646  Left_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
647 
648  // Right boundary
649  ibound=1;
650  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
651  control_nod=num_nod/2;
652  Right_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
653 
654 
655  // Set the pointer to the control node on the wall
656  num_nod= wall_mesh_pt()->nnode();
657  Wall_node_pt=wall_mesh_pt()->node_pt(Ncollapsible/2);
658 
659 
660  // Setup FSI
661  //----------
662 
663  // The velocity of the fluid nodes on the wall (fluid mesh boundary 3)
664  // is set by the wall motion -- hence the no-slip condition needs to be
665  // re-applied whenever a node update is performed for these nodes.
666  // Such tasks may be performed automatically by the auxiliary node update
667  // function specified by a function pointer:
668  ibound=3;
669  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
670  for (unsigned inod=0;inod<num_nod;inod++)
671  {
672  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->
673  set_auxiliary_node_update_fct_pt(
674  FSI_functions::apply_no_slip_on_moving_wall);
675  }
676 
677 
678 
679  // Work out which fluid dofs affect the residuals of the wall elements:
680  // We pass the boundary between the fluid and solid meshes and
681  // pointers to the meshes. The interaction boundary is boundary 3 of the
682  // 2D fluid mesh.
683  FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
684  (this,3,Bulk_mesh_pt,Wall_mesh_pt);
685 
686  // Setup equation numbering scheme
687  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
688 
689 
690 }//end of constructor
691 
692 
693 
694 
695 //====start_of_doc_solution===================================================
696 /// Doc the solution
697 //============================================================================
698 template <class ELEMENT>
700  ofstream& trace_file)
701 {
702 
703 #ifdef MACRO_ELEMENT_NODE_UPDATE
704 
705  // Doc fsi
706  if (CommandLineArgs::Argc>1)
707  {
708  FSI_functions::doc_fsi<MacroElementNodeUpdateNode>
709  (Bulk_mesh_pt,Wall_mesh_pt,doc_info);
710  }
711 
712 #else
713 
714  // Doc fsi
715  if (CommandLineArgs::Argc>1)
716  {
717  FSI_functions::doc_fsi<AlgebraicNode>(Bulk_mesh_pt,Wall_mesh_pt,doc_info);
718  }
719 
720 #endif
721 
722  ofstream some_file;
723  char filename[100];
724 
725  // Number of plot points
726  unsigned npts;
727  npts=5;
728 
729  // Output fluid solution
730  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
731  doc_info.number());
732  some_file.open(filename);
733  bulk_mesh_pt()->output(some_file,npts);
734  some_file.close();
735 
736  // Document the wall shape
737  sprintf(filename,"%s/beam%i.dat",doc_info.directory().c_str(),
738  doc_info.number());
739  some_file.open(filename);
740  wall_mesh_pt()->output(some_file,npts);
741  some_file.close();
742 
743  // Loop over all elements do dump out previous solutions
744  // (get the number of previous timesteps available from the wall
745  // timestepper)
746  unsigned nsteps=time_stepper_pt(1)->nprev_values();
747  for (unsigned t=0;t<=nsteps;t++)
748  {
749  sprintf(filename,"%s/wall%i-%i.dat",doc_info.directory().c_str(),
750  doc_info.number(),t);
751  some_file.open(filename);
752  unsigned n_elem=wall_mesh_pt()->nelement();
753  for (unsigned ielem=0;ielem<n_elem;ielem++)
754  {
755  dynamic_cast<FSIHermiteBeamElement*>(wall_mesh_pt()->element_pt(ielem))->
756  output(t,some_file,npts);
757  }
758  some_file.close();
759  } // end of output of previous solutions
760 
761 
762 
763  // Write trace file
764  trace_file << time_pt()->time() << " "
765  << Wall_node_pt->x(1) << " "
766  << Left_node_pt->value(0) << " "
767  << Right_node_pt->value(0) << " "
769  << std::endl;
770 
771 } // end_of_doc_solution
772 
773 
774 
775 
776 
777 //=====start_of_create_traction_elements======================================
778 /// Create the traction elements
779 //============================================================================
780 template <class ELEMENT>
782  const unsigned &b, Mesh* const &bulk_mesh_pt, Mesh* const &traction_mesh_pt)
783 {
784 
785  // How many bulk elements are adjacent to boundary b?
786  unsigned n_element = bulk_mesh_pt->nboundary_element(b);
787 
788  // Loop over the bulk elements adjacent to boundary b?
789  for(unsigned e=0;e<n_element;e++)
790  {
791  // Get pointer to the bulk element that is adjacent to boundary b
792  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>
793  (bulk_mesh_pt->boundary_element_pt(b,e));
794 
795  //What is the index of the face of element e along boundary
796  int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
797 
798  // Build the corresponding prescribed-traction element
799  NavierStokesTractionElement<ELEMENT>* flux_element_pt =
800  new NavierStokesTractionElement<ELEMENT>(bulk_elem_pt,face_index);
801 
802  //Add the prescribed-traction element to the surface mesh
803  traction_mesh_pt->add_element_pt(flux_element_pt);
804 
805  } //end of loop over bulk elements adjacent to boundary b
806 
807 } // end of create_traction_elements
808 
809 
810 
811 //============start_of_delete_traction_elements==============================
812 /// Delete traction elements and wipe the surface mesh
813 //=======================================================================
814 template<class ELEMENT>
816 delete_traction_elements(Mesh* const &surface_mesh_pt)
817 {
818  // How many surface elements are in the surface mesh
819  unsigned n_element = surface_mesh_pt->nelement();
820 
821  // Loop over the surface elements
822  for(unsigned e=0;e<n_element;e++)
823  {
824  // Kill surface element
825  delete surface_mesh_pt->element_pt(e);
826  }
827 
828  // Wipe the mesh
829  surface_mesh_pt->flush_element_and_node_storage();
830 
831 } // end of delete_traction_elements
832 
833 
834 
835 //====start_of_apply_initial_condition========================================
836 /// Apply initial conditions
837 //============================================================================
838 template <class ELEMENT>
840 {
841  // Check that timestepper is from the BDF family
842  if (time_stepper_pt()->type()!="BDF")
843  {
844  std::ostringstream error_stream;
845  error_stream << "Timestepper has to be from the BDF family!\n"
846  << "You have specified a timestepper from the "
847  << time_stepper_pt()->type() << " family" << std::endl;
848 
849  throw OomphLibError(error_stream.str(),
850  OOMPH_CURRENT_FUNCTION,
851  OOMPH_EXCEPTION_LOCATION);
852  }
853 
854  // Update the mesh
855  bulk_mesh_pt()->node_update();
856 
857  // Loop over the nodes to set initial guess everywhere
858  unsigned num_nod = bulk_mesh_pt()->nnode();
859  for (unsigned n=0;n<num_nod;n++)
860  {
861  // Get nodal coordinates
862  Vector<double> x(2);
863  x[0]=bulk_mesh_pt()->node_pt(n)->x(0);
864  x[1]=bulk_mesh_pt()->node_pt(n)->x(1);
865 
866  // Assign initial condition: Steady Poiseuille flow
867  bulk_mesh_pt()->node_pt(n)->set_value(0,6.0*(x[1]/Ly)*(1.0-(x[1]/Ly)));
868  bulk_mesh_pt()->node_pt(n)->set_value(1,0.0);
869  }
870 
871  // Assign initial values for an impulsive start
872  bulk_mesh_pt()->assign_initial_values_impulsive();
873  wall_mesh_pt()->assign_initial_values_impulsive();
874 
875 } // end of set_initial_condition
876 
877 
878 
879 
880 
881 //=========start_of_actions_before_adapt==================================
882 /// Actions before adapt: Wipe the mesh of prescribed traction elements
883 //========================================================================
884 template<class ELEMENT>
886 {
887  // Kill the traction elements and wipe surface mesh
888  delete_traction_elements(Applied_fluid_traction_mesh_pt);
889 
890  // Rebuild the global mesh.
891  rebuild_global_mesh();
892 
893 } // end of actions_before_adapt
894 
895 
896 
897 //==========start_of_actions_after_adapt==================================
898 /// Actions after adapt: Rebuild the mesh of prescribed traction elements
899 //========================================================================
900 template<class ELEMENT>
902 {
903  // Create prescribed-flux elements from all elements that are
904  // adjacent to boundary 5 and add them to surface mesh
905  create_traction_elements(5,Bulk_mesh_pt,Applied_fluid_traction_mesh_pt);
906 
907  // Rebuild the global mesh
908  rebuild_global_mesh();
909 
910  // Unpin all pressure dofs
911  RefineableNavierStokesEquations<2>::
912  unpin_all_pressure_dofs(Bulk_mesh_pt->element_pt());
913 
914  // Pin redundant pressure dofs
915  RefineableNavierStokesEquations<2>::
916  pin_redundant_nodal_pressures(Bulk_mesh_pt->element_pt());
917 
918  // Loop over the traction elements to pass pointer to prescribed
919  // traction function
920  unsigned n_element=Applied_fluid_traction_mesh_pt->nelement();
921  for(unsigned e=0;e<n_element;e++)
922  {
923  // Upcast from GeneralisedElement to NavierStokesTractionElement element
924  NavierStokesTractionElement<ELEMENT> *el_pt =
925  dynamic_cast<NavierStokesTractionElement<ELEMENT>*>(
926  Applied_fluid_traction_mesh_pt->element_pt(e));
927 
928  // Set the pointer to the prescribed traction function
929  el_pt->traction_fct_pt() = &Global_Physical_Variables::prescribed_traction;
930  }
931 
932  // The functions used to update the no slip boundary conditions
933  // must be set on any new nodes that have been created during the
934  // mesh adaptation process.
935  // There is no mechanism by which auxiliary update functions
936  // are copied to newly created nodes.
937  // (because, unlike boundary conditions, they don't occur exclusively
938  // at boundaries)
939 
940  // The velocity of the fluid nodes on the wall (fluid mesh boundary 3)
941  // is set by the wall motion -- hence the no-slip condition needs to be
942  // re-applied whenever a node update is performed for these nodes.
943  // Such tasks may be performed automatically by the auxiliary node update
944  // function specified by a function pointer:
945  unsigned ibound=3;
946  unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
947  for (unsigned inod=0;inod<num_nod;inod++)
948  {
949  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->
950  set_auxiliary_node_update_fct_pt(
951  FSI_functions::apply_no_slip_on_moving_wall);
952  }
953 
954  // (Re-)setup fsi: Work out which fluid dofs affect wall elements
955  // the correspondance between wall dofs and fluid elements is handled
956  // during the remeshing, but the "reverse" association must be done
957  // separately. We need to set up the interaction every time because the fluid
958  // element adjacent to a given solid element's integration point may have
959  // changed.We pass the boundary between the fluid and solid meshes and
960  // pointers to the meshes. The interaction boundary is boundary 3 of
961  // the Fluid mesh.
962  FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
963  (this,3,Bulk_mesh_pt,Wall_mesh_pt);
964 
965 
966 } // end of actions_after_adapt
967 
968 
969 
970 //============start_of_main====================================================
971 /// Driver code for a collapsible channel problem with FSI.
972 /// Presence of command line arguments indicates validation run with
973 /// coarse resolution and small number of timesteps.
974 //=============================================================================
975 int main(int argc, char* argv[])
976 {
977 
978  // Store command line arguments
979  CommandLineArgs::setup(argc,argv);
980 
981  // Reduction in resolution for validation run?
982  unsigned coarsening_factor=4;
983  if (CommandLineArgs::Argc>1)
984  {
985  coarsening_factor=4;
986  }
987 
988  // Number of elements in the domain
989  unsigned nup=20/coarsening_factor;
990  unsigned ncollapsible=40/coarsening_factor;
991  unsigned ndown=40/coarsening_factor;
992  unsigned ny=16/coarsening_factor;
993 
994  // Length of the domain
995  double lup=5.0;
996  double lcollapsible=10.0;
997  double ldown=10.0;
998  double ly=1.0;
999 
1000  // Set external pressure (on the wall stiffness scale).
1002 
1003  // Pressure on the left boundary: This is consistent with steady
1004  // Poiseuille flow
1005  Global_Physical_Variables::P_up=12.0*(lup+lcollapsible+ldown);
1006 
1007 
1008 #ifdef MACRO_ELEMENT_NODE_UPDATE
1009 
1010 #ifdef TAYLOR_HOOD
1011 
1012  // Build the problem with QTaylorHoodElements
1014  <MacroElementNodeUpdateElement<RefineableQTaylorHoodElement<2> > >
1015  problem(nup, ncollapsible, ndown, ny,
1016  lup, lcollapsible, ldown, ly);
1017 
1018 #else
1019 
1020  // Build the problem with QCrouzeixRaviartElements
1022  <MacroElementNodeUpdateElement<RefineableQCrouzeixRaviartElement<2> > >
1023  problem(nup, ncollapsible, ndown, ny,
1024  lup, lcollapsible, ldown, ly);
1025 
1026 #endif
1027 
1028 #else
1029 
1030 #ifdef TAYLOR_HOOD
1031 
1032  // Build the problem with QTaylorHoodElements
1034  <AlgebraicElement<RefineableQTaylorHoodElement<2> > >
1035  problem(nup, ncollapsible, ndown, ny,
1036  lup, lcollapsible, ldown, ly);
1037 
1038 #else
1039 
1040  // Build the problem with QCrouzeixRaviartElements
1042  <AlgebraicElement<RefineableQCrouzeixRaviartElement<2> > >
1043  problem(nup, ncollapsible, ndown, ny,
1044  lup, lcollapsible, ldown, ly);
1045 
1046 #endif
1047 
1048 #endif
1049 
1050 
1051  // Timestep. Note: Preliminary runs indicate that the period of
1052  // the oscillation is about 1 so this gives us 40 steps per period.
1053  double dt=1.0/40.0;
1054 
1055  // Initial time for the simulation
1056  double t_min=0.0;
1057 
1058  // Maximum time for simulation
1059  double t_max=3.5;
1060 
1061  // Initialise timestep
1062  problem.time_pt()->time()=t_min;
1063  problem.initialise_dt(dt);
1064 
1065  // Apply initial condition
1066  problem.set_initial_condition();
1067 
1068  //Set output directory
1069  DocInfo doc_info;
1070  doc_info.set_directory("RESLT");
1071 
1072  // Open a trace file
1073  ofstream trace_file;
1074  char filename[100];
1075  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
1076  trace_file.open(filename);
1077 
1078  // Output the initial condition
1079  problem.doc_solution(doc_info, trace_file);
1080 
1081  // Increment step number
1082  doc_info.number()++;
1083 
1084  // Find number of timesteps (reduced for validation)
1085  unsigned nstep = unsigned((t_max-t_min)/dt);
1086  if (CommandLineArgs::Argc>1)
1087  {
1088  nstep=3;
1089  }
1090 
1091 
1092  // Set targets for spatial adaptivity
1093  problem.bulk_mesh_pt()->max_permitted_error()=1.0e-3;
1094  problem.bulk_mesh_pt()->min_permitted_error()=1.0e-5;
1095 
1096  // Overwrite for validation run
1097  if (CommandLineArgs::Argc>1)
1098  {
1099  // Set targets for spatial adaptivity
1100  problem.bulk_mesh_pt()->max_permitted_error()=0.5e-2;
1101  problem.bulk_mesh_pt()->min_permitted_error()=0.5e-4;
1102  }
1103 
1104  // When performing the first timestep, we can adapt the mesh as many times
1105  // as we want because the initial condition can be re-set
1106  unsigned max_adapt=3;
1107  bool first=true;
1108 
1109  // Timestepping loop
1110  for (unsigned istep=0;istep<nstep;istep++)
1111  {
1112  // Solve the problem
1113  problem.unsteady_newton_solve(dt,max_adapt,first);
1114 
1115  // Outpt the solution
1116  problem.doc_solution(doc_info, trace_file);
1117 
1118  // Step number
1119  doc_info.number()++;
1120 
1121  // We've done the first step
1122  first=false;
1123  max_adapt=1;
1124  }
1125 
1126 
1127  // Close trace file.
1128  trace_file.close();
1129 
1130 }//end of main
1131 
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)
RefineableAlgebraicCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
void set_initial_condition()
Apply initial conditions.
MacroElementNodeUpdateRefineableCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
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...
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 ...
double P_up
Default pressure on the left boundary.
double ReSt
Womersley = Reynolds times Strouhal.
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.
virtual 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.
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.
int main(int argc, char *argv[])
void actions_before_adapt()
Actions before adapt: Wipe the mesh of prescribed traction elements.
double Fract_in_BL
Fraction of points in boundary layer.
RefineableAlgebraicCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.
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.
double H
Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.
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.
void delete_traction_elements(Mesh *const &traction_mesh_pt)
Delete prescribed traction elements from the surface mesh.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of prescribed traction elements.
MacroElementNodeUpdateRefineableCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.