A collocation IGA-BEM for 3D potential problems on unbounded domains

The framework of Isogeometric Analysis (IGA), introduced in 2005 by Thomas J. R. Hughes and coauthors, has emerged as a very attractive approach to solve differential problems LibroHughes ; Hughes_2005 . IGA bridges the gap between Computational Mechanics and Computer Aided Design (CAD) by employing a unified representation of geometric objects and physical quantities. In particular, the spline basis functions commonly used for CAD geometry parameterization are also used to span the discretization space adopted for the numerical solution. This not only eliminates the need for slow and error-prone conversion processes between the geometry and simulation models, but also enables the possibility of a widespread use of design optimization tools in simulation-based engineering.

While IGA has achieved great success within the finite element paradigm by attracting most of the IGA researchers, a remaining critical challenge is how to efficiently and exactly parameterize a volumetric domain from its boundary al2016Brep ; cohen_etal2010 ; pan2020volumetric . Three dimensional problems need trivariate NURBS solids for analysis, but CAD systems use a boundary representation where only bivariate NURBS are involved. A possible solution to this issue consists in adopting boundary element methods (BEMs) beer2020LIBRO ; BEMbook , whose isogeometric formulation is briefly referred to as IGA-BEM and has recently been succesfully applied to several classical problems, e.g., Laplace and Helmholtz equation TauRodHug ; Dolz18 , elastostatics and elasticity Simp2012 ; Wang15 ; beer_etal2017 ; nguyen2017 , potential flow gong2017 and wave-resistance problems ginnis . BEMs reformulate the original differential problem through the fundamental solution into boundary integral equations. Thus, adopting a boundary representation, they can be easily coupled with standard CAD techniques and ensure dimension reduction of the computational domain where the numerical simulation is developed. Besides these key points, they are attractive because offering an easy treatment of problems on unbounded domains, see for example the potential or acoustic problems considered in Coox ; Nguyen16 . Within this kind of problems, screen or crack problems are particularly important for applications, being in this case the domain given by the whole Euclidean space external to an open limited arc in 2D or surface in 3D.

Relying on the general NURBS parametric description of the domain boundary and on (refinements of) the related spline space for the approximation of the unknowns, in this paper we confirm the effectiveness of the IGA-BEM in the 3D setting, previously experimented in 2D ADSS1 . In particular, we focus on a collocation discretization strategy for 3D exterior potential problems. The combination of the superior smoothness of spline basis functions with the low computational cost and simplicity of collocation techniques seems to constitute an optimal basis for accurately modeling complex and computationally demanding problems, at least when smooth geometries are dealt with, as also outlined in TauRodHug . However it is worth mentioning that considering a Galerkin based variant of the developed method could be of great interest in the future, either because the method could become more robust for treating problems with low regular solutions sutradhar2008SGBEM or even from a more theoretical point of view, since most of the convergence results available in the literature refer to the Galerkin approach.
Since an effective implementation of BEMs surely requires suitable rules for the approximation of arising singular integrals, in this paper we focus on such aspect. In particular, the cubature rules for bivariate weakly singular integrals developed in Nash20 are here firstly applied to the 3D IGA-BEM context. These rules are based on spline Quasi-Interpolation (QI) and on the usage of a spline product formula Morken91 . Specifically formulated for weakly singular integrals containing a tensor-product B-spline factor in the integrand, their peculiarity consists in being applicable on its whole support. Since tensor product B-splines are here adopted to span the discretization space, this feature implies that the assembly phase of the system matrix can be done function-by function. This is a clear advantage of our approach, as already outlined in the 2D setting ACDS3 . Furthermore, as a preliminary step the multiplicative singularity extraction procedure considered in Nash20 is here improved with the aim to apply the QI operator adopted by our rules to sufficiently smooth functions. Note that in the IGA-BEM setting, according to the distance between the source point and the integration domain, not all the integrals to be dealt with are singular, some being near singular and others regular. The final accuracy of our IGA-BEM implementation gained from treating near singular integrals with the same strategy adopted for singular ones. For regular integrals, a simplified formulation of the mentioned cubature rules has been developed, relying on the same QI operator, but possibly selecting for it a different spline degree and/or knot spacing.

The structure of the paper is as follows. Section 2 recalls the indirect integral formulation of a potential problem on a 3D unbounded domain, and then it introduces the related discretization based on an isogeometric collocation BEM. Section 3 presents cubature rules based on quasi-interpolation and on the spline product formula. Section 4 reports the results obtained with a Matlab implementation of the scheme for a screen problem and for two different problems on the unbounded domain external to a torus. Finally, some conclusions are given in Section 5.

