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-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 meshes
35#include "meshes/annular_mesh.h"
36#include "meshes/triangle_mesh.h"
37
38using namespace std;
39using 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//=========================================================================
52class MyStraightLine : public GeomObject
53{
54
55public:
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
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
87private:
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//================================================================
107namespace 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 {
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//======================================================================
178template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
179class CoatedDiskProblem : public Problem
180{
181
182public:
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
201
202private:
203
204 // Complete problem setup
206
207 /// 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
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
233
234 /// Pointer to Helmholtz mesh
235 TreeBasedRefineableMeshBase* Helmholtz_mesh_pt;
236
237 /// Pointer to mesh of Helmholtz FSI flux elements
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//======================================================================
273template<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//========================================================================
706template<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//========================================================================
732template<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//========================================================================
768template<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//==========================================================
841template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
843delete_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//=======================================================================
865template<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//=======================================================================
929template<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//=======================================================================
975template<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//===========================================================================
1019template<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//========================================================================
1060template<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//==================================================================
1098template<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) ="
1130 << " density ratio (rib) ="
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//======================================================================
1256int 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
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(
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_fsi_flux_mesh_pt
Pointer to mesh of Helmholtz FSI flux 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.
ofstream Trace_file
Trace file.
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.
DocInfo Doc_info
DocInfo object for output.
unsigned Outer_boundary_id
Boundary ID of outer boundary.
void setup_interaction()
Setup interaction.
TreeBasedRefineableMeshBase * Helmholtz_mesh_pt
Pointer to Helmholtz mesh.
void create_helmholtz_ABC_elements()
Create ABC face elements.
unsigned Upper_inner_boundary_id
Boundary ID of upper inner boundary.
Mesh * FSI_traction_mesh_pt
Pointer to mesh of FSI traction elements.
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.
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.