The Singular Function Boundary Integral Method for the 2-D and 3-D Laplace Equation Problems in Mechanics, with a Boundary Singularity

In this study the Singular Function Boundary Integral Method (SFBIM) is implemented in the case of a planar elliptic boundary value problem in Mechanics, with a point boundary singularity. The method is also extended in the case of a typical problem of Solid Mechanics, concerning the Laplace equation problem in three dimensions, defined in a domain with a straight edge singularity on the surface boundary. In both the 2-D and 3-D cases, the general solution of the Laplace equation is approximated by the leading terms (which contain the singular functions) of the local asymptotic solution expansion. The singular functions are used to weight the governing equation in the Galerkin sense. For the 2-D Laplacian model problem of this study, which is defined over a domain with a re-entrant corner, the resulting discretized equations are reduced to boundary integrals by means of Green’s second identity. For the 3-D model problem of this work, the volume integrals of the discretized equations are reduced to surface integrals by implementing Gauss’ divergence theorem. The Dirichlet boundary conditions are then weakly enforced by means of Lagrange multipliers. The values of the latter are calculated together with the singular coefficients, in the 2-D case or the Edge Flux Intensity Functions (EFIFs), in the 3-D model problem, which appear in the local solution expansion. For the planar problem, the numerical results are favorably compared with the analytic solution. Especially for the extension of the method in three dimensions, the preliminary numerical results compare favorably with available postprocessed finite element results.


