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: 1176 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-11 17:57:27 +0100 (Wed, 11 May 2016) $
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 mesh
39 #include "meshes/annular_mesh.h"
40 
41 using namespace std;
42 using namespace oomph;
43 
44 //////////////////////////////////////////////////////////////////////
45 //////////////////////////////////////////////////////////////////////
46 //////////////////////////////////////////////////////////////////////
47 
48 
49 //=======start_namespace==========================================
50 /// Global variables
51 //================================================================
52 namespace Global_Parameters
53 {
54 
55  /// \short Square of wavenumber for the Helmholtz equation
56  double K_squared=10.0;
57 
58  /// \short Radius of outer boundary of Helmholtz domain
59  double Outer_radius=4.0;
60 
61  /// FSI parameter
62  double Q=0.0;
63 
64  /// Non-dim thickness of elastic coating
65  double H_coating=0.3;
66 
67  /// Poisson's ratio
68  double Nu = 0.3;
69 
70  /// The elasticity tensor for the solid
71  TimeHarmonicIsotropicElasticityTensor E(Nu);
72 
73  /// Density ratio: solid to fluid
74  double Density_ratio=0.0;
75 
76  /// Non-dim square of frequency for solid -- dependent variable!
77  double Omega_sq=0.0;
78 
79  /// Function to update dependent parameter values
81  {
83  }
84 
85  /// \short Azimuthal wavenumber for imposed displacement of coating
86  /// on inner boundary
87  unsigned N=0;
88 
89  /// \short Displacement field on inner boundary of solid
90  void solid_boundary_displacement(const Vector<double>& x,
91  Vector<std::complex<double> >& u)
92  {
93  Vector<double> normal(2);
94  double norm=sqrt(x[0]*x[0]+x[1]*x[1]);
95  double phi=atan2(x[1],x[0]);
96  normal[0]=x[0]/norm;
97  normal[1]=x[1]/norm;
98 
99  u[0]=complex<double>(normal[0]*cos(double(N)*phi),0.0);
100  u[1]=complex<double>(normal[1]*cos(double(N)*phi),0.0);
101  }
102 
103 
104  /// Output directory
105  string Directory="RESLT";
106 
107  /// Multiplier for number of elements
108  unsigned El_multiplier=1;
109 
110  /// Interface to Hankel function in maple style
111  std::complex<double> HankelH1(const double& k, const double& x)
112  {
113  Vector<std::complex<double> > h(2);
114  Vector<std::complex<double> > hp(2);
115  Hankel_functions_for_helmholtz_problem::Hankel_first(2,x,h,hp);
116 
117  if (k==0.0)
118  {
119  return h[0];
120  }
121  else if (k==1.0)
122  {
123  return h[1];
124  }
125  else
126  {
127  cout << "Never get here. k=" << k << std::endl;
128  assert(false);
129  // Dummy return
130  return std::complex<double>(1.0,1.0);
131  }
132  }
133 
134 
135  /// \short Coefficient in front of Hankel function for axisymmetric solution
136  /// of Helmholtz potential
137  std::complex<double> axisym_coefficient()
138  {
139  std::complex<double> MapleGenVar1;
140  std::complex<double> MapleGenVar2;
141  std::complex<double> MapleGenVar3;
142  std::complex<double> MapleGenVar4;
143  std::complex<double> MapleGenVar5;
144  std::complex<double> MapleGenVar6;
145  std::complex<double> MapleGenVar7;
146  std::complex<double> MapleGenVar8;
147  std::complex<double> t0;
148 
149 
150  MapleGenVar3 = (-2.0/(2.0+2.0*Nu)*1.0+2.0/(2.0+2.0*Nu)*1.0*
151 H_coating)/(Nu/(1.0+Nu)/(1.0-2.0*Nu)+2.0/(2.0+2.0*Nu)-2.0/(2.0+2.0*Nu)*
152 H_coating+1/(2.0+2.0*Nu)*H_coating*H_coating)/2.0;
153  MapleGenVar5 = -1/(Nu/(1.0+Nu)/(1.0-2.0*Nu)+2.0/(2.0+2.0*Nu)-2.0/(2.0+2.0
154 *Nu)*H_coating+1/(2.0+2.0*Nu)*H_coating*H_coating)*Q/2.0;
155  MapleGenVar6 = HankelH1(0.0,sqrt(K_squared))*(-(-2.0/(2.0+2.0*Nu)*1.0
156 +2.0/(2.0+2.0*Nu)*1.0*H_coating)/(Nu/(1.0+Nu)/(1.0-2.0*Nu)+2.0/(2.0+2.0*Nu)
157 -2.0/(2.0+2.0*Nu)*H_coating+1/(2.0+2.0*Nu)*H_coating*H_coating)/2.0+(2.0/(2.0+
158 2.0*Nu)*1.0-2.0/(2.0+2.0*Nu)*1.0*H_coating+2.0*1.0*Nu/(1.0+Nu)/(1.0
159 -2.0*Nu)-2.0*1.0*H_coating*Nu/(1.0+Nu)/(1.0-2.0*Nu))/(Nu/(1.0+Nu)/(1.0-2.0*
160 Nu)+2.0/(2.0+2.0*Nu)-2.0/(2.0+2.0*Nu)*H_coating+1/(2.0+2.0*Nu)*H_coating*
161 H_coating)/2.0)/(HankelH1(1.0,sqrt(K_squared))*sqrt(K_squared)-Q*HankelH1(0.0,
162 sqrt(K_squared))/(Nu/(1.0+Nu)/(1.0-2.0*Nu)+2.0/(2.0+2.0*Nu)-2.0/(2.0+2.0*Nu)*
163 H_coating+1/(2.0+2.0*Nu)*H_coating*H_coating)/2.0+Q*HankelH1(0.0,sqrt(K_squared
164 ))*(1.0-2.0*H_coating+H_coating*H_coating)/(Nu/(1.0+Nu)/(1.0-2.0*Nu)+2.0/(2.0+
165 2.0*Nu)-2.0/(2.0+2.0*Nu)*H_coating+1/(2.0+2.0*Nu)*H_coating*H_coating)/2.0);
166  MapleGenVar4 = MapleGenVar5*MapleGenVar6;
167  MapleGenVar2 = MapleGenVar3+MapleGenVar4;
168  MapleGenVar3 = MapleGenVar2;
169  MapleGenVar5 = -(2.0/(2.0+2.0*Nu)*1.0-2.0/(2.0+2.0*Nu)*1.0*
170 H_coating+2.0*1.0*Nu/(1.0+Nu)/(1.0-2.0*Nu)-2.0*1.0*H_coating*Nu/(1.0+Nu
171 )/(1.0-2.0*Nu))/(Nu/(1.0+Nu)/(1.0-2.0*Nu)+2.0/(2.0+2.0*Nu)-2.0/(2.0+2.0*Nu)*
172 H_coating+1/(2.0+2.0*Nu)*H_coating*H_coating)/2.0;
173  MapleGenVar7 = (1.0-2.0*H_coating+H_coating*H_coating)/(Nu/(1.0+Nu)/(1.0
174 -2.0*Nu)+2.0/(2.0+2.0*Nu)-2.0/(2.0+2.0*Nu)*H_coating+1/(2.0+2.0*Nu)*H_coating*
175 H_coating)/2.0;
176  MapleGenVar8 = Q*HankelH1(0.0,sqrt(K_squared))*(-(-2.0/(2.0+2.0*Nu)*
177 1.0+2.0/(2.0+2.0*Nu)*1.0*H_coating)/(Nu/(1.0+Nu)/(1.0-2.0*Nu)+2.0/(2.0+
178 2.0*Nu)-2.0/(2.0+2.0*Nu)*H_coating+1/(2.0+2.0*Nu)*H_coating*H_coating)/2.0+(2.0
179 /(2.0+2.0*Nu)*1.0-2.0/(2.0+2.0*Nu)*1.0*H_coating+2.0*1.0*Nu/(1.0+Nu
180 )/(1.0-2.0*Nu)-2.0*1.0*H_coating*Nu/(1.0+Nu)/(1.0-2.0*Nu))/(Nu/(1.0+Nu)/(
181 1.0-2.0*Nu)+2.0/(2.0+2.0*Nu)-2.0/(2.0+2.0*Nu)*H_coating+1/(2.0+2.0*Nu)*
182 H_coating*H_coating)/2.0)/(HankelH1(1.0,sqrt(K_squared))*sqrt(K_squared)-Q*
183 HankelH1(0.0,sqrt(K_squared))/(Nu/(1.0+Nu)/(1.0-2.0*Nu)+2.0/(2.0+2.0*Nu)-2.0/(
184 2.0+2.0*Nu)*H_coating+1/(2.0+2.0*Nu)*H_coating*H_coating)/2.0+Q*HankelH1(0.0,
185 sqrt(K_squared))*(1.0-2.0*H_coating+H_coating*H_coating)/(Nu/(1.0+Nu)/(1.0-2.0*
186 Nu)+2.0/(2.0+2.0*Nu)-2.0/(2.0+2.0*Nu)*H_coating+1/(2.0+2.0*Nu)*H_coating*
187 H_coating)/2.0);
188  MapleGenVar6 = MapleGenVar7*MapleGenVar8;
189  MapleGenVar4 = MapleGenVar5+MapleGenVar6;
190  MapleGenVar1 = MapleGenVar3+MapleGenVar4;
191  MapleGenVar2 = 1.0/HankelH1(1.0,sqrt(K_squared))/sqrt(K_squared);
192  t0 = MapleGenVar1*MapleGenVar2;
193 
194  return t0;
195  }
196 
197 
198  /// Exact solution for Helmholtz potential for axisymmetric solution
199  void exact_axisym_potential(const Vector<double>& x,
200  Vector<double>& soln)
201  {
202  std::complex<double> C=axisym_coefficient();
203  double r=sqrt(x[0]*x[0]+x[1]*x[1]);
204  soln[0]=real(HankelH1(0,sqrt(K_squared)*r)*C);
205  soln[1]=imag(HankelH1(0,sqrt(K_squared)*r)*C);
206  }
207 
208  /// Exact radiated power for axisymmetric solution
210  {
211  // Solution is independent of where it's evaluated (have checked!)
212  double r=1.0;
213 
214  // Get the coefficient for the axisymmetric potential
215  std::complex<double> C=axisym_coefficient();
216 
217  // Argument for Bessel/Hankel functions
218  double rr=sqrt(K_squared)*r;
219 
220  // Evaluate Hankel functions -- only need zero-th term
221  Vector<std::complex<double> > h(2),hp(2);
222  Hankel_functions_for_helmholtz_problem::Hankel_first(1,rr,h,hp);
223 
224  // Compute time-average radiated power
225  double power=sqrt(K_squared)*0.5*r*2.0*MathematicalConstants::Pi*
226  (imag(C*hp[0])*real(C*h[0])-real(C*hp[0])*imag(C*h[0]));
227 
228  return power;
229  }
230 
231 } //end namespace
232 
233 
234 
235 //=============begin_problem============================================
236 /// Coated disk FSI
237 //======================================================================
238 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
239 class CoatedDiskProblem : public Problem
240 {
241 
242 public:
243 
244  /// Constructor:
246 
247  /// Update function (empty)
249 
250  /// Update function (empty)
252 
253  /// Recompute gamma integral before checking Newton residuals
255  {
256  Helmholtz_outer_boundary_mesh_pt->setup_gamma();
257  }
258 
259  /// Actions before adapt: Wipe the mesh of traction elements
260  void actions_before_adapt();
261 
262  /// Actions after adapt: Rebuild the mesh of traction elements
263  void actions_after_adapt();
264 
265  /// Doc the solution
266  void doc_solution();
267 
268 private:
269 
270  /// \short Create FSI traction elements
271  void create_fsi_traction_elements();
272 
273  /// \short Create Helmholtz FSI flux elements
274  void create_helmholtz_fsi_flux_elements();
275 
276  /// Delete (face) elements in specified mesh
277  void delete_face_elements(Mesh* const & boundary_mesh_pt);
278 
279  /// \short Create DtN face elements
280  void create_helmholtz_DtN_elements();
281 
282  /// Setup interaction
283  void setup_interaction();
284 
285  /// Pointer to solid mesh
286  TreeBasedRefineableMeshBase* Solid_mesh_pt;
287 
288  /// Pointer to mesh of FSI traction elements
290 
291  /// Pointer to Helmholtz mesh
292  TreeBasedRefineableMeshBase* Helmholtz_mesh_pt;
293 
294  /// Pointer to mesh of Helmholtz FSI flux elements
296 
297  /// \short Pointer to mesh containing the DtN elements
298  HelmholtzDtNMesh<HELMHOLTZ_ELEMENT>* Helmholtz_outer_boundary_mesh_pt;
299 
300  /// DocInfo object for output
301  DocInfo Doc_info;
302 
303  /// Trace file
304  ofstream Trace_file;
305 
306 };
307 
308 
309 //===========start_of_constructor=======================================
310 /// Constructor
311 //======================================================================
312 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
314 {
315 
316  // The coating mesh is periodic
317  bool periodic=true;
318  double azimuthal_fraction_of_coating=1.0;
319 
320  // Solid mesh
321  //-----------
322  // Number of elements in azimuthal direction
323  unsigned ntheta_solid=10*Global_Parameters::El_multiplier;
324 
325  // Number of elements in radial direction
326  unsigned nr_solid=3*Global_Parameters::El_multiplier;
327 
328  // Innermost radius for solid mesh
329  double a=1.0-Global_Parameters::H_coating;
330 
331  // Build solid mesh
332  Solid_mesh_pt = new
333  RefineableTwoDAnnularMesh<ELASTICITY_ELEMENT>
334  (periodic,azimuthal_fraction_of_coating,
335  ntheta_solid,nr_solid,a,Global_Parameters::H_coating);
336 
337 
338  // Helmholtz mesh
339  //---------------
340 
341  // Number of elements in azimuthal direction in Helmholtz mesh
342  unsigned ntheta_helmholtz=11*Global_Parameters::El_multiplier;
343 
344  // Number of elements in radial direction in Helmholtz mesh
345  unsigned nr_helmholtz=3*Global_Parameters::El_multiplier;
346 
347  // Innermost radius of Helmholtz mesh
348  a=1.0;
349 
350  // Thickness of Helmholtz mesh
351  double h_thick_helmholtz=Global_Parameters::Outer_radius-a;
352 
353  // Build mesh
354  Helmholtz_mesh_pt = new
355  RefineableTwoDAnnularMesh<HELMHOLTZ_ELEMENT>
356  (periodic,azimuthal_fraction_of_coating,
357  ntheta_helmholtz,nr_helmholtz,a,h_thick_helmholtz);
358 
359 
360  // Set error estimators
361  Solid_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
362  Helmholtz_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
363 
364 
365  // Mesh containing the Helmholtz DtN
366  // elements. Specify outer radius and number of Fourier terms to be
367  // used in gamma integral
368  unsigned nfourier=20;
369  Helmholtz_outer_boundary_mesh_pt =
370  new HelmholtzDtNMesh<HELMHOLTZ_ELEMENT>(Global_Parameters::Outer_radius,
371  nfourier);
372 
373  //Assign the physical properties to the elements before any refinement
374  //Loop over the elements in the solid mesh
375  unsigned n_element=Solid_mesh_pt->nelement();
376  for(unsigned i=0;i<n_element;i++)
377  {
378  //Cast to a solid element
379  ELASTICITY_ELEMENT *el_pt =
380  dynamic_cast<ELASTICITY_ELEMENT*>(Solid_mesh_pt->element_pt(i));
381 
382  // Set the constitutive law
383  el_pt->elasticity_tensor_pt() = &Global_Parameters::E;
384 
385  // Square of non-dim frequency
386  el_pt->omega_sq_pt()= &Global_Parameters::Omega_sq;
387  }
388 
389 
390  // Same for Helmholtz mesh
391  n_element =Helmholtz_mesh_pt->nelement();
392  for(unsigned i=0;i<n_element;i++)
393  {
394  //Cast to a solid element
395  HELMHOLTZ_ELEMENT *el_pt =
396  dynamic_cast<HELMHOLTZ_ELEMENT*>(Helmholtz_mesh_pt->element_pt(i));
397 
398  //Set the pointer to square of Helmholtz wavenumber
399  el_pt->k_squared_pt() = &Global_Parameters::K_squared;
400  }
401 
402 
403  // Output meshes and their boundaries so far so we can double
404  // check the boundary enumeration
405  Solid_mesh_pt->output("solid_mesh.dat");
406  Helmholtz_mesh_pt->output("helmholtz_mesh.dat");
407  Solid_mesh_pt->output_boundaries("solid_mesh_boundary.dat");
408  Helmholtz_mesh_pt->output_boundaries("helmholtz_mesh_boundary.dat");
409 
410 
411  // Create FaceElement meshes for boundary conditions
412  //---------------------------------------------------
413 
414  // Construct the fsi traction element mesh
415  FSI_traction_mesh_pt=new Mesh;
416  create_fsi_traction_elements();
417 
418  // Construct the Helmholtz fsi flux element mesh
419  Helmholtz_fsi_flux_mesh_pt=new Mesh;
420  create_helmholtz_fsi_flux_elements();
421 
422  // Create DtN elements on outer boundary of Helmholtz mesh
423  create_helmholtz_DtN_elements();
424 
425 
426  // Combine sub meshes
427  //-------------------
428 
429  // Solid mesh is first sub-mesh
430  add_sub_mesh(Solid_mesh_pt);
431 
432  // Add traction sub-mesh
433  add_sub_mesh(FSI_traction_mesh_pt);
434 
435  // Add Helmholtz mesh
436  add_sub_mesh(Helmholtz_mesh_pt);
437 
438  // Add Helmholtz FSI flux mesh
439  add_sub_mesh(Helmholtz_fsi_flux_mesh_pt);
440 
441  // Add Helmholtz DtN mesh
442  add_sub_mesh(Helmholtz_outer_boundary_mesh_pt);
443 
444  // Build combined "global" mesh
445  build_global_mesh();
446 
447 
448  // Solid boundary conditions:
449  //---------------------------
450  // Pin displacements on innermost boundary (boundary 0) of solid mesh
451  unsigned n_node = Solid_mesh_pt->nboundary_node(0);
452  Vector<std::complex<double> > u(2);
453  Vector<double> x(2);
454  for(unsigned i=0;i<n_node;i++)
455  {
456  Node* nod_pt=Solid_mesh_pt->boundary_node_pt(0,i);
457  nod_pt->pin(0);
458  nod_pt->pin(1);
459  nod_pt->pin(2);
460  nod_pt->pin(3);
461 
462  // Assign displacements
463  x[0]=nod_pt->x(0);
464  x[1]=nod_pt->x(1);
466 
467  // Real part of x-displacement
468  nod_pt->set_value(0,u[0].real());
469 
470  // Imag part of x-displacement
471  nod_pt->set_value(1,u[1].real());
472 
473  // Real part of y-displacement
474  nod_pt->set_value(2,u[0].imag());
475 
476  //Imag part of y-displacement
477  nod_pt->set_value(3,u[1].imag());
478  }
479 
480 
481  // Setup fluid-structure interaction
482  //----------------------------------
483  setup_interaction();
484 
485  // Assign equation numbers
486  oomph_info << "Number of unknowns: " << assign_eqn_numbers() << std::endl;
487 
488  // Set output directory
489  Doc_info.set_directory(Global_Parameters::Directory);
490 
491  // Open trace file
492  char filename[100];
493  sprintf(filename,"%s/trace.dat",Doc_info.directory().c_str());
494  Trace_file.open(filename);
495 
496 
497 } //end of constructor
498 
499 
500 //=====================start_of_actions_before_adapt======================
501 /// Actions before adapt: Wipe the meshes face elements
502 //========================================================================
503 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
506 {
507  // Kill the fsi traction elements and wipe surface mesh
508  delete_face_elements(FSI_traction_mesh_pt);
509 
510  // Kill Helmholtz FSI flux elements
511  delete_face_elements(Helmholtz_fsi_flux_mesh_pt);
512 
513  // Kill Helmholtz BC elements
514  delete_face_elements(Helmholtz_outer_boundary_mesh_pt);
515 
516  // Rebuild the Problem's global mesh from its various sub-meshes
517  rebuild_global_mesh();
518 
519 }// end of actions_before_adapt
520 
521 
522 
523 //=====================start_of_actions_after_adapt=======================
524 /// Actions after adapt: Rebuild the meshes of face elements
525 //========================================================================
526 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
529 {
530  // Create fsi traction elements from all elements that are
531  // adjacent to FSI boundaries and add them to surface meshes
532  create_fsi_traction_elements();
533 
534  // Create Helmholtz fsi flux elements
535  create_helmholtz_fsi_flux_elements();
536 
537  // Create DtN elements from all elements that are
538  // adjacent to the outer boundary of Helmholtz mesh
539  create_helmholtz_DtN_elements();
540 
541  // Setup interaction
542  setup_interaction();
543 
544  // Rebuild the Problem's global mesh from its various sub-meshes
545  rebuild_global_mesh();
546 
547 }// end of actions_after_adapt
548 
549 
550 //============start_of_delete_face_elements================
551 /// Delete face elements and wipe the mesh
552 //==========================================================
553 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
555 delete_face_elements(Mesh* const & boundary_mesh_pt)
556 {
557  // How many surface elements are in the surface mesh
558  unsigned n_element = boundary_mesh_pt->nelement();
559 
560  // Loop over the surface elements
561  for(unsigned e=0;e<n_element;e++)
562  {
563  // Kill surface element
564  delete boundary_mesh_pt->element_pt(e);
565  }
566 
567  // Wipe the mesh
568  boundary_mesh_pt->flush_element_and_node_storage();
569 
570 } // end of delete_face_elements
571 
572 
573 
574 //============start_of_create_fsi_traction_elements======================
575 /// Create fsi traction elements
576 //=======================================================================
577 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
580 {
581  // We're on boundary 2 of the solid mesh
582  unsigned b=2;
583 
584  // How many bulk elements are adjacent to boundary b?
585  unsigned n_element = Solid_mesh_pt->nboundary_element(b);
586 
587  // Loop over the bulk elements adjacent to boundary b
588  for(unsigned e=0;e<n_element;e++)
589  {
590  // Get pointer to the bulk element that is adjacent to boundary b
591  ELASTICITY_ELEMENT* bulk_elem_pt = dynamic_cast<ELASTICITY_ELEMENT*>(
592  Solid_mesh_pt->boundary_element_pt(b,e));
593 
594  //Find the index of the face of element e along boundary b
595  int face_index = Solid_mesh_pt->face_index_at_boundary(b,e);
596 
597  // Create element
598  TimeHarmonicLinElastLoadedByHelmholtzPressureBCElement
599  <ELASTICITY_ELEMENT,HELMHOLTZ_ELEMENT>* el_pt=
600  new TimeHarmonicLinElastLoadedByHelmholtzPressureBCElement
601  <ELASTICITY_ELEMENT,HELMHOLTZ_ELEMENT>(bulk_elem_pt,
602  face_index);
603  // Add to mesh
604  FSI_traction_mesh_pt->add_element_pt(el_pt);
605 
606  // Associate element with bulk boundary (to allow it to access
607  // the boundary coordinates in the bulk mesh)
608  el_pt->set_boundary_number_in_bulk_mesh(b);
609 
610  // Set FSI parameter
611  el_pt->q_pt()=&Global_Parameters::Q;
612  }
613 
614 } // end of create_traction_elements
615 
616 
617 
618 
619 
620 //============start_of_create_helmholtz_fsi_flux_elements================
621 /// Create Helmholtz fsii flux elements
622 //=======================================================================
623 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
626 {
627 
628  // Attach to inner boundary of Helmholtz mesh (0)
629  unsigned b=0;
630 
631  // How many bulk elements are adjacent to boundary b?
632  unsigned n_element = Helmholtz_mesh_pt->nboundary_element(b);
633 
634  // Loop over the bulk elements adjacent to boundary b
635  for(unsigned e=0;e<n_element;e++)
636  {
637  // Get pointer to the bulk element that is adjacent to boundary b
638  HELMHOLTZ_ELEMENT* bulk_elem_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
639  Helmholtz_mesh_pt->boundary_element_pt(b,e));
640 
641  //Find the index of the face of element e along boundary b
642  int face_index = Helmholtz_mesh_pt->face_index_at_boundary(b,e);
643 
644  // Create element
645  HelmholtzFluxFromNormalDisplacementBCElement
646  <HELMHOLTZ_ELEMENT,ELASTICITY_ELEMENT>* el_pt=
647  new HelmholtzFluxFromNormalDisplacementBCElement
648  <HELMHOLTZ_ELEMENT,ELASTICITY_ELEMENT>(bulk_elem_pt,
649  face_index);
650 
651  // Add to mesh
652  Helmholtz_fsi_flux_mesh_pt->add_element_pt(el_pt);
653 
654  // Associate element with bulk boundary (to allow it to access
655  // the boundary coordinates in the bulk mesh)
656  el_pt->set_boundary_number_in_bulk_mesh(b);
657  }
658 
659 } // end of create_helmholtz_flux_elements
660 
661 
662 
663 //============start_of_create_DtN_elements==============================
664 /// Create DtN elements on the outer boundary of
665 /// the Helmholtz mesh
666 //===========================================================================
667 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
670 {
671  // We're on boundary 2 of the Helmholtz mesh
672  unsigned b=2;
673 
674  // How many bulk elements are adjacent to boundary b?
675  unsigned n_element = Helmholtz_mesh_pt->nboundary_element(b);
676 
677  // Loop over the bulk elements adjacent to boundary b?
678  for(unsigned e=0;e<n_element;e++)
679  {
680  // Get pointer to the bulk element that is adjacent to boundary b
681  HELMHOLTZ_ELEMENT* bulk_elem_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
682  Helmholtz_mesh_pt->boundary_element_pt(b,e));
683 
684  //Find the index of the face of element e along boundary b
685  int face_index = Helmholtz_mesh_pt->face_index_at_boundary(b,e);
686 
687  // Build the corresponding DtN element
688  HelmholtzDtNBoundaryElement<HELMHOLTZ_ELEMENT>* flux_element_pt = new
689  HelmholtzDtNBoundaryElement<HELMHOLTZ_ELEMENT>(bulk_elem_pt,face_index);
690 
691  // Set pointer to the mesh that contains all the boundary condition
692  // elements on this boundary
693  flux_element_pt->set_outer_boundary_mesh_pt(
694  Helmholtz_outer_boundary_mesh_pt);
695 
696  //Add the DtN element to the mesh
697  Helmholtz_outer_boundary_mesh_pt->add_element_pt(flux_element_pt);
698 
699  }
700 
701 } // end of create_helmholtz_DtN_elements
702 
703 
704 
705 //=====================start_of_setup_interaction======================
706 /// Setup interaction between two fields
707 //========================================================================
708 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
711 {
712  // Setup Helmholtz "pressure" load on traction elements
713  unsigned boundary_in_helmholtz_mesh=0;
714  Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
715  <HELMHOLTZ_ELEMENT,2>
716  (this,boundary_in_helmholtz_mesh,Helmholtz_mesh_pt,FSI_traction_mesh_pt);
717 
718  // Setup Helmholtz flux from normal displacement interaction
719  unsigned boundary_in_solid_mesh=2;
720  Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
721  <ELASTICITY_ELEMENT,2>(
722  this,boundary_in_solid_mesh,Solid_mesh_pt,Helmholtz_fsi_flux_mesh_pt);
723 }
724 
725 
726 
727 //==============start_doc===========================================
728 /// Doc the solution
729 //==================================================================
730 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
732 {
733 
734  ofstream some_file,some_file2;
735  char filename[100];
736 
737  // Number of plot points
738  unsigned n_plot=5;
739 
740  // Compute/output the radiated power
741  //----------------------------------
742  sprintf(filename,"%s/power%i.dat",Doc_info.directory().c_str(),
743  Doc_info.number());
744  some_file.open(filename);
745 
746  // Accumulate contribution from elements
747  double power=0.0;
748  unsigned nn_element=Helmholtz_outer_boundary_mesh_pt->nelement();
749  for(unsigned e=0;e<nn_element;e++)
750  {
751  HelmholtzBCElementBase<HELMHOLTZ_ELEMENT> *el_pt =
752  dynamic_cast<HelmholtzBCElementBase<HELMHOLTZ_ELEMENT>*>(
753  Helmholtz_outer_boundary_mesh_pt->element_pt(e));
754  power += el_pt->global_power_contribution(some_file);
755  }
756  some_file.close();
757  oomph_info << "Step: " << Doc_info.number()
758  << ": Q=" << Global_Parameters::Q << "\n"
759  << " k_squared=" << Global_Parameters::K_squared << "\n"
760  << " density ratio=" << Global_Parameters::Density_ratio << "\n"
761  << " omega_sq=" << Global_Parameters::Omega_sq << "\n"
762  << " Total radiated power " << power << "\n"
763  << " Axisymmetric radiated power " << "\n"
765  << std::endl;
766 
767 
768  // Write trace file
769  Trace_file << Global_Parameters::Q << " "
773  << power << " "
775  << std::endl;
776 
777  std::ostringstream case_string;
778  case_string << "TEXT X=10,Y=90, T=\"Q="
780  << ", k<sup>2</sup>="
782  << ", density ratio="
784  << ", omega_sq="
786  << "\"\n";
787 
788 
789  // Output displacement field
790  //--------------------------
791  sprintf(filename,"%s/elast_soln%i.dat",Doc_info.directory().c_str(),
792  Doc_info.number());
793  some_file.open(filename);
794  Solid_mesh_pt->output(some_file,n_plot);
795  some_file.close();
796 
797 
798  // Output fsi traction elements
799  //-----------------------------
800  sprintf(filename,"%s/traction_soln%i.dat",Doc_info.directory().c_str(),
801  Doc_info.number());
802  some_file.open(filename);
803  FSI_traction_mesh_pt->output(some_file,n_plot);
804  some_file.close();
805 
806 
807  // Output Helmholtz fsi flux elements
808  //-----------------------------------
809  sprintf(filename,"%s/flux_bc_soln%i.dat",Doc_info.directory().c_str(),
810  Doc_info.number());
811  some_file.open(filename);
812  Helmholtz_fsi_flux_mesh_pt->output(some_file,n_plot);
813  some_file.close();
814 
815 
816  // Output Helmholtz
817  //-----------------
818  sprintf(filename,"%s/helmholtz_soln%i.dat",Doc_info.directory().c_str(),
819  Doc_info.number());
820  some_file.open(filename);
821  Helmholtz_mesh_pt->output(some_file,n_plot);
822  some_file << case_string.str();
823  some_file.close();
824 
825 
826  // Output exact solution for Helmholtz
827  //------------------------------------
828  sprintf(filename,"%s/exact_helmholtz_soln%i.dat",Doc_info.directory().c_str(),
829  Doc_info.number());
830  some_file.open(filename);
831  Helmholtz_mesh_pt->output_fct(some_file,n_plot,
833  some_file.close();
834 
835 
836  cout << "Doced for Q=" << Global_Parameters::Q << " (step "
837  << Doc_info.number() << ")\n";
838 
839  // Increment label for output files
840  Doc_info.number()++;
841 
842 } //end doc
843 
844 
845 
846 //=======start_of_main==================================================
847 /// Driver for acoustic fsi problem
848 //======================================================================
849 int main(int argc, char **argv)
850 {
851 
852  // Store command line arguments
853  CommandLineArgs::setup(argc,argv);
854 
855  // Define possible command line arguments and parse the ones that
856  // were actually specified
857 
858  // Output directory
859  CommandLineArgs::specify_command_line_flag("--dir",
861 
862  // Azimuthal wavenumber of forcing
863  CommandLineArgs::specify_command_line_flag("--n",&Global_Parameters::N);
864 
865  // Minimum refinement level
866  CommandLineArgs::specify_command_line_flag("--el_multiplier",
868 
869  // Outer radius of Helmholtz domain
870  CommandLineArgs::specify_command_line_flag("--outer_radius",
872 
873  // Number of steps in parameter study
874  unsigned nstep=2;
875  CommandLineArgs::specify_command_line_flag("--nstep",&nstep);
876 
877  // Increment in FSI parameter in parameter study
878  double q_increment=5.0;
879  CommandLineArgs::specify_command_line_flag("--q_increment",&q_increment);
880 
881  // Max. number of adaptations
882  unsigned max_adapt=3;
883  CommandLineArgs::specify_command_line_flag("--max_adapt",&max_adapt);
884 
885  // Parse command line
886  CommandLineArgs::parse_and_assign();
887 
888  // Doc what has actually been specified on the command line
889  CommandLineArgs::doc_specified_flags();
890 
891  //Set up the problem
893  RefineableQHelmholtzElement<2,3> > problem;
894 
895 
896  // Initial values for parameter values
899 
900  //Parameter incrementation
901  for(unsigned i=0;i<nstep;i++)
902  {
903  // Solve the problem with Newton's method, allowing
904  // up to max_adapt mesh adaptations after every solve.
905  problem.newton_solve(max_adapt);
906 
907  // Doc solution
908  problem.doc_solution();
909 
910  // Increment FSI parameter
911  Global_Parameters::Q+=q_increment;
913  }
914 
915 } //end of main
916 
917 
918 
919 
920 
921 
922 
923 
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:
std::complex< double > axisym_coefficient()
Coefficient in front of Hankel function for axisymmetric solution of Helmholtz potential.
std::complex< double > HankelH1(const double &k, const double &x)
Interface to Hankel function in maple style.
double Outer_radius
Radius of outer boundary of Helmholtz domain.
Definition: acoustic_fsi.cc:59
void create_helmholtz_DtN_elements()
Create DtN face elements.
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
void doc_solution()
Doc the solution.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of traction elements.
double K_squared
Square of wavenumber for the Helmholtz equation.
Definition: acoustic_fsi.cc:56
Coated disk FSI.
int main(int argc, char **argv)
Driver for acoustic fsi problem.
void create_helmholtz_fsi_flux_elements()
Create Helmholtz FSI flux elements.
void setup_interaction()
Setup interaction.
double exact_axisym_radiated_power()
Exact radiated power for axisymmetric solution.
string Directory
Output directory.
TreeBasedRefineableMeshBase * Helmholtz_mesh_pt
Pointer to Helmholtz mesh.
unsigned El_multiplier
Multiplier for number of elements.
double H_coating
Non-dim thickness of elastic coating.
Definition: acoustic_fsi.cc:65
HelmholtzDtNMesh< HELMHOLTZ_ELEMENT > * Helmholtz_outer_boundary_mesh_pt
Pointer to mesh containing the DtN elements.
TimeHarmonicIsotropicElasticityTensor E(Nu)
The elasticity tensor for the solid.
TreeBasedRefineableMeshBase * Solid_mesh_pt
Pointer to solid mesh.
void delete_face_elements(Mesh *const &boundary_mesh_pt)
Delete (face) elements in specified mesh.
ofstream Trace_file
Trace file.
double Nu
Poisson's ratio.
Definition: acoustic_fsi.cc:68
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of traction elements.
void solid_boundary_displacement(const Vector< double > &x, Vector< std::complex< double > > &u)
Displacement field on inner boundary of solid.
Definition: acoustic_fsi.cc:90
void actions_before_newton_solve()
Update function (empty)
void actions_after_newton_solve()
Update function (empty)
Mesh * FSI_traction_mesh_pt
Pointer to mesh of FSI traction elements.
unsigned N
Azimuthal wavenumber for imposed displacement of coating on inner boundary.
Definition: acoustic_fsi.cc:87
Mesh * Helmholtz_fsi_flux_mesh_pt
Pointer to mesh of Helmholtz FSI flux elements.
void create_fsi_traction_elements()
Create FSI traction elements.
void exact_axisym_potential(const Vector< double > &x, Vector< double > &soln)
Exact solution for Helmholtz potential for axisymmetric solution.
DocInfo Doc_info
DocInfo object for output.
void actions_before_newton_convergence_check()
Recompute gamma integral before checking Newton residuals.