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

*A Pilot Fortran Software Library for the Solution of Laplace's Equation by the Boundary… DOI: http://dx.doi.org/10.5772/intechopen.86507*

$$a(\mathbf{p})\boldsymbol{\varrho}(\mathbf{p}) + \beta(\mathbf{p})\frac{\partial \boldsymbol{\varrho}}{\partial n\_p}(\mathbf{p}) = f(\mathbf{p})(\mathbf{p} \in \mathbb{S}).\tag{3}$$

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

$$\int\_{\mathcal{S}} \frac{\partial G(\mathbf{p}, \mathbf{q})}{\partial n\_q} \, \rho(\mathbf{q}) dS\_q + \frac{1}{2} \rho(\mathbf{q}) = \int\_{\mathcal{S}} G(\mathbf{p}, \mathbf{q}) \frac{\partial \rho(\mathbf{q})}{\partial n\_q} \, dS\_q(\mathbf{p} \in \mathcal{S}), \tag{4}$$

$$\int\_{\mathcal{S}} \frac{\partial G(\mathbf{p}, \mathbf{q})}{\partial \mathfrak{n}\_q} \, \varrho(\mathbf{q}) d\mathbf{S}\_q + \varrho(\mathbf{q}) = \int\_{\mathcal{S}} G(\mathbf{p}, \mathbf{q}) \frac{\partial \wp(\mathbf{q})}{\partial \mathfrak{n}\_q} \, dS\_{\mathfrak{q}}(\mathbf{p} \in D). \tag{5}$$

The terminology *<sup>∂</sup>* <sup>∗</sup> *<sup>∂</sup>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

$$G(\mathbf{p}, \mathbf{q}) = -\frac{1}{2\pi} \ln r \tag{6}$$

for two-dimensional Laplace problems, where *r* ¼ ∣*q* � *p*∣.

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 <sup>φ</sup>(*q*), *<sup>∂</sup><sup>φ</sup>*ð Þ*<sup>q</sup> <sup>∂</sup>nq* 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**

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

The simplest partial differential equation that is amenable to the BEM is

*∂*2 *φ*ð Þ *p ∂x*<sup>2</sup> *i*

∇2

steady-state electric potential, gravitation and groundwater flow [4–13].

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,

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

There have been a number of works on coding the boundary element method [17–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–24], this is about continuing with the development of a base library of methods and corresponding software. The

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

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

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

¼ 0 (1)

*φ* ¼ 0*:* (2)

∑ *N i*¼1

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

codes and guides can be found on the first author's website [25].

**2. The BEM and the 2D interior Laplace problem**

foundation for the 3D BEM development in later sections.

this work the direct boundary element method is developed.

*Numerical Modeling and Computer Simulation*

Laplace's equation:

ary element method [14–16].

to the library of numerical software [26].

**4**

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

$$\int\_{\Gamma} G(\mathfrak{p}, \mathfrak{q}) \zeta(\mathfrak{q}) d\mathbb{S}\_{\mathfrak{q}} = \mu(\mathfrak{p})(\mathfrak{p} \in \mathbb{S}) \tag{7}$$

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

$$\{L\zeta\}\_\Gamma(\mathbf{p}) = \mathfrak{\mu}(\mathbf{p}).\tag{8}$$

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:

$$\{M\zeta\}\_\Gamma(\mathfrak{p}) = \int\_\Gamma \frac{\partial G(\mathfrak{p}, \mathfrak{q})}{\partial n\_q} \zeta(\mathfrak{q}) d\mathfrak{S}\_q,\tag{9}$$

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

*A Pilot Fortran Software Library for the Solution of Laplace's Equation by the Boundary…*

ð Þ *<sup>p</sup> <sup>φ</sup><sup>j</sup>* ≈ ∑*<sup>n</sup>*

*p*, are the numerical values of definite integrals that together can be interpreted as

