Efficient assembly based on B-spline tailored quadrature rules for the IgA-SGBEM

TL;DR

This paper introduces B-spline tailored quadrature rules within Isogeometric Analysis for Boundary Element Methods, significantly improving computational efficiency and accuracy in assembling linear systems for 2D elliptic problems.

Contribution

It develops new weighted quadrature schemes tailored for B-splines, enhancing the assembly process in IgA-SGBEM by reducing computational cost and increasing accuracy.

Findings

New quadrature rules are exact for the test space including singular kernels.

The proposed method reduces assembly time compared to standard element-by-element approaches.

Simulations demonstrate improved accuracy and efficiency in integral evaluation.

Abstract

This paper deals with the discrete counterpart of 2D elliptic model problems rewritten in terms of Boundary Integral Equations. The study is done within the framework of Isogeometric Analysis based on B-splines. In such a context, the problem of constructing appropriate, accurate and efficient quadrature rules for the Symmetric Galerkin Boundary Element Method is here investigated. The new integration schemes, together with row assembly and sum factorization, are used to build a more efficient strategy to derive the final linear system of equations. Key ingredients are weighted quadrature rules tailored for B--splines, that are constructed to be exact in the whole test space, also with respect to the singular kernel. Several simulations are presented and discussed, showing accurate evaluation of the involved integrals and outlining the superiority of the new approach in terms of computational cost and elapsed time with respect to the standard element-by-element assembly.

Tables

Table 1. Table 1: Error E ​ R ​ R d , N h ​ [ v ​ ( t ) ] 𝐸 𝑅 subscript 𝑅 𝑑 subscript 𝑁 ℎ delimited-[] 𝑣 𝑡 ERR_{d,N_{h}}[v(t)] as defined in ( 27 ). Notice that in the case of functions B j , d ​ ( t ) subscript 𝐵 𝑗 𝑑 𝑡 B_{j,d}(t) the integrating function changes both rowwise, due to degree d 𝑑 d , and columnwise, due to N h subscript 𝑁 ℎ N_{h} .

$N_h$		$10$	$20$	$40$	$80$	$100$
$d=2$	$v\left(t\right)=B_{4,3}\left(t\right)$	$3.46\cdot {10}^{-7}$	$1.41\cdot {10}^{-7}$	$5.62\cdot {10}^{-8}$	$2.08\cdot {10}^{-8}$	$1.49\cdot {10}^{-8}$
$3$	$B_{5,4}\left(t\right)$	$2.92\cdot {10}^{-8}$	$1.40\cdot {10}^{-8}$	$7.20\cdot {10}^{-9}$	$3.75\cdot {10}^{-9}$	$3.03\cdot {10}^{-9}$
$4$	$B_{6,5}\left(t\right)$	$5.96\cdot {10}^{-10}$	$2.39\cdot {10}^{-10}$	$1.03\cdot {10}^{-10}$	$4.40\cdot {10}^{-11}$	$3.44\cdot {10}^{-11}$
$5$	$B_{7,6}\left(t\right)$	$2.58\cdot {10}^{-11}$	$9.64\cdot {10}^{-12}$	$3.92\cdot {10}^{-12}$	$1.06\cdot {10}^{-10}$	$4.67\cdot {10}^{-9}$
$d=2$	$v\left(t\right)=t^3$	$3.16\cdot {10}^{-6}$	$1.99\cdot {10}^{-7}$	$1.24\cdot {10}^{-8}$	$7.74\cdot {10}^{-10}$	$3.17\cdot {10}^{-10}$
	$t^4$	$5.21\cdot {10}^{-5}$	$3.25\cdot {10}^{-6}$	$2.04\cdot {10}^{-7}$	$1.27\cdot {10}^{-8}$	$5.22\cdot {10}^{-9}$
$d=3$	$v\left(t\right)=t^4$	$1.60\cdot {10}^{-5}$	$1.07\cdot {10}^{-6}$	$6.99\cdot {10}^{-8}$	$4.47\cdot {10}^{-9}$	$1.84\cdot {10}^{-9}$
$d=2$	$v\left(t\right)=\frac{\sqrt{1-t^2}}{t^2+25}$	$6.89\cdot {10}^{-4}$	$2.60\cdot {10}^{-4}$	$9.45\cdot {10}^{-5}$	$3.37\cdot {10}^{-5}$	$2.41\cdot {10}^{-5}$
$3$		$4.05\cdot {10}^{-4}$	$1.50\cdot {10}^{-4}$	$5.30\cdot {10}^{-5}$	$1.85\cdot {10}^{-5}$	$1.31\cdot {10}^{-5}$
$4$		$2.92\cdot {10}^{-4}$	$1.06\cdot {10}^{-4}$	$3.66\cdot {10}^{-5}$	$1.25\cdot {10}^{-5}$	$8.80\cdot {10}^{-6}$
$5$		$2.04\cdot {10}^{-4}$	$7.09\cdot {10}^{-5}$	$2.44\cdot {10}^{-5}$	$8.23\cdot {10}^{-6}$	$5.76\cdot {10}^{-6}$

Table 2. Table 2: Parabola test : comparison between the two different assembly strategy. On the top, for degree d = 2 𝑑 2 d=2 and various spacing h ℎ h . On the bottom, for constant spacing h = 1 / 5 ℎ 1 5 h=1/5 and various degrees d 𝑑 d . On the first column we report calculated spectral condition number of the matrix. On the second and third we report the calculated errors: true solution is known, thus the relative error E R subscript 𝐸 𝑅 E_{R} in L 2 superscript 𝐿 2 L^{2} norm and the absolute errors E M subscript 𝐸 𝑀 E_{M} in maximum norm are computed. Finally the elapsed time in seconds for the IgA-SGBEM matrix generation is reported. Notice that the error when d 𝑑 d increases decrees theoretically of a factor h ℎ h , and this convergence is partially maintained from the methods, being the condition number increasing. By the other side the convergence in h ℎ h is as predicted in both cases, namely of order d + 1 𝑑 1 d+1 . In all tested cases the condition number of the system is almost the same for the two assembly strategy and the overall time for assembly is strongly reduced by the new strategy.

		element-by-element				new assembly