Introduction
In the past few decades, many numerical methods have been proposed for the solution of 2-D elliptic boundary value problems in Mechanics, with a boundary singularity. The lack of adequate accuracy and poor convergence, were some of the difficulties which appeared in the attempt of many researchers to solve this kind of problems. Remedies used were special mesh-refinement schemes, multigrid methods, singular elements, p/hp finite elements and many other techniques (e.g. [1-3, 5, 24, 25]). An extensive survey of the treatment of singularities in elliptic 2-D boundary value problems is provided in [8]. Of course, in the absence of a boundary singularity, there are many other efficient numerical techniques available in the literature (e.g. [1,2,4,22,29]) to be used for the solution of problems in Mechanics.
For the 2-D problems with a boundary singularity there are more sophisticated techniques in the literature, than those mentioned above, which incorporate the form of the local asymptotic expansion. For example, in [18] a singular finite element method is developed for Stoke' s flow problems, in which special elements are employed, in the neighborhood of the singularity, with which the radial form of the local expansion is utilized, in order to resolve the convergence difficulties and improve the accuracy of the global solution.
The interest in such two-dimensional Laplace equation problems with a boundary singularity [20], is motivated by In the above asymptotic expansion, µ j are the eigenvalues and f j (θ) are the eigenfuntions of the problem and they are known. The polar coordinates (r, θ) have as a center the singular point. Singular coefficients β j are the primary unknowns in this kind of problems. They are also known as generalized stress intensity factors (GSIFs). These are determined by the boundary conditions on the boundary parts which are away from the singularity. They are of significant importance in many applications in the area of Fracture Mechanics [16]. Thus, after the calculation of the leading singular coefficients β j , the approximation of Φ(x,y) is known (by converting into Cartesian coordinates) and the approximate stress state, at any position, is found from the known expressions: σ xx =∂ 2 Φ/∂y 2 , σ yy =∂ 2 Φ/∂x 2 and σ xy =− ∂ 2 Φ/∂x∂y.
In the past two decades, the Singular Function Boundary Integral Method (SFBIM) was developed [9-14, 17, 19], which belongs to the general group of collocation methods [26] and with which the unknown singular coefficients are calculated directly. Thus, the method gives directly the approximation of the Airy stress function (for solids) or the stream function (for fluids) in planar Laplacian and biharmonic problems. Basically, for the implementation of the SFBIM to solve 2-D problems with a boundary singularity, it is required that the solution is approximated by the leading terms of expansion (1). The Dirichlet conditions, on the boundary parts, which are away from the singularity, are enforced by means of Lagrange multipliers. The method has been tested on standard elliptic problems exhibiting exponential convergence and high accuracy with respect to the number of singular functions. This behavior of the method has been theoretically proved in [28]. The results obtained in the 2-D case have encouraged the extension of the method in 3-D problems with a straight-edge singularity [6,13].
The 2-D model problem, treated in this study, is a Laplacian problem defined over a disk, with a boundary ∂Ω on which a quarter of a circle is missing and thus a boundary singularity at the centre, is created. This could be a heat transfer problem in two dimensions, but in this study it is a plane stress problem of a perfectly elastic solid, on which the interest is to find the unknown singular coefficients β j , of expansion (1). As for the 3-D case, we consider a Laplacian boundary value problem in a 3-D domain. This is a model for an elastic cylindrical body, having a typical V-notch, made of an isotropic material which obeys to Hooke's Law and which is subjected to certain physical conditions. The form of the local solution depends on the geometry in the neighbourhood of the straight-edge singularity and on the boundary conditions of the boundary parts which share the singularity. Also, the local solution is characterized by the presence of certain eigenpairs (arising from the 2-D problem) and the socalled Edge Flux Intensity Functions (EFIFs), which are the primary unknowns in the 3-D case [30]. The interest in such a problem is motivated by the need to compute generalized stress intensity functions for the V-notched solids loaded by static loads. For the solution of this class of 3-D problems, few methods have been proposed so far, such as the J-integral method [21], the B-and H-integral methods [23], and more recently the methods by Costabel et al. [7] and Yosibash et al. [32][33][34][35] and Zaltzman [36], in which the EFIFs are computed by means of a post-processing procedure in a p-version finite element scheme.
Based on [7] and [30], the solution u of the 3-D Laplace equation, in a domain with an edge singularity, may be written in terms of cylindrical coordinates (r, θ, z) as In the above expansion, α κ ∈ R and φ κ (θ,α κ ) is the known eigenpair, of the two-dimensional Laplace operator. The functions φ κ (θ,α κ ) are analytic in θ. The functions Α (ακ) (z) are the EFIFs. These are functions of z and for a large class of problems they are the primary unknowns. In Solid Mechanics [16], function u represents the first stress invariant (i.e. u=σ ii ).
In the SFBIM, the solution is approximated by the leading terms of the asymptotic expansion (1) for the 2-D case, or the asymptotic expansion (2) for the 3-D case. The leading terms of these expansions are also used to discretize the governing differential equation in the Galerkin sense. The discretized equations are reduced to boundary integrals by means of the Green's second identity, for the 2-D case, or the Gauss divergence theorem, for the 3-D case. A particular feature of the SFBIM is that the Dirichlet conditions are weakly enforced by means of Lagrange multipliers. In twodimensional problems the coefficients of the asymptotic expansion are constants; these are calculated directly by the SFBIM. As we have already mentioned above, in three dimensions the EFIFs are functions of the axial direction. In the present approach, these are approximated by polynomials the coefficients of which are primary unknowns and are calculated directly by the method.
The outline of the rest of this paper is as follows: in Section 2 we present both a 2-D Laplace equation problem, with a point boundary singularity and a 3-D Laplace equation problem with a straight-edge singularity. The asymptotic local solution expansion is also given for each one of these problems. In Section 3 the formulation of the SFBIM is presented both in two and three dimensions. Numerical results are given in Section 4 together with comparisons between the results obtained with the SFBIM and the exact solution. Especially for the extension of the method in three dimensions, the numerical results obtained are compared Pure and Applied Mathematics Journal 2016; 5 (6): 192-204 194 with post-processed finite element results available in the literature [30]. Finally, conclusions derived from this study are summarized in Section 5.

Governing Equation and Local
Solution in 2-D and 3-D

