Example problem: The Helmholtz equation – scattering problems

In this document we discuss the finite-element-based solution of the Helmholtz equation, an elliptic PDE that describes time-harmonic wave propagation problems. We start by reviewing the relevant theory and then present the solution of a simple model problem – the scattering of a planar wave from a circular cylinder.

 Acknowledgement: This tutorial and the associated driver codes were developed jointly with Tarak Kharrat (EnstaParisTech, Paris).

# Theory: The Helmholtz equation for time-harmonic scattering problems

The Helmholtz equation governs time-harmonic solutions of problems governed by the linear wave equation where is the wavespeed. Assuming that is time-harmonic, with frequency , we write the real function as where is complex-valued. This transforms (1) into the Helmholtz equation where is the wave number. Like other elliptic PDEs the Helmholtz equation admits Dirichlet, Neumann (flux) and Robin boundary conditions.

If the equation is solved in an infinite domain (e.g. in scattering problems) the solution must satisfy the so-called Sommerfeld radiation condition which in 2D has the form Mathematically, this conditions is required to ensure the uniqueness of the solution (and hence the well-posedness of the problem). In a physical context, such as a scattering problem, the condition ensures that scattering of an incoming wave only produces outgoing not incoming waves from infinity.

# Discretisation by finite elements

The discretisation of the Helmholtz equation itself only requires a trivial modification of oomph-lib's Poisson elements – we simply add the term to the residual. Since most practical applications of the Helmholtz equation involve complex-valued solutions, we provide separate storage for the real and imaginary parts of the solution – each Node therefore stores two unknowns values. By default, the real and imaginary parts are stored as values 0 and 1, respectively; see the section The enumeration of the unknowns for details.

The application of Dirichlet and Neumann boundary conditions is straightforward and follows the pattern employed for the solution of the Poisson equation:

• Dirichlet conditions are imposed by pinning the relevant nodal values and setting them to the appropriate prescribed values.
• Neumann (flux) boundary conditions are imposed via FaceElements (here the HelmholtzFluxElements). As usual we attach these to the faces of the "bulk" elements that are subject to the Neumann boundary conditions.

