Corrected trapezoidal rules for singular implicit boundary integrals

We present new higher-order quadratures for a family of boundary integral operators re-derived using the approach introduced in [Kublik, Tanushev, and Tsai - J. Comp. Phys. 247: 279-311, 2013]. In this formulation, a boundary integral over a smooth, closed hypersurface is transformed into an equivalent volume integral defined in a sufficiently thin tubular neighborhood of the surface. The volumetric formulation makes it possible to use the simple trapezoidal rule on uniform Cartesian grids and relieves the need to use parameterization for developing quadrature. Consequently, typical point singularities in a layer potential extend along the surface's normal lines. We propose new higher-order corrections to the trapezoidal rule on the grid nodes around the singularities. This correction is based on local decompositions of the singularity and is dependent on the angle of approach to the singularity relative to the surface's principal curvature directions. The proposed decomposition, combined with the volumetric formulation, leads to a special quadrature error cancellation.


Introduction
Boundary integral methods (BIMs) are employed in a wide range of applications for solving partial differential equations with conditions defined on boundaries of subregions and at infinity.In a BIM, one needs to solve a system of boundary integral equations (BIEs) involving singular integral operators acting on an unknown function defined on the boundaries and at infinity.Typical computational challenges for a boundary integral method involve developing high-order quadrature rules for the singular integrals and efficient dense matrix-vector computations for solving the resulting linear systems.Overcoming these challenges leads to highly efficient and accurate solutions for the associated partial differential equations.
We consider applications that require solving BIEs on a sequence of surfaces that are challenging to parametrize.These may include level set methods [24,26,23] and the closest point method [25,18,19] used to track evolving surfaces on a grid, particularly when the PDE solution is needed only at a small set far away from the surfaces.In such situations, it is not immediately convenient to use any classical BIM.The implicit boundary integral formulations are derived in [15], aiming at these situations.We refer to this as the IBIM approach.In [16,17], further analysis related to the closest point projection is reported.In [4], a similar formulation is derived to approximate the hypersingular integral equations arising from the Neumann problems of the Helmholtz equation.The method evaluates the limit of a family of surface integrals utilizing extrapolative averaging kernels.An IBIM is applied to compute electrostatic potential from large molecules submerged in a solvent in [30].That paper also demonstrates that IBIM, coupled with an "off-the-shelf" Fast Multipole Method, can easily be applied to solve the equations for very large molecules.Partial differential equations arising from calculus of variation problems defined on closed surfaces can be solved with high order convergence rates using similar strategies; see [5,11,21,22].
In the center of these formulations lie volume integrals with identical evaluations to the corresponding surface integrals.The volume integrals involve integration in thin tubular neighborhoods of the surfaces in the ambient space, and do not require surface parameterizations.In principle, the integrals can be approximated on a wide range of meshing.Among the existing work, these volume integrals are discretized on Cartesian grids using the trapezoidal rule and a lower order regularization of the layer singularities.
This paper presents higher-order accurate quadrature rules for the singular integrals arising from the non-parametric boundary integral formulation discussed above.In these formulations, the singularities in the integral operators concentrate along the surface normal lines; these lines generally do not lie on the grid.This feature is atypical in the more classical boundary integral formulations.Our approach is a generalization of the methods in [20]: Regular trapezoidal rule-based summation is performed over the grid nodes lying in the regions, excluding small, grid-dependent neighborhoods around the singularities.This approach is called the punctured trapezoidal rule.We derive additional corrections corresponding to the skipped grid nodes and add them to the punctured trapezoidal rule.The resulting corrected trapezoidal rule is second-order accurate with respect to the uniform grid spacing of the underlying Cartesian grid.We discover an additional benefit of the non-parametric approach -there are cancellation of errors which leads to a improved order of accuracy in practice.
The structure of the paper is as follows: in Section 2 we present an overview of how to express solutions to Laplace and Helmholtz problems using boundary integral equations, and introduce the volumetric extension setting to express surface integrals via volume integrals.In Section 3 we present singularity regularization methods for volume integrals of Section 2. In Section 4 we present a general singularity correction framework for the volume integrals of Section 2, and in Section 5 we go into details about how to apply this methods to the Laplace singular kernels.Finally, Section 6 presents numerical results for the methods of Sections 3 and 5 applied to the evaluation of Laplace potentials.

A short review of boundary integral equations
Boundary integral methods can be used to solve the Laplace and homogeneous Helmholtz equations in both bounded and unbounded domains.Given a bounded domain D ⊂ R 3 , the problems are where Ω = D for the interior problem, and Ω = R 3 \ D complement of the closure of D for the exterior problem, and n is the outward pointing normal to the surface ∂Ω =: Γ.For the exterior Helmholtz problem, in order to ensure uniqueness (see §3, Thm 3.13 in [6]), the solution must also satisfy the radiation (or Sommerfeld) condition The Helmholtz equation with λ = 0 becomes Laplace equation, and the two problems are heavily related.In three dimensions, the fundamental solutions for Laplace and Helmholtz are respectively: exp(iλ x − y ) x − y .