*<sup>p</sup>Si* � �*φ<sup>j</sup>* ≈ ∑*<sup>n</sup>*

relating the exact values of the boundary functions at the collocation points.

*MSS* þ 1 2 I � �*φ*^*<sup>S</sup>* <sup>¼</sup> *LSSv*

with the matrix components defined by ½ � *LSS ij* ¼ f g *Le <sup>Δ</sup>S*<sup>~</sup>

*Sj <sup>p</sup>Si* � �. The vectors *<sup>φ</sup>*^*<sup>S</sup>* and *<sup>v</sup>*

to) boundary data is known at the collocation points.

f g *Le <sup>Δ</sup>*<sup>~</sup>

*Sj <sup>p</sup>Di* � �*<sup>v</sup>*

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

*j*¼1

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

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

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

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

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 *pDi* in the domain *D*e . 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

f g *Me <sup>Δ</sup>*<sup>~</sup>

^*<sup>S</sup>* � *MDSφ*^*<sup>S</sup>*

*Sj <sup>p</sup>Di* � �.

^*<sup>j</sup>* � ∑ *n j*¼1

*φ*^*<sup>D</sup>* ¼ *LDSv*

*Sj <sup>p</sup>Di* � �*, M*½ � *DS ij* <sup>¼</sup> f g *Me <sup>Δ</sup>*<sup>~</sup>

f g *Le <sup>Δ</sup>S*<sup>~</sup>

*j*¼1

*j*

f g Le *<sup>Δ</sup>*<sup>~</sup> *Sj* ð Þ *p vj p*∈e*S*

*<sup>j</sup> <sup>p</sup>Si* � �*vj <sup>p</sup>Si* <sup>∈</sup>e*<sup>S</sup>*

ð Þ *p* , for example, for a specific point

� �*,* (16)

� �*,* (17)

^*S,* (18)

*Sj <sup>p</sup>Di* � �*φ*^*<sup>j</sup> <sup>p</sup>Di* <sup>∈</sup> *<sup>D</sup>*<sup>e</sup> � �*,* (20)

*,* (21)

*<sup>j</sup> <sup>p</sup>Si* � �*,*

^*<sup>S</sup>* are representative or approximate

*αiφ<sup>i</sup>* þ *βivi* ¼ *fi* for *i* ¼ 1*,* 2*,* …*n:* (19)

*Δ*~ *Sj*

*ΔS*~ *j*

rewrite Eq. (11) as the approximation:

*DOI: http://dx.doi.org/10.5772/intechopen.86507*

*M* þ 1 2 *I* � �*<sup>e</sup>* � �

where *e* is the unit function (*e* � 1). f g *Le <sup>Δ</sup>S*<sup>~</sup>

the discrete form of the *L* integral operator.

*M* þ 1 2 *I* � �*<sup>e</sup>* � �

∑ *n j*¼1

system of the form

½ � *MSS ij* ¼ f g *Me <sup>Δ</sup>*<sup>~</sup>

*<sup>φ</sup> <sup>p</sup>Di* � � <sup>¼</sup> <sup>∑</sup>

where ½ � *LDS ij* ¼ f g *Le <sup>Δ</sup>*<sup>~</sup>

written as

**7**

*n j*¼1

∑ *n j*¼1

$$\{\mathcal{M}^t \zeta\}\_\Gamma(\mathfrak{p}; \mathfrak{v}\_p) = \frac{\partial}{\partial v\_p} \int\_\Gamma G(\mathfrak{p}, \mathfrak{q}) \zeta(\mathfrak{q}) d\mathbb{S}\_\mathfrak{q},\tag{10}$$

$$\{N\zeta\}\_\Gamma(\mathfrak{p};\mathfrak{v}\_p) = \frac{\partial}{\partial v\_p} \int\_\Gamma \frac{\partial G(\mathfrak{p},\mathfrak{q})}{\partial n\_q} \zeta(\mathfrak{q}) d\mathbf{S}\_{\mathfrak{q}}.\tag{11}$$