The Planar Model Problem of the Laplace Equation
The 2-D model problem is a problem of a Laplace equation defined over a domain Ω. This domain is a very thin disk with radius R and with a re-entrant corner on its boundary, which has its tip on the center O of the disk ( Figure 1). Thus, at point O there is a boundary singularity. Boundary parts OA and OB, which share the singularity, intersect vertically at O and each one has the Neumann condition: ∂Φ/∂n=0. Circular part AB has the Dirichlet condition: Φ=f(θ). The material of the disk is isotropic and linearly elastic. Thus, basically, with the specific geometry and boundary conditions of this problem we need to find an expression for the Airy stress function Φ(r,θ), which is the solution of the Laplace equation. The values of Young's modulus and of Poisson's ratio are not necessary for the solution of the problem. It must be mentioned that this model problem belongs to the general class of problems of twodimensional plates with geometric singularities on the boundary. The planar model problem of the Laplace equation, illustrated by Figure 1 is the following: with the boundary conditions For this specific model problem the Dirichlet condition on the boundary S C has the following trigonometric expression: where R is the radius of boundary S C .
Using the separation of variables approach and the boundary conditions on S A and S B it can be easily found that the solution around the singular point is found to be Φ(r,θ)=βr µ j cos(µ j θ). Furthermore, by employing the principle of superposition we obtain the local solution expansion, which is expressed in terms of polar coordinates (r,θ), centered at O: Also, by considering the 2π-periodicity property of the trigonometric functions, singular coefficients β j can be calculated as Fourier coefficients by using Euler's formula. So, the exact solution of the 2-D model problem, defined by (3) and (4), is found to be Using the maximum principle theorem in its weak form (e.g. [15]) it can be shown that the above solution of the 2-D problem is unique.

The 3-D Model Problem of the Laplace Equation
The geometry and the boundary conditions of the 3-D model problem of the Laplace equation, are illustrated in Figure 2. The domain Ω is bounded by a cylindrical surface S C around the z-axis, a circular sector S D lying on the xyplane, a flat circular sector S E perpendicular to the z-axis and at a distance L from S D (here L=2) and two flat boundaries S A and S B (one parallel to the xz plane and the other parallel to the yz plane) intersecting vertically on the z-axis and thus creating a V-notch, which is the straight-edge singularity of this problem.
This solid body is made of an isotropic linearly elastic material and has no voids inside. For the specific geometry and boundary conditions of this model problem, the solution u, which, as we have already mentioned, normally represents the first stress invariant, is independent of any combination of values of the Young modulus E and the Poisson ratio ν. As we know, these are important parameters for the calculation of stresses and strains in linearly elastic materials which obey Equation Problems in Mechanics, with a Boundary Singularity to Hooke's law.

Fig. 2. Geometry and boundary conditions of the 3-D Laplacian model problem.
This problem is a version of a more general problem which was first suggested by Yosibash et al. [30] and is as follows: Find u such that with the following boundary conditions 0, where Naturally, the Laplace equation (7) is expressed in cylindrical co-ordinates, which are more convenient for both the analytical and numerical treatment of the problem.
In [30] it was shown that the solution u, of the threedimensional general Laplace equation problem with a governing equation ∇ 2 3-D u = 0 in Ω, such that the boundary conditions (8) hold, is obtained by augmenting the 2-D solution of the Laplace equation. Thus, according to [30], the solution in the neighborhood of the edge singularity, is given by where (r,θ,z)∈Ω are cylindrical coordinates, with the z-axis along the boundary straight-edge singularity. Functions A (ak) (z) are the EFIFs, which were described earlier. In this work they are chosen to be polynomials of degree N p and with unknown coefficients a k,j as follows: The constants c ki are given by Here, ω is the external angle defined by the flat boundaries ODEA and ODCB as shown in Figure 2; in the present work ω = 3π/2. Moreover, for the model problem of Figure 2, the first two EFIFs are known and they are of the form We must emphasize, at this stage, that the model problem of the 3-D case is selected in a way that the exact solution is known and given by a (finite) sum. Investigation of the uniqueness of this solution is similar with that of the 2-D case. Furthermore, since the only two non-zero EFIFs are polynomials, the accuracy of the EFIFs, computed by our method, can be measured by simply comparing the polynomial coefficients of the true and approximate EFIFs. In general, however, the true EFIFs are not always polynomials and the accuracy of the EFIFs computed by the numerical method, would need to be measured using some appropriate function norm.

The Singular Function Boundary
Integral Method

