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 *L* and *M* operators are required.

In general L2LC simply implements a Gaussian quadrature rule in order to determine the integral, using a higher-order rule when point *p* is close to the panel *ΔS* e*j* . 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:

final node with the first node to complete the boundary. The boundary conditions

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

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

The exact solution is *φ*ð Þ¼ *p* 10 þ 10*x*; 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

A set of nodal coordinates and each panel is determined through the two nodal

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—

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-

be interpreted as temperatures in a steady-state heat conduction problem.

shown in **Figure 1** are then applied.

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

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

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

indices for the endpoints of the panel.

**3. The BEM and 3D Laplace problems**

dimensional Laplace problems is

**9**

shown in **Table 1**.

**Figure 1.**

**Table 1.**

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


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*, *β<sup>i</sup>* 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.

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

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



#### **Table 1.**

**2.4 LIBEM2 and L2LC modules, test problem and results**

tant component module.

ð Þ *p* , f g *Me <sup>Δ</sup>*<sup>~</sup>

*Sj*

*Numerical Modeling and Computer Simulation*

tion, only *L* and *M* operators are required.

\* MAXPOINTS,NPOINT,POINTS, \* SALPHA,SBETA,SF,SINPHI,PINPHI,

\* L\_SS,M\_SSPMHALFI,L\_PS,M\_PS, \* PERM,XORY,C,workspace)

and (0.5, 0.5), and these are also illustrated in the figure.

\* LSOL,LVALID,TOLGEOM,

\* SPHI,SVEL,PPHI,

values in Eq. (19) for each panel.

**8**

ð Þ *<sup>p</sup> , M<sup>t</sup>* f g*<sup>e</sup> <sup>Δ</sup>*<sup>~</sup>

*Sj*

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

following form:

following form:

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

**L2LC.** The L2LC module computes the discrete Laplace operators for twodimensional problems. In the notation of this article, the routine computes

ð Þ *p* and f g *Ne <sup>Δ</sup>*<sup>~</sup>

the domain of integration and *p* is any point. The call to the subroutine has the

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

where P is point *p*; VECP is a unit directional vector that passes through *p*; QA

In general L2LC simply implements a Gaussian quadrature rule in order to deter-

mine the integral, using a higher-order rule when point *p* is close to the panel *ΔS*

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

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

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*, *β<sup>i</sup>* and *fi*

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

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

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

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

*Sj*

ð Þ *p* , where *ΔS*

e*<sup>j</sup>* is the panel that is

e*j* .

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

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 þ 10*x*; 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**.

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 threedimensional Laplace problems is

*Numerical Modeling and Computer Simulation*

$$G(p,q) = \frac{1}{4\pi r},\tag{22}$$

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

$$\left\{ \left( M - \frac{1}{2} I \right) \rho \right\}\_{\mathcal{S}} (\mathfrak{p}) = \{ Lv \}\_{\mathcal{S}} (\mathfrak{p}) (\mathfrak{p} \in \mathcal{S}), \tag{23}$$

$$
\rho(\mathfrak{p}) = \{M\mathfrak{q}\}\_S - \{L\upsilon\}\_S (\mathfrak{p} \in E). \tag{24}
$$

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

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

**Index Point Exact Computed (4 d.p.)** (0.0, 0.0) 0.0 �0.0013 (0.0, 0.5) �0.5 �0.4995 (0.0, �0.5) �0.5 �0.4995 (0.5, 0.0) 0.25 0.2477

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

**3.2 LBEM3 and L3LC modules, test problems and results**

panels. The L3LC subroutine is called as follows:

three-dimensional domain and is called as follows:

\* SALPHA,SBETA,SF,SINPHI,PINPHI,

\* L\_SS,M\_SSPMHALFI,L\_PS,M\_PS,

\* LSOL,LVALID,TOLGEOM,

\* LINTERIOR,MAXPOINTS,NPOINT,POINTS,

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

results are given in **Table 3**.

*The results from the axisymmetric interior problem.*

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

*The results from the axisymmetric exterior problem.*

**Table 2.**

**Table 3.**

vertices of the triangular panel.

\* SPHI,SVEL,PPHI,

**11**

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

The LBEM3 and L3LC subroutines implement the boundary element method

axisymmetric codes, the component module L3LC computes the integrals over the

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

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

The LBEM3 module solves Laplace's equation in a general interior or exterior

As with LIBEM2 and LBEMA, NODES and PANELS define the boundary. How-

ever, 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.

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

for general three-dimensional problems. As with the two-dimensional and

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

*,* (26)

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.
