Efficient matrix assembly in isogeometric analysis with hierarchical B-splines (Q2226302)

The primary motivation for this work is to enhance the applicability of isogeometric analysis (IGA), as a generalization with modern approach to the classical finite element analysis (FEA). The present paper describes an efficient method referred to as adaptive isogeometric method (AIGM), to derive the isogeometric Galerkin discretizations of an elliptic boundary value problem, with hierarchical B-splines that allow local refinement. Hierarchical B-splines that provide local refinement capabilities, have already been shown to be a promising tool for developing adaptive isogeometric methods. This is in line with the theory of adaptive methods developed for finite elements, recently well reviewed in the literature. Based on their previous work and pioneering work done in IGA, the authors propose an efficient matrix assembly approach for bivariate hierarchical B-splines as building blocks for approximating geometry of the physical domain and unknown solution fields, primarily to address the issue of high computational cost for assembling the system matrices in IGA. To the best of author's knowledge, the research output, at the time of preparation of the manuscript on this topic was yet to achieve the desired goal. In the 1st part of the paper, the basic concepts, definitions of the global index space of hierarchical B-splines (HB-splines), and some important properties of the admissible HB-splines are presented. This is in line with providing the theoretical foundation of the method. In the second part, to describe the proposed method for assembling the system matrices arising in isogeometric Galerkin discretizations of partial differential equations (PDEs), authors consider a general second-order linear elliptic boundary value problem and its weak (variational) formulation. In the framework of isogeometric analysis, the physical domain $\Omega$ is parameterized using HB-spline functions. Applying the Galerkin approach to the variational problem, isogeometric Galerkin discretizations of an elliptic boundary value problem is derived. Solving the discretized variational problem, requires the formation of several matrices and vectors. Matrix assembly, a crucial step, especially for high polynomial degrees involves three stages: spline projection, building look-up tables and assembling the matrices via sum-factorization. The reminder of the paper, describes the three stages in more detail. This is accomplished by the concept of global fitting or quasi-interpolation methods. Quasi-interpolation (QI) is another well-established methodology to construct approximations of a given function by splines. Typically it leads to lower computational costs than global fitting. The properties of various quasi-interpolants for bivariate hierarchical splines are summarized in a tabular form. The next stage of the approach is to build look-up tables. By exploiting the symmetry and compact supports of the B-splines, the look-up tables are compressed considerably. With the pre-computed compact look-up tables at hand, the sum-factorization technique to accelerate the matrix assembly procedure is employed. To realize the fast and efficient assembly, an algorithm for Stiffness Matrix Assembly is provided. The main contributions of the paper are: \begin{itemize} \item The authors combine the techniques of quasi-interpolation, building look-up tables and sum factorization, thereby introducing a fast matrix assembly algorithm. \item The complexity of the method has the order \(O(Np^3)\) under a mild assumption about mesh admissibility, where \(N\) and \(p\) denote the number of degrees of freedom and spline degree respectively. \item Finally, several experimental results are obtained to verify and validate the theoretical results and to establish the efficacy of the proposed method. \end{itemize} A detailed complexity analysis demonstrates that the computational cost for matrix assembly has order \(O(Np^3)\) under a mesh admissibility assumption; a vast improvement from \(O(Np^6)\), based on element-wise Gaussian quadrature. Several numerical tests in terms of the \(N\)-dependence and \(p\)-dependence of the computation time are provided graphically, which confirm the theoretical result of the new algorithm. These tests are conducted on four different computational domains. Each domain is parameterized by single-patch B-splines, and the corresponding iso-parametric curves are demonstrated diagrammatically. In addition, the convergence behavior and accuracy of the proposed method are also verified by applying it to solving a Poisson problem with an adaptive refinement strategy. This qualifies as a benchmark example for code correctness and verification with known analytical solution is used. The authors have also compared the adaptive refinement by the new method with the one by Gauss quadrature, which indicates that the accuracy of the method is maintained provided the spline projection is sufficiently accurate. Finally, authors' conclusions concerning new method are drawn based on the basis of numerical tests and benchmark example, Also scope for future research is discussed indicating extension of the present work to matrix-free applications. Observations and comments \begin{itemize} \item[1)] IGA has been shown to be a powerful and efficient high-order discretization method for the numerical solution of PDEs. However, recent study shows that, it suffers from a few drawbacks particularly with regard to computational cost. The algorithm presented is a step in the right direction. \item[2)] N and p dependence of the computation time play a crucial role in assessing the performance of the method. Numerical tests conducted illustrate the substantial gains in computation time without compromising the accuracy for the adaptive strategy. \item[3)] The fast assembling of stiffness and mass matrices is a key issue in Isogeometric Analysis, particularly if the spline degree is increased. The method proposed addresses the issue successfully. \item[4)] Motivation for this new approach was presented together with theoretical foundations of the method. The numerical tests conducted for four different computational domains in size and shape amply demonstrate the superiority of IGA over FEA. \item[5)] Throughout, the isoparametric philosophy and guiding principle is invoked, that is, the solution space for dependent variables is represented in terms of the same functions which represent the geometry. \item[6)] Isogeometric Analysis proves to be effective in cases where exact geometry representation plays crucial role. The geometrical error estimates are prerequisite for the proposed approach. \item[7)] In this work, several challenges and pitfalls are described and an effort is made to circumvent those supporting theoretical results with experimental findings. \item[8)] Theoretical background and practical applications must go hand in hand. Numerical experiments conducted underpin the theoretical observations. \item[9)] Furthermore, to fully utilize the advantages of the method, new systems have to be developed and tested. Fortunately, authors have been profited from the already existing efficient IGA algorithms from the literature survey. Every method inherently has certain limitations. In the opinion of the authors, the present implementation hasn't been fully optimized. \item[10)] Since the developments of IGA to eliminate the defects of the conventional FEM, several researchers have devoted to developing new TO methods and discussing their applications using IGA, rather than FEM. The present paper exemplifies this trend justifiably. \end{itemize}