$h$	$DoF$	$cond.$	$E_R$	$E_M$	time ($s$)	$cond.$	$E_R$	$E_M$	time ($s$)
$1/5$	$12$	$7.19\cdot {10}^1$	$1.54\cdot {10}^{-4}$	$5.08\cdot {10}^{-4}$	$6.80$	$8.56\cdot {10}^1$	$1.79\cdot {10}^{-4}$	$6.67\cdot {10}^{-4}$	$0.13$
$1/10$	$22$	$1.87\cdot {10}^2$	$1.65\cdot {10}^{-5}$	$5.73\cdot {10}^{-5}$	$15.96$	$1.87\cdot {10}^2$	$1.72\cdot {10}^{-5}$	$5.96\cdot {10}^{-5}$	$0.44$
$1/20$	$42$	$4.57\cdot {10}^2$	$1.96\cdot {10}^{-6}$	$6.92\cdot {10}^{-6}$	$41.75$	$4.50\cdot {10}^2$	$2.01\cdot {10}^{-6}$	$6.98\cdot {10}^{-6}$	$1.93$
$1/40$	$82$	$1.01\cdot {10}^3$	$2.43\cdot {10}^{-7}$	$8.50\cdot {10}^{-7}$	$134.29$	$9.90\cdot {10}^2$	$2.48\cdot {10}^{-7}$	$8.60\cdot {10}^{-7}$	$6.47$
$1/80$	$162$	$2.12\cdot {10}^3$	$3.03\cdot {10}^{-8}$	$1.07\cdot {10}^{-7}$	$566.12$	$2.07\cdot {10}^3$	$3.08\cdot {10}^{-8}$	$1.06\cdot {10}^{-7}$	$25.45$
$1/160$	$322$	$4.34\cdot {10}^3$	$3.90\cdot {10}^{-9}$	$1.97\cdot {10}^{-8}$	$2828.18$	$4.25\cdot {10}^3$	$3.85\cdot {10}^{-9}$	$1.34\cdot {10}^{-8}$	$98.87$
$1/320$	$642$	$8.78\cdot {10}^3$	$7.17\cdot {10}^{-10}$	$1.82\cdot {10}^{-8}$	$16659.45$	$8.59\cdot {10}^3$	$4.81\cdot {10}^{-10}$	$1.64\cdot {10}^{-9}$	$409.79$
$d$	$DoF$	$cond.$	$E_R$	$E_M$	time ($s$)	$cond.$	$E_R$	$E_M$	time ($s$)
$2$	$12$	$7.19\cdot {10}^1$	$1.54\cdot {10}^{-4}$	$5.08\cdot {10}^{-4}$	$6.80$	$8.56\cdot {10}^1$	$1.79\cdot {10}^{-4}$	$6.67\cdot {10}^{-4}$	$0.13$
$3$	$13$	$1.89\cdot {10}^2$	$3.23\cdot {10}^{-5}$	$1.69\cdot {10}^{-4}$	$11.93$	$2.14\cdot {10}^2$	$5.63\cdot {10}^{-5}$	$3.87\cdot {10}^{-4}$	$0.21$
$4$	$14$	$5.09\cdot {10}^2$	$1.80\cdot {10}^{-5}$	$8.01\cdot {10}^{-5}$	$20.44$	$5.61\cdot {10}^2$	$2.19\cdot {10}^{-5}$	$1.20\cdot {10}^{-4}$	$0.25$
$5$	$15$	$1.41\cdot {10}^3$	$5.43\cdot {10}^{-6}$	$2.06\cdot {10}^{-5}$	$34.27$	$1.65\cdot {10}^3$	$1.05\cdot {10}^{-5}$	$5.53\cdot {10}^{-5}$	$0.31$

Table 3. Table 3: S-shaped domain test: comparison between the two different assembly strategies for degree d = 3 𝑑 3 d=3 . For the explanation of the columns see the previous table. Notice that in all the tested cases the new assembly is accurate and efficient.

		element-by-element				new assembly
$h$	$DoF$	$cond.$	$E_R$	$E_M$	time ($s$)	$cond.$	$E_R$	$E_M$	time ($s$)
$1/6$	$15$	$3.39\cdot {10}^2$	$1.12\cdot {10}^{-1}$	$1.62\cdot {10}^{-1}$	$98.28$	$2.87\cdot {10}^2$	$1.14\cdot {10}^{-1}$	$2.33\cdot {10}^{-1}$	$2.17$
$1/12$	$27$	$8.20\cdot {10}^2$	$3.23\cdot {10}^{-2}$	$6.71\cdot {10}^{-2}$	$351.41$	$9.19\cdot {10}^2$	$3.30\cdot {10}^{-2}$	$9.39\cdot {10}^{-2}$	$5.97$
$1/24$	$51$	$3.00\cdot {10}^3$	$3.98\cdot {10}^{-3}$	$1.33\cdot {10}^{-2}$	$1240.44$	$3.67\cdot {10}^3$	$4.03\cdot {10}^{-3}$	$1.97\cdot {10}^{-2}$	$22.50$
$1/48$	$99$	$8.75\cdot {10}^3$	$5.79\cdot {10}^{-4}$	$2.03\cdot {10}^{-3}$	$4728.49$	$8.95\cdot {10}^3$	$6.32\cdot {10}^{-4}$	$3.11\cdot {10}^{-3}$	$88.80$
$1/96$	$195$	$2.09\cdot {10}^4$	$8.96\cdot {10}^{-5}$	$3.06\cdot {10}^{-4}$	$18615.21$	$2.15\cdot {10}^4$	$1.85\cdot {10}^{-4}$	$1.69\cdot {10}^{-3}$	$372.18$

Equations