The SFBIM in the 2-D Case
The first step of the method, in both the planar and the three-dimensional problems, is the approximation of the local solution expansion with its leading terms. Thus, for the 2-D problem of Figure 1, expansion (5) is approximated by its first N a terms: The above series can also be written, more simply, as where V j =r µj cos(µ j θ). Functions V j are called singular functions and from now on, they are going to play an important role in the development of the formulation of the method. Note that these functions satisfy the governing equation and the boundary conditions along S A and S B . Next step is to weight the governing equation by the singular functions V i , in the Galerkin sense. Thus, we obtain the first set of discretized linear equations: So, we have our first set of N a discretized equations. Next, we apply Green's second identity by considering that both the approximation of Φ and the singular functions V i satisfy the Laplace equation and we obtain ( ) The integrand of the first integral in (15) is equal to zero along boundary parts S A and S B , because there it is ∂Φ/∂n=0 and ∂V i /∂n=0. Thus, the system of discretized equations (15) is simplified further to 0, 1, 2, , For boundary part S C the Dirichlet boundary condition is imposed by employing the Langrange multipliers function which is expanded in terms of quadratic basis functions M j as follows: In the above expansion, λ j are the so-called Lagrange multipliers which are auxiliary parameters and they are additional unknowns in our problem. As we will see, they are calculated together with singular coefficients β j after solving the complete system of equations. The quadratic basis functions M j , in expansion (17), are defined as follows: For j odd and for j even In expressions (18) and (19) h s is the mesh-width, which is constant in this approach.
The normal derivative ∂V i /∂n, which appears in (16), is expressed as follows ( ) sin , cos , where ñ is the normal unit vector on boundary S C (see Figure  1). Furthermore, using (20) the normal derivative of V i is given by which is an easy expression to be used in calculations. Before we proceed further, it is more convenient to change, in all the above formulation, variable s with variable θ by considering that θ=s/R ( Figure 1). Now, there is one Equation Problems in Mechanics, with a Boundary Singularity more step before we have our complete system of linear equations. So, according to the SFBIM the Dirichlet condition on S C must be weighted with the basis functions M i . Thus, finally, the discretized linear equations (16) take the form of the following system of N a +N λ discretized equations: The integrands in (22) are not singular and integrations take place away from the singularity. If we substitute the approximation of Φ and function λ in (22), with their expressions in (13) and (17), respectively, then the system of equations (22) can further be written as The above system of linear equations can also be expressed in matrix form as follows (24) or K.b=c (25) where K is called stiffness matrix. It is a symmetric matrix but it becomes singular when N a < N λ .

The SFBIM in the 3-D Case
In the SFBIM the solution of the problem is approximated by the leading terms of the local solution expansion given by (9). Thus, in (9) we employ the first N a terms and we substitute the EFIFs with their polynomial expression given by (10). Then the approximate solution is written as follows: where N (which also appears in (9)) is an additional parameter that allows us to ensure that the approximation of u satisfies the 3-D Laplace equation, by selecting it according to the restriction N p <2N+1 (see Appendix). We note that, in principle, N could be taken to be infinity since after all, the sum would terminate after a finite number of terms due to the fact that we are differentiating a polynomial of degree N p . Following the notation used in previous applications of the SFBIM (e.g. [9][10][11][12][13][14]), the above expansion is written as follows: In the above expression the coefficients a k,j will be referred to as singular coefficients and they are the coefficients of the polynomial expression (10) of the EFIFs. The functions W k (j) are the singular functions (the same name was adopted for the 2-D approach) and in this case they have the form It is easy to verify that ∇ 2 3-D (W k (j) )=0 in domain Ω and that the W k (j) 's satisfy exactly the boundary conditions on S A and S B (see Appendix).
As in the 2-D case, we weight the governing equation by the singular functions W k (j) in the Galerkin sense. This gives the first N a (N p +1) discretized equations: where ∂Ω=S A ∪S B ∪S C ∪S D ∪S E . Since the singular functions W k (j) satisfy exactly the boundary conditions on S A and S B and considering the boundary conditions (8), the boundary integral in (30) is identically zero along S A and S B . The discretized equations (30) then become As in the 2-D approach, the Dirichlet condition on the cylindrical boundary S C , which is away from the singularity, is imposed by means of a Lagrange multiplier function λ, which is expanded in terms of basis functions M i (θ,z): where N λ is the number of the discrete Lagrange multipliers λ (i) on S C . The nodal values of λ appear as additional unknowns in the problem. The additional N λ required equations are then obtained by weighting the Dirichlet boundary condition on S C by the bilinear basis functions M i (θ,z) in the Galerkin sense. Thus, in the end, we obtain the following linear system of (N p +1)N a +N λ discretized equations: Integration in the neighborhood of the point singularities O and D must be avoided. Thus, the EFIFs should be of degree equal to one or zero, for this problem (i.e. N p =1 or N p =0), in order to force the surface integrals on S C and S D to vanish. Therefore, by choosing a value of N p less or equal to 1 and by substituting the approximation of u and of the Lagrange multiplier function λ, with their expressions in (27) and (32), respectively, equations (33) and (34) . , , 2 It should be noted that, as in the 2-D case, the integrands in the above equations are non-singular and that all integrations are carried out far from the boundaries SA and SB and the points O and D causing the edge singularity. The surface integrals in (35) and (36) Here, L z (l) and L θ (l) are the dimensions of each element (l) along directions z and θ, respectively and N E is the number of boundary elements employed on the surface S C .
The system of discretized linear equations, (35) and (36), can be written in block form as follows: where the vector A contains the unknown coefficients a k,j of the approximation of the EFIFs and the vector Λ contains the unknown (discrete) Lagrange multipliers λ (i) . Clearly, the stiffness matrix K is symmetric. As in the case of 2-D problems, the number (Ν p +1)N a of the unknown coefficients a k,j should be greater than or equal to the number of Lagrange multipliers N λ , since the matrix Κ becomes singular when (Ν p +1)N a < N λ .

