A Pilot Fortran Software Library for the Solution of Laplace’s Equation by the Boundary Element Method

Abstract

The boundary element method (BEM) is developed from the standpoint of software design. The Fortran language is used to produce a structured library for solving Laplace’s equation in various domain topologies and dimensions with generalised boundary conditions. Subroutines that compute the discrete Laplace operators, which are the core components for populating the matrices in the BEM, are developed. The main subroutines for solving Laplace’s equation in 2D, 3D and axisymmetric cases for open and closed boundaries are introduced. The methods are demonstrated on test problems.

1. Introduction

The boundary element method (BEM) has established itself as an important numerical technique for solving partial differential equations (PDEs) over the last half century [1, 2]. It distinguishes itself from competing methods, such as the finite element method (FEM) [3] in that the latter method requires a mesh of the domain, whereas the BEM only requires a mesh of the boundary (of the domain). The BEM is not as widely applicable as the FEM, particularly in that it is much more of a struggle to apply the BEM to non-linear problems. However, for problems to which the boundary element method is viable, the advantage of only requiring a boundary mesh is a significant one; the BEM is likely to be more efficient but also the relative simplicity of meshing, and the method is easier to use and is more accessible. This advantage is more notable for exterior problems; the domain is infinite, and ‘domain methods’ such as the FEM require special treatment, but for the BEM, only a (finite) boundary mesh is required. Computational methods may be combined or coupled [2].

The boundary element method is derived through the discretisation of an integral equation that is mathematically equivalent to the original partial differential equation. The essential reformulation of the PDE that underlies the BEM consists of an integral equation that is defined on the boundary of the domain and an integral that relates the boundary solution to the solution at points in the domain. The former is termed a boundary integral equation (BIE), and the BEM is often referred to as the boundary integral equation method or boundary integral method. There are two classes of boundary element method, termed the direct and indirect method. The direct method is based on Green’s second theorem, whereas the indirect method is based on describing the solution in terms of layer potentials. In this work the direct boundary element method is developed.

The simplest partial differential equation that is amenable to the BEM is Laplace’s equation:

$$\sum _{i=1}^N\frac{{\partial }^2\phi \left(\mathit{p}\right)}{\partial x_i^2}=0$$E1

where N is the dimension of the space or, more concisely,

$${\nabla }^2\phi =0.$$E2

Laplace’s equation therefore acts as a model problem for developing the BEM. Laplace’s equation also has a number of applications; steady-state heat conduction, steady-state electric potential, gravitation and groundwater flow [4, 5, 6, 7, 8, 9, 10, 11, 12, 13].

Initially, in this paper, the derivation of the direct boundary element method is introduced for the interior two-dimensional Laplace problem. The boundary element method is developed in Fortran for the 2D Laplace problem; then this is extended to axisymmetric three-dimensional problems and to both interior and exterior problems. The boundary element method can be extended to problems where the body being modelled is ‘thin’, like a screen or discontinuity, and these are also included. Test problems are applied to the codes, and the results are given for all problem classes. There are a number of studies on numerical error in the boundary element method [14, 15, 16].

There have been a number of works on coding the boundary element method [17, 18, 19]. The focus of this work is the algorithms and the software for solving Laplace problems by the BEM. As with the earlier works by the first author on Laplace and Helmholtz (acoustic) problems [20, 21, 22, 23, 24], this is about continuing with the development of a base library of methods and corresponding software. The codes and guides can be found on the first author’s website [25].

The codes have been developed in Fortran 77, but the language is just used to provide a simple template for exploring the methods and the organisation of coding. The algorithms and coding for Laplace’s equation considered in this work also provide a useful basis for the development of the BEM for other problems and add to the library of numerical software [26].

2. The BEM and the 2D interior Laplace problem

The Laplace equation provides a useful model problem for the boundary element method. The two-dimensional case is the simplest of these and is the best place to start to learn about the method. In this section the solution of Laplace’s equation in an interior domain by the direct BEM is outlined, and this also provides the foundation for the 3D BEM development in later sections.

2.1 Boundary integral equation formulation of the interior Laplace problem

Laplace’s equation (2) governs the interior domain D enclosed by a boundary S. The solution must also satisfy a boundary condition, and it is important in terms of maintaining the generality of the method that this is in a general (Robin) form:

$$\alpha \left(\mathit{p}\right)\mathit{\phi }\left(\mathit{p}\right)+\beta \left(\mathit{p}\right)\frac{\partial \phi }{\partial n_p}\left(\mathit{p}\right)=f\left(\mathit{p}\right)\left(\mathit{p}\in S\right).$$E3

In the direct BEM, Laplace’s equation is replaced by an equivalent integral equation of the form:

$${\int }_S\frac{\partial G\left(\mathit{p},\mathit{q}\right)}{\partial n_q}\phi \left(\mathit{q}\right)d{S}_q+\frac{1}{2}\phi \left(\mathit{q}\right)={\int }_{S}G\left(\mathit{p},\mathit{q}\right)\frac{\partial \phi \left(\mathit{q}\right)}{\partial n_q}d{S}_q\left(\mathit{p}\in S\right),$$E4