where *v<sup>p</sup>* 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) \wp \right\}\_{\mathbb{S}} (\mathfrak{p}) = \{ Lv \}\_{\mathbb{S}} (\mathfrak{p}) (\mathfrak{p} \in \mathbb{S}), \tag{12}$$

$$\rho(\mathbf{p}) = \{L\boldsymbol{\nu}\}\_{\mathcal{S}} - \{M\boldsymbol{\rho}\}\_{\mathcal{S}}(\mathbf{p} \in D),\tag{13}$$

where *v q* � � <sup>¼</sup> *<sup>∂</sup><sup>φ</sup>*ð Þ*<sup>q</sup> <sup>∂</sup>nq* .

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

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

$$\mathcal{S} \approx \widetilde{\mathcal{S}} = \sum\_{j=1}^{n} \Delta \widetilde{\mathcal{S}}\_{j}. \tag{14}$$

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�<sup>1</sup> collocation method involves representing the boundary function by a constant on each panel:

$$
\rho(\mathbf{p}) \approx \rho\_{\mathbf{j}}, \\
v(\mathbf{p}) \approx v\_{\mathbf{j}}(\mathbf{p} \in \Delta \tilde{\mathbf{S}}\_{\mathbf{j}}).\tag{15}
$$

*A Pilot Fortran Software Library for the Solution of Laplace's Equation by the Boundary… DOI: http://dx.doi.org/10.5772/intechopen.86507*

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\_{\hat{\jmath}}^{\mathbb{S}}} (\mathbf{p}) \rho\_{\hat{\jmath}} \approx \sum\_{j=1}^{n} \left\{ \text{Le} \right\}\_{\Delta\_{\hat{\jmath}}^{\mathbb{S}}} (\mathbf{p}) v\_{\triangleright} \left( \mathbf{p} \in \widetilde{\mathcal{S}} \right), \tag{16}$$

where *e* is the unit function (*e* � 1). f g *Le <sup>Δ</sup>S*<sup>~</sup> *j* ð Þ *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

$$\sum\_{j=1}^{n} \left\{ \left( M + \frac{1}{2} I \right) e \right\}\_{\Delta\_{\hat{\mathbb{S}}}^{\mathbb{S}}} (\mathbf{p}\_{\hat{\mathbb{S}}i}) \rho\_{\hat{j}} \approx \sum\_{j=1}^{n} \left\{ L e \right\}\_{\Delta\_{\hat{\mathbb{S}}}^{\mathbb{S}}} (\mathbf{p}\_{\hat{\mathbb{S}}i}) v\_{\hat{j}} \left( \mathbf{p}\_{\hat{\mathbb{S}}i} \in \tilde{\mathbb{S}} \right), \tag{17}$$

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:

$$\left(\mathbf{M}\_{\rm SS} + \frac{1}{2}\mathbf{I}\right)\hat{\underline{\boldsymbol{\rho}}}\_{\rm S} = L\_{\rm SS}\hat{\underline{\boldsymbol{\nu}}}\_{\rm S} \tag{18}$$

with the matrix components defined by ½ � *LSS ij* ¼ f g *Le <sup>Δ</sup>S*<sup>~</sup> *<sup>j</sup> <sup>p</sup>Si* � �*,* ½ � *MSS ij* ¼ f g *Me <sup>Δ</sup>*<sup>~</sup> *Sj <sup>p</sup>Si* � �. The vectors *<sup>φ</sup>*^*<sup>S</sup>* and *<sup>v</sup>* ^*<sup>S</sup>* 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):

$$
\beta\_i \rho\_i + \beta\_i v\_i = f\_i \text{ for } i = 1, 2, \dots \\ n. \tag{19}
$$

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:

