An unfitted RBF-FD method in a least-squares setting for elliptic PDEs on complex geometries

Radial basis function generated finite difference (RBF-FD) methods for PDEs require a set of interpolation points which conform to the computational domain $\Omega$. One of the requirements leading to approximation robustness is to place the interpolation points with a locally uniform distance around the boundary of $\Omega$. However generating interpolation points with such properties is a cumbersome problem. Instead, the interpolation points can be extended over the boundary and as such completely decoupled from the shape of $\Omega$. In this paper we present a modification to the least-squares RBF-FD method which allows the interpolation points to be placed in a box that encapsulates $\Omega$. This way, the node placement over a complex domain in 2D and 3D is greatly simplified. Numerical experiments on solving an elliptic model PDE over complex 2D geometries show that our approach is robust. Furthermore it performs better in terms of the approximation error and the runtime vs. error compared with the classic RBF-FD methods. It is also possible to use our approach in 3D, which we indicate by providing convergence results of a solution over a thoracic diaphragm.


Introduction
Most localized radial basis function (RBF) methods for computing solutions to partial differential equations (PDEs) (for example [1,2,3]) require a set of interpolation points that conforms to a computational domain Ω.Normally, a localized RBF method uses a collection of local interpolation problems over the subsets of interpolation points (stencils, patches,..) to generate the compactly supported cardinal functions which are then employed to solve a PDE.It is well known that locally non-uniform node distances between interpolation points increase the conditioning of the local interpolation problem and cause unwanted growth of the cardinal functions [4,2].While it is trivial to place uniform points in the very interior of Ω, it is on the other hand challenging to place them in the vicinity of a boundary of Ω and at the same time maintain their uniformity.This is especially difficult in three dimensions.This motivates the investigation of possibilities to decouple the interpolation points from Ω and avoid those restrictions.
The so called immersed methods (also unfitted methods) [5,6,7] which are a part of the finite element methods address the decoupling of Ω and a mesh: the boundary of Ω is enclosed in a box with a background mesh, where only the elements which have a non-empty intersection with Ω are taken to be active.Those methods tend to suffer from ill-conditioning in the presence of small and irregular cuts close to the boundary of Ω.This has for example been addressed in [8,9,10].An additional challenge is the enforcement of Dirichlet boundary conditions, which has been addressed by introducing a penalty term over the boundary elements [11,12].
Another related approach is the placement of ghost points (also fictitious points) in finite difference methods, where additional points (unknowns) are placed outside of Ω in order to enforce Figure 1: Node distributions over a butterfly.Left: a classical distribution of interpolation points over Ω for RBF-FD in the least-squares setting with two skewed stencils on the boundary.Here the black point is the stencil center and the red points are the members of the stencil.Right: a node distribution over Ω for the unfitted RBF-FD method in the least-squares setting with a less skewed stencil on the boundary.
Neumann type boundary conditions in a more accurate way.This is a widely used concept, an example can be found in [13], where for every added point, an additional equation is generated in order to maintain a square linear system of equations.In [14,15] the authors introduce ghost points for a global radial basis function (RBF) collocation method in order to decouple interpolation points from Ω. Their computational study shows that this is a feasible approach and that the error under node refinement tends to be smaller compared with the fitted method.However the study is limited to using basis functions with a global support, and the study does not provide an insight into how many ghost points to use and how that affects the stability properties.
A partition of unity based RBF method in a least-squares setting (RBF-PUM-LS) [16] enables a decoupling of the interpolation points from Ω by placing a set of overlapping patches over Ω, where every patch contains interpolation points independent of the shape of Ω.The authors provided numerical evidence that RBF-PUM-LS is an accurate and robust method to solve an elliptic model problem, but have not studied the effects of the patches that extend outside Ω.
A recently introduced RBF-FD method in a least-squares setting (RBF-FD-LS) [17] was proven to be significantly more robust compared to the same method in the collocation setting (RBF-FD-C), especially in the presence of Neumann-type boundary conditions [17].However the interpolation points are required to conform to Ω. Another study leading to a least-squares RBF-FD was introduced in [18].
In this paper we use the RBF-FD-LS method [17] and for that method introduce an approach to decouple interpolation nodes from Ω.A computational domain Ω is enclosed in a box that contains a set of regularly spaced interpolation nodes with reasonably good interpolation properties (see Figure 1 and Figure 2).The solution unknowns are determined by solving a system where every equation is an evaluation of a PDE in a point y ∈ Ω.Here every y picks the closest interpolation point which is used as a center of an approximation stencil built over a set of neighboring interpolation points in a box.The shape of such a stencil is given in Figure 1 and Figure 2.An additional strength of the presented approach is that the approximation stencils around the boundary become less skewed which is found to reduce the magnitude of the error around the boundary, especially when the stencil size is large.
The paper is organized as follows.In Section 2 we state the model problem.In Section 3 we provide a description of the unfitted RBF-FD-LS method, together with the formulas for computing local differentiation weights, the global differentiation matrices and the discretization of a model problem.In Section 4 we study linear independence of the cardinal functions as the interpolation points move away from the boundary of Ω and develop a heuristic criterion to keep the linear independence of cardinal functions unchanged.This is a necessary condition for the well-posedness of the discrete PDE problem.In Section 5 we study the relation between the discrete solution and the analytic solution by deriving a discrete error estimate without an a-priori bound on the stability norm.The 2D experiments are presented in Sections 6 and 7.In the former section we consider a butterfly domain and numerically investigate the behaviours of the error against true solution, stability norm and the condition numbers under node and polynomial degree refinements.The results are compared against RBF-FD-LS and RBF-FD-C.In the latter we use a drilled 24tooth sprocket as a computational domain and study the effects of the unfitted discretization on the spatial distribution of the error for a fixed internodal distance and several polynomial degrees.
In Section 8 we numerically study the convergence under node refinement in 3D, where we use a thoracic diaphragm geometry extracted from medical images.Lastly, Section 9 concludes the paper and offers directions for further work.