$$\left\{\begin{array}[]{ll}\Delta u=0&{\rm in}\;\Omega\,,\\ u=u_{D}&{\rm on}\;\Gamma\,,\end{array}\right.$$

$$u(\mathbf{x})=-\frac{1}{2\pi }{\int }_{\Gamma }\hat{K}(\mathbf{x},\mathbf{y})q(\mathbf{y})d{\gamma }_{\mathbf{y}}+\frac{1}{2\pi }{\int }_{\Gamma }\frac{\partial \hat{K}}{\partial {\mathbf{n}}_y}(\mathbf{x},\mathbf{y})u_{\mathcal{D}}(\mathbf{y})d{\gamma }_{\mathbf{y}},\mathbf{x}\in \Omega$$

$$\hat{K}(\mathbf{x},\mathbf{y}):=ln(r),\text{mbox}withr:=\Vert \mathbf{x}-\mathbf{y}{\Vert }_2,$$

$$-\frac{1}{2\pi }{\int }_{\Gamma }\hat{K}(\mathbf{x},\mathbf{y})q(\mathbf{y})d{\gamma }_{\mathbf{y}}=\frac{1}{2}{u}_{\mathcal{D}}(\mathbf{x})-\frac{1}{2\pi }{\int }_{\Gamma }\frac{\partial \hat{K}}{\partial {\mathbf{n}}_y}(\mathbf{x},\mathbf{y})u_{\mathcal{D}}(\mathbf{y})d{\gamma }_{\mathbf{y}},\mathbf{x}\in \Gamma ,$$

$$\mathcal{A}(q,p)=\mathcal{F}(p),\forall p\in H^{-1/2}(\Gamma )$$

$$\mathcal{A}(q,p):=-\frac{1}{2\pi }{\int }_{\Gamma }p(\mathbf{x}){\int }_{\Gamma }\hat{K}(\mathbf{x},\mathbf{y})q(\mathbf{y})d{\gamma }_{\mathbf{y}}d{\gamma }_{\mathbf{x}}$$

$${\cal F}(p):=\int_{\Gamma}p({\bf x})\,\Big{[}\frac{1}{2}\,u_{\cal D}({\bf x})-\frac{1}{2\pi}\,\int_{\Gamma}\frac{\partial\hat{K}}{\partial{\bf n}_{y}}({\bf x},{\bf y})\,u_{\cal D}({\bf y})\,d\gamma_{\bf y}\Big{]}\,d\gamma_{\bf x}$$

$$\left\{\begin{array}[]{ll}\Delta u=0&{\rm in}\;{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{2}\setminus\Gamma,\\ u=u_{D}&{\rm on}\;\Gamma\,.\end{array}\right.$$

$$-\frac{1}{2\pi }{\int }_{\Gamma }\hat{K}(\mathbf{x},\mathbf{y})\phi (\mathbf{y})d{\gamma }_{\mathbf{y}}=u_{\mathcal{D}}(\mathbf{x}),\mathbf{x}\in \Gamma ,$$

$$\mathcal{F}(\psi ):={\int }_{\Gamma }\psi (\mathbf{x})u_{\mathcal{D}}(\mathbf{x})d{\gamma }_{\mathbf{x}}.$$

$$\mathbf{f}(s):=\sum _{i=1}^N{\mathbf{Q}}_i{B}_{i,d}(s),s\in I,$$

$$\hat{S}:=span\{{\hat{B}}_1,\cdots ,{\hat{B}}_N\},$$

$${\hat{B}}_i(\mathbf{x}):=B_i({\mathbf{f}}^{-1}(\mathbf{x})),\mathbf{x}\in \Gamma ,i=1,\dots ,N.$$

$$\hat{S}:=span\{{\hat{B}}_1,\cdots ,{\hat{B}}_{N-d}\},$$

$${\hat{B}}_i(\mathbf{x}):=B_i^{(c)}({\mathbf{f}}^{-1}(\mathbf{x})),\mathbf{x}\in \Gamma ,i=1,\dots ,d,$$

$$B_{i}^{(c)}(t):=\left\{\begin{array}[]{ll}B_{i}(t)&\mbox{if }t\in I\cap[t_{i},t_{i+d+1}]\,,\cr B_{N-d+i}(t)&\mbox{if }t\in I\cap[t_{N-d+i},t_{N+i+1}]\,,\cr 0&\mbox{otherwise\,.}\end{array}\right.$$

$$A\,\mbox{\boldmath$\alpha$}=\mbox{\boldmath$\beta$}\,,$$

$${\int }_{\Gamma }{\hat{B}}_i(\mathbf{x}){\int }_{\Gamma }\hat{K}(\mathbf{x},\mathbf{y}){\hat{B}}_j(\mathbf{y})d{\gamma }_{\mathbf{y}}d{\gamma }_{\mathbf{x}},$$

$${\int }_{\Gamma }{u}_D(\mathbf{x}){\hat{B}}_i(\mathbf{x})d{\gamma }_{\mathbf{x}},{\int }_{\Gamma }{\hat{B}}_i(\mathbf{x}){\int }_{\Gamma }\frac{\partial \hat{K}}{\partial {\mathbf{n}}_y}(\mathbf{x},\mathbf{y})u_D(\mathbf{y})d{\gamma }_{\mathbf{y}}d{\gamma }_{\mathbf{x}}.$$

$$s:={\mathbf{f}}^{-1}(\mathbf{x}),t:={\mathbf{f}}^{-1}(\mathbf{y})$$

$$J(\cdot ):=\sqrt{(f_1^{\prime }(\cdot ){)}^2+(f_2^{\prime }(\cdot ){)}^2},$$

$$I_K^{(i,j)}:={\int }_I{B}_i(s)J(s){\int }_{I}K(s,t)B_j(t)J(t)dtds,$$

$$\hat{K}(\mathbf{f}(s),\mathbf{f}(t))=\ln (\Vert \mathbf{f}(s)-\mathbf{f}(t){\Vert }_2^2{)}^{1/2}=\frac{1}{2}ln\left(\frac{\Vert \mathbf{f}(s)-\mathbf{f}(t){\Vert }_2^2}{(s-t{)}^2}\right)+ln\left|s-t\right|$$

$$R(s,t):=\frac{\Vert \mathbf{f}(s)-\mathbf{f}(t){\Vert }_2^2}{(s-t{)}^2}={\left[\frac{f_1(s)-f_1(t)}{s-t}\right]}^2+{\left[\frac{f_2(s)-f_2(t)}{s-t}\right]}^2,$$

$$K(s,t)=K_1(s,t)+K_2(s,t),$$

$$K_1(s,t):=\frac{1}{2}ln\left(R(s,t)\right),K_2(s,t):=ln|s-t|.$$

$$\underset{t\to s}{\lim }R(s,t)=J^2(s),$$

$$I_K^{(i,j)}=I_{K_1}^{(i,j)}+I_{K_2}^{(i,j)}$$

$$I_{K_{\ell }}^{(i,j)}:={\int }_I{B}_i(s)J(s){\int }_I{K}_{\ell }(s,t)B_j(t)J(t)dtds,\ell =1,2.$$

$$b_1^{(i)}:={\int }_I{B}_i(s)J(s)u_D(s)ds,$$

Full text

1 Introduction

Boundary Element Methods (BEMs) are an important strategy for the numerical solution of linear partial differential equations appearing in many relevant problems in science and engineering applications, see [6, 17]. Through the fundamental solution associated to the considered differential equation, a large class of both exterior and interior elliptic Boundary Value Problems (BVPs) can be reformulated as linear Boundary Integral Equations (BIEs), reducing by one the dimension of the computational domain for the discretization, when compared with Finite Element Methods (FEMs). For all these problems BEMs can be adopted, offering substantial computational advantages over other numerical techniques. However, in order to achieve a general efficient numerical implementation, a number of issues has to be carefully addressed. One of the most important consists in the accurate approximation of weakly singular, Cauchy singular and even hyper-singular integrals over the boundary. Such integrals occur also when the Symmetric Galerkin Boundary Element Method (SGBEM) [12] is applied to some BIEs. Indeed, in this case a linear symmetric system of equations is obtained, whose coefficient matrix has entries defined as possibly singular double integrals on the boundary of the problem’s domain. Thus, considering the outlined difficulties, it is not surprising that, since the first appearances of SGBEM [20, 37, 48, 56], great effort has been devoted to the efficient and accurate computation of the related Galerkin linear system, as proved by many papers investigating this aspect (see [1, 2] and references therein). Such techniques were combined to the SGBEM scheme based on Lagrangian bases, always adopting the standard element-by-element assembly procedure, as customary when the Lagrangian basis is used. This implies that all the double integrals were split into a sum of integrals on pairs of elements where a local quadrature rule was adopted.

The advent of Isogeometric Analysis (IgA), [10, 21, 34] has brought a renewed interest in BEMs and very recently it has been combined for the first time to the SGBEM scheme [3, 4, 41]. Papers dealing with problems in acoustics [19], airfoils potential flows [38], Stokes flows [32], fluid-structure-interactions [33] can be found in the literature, using the IgA Galerkin BEM formulation. BEM formulations have been used also to construct computational domains for Galerkin-IgA [23]. In order to reduce complexity and gain efficiency, the use of collocation [51] or mixed collocation is also common, see e.g. [28, 30, 35, 36, 39, 40, 42, 43, 45, 46, 47].

Following the IgA approach, in [3] B-splines have been used to represent both the domain geometry and the approximated solution of the problem at hand, giving a significant reduction in the dimension of the discretization space required to attain a fixed accuracy with respect to the standard Lagrangian basis. In such preliminary implementation of the IgA-SGBEM scheme, the entries of the coefficient matrix and of the right-hand side of the linear system are evaluated by using the accurate quadrature rules introduced in [1] and suited for an element-by-element implementation of the SGBEM scheme. However, such quadrature rules are not capable to take any computational advantage from the higher regularity of B-splines. The obvious consequence is that the method cannot benefit further from its new isogeometric formulation. In order to avoid this drawback, in the present paper new quadrature rules tailored on B-splines are introduced. The aim is to compute efficiently as well as accurately all the occurring double, possibly singular, integrals, and this is a key point for a new, efficient assembly strategy. The considered quadrature rules are of two kinds. When no singularity appears a B–spline weighted quadrature is adopted [15, 16]. This is constructed for each B-spline basis function of the approximation space, taken as weight function, in order to be exact in a suitable B-spline space, usually a refinement of the previous one. In presence of singularities, a new idea is applied. Following the classic approach based on modified moments [15, 22, 27], and using the recurrence relations for B-splines [13], we propose a new procedure so that, also in this case, the final quadrature is exact in the test space or in a related refinement. As done in [5, 16] for the IgA-FEM case, we assembly the final matrix via row loop and sum-factorization. This allows us to get a remarkable gain in terms of computational costs, as confirmed by the experiments. The proposed strategy is ready to be used for problems equipped by Dirichlet data. For mixed problems, hyper-singular integrals can appear which require some more analysis, for example using coordinate transformation [44, 52], subtracting the singularities [25, 29, 31] or splitting the contributions in order to recover quantities that can be evaluated in closed forms [2].

The paper is organized as follows. In Section 2 the boundary integral formulation of two model problems is briefly introduced. In Section 3 some preliminaries and notations are given, which are necessary in the isogeometric setting for the representation of the domain geometry and of the discretization space; then the double integrals defining the entries of the coefficient matrix and of the right-hand side of the IgA-SGBEM linear system are defined. In Section 4, using a unified formulation, two families of quadrature rules are introduced to approximate the double integrals introduced before. The families are two, since different rules are needed to deal with regular or singular integrals, both occurring in any BEM formulation. The formulas used to deal with non singular integrals are taken from [16], and they are briefly summarized in Subsection 4.1. Novel formulas addressing the double singular integrals on fixed nodes and tailored on B-splines are introduced in Subsection 4.2, where a recursive relation for the computation of the necessary modified moments for B-spline functions is presented. These last quadrature rules are also preliminarily tested on some singular integrals whose exact solution is available. Section 5 deals with the efficient computation and assembly of the IgA-SGBEM linear system. The new B-spline–oriented IgA-SGBEM assembly strategy is tested through some numerical examples presented and discussed in Section 6. Finally, Section 7 reports the research conclusions.

2 Boundary integral model problems

In this work we focus on 2D interior or exterior Laplace model problems on planar domains, assuming boundary Cauchy data of Dirichlet type. In particular we deal with two different geometries: bounded simply connected domains $\Omega\subset\mathbf{R^{2}}$, and unbounded domains external to an open limited arc.

In the first case, denoting with $\Gamma$ the boundary of $\Omega$, assumed sufficiently regular, we deal with the boundary value problem

[TABLE]

where $u_{D}$ is the given boundary datum.

Choosing a direct approach [17], the boundary integral reformulation of (1) starts from the representation formula for the solution $u$, i.e.

[TABLE]

where ${\bf x}=(x_{1},x_{2})$, ${\bf y}=(y_{1},y_{2})$, the kernel $-\frac{1}{2\pi}\hat{K}({\bf x},{\bf y})$ represents the fundamental solution of the 2D Laplace operator, i.e.

[TABLE]

and $q:=\frac{\partial u}{\partial{\bf n}}$. Hence, the solution $u$ can be evaluated in any point of the domain, provided that we know the flux $q$ on $\Gamma$. So, with a limiting process for ${\bf x}$ tending to $\Gamma$ and using the boundary datum, we obtain the BIE:

[TABLE]

in the boundary unknown $q$. The boundary integral problem (3) can be finally set in the following weak form [55]:

given $u_{D}\in H^{1/2}(\Gamma)$, find $q\in H^{-1/2}(\Gamma)$ such that

[TABLE]

where

[TABLE]

and

[TABLE]

In the second case, still denoting with $\Gamma$ an open limited arc in the plane, we consider the following problem

[TABLE]

The BVP in (6) can model the electrostatic problem of finding the electric potential around a condenser, whose two faces are so near one another to be considered as overlapped, knowing the electric potential only on the condenser, see e.g. [9]. This example is a classic case where BEM are preferred to FEM: a problem in an infinite domain with an obstacle or a source lying on a curve [20, 50].

Choosing an indirect approach [17, 50], the BIE coming from the boundary integral reformulation of (6) reads:

[TABLE]

where the unknown density function $\varphi$ represents the jump of $q$ across $\Gamma$.

Boundary integral problems (7) can be set in a weak form similar to (4), where the bilinear form ${\cal A}(\varphi,\psi)$ is defined as in (5) and the right-hand side simplifies in:

[TABLE]

The price for the simplification of the indirect approach is that the discrete solution obtained by solving the linear system does not approximate directly the missing Cauchy data on the boundary, but just the density function appearing in the chosen integral representation formula, see e.g. [17].

3 IgA-SGBEM discretization

Now, we assume that in both the above model problems, $\Gamma$ is defined as a parametric curve with no self intersections, closed in the first case and open in the other. Thus $\Gamma$ is the image of a regular invertible function ${\bf f}:I\subset{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}\rightarrow\Gamma\subset{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{2}$ such that, setting ${\bf f}(s):=(f_{1}(s),f_{2}(s))\,,$ every point ${\bf x}=(x_{1},x_{2})\in\Gamma$ can be seen as the image of just one value $s\in I\,,$ where $I=[a\,,\,b]$ if $\Gamma$ is open and $I=[a\,,\,b),$ otherwise. In more detail, we assume that $\Gamma$ is parametrically represented by a function ${\bf f}$ defined as follows,

[TABLE]

where the ${\bf Q}_{i},\,i=1,\ldots,N\,,$ are ordered control points assigned in the plane and defining the shape of $\Gamma$, and where $\{B_{i,d}(\cdot),\,i=1,\ldots,N\}$ is the B-spline basis which spans a space $S$ of piecewise $d$–degree polynomials with respect to an assigned set $\Delta$ of distinct breakpoints in $I.$ Note that $N$ and the regularity required in $S$ at each inner breakpoint can be a priori established. Actually the definition of the B-spline basis needs the preliminary introduction of the extended knot vector $T=\{t_{1},\,\ldots,t_{d},t_{d+1},\ldots,t_{N+1},t_{N+2}\,,\ldots,t_{N+d+1}\},$ where $t_{d+1}=a\,,\,t_{N+1}=b\,,$ and $t_{1}\leq\cdots\leq t_{N+d+1}\,.$ Each knot $t_{i},\ i={d+2},\ldots,N,$ is a possibly repeated occurrence of an inner breakpoint, while the first and last $d$ knots in $T$ are auxiliary knots characterizing a specific B–spline basis. For brevity, when it is not strictly necessary, in the following we shall refer to it by omitting $d$. The regularity at a certain inner breakpoint in $\Delta$ of any function in $S$ is fully specified by prescribing the breakpoint multiplicity (integer between $1$ and $d+1$) in $T$. Since we want ${\bf f}$ at least with a continuous image, such multiplicities are required to be all $\leq d$ and one of them is set to $d$ only when $\Gamma$ has an angular point which will coincide with a control point. Extended knot vectors with multiple or periodic auxiliary knots are the more commonly adopted strategies for completing the definition of $T$, see e.g. [13]. Note that the in Computer Aided Geometric Design (CAGD) the standard way to represent closed curves relies on periodic extended knot vectors and with the last $d$ control points given by an ordered repetition of the first ones, see e.g. [26]. This clearly means that in practice the dimension of the considered spline space $S$ is $N-d$ in this case.

In the IgA context, dealing with an isoparametric approach, the same basis or an its appropriate refinement is considered to generate the discretization space $\hat{S}$ where we will search the approximate solution of weak BIEs. Actually spaces with dimension greater than $N$ (or $N-d$ in the closed case) can be clearly used for the experiments, without abandoning the IgA paradigm, since the knot insertion procedure can always be adopted to represent the boundary in a higher dimension spline space with a desired mesh spacing or with a reduced regularity at the breakpoints. We remind that such a procedure implies introducing new breakpoints or increasing the multiplicity of the existing ones, see e.g. [24].

In more detail, assuming for notational simplicity that no refinement is done, when $\Gamma$ is an open arc, we difine the space $\hat{S}$ as

[TABLE]

where

[TABLE]

When $\Gamma$ is a closed curve, the dimension of $\hat{S}$ must be equal to $N-d.$ Thus we fix $\hat{S}$ as

[TABLE]

where $\hat{B}_{i},\,i=d+1,\ldots,N-d$ are defined as in (8) while we have

[TABLE]

with

[TABLE]

In order to simplify the notation, in the sequel, when $\Gamma$ is closed, we shall omit the superscript $(c)$ to denote the first $d$ cyclic basis elements.

The algebraic reformulation of the IgA-SGBEM scheme leads to a linear system of equations, whose unknowns represents the coefficient of the BIE approximate solution w.r.t. the chosen basis, see e.g. [44]. In more detail, denoting such a solution with $\hat{\alpha}({\bf x}):=\sum_{j=1}^{N_{DoF}}\alpha_{j}\hat{B}_{j}({\bf x})\,,\,\mbox{with }N_{DoF}:=dim(\hat{S})\,,$ the resulting linear system has size $N_{DoF}$ and is referred as

[TABLE]

where $A$ is a symmetric positive definite matrix, \mbox{\boldmath\alpha}=(\alpha_{1},\cdots,\alpha_{N_{DoF}})^{T} is the unknown vector and \mbox{\boldmath\beta}\in{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{N_{DoF}} is the vector defining the right-hand side which depends on the given Cauchy data.

In particular, the entries of the matrix $A$ in (9) are, up to the coefficient $-\frac{1}{2\pi}$, double integral of the type (see (5)):

[TABLE]

while the right-hand side involves the computation, depending on the problem at hand, of the following integrals

[TABLE]

Introducing the scalar variables $s$ and $t$ defined as:

[TABLE]

and the parametric speed associated with $\Gamma,$

[TABLE]

in following subsections we rewrite integrals in (10), (11) over the parametric interval $I$.

3.1 Matrix entries

Referring to the double integral in (10), we can express it as

[TABLE]

where $K(s,t):=\hat{K}({\bf f}(s),{\bf f}(t))\,$. Let us investigate the nature of the involved kernel. Considering (2), we can write

[TABLE]

Setting

[TABLE]

we have that

[TABLE]

where

[TABLE]

Note that $K_{1}(s,t)$ can be defined also for $s=t$ extending its definition by continuity, since

[TABLE]

with $J$ defined as in (12).

This splitting of $K$ is useful to separate two contributions, the first coming from the geometry of $\Gamma$ and the last depending just from its singular nature. Actually, integral (13) can be evaluated as

[TABLE]

with

[TABLE]

Different quadrature rules are needed to compute $I_{K_{\ell}}$, $\ell=1,2$, since for $\ell=2$ weakly singular integrals can occur because of the logarithmic nature of the kernel $K_{2}$.

3.2 Right-hand side elements

Referring to the first integral in (11) we can rewrite it as:

[TABLE]

where $u_{D}(s)=u_{D}({\bf f}(s))\,$. Hence we have to evaluate a regular integral.

The second double integral in (11) can be expressed as

[TABLE]

where $u_{D}(t)=u_{D}({\bf f}(t))\,$ and

[TABLE]

We note that, when ${\bf f}\in C^{2}(I),$ no singularity occurs in (18). Indeed, with some computation, we can write

[TABLE]

that is,

[TABLE]

with $R(s,t)$ defined in $(\ref{Rst})$. Then, taking into account (15), we have

[TABLE]

On the other hand, without assuming ${\bf f}\in C^{2}(I)\,,$ the kernel $\bar{K}$ becomes strongly singular on boundaries with corners and weakly singular on Lyapunov curves, see [8, Section 7].

4 Novel quadrature rules

As pointed out in the introduction, in the present paper we explore the construction of the algebraic counterpart of the IgA-SGBEM scheme, by using two different weighted quadrature strategies for regular and singular integrals. This means that the usual conditions considered for the construction of the quadrature - the exactness requirements - are imposed with respect to a suitable chosen weight function.

The developed formulas have in common the vector of nodes - denoted by $\eta$ - but the weights will change according to the exactness requirements, consisting in imposing their exactness in a suitable spline space with the respect to the selected weight function. Several quadrature rules are necessary, all determined by solving a linear system whose coefficient matrix is in any case a B-splines collocation matrix. The existence of the rules is ensured by the non singularity of such a matrix, feature achieved by choosing an $\eta$ vector fulfilling the Schoenberg-Whitney conditions, see for example [13].

The selected exactness spline space is a possible refinement of the test space $S$, obtained by uniform subdivision of each element of $\Delta$ into $N_{ref}\geq 1$ elements, where $N_{ref}>1$ is adopted to improve the quadrature accuracy. This - possibly new - partition of the interval $I$ is denoted by $\bar{\Delta}$. The corresponding spline space can be generated by B-splines, denoted in the following by $\bar{B}_{j}$, $j=1,\dots,N_{E}$, with $N_{E}\approx N_{ref}N_{DoF}$.

Our choice of $\eta$, that fulfills the Schoenberg-Whitney conditions with respect to the extended knot vector associated with $\bar{B}_{j},j=1,\ldots,N_{E},$ is the following: in the first and last element of $\bar{\Delta}$, we take $d+2$ uniformly spaced points; on each inner element we take midpoints and breakpoints. Hence \mbox{\boldmath\eta}\in{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{N_{quad}} with $N_{quad}=2d+2N_{ref}N_{h}-1,$ where $N_{h}$ denotes the number of elements of $\Delta$ ($N_{h}\approx N_{DoF}$ when no multiple inner knot is included in the extended knot vector $T$). Note that this simple choice was used also in [16] and for BEM with splines in [7].

We detail the construction of the quadrature rules in the following two subsections, respectively related to regular and singular integrals.

4.1 B-spline weighted quadrature rule

When dealing with the external integrals in (16) and (18), with the integrals (17), or even with the inner integral in (16) for the kernel $K_{1}$, regular integrands occur, always with a B-spline factor, denoted in the sequel as $B_{i}(s)$. In this case, following [16], we consider this term as a weight and thus we obtain the quadrature:

[TABLE]

Let us detail the construction introducing some notation111Notice that we define support as the open set where the function is non-zero..

$\mathcal{J}^{(i)}:=\{j$ : $1\leq j\leq N_{E}$, $\textit{supp}(\bar{B}_{j})\cap\textit{supp}({B}_{i})\neq\emptyset\}$ ;  $N^{(i)}_{E}:=\#\mathcal{J}^{(i)}$ (number of active exactness functions);

• 2.

$\mathcal{N}^{(i)}:=\{n$: $1\leq n\leq N_{quad}$, with $\eta_{n}\in\textit{supp}(B_{i})\}$;  $N^{(i)}_{quad}:=\#\mathcal{N}^{(i)}$ (number of active quadrature nodes).

The weights $w_{n}^{(i)}$, $n\in\mathcal{N}^{(i)}$, are determined by imposing the exactness of the formula on all $\bar{B}_{j}$ with $j\in\mathcal{J}^{(i)}$, the other weights are set to 0. The local exactness requirements then read:

[TABLE]

where $\mu^{(i)}_{j}:=\int_{I}\bar{B}_{j}(s){B}_{i}(s)\,ds$ can be exactly computed since the integrand is a piecewise polynomial function. Let us note that the considered definition of $\eta$ always ensures that $N^{(i)}_{E}\leq N^{(i)}_{quad}$. Thus, the conditions in (20) lead to a possibly underdetermined linear system with maximum rank, thanks to the Withney-Schoenberg conditions fulfilled by $\eta$. Actually the system is non squared only when the support of $B_{i}$ includes (but is not limited to) the first or the last element of $\bar{\Delta}$, where more nodes are taken. In this case we have chosen to solve the system in the minimum Euclidean norm. Note that $N^{(i)}_{E}$ can be upper bounded independently from $N_{DoF}$, since $N^{(i)}_{E}\leq(1+N_{ref})\ (d+1)\,$. Consequently the size of the linear system in (20) never becomes prohibitive.

This construction has to be repeated for all the basis functions $B_{i}\,,\ i=1,\dots N_{DoF}$. A related Pseudo–code is given in Algorithm 1.

4.2 Singular weighted quadrature rule

When instead we deal with the inner integrals in (16), for $K=K_{2}$, weakly singular integrals can occur, as discussed in Section 3, hence a different formula is used. We choose to isolate the singular term, as done in [49] and introduce a quadrature:

[TABLE]

In order to fix the weight vector $\mathbf{w}^{s}:=(w_{n}^{s})_{n=1}^{N_{quad}}$, we impose exactness on the B-spline functions $\bar{B}_{j}$, $j=1,\ldots,N_{E}$, that is we require the fulfillment of the following conditions:

[TABLE]

with

[TABLE]

Then, $\mathbf{w}^{s}$ is the solution of a linear system of size $N_{E}\times N_{quad}$ with $N_{E}\leq N_{quad}$ whose coefficient matrix is a B-splines collocation matrix at the quadrature nodes. We remark that this is a -possibly underdetermined- linear system with maximum rank, thanks to the Withney-Schoenberg conditions. The right hand side \mbox{\boldmath\mu}^{s}:=(\mu_{j}(s))_{j=1,\dots,N_{E}}, of (22) is calculated by exploiting the recursive formula for B-splines, see below. Recalling that we are dealing with the inner integrals in (16), for $K=K_{2}$, and that the formulas introduced in the previous subsection are used for the outer integrals, the parameter $s$ varies among the entries of the vector $\eta$. So different weight vectors ${\bf w}^{\eta_{n}}\in{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{N_{quad}},\,n=1,\ldots,N_{quad}$, are pre-computed.

Let us derive now the recursive approach which can be use to compute the analytic expressions of the modified moments introduced in (23). Given a B-spline basis of degree $r$, defined with respect to the partition $\bar{\Delta}$ of $I$, we set

[TABLE]

Then we can write

[TABLE]

since $B_{j,d}=\bar{B}_{j}$.

Taking into account the Cox-De Boor recurrence relation of B-splines [13], assuming that $t_{j}$ is the $j$–th entry of the associated extended knot vector,

[TABLE]

we have that a consequent recurrence relation can be obtained for $I_{q}(B_{j,r},s)$:

[TABLE]

We remark that, when multiple knots are taken, if $t_{j+r}=t_{j}\;(t_{j+r+1}=t_{j+1}),$ the first (second) addend in the right-hand side of (25) must be set to zero. In order to compute the desired modified moments (24) by means of (25), we need to compute preliminarily $I_{k}(B_{i,0},s)$, for $k=0,\ldots,d$ and $i=1,\ldots,N_{E}+d$. Now, considering that

[TABLE]

we obtain that $I_{k}(B_{i,0},s)=\int_{t_{i}}^{t_{i+1}}ln(|t-s|)\,t^{k}\,dt\,,$ which clearly is a vanishing quantity if $t_{i}=t_{i+1}.$ In the opposite case, using the substitution $z=t-s$, we can write

[TABLE]

Thus, the procedure starts by computing the quantities:

[TABLE]

Note that with some basic analytic computation we can derive the following explicit expression of ${\cal I}_{i}^{j}(s)$ for the non trivial case:

[TABLE]

The described procedure for the computation of the singular modified moments can be seen as a variant of that introduced in [22], where it was defined for Cauchy singular integrals and used for developing quadrature formulas based on cubic spline interpolation.

A related Pseudo–code sketching the main steps required for the computation of ${\bf w}^{s}$ is given in Algorithm 2. Notice that this procedure can be used -for example in collocation BEM methods- with respect to a generic vector of points $\sigma$, while in our case we will have \mbox{\boldmath\sigma}=\mbox{\boldmath\eta}.

Let us remark that we are proposing a new weighted quadrature where the weight is the singular kernel $\ln|t-s|$ for a fixed value of $s$, and exactness is required on B-splines of a fixed degree. Following the analysis in [54], this requirement is the main ingredient for optimal convergence of the BIE approximate solution, because the quadrature error estimate implies consistency of the overall scheme. Then, it is important to have some indication of the quadrature error when applying the quadrature to a generic function. We consider the approximation of the integral

[TABLE]

for functions $v(t)$ such that these values can be evaluated analytically and for $s=\eta_{n}$, where $\eta_{n}$ is one of the quadrature points introduced in Section 4. Then, the value ${\mathrm{I}}^{s}[v]$ is approximated by ${\cal Q}^{s}[v]$ as in equation (21). The error is calculated as follows:

[TABLE]

where $N_{quad}$ is the number of quadrature points of the external quadrature. The complete analysis of the error of such rules is beyond the scope of the present paper, we only report here some tests that we have run in order to confirm their good behavior. We have done three kind of tests:

First, we have checked that the exactness requirements are fulfilled. This issue is mainly concerned with the condition number of the collocation matrix: we have chosen a-priori fixed quadrature points that could give bad-conditioned matrices. In all our tests we get machine precision in the numerical computation of the modified moments with the quadrature formula (21). Moreover, we have checked that the integral of the monomial (considered in the whole patch) where exactness is due -namely the monomial $t^{d}$- is computed properly. Also in these cases all tests up to degree $d=5$ give exact result up to machine precision. 2. 2.

Then, we have considered the computation of moments where exactness is not required: for quadrature rules where exactness is required on B-splines of degree $d$, we have computed approximated integrals of B-splines of degree $d+1$ and of the monomial $t^{d+1}$. Results are reported in Table 1. 3. 3.

At last, we have considered the integral:

[TABLE]

taken from reference [53]. This is an interesting benchmark, and in the cited paper the integral is approximated by a procedure that computes an auxiliary integral via symbolic computations. On such function our aim is to test monotone convergence confirmed by results available in Table 1.

5 Efficient computation and assembly of IgA-SGBEM linear system

Applying quadrature rules introduced in the previous Section, and referring at first to (16), we can compute

[TABLE]

and

[TABLE]

where the last equivalence is due to the local support of $B_{j}\,,$ so that the number of non vanishing addends does not depend on $N_{quad}$ and therefor neither on $N_{DoF}.$

For the numerical evaluation of $I_{K_{2}}^{(i,j)}$, we could distinguish the case $|i-j|>d+1$ where the singularity does not occur, so that the regular quadrature could be used as in (29). Nevertheless, in our implementation, we have used quadrature rule (30) also for regular integrals in the construction of matrix $I_{K_{2}}$, because we have checked that this change of rule not only gives no advantage in terms of computational time, but also gives the same accuracy.

For the numerical evaluation of (17), we can simply consider

[TABLE]

At last, since we assume ${\bf f}\in C^{2}(I),$ for (18) we proceed by computing

[TABLE]

where \left(\mbox{\boldmath\eta}^{GQ},{\bf w}^{GQ}\right) is a gaussian quadrature, being the integrand function regular.

5.1 Assembly: computational cost

Following what described in the previous sections, the proposed procedure for the assembly of the algebraic counterpart of the IgA-SGBEM scheme can be sketched as follows:

Fix $N_{ref}$ and a choice of $N_{quad}$ quadrature points $\eta$. 2. 2.

Construct the family of quadrature rules $\left(\eta^{(i)}_{n},w^{(i)}_{n}\right)_{n=1,\dots N^{(i)}_{quad}},i=1,\dots,N_{DoF}$ following the procedure in Algorithm 1. This is done by solving $N_{DoF}$ linear systems of dimension $N_{E}^{(i)}\leq(1+N_{ref})(d+1)$. 3. 3.

Construct the family of quadrature rules $\left(\eta_{n},w^{\eta_{m}}_{n}\right)_{n=1,\dots N_{quad}},m=1,\dots,N_{quad}$ following the procedure in Algorithm 2. This is done by solving $N_{quad}$ linear systems of dimension $N_{E}\times N_{quad}\approx N_{ref}N_{DoF}\times 2N_{ref}N_{DoF}$ where only the right hand size changes, so that the matrix can be factorized just once. We have performed in our tests the usual $LU$ factorization as implemented in Matlab. 4. 4.

Assembly the IgA-SGBEM matrix and right-hand side of interest. For computational efficiency, we propose row assembly via sum factorization. Details of such construction can be found in [5, 16]. In Algorithm 3 we present the technique in the case of (29).

Due to the savings introduced by the patch-strategy for the construction of the weighted quadratures, the number of quadrature points for each parametric variable is $2d+2N_{ref}N_{h}-1$. Thus the computational cost for the assembly of the IgA-SGBEM linear system matrix with the introduced B-spline based quadrature strategy and sum factorization is $O([2d+2N_{ref}N_{h}-1]^{2})$ function evaluations.

Remark. In [3] the matrix entries were numerically evaluated within the framework of the standard element-by-element assembly phase. This implies that, within a double nested cycle over the $N_{h}$ mesh elements of the partition $\Delta$ on the parametrization interval $I$, having indicated by $e^{(i)}_{\ell_{i}},\,\ell_{i}=1,\cdots,d+1$ the elements constituting the support of the B-spline $B_{i}$, we computed and assembled (13) as

[TABLE]

taking into account the polynomial nature of the B-splines over each element of their support and using quadrature schemes introduced in [2], suitable for standard local Lagrangian basis functions. In particular, in presence of kernel singularity, arising for double integration over couples $(e^{(i)}_{\ell_{i}},\,e^{(j)}_{\ell_{j}})$ of coincident or consecutive elements, related double integrals in (31)

were split into the sum of a regular part and a singular part: the first was treated using Gauss-Legendre rule for both inner and outer integrals, while the other term was evaluated using an interpolatory product rule, absorbing kernel singularity into the weights, for the inner integration and the Double Exponential rule (briefly, DE-rule), suitable for integrand function having weak singularities at the endpoints, for the outer integration. When the double integration occurred on couples of ”far” elements Gauss-Legendre rule for both inner and outer integrals was employed. The interested reader is referred to [2] and related references for more details on the above resumed quadrature schemes.

Hence, having denoted by $N_{G}$ the number of the Gauss-Legendre rule points, by $N_{DE}$ the number of DE-rule points and by $N_{prod}$ the number of points of product rules for singular integration, the computational cost of the element-by-element procedure for the generation of the IgA-SGBEM linear system matrix, in terms of integrand function evaluation, was $(d+1)^{2}(N_{h}^{2}\,N_{G}^{2}+3\,N_{h}\,N_{DE}\,N_{prod})\,.$

6 Numerical examples

In this section we present and discuss some numerical results obtained using both the old quadrature strategy with the element-by-element assembly and the new B-spline based quadrature with the related new assembly. The old code is developed in Fortran language while the new one is written in Matlab. Simulations have been performed on a laptop equipped by Intel Core i5 CPU (2.53 GHz, 4 Gb RAM, 64 bit OS).

In all the following experiments, the error $E_{M}$ represents the maximum of the discrete error on a uniform mesh of $500$ points in the parametric domain $I$, while $E_{R}$ is defined as the relative error in $L^{2}(I)$ norm. The details on the considered BVP are reported in the description of each numerical test.

6.1 Exterior Dirichlet problem: parabola test

At first, let us consider the Dirichlet BVP defined in (6) exterior to the arc of parabola $\Gamma=\{{\bf x}=(x_{1},x_{2})\,|\,x_{1}=t,\,x_{2}=1-t^{2},\,t\in[-1,1]$}, representable by means of quadratic B-splines related to the extended knot vector

[TABLE]

and to the control points ${\bf Q}_{i},\,i=0,\,1,\,2$, whose coordinates are collected in the following matrix:

[TABLE]

Here, the Dirichlet datum $u_{\cal D}$ is given in such a way that the solution of (7) is explicitly known and reads $\varphi({\bf x})=\sqrt{1+4x_{1}^{2}}$. The discrete counterpart of such BIE involves the computation of matrices $I_{K_{1}}$ and $I_{K_{2}}$ with, respectively, regular and weakly singular integrals, see equations (29)-(30). Also, the right-hand side is regular so that in our computations with the new assembly a B-spline weighted quadrature rule has been used. The comparison reported in Table 2, for different values of the parameter $h$, which uniformly decomposes the parameter interval $[-1,1]$, involves $C^{1}$ quadratic B-spline basis used in the old (element-by-element, with $N_{G}=N_{prod}=32,\,N_{DE}=63$) and new ($N_{ref}=1$ for both quadrature formulas) implementations of IgA-SGBEM.

Together with DoF and spectral condition numbers of the associated matrices, we show the relative error $E_{R}$ in $L^{2}$ norm, the absolute errors $E_{M}$ in maximum norm and the elapsed time for the IgA-SGBEM matrix generation. Note that the slightly different results concerning conditioning and accuracy are due to the different quadrature schemes adopted in the old and new assembly strategies. The elapsed time behaves in both cases as ${\cal O}(N_{h}^{2})$ as expected, but the superiority of the new approach is self-evident.

Similar conclusions can be deduced looking at the results in the bottom part of Table 2 where the comparison has been done for increasing values of the B-spline basis degree $d$ chosen for the approximation of the BIE solution. Here the elapsed time behaves in both cases as ${\cal O}(d^{2})$.

The same data are used in Figure 1 where the convergence of errors are plotted with respect to the number of degrees of freedom and with respect to assembly time. In the first plot it can be seen that the solution obtained with the new strategy is accurate as the one calculated before, in the second plot it can be seen that the new strategy obtains very good errors using modest times.

6.2 Interior Dirichlet problem: S-shaped closed domain test

We consider the interior Laplace problem (1) on the domain $\Omega$ shown in Figure 2, equipped by Dirichlet boundary condition, where $\Gamma:=\partial\Omega$ is described by cubic B-splines defined by the cyclic extended knot vector

[TABLE]

and control points depicted in Figure 2. Since $u_{\cal D}=-(x_{1}+x_{2})$ is chosen, the solution of (3) is explicitly known, it reads $q({\bf x})=q({\bf f}(t))=(f^{\prime}_{1}(t)-f^{\prime}_{2}(t))/\|{\bf f}^{\prime}(t)\|_{2}$ and has $C^{1}([-1,1])$ regularity.

Table 3 compares results obtained for this numerical test by the IgA-SGBEM element-by-element implementation ($N_{G}=N_{prod}=32,\,N_{DE}=63$) and the new one ($N_{ref}=2$ for both the quadrature formulas). Also in this case, the new quadrature and assembly strategy reveals much faster than the old one. Because of modest regularity of the solution the performed tests are limited to the case $d=3$. The plot of the approximate solution obtained from the new implementation with $h=1/48$ is depicted in Figure 2, where the analytical solution of (3) is also reported.

This example is more challenging from the numerical point of view, due to modest regularity, the use of a direct approach and the oscillations of the solution, see Figure 2. Moreover, we point out that the S-shaped domain, being not starred, is not good for some of the interior point-BEM methods introduced recently [18].

7 Conclusions and future work

In this paper we have presented a new strategy that gives fast assembly in IgA-SGBEM dealing with BVPs equipped by Dirichlet boundary conditions. The remaining analysis for the application to mixed problems will be covered in a forthcoming paper.

Note that the presented integration schemes can be easily used also in the context of collocation BEM, fixing the outer collocation node and using the new quadrature rules for the inner integration.

Moreover, the proposed method can be profitably used when treating non-linear transmission problems, in order to obtain the trace of the solution on the boundary, see [14].

An open issue consists in testing hierarchical spline spaces to develop an efficient adaptive isogeometric version of the scheme, see [11]. The hierarchical approach could overcome the problem of dealing with non-regular curves.

The extension to 3D problems where the boundary surface can be described just by one open patch and the multi-patch case for the description of closed surfaces is currently under study.

Acknowledgements

The support by Gruppo Nazionale per il Calcolo Scientifico (GNCS) of the Istituto Nazionale di Alta Matematica (INdAM) through “Progetti di ricerca” program is gratefully acknowledged. Giancarlo Sangalli was partially supported by the European Research Council through the FP7 ERC Consolidator Grant No.616563 ”HIGEOM”.