Semi‐analytical smoothed‐particle hydrodynamics correction factors for polynomial kernels and piecewise‐planar boundaries

This work presents a method of exactly and efficiently evaluating integrals of compactly supported (and otherwise) polynomial radial kernels near piecewise‐planar boundaries. The technique is motivated by and has an immediate application in the smoothed‐particle hydrodynamics method employing the semi‐analytical boundary formulation. The current implementation is exact and fully general, avoiding any need for symbolic or numerical integration which previous implementations rely on. A detailed derivation and three compact implementations are provided and the latter can be easily modified to fit new and existing simulation codes. Applications and good performance are demonstrated on a number of test problems, where the semi‐analytical boundary method can simulate complex boundaries using roughly half the boundary particles required in a ghost particle boundary method.

(A) (B) F I G U R E 1 Visualization of the semi-analytical integral in two dimensions. In (A), a fluid particle (green) near a domain boundary (black line) is shown. The center of the particle (black dot) is on the inside of the fluid domain. In (B), the problem is visually extended to 3D so that the particle density function can be plotted. The normalization factor ≤ 1 is the ratio of the volume under the green part of the surface (inside the fluid domain) to the total volume under the peak (the sum of the volumes under the green and red surfaces). The latter is a constant known a priori and typically set to unity by appropriately scaling the particle density profile. It is intuitively clear that even for general piecewise-planar (linear) boundaries integrating the volume under the green surface is a nontrivial problem, and a much harder one in three dimensions dimensions. Numerical methods capable of performing the brute force calculation are available, but their cost is generally way too high to be feasible. In recent years multiple authors addressed the issue of performing the underlying calculation in a practically viable fashion. The two-dimensional case was handled analytically by Ferrand et al. 7 for a single polynomial density function, the C 2 quintic Wendland kernel. 8 The three-dimensional case, for a special case of the density function, was first addressed by Mayrhofer et al., 9 later by Violeau et al., 10 and recently by Winchenbach et al. 11 These solutions are not fully general, do not handle nonconvex geometries correctly and rely on involved expressions computed via mathematics software. Additionally, Reference 9 employs a relatively complex and cumbersome geometric decomposition process, Reference 10 is incomplete, and Reference 11 is formulated for (and therefore accurate near) just infinite planar boundaries. An alternative approach was proposed by Chiron et al., 2 where the problem is reduced to a two-dimensional integral performed numerically on an unstructured planar mesh. This method correctly accounts for star-convexity but is again formulated for a single kernel (Wendland C 2 ), 8 uses a nonstandard low-order formulation, purposefully avoids handling planar regions bounded by curved edges (by partially integrating over trivial, noncontributing areas) and requires a meshing algorithm and a domain trimming process both of which are highly nontrivial. Our recent work 12 introduced a general, simple and compact method for arbitrary polynomial kernels where the original problem is reduced to a series of one-dimensional integrals over straight edges approximated numerically and the rest of the solution is fully analytic. The method is robust but still requires numerically integrating a well-behaved (but fairly cumbersome) function over a finite interval and a single user-specified parameter controlling its accuracy.
This work is a natural continuation and extension of Reference 12, with the part of the solution previously approximated numerically now evaluated analytically. In order to achieve this, the antiderivatives of interest are not formed explicitly but instead their values at a point are generated via inexpensive recurrence relations. Avoiding the numerical integration means the method no longer requires user-defined parameters (see in Reference 12 which controls the number of numerical integration nodes) to control the accuracy and computational cost. Instead, the method is fully exact, with computational cost controlled only by the number of boundaries, rather than the desired accuracy. No numerical quadrature need be imported or implemented, and thus as a whole, the implementation is much more user friendly.
The previously given algorithm 12 is appropriately updated and improved, in particular by minimizing the number of vector-vector operations (e.g., scalar products) and providing appropriate handling of particles along walls (but away from element edges, where the underlying problem becomes ill posed). We therefore provide the first fully general, computationally efficient and exact solution of the problem for arbitrary piecewise-planar boundaries and polynomial kernels. MATLAB, C++, and Python implementations of the three-dimensional version of the algorithm are provided and can easily be adapted to fit existing or new SPH codes. An additional application of the methodology in conjunction with the radial basis function interpolation method is also discussed.
It should be noted that this article does not address the formulation of the SPH governing equations under the SA boundary method, but rather focuses on the efficient and exact calculation of the required terms and ∇ . For a discussion of the full implementation of the SA method including how and ∇ are included in the governing equations, the reader is directed to Reference 7, which this article used as a basis for its own SPH implementation.
The rest of this article is organized as follows. The theoretical background is given in Section 2. The three-dimensional problem statement and its reduction to one-dimensional integration as in Reference 12 are briefly summarized in Section 2.1. Two special cases of line contours are motivated and discussed in Section 2.2 where the corresponding analytic solution for circular arcs is given as well. The essential part of this article, closed-form integration over straight lines previously performed numerically, is addressed in Section 2.3. A brief discussion on the closed-form two-dimensional solution follows in Section 2.4. The general two-and three-dimensional algorithms for computing the normalization factors and their gradients in SPH simulations are given in Section 3. Simplifications and an application of the methodology in conjunction with the radial basis function interpolation method are discussed in Section 4. Representative examples and conclusions follow in Sections 5 and 6 respectively.