The model problem
We choose to evaluate our method by solving the Poisson equation on an open and bounded domain Ω ⊂ R d (d = 2, 3) with mixed boundary conditions.
where ∂Ω 0 and ∂Ω 1 are two disjoint parts of a smooth boundary ∂Ω.The solution u is throughout the paper assumed to be smooth.In the theoretical parts of the paper we prefer to work with the following formulation of the same problem: where: The numerical solution is going to be sought using: where Ψ i (y), i = 1, .., N are the RBF-FD cardinal functions and u h (x i ) are the nodal values of the solution.

The unfitted RBF-FD method
In this section we discuss the choice of point sets that discretize the computational domain, the generation of the cardinal functions (4) using the RBF-FD method, the assembly of the evaluation and differentiation matrices and the discretization of the model problem (1).

The point sets
Two sets of computational points are distributed over Ω: • The interpolation point set X = {x i } N i=1 for generating the cardinal functions.
• The evaluation point set Y = {y i } M j=1 for sampling the PDE (1).We noted in [17] that the X-points are supposed to be distributed such that the internodal distance is as uniform as possible, since the Lebesgue constants associated with the cardinal functions then stay fairly small.On the other hand the evaluation point set does not influence the magnitude of the Lebesgue constants but is instead important for the implicit integration that occurs when solving a discretized system of equations in the least-squares sense [17].Thus the constraints for placing the evaluation points are far more forgiving.
We choose the interpolation point set such that it does not conform to the computational domain Ω (see Figure 1 and Figure 2).Throughout the paper we take X to be a hexagonal grid with spacing h.The motivation for that is a simplified point generation, and the benefit of a polynomial unisolvency on those points.The latter is important for the well-posedness of an interpolation problem over a stencil.
The evaluation point set Y conforms to Ω.We choose the interior part of Y such that there are q points in every Voronoi cell centered around each x ∈ X.The boundary points are then placed on ∂Ω with a uniform distance that corresponds to the distance between the interior points.For a visual representation see Figure 3.
With such a relation between the X-and the Y -points it follows that the cardinality of those sets very closely matches the relation: M ≈ qN , where M is the number of Y -points and N is the number of X-points.

