{'fraction_non_alphanumeric': 0.03590780333572327, 'fraction_numerical': 0.018429929829922777, 'mean_word_length': 4.5746297995020315, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New hybrid quadrature schemes for weakly singular kernels applied to isogeometric boundary elements for 3D Stokes flow

This work proposes four novel hybrid quadrature schemes for the efficient and accurate evaluation of weakly singular boundary integrals (1/r kernel) on arbitrary smooth surfaces. Such integrals appear in boundary element analysis for several partial differential equations including the Stokes equation for viscous flow and the Helmholtz equation for acoustics. The proposed quadrature schemes apply a Duffy transform-based quadrature rule to surface elements containing the singularity and classical Gaussian quadrature to the remaining elements. Two of the four schemes additionally consider a special treatment for elements near to the singularity, where refined Gaussian quadrature and a new moment-fitting quadrature rule are used. The hybrid quadrature schemes are systematically studied on flat B-spline patches and on NURBS spheres considering two different sphere discretizations: An exact single-patch sphere with degenerate control points at the poles and an approximate discretization that consist of six patches with regular elements. The efficiency of the quadrature schemes is further demonstrated in boundary element analysis for Stokes flow, where steady problems with rotating and translating curved objects are investigated in convergence studies for both, mesh and quadrature refinement. Much higher convergence rates are observed for the proposed new schemes in comparison to classical schemes.


Introduction
One of the main difficulties of the boundary element method (BEM) is the efficient approximation of singular integrals that appear in the boundary integral equation (BIE). One approach to overcome this challenge is to superimpose known solutions to the unknown fields such that singularities are removed (Cruse, 1974;Liu and Rudolphi, 1999;Liu, 2000;Klaseboer et al., 2009). This approach is referred to as nonsingular, regularized or desingularized BEM. It is applied to linear elasticity (Scott et al., 2013;Taus et al., 2019), to Stokes flow (Taus et al., 2016;Harmel et al., 2018) and to the Helmholtz equation (Simpson et al., 2014;Peake et al., 2015), among others. Klaseboer et al. (2012) further apply nonsingular BEM to fluid mechanics by considering Stokes equations for viscous flow, Laplace equation for potential flow and Helmholtz equation for free-streamline flow.
The nonsingular BEM avoids singular integrals and is thus commonly used in recent papers, but it requires additional integrals and knowledge about analytical solutions and is further disadvantageous in efficiency and implementation (Khayat and Wilton, 2005). An alternative approach considered is to use the classical BIE without regularization and to approximate the singular boundary integrals with special quadrature rules. Such an approach is considered in this paper for weakly singular boundary integrals (kernel proportional to 1/r) 1 that are required and sufficient for steady and three-dimensional BIEs of several linear problems such as Stokes flow, potential flow, elastostatics, heat conduction and acoustics. Existing quadrature rules for surface elements containing the singularity (denoted singular elements) are discussed in Sec. 1.1 and for the adjacent elements near to the singularity (denoted near singular elements) in Sec. 1.2. The application of these quadrature rules to various boundary element (BE) problems is then discussed in Sec. 1.3.

Approximation of weakly singular integrals
Considerable progress has been made with quadrature rules based on variable transformations that map the physical domain to a parent domain so that the singularity is removed through the introduction of the Jacobian of the transformation, e.g. see Schwartz (1969); Takahasi and Mori (1973); Lean and Wexler (1985); Khayat and Wilton (2005) and Cano and Moreno (2015). The Duffy transformation from a square to a triangle (Fairweather et al., 1979;Duffy, 1982) is of particular interest here. In finite element (FE) analysis for fracture mechanic, numerical quadrature based on the Duffy transformation (Duffy quadrature in short) is used to integrate singular shape function derivatives (Tracey, 1971;Stern and Becker, 1978) and in the context of the extended FE method (Laborde et al., 2005;Béchet et al., 2005;Lv et al., 2018). Duffy quadrature is further used to integrate rational bubble functions with multiple singularities in FE analysis for Stokes flow (Schneier, 2015). More recently, Tan et al. (2019) propose the Duffy-Distance transformation that includes a further mapping that has to be adjusted depending on the shape of the element. Numerical investigations show that weakly singular integrals on spheres and cylinders can be approximated mostly with high accuracy, whereas the location of the singularity and the aspect ratio of the elements have a strong negative impact. Other important quadrature rules based on variable transformations are the tanh rule (Haber, 1977) and the sinh-tanh rule (Borwein and Ye, 2006), the quadrature rule from Telles (1987) that is based on a non-linear coordinate transformation and a polar coordinate transformation that is also applicable to elements with collapsed edges (Taus et al., 2016). The application of the quadrature rules from Sec. 1.1, including Duffy quadrature, to various BE problems is discussed in more detail in Sec. 1.3.
In addition to variable transformation methods, there are many other approaches for the approximation of weakly singular boundary integrals: Guiggiani et al. (1992) propose the approximation of singular integrals by singularity subtraction, Dautray and Lions (1985) and Atkinson and Chien (1995) use piecewise polynomial functions to approximate surface and integrals, while Niu et al. (2005) propose semi-analytical integration of the singular kernel. Recently, Ochiai (2022) propose the triple-reciprocity BEM that allows a direct integration on Lagrange elements and Velázquez-Mata et al. (2022) present a method that constructs efficient quadrature rules for the approximation of singular integrals from the finite part of known integrals including the shape functions.

Approximation of weakly near singular integrals
Weakly singular integrals are approximated efficiently with adaptive quadrature rules (Lachat and Watson, 1976;Gao and Davies, 2000) by dividing near singular elements into subelements 1 r denotes the distance between field point and the singularity and applying Gauss-Legendre quadrature to the subelements. The number of subelements and the number of quadrature points for each subelement are determined individually based upon analytical (Stroud and Secrest, 1966) or numerical (Bu and Davies, 1995) criteria. Similarly as on singular elements (see Sec. 1.1), quadrature rules based on variable transformation can also be used on nearly singular elements to approximate weakly singular integrals. The quadrature rule from Telles (1987) is also applicable to near singular elements, while Duffy quadrature (Duffy, 1982) is only accurate on the singular element. Other variable transformation considered for the approximation of weakly near singular integrals include the radial variable transformation (Hayami and Matsumoto, 1994), the distance transformation (Ma and Kamiya, 2002;Xie et al., 2013) and the sinh transformation (Johnston et al., 2007;Xie et al., 2021). The quadrature scheme presented in Gong et al. (2020) combines the benefits of the sinh transformation method and adaptive methods such it is capable of integral kernels of the type 1/r, 1/r 2 and 1/r 3 on near singular elements. Many other quadrature rules for weakly singular integrals on nearly singular elements can be found in the literature as for example the parameterized Gaussian quadrature (Lutz, 1992), line integral approaches (Krishnasamy et al., 1994;Liu, 1998) and analytical (Padhi et al., 2001;Zhou et al., 2008) and semi-analytical methods (Sladek et al., 2001;Niu et al., 2005;Han et al., 2022).
An alternative approach is to construct quadrature rules by the moment fitting method, i.e. solve a variant formulation of the moment fitting equations for the quadrature weights and abscissas. First, moment fitting quadrature rules for the exact integration of polynomial functions on triangles (Lyness and Jespersen, 1975) and on quadrilaterals (Dunavant, 1985) were presented. Moment fitting quadrature rules for polynomials are formalized and generalized by Wandzurat and Xiao (2003), while a numerical algorithm for polynomials on triangles, squares and cubes is given by Xiao and Gimbutas (2010). In all these approaches the integral on the right hand side of the moment fitting equations is computed analytically. Moment fitting quadrature for discontinuous functions is proposed for crack propagation analysis with the extended FE method (Mousavi and Sukumar, 2010b;Zhang et al., 2018), where the moment fitting integrals are evaluated partially analytically and partially numerically applying the method of Lasserre (1998Lasserre ( , 1999. Joulaian et al. (2016) and Hubrich et al. (2017) manipulate the moment fitting integrals so they can be computed by Gaussian quadrature numerically. This quadrature rule approximates discontinuous functions on domains of arbitrary geometry and topology highly efficiently. Moment fitting quadrature rules are further used for the construction of highly efficient quadrature rules in FE analysis for shells (Zou et al., 2021). The present paper proposes a new Gauss-Legendre quadrature rule with adjusted weights that are determined by the moment fitting method. The proposed quadrature rule determines weakly near singular integrals exactly on plane surfaces with regular quadrilateral elements and accurately on curved surfaces. It requires only a moderate number of quadrature points and uses the same abscissas as standard Gauss-Legendre quadrature. The Gauss-Legendre quadrature rule with adjusted weights is thus simple to implement and computationally efficient.

Application to BE analysis
Solving boundary value problems using BE analysis requires to approximate the BIE on the entire surface including singular, near singular and the remaining regular elements. Numerous approaches using the quadrature rules from Sec. 1.1 and Sec. 1.2 to approximate the boundary integrals efficiently can be found in the literature, so the following enumeration of application is not exhaustive.
Duffy quadrature (Duffy, 1982) is applied to the singular element for Stokes flow problems (Varnhorn, 1989;Johnson and Scott, 1989;Barakat and Shaqfeh, 2018) and for acoustic prob-lems (Amini and Harris, 1990). In these approaches, the integrals on the remaining elements are approximated by standard Gauss-Legendre quadrature or by simple trapezoidal rules. Venås and Kvamsdal (2020) additionally consider local refinement of Gauss-Legendre quadrature on near singular elements for acoustic problems without investigating its influence in detail. The quadrature rule from Telles (1987) is combined with Gauss-Legendre quadrature for linear elastic problems considering homogeneous material (Karam and Telles, 1988;Simpson et al., 2012) and liquid inclusions (Dai et al., 2021) and for FSI problems considering Stokes Flow (Patiño and Nieto-Londoño, 2021). The adaptive integral method for near singular integrals (Gao and Davies, 2000) is used to solve Laplace's equation by isogeometric BE analysis for potential problems (Gong and Dong, 2017), thermal problems (Gong et al., 2018) and thermoelasticity problems (Gong et al., 2020), whereas the singularity subtraction method (Guiggiani et al., 1992) is applied in Jarvenpaa et al. (2006) to electromagnetics and in Peng et al. (2017) to fracture mechanics. Ata and Sahin (2018) consider the tanh-sinh quadrature rule (Borwein and Ye, 2006) for BE analysis of Stokes flow problems and Keuchel et al. (2017) apply the sinhquadrature (Johnston et al., 2007) to nearly singular integrals in the Burton-Miller formulation of the Helmholtz equation for acoustic problems. The quadrature rule based on polar coordinate transformation (Taus et al., 2016) is extended in Taus et al. (2019) to handle elements with large curvatures and aspect ratios efficiently and to improve the accuracy on near singular elements. The extended quadrature rule is applied to linear elasticity considering Gaussian quadrature on the regular elements. Giuliani et al. (2018) provide a generalized parallel implementation for adaptive, geometry aware, and high order boundary element methods using quadrature rules from Lachat and Watson (1976), Duffy (1982) and Telles (1987) without investigating their influence on the accuracy of the BE results.
The above mentioned approaches are limited in application (geometry and discretization) or lack systematic investigation. This paper, therefore, proposes four hybrid quadrature schemes which • approximate weakly singular integrals on singular, near singular and regular elements accurately, • adaptively combine a new moment fitting quadrature rule with Duffy quadrature (Duffy, 1982) and rings of Gauss-Legendre quadrature with varying quadrature density, • are systematically studied for discretization and quadrature refinement, • show fast convergence for three-dimensional BE examples of steady Stokes flow with flat and curved isogeometric surfaces, • are applicable to various discretizations (including isogeometric, Lagrange and Hermite shape functions) using quadrilateral elements of arbitrary genus zero surfaces.
The remainder of this paper is organized as follows: Sec. 2 presents an overview of the underlying BE theory for incompressible Stokes flow. Boundary quadrature of weakly singular integrals on singular and near singular elements is discussed in Sec. 3. In addition to known rules, the new moment fitting quadrature rule is also presented in this chapter. Sec. 4 presents the four new hybrid quadrature scheme and investigates their accuracy and efficiency in detail. The hybrid quadrature schemes are applied to BE analysis for three-dimensional Stokes flow problems in Sec. 5. The paper then concludes with Sec. 6.

Fluid constitution
The three-dimensional Cauchy stress tensor for an incompressible Newtonian fluid with dynamic viscosity η is given by where D := grad v + (grad v) T /2 is the symmetric part of the velocity gradient and p = −tr σ/3 is the fluid pressure in domain F. The incompressibility of the fluid is enforced by the continuity equation that is given in absence of mass sources or sinks by where v denotes the velocity.

Fluid equilibrium
The steady-state motion of a viscous 2 fluid flow is governed by the steady Stokes equation together with the continuity equation (2) and the Dirichlet and Neumann boundary conditions with velocity v, surface traction t := σ n and outward unit normal vector n. Dirichlet and Neumann boundary regions are denoted by ∂ d F and ∂ n F, respectively and the entire boundary of the fluid domain is denoted by S := ∂ n F ∪ ∂ d F. Apart from (4), a condition on pressure p is needed on S to define the flow field. 3 This work focuses on the efficient quadrature of boundary integrals for boundary element (BE) analysis of steady Stokes flow. However, all following results apply to arbitrary weakly singular boundary integrals including BE for the Helmholtz equation and linear elasticity. The equilibrium equation for linear elasticity is even mathematically equivalent to the steady Stokes equation (3). 4 The quadrature schemes presented in this work are further applicable to unsteady BE problems since the singular integrals investigated here also appear there.

Boundary representation of the fluid flow
This section discusses a representation of steady Stokes equation (3) The distance between an arbitrary field point x ∈ F and source point y is defined as r := r with r := x − y. The point load can then be written as where δ(r) denotes the Dirac delta function and b ∞ denotes a constant force vector that describes magnitude and orientation of the point load. The solutions for velocity and stress of a singularly forced Stokes flow (5) are given by the components where Latin indices i, j and k (from one to three) indicate components in Cartesian coordinates following Einstein's summation convention. The Green's functions for velocity and stress are given by wherer i is the component ofr := r/r and δ ij denotes Kronecker's delta. From (8) follows that the indices of both Green's functions are arbitrarily interchangable. 5 Further, the Green's functions are symmetric with respect to r such that G ij (r) = G ij (−r) and T ijk (r) = −T ijk (−r).
Integrating the Green's functions over the fluid domain F and applying the divergence theorem yield the important integral identities for the velocity Green's function and the stress Green's function where pv d denotes the Cauchy principal value integral.

Boundary integral equation
The governing equations for incompressible Stokes flow, (1), (2) and (3), are transformed into a boundary representation that is referred to as the boundary integral equation (BIE). The BIE allows to determine the fluid velocity at any point y ∈ F by solving only surface integrals defined on S. To derive the BIE, a general incompressible Stokes flow (v, σ) is related to a singularly forced incompressible Stokes flow (v ∞ , σ ∞ ) by use of the introduced Green's functions (7) and Lorentz reciprocal theorem (Lorentz, 1896). Integrating the resulting identity over the fluid domain and applying the divergence theorem yields the BIE where y ∈ F denotes a source point within the fluid domain. Considering the limiting behavior of the singular integrals yields the BIE for a point on the surface (y ∈ S) where ϕ denotes the solid angle of the enclosed domain. 6 Note that (11) and (12) consider exterior flow problems where an object or particle enclosed by surface S is surrounded by the infinite fluid domain F. Similar BIEs for interior and two-sided problems and their derivations can be found, see for instance Pozrikidis (1992) and Harmel (2022).
BIEs for unsteady Stokes flow and corresponding Green's functions can be found for instance in Power and Partridge (1993) and Jiang et al. (2012). Those BIEs involve additional complexities like trivariate domain integrals and Laplace or Fourier transforms.

Boundary discretization and collocation
The surface geometry S and the continuous BIE (12) are discretized into n el elements and n no nodes (control points in isogeometric analysis). The discretized BIE is then collocated at n no source points y A ∈ S, with A = 1, . . . , n no , to obtain a square boundary element (BE) system. The source points y A are thus referred to as collocation points in the following. Note that integration over the whole surface has to be performed once for each collocation point. Explicit expressions for the discretized BIE and the collocated BE system can be found in Appendix A.
This paper considers isogeometric discretizations, where the control points are not necessarily located on the surface, which makes them unsuitable for collocation. Instead, the location of the collocation points are determined by the Greville abscissae (see e.g. Greville (1964), Johnson (2005) and Auricchio et al. (2010)).

Singular nature of the Green's functions
The Green's functions (8) become singular for r approaching to zero. Source point y is therefore a singularity of both Green's functions. According to their definitions, the velocity Green's function G ij (r) is proportional to 1/r, while the stress Green's function T ijk (r) is proportional to 1/r 2 . This kind of function behavior is referred to as weakly singular and strongly singular, respectively. Since these singularities are an essential component of boundary element (BE) analysis, their asymptotic behavior for r → 0 is numerically investigated on spherical and on flat surfaces in the following. The Green's functions behavior is evaluated by a scalar invariant, the Frobenius norm, that is defined for the second order tensor G(r) = G ij (r) e i ⊗ e j by and for the third order tensor T (r) = T ijk (r) e i ⊗ e j ⊗ e k by where (13) and (14) again follow Einstein's summation convention. The second-order tensor T n is part of an integral kernel of the BIE in Sec. 2.3.2 and therefore of particular interest for BE analysis. Its asymptotic behavior is additionally investigated here by determining T n F according to (13).
A spherical surface with radius R is considered first, see Fig. 1a. Given source point y, the Green's functions behavior is investigated as field variable x approaches y along the blue curve. Fig. 1c shows that G F is proportional to 1/r (weakly singular ), while T F is proportional to 1/r 2 (strongly singular ). The G F and T F show the asymptotic behavior that is expected from the Green's functions definition (8). Fig. 1b further shows that T n F is at first only weakly singular, although T F is strongly singular. Let us take a closer look on the expression to substantiate this numerical result with analytical findings: From (8) follows that T n F ∝r · n/r 2 . For a smooth surface with continuously varying normal vector, the unit vectorsr and n become more and more orthogonal as x approaches y, i.e. r → 0. The dot product r · n thus changes linearly with respect to r (see Fig. 1b) such that the behavior of T n F for r → 0 is mildened to be only weakly singular. However, for r < 10 −8 R the behavior of T n F vs. 1/r changes suddenly: From a perfectly linear relation ( T n F ∝ 1/r) to a distorted cubic relation ( T n F ∝ 1/r 3 ). This behavior results from the inaccurate numerical computation of the dot productr ·n for very small r that is shown in Fig. 1b. Within the numerical quadrature, T n is typically not required to be evaluated for such small values of r. The discussed change in the convergence behavior of T n has therefore no effect on the approximations of the singular integrals in Sec. 2.3.2. a.
y Figure 1: Singular nature of Green's functions: a. Spherical surface (radius R) with singularity at y; b. Limiting behavior of |r · n| for r := r approaching 0; c. Limiting behavior of the Green's functions.
Second, a source point on a tilted planar surface of side length L is considered (see Fig. 2a). It can be seen in Fig. 2c that the Green's functions show the same asymptotic behavior here as for spherical surfaces: G F is weakly singular and T F is strongly singular for r → 0. However, the numerical evaluation of T n F yields an overall cubic asymptotic behavior deviating from theoretical expectations: For a source point located on a plane surface, the vectors n andr are orthogonal. Analytically their dot product thus yields n ·r = 0 and the expression T n consequently vanishes. The deviating behavior of T n F is caused by the diverging numerical result for r → 0 on plane surfaces that is shown in Fig. 2b. On flat surfaces it is therefore advisable to use the known value T n = 0 in the BE computation instead of evaluating it numerically. Nevertheless, it should be noted that moderate values of r result in very small T n F (smaller than 10 −10 for r > 10 −2 L). a.

Boundary quadrature of singularities
This section investigates singular boundary quadrature on the plane sheet from Fig. 3a that consists of 3 × 3 plane biquadratic B-spline elements with side length L e . The depicted collocation point y 0 is located on the center element which is therefore referred to as singular element (blue, surface S sing ), whereas the adjacent elements are referred to as near singular elements (green, surface S near ). The weakly singular integral that is representative of the behavior of the two integrals in the BIE (12) is approximated on S sing and on S near using various quadrature rules: Sec. 3.1 studies Duffy quadrature and classical and modified Gauss-Legendre quadrature on the singular element and Sec. 3.2 investigates classical Gauss-Legendre quadrature and a new Gauss-Legendre quadrature rule with adjusted weights on the near singular elements. The quadrature accuracy is evaluated by means of the relative quadrature error where I h denotes the numerical approximation of integral (15)  S near S sing Figure 3: Boundary quadrature of singularities: a. Singular element (blue, S sing ) with collocation point y 0 and the surrounding near singular elements (green, S near ). b. Duffy quadrature points on the singular element forñ gp = 3 (see Table 2

Boundary quadrature on singular elements
Classical Gauss-Legendre quadrature (Gauß, 1815), modified Gauss-Legendre quadrature (Heltai et al. (2014) for biquadratic NURBS) and Duffy quadrature (Fairweather et al., 1979;Duffy, 1982) are investigated in Appendix B with respect to their suitability for approximation of weakly singular integrals. The resulting advantages and deficiencies of these three quadrature rules are summarized here. Fig. 4 shows that the relative quadrature error (16) decreases with increasing number of quadrature points n qp (26) for all three quadrature rules. However, their convergence behavior is very different: The error for classical and modified Gauss-Legendre quadrature with respect to n qp yield the rather low convergence rates of 1/2 and 1, respectively. 9 Fig. 4a shows that the relative quadrature error for modified Gauss-Legendre quadrature can be reduced to the range of machine precision by increasing the number of quadrature points to n qp > 10 15 , while classical Gauss-Legendre quadrature cannot reach machine precision as it has an inherent error bound of e sing rel > e class ≈ 3 × 10 −9 .
Duffy quadrature leads to a much lower error for small n qp and, most importantly, to a significantly improved convergence behavior: Fig. 4b shows that the relative error converges to the analytical solution with e rel ∝ 1/n qp 2.5 for small n qp ≤ 18 and with e rel ∝ 1/n qp 12 for higher n qp such that is in the range of machine precision for n qp > 200. Hence, Duffy quadrature is several orders of magnitude more accurate than modified Gauss-Legendre quadrature, not to mention classical Gauss-Legendre quadrature.
Further, collocation points and classical Gauss-Legendre quadrature points may coincide, in 7 Analytical solutions: Ising = log 12 √ 2 + 17 L e and Inear = 6 log( √ 2 + 1) + log 2 √ 2 + 3 L e 8 The color scheme from Fig. 3 is used throughout the paper to visualize the quadrature schemes: Blue, green and white are used to distinguish between singular, near singular and regular elements (see Fig. 6a for a 4x4 mesh). Various other colors, including magenta for Duffy quadrature, are used to illustrate the quadrature rules applied to these elements (see Fig. 7 for the same 4x4 mesh). 9 The minimum distance between quadrature points and collocation points (B.7) is inversely proportional to the number of quadrature points, i.e. rmin ∝ n −1 qp .
a.   7), the vertical dashed line shows the smallest admissible distance ε machine ≈ 2.2 × 10 −16 . The horizontal lines show that classical Gauss-Legendre quadrature has an inherent error bound e class ≈ 3 × 10 −9 , while modified Gauss-Legendre quadrature theoretically meets e mod ≈ ε machine . b. Relative quadrature error e sing rel for modified Gauss-Legendre quadrature and Duffy quadrature vs. the number of quadrature points n qp . which case the integral kernel becomes infinite and useless. Classical Gauss-Legendre quadrature is therefore less robust than modified Gauss-Legendre quadrature and Duffy quadrature, where the coincidence of collocation points and quadrature points is generally prevented.
For classical and modified Gauss-Legendre quadrature, it is sufficient to evaluate the shape functions only once on the master element for all collocation points and elements. 10 Gauss-Legendre quadrature is therefore simple to implement. Duffy quadrature, on the other hand, needs to be modified depending on the location of the collocation points (see Fig. 3b) and it is thus more difficult to implement. The above findings are summarized in Table 1. efficiency robustness implementation classical Gauss-Legendre poor low simple modified Gauss-Legendre low high simple Duffy quadrature high high involved In conclusion, Duffy quadrature is highly accurate and thus most recommend for singular integral approximation. Modified Gauss-Legendre quadrature is a robust alternative that is simple to implement but less accurate, while classical Gauss-Legendre quadrature is not recommend for the quadrature of singular kernels since it lacks robustness and efficiency.

Boundary quadrature on near singular elements
This section deals with the approximation of singular integrals on so-called near singular elements. The collocation point is located near these elements, but not on within them (green elements in Fig. 3a). The accuracy of Gauss-Legendre quadrature on near singular elements is investigated in Sec. 3.2.1, while a new and efficient quadrature rule for these elements is presented in Sec. 3.2.2. The quadrature rules discussed in Sec. 3.1, modified Gauss-Legendre quadrature and Duffy quadrature, are only beneficial on singular elements and thus not considered here.
3.2.1 Classical Gauss-Legendre quadrature Fig. 5a shows that the relative quadrature error (16) for classical Gauss-Legendre quadrature on near singular elements decreases with increasing number of classical Gauss-Legendre quadrature points until it within machine precision for n e qp = 16 × 16 quadrature points per element (ñ qp = 16). 11 The refinement of classical Gauss-Legendre quadrature is on near singular elements much more effective than on singular elements. However, determining an appropriate number of quadrature points per singular element requires a careful trade-off between accuracy and computational effort. Four variables that refer to the number of quadrature points are used in this paper. An overview of these variable is given in Table 2 and more details can be found in Appendix B.    Table 4 mean number of quadrature points Table 2: Boundary quadrature of singularities: Variables that refer to the number of quadrature points.

Gauss-Legendre quadrature with adjusted weights
As discussed in Sec. 3.2.1, a large number of Gauss-Legendre quadrature points is required to obtain accurate approximations of singular integrals on near singular elements. This section presents a new quadrature rule that yields high accuracy by adjusting the quadrature weights for singular kernels instead of increasing the number of quadrature points. The new quadrature rule is therefore referred to as Gauss-Legendre quadrature with adjusted weights.
The main idea is to use the locations of n e qp classical Gauss-Legendre quadrature points in the parameter space, denoted by ξ i (see Appendix B.1), and determine the corresponding quadrature weights w e i = w e (ξ i ) by solving a linear moment fitting problem such that 1/r kernels are integrated exactly. Quadrature rules based on moment fitting have been introduced by Joulaian et al. (2016), Thiagarajan and Shapiro (2016) and Hubrich et al. (2017) and applied to surface integration by Hubrich and Düster (2019) and Zou et al. (2021). Here, the quadrature weights on near singular element e have to fulfill for each shape function N A with A = 1, . . . , (p + 1)(q + 1). The local surface stretch at ξ i is denoted by J(ξ i ) such that da(ξ i ) ≈ J(ξ i ) w e i . 12 The integral kernel N A /r is proportional to the discrete BE kernels (A.9), since G ∝ 1/r and T n ∝ 1/r as shown in Fig. 1c. Quadrature weights that fulfill (17) on element e thus determine the singular BE integrals on this element exactly. The right hand side of (17) is determined to machine precision by classical Gauss-Legendre quadrature with 16 × 16 quadrature points (see Fig. 5a).
From (17) follows that n e qp = (p + 1) × (q + 1) quadrature points are required to obtain the following square moment fitting system of equations for near singular element e A e w e = f e , with subscripts A = 1, . . . , n e qp and B = 1, . . . , n e qp . The system of equations (19) is solved on each near singular element for the respective weight vector w e .
The proposed Gauss-Legendre quadrature rule with adjusted weights is illustrated by applying it to the biquadratic B-spline sheet from Fig. 3a. The adjusted quadrature weights are determined by solving (19) with n e qp = 9 for each of the eight near singular elements. Fig. 5b shows the difference between adjusted and classical Gauss-Legendre weights, defined on singular element e by It can be seen that small adjustments of the weights are sufficient to obtain integral approximations on near singular elements that are exact to machine precision. Due to symmetry, only two different sets of quadrature weight values occur on near singular elements: Larger weight differences on elements sharing one edge with the singular element (element type 1) and smaller weight differences on elements sharing one point with the singular element (element type 2). The near singular elements in Fig. 5c are colored to indicate elements of type 1 (blue) and elements of type 2 (red). The element type, and thus the quadrature weights, are represented in this way repeatedly throughout the remainder of this paper. The classical Gauss-Legendre weights (w GL 1 , . . . , w GL 9 ) and the quadrature weight differences for element type 1 (∆w 1 1 , . . . , ∆w 1 9 ) and for element type 2 (∆w 2 1 , . . . , ∆w 2 9 ) are given in Table C.1.
The same weight values are applied to all elements of the same element type, but the orientation of the weights may differ. The application of Gauss quadrature with adjusted weights therefore requires some additional bookkeeping. However, this is the only additional effort compared to classical Gauss-Legendre quadrature, since both quadrature rules use the same quadrature point locations. Shape functions and their derivatives can simply be determined a priori at the known locations ξ GL in the parameter domain as discussed in Appendix A. This property makes Gauss quadrature with adjusted weights highly suitable for hybrid quadrature based on classical Gauss quadrature, as will be seen in Sec. 4. Both set of weights are derived for a very specific case here, but they yield exact (to the range of machine precision) integrals on regular, quadrilateral and plane NURBS elements of biquadratic order and arbitrary size, as will be seen in Sec. 4.1. Moreover, Gauss quadrature with adjusted weights is also beneficial on curved surface, where it improves the quadrature accuracy compared to classical Gauss quadrature (see Sec. 4.2).

Hybrid quadrature schemes
The quadrature rules for singular and near singular elements from Sec. 3 are combined here to obtain hybrid quadrature schemes capable of efficient singular integral approximation on the entire surface. The hybrid quadrature schemes are then applied to plane B-spline surfaces in Sec The numerical computations in Secs. 4 and 5 are carried out in MATLAB using IEEE 754 double precision (53 bit with machine precision ε machine ≈ 2.2 × 10 −16 ).
14 Fig. 4 shows that increasing the number of Duffy quadrature points drastically reduces the quadrature error on the singular element. However, the overall error for DG, DGr and DGw is dominated by the remaining elements where classical Gauss-Legendre quadrature is applied (see Fig. 9 for a B-spline sheet). A higher number c. hybrid Duffy-Gauss quadrature with progressive refinement (DGr) singular element: Duffy quadrature withñ qp = 2n 0 near singular elements: classical Gauss-Legendre quadrature withñ qp = 2n 0 d. hybrid Duffy-Gauss quadrature with adjusted weights (DGw) singular element: Duffy quadrature withñ qp = 2n 0 near singular elements: Gauss-Legendre quadrature with adjusted weights (ñ qp = 3) In all cases, classical Gauss-Legendre quadrature withñ qp = n 0 is used on the regular elements (see Fig. 6 for singular, near singular and regular elements on a plane sheet). The quadrature density n 0 controls the total number of quadrature points n qp for the hybrid quadrature schemes. 15 It should be noted that even though Duffy quadrature can be applied to the entire surface to evaluate the singular integrals efficiently, this only gives accurate results for the entire surface but not for individual elements, as is required in BE analysis.
As described in Sec. 1, hybrid Duffy-Gauss quadrature (DG) has been used in BE analysis for Stokes flow problems (Barakat and Shaqfeh, 2018). Progressive refinement of Gauss-Legendre quadrature is additionally considered for acoustic problems (Venås and Kvamsdal, 2020). This approach resembles hybrid Duffy-Gauss quadrature with progressive refinement (DGr), but Venås and Kvamsdal (2020) do neither describe the refinement strategy nor investigate its influence on the BE result. Hybrid Duffy-Gauss quadrature with adjusted weights (DGw) is an entirely new scheme that is presented in this paper for the first time.
Each of the four schemes presented is applicable to arbitrary surface discretizations with quadrilateral elements. However, the investigation here focus on biquadratic isogeometric discretizations: The quadrature schemes are applied to a biquadratic B-spline sheet in Sec. 4.1 and to a biquadratic NURBS sphere in Sec. 4.2.

Hybrid quadrature on a B-spline sheet
The introduced hybrid quadrature schemes are applied to a biquadratic B-spline sheet in this section. The coarsest discretization level = 1 with 4 × 4 elements is shown in Fig. 6a, while Figs. 6b and c show the subsequent levels = 2, 3 created by successive knot insertion. Black circles represent the n no collocation points on the sheet. Among those, y 0 (red face) is chosen to illustrate the introduced quadrature schemes. The quadrature rules applied to the B-spline elements are discussed in Sec. 4.1.1. The accuracy of the hybrid quadrature schemes is then compared in Sec. 4.1.2.

Elemental quadrature
The four hybrid quadrature schemes are applied to the B-spline sheets from Fig. 6. Fig. 7 shows the elemental quadrature rules used for quadrature schemes G (a), DG (b), DGr (c) and DGw (d) on the B-spline sheet from of Duffy quadrature points would increase the computational effort without increasing the quadrature error and is therefore not considered here. 15 Variables defining the number of quadrature points for an element or for the entire discretization are summarized in Table 2. More details can be found in Appendix B.
Elements: singular near sing. regular Collocation points:  In BE analysis, the integration over the whole surface has to be performed for each collocation point once (see Appendix A.2). The quadrature weights for Gauss-Legendre quadrature with adjusted weights as well as the weights and locations of the quadrature points for modified Gauss-Legendre quadrature and for Duffy quadrature depend on the location of the collocation point. The introduced quadrature schemes thus require small modifications based on the location of the collocation point. The collocation points on a biquadratic NURBS sheet (black circles in Fig. 6) can be divided into three types based on their location: • type C: at a corner point of an element • type E: at the midpoint of an elemental edge • type M: at the midpoint of an element (e.g. y 0 ) The location of the Duffy quadrature points are shown in Fig. 3b for collocation points of type M and in Fig. B.4 for collocation points of type C and E. Quadrature points for modified Gauss-Legendre quadrature on a biquadratic element are shown in Fig. B.1b. This modification is beneficial in terms of robustness and efficiency for collocation points of type E and M, while classical Gauss-Legendre quadrature (Fig. B.1a) is more efficient for type C collocation points.
Gauss-Legendre quadrature with adjusted weights has been presented and discussed in Sec. 3.2.2 for collocation points of type M. The quadrature weights for collocation points of other types are obtained in the same manner. Fig. 8 shows the element types for DGw and exemplary collocation points of type C, E and M. It can be seen that only seven sets of quadrature weights are required for the integration over all elements and collocation points. These seven sets, which are required and sufficient for regular biquadratic B-spline sheets of arbitrary refinement, are given in Tables C.1 and C.2. The largest weight differences occur by far for element type 1 and 5. Therefore, a simpler but still promising approach is to apply adjusted weight quadrature only to elements of type 1 and 5 and classical Gauss quadrature to the remaining near singular elements.

Gauss Duffy
Gauss with adjusted weights 1 2 3 4 5 6 7 Figure 8: Hybrid quadrature on a B-spline sheet: Element type (above) and the corresponding weight differences (below) for quadrature scheme DGw and various exemplary collocation points (red dots).

Quadrature accuracy
The absolute quadrature error on element e due to evaluating the singular integral 16 at collocation point y 0 is defined as where I h e denotes the numerical approximation of integral (9). Fig. 9 shows the absolute quadrature error on B-spline discretizations of refinement level = 1, 2, 3 with n 0 = 3. It can be seen that the maximum elemental error decreases with increasing and decreases from scheme a. (G) to d. (DGw). The total quadrature error for y 0 , defined as with the elemental summations I = e I e and I h = e I h e , also decreases with increasing and from scheme a. to d. as Fig. 9 shows. The only exception occurs for scheme DGw from = 1 to = 2 (see Fig. 9d). The reason for the increasing error is to be found in the most inner layer of classical Gauss-Legendre elements (third layer in total) that yields the highest elemental errors: The sheets for = 2, 3 contain the entire third layer, while the = 1 sheet only contains half of that layer. The same effect can be observed for DGw, where the overall error decreases only by factor 1.3 from = 1 to = 2. For all other cases, e 0 abs is halved with each refinement step, i.e. linear convergence with convergence rate µ = 1/2.
on the investigated B-spline. Fig. 10a shows linear convergence with rate µ = 1/2 for each quadrature scheme. However, the error for quadrature scheme DG is about two orders of magnitude lower than for scheme G, while Duffy-Gauss quadrature scheme with progressive refinement (DGr) reduces the error by two further orders of magnitude. Duffy-Gauss quadrature with adjusted weights (DGw) is even more accurate and yields the lowest error among the investigated schemes.  Fig. 10b shows the mean relative error vs. the total number of quadrature points that is given by where n A, e qp denotes the number of quadrature points for collocation point y A on element e. The figure shows the same convergence behavior as noted in Fig. 10a. It also shows that DGr uses slightly more quadrature points per refinement level than schemes DGw, DG and G. However, the gain in accuracy clearly outweighs the additional quadrature effort as Fig. 10b shows. Therefore DGr and especially DGw are the most efficient quadrature schemes for singular integral approximation on biquadratic B-spline sheets. Fig. 10 further shows that increasing the number of quadrature points or adjusting the quadrature weights in one ring of elements (DGr or DGw) leads to a beneficial jump in accuracy (c.f. DG). Using further rings, additional, albeit smaller jumps in accuracy can be expected.

Hybrid quadrature on a NURBS sphere
Next, the four hybrid quadrature schemes are applied to spherical surfaces to investigate their suitability for curved NURBS surfaces. The procedure advocated here is visualized and investigated by considering biquadratic NURBS spheres, but is applicable to any closed and open surface defined by a mapping from a planar parameter domain. It should be further noted that any genus zero surface can be mapped conformally onto a sphere (Gu et al., 2004), so this example has far reaching applications. Two different spherical discretizations with biquadratic NURBS elements are considered here: The single-patch NURBS sphere from Fig. 11a and b is investigated in Sec. 4.2.1, while the six-patch NURBS sphere from Fig. 11c and d is investigated in Sec. 4.2.2. Fig. 11e shows that the single-patch discretization is exactly spherical independent of , whereas the six-patch sphere is only approximately spherical. For the six-patch sphere, the L 2 norm of the radius error 17 shows almost quadratic convergence w.r.t. the number of control points (e L 2 R ∝ 1/n no 2 ).
a. b.
∝ 1 / n n o 2 ∝ 1 / n n o Figure 11: Hybrid quadrature on a NURBS sphere: Single-patch NURBS sphere of discretization level = 1 (a) and = 2 (b) and six-patch NURBS sphere of level = 1 (c) and = 2 (d). Collocation points are located on surface S h = S, while control points are not. e. L 2 norm of the radius error for = 1, . . . , 6.
The number of elements and control points are given in Table 3 for the single-patch and the six patch sphere, both of refinement level = 1, . . . , 6. The corresponding number of quadrature points are given in Table 4 for the four hybrid quadrature schemes with quadrature density n 0 = 3. Those numbers also apply to the numerical examples in Sec. 5, which consider spherical and spheroidal geometries.
The accuracy of the quadrature schemes is investigated by approximating the singular integrals of the BIE (12). The quadrature error for collocation point y A is defined with respect to identity (9) as 17 The L 2 norm of the radius error eR = (R h − R)/R is defined as e L 2 R := S e 2 R da/AS , where AS = 4πR 2 denotes the surface area of a sphere. and with respect to identity (10) as where r A := x − y A . The mean quadrature errors with respect to identities (9) and (10) Table 4: Hybrid quadrature on a NURBS sphere: The mean number of quadrature pointsn qp = n qp /n no (rounded) for NURBS spheres of discretization level = 1, . . . , 6 (red font: single-patch sphere, blue font: six-patch sphere) and quadrature density n 0 = 3.

Single-patch NURBS sphere
Single-patch NURBS surfaces are constructed by revolving a NURBS semicircle by a full rotation (Piegl and Tiller, 1997) and they are thus exactly spherical. Fig. 11a shows a coarse NURBS discretization of the surface (refinement level = 1) and the corresponding control points and collocation points. The sphere has constant radius (C 0 -continuity), continuous tangent vectors (C 1 -continuity) and constant curvature, making it C 2 -continuous everywhere. However, these continuity properties are in general not maintained during deformation: The single-patch NURBS sphere is interpolatory across the patch boundary (C −1 -continuity) and between the other octants of the sphere (C 0 -continuity), while it is C p−1 -continuous across the remaining element boundaries. Fig. 11b shows a NURBS discretization of refinement level = 2 obtained by knot insertion (Hughes et al., 2005) maintaining the continuity properties from = 1. weights (DGw) and exemplary collocation points located within elements. Particular attention has to be paid to collocation points within the degenerated elements located at the poles of the sphere (Fig. 12b). For these points, Gauss quadrature with adjusted weights cannot be applied to the near singular elements in a meaningful way. Therefore, refined Gauss-Legendre quadrature is applied instead to the elements adjacent to the pole (except the singular element where Duffy quadrature is applied). In contrast, the application of DGw to collocation points within other elements (Fig. 12a) follows directly from Sec. 4.1 and is thus straightforward. Collocation points on element boundaries only occur along the lines of reduced continuity between octants of the sphere (see Fig. 11a and b). DGw is applied to collocation points between two (Fig. 12c) and four octants (Fig. 12d) by treating each octant of the sphere separately as discussed in Sec. 4.1. For a collocation point at the pole (Fig. 12e), Duffy quadrature is applied to the 2 +2 elements adjacent to the collocation point.   Fig. 13 shows the absolute element-wise quadrature error (23) for hybrid quadrature schemes G (a.), DG (b.), DGr (c.) and DGw (d.) considering collocation point y A that is located within an element. 18 It can be seen that scheme G is very inaccurate on the singular element, while the elemental error decreases with increasing distance to y A . Using scheme DG instead reduces the quadrature error on the singular element drastically so that the highest errors now occur on the near singular elements. Scheme DGr reduces the error on those elements so that the singular integral is approximated on all elements accurately. Gauss-Legendre quadrature with adjusted weights is designed in Sec. 3.2.2 to determine singular integrals on near singular elements of plane and regular sheets exactly. Although integral (22) is not approximated exactly on the near singular elements of curved surfaces, Fig. 13d shows that scheme DGw yields a lower error than scheme DG for the same number of quadrature points.
The mean quadrature error (29) with respect to the BE identities (9) and (10) is shown in Fig. 14 vs. the number of total quadrature points (26) for mesh refinement level = 2 and varying quadrature refinement n 0 = 2, 4, 8, 16, 32. Both errors decrease with increasing quadrature refinement for quadrature schemes G, DG and DGr,19 where DG and especially DGr shows a much faster convergence than G: With n 0 = 32, the quadrature error with respect to the first identity (see Fig. 14a) is e SL ≈ 10 −3 for scheme G, e SL ≈ 3 · 10 −9 for scheme DG and e SL ≈ 6 · 10 −14 for scheme DGr. The error with respect to the second identity (see Fig. 14b) shows a similar behavior, albeit the difference between thee three schemes is slightly smaller. For both errors, schemes DG and DGr show a much better convergence behavior than scheme G whose convergence rate is linear.
The singular BE integrals are approximated with high efficiency using a coarse discretization and a moderate number of quadrature points as shown in the previous paragraph. However, refined discretizations are required to represent high order boundary conditions or more complex surface geometries. Another application that requires sufficiently small elements are coupled FE-BE simulations, where FE and BE analysis is conducted on a deforming surface. The accuracy of the quadrature schemes on refined discretization is thus of great interest for BE analysis and is investigated below. Fig. 15 shows the quadrature error w.r.t. identity (9) for all collocation points on one octant of a sphere of refinement level = 1 (a.), = 2 (b.), = 3 (c.) and = 4 (d.) using quadrature scheme DGr with n 0 = 3. The quadrature error decreases for the majority of the collocation points with increasing quadrature refinement. However, the collocation points next to the pole yield similarly high errors independently of . The mean quadrature error (29) e y A SL a. b. c. d. considering all collocation points is thus not monotonically decreasing as can be seen in Fig. 16: The quadrature error for DGr increases for the first refinement steps before it starts to decrease at = 3. Scheme G results in a much higher error that, however, shows linear convergence with an almost constant rate µ = 0.5. The curves for scheme DG and scheme DGw lie between those for G and DGr, with DGw being more accurate than DG. For higher , all investigated quadrature scheme show linear convergence with the same rate, where DGr is more accurate than DGw, which is more accurate than DG, which is in turn more accurate than G. a.

Six-patch NURBS sphere
This section discusses the application of the four hybrid quadrature schemes to the discretizations shown in Fig. 11c and d. These discretizations consist of six biquadratic NURBS patches and are thus referred to as six-patch NURBS spheres. In contrast to the single-patch discretziation from Sec. 4.2.2, the six-patch NURBS sphere is only approximately spherical (see Fig. 11e). The application of hybrid quadrature schemes to six-patch NURBS spheres is nevertheless promising, since there are no degenerated and small elements near the pole, which cause the largest quadrature errors on the single-patch sphere (cf. Fig. 15). It should be further noted that all collocation points lying on element boundaries also lie on patch boundaries as Fig. 11c and d shows.
The quadrature rules and the element types for hybrid Duffy-Gauss quadrature with adjusted weights (DGw) are shown in Fig. 17a and b for collocation points within patches, in Fig. 17c for collocation points at junctions of three patches and in Fig. 17d and e for collocation points at boundaries between two patches. Particular attention has to be paid to collocation points located near, but not directly at, junctions of three patches (see Fig. 17b and e). For those collocation points, the number of near singular elements is reduced by one (cf. Fig. 17a and d).

Gauss Duffy
Gauss with adjusted weights n 0 2n 0 Figure 17: Hybrid quadrature on a six-patch NURBS sphere: Quadrature rules used for hybrid Duffy-Gauss quadrature with adjusted weights (DGw) considering collocation points that are located within one patch (a. and b.), at the junction of three patches (c.) and along the common edge of two patches (c. and d.).
The mean quadrature error (29) with respect to BE identities (9) and (10) is shown in Fig. 18 vs. the number of total quadrature points (26) for mesh refinement level = 2 and varying quadrature refinement n 0 = 2, 4, 8, 16. Similar to the single-patch sphere from Sec. 4.2.1, the hybrid quadrature schemes G, DG, and DGr show decreasing quadrature errors as the number of quadrature points increases. The convergence behavior of quadrature schemes that use Duffy quadrature is even better for the six-patch sphere than for the single-patch sphere so that DG and DGr yield mean quadrature errors in the range of machine precision for n 0 = 16, respectively for n 0 = 8. a. b.   quadrature scheme DGr with n 0 = 3. It can be seen that the quadrature error decreases with increasing quadrature refinement for all collocation points, which was not the case for the singlepatch sphere (cf. Fig. 15). The mean quadrature error (29) with respect to identities (9) and (10) is shown in Fig. 20 for n 0 = 3 and varying mesh refinement = 1, . . . , 5. Both quadrature errors decrease monotonically with increasing mesh refinement for each of the four quadrature scheme. DGr provides by far the best result, followed by DGw and DG and eventually G, which is still quite inaccurate even for highly refined meshes. In contrast to e SL (Fig. 20a), DGw does not yield a significant improvement over DG for e DL (Fig. 20b). a.

Conclusions on hybrid quadrature
The suitability of the presented hybrid quadrature schemes for BE analysis has been demonstrated in the previous sections: Sec. 4.1 shows that Duffy-Gauss quadrature with progressive refinement (DGr) and especially hybrid Duffy-Gauss quadrature with adjusted weights (DGw) approximate singular integrals on flat surfaces very efficiently, while Sec. 4.2 shows that DGr is most efficient on curved surfaces.
It has been further shown that the hybrid quadrature schemes approximate the singular BE integrals very efficiently on coarse single-patch NURBS spheres (Sec. 4.2.1). These discretizations are exactly spherical and therefore recommended for BE analysis on rigid spheres or spheroids.
The application of hybrid quadrature schemes to the six-patch NURBS sphere (Sec. 4.2.2) also yields accurate results for coarse discretizations, but moreover a significantly improved convergence behavior for mesh refinement. The six-patch discretization is therefore recommended for more complex surface geometries and deforming surfaces.

Application to Stokes flow problems
The suitability of the introduced quadrature schemes for BE analysis is investigated here with three numerical examples: Sec. 5.1 and Sec. 5.2 consider spheres rotating in a viscous fluid and translating through a viscous fluid, respectively. 20 A rising bubble of non-spherical shape is considered as a third problem in Sec. 5.3. In all three problems, the velocity on the surface is given by Dirichlet boundary conditions and the traction on the surface is determined by BE analysis. The numerical traction results are compared to analytical results based on Chwang and Wu (1975) and Kong et al. (2012).

Flow caused by a rotating sphere
The first example considers a rigid sphere (surface S, Radius R) that is surrounded by a viscous fluid F (dynamic viscosity η). The sphere rotates around its center with the prescribed angular velocityω ∈ R 3 , so that the surface velocity is given by A rotation around the vertical axis is considered here, i.e.ω =ω e 3 . The BIE can be solved for the surface traction t := σn considering (30) as a Dirichlet boundary condition. Due to the pure rotational velocity, the traction field has only tangential components, while the normal traction (i.e. the pressure) is zero.
The BE traction error with the maximum traction value located along the equator of the sphere, is introduced to compare the numerical results systematically to the analytical solution. Fig. 21 shows the BE traction error on a single-patch NURBS sphere and on a six-patch NURBS sphere, both of mesh refinement level = 4, considering the introduced quadrature schemes from Sec. 4. Half of the symmetric BE meshes is hidden to improve the visibility of the results.
The accuracy of the introduced hybrid quadrature schemes is investigated systematically in two convergence studies: First, the influence of quadrature refinement is investigated by varying the quadrature density n 0 = 2, 4, 8, 16 on a fixed mesh ( = 2). Second, the influence of mesh refinement is investigated by varying the mesh refinement level = 1, . . . , 6, while n 0 = 3. Both convergence studies consider single-patch and six-patch NURBS spheres. The traction error on the surface is characterized by the L 2 norm of the traction error defined as 20 A fixed sphere in a steady rotational flow leads to a mathematically equivalent problem as a sphere rotating in a quiescent fluid. A fixed sphere in a steady transversal flow leads to a mathematically equivalent problem as a sphere translating through a quiescent fluid.  where A S = 4πR 2 denotes the surface area of S. Fig. 22 depicts e L 2 t for the first convergence study. On the single-patch sphere ( Fig. 22a and  b), the traction error decreases for each of the hybrid quadrature schemes, where hybrid Duffy-Gauss quadrature with progressive refinement (DGr) provides the best result, followed by hybrid Duffy-Gauss quadrature (DG) and eventually by Gauss-Legendre quadrature (G). Black markers denote results obtained from an iterative solver. 21 This convention is used for all examples in Sec. 5. Similar results are obtained for the six-patch sphere ( Fig. 22c and d), but the L 2 norm of the traction error converges only to e L 2 t ≈ 10 −3 for = 2 and to e L 2 t ≈ 5 · 10 −5 for = 4. Better results are prevented by the inaccuracy of the six-patch sphere discretization (see Fig. 11e). The six-patch sphere does not require an iterative solver and is thus preferable in terms of robustness and computational effort. The quadrature weights for adjusted weight quadrature are given in Table C.1 and C.2 only for n qp = 3 × 3. The results for DGw are thus not investigated for a varying quadrature density. Fig. 23 depicts the L 2 norm of the traction error for the second convergence study. It can be seen that the traction error decreases for each of the hybrid quadrature schemes and on both of the considered discretizations. As expected, scheme DGr and G provide the best and the worst results, respectively. The results of the remaining two schemes are in between, where DGw seems to offer no significant advantage in accuracy over DG here. Applying G or DG to a 21 The iterative solution method used is the preconditioned conjugate gradients method. For each individual combination of quadrature scheme, mesh refinement and quadrature density, the depicted error is chosen based on where e L 2 t,iter and e L 2 t,direct denote the L 2 norm of the traction error using an iterative or a direct solver, respectively. The iterative solver is more robust w.r.t. ill-conditioned BE matrices and yields thus better results on highly refined single-patch spheres with coinciding control points at the poles. c. d. single-patch sphere with > 1 requires an iterative solver, as indicated by the black markers in Fig. 23a. Schemes DGw and DGr, which also treat the near singular elements properly, allow to use a direct solver for = 2, 3 also, and are therefore preferable in regard to the computational time.
Unlike the single-patch sphere, the six-patch sphere (Fig. 23b) does not require an iterative solver and is therefore very robust, both with respect to mesh and quadrature refinement. The six-patch sphere further leads to a smaller error, where DGr provides by far the best result with a convergence rate of almost 1.5 (e L 2 t ∝ 1/n 1.5 no ). The kink in the curves for DG and DGw shows that applying Duffy quadrature to the singular elements without increasing the quadrature density on the near singular elements is insufficient to yield accurate results for > 3.

Flow caused by a rising sphere
The suitability of the introduced quadrature schemes for a pure rotational problem has been shown in Sec. 5.1. However, this problem yields only tangential tractions, while the normal tractions are zero. The suitability of the quadrature schemes for problems with non-zero velocity and traction components in both normal and tangential direction is investigated in the second for quadrature density n 0 = 3 on the single-patch NURBS sphere (a) and on the six-patch NURBS sphere (b), both of varying refinement level = 1, . . . , 6. The black markers depict results from an iterative solver. example: A sphere with radius R is transtaled with constant velocityv =v e 3 22 through the fluid. The Dirichlet boundary condition on the surface S is therefore given by while the surface traction is unknown. Fig. 24 shows the L 2 norm of the relative BE traction error for translating spheres of refinement level = 1, . . . , 6. The results are generally similar to those from Sec. 5.1. On the single-patch sphere, the translating sphere problem can be solved with higher accuracy than the rotating sphere problem. On the six-patch sphere, it is the other way around. It is noteworthy that the convergence rate for DGr and > 4 decreases slightly.
Treating additional rings with refined quadrature would allow to maintain a convergence rate of 1.5. However, a single refinement ring already gives an excellent gain in accuracy.

Flow caused by a rising ellipsoid
The suitability of the quadrature schemes for problems with non-spherical surface geometries is investigated in the third example: A ellipsoid with surface S rises with constant velocityv through a fluid of dynamic viscosity η. The surface of the ellipsoid is described by where e denotes the eccentricity of the ellipsoid (0 ≤ e < 1) and a is the length of its semi-major axis. 23 The velocity on S is given by the Dirichlet boundary condition (34). 22 The prescribed velocity can be chosen arbitrarily, i.e.v =vi ei. Here, a pure vertical velocity is chosen for the sake of simplicity. 23 The eccentricity of a ellipsoid is defined as e := 1 − b 2 /a 2 , where b denotes the length of the semi-minor axis. Two parameters out of the triplet (a,b,e) needs to be known to define the ellipsoidal surface. The example from Sec. 5.3 considers a ellipsoid with e = 0.75 and a = L/e 1/3 . Figure 24: Rising sphere: Traction error e L 2 t for quadrature density n 0 = 3 on the single-patch NURBS sphere (a) and on the six-patch NURBS sphere (b), both of varying refinement level = 1, . . . , 6. The black markers depict results from an iterative solver. Fig. 25 shows the BE traction on a six-patch ellipsoid of refinement level = 3 using quadrature scheme DGr with n 0 = 3. The normal and tangential components of the traction vector and its magnitude are shown in Fig. 25a, b and c, respectively. Fig. 25d shows that the relative BE traction error (35) yields e t ≈ 10 −3 along patch boundaries and an even much smaller error away from patch boundaries. The velocity and pressure field on the ellipsoid and in the surrounding fluid are shown in Fig. 26. The fluid velocity is determined by the BIE for points within the domain (11), while the surface velocity on S is given by BC (34). The pressure is determined on the ellipsoid surface by p = −t · n and within the fluid domain by the pressure BIE (A.11).  The L 2 norm (33) of the relative BE traction error (35) is shown in Fig. 27a vs. the number of nodes, and in Fig. 27b vs. the total number of quadrature points (26). The results are very similar to those of the rising sphere in Sec. 5.2, since the sphere is a special case of the ellipsoid (for e = 0): Quadrature schemes DGr and G provide the best and the worst of the BE results, respectively. Considering DGr with additional rings of refined quadrature would allow to maintain a convergence rate of 1.5. Quadrature schemes DG and DGw provide almost the same results, which are in between the results of the former. Thus, all hybrid quadrature schemes, in particular DGr, prove to be robust to changes in surface shape and are therefore suitable for various applications, including coupled FE-BE simulations.   The L 2 norm of the relative pressure error e p for quadrature density n 0 = 3 and six-patch NURBS ellipsoids of varying refinement level = 1, . . . , 6.

Conclusion and outlook
This work presents new quadrature schemes for the efficient approximation of weakly singular integrals (1/r kernel) and provides important findings on singular quadrature in BE (Boundary element) analysis. A new quadrature rule that approximates weakly singular integrals on plane and regular near singular elements exactly is presented in Sec. 3.2.2. Numerical examples have shown that this quadrature rule can also be advantageous on curved surfaces, as it improves the accuracy compared to classical Gaussian quadrature. The presented quadrature rule can be easily extended to determine singular integrals exactly on other surfaces like spheres and cylinders.
Sec. 4 presents four new hybrid quadrature schemes for BE analysis that combine different quadrature rules for singular, near singular, and regular elements. Numerical investigations on flat and curved isogeometric surfaces show that all presented schemes converge robustly for both quadrature and mesh refinement, with hybrid Duffy-Gauss quadrature with progressive refinement (DGr) being by far the most efficient scheme. The present paper considers a single ring of elements with refined quadrature, which gives already an excellent gain in accuracy. How-ever, considering additional rings would yield even better approximations and is thus worth to investigate in future works. The quadrature resolution for each ring could then be determined with respect to numerical criteria such as those of Bu and Davies (1995).
The application of the presented hybrid quadrature schemes to coupled FE-BE analysis of fluid structure interaction problems is planned in a forthcoming paper. The presented schemes also offer a wide range of other applications, since they are suitable for arbitrary weakly singular integrals and are not limited to Stokes flow or to BE analysis at all.

A Boundary element analysis
Solving the boundary integral equation (BIE) for Dirichlet, Neumann or mixed problems by use of boundary elements (BE) is discussed in this section. The spatial discretization of the surface geometry and the BIE (12) is introduced in Sec. A.1, while the subsequent collocation into a solvable system is discussed in Sec. A.2. The elemental mapping from the physical domain to a parameter domain is described in Sec. A.3 and the conceptual approximation of surface integrals by numerical quadrature is presented in Sec. A.4.

A.1 Boundary element discretization
The surface geometry S and the BIE (12) are discretized into n el finite boundary elements, numbered e = 1, . . . , n el , and n no nodal points. Element e occupies the surface domain Ω e ⊂ S h such that the surface geometry is approximated by where superscript h denotes approximated quantities. Point x ∈ Ω e is approximated by the nodal interpolation where x I denote one of the n e nodal points that define element e and N I denotes the corresponding shape functions 24 (e.g. Lagrange functions, B-splines or NURBS). Discrete elemental arrays are defined to simplify the notation: The nodal positions for element e are assembled in the discrete column vector of length n e dof := 3 n e whereas the nodal shape functions are assembled in the (3 × n e dof ) array N e = N 1 1, N 2 1, · · · , N ne 1 , (A.4) where 1 is the usual identity tensor in R 3 . Approximation (A.2) can be re-written in matrix form by use of the elemental position vectors as where v e and t e denote the nodal velocity and traction vector of length n e dof for element e. At source point y ∈ Ωē, the velocity is analogously approximated by v(y) ≈ Nē vē , (A.7) where Nē and vē are the shape function array and the nodal velocity vector for elementē. The components of the discrete vectors v e ,v e and t e are denoted by v e A , vē A and t e A , respectively, while the components of the shape function arrays N e and Nē are denoted by N e iA and Nē iA . The lowercase indices i, j and k run from 1 to 3, while the uppercase index A runs from 1 to n e dof .
The BIE (12) is then approximated by use of (A.6) and (A.7) as where summation over A is implied. The components of the discrete BE matrices G e and T e , both of size (3 × n e dof ), are given by The elemental vectors v e and t e are assembled into global knot vectors v and t of length n dof := 3 n n , where n n denotes the total number of nodal discretization points, i.e. nodes for Lagrange and control points for isogeometric discretizations. Similarly, the elemental arrays G e , T e and Nē are assembled into global global matrices G, T and N of size (3 × n dof ). The discretized boundary integral equation for a source point y on the interface S can thus be given in matrix-vector form by A continuous boundary integral representation for the pressure is given in Pozrikidis (1992) by with pressure Green's function (A.12) and the corresponding tensor for pressure field associated with the stresslet More details on theory, discretization, assembly and implementation of BIEs can be found in Harmel (2022).

A.2 Boundary element collocation
The discretized BIE is then collocated at n no source points y A ∈ S, for A = 1, . . . , n no , to obtain a square system. The source points y A are thus referred to as collocation points in the following.
For Lagrange discretizations on the one hand, the collocation points are simply selected to be the same as the nodes located on the surface. For isogeometric discretizations on the other hand, the control points are not necessarily located on the surface, which makes then unsuitable for collocation. Therefore, the locations of the collocation points are determined by the Greville abscissae (see e.g. Greville (1964), Johnson (2005) and Auricchio et al. (2010)). Fig. 11a and b depict the resulting collocation points on spheres with biquadratic NURBS discretizations.
The discretized BIE (A.10) is evaluated for each collocation point to obtain the discrete system (A.14) where N BE , T BE and G BE are square matrices of size (n dof ×n dof ). Evaluating the BIE for collocation point y A fills three rows of the matrices, indicated by vector where now A = 1, . . . , n no and N , G and T as defined in (A.7) and (A.9). The subscripts on the left hand side of (A.15) are given in MATLAB-like notation so that N BE (d 1 , :) , for example, refers to the first three rows of N BE . The BE system (A.14) is solvable for Dirichlet, Neumann and mixed problems. However, only Dirichlet problems are considered in the numerical examples in Sec. 5.
Remark: The probably most intuitive way to create the BE system is to evaluate the discretized BIE (A.10) for all collocation points successively. However, changing the order of the loops, i.e. approximating the integrals on element e for all collocation points y A and repeating this step for the remaining elements subsequently, is computationally much more efficient.

A.3 Mapping to a parameter domain
The surface S is fully characterized by the parametric description where ξ α with α = 1, 2 are curvilinear coordinates associated with a parameter domain P. The mapping (A.16) reflects the property that the surface is a two-dimensional object embedded within three-dimensional space. Each elemental surface Ω e is mapped to a quadrilateral master element in the parameter domain ξ α ∈ [−1, 1] of side length 2. Comparing approximations (A.2) and (A.5) with the parametric mapping (A.16) shows that the nodal shape functions are defined on the master element, i.e. N I = N I (ξ 1 , ξ 2 ) . (A.17) This paper considers quadrilateral elements since these can be conveniently related to the master element introduced above. The reader is referred to Sauer (2018) for more information on mapping (A.16) and the corresponding surface description in curvilinear coordinates.

A.4 Boundary quadrature
The boundary integrals are first mapped to the parameter space with mapping (A.16). Each elemental integral is defined on the master element as where k(x) denotes an arbitrary integral kernel and J a denotes the local surface stretch between surface S and parameter domain P . The right hand side of (A.18) is then approximated with numerical quadrature rules (Gauß, 1815;Golub and Welsch, 1969;Laurie, 2001) that are stated as weighted sums of function values at specified positions in the parameter space. Bivariate quadrature rules with n e qp quadrature points for the approximation of surface integrals on the master element are thus defined as where f (ξ) denotes the integral kernel on the master element, while ξ i and w(ξ i ) denote the position and weight of quadrature point i. Various bivariate quadrature rules are investigated in Sec. 3 with respect to approximation of singular integrals.

B Existing boundary quadrature rules
The accurate approximation of the singular boundary integrals from (A.14) is crucial in BE analysis. Various quadrature approaches are therefore investigated in this section with respect to their suitability for singular integral approximation: The classical Gauss-Legendre quadrature rule is considered in Sec. B.1, a modified Gauss-Legendre quadrature rule is introduced in Sec. B.2, while Sec. B.3 considers a Duffy transformation-based quadrature rule for singular integrals.

B.1 Classical Gauss-Legendre quadrature
Gauss-Legendre quadrature has been introduced by Gauß (1815) and is nowadays the most common and established quadrature rule for numerical integration. Bipolynomial kernels of orders p and q are integrated exactly by Gauss-Legendre quadrature with a minimum number of n e qp = (p + 1)/2 × (q + 1)/2 (B.1) quadrature points per element, while non-polynomials kernels are only integrated approximately. This work focus on bivariate Gauss-Legendre quadrature rules with n e qp =ñ qp ×ñ qp without loss of generality. The quadrature point locations are thus denoted by for i := j + (k − 1)ñ qp , where j and k run from 1 toñ qp . The scalar valuesξ j andξ k denote the j-th and the k-th root of the Legendre polynomial Pñ qp (ξ), respectively. They can be determined with numerical methods like Newton's method or by exploiting explicit expressions (Golub and Welsch, 1969) or tables (Laurie, 2001). The corresponding weights are given by  The application of classical Gauss-Legendre quadrature to regular, i.e. non-singular, kernels is accurate and robust (Chawla, 1968;Kambo, 1970;Leone et al., 1979). However, the numerical integration of singular kernels is much more challenging. The weakly singular integral  decreases for classical Gauss-Legendre quadrature (blue line) with increasing number of quadrature points, i.e.ñ qp = 2 1 , 2 2 , . . . , 2 17 , where I h A denotes the numerical approximation of integral (B.5). However, the convergence rate is only 1/2, so the highest number of quadrature points n e qp = 2 17 × 2 17 ≈ 1.72 × 10 10 still results in e rel ≈ 3 × 10 −6 . Fig. B.2b depicts the mean quadrature error vs. the minimum distance between quadrature points and collocation points that is defined by where x(ξ i ) denote the position on S for quadrature point i. The vertical dashed line shows ε machine ≈ 2.2 × 10 −16 , which is the smallest computationally admissible value for r min . Extending the curve by extrapolating the simulation data shows that the mean error of classical Gauss-Legendre quadrature is bounded by e rel > e class ≈ 3 × 10 −9 , even with unlimited memory capacity. A huge number of n qp > 10 16 quadrature points would be required to obtain such a small error. It should be further noted that classical Gauss-Legendre quadrature does in general not prevent the coincidence of collocation points and quadrature points. For instance, using an odd number forñ qp results in a quadrature point located at the middle of the element, where in many cases a collocation point is also located. The approximation of the integral I 3 (B.5) is infinite for these cases and thus useless. 26 In conclusion, three main drawbacks for singular BE integral approximation with classical Gauss-Legendre quadrature can be identified: The quadrature rule -is very inefficient with respect to the number of quadrature points, -has a lower error bound of e class ≈ 3 × 10 −9 , -does not prevent the coincidence of collocation points and quadrature points.
On the other hand, classical Gauss-Legendre quadrature is + simple to implement and computationally efficient. 27 The disadvantages in accuracy and robustness clearly outweigh the advantageous implementation, making the classical Gauss-Legendre quadrature unsuitable for the approximation of singular BE integrals.

B.2 Modified Gauss-Legendre quadrature
The classical Gauss-Legendre rule can be modified to overcome the discussed deficiencies, while maintaining the efficient equal treatment of all elements. This modification is briefly introduced by Heltai et al. (2014) for biquadratic NURBS and is generalized here for discretizations with Lagrange or isogeometric basis functions of arbitrary order. The key idea of the modification is to split the surface elements into virtual sub-elements such that collocation points are exclusively located at the corners of sub-elements. Classical Gauss-Legendre quadrature is then 26 The coincidence of collocation points and quadrature points potentially also occurs for evenñqp. The collocation point locations depend on the order of the shape functions used for discretization. 27 The shape function values and the corresponding weights are determined once at the quadrature points on the master element. Isogeometric shape function values are determined from universal Bernstein polynomials and element specific Bézier extraction operators (Borden et al., 2011). The introduced modified Gauss-Legendre quadrature is used to approximate the singular integral (B.5) for the same collocation points y 0 1 , y 0 2 and y 0 3 as in Sec. B.1 to compare the accuracy of both quadrature rules. The mean quadrature error for modified Gauss-Legendre quadrature is illustrated by red lines in Fig. B.2. It yields convergence rates of 1 with respect to the number of quadrature points (a) and with respect to the minimum distance between collocation points and quadrature points (b). It can be seen that modified Gauss-Legendre quadrature shows a twice higher convergence rate than the classical counterpart. The mean quadrature error for modified Gauss-Legendre quadrature can thus be reduced to the range of machine precision by increasing the number of quadrature points to n e qp > 10 15 .
Summarizing, the modified Gauss-Legendre quadrature rule eliminates the first and second drawback of the classical one. Also the third drawback has been addressed, as the required number of quadrature points to achieve a certain precision has been reduced due to the better convergence rate. However, a convergence rate of 1 still has additional potential for improvement. Even the modified quadrature rule still requires at least n e qp = 1,000 to obtain e rel ≤ 10 −3 . Besides addressing the drawbacks, the modified Gauss-Legendre quadrature also maintains the desirable property from the classical Gauss-Legendre quadrature: It is sufficient to evaluate the shape functions once on the master element and reuse the values for all elements independently of the collocation point position.
In conclusion, the presented modified Gauss-Legendre quadrature • shows moderate efficiency with respect to the number of quadrature points, + leads to a vanishing error with progressing quadrature refinement, + prevents the coincidence of collocation points and quadrature points, + is simple to implement and computationally efficient.
The numerical investigation of the modified Gauss-Legendre quadrature shows a strong robustness and a moderate accuracy making it suitable for the approximation of singular integrals. A special quadrature rule for singular functions is presented in Sec. B.3 to achieve even more accurate integral approximations.

B.3 Duffy transformation-based quadrature
The introduced Gauss-Legendre quadrature rule comes along with significant accuracy shortcomings as discussed in the previous two sections. The quadrature approach from Fairweather et al. (1979) and Duffy (1982) is considered here to overcome this shortcomings. This approach exploits Duffy transformation from a triangle to a square to remove corner singularities of the type 1/r. The quadrature rule is referred to as Duffy transformation-based quadrature or as Duffy quadrature in short form. While Duffy (1982) and Mousavi and Sukumar (2010a) use Duffy quadrature for integral approximation on triangular elements, two Duffy triangles are combined here to define the quadrature rule on the master element. The Duffy quadrature is therefore applicable to arbitrary quadrilateral elements. The singular integral (B.5) is approximated using Duffy quadrature with an increasing number of quadrature points to compare its accuracy to Gauss-Legendre quadrature. Fig. 4b shows that a coarse Duffy quadrature with n e qp = 2 yields a similar accuracy as a refined modified Gauss-Legendre quadrature with n e qp = 25. The Duffy quadrature also shows a much better convergence behavior: The relative error converges to the analytical solution with e rel ∝ 1/n e qp 2.5 for small n e qp ≤ 18 and with e rel ∝ 1/n e qp 12 for higher n e qp such that is in the range of machine precision for n e qp > 200.
In conclusion, the introduced Duffy transformation-based quadrature + shows outstanding efficiency with respect to the number of quadrature points, + provides a perfect approximation already for moderate quadrature refinement, + prevents the coincidence of collocation points and quadrature points.
Duffy quadrature is several orders of magnitude more accurate than modified Gauss-Legendre quadrature, not to mention classical Gauss-Legendre quadrature, and is thus highly suitable for the approximation of strongly singular integrals. However, the Duffy quadrature presented is only directly applicable to corner collocation points and needs to be modified for other locations (see Fig. B.4). Since quadrature points and weights depend on the particular collocation point, Duffy quadrature -is more complex to implement than Gauss-Legendre quadrature.
Three different quadrature approaches are investigated in Sec. 3 with respect to their suitability for singular integral approximation. First, classical Gauss-Legendre quadrature is not robust and very inefficient and is thus not recommend for the quadrature of singular kernels. Second, modified Gauss-Legendre quadrature is very robust, simple to implement and yields moderate accuracy. As for classical Gauss-Legendre quadrature, it is sufficient to evaluate the shape functions only once on the master element for all collocation points and elements. Third, Duffy quadrature approximates singular integrals very efficiently, but the quadrature points and a.
b. c. weights depend on the location of the collocation point. However, this disadvantage is clearly outweighed by the enormous gain in accuracy. In conclusion, two investigated methods are recommend for singular integral approximation: Highly accurate and efficient Duffy quadrature and modified Gauss-Legendre quadrature that is robust and simple to implement.