unstructured_two_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 for fsi on unstructured mesh
27
28//Generic includes
29#include "generic.h"
30#include "navier_stokes.h"
31#include "solid.h"
32#include "constitutive.h"
33
34// The meshes
35#include "meshes/triangle_mesh.h"
36
37using namespace std;
38using namespace oomph;
39
40//==start_of_namespace==============================
41/// Namespace for physical parameters
42//==================================================
44{
45 /// Reynolds number
46 double Re=0.0;
47
48 /// FSI parameter
49 double Q=0.0;
50
51 /// Non-dim gravity
52 double Gravity=0.0;
53
54 /// Non-dimensional gravity as body force
55 void gravity(const double& time,
56 const Vector<double> &xi,
57 Vector<double> &b)
58 {
59 b[0]=0.0;
60 b[1]=-Gravity;
61 }
62
63 /// Pseudo-solid Poisson ratio
64 double Nu=0.3;
65
66 /// Constitutive law for the solid (and pseudo-solid) mechanics
67 ConstitutiveLaw *Constitutive_law_pt=0;
68
69 /// Boolean to identify if node is on fsi boundary
70 bool is_on_fsi_boundary(Node* nod_pt)
71 {
72 if (
73 (
74 // Is it a boundary node?
75 dynamic_cast<BoundaryNodeBase*>(nod_pt)!=0)&&
76 (
77 // Horizontal extent of main immersed obstacle
78 ( (nod_pt->x(0)>1.6)&&(nod_pt->x(0)<4.75)&&
79 // Vertical extent of main immersed obstacle
80 (nod_pt->x(1)>0.1125)&&(nod_pt->x(1)<2.8) ) ||
81 // Two nodes on the bottom wall are below y=0.3
82 ( (nod_pt->x(1)<0.3)&&
83 // ...and bracketed in these two x-ranges
84 ( ( (nod_pt->x(0)>3.0)&&(nod_pt->x(0)<3.1) ) ||
85 ( (nod_pt->x(0)<4.6)&&(nod_pt->x(0)>4.5) ) )
86 )
87 )
88 )
89 {
90 return true;
91 }
92 else
93 {
94 return false;
95 }
96 }
97
98
99} // end_of_namespace
100
101
102
103/// /////////////////////////////////////////////////////////////////////
104/// /////////////////////////////////////////////////////////////////////
105/// /////////////////////////////////////////////////////////////////////
106
107
108//========start_solid_mesh=================================================
109/// Triangle-based mesh upgraded to become a solid mesh
110//=========================================================================
111template<class ELEMENT>
112class MySolidTriangleMesh : public virtual TriangleMesh<ELEMENT>,
113 public virtual SolidMesh
114{
115
116public:
117
118 /// Constructor:
119 MySolidTriangleMesh(const std::string& node_file_name,
120 const std::string& element_file_name,
121 const std::string& poly_file_name,
122 TimeStepper* time_stepper_pt=
123 &Mesh::Default_TimeStepper) :
124 TriangleMesh<ELEMENT>(node_file_name,element_file_name,
125 poly_file_name, time_stepper_pt)
126 {
127 //Assign the Lagrangian coordinates
128 set_lagrangian_nodal_coordinates();
129
130 // Identify special boundaries
131 set_nboundary(3);
132
133 unsigned n_node=this->nnode();
134 for (unsigned j=0;j<n_node;j++)
135 {
136 Node* nod_pt=this->node_pt(j);
137
138 // Boundary 1 is lower boundary
139 if (nod_pt->x(1)<0.15)
140 {
141 this->remove_boundary_node(0,nod_pt);
142 this->add_boundary_node(1,nod_pt);
143 }
144
145 // Boundary 2 is FSI interface
147 {
148 this->remove_boundary_node(0,nod_pt);
149 this->add_boundary_node(2,nod_pt);
150 }
151 }// done boundary assignment
152
153 // Identify the elements next to the newly created boundaries
154 TriangleMesh<ELEMENT>::setup_boundary_element_info();
155
156 // Setup boundary coordinates for boundaries 1 and 2
157 this->template setup_boundary_coordinates<ELEMENT>(1);
158 this->template setup_boundary_coordinates<ELEMENT>(2);
159 }
160
161 /// Empty Destructor
163
164};
165
166
167
168/// ////////////////////////////////////////////////////////////////////////
169/// ////////////////////////////////////////////////////////////////////////
170/// ////////////////////////////////////////////////////////////////////////
171
172
173
174//======================start_fluid_mesh===================================
175/// Triangle-based mesh upgraded to become a pseudo-solid mesh
176//=========================================================================
177template<class ELEMENT>
178class FluidTriangleMesh : public virtual TriangleMesh<ELEMENT>,
179 public virtual SolidMesh
180{
181
182public:
183
184 /// Constructor
185 FluidTriangleMesh(const std::string& node_file_name,
186 const std::string& element_file_name,
187 const std::string& poly_file_name,
188 TimeStepper* time_stepper_pt=
189 &Mesh::Default_TimeStepper) :
190 TriangleMesh<ELEMENT>(node_file_name,element_file_name,
191 poly_file_name, time_stepper_pt)
192 {
193 //Assign the Lagrangian coordinates
194 set_lagrangian_nodal_coordinates();
195
196 // Identify special boundaries
197 this->set_nboundary(4);
198
199 unsigned n_node=this->nnode();
200 for (unsigned j=0;j<n_node;j++)
201 {
202 Node* nod_pt=this->node_pt(j);
203
204 // Boundary 1 is left (inflow) boundary
205 if (nod_pt->x(0)<0.226)
206 {
207 this->remove_boundary_node(0,nod_pt);
208 this->add_boundary_node(1,nod_pt);
209
210 // Add overlapping nodes back to boundary 0
211 if (nod_pt->x(1)<0.2) this->add_boundary_node(0,nod_pt);
212 if (nod_pt->x(1)>4.06) this->add_boundary_node(0,nod_pt);
213 }
214
215 // Boundary 2 is right (outflow) boundary
216 if (nod_pt->x(0)>8.28)
217 {
218 this->remove_boundary_node(0,nod_pt);
219 this->add_boundary_node(2,nod_pt);
220
221 // Add overlapping nodes back to boundary 0
222 if (nod_pt->x(1)<0.2) this->add_boundary_node(0,nod_pt);
223 if (nod_pt->x(1)>4.06) this->add_boundary_node(0,nod_pt);
224
225 }
226
227 // Boundary 3 is FSI boundary
229 {
230 this->remove_boundary_node(0,nod_pt);
231 this->add_boundary_node(3,nod_pt);
232
233 //If it's below y=0.2 it's also on boundary 0 so stick it back on
234 if (nod_pt->x(1)<0.2) this->add_boundary_node(0,nod_pt);
235 }
236 }
237 TriangleMesh<ELEMENT>::setup_boundary_element_info();
238
239 // Open a file to doc the FaceElements that are used to
240 // create the boundary coordinates. The elements must
241 // form a simply-connected line. This may not work
242 // if the mesh is too coarse so that, e.g. an element
243 // that goes through the interior has endpoints that are
244 // both located on the same boundary. Outputting the
245 // FaceElements can help identify such cases.
246 ofstream some_file("RESLT/boundary_generation_test.dat");
247
248 // Setup boundary coordinates for boundaries 1, 2 and 3
249 this->template setup_boundary_coordinates<ELEMENT>(1);
250 this->template setup_boundary_coordinates<ELEMENT>(2);
251 this->template setup_boundary_coordinates<ELEMENT>(3,some_file);
252
253 // Close it again
254 some_file.close();
255 }
256
257 /// Empty Destructor
258 virtual ~FluidTriangleMesh() { }
259
260};
261
262
263
264/// ////////////////////////////////////////////////////////////////////////
265/// ////////////////////////////////////////////////////////////////////////
266/// ////////////////////////////////////////////////////////////////////////
267
268
269//==start_of_problem_class============================================
270/// Unstructured FSI Problem
271//====================================================================
272template<class FLUID_ELEMENT, class SOLID_ELEMENT>
273class UnstructuredFSIProblem : public Problem
274{
275
276public:
277
278 /// Constructor
280
281 /// Destructor (empty)
283
284 /// Access function for the fluid mesh
286 {
287 return Fluid_mesh_pt;
288 }
289
290 /// Access function for the solid mesh
292 {
293 return Solid_mesh_pt;
294 }
295
296 /// Doc the solution
297 void doc_solution(DocInfo& doc_info);
298
299private:
300
301 /// Create FSI traction elements
302 void create_fsi_traction_elements();
303
304 /// Create elements that enforce prescribed boundary motion
305 /// for the pseudo-solid fluid mesh by Lagrange multipliers
306 void create_lagrange_multiplier_elements();
307
308 /// Sanity check: Doc boundary coordinates from solid side
309 void doc_solid_boundary_coordinates();
310
311 /// Fluid mesh
313
314 /// Solid mesh
316
317 /// Pointers to mesh of Lagrange multiplier elements
319
320 /// Vector of pointers to mesh of FSI traction elements
322
323 /// GeomObject incarnation of fsi boundary in solid mesh
324 MeshAsGeomObject*
326
327};
328
329
330
331//==start_of_constructor==================================================
332/// Constructor for unstructured FSI problem.
333//========================================================================
334template<class FLUID_ELEMENT, class SOLID_ELEMENT>
336{
337
338 // Fluid mesh
339 //-----------
340
341 //Create fluid mesh
342 string fluid_node_file_name="fluid.fig.1.node";
343 string fluid_element_file_name="fluid.fig.1.ele";
344 string fluid_poly_file_name="fluid.fig.1.poly";
345 Fluid_mesh_pt = new FluidTriangleMesh<FLUID_ELEMENT>(fluid_node_file_name,
346 fluid_element_file_name,
347 fluid_poly_file_name);
348
349 // Doc pinned solid nodes
350 std::ofstream pseudo_solid_bc_file("pinned_pseudo_solid_nodes.dat");
351
352 // Set the boundary conditions for fluid problem: All nodes are
353 // free by default -- just pin the ones that have Dirichlet conditions
354 // here.
355 unsigned nbound=Fluid_mesh_pt->nboundary();
356 for(unsigned ibound=0;ibound<nbound;ibound++)
357 {
358 unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
359 for (unsigned inod=0;inod<num_nod;inod++)
360 {
361 // Pin velocity everywhere apart from outlet where we
362 // have parallel outflow
363 if (ibound!=2)
364 {
365 Fluid_mesh_pt->boundary_node_pt(ibound,inod)->pin(0);
366 }
367 Fluid_mesh_pt->boundary_node_pt(ibound,inod)->pin(1);
368
369 // Pin pseudo-solid positions everywhere apart from boundary 3,
370 // the fsi boundary
371 if ((ibound==0)||(ibound==1)||(ibound==2))
372 {
373 for(unsigned i=0;i<2;i++)
374 {
375 // Pin the node
376 SolidNode* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
377 nod_pt->pin_position(i);
378
379 // Doc it as pinned
380 pseudo_solid_bc_file << nod_pt->x(i) << " ";
381 }
382 pseudo_solid_bc_file << std::endl;
383 }
384 }
385 } // end loop over boundaries
386
387 // Close
388 pseudo_solid_bc_file.close();
389
390
391 // Complete the build of the fluid elements so they are fully functional
392 unsigned n_element = Fluid_mesh_pt->nelement();
393 for(unsigned e=0;e<n_element;e++)
394 {
395 // Upcast from GeneralisedElement to the present element
396 FLUID_ELEMENT* el_pt =
397 dynamic_cast<FLUID_ELEMENT*>(Fluid_mesh_pt->element_pt(e));
398
399 //Set the Reynolds number
400 el_pt->re_pt() = &Global_Parameters::Re;
401
402 // Set the constitutive law for pseudo-elastic mesh deformation
403 el_pt->constitutive_law_pt() =
405
406 } // end loop over elements
407
408
409 // Apply fluid boundary conditions: Poiseuille at inflow
410
411 // Find max. and min y-coordinate at inflow
412 unsigned ibound=1;
413 //Initialise both to the y-coordinate of the first boundary node
414 double y_min=fluid_mesh_pt()->boundary_node_pt(ibound,0)->x(1);;
415 double y_max=y_min;
416
417 //Loop over the rest of the boundary nodes
418 unsigned num_nod= fluid_mesh_pt()->nboundary_node(ibound);
419 for (unsigned inod=1;inod<num_nod;inod++)
420 {
421 double y=fluid_mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
422 if (y>y_max)
423 {
424 y_max=y;
425 }
426 if (y<y_min)
427 {
428 y_min=y;
429 }
430 }
431 double y_mid=0.5*(y_min+y_max);
432
433 // Loop over all boundaries
434 const unsigned n_boundary = fluid_mesh_pt()->nboundary();
435 for (unsigned ibound=0;ibound<n_boundary;ibound++)
436 {
437 const unsigned num_nod= fluid_mesh_pt()->nboundary_node(ibound);
438 for (unsigned inod=0;inod<num_nod;inod++)
439 {
440 // Parabolic inflow at the left boundary (boundary 1)
441 if (ibound==1)
442 {
443 double y=fluid_mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
444 double veloc=1.5/(y_max-y_min)*
445 (y-y_min)*(y_max-y)/((y_mid-y_min)*(y_max-y_mid));
446 fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,veloc);
447 fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
448 }
449 // Zero flow elsewhere
450 else
451 {
452 fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,0.0);
453 fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
454 }
455 }
456 } // end Poiseuille
457
458
459 // Solid mesh
460 //-----------
461
462 //Create solid mesh
463 string solid_node_file_name="solid.fig.1.node";
464 string solid_element_file_name="solid.fig.1.ele";
465 string solid_poly_file_name="solid.fig.1.poly";
466 Solid_mesh_pt = new MySolidTriangleMesh<SOLID_ELEMENT>(solid_node_file_name,
467 solid_element_file_name,
468 solid_poly_file_name);
469
470
471 // Complete the build of all solid elements so they are fully functional
472 n_element = Solid_mesh_pt->nelement();
473 for(unsigned i=0;i<n_element;i++)
474 {
475 //Cast to a solid element
476 SOLID_ELEMENT *el_pt =
477 dynamic_cast<SOLID_ELEMENT*>(Solid_mesh_pt->element_pt(i));
478
479 // Set the constitutive law
480 el_pt->constitutive_law_pt() =
482
483 //Set the body force
484 el_pt->body_force_fct_pt() = Global_Parameters::gravity;
485 }
486
487 // Pin both positions at lower boundary of solid mesh (boundary 1)
488 ibound=1;
489 num_nod=Solid_mesh_pt->nboundary_node(ibound);
490 for (unsigned inod=0;inod<num_nod;inod++)
491 {
492 Solid_mesh_pt->boundary_node_pt(ibound,inod)->pin_position(0);
493 Solid_mesh_pt->boundary_node_pt(ibound,inod)->pin_position(1);
494 }
495
496
497 // Create FSI Traction elements
498 //-----------------------------
499
500 // Now construct the (empty) traction element mesh
501 Traction_mesh_pt=new SolidMesh;
502
503 // Build the FSI traction elements and add them to the traction mesh
504 create_fsi_traction_elements();
505
506 // Create Lagrange multiplier mesh for boundary motion
507 //----------------------------------------------------
508
509 // Construct the mesh of elements that enforce prescribed boundary motion
510 // of pseudo-solid fluid mesh by Lagrange multipliers
511 Lagrange_multiplier_mesh_pt=new SolidMesh;
512 create_lagrange_multiplier_elements();
513
514
515 // Combine meshes
516 //---------------
517
518 // Add sub meshes
519 add_sub_mesh(Fluid_mesh_pt);
520 add_sub_mesh(Solid_mesh_pt);
521 add_sub_mesh(Traction_mesh_pt);
522 add_sub_mesh(Lagrange_multiplier_mesh_pt);
523
524 // Build global mesh
525 build_global_mesh();
526
527
528 // Setup FSI
529 //----------
530
531 // Document the boundary coordinate along the FSI interface
532 // of the fluid mesh during call to
533 // setup_fluid_load_info_for_solid_elements()
534 Multi_domain_functions::Doc_boundary_coordinate_file.open(
535 "fluid_boundary_test.dat");
536
537 // Work out which fluid dofs affect the residuals of the wall elements:
538 // We pass the boundary between the fluid and solid meshes and
539 // pointers to the meshes. The interaction boundary is boundary 3
540 // of the 2D fluid mesh.
541 FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
542 (this,3,Fluid_mesh_pt,Traction_mesh_pt);
543
544 // Close the doc file
545 Multi_domain_functions::Doc_boundary_coordinate_file.close();
546
547 // Sanity check: Doc boundary coordinates from solid side
548 doc_solid_boundary_coordinates();
549
550 // Setup equation numbering scheme
551 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
552
553} // end_of_constructor
554
555
556
557//============start_doc_solid_zeta=======================================
558/// Doc boundary coordinates in solid and plot GeomObject representation
559/// of FSI boundary.
560//=======================================================================
561template<class FLUID_ELEMENT,class SOLID_ELEMENT>
564{
565
566 // Doc boundary coordinates for fsi boundary in solid mesh
567 std::ofstream the_file("solid_boundary_test.dat");
568
569 // Initialise max/min boundary coordinate
570 double zeta_min= 1.0e40;
571 double zeta_max=-1.0e40;
572
573 // Loop over FSI traction elements
574 unsigned n_face_element = Traction_mesh_pt->nelement();
575 for(unsigned e=0;e<n_face_element;e++)
576 {
577 //Cast the element pointer
578 FSISolidTractionElement<SOLID_ELEMENT,2>* el_pt=
579 dynamic_cast<FSISolidTractionElement<SOLID_ELEMENT,2>*>
580 (Traction_mesh_pt->element_pt(e));
581
582 // Doc boundary coordinate
583 Vector<double> s(1);
584 Vector<double> zeta(1);
585 Vector<double> x(2);
586 unsigned n_plot=5;
587 the_file << el_pt->tecplot_zone_string(n_plot);
588
589 // Loop over plot points
590 unsigned num_plot_points=el_pt->nplot_points(n_plot);
591 for (unsigned iplot=0;iplot<num_plot_points;iplot++)
592 {
593 // Get local coordinates of plot point
594 el_pt->get_s_plot(iplot,n_plot,s);
595 el_pt->interpolated_zeta(s,zeta);
596 el_pt->interpolated_x(s,x);
597 for (unsigned i=0;i<2;i++)
598 {
599 the_file << x[i] << " ";
600 }
601 the_file << zeta[0] << " ";
602
603 // Update max/min boundary coordinate
604 if (zeta[0]<zeta_min) zeta_min=zeta[0];
605 if (zeta[0]>zeta_max) zeta_max=zeta[0];
606
607 the_file << std::endl;
608 }
609 }
610
611 // Close doc file
612 the_file.close();
613
614
615 // Doc compound GeomObject
616 the_file.open("fsi_geom_object.dat");
617 unsigned nplot=10000;
618 Vector<double> zeta(1);
619 Vector<double> r(2);
620 for (unsigned i=0;i<nplot;i++)
621 {
622 zeta[0]=zeta_min+(zeta_max-zeta_min)*double(i)/double(nplot-1);
623 Solid_fsi_boundary_pt->position(zeta,r);
624 the_file << r[0] << " " << r[1] << std::endl;
625 }
626 the_file.close();
627
628} //end doc_solid_zeta
629
630
631
632
633//============start_of_create_traction_elements==========================
634/// Create FSI traction elements
635//=======================================================================
636template<class FLUID_ELEMENT,class SOLID_ELEMENT>
639{
640 // Traction elements are located on boundary 2 of solid bulk mesh
641 unsigned b=2;
642
643 // How many bulk elements are adjacent to boundary b?
644 unsigned n_element = solid_mesh_pt()->nboundary_element(b);
645
646 // Loop over the bulk elements adjacent to boundary b
647 for(unsigned e=0;e<n_element;e++)
648 {
649 // Get pointer to the bulk element that is adjacent to boundary b
650 SOLID_ELEMENT* bulk_elem_pt = dynamic_cast<SOLID_ELEMENT*>(
651 solid_mesh_pt()->boundary_element_pt(b,e));
652
653 //What is the index of the face of the element e along boundary b
654 int face_index = solid_mesh_pt()->face_index_at_boundary(b,e);
655
656 // Create new element
657 FSISolidTractionElement<SOLID_ELEMENT,2>* el_pt=
658 new FSISolidTractionElement<SOLID_ELEMENT,2>(bulk_elem_pt,face_index);
659
660 // Add it to the mesh
661 Traction_mesh_pt->add_element_pt(el_pt);
662
663 // Specify boundary number
664 el_pt->set_boundary_number_in_bulk_mesh(b);
665
666 // Function that specifies the load ratios
667 el_pt->q_pt() = &Global_Parameters::Q;
668 }
669
670 } // end of create_traction_elements
671
672
673
674//============start_of_create_lagrange_multiplier_elements===============
675/// Create elements that impose the prescribed boundary displacement
676/// for the pseudo-solid fluid mesh
677//=======================================================================
678template<class FLUID_ELEMENT, class SOLID_ELEMENT>
681{
682
683 // Create GeomObject incarnation of fsi boundary in solid mesh
684 Solid_fsi_boundary_pt=
685 new MeshAsGeomObject
686 (Traction_mesh_pt);
687
688 // Lagrange multiplier elements are located on boundary 3 of the fluid mesh
689 unsigned b=3;
690
691 // How many bulk fluid elements are adjacent to boundary b?
692 unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
693
694 // Loop over the bulk fluid elements adjacent to boundary b?
695 for(unsigned e=0;e<n_element;e++)
696 {
697 // Get pointer to the bulk fluid element that is adjacent to boundary b
698 FLUID_ELEMENT* bulk_elem_pt = dynamic_cast<FLUID_ELEMENT*>(
699 Fluid_mesh_pt->boundary_element_pt(b,e));
700
701 //Find the index of the face of element e along boundary b
702 int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
703
704 // Create new element
705 ImposeDisplacementByLagrangeMultiplierElement<FLUID_ELEMENT>* el_pt =
706 new ImposeDisplacementByLagrangeMultiplierElement<FLUID_ELEMENT>(
707 bulk_elem_pt,face_index);
708
709 // Add it to the mesh
710 Lagrange_multiplier_mesh_pt->add_element_pt(el_pt);
711
712 // Set the GeomObject that defines the boundary shape and set
713 // which bulk boundary we are attached to (needed to extract
714 // the boundary coordinate from the bulk nodes)
715 el_pt->set_boundary_shape_geom_object_pt(Solid_fsi_boundary_pt,b);
716
717 // Loop over the nodes to apply boundary conditions
718 unsigned nnod=el_pt->nnode();
719 for (unsigned j=0;j<nnod;j++)
720 {
721 Node* nod_pt = el_pt->node_pt(j);
722
723 // Is the node also on boundary 0?
724 if (nod_pt->is_on_boundary(0))
725 {
726 // How many nodal values were used by the "bulk" element
727 // that originally created this node?
728 unsigned n_bulk_value=el_pt->nbulk_value(j);
729
730 // The remaining ones are Lagrange multipliers and we pin them.
731 unsigned nval=nod_pt->nvalue();
732 for (unsigned j=n_bulk_value;j<nval;j++)
733 {
734 nod_pt->pin(j);
735 }
736 }
737 }
738 }
739
740} // end of create_lagrange_multiplier_elements
741
742
743//==start_of_doc_solution=================================================
744/// Doc the solution
745//========================================================================
746template<class FLUID_ELEMENT, class SOLID_ELEMENT>
748doc_solution(DocInfo& doc_info)
749{
750 ofstream some_file;
751 char filename[100];
752
753 // Number of plot points
754 unsigned npts;
755 npts=5;
756
757 // Output fluid solution
758 sprintf(filename,"%s/fluid_soln%i.dat",doc_info.directory().c_str(),
759 doc_info.number());
760 some_file.open(filename);
761 Fluid_mesh_pt->output(some_file,npts);
762 some_file.close();
763
764
765 // Output solid solution
766 sprintf(filename,"%s/solid_soln%i.dat",doc_info.directory().c_str(),
767 doc_info.number());
768 some_file.open(filename);
769 Solid_mesh_pt->output(some_file,npts);
770 some_file.close();
771
772
773} // end_of_doc_solution
774
775
776
777
778/// /////////////////////////////////////////////////////////////////////
779/// /////////////////////////////////////////////////////////////////////
780/// /////////////////////////////////////////////////////////////////////
781
782
783//==start_of_main======================================================
784/// Driver for unstructured fsi problem
785//=====================================================================
786int main()
787{
788 // Label for output
789 DocInfo doc_info;
790
791 // Set output directory
792 doc_info.set_directory("RESLT");
793
794 //Create the constitutive law
795 Global_Parameters::Constitutive_law_pt = new GeneralisedHookean(
797
798 // Build the problem with triangular Taylor Hood for fluid and solid
800 PseudoSolidNodeUpdateElement<TTaylorHoodElement<2>, TPVDElement<2,3> >,
801 TPVDElement<2,3> > problem;
802
803 // Output boundaries
804 problem.fluid_mesh_pt()->output_boundaries("RESLT/fluid_boundaries.dat");
805 problem.solid_mesh_pt()->output_boundaries("RESLT/solid_boundaries.dat");
806
807 // Output the initial guess for the solution
808 problem.doc_solution(doc_info);
809 doc_info.number()++;
810
811 // Parameter study
813 double q_increment=1.0e-6;
814
815 // Solve the problem at zero Re and Q
816 problem.newton_solve();
817
818 // Output the solution
819 problem.doc_solution(doc_info);
820 doc_info.number()++;
821
822 // Bump up Re
824
825 // Now do proper parameter study with increase in Q
826 unsigned nstep=2; // 10;
827 for (unsigned i=0;i<nstep;i++)
828 {
829 // Solve the problem
830 problem.newton_solve();
831
832 // Output the solution
833 problem.doc_solution(doc_info);
834 doc_info.number()++;
835
836 // Bump up Q
837 Global_Parameters::Q+=q_increment;
838 }
839
840
841} // end_of_main
842
843
844
845
846
847
848
849
850
851
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
FluidTriangleMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor.
virtual ~FluidTriangleMesh()
Empty Destructor.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
MySolidTriangleMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor:
virtual ~MySolidTriangleMesh()
Empty Destructor.
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
void doc_solution(DocInfo &doc_info)
Doc the solution.
FluidTriangleMesh< FLUID_ELEMENT > * Fluid_mesh_pt
Fluid mesh.
MySolidTriangleMesh< SOLID_ELEMENT > *& solid_mesh_pt()
Access function for the solid mesh.
void doc_solid_boundary_coordinates()
Sanity check: Doc boundary coordinates from solid side.
void create_lagrange_multiplier_elements()
Create elements that enforce prescribed boundary motion for the pseudo-solid fluid mesh by Lagrange m...
FluidTriangleMesh< FLUID_ELEMENT > *& fluid_mesh_pt()
Access function for the fluid mesh.
void create_fsi_traction_elements()
Create FSI traction elements.
~UnstructuredFSIProblem()
Destructor (empty)
MeshAsGeomObject * Solid_fsi_boundary_pt
GeomObject incarnation of fsi boundary in solid mesh.
SolidMesh * Lagrange_multiplier_mesh_pt
Pointers to mesh of Lagrange multiplier elements.
SolidMesh * Traction_mesh_pt
Vector of pointers to mesh of FSI traction elements.
MySolidTriangleMesh< SOLID_ELEMENT > * Solid_mesh_pt
Solid mesh.
Namespace for physical parameters.
void gravity(const double &time, const Vector< double > &xi, Vector< double > &b)
Non-dimensional gravity as body force.
double Nu
Pseudo-solid Poisson ratio.
double Gravity
Non-dim gravity.
bool is_on_fsi_boundary(Node *nod_pt)
Boolean to identify if node is on fsi boundary.
double Q
FSI parameter.
double Re
Reynolds number.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law for the solid (and pseudo-solid) mechanics.
int main()
///////////////////////////////////////////////////////////////////// ///////////////////////////////...