fourier_decomposed_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 #include <complex>
32 #include <cmath>
33 
34 //Oomph-lib includes
35 #include "generic.h"
36 
37 //The Helmholtz equation
38 #include "fourier_decomposed_helmholtz.h"
39 
40 //The Elasticity equation
41 #include "time_harmonic_fourier_decomposed_linear_elasticity.h"
42 
43 // The interaction elements
44 #include "multi_physics.h"
45 
46 // The mesh
47 #include "meshes/annular_mesh.h"
48 
49 // Get the Bessel functions
50 #include "oomph_crbond_bessel.h"
51 
52 using namespace oomph;
53 using namespace std;
54 
55 
56 //////////////////////////////////////////////////////////////////////
57 //////////////////////////////////////////////////////////////////////
58 //////////////////////////////////////////////////////////////////////
59 
60 
61 //=======start_of_namespace==========================================
62 /// Global variables
63 //===================================================================
64 namespace Global_Parameters
65 {
66 
67  /// \short Square of wavenumber for the Helmholtz equation
68  double K_squared=10.0;
69 
70  /// \short Radius of outer boundary of Helmholtz domain
71  double Outer_radius=2.0;
72 
73  /// FSI parameter
74  double Q=10.0;
75 
76  /// Non-dim thickness of elastic coating
77  double H_coating=0.2;
78 
79  /// Define azimuthal Fourier wavenumber
81 
82  /// Poisson's ratio Nu
83  std::complex<double> Nu(std::complex<double>(0.3,0.0));
84 
85  /// Non-dim square of frequency for solid -- dependent variable!
86  std::complex<double> Omega_sq(std::complex<double>(100.0,0.0));
87 
88  /// Density ratio: solid to fluid
89  double Density_ratio=1.0;
90 
91  /// Function to update dependent parameter values
93  {
95  }
96 
97  /// \short Wavenumber "zenith"-variation of imposed displacement of coating
98  /// on inner boundary
99  unsigned M=4;
100 
101  /// \short Displacement field on inner boundary of solid
102  void solid_boundary_displacement(const Vector<double>& x,
103  Vector<std::complex<double> >& u)
104  {
105  Vector<double> normal(2);
106  double norm=sqrt(x[0]*x[0]+x[1]*x[1]);
107  double theta=atan2(x[1],x[0]);
108  normal[0]=x[0]/norm;
109  normal[1]=x[1]/norm;
110 
111  u[0]=complex<double>(normal[0]*cos(double(M)*theta),0.0);
112  u[1]=complex<double>(normal[1]*cos(double(M)*theta),0.0);
113  }
114 
115 
116  /// Output directory
117  string Directory="RESLT";
118 
119  /// Multiplier for number of elements
120  unsigned El_multiplier=1;
121 
122 } //end_of_namespace
123 
124 
125 
126 //=============start_of_problem_class===================================
127 /// Coated sphere FSI
128 //======================================================================
129 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
130 class CoatedSphereProblem : public Problem
131 {
132 
133 public:
134 
135  /// Constructor:
137 
138  /// \short Update function (empty)
140 
141  /// Update function (empty)
143 
144  /// Recompute gamma integral before checking Newton residuals
146  {
147  Helmholtz_DtN_mesh_pt->setup_gamma();
148  }
149 
150  /// Doc the solution
151  void doc_solution(DocInfo& doc_info);
152 
153 private:
154 
155  /// \short Create FSI traction elements
156  void create_fsi_traction_elements();
157 
158  /// \short Create Helmholtz FSI flux elements
159  void create_helmholtz_fsi_flux_elements();
160 
161  /// Setup interaction
162  void setup_interaction();
163 
164  /// \short Create DtN elements on outer boundary
165  void create_helmholtz_DtN_elements();
166 
167  /// Pointer to solid mesh
168  TwoDAnnularMesh<ELASTICITY_ELEMENT>* Solid_mesh_pt;
169 
170  /// Pointer to mesh of FSI traction elements
172 
173  /// Pointer to Helmholtz mesh
174  TwoDAnnularMesh<HELMHOLTZ_ELEMENT>* Helmholtz_mesh_pt;
175 
176  /// Pointer to mesh of Helmholtz FSI flux elements
178 
179  /// \short Pointer to mesh containing the DtN elements
180  FourierDecomposedHelmholtzDtNMesh<HELMHOLTZ_ELEMENT>* Helmholtz_DtN_mesh_pt;
181 
182  /// Trace file
183  ofstream Trace_file;
184 
185 };// end_of_problem_class
186 
187 
188 //===========start_of_constructor=======================================
189 /// Constructor:
190 //======================================================================
191 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
194 {
195 
196  // Parameters for meshes
197  bool periodic=false;
198  double azimuthal_fraction_of_coating=0.5;
199  double phi=0.5*MathematicalConstants::Pi;
200 
201  // Solid mesh
202  //-----------
203  // Number of elements in azimuthal direction
204  unsigned ntheta_solid=10*Global_Parameters::El_multiplier;
205 
206  // Number of elements in radial direction
207  unsigned nr_solid=3*Global_Parameters::El_multiplier;
208 
209  // Innermost radius for solid mesh
210  double a=1.0-Global_Parameters::H_coating;
211 
212  // Build solid mesh
213  Solid_mesh_pt = new
214  TwoDAnnularMesh<ELASTICITY_ELEMENT>(periodic,azimuthal_fraction_of_coating,
215  ntheta_solid,nr_solid,a,
217 
218 
219  // Helmholtz mesh
220  //---------------
221 
222  // Number of elements in azimuthal direction in Helmholtz mesh
223  unsigned ntheta_helmholtz=11*Global_Parameters::El_multiplier;
224 
225  // Number of elements in radial direction in Helmholtz mesh
226  unsigned nr_helmholtz=3*Global_Parameters::El_multiplier;
227 
228  // Innermost radius of Helmholtz mesh
229  a=1.0;
230 
231  // Thickness of Helmholtz mesh
232  double h_thick_helmholtz=Global_Parameters::Outer_radius-a;
233 
234  // Build mesh
235  Helmholtz_mesh_pt = new TwoDAnnularMesh<HELMHOLTZ_ELEMENT>
236  (periodic,azimuthal_fraction_of_coating,
237  ntheta_helmholtz,nr_helmholtz,a,h_thick_helmholtz,phi);
238 
239 
240  // Create mesh for DtN elements on outer boundary
241  unsigned nfourier=20;
242  Helmholtz_DtN_mesh_pt=
243  new FourierDecomposedHelmholtzDtNMesh<HELMHOLTZ_ELEMENT>(
245 
246  // Complete the solid problem setup to make the elements fully functional
247  unsigned nel=Solid_mesh_pt->nelement();
248  for (unsigned e=0;e<nel;e++)
249  {
250  // Cast to a bulk element
251  ELASTICITY_ELEMENT* el_pt=dynamic_cast<ELASTICITY_ELEMENT*>(
252  Solid_mesh_pt->element_pt(e));
253 
254  // Set the pointer to Fourier wavenumber
255  el_pt->fourier_wavenumber_pt() = &Global_Parameters::Fourier_wavenumber;
256 
257  // Set the pointer to Poisson's ratio
258  el_pt->nu_pt() = &Global_Parameters::Nu;
259 
260  // Set the pointer to square of the angular frequency
261  el_pt->omega_sq_pt() = &Global_Parameters::Omega_sq;
262  }
263 
264  // Complete the build of all Helmholtz elements so they are fully functional
265  unsigned n_element = Helmholtz_mesh_pt->nelement();
266  for(unsigned i=0;i<n_element;i++)
267  {
268  // Upcast from GeneralsedElement to the present element
269  HELMHOLTZ_ELEMENT *el_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
270  Helmholtz_mesh_pt->element_pt(i));
271 
272  //Set the k_squared pointer
273  el_pt->k_squared_pt()=&Global_Parameters::K_squared;
274 
275  // Set pointer to Fourier wave number
276  el_pt->fourier_wavenumber_pt()=&Global_Parameters::Fourier_wavenumber;
277  }
278 
279  // Output meshes and their boundaries so far so we can double
280  // check the boundary enumeration
281  Solid_mesh_pt->output("solid_mesh.dat");
282  Helmholtz_mesh_pt->output("helmholtz_mesh.dat");
283  Solid_mesh_pt->output_boundaries("solid_mesh_boundary.dat");
284  Helmholtz_mesh_pt->output_boundaries("helmholtz_mesh_boundary.dat");
285 
286  // Create FaceElement meshes for boundary conditions
287  //--------------------------------------------------
288 
289  // Construct the fsi traction element mesh
290  FSI_traction_mesh_pt=new Mesh;
291  create_fsi_traction_elements();
292 
293  // Construct the Helmholtz fsi flux element mesh
294  Helmholtz_fsi_flux_mesh_pt=new Mesh;
295  create_helmholtz_fsi_flux_elements();
296 
297  // Create DtN elements
298  create_helmholtz_DtN_elements();
299 
300 
301  // Combine sub meshes
302  //-------------------
303  add_sub_mesh(Solid_mesh_pt);
304  add_sub_mesh(FSI_traction_mesh_pt);
305  add_sub_mesh(Helmholtz_mesh_pt);
306  add_sub_mesh(Helmholtz_fsi_flux_mesh_pt);
307  add_sub_mesh(Helmholtz_DtN_mesh_pt);
308 
309  // Build the Problem's global mesh from its various sub-meshes
310  build_global_mesh();
311 
312 
313  // Solid boundary conditions:
314  //---------------------------
315 
316  // Pin the solid inner boundary (boundary 0) in all directions
317  unsigned b=0;
318  unsigned n_node = Solid_mesh_pt->nboundary_node(b);
319 
320  Vector<std::complex<double> > u(2);
321  Vector<double> x(2);
322 
323  //Loop over the nodes to pin and assign boundary displacements on
324  //solid boundary
325  for(unsigned i=0;i<n_node;i++)
326  {
327  Node* nod_pt=Solid_mesh_pt->boundary_node_pt(b,i);
328  nod_pt->pin(0);
329  nod_pt->pin(1);
330  nod_pt->pin(2);
331  nod_pt->pin(3);
332  nod_pt->pin(4);
333  nod_pt->pin(5);
334 
335  // Assign prescribed displacements
336  x[0]=nod_pt->x(0);
337  x[1]=nod_pt->x(1);
339 
340  // Real part of radial displacement
341  nod_pt->set_value(0,u[0].real());
342  // Real part of axial displacement
343  nod_pt->set_value(1,u[1].real());
344  // Real part of azimuthal displacement
345  nod_pt->set_value(2,0.0);
346  // Imag part of radial displacement
347  nod_pt->set_value(3,u[0].imag());
348  // Imag part of axial displacement
349  nod_pt->set_value(4,u[1].imag());
350  // Imag part of azimuthal displacement
351  nod_pt->set_value(5,0.0);
352  }
353 
354  // Vertical Symmetry boundary (r=0 and z<0)
355  {
356  unsigned ibound=1;
357  {
358  unsigned num_nod= Solid_mesh_pt->nboundary_node(ibound);
359  for (unsigned inod=0;inod<num_nod;inod++)
360  {
361  // Get pointer to node
362  Node* nod_pt=Solid_mesh_pt->boundary_node_pt(ibound,inod);
363 
364  // Pin radial displacement (u_0 (real) and u_3 (imag))
365  nod_pt->pin(0);
366  nod_pt->set_value(0,0.0);
367  nod_pt->pin(3);
368  nod_pt->set_value(3,0.0);
369 
370  // Pin azimuthal displacement (u_2 (real) and u_5 (imag))
371  nod_pt->pin(2);
372  nod_pt->set_value(2,0.0);
373  nod_pt->pin(5);
374  nod_pt->set_value(5,0.0);
375  }
376  }
377  }
378 
379 
380  // Vertical Symmetry boundary (r=0 and z>0)
381  {
382  unsigned ibound=3;
383  {
384  unsigned num_nod= Solid_mesh_pt->nboundary_node(ibound);
385  for (unsigned inod=0;inod<num_nod;inod++)
386  {
387  // Get pointer to node
388  Node* nod_pt=Solid_mesh_pt->boundary_node_pt(ibound,inod);
389 
390  // Pin radial displacement (u_0 (real) and u_3 (imag))
391  nod_pt->pin(0);
392  nod_pt->set_value(0,0.0);
393  nod_pt->pin(3);
394  nod_pt->set_value(3,0.0);
395 
396  // Pin azimuthal displacement (u_2 (real) and u_5 (imag))
397  nod_pt->pin(2);
398  nod_pt->set_value(2,0.0);
399  nod_pt->pin(5);
400  nod_pt->set_value(5,0.0);
401  }
402  }
403  } // done sym bc
404 
405  // Setup fluid-structure interaction
406  //----------------------------------
407  setup_interaction();
408 
409  // Open trace file
410  char filename[100];
411  sprintf(filename,"%s/trace.dat",Global_Parameters::Directory.c_str());
412  Trace_file.open(filename);
413 
414  // Setup equation numbering scheme
415  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
416 
417 }//end_of_constructor
418 
419 
420 
421 //============start_of_create_outer_bc_elements===========================
422 /// Create BC elements on outer boundary
423 //========================================================================
424 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
427 {
428  // Outer boundary is boundary 2:
429  unsigned b=2;
430 
431  // Loop over the bulk elements adjacent to boundary b?
432  unsigned n_element = Helmholtz_mesh_pt->nboundary_element(b);
433  for(unsigned e=0;e<n_element;e++)
434  {
435  // Get pointer to the bulk element that is adjacent to boundary b
436  HELMHOLTZ_ELEMENT* bulk_elem_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
437  Helmholtz_mesh_pt->boundary_element_pt(b,e));
438 
439  //Find the index of the face of element e along boundary b
440  int face_index = Helmholtz_mesh_pt->face_index_at_boundary(b,e);
441 
442  // Build the corresponding DtN element
443  FourierDecomposedHelmholtzDtNBoundaryElement<HELMHOLTZ_ELEMENT>*
444  flux_element_pt = new
445  FourierDecomposedHelmholtzDtNBoundaryElement<HELMHOLTZ_ELEMENT>
446  (bulk_elem_pt,face_index);
447 
448  //Add the flux boundary element to the helmholtz_DtN_mesh
449  Helmholtz_DtN_mesh_pt->add_element_pt(flux_element_pt);
450 
451  // Set pointer to the mesh that contains all the boundary condition
452  // elements on this boundary
453  flux_element_pt->set_outer_boundary_mesh_pt(Helmholtz_DtN_mesh_pt);
454  }
455 
456 } // end_of_create_outer_bc_elements
457 
458 
459 
460 
461 
462 //=====================start_of_setup_interaction======================
463 /// Setup interaction between two fields
464 //========================================================================
465 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
468 {
469 
470  // Setup Helmholtz "pressure" load on traction elements
471  unsigned boundary_in_helmholtz_mesh=0;
472 
473  // Doc boundary coordinate for Helmholtz
474  ofstream the_file;
475  the_file.open("boundary_coordinate_hh.dat");
476  Helmholtz_mesh_pt->Mesh::template doc_boundary_coordinates<HELMHOLTZ_ELEMENT>
477  (boundary_in_helmholtz_mesh, the_file);
478  the_file.close();
479 
480  // Setup interaction
481  Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
482  <HELMHOLTZ_ELEMENT,2>
483  (this,boundary_in_helmholtz_mesh,Helmholtz_mesh_pt,FSI_traction_mesh_pt);
484 
485  // Setup Helmholtz flux from normal displacement interaction
486  unsigned boundary_in_solid_mesh=2;
487 
488  // Doc boundary coordinate for solid mesh
489  the_file.open("boundary_coordinate_solid.dat");
490  Solid_mesh_pt->Mesh::template doc_boundary_coordinates<ELASTICITY_ELEMENT>
491  (boundary_in_solid_mesh, the_file);
492  the_file.close();
493 
494  // Setup interaction
495  Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh
496  <ELASTICITY_ELEMENT,2>(
497  this,boundary_in_solid_mesh,Solid_mesh_pt,Helmholtz_fsi_flux_mesh_pt);
498 
499 }// end_of_setup_interaction
500 
501 
502 
503 
504 
505 //============start_of_create_fsi_traction_elements======================
506 /// Create fsi traction elements
507 //=======================================================================
508 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
511 {
512  // We're on boundary 2 of the solid mesh
513  unsigned b=2;
514 
515  // How many bulk elements are adjacent to boundary b?
516  unsigned n_element = Solid_mesh_pt->nboundary_element(b);
517 
518  // Loop over the bulk elements adjacent to boundary b
519  for(unsigned e=0;e<n_element;e++)
520  {
521  // Get pointer to the bulk element that is adjacent to boundary b
522  ELASTICITY_ELEMENT* bulk_elem_pt = dynamic_cast<ELASTICITY_ELEMENT*>(
523  Solid_mesh_pt->boundary_element_pt(b,e));
524 
525  //Find the index of the face of element e along boundary b
526  int face_index = Solid_mesh_pt->face_index_at_boundary(b,e);
527 
528  // Create element
529  FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement
530  <ELASTICITY_ELEMENT,HELMHOLTZ_ELEMENT>* el_pt=
531  new FourierDecomposedTimeHarmonicLinElastLoadedByHelmholtzPressureBCElement
532  <ELASTICITY_ELEMENT,HELMHOLTZ_ELEMENT>(bulk_elem_pt,
533  face_index);
534  // Add to mesh
535  FSI_traction_mesh_pt->add_element_pt(el_pt);
536 
537  // Associate element with bulk boundary (to allow it to access
538  // the boundary coordinates in the bulk mesh)
539  el_pt->set_boundary_number_in_bulk_mesh(b);
540 
541  // Set FSI parameter
542  el_pt->q_pt()=&Global_Parameters::Q;
543  }
544 
545 } // end_of_create_fsi_traction_elements
546 
547 
548 
549 //============start_of_create_helmholtz_fsi_flux_elements================
550 /// Create Helmholtz fsi flux elements
551 //=======================================================================
552 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
555 {
556 
557  // Attach to inner boundary of Helmholtz mesh (0)
558  unsigned b=0;
559 
560  // How many bulk elements are adjacent to boundary b?
561  unsigned n_element = Helmholtz_mesh_pt->nboundary_element(b);
562 
563  // Loop over the bulk elements adjacent to boundary b
564  for(unsigned e=0;e<n_element;e++)
565  {
566  // Get pointer to the bulk element that is adjacent to boundary b
567  HELMHOLTZ_ELEMENT* bulk_elem_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
568  Helmholtz_mesh_pt->boundary_element_pt(b,e));
569 
570  //Find the index of the face of element e along boundary b
571  int face_index = Helmholtz_mesh_pt->face_index_at_boundary(b,e);
572 
573  // Create element
574  FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement
575  <HELMHOLTZ_ELEMENT,ELASTICITY_ELEMENT>* el_pt=
576  new FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement
577  <HELMHOLTZ_ELEMENT,ELASTICITY_ELEMENT>(bulk_elem_pt,
578  face_index);
579 
580  // Add to mesh
581  Helmholtz_fsi_flux_mesh_pt->add_element_pt(el_pt);
582 
583  // Associate element with bulk boundary (to allow it to access
584  // the boundary coordinates in the bulk mesh)
585  el_pt->set_boundary_number_in_bulk_mesh(b);
586  }
587 
588 } // end_of_create_helmholtz_fsi_flux_elements
589 
590 
591 
592 //==============start_of_doc_solution===============================
593 /// Doc the solution
594 //==================================================================
595 template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
597 doc_solution(DocInfo& doc_info)
598 {
599 
600  // Doc parameters
601  oomph_info << "Writing result for step " << doc_info.number()
602  << ". Parameters: "<< std::endl;
603  oomph_info << "Fourier mode number : N = "
605  oomph_info << "FSI parameter : Q = " << Global_Parameters::Q << std::endl;
606  oomph_info << "Fluid outer radius : R = " << Global_Parameters::Outer_radius
607  << std::endl;
608  oomph_info << "Fluid wavenumber : k^2 = " << Global_Parameters::K_squared
609  << std::endl;
610  oomph_info << "Solid wavenumber : Omega_sq = " << Global_Parameters::Omega_sq
611  << std::endl << std::endl;
612 
613 
614  ofstream some_file,some_file2;
615  char filename[100];
616 
617  // Number of plot points
618  unsigned n_plot=5;
619 
620  // Compute/output the radiated power
621  //----------------------------------
622  sprintf(filename,"%s/power%i.dat",doc_info.directory().c_str(),
623  doc_info.number());
624  some_file.open(filename);
625 
626  // Accumulate contribution from elements
627  double power=0.0;
628  unsigned nn_element=Helmholtz_DtN_mesh_pt->nelement();
629  for(unsigned e=0;e<nn_element;e++)
630  {
631  FourierDecomposedHelmholtzBCElementBase<HELMHOLTZ_ELEMENT> *el_pt =
632  dynamic_cast<FourierDecomposedHelmholtzBCElementBase<HELMHOLTZ_ELEMENT>*>(
633  Helmholtz_DtN_mesh_pt->element_pt(e));
634  power += el_pt->global_power_contribution(some_file);
635  }
636  some_file.close();
637  oomph_info << "Radiated power: " << power << std::endl;
638 
639  // Output displacement field
640  //--------------------------
641  sprintf(filename,"%s/elast_soln%i.dat",doc_info.directory().c_str(),
642  doc_info.number());
643  some_file.open(filename);
644  Solid_mesh_pt->output(some_file,n_plot);
645  some_file.close();
646 
647  // Output Helmholtz
648  //-----------------
649  sprintf(filename,"%s/helmholtz_soln%i.dat",doc_info.directory().c_str(),
650  doc_info.number());
651  some_file.open(filename);
652  Helmholtz_mesh_pt->output(some_file,n_plot);
653  some_file.close();
654 
655 
656  // Output fsi traction elements
657  //-----------------------------
658  sprintf(filename,"%s/fsi_traction_soln%i.dat",doc_info.directory().c_str(),
659  doc_info.number());
660  some_file.open(filename);
661  FSI_traction_mesh_pt->output(some_file,n_plot);
662  some_file.close();
663 
664 
665  // Output Helmholtz fsi flux elements
666  //-----------------------------------
667  sprintf(filename,"%s/fsi_flux_bc_soln%i.dat",doc_info.directory().c_str(),
668  doc_info.number());
669  some_file.open(filename);
670  Helmholtz_fsi_flux_mesh_pt->output(some_file,n_plot);
671  some_file.close();
672 
673  // Write trace file
674  Trace_file << Global_Parameters::Q << " "
677  << Global_Parameters::Omega_sq.real() << " "
678  << power << " "
679  << std::endl;
680 
681  // Bump up counter
682  doc_info.number()++;
683 
684 } //end_of_doc_solution
685 
686 
687 
688 //=======start_of_main==================================================
689 /// Driver for coated sphere loaded by lineared fluid loading
690 //======================================================================
691 int main(int argc, char **argv)
692 {
693 
694  // Store command line arguments
695  CommandLineArgs::setup(argc,argv);
696 
697  // Define possible command line arguments and parse the ones that
698  // were actually specified
699 
700  // Output directory
701  CommandLineArgs::specify_command_line_flag("--dir",
703 
704  // Parameter for the Helmholtz equation
705  CommandLineArgs::specify_command_line_flag("--k_squared",
707 
708  // Initial value of Q
709  CommandLineArgs::specify_command_line_flag("--q_initial",
711 
712  // Number of steps in parameter study
713  unsigned nstep=2;
714  CommandLineArgs::specify_command_line_flag("--nstep",&nstep);
715 
716  // Increment in FSI parameter in parameter study
717  double q_increment=5.0;
718  CommandLineArgs::specify_command_line_flag("--q_increment",&q_increment);
719 
720 
721  // Wavenumber "zenith"-variation of imposed displacement of coating
722  // on inner boundary
723  CommandLineArgs::specify_command_line_flag("--M",
725 
726  // Multiplier for number of elements
727  CommandLineArgs::specify_command_line_flag("--el_multiplier",
729 
730  // Parse command line
731  CommandLineArgs::parse_and_assign();
732 
733  // Doc what has actually been specified on the command line
734  CommandLineArgs::doc_specified_flags();
735 
736  // Update dependent parameters values
738 
739  // Set up doc info
740  DocInfo doc_info;
741 
742  // Set output directory
743  doc_info.set_directory(Global_Parameters::Directory);
744 
745  // Set up the problem
747  QFourierDecomposedHelmholtzElement<3> > problem;
748 
749  //Parameter incrementation
750  for(unsigned i=0;i<nstep;i++)
751  {
752  // Solve the problem with Newton's method
753  problem.newton_solve();
754 
755  // Doc solution
756  problem.doc_solution(doc_info);
757 
758  // Increment FSI parameter
759  Global_Parameters::Q+=q_increment;
761  }
762 
763 } //end_of_main
764 
765 
766 
767 
768 
769 
770 
771 
void create_fsi_traction_elements()
Create FSI traction elements.
double Outer_radius
Radius of outer boundary of Helmholtz domain.
FourierDecomposedHelmholtzDtNMesh< HELMHOLTZ_ELEMENT > * Helmholtz_DtN_mesh_pt
Pointer to mesh containing the DtN elements.
void update_parameter_values()
Function to update dependent parameter values.
int Fourier_wavenumber
Define azimuthal Fourier wavenumber.
double Density_ratio
Density ratio: solid to fluid.
void create_helmholtz_DtN_elements()
Create DtN elements on outer boundary.
void actions_before_newton_solve()
Update function (empty)
std::complex< double > Omega_sq(std::complex< double >(100.0, 0.0))
Non-dim square of frequency for solid – dependent variable!
TwoDAnnularMesh< HELMHOLTZ_ELEMENT > * Helmholtz_mesh_pt
Pointer to Helmholtz mesh.
Mesh * Helmholtz_fsi_flux_mesh_pt
Pointer to mesh of Helmholtz FSI flux elements.
void actions_before_newton_convergence_check()
Recompute gamma integral before checking Newton residuals.
std::complex< double > Nu(std::complex< double >(0.3, 0.0))
Poisson's ratio Nu.
double K_squared
Square of wavenumber for the Helmholtz equation.
int main(int argc, char **argv)
Driver for coated sphere loaded by lineared fluid loading.
void setup_interaction()
Setup interaction.
void actions_after_newton_solve()
Update function (empty)
string Directory
Output directory.
unsigned El_multiplier
Multiplier for number of elements.
double H_coating
Non-dim thickness of elastic coating.
TwoDAnnularMesh< ELASTICITY_ELEMENT > * Solid_mesh_pt
Pointer to solid mesh.
void solid_boundary_displacement(const Vector< double > &x, Vector< std::complex< double > > &u)
Displacement field on inner boundary of solid.
void create_helmholtz_fsi_flux_elements()
Create Helmholtz FSI flux elements.
unsigned M
Wavenumber "zenith"-variation of imposed displacement of coating on inner boundary.
Mesh * FSI_traction_mesh_pt
Pointer to mesh of FSI traction elements.
void doc_solution(DocInfo &doc_info)
Doc the solution.