$${\int }_S\frac{\partial G\left(\mathit{p},\mathit{q}\right)}{\partial n_q}\phi \left(\mathit{q}\right)d{S}_q+\phi \left(\mathit{q}\right)={\int }_{S}G\left(\mathit{p},\mathit{q}\right)\frac{\partial \phi \left(\mathit{q}\right)}{\partial n_q}d{S}_q\left(\mathit{p}\in D\right).$$E5

The terminology $\frac{\partial \ast }{\partial n_q}$represents the partial derivative of the function* with respect to the unit outward normal at point qon the boundary. The function Gis known as Green’s function. Physically, G(p,q)represents the effect observed at point pof a unit source at point q. For the Laplace equation, Green’s function is denoted by G and is defined as

$$G\left(\mathbf{p},\mathbf{q}\right)=-\frac{1}{2\pi }lnr$$E6

for two-dimensional Laplace problems, where $r=\mid \mathit{q}-\mathit{p}\mid$.

Integral Eqs. (4) and (5) can be derived from the Laplace equation by applying Green’s second theorem. The power of the formulation lies in the fact that Eq. (4) relates the potential φand its derivative on the boundary alone; no reference is made to φ at points in the domain in this particular boundary integral equation. In a typical boundary value problem, we may be given φ(q), $\frac{\partial \phi \left(\mathit{q}\right)}{\partial n_q}$or a combination of such data on $S$. The boundary integral equation is a means of determining the unknown boundary function(s), followed by the domain solution from the given boundary data.

2.2 Operator notation

Operator notation is a useful shorthand in writing integral equations. Moreover, it will be shown that it is a very powerful notation in that it clearly demonstrates the connection between the integral equation and the linear system of equations that results from its discretisation.

Integral equations can always be written in terms of integral operators. For example, if ζis a function defined on a (closed or open) boundary Г, then applying the following operation to ζfor all points pon Г

$${\int }_{\mathrm{Г}}G\left(\mathit{p},\mathit{q}\right)\zeta \left(\mathit{q}\right)d{S}_q=\mu \left(\mathit{p}\right)\left(\mathit{p}\in S\right)$$E7

gives a function μ. This may be viewed as the application of an operator to the function ζ to return the function μ. More simply we may write

$${\left\{\mathit{L\zeta }\right\}}_{\mathrm{Г}}\left(\mathbf{p}\right)=\mathrm{\mu }\left(\mathbf{p}\right).$$E8

In Eq. (8)Lrepresents the integral operator, and the subscript (Г) refers to the domain of integration. Гis used as a variable, representing either a whole boundary or a part of the boundary. The other three important Laplace integral operators are defined as follows:

$${\left\{\mathit{M\zeta }\right\}}_{\mathrm{Г}}\left(\mathit{p}\right)={\int }_{\mathrm{Г}}\frac{\partial G\left(\mathit{p},\mathit{q}\right)}{\partial n_q}\zeta \left(\mathit{q}\right)d{S}_q,$$E9

$${\left\{M^t\zeta \right\}}_{\mathrm{Г}}\left(\mathit{p};{\mathit{v}}_p\right)=\frac{\partial }{\partial v_p}{\int }_{\mathrm{Г}}G\left(\mathit{p},\mathit{q}\right)\zeta \left(\mathit{q}\right)d{S}_q,$$E10

$${\left\{\mathit{N\zeta }\right\}}_{\mathrm{Г}}\left(\mathit{p};{\mathit{v}}_p\right)=\frac{\partial }{\partial v_p}{\int }_{\mathrm{Г}}\frac{\partial G\left(\mathit{p},\mathit{q}\right)}{\partial n_q}\zeta \left(\mathit{q}\right)d{S}_q,$$E11

where ${\mathit{v}}_p$is any unit vector. In operator notation of the previous subsection, the integral equation formulation (3) can be written in the following form:

$${\left\{\left(M+\frac{1}{2}I\right)\phi \right\}}_S\left(\mathit{p}\right)={\left\{\mathit{Lv}\right\}}_S\left(\mathit{p}\right)\left(\mathit{p}\in S\right),$$E12

$$\phi \left(\mathit{p}\right)={\left\{\mathit{Lv}\right\}}_S-{\left\{\mathit{M\phi }\right\}}_S\left(\mathit{p}\in D\right),$$E13

where $v\left(\mathit{q}\right)=\frac{\partial \phi \left(\mathit{q}\right)}{\partial n_q}$.

2.3 Direct boundary element method

For the direct boundary element method solution of the interior Laplace problem, that is, developed in this section, the initial stage involves solving boundary integral Eq. (4), returning (approximations to) both φ and ∂φ/∂n on S. The second stage of the BEM involved finding the solution at any chosen points in the domain D. The most straightforward method for solving integral equations like Eq. (4) is that of collocation. Collocation may be applied in a remarkably elementary form, which is termed C⁻¹ collocation in this text since it is derived by approximating the boundary functions by a constant on each panel. In this subsection the C⁻¹ collocation method is briefly outlined.

To begin with, the boundary Sis assumed to be expressed as a set of panels:

