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-2022 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
39using namespace std;
40using namespace oomph;
41
42
43
44//==========start_solid_mesh===============================================
45/// Tetgen-based mesh upgraded to become a solid mesh
46//=========================================================================
47template<class ELEMENT>
48class MySolidTetgenMesh : public virtual TetgenMesh<ELEMENT>,
49 public virtual SolidMesh
50{
51
52public:
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//=========================================================================
97template<class ELEMENT>
98class FluidTetMesh : public virtual TetgenMesh<ELEMENT>,
99 public virtual SolidMesh
100{
101
102public:
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//====================================================================
211template<class FLUID_ELEMENT, class SOLID_ELEMENT>
212class UnstructuredFSIProblem : public Problem
213{
214
215public:
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
237private:
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
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//========================================================================
317template<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//=======================================================================
703template<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//=======================================================================
750template<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//=======================================================================
804template<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//=======================================================================
876template<class FLUID_ELEMENT,class SOLID_ELEMENT>
878doc_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//========================================================================
933template<class FLUID_ELEMENT, class SOLID_ELEMENT>
935doc_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//========================================================================
1002int 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.