The imposition of the Sommerfeld radiation condition for problems in infinite domains is slightly more complicated. In the following discussion we will restrict ourselves to two dimensions and assume that the infinite domain is truncated at a circular artificial boundary of radius [This assumption is also made in the implementation of oomph-lib's FaceElements that allow the (approximate) imposition of the Sommerfeld radiation condition. The methodology can easily be modified to deal with other geometries but this has not been done yet – any volunteers?] All methods exploit the fact that the relevant solution of the Helmholtz equation can be written in polar coordinates as where the are suitable coefficients and is the -th-order Hankel function of the first kind.

## Approximate/absorbing boundary conditions (ABCs)

It is possible to derive approximate versions of the Sommerfeld radiation condition in which the normal derivative of the solution on the artificial boundary is related to its value and possibly its tangential derivatives. Such boundary conditions (sometimes referred to as approximate or absorbing boundary conditions – ABCs) are typically derived from asymptotic expansions of the solution at large distances from the origin and become more accurate the larger the radius of the artificial boundary is. Higher accuracy can therefore only be achieved by increasing the size of the computational domain, with an associated increase in computational cost.

oomph-lib provides an implementation of the following three boundary conditions (all taken from J. J. Shirron & I. Babuska's paper "A comparison of approximate boundary conditions and infinite element methods for exterior Helmholtz problems", Computer Methods in Applied Mechanics and Engineering 164 121-139 (1998), in which the authors compare the accuracy of these and many other approximate boundary conditions).

• Feng's first order ABC: (This is identical to the first-order Bayliss and Turkel boundary condition).
• Feng's second order ABC: • Feng's third order ABC: All three boundary conditions are implemented in the class HelmholtzAbsorbingBCElement. The order of the approximation can be set via the member function HelmholtzAbsorbingBCElement::abc_order(). All three boundary conditions are local (relating the function to its normal derivative) and do therefore not change the sparsity of the resulting finite element equations.

## The Dirichlet-to-Neumann mapping (DtN)

Using (4), it is easy to show (see, e.g., J. Jin "The Finite Element Method in Electromagnetics (second edition)", Wiley (2002) p. 501ff – but note that Jin assumes that the potential varies like rather than as assumed here) that the normal (radial) derivative, on the artificial boundary is given by where Equation (5) again provides a condition on the normal derivative of the solution along the artificial boundary and is implemented in the HelmholtzDtNBoundaryElement class. Since depends on the solution everywhere along the artificial boundary (see (6)), the application of the boundary condition (5) introduces a non-local coupling between all the degrees of freedom located on that boundary. This is handled by classifying the unknowns that affect but are not associated with the element's own nodes as external Data.

To facilitate the setup of the interaction between the HelmholtzDtNBoundaryElements, oomph-lib provides the class HelmholtzDtNMesh which provides storage for (the pointers to) the HelmholtzDtNBoundaryElements that discretise the artificial boundary. The member function HelmholtzDtNMesh::setup_gamma() pre-computes the values required for the imposition of equation (5). The radius of the artificial boundary and the (finite) number of (Fourier) terms used in the sum in (6) are specified as arguments to the constructor of the HelmholtzDtNMesh.

NOTE: Since depends on the solution, it must be recomputed whenever the unknowns are updated during the Newton iteration. This is best done by adding a call to HelmholtzDtNMesh::setup_gamma() to Problem::actions_before_newton_convergence_check(). [If Helmholtz's equation is solved in isolation (or within a coupled, but linear problem), Newton's method will converge in one iteration. In such cases the unnecessary recomputation of after the one-and-only Newton iteration can be suppressed by setting Problem::Problem_is_nonlinear to false.]

# A specific example: Scattering of an acoustic wave from a sound-hard obstacle

We will now demonstrate the methodology for a specific example: the scattering of sound waves in an acoustic medium of density and bulk modulus . Assuming that an incoming sound wave impacts a rigid, impermeable obstacle as shown in this sketch, Scattering of an incoming wave from a sound-hard obstacle -- the scatterer.

we wish to find the wave field that is scattered from the body.

For this purpose we denote the time-dependent displacement of the fluid particle in the acoustic medium by and introduce a displacement potential such that (As usual we employ asterisks to distinguish dimensional quantities from their non-dimensional equivalents, to be introduced below.) It is easy to show that satisfies the linear wave equation (1) with wave speed .

Since the surface of the scatterer is impenetrable, the normal displacement of the fluid has to vanish on and the boundary condition for the displacement potential becomes We non-dimensionalise all lengths and displacements on some problem-dependent lengthscale (e.g. the radius of the scatterer), non-dimensionalise the potential as and scale time on the period of the oscillation, The governing equation then becomes where the square of the wavenumber is given by Assuming that the incoming wave (already satisfying (8)) is described by a (known) non-dimensional displacement potential of the form we write the total potential as where represents the displacement potential associated with the scattered field which must satisfy (2). The boundary condition (7) then becomes a Neumann (flux) boundary condition for the scattered field, For the special case of the incoming wave being a planar wave, propagating along the x-axis, the incoming field can be written in polar coordinates as where is the Bessel function of the first kind of order . The exact solution for the scattering of such a wave from a circular disk is given by the series where we have chosen the disk's radius, , as the lengthscale by setting . In the above expression, denotes the Hankel function of the first kind of order and the prime denotes differentiation with respect to the function's argument.

A quantity that is of particular interest in wave propagation problems is the time-average of the power radiated by the scatterer, In the context of an acoustic wave, the total instantaneous power, radiated over a closed boundary is where the pressure is related to the displacement potential via The non-dimensional time-averaged radiated power can be expressed in terms of the complex potential as # Results

The figure below shows an animation of the displacement potential for scattering from a circular disk for a non-dimensional wavenumber of over one period of the oscillation. The simulation was performed in an annular computational domain, bounded by the outer surface the (unit) disk and an artificial outer boundary of non-dimensional radius The Sommerfeld radiation condition was imposed using the DtN mapping and the simulation was performed with spatial adaptivity (note the non-uniform refinement).

The "carpet plot" compares the exact (green) and computed (red) solutions for the displacement potential. The colours in the contour plot at the bottom of the figure provide an alternative visualisation of the magnitude of the scattered field. The displacement potential associated with the scattered wave, animated over one period of the oscillation.

# The numerical solution

## The global namespace

As usual, we define the problem parameters in a global namespace. The main physical parameter is the (square of the) wave number, . N_fourier is the number of (Fourier) terms to be used in evaluation of the series in equations (6) and (10). The remaining parameters determine how the Sommerfeld radiation condition is applied.

//===== start_of_namespace=============================================

The function get_exact_u returns the exact solution for the scattering problem. We will use this function for the validation of our results.

/// Exact solution for scattered field
/// (vector returns real and impaginary parts).
void get_exact_u(const Vector<double>& x, Vector<double>& u)
{
// Switch to polar coordinates
double r;
r=sqrt(x*x+x*x);
double theta;
theta=atan2(x,x);
// Argument for Bessel/Hankel functions
double rr=sqrt(K_squared)*r;
// Evaluate Bessel/Hankel functions
complex <double > u_ex(0.0,0.0);
Vector<double> jn(N_fourier+1), yn(N_fourier+1),
jnp(N_fourier+1), ynp(N_fourier+1);
Vector<double> jn_a(N_fourier+1),yn_a(N_fourier+1),
jnp_a(N_fourier+1), ynp_a(N_fourier+1);
Vector<complex<double> > h(N_fourier+1),h_a(N_fourier+1),
hp(N_fourier+1), hp_a(N_fourier+1);
// We want to compute N_fourier terms but the function
// may return fewer than that.
int n_actual=0;
CRBond_Bessel::bessjyna(N_fourier,sqrt(K_squared),n_actual,
&jn_a,&yn_a,
&jnp_a,&ynp_a);
// Shout if things went wrong
#ifdef PARANOID
if (n_actual!=int(N_fourier))
{
std::ostringstream error_stream;
error_stream << "CRBond_Bessel::bessjyna() only computed "
<< n_actual << " rather than " << N_fourier
<< " Bessel functions.\n";
throw OomphLibError(error_stream.str(),
OOMPH_CURRENT_FUNCTION,
OOMPH_EXCEPTION_LOCATION);
}
#endif
// Evaluate Hankel at actual radius
Hankel_functions_for_helmholtz_problem::Hankel_first(N_fourier,rr,h,hp);
// Evaluate Hankel at inner (unit) radius
Hankel_functions_for_helmholtz_problem::Hankel_first(N_fourier
,sqrt(K_squared),
h_a,hp_a);
// Compute the sum: Separate the computation of the negative
// and positive terms
for (unsigned i=0;i<N_fourier;i++)
{
u_ex-=pow(I,i)*h[i]*((jnp_a[i])/hp_a[i])*pow(exp(I*theta),i);
}
for (unsigned i=1;i<N_fourier;i++)
{
u_ex-=pow(I,i)*h[i]*((jnp_a[i])/hp_a[i])*pow(exp(-I*theta),i);
}
// Get the real & imaginary part of the result
u=real(u_ex);
u=imag(u_ex);
}// end of get_exact_u
std::complex< double > I(0.0, 1.0)
Imaginary unit.
double K_squared
Square of the wavenumber.
Definition: scattering.cc:64
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution for scattered field (vector returns real and impaginary parts).
Definition: scattering.cc:86
unsigned N_fourier
Number of terms used in the computation of the exact solution.
Definition: scattering.cc:68

Next we provide a function that computes the prescribed flux (normal derivative) of the solution, , evaluated on the surface of the unit disk.

/// Flux (normal derivative) on the unit disk
/// for a planar incoming wave
void prescribed_incoming_flux(const Vector<double>& x,
complex<double>& flux)
{
// Switch to polar coordinates
double r;
r=sqrt(x*x+x*x);
double theta;
theta=atan2(x,x);
// Argument of the Bessel/Hankel fcts
double rr=sqrt(K_squared)*r;
// Compute Bessel/Hankel functions
Vector<double> jn(N_fourier+1), yn(N_fourier+1),
jnp(N_fourier+1), ynp(N_fourier+1);
// We want to compute N_fourier terms but the function
// may return fewer than that.
int n_actual=0;
CRBond_Bessel::bessjyna(N_fourier,rr,n_actual,&jn,&yn,
&jnp,&ynp);
// Shout if things went wrong...
#ifdef PARANOID
if (n_actual!=int(N_fourier))
{
std::ostringstream error_stream;
error_stream << "CRBond_Bessel::bessjyna() only computed "
<< n_actual << " rather than " << N_fourier
<< " Bessel functions.\n";
throw OomphLibError(error_stream.str(),
OOMPH_CURRENT_FUNCTION,
OOMPH_EXCEPTION_LOCATION);
}
#endif
// Compute the sum: Separate the computation of the negative and
// positive terms
flux=std::complex<double>(0.0,0.0);
for (unsigned i=0;i<N_fourier;i++)
{
flux+=pow(I,i)*(sqrt(K_squared))*pow(exp(I*theta),i)*jnp[i];
}
for (unsigned i=1;i<N_fourier;i++)
{
flux+=pow(I,i)*(sqrt(K_squared))*pow(exp(-I*theta),i)*jnp[i];
}
}// end of prescribed_incoming_flux
} // end of namespace
void prescribed_incoming_flux(const Vector< double > &x, complex< double > &flux)
Flux (normal derivative) on the unit disk for a planar incoming wave.
Definition: scattering.cc:156

## The driver code

The driver code is very straightforward. We parse the command line to determine which boundary condition to use and set the flags in the global namespace accordingly.

//==========start_of_main=================================================
/// Solve 2D Helmholtz problem for scattering of a planar wave from a
/// unit disk
//========================================================================
int main(int argc, char **argv)
{
// Store command line arguments
CommandLineArgs::setup(argc,argv);
// Define case to be run
unsigned i_case=0;
CommandLineArgs::specify_command_line_flag("--case",&i_case);
// Parse command line
CommandLineArgs::parse_and_assign();
// Doc what has actually been specified on the command line
CommandLineArgs::doc_specified_flags();
// Now set flags accordingly
switch(i_case)
{
case 0:
break;
case 1:
break;
case 2:
break;
case 3:
break;
}
unsigned ABC_order
Flag to choose wich order to use.
Definition: scattering.cc:76
bool DtN_BC
Flag to choose the Dirichlet to Neumann BC or ABC BC.
Definition: scattering.cc:72
int main(int argc, char **argv)
Solve 2D Helmholtz problem for scattering of a planar wave from a unit disk.
Definition: scattering.cc:709

Next we build the problem, either with or without enabling spatial adaptivity and define the output directory.

//Set up the problem
//------------------
//Set up the problem with 2D nine-node elements from the
//QHelmholtzElement family.
problem;
#else
//Set up the problem with 2D nine-node elements from the
//QHelmholtzElement family.
problem;
#endif
// Create label for output
//------------------------
DocInfo doc_info;
// Set output directory
doc_info.set_directory("RESLT");
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: scattering.cc:223

Finally, we solve the problem and document the results.

// Solve the problem with Newton's method, allowing
#else
// Solve the problem with Newton's method
problem.newton_solve();
#endif
//Output solution
problem.doc_solution(doc_info);
} //end of main
void doc_solution(DocInfo &doc_info)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to.
Definition: scattering.cc:516

## The problem class

The problem class is very similar to that employed for the adaptive solution of the 2D Poisson equation with flux boundary conditions. The only difference is that we provide two separate meshes of FaceElements: one for the inner boundary where the HelmholtzFluxElements apply the Neumann condition (9), and one for the outer boundary where we apply the (approximate) Sommerfeld radiation condition. As discussed in section The Dirichlet-to-Neumann mapping (DtN) , we use the function actions_before_newton_convergence_check() to recompute the integral whenever the unknowns are updated during the Newton iteration.

//========= start_of_problem_class=====================================
/// Problem class to compute scattering of planar wave from unit disk
//=====================================================================
template<class ELEMENT>
class ScatteringProblem : public Problem
{
public:
/// Constructor
/// Destructor (empty)
/// Doc the solution. DocInfo object stores flags/labels for where the
/// output gets written to
void doc_solution(DocInfo& doc_info);
/// Update the problem specs before solve (empty)
/// Update the problem specs after solve (empty)
/// Recompute gamma integral before checking Newton residuals
{
{
}
}
/// Actions before adapt: Wipe the mesh of prescribed flux elements
/// Actions after adapt: Rebuild the mesh of prescribed flux elements
/// Create BC elements on boundary b of the Mesh pointed
/// to by bulk_mesh_pt and add them to the specified survace Mesh
const unsigned &b, Mesh* const &bulk_mesh_pt,
Mesh* const & helmholtz_outer_boundary_mesh_pt);
/// Create Helmholtz flux elements on boundary b of the Mesh pointed
/// to by bulk_mesh_pt and add them to the specified surface Mesh
void create_flux_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
Mesh* const & helmholtz_inner_boundary_mesh_pt);
/// Delete boundary face elements and wipe the surface mesh
void delete_face_elements( Mesh* const & boundary_mesh_pt);
/// Set pointer to prescribed-flux function for all
/// elements in the surface mesh on the surface of the unit disk
/// Set up boundary condition elements on outer boundary
/// Pointer to the "bulk" mesh
RefineableTwoDAnnularMesh<ELEMENT>* Bulk_mesh_pt;
#else
/// Pointer to the "bulk" mesh
TwoDAnnularMesh<ELEMENT>* Bulk_mesh_pt;
#endif
/// Pointer to mesh containing the DtN (or ABC) boundary
/// condition elements
HelmholtzDtNMesh<ELEMENT>* Helmholtz_outer_boundary_mesh_pt;
/// Pointer to the mesh containing
/// the Helmholtz inner boundary condition elements
}; // end of problem class
~ScatteringProblem()
Destructor (empty)
Definition: scattering.cc:231
Actions after adapt: Rebuild the mesh of prescribed flux elements.
Definition: scattering.cc:429
void create_outer_bc_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &helmholtz_outer_boundary_mesh_pt)
Create BC elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to the specified...
Definition: scattering.cc:645
Mesh * Helmholtz_inner_boundary_mesh_pt
Pointer to the mesh containing the Helmholtz inner boundary condition elements.
Definition: scattering.cc:297
void setup_outer_boundary()
Set up boundary condition elements on outer boundary.
Definition: scattering.cc:451
RefineableTwoDAnnularMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
Definition: scattering.cc:282
void actions_before_newton_solve()
Update the problem specs before solve (empty)
Definition: scattering.cc:238
ScatteringProblem()
Constructor.
Definition: scattering.cc:308
HelmholtzDtNMesh< ELEMENT > * Helmholtz_outer_boundary_mesh_pt
Pointer to mesh containing the DtN (or ABC) boundary condition elements.
Definition: scattering.cc:293
void delete_face_elements(Mesh *const &boundary_mesh_pt)
Delete boundary face elements and wipe the surface mesh.
Definition: scattering.cc:688
void actions_after_newton_solve()
Update the problem specs after solve (empty)
Definition: scattering.cc:241
Actions before adapt: Wipe the mesh of prescribed flux elements.
Definition: scattering.cc:413
void actions_before_newton_convergence_check()
Recompute gamma integral before checking Newton residuals.
Definition: scattering.cc:244
void set_prescribed_incoming_flux_pt()
Set pointer to prescribed-flux function for all elements in the surface mesh on the surface of the un...
Definition: scattering.cc:493
void create_flux_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &helmholtz_inner_boundary_mesh_pt)
Create Helmholtz flux elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to t...
Definition: scattering.cc:611