$$S\approx \stackrel{\sim }{S}=\sum _{j=1}^n\mathrm{∆}\stackrel{\sim }{S_j}.$$E14

Usually the panels have a characteristic form and cannot represent a given boundary exactly. For example, a two-dimensional boundary can be approximated by a set of straight lines. In order to complete the discretisation of the integral equations, the boundary functions also need to be approximated on each panel. In this work, it is the characteristics of the panel and the representation of the boundary function on the panel that together define the element in the boundary element method. By representing the boundary functions by a characteristic form on each panel, the boundary integral equations can be written as a linear system of equations of the form introduced earlier.

The term element refers not only to the form of ${\mathit{\Delta S}}_j$but also to the method of representing the boundary functions on ${\mathit{\Delta S}}_j$. The C⁻¹ collocation method involves representing the boundary function by a constant on each panel:

$$\phi \left(\mathit{p}\right)\approx {\phi }_j,v\left(\mathit{p}\right)\approx v_j\left(\mathit{p}\in \Delta {\stackrel{\sim }{S}}_j\right).$$E15

The substitution of representations of this form for the boundary functions in the integral equation reduces it to discrete form. The simplifications allow us to rewrite Eq. (11) as the approximation:

$$\sum _{j=1}^n{\left\{\left(M+\frac{1}{2}I\right)e\right\}}_{\Delta \widetilde{S_j}}\left(\mathit{p}\right){\phi }_j\approx \sum _{j=1}^n{\left\{\mathrm{Le}\right\}}_{\Delta \widetilde{S_j}}\left(\mathit{p}\right)v_j\left(\mathit{p}\in \stackrel{\sim }{S}\right),$$E16

where eis the unit function (e$\equiv$1). ${\left\{\mathit{Le}\right\}}_{\Delta \widetilde{S_j}}\left(\mathit{p}\right)$, for example, for a specific point p, are the numerical values of definite integrals that together can be interpreted as the discrete form of the Lintegral operator.

The constant approximation is taken to be the value of the boundary functions at the representative central point (the collocation point) on each panel. By finding the discrete forms of the relevant integral operators for all the collocation points, a system of the form

$$\sum _{j=1}^n{\left\{\left(M+\frac{1}{2}I\right)e\right\}}_{\Delta \widetilde{S_j}}\left({\mathit{p}}_{\mathit{Si}}\right){\phi }_j\approx \sum _{j=1}^n{\left\{\mathit{Le}\right\}}_{\Delta \widetilde{S_j}}\left({\mathit{p}}_{\mathit{Si}}\right)v_j\left({\mathit{p}}_{\mathit{Si}}\in \stackrel{\sim }{S}\right),$$E17

for i =1, 2 and nis obtained by putting $\mathit{p}={\mathit{p}}_{\mathit{Si}}$in the previous approximation. Note that because of the approximation of the boundary functions (and also the boundary approximation, if applicable), the discrete equivalent of Eq. (12) is an approximation relating the exact values of the boundary functions at the collocation points.

This system of approximations can now be written in the matrix-vector form:

$$\left(M_{\mathit{SS}}+\frac{1}{2}\mathrm{I}\right){\underset{¯}{\widehat{\phi }}}_S=L_{\mathit{SS}}{\underset{¯}{\widehat{v}}}_S,$$E18

with the matrix components defined by ${\left[L_{\mathit{SS}}\right]}_{\mathit{ij}}={\left\{\mathit{Le}\right\}}_{\Delta \widetilde{S_j}}\left({\mathit{p}}_{\mathit{Si}}\right),{\left[M_{\mathit{SS}}\right]}_{\mathit{ij}}={\left\{\mathit{Me}\right\}}_{\Delta \widetilde{S_j}}\left({\mathit{p}}_{\mathit{Si}}\right)$. The vectors ${\underset{¯}{\widehat{\phi }}}_S$and ${\underset{¯}{\widehat{v}}}_S$are representative or approximate values of φ and $v$at the collocation points. In the first stage of the boundary element, the system (18) is solved alongside the discrete form of the boundary condition (3):

$${\alpha }_i{\phi }_i+{\beta }_i{v}_i=f_i\text{for}i=1,2,\dots n.$$E19

The discrete forms are definite integrals that need to be computed usually by numerical integration. For the solution of Eqs. (18) and (19), the (approximation to) boundary data is known at the collocation points.

Once the (approximations to) functions on the boundary are known, after completing the initial stage of the boundary element method, the domain solution can be found. In the case of the interior Laplace problem, Eq. (13) will yield the domain solution. Similarly, the discrete equivalent of Eq. (11) may be derived:

$$\phi \left({\mathit{p}}_{\mathit{Di}}\right)=\sum _{j=1}^n{\left\{\mathit{Le}\right\}}_{\Delta \widetilde{S_j}}\left({\mathit{p}}_{\mathit{Di}}\right){\widehat{v}}_j-\sum _{j=1}^n{\left\{\mathit{Me}\right\}}_{\Delta \widetilde{S_j}}\left({\mathit{p}}_{\mathit{Di}}\right){\widehat{\phi }}_j\left({\mathit{p}}_{\mathit{Di}}\in \stackrel{\sim }{D}\right),$$E20

