Open access peer-reviewed chapter

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

Written By

Stephen Kirkup and Javad Yazdani

Submitted: 15 November 2018 Reviewed: 24 April 2019 Published: 06 May 2020

DOI: 10.5772/intechopen.86507

From the Edited Volume

Numerical Modeling and Computer Simulation

Edited by Dragan M. Cvetković and Gunvant A. Birajdar

Chapter metrics overview

948 Chapter Downloads

View Full Metrics

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:

i=1N2φpxi2=0E1

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

2φ=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].

Advertisement

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:

αpφp+βpφnpp=fppS.E3

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

SGpqnqφqdSq+12φq=SGpqφqnqdSqpS,E4
SGpqnqφqdSq+φq=SGpqφqnqdSqpD.E5

The terminology nq represents the partial derivative of the function* with respect to the unit outward normal at point 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

Gpq=12πlnrE6

for two-dimensional Laplace problems, where r=qp.

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), φqnq 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 p on Г

ГGpqζqdSq=μppSE7

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

Гp=μp.E8

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:

Гp=ГGpqnqζqdSq,E9
MtζГpvp=vpГGpqζqdSq,E10
Гpvp=vpГGpqnqζqdSq,E11

where vp is any unit vector. In operator notation of the previous subsection, the integral equation formulation (3) can be written in the following form:

M+12IφSp=LvSppS,E12
φp=LvSSpD,E13

where vq=φqnq.

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:

SS=j=1nSj.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 ΔSj but also to the method of representing the boundary functions on ΔSj. The C−1 collocation method involves representing the boundary function by a constant on each panel:

φpφj,vpvjpΔSj.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:

j=1nM+12IeΔSj˜pφjj=1nLeΔSj˜pvjpS,E16

where e is the unit function (e 1). LeΔSj˜p, 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 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

j=1nM+12IeΔSj˜pSiφjj=1nLeΔSj˜pSivjpSiS,E17

for i = 1, 2 and n is obtained by putting p=pSi 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:

MSS+12Iφ̂¯S=LSSv̂¯S,E18

with the matrix components defined by LSSij=LeΔSj˜pSi,MSSij=MeΔSj˜pSi. The vectors φ̂¯S and 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):

αiφi+βivi=fifori=1,2,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:

φpDi=j=1nLeΔSj˜pDiv̂jj=1nMeΔSj˜pDiφ̂jpDiD,E20

for each point pDi in the domain D. Let the solution be sought at m domain points pDi for i=1,2,m, and then the equation above, for all the domain points, is written as

φ̂¯D=LDSv̂¯SMDSφ̂¯S,E21

where LDSij=LeΔSj˜pDi,MDSij=MeΔSj˜pDi.

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 LeΔSj˜p, MeΔSj˜p,MteΔSj˜p and NeΔSj˜p, where ΔSj is the panel that is the domain of integration and 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 p; VECP is a unit directional vector that passes through 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 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 Land 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 p is close to the panel ΔSj. However, when point 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 αi, βi and fi 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 φp=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 φ=12.5 at the two points on the left, φ=17.5 at the two points on the right and φ=15.0 at the central point. The exact and computed results are shown in Table 1.

IndexPointExactComputed (5 d.p.)
1(0.25, 0.25)12.512.49568
2(0.75, 0.25)17.517.50432
3(0.25, 0.75)12.512.49568
4(0.75, 0.75)17.517.50432
5(0.5, 0.5)15.015.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.

Advertisement

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

Gpq=14πr,E22

where r is the distance between points p and 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

M12IφSp=LvSppS,E23
φp=SLvSpE.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 rz 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 rz 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:

φ=r22z2,E25

which is easily shown to be a solution of Laplace’s equation by writing r2=x2+y2. 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.