Solutions as layer potentials
A solution u to ∆u + λ 2 u = 0 in Ω = R 3 \ D, can be expressed, for x ∈ Ω, as (Single-Layer potential) u(x) =S[α](x) := ∂Ω G λ (x, y)α(y)dσ y , The functions G λ and ∂G λ ∂ny are called single-layer (SL) and double-layer (DL) kernels respectively.The functions α, β, ζ, are called single-layer density, double-layer density, and combined-layer density.Their expressions are derived by the application of the Green's identities and the properties of the solutions to interior and exterior Helmholtz problems (see §2.4, 2.5, 3.2 and 3.4 in [6]): for y ∈ ∂Ω where v SL , v DL , v CL are solutions to ∆v + λ 2 v = 0 in Ω c = D with boundary conditions: The uniqueness of v CL requires ξ = 0.In practice, the value for ξ is usually tuned to improve the properties of the numerical methods.
When expressing the solution to the problem as a layer potential, the density is unknown.Using the boundary conditions, we can find BIEs to which the solution is the density.In the Dirichlet problem with boundary conditions u = f on ∂Ω the solution can be expressed as either a single-layer potential or a double-layer potential: with minus for the interior and plus for the exterior problem.Note that in general the double-layer formulation is preferable as it involves the solution of an integral equation of the second kind: this leads to non-singular matrices when discretizing the integral operators with Nyström methods.
For the Neumann problem with boundary conditions ∂u ∂n = g on ∂Ω, the singlelayer formulation is preferable as it avoids the appearance of a hypersingular kernel, and the boundary integral equations to solve is with plus for the interior and minus for the exterior problem.
In addition to the single-and double-layer potentials, we also consider the potential appearing in (7), called the double-layer conjugate potential: Its treatment is going to be analogous to the treatment of the double-layer potential.The function ∂G λ ∂nx is called double-layer conjugate (DLC) kernel.It is important to note that these boundary integral equations (5 -7) are valid for both Laplace and Helmholtz equation, and the only thing that changes in the formulae is the parameter λ (the wavenumber), and consequently the kernel G λ .
If the function F is smooth, the freedom of choice of the position of the nodes and of the quadrature rule makes it easy to attain high accuracy.However, if the function F is singular, for example in the origin, then using standard quadrature rules results in a great loss of accuracy.To remedy this loss, several classes of methods have been developed, for different kinds of singular integrands.
An important class of methods to handle these singular integrands, the methods of singularity subtraction, approaches the problem by locally approximating the surface and evaluating the integral analytically, and then adding a correction term dependent on the surface approximation.If the kernels are similar to the ones for which the analytical results exist, singularity subtraction is applied and those same results are used [8,7].
Another class of methods, the methods of singularity cancellation, use a change of variables to put themselves in a setting where it is possible to split the integral in a smooth one and a singular one which is only defined close to the singularity points, and can be computed to high accuracy using exact local parametrization and a suitable quadrature rule, e.g.trapezoidal rule in polar coordinates [3,29].
A third class of methods, the methods of singularity regularization, relies on regularizing the kernel so that rules for smooth integrands can be applied, and then adding corrections to account for the different kernel based on analytical results, or on Richardson extrapolation [2,10,15].
A newer class of methods, called quadrature by expansion (QBX), handles the problem by expanding and treating the kernel (and the corresponding layer potential) away from the surface target point, e.g. with Taylor or spherical harmonics, consequently working with smooth integrands, and then evaluate the results back on the surface.It relies on the smoothness of the expansion terms because the new target point is not on the surface, and on the convergence of the expansion in the surface target point [14].
Finally, the methods of singularity correction aim to develop specialized quadrature rules to deal with families of singular integrands by modifying the weights of an existing quadrature rule, often trapezoidal rule, close to the singularity point [13,20,28].Marin, Tornberg and Runborg [20] developed corrections to the trapezoidal rule for singularities of the kind |x| γ , γ ∈ (−1, 0), in one dimension and x −1 in two dimensions, proved convergence order, and found an analytic expression for the weights in one dimension.Wu and Martinsson expanded these results to log |x| in one dimension [28] and found analytic expression for weights in one and two dimen-sions [27].
The majority of above-mentioned methods require explicit knowledge of the parametrization, and the possibility to choose the position of the nodes around the singularity; moreover, often the singularity point lies in one of the nodes.
In our setting however, the position of the nodes on the surface is going to be determined by the projections of the nodes in the volume onto the surface, which cannot be assumed to have any particular structure (see Figure 8).Moreover, because of how the integrand is extended from the surface to the volume, instead of a single singularity point on the surface, we will have the singularity lying along a straight line in three dimensions.
We will approach this problem then by splitting the three-dimensional trapezoidal rule into the weighted sum of all the two-dimensional trapezoidal rules on each two-dimensional grid and correct each one separately.

Volumetric extensions of the layer integrals
Let Ω ⊂ R n be a bounded open set with C 2 boundaries, and ∂Ω =: Γ.We shall refer to Γ as the surface.Let f be a function defined on Γ (or R n ).In this Section we present an approach for extending a boundary integral to a volumetric integral around the surface.Instead of parameterizations, this approach relies on the Euclidean distance to the surface, and its derivatives.More precisely, we define the signed distance function and the closest point projection If there is more than one global minimum, we pick one randomly from the set.Let C Γ denote the set of points in R n which are equidistant to at least two distinct points on Γ.The reach τ is defined as inf x∈Γ,y∈C Γ x − y .Clearly, τ is restricted by the local geometry (the curvatures) and the global structure of Γ (the Euclidean and geodesic distances between any two points on Γ).
In this paper, we assume that Γ is C 2 and has a non-zero reach.Let T ε denote the set of points of distance at most ε from Γ: Then, for ε < τ , P Γ is a diffeomorphism between the level sets of d Γ and We define the extension (or restriction) of the integrand f by As in [15,16], we can then rewrite the surface integral (9), for any η ∈ [−ε, +ε], as where J η (x ) is the Jacobian of the transformation from Γ to the level set Γ η := {x ∈ R n : d Γ (x) = η}.In R 3 , the Jacobian J η (x ) is a quadratic polynomial in η: where H(x ) and G(x ) are respectively the mean and Gaussian curvatures of Γ η at x , and σ 1 σ 2 (P Γ x ) is the product of the first two singular values of the Jacobian matrix of P Γ evaluated at x .See [16] for more detail.
To extend (9) to a volumetric integral, we now average the integral on the right hand side in ( 14) over η ranging from −ε to ε, using Applying the coarea formula, we have Thus from the surface integral, we derive a volume integral with the same evaluation: where Our primary focus is when f (x) is replaced by a function K(x, y)ζ(y), with K : R n × R n → R and K(x, y) singular for x = y (corresponding to the layer potentials reviewed in the previous Section): When a function g : Γ → R is given, we may form an integral equation for the unknown density ζ.For example, in the case of the double-layer potential (6), the equation is Suppose that for any x, we are interested in evaluating K(x, y) at the point on Γ that is closest to y.This can be done by Hence we refer to K(x, y) as the restriction of K.If K(x, y) is singular for x = y, then K(P Γ x, y) is singular on the set i.e. for a fixed x * ∈ Γ, K(x * , y) is singular along the normal line passing through x * , while K(x * , y) is singular in a point.In Figure 1 the singular behavior of K(x * , y) along the normal is illustrated.
In conclusion, instead of approximating (16), we approximate for functions ρ that are integrable in T ε .
Corresponding to (17), we have the equivalent implicit boundary integral equation The solution ρ will coincide with the constant extension along the normals of ζ.To see this, we write the two equations: The first equation corresponds to (17) where the integral has been rewritten and the target point x ∈ T ε is projected onto Γ.The second is (20) where the equation is imposed for ρ function defined in T ε .If we take the difference of these two equations, we find The kernel of the operator on the left-hand side coincides with the kernel of the original operator, so whenever the solution is unique, ρ(x) = ζ(P Γ x) for any x ∈ T ε .
In this paper, we will concentrate on developing numerical quadratures for the extended singular integral operator J[ρ](x) for x ∈ Γ (equivalently, J[ρ](P Γ x) for x ∈ T ε ).The quadrature rules will be constructed based on the trapezoidal rule for the grid nodes T h ε := T ε ∩ hZ 3 , which corresponds to the portion of the uniform Cartesian grid hZ 3 within T ε .Since the integrand in (19) is singular for x ∈ Γ, the trapezoidal rule should be corrected near x for faster convergence.Correction will be defined by summing the judiciously derived weights over a set of grid nodes denoted by N h (x).The sum will be denoted by R h (x).Ultimately, the quadrature for J[ρ](P Γ x) will involve the regular Riemann sum of the integrand in T h ε \ N h (x), and the correction R h (x) in N h (x): The contribution of this paper is a high order, trapezoidal rule-based, quadrature rule for J via J . Figure 8 demonstrates a typical configurations the points y m in the summation for a torus.
In the following two Sections, we will see how to build the correction term R h (x) using two different approaches: a function regularization independent of trapezoidal rule (Section 3), and the corrected trapezoidal rule (Section 4).Each approach will determine the set N h (x) differently.