Numerical Results for the 2-D Problem
Expansion (17) indicates clearly that the Lagrange multiplier function λ, which is used to impose the Dirichlet boundary condition on S C , has specific numerical values λ (j) at positions j which are chosen according to a selected subdivision of S C into boundary elements. For the 2-D model problem (Figure 1), boundary S C is subdivided into N E elements of equal size. Therefore, the number of Lagrange multipliers is N λ =2N E +1.
The integrals in (23) are calculated numerically by subdividing each quadratic element into 10 subintervals and using a 15-point Gauss-Legendre quadrature over each subinterval. Since the stiffness matrix is symmetric only the elements which are on or above its principal diagonal are calculated. Now, in order to obtain the optimum combination of N a and N λ , several series of runs have been made. Previous applications of the method in 2-D problems [9][10][11][12][13][14], indicated that N λ should be large enough. This means that we must have an adequate number of elements on the boundary which will also help in having a better accuracy in numerical integration. However, N λ should be smaller than or equal to N a in order to avoid ill-conditioning of the stiffness matrix. In general, very high values of N a are also avoided because although double precision is used in the codes, the computer accuracy cannot handle the contributions of the higher-order singular functions, which become very small for r < 1 or very Equation Problems in Mechanics, with a Boundary Singularity large for r > 1. For the specific 2-D problem, the radius of the disk was chosen to be R=1. Hence, we do not have the phenomena described above. Thus, N λ was varied from 3 to 15 and N a was varied from 13 to 30. After performing calculations within this range of combinations of values for N a and N λ , it was observed that convergence occurs when N a =N λ =13 (optimum choice). Table 1 presents the convergence in the values of the four leading coefficients, with respect to the number of Lagrange multipliers N λ and for N a =13. As in previous implementations of the method [9][10][11][12][13][14], one may observe that the values of singular coefficients converge rapidly with N λ . In fact, in [28] a theoretical analysis of the method proved algebraic convergence in N λ . Also, very accurate estimates are obtained. Clearly, since the exact values of the leading singular coefficients of expansion (5) are β 1 =0, β 2 =2/5=0.40000000, β 3 =4/7≈0.571428571428571, β 4 =0 and β j =0, for j≥5, we can see, from Table 1, how fast convergence is achieved at N λ =13.
The values of the leading singular coefficients calculated for N λ =13 and various values of N a are shown in Table 2. Extremely fast convergence with respect to N a is observed (much faster than in previous implementations) and extremely accurate estimates are obtained. In fact, the leading singular coefficients converge from the very beginning. For N a >13 the solution starts to deteriorate. This is because the size of the system increases and thus the number of numerical calculations is increased together with the numerical errors which, of course, are not avoidable at the calculation of the matrix elements or at the numerical inversion of the stiffness matrix of the model problem. For the optimum combination N a =N λ =13, the plot of λ as a function of variable θ, is shown in Figure 3. As expected, it has the form of a cosine trigonometric function because on S C we have ∂V i /∂r=µ i r µi-1 cos(µ i θ). In this graph the curve is smooth and free from oscillations.
In Table 3 the converged values of the singular coefficients, calculated for the optimal choices N a =N λ =13 and R=1, are presented. The values of this table indicate clearly that the contribution of the higher order terms, vanishes immediately for this 2-D model problem and that the approximate values are, practically, equal to the exact values. Figure 4 shows the plot of the errors which appear in the calculated values of the leading singular coefficients β 2 and β 3 , as N λ varies, when N a =13. Here, the error is defined as the absolute value of the difference between the exact and the approximate solution of coefficients β j : This is not the only definition of error. Any suitable norm can also be used to define the error in numerical calculations. However, in the present problem the absolute value of the difference between the exact value and the approximate value of β j were much adequate in our present research. Both curves in Figure 4 indicate an algebraic convergence with N λ . This behavior will be explained theoretically below. The above behavior of the present technique was also demonstrated in previous implementations of the method, in both 2-D and 3-D problems [9][10][11][12][13][14]. According to reference [28], if λ ∈ H k (∂Ω AB ), for some k≥1 and λ h is the approximation of the Lagrange multiplier function with h being always the mesh-width, then there exist positive constants C and γ∈(0,1), independent of N α and h such that where m=min{k, p+1} and H k (∂Ω AB ), k∈N, is the usual Sobolev space on boundary part AB, which contains functions that have k generalized derivatives in the space of the square integrable functions L 2 (Ω). Also, in [28] it was shown that Inequalities (39) and (40) clearly indicate exponential convergence of the method with respect to the number of singular functions N a and an algebraic convergence with the number of Lagrange multipliers N λ . This fast convergence is visible in the results presented in Tables 1 and 2 and in the plot of Figure 4. In this study it is k=2 for the Sobolev space used. Also, quadratic basis functions were chosen (p=2). Now, if the two errors in (39) are forced to be equal to each other then we obtain 2 1/ 2 2 1/ 2 2 3 2 1 Hence using (41) and the "optimal" pair N a =N λ =13 we find that γ≈0.77≈0.80 ∈ (0,1). Therefore, with this value of γ and by prescribing parameter N a , one can find N λ from (41) and thus can have another pair of values. Using these values convergence of the method can be achieved.

