unstructured_acoustic_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// Version 1.0; svn revision $LastChangedRevision: 1282 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Driver for Helmholtz/TimeHarmonicTimeHarmonicLinElast coupling
31 
32 //Oomph-lib includes
33 #include "generic.h"
34 #include "helmholtz.h"
35 #include "time_harmonic_linear_elasticity.h"
36 #include "multi_physics.h"
37 
38 //The meshes
39 #include "meshes/annular_mesh.h"
40 #include "meshes/triangle_mesh.h"
41 
42 using namespace std;
43 using namespace oomph;
44 
45 ///////////////////////////////////////////////////////////////////////
46 ///////////////////////////////////////////////////////////////////////
47 // Straight line as geometric object
48 ///////////////////////////////////////////////////////////////////////
49 ///////////////////////////////////////////////////////////////////////
50 
51 
52 
53 //=========================================================================
54 /// Straight 1D line in 2D space
55 //=========================================================================
56 class MyStraightLine : public GeomObject
57 {
58 
59 public:
60 
61  /// Constructor: Pass start and end points
62  MyStraightLine(const Vector<double>& r_start, const Vector<double>& r_end)
63  : GeomObject(1,2), R_start(r_start), R_end(r_end)
64  { }
65 
66 
67  /// Broken copy constructor
69  {
70  BrokenCopy::broken_copy("MyStraightLine");
71  }
72 
73  /// Broken assignment operator
74  void operator=(const MyStraightLine&)
75  {
76  BrokenCopy::broken_assign("MyStraightLine");
77  }
78 
79 
80  /// Destructor: Empty
82 
83  /// \short Position Vector at Lagrangian coordinate zeta
84  void position(const Vector<double>& zeta, Vector<double>& r) const
85  {
86  // Position Vector
87  r[0] = R_start[0]+(R_end[0]-R_start[0])*zeta[0];
88  r[1] = R_start[1]+(R_end[1]-R_start[1])*zeta[0];
89  }
90 
91 private:
92 
93  /// Start point of line
94  Vector<double> R_start;
95 
96  /// End point of line
97  Vector<double> R_end;
98 
99 };
100 
101 
102 
103 //////////////////////////////////////////////////////////////////////
104 //////////////////////////////////////////////////////////////////////
105 //////////////////////////////////////////////////////////////////////
106 
107 
108 //=======start_namespace==========================================
109 /// Global variables
110 //================================================================
111 namespace Global_Parameters
112 {
113 
114  /// \short Square of wavenumber for the Helmholtz equation
115  double K_squared=10.0;
116 
117  /// \short Radius of outer boundary of Helmholtz domain
118  double Outer_radius=4.0;
119 
120  /// Order of absorbing/appproximate boundary condition
121  unsigned ABC_order=3;
122 
123  /// FSI parameter
124  double Q=0.0;
125 
126  /// Non-dim thickness of elastic coating
127  double H_coating=0.1;
128 
129  /// Poisson's ratio
130  double Nu = 0.3;
131 
132  /// Square of non-dim frequency for the two regions -- dependent variable!
133  Vector<double> Omega_sq(2,0.0);
134 
135  /// Density ratio for the two regions: solid to fluid
136  Vector<double> Density_ratio(2,0.1);
137 
138  /// The elasticity tensors for the two regions
139  Vector<TimeHarmonicIsotropicElasticityTensor*> E_pt;
140 
141  /// Function to update dependent parameter values
143  {
144  Omega_sq[0]=Density_ratio[0]*Q;
145  Omega_sq[1]=Density_ratio[1]*Q;
146  }
147 
148  /// Uniform pressure
149  double P = 0.1;
150 
151  /// Peakiness parameter for pressure load
152  double Alpha=200.0;
153 
154  /// Pressure load (real and imag part)
155  void pressure_load(const Vector<double> &x,
156  const Vector<double> &n,
157  Vector<std::complex<double> >&traction)
158  {
159  double phi=atan2(x[1],x[0]);
160  double magnitude=exp(-Alpha*pow(phi-0.25*MathematicalConstants::Pi,2));
161 
162  unsigned dim = traction.size();
163  for(unsigned i=0;i<dim;i++)
164  {
165  traction[i] = complex<double>(-magnitude*P*n[i],magnitude*P*n[i]);
166  }
167  } // end_of_pressure_load
168 
169  // Output directory
170  string Directory="RESLT";
171 
172  // Multiplier for number of elements
173  unsigned El_multiplier=1;
174 
175 } //end namespace
176 
177 
178 
179 //=============begin_problem============================================
180 /// Coated disk FSI
181 //======================================================================
182 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
183 class CoatedDiskProblem : public Problem
184 {
185 
186 public:
187 
188  /// Constructor:
190 
191  /// Update function (empty)
193 
194  /// Update function (empty)
196 
197  /// Actions before adapt: Wipe the face meshes
198  void actions_before_adapt();
199 
200  /// Actions after adapt: Rebuild the face meshes
201  void actions_after_adapt();
202 
203  /// Doc the solution
204  void doc_solution();
205 
206 private:
207 
208  // Complete problem setup
209  void complete_problem_setup();
210 
211  /// \short Create solid traction elements
212  void create_solid_traction_elements();
213 
214  /// \short Create FSI traction elements
215  void create_fsi_traction_elements();
216 
217  /// \short Create Helmholtz FSI flux elements
218  void create_helmholtz_fsi_flux_elements();
219 
220  /// Delete (face) elements in specified mesh
221  void delete_face_elements(Mesh* const & boundary_mesh_pt);
222 
223  /// \short Create ABC face elements
224  void create_helmholtz_ABC_elements();
225 
226  /// Setup interaction
227  void setup_interaction();
228 
229  /// Pointer to refineable solid mesh
230  RefineableTriangleMesh<ELASTICITY_ELEMENT>* Solid_mesh_pt;
231 
232  /// Pointer to mesh of solid traction elements
234 
235  /// Pointer to mesh of FSI traction elements
236  Mesh* FSI_traction_mesh_pt;
237 
238  /// Pointer to Helmholtz mesh
239  TreeBasedRefineableMeshBase* Helmholtz_mesh_pt;
240 
241  /// Pointer to mesh of Helmholtz FSI flux elements
242  Mesh* Helmholtz_fsi_flux_mesh_pt;
243 
244  /// \short Pointer to mesh containing the ABC elements
246 
247  /// Boundary ID of upper symmetry boundary
249 
250  /// Boundary ID of lower symmetry boundary
252 
253  /// Boundary ID of upper inner boundary
255 
256  /// Boundary ID of lower inner boundary
258 
259  /// Boundary ID of outer boundary
261 
262  // Boundary ID of rib divider
264 
265  /// DocInfo object for output
266  DocInfo Doc_info;
267 
268  /// Trace file
269  ofstream Trace_file;
270 
271 };
272 
273 
274 //===========start_of_constructor=======================================
275 /// Constructor
276 //======================================================================
277 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
279 {
280  // To suppress boundary warnings (to avoid a lot of warnings) uncomment
281  // the following code:
282  //UnstructuredTwoDMeshGeometryBase::
283  // Suppress_warning_about_regions_and_boundaries=true;
284 
285  // Solid mesh
286  //-----------
287 
288  // Start and end coordinates
289  Vector<double> r_start(2);
290  Vector<double> r_end(2);
291 
292  // Outer radius of hull
293  double r_outer = 1.0;
294 
295  // Inner radius of hull
296  double r_inner = r_outer-Global_Parameters::H_coating;
297 
298  // Thickness of rib
299  double rib_thick=0.05;
300 
301  // Depth of rib
302  double rib_depth=0.2;
303 
304  // Total width of T
305  double t_width=0.2;
306 
307  // Thickness of T
308  double t_thick=0.05;
309 
310  // Half-opening angle of rib
311  double half_phi_rib=asin(0.5*rib_thick/r_inner);
312 
313  // Pointer to the closed curve that defines the outer boundary
314  TriangleMeshClosedCurve* closed_curve_pt=0;
315 
316  // Provide storage for pointers to the parts of the curvilinear boundary
317  Vector<TriangleMeshCurveSection*> curvilinear_boundary_pt;
318 
319  // Outer boundary
320  //---------------
321  Ellipse* outer_boundary_circle_pt = new Ellipse(r_outer,r_outer);
322  double zeta_start=-0.5*MathematicalConstants::Pi;
323  double zeta_end=0.5*MathematicalConstants::Pi;
324  unsigned nsegment=50;
325  unsigned boundary_id=curvilinear_boundary_pt.size();
326  curvilinear_boundary_pt.push_back(
327  new TriangleMeshCurviLine(
328  outer_boundary_circle_pt,zeta_start,zeta_end,nsegment,boundary_id));
329 
330  // Remember it
331  Outer_boundary_id=boundary_id;
332 
333 
334  // Upper straight line segment on symmetry axis
335  //---------------------------------------------
336  r_start[0]=0.0;
337  r_start[1]=r_outer;
338  r_end[0]=0.0;
339  r_end[1]=r_inner;
340  MyStraightLine* upper_sym_pt = new MyStraightLine(r_start,r_end);
341  zeta_start=0.0;
342  zeta_end=1.0;
343  nsegment=1;
344  boundary_id=curvilinear_boundary_pt.size();
345  curvilinear_boundary_pt.push_back(
346  new TriangleMeshCurviLine(
347  upper_sym_pt,zeta_start,zeta_end,nsegment,boundary_id));
348 
349  // Remember it
350  Upper_symmetry_boundary_id=boundary_id;
351 
352  // Upper part of inner boundary
353  //-----------------------------
354  Ellipse* upper_inner_boundary_pt =
355  new Ellipse(r_inner,r_inner);
356  zeta_start=0.5*MathematicalConstants::Pi;
357  zeta_end=half_phi_rib;
358  nsegment=20;
359  boundary_id=curvilinear_boundary_pt.size();
360  curvilinear_boundary_pt.push_back(
361  new TriangleMeshCurviLine(
362  upper_inner_boundary_pt,
363  zeta_start,zeta_end,nsegment,boundary_id));
364 
365  // Remember it
366  Upper_inner_boundary_id=boundary_id;
367 
368  // Data associated with rib
369  MyStraightLine* upper_inward_rib_pt=0;
370  MyStraightLine* lower_inward_rib_pt=0;
371  TriangleMeshCurviLine* upper_inward_rib_curviline_pt=0;
372  Vector<TriangleMeshOpenCurve*> inner_boundary_pt;
373  TriangleMeshCurviLine* lower_inward_rib_curviline_pt=0;
374  Vector<double> rib_center(2);
375 
376  // Upper half of inward rib
377  //-------------------------
378  r_start[0]=r_inner*cos(half_phi_rib);
379  r_start[1]=r_inner*sin(half_phi_rib);
380  r_end[0]=r_start[0]-rib_depth;
381  r_end[1]=r_start[1];
382  upper_inward_rib_pt = new MyStraightLine(r_start,r_end);
383  zeta_start=0.0;
384  zeta_end=1.0;
385  nsegment=1;
386  boundary_id=curvilinear_boundary_pt.size();
387  upper_inward_rib_curviline_pt=
388  new TriangleMeshCurviLine(
389  upper_inward_rib_pt,zeta_start,zeta_end,nsegment,boundary_id);
390  curvilinear_boundary_pt.push_back(upper_inward_rib_curviline_pt);
391 
392  // Vertical upper bit of T
393  //------------------------
394  r_start[0]=r_end[0];
395  r_start[1]=r_end[1];
396  r_end[0]=r_start[0];
397  r_end[1]=r_start[1]+0.5*(t_width-rib_thick);
398  MyStraightLine* vertical_upper_t_rib_pt = new MyStraightLine(r_start,r_end);
399  zeta_start=0.0;
400  zeta_end=1.0;
401  nsegment=1;
402  boundary_id=curvilinear_boundary_pt.size();
403  curvilinear_boundary_pt.push_back(
404  new TriangleMeshCurviLine(
405  vertical_upper_t_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
406 
407 
408  // Horizontal upper bit of T
409  //-----------==-------------
410  r_start[0]=r_end[0];
411  r_start[1]=r_end[1];
412  r_end[0]=r_start[0]-t_thick;
413  r_end[1]=r_start[1];
414  MyStraightLine* horizontal_upper_t_rib_pt = new MyStraightLine(r_start,r_end);
415  zeta_start=0.0;
416  zeta_end=1.0;
417  nsegment=1;
418  boundary_id=curvilinear_boundary_pt.size();
419  curvilinear_boundary_pt.push_back(
420  new TriangleMeshCurviLine(
421  horizontal_upper_t_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
422 
423  // Vertical end of rib end
424  //------------------------
425  r_start[0]=r_end[0];
426  r_start[1]=r_end[1];
427  r_end[0]=r_start[0];
428  r_end[1]=-r_start[1];
429  MyStraightLine* inner_vertical_rib_pt = new MyStraightLine(r_start,r_end);
430  zeta_start=0.0;
431  zeta_end=1.0;
432  nsegment=1;
433  boundary_id=curvilinear_boundary_pt.size();
434  curvilinear_boundary_pt.push_back(
435  new TriangleMeshCurviLine(
436  inner_vertical_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
437 
438 
439  // Horizontal lower bit of T
440  //-----------==-------------
441  r_start[0]=r_end[0];
442  r_start[1]=r_end[1];
443  r_end[0]=r_start[0]+t_thick;
444  r_end[1]=r_start[1];
445  MyStraightLine* horizontal_lower_t_rib_pt = new MyStraightLine(r_start,r_end);
446  zeta_start=0.0;
447  zeta_end=1.0;
448  nsegment=1;
449  boundary_id=curvilinear_boundary_pt.size();
450  curvilinear_boundary_pt.push_back(
451  new TriangleMeshCurviLine(
452  horizontal_lower_t_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
453 
454 
455  // Vertical lower bit of T
456  //------------------------
457  r_start[0]=r_end[0];
458  r_start[1]=r_end[1];
459  r_end[0]=r_start[0];
460  r_end[1]=r_start[1]+0.5*(t_width-rib_thick);
461  MyStraightLine* vertical_lower_t_rib_pt = new MyStraightLine(r_start,r_end);
462  zeta_start=0.0;
463  zeta_end=1.0;
464  nsegment=1;
465  boundary_id=curvilinear_boundary_pt.size();
466  curvilinear_boundary_pt.push_back(
467  new TriangleMeshCurviLine(
468  vertical_lower_t_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
469 
470 
471  // Lower half of inward rib
472  //-------------------------
473  r_end[0]=r_inner*cos(half_phi_rib);
474  r_end[1]=-r_inner*sin(half_phi_rib);
475  r_start[0]=r_end[0]-rib_depth;
476  r_start[1]=r_end[1];
477  lower_inward_rib_pt = new MyStraightLine(r_start,r_end);
478  zeta_start=0.0;
479  zeta_end=1.0;
480  nsegment=1;
481  boundary_id=curvilinear_boundary_pt.size();
482  lower_inward_rib_curviline_pt=
483  new TriangleMeshCurviLine(
484  lower_inward_rib_pt,zeta_start,zeta_end,nsegment,boundary_id);
485  curvilinear_boundary_pt.push_back(lower_inward_rib_curviline_pt);
486 
487 
488  // Lower part of inner boundary
489  //-----------------------------
490  Ellipse* lower_inner_boundary_circle_pt = new Ellipse(r_inner,r_inner);
491  zeta_start=-half_phi_rib;
492  zeta_end=-0.5*MathematicalConstants::Pi;
493  nsegment=20;
494  boundary_id=curvilinear_boundary_pt.size();
495  curvilinear_boundary_pt.push_back(
496  new TriangleMeshCurviLine(
497  lower_inner_boundary_circle_pt,zeta_start,zeta_end,nsegment,boundary_id));
498 
499  // Remember it
500  Lower_inner_boundary_id=boundary_id;
501 
502  // Lower straight line segment on symmetry axis
503  //---------------------------------------------
504  r_start[0]=0.0;
505  r_start[1]=-r_inner;
506  r_end[0]=0.0;
507  r_end[1]=-r_outer;
508  MyStraightLine* lower_sym_pt = new MyStraightLine(r_start,r_end);
509  zeta_start=0.0;
510  zeta_end=1.0;
511  nsegment=1;
512  boundary_id=curvilinear_boundary_pt.size();
513  curvilinear_boundary_pt.push_back(
514  new TriangleMeshCurviLine(
515  lower_sym_pt,zeta_start,zeta_end,nsegment,boundary_id));
516 
517  // Remember it
518  Lower_symmetry_boundary_id=boundary_id;
519 
520  // Combine to curvilinear boundary
521  //--------------------------------
522  closed_curve_pt=
523  new TriangleMeshClosedCurve(curvilinear_boundary_pt);
524 
525  // Vertical dividing line across base of T-rib
526  //--------------------------------------------
527  Vector<TriangleMeshCurveSection*> internal_polyline_pt(1);
528  r_start[0]=r_inner*cos(half_phi_rib);
529  r_start[1]=r_inner*sin(half_phi_rib);
530  r_end[0]=r_inner*cos(half_phi_rib);
531  r_end[1]=-r_inner*sin(half_phi_rib);
532 
533  Vector<Vector<double> > boundary_vertices(2);
534  boundary_vertices[0]=r_start;
535  boundary_vertices[1]=r_end;
536  boundary_id=100;
537  TriangleMeshPolyLine* rib_divider_pt=
538  new TriangleMeshPolyLine(boundary_vertices,boundary_id);
539  internal_polyline_pt[0]=rib_divider_pt;
540 
541  // Remember it
542  Rib_divider_boundary_id=boundary_id;
543 
544  // Make connection
545  double s_connect=0.0;
546  internal_polyline_pt[0]->connect_initial_vertex_to_curviline(
547  upper_inward_rib_curviline_pt,s_connect);
548 
549  // Make connection
550  s_connect=1.0;
551  internal_polyline_pt[0]->connect_final_vertex_to_curviline(
552  lower_inward_rib_curviline_pt,s_connect);
553 
554  // Create open curve that defines internal bondary
555  inner_boundary_pt.push_back(new TriangleMeshOpenCurve(internal_polyline_pt));
556 
557  // Define coordinates of a point inside the rib
558  rib_center[0]=r_inner-rib_depth;
559  rib_center[1]=0.0;
560 
561 
562  // Now build the mesh
563  //===================
564 
565  // Use the TriangleMeshParameters object for helping on the manage of the
566  // TriangleMesh parameters. The only parameter that needs to take is the
567  // outer boundary.
568  TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
569 
570  // Target area
571  triangle_mesh_parameters.element_area()=0.2;
572 
573  // Specify the internal open boundary
574  triangle_mesh_parameters.internal_open_curves_pt()=inner_boundary_pt;
575 
576  // Define the region
577  triangle_mesh_parameters.add_region_coordinates(1,rib_center);
578 
579  // Build the mesh
580  Solid_mesh_pt=new
581  RefineableTriangleMesh<ELASTICITY_ELEMENT>(triangle_mesh_parameters);
582 
583  // Helmholtz mesh
584  //---------------
585 
586  // Number of elements in azimuthal direction in Helmholtz mesh
587  unsigned ntheta_helmholtz=11*Global_Parameters::El_multiplier;
588 
589  // Number of elements in radial direction in Helmholtz mesh
590  unsigned nr_helmholtz=3*Global_Parameters::El_multiplier;
591 
592  // Innermost radius of Helmholtz mesh
593  double a=1.0;
594 
595  // Thickness of Helmholtz mesh
596  double h_thick_helmholtz=Global_Parameters::Outer_radius-a;
597 
598  // Build mesh
599  bool periodic=false;
600  double azimuthal_fraction=0.5;
601  double phi=MathematicalConstants::Pi/2.0;
602  Helmholtz_mesh_pt = new
603  RefineableTwoDAnnularMesh<HELMHOLTZ_ELEMENT>
604  (periodic,azimuthal_fraction,
605  ntheta_helmholtz,nr_helmholtz,a,h_thick_helmholtz,phi);
606 
607  // Set error estimators
608  Solid_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
609  Helmholtz_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
610 
611  // Mesh containing the Helmholtz ABC
612  // elements.
613  Helmholtz_outer_boundary_mesh_pt = new Mesh;
614 
615  // Create elasticity tensors
616  Global_Parameters::E_pt.resize(2);
617  Global_Parameters::E_pt[0]=new TimeHarmonicIsotropicElasticityTensor(
619  Global_Parameters::E_pt[1]=new TimeHarmonicIsotropicElasticityTensor(
621 
622  // Complete build of solid elements
623  complete_problem_setup();
624 
625  // Same for Helmholtz mesh
626  unsigned n_element =Helmholtz_mesh_pt->nelement();
627  for(unsigned i=0;i<n_element;i++)
628  {
629  //Cast to a solid element
630  HELMHOLTZ_ELEMENT *el_pt =
631  dynamic_cast<HELMHOLTZ_ELEMENT*>(Helmholtz_mesh_pt->element_pt(i));
632 
633  //Set the pointer to square of Helmholtz wavenumber
634  el_pt->k_squared_pt() = &Global_Parameters::K_squared;
635  }
636 
637  // Output meshes and their boundaries so far so we can double
638  // check the boundary enumeration
639  Solid_mesh_pt->output("solid_mesh.dat");
640  Helmholtz_mesh_pt->output("helmholtz_mesh.dat");
641  Solid_mesh_pt->output_boundaries("solid_mesh_boundary.dat");
642  Helmholtz_mesh_pt->output_boundaries("helmholtz_mesh_boundary.dat");
643 
644  // Create FaceElement meshes for boundary conditions
645  //---------------------------------------------------
646 
647  // Construct the solid traction element mesh
648  Solid_traction_mesh_pt=new Mesh;
649  create_solid_traction_elements();
650 
651  // Construct the fsi traction element mesh
652  FSI_traction_mesh_pt=new Mesh;
653  create_fsi_traction_elements();
654 
655  // Construct the Helmholtz fsi flux element mesh
656  Helmholtz_fsi_flux_mesh_pt=new Mesh;
657  create_helmholtz_fsi_flux_elements();
658 
659  // Create ABC elements on outer boundary of Helmholtz mesh
660  create_helmholtz_ABC_elements();
661 
662 
663  // Combine sub meshes
664  //-------------------
665 
666  // Solid mesh is first sub-mesh
667  add_sub_mesh(Solid_mesh_pt);
668 
669  // Add solid traction sub-mesh
670  add_sub_mesh(Solid_traction_mesh_pt);
671 
672  // Add FSI traction sub-mesh
673  add_sub_mesh(FSI_traction_mesh_pt);
674 
675  // Add Helmholtz mesh
676  add_sub_mesh(Helmholtz_mesh_pt);
677 
678  // Add Helmholtz FSI flux mesh
679  add_sub_mesh(Helmholtz_fsi_flux_mesh_pt);
680 
681  // Add Helmholtz ABC mesh
682  add_sub_mesh(Helmholtz_outer_boundary_mesh_pt);
683 
684  // Build combined "global" mesh
685  build_global_mesh();
686 
687 
688  // Setup fluid-structure interaction
689  //----------------------------------
690  setup_interaction();
691 
692  // Assign equation numbers
693  oomph_info << "Number of unknowns: " << assign_eqn_numbers() << std::endl;
694 
695  // Set output directory
696  Doc_info.set_directory(Global_Parameters::Directory);
697 
698  // Open trace file
699  char filename[100];
700  sprintf(filename,"%s/trace.dat",Doc_info.directory().c_str());
701  Trace_file.open(filename);
702 
703 
704 } //end of constructor
705 
706 
707 //=====================start_of_actions_before_adapt======================
708 /// Actions before adapt: Wipe the meshes face elements
709 //========================================================================
710 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
713 {
714  // Kill the solid traction elements and wipe surface mesh
715  delete_face_elements(Solid_traction_mesh_pt);
716 
717  // Kill the fsi traction elements and wipe surface mesh
718  delete_face_elements(FSI_traction_mesh_pt);
719 
720  // Kill Helmholtz FSI flux elements
721  delete_face_elements(Helmholtz_fsi_flux_mesh_pt);
722 
723  // Kill Helmholtz BC elements
724  delete_face_elements(Helmholtz_outer_boundary_mesh_pt);
725 
726  // Rebuild the Problem's global mesh from its various sub-meshes
727  rebuild_global_mesh();
728 
729 }// end of actions_before_adapt
730 
731 
732 
733 //=====================start_of_actions_after_adapt=======================
734 /// Actions after adapt: Rebuild the meshes of face elements
735 //========================================================================
736 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
739 {
740  // Complete problem setup
741  complete_problem_setup();
742 
743  // Construct the solid traction elements
744  create_solid_traction_elements();
745 
746  // Create fsi traction elements from all elements that are
747  // adjacent to FSI boundaries and add them to surface meshes
748  create_fsi_traction_elements();
749 
750  // Create Helmholtz fsi flux elements
751  create_helmholtz_fsi_flux_elements();
752 
753  // Create ABC elements from all elements that are
754  // adjacent to the outer boundary of Helmholtz mesh
755  create_helmholtz_ABC_elements();
756 
757  // Setup interaction
758  setup_interaction();
759 
760  // Rebuild the Problem's global mesh from its various sub-meshes
761  rebuild_global_mesh();
762 
763 }// end of actions_after_adapt
764 
765 
766 
767 
768 //=====================start_of_complete_problem_setup=======================
769 /// Complete problem setup: Apply boundary conditions and set
770 /// physical properties
771 //========================================================================
772 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
775 {
776 
777  // Solid boundary conditions:
778  //---------------------------
779  // Pin real and imag part of horizontal displacement components
780  //-------------------------------------------------------------
781  // on vertical boundaries
782  //-----------------------
783  {
784  //Loop over the nodes to pin and assign boundary displacements on
785  //solid boundary
786  unsigned n_node = Solid_mesh_pt->nboundary_node(Upper_symmetry_boundary_id);
787  for(unsigned i=0;i<n_node;i++)
788  {
789  Node* nod_pt=Solid_mesh_pt->boundary_node_pt(Upper_symmetry_boundary_id,i);
790 
791  // Real part of x-displacement
792  nod_pt->pin(0);
793  nod_pt->set_value(0,0.0);
794 
795  // Imag part of x-displacement
796  nod_pt->pin(2);
797  nod_pt->set_value(2,0.0);
798  }
799  }
800  {
801  //Loop over the nodes to pin and assign boundary displacements on
802  //solid boundary
803  unsigned n_node = Solid_mesh_pt->nboundary_node(Lower_symmetry_boundary_id);
804  for(unsigned i=0;i<n_node;i++)
805  {
806  Node* nod_pt=Solid_mesh_pt->boundary_node_pt(Lower_symmetry_boundary_id,i);
807 
808  // Real part of x-displacement
809  nod_pt->pin(0);
810  nod_pt->set_value(0,0.0);
811 
812  // Imag part of x-displacement
813  nod_pt->pin(2);
814  nod_pt->set_value(2,0.0);
815  }
816  }
817 
818 
819 
820  //Assign the physical properties to the elements
821  //----------------------------------------------
822  unsigned nreg=Solid_mesh_pt->nregion();
823  for (unsigned r=0;r<nreg;r++)
824  {
825  unsigned nel=Solid_mesh_pt->nregion_element(r);
826  for (unsigned e=0;e<nel;e++)
827  {
828  //Cast to a solid element
829  ELASTICITY_ELEMENT *el_pt =
830  dynamic_cast<ELASTICITY_ELEMENT*>(Solid_mesh_pt->
831  region_element_pt(r,e));
832 
833  // Set the constitutive law
834  el_pt->elasticity_tensor_pt() = Global_Parameters::E_pt[r];
835 
836  // Square of non-dim frequency
837  el_pt->omega_sq_pt()= &Global_Parameters::Omega_sq[r];
838  }
839  }
840 }
841 
842 //============start_of_delete_face_elements================
843 /// Delete face elements and wipe the mesh
844 //==========================================================
845 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
847 delete_face_elements(Mesh* const & boundary_mesh_pt)
848 {
849  // How many surface elements are in the surface mesh
850  unsigned n_element = boundary_mesh_pt->nelement();
851 
852  // Loop over the surface elements
853  for(unsigned e=0;e<n_element;e++)
854  {
855  // Kill surface element
856  delete boundary_mesh_pt->element_pt(e);
857  }
858 
859  // Wipe the mesh
860  boundary_mesh_pt->flush_element_and_node_storage();
861 
862 } // end of delete_face_elements
863 
864 
865 
866 //============start_of_create_solid_traction_elements====================
867 /// Create solid traction elements
868 //=======================================================================
869 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
872 {
873  // Loop over pressure loaded boundaries
874  unsigned b=0;
875  unsigned nb=3;
876  for (unsigned i=0;i<nb;i++)
877  {
878  switch(i)
879  {
880  case 0:
881  b=Upper_inner_boundary_id;
882  break;
883 
884  case 1:
885  b=Lower_inner_boundary_id;
886  break;
887 
888  case 2:
889  b=Rib_divider_boundary_id;
890  break;
891  }
892 
893  // We're attaching face elements to region 0
894  unsigned r=0;
895 
896  // How many bulk elements are adjacent to boundary b?
897  unsigned n_element = Solid_mesh_pt->nboundary_element_in_region(b,r);
898 
899  // Loop over the bulk elements adjacent to boundary b
900  for(unsigned e=0;e<n_element;e++)
901  {
902  // Get pointer to the bulk element that is adjacent to boundary b
903  ELASTICITY_ELEMENT* bulk_elem_pt = dynamic_cast<ELASTICITY_ELEMENT*>(
904  Solid_mesh_pt->boundary_element_in_region_pt(b,r,e));
905 
906  //Find the index of the face of element e along boundary b
907  int face_index = Solid_mesh_pt->face_index_at_boundary_in_region(b,r,e);
908 
909  // Create element
910  TimeHarmonicLinearElasticityTractionElement<ELASTICITY_ELEMENT>* el_pt=
911  new TimeHarmonicLinearElasticityTractionElement<ELASTICITY_ELEMENT>
912  (bulk_elem_pt,face_index);
913 
914  // Add to mesh
915  Solid_traction_mesh_pt->add_element_pt(el_pt);
916 
917  // Associate element with bulk boundary (to allow it to access
918  // the boundary coordinates in the bulk mesh)
919  el_pt->set_boundary_number_in_bulk_mesh(b);
920 
921  //Set the traction function
922  el_pt->traction_fct_pt() = Global_Parameters::pressure_load;
923  }
924  }
925 } // end of create_traction_elements
926 
927 
928 
929 
930 //============start_of_create_fsi_traction_elements======================
931 /// Create fsi traction elements
932 //=======================================================================
933 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
936 {
937  // We're on the outer boundary of the solid mesh
938  unsigned b=Outer_boundary_id;
939 
940  // How many bulk elements are adjacent to boundary b?
941  unsigned n_element = Solid_mesh_pt->nboundary_element(b);
942 
943  // Loop over the bulk elements adjacent to boundary b
944  for(unsigned e=0;e<n_element;e++)
945  {
946  // Get pointer to the bulk element that is adjacent to boundary b
947  ELASTICITY_ELEMENT* bulk_elem_pt = dynamic_cast<ELASTICITY_ELEMENT*>(
948  Solid_mesh_pt->boundary_element_pt(b,e));
949 
950  //Find the index of the face of element e along boundary b
951  int face_index = Solid_mesh_pt->face_index_at_boundary(b,e);
952 
953  // Create element
954  TimeHarmonicLinElastLoadedByHelmholtzPressureBCElement
955  <ELASTICITY_ELEMENT,HELMHOLTZ_ELEMENT>* el_pt=
956  new TimeHarmonicLinElastLoadedByHelmholtzPressureBCElement
957  <ELASTICITY_ELEMENT,HELMHOLTZ_ELEMENT>(bulk_elem_pt,
958  face_index);
959  // Add to mesh
960  FSI_traction_mesh_pt->add_element_pt(el_pt);
961 
962  // Associate element with bulk boundary (to allow it to access
963  // the boundary coordinates in the bulk mesh)
964  el_pt->set_boundary_number_in_bulk_mesh(b);
965 
966  // Set FSI parameter
967  el_pt->q_pt()=&Global_Parameters::Q;
968  }
969 
970 } // end of create_traction_elements
971 
972 
973 
974 
975 
976 //============start_of_create_helmholtz_fsi_flux_elements================
977 /// Create Helmholtz fsi flux elements
978 //=======================================================================
979 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
982 {
983 
984  // Attach to inner boundary of Helmholtz mesh (0)
985  unsigned b=0;
986 
987  // How many bulk elements are adjacent to boundary b?
988  unsigned n_element = Helmholtz_mesh_pt->nboundary_element(b);
989 
990  // Loop over the bulk elements adjacent to boundary b
991  for(unsigned e=0;e<n_element;e++)
992  {
993  // Get pointer to the bulk element that is adjacent to boundary b
994  HELMHOLTZ_ELEMENT* bulk_elem_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
995  Helmholtz_mesh_pt->boundary_element_pt(b,e));
996 
997  //Find the index of the face of element e along boundary b
998  int face_index = Helmholtz_mesh_pt->face_index_at_boundary(b,e);
999 
1000  // Create element
1001  HelmholtzFluxFromNormalDisplacementBCElement
1002  <HELMHOLTZ_ELEMENT,ELASTICITY_ELEMENT>* el_pt=
1003  new HelmholtzFluxFromNormalDisplacementBCElement
1004  <HELMHOLTZ_ELEMENT,ELASTICITY_ELEMENT>(bulk_elem_pt,
1005  face_index);
1006 
1007  // Add to mesh
1008  Helmholtz_fsi_flux_mesh_pt->add_element_pt(el_pt);
1009 
1010  // Associate element with bulk boundary (to allow it to access
1011  // the boundary coordinates in the bulk mesh)
1012  el_pt->set_boundary_number_in_bulk_mesh(b);
1013  }
1014 
1015 } // end of create_helmholtz_flux_elements
1016 
1017 
1018 
1019 //============start_of_create_ABC_elements==============================
1020 /// Create ABC elements on the outer boundary of
1021 /// the Helmholtz mesh
1022 //===========================================================================
1023 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
1026 {
1027  // We're on boundary 2 of the Helmholtz mesh
1028  unsigned b=2;
1029 
1030  // How many bulk elements are adjacent to boundary b?
1031  unsigned n_element = Helmholtz_mesh_pt->nboundary_element(b);
1032 
1033  // Loop over the bulk elements adjacent to boundary b?
1034  for(unsigned e=0;e<n_element;e++)
1035  {
1036  // Get pointer to the bulk element that is adjacent to boundary b
1037  HELMHOLTZ_ELEMENT* bulk_elem_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
1038  Helmholtz_mesh_pt->boundary_element_pt(b,e));
1039 
1040  //Find the index of the face of element e along boundary b
1041  int face_index = Helmholtz_mesh_pt->face_index_at_boundary(b,e);
1042 
1043  // Build the corresponding ABC element
1044  HelmholtzAbsorbingBCElement<HELMHOLTZ_ELEMENT>* flux_element_pt = new
1045  HelmholtzAbsorbingBCElement<HELMHOLTZ_ELEMENT>(bulk_elem_pt,face_index);
1046 
1047  // Set pointer to outer radius of artificial boundary
1048  flux_element_pt->outer_radius_pt()=&Global_Parameters::Outer_radius;
1049 
1050  // Set order of absorbing boundary condition
1051  flux_element_pt->abc_order_pt()=&Global_Parameters::ABC_order;
1052 
1053  //Add the flux boundary element to the helmholtz_outer_boundary_mesh
1054  Helmholtz_outer_boundary_mesh_pt->add_element_pt(flux_element_pt);
1055  }
1056 
1057 } // end of create_helmholtz_ABC_elements
1058 
1059 
1060 
1061 //=====================start_of_setup_interaction======================
1062 /// Setup interaction between two fields
1063 //========================================================================
1064 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
1067 {
1068 
1069  // Setup Helmholtz "pressure" load on traction elements
1070  unsigned boundary_in_helmholtz_mesh=0;
1071 
1072  // Doc boundary coordinate for Helmholtz
1073  ofstream the_file;
1074  the_file.open("boundary_coordinate_hh.dat");
1075  Helmholtz_mesh_pt->Mesh::doc_boundary_coordinates<HELMHOLTZ_ELEMENT>
1076  (boundary_in_helmholtz_mesh, the_file);
1077  the_file.close();
1078 
1079  Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
1080  <HELMHOLTZ_ELEMENT,2>
1081  (this,boundary_in_helmholtz_mesh,Helmholtz_mesh_pt,FSI_traction_mesh_pt);
1082 
1083  // Setup Helmholtz flux from normal displacement interaction
1084  unsigned boundary_in_solid_mesh=Outer_boundary_id;
1085 
1086  // Doc boundary coordinate for solid mesh
1087  the_file.open("boundary_coordinate_solid.dat");
1088  Solid_mesh_pt->Mesh::template doc_boundary_coordinates<ELASTICITY_ELEMENT>
1089  (boundary_in_solid_mesh, the_file);
1090  the_file.close();
1091 
1092  Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
1093  <ELASTICITY_ELEMENT,2>(
1094  this,boundary_in_solid_mesh,Solid_mesh_pt,Helmholtz_fsi_flux_mesh_pt);
1095 }
1096 
1097 
1098 
1099 //==============start_doc===========================================
1100 /// Doc the solution
1101 //==================================================================
1102 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
1104 {
1105 
1106  ofstream some_file,some_file2;
1107  char filename[100];
1108 
1109  // Number of plot points
1110  unsigned n_plot=5;
1111 
1112  // Compute/output the radiated power
1113  //----------------------------------
1114  sprintf(filename,"%s/power%i.dat",Doc_info.directory().c_str(),
1115  Doc_info.number());
1116  some_file.open(filename);
1117 
1118  // Accumulate contribution from elements
1119  double power=0.0;
1120  unsigned nn_element=Helmholtz_outer_boundary_mesh_pt->nelement();
1121  for(unsigned e=0;e<nn_element;e++)
1122  {
1123  HelmholtzBCElementBase<HELMHOLTZ_ELEMENT> *el_pt =
1124  dynamic_cast<HelmholtzBCElementBase<HELMHOLTZ_ELEMENT>*>(
1125  Helmholtz_outer_boundary_mesh_pt->element_pt(e));
1126  power += el_pt->global_power_contribution(some_file);
1127  }
1128  some_file.close();
1129  oomph_info << "Step: " << Doc_info.number()
1130  << ": Q=" << Global_Parameters::Q << "\n"
1131  << " k_squared=" << Global_Parameters::K_squared << "\n"
1132  << " density ratio (annulus) ="
1133  << Global_Parameters::Density_ratio[0] << "\n"
1134  << " density ratio (rib) ="
1135  << Global_Parameters::Density_ratio[1] << "\n"
1136  << " omega_sq (annulus)=" << Global_Parameters::Omega_sq[0] << "\n"
1137  << " omega_sq (rib )=" << Global_Parameters::Omega_sq[1] << "\n"
1138  << " Total radiated power " << power << "\n"
1139  << std::endl;
1140 
1141 
1142  // Write trace file
1143  Trace_file << Global_Parameters::Q << " "
1147  << Global_Parameters::Omega_sq[0] << " "
1148  << Global_Parameters::Omega_sq[1] << " "
1149  << power << " "
1150  << std::endl;
1151 
1152  std::ostringstream case_string;
1153  case_string << "TEXT X=10,Y=90, T=\"Q="
1155  << ", k<sup>2</sup> = "
1157  << ", density ratio = "
1158  << Global_Parameters::Density_ratio[0] << " and "
1160  << ", omega_sq = "
1161  << Global_Parameters::Omega_sq[0] << " and "
1162  << Global_Parameters::Omega_sq[1] << " "
1163  << "\"\n";
1164 
1165 
1166  // Output displacement field
1167  //--------------------------
1168  sprintf(filename,"%s/elast_soln%i.dat",Doc_info.directory().c_str(),
1169  Doc_info.number());
1170  some_file.open(filename);
1171  Solid_mesh_pt->output(some_file,n_plot);
1172  some_file.close();
1173 
1174  // Output solid traction elements
1175  //--------------------------------
1176  sprintf(filename,"%s/solid_traction_soln%i.dat",Doc_info.directory().c_str(),
1177  Doc_info.number());
1178  some_file.open(filename);
1179  Solid_traction_mesh_pt->output(some_file,n_plot);
1180  some_file.close();
1181 
1182  // Output fsi traction elements
1183  //-----------------------------
1184  sprintf(filename,"%s/traction_soln%i.dat",Doc_info.directory().c_str(),
1185  Doc_info.number());
1186  some_file.open(filename);
1187  FSI_traction_mesh_pt->output(some_file,n_plot);
1188  some_file.close();
1189 
1190 
1191  // Output Helmholtz fsi flux elements
1192  //-----------------------------------
1193  sprintf(filename,"%s/flux_bc_soln%i.dat",Doc_info.directory().c_str(),
1194  Doc_info.number());
1195  some_file.open(filename);
1196  Helmholtz_fsi_flux_mesh_pt->output(some_file,n_plot);
1197  some_file.close();
1198 
1199 
1200  // Output Helmholtz
1201  //-----------------
1202  sprintf(filename,"%s/helmholtz_soln%i.dat",Doc_info.directory().c_str(),
1203  Doc_info.number());
1204  some_file.open(filename);
1205  Helmholtz_mesh_pt->output(some_file,n_plot);
1206  some_file << case_string.str();
1207  some_file.close();
1208 
1209  // Output regions of solid mesh
1210  //-----------------------------
1211  unsigned nreg=Solid_mesh_pt->nregion();
1212  for (unsigned r=0;r<nreg;r++)
1213  {
1214  sprintf(filename,"%s/region%i_%i.dat",Doc_info.directory().c_str(),
1215  r,Doc_info.number());
1216  some_file.open(filename);
1217  unsigned nel=Solid_mesh_pt->nregion_element(r);
1218  for (unsigned e=0;e<nel;e++)
1219  {
1220  FiniteElement* el_pt=Solid_mesh_pt->region_element_pt(r,e);
1221  el_pt->output(some_file,n_plot);
1222  }
1223  some_file.close();
1224  }
1225 
1226 
1227  // Do animation of Helmholtz solution
1228  //-----------------------------------
1229  unsigned nstep=40;
1230  for (unsigned i=0;i<nstep;i++)
1231  {
1232  sprintf(filename,"%s/helmholtz_animation%i_frame%i.dat",
1233  Doc_info.directory().c_str(),
1234  Doc_info.number(),i);
1235  some_file.open(filename);
1236  double phi=2.0*MathematicalConstants::Pi*double(i)/double(nstep-1);
1237  unsigned nelem=Helmholtz_mesh_pt->nelement();
1238  for (unsigned e=0;e<nelem;e++)
1239  {
1240  HELMHOLTZ_ELEMENT* el_pt=dynamic_cast<HELMHOLTZ_ELEMENT*>(
1241  Helmholtz_mesh_pt->element_pt(e));
1242  el_pt->output_real(some_file,phi,n_plot);
1243  }
1244  some_file.close();
1245  }
1246 
1247  cout << "Doced for Q=" << Global_Parameters::Q << " (step "
1248  << Doc_info.number() << ")\n";
1249 
1250 // Increment label for output files
1251  Doc_info.number()++;
1252 
1253 } //end doc
1254 
1255 
1256 
1257 //=======start_of_main==================================================
1258 /// Driver for acoustic fsi problem
1259 //======================================================================
1260 int main(int argc, char **argv)
1261 {
1262 
1263  // Store command line arguments
1264  CommandLineArgs::setup(argc,argv);
1265 
1266  // Define possible command line arguments and parse the ones that
1267  // were actually specified
1268 
1269  // Output directory
1270  CommandLineArgs::specify_command_line_flag("--dir",
1272 
1273  // Peakiness parameter for loading
1274  CommandLineArgs::specify_command_line_flag("--alpha",
1276 
1277  // Multiplier for number of elements in solid mesh
1278  CommandLineArgs::specify_command_line_flag("--el_multiplier",
1280 
1281  // Outer radius of Helmholtz domain
1282  CommandLineArgs::specify_command_line_flag("--outer_radius",
1284 
1285  // Validaton run?
1286  CommandLineArgs::specify_command_line_flag("--validation");
1287 
1288  // Max. number of adaptations
1289  unsigned max_adapt=3;
1290  CommandLineArgs::specify_command_line_flag("--max_adapt",&max_adapt);
1291 
1292  // Parse command line
1293  CommandLineArgs::parse_and_assign();
1294 
1295  // Doc what has actually been specified on the command line
1296  CommandLineArgs::doc_specified_flags();
1297 
1298  //Set up the problem
1299  CoatedDiskProblem<ProjectableTimeHarmonicLinearElasticityElement
1300  <TTimeHarmonicLinearElasticityElement<2,3> >,
1301  RefineableQHelmholtzElement<2,3> > problem;
1302 
1303 
1304  // Set values for parameter values
1305  Global_Parameters::Q=5.0;
1309 
1310  //Parameter study
1311  unsigned nstep=3;
1312  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
1313  {
1314  nstep=1;
1315  max_adapt=2;
1316  }
1317  for(unsigned i=0;i<nstep;i++)
1318  {
1319  // Solve the problem with Newton's method, allowing
1320  // up to max_adapt mesh adaptations after every solve.
1321  problem.newton_solve(max_adapt);
1322 
1323  // Doc solution
1324  problem.doc_solution();
1325 
1326  // Make rib a lot heavier but keep its stiffness
1327  if (i==0)
1328  {
1329  Global_Parameters::E_pt[1]->update_constitutive_parameters(
1330  Global_Parameters::Nu,1.0);
1333  }
1334 
1335  // Make rib very soft and inertia-less
1336  if (i==1)
1337  {
1338  Global_Parameters::E_pt[1]->update_constitutive_parameters(
1339  Global_Parameters::Nu,1.0e-16);
1342  }
1343 
1344 
1345  }
1346 
1347 } //end of main
1348 
1349 
1350 
1351 
1352 
1353 
1354 
1355 
double Omega_sq
Non-dim square of frequency for solid – dependent variable!
Definition: acoustic_fsi.cc:77
double Q
FSI parameter.
Definition: acoustic_fsi.cc:62
CoatedDiskProblem()
Constructor:
Vector< double > Omega_sq(2, 0.0)
Square of non-dim frequency for the two regions – dependent variable!
RefineableTriangleMesh< ELASTICITY_ELEMENT > * Solid_mesh_pt
Pointer to refineable solid mesh.
double Outer_radius
Radius of outer boundary of Helmholtz domain.
Definition: acoustic_fsi.cc:59
int main(int argc, char **argv)
Driver for acoustic fsi problem.
Vector< double > Density_ratio(2, 0.1)
Density ratio for the two regions: solid to fluid.
void update_parameter_values()
Function to update dependent parameter values.
Definition: acoustic_fsi.cc:80
double Density_ratio
Density ratio: solid to fluid.
Definition: acoustic_fsi.cc:74
unsigned Upper_symmetry_boundary_id
Boundary ID of upper symmetry boundary.
void doc_solution()
Doc the solution.
void create_solid_traction_elements()
Create solid traction elements.
void operator=(const MyStraightLine &)
Broken assignment operator.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of traction elements.
Vector< TimeHarmonicIsotropicElasticityTensor * > E_pt
The elasticity tensors for the two regions.
unsigned Upper_inner_boundary_id
Boundary ID of upper inner boundary.
Vector< double > R_start
Start point of line.
MyStraightLine(const Vector< double > &r_start, const Vector< double > &r_end)
Constructor: Pass start and end points.
double K_squared
Square of wavenumber for the Helmholtz equation.
Definition: acoustic_fsi.cc:56
Coated disk FSI.
~MyStraightLine()
Destructor: Empty.
void create_helmholtz_fsi_flux_elements()
Create Helmholtz FSI flux elements.
void pressure_load(const Vector< double > &x, const Vector< double > &n, Vector< std::complex< double > > &traction)
Pressure load (real and imag part)
void setup_interaction()
Setup interaction.
Mesh * Solid_traction_mesh_pt
Pointer to mesh of solid traction elements.
string Directory
Output directory.
unsigned Lower_inner_boundary_id
Boundary ID of lower inner boundary.
unsigned El_multiplier
Multiplier for number of elements.
double H_coating
Non-dim thickness of elastic coating.
Definition: acoustic_fsi.cc:65
Vector< double > R_end
End point of line.
Straight 1D line in 2D space.
void delete_face_elements(Mesh *const &boundary_mesh_pt)
Delete (face) elements in specified mesh.
double Nu
Poisson's ratio.
Definition: acoustic_fsi.cc:68
void create_helmholtz_ABC_elements()
Create ABC face elements.
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of traction elements.
unsigned ABC_order
Order of absorbing/appproximate boundary condition.
MyStraightLine(const MyStraightLine &dummy)
Broken copy constructor.
void actions_before_newton_solve()
Update function (empty)
double Alpha
Peakiness parameter for pressure load.
void actions_after_newton_solve()
Update function (empty)
double P
Uniform pressure.
unsigned Outer_boundary_id
Boundary ID of outer boundary.
void create_fsi_traction_elements()
Create FSI traction elements.
void position(const Vector< double > &zeta, Vector< double > &r) const
Position Vector at Lagrangian coordinate zeta.
Mesh * Helmholtz_outer_boundary_mesh_pt
Pointer to mesh containing the ABC elements.
unsigned Lower_symmetry_boundary_id
Boundary ID of lower symmetry boundary.