Correction via regularization of singularity
In this Section we present an approach that locally regularizes a singular kernel before the discretization.In this approach, a special Lipschitz continuous function, Ψ, is used to replace the kernel in a neighborhood around the kernel's singularity.The integral with the regularized kernel is then extended following (19).Again, the resulting implicit boundary integral can be discretized on different meshings.When the trapezoidal rule is applied to J involving the locally regularized K, we find an expression of the kind (21) where the term R h (x) involves a local sum of the integrand with K replaced by Ψ.

A localized regularization approach
We consider regularization of K(x, y) constructed in the following fashion: where Ψ Γ,r 0 thus is a function which substitutes K close to the singularity point.We choose this as a simple function (constant, or linear in x − y ) which approximates K(x, • ) weakly for C 1 functions on the r 0 neighborhood of x, such that Here, B r 0 (x) is the ball with radius r 0 , centered at x. Applying the trapezoidal rule to the integral (19) with the regularized kernel ( 22), we get a correction to the trapezoidal rule of the form (21): = h 3 where: K, K reg r 0 , and Ψ Γ,r 0 are the restrictions (18) of K, K reg r 0 , and Ψ Γ,r 0 respectively.The formula (21) holds with: In order to determine Ψ Γ,r 0 we want it to satisfy (23) for ρ ≡ 1 with an error at most O(r 2 0 ).However, given the lack of an explicit parametrization of the surface, we approximate Γ in the integrals in ( 23) by a suitable paraboloid, Γx , defined from the principal curvatures of Γ at x (as shown in [15]).The domain Γ ∩ B r 0 (x) is furthermore replaced by a neighborhood M(x, r 0 ) ≈ Γx ∩ B r 0 (x).Eventually, we seek Ψ Γ,ro satisfying If v is a Lipschitz continuous function on Γ, we can write for the kernels we are interested in.This approach is used in [15] and [4].
The rest of this Section will be now dedicated to showing results of this approach for the Laplace and Helmholtz double-layer kernels.

Application to the Laplace and Helmholtz double-layer kernels
Consider the double-layer kernel K(x, y) = ∂G 0 ∂n y (x, y) for Laplace.In [15] the The constant C Γ,r 0 represents the average of the integrand on the set, and Ψ Γ,r 0 (x, y) = C Γ,r 0 regularizes the double-layer kernel: The expression for C Γ,r 0 found in this setting, dependent on the principal curvatures and r 0 , is: The calculations and details about the setting together with the exact definition of M(x, r 0 ) can be found in the A.2.When treating the Helmholtz double-layer kernel, we can apply this constant regularization approach to the additional term which differentiates it from the Laplace double-layer kernel: the gradient of the Helmholtz fundamental solution in three dimensions is hence, the double-layer kernel for Helmholtz takes the form: In the expression above the factor exp(iλ x − y ) is a Lipschitz continuous function in y, and we know how to deal numerically with the Laplace double-layer kernel ∂G 0 ∂ny (using regularizations (26) or (30), or the corrected trapezoidal rule which will be the focus of Section 4), so we focus only on the secondary kernel which, if κ 1 = κ 2 , is undetermined in x = y, as the value depends on the direction of approach.The maximum and minimum limit values are the ones found traveling along the principal directions, equal to κ 1 2 and κ 2 2 respectively.
Then we wish to find a constant function Ψ Γ,r 0 (x, y) ≡ CΓ,r 0 such that the integral around the singularity point is approximated well, This requirement gives It is interesting to notice that the first term in CΓ,r 0 is the average of the maximum and minimum limits of the integrand.The function (x−y)•ny x−y 2 , for x, y ∈ Γ, can then be regularized using (28):

New regularization with linear function (cappuccio)
Here, we consider regularizing with a class of function that are linear with respect to the distance to the singularity.We construct the hood-like (cappuccio in Italian) function This condition imposes one constraint; the second constraint we impose is that which means that Ψ L Γ,r 0 takes as outermost value the average of the kernel on the boundary of M(x, r 0 ).
The motivation is that the discontinuity in the regularized kernels can be significantly smaller, for small r 0 , than the ones regularized by constants.Consequently, the quadrature errors can be smaller.See Figures 2 and 3 for a comparison.In particular, the third columns in Figures 2 and 3 show that the existing and proposed regularizations will lead to discontinuity between the regularized and original functions if the direction of approach to the target point is not taken into account.Consequently we believe that a possible future work is to develop a continuous regularization by including dependence on the principal directions τ 1 , τ 2 in addition to the curvatures of the surface, e.g.
We get the following expression for the regularized kernel,  where Note that Ψ L Γ,r 0 scales as 1/r 0 for small r 0 .In Section 6, we shall present some numerical convergence studies of the approaches mentioned in this Section.

The corrected trapezoidal rules
We have shown in Section 2.3 that the singular integrals of interest can be characterized by their singular behavior along lines in R 3 .We view the trapezoidal rule on a three-dimensional uniform Cartesian grid as the sum over the trapezoidal rules applied to the two-dimensional uniform grids.On each two-dimensional grid, the case is reduced to correction of trapezoidal rule for functions that are singular only at a single point.However this point is typically not lying on any grid nodes.
We will first present the trapezoidal rule and the existing methods for correcting it to achieve higher order convergence rates.We will then present our generalization of these works.