for each point ${\mathit{p}}_{\mathit{Di}}$in the domain $\stackrel{\sim }{D}$. Let the solution be sought at mdomain points ${\mathit{p}}_{\mathit{Di}}$for $i=1,2,\dots m$, and then the equation above, for all the domain points, is written as

$${\underset{¯}{\widehat{\phi }}}_D=L_{\mathit{DS}}{\underset{¯}{\widehat{v}}}_S-M_{\mathit{DS}}{\underset{¯}{\widehat{\phi }}}_S,$$E21

where ${\left[L_{\mathit{DS}}\right]}_{\mathit{ij}}={\left\{\mathit{Le}\right\}}_{\Delta \widetilde{S_j}}\left({\mathit{p}}_{\mathit{Di}}\right),{\left[M_{\mathit{DS}}\right]}_{\mathit{ij}}={\left\{\mathit{Me}\right\}}_{\Delta \widetilde{S_j}}\left({\mathit{p}}_{\mathit{Di}}\right)$.

2.4 LIBEM2 and L2LC modules, test problem and results

This section includes an outline of the Laplace interior BEM 2D (LIBEM2) and Laplace 2D linear boundary approximation constant element (L2LC) modules, and they are demonstrated by means of a test problem. The LIBEM2 module solves Laplace’s equation in an interior two-dimensional domain. L2LC is the most important component module.

L2LC.The L2LC module computes the discrete Laplace operators for two-dimensional problems. In the notation of this article, the routine computes ${\left\{\mathit{Le}\right\}}_{\Delta \widetilde{S_j}}\left(\mathit{p}\right)$, ${\left\{\mathit{Me}\right\}}_{\Delta \widetilde{S_j}}\left(\mathit{p}\right),{\left\{M^{t}e\right\}}_{\Delta \widetilde{S_j}}\left(\mathit{p}\right)$and ${\left\{\mathit{Ne}\right\}}_{\Delta \widetilde{S_j}}\left(\mathit{p}\right)$, where $\Delta \stackrel{\sim }{S_j}$is the panel that is the domain of integration and $\mathit{p}$is any point. The call to the subroutine has the following form:

SUBROUTINE L2LC(P,VECP,QA,QB,LPONEL, LVALID,EGEOM,LFAIL,

  * NEEDL,NEEDM,NEEDMT,NEEDN,L0,M0,M0T,N0),

where P is point $\mathit{p}$; VECP is a unit directional vector that passes through $\mathit{p}$; QA and QB are the points, either side of the panel and hence defining the panel; LPONEL is a logical switch that declares whether $\mathit{p}$is on the panel; and NEEDL is a logical switch that states whether the discrete $L$operator is required and similar to the other operators. The computed values for the integrals are output in L0, M0, M0T and N0. For the straightforward direct BEM, developed in the previous section, only $L$and $M$operators are required.

In general L2LC simply implements a Gaussian quadrature rule in order to determine the integral, using a higher-order rule when point $\mathit{p}$is close to the panel $\Delta \stackrel{\sim }{S_j}$. However, when point $\mathit{p}$is on the panel, then an exact integration is used [21, 22].

LIBEM2.The LIBEM2 module solves the interior Laplace problem and has the following form:

LIBEM2(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

  * MAXPOINTS,NPOINT,POINTS,

  * SALPHA,SBETA,SF,SINPHI,PINPHI,

  * LSOL,LVALID,TOLGEOM,

  * SPHI,SVEL,PPHI,

  * L_SS,M_SSPMHALFI,L_PS,M_PS,

  * PERM,XORY,C,workspace)

The boundary is set up through listing a set of nodal coordinates, and each panel is determined through the two nodal indices for the endpoints of the panel. The nodal coordinates are input through NODES and the panel information through PANELS. The nodes are oriented clockwise on each panel for an outer boundary and anticlockwise for any inner boundary. Usually, a solution in the domain is sought, and for this a set of (interior) domain points are set in POINTS. The boundary condition is set with the parameters SALPHA, SBETA and SF, setting ${\alpha }_i$, ${\beta }_i$and $f_i$values in Eq. (19) for each panel.

Test problem.The test problem is that of solving Laplace’s equation on a unit square with the boundary conditions defined as shown in Figure 1. The solution is sought at the five interior points (0.25, 0.25), (0.75, 0.25), (0.25, 0.75), (0.75, 0.75) and (0.5, 0.5), and these are also illustrated in the figure.

Figure 1.

The test problem of the unit square domain with boundary conditions.

The test problem is set up in the file LIBEM2_T.FOR. The boundary is defined by 32 nodes and panels. The nodes are indexed, starting with 1: (0.0, 0.0), 2: (0.0, 0.125), 3: (0.0, 0.25) and continue clockwise around the boundary until the final node 32: (0.125, 0.0). The panels are similarly set up in the clockwise sense with panel 1:1–2 (panel 1 links node 1 with node 2) and 2:2–3 until the final panel 32:32–1, linking the final node with the first node to complete the boundary. The boundary conditions shown in Figure 1 are then applied.

