IGA-Energetic BEM: An Effective Tool for the Numerical Solution of Wave Propagation Problems in Space-Time Domain

: The Energetic Boundary Element Method (BEM) is a recent discretization technique for the numerical solution of wave propagation problems, inside bounded domains or outside bounded obstacles. The differential model problem is converted into a Boundary Integral Equation (BIE) in the time domain, which is then written into an energy-dependent weak form successively discretized by a Galerkin-type approach. Taking into account the space-time model problem of 2D soft-scattering of acoustic waves by obstacles described by open arcs by B-spline (or NURBS) parametrizations, the aim of this paper is to introduce the powerful Isogeometric Analysis (IGA) approach into Energetic BEM for what concerns discretization in space variables. The same computational beneﬁts already observed for IGA-BEM in the case of elliptic (ı.e., static) problems, is emphasized here because it is gained at every step of the time-marching procedure. Numerical issues for an efﬁcient integration of weakly singular kernels, related to the fundamental solution of the wave operator and dependent on the propagation wavefront, will be described. Effective numerical results will be given and discussed, showing, from a numerical point of view, convergence and accuracy of the proposed method, as well as the superiority of IGA-Energetic BEM compared to the standard version of the method, which employs classical Lagrangian basis functions.


Introduction
The Isogeometric Analysis (IGA) paradigm, introduced by Hughes and collaborators in the seminal paper [1] in the context of the Finite Element Method (FEM), has become more and more popular in the Boundary Element Method (BEM) community [2]. In fact, from one side, BEM needs to work only on the boundary of the problem domain, hence avoiding domain meshing, and, from the other, IGA uses the same B-spline or NURBS basis, as adopted in CAD systems, to describe both the boundary and the approximate solution of the problem at hand, drastically reducing the number of degrees of freedom (DoFs) and giving surprising computational advantages with regard to the classical Lagrangian basis. In the last 10 years, IGA-BEM has been fruitfully applied to a variety of elliptic problems reformulated in terms of Boundary Integral Equations (BIEs), such as the Laplace equation [3][4][5][6][7][8][9][10][11][12][13], Helmholtz equation [9,[13][14][15][16], elastostatics, and crack problems [17][18][19][20][21]; on the other hand, little attention has been devoted to the numerical solution of differential problems in time domains (see [22][23][24] as very recent works in this context).
In particular, for what concerns hyperbolic problems, the analysis of the propagation or the scattering of acoustic or elastic waves is extremely important in various sciences and application fields, from engineering to geology, and from seismic risk assessment to the design of electronic devices [25][26][27][28][29][30]. Moreover, avoiding the frequency domain, simulations in the time domain allow to study the phenomenon as it evolves.
For this kind of problem, an accurate and stable space-time BEM, namely Energetic BEM, has been recently introduced by the first author and her collaborators, and applied to the wave equation inside bounded domains or outside bounded obstacles, with or without damping terms, proving its robustness, accuracy and long-time stability, both from theoretical and numerical points of view (see e.g., [31]). The differential initial-boundary value problem is converted into a space-time BIE, which is then written into an energydependent weak form successively discretized by a Galerkin-type approach.
This work takes into account the space-time model problem of 2D soft-scattering of acoustic waves by obstacles described by open arcs with B-spline (or NURBS) parametrizations. It aims at introducing the IGA paradigm into energetic BEM, for what concerns the discretization in space variables, showing the same benefits already observed in the context of IGA-BEM for elliptic problems [5], especially the superiority of the IGA approach versus the standard Energetic BEM. The latter is based on the classical local Lagrangian basis for the approximation of the problem solution. The computational gain already noted for elliptic problems is more emphasized here, because it is obtained at every step of the time-marching procedure.
The paper is structured as follows: at first, the model problem and its energetic boundary integral weak formulation will be introduced, then the IGA-Energetic BEM will be outlined. Successively, some important issues involved in the numerical evaluation of the entries of the final discretization linear system matrix will be taken into account. In fact, these entries are double space-time integrals involving the weakly singular fundamental solution of the 2D wave operator, which further depends on the propagation wavefront: hence, the related integration domain changes at every time-step and is reduced until the wavefront completely covers it. At last, several effective numerical results will be presented and discussed, showing, from a numerical point of view, convergence and accuracy of IGA-Energetic BEM, as well as its optimal computational performance in terms of DoFs. A conclusive section ends the paper.

Model Problem on Curved Boundary and Its Energetic Boundary Integral Weak Formulation
Let us consider the following differential problem for the wave equation, exterior to a smooth and open arc Γ ⊂ R 2 , equipped by trivial initial conditions and the Dirichlet boundary datum, which models a soft-scattering phenomenon: u(x, 0) = u t (x, 0) = 0, u(x, t) = g(x, t), where c is the propagation velocity of a perturbation inside the domain, T is the final time instant of analysis, and the boundary datum g(x, t) represents the time history of the excitation field over Γ. Let us suppose that Γ is defined as a parametric curve with no self-intersection. Thus, Γ is the image of a regular invertible function C : [a, b] ⊂ R −→ Γ ⊂ R 2 , such that setting C(s) := (C 1 (s), C 2 (s)), every point x = (x 1 , x 2 ) ∈ Γ can be seen as the image of just one value s ∈ [a, b].
In order to rewrite problem (1)-(3) as a BIE, we need to start from the single-layer representation formula for the solutions of the Equation (1) [32]: u(x, t) = Γ t 0 G(r, t − τ)ψ(y, τ) dτ dγ y , x ∈ R 2 \ Γ, t ∈ (0, T], Mathematics 2022, 10, 334 3 of 30 written in terms of the forward fundamental solution of the 2D wave operator where r = x − y 2 and H(·) is the Heaviside step-function, and in terms of the unknown density function ψ(x, t), which represents the time history of the jump of the normal derivative of u across Γ. From (4), it is clear that the differential problem solution is recovered at any point outside the obstacle and at any time instant, provided the density function is found. To this aim, performing a limiting process for x tending to Γ in (4) and using the assigned Dirichlet boundary condition (3), the weakly singular space-time BIE Γ t 0 G(r, t − τ)ψ(y, τ) dτ dγ y = g(x, t), in the unknown ψ(x, t) can be obtained and written in the compact notation where V is the so-called single-layer space-time integral operator. Now, considering the well-known conservation law satisfied by the (real valued) solutions of the d'Alembert equation: integrating (1) over [0, T] × R 2 \ Γ , using integration by parts in space and taking into account that u and u t vanish for t = 0, one obtains the energy of the solution u at the final time of analysis T, defined by E (u, T) := 1 2 R 2 \Γ ∇ x u(x, T) 2 2 + 1 c 2 u 2 t (x, T) dγ x (9) which can be rewritten as The quadratic form appearing in the last term of (10) leads to a natural space-time weak formulation of the corresponding BIE (7) with robust theoretical properties, which allow to deduce stability and convergence of the related Galerkin approximate solution [31]. Hence, projecting (7), derived with regard to time, by means of test functions φ belonging to the same functional space of the unknown density ψ, the energetic weak problem finally reads: find ψ ∈ L 2 ([0, T]; H −1/2 (Γ)) such that, ∀φ ∈ L 2 ([0, T]; H −1/2 (Γ)), Let us note that Equation (11) can be equivalently written as For the interested reader, the theoretical analysis of the bilinear form in the left-hand side of (11) was carried out in the Ref. [31] where, under suitable assumptions, coercivity was proved with some technicalities. Now, since Γ has a parametric representation, the energetic weak problem (12) can be rewritten in the parameter space as (13) where J(·) is the parametric speed of the boundary Γ, that is, The left-hand side of (13) can be rewritten more extensively as b a J(s) with r = C(s) − C(z) 2 . Further, in the IGA setting, Γ is given by a parametric representation based on B-Splines or NURBS. In the following, we will summarize for brevity only the basic ideas of B-Splines construction. More precisely, the definition of a B-spline basis used in the IGA approach is briefly recalled, and the Cox-De Boor recurrence relation is then presented. Given a partition of the parametrization interval [a, b] ⊂ R, it follows that [a, b] = n−1 l=0 e l is made up by n subintervals e l = [σ l , σ l+1 ]; a polynomial spline space B of order k and maximal regularity on such partition is composed by piecewise polynomial functions of degree d = k − 1 with assigned regularity C k−2 at the breakpoints σ i , i = 1, . . . , n − 1. In this way, B is a subset of C k−2 [a, b]. It is quite easy to verify that the dimension of such space is The easiest way to define in B a B-spline basis B i,k (t), i = 0, . . . , N, with N + 1 = dim(B), is based on the usage of a recursion formula and can be described through two easy steps [33]. The first step consists in associating to B an extended knot vector Z = {ζ 0 , · · · , ζ N+k } whose elements constitute a non-decreasing sequence of abscissas, where {ζ k−1 , · · · , ζ N+1 } coincides with those given in (16). The remaining knots in Z, {ζ 0 , · · · , ζ k−2 } and {ζ N+2 , · · · , ζ N+k } form two sets of k − 1 knots called auxiliary left and right knots which, under the standard assumption of selecting an open extended knot vector, are defined as ζ 0 = · · · = ζ k−2 = ζ k−1 = a and ζ N+1 = ζ N+2 = · · · = ζ N+k= = b.
In the second step, the basis is defined by using the following recursion [33]: with (fractions with zero denominator are considered vanishing). Note that from the above recursive definition, it is easy to verify the nonnegativity of B-splines and that the support of B i,k is the subinterval [ζ i , ζ i+k ]. From now on, if not otherwise specified, we will suppose that the obstacle Γ is given by a B-spline parametric representation of order k over the interval [a, b], that is, where P i , i = 0, · · · , N are suitably given control points in R 2 .

IGA-Energetic BEM Discretization
For the approximation of the unknown boundary density ψ in the space variable, the functional background compels to fix space basis functions belonging to L 2 (Γ), but for smooth obstacles, one can require at least C 0 (Γ) regularity. In the IGA paradigm, boundary element basis functions will be defined by means of the B-spline basis related to the partition ∆, as described in the previous section; that is, w j (x) := B j,k (C −1 (x)), j = 0, . . . , N. The total number of DoFs in the space variable will be N DoF = N + 1. Note that the partition ∆ induces on the obstacle Γ a boundary mesh constituted by n non-overlapping curvilinear elements {I 0 , · · · , I n−1 }, I l = C(e l ), such that n−1 i=0 I i = Γ. For time discretization, a uniform decomposition of the time interval [0, T] with timestep ∆ t = T/N ∆t , N ∆t ∈ N + , generated by the N ∆t + 1 instants t k = k ∆ t, k = 0, · · · , N ∆t , is considered, and the time mesh is equipped by piece-wise constant shape functions. Note that, for this particular choice which will allow for analytical double-integration in time variables, shape functions, denoted by v k (t), k = 0, . . . , N ∆t − 1, will be defined as Hence, the approximate solution of the boundary integral problem at hand will be expressed as Substitutingψ(x, t) in (12) and considering test functions φ ih (x, t), for i = 0, . . . , N and h = 0, . . . , the Energetic BEM discretization coming from energetic weak formulation (12) produces the linear system of order N DoF · N ∆t , where matrix A has a block lower triangular Toeplitz structure, owing respectively to the choice (21) and the fact that kernel G in (5) depends on the difference between t and τ. Note that higher-order time basis functions would produce, instead, a block Hessenberg type matrix structure.
If we indicate with A ( ) the block obtained when ∆ hk = t h − t k = ∆t, = 0, . . . , N ∆t − 1, each block has dimension N DoF and the solution is obtained by block forward substitution, that is, at every time instant t , one solves a reduced linear system of the type: where α ( ) = α ( ) j and β ( ) = β ( ) j , j = 0, . . . , N. Procedure (25) is a time marching technique, where the only matrix to be factorized once and for all is the symmetric nonsingular A (0) diagonal block, while all the other blocks are used to update the right-hand side at every time-step. Owing to this procedure, one can construct and store only the blocks A (0) , · · · , A (N ∆t −1) with a considerable reduction of computational cost and memory requirement. Now, remembering (15) and (21), matrix entries in blocks A ( ) are expressed by the sum with where the outer time integral in (15) has already been analytically performed and the notation w l (s) = B l,k (s), w m (z) = B m,k (z) has been adopted. With a similar procedure, the generic element of the right-hand side is expressed as After a further analytical integration in time, the matrix entries in blocks A ( ) are of the form and H(c ∆ hk − r) models the wavefront propagation.
Using the standard element-by-element technique, the evaluation of every doubleintegral in (28) is reduced to the assembling of local contributions of the type where e (l) i , e (m) j belong to (or coincide with) the support of w l , w m , respectively. Further, in order to simplify the notation, in the sequel, the upper index will be dropped from e where and Note that, for what concerns the logarithmic term, the above splitting has already been adopted in the Ref. [7]. Setting Note that (33) can also be defined for z = s, extending its definition by continuity, since lim z→s R(s, z) = J 2 (s) .
The splitting of (29) is useful to separate three contributions, the first and second coming from the geometry of Γ and the last depending only on the weakly singular nature of the kernel. Actually, integral (30) can be evaluated as the sum of with For δ = 3, weakly singular integrals can occur because of the nature of the kernel V 3 . Therefore, efficient evaluation of double integrals (38) is particularly required when e i ≡ e j or e j = e i±1 , where r can vanish. Further, the square root function c 2 ∆ 2 hk − r 2 in (32) gives rise to unbounded space derivatives when c ∆ hk = r, lowering the efficiency of the standard Gauss-Legendre rule, and this has to be properly managed. Hence, different quadrature rules are needed to compute I (i,j) V δ , δ = 1, 2, 3. The adopted schemes will be briefly recalled at the end of the next section. The most challenging issue, however, is the possible reduction and splitting of the integration domain, due to the presence of the Heaviside step-function in (38), modeling wavefront propagation, which will be taken into account in the following section.

Reduction and Splitting of the Integration Domain in Space Variables
Due to the presence of Heaviside function H(c∆ hk − r) in (38), a reduction of the integration domain can occur and a careful subdivision of internal and external integration intervals is needed. In fact, the derivative of the external integrand function in the s variable can show step discontinuity in some points of the outer integration interval, as already studied in the Ref. [34] for the simplest case of the rectilinear scatterer. In absence of this careful subdivision, a very large number of quadrature nodes is needed for the outer numerical integration, in order to achieve single-precision accuracy. Automatic detection of possible subdivision points, both for the inner and the outer variables of integration, therefore becomes the key point for an accurate and effective evaluation of the double integrals (38).
While in the case of C = (Id, 0) these points are expressed in closed form [34], in the general case, they have to be found as roots of the non-linear equation either in the outer integration variable s or in the inner integration variable z. More precisely, once intervals of outer and inner integration parametric intervals e i and e j are fixed, one needs to find where the function related to the argument of the Heaviside function, is less than or equal to zero. The search for subdivision points has been differently implemented depending on the mutual position between inner and outer integration intervals, distinguishing three different geometrical situations: coincident elements, that is, i = j, contiguous elements, that is, i = j + 1 and i = j − 1, or separate elements. In the next subsections, for the sake of brevity, we will resume the first two scenarios, which are the most important due to the appearance of kernel singularity.
Assuming the discretization of the parametric domain is fine enough and the curve C is at least C 1 , we can approximate curvilinear elements I i by an expansion of the curve C on the corresponding parametric element e i ; for the approximated elements, possible subdivision points can be determined in closed form, as detailed in the subsequent sections. The subdivision points determined in this way are then used as good initial guesses in a Newton-type solver to get closer at double-precision accuracy to the actual roots of the F(s, z) function (40).
In the following, to simplify the notation, the wave propagation speed is set as c = 1.

Double Integration over Coincident Elements
In the case where the double integrals (38) are computed on coincident elements e i = e j , the corresponding curvilinear element I i on the scatterer Γ is approximated by a first-order expansion around the left endpoint σ i of e i . For s, z ∈ e i , we can write with σ, ζ ∈ [0, l i ], l i being the length of e i ; the positions of points C(s), C(z) on I i are then approximated by the positions of the corresponding points on the segment tangent to Γ at C(σ i ): the distance between points on the curve is then approximated by the distance along the tangent segment, that is, and, consequently, the nonlinear Equation (39) is approximated in this geometric setting by The double-integration domain can hence be approximated by the intersection of the square e i × e i with the domain where the Heaviside function is not trivial. Having set the inner integration extremes the description of the integration domain with regard to the outer integration variable s changes where the min and max in (45) switch branch, or, equivalently, where z 1 (s) = σ i and z 2 (s) = σ i+1 , that is, at As already explained, the outer integration interval needs to be split in correspondence of these points, whenever they lie in the considered element e i .
Let us remark that reduction and subdivision points determined by means of (45) and (46), respectively, are approximations of the exact roots of F(s, z) in (40) and might need to be refined in order to get the desired accuracy.
Hence, starting with the search of possible subdivision points in the outer integration variable s, let s (0) * be the approximation of a subdivision point of the outer integration variable, obtained by (46), associated to an equation of the form z * (s) = σ * , where z * (s) is one of the inner integration extremes and σ * is a suitable endpoint of e i (the symbol * denotes any index which makes the present notation consistent with the one employed above). In order to test if the approximation is good enough, we compute the residual F(s (0) * , σ * ) -if it is smaller than a fixed tolerance , the approximation is accepted; otherwise, the following Newton iteration is applied: until the residual stopping criterion F(s The proposed procedure, starting from good initial guesses, produced fast converging sequences of approximations in all performed numerical experiments. In Table 1, for example, we report a sequence of approximations of the outer integral splitting points for the parametric element e i = [0, 0.1] on a semi-circular curve parametrized by C(s) = (cos(s), sin(s)), having fixed ∆ hk = 0.05. Table 1. Initial guess (k = 0) and refinements of external integration domain splitting points along with the corresponding residual values for a semi-circular curve parametrized by C(s) = (cos(s), sin(s)), In Table 2, we see that on the parabolic curve C(s) = (s, s 2 ) (with non-unitary speed), using a coarser parametric element e i = [0, 0.2] and having fixed ∆ hk = 0.05, the initial approximations obtained by linearization (46) have a higher residual compared to the previous example; nonetheless, the residual goes to zero in double-precision accuracy within the third iteration of (47). Table 2. Initial guess (k = 0) and refinements of external integration domain splitting points along with the corresponding residual values for a parabolic curve parametrized by C(s) = (s, s 2 ), After all subdivision points in the s variable have been determined within the prescribed accuracy, the external integration domain is decomposed into subintervals [s 1 ,s 2 ] ⊆ e i withs 1 ≤s 2 , wheres 1 ,s 2 are either extremes of e i or subdivision points for the s variable. If [s 1 ,s 2 ] is not empty, at each selected outer quadrature node s ∈ [s 1 ,s 2 ], approximated extremes for the inner integral z (0) * (s) are computed from (45). For those that lie within e i , the residual F s, z (0) * (s) is computed: if it is smaller than , the approximation is accepted; otherwise, the following iteration is applied: until the residual stopping criterion F(s, z (k+1) * (s)) ≤ is satisfied. The inner integration endpointsz 1 (s) <z 2 (s) are hence set up after the last iteration. A sequence of approximations of the inner integration extremes is reported in Table 3, in relation to the parabolic curve C(s) = (s, s 2 ) and the parametric element e i = [0, 0.2] as in the previous Table, having fixed ∆ hk = 0.05 and the outer variable s = 0.1. Again, starting from an initial guess obtained by (45), the sequence converges in double precision accuracy after three iterations of (48).
The overall procedure allows to express the double integral (38) as the sum of contributions of the form accurately evaluated by means of quadrature schemes described in Section 4.3. Table 3. Initial guess (k = 0) and refinements of inner integration domain extremes, related to s = 0.1, and corresponding residual values for a parabolic curve parametrized by C(s) = (s, 0.0505572313270522 0.148523936468889 9.22 · 10 −12 5.67 · 10 −10 3 0.0505572314189283 0.148523931028912 0 · 10 0 3.21 · 10 −17

Double Integration over Contiguous Elements
When double integrals (38) are computed on contiguous element e i , e j , with j = i + 1, we consider the same linear approximation as in Section 4.1 for the curvilinear elements I i , I j ⊂ Γ, this time around the common endpoint σ i+1 = σ j .
Let s ∈ e i and z ∈ e j be expressed as l i , l j being the lengths of e i and e j , respectively; the positions of points C(s) ∈ I i and C(z) ∈ I j on Γ are then approximated by the positions of points on the common tangent segment to Γ at C(σ j ), expressed by and the distance between points on the curve is then approximated by the distance along the tangent segment, that is, Consequently, the subset of the rectangle e i × e j where the Heaviside function is non-zero is approximated by the inner integration upper extreme is therefore set to In the points where z(s) intersects either the lower or the upper horizontal edge of the rectangle e i × e j , we find subdivision points for the outer integration variable respectively. Starting from approximations of reduction and subdivision points in Equations (53) and (54), Newton iterations are applied, similarly to the procedure described in Section 4.1, to get closer to the exact roots of the function F(s, z) in (40).
Finally, after reduction and subdivision points have been determined within a prescribed tolerance, double integrals (38) are expressed, using a notation similar to that one employed in (49), as sums of contributions of the form wheres 2 ,s 1 coincide either with the endpoints of e i or with possible subdivision points for the outer integration variable s lying within e i . In the case where the two contiguous elements are in the opposite order, that is, e j = e i−1 , the same procedure and formulas can be used by simply swapping the indices of both of the elements and shape functions in Equation (38).

Quadrature Schemes
We briefly explain below the different quadrature schemes adopted according to regular and singular integrals, which are needed to numerically evaluate the double integral (38). Once computed, if required by the Heaviside function, the appropriate splitting of the outer and inner integration intervals, as discussed above, the integration (38), for δ = 1, 2, 3, can be split in contributions of the type whereē i ⊆ e i ,ē j ⊆ e j in relation to the possible performed decomposition. For δ = 1, looking at (32), a numerical issue arises in relation to the integrand function log c ∆ hk + c 2 ∆ 2 hk − r 2 , even if it is not singular for r → 0. In fact, the square root argument is always positive, but it can assume very small values, and in the limit for the argument tending to zero, the derivative of the square root with respect to the inner variable of integration z becomes unbounded. To overcome this difficulty, the inner integration in (56) has been performed considering the regularization procedure in the Ref. [35], which suitably pushes the Gaussian nodes towards the endpoints of the intervalē j and modifies the Gaussian weights in order to regularize integrand functions with mild boundary singularities (the higher certain positive integer parameters p, q, the smoother the regularization). We will indicate nodes and weights of this quadrature rule with {z The outer integral in (56) is numerically evaluated by a classical Gauss-Legendre rule, with nodes and weights {s i G , ω i G } N G i G =1 , finally having: For δ = 2, which involves the double-integration of the regular kernel (33), we operate with the product of two Gauss-Legendre rules, yielding To conclude, for δ = 3, the numerical treatment of a weakly singular kernel in (34) on coincident and contiguous elements has been operated through the quadrature schemes proposed in the Ref. [36], which are widely used for weakly singular kernels in the context of Galerkin BEM applied to elliptic problems. The adopted quadrature rule for the inner integration is based on classical Gauss-Legendre nodes, with Gaussian weights that are suitably modified in order to absorb the logarithmic singularity and depend on the outer integration variable. The new weights are indicated by ω log j G (s). The outer integration is performed by the regularization technique [35], finally yielding Otherwise, for the evaluation of (38) on separate elements when δ = 3, we have operated with the product of two Gauss-Legendre rules.

Numerical Results
In this section, several numerical results obtained by the IGA-Energetic BEM applied to the analysis of 2D soft-scattering problems by open arcs having B-spline or NURBS parametric representation are presented and discussed, showing, from a numerical point of view, convergence and accuracy of the proposed method, as well as its optimal computational performance in terms of DoFs.
Firstly, we consider a straight obstacle, setting C = (Id, 0), and treat the problem (1)-(3) already studied in the Ref. [34], where the error was analyzed only with regard to the space variable at a fixed time instant. Subsequently, we discuss the case where the obstacle is represented by a semicircumference. Then, we consider the problem (1)-(3) with a parabola arc as a scatterer, with the aim of testing the IGA-Energetic BEM approach when the parametric speed of Γ is different from 1. At last, the scatter is chosen in such a way that it changes its concavity.
For the segment and parabola arc, we employed a B-Spline parametric representation, while for the semi-circumference and the last obstacle, we employed a quadratic NURBS representation. Numerical results reported below are related to simulations based on keeping the space mesh fixed and increasing the B-splines degree, instead of refining the space mesh keeping the B-spline order fixed, determined by Γ representation: this choice better highlights the benefits of using an IGA instead of a standard Energetic BEM approach, which employs the C 0 Lagrangian basis functions to represent the approximate solution, for what concerns the computational saving in terms of DoFs. In any case, it is worth noting that the boundary curve is exactly expressed in all the simulations, because its representation can be obtained by a standard degree elevation procedure [37].
In each simulation, we set the wave propagation velocity c = 1 and we consider uniform decompositions of the parametrization interval [a, b], setting ∆x = l i , ∀ i. This choice, in general, does not imply a uniform decomposition on Γ, as in the case of a parabola arc.
The numerical study of the approximate solution will be done in squared space-time energy norm, which, remembering (10), for ψ on Σ T and in the proper functional space, is defined by: For the considered numerical tests, to the authors' knowledge, there are no exact spacetime solutions in closed form. Nevertheless, by using a boundary datum which becomes independent of time, the IGA-Energetic BEM solution stabilizes to a time-independent function. Therefore, the transient solutionψ evaluated on Γ at a fixed (large) time instant can be compared with the analytical solution ψ ∞ , if available, of the BIE related to the stationary differential problem with g ∞ (x) = lim t→∞ g(x, t) chosen in such a way that the Sommerfeld radiation condition (63) is satisfied [38].

Straight Scatterer
Let us consider the model problem (1)-(3), where the obstacle Γ is described by a B-Spline parametric representation of order 2 over the interval [−0.5, 0.5], where P 0 = [−0.5 0] and P 1 = [0.5 0] are the two control points, and Z = [−0.5 − 0.5 0.5 0.5] is the extended knot vector. We fix the Dirichlet boundary datum given by This is due to the fact that the boundary datum becomes independent of time, so a convergence of the transient solution ψ on Γ to a stationary solution can be expected. The analytical solution ψ ∞ of the stationary BIE on Γ related to the limit problem (61)-(63), with g ∞ (x) = −x, is known in closed form and reads This analytical solution clearly presents mild singularities at the endpoints of Γ.
In Table 4, the relative error ψ (·, 10) between the approximate transient solution at time instant T = 10 and the limit stationary solution is reported, varying the degree of the B-spline basis, together with the space DoFs. In the same  [34]). This result is the same as what has been known since the Ref. [39] for elliptic problems. The comparison between the two error decays is reported in Figure 2 with regard to DoFs. The plot shows that the IGA approach, despite having nearly the same speed of convergence of the standard approach with regard to d, presents a higher rate with regard to DoFs. In fact, for the Lagrangian basis, it holds that d N DoF /n, while for the B-spline basis, we have d = N DoF − n. In the same figure, the decays O(N −1 DoF ) and O((N DoF − n) −1 ) are also shown. This behavior highlights and confirms the same computational benefits as observed for IGA-BEM applied to elliptic problems in the Ref. [5]. Table 4. Relative error ψ (·, 10)   For completeness, in Figure 3 (left), the whole time-history of the relative error in H −1/2 (Γ) norm is reported for the considered B-spline degrees: at the beginning of the simulations, the errors are nearly the same for all the considered B-spline bases, but the curves soon become separated and stagnate at different levels; in this case, the quality of the approximation is bounded by the ability of the discretization space to catch the singularity of the solution at the endpoints. Increasing the B-spline degrees, the space approximation becomes more capable of representing the solution behavior, so better results are achieved, even if the relaxation to equilibrium requires more time. The same remarks can be made for Figure 3 (right), related to the standard Energetic BEM.

177
In Table 5  In Figure 4, we show the approximate solutionsψ(x, 10), obtained by IGA-Energetic BEM for different values of degrees d of the B-spline basis, at the final instant of analysis T = 10 on Γ. Additionally, ψ ∞ (x) has been reported in Figure 4 and, as one can observe, the approximate solutionsψ(x, 10) tend to overlap the curve related to the stationary solution ψ ∞ (x), diminishing their oscillatory behavior near the endpoints of Γ for growing d.  Of course, the complete analysis of numerical results has to be done in the whole space-time domain, but unfortunately, the analytical space-time solution is not known in closed form. Nevertheless, we can proceed with the following numerical study, employing the value of space-time energy norm of the approximated solution as a mean to assess the convergence of the method for finer discretization, both in space, by p-refinement, and, in time, by h-refinement.
In Table 5 the space-time energy norm squared (60) is reported for both IGA-and standard Energetic BEM, for growing degree d of the B-spline basis and diminishing timestep defined as ∆t = 2 −d ∆x, having fixed ∆x = 0.1, together with the space DoFs. As one can observe, these values tend to stabilize towards a limit value conceived as the spacetime energy norm of the exact solution of the BIE wave problem, but the convergence of IGA-Energetic BEM is much faster with regard to DoFs, as already shown above for a fixed time instant. This behavior is perfectly visible in Figure 5.  The analysis of the space-time energy norm for the approximate solution given by IGA-Energetic BEM is also reported in Figure 6, versus both increasing B-spline degrees and refined time-steps, highlighting the contribution of these two parameters to the convergence to the reference value. Each curve in the left plot corresponds to a fixed time-step, and the B-splines degree is increased along each curve; moving along the vertical axis, we see that for every fixed B-spline degree, the values of the energy accumulate around a limiting value. Moving horizontally, we instead see that the energies of the approximated solutions stabilize at a value which depends on the time-step. Finally, the two contributions can be combined by extracting the diagonal sequence (black line), reported also in Table 5 and Figure 5. The plot on the right reports the same data, where for each fixed B-spline degree, we can observe that the values of the energy of the approximated solutions converge under time-step refinement.

Semi-Circular Scatterer
Let us consider the model problem (1)-(3), where the obstacle is a semi-circular arc of unitary radius: This curve does not admit an exact B-spline representation, but can be parametrized using NURBS (see [40]) over the interval [0, π] as The model problem is equipped by the Dirichlet boundary datum given by where f is given in (65) and Looking at the time history of the approximated boundary density computed at α = π 8 , reported in Figure 7, we see again that after the initial transient phase, the approximated solution stabilizes at a constant value, which is the value of the stationary density ψ ∞ at the same point of the boundary. The stationary solution of the BIE on the semi-circumference Γ related to the limit problem (61)-(63) is known to be ψ ∞ (α) = cos(α).
Since the limiting stationary solution (70) is regular, we can expect the space dependence to be easier to approximate than in the previous example, even with low-degree shape functions and coarse meshes over Γ. In fact, having fixed ∆α = π 10 , ∆t = 0.3 in Figure 8, we see that the computed density at final time T = 1200 approximates the analytic stationary density very closely. In this case, the stationary density ψ ∞ being smooth, the approximations for degrees ≥ 2 are all overlapped. Moreover, since in this example the stationary solution is defined and smooth everywhere, we can compute, at every time instant t, the relative error ψ (·, t) − ψ ∞ (·) H 0 (Γ) / ψ ∞ (·) H 0 (Γ) in the H 0 (Γ) norm between the approximate transient solution and the limit stationary solution over the obstacle Γ. In Figure 9, the time history of the relative error obtained by IGA (left) and standard (right) Energetic BEM is shown, and we observe that at the beginning of the simulations, the errors are practically the same, meaning that in the transient phase, the dominating error is due to the time variable. Then, low-order shape functions stagnate first, while the curves related to shape functions of higher degree proceed whilst packed together. At large enough time instants, the curves separate one by one and stagnate as the time-dependent part of the error fades away, revealing the space discretization error. In particular, in Figure 9 (left) we see that the curve related to d = 4 is starting to stagnate at the end of the simulations, meaning that the equilibrium has almost been reached, while B-splines of degree 5 and 6 are still equilibrating and have the same relative error. Using Lagrangian shape functions instead, it appears from Figure 9 (right) that the approximation of degree 4 is still improving at the end of the simulation and would give better results on a longer time interval. Let us stress the fact that this example is different from the previous one on the straight scatterer where the space component of the error was almost immediately dominant.
In Table 6, the above-defined relative error at the final time-instant of analysis T = 1200 is reported, varying the degree of the B-spline basis, together with the space DoFs. In the same Table, the analogous results obtained by standard Energetic BEM, employing C 0 Lagrangian basis functions to represent the approximate solution, are reported. As one can see, for the selected T, the errors stagnate for degrees higher than 4, even if they reach the single precision accuracy.
To investigate the influence of the time discretization on the convergence of the transient solution, we repeat the same simulation on [0, T] = [0, 1200] decomposed in 8000 time-steps of size ∆t = 0.15. Looking at Figure 10, we can observe, as expected, a finer resolution, capturing more minute details in the initial transient phase, but the limit values of the error, coming from the space discretization are, as expected, the same as in Figure 9.
To conclude, we present the numerical results related to the space-time energy norm squared (60), reported in Table 7.
10 0 10 1 10 2 10 3 10 −7        In Figure 11 (right) we see that the energy values of approximations obtained with finer time-steps stabilize at a fixed value, indicating convergence by time-step refinement. We see instead that all degrees give almost the same energy values. Similarly, in Figure 11 (left) the main variation to the 'diagonal' sequence of refinements with discretization parameters is given by time-step refinement. This can be explained observing that the Dirichlet datum (68) becomes independent of time very quickly, hence the space-time energy norm of the solution (60) only accounts for the first time-instants during the initial transient phase. Since, at the beginning of the simulation shape functions of all degrees give almost the same results, as visible from the previous errors plots, it is not surprising that degree elevation has almost no effect on the space-time energy in this example.

Wave Scattering by Parabola Arc
Let us consider the model problem (1)- (3), where the obstacle Γ is described by a B-Spline parametric representation of order 3 over the interval [−1, 1]: with control points P i , i = 0, 1, 2 whose coordinates are given in the i-th column of the following matrix P: and where Z = [−1 , −1 , −1 , 1 , 1 , 1] is the extended knot vector. The arc is shown in Figure 12, together with its control points. We fix the Dirichlet boundary datum given by where f has already been given in (65) and in such a way that the analytical solution ψ ∞ of the stationary BIE on Γ was related to the limit problem (61)-(63), with g ∞ (x) = q(x), it is known in closed form and reads At first, we consider the time interval [0, T] = [0, 100], fixing the time-step ∆t = 0.2 and a uniform decomposition of the parametrization interval [−1, 1] constituted by 10 elements equipped by quadratic B-Splines. Note that the mesh lifted on Γ by the considered uniform partition of the parametrization interval is not uniform on the scatterer. Figure 13 presents the time history of the BIE numerical solutionψ(x, t), having set x = −1, zooming in on the time interval [0, 60] to show the initial oscillations, and simi-larly to the previous cases, the IGA-Energetic BEM solution stabilizes at a constant value, showing no long-time instability. Hence, convergence of the transient solutionψ on Γ to ψ ∞ can be expected. In Table 8, the relative error ψ (·, 100) − ψ ∞ (·) H 0 (Γ) / ψ ∞ (·) H 0 (Γ) between the approximate transient solution at the final time instant of analysis T = 100 and the limited stationary solution is reported, varying the degree of the B-spline basis, together with the space DoFs. This norm is applicable due to the regularity of the solution. In the same Table, the analogous results obtained by standard Energetic BEM, employing C 0 Lagrangian basis functions to represent the approximate solution, are reported: as one can see, for the selected T the errors stagnate for degrees higher than 2. Nevertheless, decreasing errors behave nearly as O(∆x d+1 ) but, due to the already reported different representations of d in terms of space DoFs, IGA approach also has, in this case, a higher rate with regard to DoFs, as shown in Figure 14. This behavior highlights and confirms the same computational benefits, as observed for IGA-BEM applied to elliptic problems in the Ref. [5]. Table 8. Relative error ψ (·, 100) − ψ ∞ (·) H 0 (Γ) / ψ ∞ (·) H 0 (Γ) varying d.  For completeness, in Figure 15 (left) the whole time-history of the relative error in energy norm is reported for the considered B-spline degrees. In this example, the errors related to all the considered B-splines bases remain nearly the same for long times, and the curves become separated when the error in time, for growing time, becomes smaller than the error due to the chosen B-spline degree and visible at different levels-the smaller this space error, the longer the time needed for the separation of the curve from the others. Thus, better results are achieved for higher degrees, even if the relaxation to equilibrium requires more time. The same plot related to standard Energetic BEM can be seen in Figure 15 (right), where curves related to Lagrangian basis degrees greater than 2 are not yet separated at T = 100. This is the reason for stagnation in Table 8.  x ψ(x,  In Figure 16 we show the approximate solutionsψ(x, 100), obtained by IGA-Energetic BEM for growing values of degree d of the B-spline basis, at the final instant of analysis, T = 100. Additionally, ψ ∞ (x) has been reported in Figure 16 and, as one can observe, all the curves related to the approximate solutions overlap the curve related to the stationary one.    x ψ(x,  In the previous simulation, most of the approximated solutions did not reach equilibrium within the final time, as prominently suggested by the error curves in Figure 15. We therefore extend the time interval to [0, T] = [0, 300], keeping the same discretization parameters. Figure 17 (left) shows that by time T = 300, error curves for B-spline approximations are almost flat, meaning that the equilibrium has been reached. From this time instant, the dominating error comes from the space discretization. Using Lagrangian shape functions, we can see in Figure 17 (right) the separation of the curve related to d = 3 from the pack of higher-degree shape functions; the overlap of curves related to d = 4, 5, 6 at final time T = 300 suggests that here, the dominating contribution still comes from the error in time, and better approximations of the stationary solution can be achieved with even longer simulations.   In Table 9, the relative error ψ (·, 300) − ψ ∞ (·) H 0 (Γ) / ψ ∞ (·) H 0 (Γ) between the approximate In Table 9, the relative error ψ (·, 300) − ψ ∞ (·) H 0 (Γ) / ψ ∞ (·) H 0 (Γ) between the approximate transient solution at the final time instant of analysis T = 300 and the limit stationary solution is reported, varying the degree of the B-spline basis, together with the space DoFs. In the same Table, the analogous results obtained by standard Energetic BEM, employing C 0 Lagrangian basis functions to represent the approximate solution, are reported: as one can see, now for the selected T, the errors stagnate for degrees higher than 3. As before, decreasing errors almost behave like O(∆x d+1 ). We now refine the space mesh to 20 elements of size ∆x = 0.1, keeping the time discretization parameter fixed to ∆t = 0.2 for the time interval [0, T] = [0, 300].
B-splines built on this refined space mesh allow for a much more accurate approximation of the stationary solution; as we can see in Figure 18 (left), the equilibrium values of the errors are significantly lower compared to those obtained on the coarse mesh. Moreover, in Figure 17 (left) all curves are flat after T = 200, leaving no room for improvement on longer simulations, while on the finer space mesh, we instead see that curves for d = 4, 5, 6 still overlap at T = 300, meaning that the error from the space discretization will reach a lower value at equilibrium reached for longer times. Using Lagrangian shape functions, we get very similar results, visible in Figure 18 (right). Looking at Figure 19, we see that, for the degrees that reached equilibrium within final time, the convergence rate is close to the expected O ∆x d+1 .     As done in the previous numerical example, we show results in the whole space-time domain, even if, as before, the analytical space-time solution is not known in closed form. In Table 10, the space-time energy norm squared (60) is reported for both IGA-and standard Energetic BEM, for a growing degree d of the B-spline basis and diminishing time-step defined as ∆t = 2 −d ∆x, having fixed ∆x = 0.1, together with the space DoFs. As one can observe, these values tend to stabilize towards a limit value conceived as the space-time energy norm of the exact solution of the BIE wave problem, but the convergence of IGA-Energetic BEM is much faster with regard to DoFs, as already shown above for a fixed time instant. This behavior is perfectly visible in Figure 20.     To conclude, in Figure 21 (right) we can observe the expected convergence of the space-time energy norm with regard to growing B-spline degree and time-step refinement (right). As already noted in the previous example, degree elevation makes hardly any difference, because for the given Dirichlet datum this type of norm takes into account the initial transient phase of the simulation, where we already observed that the error curves are packed together; the main variations come from time-step refinement.     At last, let us consider the obstacle Γ given by This scatterer can be described by a NURBS parametric representation of order 3 over the interval [0, 2 π], i.e.,

S-Shaped Scatterer
At last, let us consider the obstacle Γ given by This scatterer can be described by a NURBS parametric representation of order 3 over the interval [0, 2 π], that is, where B i,2 , i = 0, . . . , 4 are the quadratic B-spline basis functions used in the construction of NURBS, P i , i = 0, . . . , 4 are the control points, whose coordinates are given in the i-th column of matrix P: and Z = [0 , 0 , 0 , π/2 , π/2 , π , π , 3π/2 , 3π/2 , 2π , 2π , 2π] is the related extended knot vector. The scatterer is shown in Figure 22, together with its control points.
Let us consider the model problem (1)-(3), fixing Dirichlet boundary datum given by where f has already been given in (65). In this case, neither the analytical space-time solution nor the limit stationary solution in closed form, as in the previous cases, are known. In Figure 23, the Energetic BEM approximate solution, obtained at T = 60 on Γ using quadratic B-splines and fixing ∆x = π/10 and ∆t = 0.3, is shown. Furthermore, in Table 11 the space-time energy norm squared (60) is reported for both IGA-and standard Energetic BEM, for growing degree d of the B-spline basis and diminishing time-step defined as ∆t = 2 −d ∆x, having fixed ∆x = π/30, together with the space DoFs. As one can observe, these values tend to stabilize towards a limit value representing the space-time energy norm of the exact solution of the BIE wave problem, but the convergence of IGA-Energetic BEM is much faster with regard to DoFs. This behavior is perfectly visible in Figure 24.   To conclude, in Figure 25 (right) we can observe the expected convergence of the space-time energy norm with regard to growing B-spline degree and time-step refinement (right). As already noted in the previous examples, the main variations come from time-step refinement, due to the fact that for the given Dirichlet datum this norm takes into account the initial transient phase of the simulations.

Conclusions
In this work, the numerical investigation of a B-spline-based IGA approach to the so-called Energetic BEM, here applied to the numerical solution of 2D wave propagation exterior problems equipped by a Dirichlet boundary condition (i.e., acoustic soft-scattering), was carried out. The paper completes and enriches the preliminary work [34]. In particular, the simulations computed on curvilinear obstacles Γ confirmed the same computational advantages already observed in IGA-BEM for elliptic problems, with regard to standard Energetic BEM based on Lagrangian basis functions for the representation of the approximate solution for what concerns space variables. In comparison with standard Energetic BEM, IGA-Energetic BEM appears really competitive in reducing the number of DoFs, and this gain is obtained at every step of the time-marching procedure. Further investigations are needed to apply the so-called basis-function-by-basis-function approach instead of the classical element-by-element technique, principally for what concerns the adoption of suitable quadrature rules with nodes not on the element, but directly on the whole B-spline support, as done in the Refs. [12,41] for IGA-BEM applied to elliptic problems. Currently, the investigation is focused on the IGA implementation of Energetic BEM involving the remaining double layer and hypersingular integral operators, as well as on theoretical analysis of convergence and accuracy. Future research directions will also include the extension of IGA-Energetic BEM to efficiently solve elastodynamics problems.

Data Availability Statement:
The data presented in this study are available on request to the authors.