The punctured trapezoidal rules
Let f : R n → R be a compactly supported smooth function.We are interested in approximating its integral R n f (x)dx by utilizing values of f on the uniform grid hZ n .By the compact support, the trapezoidal rule applied to f becomes the following simple Riemann sum: When f is compactly supported the order of accuracy of such summations depends on the regularity of f : if f ∈ C p , the error is O(h p ) (see §25.4.3 in [1] and §5.1 in [12]); the trapezoidal rule enjoys spectral accuracy if f ∈ C ∞ .When f is continuous in R n \ {x 0 } and singular at x 0 , where R n f (x)dx exists as a Cauchy principal value, it is natural to modify the trapezoidal rule by skipping the summation over the grid nodes within certain distance to x 0 : where N h (x 0 ) determines which grid nodes we remove.We will call (32) the punctured trapezoidal rule when N h includes only a single grid node; in other words, The punctured trapezoidal rule converges, but with lower order rates at best, even though f may be C ∞ in the punctured domains.For example, in one dimension for f (x) = log |x|, the order of convergence is sublinear O(h p ), p < 1, and in two dimensions for f (x) = x −1 the order is 1.The large decrease in order is exactly the property we would like to address with the correction technique.
The idea is to add a correction term to (32), which makes up for the integral over N h .In the following, we describe an approach for defining such corrections in detail.

Corrections for the trapezoidal rule
From this point forward, we will assume the function f can be factored into the following form where s represents an integrable function, singular in the origin, and v represents a smooth compactly supported function in R n .In this Section we discuss a general approach to developing high order quadratures for the integration of such type of functions.In Section 5 we will provide specific choices of s for use with single-and double-layer kernels arising from the Laplace or Helmholtz operator.The trapezoidal rule is a sum of the function values on the grid, where all values have the same weight h n .Improving the order of accuracy of the trapezoidal rule by modifying the weights close to the singularity point has been an approach studied and applied successfully with different kinds of singular behaviors and in different dimensions.See for example [13,20,28].
The following is a brief presentation of the one-and two-dimensional corrections found in [20], where x 0 is always assumed to be the origin.
The starting point is the punctured trapezoidal rule in one dimension.When s(x) = |x| γ for −1 < γ < 0, an error expansion of the following type can be derived, The goal is to find the constant ω, which is independent of v (but depends on γ), and use it to correct the rule as While T 0 h is of order 1 + γ, the method Q 0 h is of order 3 + γ.Note that Q 0 h only modifies the original trapezoidal rule in one point; the value in the singular point is replaced by the value of the smooth part v(0), weighted by ω and a suitable power of h.In general ω is a functional of the singular function s and we write ω = ω[s].
In order to find ω[s] we define ω(h) as the actual error of T 0 h for a fixed h and a smooth test function g with g(0) = 1, scaled by h 1+γ .More precisely, From the error expansion above, since g(0) = 1, we see that Hence, the weight ω(h) converges to ω[s] = 0 for h → 0 + , independent of g.This is a crucial property, which makes it possible to compute, store and reuse ω for the integration of any integrand of the kind (34).In order to do that one needs to be able to accurately compute the integral s(x)g(x)dx containing the test function.As this computation is only needed for one function g, it can be done either by analytical means or adaptive high order numerical integration.If the test function g is chosen more flat at the singularity point, such that 0 = g (0) = g (0) = . . .one can show that the convergence will be faster, which makes the numerical computations easier.In fact, g would be ideally the constant function g(x) = 1 but in order to avoid dealing with the boundary conditions of trapezoidal rule and keep the expression of T 0 h equal to the Riemann sum with the exclusion of a single node, g is taken compactly supported.
Higher order corrections are also possible, where more terms in the error expansion are cancelled.The weights must then be modified in more points close to the singularity.The condition (35) can be interpreted as requiring that T 0 h , corrected with the weight ω(h), integrate s(x)g(x) exactly.When multiple weights are used, the weights are similarly defined by requiring that the modified method integrates not only s(x)g(x) exactly but also s(x)g(x)x, s(x)g(x)x 2 , . . . .A set of h-dependent weights are then obtained, which converge as h → 0 + .One can also apply the same idea to other singularities.In [28] this was done for s(x) = − log |x|.Then the factor h 1+γ must be replaced by an expression a(h) = h(2 − log h) and a second order method is obtained In two dimensions, similar to one dimension, for functions (34) with s(x) = x −1 , the corrected trapezoidal rule is defined as where ω[s] is calculated as the limit: for a test function g with g(0) = 1.In [20] it was proven that the corrected method for s(x) = x −1 is third order accurate,

Corrections for singularity unaligned to the grid
To prepare for the proposed quadrature rules for implicit boundary integrals, we first generalize the approach presented in Section 4.2 to the case when the singularity does not lie on a grid node.We consider two dimensions and retain the assumption that f can be factorized as s(x − x 0 )v(x) where s has a singularity and v is smooth and compactly supported.However, the singularity is now in a point x 0 which may not be part of the grid.We let x ∆ be the grid node closest to x 0 , or one of the closest in the case that it may not be unique, such that as shown in Figure 4.When (α, β) = 0 the usual trapezoidal rule is well-defined also for the singular function and the same type of error expansion holds as for the punctured trapezoidal rule in the previous Section.However, the error constant is not uniform and blows up as (α, β) → 0. We therefore use the punctured trapezoidal rule also for unaligned grids as the base method for correction.The singular functions considered in this paper are of the form f (x) = s(x − x 0 )v(x) where |s(x)| ∼ ||x|| −1 .For those functions the same scaling in h as s(x) = ||x|| −1 is appropriate and we define the single-correction trapezoidal rule for unaligned grids in two dimensions as The weight is given as the limit of the sequence: where as before ω δ [s; α, β] is defined using a smooth compactly supported test function g with g(0) = 1, In the last step we shifted the exact integral by x 0 and the trapezoidal rule by x ∆ , the closest node of the grid δZ 2 , to show that the weight, in addition to s, only depends on the difference x ∆ − x 0 , i.e. on α and β, not on x 0 itself.The expression converges quickly when the stepsize δ is halved, and the accurate computation of ω is possible without needing specialized quadratures for singular integrands.In this paper, ω[s; α, β] will be computed offline and tabulated for a suitable set of (α, β), and for relevant functions s; for values outside of the tabulation, we will use interpolation.In the next Section, we will discuss a few specific cases involving layer kernels, and we shall then present more details about the approximation of ω via tabulation and interpolation.
When α = β = 0 we get the same weights as in the aligned case.In particular, for s(x) = x −1 the limit (39) will find the same weight as the one found in [20].For the more general kernels considered in this paper and with unaligned grids, in numerical experiments we observe an error expansion of the type where F 1 is a smooth function of α, β.Moreover, F 1 (0, 0) = 0 = F 2 (0, 0), which means that we have the same third order error (37) as for s(x) = x −1 when the grid is aligned.The properties of F 1 will be further explored in Section 6.2.Proving these error results rigorously is in program for future research.

