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