METHOD
The three-dimensional problem and its reduction to one-dimensional integration for planar boundary elements 12 are summarized first. Two special cases are then motivated and introduced and closed-form expressions for both are given. In particular, integration along straight edges previously handled numerically is tackled analytically in Section 2.3. A short discussion on the two-dimensional version of the problem follows.

Reduction to 1D
In this section we summarize the problem and the derivation given in Reference 12. Let Ω be a compact subspace of R 3 with boundary Ω. The aim is to evaluate where n = n(x) is the unit outward normal to Ω and W = W(||x − x a ||) is a smooth scalar function dependent only on the distance between the argument x and some fixed point x a ∈ R 3 . The gradient ∇ in the above is taken with respect to x a , since is implicitly a function of the latter; however, in further calculation we treat x a as constant and the definition (1) of ∇ is used directly. We write so that W = W(r) is a so-called radial function (spherically symmetric around x a ). By applying Gauss' theorem and under mild conditions on W it can be shown 12 that so that both and ∇ are now expressed as surface integrals. Next, the case of piecewise-planar Ω is investigated. Let Ω j denote a planar portion of Ω where and Ω j are disjoint (they can only share edges). Clearly, where j runs over all boundary elements Ω j and j ∕∇ j are their respective contributions (surface integrals performed over Ω j ); it therefore suffices to consider a single planar portion of the boundary (a single j ∕∇ j pair). The plane of element Ω j is written P and x P denotes the projection of x a onto P. For any x ∈ P we have since n is constant along a planar boundary. Here is the signed distance between x a and its projection onto P. It also follows that for all x ∈ P and by defining = ||x − x P || the contributions j ∕∇ j as in (1), (3) can be written as By then applying Stokes' theorem the above can be integrated up once more, 12 resulting in where C j is the contour of Ω j oriented anticlockwise w.r.t. the (constant) element outward normal n, x = x(t) is an arbitrary parameterization of C j and ds, Finally, it is observed that for any polynomial W, written symbolically as the functions H (2,3) W can be given in closed form, 12 namely, ∑ a n r n = W(r) → H (3) with the expression for H (2) W being the same without the term (n + 3) in each denominator in (12). Only polynomial W will be tackled in the rest of this work; for a generic W we recommend using a polynomial approximation so that the methodology presented here can be applied.

Special cases
Two types of contours C j are of particular interest in this work. The first one is a straight line which enables integrating over domains bounded by polygonal (e.g., triangular) meshes, the latter being the typical way of representing geometries in numerical models. Second, we are also interested in (a part of) C j being an arc of a circle in the plane P of the element Ω j centered at x P . This shape occurs naturally in the SPH method, 12 the latter being the main motivation and application for this work. For any such circle = const. and since is also constant in the plane P, along such C j both H (2,3) W are constant. It has already been pointed out 12 that if C j is an arc with radius ⟂ spanning the angle in the plane P around x P then for k = 2, 3 with the sign of indicating the direction of integration around the normal n. Contour contributions of this type are therefore fully analytic for polynomial W. Previously, 12 the integral (9) along straight edges was suggested to be performed numerically. The reason is a nontrivial dependence of H (2,3) W on and , which together with itself being a nonlinear function of the natural parameterization makes it difficult to generate closed-form results. In the next section we address this problem by providing a procedure for generating closed-form antiderivatives of H (2,3) W and a method for evaluating them without constructing the analytic expressions explicitly.

