Fast formation and assembly of finite element matrices with application to isogeometric linear elasticity

Recently, a new formation and assembly strategy was proposed in [1], which resulted in significant speedups in the formation and assembly time of the Galerkin mass matrix in isogeometric analysis. The strategy relies on three key ingredients: (1) assembly row by row, instead of element by element; and an efficient formation strategy based on (2) sum factorization and (3) weighted quadrature, that is applied to each specific row of the matrix. Compared to traditional element procedures applied to three dimensional problems, the computational complexity is lowered from O(p(9)) per degree of freedom to O(p(4)), where p is the order of polynomials. This is close to the theoretical minimum of O(p(3)), attained by, for example, collocation. Consequently, this type of formation and assembly scales favorably with polynomial degree, which opens the way for high order isogeometric analysis employing k-refinement, that is, use of maximally smooth, higher order splines. In this work we discuss various important details for the practical implementation of the weighted quadrature formation strategy proposed in [1]. Specifically, we extend the weighted quadrature scheme to accurately integrate the elements of the stiffness matrix in linear elasticity and propose a means of distributing quadrature points for non-uniform, mixed continuity, spline spaces. Furthermore, we discuss efficient access and assignment into the prevalent sparse matrix data structures, namely, Compressed Sparse Row (CSR) and Compressed Sparse Column (CSC). In particular, row-by-row or column-by-column assembly allows matrix rows or columns, respectively, to be formed contiguously in the storage order of the sparse matrix, thereby minimizing the memory overhead and eliminating the addition assignment operation on sparse matrices. Several three-dimensional benchmark problems illustrate the efficiency and efficacy of the proposed formation and assembly technique applied to isogeometric linear elasticity. We show that the accuracy of full Gauss quadrature is maintained while the computational burden of forming the matrix equations is significantly reduced.