Fast formation and assembly of isogeometric Galerkin matrices for trimmed patches
FFast formation and assembly of isogeometricGalerkin matrices for trimmed patches
Abstract
This work explores the application of the fast assembly and formation strategyfrom [8, 17] to trimmed bi-variate parameter spaces. Two concepts for the treatmentof basis functions cut by the trimming curve are investigated: one employs a hy-brid Gauss-point-based approach, and the other computes discontinuous weightedquadrature rules. The concepts’ accuracy and efficiency are examined for the forma-tion of mass matrices and their application to L -projection. Significant speed-upscompared to standard element by element finite element formation are observed.There is no clear preference between the concepts proposed. While the discontinu-ous weighted scheme scales favorably with the degree of the basis, it also requiresadditional effort for computing the quadrature weights. The hybrid Gauss approachdoes not have this overhead, which is determined by the complexity of the trimmingcurve. Hence, it is well-suited for moderate degrees, whereas discontinuous-weighted-quadrature has potential for high degrees, in particular, if the related weights arecomputed in parallel. Isogeometric analysis (IGA) has been introduced to overcome the profound inefficienciesin the conventional interaction of CAD and simulation tools [12, 19]. The ground-breaking idea of IGA is to perform numerical simulations using CAD technology suchas (non-uniform rational) B-splines. During the last decade, it has been demonstratedthat this paradigm provides several computational benefits. Indeed, IGA outperforms a r X i v : . [ c s . C E ] J a n raditional simulations in most academic benchmarks [11, 10, 22, 28] and is nowadaysgenerally recognized as a powerful alternative to the conventional analysis methodology.The straightforward utilization of high-degree and high-continuity basis functions isan outstanding benefit of IGA. Together with the concept of k -refinement, the resultinganalysis features high robustness and accuracy w.r.t. the degrees of freedom employed,see, e.g. [4, 13, 14]. This ability to perform high-order accurate simulations is, however,somewhat limited due to the computational cost, because the formation and assemblyof the system matrix gets more involved with increasing polynomial degree p . Thestate-of-the-art at the core of standard finite element codes is an element-wise assembly.IGA adopted this concept, and the introduction of Bézier extraction [5, 6] has provideda means to map between a smooth spline basis and an element-based representation.Being compatible with conventional analysis routines has played an essential role forthe acceptance and dispersion of IGA in the numerical analysis community and alloweda simple integration of this new paradigm into existing simulation software. Yet, thecorresponding cost for setting up the system matrix for a C p − -continuous d -dimensionaltensor product B-spline basis is O (cid:16) p d (cid:17) per degree of freedom [8]. Consequently, high-order analysis is doomed to be computationally expensive when a conventional matrixformation strategy is employed.Reduced quadrature rules [3, 16, 20, 21, 36] improve the efficiency of the numericalintegration on the element-level. However, the element-wise assembly by itself limits thecomputation cost to O (cid:16) p d (cid:17) [8]. Sum factorization [1, 35] approaches this thresholdby rearranging the computations to exploit the tensor product structure of the B-splinebasis. The resulting cost using element-wise Gaus quadrature is O (cid:16) p d +1 (cid:17) . As shown in[7], sum factorization with Gauss quadrature can even yield a computational complex-ity of O (cid:16) p d +2 (cid:17) , when used globally and not on an element level. Weighted quadrature introduced in [8] is a new integration technique that also drops the element perspective.A quadrature rule is set up for each test function by incorporating the test functioninto the quadrature weights. The outstanding feature of the resulting quadrature is itsindependence on the degree p for spline with maximal smoothness, i.e., C p − . The appli-cation of weighted quadrature to non-uniform spline with mixed continuity is addressedin [17]. Furthermore, the integration of sum factorization and weighted quadrature intoa row-based assembly strategy results in the fast formation and assembly approach de-tailed in [8, 17]. These three ingredients – row assembly, sum factorization, and weightedquadrature – reduce the computational cost to O (cid:16) p d +1 (cid:17) .The treatment of trimmed patches has been denoted as an open challenge for thisfast formation and assembly approach at the recent INdAM Workshop on GeometricChallenges in Isogeometric Analysis. Trimming is an essential concept for representingcomplex geometries with tensor product B-splines. The main component involved is the trimming curve , which is specified in the parameter space and restricts the visible partof the spline object to a subregion. The resulting trimmed space consists of interior,exterior, and cut basis functions. The latter introduces several computational challengesas detailed in [29]. For example, the realization of a high-order accurate numerical inte-2ration scheme is far from trivial even in the classical finite element setting, see e.g. [15].Moreover, cut basis functions do not follow the tensor product structure of the basisanymore. In order words, trimming violates an essential property for the fast formationprocess. Hence, the question arises if a significant reduction of the computational costsis restricted to the non-trimmed case.This paper provides extensions to the fast formation and assembly approach presentedin [8, 17] that allow the analysis of trimmed domains. In particular, the integrationover cut basis functions is investigated using either a (i) Gauss quadrature or a (ii)weighted quadrature approach. Their performance and accuracy are compared for themass matrix formation of trimmed bi-variate spaces with different complexity. Bothconcepts maintain the optimal approximation order and significantly improve efficiencycompared to a standard element based Gaussian assembly. From a conceptional pointof view, the extension to tri-variate spaces is straightforward, but the implementationgets more involved due to the increased topological complexity of the cut elements.The paper is structured as follows: Sect. 2.1 outlines the fast formation and trimmingconcepts and highlights the contradiction of their underlying ideas. The proposed ap-proaches to overcome this barrier are presented in Sect. 3 and they are compared bynumerical experiments in Sect. 4. The paper closes with concluding remarks. This section provides short introductions to weighted quadrature, sum factorization,and trimmed patches, which are the preliminaries for the proposed approach detailedsubsequently. First, the main aspects presented in [8, 17] are recapitulated. Then, thefundamentals of trimmed spaces and the challenge of deriving a fast formation techniquefor them is outlined.
Let us recapitulate some essential properties of the basis functions used. A
B-spline B i,p is described by piecewise polynomial segments of degree p . The continuity betweenthem is specified by the knot vector Ξ , which is a non-decreasing sequence of parametriccoordinates ξ j (cid:54) ξ j +1 . The values of these knots ξ j define the location where adjacentsegments join. The continuity at these breakpoints is C p − m , with m denoting the mul-tiplicity of the corresponding knot value, i.e., ξ j = ξ j +1 = · · · = ξ j + m − . The knotspan s refers to the half-open interval [ ξ s , ξ s +1 ), and if its size is non-zero, it marks anelement. Furthermore, the knot vector Ξ defines an entire set of linearly independent B-splines { B i,p } ni =0 on the parametric domain Ω. Each B i,p has local support, supp { B i,p } ,specified by the knots { ξ i , . . . , ξ i + p +1 } , and each knot span s contains p + 1 non-zeroB-splines. Bi-variate basis functions are obtained by computing the tensor product ofuni-variate B-splines B i ,p and B i ,p of degrees p and p , which are defined by separateknot vectors Ξ and Ξ for the parametric directions ξ and ξ , respectively. This can3e generally expressed as B i ( ξ ) = B i ,...,i d ( ξ , . . . , ξ d ) = d Y k =1 B i k ( ξ k ) . (1)In this paper, the focus lies on weighted quadrature for the formation of mass matrices.Hence, the aim is to derive an efficient and accurate evaluation of the following uni-variateintegral Z Ω B i ( ξ ) B j ( ξ ) c ( ξ ) dξ. (2) B i ( ξ ) and B j ( ξ ) are test and trial functions and c ( ξ ) is determined by the geometrymapping. In general, numerical quadrature rules Q are designed to be exact for the casethat c ( ξ ) = 1, and they provide quadrature weights w k , which allow the expression ofequation (2) as a sum over corresponding quadrature points x k , i.e., Q = X k B i ( x k ) B j ( x k ) w k := Z Ω B i ( ξ ) B j ( ξ ) dξ. (3)The novelty of weighted quadrature is that a quadrature rule Q i is designed for eachtest function B i by incorporating the test function into the quadrature weights Q i = X k B j ( x k ) w k,i := Z Ω B j ( ξ ) ( B i ( ξ ) dξ ) . (4)The computation of the weights w k,i requires an adequate distribution of the quadraturepoints x k . Here, the procedure proposed in [17] is employed, and the reader is referredto this publication for details. It determines the minimal number of quadrature pointsfor each element of a non-uniform knot vectors with mixed continuity. Subsequently, auniform distribution within the elements’ interior is employed so that no x k coincideswith a knot ξ i .Once the position of the quadrature points x k is fixed, the corresponding weights canbe computed by the following system of equations Q i ( B j ) = X k ∈Q i B j ( x k ) w k,i := Z Ω B j ( ξ ) ( B i ( ξ ) dξ )... ... Q i ( B j n ) = X k ∈Q i B j n ( x k ) w k,i := Z Ω B j n ( ξ ) ( B i ( ξ ) dξ ) (5)where j , . . . , j n are the indices of all trial functions whose support overlaps with theone of the current test function, supp { B i } , and the index set Q i refers to all quadraturepoints that lie within supp { B i } . Solving (5) by QR-factorization yields w k,i for ∀ k ∈ Q i ,4nd the remaining quadrature weights are set to zero. The solution of (5) maybe notunique, and the system of equations can be weighted to improve positivity and bound-edness of the weights as detailed in [17]. In the case of multi-variate tensor productbasis functions, quadrature rules Q i , . . . , Q i d are computed for each parametric direc-tion. Fig. 1 shows the quadrature layout for two bases of degree 2 and 6, respectively.Note that the number of quadrature points for the inner elements does not change. (a) Bi-degree p = 2 (b) Bi-degree p = 6 Figure 1: Weighted quadrature points for two bi-variate bases with the same number ofelements but with a different polynomial degree.
Sum factorization takes advantage of the tensor product structure of the quadraturerule and the spline basis by a reordering of the numerical operations such that onlyuni-variate quadrature rules are required. Consider the computation of an entry of amass matrix for a bi-variate basis m i , j = Z Ω B i ( ξ ) B j ( ξ ) c ( ξ ) d ξ . (6)Starting from (1), equation (6) can be expressed as a recursion of uni-variate integrals m i , j = Z Ω B i ( ξ ) B j ( ξ ) × Z Ω B i ( ξ ) B j ( ξ ) c ( ξ , ξ ) dξ dξ . (7)Consequently, it is straightforward to apply the weighted quadrature rules Q i , Q i asso-ciated with a bi-variate test function B i . Finally, it is noted that the equations have tobe pulled back to the parametric domain [1, 17] to exploit the structure of equation (7).5 .2 Trimmed patches Tensor product B-splines allow full control over the continuity and degree of the basis, butthey possess an intrinsic four-sided topology limiting their ability to represent arbitrarydomains. Trimming provides a remedy to this restriction by defining the valid area Ω v independent from the basis’ structure. In particular, curves in the parametric spacespecify Ω v , and the curve’s orientation determines the interior and exterior domain.Usually, these trimming curves C t are represented by B-splines as well [29], but thischoice is more or less free.The presence of C t divides the basis into three different functions types based on theoverlap of the support with the valid domain, i.e., S v i := supp { B i } ∩ Ω v . That is, aB-spline B i is classified as:• Exterior if S v i = ∅ ,• Interior if S v i = supp { B i } ,• Cut if 0 < |S v i | < | supp { B i }| ,where |·| denotes the Lebesgue measure in R d . Fig. 2 illustrated these different classesfor a bi-variate cubic basis. In the context of analysis, exterior basis functions can be ne-glected from the system of equations, and interior ones can be treated as usual. Cut basisfunctions, however, induce profound numerical challenges regarding the application ofboundary conditions, the conditioning of system matrices, and accurate integration. Ad-dressing these aspects is far from trivial, even in the conventional low-order finite element Ω v C t supp { B , } supp { B , } supp { B , } Figure 2: Trimmed cubic bi-variate basis with the trimming curve C t specifying the validdomain Ω v (gray). The resulting B-splines types are interior (green), cut (red),or exterior (yellow) based on the overlap of the support, supp { B i ,i } , with Ω v .6etting. The loss of the tensor product structure and the requirement of higher-orderaccuracy complicates the situation further for the fast formation outlined in Sect. 2.1.In the following, the focus is on integrating cut basis functions. The application ofboundary conditions is not an issue in this paper because only the formation of massmatrices and L -projection is considered. The extended B-spline concept is employedto address the conditioning aspect. This procedure is independent of the assembly andformation process, and hence, it is not described here. The interested reader is referredto [18, 31, 30]. A trimming curve C t introduces an arbitrarily located jump discontinuity within theparameter space. Simply integrating over this interface or neglecting quadrature pointsthat lie outside of the valid domain Ω v will evidently lead to incorrect results. Conse-quently, the numerical integration scheme has to account for these arbitrarily locateddiscontinuities. In the case of weighted quadrature, this circumstance affects not onlythe correct representation of the integration domain S v i of a cut B-spline B i but also thecomputation of its quadrature rules Q i and Q i . Moreover, sum factorization cannotbe applied because S v i does not follow a tensor product structure in general.However, we can split the domain S v i into a regular part S r i , which follows the tensorproduct structure (at least on the element-level), and a trimmed part S t i , which consistsof all elements cut by the trimming curve. The integral over a cut basis function can bewritten as Z S v i B i ( ξ ) d ξ = Z S r i B i ( ξ ) d ξ + Z S t i B i ( ξ ) d ξ . (8)In the following, the numerical integration of S t i employs a standard element-wise as-sembly procedure and Sect. 3.1 lists various concepts for the distribution of quadraturepoints. Furthermore, this work investigates two options for the treatment of the re-maining regular part S r i : In Sect. 3.2, it is integrated using Gaussian quadrature, whileSect. 3.3 derives a weighted quadrature rule that acknowledges the presence of the trim-ming curve within the support. The fast formation and assembly approach opens the path towards efficient higher-order simulation. Hence, it is of utmost importance that the analysis, and therefore theintegration of cut elements, is performed with higher-order accuracy. Luckily, there is asubstantial body of literature on this topic because integrating over elements cut by anarbitrary interface is a canonical problem in various analysis schemes such as fictitiousdomain methods, extended finite element approaches, and the simulation with trimmedspline geometries. The proposed concepts can be broadly divided into strategies that(i) set up tailored integration rules or (ii) decompose cut elements into sub-elements,7hich then employ standard quadrature rules. For the latter type, several works [2, 9,15, 23, 25, 24, 27] have demonstrated that the introduction of curved sub-elements withhigh degree, whose edges or faces capture the interface, allows for higher-order accurateintegration with a modest number of integration points. However, generating these sub-elements is a challenge on its own, and its complexity is determined by the dimension ofthe domain and the shape of the interface. An alternative is the construction of tailoredintegration rules for a cut element, see e.g. [32, 33, 34]. These strategies usually obtainintegration weights by solving a system of moment-fitting equations. Utilizing Lasserre’stheorems [26], integrals over the boundaries of the cut elements provide the right-handside required. Müller et al. [33] showed that this concept achieves higher-order accurateintegration over curved three-dimensional domains.Herein, the decomposition-based approach detailed in [15, 37] is utilized. It representsthe valid domain Ω v in the parameter space by a set of higher-order Lagrange elements,which are aligned with the parameter space as well as the trimming curve C t . Further-more, they have the same polynomial degree as the underlying B-spline basis so that therepresentation quality of C t matches the approximation power sought for the analysis.The mappings of these Lagrange elements distribute Gauss quadrature points withincut elements, as shown in Fig. 3. It is worth noting that this particular choice for theintegration of cut elements is not essential for the ideas of this work. Any integrationscheme could be employed, as long as it allows for higher-order accurate integration oftrimmed domains. Therefore, a detailed discussion on this part of the overall integrationof cut B-splines is omitted. (a) Bi-degree p = 2 (b) Bi-degree p = 6 Figure 3: Higher-order accurate distribution of Gaussian quadrature points within cutelements for different degree. 8 .2 Hybrid Gauss approach
The perhaps most straightforward way to integrate over the regular portion S r i of a cutbasis function is to employ standard Gauss quadrature using a conventional element-wise assembly. In contrast to the cut elements of S t i , the evaluations during the for-mation can exploit the tensor product structure. Hence, it is beneficial to implementseparate routines for S r i and S t i . As a result, a cut B-spline is treated by a hybrid Gauss-quadrature-based formation concept. Regarding the entire trimmed space, Gauss pointsare restricted to the vicinity of the trimming curve. Yet, they propagate further into theinterior when the degree increases since the supports of the trimmed B-splines increase.Consider the example shown in Fig. 4, all elements are affected in the degree 6 case.Nevertheless, the portion of the computational cost for setting up the entire system ofequations decreases with the fineness of the parameter space. (a) Degree 2 (b) Degree 6 Figure 4: Gaussian quadrature for the entire valid support S v of all cut basis functions.The support of cut basis functions is highlighted in red. Remark . Since the elements of S r i follow the tensor product structure, sum factorizationwith Gauss points x k may be used by setting the precomputed coefficient c ( x k ) ofequation (6) to zero if x k / ∈ S r i . This option is, however, not considered here. The domain splitting (8) allows us to extract those elements that follow the tensor prod-uct structure. Hence, we can apply sum factorization there, as noted in the previous sec-tion. In the context of weighted quadrature, however, we cannot set the coefficients c ( x k )in equation (7) to zero if the quadrature point x k / ∈ S r i . In contrast to Gauss quadra-ture, the weighted rules take advantage of the continuity between elements. Thus, the9ero coefficients introduce artificial discontinuities within the domain of the quadrature,which leads to an error in the evaluation of the integral.The construction of discontinuous weighted quadrature is proposed to resolve thisproblem. The idea is to incorporate the information on the location ξ disc of these artificialdiscontinuities into the quadrature weights. As a result, the numerical integration onone side of ξ disc can be treated independently from the other. Each ξ disc marks a knotadjacent to a cut element. Thus, the trimming curve determines the number of ξ disc required. The tools for deriving related quadrature rules are the fact that quadrature isa linear operator and the properties of knot insertion.Knot insertion denotes the refinement of a B-spline object by adding knots ˜ ξ into itsknot vector Ξ . This procedure leads to nested spline spaces, and the subdivision matrix S : R n +1 R ˜ n +1 encodes the coefficients of the initial coarse representation and therefined one. Considering quadrature weights, we obtain the relation˜ w i = n X j =0 S ij w j for i = 0 , . . . , ˜ n (9)where ˜ w i refers to quadrature weights computed for the refined basis functions. If onlyone knot is inserted, i.e., ˜ Ξ = Ξ ∪ ˜ ξ , where ˜ ξ ∈ [ ξ s , ξ s +1 ), the non-zero entries of S aredetermined by ( S ( k, k −
1) = 1 − α k S ( k, k ) = α k α k = k (cid:54) s − p ˜ ξ − ξ k ξ k + p − ξ k s − p + 1 (cid:54) k (cid:54) s k (cid:62) s + 1 (10)Multiple knots can be inserted by repeating this process, and the multiplication of theindividual single-knot matrices yields the overall subdivision matrix.Using the relations provided by S , the construction of discontinuities weighted quadra-ture rules for each artificial discontinuities ξ disc goes as follows:1. Knot insertion at ξ disc so that the uni-variate basis become C − continuous there,and storing of the corresponding S .2. Determination of the minimal number of quadrature points n min for the refinedbasis [17], and addition of new nested quadrature points if necessary.3. Computation of ˜ w for basis functions of the refined basis using equation (5).4. Multiplication of ˜ w by S T to obtain w for the initial basis.Fig. 5 illustrates these steps for a single cut basis function. Using nested quadraturepoints in step 2 of the construction allows the reuse of all coefficients set up for theweighted quadrature for interior B-splines during the integration process. The resultingweights w of the initial and nested quadrature points are obtained by a linear combinationof weights ˜ w that account for the discontinuity at ξ disc . Thus, we can set all coefficientson one side of ξ disc to zero without affecting the integral on the other side.10 , a Cut B-spline andits WQ ruleΩ v t ˜ B , ˜ B , ˜ B , ˜ B , b Refined basis andadditional nestedquadrature points0000 1 2 3333 4 5555 ξ disc B , c Initial cut B-splineand its DWQ rule0000 1 2 3 4 5555 S r3 S t3 B , d Quadrature forthe entire trimmedsupport S v3 Gauss pointDWQ point S r3 S t3 Figure 5: Discontinuous weighted quadrature (DWQ) for the cubic B-spline B , cut atthe trimming position t . (a) Conventional weighted quadrature (WQ) points(black dots) for B , . (b) Quadrature layout with additional quadrature points(white) for WQ of the refined discontinuous ˜ B j, associated to B , . (c) Linearcombination of the refined discontinuous WQ rules to obtain the DWQ for B , . In (a,c), the points’ height indicates the related weight value. Remark . The nested quadrature points are uniformly distributed between the initialones. The number of added points per interval is chosen such that the distance betweenbetween them gets as large as possible.In the context of sum factorization, only discontinuous weighted quadrature pointswithin S r i are of interest. In a general trimming situation, these points do not cover allelements of S r i . The remaining parts are integrated by Gaussian quadrature. Considerthe example shown in Fig. 6. Note that the artificial discontinuity restricts the Gausspoints to the vicinity of the trimming curve, independent of the degree. Remark . Every artificial discontinuity ξ disc increases the number of required quadra-ture points in the elements adjacent to ξ disc . Hence, the advantage of providing a dis-continuous weighted quadrature rule for a cut basis functions decreases with the numberof ξ disc . In the present implementation, an individual quadrature is constructed only ifnot more than one ξ disc for each parametric direction is present in the basis function’ssupport. 11 disc2 DWQGaussianquadrature
Figure 6: All quadrature points of a cut B-spline of bi-degree 6: above the artificial dis-continuity ξ disc2 discontinuous weighted quadrature (DWQ) is employed, whileGauss points are used below. The white dots mark the nested quadraturepoints added during the construction. The numerical experiments focus on the formation of the mass matrix and L -projectionfor trimmed bi-variate B-splines. Fig. 7 illustrates the three investigated trimming casesof different complexity.In the following, the standard element-wise Gauss procedure provides the referencesolutions, and the fast formation and assembly strategies considered are:• Weighted quadrature (Sect. 2.1) for the regular support S r of cut and interior B-splines.• Hybrid Gauss (Sect. 3.2) for the regular support S r of cut B-splines and weightedquadrature for interior B-splines.• Discontinuous weighted quadrature (Sect. 3.3) for the regular support S r of cutB-splines and weighted quadrature for interior B-splines.Cut elements are integrated by the same element-wise higher-order accurate procedure(Sect. 3.1) in all cases. The related decomposition by Lagrange elements for the distri-bution of the integration point goes up to degree 6. Thus, this degree marks the upperthreshold for the following numerical experiments.12 a) Line (b) Circle (c) Curvy corner cut Figure 7: The different trimming cases of the numerical examples.
The first numerical example addresses the accuracy of the formation of the mass ma-trix. Therefore, the line-trim-case (Fig. 7a) is considered for a 10 by 10 element basis.Standard Gauss provides the reference matrix, and the deviation of the other matricesbuild by the fast concepts is measured in the Euclidean norm. Table 1 summarized theresults for various degrees. The results of the other trimming cases are omitted becausethey lead to the same findings and do not provide further insights.Table 1: The error of the computed mass matrix for the line-trim-case (Fig. 7a) usingdifferent fast formation and assembly concepts.Degree Hybrid Gauss Disc. weightedquadrature Weightedquadrature1 4.68883379e-17 4.69071638e-17 3.33078900e-032 2.18894391e-17 2.19723138e-17 3.38923479e-043 4.98024622e-17 4.98024919e-17 3.03838755e-044 3.80411220e-16 3.84739261e-16 1.78250728e-045 1.01851172e-15 1.01851512e-15 1.59974261e-046 2.00575179e-15 2.00590425e-15 1.28293909e-04It is apparent that the direct application of weighted quadrature to cut basis functionsis not able to derive the correct mass matrix, whereas the proposed concepts are accurateup to machine precision. Thus, the remaining numerical examples will only employ thehybrid Gauss and discontinuous weighted quadrature approaches. L -projection Here, the accuracy and efficiency of the hybrid Gauss approach and the discontinuousweighted quadrature are investigated for all three trimming cases shown in Fig. 7. One13ngredient for the quality of the results is the conditioning of the mass matrix, which maysuffer due to the presents of trimmed basis functions. The extended B-spline conceptis employed for all simulations to guarantee well-conditioned mass matrices. For detailson this approach, the interested reader is referred to [18, 30, 31, 37]. L -projection is used to fit the trimmed patch to the target function f ( x, y ) = sin(2 x ) cos(3 y ),and the quality of the resulting approximation is measured in the relative L -error norm k (cid:15) rel k L . The convergence results for the line-, circle-, and curvy corner cut trimmingcases are illustrated in Figs. 8 to 10, respectively. In all cases, optimal convergence ratesand excellent agreement to the reference solutions can be observed for both approachesproposed. − . − − . − − − − h − k (cid:15) r e l k L Hybrid Gauss − . − − . − − − − h − Disc. weighted quadrature
Standard Gauss p = 1 Standard Gauss p = 2 Standard Gauss p = 3Standard Gauss p = 4 Standard Gauss p = 5 Standard Gauss p = 6Fast-Ass. p = 1 Fast-Ass. p = 2 Fast-Ass. p = 3Fast-Ass. p = 4 Fast-Ass. p = 5 Fast-Ass. p = 6 Figure 8: Line convergence study14 − − . − − − − h − k (cid:15) r e l k L Hybrid Gauss − − . − − − − h − Disc. weighted quadrature
Standard Gauss p = 1 Standard Gauss p = 2 Standard Gauss p = 3Standard Gauss p = 4 Standard Gauss p = 5 Standard Gauss p = 6Fast-Ass. p = 1 Fast-Ass. p = 2 Fast-Ass. p = 3Fast-Ass. p = 4 Fast-Ass. p = 5 Fast-Ass. p = 6 Figure 9: Circle convergence study − . − − . − − − − h − k (cid:15) r e l k L Hybrid Gauss − . − − . − − − − h − Disc. weighted quadrature
Standard Gauss p = 1 Standard Gauss p = 2 Standard Gauss p = 3Standard Gauss p = 4 Standard Gauss p = 5 Standard Gauss p = 6Fast-Ass. p = 1 Fast-Ass. p = 2 Fast-Ass. p = 3Fast-Ass. p = 4 Fast-Ass. p = 5 Fast-Ass. p = 6 Figure 10: Curvy corner cut convergence study15 .2.2 Efficiency comparison
The assessment of the efficiency of the different formation and assembly strategies isbased on the timings of the finest discretizations used in the convergence studies. All rou-tines of the proposed concepts have been implemented in an in-house MATLAB® code,which does not utilize any parallelization capabilities. Total timings have been obtainedon an Intel® Core™ i7-8700 3.2 GHz processor with 16 GB RAM. The timings havebeen measured with MATLAB’s tic-toc command.Fig. 11 illustrates the total formation and assembly time obtained. The standardGauss formation is the most expensive concept, as expected. The performance of theother approaches, however, is a bit astonishing since the hybrid Gauss approach is moreefficient in most situations. In order to gain deeper insight, the various components ofthe total formation and assembly time are compared for each trimming case in Fig. 12. p t i m e s Hybrid Gauss Disc. weighted quadrature Standard Gauss(a) Line2 4 610 p t i m e s (b) Circle 2 4 610 p t i m e s (c) Curvy corner cut2 4 610 p t i m e s Hybrid Gauss Disc. weighted quadrature Standard Gauss
Figure 11: Total formation and assembly time including (i) the computation of quadra-ture weights and (ii) the integration of the entire valid support of cut and inte-rior basis functions. Each sub-figure corresponds to a trimming case (Fig. 7),while the different graphs refer to the assembly and formation strategy used.16 − − t i m e s Hybrid Gauss − − Disc. weighted quadrature − − t i m e s Hybrid Gauss − − Disc. weighted quadrature − − p t i m e s Hybrid Gauss − − p Disc. weighted quadrature
Quad. weights Interior B-splines S r of cut B-splines Cut elements a Lineb Circlec Curvy corner cutFigure 12: The various components of the total formation and assembly time for eachtrimming case (Fig. 7). The graphs refer to the set-up of the quadratureweights (Quad. weights), as well as the integration of the regular (InteriorB-splines) and trimmed ( S r of cut B-splines, Cut elements) B-splines.17hese components are (i) the computation of the (normal and discontinuous) weightedquadrature rules, (ii) the formation of interior B-splines, the integration over (iii) the reg-ular part S r of cut B-splines and (iv) the remaining portion S t covered by cut elements.The components (ii) and (iv) employ the same routines in both approaches. Hence,there are only minor deviations due to the usual time irregularities between two runs ofthe same code. Looking at the timings related to the regular part S r of cut B-splinesreveals that the discontinuous weighted quadrature scales better w.r.t. the degree. Thecomputation of the corresponding quadrature weights, however, is more expensive dueto the additional expenses for setting up the discontinuous weighted quadrature rules forcut B-splines. This overhead varies with the trimming situation, which determines thenumber of artificial discontinuous introduced, and therefore, the number of additionalquadrature rule computations. It has been demonstrated that the fast assembly and formation strategy presented in[8, 17] cannot be applied directly to trimmed spline spaces. This work examines conceptsto overcome this limitation. In particular, B-splines B i cut by the trimming curve C t require an adaptation of the fast procedure. Therefore, their support is divided into aregular part S r i and a trimmed part S t i . The latter consists of the elements intersected by C t . These elements do not follow the tensor product structure of the underlying splinebasis and require tailored integration rules. There are different integration techniquesfor this task available in the literature. This work employs the approach presented in[15, 37], but any scheme that allows higher-order accurate integration (favorably in anefficient manner) may be used. This property, however, is essential. Otherwise, theintegration of the cut elements outweighs the benefits of applying fast assembly andformation of the remaining basis.The formation of the remaining contributions of a cut B-spline B i involves the integra-tion over its S r i , i.e., all non-cut elements of B i enclosed by C t . This paper presents twoconcepts: The first one is a hybrid Gauss approach which employs standard Gaussianquadrature. Hence, the formation of a cut B-spline utilizes standard Gauss quadrature inall elements, but the evaluations for S r i can be implemented more efficiently. The secondconcept introduces tailored weighted quadrature rules for a cut B-spline’s S r i . These dis-continuous weighted quadratures define artificial discontinuous ξ disc between elements forthe construction of the quadrature weights. The trimming curve C t determines the lo-cations of ξ disc . To be precise, they mark the knots adjacent to cut elements. Comparedto normal weighted quadrature, additional quadrature points are needed next to ξ disc tofulfill the exactness condition. The resulting discontinuous quadrature rules allow thatthe coefficients on one side of ξ disc can be set to zero without affecting the contributionsof the coefficients on the other side. This property enables a straightforward applicationof sum factorization to a subregion of a support. Finally, it is emphasized that interiorB-splines, i.e., those not cut by C t , are treated by normal weighted quadrature and sumfactorization in both approaches. 18he approximation quality and efficiency of the proposed approaches have been in-vestigated by the formation of mass matrices and L -projection for bi-variate trimmedspline spaces. It has been shown that the hybrid Gauss approach and discontinuousweighted quadrature yield the same accuracy as the reference solutions obtained bystandard element-wise formation using Gauss quadrature. Both schemes show signifi-cant speedups compared to the reference computations. Using discontinuous weightedquadrature, computing the contribution corresponding to S r i scales better with the de-gree. However, the computational effort for setting up the additional quadrature rulesdominates the total formation and assembly time, especially for complex trimming caseswith multiple ξ disc . Overall, the hybrid Gauss approach turned out as the more efficienttechnique for the majority of the test cases considered. From the numerical experiments,the following conclusions may be drawn: (i) The hybrid Gauss scheme is sufficient formoderate degrees; (ii) discontinuous weighted quadrature rules are beneficial for highdegrees such as p (cid:62)
5, but (iii) the computation of their weights must be performed inparallel.Future work will focus on this parallelization aspect for the computation of the dis-continuous weighted quadrature rules as well as the overall fast assembly and formationscheme. Furthermore, the efficiency of the current hybrid Gauss can be enhanced byutilizing sum factorization. The performance in the case of three-dimensional problemsis another topic worth exploring. From a conceptional point of view, the extension ofthe concepts presented to another parametric dimension is uncomplicated. Consideringdiscontinuous weighted quadrature, adding another dimension merely affects the deter-mination of proper artificial discontinuities. However, the implementation is much moreinvolved due to the increased complexity of the trimming situations. Higher-order ac-curate integration of trimmed tri-variate splines is already a challenge in the context oftraditional finite element setting. Solving this challenge is a prerequisite for applyingfast assembly and formation.
References [1] Antolin, P.; Buffa, A.; Calabrò, F.; Martinelli, M.; Sangalli, G.: Efficient matrixcomputation for tensor-product isogeometric analysis: The use of sum factorization,
Computer Methods in Applied Mechanics and Engineering , 285:817–828, 2015, doi:10.1016/j.cma.2014.12.013.[2] Antolin, P.; Buffa, A.; Martinelli, M.: Isogeometric analysis on V-reps: First results,
Computer Methods in Applied Mechanics and Engineering , 355:976–1002, 2019, doi:10.1016/j.cma.2019.07.015.[3] Auricchio, F.; Calabrò, F.; Hughes, T.J.R.; Reali, A.; Sangalli, G.: A simple algo-rithm for obtaining nearly optimal quadrature rules for nurbs-based isogeometricanalysis,
Computer Methods in Applied Mechanics and Engineering , 249–252:15–27,2012, doi: 10.1016/j.cma.2012.04.014. 194] Beer, G.; Marussig, B.; Duenser, C.
Basis Functions, B-splines , pages 35–71.Springer International Publishing, Cham, 2020. ISBN 978-3-030-23339-6, doi:10.1007/978-3-030-23339-6_3.[5] Borden, M.J.; Scott, M.A.; Evans, J.A.; Hughes, T.J.R.: Isogeometric finite elementdata structures based on Bézier extraction of NURBS,
International Journal forNumerical Methods in Engineering , 87(1-5):15–47, 2011, doi: 10.1002/nme.2968.[6] de Borst, R.; Chen, L.: The role of Bézier extraction in adaptive isogeometricanalysis: Local refinement and hierarchical refinement,
International Journal forNumerical Methods in Engineering , 2017, doi: 10.1002/nme.5696.[7] Bressan, A.; Takacs, S.: Sum factorization techniques in isogeometric analysis,
Computer Methods in Applied Mechanics and Engineering , 352:437–460, 2019, doi:10.1016/j.cma.2019.04.031.[8] Calabrò, F.; Sangalli, G.; Tani, M.: Fast formation of isogeometric Galerkin ma-trices by weighted quadrature,
Computer Methods in Applied Mechanics and Engi-neering , 316:606–622, 2017, doi: 10.1016/j.cma.2016.09.013.[9] Cheng, K.W.; Fries, T.-P.: Higher-order XFEM for curved strong and weak dis-continuities,
International Journal for Numerical Methods in Engineering , 82(5):564–590, 2010, doi: 10.1002/nme.2768.[10] Cottrell, J.A.; Reali, A.; Bazilevs, Y.; Hughes, T.J.R.: Isogeometric analysis ofstructural vibrations,
Computer Methods in Applied Mechanics and Engineering ,195(41–43):5257–5296, 2006, doi: 10.1016/j.cma.2005.09.027.[11] Cottrell, J.A.; Hughes, T.J.R.; Reali, A.: Studies of refinement and continuityin isogeometric structural analysis,
Computer Methods in Applied Mechanics andEngineering , 196(41–44):4160–4183, 2007, doi: 10.1016/j.cma.2007.04.007.[12] Cottrell, J.A.; Hughes, T.J.R.; Bazilevs, Y.:
Isogeometric Analysis: Toward Inte-gration of CAD and FEA . John Wiley & Sons, Chichester, England, 2009. ISBN9780470749098.[13] Da Veiga, L.B.; Buffa, A.; Rivas, J.; Sangalli, G.: Some estimates for h – p – k -refinement in isogeometric analysis, Numerische Mathematik , 118(2):271–305, 2011,doi: 10.1007/s00211-010-0338-z.[14] Evans, J.A.; Bazilevs, Y.; Babuška, I.; Hughes, T.J.R.: n -widths, sup–infs, andoptimality ratios for the k -version of the isogeometric finite element method, Com-puter Methods in Applied Mechanics and Engineering , 198(21-26):1726–1741, 2009,doi: 10.1016/j.cma.2009.01.021.[15] Fries, T.-P.; Omerović, S.: Higher-order accurate integration of implicit geometries,
International Journal for Numerical Methods in Engineering , 106(5):323–371, 2016,doi: 10.1002/nme.5121. 2016] Hiemstra, R.R.; Calabrò, F.; Schillinger, D.; Hughes, T.J.R.: Optimal and reducedquadrature rules for tensor product and hierarchically refined splines in isogeometricanalysis,
Computer Methods in Applied Mechanics and Engineering , 316:966–1004,2017, doi: 10.1016/j.cma.2016.10.049.[17] Hiemstra, R.R.; Sangalli, G.; Tani, M.; Calabrò, F.; Hughes, T.J.R.: Fast forma-tion and assembly of finite element matrices with application to isogeometric linearelasticity,
Computer Methods in Applied Mechanics and Engineering , 355:234–260,2019, doi: 10.1016/j.cma.2019.06.020.[18] Höllig, K.:
Finite Element Methods with B-Splines , volume 26 of
Frontiers in Ap-plied Mathematics . SIAM, 2003.[19] Hughes, T.J.R.; Cottrell, J.A.; Bazilevs, Y.: Isogeometric analysis: CAD, finiteelements, NURBS, exact geometry and mesh refinement,
Computer Methods inApplied Mechanics and Engineering , 194(39–41):4135–4195, 2005, doi: 10.1016/j.cma.2004.10.008.[20] Hughes, T.J.R.; Reali, A.; Sangalli, G.: Efficient quadrature for nurbs-based iso-geometric analysis,
Computer Methods in Applied Mechanics and Engineering , 199(5-8):301–313, 2010, doi: 10.1016/j.cma.2008.12.004.[21] Johannessen, K.A.: Optimal quadrature for univariate and tensor product splines,
Computer Methods in Applied Mechanics and Engineering , 316:84–99, 2017, doi:10.1016/j.cma.2016.04.030.[22] Kiendl, J.; Bletzinger, K.-U.; Linhard, J.; Wüchner, R.: Isogeometric shell anal-ysis with Kirchhoff-Love elements,
Computer Methods in Applied Mechanics andEngineering , 198(49–52):3902–3914, 2009, doi: 10.1016/j.cma.2009.08.013.[23] Kudela, L. Highly accurate subcell integration in the context of the finite cellmethod. Master’s thesis, Technical University Munich, 2013.[24] Kudela, L.; Zander, N.; Bog, T.; Kollmannsberger, S.; Rank, E.: Efficient andaccurate numerical quadrature for immersed boundary methods,
Advanced Mod-eling and Simulation in Engineering Sciences , 2(1):1–22, 2015, doi: 10.1186/s40323-015-0031-y.[25] Kudela, L.; Zander, N.; Kollmannsberger, S.; Rank, E.: Smart octrees: Accuratelyintegrating discontinuous functions in 3D,
Computer Methods in Applied Mechanicsand Engineering , 306:406–426, 2016, doi: 10.1016/j.cma.2016.04.006.[26] Lasserre, J.: Integration on a convex polytope,
Proceedings of the American Math-ematical Society , 126(8):2433–2441, 1998, doi: 10.1090/S0002-9939-98-04454-2.[27] Legay, A.; Wang, H.W.; Belytschko, T.: Strong and weak arbitrary discontinuitiesin spectral finite elements,
International Journal for Numerical Methods in Engi-neering , 64(8):991–1008, 2005, doi: 10.1002/nme.1388.2128] Lipton, S.; Evans, J.A.; Bazilevs, Y.; Elguedj, T.; Hughes, T.J.R.: Robustnessof isogeometric structural discretizations under severe mesh distortion,
ComputerMethods in Applied Mechanics and Engineering , 199(5–8):357–373, 2010, doi: 10.1016/j.cma.2009.01.022.[29] Marussig, B.; Hughes, T.J.R.: A review of trimming in isogeometric analysis: Chal-lenges, data exchange and simulation aspects,
Archives of Computational Methodsin Engineering , 25(4):1059–1127, 2018, doi: 10.1007/s11831-017-9220-9.[30] Marussig, B.; Zechner, J.; Beer, G.; Fries, T.-P.: Stable isogeometric analysis oftrimmed geometries,
Computer Methods in Applied Mechanics and Engineering ,316:497–521, 2016, doi: 10.1016/j.cma.2016.07.040.[31] Marussig, B.; Hiemstra, R.; Hughes, T.J.R.: Improved conditioning of isogeometricanalysis matrices for trimmed geometries,
Computer Methods in Applied Mechanicsand Engineering , 334:79–110, 2018, doi: 10.1016/j.cma.2018.01.052.[32] Mousavi, S.; Sukumar, N.: Numerical integration of polynomials and discontinuousfunctions on irregular convex polygons and polyhedrons,
Computational Mechanics ,47(5):535–554, 2011, doi: 10.1007/s00466-010-0562-5.[33] Müller, B.; Kummer, F.; Oberlack, M.: Highly accurate surface and volume inte-gration on implicit domains by means of moment-fitting,
International Journal forNumerical Methods in Engineering , 96(8):512–528, 2013, doi: 10.1002/nme.4569.[34] Nagy, A.P.; Benson, D.J.: On the numerical integration of trimmed isogeometricelements,
Computer Methods in Applied Mechanics and Engineering , 284:165–185,2015, doi: 10.1016/j.cma.2014.08.002.[35] Orszag, S.A.: Spectral methods for problems in complex geometries,
Journal ofComputational Physics , 37(1):70–92, 1980, doi: 10.1016/0021-9991(80)90005-4.[36] Schillinger, D.; Hossain, S.J.; Hughes, T.J.R.: Reduced Bézier element quadraturerules for quadratic and cubic splines in isogeometric analysis,
Computer Methodsin Applied Mechanics and Engineering , 277:1–45, 2014, doi: 10.1016/j.cma.2014.04.008.[37] Schöllhammer, D.; Marussig, B.; Fries, T.-P.: A consistent higher-order isogeometricshell formulation, arXiv preprint arXiv:2012.11975arXiv preprint arXiv:2012.11975