unstructured_three_d_fsi.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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 //LIC//
8 //LIC// This library is free software; you can redistribute it and/or
9 //LIC// modify it under the terms of the GNU Lesser General Public
10 //LIC// License as published by the Free Software Foundation; either
11 //LIC// version 2.1 of the License, or (at your option) any later version.
12 //LIC//
13 //LIC// This library is distributed in the hope that it will be useful,
14 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 //LIC// Lesser General Public License for more details.
17 //LIC//
18 //LIC// You should have received a copy of the GNU Lesser General Public
19 //LIC// License along with this library; if not, write to the Free Software
20 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 //LIC// 02110-1301 USA.
22 //LIC//
23 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 //LIC//
25 //LIC//====================================================================
26 // Driver code for a simple unstructured fsi problem using meshes
27 // generated from input files generated by the 3d mesh generator
28 // tetgen.
29 
30 //Generic libraries
31 #include "generic.h"
32 #include "solid.h"
33 #include "constitutive.h"
34 #include "navier_stokes.h"
35 
36 // Get the mesh
37 #include "meshes/tetgen_mesh.h"
38 
39 using namespace std;
40 using namespace oomph;
41 
42 
43 
44 //==========start_solid_mesh===============================================
45 /// Tetgen-based mesh upgraded to become a solid mesh
46 //=========================================================================
47 template<class ELEMENT>
48 class MySolidTetgenMesh : public virtual TetgenMesh<ELEMENT>,
49  public virtual SolidMesh
50 {
51 
52 public:
53 
54  /// Constructor:
55  MySolidTetgenMesh(const std::string& node_file_name,
56  const std::string& element_file_name,
57  const std::string& face_file_name,
58  TimeStepper* time_stepper_pt=
59  &Mesh::Default_TimeStepper) :
60  TetgenMesh<ELEMENT>(node_file_name, element_file_name,
61  face_file_name, time_stepper_pt)
62  {
63  //Assign the Lagrangian coordinates
64  set_lagrangian_nodal_coordinates();
65 
66  // Find elements next to boundaries
67  setup_boundary_element_info();
68 
69  // Setup boundary coordinates for all boundaries
70  char filename[100];
71  ofstream some_file;
72  unsigned nb=this->nboundary();
73  for (unsigned b=0;b<nb;b++)
74  {
75  sprintf(filename,"RESLT/solid_boundary_test%i.dat",b);
76  some_file.open(filename);
77  this->template setup_boundary_coordinates<ELEMENT>(b,some_file);
78  some_file.close();
79  }
80 
81  }
82 
83  /// Empty Destructor
84  virtual ~MySolidTetgenMesh() { }
85 
86 };
87 
88 /// ////////////////////////////////////////////////////////////////////
89 /// ////////////////////////////////////////////////////////////////////
90 /// ////////////////////////////////////////////////////////////////////
91 
92 
93 
94 //==============start_fluid_mesh===========================================
95 /// Tetgen-based mesh upgraded to become a (pseudo-) solid mesh
96 //=========================================================================
97 template<class ELEMENT>
98 class FluidTetMesh : public virtual TetgenMesh<ELEMENT>,
99  public virtual SolidMesh
100 {
101 
102 public:
103 
104  /// Constructor:
105  FluidTetMesh(const std::string& node_file_name,
106  const std::string& element_file_name,
107  const std::string& face_file_name,
108  const bool& split_corner_elements,
109  TimeStepper* time_stepper_pt=
110  &Mesh::Default_TimeStepper) :
111  TetgenMesh<ELEMENT>(node_file_name, element_file_name,
112  face_file_name, split_corner_elements,
113  time_stepper_pt)
114  {
115  //Assign the Lagrangian coordinates
116  set_lagrangian_nodal_coordinates();
117 
118  // Find out elements next to boundary
119  setup_boundary_element_info();
120 
121  // Setup boundary coordinates for boundary.
122  // To be consistent with the boundary coordinates generated
123  // in the solid, we switch the direction of the normal.
124  // (Both meshes are generated from the same polygonal facets
125  // at the FSI interface).
126  bool switch_normal=true;
127 
128  // Setup boundary coordinates for all boundaries
129  char filename[100];
130  ofstream some_file;
131  unsigned nb=this->nboundary();
132  for (unsigned b=0;b<nb;b++)
133  {
134  sprintf(filename,"RESLT/fluid_boundary_test%i.dat",b);
135  some_file.open(filename);
136  this->template setup_boundary_coordinates<ELEMENT>(b,switch_normal,some_file);
137  some_file.close();
138  }
139 
140  }
141 
142  /// Empty Destructor
143  virtual ~FluidTetMesh() { }
144 
145 };
146 
147 
148 /// ///////////////////////////////////////////////////////////////
149 /// ///////////////////////////////////////////////////////////////
150 /// ///////////////////////////////////////////////////////////////
151 
152 
153 //=======start_of_namespace==========================================
154 /// Global variables
155 //================================================================
157 {
158 
159  /// Default Reynolds number
160  double Re=100.0;
161 
162  /// Default FSI parameter
163  double Q=0.0;
164 
165  /// Pointer to constitutive law
166  ConstitutiveLaw* Constitutive_law_pt=0;
167 
168  /// Poisson's ratio for generalised Hookean constitutive equation
169  double Nu=0.3;
170 
171  /// Fluid pressure on inflow boundary
172  double P_in=0.5;
173 
174  /// Applied traction on fluid at the inflow boundary
175  void prescribed_inflow_traction(const double& t,
176  const Vector<double>& x,
177  const Vector<double>& n,
178  Vector<double>& traction)
179  {
180  traction[0]=0.0;
181  traction[1]=0.0;
182  traction[2]=P_in;
183  }
184 
185 
186  /// Fluid pressure on outflow boundary
187  double P_out=-0.5;
188 
189  /// Applied traction on fluid at the inflow boundary
190  void prescribed_outflow_traction(const double& t,
191  const Vector<double>& x,
192  const Vector<double>& n,
193  Vector<double>& traction)
194  {
195  traction[0]=0.0;
196  traction[1]=0.0;
197  traction[2]=-P_out;
198  }
199 
200 
201 } //end_of_namespace
202 
203 
204 
205 
206 
207 
208 //===============start_of_problem_class===============================
209 /// Unstructured 3D FSI problem
210 //====================================================================
211 template<class FLUID_ELEMENT, class SOLID_ELEMENT>
212 class UnstructuredFSIProblem : public Problem
213 {
214 
215 public:
216 
217  /// Constructor:
219 
220  /// Destructor (empty)
222 
223  /// Doc the solution
224  void doc_solution(DocInfo& doc_info);
225 
226  /// Create fluid traction elements at inflow
227  void create_fluid_traction_elements();
228 
229  /// Create FSI traction elements
230  void create_fsi_traction_elements();
231 
232  /// Create elements that enforce prescribed boundary motion
233  /// for the pseudo-solid fluid mesh by Lagrange multipliers
234  void create_lagrange_multiplier_elements();
235 
236 
237 private:
238 
239  /// Sanity check: Doc boundary coordinates on i-th solid FSI interface
240  void doc_solid_boundary_coordinates(const unsigned& i);
241 
242  /// Return total number of mesh boundaries that make up the inflow
243  /// boundary
245  {return Inflow_boundary_id.size();}
246 
247  /// Return total number of mesh boundaries that make up the outflow
248  /// boundary
250  {return Outflow_boundary_id.size();}
251 
252  /// Return total number of mesh boundaries that make up the
253  /// in- and outflow boundaries where a traction has to be applied
255  {return Inflow_boundary_id.size()+Outflow_boundary_id.size();}
256 
257  /// Return total number of mesh boundaries in the solid mesh that
258  /// make up the FSI interface
260  {return Solid_fsi_boundary_id.size();}
261 
262  /// Return total number of mesh boundaries in the fluid mesh that
263  /// make up the FSI interface
265  {return Fluid_fsi_boundary_id.size();}
266 
267  /// Return total number of mesh boundaries in the solid mesh
268  /// where the position is pinned.
270  {return Pinned_solid_boundary_id.size();}
271  //end npinned_solid_boundary
272 
273 
274  /// Bulk solid mesh
276 
277  /// Meshes of FSI traction elements
278  Vector<SolidMesh*> Solid_fsi_traction_mesh_pt;
279 
280  /// Bulk fluid mesh
282 
283  /// Meshes of fluid traction elements that apply pressure at in/outflow
284  Vector<Mesh*> Fluid_traction_mesh_pt;
285 
286  /// Meshes of Lagrange multiplier elements
287  Vector<SolidMesh*> Lagrange_multiplier_mesh_pt;
288 
289  /// GeomObject incarnations of the FSI boundary in the solid mesh
290  Vector<MeshAsGeomObject*>
292 
293  /// IDs of solid mesh boundaries where displacements are pinned
294  Vector<unsigned> Pinned_solid_boundary_id;
295 
296  /// IDs of solid mesh boundaries which make up the FSI interface
297  Vector<unsigned> Solid_fsi_boundary_id;
298 
299  /// IDs of fluid mesh boundaries along which inflow boundary conditions
300  /// are applied
301  Vector<unsigned> Inflow_boundary_id;
302 
303  /// IDs of fluid mesh boundaries along which inflow boundary conditions
304  /// are applied
305  Vector<unsigned> Outflow_boundary_id;
306 
307  /// IDs of fluid mesh boundaries which make up the FSI interface
308  Vector<unsigned> Fluid_fsi_boundary_id;
309 
310 };
311 
312 
313 
314 //==========start_of_constructor==========================================
315 /// Constructor for unstructured 3D FSI problem
316 //========================================================================
317 template<class FLUID_ELEMENT, class SOLID_ELEMENT>
319 {
320  // Define fluid mesh and its distinguished boundaries
321  //---------------------------------------------------
322 
323  //Create fluid bulk mesh, sub-dividing "corner" elements
324  string node_file_name="fsi_bifurcation_fluid.1.node";
325  string element_file_name="fsi_bifurcation_fluid.1.ele";
326  string face_file_name="fsi_bifurcation_fluid.1.face";
327  bool split_corner_elements=true;
328  Fluid_mesh_pt = new FluidTetMesh<FLUID_ELEMENT>(node_file_name,
329  element_file_name,
330  face_file_name,
331  split_corner_elements);
332 
333 
334  // The following corresponds to the boundaries as specified by
335  // facets in the tetgen input:
336 
337  // Fluid mesh has one inflow boundary: Boundary 0
338  Inflow_boundary_id.resize(1);
339  Inflow_boundary_id[0]=0;
340 
341  // Fluid mesh has two outflow boundaries: Boundaries 1 and 2
342  Outflow_boundary_id.resize(2);
343  Outflow_boundary_id[0]=1;
344  Outflow_boundary_id[1]=2;
345 
346  // The remaining fluid boundaries are FSI boundaries.
347  // Note that their order (as indexed in this vector, not
348  // their actual numbers) have to match those in the corresponding
349  // lookup scheme for the solid.
350  Fluid_fsi_boundary_id.resize(12);
351  for (unsigned i=0;i<12;i++)
352  {
353  Fluid_fsi_boundary_id[i]=i+3;
354  }
355 
356 
357  // Define solid mesh and its distinguished boundaries
358  //---------------------------------------------------
359 
360  //Create solid bulk mesh
361  node_file_name="fsi_bifurcation_solid.1.node";
362  element_file_name="fsi_bifurcation_solid.1.ele";
363  face_file_name="fsi_bifurcation_solid.1.face";
364  Solid_mesh_pt = new MySolidTetgenMesh<SOLID_ELEMENT>(node_file_name,
365  element_file_name,
366  face_file_name);
367 
368  // The following corresponds to the boundaries as specified by
369  // facets in the tetgen input:
370 
371  /// IDs of solid mesh boundaries where displacements are pinned
372  Pinned_solid_boundary_id.resize(3);
373  Pinned_solid_boundary_id[0]=0;
374  Pinned_solid_boundary_id[1]=1;
375  Pinned_solid_boundary_id[2]=2;
376 
377  // The solid and fluid fsi boundaries are numbered int he same way.
378  Solid_fsi_boundary_id.resize(12);
379  for (unsigned i=0;i<12;i++)
380  {
381  Solid_fsi_boundary_id[i]=i+3;
382  }
383 
384 
385 
386  // Create (empty) meshes of fluid traction elements at inflow/outflow
387  //-----------------------------------------------------------
388 
389  // Create the meshes
390  unsigned n=nfluid_traction_boundary();
391  Fluid_traction_mesh_pt.resize(n);
392  for (unsigned i=0;i<n;i++)
393  {
394  Fluid_traction_mesh_pt[i]=new Mesh;
395  }
396 
397  // Populate them with elements
398  create_fluid_traction_elements();
399 
400 
401 // Create FSI Traction elements
402 //-----------------------------
403 
404 // Create (empty) meshes of FSI traction elements
405  n=nsolid_fsi_boundary();
406  Solid_fsi_traction_mesh_pt.resize(n);
407  for (unsigned i=0;i<n;i++)
408  {
409  Solid_fsi_traction_mesh_pt[i]=new SolidMesh;
410  }
411 
412  // Build the FSI traction elements
413  create_fsi_traction_elements();
414 
415 
416  // Create Lagrange multiplier mesh for boundary motion of fluid mesh
417  //------------------------------------------------------------------
418 
419  // Construct the mesh of elements that enforce prescribed boundary motion
420  // of pseudo-solid fluid mesh by Lagrange multipliers
421  n=nfluid_fsi_boundary();
422  Lagrange_multiplier_mesh_pt.resize(n);
423  for (unsigned i=0;i<n;i++)
424  {
425  Lagrange_multiplier_mesh_pt[i]=new SolidMesh;
426  }
427 
428  // Create elements
429  create_lagrange_multiplier_elements();
430 
431 
432  // Combine the lot
433  //----------------
434 
435  // Add sub meshes:
436 
437  // The solid bulk mesh
438  add_sub_mesh(Solid_mesh_pt);
439 
440  // Fluid bulk mesh
441  add_sub_mesh(Fluid_mesh_pt);
442 
443  // The fluid traction meshes
444  n=nfluid_traction_boundary();
445  for (unsigned i=0;i<n;i++)
446  {
447  add_sub_mesh(Fluid_traction_mesh_pt[i]);
448  }
449 
450  // The solid fsi traction meshes
451  n=nsolid_fsi_boundary();
452  for (unsigned i=0;i<n;i++)
453  {
454  add_sub_mesh(Solid_fsi_traction_mesh_pt[i]);
455  }
456 
457  // The Lagrange multiplier meshes for the fluid
458  n=nfluid_fsi_boundary();
459  for (unsigned i=0;i<n;i++)
460  {
461  add_sub_mesh(Lagrange_multiplier_mesh_pt[i]);
462  }
463 
464  // Build global mesh
465  build_global_mesh();
466 
467 
468 
469 
470  // Apply BCs for fluid
471  //--------------------
472 
473  // Doc position of pinned pseudo solid nodes
474  std::ofstream pseudo_solid_bc_file("RESLT/pinned_pseudo_solid_nodes.dat");
475 
476  // Loop over inflow/outflow boundaries to impose parallel flow
477  // and pin pseudo-solid displacements
478  for (unsigned in_out=0;in_out<2;in_out++)
479  {
480  // Loop over in/outflow boundaries
481  unsigned n=nfluid_inflow_traction_boundary();
482  if (in_out==1) n=nfluid_outflow_traction_boundary();
483  for (unsigned i=0;i<n;i++)
484  {
485 
486  // Get boundary ID
487  unsigned b=0;
488  if (in_out==0)
489  {
490  b=Inflow_boundary_id[i];
491  }
492  else
493  {
494  b=Outflow_boundary_id[i];
495  }
496 
497  // Number of nodes on that boundary
498  unsigned num_nod=Fluid_mesh_pt->nboundary_node(b);
499  for (unsigned inod=0;inod<num_nod;inod++)
500  {
501  // Get the node
502  SolidNode* nod_pt=Fluid_mesh_pt->boundary_node_pt(b,inod);
503 
504  // Pin transverse velocities
505  nod_pt->pin(0);
506  nod_pt->pin(1);
507 
508  // Pin the nodal (pseudo-solid) displacements
509  for(unsigned i=0;i<3;i++)
510  {
511  nod_pt->pin_position(i);
512 
513  // Doc it as pinned
514  pseudo_solid_bc_file << nod_pt->x(i) << " ";
515  }
516  }
517  }
518  }
519 
520  // Close
521  pseudo_solid_bc_file.close();
522 
523  // Doc bcs for Lagrange multipliers
524  ofstream pinned_file("RESLT/pinned_lagrange_multiplier_nodes.dat");
525 
526  // Loop over all fluid mesh boundaries and pin velocities
527  // of nodes that haven't been dealt with yet
528  unsigned nbound=nfluid_fsi_boundary();
529  for(unsigned i=0;i<nbound;i++)
530  {
531  //Get the mesh boundary
532  unsigned b = Fluid_fsi_boundary_id[i];
533 
534  unsigned num_nod=Fluid_mesh_pt->nboundary_node(b);
535  for (unsigned inod=0;inod<num_nod;inod++)
536  {
537  // Get node
538  Node* nod_pt= Fluid_mesh_pt->boundary_node_pt(b,inod);
539 
540  // Pin all velocities
541  nod_pt->pin(0);
542  nod_pt->pin(1);
543  nod_pt->pin(2);
544 
545  // Find out whether node is also on in/outflow
546  bool is_in_or_outflow_node=false;
547  unsigned n=nfluid_inflow_traction_boundary();
548  for (unsigned k=0;k<n;k++)
549  {
550  if (nod_pt->is_on_boundary(Inflow_boundary_id[k]))
551  {
552  is_in_or_outflow_node=true;
553  break;
554  }
555  }
556  if (!is_in_or_outflow_node)
557  {
558  unsigned n=nfluid_outflow_traction_boundary();
559  for (unsigned k=0;k<n;k++)
560  {
561  if (nod_pt->is_on_boundary(Outflow_boundary_id[k]))
562  {
563  is_in_or_outflow_node=true;
564  break;
565  }
566  }
567  }
568 
569  // Pin the Lagrange multipliers on the out/in-flow boundaries
570  if (is_in_or_outflow_node)
571  {
572  //Cast to a boundary node
573  BoundaryNode<SolidNode> *bnod_pt =
574  dynamic_cast<BoundaryNode<SolidNode>*>
575  ( Fluid_mesh_pt->boundary_node_pt(b,inod) );
576 
577  // Loop over the Lagrange multipliers
578  for (unsigned l=0;l<3;l++)
579  {
580  // Pin the Lagrange multipliers that impose the displacement
581  // because the positon of the fluid nodes at the in/outflow
582  // is already determined.
583  nod_pt->pin
584  (bnod_pt->index_of_first_value_assigned_by_face_element()+l);
585  }
586 
587  // Doc that we've pinned the Lagrange multipliers at this node
588  pinned_file << nod_pt->x(0) << " "
589  << nod_pt->x(1) << " "
590  << nod_pt->x(2) << endl;
591  }
592  }
593 
594  } // done no slip on fsi boundary
595 
596  // Done
597  pinned_file.close();
598 
599  // Complete the build of the fluid elements so they are fully functional
600  //----------------------------------------------------------------------
601  unsigned n_element = Fluid_mesh_pt->nelement();
602  for(unsigned e=0;e<n_element;e++)
603  {
604  // Upcast from GeneralisedElement to the present element
605  FLUID_ELEMENT* el_pt =
606  dynamic_cast<FLUID_ELEMENT*>(Fluid_mesh_pt->element_pt(e));
607 
608  //Set the Reynolds number
609  el_pt->re_pt() = &Global_Parameters::Re;
610 
611  // Set the constitutive law for pseudo-elastic mesh deformation
612  el_pt->constitutive_law_pt() =
614 
615  } // end loop over elements
616 
617 
618 
619  // Apply BCs for solid
620  //--------------------
621 
622  // Doc pinned solid nodes
623  std::ofstream bc_file("RESLT/pinned_solid_nodes.dat");
624 
625  // Pin positions at inflow boundary (boundaries 0 and 1)
626  n=npinned_solid_boundary();
627  for (unsigned i=0;i<n;i++)
628  {
629  // Get boundary ID
630  unsigned b=Pinned_solid_boundary_id[i];
631  unsigned num_nod= Solid_mesh_pt->nboundary_node(b);
632  for (unsigned inod=0;inod<num_nod;inod++)
633  {
634  // Get node
635  SolidNode* nod_pt=Solid_mesh_pt->boundary_node_pt(b,inod);
636 
637  // Pin all directions
638  for (unsigned i=0;i<3;i++)
639  {
640  nod_pt->pin_position(i);
641 
642  // ...and doc it as pinned
643  bc_file << nod_pt->x(i) << " ";
644  }
645 
646  bc_file << std::endl;
647  }
648  }
649  bc_file.close();
650 
651 
652 
653  // Complete the build of Solid elements so they are fully functional
654  //----------------------------------------------------------------
655  n_element = Solid_mesh_pt->nelement();
656  for(unsigned i=0;i<n_element;i++)
657  {
658  //Cast to a solid element
659  SOLID_ELEMENT *el_pt = dynamic_cast<SOLID_ELEMENT*>(
660  Solid_mesh_pt->element_pt(i));
661 
662  // Set the constitutive law
663  el_pt->constitutive_law_pt() =
665  }
666 
667 
668  // Setup FSI
669  //----------
670 
671  // Work out which fluid dofs affect the residuals of the wall elements:
672  // We pass the boundary between the fluid and solid meshes and
673  // pointers to the meshes.
674  n=nsolid_fsi_boundary();
675  for (unsigned i=0;i<n;i++)
676  {
677  // Sanity check: Doc boundary coordinates from solid side
678  doc_solid_boundary_coordinates(i);
679 
680  //Doc boundary coordinates in fluid
681  char filename[100];
682  sprintf(filename,"RESLT/fluid_boundary_coordinates%i.dat",i);
683  Multi_domain_functions::Doc_boundary_coordinate_file.open(filename);
684 
685  // Setup FSI: Pass ID of fluid FSI boundary and associated
686  // mesh of solid fsi traction elements.
687  FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,3>
688  (this,Fluid_fsi_boundary_id[i],Fluid_mesh_pt,Solid_fsi_traction_mesh_pt[i]);
689 
690  // Close the doc file
691  Multi_domain_functions::Doc_boundary_coordinate_file.close();
692  }
693 
694  // Setup equation numbering scheme
695  std::cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
696 
697 }
698 
699 
700 //============start_of_create_fsi_traction_elements======================
701 /// Create FSI traction elements
702 //=======================================================================
703 template<class FLUID_ELEMENT,class SOLID_ELEMENT>
706 {
707 
708  // Loop over FSI boundaries in solid
709  unsigned n=nsolid_fsi_boundary();
710  for (unsigned i=0;i<n;i++)
711  {
712  // Get boundary ID
713  unsigned b=Solid_fsi_boundary_id[i];
714 
715  // How many bulk elements are adjacent to boundary b?
716  unsigned n_element = Solid_mesh_pt->nboundary_element(b);
717 
718  // Loop over the bulk elements adjacent to boundary b
719  for(unsigned e=0;e<n_element;e++)
720  {
721  // Get pointer to the bulk element that is adjacent to boundary b
722  SOLID_ELEMENT* bulk_elem_pt = dynamic_cast<SOLID_ELEMENT*>(
723  Solid_mesh_pt->boundary_element_pt(b,e));
724 
725  //What is the index of the face of the element e along boundary b
726  int face_index = Solid_mesh_pt->face_index_at_boundary(b,e);
727 
728  // Create new element
729  FSISolidTractionElement<SOLID_ELEMENT,3>* el_pt=
730  new FSISolidTractionElement<SOLID_ELEMENT,3>(bulk_elem_pt,face_index);
731 
732  // Add it to the mesh
733  Solid_fsi_traction_mesh_pt[i]->add_element_pt(el_pt);
734 
735  // Specify boundary number
736  el_pt->set_boundary_number_in_bulk_mesh(b);
737 
738  // Function that specifies the load ratios
739  el_pt->q_pt() = &Global_Parameters::Q;
740  }
741  }
742 
743 } // end of create_fsi_traction_elements
744 
745 
746 //============start_of_create_lagrange_multiplier_elements===============
747 /// Create elements that impose the prescribed boundary displacement
748 /// for the pseudo-solid fluid mesh
749 //=======================================================================
750 template<class FLUID_ELEMENT, class SOLID_ELEMENT>
753 {
754  // Make space
755  unsigned n=nfluid_fsi_boundary();
756  Solid_fsi_boundary_pt.resize(n);
757 
758  // Loop over FSI interfaces in fluid
759  for (unsigned i=0;i<n;i++)
760  {
761  // Get boundary ID
762  unsigned b=Fluid_fsi_boundary_id[i];
763 
764  // Create GeomObject incarnation of fsi boundary in solid mesh
765  Solid_fsi_boundary_pt[i]=
766  new MeshAsGeomObject
767  (Solid_fsi_traction_mesh_pt[i]);
768 
769  // How many bulk fluid elements are adjacent to boundary b?
770  unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
771 
772  // Loop over the bulk fluid elements adjacent to boundary b?
773  for(unsigned e=0;e<n_element;e++)
774  {
775  // Get pointer to the bulk fluid element that is adjacent to boundary b
776  FLUID_ELEMENT* bulk_elem_pt = dynamic_cast<FLUID_ELEMENT*>(
777  Fluid_mesh_pt->boundary_element_pt(b,e));
778 
779  //Find the index of the face of element e along boundary b
780  int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
781 
782  // Create new element
783  ImposeDisplacementByLagrangeMultiplierElement<FLUID_ELEMENT>* el_pt =
784  new ImposeDisplacementByLagrangeMultiplierElement<FLUID_ELEMENT>(
785  bulk_elem_pt,face_index);
786 
787  // Add it to the mesh
788  Lagrange_multiplier_mesh_pt[i]->add_element_pt(el_pt);
789 
790  // Set the GeomObject that defines the boundary shape and set
791  // which bulk boundary we are attached to (needed to extract
792  // the boundary coordinate from the bulk nodes)
793  el_pt->set_boundary_shape_geom_object_pt(Solid_fsi_boundary_pt[i],b);
794  }
795  }
796 
797 } // end of create_lagrange_multiplier_elements
798 
799 
800 
801 //============start_of_fluid_traction_elements==============================
802 /// Create fluid traction elements
803 //=======================================================================
804 template<class FLUID_ELEMENT,class SOLID_ELEMENT>
807 {
808 
809  // Counter for number of fluid traction meshes
810  unsigned count=0;
811 
812  // Loop over inflow/outflow boundaries
813  for (unsigned in_out=0;in_out<2;in_out++)
814  {
815  // Loop over boundaries with fluid traction elements
816  unsigned n=nfluid_inflow_traction_boundary();
817  if (in_out==1) n=nfluid_outflow_traction_boundary();
818  for (unsigned i=0;i<n;i++)
819  {
820 
821  // Get boundary ID
822  unsigned b=0;
823  if (in_out==0)
824  {
825  b=Inflow_boundary_id[i];
826  }
827  else
828  {
829  b=Outflow_boundary_id[i];
830  }
831 
832  // How many bulk elements are adjacent to boundary b?
833  unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
834 
835  // Loop over the bulk elements adjacent to boundary b
836  for(unsigned e=0;e<n_element;e++)
837  {
838  // Get pointer to the bulk element that is adjacent to boundary b
839  FLUID_ELEMENT* bulk_elem_pt = dynamic_cast<FLUID_ELEMENT*>(
840  Fluid_mesh_pt->boundary_element_pt(b,e));
841 
842  //What is the index of the face of the element e along boundary b
843  int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
844 
845  // Create new element
846  NavierStokesTractionElement<FLUID_ELEMENT>* el_pt=
847  new NavierStokesTractionElement<FLUID_ELEMENT>(bulk_elem_pt,
848  face_index);
849 
850  // Add it to the mesh
851  Fluid_traction_mesh_pt[count]->add_element_pt(el_pt);
852 
853  // Set the pointer to the prescribed traction function
854  if (in_out==0)
855  {
856  el_pt->traction_fct_pt() =
858  }
859  else
860  {
861  el_pt->traction_fct_pt() =
863  }
864  }
865  // Bump up counter
866  count++;
867  }
868  }
869 
870  } // end of create_traction_elements
871 
872 
873 //============start_doc_solid_zeta=======================================
874 /// Doc boundary coordinates of i-th solid FSI boundary.
875 //=======================================================================
876 template<class FLUID_ELEMENT,class SOLID_ELEMENT>
878 doc_solid_boundary_coordinates(const unsigned& i)
879 {
880 
881  //Doc boundary coordinates in fluid
882  char filename[100];
883  sprintf(filename,"RESLT/solid_boundary_coordinates%i.dat",i);
884  std::ofstream the_file(filename);
885 
886  // Loop over traction elements
887  unsigned n_face_element = Solid_fsi_traction_mesh_pt[i]->nelement();
888  for(unsigned e=0;e<n_face_element;e++)
889  {
890  //Cast the element pointer
891  FSISolidTractionElement<SOLID_ELEMENT,3>* el_pt=
892  dynamic_cast<FSISolidTractionElement<SOLID_ELEMENT,3>*>
893  (Solid_fsi_traction_mesh_pt[i]->element_pt(e));
894 
895  // Doc boundary coordinate
896  Vector<double> s(2);
897  Vector<double> zeta(2);
898  Vector<double> x(3);
899  unsigned n_plot=3;
900  the_file << el_pt->tecplot_zone_string(n_plot);
901 
902  // Loop over plot points
903  unsigned num_plot_points=el_pt->nplot_points(n_plot);
904  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
905  {
906  // Get local coordinates of plot point
907  el_pt->get_s_plot(iplot,n_plot,s);
908  el_pt->interpolated_zeta(s,zeta);
909  el_pt->interpolated_x(s,x);
910  for (unsigned i=0;i<3;i++)
911  {
912  the_file << x[i] << " ";
913  }
914  for (unsigned i=0;i<2;i++)
915  {
916  the_file << zeta[i] << " ";
917  }
918 
919  the_file << std::endl;
920  }
921  el_pt->write_tecplot_zone_footer(the_file,n_plot);
922  }
923 
924  // Close doc file
925  the_file.close();
926 
927 } // end doc_solid_zeta
928 
929 
930 //========start_of_doc_solution===========================================
931 /// Doc the solution
932 //========================================================================
933 template<class FLUID_ELEMENT, class SOLID_ELEMENT>
935 doc_solution(DocInfo& doc_info)
936 {
937 
938  ofstream some_file;
939  char filename[100];
940 
941  // Number of plot points
942  unsigned npts;
943  npts=5;
944 
945  // Output solid boundaries
946  //------------------------
947  sprintf(filename,"%s/solid_boundaries%i.dat",doc_info.directory().c_str(),
948  doc_info.number());
949  some_file.open(filename);
950  Solid_mesh_pt->output_boundaries(some_file);
951  some_file.close();
952 
953 
954  // Output solid solution
955  //-----------------------
956  sprintf(filename,"%s/solid_soln%i.dat",doc_info.directory().c_str(),
957  doc_info.number());
958  some_file.open(filename);
959  Solid_mesh_pt->output(some_file,npts);
960  some_file.close();
961 
962 
963  // Output fluid boundaries
964  //------------------------
965  sprintf(filename,"%s/fluid_boundaries%i.dat",doc_info.directory().c_str(),
966  doc_info.number());
967  some_file.open(filename);
968  Fluid_mesh_pt->output_boundaries(some_file);
969  some_file.close();
970 
971 
972  // Output fluid solution
973  //-----------------------
974  sprintf(filename,"%s/fluid_soln%i.dat",doc_info.directory().c_str(),
975  doc_info.number());
976  some_file.open(filename);
977  Fluid_mesh_pt->output(some_file,npts);
978  some_file.close();
979 
980 
981  // Output fsi traction
982  //--------------------
983  sprintf(filename,"%s/fsi_traction%i.dat",doc_info.directory().c_str(),
984  doc_info.number());
985  some_file.open(filename);
986  unsigned n=nsolid_fsi_boundary();
987  for (unsigned i=0;i<n;i++)
988  {
989  Solid_fsi_traction_mesh_pt[i]->output(some_file,npts);
990  }
991  some_file.close();
992 
993 } // end_of_doc
994 
995 
996 
997 
998 
999 //========================= start_of_main=================================
1000 /// Demonstrate how to solve an unstructured 3D FSI problem
1001 //========================================================================
1002 int main(int argc, char **argv)
1003 {
1004  // Label for output
1005  DocInfo doc_info;
1006 
1007  // Output directory
1008  doc_info.set_directory("RESLT");
1009 
1010  // Create generalised Hookean constitutive equations
1012  new GeneralisedHookean(&Global_Parameters::Nu);
1013 
1014  //Set up the problem
1016  PseudoSolidNodeUpdateElement<TTaylorHoodElement<3>, TPVDElement<3,3> >,
1017  TPVDElement<3,3> > problem;
1018 
1019  //Output initial configuration
1020  problem.doc_solution(doc_info);
1021  doc_info.number()++;
1022 
1023  // Parameter study
1024  unsigned nstep=2;
1025 
1026  // Increment in FSI parameter
1027  double q_increment=5.0e-2;
1028 
1029  for (unsigned istep=0;istep<nstep;istep++)
1030  {
1031  // Solve the problem
1032  problem.newton_solve();
1033 
1034  //Output solution
1035  problem.doc_solution(doc_info);
1036  doc_info.number()++;
1037 
1038  // Bump up FSI parameter
1039  Global_Parameters::Q+=q_increment;
1040  }
1041 
1042 } // end_of_main
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
FluidTetMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &face_file_name, const bool &split_corner_elements, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor:
virtual ~FluidTetMesh()
Empty Destructor.
Tetgen-based mesh upgraded to become a solid mesh.
MySolidTetgenMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &face_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor:
virtual ~MySolidTetgenMesh()
Empty Destructor.
Unstructured 3D FSI problem.
void doc_solution(DocInfo &doc_info)
Doc the solution.
Vector< Mesh * > Fluid_traction_mesh_pt
Meshes of fluid traction elements that apply pressure at in/outflow.
unsigned npinned_solid_boundary()
Return total number of mesh boundaries in the solid mesh where the position is pinned.
Vector< unsigned > Solid_fsi_boundary_id
IDs of solid mesh boundaries which make up the FSI interface.
Vector< MeshAsGeomObject * > Solid_fsi_boundary_pt
GeomObject incarnations of the FSI boundary in the solid mesh.
Vector< SolidMesh * > Lagrange_multiplier_mesh_pt
Meshes of Lagrange multiplier elements.
unsigned nsolid_fsi_boundary()
Return total number of mesh boundaries in the solid mesh that make up the FSI interface.
unsigned nfluid_traction_boundary()
Return total number of mesh boundaries that make up the in- and outflow boundaries where a traction h...
MySolidTetgenMesh< SOLID_ELEMENT > * Solid_mesh_pt
Bulk solid mesh.
Vector< SolidMesh * > Solid_fsi_traction_mesh_pt
Meshes of FSI traction elements.
void create_lagrange_multiplier_elements()
Create elements that enforce prescribed boundary motion for the pseudo-solid fluid mesh by Lagrange m...
unsigned nfluid_inflow_traction_boundary()
Return total number of mesh boundaries that make up the inflow boundary.
Vector< unsigned > Inflow_boundary_id
IDs of fluid mesh boundaries along which inflow boundary conditions are applied.
Vector< unsigned > Pinned_solid_boundary_id
IDs of solid mesh boundaries where displacements are pinned.
FluidTetMesh< FLUID_ELEMENT > * Fluid_mesh_pt
Bulk fluid mesh.
Vector< unsigned > Fluid_fsi_boundary_id
IDs of fluid mesh boundaries which make up the FSI interface.
unsigned nfluid_fsi_boundary()
Return total number of mesh boundaries in the fluid mesh that make up the FSI interface.
void doc_solid_boundary_coordinates(const unsigned &i)
Sanity check: Doc boundary coordinates on i-th solid FSI interface.
void create_fsi_traction_elements()
Create FSI traction elements.
~UnstructuredFSIProblem()
Destructor (empty)
Vector< unsigned > Outflow_boundary_id
IDs of fluid mesh boundaries along which inflow boundary conditions are applied.
unsigned nfluid_outflow_traction_boundary()
Return total number of mesh boundaries that make up the outflow boundary.
void create_fluid_traction_elements()
Create fluid traction elements at inflow.
/////////////////////////////////////////////////////////////// /////////////////////////////////////...
double P_in
Fluid pressure on inflow boundary.
double Nu
Poisson's ratio for generalised Hookean constitutive equation.
double Q
Default FSI parameter.
void prescribed_outflow_traction(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Applied traction on fluid at the inflow boundary.
double Re
Default Reynolds number.
double P_out
Fluid pressure on outflow boundary.
ConstitutiveLaw * Constitutive_law_pt
Pointer to constitutive law.
void prescribed_inflow_traction(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Applied traction on fluid at the inflow boundary.
int main(int argc, char **argv)
Demonstrate how to solve an unstructured 3D FSI problem.