Programming finite element methods with weighted B-splines

Computers and Mathematics with Applications ( ) – Contents lists available at ScienceDirect Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa Programming finite element methods with weighted B-splines Klaus Höllig, Jörg Hörner ∗ IMNG, Universität Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany article info Article history: Available online xxxx Keywords: B-splines Finite element method Numerical integration abstract We discuss the main components of the recently developed FEMB program package, which implements finite element methods with weighted B-splines for basic linear elliptic boundary value problems in two and three dimensions. A three-dimensional implementation without topological restrictions has not been available before. We describe in particular the mathematical background for the recursive quadrature/cubature over boundary cells and explain how to utilize the regular data structure of uniform B-splines efficiently. Considering the Lamé–Navier equations of linear elasticity as a typical example, we illustrate the performance of the main FEMB routines. The numerical tests confirm that, due to the new integration routines, the weighted B-spline solvers have become considerably more efficient. © 2015 Elsevier Ltd. All rights reserved. 1. Introduction B-splines have become standard tools in approximation, computer graphics, design, and manufacturing. Recently, they have also been used to construct basis functions for finite element methods. The resulting techniques combine the computational efficiency of regular grids with the geometric flexibility of classical finite elements. Some key advantages are the free choice of order and smoothness, a simple data structure with one parameter per grid point, and the exact representation of boundary conditions. Over the past decade, essentially two B-spline based techniques have been developed: finite element methods with web-splines (weighted extended B-splines) [1] and isogeometric analysis [2], respectively. The books [3,4] provide a comprehensive treatment of the relevant theory. Which technique is best suited for a particular problem depends to some extent on the topological form of the simulation domain and its representation. Weighted methods are a good choice for problems with natural (Neumann) boundary conditions or if the part of the boundary, where essential (Dirichlet) boundary conditions are prescribed, has a convenient implicit description. Isogeometric methods can handle domains well which are parametrized over rectangles or cuboids or which can be expressed as union of few such parametrizations, e.g., NURBS representations for CAD/CAM applications. Fig. 1 shows two typical examples. There are also problems where a combination of both techniques might be useful [5]. In this article we discuss the implementation of weighted finite element methods based on the FEMB program package, which has been recently developed and is available on the website [6] (latest update FEMB1.0.1, November 2014). The first two sections are devoted to the construction of finite element bases with weighted B-splines, mainly providing the necessary background for the description of our algorithms. Then, we briefly describe in Section 4 the scope and structure of the FEMB routines. The key algorithmic components, numerical integration and solution of Ritz–Galerkin systems, are discussed in ∗ Corresponding author. E-mail addresses: Klaus.Hoellig@mathematik.uni-stuttgart.de (K. Höllig), Joerg.Hoerner@mathematik.uni-stuttgart.de (J. Hörner). http://dx.doi.org/10.1016/j.camwa.2015.02.019 0898-1221/© 2015 Elsevier Ltd. All rights reserved. 2 K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – Fig. 1. Domains suited for weighted (left) and isogeometric (right) finite element methods with B-splines. Fig. 2. Uniform B-spline. Sections 5 and 6, respectively. In particular, we establish the mathematical justification for the recursive boundary cell quadrature/cubature. Our original implementation, made available in 2012, was partly based on heuristic considerations and required minor modifications to fit the theory of Section 5. Finally, in Section 7, we illustrate the performance of weighted B-spline approximations for a model problem in linear elasticity. To a large extent due to the new integration algorithm, the FEMB package now provides an efficient implementation also for three-dimensional applications. It replaces the software, developed along with the textbook [4] which was limited to two-dimensional boundary value problems. 2. B-Spline bases There are several different generalizations of univariate B-splines, each giving rise to a very rich mathematical theory. From a computational point of view, and in particular for finite element schemes, the simplest multivariate concept, based on forming tensor products, is adequate [7–9]. Definition 2.1 (Uniform B-Spline). A d-variate uniform B-spline of degree n > 0 and grid width h is a product of scaled translates of univariate cardinal B-splines: bk (x) = d  ν=1 b(xν /h − kν ), with k = (k1 , . . . , kd ) ∈ Zd , x = (x1 , . . . , xd ) ∈ Rd , and b the univariate B-spline with knots 0, . . . , n + 1. As is illustrated in Fig. 2, bk is a nonnegative, (n − 1)-times continuously differentiable bell-shaped function with support supp bk = kh + [0, n + 1]d h. On each grid cell Q = ℓh + (0, 1)d h it coincides with a d-variate polynomial of degree n in each coordinate. The center of the support (red dot) is often used to identify the B-spline on the grid. Definition 2.2 (Splines). The spline space for a domain D is spanned by all relevant B-splines, i.e., those B-splines which have some support in D: B=  bk . k∼D For computations it is convenient to keep also irrelevant indices, i.e., to work with linear combinations  k∈K uk bk , where K is the smallest rectangular array containing all relevant k, and the coefficients uk with k ≁ D, are set to zero. K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – 3 Fig. 3. Splines on a bounded domain. Fig. 4. Representing essential boundary conditions with R-functions. The definition is illustrated for biquadratic B-splines in Fig. 3, where relevant (some support in D) and irrelevant (no support in D) B-splines are marked with dots and circles at the centers of their supports, respectively. Clearly, using scaled translates of a single basis function leads to an exceptionally simple data structure and has obvious computational advantages. In particular, no index lists are necessary as in conventional finite element programming, and many expressions can be precomputed. Splines provide suitable finite element subspaces for problems with natural boundary conditions. Whenever the weak formulation of an elliptic problem is based on the Sobolev spaces H k (D) without any constraints on the boundary values, just restricting uniform B-splines to the domain D yields admissible finite elements. A simple example is the minimization of 1 2  D   |grad u|2 + u2 −  fu D over H 1 (D). Essential boundary conditions, however, must be incorporated into the spline space. This can be done via suitable weight functions, as will be discussed in the next section. A slight drawback of the spline space B is that, in general, the B-splines do not provide a stable basis. This is due to Bsplines near the boundary with very small support in D. Such B-splines can be adjoined to neighboring B-splines with an extension procedure suggested by Höllig, Reif, and Wipper in [1]. The resulting extended B-splines are uniformly stable with respect to the grid width h. For multivariate B-spline bases, the importance of stability is often overestimated. We have found that stabilization is not necessary in order to obtain accurate approximations. Moreover, simple preconditioning procedures appear to be sufficient to solve the ill-conditioned finite element systems with acceptable accuracy. An application, for which stability seems indispensable, is multigrid methods [10]. So far it is still an open problem whether the extension procedure, required for theoretical purposes, is essential also for an optimal implementation. 3. Weight functions The idea of using weight functions to represent essential boundary conditions is due to Kantorowitsch and Krylow [11] and has been intensively investigated by Rvachev et al. (cf., e.g., [12]). Fig. 4 illustrates this method with an elementary example. The B-shaped domain D, shown in the right part of Fig. 4, is obtained as the union of two discs intersected with a half space: D = (D− ∪ D+ ) ∩ H with D± : w± (x) = r 2 − (x1 ∓ 1)2 − x22 > 0, H : w0 (x) = x2 > 0. 4 K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – Table 1 Standard R-function system. Set operation Complement: D Corresponding R-function c Union: D1 ∪ D2 Intersection: D1 ∩ D2 rc (w) = −w  w2 + w22  1 r∩ (w1 , w2 ) = w1 + w2 − w12 + w22 r∪ (w1 , w2 ) = w1 + w2 + If homogeneous Dirichlet boundary conditions are prescribed at the lower (flat) boundary, multiplying the B-splines bk by w0 yields an admissible finite element subspace: w0 B =  w0 bk . k∼D To approximate functions u which vanish on the entire boundary of D, one is tempted to use the elements w0 w− w+ bk . However, while these products obviously vanish on ∂ D, the zeros in the interior of the domain cause a severe loss of accuracy. Only functions u with the same zero pattern can be well approximated. Rvachev’s R-function method overcomes this problem in an elegant fashion. A weight function w , which vanishes on ∂ D and is positive on D, is constructed by combining the elementary weight functions according to the Boolean operations leading to the B-shaped domain. In its simplest form, the method employs the R-function system of Table 1. This yields the weight function wunion = w− + w+ +  2 2 w− + w+ for the union of the two discs, and w = wunion + w0 −  2 wunion + w02 is a positive weight function for Dirichlet boundary conditions on ∂ D. This construction is illustrated in Fig. 4, showing the weight functions wunion (left) and w (right). The simple motivating example leads to the following definition. Definition 3.1 (Weight Function). A weight function w associated with a part Γ of the boundary of the simulation domain D has bounded gradient, is positive on D \ Γ , and vanishes on Γ . The weighted spline space wB =  w bk k∼D provides admissible finite elements for problems with essential homogeneous boundary conditions on Γ . Several techniques for constructing weight functions are available (cf., e.g., sections 4.3 and 8.3 in [4]). As illustrated above, the R-function calculus provides a convenient tool for defining weight functions for objects arising in constructive solid geometry [12]. For general domains with smooth free form boundaries the distance function w(x) = dist(x, ∂ D) is a possible weight function. To avoid discontinuous gradients one can blend w with a plateau of fixed height outside of a small boundary strip. This has the additional advantage that w B = B on the major portion of the domain. The evaluation and differentiation of distance functions can be time-consuming since one has to resort to numerical procedures. Hence, one might want to replace w by a spline approximation. The work by Bajaj et al. on A-splines (cf., e.g., [13]) could be a starting point for using such implicit spline representations in the context of finite element methods. 4. FEMB program package The FEMB routines [6] implement finite element methods with weighted B-splines for basic elliptic model problems in two and three dimensions. The programs are intended for educational purposes with an emphasis on simplicity and readability rather than generality and professional standards. Each problem considered has a variational formulation a(u, v) = ⟨f , v⟩, v ∈ H, where a is an elliptic bilinear form and ⟨f , ·⟩ is a bounded linear functional on a Hilbert space H . These familiar equations admit a standard Ritz–Galerkin approximation. Hence, the algorithmic steps are familiar from conventional finite element codes: (A1) specify weight functions for the domain and essential boundary conditions; (A2) compute integration parameters; K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – 5 Fig. 5. Deformation of a rotating solid of revolution. (A3) assemble the Ritz–Galerkin system; (A4) solve the Ritz–Galerkin system; (A5) evaluate the solution and the residual. As an illustration of these program components we consider one of the FEMB demos, the program demo_elasticity_3d. It computes the deformation u(x) of a rotating elastic solid of revolution, as depicted in Fig. 5. This application will serve as a basic example throughout the remainder of this article. The symmetry of the domain is not exploited by the B-spline discretization (e.g., by using polar coordinates). Hence, there are no essential simplifications compared to general implicit free-form boundary representations so that the example provides a legitimate test case (other examples, suitable for assessing the accuracy of the FEMB algorithms, are the programs example_bvp_2d/3d_convergence). According to the physical model, Lu = −div σ (u) = f in D, (u1 , u2 , u3 ) = (0, 0, 0) on Γ , with σ the stress tensor, (f1 , f2 , f3 ) the body force (equal to the centrifugal force in the example considered), D the volume occupied by the solid, and Γ the boundary of the cylindrical axis of rotation. The weak form of this elliptic boundary value problem is a(u, v) =  D σ (u) : ε(v) =  fv D (1) ∀v ∈ H , where ε is the strain tensor and H consists of all displacements which vanish on Γ and have square integrable gradients on D (cf. [4, Section 6.5]). The Ritz–Galerkin approximation is defined in the usual way. Using B-splines, we have to incorporate the essential boundary conditions, as described in the previous sections. With the weight function w(x1 , x2 , x3 ) = (x1 − 1/2)2 + (x2 − 1/2)2 − r 2 representing the cylindrical boundary (w = 0 on Γ ) where u = (0, 0, 0), the finite element approximation w uh has the form w (uh1 , uh2 , uh3 ) = w  k∈K (uk,1 , uk,2 , uk,3 ) bk = w 3  uk,ν bk eν , k∈K ν=1 with e1 , e2 , e3 the unit vectors in R3 . Since we keep irrelevant B-splines (with zero coefficients), the index set K is a rectangular array. For example, if D ⊆ (0, 1)3 , as is assumed in the FEMB package, then K = {−n, . . . , 1/h − 1}3 . The analysis of the approximation properties of weighted approximations w uh is not straightforward. In general, polynomials (in particular linear functions) cannot be represented exactly. Nevertheless, weighted spline spaces possess the optimal approximation order (polynomial degree plus one) if both the function u to be approximated and the boundary of the domain D are smooth (cf. [4, Sections 5.4 and 5.6]). In case of singularities of u, caused e.g. by edges and corners of the simulation domain, the approximation order of (uniform) weighted elements usually deteriorates. One has to resort to hierarchical refinement to obtain high accuracy solutions with a moderate number of parameters. Let us describe the program demo_elasticity_3d in some more detail. In addition to the weight function w for representing the essential boundary condition, the user  has to specify the profile R of the solid. The choice corresponding to Fig. 5 is r = 1/5 in the definition of w and R(x3 ) = 1 + x23 /3. This leads to a second weight function   wD (x) = w(x) R(x3 )2 − (x1 − 1/2)2 − (x2 − 1/2)2 , (2) 3 describing the domain in implicit form: D = {x ∈ (0, 1) : wD (x) > 0}. Equivalently, D : (1/5)2 < (x1 − 1/2)2 + (x2 − 1/2)2 < (1 + x23 )/9, 0 < x3 < 1. The bounds on x3 need not to be taken into account in the definition of wD since they conform to the B-spline grid. (3) 6 K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – The definitions of w and wD complete step (A1). Steps (A2–A5) are implemented in separate Matlab1 functions. The main input variables are the grid width h = 1/H with H ∈ N, the degree n, the force function f , and the material constants (Young modulus E and Poisson ratio ν ). A variable PAR, which can be set automatically with the statement PAR = set_par(H,n), contains various algorithmic parameters. The program int = integrate_3d(w_D,H,PAR) determines Gauß points and weights for each grid cell Q in (0, 1)3 (step (A2)). This is nontrivial only if Q intersects the domain boundary and will be explained in detail in the next section. The program [G,F] = assemble_3de(E,nu,f,w,H,n,int,PAR) assembles the Ritz–Galerkin system which is solved with the program U = solve_3de(G,F,PAR) (steps (A3, A4)). Finally, the program [wuXYZ,X,Y,Z,R,e] = evaluate_3de(E,nu,f,wD,U,H,n,PAR) evaluates the finite element approximation w uh and the residual R = L(w uh ) − f on a regular XYZ-grid and computes the relative error in the L2 -norm: e = ∥R∥2 /∥f ∥2 . While f ̸= 0 for the model problems considered in the FEMB program package, for problems without a force term, a relative error could be computed by dividing by the norm of the numerical approximation or, alternatively, by the norm of boundary data (e.g., inhomogeneous Neumann boundary values). In contrast to standard finite element approximations, computing the residual is possible, since for n > 1 the weighted B-spline basis functions are continuously differentiable. Clearly, substituting back into the partial differential equation Lu = f provides a very reliable assessment of the accuracy of w uh . So far, the FEMB programs employ only uniform grids. Clearly, it would be natural to use the local size of the residual to control adaptive refinement. Hierarchical bases naturally fit into the B-spline framework (cf. [4, Sections 4.5] and [14]). Perhaps it is also possible to use weighted B-splines of different degrees. Both topics are subject of future research guided by the comprehensive theory available for classical mesh-based finite elements. 5. Numerical integration The assembly of the Ritz–Galerkin system involves integrals of the form IQ p =  (4) p, Q ∩D where the region of integration is the portion of the grid cell Q = ℓ/H + (0, 1)d /H in the domain D and the integrands p are smooth functions, computed from the finite element basis. If Q does not intersect the boundary of D, i.e., if Q ∩ D = Q , tensor product Gauß formulas yield very accurate approximations to IQ p. If m nodes are used in each coordinate direction, the errors of Gauß approximations to p over a grid cell Q and over a domain D, which is a union of grid cells, are of order eQ = O(H −2m−d ), eD = O(H −2m ), (5) respectively. Boundary cells, however, require a more subtle procedure because of a possibly rather complicated form of the set Q ∩ D (see Fig. 7). Extending the integrand p by 0 on the complement of D in Q and applying Gauß approximations in a straightforward manner would result in very poor accuracy. Before describing the idea used in the subroutines integrate_2d/3d of the FEMB package, we review some facts about univariate Gauß formulas, adapted to our application (cf. [15,16] for a comprehensive treatment of univariate and multivariate numerical integration). R 1 Matlab⃝ is a registered trademark of The MathWorks, Inc., Natick, MA, USA K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) 7 – Theorem 5.1 (Gauß Formulas for Smooth Integrands). The m-point Gauß formula m  sm hϕ = h γαm ϕ(hταm ) α=1 h ϕ is exact for polynomials of degree <2m. More precisely,  h  2m+1  cm max |ϕ (2m) (t )| ϕ − sm h ϕ ≤ h t ∈[0,h] 0  2 √ π /4 (e /64)m m−2m−1/2 (1 + o(1)) as m → ∞. with cm = for approximating     0 Proof. We combine the classical representation of the error [17],  h2m+1 (m!)4 h ϕ − sm hϕ = 0 ϕ (2m) (τ ) for some τ ∈ [0, h], (2m + 1) ((2m)!)3 with Stirling’s formula [18], m! = √ 2π m (m/e)m (1 + o(1)) as m → ∞. This gives cm = (2π m)2 m4m e−4m (1 + o(1))4 . (2m + 1) (4π m)3/2 (2m)6m e−6m (1 + o(1))3 Noting that 2m + 1 = 2m(1 + o(1)) and (1 + o(1))k = 1 + o(1) for k ∈ Z we obtain the expression for cm after simplification.  The error estimate confirms the extremely rapid convergence of Gauß formulas for smooth integrands. This is also the case for piecewise smooth integrands if we apply the Gauß approximations separately on the appropriate subintervals. For example, if ϕ has a discontinuity at ξ = h/3, then  h k γkm,ξ ϕ(hτkm,ξ ) = h 6 h h h 6 3 3 ϕ(hs1− ) + ϕ(hs1+ ) + ϕ(hs2− ) + ϕ(hs2+ ) √ √ with s1± = (3 ± 3)/18 and s2± = (6 ± 3)/9 is the piecewise two-point Gauß approximation of order 2m = 4 for we ignore discontinuities, Gauß formulas still converge. However, the error decreases, in general, very slowly. h 0 ϕ . If Theorem 5.2 (Gauß Formulas for Piecewise Continuous Integrands). If ϕ is piecewise continuous with a finite number of break points ξi , then     h 0 ϕ− sm h as m → ∞.   ϕ  = h o(1) Proof. Scaling the interval of integration, we may assume that h = 1. Let us first assume that ϕ is continuous. In this case the convergence of sm 1 is well known. The argument is based on the Weierstrass approximation theorem [19] which states that any continuous function can be approximated arbitrarily well by polynomials in the maximum norm. We split the error in the form  1 ϕ− 0 sm 1 ϕ=  1 ϕ− 0  1 pℓ 0  +  1 pℓ − 0 sm 1 ϕ  =: E1 + E2 , where pℓ is a polynomial of degree ≤ ℓ which differs from ϕ in absolute value by less than ε/2. Consequently, |E1 | < ε/2. 1 m For m > ℓ/2, 0 pℓ = sm 1 pℓ by the exactness of Gauß formulas. Hence, since the weights γα are positive and sum to one, ( p − ϕ)| < ε/ 2, completing the proof for continuous ϕ . we also have |E2 | = |sm ℓ 1 Turning to the general case, we may assume by linearity that the integrand ϕ has at most a single discontinuity at a point ξ ∈ (0, 1). As is illustrated in Fig. 6, we split the error in the form  1 0 ϕ − sm 1ϕ =  0 1 ϕ−  1 0 =: E1 + E2 + E3 ,   ϕδ + 0 1    m + sm ϕδ − sm ϕ 1 ϕδ − s1 ϕ 1 δ 8 K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – Fig. 6. Gauß approximation for discontinuous integrands. where ϕδ is a continuous approximation of ϕ which differs from ϕ only in the interval [ξ , ξ + δ] and has no smaller infimum and no larger supremum than ϕ . We claim that |Ek | < ε/3 for sufficiently small δ(ϕ) and sufficiently large m(δ). For the first term we have |E1 | ≤ 2C δ with C = supt ∈[0,1] |ϕ(t )|. To bound the third term, let ψδ be a continuous majorant of |ϕ − ϕδ | with support [ξ − δ, ξ + 2δ] and |ψδ | ≤ 2C . Then    1     m m m m m  γα ψδ (τα ) = s1 ψδ ≤ 2 ψδ ≤ 2 (3δ) 2C γ (ϕδ − ϕ)(τα ) ≤ |E3 | =    α α 0 α 1 for sufficiently large m since sm 1 ψδ → 0 ψδ in view of the continuity of ψδ . Hence, choosing δ < (ε/3)/(12C ), the terms E3 and E1 are <ε/3 in absolute value. The same bound is valid for |E2 | for large enough m by the continuity of ϕδ .  Based on the above theorems, we obtain an algorithm for computing multiple integrals with piecewise smooth integrands. We apply Fubini’s theorem with a recursive approximation of the resulting univariate integrals. Identifying most of the discontinuities (break points), which can occur also in derivatives of the integrand, will, in general, lead to rapid convergence. Algorithm 5.1 (Gauß Integration of Piecewise Smooth Functions). An approximation s of mated error <ε hd is computed via a recursive application of univariate Gauß formulas: h 0 ··· h 0 ϕ(x) dx1 · · · dxd with esti- function s = integrate(ϕ, d, h, ε) determine the break points ξi of h h ψ(t ) = 0 · · · 0 ϕ(x1 , . . . , xd−1 , t ) dx1 · · · dxd−1 for m = mmin : mmax m m define weights hγα,ξ and nodes hτα,ξ for a piecewise Gauß formula, based on the partition determined by ξ for all α if d > 1 m m ψ̃α,ξ = integrate(ϕ(·, hτα,ξ ), d − 1, h, ε/2) else m m ψ̃α,ξ = ϕ(hτα,ξ ) end end  m m s̃m = h α γα,ξ ψ̃α,ξ if m > mmin and |s̃m − s̃m−1 | < ε hd /2 break end end s = s̃m m The recursive strategy (ψα,ξ are integrals of the same type; ϕ(·, t ) denotes the function (x1 , . . . , xd−1 ) → ϕ(x1 , . . . , xd−1 , t )) eliminates much of the geometric complexity of the integration problem. For the application to finite element matrix assembly, ϕ = χ p, where, according to (4), χ is the characteristic function of the domain D and p is a smooth expression involving the weighted B-spline basis functions. With the aid of the algorithm we obtain a formula of the form IQ p =  Q ∩D p=  Q χp ≈  Q Q cαQ p(xα,1 , . . . , xα,d ) α after weights and coordinates of evaluation points have been combined and stored appropriately. (6) K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – 9 Fig. 7. Integration over boundary cells. Theorem 5.2 guarantees the convergence of the algorithm if no upper bound on m is imposed, which is, of course, necessary for practical purposes. Since, in general, most break points will be found, the iteration usually terminates instantly or after very few steps in view of the high accuracy predicted by Theorem 5.1. As a consequence, the number of Gauß points for boundary cells is not excessively large. Let us comment on heuristic error estimation in some more detail. In the for-loop convergence is tested by comparing two consecutive Gauß approximations (degrees m and m − 1). As is customary for numerical iterations, the higher degree is assumed to represent the exact value of the integral: I =  h 0 ···  h 0 ϕ(x) dx1 · · · dxd = sm , where sm denotes the exact piecewise Gauß approximation based on true values of ψ . Taking the errors in the computation m of ψα,ξ into account, the error of s̃m−1 upon termination can be bounded by |I − s̃m−1 | ≤ |I − sm | + |sm − s̃m | + |s̃m − s̃m−1 |  m m m γα,ξ |ψα,ξ − ψ̃α,ξ | + ε hd /2 < εhd < 0+h α since the Gauß weights γα,ξ sum to one and ψ was supposed to be computed with accuracy ε hd−1 /2. Clearly, while the error m of s̃m−1 is estimated, s̃m is returnedas final approximation by the algorithm. For the computation of IQ p = Q χ p with the procedure of the algorithm, the location of the discontinuities ξi , caused by the characteristic function χ , is of considerably more importance than the particular form of the smooth part p of the integrand. Therefore, the programs integrate_2d/3d of the FEMB package determine the Gauß parameters (c Q , xQ ) of the approximation (6) with ϕ = χ as test function in the recursive algorithm (⇔ p = 1) postponing the computation of the actual finite element integrals. The generated parameters are then used in the assembly of the Ritz–Galerkin matrix for weighted B-spline integrands p over the sets D ∩ Q without further testing the accuracy of the results. It remains to discuss how the break points, which are crucial for the algorithm, can be determined with the aid of the implicit representation of the simulation domain: D = {x ∈ (0, 1)3 : wD (x) > 0}. A typical situation is illustrated in Fig. 7. As is usually the case, the break points ξi correspond to intersections with the edges of the cube (d = 3), the sides of the square (d = 2, first level of recursion), or the line segment (d = 1, second level of recursion) under consideration: d=3: d=2: d=1: wD (x1 , x2 , ξi ) = 0 ∧ (x1 , x2 ) ∈ {0, h}2 wD (x1 , ξi , hτ ) = 0 ∧ x1 ∈ {0, h} wD (ξi , hτ ′ , hτ ) = 0 where hτ , hτ ′ are nodes of piecewise Gauß formulas, and we have assumed that the lower left corner of Q is the origin. In all cases the break points are zeros of univariate functions. Their computation is facilitated by using polynomial approximations of the weight function wD and the Matlab function roots. This has the additional advantage that several useful expressions can be precomputed (files parameters_2d/3d.mat in the FEMB package). In the example in Fig. 7, there is only one break point ξ1 in the x3 -direction. Accordingly, if a two-point piecewise 2 Gauß formula is used, 4 nodes hτα,ξ are placed on the interval [0, h] of the x3 -axis. For each of these nodes the x1 - and x2 coordinates of the three-dimensional Gauß points are determined by a recursive application of Algorithm 5.1 with d = 2. This is illustrated on the right part of the figure for x3 = hτ12,ξ . For the two-dimensional cross section, break points are now sought in the x2 -direction (x1 ∈ {0, h}). Again, there exists only one break point for the example under consideration, marked with a large tick on the x2 -axis. The small ticks represent the nodes of the corresponding piecewise Gauß formula. For each of these nodes, a final recursive call of the integration algorithm with d = 1 yields the x1 -coordinates of the three-dimensional Q points xβ marked with dots in the cross section on the right of the figure. 10 K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – Fig. 8. Test domains D1 and D2 . Of course, the intersection pattern can be more complicated which leads to more break points ξi . Fig. 7 is kept simple in order to focus on the key strategy of Algorithm 5.1. For further illustrations of the point selection mechanism see the programs demo_integrate_2d/3d in the FEMB package. Does Algorithm 5.1 work sufficiently well? As mentioned before, the integration strategy for boundary cells, described above, is implemented in the Matlab function int = integrate_3d(w_D,H,PAR) with the structure PAR containing the algorithmic parameters ε , mmin , and mmax . As a first illustration of the performance of Algorithm 5.1 we use the Gauß points and weights generated by this program (for interior grid cells the standard parameters based on mmin nodes are returned) to compute  f, D (0, 1)3 ⊇ D : wD (x) > 0. Already for this rather elementary application, the principle difficulty of integrating over boundary cells is present. In accordance with the standard accuracy of Gauß formulas we set ε = H −2mmin . The upper bound mmax is chosen sufficiently large in order not to affect the termination criterion in Algorithm 5.1. This leads to the following program which we include for convenience of the reader who wishes to experiment with our integrator. function s = demo_integrate(w_D,f,H,m_min,m_max) % s ~ int_{w_D>0, (x,y,z) in (0,1)^3} f % initialize algorithmic parameters PAR.gauss = m_min; PAR.steps = m_max-m_min; PAR.tol_int = H^(-2*PAR.gauss); PAR.fit = PAR.gauss+1; PAR.tol_zero = 100*eps; PAR.info = 2; % compute Gauss parameters int = integrate_3d(w,H,PAR); % sum contributions of Gauss rules from grid cells s = 0; for m1=1:H, for m2=1:H, for m3=1:H x = int.p{m1,m2,m3}(:,1); y = int.p{m1,m2,m3}(:,2); z = int.p{m1,m2,m3}(:,3); s = s + int.c{m1,m2,m3}’*f(x,y,z); end, end, end For example, the program segment w_D = @(x,y,z) 1-x.^2-y.^2-z.^2 f = @(x,y,z) = x.^2+y.^2+z.^2 s = demo_integrate(w_D,f,5,3,4) results in s ≈ π /10. As test domains and test function we choose D1 : spade-shaped domain (smooth) D2 : intersection of spheres (non-smooth) f (x, y, z ) = sin 27xy2 z 3 (cf. Fig. 8).   (smooth) K. Höllig, J. Hörner / Computers and Mathematics with Applications (  Fig. 9. Errors of mmin = 1, 2, 3, 4. D1 f (left) and  D2 ) – 11 f (right) as a function of the grid width 1/H; the symbols ×, ⃝, △,  correspond to the Gauß approximations with Fig. 9 illustrates the accuracy of the integration scheme. We have plotted the logarithm of the error e as a function of H for mmin = 1, 2, 3, 4. As expected, the smoothness of D is crucial for convergence. For the smooth domain, a higher order results in a better convergence rate, leading to a relative error smaller than 10−10 for H = 24 and mmin = 4. For the non-smooth domain the integrals still converge, but at nearly the same rate. It should be noted that the programs integrate_2d/3d are not optimized for integrating a single function. Instead, with regard to finite element matrix assembly, many smooth functions are integrated on each grid cell using the same weights and nodes. We now turn to the application to weighted finite elements. Continuing to use the elasticity problem (1) as a test case, with the set D defined in (3) and shown on the left of Fig. 5, we first verify empirically that the boundary integration does not significantly fall behind the order of the Gauß formulas. To this end we consider the integral ID =  dx1 dx2 dx3 D  √ 9π = 5π arctan(5) − √ arctan( 10)  2 2 (x1 − 1/2)2 + (x2 − 1/2)2 + x3 10 as another test case. As before, we set ε = H −2mmin and choose mmax sufficiently large. Assuming that the error of the numerical approximation of ID satisfies eD ≈ cH −α , we plotted the values of the exponent α , estimated from consecutive approximations with different grid widths 1/H, in the left diagram of Fig. 10, using H = 3, 6, . . . , 36 and mmin = 1, 2, 3, 4. The results agree reasonably well with the standard rate α = 2mmin , confirming that the boundary integration strategy achieves the desired accuracy. Actually, the integration over boundary cells is often more accurate than over inner cells, since at least mmin + 1 Gauß nodes are used in each coordinate direction and subdivisions are performed. This explains the deterioration of the rate in the first few steps (H = 6, 9, 12); initially, the percentage of boundary cells is large and their, potentially by several orders superior, error behavior dominates. For smaller grid width 1/H, the majority of grid cells Q does not intersect the boundary, and the order 2mmin of the minimal Gauß rules, used on these cells, will eventually prevail. The for-loop in Algorithm 5.1 could lead to a substantial increase of the number of Gauß points. However, in the example considered, this is not the case. The average value of m upon termination of the for-loop is in all cases less than mmin + 2. As a second, more significant test, the right diagram in Fig. 10 shows a logarithmic plot (base 10) of the relative error of the Ritz–Galerkin matrices, computed with the program assemble_3de based on the points and weights generated by integrate_3d. In order to maintain optimal accuracy of the numerical approximation uh , the order 2mmin of Gauß integration has to be adapted to the degree n of the finite element basis functions. For conventional finite elements, the optimal (minimal) order was determined by Strang in his classical analysis of variational crimes [20]. An analogue of his results for weighted B-splines does not yet exist. However, our experience with the FEMB programs has shown that mmin = n + 1 is an adequate choice. The tolerance ε is set to H −2mmin −1 , one order higher than the standard Gauß rate, as an additional precaution. With these choices, we generated the 8-dimensional arrays G for degrees n = 1, 2, 3, 4 and H = 6, 9, . . . , 24. The relative error of the vector of entries eG = ∥1G(:)∥∞ /∥Gexact (:)∥∞ was estimated using higher Gauß orders. The right diagram in Fig. 10 confirms the very good accuracy of the numerical integration procedure; eG ranges from 2.56 e-05 (degree 1, grid width 1/6) to 4.26 e-12 (n = 4, h = 1/24). These are only a few examples. The finite element approximations considered in Section 7 also serve as illustrations; every algorithm can only be as good as its weakest component. The reader can easily generate his own examples by modifying the demos provided by the FEMB package. It should be noted that the convergence of the integration procedure can deteriorate if ∂ D has edges and corners in the interior of grid cells. For the resulting complicated boundary configurations Algorithm 5.1 will not find all break points. While convergence is guaranteed by Theorem 5.2, the tolerance might not be met since the maximal Gauß order is limited. Such cases are asymptotically rare, and additional effort could be spent to resolve critical situations more accurately without 12 K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – Fig. 10. Accuracy of Algorithm 5.1 for computing integrals over D (left) and the entries of the Ritz–Galerkin matrix (right); the symbols ×, ⃝, △,  correspond to the Gauß orders mmin = 1, 2, 3, 4 (left) and the B-spline degrees n = 1, 2, 3, 4 (right). significantly increasing the overall computing time. In any case, as is to be expected for a complicated algorithm, a definite a priori assessment cannot be made. However, as mentioned before, for B-spline elements we can compute the pointwise residual of the partial differential equation. This provides an a posteriori error measure without having to resort to numerical error estimation. If the result is not satisfactory, the parameters ε , mmin , and mmax can be adjusted and the finite element approximation recomputed. 6. Ritz–Galerkin system The assembly of Ritz–Galerkin systems is straightforward. However, the regular data structure of uniform B-splines can be exploited, simplifying the standard procedures for mesh-based elements. For the model problem (1), which we continue to use as a typical example, the finite element equations have the form 3  k′ ∈K ν ′ =1 a(w bk eν , w bk′ eν ′ ) uk′ ,ν ′ =  fν w bk , D k ∈ K , ν ∈ {1, 2, 3}, where k = (k1 , k2 , k3 ) and K is a rectangular array, as explained in Section 4. It would be a crime against uniform B-splines to label the unknowns uk,ν with a single index and to resort to standard linear solvers. Instead, it is imperative to preserve the regular data structure leading to an exceptionally simple code. If D is a subset of (0, 1)3 , there are H 3 (H = 1/h) grid cells and each of the three components uhν of the finite element approximation is a linear combination of (H + n)3 weighted B-splines w bk . Hence, it is natural to store the unknowns and the right hand side in arrays U and F of dimension (H + n)3 × 3. Accordingly, the Ritz–Galerkin matrix is represented by an 8-foldly indexed array G. Its dimension can be reduced by observing that each B-spline bk shares common support with exactly (2n + 1)3 B-splines bk′ . Hence, we can replace k′ by an index s = k′ − k ∈ S := {−n, . . . , n}3 , representing the offset between the two supports. This banded storage scheme yields an array G of dimension (H + n)3 × (2n + 1)3 × 32 . With these definitions, the Ritz–Galerkin system has the form 3  s∈S ν ′ =1 G(k1 , k2 , k3 , s1 , s2 , s3 , ν, ν ′ )U (k1 + s1 , k2 + s2 , k3 + s3 , ν ′ ) = F (k1 , k2 , k3 , ν), (7) where k ∈ K , ν ∈ {1, 2, 3}. For simplicity of exposition and consistency with the theory, we allowed negative array indices and use the convention that array elements are set to 0 if the indices are out of range. The necessary modifications for Matlab code are trivial, but can be a programmer’s nightmare. After these preliminaries, we can describe the assembly procedure for the Ritz–Galerkin matrix. For each grid cell Q = ℓ/H + (0, 1)3 /H and each pair of B-splines bk , bk′ , kα , k′α ∈ {ℓα − n, . . . , ℓα }, α = 1, 2, 3, with supports overlapping Q we add the contributions of the bilinear form  Q ∩D σ (w bk eν ) : ε(w bk′ eν ′ ), ν, ν ′ = 1, 2, 3, to the appropriate entry of the array G. K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – 13 Algorithm 6.1 (Matrix Assembly for Linear Elasticity). The entries of the Ritz–Galerkin array G are computed by adding the contributions to the finite element integrals from each grid cell: initialize G with zeros for ℓ1 , ℓ2 , ℓ3 = 0 : H − 1 Q = ℓ/H + (0, 1)3 /H for k1 , k′1 = ℓ1 − n : ℓ1 , k2 , k′2 = ℓ2 − n : ℓ2 , k3 , k′3 = ℓ3 − n : ℓ3 s = k′ − k for ν, ν ′ = 1 : 3 add IQ p with p = σ (w bk eν ) : ε(w bk′ eν ′ ) to G(k1 , k2 , k3 , s1 , s2 , s3 , ν, ν ′ ) end end end set diagonal entries G(k1 , k2 , k3 , 0, 0, 0, ν, ν) for irrelevant k to 1 If a subroutine for evaluating the integrand p is supplied, there is only a single nontrivial statement (a weighted sum of function values according to the Gauß approximation (6)) inside nested loops of lengths H 3 , (n + 1)2 · (n + 1)2 · (n + 1)2 , 32 , respectively. Obviously, the B-spline matrix assembly is ideally suited for vectorization, as was demonstrated for Poisson’s equation in [21]. The modification of G in the last program statement assures that the Ritz–Galerkin system is non-singular despite irrelevant B-splines. Since the entries of the right hand side are zero for irrelevant k, irrelevant entries of the solution U are 0, in accordance with our convention (cf. Definition 2.2). The right hand side of the Ritz–Galerkin system, represented by an array F of dimension (H + n)3 × 3, is assembled in a similar fashion. Both arrays, G and F , are returned by the program assemble_3de. Having refrained from converting the Ritz–Galerkin system to standard format, we must adapt iterative solvers to the Bspline data structure. The FEMB package uses a preconditioned conjugate gradient iteration, implemented for the elasticity problem in the program solve_3de. Examining what is required, we need routines for scalar multiplication, addition and subtraction, scalar product, and matrix/vector multiplication. Only the last operation, i.e., the computation of the left side of (7) is not completely straightforward. Generalizing the familiar procedure for multiplication with banded matrices leads to the following algorithm. Algorithm 6.2 (Multiplication by the Ritz–Galerkin Array G). The product F of G with an array of unknowns U is computed by adding the contributions from the matrix diagonals corresponding to the offsets s: initialize F with zeros for s ∈ S for ν = 1 : 3 for k ∈ K update F (k1 , k2 , k3 , ν) by adding 3 ′ ′ ν ′ =1 G(k1 , k2 , k3 , s1 , s2 , s3 , ν, ν )U (k1 + s1 , k2 + s2 , k3 + s3 , ν ) where undefined entries of U are set to 0 end end end There are more sophisticated linear solvers. However, Fig. 13 confirms that the overall solution time is, at least up to degree 3, dominated by the programs integrate_3d and assemble_3de. Hence, an elementary preconditioned conjugate gradient algorithm is adequate for the problems under consideration. 7. Convergence rates and computing times All statements made about the components of the finite element algorithm are of relatively minor importance compared to the assessment of the accuracy of the approximation uh and the efficiency of its computation. We consider these criteria for the model problem (1) which served as a test case already in Section 5. In all computations we use the default values for the integration parameters: mmin = n + 1, ε = H −2mmin −1 . The top two diagrams in Fig. 11 show the logarithm of the L2 - and H1 -errors as a function of the grid width h = 1/H for H = 6, 9, . . . , 24 and degrees n = 1, 2, 3. As usual, the errors are estimated via comparison with solutions of higher degree (which restricts the size of admissible H due to computer limitations). 14 K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – Fig. 11. Logarithmic (base 10) L2 - (top left) and H1 -errors (top right) as well as the corresponding estimated convergence rates for the finite element approximation of the elasticity equations; the markers ×, ⃝, △ correspond to the B-spline degrees n = 1, 2, 3. According to the standard finite element theory, we expect that eL2 ≈ cH −n−1 , eH1 ≈ cH −n for smooth solutions. These approximation orders can deteriorate due to singularities caused, e.g., by edges and corners in the domain. To minimize such effects in the relevant low order partial derivatives, we choose a different profile R for the outer boundary of the rotating solid in (3); the outer radius R is replaced by R(x3 ) = (1 + 3 sin2 (π x3 ))/5 in definition (2) of wD . Then, the rates in the bottom diagrams in Fig. 11, estimated from consecutive errors, are in good agreement with the theoretical predictions for degree n = 1, 2. For cubic B-splines, the error behavior is slightly irregular, which suggests that the grid width has not yet reached the asymptotic range. As mentioned before, in view of the smoothness of the weighted B-spline basis functions, we can substitute the finite element approximation into the partial differential equation for degree n > 1. Fig. 12 shows the relative residual for grid widths 1/H = 1/6, . . . , 1/24. While the residuals decay, estimates of the rates exhibit an irregular behavior (e.g., the residual norm increases from 7.83 e-03 to 7.08 e-02 for degree 4 when H is changed from 18 to 21) which we attribute to numerical instabilities. Approximation theory suggests that eR ≈ cH −n+1 . However, to our knowledge, this has not been rigorously proved for finite element approximations. Since several numerical approximations are involved (integration, conjugate gradient iteration, rate and error estimation) and termination criteria rely on asymptotic considerations, the cause of the oscillatory behavior for degree 4 is hard to tell. We think that the most likely reason are instabilities due to huge condition numbers of the Ritz–Galerkin matrix G, in particular for small grid width and high B-spline degree. The extension procedure, introduced in [1] reduces cond G to the optimal bound O(H 2 ). However, this technique has been implemented so far only in two dimensions, in a software package available at the time of the publication of [4]. As a final illustration, Fig. 13 compares the computing times for the main programs of the FEMB elasticity solver: integrate_3d, assemble_3de, solve_3de. On a 3.3 GHz dual core computer with 8 GB memory the total cpu-time ranges from 18 (n = 1, H = 6) to 106 154 (n = 4, H = 24) seconds. This is quite acceptable in view of the size of the problem. For example, for degree n = 4 and grid width 1/H = 1/24 • Gauß parameters for 243 = 13 824 grid cells have to be determined; • (24 + 4)3 · 93 · 32 = 144 027 072 nonzero entries of the Ritz–Galerkin matrix have to be generated; • a linear system with (24 + 4)3 · 3 = 65 856 unknowns has to be solved. Returning to the statement made at the beginning of this section, the statistics for the example confirm that the use of B-splines pays off for smooth solutions. Very accurate finite element approximations can be obtained already with relatively low-dimensional finite element subspaces. K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – 15 Fig. 12. Logarithm (base 10) of the relative residual eR = ∥Lu − f ∥L2 /∥f ∥L2 for degree n = 2, 3, 4 (markers ⃝, △, ). Fig. 13. Total cpu-times s (in seconds, sum for active cores) for integration (left), system assembly (middle), and conjugate gradient iteration (right). 8. Conclusion Using the equations of linear elasticity as a model problem, we have described the implementation of finite element methods with weighted B-splines using the FEMB program package. A crucial component, the integration over boundary cells, has been discussed in detail, providing the mathematical background for our new algorithm. Moreover, we have explained how to use the regular data structure of uniform B-splines efficiently. Numerical examples illustrate the performance of typical FEMB programs. We note in particular that the new integration routines have considerably improved the efficiency of three-dimensional simulations for weighted approximations; previously a general purpose implementation has been available only for twodimensional applications (cf. [4, Chapter 6]). Summarizing, for basic elliptic model problems, such as the Lamé–Navier system, the FEMB package provides fully automatic highly accurate finite element solvers for arbitrary simulation domains, described in implicit form. There are a number of topics for future research. While the numerical integration algorithm yields satisfactory results in many cases, there is room for improvement. For example, for topologically difficult cases, a combination with subdivision techniques or the utilization of precomputed information seems promising. The latter strategy has already been successfully implemented for linear B-splines on regular triangulations in [21]. Moreover, computing the break points ξi directly, without resorting to polynomial interpolation, is more efficient for special boundary representations such as R-functions, but has not been implemented yet. For the FEMB package an implicit domain description is essential. Converting parametric and polygonal geometry representations to this format leads to interesting research projects. Spline approximation methods will play a crucial role. Optimizing the qualitative properties of weight functions (smoothness, ease of computation) is a closely related issue. Very important is the development of hierarchical techniques for handling problems with singularities. In conjunction with subdivision rules for B-splines, the FEMB package can be used to implement such methods without much additional programming effort [14]. Another area, which is subject of our current research, is free boundary problems and shape optimization. While conventional finite element methods often require remeshing for such applications, weight functions can represent variable simulation domains by varying appropriate parameters. Finally, the regular B-spline data structure is ideally suited for parallelization and vectorization. So far, for weighted approximations, this has been exploited only for a very simple case [21]: a weighted piecewise linear finite element solver for Poisson’s equation. These are just samples of a longer list serving as an advertisement for B-splines, our favorite research topic! 16 K. Höllig, J. Hörner / Computers and Mathematics with Applications ( ) – References [1] K. Höllig, U. Reif, J. Wipper, Weighted extended B-spline approximation of Dirichlet problems, SIAM J. Numer. Anal. 39 (2001) 442–462. [2] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [3] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, Wiley, 2009. [4] K. Höllig, Finite Element Methods with B-Splines, SIAM, 2003. [5] K. Höllig, J. Hörner, A. Hoffacker, Finite element analysis with B-splines: weighted and isogeometric methods, in: J.D. Boissonnat, et al. (Eds.), Curves and Surfaces 2011, in: LNCS, vol. 6920, Springer, 2012, pp. 330–350. [6] K. Höllig, J. Hörner, Finite Element Methods with B-Splines: Supplementary Material. http://www.siam.org/books/fr26/, website, 2012. [7] C. de Boor, A Practical Guide to Splines, Springer, 1978. [8] K. Höllig, J. Hörner, Approximation and Modeling with B-Splines, SIAM, 2013. [9] K. Höllig, J. Hörner, Approximation and Modeling with B-Splines: Supplementary Material, http://www.siam.org/books/ot132/, website, 2014. [10] K. Höllig, U. Reif, J. Wipper, Multigrid methods with WEB-splines, Numer. Math. 91 (2002) 237–256. [11] L.W. Kantorowitsch, W.I. Krylow, Näherungsmethoden der Höheren Analysis, VEB Deutscher Verlag der Wissenschaften, 1956. [12] V.L. Rvachev, T.I. Sheiko, R-functions in boundary value problems in mechanics, Appl. Mech. Rev. 48 (1995) 151–188. [13] C. Bajaj, R. Bettadapura, N. Lei, A. Mollere, C. Peng, A. Rand, Constructing A-spline weight functions for Stable WEB-spline finite element methods, in: Proc. ACM Symposium on Solid and Physical Modeling, 2010, pp. 153–158. [14] C. Apprich, K. Höllig, J. Hörner, A. Keller, E. Nava Yazdani, Hierarchical finite element approximation with B-splines, (submitted for publication). [15] H. Brass, K. Petras, Quadrature Theory: The Theory of Numerical Integration on a Compact Interval, in: Mathematical Surveys and Monographs, vol. 178, AMS, 2011. [16] A.H. Stroud, Approximate Calculation of Multiple Integrals, Prentice-Hall, Englewood Cliffs, NJ, 1971. [17] D. Kahaner, C. Moler, S. Nash, Numerical Methods and Software, Prentice-Hall, 1989. [18] D. Romik, Stirling’s approximation for n!: the ultimate short proof? Amer. Math. Monthly 107 (6) (2000) 556–557. [19] K. Weierstrass, Über die analytische Darstellbarkeit sogenannter willkürlicher Functionen einer reellen Veränderlichen, Sitzungsberichte der Preußischen Akademie der Wissenschaften zu Berlin, 1885, I 633–639 and II 789–805. [20] G. Strang, G. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ, 1973. [21] K. Höllig, J. Hörner, M. Pfeil, Parallel finite element methods with weighted linear B-splines, in: Nagel, Kröner, Resch (Eds.), High Performance Computing, Springer, 2008, pp. 667–676.