As well as for interior problems, even for general exterior problems different integral formulations are available in the literature, coming from different simple or double layer representation formulas. BEMs derived from an integral equation which directly involves the unknown Cauchy datum are called (simple or double layer) “direct” methods. As an alternative, (simple or double layer) “indirect” approaches can also be adopted. In such case the unknown appearing in the integral equation and in the associated representation formula is a different function, usually referred to as “density” function. In this paper we adopt a BEM indirect formulation, since it is the only one applicable also to screen/crack problems, costabel1986principles . In particular we focus on such formulation for 3D Dirichlet exterior potential problems,

where $g$ is the given Dirichlet datum and $\Gamma:=\partial\Omega,$ with $\Omega$ denoting a bounded domain in ${\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{3}$. Note that for screen problems $\Omega$ has to be replaced with $\Gamma,$ where $\Gamma$ is a given open and limited surface. Observe also that, since we deal with unbounded domains, it is necessary to specify the required behavior of the harmonic solution at infinity to ensure uniqueness BEMbook . The boundary integral formulation of (1) is given by the following BIE,

where the unknown density function $\psi$ represents the jump across the boundary of the normal derivative of $u$ (in the outward direction from ${\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{3}\setminus\overline{\Omega}$ to $\Omega$) and the kernel $U$ is the fundamental solution of the Laplace equation in 3D,

The operator $V$ defined in (2) is an elliptic isomorphism between its functional domain $H^{-\frac{1}{2}}(\Gamma)$ and its codomain $H^{\frac{1}{2}}(\Gamma),$ while the BIE also defined in (2) is a Fredholm integral equation of the first kind with weakly singular kernel.

Once the density function $\psi$ has been (approximately) computed, the (approximated) expression of $u$ at any point in ${\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{3}\setminus\overline{\Omega}$ can be derived by relying on its simple layer representation formula,

Note that this formula automatically ensures the required decay behavior of $u$ at infinity if $\psi$ is continuous in $\Gamma$ BEMbook , and also that its numerical approximation can be challenging only when ${\bf x}$ is very close to $\Gamma.$

To compute a numerical solution of the BIE in (2) we adopt collocation, since this is probably the easiest discretization method suited for smooth problems (i.e., problems with a smooth $\Gamma$ and sufficiently smooth Dirichlet datum). However, we observe that the related theory is not yet fully developed in particular for the 3D case; for example in Sections 9.2 and 9.3 of Atkinson_2009 the available theoretical results can be applied just for the direct approach which is not usable for screen problems, as we already remarked. In general, after having numerically solved the considered BIE with a BEM approach, the representation formula associated with it has to be considered in order to compute values of the approximated solution inside the domain. However there are several applications, e.g. the simulation for screen or crack problems, which do not require such post-processing step, since their physically relevant quantity is the unknown Cauchy datum on the boundary itself, see for example costabel1986principles .

We assume that the geometry is CAD-generated, that is $\overline{\Gamma}$ has a parametric representation, $\overline{\Gamma}=\mbox{Image}({\bf F})\,,$ where ${\bf F}:R\rightarrow{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{3}$ is a one-to-one mapping¹¹1With the exception of the boundary of the parametric domain in case of a closed $\Gamma$ or with cylinder-like shape. defined on the parametric rectangular domain $R:=[a_{1}\,,\,b_{1}]\times[a_{2}\,,\,b_{2}]$ in the following standard NURBS form,

where ${\bf t}:=(t_{1},t_{2})$. Each set of functions $B^{d_{j}}_{i,T_{j}},i=0,\ldots,n_{j},j=1,2$ is made up of univariate B-spline functions defined in $[a_{j}\,,\,b_{j}],$ of degree $d_{j}$ and with extended knot vector $T_{j}$. We assume that the mapping ${\bf F}$ is at least differentiable. The set $\{{\bf Q}_{i,j}\in{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{3},i=0,\ldots,n_{1},j=0,\ldots,n_{2}\}$ defines the net of control points²²2In order to avoid degenerate nets we assume the control net composed by distinct points. and each $w_{i,j}$ is a positive weight used to obtain additional shape control (obviously, if the weights are all equal to a common constant value, ${\bf F}$ becomes just a vector tensor-product polynomial spline represented in B-spline form). This surface representation is often adopted in CAD environments to design free-form surfaces and it has the facility of involving just one patch, that is a unique continuous vector function ${\bf F}$ is used to define the whole considered surface. On this concern we note that clearly this representation can be not flexible enough for complex geometries characterized for example by linkage curves with just geometric continuity. Note also that, even if the NURBS form was originally introduced for guaranteeing exact representation of conic sections and quadric patches, it cannot be used for example for representing closed surfaces with zero genus without singularities in the parameterization. For example in reference TauRodHug the one-patch NURBS parameterization of the sphere is singular at the poles and this has required the development of special cubature formulas to obtain accurate results. In order to deal with surface representations more flexible and more suitable for the analysis, multi-patch NURBS representation could be considered, as well as different kinds of spline spaces, e.g., multi-degree polar splines TSHH20 .

Focusing on the basic standard single patch NURBS representation and referring for brevity only to the case of clamped knot vectors $T_{j},j=1,2$, let us define ${\cal T}_{j}$ as a possible refinement of $T_{j},j=1,2,$ obtained by inserting new breakpoints in $[a_{j}\,,\,b_{j}],$ or by increasing the multiplicity of a knot already included in $T_{j}.$ Denoting with $h_{j}$ the maximal distance between successive knots in ${\cal T}_{j},j=1,2$ and with $S_{h_{j}}$ the univariate spline space of dimension $m_{j}+1\,,m_{j}\geq n_{j},$ spanned by the B-splines of degree $d_{j}$ associated with ${\cal T}_{j},$ we can define the tensor-product spline space ${\cal S}^{h}:=S_{h_{1}}\otimes S_{h_{2}}$ of dimension $D:=(m_{1}+1)(m_{2}+1),$ where $h:=\max\{h_{1},\,h_{2}\}.$ Note that in the following we simplify the notation by denoting as $B_{j}$ the bivariate B-spline $B^{d_{1}}_{i_{1},{\cal T}_{1}}B^{d_{2}}_{i_{2},{\cal T}_{2}}$ generating ${\cal S}^{h},$ where $j=i_{2}(m_{1}+1)+i_{1}.$ Thus, the density function $\psi$ is approximated in a finite dimensional space $\hat{\cal S}^{h}$ derived from ${\cal S}^{h}$ by lifting all the B-splines to the physical space

Using collocation, the coefficient vector $\mbox{\boldmath$\alpha$}_{h}=(\alpha_{0},\cdots,\alpha_{D-1})^{T}\in{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{D}$ associated with the IGA-BEM discrete solution $\psi_{h}\in\hat{\cal S}^{h},\ \psi_{h}({\bf x})=\sum_{j=0}^{D-1}\alpha_{j}(B_{j}\circ{\bf F}^{-1})({\bf x})\,,$ is determined by solving the following linear system with $D$ equations,

where

and $({\bf g}_{h})_{i}:=g({\bf x}_{i})\,.$

The points ${\bf x}_{i},i=0,\ldots,D-1,$ are $D$ distinct collocation points belonging to $\Gamma$ defined as the image through ${\bf F}$ of the Cartesian product of two sets of distinct abscissas, defined in $[a_{1},\,b_{1}]$ and in $[a_{2},\,b_{2}],$ respectively, as for example Greville abscissas (or suitable variants, see for example Wang15 , to avoid collocation on the boundary of $\Gamma$). Thus, for each collocation point ${\bf x}_{i}\in\Gamma,i=1,\ldots,D,$ there exists ${\bf s}_{i}\in R$ such that ${\bf x}_{i}={\bf F}({\bf s}_{i}),$ with ${\bf s}_{i}\neq{\bf s}_{j},i\neq j.$

In order to be able to compute the discrete solution efficiently and accurately, it is important to have an efficient and stable method for solving the linear system in (6) but even more to consider suitable cubature rules for weakly singular integrals to be used in the assembly phase. In this paper we focus on this second issue, extending to cubature the quadrature rules previously developed for analogous univariate weakly singular integrals. As already mentioned in the introduction, the use of such rules (which are specific for integrals containing an explicit B-spline factor in the integrand) allows us to avoid the element-by-element assembly strategy, since they are directly applied on the whole support of each $B_{j}.$

In this section we present the cubature rules introduced in Nash20 and here firstly applied in the 3D IGA-BEM context for the numerical approximation of all the integrals on the right hand side of (8). Such rules rely on spline quasi-interpolation and are an extension to the bivariate case of those firstly studied in the univariate setting CFSS18 and successfully tested also working in adaptive spline spaces FGKSS_2019 .

We first observe that the integrals defining the entries of $A_{h}$ are weakly singular only if the source point ${\bf s}_{i}$ belongs to ${\overline{R}}_{j}.$ If this is not the case, they are nearly singular or regular, depending on the distance between ${\bf s}_{i}$ and the integration domain $R_{j},$ clearly taking also into account if $\Gamma$ is a closed surface. This initial classification is important because nearly singular integrals need as much care as weakly singular ones. In our implementation we adopt for them the same strategy used for weakly singular integrals, while a simplified version of the developed rules is used to approximate regular integrals, possibly using also higher QI spline degree in such a case.

Referring first of all to singular and nearly singular integrals, we observe that, as well as in the 2D case ACDS3 ; FGKSS_2019 , a preliminary singularity extraction procedure is necessary in order to increase the performance of our cubature rule. In Nash20 an analysis has been developed to show that two singularity extraction procedures, based on a Taylor expansion of ${\bf F}({\bf t})$ about ${\bf s}$, are possible, one using an additive splitting of $U$ and the other of multiplicative type. Considering that the first splitting requires the numerical computation of two integrals instead of one, here we have chosen to rely on the multiplicative technique. In this case the idea is to rewrite the entries of the matrix $A_{h}$ given in (8) in the following equivalent form,

where

and $K_{\bf F}$ is a suitably locally defined weakly singular kernel such that the function $\rho_{{\bf s},{\bf F}}$ is at least continuous. In particular in Nash20 we considered the following definition of such reference kernel,

with $P_{{\bf s},{\bf F}}$ denoting the bivariate homogeneous quadratic polynomial vanishing with its first derivatives at ${\bf t}={\bf s}$ and obtained by replacing ${\bf F}({\bf t})$ with its first degree Taylor expansion about ${\bf s}$ in $\|{\bf F}({\bf t})-{\bf F}({\bf s})\|_{2}^{2}$. Note anyway that for a general surface this definition just ensures continuity to $\rho_{{\bf s},{\bf F}}$ at ${\bf t}={\bf s}$. However, if ${\bf F}$ is sufficiently smooth at the singular point ${\bf s}$, it is possible to use higher order expansions for ${\bf F}$, in order to ensure higher regularity to $\rho_{{\bf s},{\bf F}}$. In particular here we consider an expansion with up to $3$ terms, which ensures up to $C^{2}$ smoothness to $\rho_{{\bf s},{\bf F}}$. Details about this improved singularity extraction procedure will be presented in a forthcoming paper³³3T. Kanduč. Isoparametric singularity extraction technique for 3D potential problems in BEM, arXiv:2203.11538., since they form a technical part necessary for an exhaustive explanation but too long for being here reported. We just mention that in this case $K_{\bf F}$ is defined with the following expression which reduces to (10) for $n=1,$

where $P_{3\ell-3}^{({\bf s},{\bf F})}$ are appropriate homogeneous polynomials of total degree $3\ell-3\,,\ell=1,\ldots,n.$ Concerning this improved singularity extraction procedure, finally we add that, when $n>1$ and the computational mesh is not fine enough, a correction term $\eta\|{\bf s}-{\bf t}\|_{2}^{2}$ for some suitably small $\eta>0$ has to be added on the right hand side of (11) to avoid the introduction of new singularities in $\rho_{{\bf s},{\bf F}}$.

After this preliminary rewriting of the entries of the matrix $A_{h},$ each of them is approximated as follows,

where $\Delta_{j}$ is a uniform partition of $R_{j}$ and $Q_{{\bf p},\Delta_{j}}$ denotes a tensor-product QI operator which approximates a bivariate function with a tensor-product spline of bi-degree ${\bf p}=(p_{1},p_{2})$ with respect to the partition $\Delta_{j}.$

In order to compute the integral on the right hand side of (12), the spline $Q_{{\bf p},\Delta_{j}}\left(J\ \rho_{{\bf s}_{i},{\bf F}}\right)$ is multiplied by the B-spline $B_{j},$ using for this aim the extension to the tensor-product setting of the Mørken formula Morken91 for spline product. As a result, the whole product $\left(J\ B_{j}\ \rho_{{\bf s}_{i},{\bf F}}\right)$ is approximated with a local tensor-product spline function defined in $R_{j},$ having bi-degree ${\bf p}+{\bf d},$ and breakpoints in $R_{j}.$ We observe that, in order to implement the so derived cubature rules, the evaluation of the following modified moments is firstly necessary,

with $B$ denoting any tensor-product B-spline. This computation is done by relying on the recursion formula for B-splines and on the preliminary analytic computation of element integrals whose integrand is the product of a monomial times $K_{\bf F},$ procedure possible for a kernel $K_{\bf F}$ defined as in (10) or more generally as in (11).

When the integral defining the $(i,j)$–th entry of $A_{h}$ is classified as a regular one, we do not need the preliminary computation of the modified moments and the strategy is simplified by directly setting

where ${\bf q}=(q_{1}\,,\,q_{2})$ denotes a bi-degree for the local splines possibly different from ${\bf p}$ and $\Lambda_{j}$ is a uniform partition of $R_{j}$, also possibly different from $\Delta_{j}.$

Concerning the QI spline operator used by the proposed cubature rules, for our experiments we have always considered the uniform tensor-product formulation of the derivative free variant of the Hermite QI operator introduced in MSbit09 , see for example BGMS16 for an introduction to its formulation in the bivariate setting and for an application in the context of adaptive bivariate spline approximation. This univariate spline QI operator has approximation order $p+1$ and in the context of numerical integration it has special interest because its application to both regular and weakly singular quadrature is superconvergent when even spline degrees are considered on uniform knot distributions MSJcam12 ; CFSS18 . Thus for smooth functions $Q_{{\bf p},\Delta_{j}}$ has approximation order $\min\{p_{1},p_{2}\}+1$ and it inherits the superconvergence feature when analogously applied for cubature. Observe however that, in our approach, choosing values of $p_{k},k=1,2,$ greater than $2$ is not useful, because of the limited regularity of the function $\rho_{{\bf s}_{i},{\bf F}}.$ It can be profitable instead to select $q_{k}>2,k=1,2$, if ${\bf F}$ is smooth in $R_{j}.$ Considering the assumed analytic expression of ${\bf F},$ this is surely true when the interior of $R_{j}$ does not contain any straight line $t_{j}=\tau_{j},$ with $\tau_{j}$ geometric knot in the $j$–th parametric direction, that is $\tau_{j}\in T_{j}\,,j=1,2,$ where $T_{j},j=1,2$ are the extended knot vectors used to define ${\bf F}.$

In this section we test the uniform formulation of the proposed IGA boundary element model for the numerical solution of two Laplace problems on two different 3D unbounded domains. Furthermore a nonuniform implementation of the scheme is also tested for solving a more difficult Laplace problem exterior to the first considered domain. In all the examples the improved singularity extraction procedure is used within the developed singular cubature rules, always setting the parameter $\eta$ (see Section 3) equal to $1.$ Concerning the discretization, in all the presented experiments approximate densities are constructed on a given initial mesh and also on successive finer ones obtained from a dyadic $h$-refinement procedure which inserts a simple knot at the midpoint between every pair of consecutive breakpoints in each parametric direction.
For the first two considered examples the density $\psi$ is known, while the Dirichlet datum $g$ is not available in closed form. Thus, the entries of the right-hand side in (6) are computed numerically,

by applying the same approach adopted to compute the system matrix. For these examples the accuracy of the scheme is measured by the $L^{2}(\Gamma)$ norm of the absolute density error $(\psi_{h}-\psi),$ in particular considering its dependency on the mesh size $h.$ The third exterior Laplace problem is based on Test $5.3$ of Dolz18 where the more general Helmholtz equation is solved. Since in this case the analytic expression of the exact potential $u$ is known, the Dirichlet datum is given and the accuracy in the potential approximation $u_{h}\approx u$ can be measured, provided that the representation formula given in (4) is previously applied to $\psi_{h}$ to define $u_{h}$ in the considered unbounded domain ${\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{3}\setminus\overline{\Omega}.$

.

In this paper a collocation IGA-BEM scheme relying on the indirect integral formulation of Laplace problems is proposed and tested for 3D exterior problems, obtaining the expected convergence order in the density approximation. Cubature rules based on spline quasi-interpolation and on a spline product formula are applied for the numerical approximation of both weakly and nearly singular bivariate integrals appearing in the 3D IGA-BEM scheme. They are combined with a robust multiplicative singularity extraction procedure to increase the accuracy of the rules. A variant of such rules is proposed for the regular integrals to be also dealt with in the assembly phase. This can be efficiently implemented in a function-by-function form thanks to the peculiarity of the considered cubature.