Closed-form integration along straight edges
Consider a straight line segment C (i) j of C j in the plane P spanned by z 1,2 in this order, so that the unit direction along the segment isê The projection x l of x P onto the line l passing through z 1,2 is then given by and in the following we employ the arc length parameterization of C (i) where Only the computation of j is addressed explicitly in algebraic detail; the derivation for ∇ j is fully analogous and the corresponding results are stated when appropriate. With the above parameterization the contribution to j (9) along C (i) j , written (i) j , becomes being the signed distance between x P and x l ( Figure 2) and trivially from Pythagoras' theorem We point out that both (17) and (19) can be evaluated without constructing x P and x l explicitly, but with employing just the center of the radial function and the end points of the edge; such treatment makes implementations more efficient.
The case , ≠ 0 in (12) is tackled first. By direct substitution to (18) and with the change of variables t = | |s it follows that Next, for > 0 we define (i.e., an indefinite integral/antiderivative w.r.t. x) so that with the above notation (and analogously for ∇ (i) j ) F I G U R E 2 A spherical isosurface of a radial function (green) in the vicinity of a straight edge (black) in a given plane P (gray mesh).
The signed distance between the center of the function x a and its projection x P onto the plane is plotted in light blue. The signed distance between x P and its projection x l onto the line l containing the edge is plotted in red. For x along the edge we have ||x − x a || 2 = 2 + 2 = 2 + 2 + ||x − x l || 2 = 2 + 2 + t 2 . In our notation t is the appropriately oriented arc length parameterization of l with t = 0 corresponding to x l where in the limit → 0 which follows easily from tan −1 (x) ∼ x for |x| ≪ 1.
It is important to note that while closed-form expressions for I of interest can now be constructed explicitly, they will generally be cumbersome and involved. A different, in particular a much more efficient approach is highly recommended in implementations. First, (27) is applied "pointwise" in the following way. For a general polynomial W (11) and given x, 2 (e.g., as in (23)), (28) are evaluated if W has any even degree terms and (29) are evaluated if W has any odd degree terms. Any other required I (with the maximum dictated by the degree of the polynomial W) are evaluated by repeatedly applying (27) to those initial values while treating x and as constant. In this way the series of values is generated in a sequence of simple and cheap steps only requiring a small number of multiplications. Conversely, explicitly constructing I starting from (29), say, would make each subsequent function increasingly involved and result in a much more computationally expensive evaluation (cf. References 9,10).
Second, note that iterating the first equation in (27) for positive integer k yields I k = I 1 + (J 1 + … + J k−1 ), where in the above the dependence on x, 2 (treated as constants) is implicit for simplicity of notation. The above allows for sequentially computing a series of values of I by only iterating the self-contained recurrence relation for J in (27) and storing a single value of J in memory. Note that the case = 0 in (12) is relevant only if integration is performed numerically (so that 0∕0 does not occur in implementations). We therefore only need to consider the last remaining option = 0. In this case (i) j ≡ 0 since multiplies the first integral in (8) (more generally j ≡ 0 for = 0) and with the same parameterization and notation direct substitution yields In the original definition (24) of J the second argument was nonnegative which is not guaranteed in the above. However, J defined as the integral in (24) is well defined for the second argument ≥ −1, which makes the notation in (32) appropriate. In this case only the self-contained recurrence for J in (27) (second line) is used, together with J 0.5 in (29) and the limits where now the recurrence for even n starts from = 0 (not = 1 as before).

Two dimensions
The two-dimensional version of the problem can be tackled in a similar fashion. The two integrals ∕∇ are defined as in (1) but in two dimensions. It has already been discussed 12 that where n = n(x) is the outward normal to Ω, x = x(t) is an arbitrary parameterization of Ω, = ||x − x a ||, P is the two-dimensional + 2 rotation matrix and For polynomial W the above can be given explicitly For piecewise-linear boundaries in two dimensions it suffices to consider a single line element Ω j for general piecewise-planar geometries, including in the SPH method (i.e., there is no equivalent of integrating over circular arcs). As in three dimensions, the line segment of interest is spanned by z 1,2 in this order, the unit direction is as in (14) and x P → x a in the definition (15) of x l . The parameterization (16) and the range (17) are analogous. By direct substitution to (34) it follows that where n is constant along linear segments and is the signed distance between x l and x a . With polynomial W (36) and our previous notation and the evaluation process is the same as in the special case = 0 in three dimensions.