Corrected trapezoidal rules for implicit boundary integrals
We describe our approach in developing corrected trapezoidal rules for the family of integrals defined in (19): Without loss of generality, we consider, as the target point, x * = (x * , y * , z * ) ∈ Γ, where the surface normal at x * is n = (n 1 , n 2 , 1): from this point forward we will only consider this case, and if the normal direction points instead more towards the x or y directions, we can apply a change of coordinates and proceed with the same reasoning.To simplify notation we let f be the integrand Furthermore, since the restricted kernel K(x * , y) = K(x * + t n, y) for all t ∈ R, the integrand f is singular along this line.
The plan is then to construct a quadrature for (42) "plane-by-plane" on the grid hZ 3 .First the standard trapezoidal rule is used in the z-direction, Then, the corrected trapezoidal rule is used to compute the integrals on each plane, See Figure 5 for an illustration.As the z-component of n is 1, f is singular in one point only when restricted to the planes.We can therefore use the quadrature rules described above.
We let ȳ = (x, y) denote a point in the xy-plane and introduce the projection onto this plane ȳ = πy = π(x, y, z) = (x, y) .
For a fixed z, the singular point ȳ0 of f ( • , • , z) is then given by The corresponding closest grid node is denoted by ȳ∆ (z) and its shift parameters α = α(z), β = β(z).The factorization of f ( • , • , z) will be of the kind where s is smooth in the second argument, which ensures that the partial integral f dxdy is smooth in z, justifying the use of the standard trapezoidal rule in this variable.The functions s and v correspond to factorizations specific to the kernel K and the geometry of Γ.They will be discussed in detail in the next Section.
Figure 5: Intersections of the line in three dimensions Intersection of the line passing through x * with direction n with the planes {z = z k }: on every plane, the intersection will be y 0 (z k ) (orange circle), and the closest grid node will be y ∆ (z k ) (red square).The parameters which characterize the position of y 0 (z k ) with respect to the grid hZ 2 are (α With this notation we can now give the precise form of the corrected method After the discretization, z k = kh, and we can write the full method as where Then we can see that this method is of the form (21) with N h (x) as in (48) and 5 Factorization of the Laplace kernels and the resulting quadratures In the previous Section 4. 3 we have seen the corrected trapezoidal rule (47) for the family of integrals of the kind (42).We have however not gone into detail about the form of the singular functions and the consequent splitting (46) involved, as they depend on the functions (43) and the corresponding kernels where K is one of the Laplace layer kernels: In this section, we will derive the proposed quadrature rules for the above kernels.In the same way the surface integral of the single-layer potential (2) in the non-parametric setting becomes the volume integral (19), the double-layer (3) and double-layer conjugate (8) potentials are extended to volume integrals in the tubular neighborhood T ε ⊂ R 3 .Between the DL and DLC kernels we will only consider the DLC kernel, because the singularity behavior is identical.
In Section 5.1 we present our approach to define the factorization f = s v in (46).s describes K close to the singularity point and it will be written as the product of S , where S takes a very simple form that is easy to work with.After introducing the function S, we will explain how to find the corresponding function in Section 5.1.1,and then show the full expression of (47) in Section 5.1.2.Finally, in Section 5.2, we will describe how the weights ω[s; α, β] are defined from S and computed.A brief outline of the process can be found in Table 1.

Correction formula for the three kernels
Let x * ∈ Γ be a target point with normal to the surface n = (n 1 , n 2 , n 3 ), n 3 = 0, n = 1.We want to apply the three-dimensional second order correction formula (47) to the layer potential R 3 K(x * , y)ρ(y)δ Γ,ε (y)dy with one of the three layer kernels (50), for example the single-layer kernel K(x, y) = G 0 (x, y) = (4π x−y ) −1 .The starting point is the following singular function: which represents the reciprocal of the distance from a point r to the line with direction n passing through the origin: {t n : t ∈ R}.This choice is equivalent to approximating the distance on the denominator of the three kernels (50) as x * − P Γ y ≈ x * − y T M (x * ) , where y T M (x * ) is the projection of y onto the tangent plane to Γ at x * .The reciprocal of the distance x * −y T M (x * ) is conveniently given by S(x * − y, n).
With S, we write formula (46) as Table 1: Key ingredients in Section 5.
We write ṽ as This defines the new function v which is smooth; hence f = S • • v, and the weights for the corrected trapezoidal rule will be therefore derived for the singular function

Derivation of the formula for the new factor
Due to the closest point projection in K, the limit function depends on the intrinsic geometry of Γ (the principal curvatures) as well as the orientation and distance of y 0 (z) to Γ (the signed distance η = η(z)).For different values of z but fixed direction q, P Γ may map the lines (tq, 0) + y 0 (z) to different curves on Γ.The speed at which these curves, P Γ ((tq, 0) + y 0 (z)), pass through the target point x * may vary, depending on the curvature of the curve.
To calculate explicitly the limit that defines , we will replace the projection P Γ by P Γ, a high order local approximation of the closest point projection to the osculating paraboloid at x * .In the following, we will present the derivation of an explicit formula for , based on and building up from the simplest case.In the derivations, we let τ 1 , τ 2 , n be the orthonormal basis of R 3 composed of the principal directions τ i , with corresponding principal curvature κ i , ordered such that τ 1 × τ 2 = n.We write a point in this basis as (a, b, c) T := aτ 1 + bτ 2 + c n .
We first assume the z planes are parallel to the tangent plane T M (x * ), i.e. n = e 3 .Thus we want to find the limit For convenience, we translate the problem so that x * is in the origin, and work only in the τ 1 , τ 2 , n basis.See a representation of this setting in the left plot of Figure 7.The paraboloid Γ will then be In a sufficiently small neighborhood of the origin, a point (x, y, η) T and its closest point i.e. the vector pointing to the closest point on Γ should be normal to the surface, with magnitude equal to the distance to the surface.Along (tp 1 , tp 2 , η) T , we have: For example if p = (1, 0) the limit is taken along the τ 1 direction, the projected point will travel along the curve corresponding to the first principal direction, and the limit value for the double-layer conjugate kernel will be: In general, we have We now consider the more general case in which n = e 3 .Again we consider the target point to be in the origin: x * = 0. We define the plane Π z := {(x, y, η(z)) T : x, y ∈ R} parallel to the tangent plane T M (x * ) at distance η(z) := d Γ (y 0 (z)).Fixed z ∈ R, the projection takes a point x to the intersection of the line {x + tn : t ∈ R} and the plane Π z .Let q = (q 1 , q 2 ) ∈ S 1 , and let (tq 1 , tq 2 , 0) + y 0 (z) be a point on the z plane.To find the limit, we first consider the projection of the line (tq 1 , tq 2 , 0) + y 0 (z) onto the plane Π z , and then apply the previous formula; in the expression of (58) we consequently have (t p 1 (q 1 , q 2 ), t p 2 (q 1 , q 2 ), η(z)) T = P Πz ((tq 1 , tq 2 , 0) + y 0 (z)) instead of (tq 1 , tq 2 , 0) + y 0 (z).This change does not affect the limit expression because of the property for small values of t, which expresses how the orientation of the plane z with respect to the basis τ 1 , τ 2 , n does not affect significantly the projection of points close to the singularity point y 0 (z).
The parameters a, b, c, θ 0 relate τ 1 , τ 2 to the standard R 3 basis {e i } 3 i=1 .They are given by (a, b, c) = (sin θ cos ξ, sin θ sin ξ, cos θ) , where θ is such that cos θ = n 3 , i.e. it is the second spherical coordinate of n; ξ is such that tan ξ = e T 3 τ 2 /e T 3 τ 1 , hence it is the angle between τ 1 and the projection The surface is approximated around the target point with a paraboloid defined by the surface's principal curvatures and directions.Points on the circles on each plane are mapped to the closest points on the paraboloid Γ (instead of Γ) for calculation of the limit defined in (55).Right plot: n = z; a circle (blue) drawn on the z plane (yellow plane) around the singular point y 0 (z) becomes an ellipse (red) when projected on the plane Π z (blue plane) parallel to the tangent plane of Γ in x * .The angle ψ between e 1 and a given direction (red direction) on the plane z will correspond to the angle ψ between τ 1 and the projected direction (yellow direction) on the plane Π z .
of e 3 on the plane τ 1 , τ 2 ; and θ 0 is such that tan θ 0 = e T 1 τ 2 /e T 1 τ 1 , meaning it is the angle between τ 1 and the projection of e 1 on the plane τ 1 , τ 2 .
Given a unit vector (q 1 , q 2 ) on the z plane, formulae (59-60) map it to the unit vector p 1 τ 1 + p 2 τ 2 .The unit direction (p 1 , p 2 ) on the plane Π z , (p 1 , p 2 , 0) T , is the corresponding direction in which the limit (in the definition (55) of (q; z)) will be evaluated.Recall that the formula for the limit for any given direction is already derived in (58).For convenience, we write the unit vector (q 1 , q 2 ) = (cos ψ, sin ψ), and correspondingly we will treat p 1 , p 2 , defined in (59-61), as functions of ψ.
Thus we write the formulae for for the double-layer conjugate, double-layer, and single-layer kernels as: Again, η(z) = d Γ (y 0 (z)) is the signed distance of y 0 (z) to the surface.
The center and right column in Figure 6 illustrate how the formulae found approximate the behaviour f /S.The center column plots the function in (62-63) found for the three Laplace layer kernels, while the right column shows the difference between the function found and the values f /S.

The quadrature formulae
With defined above, we will work with the following singular function-smooth function factorization: with where the weight for the corrected trapezoidal rule applied to f is computed with s = S • .The function s completely captures the asymptotic behavior of the given kernel K in ȳ0 (z) and we can apply (38) to f (ȳ, z): As long as is non-zero, the function v is well-defined through (64), away from ȳ0 (z) and by continuity at ȳ0 (z).As can be seen from ( 62) and (63), this is always the case for the single-layer kernel but for the double-layer kernels in general only if κ 1 and κ 2 have the same sign.Even though is zero only at isolated points, we cannot apply formula (65) as it is numerically problematic.We need a different approach which works as follows.
We use the discontinuity subtraction from K: first, to shorten the formulae below, we define the smooth function We then replace the splitting in (64) by Then v is well-defined everywhere, bounded and continuous around ȳ0 , by construc-tion of via the limit (53).We therefore rewrite (65) as where inside the braces only the second term remains if (α(z), β(z)) = (0, 0).Our corrected trapezoidal rule (66) for the implicit boundary integral takes the form: where The general correction form (21) of this method is then valid for N h (x) as in (48) and

Approximation and tabulation of the weights
We approximate the singular functions using a Fourier interpolation, then tabulate the weights for the simpler singular terms of the expansion, and compose the general weights for any behavior needed.
Given a π-periodic function ( • ; z), such as (62) or (63), we wish to compute the weight ω[s( • , z); α, β], where s comes from the factorization (64): Both the factors in this expression can be seen as functions of the angle of approach, ψ, to the singular point 0: Given the dependence on z is only present in through η(z), we write ω[s; α, β] We use Fourier interpolation to approximate this function with a trigonometric polynomial: Since the weight ω[s; α, β] is a linear functional of the singular function, Therefore, we can precompute the weights for the basis functions and for certain values of α and β: For values outside of the precomputed tables, we interpolate.
In the evaluation of ω cos(2jψ(ȳ)) ȳ ; α, β using formulae (40) and (39), we need the value of the integral: where δ 0j is the Kronecker delta; R is set to 1.9 as the integrand is essentially zero at double precision.The integral R 0 exp(−r 8 )dr is computed once with high precision using common integration libraries and reused for all instances as a constant.
It is impractical to tabulate precomputations of the coefficients {c j } N j=0 , {d j } N j=1 as they depend on too many variables.Instead, we compute c j and d j on the fly, by solving the square linear system Because of the smoothness of the π-periodic functions we deal with, we can use a relatively small N in order to accurately approximate the weights.
In the next Section, we present numerical convergence studies using weights computed with N = 22.The small linear system can be inverted efficiently, e.g. using FFT, with a negligible computational time.
In the convergence studies, an array of weights of dimensions (45,101,101) has been precomputed, with 101 values for α and β in − 1 2 , 1 2 and 45 for the Fourier series with N = 22.Biquintic interpolation is used to approximate the weights for given (α, β) outside of the precomputed values.
In the case ≡ 1, {c j } N j=0 and {d j } N j=1 depend only on θ and φ.In that case, they can also be precomputed, stored, and used in an interpolation process when needed.

Numerical Examples
We demonstrate the convergence and accuracy of the proposed quadrature rules by evaluating the double layer potential with constant density on the surface Γ ⊂ R 3 .We demonstrate the numerical errors computed by the proposed corrected trapezoidal rule for approximating The value of I is known explicitly to be −1/2 for any x * ∈ Γ.So, we report for several randomly chosen x * ∈ Γ.We will compare the results for the four different quadrature rules, including the two new quadrature rules Q L IBIM defined by the regularization (30) and Q h defined in (67).
The integral is first extended to the tubular neighborhood of T ε , as in ( 16)-( 18), using the compactly supported C ∞ averaging function (73) here a ≈ 7.51393 normalizes the integral R φ(x)dx to 1.The surfaces chosen for the tests are a sphere and a torus, centered in a random point in 3D and rotated with random angles along the x-, y-and z-axes.
The sphere is characterized by center C = (0.5475547095598521, 0.6864792402110276, 0.3502726366462485) • 10 -1 and radius R = 0.7.The torus is described by the following parametrization where R 1 = 0.7, R 2 = 0.2, C is the same as the sphere, and Of course to test our algorithms, we retain no information about the parameterizations.The test sphere and torus are represented only by d Γ and P Γ on the given grid.Figure 8 shows the torus that we use and the points used in the quadrature rule in a configuration.The Jacobian J Γ is approximated using a fourth-order centered differencing of P Γ on the grid, see [16].72) on 50 random target points on a sphere.Errors for the punctured trapezoidal rule (32) (black crosses) and the three considered methods in the evaluation of the double-layer potential: Q C constant regularization (26) (blue upward triangles); Q L cappuccio regularization (30) (red downward triangles); Q h corrected trapezoidal rule (67) (magenta circles).The three plots reflect three different settings for the tubular neighborhood width: left plot ε = 0.1; center plot ε ∼ h 0.7 ; right plot ε ∼ h 0.8 .