The RBF-FD trial space
Let Ω i be a subdomain holding a collection of points X Ωi = {x j } n j=1 that are a subset of the interpolation points placed on top of the computational domain Ω.Then every stencil is defined by a tuple (Ω i , x i ), where x i ∈ X Ωi is the stencil center point (see Figure 1 for a visual representation of a stencil).The solution u h (x) over a stencil is spanned by a combination of cubic polyharmonic splines φ l (x) = ||x − x l || 3 and multivariate monomials {p k } m k=1 of degree m [19,20,21]: subject to where n is the stencil size, c l are the interpolation coefficients and β k are the Lagrange multipliers.
The coefficients c l and β k from ( 5) are computed by requiring the interpolation conditions u h (X Ωi ) = u(X Ωi ), where u(X Ωi ) := u are the stencil data coefficients.This results in a square system of equations: Here A jl = φ l (x j ) for indices j, l = 1, .., n and P jk = p k (x j ) for index k = 1, .., m and u j = u(x j ).The stencil-based solution in any point x ∈ Ω i is expressed by reusing the computed coefficients c and β from (6) within the linear combination (5): where Ψ 1 (x), .., Ψ n (x) are the stencil-based cardinal functions with the Kronecker delta property and w are the local stencil weights.The next step is to use the formulation from (7) to represent the solution over the whole computational domain Ω.First the coefficients in ( 6) are computed for every stencil (Ω i , x i ).After that each evaluation point y ∈ Y is associated with an index of the closest stencil center ρ(y) defined as: Finally ( 7) is evaluated for every y with the matrices A and P based on the closest stencil criterion.The solution at any y ∈ Ω is then formally written as: where Ψ i are the global cardinal functions, w ρ(y) j are the local weights over the stencil centered at x ρ(y) .Furthermore the operator Γ(j, ρ(y)) : I[1, n] → I[1, N ] is an index mapping from the j − th local weight of a stencil with index ρ(y) to its global equivalent.Applying a differentiation operator L on ( 9) we obtain a representation of L in y:

Evaluation and differentiation matrices
Setting y = Y in equation ( 9) we arrive at the discrete representation of a solution in the Y -points: where E h (Y, X) is a rectangular evaluation matrix interpolating u h (X) from X to Y .Its components are: The discrete representation of a differential operator L is obtained by setting y = Y in (10): where D L h is a rectangular differentiation matrix with components (D L h ) ik = LΨ i (y k ).

The unfitted discretization of a PDE
The model problem ( 3) is discretized with the RBF-FD operators E h and D L h given in ( 11) and ( 12) respectively.The result is the semi-discrete matrix D h (y, X) and the semi-discrete vector F (y), where: Here β 2 , β 0 , β 1 are the scalings of the PDE and the boundary conditions.Setting y = Y we obtain a rectangular linear system of size M × N : For that system choose the scalings: where h is the average node distance in the node set X and where M 2 , M 0 and M 1 are the number of evaluation points placed over Ω, ∂Ω 0 and ∂Ω 1 respectively.The motivation to use this scaling is two fold.Firstly, the factors which are discrete inner products, to continuous inner products plus a first order integration error [17].Secondly, the factor h −1 is used to impose the Dirichlet condition in a weak sense such that the matrix D h is nonsingular: this is a classical approach in those finite element methods which use a solution space that does not exactly satisfy the Dirichlet condition.In order to prove uniqueness of the solution in such a setup, a mesh dependent parameter h −1 has to be introduced via inverse inequalities, which is at the end multiplying the added Dirichlet penalty term [12].We also note that we are not able to impose a Dirichlet condition exactly in an efficient way due to using X-points that are unfitted with respect to the boundary.
The numerical solution is obtained by solving (14) for the solution coefficients u h (X) and then interpolating this data onto the evaluation points Y .This can be written as: where h is a pseudoinverse in practice computed using the QR decomposition.
Once u h (X) is computed, the residual is given by: The least-squares residual is by definition orthogonal to the column-space of D h (Y, X), which implies the relation:

Linear independence of cardinal functions
In this section we address the difficulties related to the linear independence of the cardinal functions which arise when using the unfitted discretization.We provide a criterion upon which a certain amount of the interpolation points that extend outside of the computational domain is removed.A similar study, but in a context of the isogeometric finite element method is performed in [10].
The PDE matrix D h from (13) has to be constructed using a family of linearly independent cardinal functions on Ω.This is a necessary requirement for D h to have nonzero singular values.Since every column of the matrix E h from (11) contains a cardinal function in its undifferentiated form, we investigate the linear independence of those columns (sampled cardinal functions) by computing the smallest singular values σ min (E h ).When at least one singular value is zero we have that the columns have a nonzero nullspace and thus the sampled family of cardinal functions is linearly dependent.
Whereas the RBF-FD cardinal functions are indeed linearly independent when the interpolation points conform to Ω and the interpolation points X are a subset of the evaluation points Y [17], it is important to check whether this is true in the unfitted case as well.
A cardinal function Ψ k (y) has a Kronecker delta property in the X-points: which guarantees linear independence as long as X ⊆ Y and X ⊆ Ω since in this case, there is always at least one point in Ω for every cardinal function (e.g.x k for Ψ k ) where Ψ k is one, but all other Ψ j for indices k = j are 0. A problem when using the unfitted discretization can occur due to the compact support of Ψ k .When an external x k is placed such that its corresponding Ψ k vanishes before it reaches the interior of Ω, then Ψ k (y) = 0 for every y ∈ Ω and the basis function becomes linearly dependent (on Ω) with all others.
In Figure 4 we can see a one-dimensional setup, where in the top plot, the X-points are placed outside of Ω such that the left-most cardinal function Ψ left does not have a support (red ellipse) in Ω, which results in a singular E h (σ min = 0).As the support of Ψ left enters Ω, then E h becomes non-singular (σ min = 5.2 • 10 −5 ), and when the support is fully contained inside Ω then the smallest singular value gets considerably larger (σ min = 3.1).
The same setup is used in Figure 5, where for different polynomial degrees p, σ min (E h ) is computed as a function of: • the approximate area under Ψ left inside Ω, • the percentage of the compact support of Ψ left inside Ω, • the percentage of stencil points of the left-most stencil inside Ω.From Figure 5 we can see that as long as the percentage of support inside Ω is larger than 0, σ min is also larger than 0. When the percentage of support is gradually increased σ min is also increased, approximately in the same way as the area of Ψ left inside Ω.In the third plot we can see that when the percentage of the stencil points inside Ω is more than 50%, E h is always non-singular given that the cardinal functions (columns of E h ) are well sampled.We note that any cardinal function Ψ k centered in x k is not genereated by a single stencil, but by several stencils which have x k as a neighboring point.In this sense the compact support of Ψ k is decoupled from the support of one stencil centered around x k .This is the reason why the measurements of σ min do not coincide when considered as a function of a compact support inside Ω and stencil points inside Ω.
Figure 6 shows analogous results to Figure 5 but for a two-dimensional case, where the computational domain is a square of size [−1, 1] × [−1, 1].Similar results can be observed as for the one-dimensional case in terms of the area under Ψ left and the percentage of its support inside Ω, however, the percentage of stencil points inside Ω is allowed to be smaller in the 2-dimensional case.This indicates that the relation between the support of a stencil and the compact support of a cardinal function is in this case tighter.
In our experience the smallest percentage of stencil points inside Ω is the criterion which is the easiest to implement prior to computing RBF-FD differentiation matrices and it is therefore our choice for all further experiments.We remove some of the initial interpolation points (placed in the box around Ω) such that at least 50% of points of each stencil are contained inside Ω.Once we decide on a stencil size n and the initial X-points and Y -points are computed, the criterion can be used by invoking one command in Matlab.X = X(unique(knnsearch(X,Y,'k',ceil(0.5*n)),:);

Analysis of the approximation error under node refinement
In this section we develop an understanding of the behavior of the error e between the true solution u(Y ) and the approximate solution u h (Y ), both restricted to the evaluation points.

Preliminaries I: A norm for measuring the error
Our choice of a norm that measures the error e(Y ) = u(Y ) − u h (Y ) is given by: The 2 norm is a good choice for the discrete least-squares problems since it is up to O(h y ) (h y is the average spacing between the Y -points) equivalent to the L 2 norm which is a natural norm for the continuous least-squares problem [17].In this sense the stability properties can carry over to the discrete formulation, under some assumptions [17].Furthermore we have the relation e h (Y ) 2 ≤ e h (Y ) ∞ due to: which is later used to bound the consistency terms.

Preliminaries II: Consistency estimates
Evaluate a function û : Ω → R in the X-points and denote this data by û(X).Then we can, for any y ∈ Ω, where Ω contains Ω, construct a finite dimensional representation u p (y) that interpolates the data û(X): where Ψ i (y) are the RBF-FD cardinal functions defined in (9).The equivalent matrix-vector formulation is: The interpolation error is estimated by [20,17]: where h is the internodal distance of the X-points, and |û| , where D p+1 is a partial derivative of degree p + 1.An application of a PDE operator D from ( 3) to (20) gives: The (semi-discrete) matrix-vector equivalent is: There are three parts involved in D, see (3), the Laplacian, the normal derivative and the Dirichlet condition.Each of them has its own consistency estimate depending on the order of the derivative.The overall consistency [17] is bounded by:

Preliminaries III: Extension of the true solution
Generally speaking the error estimation for the unfitted RBF-FD requires a special treatment since an intermediate numerical solution u h (X), see (15), also lives in the exterior of the computational domain.In order to be able to compare the solution with the true solution we have to define an extension of the true solution on some extended open domain.Let u = u(y) be a true unique solution of the PDE problem (1), where y ∈ Ω (the closure of Ω) and let Ω = ( Ω + Ω) ⊂ R 2 be an extended domain.The extended smooth solution û = û(ŷ), where ŷ ∈ ( Ω + Ω) is then defined such that: This definition makes it possible to bound the error in terms of the stability and consistency terms, where the latter then depends on the size of the partial derivative of degree p + 1.There exists an extension lemma for smooth functions that we can use, which is stated below.The lemma from above is directly applicable to extending u : Ω → R, since Ω is a closed domain and since our u is a smooth function.We can therefore conclude that the extension defined in (25) exists, and that a discrete implication is the relation: where Y is an evaluation point set which conforms to Ω.

The PDE error estimate
Now we estimate the error e = u h (Y ) − u(Y ) between the numerical solution u h (Y ) and the true solution u(Y ): where we first added and subtracted u p (Y ), then used the triangle inequality to make a split into the PDE error u h (Y ) − u p (Y ) 2 and the interpolation error u(Y ) − u p (Y ) 2 .The latter is trivially bounded by (22).The term u h (Y ) − u p (Y ) 2 remains to be estimated.We first use the definition of u h (Y ) from (11) and the definition of u p (Y ) from ( 21) and then multiply with D + h D h = I where D h is defined in (14).16) and the fact that D + h r h (Y ) = 0 from (17) we then obtain:

Using the relation D
After that we use that by the norm relation from (18) and use the Cauchy-Schwarz inequality for the matrices D + h and E h consequtively to arrive at: It now remains to insert the estimate (28) into (27) and then also combine this with consistency estimates (22), (24) to arrive at the final error estimate: The term 1 is what we call the stability norm, which should remain constant so that the error overall decays with at least order p − 1.We numerically test that in the following section.

Detailed numerical experiments on a 2D butterfly domain
In this section we perform computational experiments to further explore the numerical properties of the unfitted RBF-FD method when solving (1) and compare them to the classical RBF-FD method in the least-squares (RBF-FD-LS) and collocation settings (RBF-FD-C).The involved parameters are (the internodal distance of the X-points), p (the polynomial degree used to form the interpolant over a stencil) and q (the oversampling parameter).More precisely, we define h as the average distance between all pairs of the neighboring interpolation points: The relation between the stencil size n and the polynomial degree p is [21]: Throughout the section we compute the relative error e as: where u h (Y ) and u(Y ) are the numerical and exact solutions sampled in the Y -points.All of the computations are performed in Matlab on a laptop with an Intel i7-7500U processor and 16 Gb of RAM.

Domain Ω
The boundary of a computational domain has a butterfly-like shape (see Figure 1) and is prescribed using a polar function: where r is a radial coordinate and θ ∈ [0, 2π] the angle.

