An accurate treatment of non-homogeneous boundary conditions for development of the BEM

An accurate treatment of non-homogeneous boundary conditions for development of the BEM A. Romeroa,∗, P. Galv´ına, A. Tadeub,c aEscuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla, Camino de los Descubrimientos, 41092 Sevilla, Spain bITeCons - Institute for Research and Technological Development in Construction, Energy, Environment and Sustainability, Rua Pedro Hispano s/ n., 3030-289 Coimbra, Portugal cADAI - LAETA, Department of Civil Engineering, University of Coimbra, P´olo II, Rua Lu´ıs Reis Santos, 3030-788 Coimbra, Portugal Abstract This paper proposes an enhancement of the treatment of non-homogeneous boundary conditions to improve the boundary element method (BEM) formulation. The standard formulation is modified by introducing the boundary conditions in the integral kernels. The boundary conditions are implicitly defined through known parameters depending on the geometry, rather than by prescribing nodal values as is done in the standard formulation. The main advantage of this procedure is that the right-hand side of the system of equations is integrated taking the exact distribution of loads into account. This approach is implemented in the B´ezierBernstein space to yield a geometry-independent field approximation. We use the B´ezier-Bernstein form of a polynomial as an approximation basis to represent both geometry and field variables. The application of the proposed method covers the resolution of complex boundary value problems as optimization with uncertain data, material modelling with graded impedance, and the definition of general boundary constraints. The performance of the proposed method is shown by solving the Helmholtz equation in two dimensions. The proposed method is numerically compared to the standard BEM formulation in two benchmark problems. Finally, the application of complex impedance boundary conditions is analysed in a numerical example. Keywords: Boundary conditions; Element approximation; Computer-aided design; B´ezier curve; Bernstein polynomial; Acoustic impedance. 1. Introduction Of the wide range of numerical methods, the most popular are the finite element method (FEM), the finite difference method (FDM) and the boundary element method (BEM). However, the BEM has a distinctive role as a boundary method, as opposed to FEM and FDM, as noted by Cheng and Cheng [1]. Owing to its ∗Corresponding author Email address: [email protected] (A. Romero) Preprint submitted to Engineering Analysis with Boundary Elements May 5, 2020 higher accuracy and efficiency, the BEM has been widely used in many engineering fields such as acoustics scattering, fracture mechanics, and soil wave propagation. The BEM is a collocation method that uses the fundamental solution to a problem as a weighted function [2]. The starting point for the formulation of the BEM is the boundary integral equation (BIE), which can be deduced from the weighted residual method, Green’s third identity or Betti’s reciprocal theorem. The boundary is then discretised into elements where the unknown values of the field variables coincide with the nodal positions. The field variables are usually approximated by nodal values that have a clear physical meaning. Typically, the element approximation involves a linear combination of nodal values and interpolation functions. Once the element approximation is introduced, we find that the BIE cannot be satisfied exactly unless the approximation of the field variables is exact [3]. This induces an error that determines the accuracy of the method, as will be discussed later. The concept of residual is introduced to represent the error by non-achievement of the integral representation. The collocation method is then used to require that the residual vanishes at a set of points called collocation points [4]. These points are usually chosen to be coincident with nodes and this procedure then yields to a nonsymmetric linear system of equations that relates the nodal variables. It proceeds by applying the necessary constraints for the solution of the boundary value problem. Typically, the system of equations is rearranged by moving the columns of the influence matrix according to the boundary conditions. The right-hand side of the system of equations is computed by multiplying those columns by the nodal values of the field variables. Some recent applications of the BEM regarding the reliability analysis coincided with the importance of achieving accurate and efficient numerical solutions [5–10]. Accordingly, Vable [11] discussed the various sources of errors in BEM that were classified as: i) formulation error, ii) interpolation error, iii) integration error, iv) continuity error, v) collocation error, vi) matrix conditioning error, and vii) mesh error. The interpolation error is one of the most important sources of error in BEM analysis [12] because of a lack of precision in the geometry discretisation and the field variables approximation. The loss of accuracy in the model discretisation is mainly given by the difficulty to represent the field variables through a polynomial approximation. Hence, many authors have conducted improvements related to geometry representation and field approximation. These studies covered the parametric representation of the geometry [13, 14], the isogeometric approach [15, 16], and the spectral formulations [17–20]. These publications have led to progress in the formulation of the BEM. This article presents an improvement in the treatment of boundary conditions based on the idea of developing a reliable method for engineering analysis. The novelty of the proposed method lies in the modification of BEM integral kernels to include the actual distribution of the boundary conditions instead of the polynomial approximation used in the standard BEM formulation. Therefore, these terms of the integral representation related to known values of the field variables are satisfied exactly since the exact 2 distribution given by boundary conditions is considered. The proposed method has been implemented in the BEM formulation based on the B´ezier-Bernstein space [21] to consider the exact boundary geometry. Also, this formulation allows the use of arbitrary highorder elements. Thus, the proposed method presents some enhancements according to Reference [11]: i) low interpolation error due to the use of high-order elements, ii) the collocation error is reduced because the boundary conditions are implicit at the integration kernels, and iii) no mesh error since exact CAD model geometry is considered. The outline of the paper is as follows. In Section 2, the definition of boundary conditions in the BEM is presented and the B´ezier-Bernstein space is briefly described. The proposed method is verified in Section 3 from two benchmark problems with a known analytical solution. A numerical example is presented and investigated in Section 4. Finally, the main results drawn from this research are summarised in the Conclusion section. 2. Numerical model In this work, we consider the Helmholtz equation, without loss of generality, to present the proposed formulation. The governing equation in the bounded domain Ω with the piecewise smooth boundary Γ defined by its normal outward n(x) is: ∇2u(x)+κ2u(x) = 0, x∈Ω, (1) where u(x) is the velocity potential and κ > 0 denotes the wavenumber. The problem definition is completed by setting the boundary conditions for on Γ = ∂Ω. Dirichlet and Neumann boundary conditions are defined as u(x) = u(x) and q(x) = ∂u(x)/∂n(x) = q(x), respectively. Also, it is possible to adopt the Robin boundary condition to gather the three types of boundary conditions in one expression: α(x)u(x) + β(x)q(x) = γ(x), x∈Γ, (2) where α(x), β(x) and γ(x) are known parameters along the boundary that allow the following definitions: β(x) = 0: u(x) = u(x) = γ(x)/α(x) Dirichlet boundary condition (3) α(x) = 0: q(x) = q(x) = γ(x)/β(x) Neumann boundary condition (4) α(x)6= 0, β(x)6= 0: q(x) = q(x)=(γ(x)−α(x)u(x))/β(x) Robin boundary condition (5) We can see from the above definitions that the Neumann condition is a particular case of the Robin type when α(x) = 0. Therefore, it is not necessary to distinguish between the two types of boundary conditions in the problem definition. Then, the boundary is defined as Γ = ΓuSΓqaccording to the following definitions: 3 u(x) = u(x) = γ(x)/α(x) , x∈Γu, (6) q(x) = q(x) = (γ(x)−α(x)u(x))/β(x) , x∈Γq. (7) It is assumed that κ2is not an eigenvalue of ∇2u(x) + κ2u(x) = 0 subject to the homogeneous form of the imposed boundary conditions. 2.1. Boundary element formulation The basic integral equation is the starting point as we introduced before, which can be written as follows [22]: c(ξ)u(ξ) = ZΓ∂u(x) ∂n(x)Ψ(x, ξ)−u(x)∂Ψ(x, ξ) ∂n(x)dΓ(x) , (8) where ξis the collocation point, q(x) = ∂u(x)/∂n(x) is the normal flux, Ψ(x, ξ) is the fundamental solution (the Hankel function H(1) 0(κ|x−ξ|)) at point xdue to a point source located at ξ, and the integral-free term c(ξ) depends only on the boundary geometry at the collocation point ξ. The dependence on the wavenumber κhas been omitted intentionally in the above equation for simplicity of notation. The basic integral representation implies that the potential u(ξ) is obtained by integrating the field variables u(x) and q(x) over the boundary Γ, weighted by the fundamental solution Ψ(x, ξ). Typically, the accuracy in the computation of u(ξ) depends on the approximation of the field variables at the boundary, and it is also conditioned by the precision of the integration scheme to solve Equation (8). The accuracy in the computation of u(ξ) will only rely on the precision of the quadrature rule when the field variables represent the actual solution at the boundary. This topic is discussed below. Following the scheme depicted in section 1, the boundary is discretised into Nelements with Γ = SN j=1 Γj. Thus, Equation (8) is rewritten as follows: c(ξ)u(ξ) = N X j=1 ZΓj∂u(x) ∂n(x)Ψ(x, ξ)−u(x)∂Ψ(x, ξ) ∂n(x)dΓ . (9) Equation (9) represents the discretised basic integral representation for u(ξ), which is computed as a piecewise integration of the field variables over the boundary. The field variables within an element Γjare interpolated from the nodal values ukusing element shape functions φk(x) of order p: u(x) = p X i=0 φi(x)ui=φ(x)ue. (10) Then, the element approximation is substituted into Equation (9) to yield the following expression: c(ξ)u(ξ) = N X j=1 ZΓj φ(x)Ψ(x, ξ)dΓ∂ue ∂n−ZΓj φ(x)∂Ψ(x, ξ) ∂n(x)dΓue. (11) 4 In this case, the basic integral Equation (11) differs from Equation (9) because the potential u(ξ) is given by integration of field nodal values interpolated by the element shape functions. The accuracy is constrained by the element approximation and should be lower than in Equation (9), even if uetakes the exact value of the field variable u(x) at nodal positions unless the element approximation of the field variables is exact. Therefore, the element approximation implies a loss of accuracy in the evaluation of the integral representation. Finally, the collocation method allows the definition of a system of equations that relates nodal values uand qover the boundary: Hu =Gq , (12) where Hand Gare the fully non-symmetrical boundary element system matrices. After boundary condition have been prescribed, Equation (12) is rewritten as AX =b[23], where Ais the matrix of coefficients, X collects the unknown nodal values, and bis the right-hand side. Our concern begins with the approximation of the right-hand side of the system of equations. This term is computed from the rearrangement of Equation (11) according to the boundary conditions, as: b(ξ) = X Γj∈ΓqZΓj φ(x)Ψ(x, ξ)dΓ∂ue ∂n−X Γj∈ΓuZΓj φ(x)∂Ψ(x, ξ) ∂n(x)dΓue, (13) where Γuand Γqdenote those parts of the boundary where Dirichlet and Robin (or Neumann) conditions are respectively defined. In this equation, the element approximation is used to interpolate the field variable u(x), and its derivative, through the nodal value ueand the element shape functions φ(x). Following this procedure, the boundary conditions are prescribed as known nodal values, but the approximation of the field variables through element shape function implies a loss of accuracy in the collocation method, as we mentioned above. Therefore, this approach does not seem to be sufficiently justified for the application of boundary conditions, given that the distribution of field variables is fully known. The novelty of the proposed method derives from this lack of accuracy. Thus, the actual distribution of u(x) is included in the integration kernels of Equation (13). Moreover, the definition of boundary conditions given by Equation (2) is used to develop a more general method. The basic integral representation of u given by Equation (9) is modified as follows, according to the boundary conditions given by Equation (2): c(ξ)u(ξ) = X Γj∈ΓuZΓj ∂u(x) ∂n(x)Ψ(x, ξ)dΓ−X Γj∈ΓuZΓj γ(x) α(x) ∂Ψ(x, ξ) ∂n(x)dΓ −X Γj∈ΓqZΓjα(x) β(x)Ψ(x, ξ) + ∂Ψ(x, ξ) ∂n(x)u(x)dΓ + X Γj∈ΓqZΓj γ(x) β(x)Ψ(x, ξ)dΓ. (14) Once the element approximation given by Equation (10) is introduced in the above equation, it becomes: 5 c(ξ)u(ξ)−X Γj∈ΓuZΓj φ(x)Ψ(x, ξ)dΓ∂ue ∂n +X Γj∈ΓqZΓjα(x) β(x)Ψ(x, ξ) + ∂Ψ(x, ξ) ∂n(x)φ(x)dΓue =X Γj∈ΓqZΓj γ(x) β(x)Ψ(x, ξ)dΓ−X Γj∈ΓuZΓj γ(x) α(x) ∂Ψ(x, ξ) ∂n(x)dΓ. (15) Equation (15) has been rearranged according to the boundary conditions. It should be noticed that the righthand side of this equation does not depend on the element approximation and, therefore, higher accuracy than a standard formulation is expected. After all known u(ξ) for ξbelonging to Γuare passed to the right-hand side, the following system of equations is obtained: AX =b. (16) The solution for the system of equations gives unknown nodal values of u(x) and q(x) at the boundary. Afterwards, the potential u(ξ) at internal point ξin the domain Ω can be computed from Equation (15), making the free term c(ξ) = 1. The solution at internal points also benefits from the proposed treatment of the boundary conditions. This approach has the following advantages over the standard BEM formulation: i) the computation of the right-hand side is independent of the element approximation; and ii) the Robin boundary condition is implicitly considered in the element matrix through parameters α(x), β(x) and γ(x) keeping the spatial information, instead of a relationship between nodal variables as is done in the standard BEM formulation. Therefore, the proposed method has high accuracy as will be shown in the next section. Although the proposed method is valid for any BEM formulation, it is implemented in the BEM formulation based on the B´ezier-Bernstein space [21] to show its capabilities. The next section summarises the main ideas for the method. 2.2. The BEM formulation in the B´ezier-Bernstein space The B´ezier-Bernstein space is used to describe the exact boundary geometry and to approximate the field variables. It is based on the application of polynomials in Bernstein form, that grows with the development of B´ezier curves rn(t) in computer-aided design. The B´ezier curve is defined over the interval t∈[0,1] as: rn(t) = n X k=0 bkBn k(t), (17) where bkare the control points used to approximate the geometry and nthe curve degree. The de Casteljau algorithm is often used for evaluating and splitting a B´ezier curve rn(t) at a given point t[24]. Although 6 the de Casteljau algorithm allows an easy evaluation of a B´ezier curve, it is computationally expensive. An efficient curve computation is achieved using the polar form (or blossom) of a B´ezier curve rn(t) [25], which defines a multi-affine transformation satisfying: bk=R(0,...,0 | {z } n−k ,1,...,1 | {z } k ), (18) where R(t1, . . . , tn) is computed as: R(t1, . . . , tn) = X I∩J=∅ I∪J={1,2,...,n}Y i∈I (1 −ti)Y j∈J tjb|J|. (19) Thus, a polynomial in Bernstein form can be formulated in the polar form, substituting Equation (18) into Equation (17) as follows: rn(t) = n X k=0 R(0,...,0 | {z } n−k ,1,...,1 | {z } k )Bn k(t) = R(t, . . . , t). (20) The B´ezier-Bernstein space is used to describe the exact element geometry as Γj(x) = rj n(t). Hence, the element integrals in Equation (15) are rewritten in the univariate basis t∈[0,1] as [21, 26]: ZΓj f(x, ξ)dΓ = Z1 0 f(x(t), ξ) drj n(t) dt dt, (21) where f(x, ξ) represents the integration kernels. Moreover, Equation (21) can be transformed into the integration interval [−1,1] to employ a Gauss-Legendre quadrature. Moreover, the proposed method employs the Lagrange interpolant relative to the Bernstein basis for the field variable approximation to an element [27]. The field approximation given by Equation (10) interpolates n+ 1 nodal values through the element shape functions Pi n=φiof order n, for i= 0, . . . , n. The proposed element is defined by the nodal positions tjin the univariate basis. The Lagrange interpolant Pi nderived from the Bernstein basis must fulfil the following condition at element nodes tj: Pi n(tj) = φi(tj) = n X k=0 ckBn k(tj) = δij, j = 0, . . . , n, (22) where, ckare the control points used to define the Lagrange polynomial Pi n. This condition is commonly expressed as a linear system of equations through the Bernstein-Vandermonde matrix Aij =Bn i(tj) as: 7               Bn 0(t0)Bn 1(t0). . . Bn k(t0). . . Bn n(t0) Bn 0(t1)Bn 1(t1). . . Bn k(t1). . . Bn n(t1) . . . Bn 0(ti)Bn 1(ti). . . Bn k(ti). . . Bn n(ti) . . . Bn 0(tn)Bn 1(tn). . . Bn k(tn). . . Bn n(tn)                                          ci 0 ci 1 . . . ci k . . . ci n                            =                            0 0 . . . 1 . . . 0                            .(23) Thus, the element shape function φiis defined by control points obtained from the solution of (23). Then, the field approximation given by Equation (10) in the univariate basis tbecomes: u(t) = p X i=0 φi(t)ui= p X i=0 (n X k=0 ck kBn k(t))ui= p X i=0 Ri(t, . . . , t)ui, (24) where the evaluation of the element shape function φi(t) also benefits from the computational advantages of using the polar form Ri(t1, . . . , tn) according to Equation (19). Once the geometry and the field approximation given by Equations (20) and (24) are introduced in Equation (15), the boundary integrals are computed using a standard Gauss-Legendre quadrature with p+ 1 integration points whenever the collocation point is sufficiently distant from the integration element. Otherwise, the solution of singular or weakly singular integrals is numerically computed using a smoothing transformation using a Gauss-Legendre quadrature [28]. Figure 1 shows a scheme for the treatment of singular and weakly singular integrals. This figure represents a collocation point ξand an integration element Γj. The parametric coordinate t(ξ) is found from ξas the point that minimizes the distance r(ξ) to the integration element. t(ξ) coincides with the coordinate of this node if the collocation point belongs to the integration element. After point t(ξ) has been identified, the element integration is subdivided into two intervals. The element integrals are numerically solved by a smoothing transformation ϕl,r(s) [28]. We chose a critical radius r∗between the collocation point and the integration element for the identification of singular or nearly singular integrals. Otherwise, the collocation point is sufficiently distant from the integration element and the resulting integrals are regular. The asymptotic behaviour of the fundamental solution is accounted for by selecting the critical radius r∗when the integral becomes singular [21]. 3. Numerical verification In this section, we analyse the performance of the proposed method for solving the Helmholtz equation in a square boundary Ω := [−1,1]×[−1,1] at a high wavenumber κ= 100 rad/m. Four linear B´ezier patches were used to define the boundary geometry: 8 ξ t(ξ) Γj t r(ξ) ϕr(s) ϕl(s) Figure 1: Treatment of singular and weakly singular integral. Γ1:= [−1,1] ×[−1,−1] , Γ2:= [1,1] ×[−1,1] , Γ3:= [−1,1] ×[1,1] , Γ4:= [−1,−1] ×[−1,1] . (25) The proposed method was then tested in two benchmark problems. Different boundary conditions were chosen such that the exact solution satisfies uI(x) = exp(ικd·x), where the polarised direction was set to d= [1,1], and the unit imaginary number was denoted by the Greek letter ιto prevent confusion with some subscripts used in the paper. Numerical results were compared with a reference solution using the l2scaled error 2to assess the accuracy. A convergence investigation was carried out for several element lengths hwith successive p−enrichment. Three different discretisation schemes were tested with element lengths given by κh = 7.5, κh = 3 and κh = 1. The element shape functions were obtained from the interpolation functions defined at Chebyshev points of the first kind [21]. The element order was increased until convergence was reached. For this purpose, we considered that the problem solution was properly approximated if the error satisfied log(2(h, p − 1)/2(h, p + 1)) ≤1. The accuracy of the proposed method was compared with a standard BEM formulation (11). The type of element and the integration scheme was the same in both methods. Only the treatment of boundary conditions was modified. 3.1. Example 1 In this example, Dirichlet and Neumann conditions were prescribed on the boundaries Γ1,3and Γ2,4, respectively. The boundary conditions were prescribed as follows, according to Equation (2): 9 (a) A0= 0.2λ(b) A0= 0.4λ (c) A0= 0.6λ(d) A0= 0.8λ Figure 8: Imaginary part of the pressure field at the boundary computed from variable impedance (blue line) and fixed impedance (red line) models and over the domain (colour map) for different roughness amplitude A0. The characteristic roughness period is T0= 0.4 m. 16 (a) A0= 0.2λ(b) A0= 0.4λ (c) A0= 0.6λ(d) A0= 0.8λ Figure 9: Imaginary part of the pressure field at the boundary computed from variable impedance (blue line) and fixed impedance (red line) models and over the domain (colour map) for different roughness amplitude A0. The characteristic roughness period is T0= 0.05 m. 17 5. Conclusions This work has proposed a new treatment of boundary conditions to improve the accuracy of the BEM. The boundary conditions were defined as the type of Dirichlet, Neumann and Robin conditions, using known parameters along the boundary that were included in the integral kernels. This procedure avoided the element approximation and allowed the domain geometry and boundary conditions to be considered exactly in the computation of the right-hand side of the BEM system of equations, as opposed to the standard formulation, which interpolates the boundary conditions from nodal values using the element shape functions. The proposed method was implemented in the BEM formulation based on the B´ezier-Bernstein space, which allowed an independent approximation of the geometry and the field variables. However, the method can be extended to other BEM formulations. The accuracy of the proposed method was demonstrated in two benchmark problems. The results have shown that: •The discretisation of the basic integral equation defines a piecewise integration when field variables are known. The best solution for an internal domain point is obtained if the field variables achieve the exact solution at the boundary. This means the accuracy depends only on the boundary discretisation and the quadrature rule to compute element integrals. The element approximation of the field variables induces a loss of accuracy related to the element order. •We have found that the way the BEM prescribes the boundary conditions could be inaccurate. Typically, boundary conditions are directly incorporated into the system of equations, which is obtained from element integration using shape functions. Therefore, the right-hand side has a base error related to the approximation of distributed loads by the element shape functions. The element order improves the precision in the numerical results. •Although the loss of accuracy can be reduced with a hp-refinement, the proposed method yielded better results in the domain than the standard BEM formulation. A numerical example has been proposed to show the robustness of the proposed method to model complex boundary conditions that are dependent on the geometry. The BEM formulation based on the B´ezierBernstein space allowed an exact definition of the geometry regardless of the field variable approximation which, in addition to the proposed treatment of the boundary condition, allowed the accurate definition of variable impedance in a duct with a rough surface. Acknowledgements The research work presented was supported by the Spanish Ministry of Economy and Competitiveness (Ministerio de Econom´ıa y Competitividad) through research project BIA2016-75042-C2-1-R and by the 18 Spanish Ministry of Science, Innovation and Universities (Ministerio de Ciencia, Innovaci´on y Universidades) through the scholarship Salvador de Madariaga scholarship, Reference PRX19/00298. We are also grateful for support from project POCI-01-0247-FEDER-017759, funded by Portugal 2020 through the Operational Programme for Competitiveness Factors (COMPETE 2020). Financial support is gratefully acknowledged. The authors also wish to acknowledge the support provided by the Andalusian Scientific Computing Centre (CICA). References [1] A.H.-D. Cheng, D.T. Cheng, Heritage and early history of the boundary element method, Engineering Analysis with Boundary Elements 29 (3) (2005) 268–302. [2] C. Brebbia, J. Dom´ınguez, Boundary Elements An Introductory Course, 1st Edition, WIT Press, 1992. [3] M. Dehghan, H. Hosseinzadeh, Improvement of the accuracy in boundary element method based on high-order discretization, Computers and Mathematics with Applications 62 (12) (2011) 4461–4471. [4] M. Schanz, T. R¨uberg, L. Kielhorn, Time Domain BEM: Numerical Aspects of Collocation and Galerkin Formulations, Springer Netherlands, Dordrecht, 2009, pp. 415–432. [5] E.D. Leonel, A. Chateauneuf, W.S. Venturini, P. Bressolette, Coupled reliability and boundary element model for probabilistic fatigue life assessment in mixed mode crack propagation, International Journal of Fatigue 32 (11) (2010) 1823–1834. [6] C. Su, S. Zhao, H. Ma, Reliability analysis of plane elasticity problems by stochastic spline fictitious boundary element method, Engineering Analysis with Boundary Elements 36 (2) (2012) 118–124. [7] C. Su, J. Xu, Reliability analysis of Reissner plate bending problems by stochastic spline fictitious boundary element method, Engineering Analysis with Boundary Elements 51 (2015) 37–43. [8] M.S. Chowdhury, C. Song, W. Gao, C. Wang, Reliability analysis of homogeneous and bimaterial cracked structures by the scaled boundary finite element method and a hybrid random-interval model, Structural Safety 59 (2016) 53–66. [9] A. Johari, A. Heydari, Reliability analysis of seepage using an applicable procedure based on stochastic scaled boundary finite element method, Engineering Analysis with Boundary Elements 94 (2018) 44–59. [10] L. Morse, Z. Sharif Khodaei, M.H. Aliabadi, A multi-fidelity boundary element method for structural reliability analysis with higher-order sensitivities, Engineering Analysis with Boundary Elements 104 (2019) 183–196. [11] M. Vable, Controlling errors in the process of automating boundary element method analysis, Engineering Analysis with Boundary Elements 26 (5) (2002) 405–415. [12] Z. Zhao, X. Wang, Error estimation and h adaptive boundary elements, Engineering Analysis with Boundary Elements 23 (10) (1999) 793–803. [13] E. Zieniuk, B´ezier curves in the modification of boundary integral equations (BIE) for potential boundary-values problems, International Journal of Solids and Structures 40 (9) (2003) 2301–2320. [14] E. Zieniuk, K. Szersze´n, Nonelement boundary representation with B´ezier surface patches for 3D linear elasticity problems in parametric integral equation system (PIES) and its solving using Lagrange polynomials, Numerical Methods for Partial Differential Equations 34 (1) (2018) 51–79. [15] R.N. Simpson, S.P.A. Bordas, J. Trevelyan, T. Rabczuk, A two-dimensional Isogeometric Boundary Element Method for elastostatic analysis, Computer Methods in Applied Mechanics and Engineering 209-212 (2012) 87–100. [16] M.A. Scott, R.N. Simpson, J.A. Evans, S. Lipton, S.P.A. Bordas, T.J.R. Hughes, T.W. Sederberg, Isogeometric boundary element analysis using unstructured T-splines, Computer Methods in Applied Mechanics and Engineering 254 (2013) 197–221. 19 [17] A.T. Patera, A spectral element method for fluid dynamics: Laminar flow in a channel expansion, Journal of Computational Physics 54 (3) (1984) 468–488. [18] B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 1996. [19] U. Iemma, V. Marchese, L. Morino, High-order bem for potential transonic flows, Computational Mechanics 21 (3) (1998) 243–252. [20] M. L¨ohndorf, J.M. Melenk, Wavenumber-explicitHp-BEM for high frequency scattering, SIAM Journal on Numerical Analysis 49 (6) (2011) 2340–2363. [21] A. Romero, P. Galv´ın, J. C´amara-Molina, A. Tadeu, On the formulation of a BEM in the B´ezier–Bernstein space for the solution of Helmholtz equation, Applied Mathematical Modelling 74 (2019) 301 – 319. [22] T.W. Wu, Boundary Element Acoustics Fundamentals and Computer Codes, Advances in Boundary Elements, WIT Press, 2000. [23] J. Dom´ınguez, Boundary elements in dynamics, Computational Mechanics Publications and Elsevier Aplied Science, Southampton, 1993. [24] R.T. Farouki, The Bernstein polynomial basis: A centennial retrospective, Computer Aided Geometric Design 29 (6) (2012) 379–419. [25] L. Ramshaw, Blossoming: A connect-the-dots approach to splines, Digital Equipment Corporation SRC Report No. 19. [26] H.-P. Seidel, An introduction to Polar Forms, IEEE Computer Graphics and Applications 13 (1) (1993) 38–46. [27] R.T. Farouki, T.N.T. Goodman, T. Sauer, Construction of orthogonal bases for polynomials in Bernstein form on triangular and simplex domains, Computer Aided Geometric Design 20 (4) (2003) 209–230. [28] G. Monegato, L. Scuderi, Numerical integration of functions with boundary singularities, Journal of Computational and Applied Mathematics 112 (1-2) (1999) 201–214. [29] I. Babuska, F. Ihlenburg, E.T. Paik, S.A. Sauter, A Generalized Finite Element Method for solving the Helmholtz equation in two dimensions with minimal pollution, Computer Methods in Applied Mechanics and Engineering 128 (3-4) (1995) 325–359. [30] J. Biermann, O. von Estorff, S. Petersen, C. Wenterodt, Higher order finite and infinite elements for the solution of Helmholtz problems, Computer Methods in Applied Mechanics and Engineering 198 (13–14) (2009) 1171 – 1188. 20