**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 *Γ* and denote them by '+' and '�'. We then introduce the

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

**Figure 2.** *A hemispherical shell.*

The interior test problem is that of a unit sphere approximated by 36 triangular

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

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

The exterior test problem is that of a unit sphere approximated by 36 triangular

*φ* ¼ *x* þ *y* þ *z:* (27)

*,* (28)

panels. The exact solution that is applied as a Dirichlet boundary condition is

panels, as in the previous test. The exact solution that is applied as a Dirichlet

*<sup>φ</sup>* <sup>¼</sup> <sup>1</sup> *r*

where *r* is the distance from 0ð Þ *;* 0*;* 0*:*5 *:* The results at four exterior points are

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

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

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

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

modelling of the steady-state electric field in a capacitor in Kirkup [9].

**4.1 Integral equations and boundary element equations for thin shells**

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

*The results from the three-dimensional interior problem.*

*Numerical Modeling and Computer Simulation*

*The results from the three-dimensional exterior problem.*

boundary condition is

**Table 4.**

**Table 5.**

given in **Table 5**.

**12**

or discontinuities in the potential field.

quantities of difference and average of the potential and its normal derivative across the surface:

$$\delta(\mathbf{p}) = \varrho(\mathbf{p}\_+) - \varrho(\mathbf{p}\_+) (\mathbf{p} \in \Gamma),\tag{29}$$

$$\Phi(\mathbf{p}) = \frac{1}{2} \left( \wp(\mathbf{p}\_+) + \wp(\mathbf{p}\_+) \right) (\mathbf{p} \in \varGamma), \tag{30}$$

$$\nu(\mathfrak{p}) = \nu(\mathfrak{p}\_+) + \nu(\mathfrak{p}\_+)(\mathfrak{p} \in \Gamma),\tag{31}$$

$$V(\mathbf{p}) = \frac{1}{2} \left( v(\mathbf{p}\_+) - v(\mathbf{p}\_+) \right) (\mathbf{p} \in \varGamma). \tag{32}$$

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

$$\rho(\mathbf{p}) = \{\mathbf{M}\delta\}\_{\varGamma}(\mathbf{p}) - \{L\nu\}\_{\varGamma}(\mathbf{p})(\mathbf{p} \in E),\tag{33}$$

$$\Phi(\mathbf{p}) = \{M\delta\}\_{\Gamma}(\mathbf{p}) - \{L\nu\}\_{\Gamma}(\mathbf{p})(\mathbf{p} \in \Gamma),\tag{34}$$

$$V(\mathbf{p}) = \{N\delta\}\_{\varGamma}(\mathbf{p}) - \{M^t\nu\}\_{\varGamma}(\mathbf{p})(\mathbf{p} \in \varGamma). \tag{35}$$

The boundary condition may be expressed in the following form:

$$a(\mathbf{p})\delta(\mathbf{p}) + \beta(\mathbf{p})\nu(\mathbf{p}) = f(\mathbf{p})(\mathbf{p} \in \Gamma),\tag{36}$$

$$A(\mathbf{p})\Phi(\mathbf{p}) + \beta(\mathbf{p})V(\mathbf{p}) = F(\mathbf{p})(\mathbf{p} \in \Gamma). \tag{37}$$

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

$$
\underline{\hat{\varrho}}\_{E} = \mathbf{M}\_{EI}\hat{\underline{\delta}}\_{I} - L\_{EI}\hat{\underline{\nu}}\_{I^{\gamma}} \tag{38}
$$

$$
\underline{\hat{\Phi}}\_{\Gamma} = \mathbf{M}\_{\Gamma\Gamma} \hat{\underline{\hat{\delta}}}\_{\Gamma} - L\_{\Gamma\Gamma} \hat{\underline{\hat{\nu}}}\_{\Gamma^{\flat}} \tag{39}
$$

$$
\widehat{\underline{V}}\_{\Gamma} = \mathsf{N}\_{\Gamma\Gamma} \widehat{\underline{\delta}}\_{\Gamma} - \mathsf{M}\_{\Gamma\Gamma}^{t} \widehat{\underline{\nu}}\_{\Gamma}. \tag{40}
$$

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


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, *<sup>Φ</sup>*^ *<sup>Γ</sup>*; VELDIF, *<sup>ν</sup>*^*Γ*; VELAV, *<sup>V</sup>*^ *<sup>Γ</sup>*; and PPHI, *<sup>φ</sup>*^*E*.

**5. Conclusions**

*The results from the three-dimensional shell problem.*

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

**Table 7.**

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

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

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

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

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

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 n*<sup>2</sup> ð Þ time and memory. Solving the linear system by a direct method, like LU factorisation used in this work, takes *O n*<sup>3</sup> ð Þ 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

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

library of codes and the way they are linked are set out in Appendix.

this work the validation is developed through the VGEOM\* modules.

\*INPHI and \*INVEL parameters), but these have not been tested.

matrices may also become an issue.

**A. Appendix**

**15**

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