Solution functions
We pick two solution functions to compute the right-hand-sides of (1).Those are: where u 1 is the Franke function, a commonly used infinitely smooth test function for benchmarking multivariate approximations.Function u 2 is a truncated series of an infinitely smooth function that is at the same time not analytic: we refer to u 2 as the Non-analytic function.Both functions over the butterfly domain are displayed in Figure 7.
skewed stencils are expected to significantly contribute towards a smaller PDE error.This can be accounted to smaller Lebesgue constants which enable a better approximation error.

Error as a function of runtime with Dirichlet and Neumann conditions
In the subsections above we confirmed that the approximation error is smaller for the unfitted RBF-FD-LS method compared with RBF-FD-LS and RBF-FD-C.The computational time for the unfitted RBF-FD-LS method is expected to be slightly larger than for RBF-FD-LS and RBF-FD-C due to the additional degrees of freedom that extend over the boundary of Ω.Those give rise to an increased number of columns in the matrices D h and E h from (15).The question that we address in this subsection is whether the error is small enough to compensate for the larger computational cost.In Figure 12 we can observe that in the Franke case error vs. runtime is comparable for all three settings when p = 2 and p = 4.When p = 6 the ratio is in favor of the unfitted RBF-FD method.In Figure 13 we see that in the Non-analytic case the ratio is in favor of the unfitted RBF-FD method for p = 4 and p = 6, while it is tied between the three settings for p = 2.
The conclusion is that the unfitted RBF-FD method does not lag behind in error vs. runtime compared to the two fitted settings, despite a larger amount of degrees of freedom.Moreover, in cases when stencils are large and the approximated solution oscillatory around the boundaries (the Non-analytic case) it performs significantly better.

