The results from the two-dimensional interior problem.

## 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.

### Keywords

- boundary element method
- Laplace’s equation
- Fortran

## 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:

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

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:

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

The terminology **q** on the boundary. The function *G* is known as Green’s function. Physically, *G(p, q)* represents the effect observed at point **p** of a unit source at point **q**. For the Laplace equation, Green’s function is denoted by G and is defined as

for two-dimensional Laplace problems, where

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**),

### 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 **p** on *Г*

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

In Eq. (8)*L* represents 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:

where

where

### 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^{−1} collocation in this text since it is derived by approximating the boundary functions by a constant on each panel. In this subsection the C^{−1} collocation method is briefly outlined.

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

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 ^{−1} collocation method involves representing the boundary function by a constant on each panel:

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:

where *e* is the unit function (*e***p**, are the numerical values of definite integrals that together can be interpreted as the discrete form of the *L* integral 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

for *i =* 1, 2 and *n* is obtained by putting

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

with the matrix components defined by

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:

for each point *m* domain points

where

### 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

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

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

where P is point

In general L2LC simply implements a Gaussian quadrature rule in order to determine the integral, using a higher-order rule when point

**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

**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.

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

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 |

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

where

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

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

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

which is easily shown to be a solution of Laplace’s equation by writing

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 |

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

where

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 |

### 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

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 |

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

where

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 |

## 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.

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

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

The boundary condition may be expressed in the following form:

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

### 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

The test problem is in file LSEMA_T. It consists of two circular coaxial parallel plates in the

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 |

### 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

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 |

## 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

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 |

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.