SMOOTHED-PARTICLE HYDRODYNAMICS
The closed-form integration method given in the previous section allows for updating the two-and three-dimensional algorithms given in Reference 12 in order to make them fully analytic. In the SPH method the radial function W satisfies the additional conditions where Ω a is the particle support, for example, a sphere (circle in 2D) with a known radius R centered at x a . By exploiting the above it follows that for piecewise-planar boundaries and for with the index j running over all boundary elements Ω j such that Ω j ∩ Ω a is nontrivial; detection of these is an essential part of SPH codes and is not discussed in further detail. The contributions j ∕∇ j in (41) are integrated over Ω j ∩ Ω a (or their contours in 3D), A 0 is the area of the unit sphere (4 in three dimensions, 2 in two dimensions) and A j is the area of the projection of Ω j ∩ Ω a onto the unit sphere (circle in 2D) centered at x a . In a more general case, when in (40) the first equation is replaced with the first equation in (41) is accordingly replaced with and we note that for polynomial kernels ∑ a n r n = W(r) → M = 4 R 3 ∑ a n R n n + 3 .
Formulation (41) splits the evaluation of into the integral over the whole support M (usually M = 1 in SPH codes), addition of piecewise-planar element contributions j and correcting for the region of Ω a obscured by Ω j ∩ Ω a . The sign factor sgn( j ) guarantees the correct result in cases when the domain of integration is not star-convex w.r.t. x a and here is defined at zero via sgn(0) = 1. Evaluation of ∇ in (41) relies only on direct contributions from boundary elements due to the natural second condition on W in (40). In the rest of the section we discuss computing j , A j and ∇ j for a single piecewise-planar (linear) element Ω j in the context of implementation in SPH simulators.