## The problem constructor

We start by building the bulk mesh, using the refineable or non-refineable version of the TwoDAnnularMesh, depending on the macro ADAPTIVE. (The error tolerances for the adaptive version are chosen such that the mesh is refined non-uniformly – with the default tolerances, oomph-lib's automatic mesh adaptation procedure refine the mesh uniformly.)

//=======start_of_constructor=============================================
/// Constructor for Helmholtz problem
//========================================================================
template<class ELEMENT>
{
// Setup "bulk" mesh
// # of elements in theta
unsigned n_theta=15;
// # of elements in radius
unsigned n_r=5;
double a=1.0;
// Thickness of annular computational domain
double h=0.5;
// Mesh is periodic
bool periodic=true;
// Full circle
double azimuthal_fraction=1.0;
// Build "bulk" mesh
Bulk_mesh_pt=
new RefineableTwoDAnnularMesh<ELEMENT>(periodic,
azimuthal_fraction,n_theta,n_r,a,h);
// Create/set error estimator
Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
// Choose error tolerances to force some uniform refinement
Bulk_mesh_pt->min_permitted_error()=0.004;
Bulk_mesh_pt->max_permitted_error()=0.01;
#else
// Build "bulk" mesh
Bulk_mesh_pt=
new TwoDAnnularMesh<ELEMENT>(periodic,
azimuthal_fraction,n_theta,n_r,a,h);
#endif
Radius of outer boundary (must be a circle!)
Definition: scattering.cc:79

Next we create the two (empty) meshes for the FaceElements,

// Pointer to mesh containing the Helmholtz outer boundary condition
// elements. Specify outer radius and number of Fourier terms to be
// used in gamma integral
Helmholtz_outer_boundary_mesh_pt =
new HelmholtzDtNMesh<ELEMENT>(a+h,GlobalParameters::N_fourier);
// Pointer to mesh containing the Helmholtz inner boundary condition
Helmholtz_inner_boundary_mesh_pt = new Mesh;

and populate them using the functions create_flux_elements(...) and create_outer_bc_elements(...).

// Create prescribed-flux elements from all elements that are
// adjacent to the inner boundary , but add them to a separate mesh.
create_flux_elements(0,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);
// Create outer boundary elements from all elements that are
// adjacent to the outer boundary , but add them to a separate mesh.
create_outer_bc_elements(2,Bulk_mesh_pt,Helmholtz_outer_boundary_mesh_pt);

We add the various (sub-)meshes to the problem and build the global mesh

// Add the several sub meshes to the problem
// Build the Problem's global mesh from its various sub-meshes
build_global_mesh();

Finally, we complete the build of the various elements by by passing pointers to the relevant quantities to them, and assign the equation numbers.

// Complete the build of all elements so they are fully functional
// Loop over the Helmholtz bulk elements to set up element-specific
// things that cannot be handled by constructor: Pass pointer to
// wave number squared
unsigned n_element = Bulk_mesh_pt->nelement();
for(unsigned e=0;e<n_element;e++)
{
// Upcast from GeneralisedElement to Helmholtz bulk element
ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
//Set the k_squared pointer
el_pt->k_squared_pt() = &GlobalParameters::K_squared;
}
// Set up elements on outer boundary
setup_outer_boundary();
// Set pointer to prescribed flux function for flux elements
set_prescribed_incoming_flux_pt();
// Setup equation numbering scheme
cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
} // end of constructor

The problem is now ready to be solved.

The mesh adaptation is driven by the error estimates for the bulk elements. The various FaceElements must therefore be removed from the global mesh before the adaptation takes place. We do this by calling the function delete_flux_elements(...) for the two face meshes, before rebuilding the Problem's global mesh.

/// Actions before adapt: Wipe the mesh of face elements
//========================================================================
template<class ELEMENT>
{
// Kill the flux elements and wipe the boundary meshs
delete_face_elements(Helmholtz_outer_boundary_mesh_pt);
delete_face_elements(Helmholtz_inner_boundary_mesh_pt);
// Rebuild the Problem's global mesh from its various sub-meshes
rebuild_global_mesh();

After the (bulk-)mesh has been adapted, the flux elements must be re-attached. This is done by calling the functions create_flux_elements(...) and create_outer_bc_elements, followed by rebuilding the Problem's global mesh. Finally, we complete the build of the FaceElements by calling the functions setup_outer_boundary() and set_prescribed_incoming_flux_pt().

/// Actions after adapt: Rebuild the face element meshes
//========================================================================
template<class ELEMENT>
{
// Create prescribed-flux elements and BC elements
// from all elements that are adjacent to the boundaries and add them to
// Helmholtz_boundary_meshes
create_outer_bc_elements(2,Bulk_mesh_pt,Helmholtz_outer_boundary_mesh_pt);
create_flux_elements(0,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);
// Rebuild the Problem's global mesh from its various sub-meshes
rebuild_global_mesh();
// Set pointer to prescribed flux function and DtN mesh
setup_outer_boundary();
set_prescribed_incoming_flux_pt();

## Delete flux elements

The helper function delete_face_elements() is used to delete all FaceElements in a given surface mesh before the mesh adaptation.

//============start_of_delete_face_elements================
/// Delete face elements and wipe the boundary mesh
//==========================================================
template<class ELEMENT>
delete_face_elements(Mesh* const & boundary_mesh_pt)
{
// Loop over the surface elements
unsigned n_element = boundary_mesh_pt->nelement();
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_outer_face_elements

## Creating the face elements

The functions create_flux_elements(...) and create_outer_bc_elements(...) create the FaceElements required to apply the boundary conditions on the inner and outer boundaries of the annular computational domain. They both loop over the bulk elements that are adjacent to the appropriate mesh boundary and attach the required FaceElements to their faces. The newly created FaceElements are then added to the appropriate mesh.

//============start_of_create_outer_bc_elements==============================
/// Create outer BC elements on the b-th boundary of
/// the Mesh object pointed to by bulk_mesh_pt and add the elements
/// to the Mesh object pointed to by helmholtz_outer_boundary_mesh_pt.
//===========================================================================
template<class ELEMENT>
create_outer_bc_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
Mesh* const & helmholtz_outer_boundary_mesh_pt)
{
// Loop over the bulk elements adjacent to boundary b?
unsigned n_element = bulk_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
ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
bulk_mesh_pt->boundary_element_pt(b,e));
//Find the index of the face of element e along boundary b
int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
// Build the corresponding outer flux element
// Dirichlet to Neumann boundary conditon
{
HelmholtzDtNBoundaryElement<ELEMENT>* flux_element_pt = new
HelmholtzDtNBoundaryElement<ELEMENT>(bulk_elem_pt,face_index);
//Add the flux boundary element to the helmholtz_outer_boundary_mesh
}
// ABCs BC
else
{
HelmholtzAbsorbingBCElement<ELEMENT>* flux_element_pt = new
HelmholtzAbsorbingBCElement<ELEMENT>(bulk_elem_pt,face_index);
//Add the flux boundary element to the helmholtz_outer_boundary_mesh
}
} //end of loop over bulk elements adjacent to boundary b
} // end of create_outer_bc_elements

(We omit the listing of the function create_flux_elements(...) because it is very similar. Feel free to inspect in the source code.)

## Post-processing

The post-processing function doc_solution(...) computes and outputs the total radiated power, and plots the computed and exact solutions (real and complex parts).

//=====================start_of_doc=======================================
/// Doc the solution: doc_info contains labels/output directory etc.
//========================================================================
template<class ELEMENT>
doc_info)
{
ofstream some_file,some_file2;
char filename;
// Number of plot points
unsigned npts;
npts=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_outer_boundary_mesh_pt->nelement();
for(unsigned e=0;e<nn_element;e++)
{
HelmholtzBCElementBase<ELEMENT> *el_pt =
dynamic_cast< HelmholtzBCElementBase<ELEMENT>*>(
Helmholtz_outer_boundary_mesh_pt->element_pt(e));
power += el_pt->global_power_contribution(some_file);
}
some_file.close();
oomph_info << "Total radiated power: " << power << std::endl;
// Output solution
//-----------------
sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
Bulk_mesh_pt->output(some_file,npts);
some_file.close();
// Output exact solution
//----------------------
sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
Bulk_mesh_pt->output_fct(some_file,npts,GlobalParameters::get_exact_u);
some_file.close();
double error,norm;
sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
Bulk_mesh_pt->compute_error(some_file,GlobalParameters::get_exact_u,
error,norm);
some_file.close();
// Doc L2 error and norm of solution
oomph_info << "\nNorm of error : " << sqrt(error) << std::endl;
oomph_info << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;

Finally, we create the data required to produce an animation of the actual (real) potential at 40 instants during a period of the oscillation.

// Do animation of Helmholtz solution
//-----------------------------------
unsigned nstep=40;
for (unsigned i=0;i<nstep;i++)
{
sprintf(filename,"%s/helmholtz_animation%i_frame%i.dat",
doc_info.directory().c_str(),
doc_info.number(),i);
some_file.open(filename);
sprintf(filename,"%s/exact_helmholtz_animation%i_frame%i.dat",
doc_info.directory().c_str(),
doc_info.number(),i);
some_file2.open(filename);
double phi=2.0*MathematicalConstants::Pi*double(i)/double(nstep-1);
unsigned nelem=Bulk_mesh_pt->nelement();
for (unsigned e=0;e<nelem;e++)
{
ELEMENT* el_pt=dynamic_cast<ELEMENT*>(
Bulk_mesh_pt->element_pt(e));
el_pt->output_real(some_file,phi,npts);
el_pt->output_real_fct(some_file2,phi,npts,
}
some_file.close();
some_file2.close();
}
} // end of doc

## The enumeration of the unknowns

As discussed in the introduction, most practically relevant solutions of the Helmholtz equation are complex valued. Since oomph-lib's solvers only deal with real (double precision) unknowns, the equations are separated into their real and imaginary parts. In the implementation of the Helmholtz elements, we store the real and imaginary parts of the solution as two separate values at each node. By default, the real and imaginary parts are accessible via Node::value(0) and Node::value(1). However, to facilitate the use of the elements in multi-physics problems we avoid accessing the unknowns directly in this manner but provide the virtual function

std::complex<unsigned> HelmholtzEquations<DIM>::u_index_helmholtz()

which returns a complex number made of the two unsigneds that indicate which nodal value represents the real and imaginary parts of the solution. This function may be overloaded in combined multi-physics elements in which a Helmholtz element is combined (by multiple inheritance) with another element, using the strategy described in the Boussinesq convection tutorial.

## Exercises

### Exploiting linearity

Confirm that the (costly) re-computation of the integral in actions_before_newton_convergence_check() after the first (and only) linear solve in the Newton iteration can be avoided by declaring the problem to be linear.

### The accuracy of the boundary condition elements

Explore the accuracy (and computational cost) of the various FaceElements that apply the Sommmerfeld radiation condition. In particular, confirm that the accuracy of the DtN boundary condition is (nearly) independent of the radius of the artificial outer boundary, whereas the accuracy of the ABC boundary condition can only be improved by increasing the size of the computational domain.