The exact solution is $\phi \left(\mathit{p}\right)=10+10x$; this is clearly a solution of Laplace’s equation and satisfies the boundary conditions. The exact solution at the interior points is therefore $\phi =12.5$at the two points on the left, $\phi =17.5$at the two points on the right and $\phi =15.0$at the central point. The exact and computed results are shown in Table 1.

Index	Point	Exact	Computed (5 d.p.)
1	(0.25, 0.25)	12.5	12.49568
2	(0.75, 0.25)	17.5	17.50432
3	(0.25, 0.75)	12.5	12.49568
4	(0.75, 0.75)	17.5	17.50432
5	(0.5, 0.5)	15.0	15.00000

Table 1.

The results from the two-dimensional interior problem.

A set of nodal coordinates and each panel is determined through the two nodal indices for the endpoints of the panel.

The results from this test problem are also intuitively correct. With the left and right sides of the square at different potentials, it is common sense to expect the potential in the middle to be halfway between etc. The potentials can—most simply—be interpreted as temperatures in a steady-state heat conduction problem.

3. The BEM and 3D Laplace problems

In this section, the boundary element method—introduced for two-dimensional problems in the previous section—is extended to include three-dimensional problems in this section. In this section, the three-dimensional boundary may be general, but the special case of axisymmetric problems is also developed in the modules LBEM3 and LBEMA. The modules can solve interior and exterior Laplace problems.

For interior three-dimensional problems, the basic integral formulation is the same as for 2D problems (12) and (13), except that Green’s function for three-dimensional Laplace problems is

$$G\left(\mathit{p},\mathit{q}\right)=\frac{1}{4\mathit{\pi r}},$$E22

where $r$is the distance between points $\mathit{p}$and $\mathit{q}$and the integrals are over surfaces rather than lines. The equations for the exterior problem are the same as for the interior problem, but for some changes of sign

$${\left\{\left(M-\frac{1}{2}I\right)\phi \right\}}_S\left(\mathit{p}\right)={\left\{\mathit{Lv}\right\}}_S\left(\mathit{p}\right)\left(\mathit{p}\in S\right),$$E23

$$\phi \left(\mathit{p}\right)={\left\{\mathit{M\phi }\right\}}_S-{\left\{\mathit{Lv}\right\}}_S\left(\mathit{p}\in E\right).$$E24

For general three-dimensional problems, the simplest elements are triangular panels, and for axisymmetric problems, they are lateral sections of a cone, with surface functions approximated by a constant on each panel.

3.1 LBEMA and L3ALC modules, test problems and results

Let us start on the introduction of three-dimensional problems with the axisymmetric codes. These codes are used in a very similar manner. As with the two-dimensional problem, the component module L3ALC computes the integrals over the panels and is called as follows:

SUBROUTINE L3ALC(P,VECP,QA,QB,LPONEL,LVALID,EGEOM,LFAIL,

  * NEEDL,NEEDM,NEEDMT,NEEDN,DISL,DISM,DISMT,DISN).

For axisymmetric problems, the surface is defined by conical panels, which are defined by piecewise straight lines along the generator. The parameters follow a similar pattern as L2LC, except the points and vectors are in cylindrical $\left(r,z\right)$coordinates. QA and QB are the two points either side of the panel on the generator.

The LBEMA subroutine computes the solution of the Laplace equation by the direct boundary element method and has the following form:

LBEMA(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

  * LINTERIOR,MAXPOINTS,NPOINT,POINTS,

  * SALPHA,SBETA,SF,SINPHI,PINPHI,

  * LSOL,LVALID,TOLGEOM,

  * SPHI,SVEL,PPHI,

  * L_SS,M_SSPMHALFI,L_PS,M_PS,

  * PERM,XORY,C,workspace)

In LBEMA, NODES lists the $\left(r,z\right)$coordinates of the nodes on the generator of the surface, and PANELS states the two nodes that together define each panel. LINTERIOR is a logical input, which is set to TRUE if an interior problem is to be solved and FALSE for an interior problem.

The interior test problem is in file LBEMA_IT. The test problem is the unit sphere with the exact solution:

$$\phi =r^2-2z^2,$$E25

which is easily shown to be a solution of Laplace’s equation by writing $r^2=x^2+y^2.$A Dirichlet boundary condition is applied; the solution is sought at four interior points, and the results for 18 elements are given in Table 2.

Index	Point	Exact	Computed (4 d.p.)
1	(0.0, 0.0)	0.0	−0.0013
2	(0.0, 0.5)	−0.5	−0.4995
3	(0.0, −0.5)	−0.5	−0.4995
4	(0.5, 0.0)	0.25	0.2477

Table 2.

The results from the axisymmetric interior problem.

The exterior test problem is in file LBEMA_ET. The test problem is the unit sphere (approximated by 18 elements) with the exact solution:

$$\mathrm{\phi }=\frac{1}{\mathrm{r}},$$E26

where $r$is the distance from the origin. $\phi$is a solution of Laplace’s equation as it is a simple multiplication of Green’s function (22). A Dirichlet boundary condition is applied to the upper hemisphere, and a Neumann boundary condition is applied on the lower hemisphere. The solution is sought at four interior points, and the results are given in Table 3.