Convergence studies
We compare the numerical orders of convergence for the quadrature rules discussed in this paper.The quadratures are defined on the grid nodes in T ε .The parameter ε, which describes the width of the tubular neighborhood, comes into play through the function δ Γ,ε , and the factor h/ε determines the number of grid nodes in a cross section of the tubular neighborhood T h ε .It consequently determines how well the integrand is resolved.The errors for the proposed quadratures are formally O h ε p , p ≥ 2, as h → 0. Thus, for fixed ε = O(1), we see the order of convergence resembling p.If we choose ε ∼ h 1/q , we will formally have the errors scale as O(h p(1−1/q) ).If ε ∼ h, then the method will not converge formally, but in the range of the grid resolution considered in practice, the method may yield results with acceptable accuracy.
In Figures 9, 10, and 11 the errors are shown as function of the Cartesian grid's spacing h; the first figure shows the errors for the sphere, while the others show them for the tilted torus.We show also how the errors scale for ε ∼ h 0 (left plot in Figure 9, and Figure 10), and ε ∼ h α for different αs (other plots in Figure 9, and Figure 11).
In Figure 10 we show the error curves for three target points; affected by their relative positions to the grid and the surface, the errors at some target point is larger than at others.Furthermore, the errors at each target point oscillate as one varies h.
We now show that the accuracy of the weights used in the tests is sufficient for the discretization used.In the previous tests, the number of terms in the Fourier expansion was 2N + 1 with N = 22, and each term was tabulated in α, β with N α,β = 101 values each.We repeated the same test as in Figure 10  The table suggests that for this range of parameters, the error from the correction of trapezoidal rule is dominating.
For ε ∼ O(1), the formal order of convergence for the proposed quadrature is O(h 2 ), but we have observed a rate of O(h 2.5 ).In the next subsection, we present a property of our quadrature that we believe leads to the increase in accuracy.