Stability norms and condition numbers under node refinement
Here we measure the condition number of the PDE matrix D h given in (14).In addition we also measure the stability norm ||E h || 2 ||D + h || 2 which can be understood as the well-posedness constant that multiplies the consistency term in the error estimate (29).The condition number of (for example) a matrix D h is computed by: and the stability norm involved in the error estimate (29) by:  In Figure 14 we see that the stability norm of the unfitted RBF-FD-LS method is constant for p = 2.For p = 4 and p = 6 it is partially decaying until it levels out at a value slightly larger than the stability norm of fitted RBF-FD-LS.This effect is in line with the behavior of the exterior interpolation points, which move closer to the boundary as 1 h gets larger.In this way the node layout of the unfitted RBF-FD-LS method is getting increasingly similar to the node layout of the classic RBF-FD-LS method.
The condition numbers specific to the collocation, the least-squares and the unfitted least-squares setups are referred to by κ C , κ LS and κ U-LS respectively.In Figure 15 we see that when p = 2 the condition numbers κ LS and κ U-LS behave as h −2 while κ C does not follow any specific pattern (already observed in [17]).Growth with h −2 is expected, since it is known that the conditioning of any discretization involving the Laplacian operator scales with at least that rate.In the p = 4 and p = 6 cases we see a different behavior of κ U-LS compared with κ LS : here κ U-LS is at first large and remains constant until it coincides with κ LS and starts growing with approximately h −2 .This effect is conceptually very similar to the behavior of the stability norm of the unfitted RBF-FD-LS method which approaches the stability norm of RBF-FD-LS as h → 0.

Approximation properties as the polynomial degree is increased
This test gives an insight into the approximation error when the stencil size is increased together with the polynomial degree, while the internodal distance h is fixed.We choose to work with two fixed values of h, namely h = 0.05 and 0.015, where the former corresponds to a case where the approximated solution is unresolved and the latter when the approximated solution is well resolved.Results for the Franke and Non-analytic functions in the role of exact solutions are given in Figure 16 and Figure 17 respectively.In both figures the unfitted RBF-FD-LS method is superior in error to RBF-FD-LS and RBF-FD-C, especially when p is large.This is expected: when the stencil size is increasing, the stencils on the boundary get increasingly more skewed in both fitted setups, while in the unfitted RBF-FD-LS method the stencils remain fairly unskewed.See Figure 1

Experiments on a 2D drilled sprocket domain with interior boundary conditions
In this section we consider the drilled 24-tooth sprocket from Figure 2 as our computational domain.We solve (1), but in addition to the exterior mixed boundary conditions we also introduce

Point sets
The interpolation and evaluation points are -conceptually speaking -constructed in the same way as in the two-dimensional cases.We start by using a surface point-cloud of the diaphragm [23] in place of the boundary evaluation points, which are placed in a three-dimensional box that contains interpolation points computed using an algorithm from [24].Then we enforce q evaluation points (again computed by an algorithm from [24]) around each interpolation point and after that remove those evaluation points which are placed outside of the diaphragm.At last we remove the interpolation points according to the criterion given in Section 4.An instance of the resulting two point sets can be observed in the left image of Figure 22.

Convergence under node refinement
We test the convergence of the solution under node refinement for different polynomial degrees used to form the stencil approximations.The oversampling parameter is always fixed to q = 10.The spatial distribution of the magnitude of the relative error in logarithmic scale when h = 0.066  and p = 5 is given in Figure 22.The convergence results can be observed in Figure 23.We can see that the numerical solution converges for all polynomial degrees p with at least O(h p−1 ).This is an expected result according to the error estimate (29) and also according to the numerical experiments previously made for the two-dimensional cases in Section 6 and Section 7.

Final remarks
In this paper we presented the unfitted version of the RBF-FD method in the least-squares setting and in this way simplified handling of complex 2D and 3D geometries.We developed a criterion that numerically establishes the linear independence of the cardinal functions when the interpolation points are placed in the exterior of Ω.
Next, we numerically verified on a two-dimensional butterfly domain that the presented method is stable and that the stability properties closely follow the properties of the fitted RBF-FD-LS method.Moreover, the experiments confirmed that the error under node refinement decays as O(h p−1 ), which was also indicated in our theoretical work.The error was in the majority of the cases found to be smaller compared with the error RBF-FD-LS and RBF-FD-C.
Through an example of a drilled sprocket we numerically demonstrated that the less skewed boundary stencils allow the unfitted RBF-FD-LS method to outperform RBF-FD-LS and RBF-FD-C in the sense of the spatial error distribution, when the discretization is built upon large stencils.
Lastly, using a thoracic diaphragm of a human being as a computational domain, we demonstrated that the unfitted RBF-FD-LS method can also be used in a three dimensional setup, where the observed convergence trends follow O(h p−1 ), in the same way as in the two-dimensional cases.
Future work includes using the unfitted RBF-FD-LS method for other elliptic PDEs and investigating a formulation for time-dependent PDEs.