Index	Point	Exact (4 d.p.)	Computed (4 d.p.)
1	(0.0, 2.0)	0.5	0.4986
2	(1.0, 1.0)	0.7071	0.7051
3	(0.0, 100.0)	0.0100	0.0100

Table 3.

The results from the axisymmetric exterior problem.

3.2 LBEM3 and L3LC modules, test problems and results

The LBEM3 and L3LC subroutines implement the boundary element method for general three-dimensional problems. As with the two-dimensional and axisymmetric codes, the component module L3LC computes the integrals over the panels. The L3LC subroutine is called as follows:

SUBROUTINE L3LC(P,VECP,QA,QB,QC,LPONEL,LVALID,EGEOM,LFAIL,

 * NEEDL,NEEDM,NEEDMT,NEEDN,DISL,DISM,DISMT,DISN)

The parameters follow a similar purpose as did in the L2LC, except that the points and vectors have three values. QA, QB and QC are the coordinates of the vertices of the triangular panel.

The LBEM3 module solves Laplace’s equation in a general interior or exterior three-dimensional domain and is called as follows:

LBEM3(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

 * LINTERIOR,MAXPOINTS,NPOINT,POINTS,

 * SALPHA,SBETA,SF,SINPHI,PINPHI,

 * LSOL,LVALID,TOLGEOM,

 * SPHI,SVEL,PPHI,

 * L_SS,M_SSPMHALFI,L_PS,M_PS,

 * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3)

As with LIBEM2 and LBEMA, NODES and PANELS define the boundary. However, in this case, NODES lists the three coordinates of each surface node, and PANELS lists the three nodal indices that make up each triangular panel.

The interior test problem is that of a unit sphere approximated by 36 triangular panels. The exact solution that is applied as a Dirichlet boundary condition is

$$\phi =x+y+z.$$E27

The results at four interior points are given in Table 4.

Index	Point	Exact	Computed (4 d.p.)
1	(0.5, 0.0, 0.0)	0.5	0.4772
2	(0.0, 0.5, 0.0)	0.5	0.4836
3	(0.0, 0.0, 0.5)	0.5	0.4817
4	(0.1, 0.2, 0.3)	0.6	0.5802

Table 4.

The results from the three-dimensional interior problem.

The exterior test problem is that of a unit sphere approximated by 36 triangular panels, as in the previous test. The exact solution that is applied as a Dirichlet boundary condition is

$$\phi =\frac{1}{r},$$E28

where $r$is the distance from $\left(0,0,0.5\right).$The results at four exterior points are given in Table 5.

Index	Point	Exact (4 d.p.)	Computed (4 d.p.)
1	(2.0, 0.0, 0.0)	0.4851	0.4969
2	(0.0, 4.0, 0.0)	0.2481	0.2536
3	(0.0, 0.0, 8.0)	0.1333	0.1360
4	(2.0, 2.0, 2.0)	0.3123	0.3189

Table 5.

The results from the three-dimensional exterior problem.

4. The solution of the 3D Laplace equation around a thin shell

Let us now consider the integral equation formulation for thin shells. An illustration of a typical problem of a hollow hemispherical cap is illustrated in Figure 2. In the traditional boundary element method, the boundaries are closed. This analysis and software design extends the boundary element method to open boundaries or discontinuities in the potential field.

Figure 2.

A hemispherical shell.

In this section, the integral equations that are a reformulation of Laplace’s equation surrounding a thin shell are stated. Fortran codes that implement the boundary element method for axisymmetric and general three-dimensional problems are outlined in this section and demonstrated on simple test problems, similar to the modelling of the steady-state electric field in a capacitor in Kirkup [9].

4.1 Integral equations and boundary element equations for thin shells

Following the work of Warham [27], the first step is to designate an ‘upper’ and ‘lower’ surface of a shell $\Gamma$and denote them by ‘+’ and ‘−’. We then introduce the quantities of difference and average of the potential and its normal derivative across the surface:

$$\delta \left(\mathit{p}\right)=\phi \left({\mathit{p}}_{+}\right)-\phi \left({\mathit{p}}_{+}\right)\left(\mathit{p}\in \Gamma \right),$$E29

$$\Phi \left(\mathit{p}\right)=\frac{1}{2}\left(\phi \left({\mathit{p}}_{+}\right)+\phi \left({\mathit{p}}_{+}\right)\right)\left(\mathit{p}\in \Gamma \right),$$E30

$$\nu \left(\mathit{p}\right)=v\left({\mathit{p}}_{+}\right)+v\left({\mathit{p}}_{+}\right)\left(\mathit{p}\in \Gamma \right),$$E31

$$V\left(\mathit{p}\right)=\frac{1}{2}\left(v\left({\mathit{p}}_{+}\right)-v\left({\mathit{p}}_{+}\right)\right)\left(\mathit{p}\in \Gamma \right).$$E32

The integral equation formulations for the Laplace equation in the exterior domain can now be written using the operator notation introduced earlier:

$$\phi \left(\mathit{p}\right)={\left\{\mathit{M\delta }\right\}}_{\Gamma }\left(\mathit{p}\right)-{\left\{\mathit{L\nu }\right\}}_{\Gamma }\left(\mathit{p}\right)\left(\mathit{p}\in E\right),$$E33

