8  Subroutine H3LC

In this section the Fortran subroutine H3LC is described. The subroutine computes the discrete form of the three-dimensional Helmholtz integral operators. Details of the background methods employed by the subroutines are given. The subroutine is applied to the Burton and Miller equation for the test problem where the boundary is a cube and results are given.

8.1  Background Methods

The regular integrals that arise are approximated by a quadrature rule defined on a triangle. Laursen and Gellert (1978) contains a selection of Gauss-Legendre quadrature rules for the triangle. The non-regular integrals that arise in the formulae (41) and (42) are computed by the following methods. See Jaswon and Symm [11], Terai [34], Banerjee and Butterfield [4] and Kirkup [14] for the background to these methods.

The M0 and M0t operators have regular kernels, hence the aim is to find expressions for:

{ L0 ~
e
 
 }DG(p) =
ó
õ
DG 
  G0(p,q)  dSq  ,
(54)

{ N0 ~
e
 
 }DG(p; vp) =
vp
ó
õ
DG 
   G0
nq
 (p,q)  dSq   ,
(55)
where DG is a planar triangular element, p Î DG (though not on an edge or corner of the element). Let R(q) be the distance from p to the edge of the element for q Î [0, 2 p].

The integrals (54) and (55) may be written in the form:

{ L0 ~
e
 
 }DG (p) = 1
4 p
 ó
õ 		2 p

0 
 R(q) d q ,

{ N0 ~
e
 
 }DG (p; vp) = - 1
4 p
 ó
õ 		2 p

0 
 1
R(q)
 d q .
In order to evaluate the integrals, the triangular element DG is divided into three D1, D2 and D3 by joining the point p to the vertices. After some elementary analysis, we obtain

{L0 ~
e
 
 }DS (p) =
å
 D1, D2, D3 
 1
4 p
 R(56) sinB( logtan( B+A
2
 ) - logtan B
2
 )   and

{N0 ~
e
 
 }DS (p; vp) =
å
 D1, D2, D3 
 1
4 p
 cos(B+A) - cosB
R(0) sinB
 .

8.2  Test Program and Results

The test problem is in the file H3LC.FOR . The test problem consists of a boundary of a cube with vertices (1,1,1), (1,1,-1), (1,-1,1), (1,-1,-1), (-1,1,1), (-1,1,-1), (-1,-1,1) and (-1,-1,-1). The exterior solution and the Neumann boundary condition are determined by a point source at the centre of the cube. The cube is divided into 24 uniform triangular elements. Seven point Gaussian quadrature rules are applied to all but the computation of the diagonal components of the matrices, where a 3×7 point rule is used. Solutions are sought on the boundary and at the point (0,0,2) at wavenumbers k = 0, 1, i, 1+i. The values of a and b in the Burton and Miller equation (41) are set as a = 1, b = 1 for all four values of k.

The functions FNSQRT and FNEXP simply call the corresponding standard Fortran functions. Since the solution is similar at each collocation point, a comparison between one exact and one computed solution is given at one surface point and the exterior point. The results are given in table 2.

Table 2: Solution for the cube.
k a and b Point Exact Solution Computed Solution
0.0 a = 1, b = i Coll. Pt. 0.832050 0.829128 + i0.119394
0.0 a = 1, b = i (0,0,2) 0.500000 0.535075 - i0.000003
1.0 a = 1, b = i Coll. Pt. 0.300064 + i0.776060 0.163472 + i0.840514
1.0 a = 1, b = i (0,0,2) -0.208073 + i0.454649 -0.242757 + i0.516164
1.0i a = 1, b = i Coll. Pt. 0.250145 0.252107 + i0.048411
1.0i a = 1, b = i (0,0,2) 0.067668 0.068024 + i0.002300
1.0 + 1.0i a = 1, b = i Coll. Pt. 0.090211 + i0.233313 0.049252 + i0.246489
1.0 + 1.0i a = i, b = i (0,0,2) -0.002816 + i0.061530 -0.003732 + i0.061234