Two dimensions
Let a linear boundary segment spanned by x 1,2 in this order be given. Generally, the segment may be (partially) outside Ω a , where W is identically zero. These regions are trivial and do not need to be integrated over and will be accounted for as previously. 12 First, the unit tangent is given byê and then we compute the parameterization limits t 1 and t 2 in (16) corresponding to x 1,2 , namely Next, the intersections of the line passing through x 1,2 and Ω a are investigated. With the parameterization (16) the intersection condition R 2 = ||x(t) − x a || 2 can be rewritten as which describes a segment contributing in a nontrivial way only if 2 < R 2 (the contributions j , A j and ∇ j are all zero otherwise). Under this assumption we then truncate the parameterization [t 1 , and proceed only if after the truncation t 2 > t 1 -otherwise the segment is outside of the particle support and all its contributions are zero. With the definition of in (47) and after the truncation (48) the element contributions j ∕∇ j for polynomial W are as in (39), and the area of the projection of the truncated segment onto the unit circle centered at x a is given by with only the case = t 1 t 2 = 0 not well defined given the current information. The latter corresponds to the fluid particle being precisely at a corner of Ω and requires information about both segments connected by the corner. A simple safety mechanism dealing with this issue was described in Reference 12, where the particle is moved away from the corner by a tiny amount in the negative mean direction of the surrounding normals. However, in practice we would not expect x a ∈ Ω to occur in simulations (particles stay away from the walls). An efficient implementation of the proposed method for polynomial W is summarized in Algorithm 1 below. The implementation is valid for all x a inside Ω and along Ω with the exception of corner points (detectable via = t 1 t 2 = 0).

Algorithm 1.
Pseudocode for the two-dimensional algorithm for polynomial W with the support radius R Data: Appropriately oriented end points of a line segment Ω j , x 1,2 and the outward unit normal n Result: Element contributions to ∕ , j ∕ j , and the area A j to be used in (41)/(43) Evaluate j ∕ j as in (39)

Three dimensions
Let Ω j be a convex polygon in the plane P with corners x 1, … ,n oriented counter-clockwise around the outward normal n. Due to the finite size of the spherical support Ω a the general nontrivial intersection Ω j ∩ Ω a will be a planar element with its edge described by a combination of straight lines and circular arcs.
First, is computed from where x 1 can be replaced by any other corner of the polygon. The element contributes in a nontrivial way only if R 2 − 2 > 0, which is assumed, and is skipped otherwise. The intersection Ω a ∩ P is a circle centered at x P with radius Next, we consider a single edge l i of the polygon, spanned by x i,i+1 in this order (x n+1 = x 1 ). For this edge we computê and Again, the segment has a nontrivial contribution only if 2 ⟂ > 2 . Also note that x P is on the left of l i in plane P if ≥ 0 and vice versa. For a contributing segment we evaluate and truncate while recording if any of t 1,2 was truncated from either above or below by the min/max functions. As in two dimensions, the segment contributes only if t 2 > t 1 and is skipped otherwise. The contributions to j ∕∇ j from the edge l i , written (23) with the truncated parameterization (54) and , as in (50), (53). For any contributing line segment we also calculate the arc length (i) of its projection onto the unit circle in the plane P centered at x P , given by the same expression as in (49) but with the values of , t 1,2 as in (53) and (54).
After all sides l 1, … ,n are considered, we initialize where the index i * runs only over contributing segments (i.e., those with 2 + 2 < R 2 ) and (i * ) is the value of as defined in (53) along the (contributing) edge l i * . The sign function in the above is again defined at zero as sgn(0) = 1. The above values of j ∕∇ j are not yet exact, since they miss contributions from potential circular arcs along C j . First, if and only if along all contributing straight edges ≥ 0 then is updated according to which indicates that x P ∈ Ω j (if Ω j is a convex polygon). In order to account for circular arcs we finally update which correctly evaluates the contributions j ∕∇ j to be used in (41). The last remaining part of the solution is the evaluation of the projection area A j of Ω (a) j = Ω j ∩ Ω a onto the unit sphere centered at x a . The general shape of Ω (a) j consists of a convex polygon and a number of sections of the circle O P with radius ⟂ in the plane P centered at x P . The corners of the polygon can be constructed from the parameterizations (54) along the contributing edges and sections of the circle can be detected and characterized by investigating if (54) was truncated, as mentioned previously. The projection area of the smaller segment of O P bounded by the chord with length 2d is given by which is equivalent to the formula given in Reference 12 with the additional handling of = 0, d = ⟂ . The above can be expanded to provide an expression for the projection area of either the larger or the smaller section of O P , where the chord is spanned by y 1,2 in this order with positive if the chord goes around x P in the anticlockwise fashion (corresponding to the smaller area) and vice versa.
In the above we can either again use sgn(0) = 1 or the typical convention sgn(0) = 0. Finally, let be an array of projections of the corners of an arbitrary convex planar polygon onto the unit sphere centered at x a (||c k || = 1 for all k = 1, … , m). The ordering of c k should correspond to an arbitrary consecutive ordering of the corners of the original polygon. The area of the spherical polygon spanned by c k (either equal to A j or directly contributing to it) is given by with the above equivalent to the expression provided and discussed in further detail in Reference 12. Note that in implementations it is preferred to evaluate just p 1, … ,m first and then use them in the above formula in order to minimize the computational expense (the naive implementation would evaluate each p k twice). Also observe that necessarily m ≤ 2n where n is the number of sides/corners of the boundary element Ω j ; the latter bound may be important in implementations if the amount of memory allocated to C has to be specified a priori. Many of the terms involved in the calculations presented, particularly (61), produce singularities for certain inputs. While these inputs should in practice never be encountered, our experience shows that finite machine precision can often lead to situations where the function evaluations do become singular. This is because the formulae presented involve repeated floating point operations which can produce results where the values are indistinguishable to machine precision. This is often the case when an element corner is relatively close (within 10 −10 ) to Ω a . In order to avoid floating point errors, in (61) we first explicitly limit to ensure that the argument of √ … in the denominator of each term in (61) is nonnegative. Second, each term in the sum is modified according to where is a small positive user-specified tolerance on the order of machine precision; in our implementations = 10 −14 is used. The above modifications directly prohibit singularities from occurring and guarantee that the argument of cos −1 is in the appropriate range [−1, 1]. Similar but simpler safeties are used in (58) for the arguments of sin −1 and cos −1 .
Since both of them are positive and ⟂ , R 2 − d 2 > 0 are already guaranteed when A(d) is evaluated, each argument is only limited from the top via min(1, ⋅). An explicit mechanism for avoiding singularities in (58) from division by zero can also be used (cf. max( , ⋅) in (63)). The above methodology is visualized in Figure 3. The main view and the three projection angles (a − c) are further supplemented by an interactive U3D model * .
An efficient and fully general implementation of our approach for polynomial W is summarized in Algorithm 2 below. MATLAB, C++ and Python implementations of the general three-dimensional algorithm together with examples corresponding to Figure 3 are provided in Appendices A, B, and C, respectively. Similarly to 2D, the method handles arbitrary x a within Ω and along Ω with the exception of regions of nonsmoothness, namely the element edges and corners, that is, all cases when = 0 and t 1 t 2 ≤ 0. In those cases the projection area of the element A j onto the unit sphere is not well defined and evaluating the appropriate contribution requires information on neighboring elements. The previously mentioned safety mechanism 12 applies in such scenarios. * 3D content activates on click. The U3D format is supported by Acrobat Reader (version 7 and above). For more information visit https://helpx.adobe. com/acrobat/using/enable-3d-content-pdf.html.

F I G U R E 3
Visualization of the three-dimensional problem where a particle (green) with radius R = 1 (so that Ω a is a unit sphere) interacts with a quadrilateral boundary element Ω j . The contour C j of Ω j ∩ Ω a is composed of three straight lines (yellow, solid) and two arcs (red) with their end points connected by yellow dashed lines (also see views (a − c) for more detail). The projection area A j (light blue) is the sum of the projection areas of the pentagon bounded by all yellow segments (dark blue contour, cf. A [C] in (61)) and two sections of the circle bounded by red arcs and dashed yellow lines (here two spherical regions between red and blue arcs, A (s) in (59)) The correct handling of , = 0 in implementations is crucial, even if the particle position x a is always assumed to be away from walls. Indeed, for x a along the line containing a straight edge of an element within the particle support (but away from the edge itself and Ω generally) one has = = 0 but can be strictly on the inside of Ω (only if Ω is not convex).
Equations (41), (43), with contributions j , ∇ j , A j as defined in this work, produce correct results for x a ∈ Ω (with the exceptions given before). For x a strictly outside Ω, and without any change in the implementation, the resulting value of ∇ will still be computed correctly. However, in such cases the computed and true values of , written C , T respectively, satisfy and C > M if x a ∉ Ω. The latter condition is unlikely to occur and generally not desired in simulations. Still, it may potentially occur and be hard to fully eliminate in extreme cases such as very high-velocity impacts. A simple general solution ensuring the correct value of is setting after the calculation (41)/(43) involving all contributing boundary elements.

SIMPLIFIED ALGORITHMS
The source of singularities (and the main source of complexity in three dimensions) in Algorithms 1 and 2 is evaluation of the projection area A j . The latter is necessary for the evaluation of (Equation 41) and is essentially an indirect way of integrating over the spherical region Ω a ⧵ Ω.

Algorithm 2.
Pseudocode for the three-dimensional algorithm for polynomial W with the support radius R. All of j , A j and ∇ j are computed in the same loop. An additional switch can be added after line 15, so that t 1 t 2 ≤ 0 = raises an error Data: Appropriately oriented corners of a convex polygon Ω j , x 1,…,n and the outward unit normal n Result: Element contributions to ∕ , j ∕ j , and the area A j to be used in (41)/(43) Evaluateê and as in (52)  Update j as in (57);

end
However, Ferrand et al. 7 proposed a formulation where only ∇ is needed and an ODE in time for is solved numerically. With their approach the calculation of A j can be omitted and our method simplifies significantly in three dimensions as summarized in Algorithm 3 below. Observe that the same approach in two dimensions (Algorithm 1) leaves the corresponding algorithm practically unchanged since the calculation of A j in two dimensions (49) is evaluated in a single inexpensive step. In both two and three dimensions the simplified algorithm is valid for all x a , including along the edges and at the corners of boundary elements. At the same time, direct calculation of is very desirable in implementations and naturally leads to much simpler and more stable numerical schemes; we therefore recommend the full method to be used in practice.
Finally, we note that the methodology derived in this work can be used for the purpose of more generic efficient numerical integration over two-and tree-dimensional piecewise-planar domains. In order to do so, a generic function f can be approximated via the Radial Basis Function (RBF) interpolation method with a polynomial basis function of infinite support (e.g., W(r) = r 3 , the cubic spline), written as where w k are constants andx k are a prespecified collection of RBF centers. Typically,x k are specified a priori and correspond to points at which the value of f is known. The set of weights w k is then constructed by solving a N w × N w linear system equivalent to enforcing the agreement between the RBF surrogate and the original function. The details of the RBF interpolation method are not crucial for this work and for further information the reader is kindly referred to Reference 13.
Once the surrogate (66) is given, we approximate . Each of the integrals on the right-hand side of Equation (67) can then be evaluated via the methodology proposed. If W has infinite support then the decomposition (5) applies, and [f ] in (67) can be evaluated in a single loop over the elements of Ω j , but at each element a series N w contributions from N w centersx k are evaluated. Note that the truncation of t 1,2 and the contribution of circular arcs (controlled via in our notation) are no longer applicable, further simplifying the implementation. The main advantage of the approach is that only the basic boundary mesh information can be used to evaluate three-dimensional integrals without complex meshing of the interior of the domain or any degree of approximation other than constructing the RBF surrogate for f . Moreover, the method resolution (here the number of RBF centers N w ) can be adjusted just w.r.t. the complexity of f rather than the domain of integration Ω. Pseudocode for the appropriately simplified three-dimensional method is given in Algorithm 4. We note that a very similar approach can be used to integrate scalar functions over (open and closed) piecewise-planar surfaces in three dimensions. The method would replace j with ∇ j in Algorithm 4 but with the multiplication by −n omitted in both expressions for ∇ j in (23) and (32).

RESULTS
In this section example applications of the methodology proposed in this work are presented, both in the context of SPH simulations and numerical integration with the aid of the RBF interpolation method.

Smoothed-particle hydrodynamics
Generality and accuracy of the semi-analytical formulation have already been derived and established, and the discussed method of computing the correction factors ∕∇ is exact (maximally accurate). The crucial remaining issue is the numerical efficiency, which if low, may provide an incentive to still use more widely spread and less accurate heuristics.
In order to quantify the performance of the approach, the supplied generic C++ implementation (Appendix B) was directly ported to an existing SPH code, originally employing the (most common) fluid extension (aka ghost particle) boundary method. 14 Both the original and modified (semi-analytical) simulators share the same base code and differ just in boundary handling, ensuring a fair comparison in the context of interest. The SA method implemented mostly followed the methodology of Ferrand et al. 7 with some differences such as the time integration scheme and boundary pressure calculation in order to be as close as possible to the existing ghost particle implementation. Numerical efficiency was quantified on a typical SPH benchmark test case visualized in Figure 4.
Five fluid resolutions between 3000 and 100,000 particles were investigated and the C 2 quintic Wendland kernel 8 was used. A direct comparison of simulation times for both formulations is given in Figure 5. At each resolution two or three independent simulations are run for both boundary methods in order to provide a basic statistical characterization ( Figure 5, bottom plot).
The semi-analytical code is practically as fast as the original (fluid extension). On average, in our test it runs slightly faster at low (<10,000) and high (>90,000) resolutions and slightly slower at the intermediate ones. However, in all cases the difference is practically negligible-the relative speed of both methods is within ±3%, making them practically equal in this respect. It is noted that this level of performance was achieved with minimal modifications of the source supplied in Appendix B. The only difference between the provided and our code is the use of an already in place framework for handling vectors (including dot and cross products). Further increase in speed and efficiency are easily achievable; the latter two were not the only design criteria for the implementations provided, with generality and easy portability also playing a crucial role. If better performance is desired, we first recommend reducing or eliminating the need for dynamic memory allocation and precomputing and hard-wiring (i.e., fixing at compile time) the coefficients a n (n+2)(n+3) , and so on instead of passing a n 's for each element-particle pair at run time.
Finally, the generality and robustness of the semi-analytical framework is illustrated on a coastal defense example ( Figure 6, on-click interactive U3D model). Here the simple obstacle from the previous test is replaced with a breakwater comprising of five irregularly spaced pieces of rock armor fixed in place. Note the complex shape of the domain in the vicinity of the defense structure.
Geometric configurations such as the one presented in Figure 6 are often too complex to be handled efficiently by typical techniques such as the fluid extension method. For example, the discretization of a single breakwater such as the one in Figure 6 can be seen in Figure 7. The SA method at most discretizes just the surface of the walls, while the fluid extension method requires ghost particles to fill the boundaries up to a fixed distance, which is difficult and inefficient for such complex and high-resolution geometries. F I G U R E 7 Discretization of a single breakwater by different boundary methods. Image (A) is the stl geometry file used to initialize the method. Image (B) presents the discretization used in the SA method, which places ghost particles along the surface of the geometry. In this case a total of 4223 particles are required. Image (C) presents the discretization used in a typical fluid extension method, which fills part of the geometry with ghost particles. In this case, a total of 8016 particles are required, almost double that of the SA implementation In the example in Figure 6, the fixed ghost particle method requires a total of ∼350,000 particles to fully model the boundaries, while the SA method at most requires only ∼185,000 particles. The boundary particles in the SA method are required only if placing a single layer of boundary particles to improve the calculation of gradients and/or shear forces near the wall. For example, the SA formulation of Ferrand et al. 7 requires the surfaces of the walls to be discretized by boundary particles in order to improve calculations of gradients near the boundary. However, the formulation of Kulasegaram et al. 5 does not explicitly require any boundary particles. The implementation in this work is agnostic to the presence of boundary particles and can be equivalently employed by formulations with or without them. In this work, and in general, a Ferrand style formulation is used, but if there is a focus on speed and minimizing the number of particles, then a Kulasegaram style formulation can instead be used. Special note is given to Ferrand et al.'s paper which does present both formulations. 7 The full setup process is minimal and comprises of just providing a generic triangulated mesh file (cf. Figure 6) and specifying the initial body of fluid. A number of illustrative simulation snapshots are given in Figure 8. While not presented here, these results are comparable in overall dynamics to those produced by a typical ghost particle method. Other literature addresses direct comparison of the two methods. 7

Volume and surface integration
For demonstration purposes we consider the Rosenbrock (banana) function, given by over the domain Ω 0 = [−1, 1] 3 . As suggested in Section 4, an RBF interpolant of (68) was first constructed, with the samples taken over a regular 8 × 8 × 8 grid (Figure 9) and the cubic spline kernel function W(r) = r 3 . The relative L 2 error between f and its interpolant over Ω 0 (estimated numerically via dense regular sampling) is under 10 −4 . After f is approximated as a weighted sum of polynomial radial functions, integrating it over any structure (surface/volume) described by piecewise-planar polygonal elements is exact and cheap. The only source of errors is the inaccuracy of the interpolant, which can be decreased easily to arbitrarily low values via more dense sampling. In Figure 10 two examples of numerical integration of the interpolant via the method proposed are given. Note that the area and volume of any shape can be evaluated in the special case W(r) = 1 (a 0 = 1, a n≥1 = 0) centered at an arbitrary x a . F I G U R E 8 Visualization of the sea defense example. As in the previous problem, the body of water is released from rest at the back end of the container and splashes against the now highly nontrivial obstacle. The complex interaction between the fluid and five pieces of rock armor is handled automatically by the semi-analytical approach, which only requires the basic mesh information (here a triangulated geometry in the .stl format)

CONCLUSIONS
A general closed-form method for exactly evaluating integrals of polynomial radial kernels in two and three dimensions was given. The approach does not utilize any degree of numerical approximation and avoids explicitly constructing involved antiderivative expressions by employing pointwise recurrence relations. Crucially, arbitrary nonconvex domains are handled correctly by the method. The original motivation for the technique is the semi-analytical formulation of the SPH method, where the calculation discussed is responsible for accurate and general modeling of wall boundary conditions. Alternative applications, namely volume and surface integration by employing polynomial RBF interpolants, were discussed and examples of all proposed applications were given. Fully self-contained implementations of the general three-dimensional method were supplied in three popular programming languages (Appendices A-C) and can be easily adapted to work with new or existing SPH simulators (via Equation 41). We remind the reader that the (cheap) correction (65) can be used in implementations for full generality.
Our approach is already fully analytic, hence further accuracy improvements are not possible. Conceptually and implementationwise, the method is designed to be compact and efficient; involved calculations are performed only in cases where an edge of the boundary is within particle support. Source codes provided can be easily optimized further, for example, by employing precomputations or clever memory handling. Optimizations of this kind are implementation-specific and beyond the scope of this work, but have significant potential in particular in applications where execution speed is crucial. Still, even without such tweaks a direct port of our generic C++ source (Appendix B) was demonstrated to be practically as fast as the most common alternative (fluid extension) when coupled with the same base SPH code; the generality, stability and accuracy of the semi-analytical framework has already been achieved without a performance drop.
This work opens the possibility for the semi-analytical formulation to finally become a widely available and easily accessible standard. With the method and source codes provided, all semi-analytical formulations become vastly easier and more natural to implement, since the correction terms ∕∇ in the governing equations can now generally be evaluated directly, cheaply, and exactly (i.e., can be treated as "generally known"). While there is still room for improvement on the semi-analytical formulation itself, including how the governing equations are formulated, in this article we have effectively solved the issue of calculation of ∕∇ which had previously been the main issue holding the method back. Additionally, our method can be completely insensitive to the boundary element shape and size. While some semi-analytical SPH formulations prefer to include a layer of boundary particles and rely on a fine mesh for accurate estimates of ∕∇ , this is not strictly necessary. For example, the formulation of Kulasegaram 5 is independent of the existence of boundary particles, and the method presented here can be used either with or without boundary particles. In the case where none are used, the improvements in speed of the method can be even greater, with no loss in accuracy beyond the accuracy of the initial SPH formulation. This would allow any mesh suitable for the geometry of interest to be used independently of the fluid resolution. Typical methods require the smallest boundary element to have an area greater than dx × dx, but this need not be the case with the semi-analytical approach.
This article solves the general case of polynomial radial kernels and piecewise-planar boundaries comprising of arbitrary convex polygons. Polynomial kernels are by far the most common type used in practice and can be used as approximations of arbitrary (sufficiently well behaved) nonpolynomial alternatives, so limiting the investigation to this class of functions is not a major issue. However, a lot of work remains regarding the evaluation of correction terms near nonplanar boundaries. One of the main strengths of the semi-analytical formulation is that it applies to arbitrary boundaries as long as the corrections can be evaluated. This principle has been discussed and employed in our recent work 15 with promising results but so far only the simplest nonplanar shape, a sphere, has been handled as a proof of concept. Further development of methods of evaluating semi-analytical corrections near nonplanar walls remains an active area of our research.

ACKNOWLEDGMENT
The authors thank Professor Stephen Neethling (Department of Earth Science & Engineering, Imperial College London) for providing the SPH code upon which the current methods were implemented.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.