The spatially-adaptive solution of azimuthally Fourier-decomposed time-harmonic 3D acoustic fluid-structure interaction problems on unstructured meshes

In this document we discuss the spatially-adaptive solution of 3D time-harmonic acoustic fluid-structure interaction problems in cylindrical polar coordinates on unstructured meshes.

The driver code is very similar to the one presented in another tutorial and we only discuss the changes necessary to deal with the generation of the adaptive, unstructured meshes and the assignment of different material properties to different parts of the solid domain.

# A test problem

The sketch below shows the problem setup: An elastic sphere which is reinforced with an azimuthal T-rib is immersed in a compressible fluid and subjected to a time-periodic pressure load of magnitude

(where is the zenith angle) along its inner surface. The parameter controls the "sharpness" of the pressure load. For we obtain a uniform, spherically symmetric load; the sketch below shows the pressure distribution (red vectors indicating the traction) for

Sketch of the problem setup.

# Results

The figure below shows an animation of the structure's time-harmonic oscillation. The blue shaded region shows the shape of the oscillating structure while the pink region shows its undeformed configuration. The left half of the plot is used to show the (mirror image of the) adaptive unstructured mesh on which the displacement field was computed:

Animation of the time-harmonic deformation.

Here is a plot of the corresponding fluid displacement potential, a measure of the fluid pressure:

The fluid displacement potential, a measure of the fluid pressure. Elevation: real part; contours: imaginary part.

This looks very pretty and shows that we can solve acoustic FSI problems in non-trivial geometries but should you believe the results? Here's an attempt to convince you: If we make the rib much softer than the sphere and set its inertia to zero, the rib will not offer much structural resistance and the sphere will deform as if the rib was not present. If we then set we apply a spherically symmetric forcing onto the structure and would expect the resulting displacement field (at least in the sphere) to be spherically symmetric, too.

The animation of the displacement field for this case, shown below, shows that this is indeed the case:

Animation of the time-harmonic deformation for a uniform pressure load and a very soft and inertia-less rib.

Here is a plot of the corresponding fluid displacement potential, a measure of the fluid pressure:

The fluid displacement potential, a measure of the fluid pressure for a uniform pressure load and a very soft and inertia-less rib. Elevation: real part; contours: imaginary part.

# The numerical solution

The driver code for this problem is very similar to the one discussed in another tutorial. Running sdiff on the two driver codes

demo_drivers/interaction/fourier_decomposed_acoustic_fsi/fourier_decomposed_acoustic_fsi.cc

and

demo_drivers/interaction/fourier_decomposed_acoustic_fsi/unstructured_fourier_decomposed_acoustic_fsi.cc

shows you the differences, the most important of which are:

• The provision of multiple non-dimensional Young's moduli and frequency parameters for the two different regions (the rib and the sphere). Recall that the non-dimensional Young's modulus specified via TimeHarmonicFourierDecomposedLinearElasticityEquations::youngs_modulus_pt() represents the ratio of the material's actual Young's modulus to the Young's modulus used in the non-dimensionalisation of the equations.

• The change of forcing from a prescribed time-harmonic displacement to a pressure load on the inside boundary – this requires yet another mesh of FaceElements.

• The provision of the actions_before/after_adapt() functions and a helper function complete_problem_setup() which rebuilds the elements (by passing the problem parameters to the elements) following the unstructured mesh adaptation. (The need/rationale for such a function is discussed in another tutorial.)

• The mesh generation – the specification of the curvilinear boundaries and the geometry of the rib is somewhat tedious. We refer to another tutorial for a discussion of how to define the internal mesh boundary that separates the two regions (the rib and the sphere) so that we can assign different material properties to them.

.
All of this is reasonably straightforward and provides a powerful code that automatically adapts both meshes while respecting the curvilinear boundaries of the domain. Have a look through the driver code and play with it.

# Code listing

