PaperlayerView exact revision on arXiv

Source-linked AI summary

Fast formation of isogeometric Galerkin matrices by weighted quadrature

Francesco Calabrò, Giancarlo Sangalli, Mattia Tani

arXiv:1605.01238v2math.NAcs.NA

TL;DR

Isogeometric matrix formation is computationally costly, while existing algorithms grow faster with spline degree than the scale of the matrix’s nonzero entries. The paper uses row-wise weighted quadrature and tensor-product structure to form matrices, reporting substantial speedups and memory-dominated runtimes at sufficiently high degree.

  • Problem

    Current isogeometric assembly algorithms grow faster with spline degree than the O(NDOF p^d) scale of the matrix’s nonzero entries, while matrix formation can dominate linear-solver cost.

    Source: §1 · pp. 2–3

  • Method

    The algorithm loops over matrix rows, applies weighted quadrature, and exploits the tensor-product structure of the basis functions.

    Source: §1 · p. 4, §8 · p. 25

  • Results

    About 62 hours with element-loop SGQ versus 27 seconds with row-loop WQ formed a 203-element-grid mass matrix at p = 10; at sufficiently high p, memory operations dominated measured formation time.

    Source: §1 · pp. 5–6

  • Takeaways & Limitations

    The reported assembly-time savings make row-wise weighted quadrature a faster option than element-loop SGQ in the tested mass-matrix benchmarks.

    Source: §1 · p. 5

    Full mathematical analysis remains future work, and the paper provides numerical optimal-convergence evidence only on a simple one-dimensional benchmark.

    Source: §1 · p. 6, §8 · p. 27

Abstract

from arXiv · show

In this paper we propose an algorithm for the formation of matrices of isogeometric Galerkin methods. The algorithm is based on three ideas. The first is that we perform the external loop over the rows of the matrix. The second is that we calculate the row entries by weighted quadrature. The third is that we exploit the (local) tensor product structure of the basis functions. While all ingredients have a fundamental role for computational efficiency, the major conceptual change of paradigm with respect to the standard implementation is the idea of using weighted quadrature: the test function is incorporated in the integration weight while the trial function, the geometry parametrization and the PDEs coefficients form the integrand function. This approach is very effective in reducing the computational cost, while maintaining the optimal order of approximation of the method. Analysis of the cost is confirmed by numerical testing, where we show that, for $p$ large enough, the time required by the floating point operations is less than the time spent in unavoidable memory operations (the sparse matrix allocation and memory write). The proposed algorithm allows significant time saving when assembling isogeometric Galerkin matrices for all the degrees of the test spline space and paves the way for a use of high-degree $k$-refinement in isogeometric analysis.

1. Introduction

Isogeometric Galerkin matrix formation is costly at high spline degrees, despite the accuracy benefits of smooth splines. The paper proposes row-wise weighted quadrature to reduce assembly cost while retaining the method’s convergence properties.

1. Introduction

  • High-degree, high-regularity splines can deliver higher accuracy per degree of freedom than C0 finite elements, but their computational cost makes quadratic or cubic splines the practical preference in available codes.System-matrix formation can dominate the cost of solving linear PDE problems, while current isogeometric assembly algorithms grow faster than the O(NDOF p^d) nonzero-entry scale.Source: §1 · pp. 2–3
  • Most isogeometric codes use element-wise standard Gaussian quadrature, costing O(NDOF p^(3d)) FLOPs; for the k-method, NEL≈NDOF, while sum-factorization exploits tensor-product structure at O(NDOF p^(2d+1)) FLOPs.Source: §1 · pp. 3–4
  • Prior alternatives interpolate geometry and coefficient factors with precomputed integrals or approximate stiffness matrices by sums of Kronecker products.The reported costs are O(NDOF p^(2d)) for the former and O(NDOF R p^d) FLOPs for the latter.Source: §1 · p. 4
  • The proposed algorithm replaces element-wise assembly with a matrix-row loop, using a row-specific weighted quadrature rule and the tensor-product basis structure.The test function enters the integration weight; the trial function and coefficient factors form the integrand.Source: §1 · pp. 4–5
  • At maximum regularity, weighted quadrature uses two points per direction away from boundaries and p+1 points on boundary knot-spans; sum-factorization gives O(NDOF p^(d+1)) FLOPs.Its exactness conditions are linear in the weights, local to the test function’s support, and compatible with nonuniform knot vectors.Source: §1 · p. 5
  • About 62 hours with element-loop SGQ versus 27 seconds with row-loop WQ formed the mass matrix on a 203-element grid at p = 10; WQ was more than four times faster for quadratics.For sufficiently high p, measured growth was at most O(NDOF p^d), and memory allocation and writes dominated formation time.Source: §1 · pp. 5–6
  • The paper postpones full convergence analysis, offering numerical evidence of optimal-order convergence on a simple one-dimensional benchmark instead.It also leaves parallelization and complete benchmarking beyond the paper’s scope.Source: §1 · p. 6