Figure 2 :
Figure2: Node distributions over a drilled 24-tooth sprocket.Left: a classical distribution of interpolation points over Ω for RBF-FD in the least-squares setting with two skewed stencils on the boundary.Here the black point is the stencil center and the red points are the members of the stencil.Right: a node distribution over Ω for the unfitted RBF-FD method in the least-squares setting with a less skewed stencil on the boundary.

Figure 3 :
Figure 3: The image on the left displays a butterfly domain with evaluation points (smaller red markers) in the Voronoi cells (grey lines) centered around every interpolation point (larger blue markers).The image on the right is a closer view over the boundary where q = 8 points are placed to every Voronoi cell.

Figure 4 :
Figure 4: Setup for measuring the smallest singular value of the evaluation matrix E h as the size of the computational domain Ω grows.The green lines represent three cardinal functions.The encircled red points represent the support of the two outmost cardinal functions.

Figure 5 :Figure 6 :
Figure 5: One-dimensional case: The smallest singular value as a function of (i) the approximate area of the outmost cardinal function Ψoutmost that penetrates inside Ω, (ii) the support of Ψoutmost inside Ω, (iii) the percentage of stencil points inside Ω, which are a part of the outmost stencil.

Lemma 1 ([ 22 ,
Lemma 4.1]).Suppose M is a smooth manifold with or without boundary, A ⊆ M is a closed subset, and f : A → R k is a smooth function.For any open subset U containing A, there exists a smooth function f : M → R k such that f | A = f and supp f ⊆ U .

u 1 u 2 Figure 7 :
Figure 7: Solution functions on a butterfly domain.Franke function over a butterfly (u 1 ), and the Non-analytic function over a butterfly (u 2 ).The black outward normals indicate the locations of the Neumann condition.The Dirichlet condition is enforced at locations where there are no normals displayed over the boundary.

Figure 13 :
Figure 13: Error versus runtime for three polynomial degrees p when the Non-analytic function is used as the exact solution.

Figure 14 :
Figure 14: Stability norm as a function of the inverse internodal distance 1/h for three different polynomial degrees p.

Figure 15 :
Figure 15: Conditioning of the PDE matrix D h as a function of the inverse internodal distance 1/h for three different polynomial degrees p.

Figure 16 :
Figure16: Error as a function of the polynomial degree p for an under-resolved case (h = 0.05) and a well-resolved case (h = 0.015) when the Franke function is chosen as the exact solution.

Figure 20 :
Figure 20: A comparison of the error distribution (logarithmic scale) over a drilled 24-tooth sprocket .The outward normals indicate the parts of ∂Ω where the Neumann condition is imposed.The internodal distance is set to h = 0.007, the oversampling parameter to q = 5 and the polynomial degree to p = 2 (first row), p = 4 (second row), and p = 6 (third row).

Figure 21 :
Figure 21: The solution function u 4 = sin(6πxyz) over a diaphragm viewed from four different angles.

Figure 22 :
Figure 22: Both pictures show a diaphragm which is sliced for visualization purposes only.Left: interpolation points (blue markers) are placed over the evaluation points (smaller red markers).The evaluation points discretize the diaphragm.Right: The spatial distribution of the magnitude of the relative error in the logarithmic scale when h = 0.066, p = 5 and q = 10.

Figure 23 :
Figure 23: Error as a function of the inverse internodal distance 1/h when the oversampling parameter is set to q = 10, for polynomial degrees p = 2, p = 4, p = 6 (left image) and polynomial degrees p = 3, p = 5 (right image).