Order increase from error cancellation
As discussed in Section 4.3, our corrected trapezoidal rule is applied on every plane in T ε (see Figure 5), and the error for the whole integral is a sum of the quadrature errors on each relevant plane.Tilted torus with 50 random target points; tubular neighborhood width ε dependent on h.Mean error for the three considered methods in the evaluation of the double-layer potential.Top plots: Q C IBIM constant regularization (blue upward triangles); Q L IBIM linear regularization, cappuccio (red downward triangles); Q h corrected trapezoidal rule (magenta circles).From left to right, the four plots represent: ε ∼ h 0.5 , ε ∼ h 0.75 , ε ∼ h 0.9 , ε = 5h.
For a fixed target point the singular line intersects each plane at a different location relative to the grid.The relative positions are given by the shift parameters α and β, which depend on both the plane's z-coordinate and the grid spacing h.More precisely, let ȳ0 (z) = (x 0 (z), y 0 (z)).Then where {x} denotes the fractional part of x.Defining the 1-periodic function r(x) = {x + 1/2} − 1/2 and using the expression for the singular line ȳ0 (z) in (45) where n = (n 1 , n 2 , 1), we can write This shows that the relative location (α, β) may vary rapidly between planes both in h and z when h is small, and since r is discontinuous, the variation is non-smooth.Let E(z, h) be the quadrature error for one plane.Based on (41) we can express it as where F 1 is a smooth function of (α, β, z).In Figure 13 we can see the function F 1 (α, β, z) for a specific value of z.The total error for the three-dimensional integral, is then where, noting that z The coefficient D(h) is thus the mean of the error coefficients on the different planes.Since F 1 is smooth, and evaluated in a compact set, D is therefore bounded in h.However, α and β are non-smooth in h and underresolved in the second argument in the sum.Therefore, D is not a smooth function of h.This accounts for the irregular convergence plots.See for instance the left subplot in Figure 12 or the right one in Figure 10.
The analysis above would predict second order accuracy for the method, when ε is independent of h.However, in practical computations we typically observe the higher order convergence rate 2.5.We believe this can be explained by a further property of F 1 .Looking at Figure 13, we may notice a skew-symmetry in F 1 (α, β, z) for fixed z.It is reasonable to expect that the average value of F 1 over α and β is much smaller in module compared to the maximum error.In fact, we conjecture that the average value of F 1 for fixed z is zero, In the sum of F 1 (α(h, k), β(h, k), z k ), defining D(h) the first two arguments (α, β) vary much faster than the third (z k ).If the sequence k → (α(h, k), β(h, k)) has some ergodic property the sum will behave similar to the full integral, which would be zero A precise analysis of this effect is beyond the scope of this article.Here we just show in Figure 12 an example of how F 1 (α(h, k), β(h, k), z k ) and (α(h, k), β(h, k)) may vary for two different h that are very close to each other.