2. Weighted quadrature

Weighted quadrature (WQ) incorporates the test function into the integration weight, while the trial function and other integrand factors remain in the integrand. Its construction uses a number of quadrature points independent of spline degree, with positive weights providing numerical evidence of stability.

2. Weighted quadrature

  • WQ incorporates the test function into the integration weight, leaving the trial function and other factors in the integrand.For p = 1, the paper starts from composite Cavalieri–Simpson quadrature, exact under the stated setup.Source: §2 · p. 7
  • For degree-p splines, 2p + 1 active points per test function determine weights, while the total number of quadrature points is independent of p.The selected points are knots and knot-span midpoints; the general rule can be computed numerically for arbitrary degree and non-uniform spacing.Source: §2 · pp. 8–9
  • For degrees p = 2 through p = 5, the plotted weights cluster around the basis-function values, providing numerical evidence of positivity and stability.Source: §2 · p. 10, Figure 2 · p. 10
  • For non-constant c, the weighted-quadrature evaluation is generally an approximation, and the symmetry of the integral is not preserved.Source: §2 · p. 11
  • A detailed mathematical analysis of weighted-quadrature error is left for future work.Source: §2 · p. 11

3. Integral arising in isogeometric Galerkin methods

The section formulates isogeometric Galerkin matrix integrals on a spline-parameterized domain using tensor-product B-spline bases. It also describes the factors and sparsity relevant to those integrals.

3. Integral arising in isogeometric Galerkin methods

  • The basis functions are tensor products of univariate B-splines, and the physical basis is obtained by composing the parametric basis with the inverse geometry map.Integrals are computed by change of variable.Source: §3 · p. 12
  • After the change of variable, the matrix integrands involve products of parametric basis functions and the geometry-map Jacobian.Source: §3 · pp. 12–13
  • In more general cases, the factor c includes the equation coefficient and, for NURBS functions, the polynomial denominator.Source: §3 · p. 13
  • The number of nonzero matrix entries is O(NDOF p^d); for d = 1 with maximal regularity, it is (2p + 1)NDOF − p(p + 1).The number depends on degree p and required regularity.Source: §3 · p. 13

4. Computation of the WQ rules

The WQ rules are constructed by incorporating each test basis function into the integration weight. The rules use specified quadrature points and exactness requirements, with weights chosen to satisfy them.

4. Computation of the WQ rules

  • The test B-spline factors are isolated in each univariate integral and used as the weight for a quadrature rule for each index il.Source: §4 · p. 14
  • The rules must satisfy exactness requirements that recover the matrix entry when c is constant; stiffness terms also require derivative terms to be exact.Source: §4 · p. 14
  • The rule for each index is determined by solving a linear system of equations for the unknown weights.Source: §4 · p. 15
  • For maximum regularity, the global grid uses knots and internal-span midpoints, plus p + 1 equally spaced points on boundary spans.The dominant grid size is approximately 2^d NEL.Source: §4 · p. 17
  • Given the points, weights are chosen to meet the exactness conditions; when they are not unique, a minimum-norm condition selects them.The paper states that suitable weights exist for its chosen points by the Schoenberg-Whitney interpolation theorem [18, Theorem 2 ¶XIII].Source: §4 · p. 17

5. Pseudo-codes and computational cost

The operation-count analysis assumes equal element and degree-of-freedom counts across directions, maximum regularity, and many more elements than the degree. Under the chosen quadrature points, each univariate rule uses at most 2p + 1 points.

5. Pseudo-codes and computational cost

  • The analysis assumes the same element and degree-of-freedom counts in every direction and maximum regularity r = p − 1.Source: §5 · p. 17
  • It also assumes nEL ≫ p, which implies nEL ≈ nDOF.Source: §5 · p. 17
  • With the chosen quadrature points, each rule has at most 2p + 1 points.Source: §5 · p. 18

14 end

The cost analysis accounts for basis-function evaluations, exact integrals, coefficient evaluations, and construction of the univariate WQ rules. It uses element-wise Gaussian quadrature for exact integrals of B-spline products.

14 end

  • Because closed forms for exact integrals of B-spline products are available only in particular cases, the analysis uses element-wise Gaussian quadrature.The exact-integral literature is noted in.Source: §14 · pp. 18–19
  • Evaluating the basis functions at the quadrature points costs O(p^4 nDOF) FLOPs.Source: §14 · p. 19
  • Evaluating the multivariate coefficient functions at all quadrature points has cost dependent on how those coefficients are evaluated.If coefficients use O(p^d) linear combinations of B-spline values or derivatives, the stated cost is C(d)p^d per point.Source: §14 · p. 19
  • The univariate WQ rule construction solves possibly underdetermined linear systems of dimension approximately (2p)^2, with total cost proportional to nDOF.The text describes this cost as negligible.Source: §14 · p. 19