Here's a listing of the complete driver code:

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented,
//LIC// multi-physics finite-element library, available
//LIC// at http://www.oomph-lib.org.
//LIC//
//LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
//LIC//
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC//
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
//LIC// Lesser General Public License for more details.
//LIC//
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301 USA.
//LIC//
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC//
//LIC//====================================================================
// Driver for Helmholtz/TimeHarmonicTimeHarmonicLinElast coupling
#include <complex>
#include <cmath>
//Oomph-lib includes
#include "generic.h"
//The Helmholtz equation
#include "fourier_decomposed_helmholtz.h"
//The Elasticity equation
#include "time_harmonic_fourier_decomposed_linear_elasticity.h"
// The interaction elements
#include "multi_physics.h"
// The meshes
#include "meshes/annular_mesh.h"
#include "meshes/triangle_mesh.h"
// Get the Bessel functions
#include "oomph_crbond_bessel.h"
using namespace oomph;
using namespace std;
/// ////////////////////////////////////////////////////////////////////
/// ////////////////////////////////////////////////////////////////////
// Straight line as geometric object
/// ////////////////////////////////////////////////////////////////////
/// ////////////////////////////////////////////////////////////////////
//=========================================================================
/// Straight 1D line in 2D space
//=========================================================================
class MyStraightLine : public GeomObject
{
public:
/// Constructor: Pass start and end points
MyStraightLine(const Vector<double>& r_start, const Vector<double>& r_end)
: GeomObject(1,2), R_start(r_start), R_end(r_end)
{ }
/// Broken copy constructor
{
BrokenCopy::broken_copy("MyStraightLine");
}
/// Destructor: Empty
/// Position Vector at Lagrangian coordinate zeta
void position(const Vector<double>& zeta, Vector<double>& r) const
{
// Position Vector
r[0] = R_start[0]+(R_end[0]-R_start[0])*zeta[0];
r[1] = R_start[1]+(R_end[1]-R_start[1])*zeta[0];
}
private:
/// Start point of line
Vector<double> R_start;
/// End point of line
Vector<double> R_end;
};
/// ///////////////////////////////////////////////////////////////////
/// ///////////////////////////////////////////////////////////////////
/// ///////////////////////////////////////////////////////////////////
//=======start_of_namespace==========================================
/// Global variables
//===================================================================
{
/// Square of wavenumber for the Helmholtz equation
double K_squared=10.0;
/// Radius of outer boundary of Helmholtz domain
/// FSI parameter
double Q=10.0;
/// Non-dim thickness of elastic coating
double H_coating=0.1;
/// Define azimuthal Fourier wavenumber
/// Poisson's ratio Nu
std::complex<double> Nu(std::complex<double>(0.3,0.0));
/// Define the non-dimensional Young's modulus
Vector<std::complex<double> > E(2,std::complex<double>(1.0,0.0));
/// Non-dim square of frequency for solid -- dependent variable!
Vector<std::complex<double> > Omega_sq(2,std::complex<double>(100.0,0.0));
/// Density ratio: solid to fluid
Vector<double> Density_ratio(2,1.0);
/// Function to update dependent parameter values
{
}
/// Uniform pressure
double P = 0.1;
/// Peakiness parameter for pressure load
double Alpha=0.0;
/// Pressure load (real and imag part)
const Vector<double> &n,
Vector<std::complex<double> >&traction)
{
double phi=atan2(x[1],x[0]);
double magnitude=exp(-Alpha*pow(phi-0.25*MathematicalConstants::Pi,2));
unsigned dim = 2;
for(unsigned i=0;i<dim;i++)
{
traction[i] = complex<double>(-magnitude*P*n[i],magnitude*P*n[i]);
}
/// Output directory
string Directory="RESLT";
} //end_of_namespace
//=============start_of_problem_class===================================
/// Coated sphere FSI
//======================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
class CoatedSphereProblem : public Problem
{
public:
/// Constructor:
/// Update function (empty)
/// Update function (empty)
/// Actions before adapt: Wipe the face meshes
/// Actions after adapt: Rebuild the face meshes
/// Recompute gamma integral before checking Newton residuals
{
Helmholtz_DtN_mesh_pt->setup_gamma();
}
/// Doc the solution
void doc_solution(DocInfo& doc_info);
private:
/// Create FSI traction elements
/// Create Helmholtz FSI flux elements
/// Setup interaction
/// Create DtN elements on outer boundary
/// Create solid traction elements
/// Delete (face) elements in specified mesh
void delete_face_elements(Mesh* const & boundary_mesh_pt);
// Complete problem setup
/// Boundary ID of upper symmetry boundary
/// Boundary ID of lower symmetry boundary
/// Boundary ID of upper inner boundary
/// Boundary ID of lower inner boundary
/// Boundary ID of outer boundary
/// Boundary ID of rib divider
/// Boundary ID of outer boundary in Helmholtz mesh
/// Boundary ID of inner boundary in Helmholtz mesh
/// Boundary ID of upper boundary in Helmholtz mesh
/// Boundary ID of lower boundary in Helmholtz mesh
/// Pointer to solid mesh
RefineableTriangleMesh<ELASTICITY_ELEMENT>* Solid_mesh_pt;
#else
/// Pointer to solid mesh
TriangleMesh<ELASTICITY_ELEMENT>* Solid_mesh_pt;
#endif
/// Pointer to mesh of solid traction elements
/// Pointer to mesh of FSI traction elements
/// Pointer to Helmholtz mesh
RefineableTriangleMesh<HELMHOLTZ_ELEMENT>* Helmholtz_mesh_pt;
#else
/// Pointer to Helmholtz mesh
TriangleMesh<HELMHOLTZ_ELEMENT>* Helmholtz_mesh_pt;
#endif
/// Pointer to mesh of Helmholtz FSI flux elements
/// Pointer to mesh containing the DtN elements
FourierDecomposedHelmholtzDtNMesh<HELMHOLTZ_ELEMENT>* Helmholtz_DtN_mesh_pt;
/// Trace file
ofstream Trace_file;
};// end_of_problem_class
//===========start_of_constructor=======================================
/// Constructor:
//======================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
{
// Solid mesh
//-----------
{
// Start and end coordinates
Vector<double> r_start(2);
Vector<double> r_end(2);
double r_outer = 1.0;
double r_inner = r_outer-Global_Parameters::H_coating;
// Thickness of rib
double rib_thick=0.05;
// Depth of rib
double rib_depth=0.2;
// Total width of T
double t_width=0.2;
// Thickness of T
double t_thick=0.05;
// Half-opening angle of rib
double half_phi_rib=asin(0.5*rib_thick/r_inner);
// Pointer to the closed curve that defines the outer boundary
TriangleMeshClosedCurve* closed_curve_pt=0;
// Provide storage for pointers to the parts of the curvilinear boundary
Vector<TriangleMeshCurveSection*> curvilinear_boundary_pt;
// Outer boundary
//---------------
Ellipse* outer_boundary_circle_pt = new Ellipse(r_outer,r_outer);
double zeta_start=-0.5*MathematicalConstants::Pi;
double zeta_end=0.5*MathematicalConstants::Pi;
unsigned nsegment=50;
unsigned boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
outer_boundary_circle_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Remember it
Outer_boundary_id=boundary_id;
// Upper straight line segment on symmetry axis
//---------------------------------------------
r_start[0]=0.0;
r_start[1]=r_outer;
r_end[0]=0.0;
r_end[1]=r_inner;
MyStraightLine* upper_sym_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
upper_sym_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Remember it
Upper_symmetry_boundary_id=boundary_id;
// Upper part of inner boundary
//-----------------------------
Ellipse* upper_inner_boundary_pt =
new Ellipse(r_inner,r_inner);
zeta_start=0.5*MathematicalConstants::Pi;
zeta_end=half_phi_rib;
nsegment=20;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
upper_inner_boundary_pt,
zeta_start,zeta_end,nsegment,boundary_id));
// Remember it
Upper_inner_boundary_id=boundary_id;
// Data associated with rib
MyStraightLine* upper_inward_rib_pt=0;
MyStraightLine* lower_inward_rib_pt=0;
TriangleMeshCurviLine* upper_inward_rib_curviline_pt=0;
Vector<TriangleMeshOpenCurve*> inner_boundary_pt;
TriangleMeshCurviLine* lower_inward_rib_curviline_pt=0;
Vector<double> rib_center(2);
// Upper half of inward rib
//-------------------------
r_start[0]=r_inner*cos(half_phi_rib);
r_start[1]=r_inner*sin(half_phi_rib);
r_end[0]=r_start[0]-rib_depth;
r_end[1]=r_start[1];
upper_inward_rib_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
upper_inward_rib_curviline_pt=
new TriangleMeshCurviLine(
upper_inward_rib_pt,zeta_start,zeta_end,nsegment,boundary_id);
curvilinear_boundary_pt.push_back(upper_inward_rib_curviline_pt);
// Vertical upper bit of T
//------------------------
r_start[0]=r_end[0];
r_start[1]=r_end[1];
r_end[0]=r_start[0];
r_end[1]=r_start[1]+0.5*(t_width-rib_thick);
MyStraightLine* vertical_upper_t_rib_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
vertical_upper_t_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Horizontal upper bit of T
//-----------==-------------
r_start[0]=r_end[0];
r_start[1]=r_end[1];
r_end[0]=r_start[0]-t_thick;
r_end[1]=r_start[1];
MyStraightLine* horizontal_upper_t_rib_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
horizontal_upper_t_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Vertical end of rib end
//------------------------
r_start[0]=r_end[0];
r_start[1]=r_end[1];
r_end[0]=r_start[0];
r_end[1]=-r_start[1];
MyStraightLine* inner_vertical_rib_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
inner_vertical_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Horizontal lower bit of T
//-----------==-------------
r_start[0]=r_end[0];
r_start[1]=r_end[1];
r_end[0]=r_start[0]+t_thick;
r_end[1]=r_start[1];
MyStraightLine* horizontal_lower_t_rib_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
horizontal_lower_t_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Vertical lower bit of T
//------------------------
r_start[0]=r_end[0];
r_start[1]=r_end[1];
r_end[0]=r_start[0];
r_end[1]=r_start[1]+0.5*(t_width-rib_thick);
MyStraightLine* vertical_lower_t_rib_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
vertical_lower_t_rib_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Lower half of inward rib
//-------------------------
r_end[0]=r_inner*cos(half_phi_rib);
r_end[1]=-r_inner*sin(half_phi_rib);
r_start[0]=r_end[0]-rib_depth;
r_start[1]=r_end[1];
lower_inward_rib_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
lower_inward_rib_curviline_pt=
new TriangleMeshCurviLine(
lower_inward_rib_pt,zeta_start,zeta_end,nsegment,boundary_id);
curvilinear_boundary_pt.push_back(lower_inward_rib_curviline_pt);
// Lower part of inner boundary
//-----------------------------
Ellipse* lower_inner_boundary_circle_pt = new Ellipse(r_inner,r_inner);
zeta_start=-half_phi_rib;
zeta_end=-0.5*MathematicalConstants::Pi;
nsegment=20;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
lower_inner_boundary_circle_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Remember it
Lower_inner_boundary_id=boundary_id;
// Lower straight line segment on symmetry axis
//---------------------------------------------
r_start[0]=0.0;
r_start[1]=-r_inner;
r_end[0]=0.0;
r_end[1]=-r_outer;
MyStraightLine* lower_sym_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
lower_sym_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Remember it
Lower_symmetry_boundary_id=boundary_id;
// Combine to curvilinear boundary
//--------------------------------
closed_curve_pt=
new TriangleMeshClosedCurve(curvilinear_boundary_pt);
// Vertical dividing line across base of T-rib
//--------------------------------------------
Vector<TriangleMeshCurveSection*> internal_polyline_pt(1);
r_start[0]=r_inner*cos(half_phi_rib);
r_start[1]=r_inner*sin(half_phi_rib);
r_end[0]=r_inner*cos(half_phi_rib);
r_end[1]=-r_inner*sin(half_phi_rib);
Vector<Vector<double> > boundary_vertices(2);
boundary_vertices[0]=r_start;
boundary_vertices[1]=r_end;
boundary_id=100;
TriangleMeshPolyLine* rib_divider_pt=
new TriangleMeshPolyLine(boundary_vertices,boundary_id);
internal_polyline_pt[0]=rib_divider_pt;
// Remember it
Rib_divider_boundary_id=boundary_id;
// Make connection
double s_connect=0.0;
internal_polyline_pt[0]->connect_initial_vertex_to_curviline(
upper_inward_rib_curviline_pt,s_connect);
// Make connection
s_connect=1.0;
internal_polyline_pt[0]->connect_final_vertex_to_curviline(
lower_inward_rib_curviline_pt,s_connect);
// Create open curve that defines internal bondary
inner_boundary_pt.push_back(new TriangleMeshOpenCurve(internal_polyline_pt));
// Define coordinates of a point inside the rib
rib_center[0]=r_inner-rib_depth;
rib_center[1]=0.0;
// Now build the mesh
//===================
// Use the TriangleMeshParameters object for helping on the manage of the
// TriangleMesh parameters. The only parameter that needs to take is the
// outer boundary.
TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
// Target area
triangle_mesh_parameters.element_area()=0.2;
// Specify the internal open boundary
triangle_mesh_parameters.internal_open_curves_pt()=inner_boundary_pt;
// Define the region
// Build the mesh
Solid_mesh_pt=new
RefineableTriangleMesh<ELASTICITY_ELEMENT>(triangle_mesh_parameters);
#else
// Build the mesh
Solid_mesh_pt=new
TriangleMesh<ELASTICITY_ELEMENT>(triangle_mesh_parameters);
#endif
}
// Helmholtz mesh
//---------------
{
// Start and end coordinates
Vector<double> r_start(2);
Vector<double> r_end(2);
// Inner radius of helmholtz region
double r_inner = 1.0;
// Outer radius of Helmholtz region
// Pointer to the closed curve that defines the outer boundary
TriangleMeshClosedCurve* closed_curve_pt=0;
// Provide storage for pointers to the parts of the curvilinear boundary
Vector<TriangleMeshCurveSection*> curvilinear_boundary_pt;
// Outer boundary
//---------------
Ellipse* outer_boundary_circle_pt = new Ellipse(r_outer,r_outer);
double zeta_start=-0.5*MathematicalConstants::Pi;
double zeta_end=0.5*MathematicalConstants::Pi;
unsigned nsegment=50;
unsigned boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
outer_boundary_circle_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Remember it
HH_outer_boundary_id=boundary_id;
// Upper straight line segment on symmetry axis
//---------------------------------------------
r_start[0]=0.0;
r_start[1]=r_outer;
r_end[0]=0.0;
r_end[1]=r_inner;
MyStraightLine* upper_sym_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
upper_sym_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Remember it
HH_upper_symmetry_boundary_id=boundary_id;
// Inner boundary
//---------------
Ellipse* upper_inner_boundary_pt =
new Ellipse(r_inner,r_inner);
zeta_start=0.5*MathematicalConstants::Pi;
zeta_end=-0.5*MathematicalConstants::Pi;
nsegment=40;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
upper_inner_boundary_pt,
zeta_start,zeta_end,nsegment,boundary_id));
// Remember it
HH_inner_boundary_id=boundary_id;
// Lower straight line segment on symmetry axis
//---------------------------------------------
r_start[0]=0.0;
r_start[1]=-r_inner;
r_end[0]=0.0;
r_end[1]=-r_outer;
MyStraightLine* lower_sym_pt = new MyStraightLine(r_start,r_end);
zeta_start=0.0;
zeta_end=1.0;
nsegment=1;
boundary_id=curvilinear_boundary_pt.size();
curvilinear_boundary_pt.push_back(
new TriangleMeshCurviLine(
lower_sym_pt,zeta_start,zeta_end,nsegment,boundary_id));
// Remember it
HH_lower_symmetry_boundary_id=boundary_id;
// Combine to curvilinear boundary
//--------------------------------
closed_curve_pt=
new TriangleMeshClosedCurve(curvilinear_boundary_pt);
// Now build the mesh
//===================
// Use the TriangleMeshParameters object for helping on the manage of the
// TriangleMesh parameters. The only parameter that needs to take is the
// outer boundary.
TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
// Target area
triangle_mesh_parameters.element_area()=0.2;
// Build the mesh
Helmholtz_mesh_pt=new
RefineableTriangleMesh<HELMHOLTZ_ELEMENT>(triangle_mesh_parameters);
#else
// Build the mesh
Helmholtz_mesh_pt=new
TriangleMesh<HELMHOLTZ_ELEMENT>(triangle_mesh_parameters);
#endif
}
// Create mesh for DtN elements on outer boundary
unsigned nfourier=20;
Helmholtz_DtN_mesh_pt=
new FourierDecomposedHelmholtzDtNMesh<HELMHOLTZ_ELEMENT>(
// Set error estimators
Solid_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
Helmholtz_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
#endif
// Output meshes and their boundaries so far so we can double
// check the boundary enumeration
Solid_mesh_pt->output("solid_mesh.dat");
Helmholtz_mesh_pt->output("helmholtz_mesh.dat");
Solid_mesh_pt->output_boundaries("solid_mesh_boundary.dat");
Helmholtz_mesh_pt->output_boundaries("helmholtz_mesh_boundary.dat");
// Create FaceElement meshes for boundary conditions
//--------------------------------------------------
// Construct the solid traction element mesh
Solid_traction_mesh_pt=new Mesh;
create_solid_traction_elements();
// Construct the fsi traction element mesh
FSI_traction_mesh_pt=new Mesh;
create_fsi_traction_elements();
// Construct the Helmholtz fsi flux element mesh
Helmholtz_fsi_flux_mesh_pt=new Mesh;
create_helmholtz_fsi_flux_elements();
// Create DtN elements
create_helmholtz_DtN_elements();
// Combine sub meshes
//-------------------
// Build the Problem's global mesh from its various sub-meshes
build_global_mesh();
// Complete problem setup
complete_problem_setup();
// Setup fluid-structure interaction
//----------------------------------
setup_interaction();
// Open trace file
char filename[100];
sprintf(filename,"%s/trace.dat",Global_Parameters::Directory.c_str());
Trace_file.open(filename);
// Setup equation numbering scheme
cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
}//end_of_constructor
/// Actions before adapt: Wipe the meshes face elements
//========================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
{
// Kill the solid traction elements and wipe surface mesh
delete_face_elements(Solid_traction_mesh_pt);
// Kill the fsi traction elements and wipe surface mesh
delete_face_elements(FSI_traction_mesh_pt);
// Kill Helmholtz FSI flux elements
delete_face_elements(Helmholtz_fsi_flux_mesh_pt);
// Kill Helmholtz BC elements
delete_face_elements(Helmholtz_DtN_mesh_pt);
// Rebuild the Problem's global mesh from its various sub-meshes
rebuild_global_mesh();
/// Actions after adapt: Rebuild the meshes of face elements
//========================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
{
// Complete problem setup
complete_problem_setup();
// Construct the solid traction elements
create_solid_traction_elements();
// Create fsi traction elements from all elements that are
create_fsi_traction_elements();
// Create Helmholtz fsi flux elements
create_helmholtz_fsi_flux_elements();
// Create DtN elements from all elements that are
// adjacent to the outer boundary of Helmholtz mesh
create_helmholtz_DtN_elements();
// Setup interaction
setup_interaction();
// Rebuild the Problem's global mesh from its various sub-meshes
rebuild_global_mesh();
/// Complete problem setup: Apply boundary conditions and set
/// physical properties
//========================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
{
// Solid boundary conditions:
//---------------------------
// Pin real and imag part of horizontal and azimuthal displacement components
//---------------------------------------------------------------------------
// on vertical boundaries
//-----------------------
{
//Loop over the nodes to pin and assign boundary displacements on
//solid boundary
unsigned n_node = Solid_mesh_pt->nboundary_node(Upper_symmetry_boundary_id);
for(unsigned i=0;i<n_node;i++)
{
Node* nod_pt=Solid_mesh_pt->boundary_node_pt(Upper_symmetry_boundary_id,i);
// Real part of r-displacement
nod_pt->pin(0);
nod_pt->set_value(0,0.0);
// Imag part of r-displacement
nod_pt->pin(3);
nod_pt->set_value(3,0.0);
// Real part of phi-displacement
nod_pt->pin(2);
nod_pt->set_value(2,0.0);
// Imag part of phi-displacement
nod_pt->pin(5);
nod_pt->set_value(5,0.0);
}
}
{
//Loop over the nodes to pin and assign boundary displacements on
//solid boundary
unsigned n_node = Solid_mesh_pt->nboundary_node(Lower_symmetry_boundary_id);
for(unsigned i=0;i<n_node;i++)
{
Node* nod_pt=Solid_mesh_pt->boundary_node_pt(Lower_symmetry_boundary_id,i);
// Real part of r-displacement
nod_pt->pin(0);
nod_pt->set_value(0,0.0);
// Imag part of r-displacement
nod_pt->pin(3);
nod_pt->set_value(3,0.0);
// Real part of phi-displacement
nod_pt->pin(2);
nod_pt->set_value(2,0.0);
// Imag part of phi-displacement
nod_pt->pin(5);
nod_pt->set_value(5,0.0);
}
}
//Assign the physical properties to the elements
//----------------------------------------------
unsigned nreg=Solid_mesh_pt->nregion();
for (unsigned r=0;r<nreg;r++)
{
unsigned nel=Solid_mesh_pt->nregion_element(r);
for (unsigned e=0;e<nel;e++)
{
//Cast to a solid element
ELASTICITY_ELEMENT *el_pt =
dynamic_cast<ELASTICITY_ELEMENT*>(Solid_mesh_pt->
region_element_pt(r,e));
// Set the pointer to Fourier wavenumber
el_pt->fourier_wavenumber_pt() = &Global_Parameters::Fourier_wavenumber;
// Set the pointer to Poisson's ratio
el_pt->nu_pt() = &Global_Parameters::Nu;
// Square of non-dim frequency
el_pt->omega_sq_pt()= &Global_Parameters::Omega_sq[r];
// Set the pointer to non-dim Young's modulus
el_pt->youngs_modulus_pt() = &Global_Parameters::E[r];
}
}
// Complete the build of all Helmholtz elements so they are fully functional
unsigned n_element = Helmholtz_mesh_pt->nelement();
for(unsigned i=0;i<n_element;i++)
{
// Upcast from GeneralsedElement to the present element
HELMHOLTZ_ELEMENT *el_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
Helmholtz_mesh_pt->element_pt(i));
//Set the k_squared pointer
el_pt->k_squared_pt()=&Global_Parameters::K_squared;
// Set pointer to Fourier wave number
el_pt->fourier_wavenumber_pt()=&Global_Parameters::Fourier_wavenumber;
}
}
//============start_of_delete_face_elements================
/// Delete face elements and wipe the mesh
//==========================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
delete_face_elements(Mesh* const & boundary_mesh_pt)
{
// How many surface elements are in the surface mesh
unsigned n_element = boundary_mesh_pt->nelement();
// Loop over the surface elements
for(unsigned e=0;e<n_element;e++)
{
// Kill surface element
delete boundary_mesh_pt->element_pt(e);
}
// Wipe the mesh
boundary_mesh_pt->flush_element_and_node_storage();
} // end of delete_face_elements
//============start_of_create_outer_bc_elements===========================
/// Create BC elements on outer boundary
//========================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
{
// Outer boundary:
unsigned b=HH_outer_boundary_id;
// Loop over the bulk elements adjacent to boundary b?
unsigned n_element = Helmholtz_mesh_pt->nboundary_element(b);
for(unsigned e=0;e<n_element;e++)
{
// Get pointer to the bulk element that is adjacent to boundary b
HELMHOLTZ_ELEMENT* bulk_elem_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
Helmholtz_mesh_pt->boundary_element_pt(b,e));
//Find the index of the face of element e along boundary b
int face_index = Helmholtz_mesh_pt->face_index_at_boundary(b,e);
// Build the corresponding DtN element
FourierDecomposedHelmholtzDtNBoundaryElement<HELMHOLTZ_ELEMENT>*
flux_element_pt = new
FourierDecomposedHelmholtzDtNBoundaryElement<HELMHOLTZ_ELEMENT>
(bulk_elem_pt,face_index);
//Add the flux boundary element to the helmholtz_DtN_mesh
// Set pointer to the mesh that contains all the boundary condition
// elements on this boundary
flux_element_pt->set_outer_boundary_mesh_pt(Helmholtz_DtN_mesh_pt);
}
} // end_of_create_outer_bc_elements
//=====================start_of_setup_interaction======================
/// Setup interaction between two fields
//========================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
{
// Setup Helmholtz "pressure" load on traction elements
unsigned boundary_in_helmholtz_mesh=HH_inner_boundary_id;
// Doc boundary coordinate for Helmholtz
ofstream the_file;
the_file.open("boundary_coordinate_hh.dat");
Helmholtz_mesh_pt->Mesh::template doc_boundary_coordinates<HELMHOLTZ_ELEMENT>
(boundary_in_helmholtz_mesh, the_file);
the_file.close();
// Setup interaction
<HELMHOLTZ_ELEMENT,2>
(this,boundary_in_helmholtz_mesh,Helmholtz_mesh_pt,FSI_traction_mesh_pt);
// Setup Helmholtz flux from normal displacement interaction
unsigned boundary_in_solid_mesh=Outer_boundary_id;
// Doc boundary coordinate for solid mesh
the_file.open("boundary_coordinate_solid.dat");
Solid_mesh_pt->Mesh::template doc_boundary_coordinates<ELASTICITY_ELEMENT>
(boundary_in_solid_mesh, the_file);
the_file.close();
// Setup interaction
<ELASTICITY_ELEMENT,2>(
this,boundary_in_solid_mesh,Solid_mesh_pt,Helmholtz_fsi_flux_mesh_pt);
}// end_of_setup_interaction
//============start_of_create_fsi_traction_elements======================
/// Create fsi traction elements
//=======================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
{
// We're on outer boundary of the solid mesh
unsigned b=Outer_boundary_id;
// How many bulk elements are adjacent to boundary b?
unsigned n_element = Solid_mesh_pt->nboundary_element(b);
// Loop over the bulk elements adjacent to boundary b
for(unsigned e=0;e<n_element;e++)
{
// Get pointer to the bulk element that is adjacent to boundary b
ELASTICITY_ELEMENT* bulk_elem_pt = dynamic_cast<ELASTICITY_ELEMENT*>(
Solid_mesh_pt->boundary_element_pt(b,e));
//Find the index of the face of element e along boundary b
int face_index = Solid_mesh_pt->face_index_at_boundary(b,e);
// Create element
<ELASTICITY_ELEMENT,HELMHOLTZ_ELEMENT>* el_pt=
<ELASTICITY_ELEMENT,HELMHOLTZ_ELEMENT>(bulk_elem_pt,
face_index);
// Associate element with bulk boundary (to allow it to access
// the boundary coordinates in the bulk mesh)
el_pt->set_boundary_number_in_bulk_mesh(b);
// Set FSI parameter
el_pt->q_pt()=&Global_Parameters::Q;
}
} // end_of_create_fsi_traction_elements
//============start_of_create_helmholtz_fsi_flux_elements================
/// Create Helmholtz fsi flux elements
//=======================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
{
// Attach to inner boundary of Helmholtz mesh
unsigned b=HH_inner_boundary_id;
// How many bulk elements are adjacent to boundary b?
unsigned n_element = Helmholtz_mesh_pt->nboundary_element(b);
// Loop over the bulk elements adjacent to boundary b
for(unsigned e=0;e<n_element;e++)
{
// Get pointer to the bulk element that is adjacent to boundary b
HELMHOLTZ_ELEMENT* bulk_elem_pt = dynamic_cast<HELMHOLTZ_ELEMENT*>(
Helmholtz_mesh_pt->boundary_element_pt(b,e));
//Find the index of the face of element e along boundary b
int face_index = Helmholtz_mesh_pt->face_index_at_boundary(b,e);
// Create element
FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement
<HELMHOLTZ_ELEMENT,ELASTICITY_ELEMENT>* el_pt=
new FourierDecomposedHelmholtzFluxFromNormalDisplacementBCElement
<HELMHOLTZ_ELEMENT,ELASTICITY_ELEMENT>(bulk_elem_pt,
face_index);
// Associate element with bulk boundary (to allow it to access
// the boundary coordinates in the bulk mesh)
el_pt->set_boundary_number_in_bulk_mesh(b);
}
} // end_of_create_helmholtz_fsi_flux_elements
//============start_of_create_solid_traction_elements====================
/// Create solid traction elements
//=======================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
{
// Loop over pressure loaded boundaries
unsigned b=0;
unsigned nb=3;
for (unsigned i=0;i<nb;i++)
{
switch(i)
{
case 0:
b=Upper_inner_boundary_id;
break;
case 1:
b=Lower_inner_boundary_id;
break;
case 2:
b=Rib_divider_boundary_id;
break;
}
// We're attaching face elements to region 0
unsigned r=0;
// How many bulk elements are adjacent to boundary b?
unsigned n_element = Solid_mesh_pt->nboundary_element_in_region(b,r);
// Loop over the bulk elements adjacent to boundary b
for(unsigned e=0;e<n_element;e++)
{
// Get pointer to the bulk element that is adjacent to boundary b
ELASTICITY_ELEMENT* bulk_elem_pt = dynamic_cast<ELASTICITY_ELEMENT*>(
Solid_mesh_pt->boundary_element_in_region_pt(b,r,e));
//Find the index of the face of element e along boundary b
int face_index = Solid_mesh_pt->face_index_at_boundary_in_region(b,r,e);
// Create element
TimeHarmonicFourierDecomposedLinearElasticityTractionElement
<ELASTICITY_ELEMENT>* el_pt=
new TimeHarmonicFourierDecomposedLinearElasticityTractionElement
<ELASTICITY_ELEMENT>(bulk_elem_pt,face_index);
// Associate element with bulk boundary (to allow it to access
// the boundary coordinates in the bulk mesh)
el_pt->set_boundary_number_in_bulk_mesh(b);
//Set the traction function
}
}
} // end of create_traction_elements
//==============start_of_doc_solution===============================
/// Doc the solution
//==================================================================
template<class ELASTICITY_ELEMENT, class HELMHOLTZ_ELEMENT>
doc_solution(DocInfo& doc_info)
{
// Doc parameters
oomph_info << "Writing result for step " << doc_info.number()
<< ". Parameters: "<< std::endl;
oomph_info << "Fourier mode number : N = "
oomph_info << "FSI parameter : Q = " << Global_Parameters::Q << std::endl;
<< std::endl;
oomph_info << "Fluid wavenumber : k^2 = " << Global_Parameters::K_squared
<< std::endl;
oomph_info << "Solid wavenumber : Omega_sq = "
<< Global_Parameters::Omega_sq[0] << std::endl;
oomph_info << "Solid wavenumber : Omega_sq = "
<< std::endl << std::endl;
ofstream some_file,some_file2;
char filename[100];
// Number of plot points
unsigned n_plot=5;
//----------------------------------
sprintf(filename,"%s/power%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
// Accumulate contribution from elements
double power=0.0;
unsigned nn_element=Helmholtz_DtN_mesh_pt->nelement();
for(unsigned e=0;e<nn_element;e++)
{
FourierDecomposedHelmholtzBCElementBase<HELMHOLTZ_ELEMENT> *el_pt =
dynamic_cast<FourierDecomposedHelmholtzBCElementBase<HELMHOLTZ_ELEMENT>*>(
Helmholtz_DtN_mesh_pt->element_pt(e));
power += el_pt->global_power_contribution(some_file);
}
some_file.close();
oomph_info << "Radiated power: " << power << std::endl;
// Output displacement field
//--------------------------
sprintf(filename,"%s/elast_soln%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
Solid_mesh_pt->output(some_file,n_plot);
some_file.close();
// Output Helmholtz
//-----------------
sprintf(filename,"%s/helmholtz_soln%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
Helmholtz_mesh_pt->output(some_file,n_plot);
some_file.close();
// Output fsi traction elements
//-----------------------------
sprintf(filename,"%s/fsi_traction_soln%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
FSI_traction_mesh_pt->output(some_file,n_plot);
some_file.close();
// Output Helmholtz fsi flux elements
//-----------------------------------
sprintf(filename,"%s/fsi_flux_bc_soln%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
Helmholtz_fsi_flux_mesh_pt->output(some_file,n_plot);
some_file.close();
// Write trace file
Trace_file << Global_Parameters::Q << " "
<< Global_Parameters::Omega_sq[0].real() << " "
<< Global_Parameters::Omega_sq[1].real() << " "
<< power << " "
<< std::endl;
// Bump up counter
doc_info.number()++;
} //end_of_doc_solution
//=======start_of_main==================================================
//======================================================================
int main(int argc, char **argv)
{
// Store command line arguments
CommandLineArgs::setup(argc,argv);
// Define possible command line arguments and parse the ones that
// were actually specified
// Output directory
CommandLineArgs::specify_command_line_flag("--dir",
// Parameter controlling the sharpness of the pressure load
CommandLineArgs::specify_command_line_flag("--alpha",
// Parse command line
CommandLineArgs::parse_and_assign();
// Doc what has actually been specified on the command line
CommandLineArgs::doc_specified_flags();
// Set values for parameter values
// Update dependent parameters values
// Set up doc info
DocInfo doc_info;
// Set output directory
doc_info.set_directory(Global_Parameters::Directory);
// Set up the problem
ProjectableTimeHarmonicFourierDecomposedLinearElasticityElement
<TTimeHarmonicFourierDecomposedLinearElasticityElement<3> >,
ProjectableFourierDecomposedHelmholtzElement
<TFourierDecomposedHelmholtzElement<3> > > problem;
#else
// Set up the problem
TFourierDecomposedHelmholtzElement<3> > problem;
#endif
//Parameter incrementation
unsigned nstep=3;
for(unsigned i=0;i<nstep;i++)
{
// Solve the problem with Newton's method, allowing
#else
// Solve the problem with Newton's method
problem.newton_solve();
#endif
// Doc solution
problem.doc_solution(doc_info);
// Make rib a lot heavier but keep its stiffness
if (i==0)
{
}
// Make rib very soft and inertia-less
if (i==1)
{
}
}
} //end_of_main
unsigned HH_upper_symmetry_boundary_id
Boundary ID of upper boundary in Helmholtz mesh.
Actions before adapt: Wipe the face meshes.
void delete_face_elements(Mesh *const &boundary_mesh_pt)
Delete (face) elements in specified mesh.
void create_fsi_traction_elements()
Create FSI traction elements.
unsigned HH_inner_boundary_id
Boundary ID of inner boundary in Helmholtz mesh.
void actions_before_newton_convergence_check()
Recompute gamma integral before checking Newton residuals.
Mesh * FSI_traction_mesh_pt
Pointer to mesh of FSI traction elements.
TwoDAnnularMesh< ELASTICITY_ELEMENT > * Solid_mesh_pt
Pointer to solid mesh.
unsigned Rib_divider_boundary_id
Boundary ID of rib divider.
TwoDAnnularMesh< HELMHOLTZ_ELEMENT > * Helmholtz_mesh_pt
Pointer to Helmholtz mesh.
void complete_problem_setup()
Complete problem setup: Apply boundary conditions and set physical properties.
unsigned HH_outer_boundary_id
Boundary ID of outer boundary in Helmholtz mesh.
Actions after adapt: Rebuild the face meshes.
void create_helmholtz_DtN_elements()
Create DtN elements on outer boundary.
unsigned HH_lower_symmetry_boundary_id
Boundary ID of lower boundary in Helmholtz mesh.
void create_solid_traction_elements()
Create solid traction elements.
void create_helmholtz_fsi_flux_elements()
Create Helmholtz FSI flux elements.
unsigned Lower_symmetry_boundary_id
Boundary ID of lower symmetry boundary.
void actions_before_newton_solve()
Update function (empty)
Mesh * Solid_traction_mesh_pt
Pointer to mesh of solid traction elements.
void doc_solution(DocInfo &doc_info)
Doc the solution.
FourierDecomposedHelmholtzDtNMesh< HELMHOLTZ_ELEMENT > * Helmholtz_DtN_mesh_pt
Pointer to mesh containing the DtN elements.
unsigned Upper_inner_boundary_id
Boundary ID of upper inner boundary.
unsigned Outer_boundary_id
Boundary ID of outer boundary.
unsigned Lower_inner_boundary_id
Boundary ID of lower inner boundary.
unsigned Upper_symmetry_boundary_id
Boundary ID of upper symmetry boundary.
void actions_after_newton_solve()
Update function (empty)
void setup_interaction()
Setup interaction.
Mesh * Helmholtz_fsi_flux_mesh_pt
Pointer to mesh of Helmholtz FSI flux elements.
////////////////////////////////////////////////////////////////////
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()
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.
int main(int argc, char **argv)
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
void pressure_load(const Vector< double > &x, const Vector< double > &n, Vector< std::complex< double > > &traction)
Pressure load (real and imag part)
string Directory
Output directory.
double P
Uniform pressure.
std::complex< double > Nu(std::complex< double >(0.3, 0.0))
Poisson's ratio Nu.
double Density_ratio
Density ratio: solid to fluid.
std::complex< double > Omega_sq(std::complex< double >(100.0, 0.0))
Non-dim square of frequency for solid – dependent variable!
Vector< std::complex< double > > E(2, std::complex< double >(1.0, 0.0))
Define the non-dimensional Young's modulus.
Radius of outer boundary of Helmholtz domain.
double K_squared
Square of wavenumber for the Helmholtz equation.
void update_parameter_values()
Function to update dependent parameter values.
int Fourier_wavenumber
Define azimuthal Fourier wavenumber.
double H_coating
Non-dim thickness of elastic coating.
double Alpha