$$\Phi \left(\mathit{p}\right)={\left\{\mathit{M\delta }\right\}}_{\Gamma }\left(\mathit{p}\right)-{\left\{\mathit{L\nu }\right\}}_{\Gamma }\left(\mathit{p}\right)\left(\mathit{p}\in \Gamma \right),$$E34

$$V\left(\mathit{p}\right)={\left\{\mathit{N\delta }\right\}}_{\Gamma }\left(\mathit{p}\right)-{\left\{M^t\nu \right\}}_{\Gamma }\left(\mathit{p}\right)\left(\mathit{p}\in \Gamma \right).$$E35

The boundary condition may be expressed in the following form:

$$\alpha \left(\mathit{p}\right)\mathit{\delta }\left(\mathit{p}\right)+\beta \left(\mathit{p}\right)\nu \left(\mathit{p}\right)=f\left(\mathit{p}\right)\left(\mathit{p}\in \Gamma \right),$$E36

$$A\left(\mathit{p}\right)\Phi \left(\mathit{p}\right)+\beta \left(\mathit{p}\right)V\left(\mathit{p}\right)=F\left(\mathit{p}\right)\left(\mathit{p}\in \Gamma \right).$$E37

The discrete equivalents of Eq. (21) are as follows:

$${\underset{¯}{\widehat{\phi }}}_E=M_{\mathit{E\Gamma }}{\underset{¯}{\widehat{\delta }}}_{\Gamma }-L_{\mathit{E\Gamma }}{\underset{¯}{\widehat{\nu }}}_{\Gamma },$$E38

$${\underset{¯}{\widehat{\Phi }}}_{\Gamma }=M_{\mathit{\Gamma \Gamma }}{\underset{¯}{\widehat{\delta }}}_{\Gamma }-L_{\mathit{\Gamma \Gamma }}{\underset{¯}{\widehat{\nu }}}_{\Gamma },$$E39

$${\underset{¯}{\widehat{V}}}_{\Gamma }=N_{\mathit{\Gamma \Gamma }}{\underset{¯}{\widehat{\delta }}}_{\Gamma }-M_{\mathit{\Gamma \Gamma }}^t{\underset{¯}{\widehat{\nu }}}_{\Gamma }.$$E40

4.2 LSEMA module, test problem and results

The LSEMA subroutine computes the solution of Laplace’s equation surrounding thin shells or discontinuities. As with the LBEMA, the subroutine relies on L3ALC to compute the matrix components in the systems (38)–(40). In this subsection, the LSEMA routine is demonstrated through solving a test problem.

The module LSEMA has the form:

LSEMA(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

 * MAXPOINTS,NPOINT,POINTS,

 * HA,HB,HF,HAA,HBB,HFF,

 * HIPHI,HIVEL,PINPHI,

 * LSOL,LVALID,TOLGEOM,

 * PHIDIF,PHIAV,VELDIF,VELAV,PPHI,

 * AMAT,BMAT,L_EH, M_EH,

 * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3).

The LSEMA parameters are similar to the LBEMA ones. However the expressions of the boundary condition and the boundary function are different.

HA stores the values of $\alpha$on the shell panels, similarly HB, $\beta$; HAA, $A$; and HBB, $B.$The main output from the subroutine is PHIDIF that corresponds to ${\underset{¯}{\widehat{\delta }}}_{\Gamma }$; PHIAV, ${\underset{¯}{\widehat{\Phi }}}_{\Gamma }$; VELDIF, ${\underset{¯}{\widehat{\nu }}}_{\Gamma }$; VELAV, ${\underset{¯}{\widehat{V}}}_{\Gamma }$; and PPHI, ${\underset{¯}{\widehat{\phi }}}_E$.

The test problem is in file LSEMA_T. It consists of two circular coaxial parallel plates in the $r,\theta$plane, of radius 1.0 and a distance of 0.1 apart in the planes where $z=0.0$and $z=0.1.$A Dirichlet boundary condition is applied to both plates. On the plate at $z=0.0$, the potential of 0.0 is applied, and a potential ($\mathit{\delta }$ = 0, $\Phi =0$) of 1.0 is applied on the other plate ($\mathit{\delta }$ = 0, $\Phi =1$). A complete analytic solution is not available. However in the central region between the plates, a simple gradient of potential is intuitive, as discussed. The results from the test problem are listed in Table 6.

Index	Point	Expected (4 d.p.)	Computed (4 d.p.)
1	(0.0, 0.025)	0.25	0.2495
2	(0.0, 0.05)	0.5	0.5000
3	(0.0, 0.075)	0.75	0.7506

Table 6.

The results from the axisymmetric shell problem.

4.3 LSEM3 module, test problem and results

The LSEM3 module solves Laplace’s equation exterior to a thin shell in three dimensions. The subroutine call has the following form:

LSEM3(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

 * MAXPOINTS,NPOINT,POINTS,

 * HA,HB,HF,HAA,HBB,HFF,

 * HINPHI,HINVEL,PINPHI,

 * LSOL,LVALID,TOLGEOM,

 * PHIDIF,PHIAV,VELDIF,VELAV,PPHI,

 * AMAT,BMAT,L_EH,M_EH,

 * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3)