6. Formation of the mass matrix

The mass-matrix algorithm computes entries row by row using weighted quadrature, then expresses the nested sums through matrix-tensor products and sum factorization. This structure gives a cost of O(NDOF p^(d+1)) FLOPs.

6. Formation of the mass matrix

  • For each matrix row, the algorithm computes the mass-matrix entries using the quadrature rules; analogous formulas and algorithms apply to the stiffness matrix.The row-wise computation is summarized in Algorithm 3.Source: §6 · p. 19
  • The nested summations are rewritten as matrix-tensor products, enabling a sum-factorization computation of each row.The m-mode product replaces tensor summations with operations that can be implemented as matrix-matrix products.Source: §6 · pp. 19–21
  • With quadrature-point counts of O(p), computing one row costs O(p^(d+1)) FLOPs, for a total cost of O(NDOF p^(d+1)).The row computation includes all nonzero entries in that row.Source: §6 · p. 21
  • In MATLAB, each matrix-tensor product is computed as a matrix-matrix product, a BLAS level 3 operation that typically runs efficiently on modern computers.Tensorizing the operations helps avoid interpreted-language loops.Source: §6 · p. 21

7. Numerical tests

The numerical tests examine solution convergence in 1D and mass-matrix formation in 3D. In 1D, weighted quadrature retains optimal convergence rates, with a small accuracy difference from standard Gaussian quadrature for some even degrees.

7. Numerical tests

  • The tests cover 1D solution convergence and 3D mass-matrix formation, with the 3D benchmarks run sequentially on one core.The authors note that the row-loop algorithm may be better suited to parallel implementation, but benchmark sequential execution.Source: §7 · pp. 21–22

7.1. Convergence of approximate solution in 1D

  • The proposed procedure preserves the Galerkin method's optimal convergence rates in all tested degrees for the 1D problem.The convergence history reports L∞, L2, and H1 errors; the figure compares WQ with element-wise SGQ.Source: §7.1 · pp. 22–23, Figure 4 · p. 23
  • For even degrees greater than 2, SGQ is slightly more accurate than WQ in the L2 and L∞ norms, while WQ's loss of symmetry did not deteriorate convergence order.Source: §7.1 · p. 23

in energy and lower-order norms in our numerical benchmarking, see Figure

The mass-matrix timing tests compare row-loop weighted quadrature with GeoPDEs' element-loop standard Gaussian quadrature. WQ's advantage grows with degree, and at higher degrees sparse-matrix allocation and writing dominate its runtime.

7.2. Time for the formation of matrices

  • At p=10, GeoPDEs takes more than 62 hours versus 27 seconds for row-loop WQ mass-matrix formation on a 203-element grid.The comparison is with the optimized, SGQ-based GeoPDEs MATLAB library.Source: §7.2 · p. 24, Figure 5 · p. 24, in energy and lower-order norms in our… · p. 24
  • For p>5, MATLAB's sparse call dominates total runtime, while matrix-tensor product time grows more mildly than the theoretical O(NDOF p^4) FLOP estimate.The sparse-call time is proportional to p^3, matching the O(NDOF p^3) nonzero count.Source: §7.2 · p. 25
  • The measured computation time depends linearly on NDOF.Source: §7.2 · p. 25

8. Conclusions

The proposed method combines row-wise assembly, weighted quadrature, and tensor-product sum factorization. Its tests report substantial gains over standard approaches, particularly at high degree, while identifying non-symmetry and several open development directions.

8. Conclusions

  • The method combines a row loop, weighted quadrature to reduce quadrature points, and sum factorization that exploits the tensor-product B-spline basis.The sum-factorization implementation follows.Source: §8 · p. 25
  • At p=6, WQ forms a mass matrix in seconds rather than the hours required by GeoPDEs, with runtime dominated by MATLAB's unavoidable sparse call.The comparison uses GeoPDEs 3.0, which has an optimized standard design.Source: §8 · p. 26
  • WQ does not preserve symmetry for symmetric differential operators, although the authors expect it to retain the other main properties of Galerkin formulations.Source: §8 · p. 27
  • Future work includes full mathematical analysis, implementation within GeoPDEs, and extension to non-tensor-product spaces such as T-splines and hierarchical splines.The authors expect some advantages to be maintained for those spaces.Source: §8 · p. 27
Loading 1605.01238v2…