$$\log \left( \mathbf{p}\_{Di} \right) = \sum\_{j=1}^{n} \left\{ L\mathcal{c} \right\}\_{\Delta \tilde{\mathbb{S}}\_{j}} \left( \mathbf{p}\_{Di} \right) \hat{\boldsymbol{\upsilon}}\_{j} - \sum\_{j=1}^{n} \left\{ M\mathcal{c} \right\}\_{\Delta \tilde{\mathbb{S}}\_{j}} \left( \mathbf{p}\_{Di} \right) \hat{\boldsymbol{\rho}}\_{j} \left( \mathbf{p}\_{Di} \in \tilde{\mathbf{D}} \right), \tag{20}$$

for each point *pDi* in the domain *D*e . 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

$$
\underline{\hat{\varrho}}\_{D} = L\_{DS} \underline{\hat{\nu}}\_{S} - M\_{DS} \underline{\hat{\varrho}}\_{S} \tag{21}
$$

$$\text{where } [L\_{DS}]\_{\vec{\boldsymbol{\eta}}} = \{ L\boldsymbol{\epsilon} \}\_{\Delta\_{\vec{\boldsymbol{\eta}}}^{\vec{\boldsymbol{\epsilon}}}} (\boldsymbol{\mathsf{p}}\_{Di})\_{\boldsymbol{\epsilon}} \, [\boldsymbol{M}\_{DS}]\_{\vec{\boldsymbol{\eta}}} = \{ \boldsymbol{M}\boldsymbol{\epsilon} \}\_{\Delta\_{\vec{\boldsymbol{\eta}}}^{\vec{\boldsymbol{\epsilon}}}} (\boldsymbol{\mathsf{p}}\_{Di})\_{\boldsymbol{\epsilon}}.$$

or a part of the boundary. The other three important Laplace integral operators are

*∂G p; q* � � *∂nq*

*ζ q*

*G p; q* � �*ζ q*

*ζ q*

*φ*ð Þ¼ *p* f g *Lv <sup>S</sup>* � f g *Mφ <sup>S</sup>*ð Þ *p*∈ *D ,* (13)

ð Þ¼ *p* f g *Lv <sup>S</sup>*ð Þ *p* ð Þ *p*∈*S ,* (12)

e*j:* (14)

*:* (15)

*∂G p; q* � � *∂nq*

� �*dSq,* (9)

� �*dSq,* (10)

� �*dSq,* (11)

ð Г

> *∂vp* ð Г

*∂vp* ð Г

where *v<sup>p</sup>* is any unit vector. In operator notation of the previous subsection, the

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�<sup>1</sup> collocation in this text since it is derived by approximating the boundary functions by a constant on each panel. In this subsection the C�<sup>1</sup> colloca-

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

*n j*¼1 Δ*S*

*S* ≈ e*S* ¼ ∑

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

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�<sup>1</sup> collocation method involves

*, v*ð Þ *p* ≈ *vj p*∈*Δ*e*Sj*

� �

representing the boundary function by a constant on each panel:

*φ*ð Þ *p* ≈ *φ<sup>j</sup>*

� � <sup>¼</sup> *<sup>∂</sup>*

� � <sup>¼</sup> *<sup>∂</sup>*

integral equation formulation (3) can be written in the following form:

*φ*

*S*

f g *Mζ* <sup>Г</sup>ð Þ¼ *p*

*Mt* f g*<sup>ζ</sup>* <sup>Г</sup> *<sup>p</sup>*; *<sup>v</sup><sup>p</sup>*

f g *Nζ* <sup>Г</sup> *p*; *v<sup>p</sup>*

*M* þ 1 2 *I* � �

� �

defined as follows:

*Numerical Modeling and Computer Simulation*

where *v q*

� � <sup>¼</sup> *<sup>∂</sup><sup>φ</sup>*ð Þ*<sup>q</sup> <sup>∂</sup>nq* .

tion method is briefly outlined.

tions of the form introduced earlier.

**6**

**2.3 Direct boundary element method**