IndexPointExactComputed (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.250.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:

φ=1r,E26

where r is the distance from the origin. φ 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.

IndexPointExact (4 d.p.)Computed (4 d.p.)
1(0.0, 2.0)0.50.4986
2(1.0, 1.0)0.70710.7051
3(0.0, 100.0)0.01000.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

φ=x+y+z.E27

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

IndexPointExactComputed (4 d.p.)
1(0.5, 0.0, 0.0)0.50.4772
2(0.0, 0.5, 0.0)0.50.4836
3(0.0, 0.0, 0.5)0.50.4817
4(0.1, 0.2, 0.3)0.60.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

φ=1r,E28

where ris the distance from 000.5. The results at four exterior points are given in Table 5.

IndexPointExact (4 d.p.)Computed (4 d.p.)
1(2.0, 0.0, 0.0)0.48510.4969
2(0.0, 4.0, 0.0)0.24810.2536
3(0.0, 0.0, 8.0)0.13330.1360
4(2.0, 2.0, 2.0)0.31230.3189

Table 5.

The results from the three-dimensional exterior problem.

Advertisement

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 Γ and denote them by ‘+’ and ‘−’. We then introduce the quantities of difference and average of the potential and its normal derivative across the surface:

δp=φp+φp+pΓ,E29
Φp=12φp++φp+pΓ,E30
νp=vp++vp+pΓ,E31
Vp=12vp+vp+pΓ.E32

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

φp=ΓpΓppE,E33
Φp=ΓpΓppΓ,E34
Vp=ΓpMtνΓppΓ.E35

The boundary condition may be expressed in the following form:

αpδp+βpνp=fppΓ,E36
ApΦp+βpVp=FppΓ.E37

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

φ̂¯E=Mδ̂¯ΓLν̂¯Γ,E38
Φ̂¯Γ=MΓΓδ̂¯ΓLΓΓν̂¯Γ,E39
V̂¯Γ=NΓΓδ̂¯ΓMΓΓtν̂¯Γ.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 α on the shell panels, similarly HB, β; HAA, A; and HBB, B. The main output from the subroutine is PHIDIF that corresponds to δ̂¯Γ; PHIAV, Φ̂¯Γ; VELDIF, ν̂¯Γ; VELAV, V̂¯Γ; and PPHI, φ̂¯E.

The test problem is in file LSEMA_T. It consists of two circular coaxial parallel plates in the r,θ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 (δ = 0, Φ=0) of 1.0 is applied on the other plate (δ = 0, Φ=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.

IndexPointExpected (4 d.p.)Computed (4 d.p.)
1(0.0, 0.025)0.250.2495
2(0.0, 0.05)0.50.5000
3(0.0, 0.075)0.750.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 xy 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.

IndexPointExpected (4 d.p.)Computed (4 d.p.)
1(0.5, 0.5, 0.1)0.10.0962
2(0.5, 0.5, 0.3)0.30.02994
3(0.5, 0.5, 0.5)0.50.5000
4(0.5, 0.5, 0.7)0.70.7006
5(0.5, 0.5, 0.9)0.90.9041

Table 7.

The results from the three-dimensional shell problem.

Advertisement

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 On2 time and memory. Solving the linear system by a direct method, like LU factorisation used in this work, takes On3 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].

Advertisement

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/codePurpose of moduleLIBEM2LBEMALBEM3LSEMALSEM3
L2LCComputes the discrete Laplace operators (2D)X
L3ALCComputes the discrete Laplace operators (axisym)XX
L3LCComputes the discrete Laplace operators (3D)XX
GLS2Solves a generalised linear system of equationsXXXXX
LUFACCarries out LU factorisation of the matrixXXXXX
LUFBSUBCarries out forward and back substitutionXXXXX
GEOM2DGeometrical operations (2D)XXX
GEOM3DGeometrical operations (3D)XXXX
GLRULESGauss-Legendre quadrature rulesXXX
GLT77-point Gaussian quadrature rule for triangleXX
GLT2525-point Gaussian quadrature rule for triangleXX
VGEOM2Verifies the geometry (2D)X
VGEOMAVerifies the geometry (axisym)XX
VGEOM3Verifies the geometry (3D)XX
VG2LCVerifies the use of the L2LC moduleX
L3ALCCCopy of L3ALC (to fake recursion)XX

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.

References

  1. 1. Wrobel LC, Aliabadi MH. The Boundary Element Method, Applications in Thermo-Fluids and Acoustics. Chichester, UK: John Wiley & Sons; 2002
  2. 2. Li ZC. Combined Methods for Elliptic Equations with Singularities, Interfaces and Infinities. Vol. 444. Boston, London: Kluwer Academic Publishers; 1998
  3. 3. Lindgren LE. From Weighted Residual Methods to Finite Element Methods. Available from: https://www.ltu.se/cms_fs/1.47275!/mwr_galerkin_fem.pdf [Accessed: 10 April 2019]
  4. 4. Gato C. Detonation-driven fracture in thin shell structures: Numerical studies. Applied Mathematical Modelling. 2010;34:3741-3753
  5. 5. Semwogerere T. Application analysis of the curved elements to a two-dimensional Laplace problem using the boundary element method. International Journal of Applied Physics and Mathematics. 2015;5:167-176
  6. 6. Seydou F, Seppanen T, Ramahi O. Numerical solution of 3D laplace and helmholtz equations for parallel scatterers. In: Proceedings of the Antennas and Propagation Society International Symposium; Washington, DC, USA: IEEE; 2005. pp. 6-9
  7. 7. Shahbazi M, Mansourzadeh S, Pishevar AR. Hydrodynamic analysis of autonomous underwater vehicle (AUV) flow through boundary element method and computing added-mass coefficients. International Journal of Artificial Intelligence and Mechatronics. 2015;3:212-217
  8. 8. Weber C. Electron Energy Loss Spectroscopy with Plasmonic Nanoparticles. Graz, Austria: Karl-Franzens-Universität Graz; 2011
  9. 9. Kirkup SM. DC capacitor simulation by the boundary element method. Communications in Numerical Methods in Engineering. 2007;23:855-869
  10. 10. Cunderlik R, Mikula K, Mikula K. A comparison of the variational solution to the Neumann geodetic boundary value problem with the geopotential model. Contributions to Geophysics and Geodesy. 2004;34:209-225
  11. 11. Sklyar O, Kueng A, Kranz C, Mizaikoff B, Lugstein A, Bertagnolli E, et al. Numerical simulation of scanning electrochemical microscopy experiments with frame-shaped integrated atomic force microscopy: SECM probes using the boundary element method. Analytical Chemistry. 2005;77:764-771
  12. 12. Nkurunziza D, Kakuba G, Mango JM, Rugeihyamu SE, Muyinda N. Boundary element method of modelling steady state groundwater flow. Applied Mathematical Sciences. 2014;8:8051-8078
  13. 13. Hashemi AR, Pishevar AR, Valipouri A, Pǎrǎu EI. Numerical and experimental study on the steady cone-jet mode of electro-centrifugal spinning. Physics of Fluids. 2018;30:114108
  14. 14. Kirkup SM, Henwood DJ. An empirical error analysis of the boundary element method applied to Laplace’s equation. Applied Mathematical Modelling. 1994;18:32-38
  15. 15. Hu H. An analysis of double Laplace equations on a concave domain. An Anal Double Laplace Equations a Concave Domain. 2018;6:304-314
  16. 16. Menin OH, Rolnik V. Relation between accuracy and computational time for boundary element method applied to Laplace equation. Journal of Computational Interdisciplinary Sciences. 2013;4:1-6
  17. 17. Fiala P, Rucz P. NiHu: An open source C++ BEM library. Advances in Engineering Software. 2014;75:101-112
  18. 18. Wieleba P, Sikora J. Open source BEM library. Advances in Engineering Software. 2009;40:564-569
  19. 19. Kumara PK. A MATLAB code for three dimensional linear elastostatics using constant boundary elements. International Journal of Advanced Engineering. 2012;2:9-20
  20. 20. Kirkup SM, Yazdani J. A gentle introduction to the boundary element method in Matlab/freemat. In: Proceedings of the WSEAS MAMECTIS; Corfu, Greece; 2008
  21. 21. Kirkup SM. Fortran codes for computing the discrete Helmholtz integral operators. Advances in Computational Mathematics. 1998;9:391-404
  22. 22. Kirkup SM. The Boundary Element Method in Acoustics. Hebden Bridge: Integrated Sound Software; 1998
  23. 23. Kirkup SM. The boundary element method in excel for teaching vector calculus and simulation. World Academy of Science, Engineering and Technology, International Science Index 144, International Journal of Social, Behavioral, Educational, Economic, Business and Industrial Engineering. 2019;12:1605-1613
  24. 24. Kirkup SM. The boundary element method in acoustics: A survey. Applied Sciences. 2019;9:1642-1690
  25. 25. Kirkup SM. Boundary Element Method. Available from: www.boundary-element-method.com [Accessed: 11 January 2018]
  26. 26. Shaukat K, Tahir F, Iqbal U, Amjad S. A comparative study of numerical analysis packages. International Journal of Computer Theory and Engineering. 2019;10:67-72
  27. 27. Warham AGP. The Helmholtz Integral Equation for a Thin Shell. London, UK: National Physical Laboratory; 1988
  28. 28. Kirkup SM. The boundary and shell element method. Applied Mathematical Modelling. 1994;18:418-422
  29. 29. Kirkup SM. The computational modelling of acoustic shields by the boundary and shell element method. Computers and Structures. 1991;40:1177-1183
  30. 30. Kirkup SM. Solution of discontinuous interior Helmholtz problems by the boundary and shell element method. Computer Methods in Applied Mechanics and Engineering. 1997;140:393-404
  31. 31. Kirkup SM. Solving the linear systems of equations in the generalized direct boundary element method. In: Proceedings of the 1st International Conference on Numerical Modelling in Engineering; Ghent, Belgium; 2018

Written By

Stephen Kirkup and Javad Yazdani

Submitted: 15 November 2018 Reviewed: 24 April 2019 Published: 06 May 2020