The definition of the important parameters can be found from the previous notes on LBEM3 and LSEMA. The test problem is in the file LSEM3_T, and it is similar to the test problem for LSEMA. This time, the open boundaries are two unit square plates of in $x-y$planes. The two squares are 0.1 apart: one is at a potential of zero and the other is at a potential of one. The squares are each divided into 32 panels. The results at points between the squares, along a central axis, are shown in Table 7.

Index	Point	Expected (4 d.p.)	Computed (4 d.p.)
1	(0.5, 0.5, 0.1)	0.1	0.0962
2	(0.5, 0.5, 0.3)	0.3	0.02994
3	(0.5, 0.5, 0.5)	0.5	0.5000
4	(0.5, 0.5, 0.7)	0.7	0.7006
5	(0.5, 0.5, 0.9)	0.9	0.9041

Table 7.

The results from the three-dimensional shell problem.

5. Conclusions

In this paper a design of a software library has been set out and implemented in Fortran. In taking a ‘library’ approach, components can be developed that can be shared. There is, therefore, an overall reduction in coding, in line with good software engineering practice. For the three-dimensional problems, it is shown how exterior problems can be solved with the same code as interior problems. It is also shown how the core discrete operator components can be reused for codes solving problems in the same dimensional space. The method for solving the linear system of equations can also often be shared, as with LU factorisation, applied in these codes. A test problem has been developed in order to demonstrate each code. The library of codes and the way they are linked are set out in Appendix.

There are several areas for further development. It is good for software engineering also to widen participation to provide strong validation in the BEM, so that errors, for example, in the boundary mesh are noted before executing the BEM. In this work the validation is developed through the VGEOM* modules.

In this paper, the BEM codes have been applied to a set of simple test problems. It would be useful if a standard library of test problems emerged, so that all existing and future codes can be benchmarked against the same tests, with information such as error and processing time. More complex geometries—such as multiple surfaces in exterior problems or cavities in the domain for interior problems—would benefit from standard test problems. The codes are also adaptable to problems in which there is an existing field that the boundary and boundary conditions modify (via the *INPHI and *INVEL parameters), but these have not been tested.

Central to the efficiency of the method, as the number of elements increases, is the method for solving the linear system of equations and the method of storing the matrices. Computing the matrices in the BEM takes $O\left(n^2\right)$time and memory. Solving the linear system by a direct method, like LU factorisation used in this work, takes $O\left(n^3\right)$time. Hence, in order to scale up the method, LU factorisation needs to be replaced by an interative method, and methods of storing and computing the matrices may also become an issue.

In the software engineering approach in this work, a generalised form of the boundary condition is also operational, and interior and exterior problems in 3D are dealt with in the same code. Further generality may be achieved by forming a hybrid of the method that allows both open and closed surfaces [28, 29, 30].

The main codes for solving Laplace problems by the boundary element method in this work are LIBEM2 for the two-dimensional problem interior to a closed boundary, LBEM3 for the general three-dimensional problem interior or exterior to a closed boundary, LBEMA for the axisymmetric three-dimensional problem interior or exterior to a closed boundary, LSEM3 for the general three-dimensional problem exterior to an open boundary and LSEMA for the axisymmetric three-dimensional problem exterior to an open boundary. The linkage between these and the supporting codes in the library is shown in Table 8.

File/code	Purpose of module	LIBEM2	LBEMA	LBEM3	LSEMA	LSEM3
L2LC	Computes the discrete Laplace operators (2D)	X
L3ALC	Computes the discrete Laplace operators (axisym)		X		X
L3LC	Computes the discrete Laplace operators (3D)			X		X
GLS2	Solves a generalised linear system of equations	X	X	X	X	X
LUFAC	Carries out LU factorisation of the matrix	X	X	X	X	X
LUFBSUB	Carries out forward and back substitution	X	X	X	X	X
GEOM2D	Geometrical operations (2D)	X	X		X
GEOM3D	Geometrical operations (3D)		X	X	X	X
GLRULES	Gauss-Legendre quadrature rules	X	X		X
GLT7	7-point Gaussian quadrature rule for triangle			X		X
GLT25	25-point Gaussian quadrature rule for triangle			X		X
VGEOM2	Verifies the geometry (2D)	X
VGEOMA	Verifies the geometry (axisym)		X		X
VGEOM3	Verifies the geometry (3D)			X		X
VG2LC	Verifies the use of the L2LC module	X
L3ALCC	Copy of L3ALC (to fake recursion)		X		X

Table 8.

The main codes and supporting library.

The main subroutines have the control parameters LSOL, LVALID and TOLGEOM. LSOL is set to TRUE if the full solution is sought and FALSE if the linear system is the output. LVALID is set to TRUE if validation is required and FALSE if it is not. TOLGEOM sets the geometrical tolerance.

The GLS algorithm in file GLS2 carries out a column-swapping method [31] in order to prepare the linear system for solution by a standard method. The standard method in this work is LU factorization and back substitution in files LUFAC and LUFBSUBS.