Numerical Results for the 3-D Problem
In order to implement the SFBIM for the 3-D model problem, the boundary part S C is subdivided into standard 3-D boundary elements (having dimensions along z and θ directions only) as shown in Figure 5. Specifically, N E =N z ×N θ elements are employed. Figure 5 shows the division of cylindrical boundary S C into elements. As it is indicated by equations (35) and (36), calculations of the integrals, finally take place in a 2-D space ( Figure 5). On the same figure the bilinear basis functions are also presented. The total number of Lagrange multipliers is N λ =(N z +1)×(N θ +1). The surface integrals are estimated using a 9-point Gaussian quadrature rule over each element.
For all computations presented we take N p =1 in (26), i.e. the EFIFs are approximated by polynomials of linear form. As explained earlier, this choice is made in order to force boundary integrals on S D and S E to vanish. In this way the integration in the vicinity of point singularities O and D is avoided. Now, with this choice of N p , the parameter N in (26) is determined to be equal to 1; choosing a value N > 1 will not yield any additional terms in (26). Systematic runs have Equation Problems in Mechanics, with a Boundary Singularity been carried out for the model problem in order to study the effects of N a and N λ and of the other parameters, on the numerical results.
The runs for different values of R included various combinations between N a and N λ in an attempt to find the "ideal'' combination of these parameters, for which the method converges. Calculations were made with different combinations of N z and N θ but the pair of N z =2 and N θ =2 indicated (in combination with a suitable value of N a ) convergence of the method. Table 4 contains the values of the first two singular coefficients of A (a1) (z) (i.e. of a 1,1 and a 1,2 ) obtained with N λ =9, R =0.01 and for different values of N a . Note that since N p +1=2, the step increment for N a p =(N p +1)N a in all the trials, is equal to 2. One may immediately observe very high rate of convergence with respect to N a and great accuracy in the values obtained, even up to the 14 th decimal digit. Also, Table 5 shows the values of the first two EFIFs at the points z = 0.5 and 1, respectively, (i.e. of A (a1) (0.5) and A (a2) (1.0)) calculated for the same values of N λ and R and for various values of N a .
The values of z were selected so that a comparison can be made between the SFBIM and the energy projection method presented in [30]. We may observe, again, very fast convergence with respect to N a and high accuracy in the values obtained. In both Tables 4 and 5 convergence corresponds to the "optimal" combination of values N a =7, N a p =(N p +1)N a =14 and N λ = 9. The value of the radius of the domain, for our 3-D model problem, is taken as R = 0.01. A model problem of the same geometry and the same dimensions was also adopted in [30], where a special finite element technique was implemented in order to solve it.  Table 6 contains the converged approximate values of the singular coefficients obtained for R = 0.01 and for the "optimal" combination of N a and N λ . One may observe that the converged values of the coefficients, obtained with the SFBIM, are practically the same with those contained in the exact form of the EFIFs which is known for this model problem as already explained above. Finally, Table 7 compares the values obtained by the SFBIM for R=0.01 and z=1.0 with the results given by the energy projection method [30], for the same values of R and z. The exact solution, for our 3-D model problem, is also presented on the same table. Clearly, the SFBIM gives results which are comparable with the exact values.
In fact, for the accuracy considered in this work, the approximate values obtained with the SFBIM, coincide with the exact values. Also, it is obvious that the SFBIM, is significantly more accurate than the finite element technique which is implemented in [30]. This is because the SFBIM is a direct method and no post-processing is required. Also, the numerical error and the CPU time in the SFBIM, are much smaller than those observed in other numerical techniques, for this type of boundary value problems in Applied Mathematics. In Figure 6 there is a graph presenting with "squares" the approximation of the first EFIF A (a1) (z)=A 1 (z) obtained at convergence, for which it is N a p =N a (N p +1)=14 and N λ =9. On the same graph and with a solid line, the analytic form of the first EFIF A (a1) (z) is also presented. One may observe that, practically, the approximation of the polynomial A (a1) (z) coincides with the exact solution. Finally, the graph of Figure  7 shows how the error changes its value as N a p varies. In this case, the error is defined by using the infinity norm as follows:  At N a p =N a (N p +1)=14, the error is minimized and it is of the order of 10 -16 . This value of error corresponds to the optimum combination of the parameters (i.e. for N a =7, N p =1, N=1 and N λ =9). After a significant number of runs it was found that this is the minimum of the Ε vs. N a p graph for all values of N λ . A general conclusion derived from these graphs is that the method requires a small number of boundary elements and that it exhibits very fast convergence.

Conclusions
The Singular Function Boundary Integral Method (SFBIM) has been formulated for both a 2-D elliptic boundary value problem with a point boundary singularity and a 3-D Laplace equation problem with a straight-edge boundary singularity. With this method the singular coefficients and the Edge Flux Inensity Functions (EFIFs), are the primary unknowns in the 2-D and the 3-D cases, respectively. The latter are approximated by polynomials, the coefficients of which are unknowns in the formulation.
Both the singular coefficients, of the two-dimensional local solution and the coefficients of the EFIFs, in the 3-D case, are calculated directly and not by post-processing of the solution. The implementation of the SFBIM to both a 2-D and a 3-D Laplacian model problems, yielded highly accurate results for the singular coefficients and the EFIFs and exhibited fast convergence, as in other two-dimensional applications [9][10][11][12] of the method. For the planar problem, the numerical results are favorably compared with the analytic solution. The numerical results for the 3-D case compare favorably with both the exact solution and those of the energy projection method [30].
Currently a particular version of the method is being investigated for other three-dimensional Laplace equation problems. This approach will aim at the better treatment of the inner sum in Eq. (26). Also, a hybrid method is being developed in which the advantages of both the SFBIM and the finite element techniques are being exploited for problems in Mechanics with boundary singularities.