A.1 Relating curvatures and the principal directions on parallel surfaces
Let Ω ⊂ R 3 be a bounded domain, and let d Γ and P Γ be the signed distance function and the closest point projection defined in ( 10) and ( 11) in Section 2.3.We assume that the distance function is twice continuously differentiable in the tubular neighborhood The derivation of the proposed quadratures relies heavily on the knowledge of geometrical information of the surface Γ, through that of the level sets of d Γ .In this Section, we relate the principal curvatures and the corresponding directions on different parallel surfaces Let z be an arbitrary point in T ε , and η = d Γ (z).The curvature information of Γ η at z can be retrieved from the Hessian of d Γ .Through eigenvalue decomposition, we have where κ1 and κ2 are the principal curvatures of Γ η at z and τ 1 , τ 2 the corresponding principal directions.One can derive easily that the following formula, relating the principle curvatures κ i of Γ at P Γ z and κi of Γ η at z: See for example [9] ( §14.6 Appendix: Boundary Curvatures and Distance Function).
The principal directions will remain the same: Lemma 1.Let Γ be a C 2 surface, z ∈ Γ; let Γ η be a parallel surface, and z ∈ Γ η such that z = P Γ z.The principal directions at z coincide with the principal directions at z.
Proof.The tangent plane T M η (z) for Γ η at z is parallel to the tangent plane T M 0 (z) for Γ at z. Let (b 1 , b 2 ) be an orthonormal basis for the plane T M 0 (z), and v = b 1 cos θ + b 2 sin θ a unit vector.We can consider the plane H v passing though z and parallel to the normal vector, and H v ∩ Γ will locally be the support of the regular curve γ θ (s), which is the normal section of Γ at z. Corresponding to the normal section we can calculate the normal curvature κ(θ) of Γ at z along v.
The maximum and minimum values of this function are going to be again θ 1 and θ 2 , as then the values θ i , i = 1, 2 are extrema also for this case, and the second derivatives have the same sign as the ones on Γ. Consequently the angles at which maximum and minimum are attained are the same, and the principal directions on the two parallel surfaces coincide.

A.2 Calculation of the regularizations of the double-layer kernels
Given a target point x ∈ Γ, and r 0 > 0, let M(x, r 0 ) ⊂ Γ be a neighborhood of x dependent on the parameter r 0 ; we will define it more clearly later.We want to find Ψ Γ,r 0 such that Ψ Γ,r 0 (x, y)dσ y .
The paraboloid is the surface Γ defined for points close to the origin with coordinates (x, y, z(x, y)) with z(x, y) := 1 2 (κ 1 x 2 + κ 2 y 2 ).The paraboloid Γ approximates the surface Γ with errors of the third order: O(x 3 , y 3 , x 2 y, xy 2 ).The Jacobian J(x, y) is going to be the norm of the normal vector to the surface J(x, y) = 1 + ∂z ∂x (x, y) By using this approximation of the surface and considering as the neighborhood M(x, r 0 ) the set M r 0 := P(x, y) : x 2 + y 2 ≤ r 0 we can rewrite the integral as x 2 + y 2 + 1 4 (κ 1 x 2 + κ 2 y 2 ) In the article [15] the function Ψ Γ,r 0 (x, y) = C Γ,r 0 is defined as a constant with respect to x and y: Then Finally, ∂G 0 ∂ny (x − y), x, y ∈ Γ can then be regularized as: (76) The same reasoning can be applied to the double-layer conjugate kernel, where in the previous calculations the expression of F is F (x, y) := 1 8π κ 1 x 2 + κ 2 y 2 x 2 + y 2 + 1 4 (κ 1 x 2 + κ 2 y 2 ) 2 where C HL Γ,r 0 = 2 )r 2 0 + O(r 4 0 ) .

Figure 1 :
Figure 1: Kernel restriction Visualization of the restriction of K(x * , y) = x * − y −1 to the unit circle and of the singular properties of K. Since K(x * , y) is singular at y * , K(x * , y) is singular along the line P Γ y = x * .One observes that the gradient of K(x * , y) −1 (and thus K(x * , y)) is orthogonal to the normal of the interface.

Figure 2 :
Figure 2: Singular behavior -1 Curvatures with same sign: same (first and second row) and different values (third and fourth row).Left column: double-layer kernel with target point in the origin, without any regularization.Center column: constant regularization.Right column: linear regularization.The second and fourth rows show the constant (center) and cappuccio (right) regularizations with r 0 halved compared to the first and third rows respectively.

Figure 3 :
Figure 3: Singular behavior -2 Curvatures with different sign.Left column: double-layer kernel with target point in the origin, without any regularization.Center column: constant regularization.Right column: linear regularization.The second row shows the constant (center) and cappuccio (right) regularizations with r 0 halved compared to the first row.

Figure 4 :
Figure 4: Singularity unaligned to the grid Position of the singularity point x 0 relative to the closest grid node x ∆ ; the parameters α, β are used to characterize its position relative to the grid.

Figure 6 :
Figure 6: Discontinuous behavior Left column: kernel multiplied by the singular function (51).Center column: function analytically found by taking the limit (53).Right column: plot of the function used in the additive splitting f /S − .

Figure 8 :
Figure 8: Torus Left: the torus used in the tests.Right: the torus and the projections of the Cartesian grid nodes inside the tubular neighborhood T ε .The projected nodes serve as the quadrature nodes.
is the composition of the three rotation matrices.The terms Q x (a), Q y (a), and Q z (a) are the matrices corresponding to a rotation by an angle a around the x, y, and z axes respectively.The parameters used for the rotations were:a = 0.2440241225550843 b = 0.7454097947651017 c = 0.2219760487439292 • 10 1

Figure 9 :
Figure 9: Sphere testsAveraged E 1 (h) errors (72) on 50 random target points on a sphere.Errors for the punctured trapezoidal rule (32) (black crosses) and the three considered methods in the evaluation of the double-layer potential: Q C constant regularization (26) (blue upward triangles); Q L cappuccio regularization (30) (red downward triangles); Q h corrected trapezoidal rule (67) (magenta circles).The three plots reflect three different settings for the tubular neighborhood width: left plot ε = 0.1; center plot ε ∼ h 0.7 ; right plot ε ∼ h 0.8 .

Figure 12 :
Figure 12: Error behavior and nonsmooth convergence Single target point on a "flat" torus.In the left half of the figure, it is shown the convergence behavior of the error for the double-layer potential, with two discretizations highlighted (upward and downward triangles).In the right half of the figure, top plots show the error behavior for the discretization h ≈ 4.57 • 10 −3 (downward triangle in the left figure) on the corrected planes as function of η(z) = n(z − z * ) with n = n 2 1 + n 2 2 + 1, which corresponds to the error 3.61 • 10 −8 .Bottom plots show the error behavior for the discretization h ≈ 4.66 • 10 −3 (upward triangle in the left figure) on the corrected planes as function of η(z) which corresponds to the error 5.03 • 10 −6 .The left figures show the signed error as function of η(z); the right plot shows the distribution of the (α(h, k), β(h, k)) values on the different planes, where the color represents the value of the error.From the mean µ and variance σ printed on top of the left plots we can see that the top case has mean much smaller than the bottom case, and the variance is half.This explains the much smaller error (downward triangle) compared to the other (upward triangle).

C r 2 F
(r 0 cos θ, r 0 sin θ)dθ then the conditions imposed form the following linear system:a 0 + a 1 = φ 0 a 0 r 0 C r 2 + a 1 C Γ = C DL Ffrom which we find:
for the smallest h ≈0.00437, first decreasing N to 11, and then decreasing N α,β to 51.The corresponding results are in the following table: