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