Computational Costs of Multi-Frontal Direct Solvers with Analysis-Suitable T-Splines

: In this paper, we consider the computational cost of a multi-frontal direct solver used for the factorization of matrices resulting from a discretization of isogeometric analysis with T-splines, and analysis-suitable T-splines. We start from model projection or model heat transfer problems discretized over two-dimensional meshes, either uniformly reﬁned or reﬁned towards a point or an edge. These grids preserve several symmetries and they are the building blocks of more complicated grids constructed during adaptive isotropic reﬁnement procedures. A large class of computational problems construct meshes reﬁned towards point or edge singularities. We propose an ordering that permutes the matrix in a way that the computational cost of a multi-frontal solver executed on adaptive grids is linear. We show that analysis-suitable T-splines with our ordering, besides having other well-known advantages, also signiﬁcantly reduce the computational cost of factorization with the multi-frontal direct solver. Namely, the factorization with N T-splines of order p over meshes reﬁned to a point has a linear O( Np 4 ) cost, and the factorization with T-splines on meshes reﬁned to an edge has O( N2 p p 2 ) cost. We compare the execution time of the multi-frontal solver with our ordering to the Approximate Minimum Degree (AMD) and Cuthill–McKee orderings available in Octave.


Introduction
Higher-order and continuity basis functions, such as B-splines, Non-Uniform Rational B-splines (NURBS) [1], and T-splines [2], are used in computer-aided design (CAD)/computer-aided engineering (CAE) systems for modeling of the geometry of engineering design.The idea of the isogeometric analysis (IGA) [1] is to use such higher-order and continuity basis functions for solving different engineering problems with the finite element method.Besides B-splines, NURBS, and T-splines, some other families of functions were introduced recently in the context of IGA [3][4][5][6][7][8][9][10].
Multi-frontal solvers are popular tools for the factorization of sparse matrices [11][12][13].The advantage of multi-frontal solvers is that they provide an accurate numerical solution for any computational problem, symmetric, unsymmetric, positive definite, non-positive definite, double precision, or complex.The price to pay is the computational cost, usually higher than that for iterative solvers.This computational cost in the context of isogeometric analysis with higher-order B-spline basis functions has been analyzed in [14], and in the context of B-splines with C 0 separators, which are equivalent to Lagrange polynomials, in [15].The B-spline basis functions are defined over regular patches of elements, and the computational cost of factorization for a mesh with N B-spline basis functions of order p can be estimated as O(N 2 ) for B-splines with C 0 continuity [15] and as O(N 2 p 3 ) for B-splines with C p−1 continuity for three-dimensional grids [14].It has also been shown in [16] that local reduction in B-spline continuity on uniform grids allows reducing the constant in front of the computational cost.
However, these results refer to the uniform patch of elements, and the computational cost of the multi-frontal solver of IGA on adaptive grids has not been analyzed so far.For extension of IGA into a non-uniform adaptive grid, we need to move from B-spline basis functions to one of their extensions.We selected the T-spline basis functions [2] since they are most popular for adaptive computations with higher-order and continuity basis functions.Still, there are several other options for discretization available [3][4][5][6][7][8][9][10], and we may analyze the computational costs of these alternative basis functions in possible future work.
The standard finite element method uses the basis functions of C 0 continuity, and the computational costs of the multi-frontal solver for adaptive grids with standard FEM have been analyzed in [15].In this paper, for the first time, we will analyze the computational cost of multi-frontal solvers applied to isogeometric analysis with T-spline basis functions on adaptive grids.
Multi-frontal solvers construct orderings and elimination trees based on a matrix's sparsity pattern resulting from discretization over the computational mesh.They utilize different algorithms for the construction of the orderings.There are several possible ordering algorithms focused on analyzing of the sparsity pattern of the matrix (e.g., Approximate Minimum Degree (AMD) [17], Paderborn ORDering tools (PORD) [18], and nested dissections [19]).It has been shown that the nested dissection algorithm is optimal for computations with uniform grids [20].There are also some ordering algorithms analyzing the structure of the computational mesh [21][22][23][24] instead of the structure of the sparse matrix.
This paper focuses on the T-spline basis functions, which can be defined on adaptive grids.We start from recalling the definition of B-spline basis functions, where B i,0 (x) = 1 for x in [x i ,x i+1 ], and B i,p (x) = (x − x i )/(x i+p − x i )B i,p−1 (x) + (x i+p+1 − x)/(x i+p+1 − x i+1 )B i+1,p−1 (x) for x in [x i ,x i+p+1 ].If the knot points are repeated in the denominator, this term is cancelled (since division by zero is not allowed).The B-splines of order p span over p + 1 elements and they are defined with p + 2 knot points, which are introduced into the recursive formula for B-spline basis functions.Thus, the B-spline basis in 1D is defined with a knot vector, e.g., [0 0 1 2 2] defines linear C 0 B-splines, [0 0 0 1 2 2 2] defines quadratic C 1 B-splines, and [0 0 0 1 1 2 2 2] defines quadratic C 0 B-splines.
Let us follow the state-of-the-art process of constructing the sparse matrix, ordering, and the elimination tree for an exemplary two-dimensional mesh with C 0 B-spline basis functions.The B-spline basis functions are usually defined with a tensor product of one dimensional B-splines B ij;p (x,y) = B i;p (x)B j,p (y), described by knot vectors [1].The basis functions obtained by a tensor product of the 1D B-splines defined by knot vectors [0 0 0 1 1 2 2 2] and [0 0 0 1 1 2 2 2] are presented in Figure 1.
Symmetry 2020, 12, x FOR PEER REVIEW 2 of 20 a non-uniform adaptive grid, we need to move from B-spline basis functions to one of their extensions.We selected the T-splines basis functions [2], since they are most popular for adaptive computations with higher-order and continuity basis functions.Still, there are several other options for discretization available [3,4,5,6,7,8,9,10], and we may analyze the computational costs of these alternative basis functions in possible future work.
The standard finite element method uses the basis functions of C 0 continuity, and the computational costs of the multi-frontal solver for adaptive grids with standard FEM have been analyzed in [15].In this paper, for the first time, we will analyze the computational cost of multi-frontal solvers applied to isogeometric analysis with T-spline basis functions on adaptive grids.
The multi-frontal solvers construct the orderings and the elimination trees based on the matrix's sparsity pattern resulting from discretization over the computational mesh.They utilize different algorithms for the construction of the orderings.There are several possible ordering algorithms focused on analyzing of the sparsity pattern of the matrix (e.g., AMD [17], PORD [18], nested-dissections [19]).It has been shown that the nested-dissection algorithm is optimal for computations with uniform grids [20].There are also some ordering algorithms analyzing the structure of the computational mesh [21,22,23,24], instead of the structure of the sparse matrix.
This paper focuses on the T-spline basis functions, which can be defined on adaptive grids.
Let us follow the state-of-the-art process of constructing the sparse matrix, ordering, and the elimination tree for an exemplary two-dimensional mesh with C 0 B-spline basis functions.The B-spline basis functions are usually defined with a tensor product of one dimensional B-splines Bij;p(x,y)=Bi;p(x)Bj,p(y), described by knot-vectors [1].The basis functions obtained by a tensor product of the 1D B-splines defined by knot vectors [0 0 0 1 1 2 2 2] and [0 0 0 1 1 2 2 2] are presented in Figure 1.We enumerate the B-spline basis functions from 1 to 25, since we have five one-dimensional B-splines in each direction, and the tensor product results in 5•5 = 25 basis functions.We construct a matrix (see Figure 2) where rows and columns correspond to the B-splines, and the non-zero values represent overlaps of spans of the B-splines.The particular values depend on the problem being discretized with B-splines over the computational mesh-e.g., the projection problem involves integrals with the multiplications of the pairs of B-splines B i;p (x)B j,p (y)B k;p (x)B l,p (y)dxdy, and the heat transfer problem involves integrals with the multiplications of the gradients of the pairs B-splines ∇(B i;p (x)B j,p (y)) • ∇(B k;p (x)B l,p (y))dxdy.We enumerate the B-spline basis functions from 1 to 25, since we have five one-dimensional B-splines in each direction, and the tensor product results in 55 = 25 basis functions.We construct a matrix (see Figure 2) where rows and columns correspond to the B-splines, and the non-zero values represent overlaps of spans of the B-splines.The particular values depend on the problem being discretized with B-splines over the computational mesh-e.g., the projection problem involves integrals with the multiplications of the pairs of B-splines Bi;p(x)Bj,p(y)Bk;p(x)Bl,p(y)dxdy, and the heat transfer problem involves integrals with the multiplications of the gradients of the pairs B-splines (Bi;p(x)Bj,p(y))(Bk;p(x)Bl,p(y))dxdy.The nested dissections ordering algorithm [19] constructs the elimination graph (see Figure 2) representing the matrix's sparsity pattern.The nodes in the graph represent the B-spline basis functions.They also represent the rows in the matrix.The edges represent the overlaps of the B-splines and the non-zero values in the columns of the matrix.
To describe this state-of-the-art nested dissection algorithm, we need to recall the Breadth-First Search (BFS) algorithm and introduce the concept of a separator and a peripheral node.The nested dissections ordering algorithm [19] constructs the elimination graph (see Figure 2) representing the matrix's sparsity pattern.The nodes in the graph represent the B-spline basis functions.They also represent the rows in the matrix.The edges represent the overlaps of the B-splines and the non-zero values in the columns of the matrix.
To describe this state-of-the-art nested dissection algorithm, we need to recall the Breadth-First Search (BFS) algorithm and introduce the concept of a separator and a peripheral node.
The Breadth-First Search algorithm (BFS) is used for visiting nodes of the graph; it starts at some arbitrary node of a graph and explores the neighboring nodes first before moving to the next layer of neighbors.The execution of the BFS algorithm can identify levels with a group of nodes.In the example presented in Figure 2, the BFS algorithm was executed from node 23, and it has found three levels: the first level with node 23, the second level with all the neighbors of node 23, and the third level with all neighbors of the nodes from the second level, which were not visited before.
A separator it is a set of nodes from the middle level having neighbors from the next level.In this example, the separator consists of nodes 11,12,13,14, and 15.
A peripheral node is a node with a higher number of levels obtained from the BFS algorithm's execution, starting from this node, and the shortest separator resulting from execution of the BFS algorithm (starting from this node).In the example presented in Figure 2, it is node 23 (but some other boundary nodes may also qualify, e.g., 3,11, or 15).
The nested dissections algorithm recursively finds separators.It first finds a peripheral node and a separator; then, it partitions the graph into two subgraphs according to the separator and it finds new separators for subgraphs in a recursive way.This is illustrated in Figure 3.
On the base of graph partition and separators, the elimination tree can be identified by listing the separators in a tree structure (see Figure 3).The elimination of unknowns is then performed in a multi-frontal manner.The rows are eliminated in the order resulting from browsing of the elimination tree from leaves up to the root.In our example, the order of elimination is:  (19,20,24, and 25).
A peripheral node is a node with a higher number of levels obtained from the BFS algorithm's execution, starting from this node, and the shortest separator resulting from execution of the BFS algorithm (starting from this node).In the example presented in Figure 2, it is node 23 (but some other boundary nodes may also qualify, e.g., 3,11, or 15) The nested dissections algorithm recursively finds separators.It first finds a peripheral node and a separator; then, it partitions the graph into two subgraphs according to the separator and it finds new separators for subgraphs in a recursive way.This is illustrated in Figure 3.
On the base of graph partition and separators, the elimination tree can be identified by listing the separators in a tree structure (see Figure 3).The elimination of unknowns is then performed in a multi-frontal manner.The rows are eliminated in the order resulting from browsing of the elimination tree from leaves up to the root.In our example, the order of elimination is:  In the next level, we merge the frontal matrices according to the structure of the elimination tree and we perform further eliminations:

•
Merge the first and the second frontal matrices into a fifth frontal matrix, and eliminate rows 3 and 8.This time, the frontal matrix has a size of seven (it contains rows 3,8,11,12,13,14, and 15) and this step eliminates two rows (3 and 8).

•
Merge the third and fourth frontal matrices into a sixth frontal matrix and eliminate rows 18 and 23.This time, the frontal matrix has a size of seven (it contains rows 18,23,11,12,13,14, and 15) and this step eliminates two rows (18 and 23).

•
Merge the fifth and the sixth frontal matrices and eliminate all the rows.This frontal matrix has a size of five (it contains rows 11,12,13,14, and 15) and this step eliminates all the rows.
This procedure is implemented in multi-frontal solvers, e.g., those available through the PETSc interface [25], such as MUMPS solver [11][12][13], Scalapack [26], PaStiX [27], and SuperLU [28].To find the computational cost of the elimination process, we can look at the spans of B-spline basis functions and the way they overlap with other B-spline functions.We can identify the dimensions of the frontal matrices (called a) and the number of eliminated rows (B-splines basis functions, called b).The cost of elimination of b rows from the matrix of size a is O(ab 2 ) [29].Thus, the procedure's total cost is Σi = 1, . . ., s a i b i 2 , where s is the number of steps (number of layers in the elimination tree representing the separators), and we need to identify the number of steps and the dimensions of the frontal matrices and the eliminated rows.T-spline basis functions are a generalization of B-splines into non-uniform grids.From the point of view of the multi-frontal solver's computational cost, the exact formula of the basis functions is not important, only the support of the function.The knot vectors for T-splines span over the adaptive grids.For the simplicity of presentation, we restricted our example to even T-splines.The centers of such even T-splines are located over the center of an element.We have one knot vector span in the horizontal direction, where the knot points are located at the cross-sections with mesh edges.We also have another knot vector span in the vertical direction, defined in a similar way.We have one T-spline associated with one mesh element, where its center is located.Thus, the T-splines can be identified with mesh elements.In Figure 4, we present the T-spline span over the mesh refined towards the corner.


Merge the fifth and the sixth frontal matrices and eliminate all the rows.This frontal matrix has a size of five (it contains rows 11,12,13,14, and 15) and this step eliminates all the rows.This procedure is implemented in multi-frontal solvers, e.g., those available through the PETSc interface [25], such as MUMPS solver [11][12][13], Scalapack [26], PaStiX [27], and SuperLU [28].To find the computational cost of the elimination process, we can look at the spans of B-spline basis functions and the way they overlap with other B-spline functions.We can identify the dimensions of the frontal matrices (called a) and the number of eliminated rows (B-splines basis functions, called b).The cost of elimination of b rows from the matrix of size a is O(ab 2 ) [29].Thus, the procedure's total cost is Σi = 1, …, s aibi 2 , where s is the number of steps (number of layers in the elimination tree representing the separators), and we need to identify the number of steps and the dimensions of the frontal matrices and the eliminated rows.
T-spline basis functions are a generalization of B-splines into non-uniform grids.From the point of view of the multi-frontal solver's computational cost, the exact formula of the basis functions is not important, only the support of the function.The knot vectors for T-splines span over the adaptive grids.For the simplicity of presentation, we restricted our example to even T-splines.The centers of such even T-splines are located over the center of an element.We have one knot vector span in the horizontal direction, where the knot points are located at the cross-sections with mesh edges.We also have another knot vector span in the vertical direction, defined in a similar way.We have one T-spline associated with one mesh element, where its center is located.Thus, the T-splines can be identified with mesh elements.In Figure 4, we present the T-spline span over the mesh refined towards the corner.T-splines are usually defined over so-called analysis-suitable grids.They were proposed to deal with the approximation properties of T-splines [30], but they were not analyzed in the context of the solvers' computational cost, which is the topic of this paper.Analysis-suitable grids are obtained by adding so-called T-junction extensions [29].The edges that end up with hanging nodes (nodes that belong to three elements, located on the edge between two small elements and one big element) are extended into p/2 blocks, where p denotes the order of T-splines' span over the analysis-suitable T-splines are usually defined over so-called analysis-suitable grids.They were proposed to deal with the approximation properties of T-splines [30], but they were not analyzed in the context of the solvers' computational cost, which is the topic of this paper.Analysis-suitable grids are obtained by adding so-called T-junction extensions [29].The edges that end up with hanging nodes (nodes that belong to three elements, located on the edge between two small elements and one big element) are extended into p/2 blocks, where p denotes the order of T-splines' span over the analysis-suitable mesh.The analysis-suitable grid obtained by the extension of the grid refined towards a point and the resulting T-splines are presented in Figure 4.The T-splines spanning over the analysis-suitable mesh are called the analysis-suitable T-splines.
The structure of this paper is as follows.In Section 2, we present the motivation for this research.In Section 3, we analyze the computational cost of T-splines on uniform grids.In this case, the T-splines are equivalent to B-splines (in the sense of supports), and the analysis-suitable grid is the same as the original uniform grid.Next, in Section 4, we estimate the computational cost of T-splines on an adaptive grid, in Section 4.1 on grids refined towards a point, and in Section 4.2 on grids refined towards an edge.In the following Section 5, we estimate the computational cost of the analysis-suitable T-splines defined on grids refined towards a point, in Section 5.1.and defined on grids refined towards an edge. in Section 5.2.Section 6 is devoted to numerical verification of the proposed ordering algorithms.We conclude the paper in Section 7.

Motivation
The computational cost of a multi-frontal solver varies from N to N 3 , depending on the computational mesh's sparsity pattern or the computational mesh structure and basis functions used for discretization.The general formula for any mesh and basis functions is not known.It has been estimated for uniform and some adaptive grids and selected basis functions.
The computational costs for the multi-frontal solver executed on uniform meshes are known for discretization with the linear finite element method (FEM) [20].The cost has also been estimated for the finite element method with hierarchical Lagrange basis functions of order p [31] and the isogeometric finite element method with B-spline basis functions [14].Namely, for the classical FEM on 2D uniform grids, it is O(N 1.5 ), and for 3D uniform grids, it is O(N 2 ).For IGA discretization, we have an extra term related to the polynomial order, and for B-spline-based discretization with uniform grids, it is O(N 2 p 3 ) for 3D grids.This cost has been not estimated for 2D grids yet.For uniform grids, the T-splines and B-splines and analysis-suitable T-splines are equivalent by definition.The computational costs for the multi-frontal solver executed on adaptive grids are known for hierarchical Lagrange polynomials for different singularities [15,32].They are not estimated for T-spline and analysis-suitable T-spline basis functions.In [32], there is an asymptotic estimation of the computational cost for refined grids based on Lagrange polynomials.It is estimated as O(N max(3(q−1)/q,1) ), where q is the type of singularity, and this result does not depend on the dimension of the mesh.However, this computational cost is not known for T-spline and analysis-suitable T-spline basis functions for any kind of grids.In general, for discretizations with T-splines, we expect these costs to be higher than those for standard FEM and IGA, and for analysis-suitable T-splines, we expect them to have an extra factor related to the polynomial orders.In this paper, we derive these costs and factors.The optimization of computational cost of IGA solvers is important, since IGA has multiple applications to phase field modeling with either Cahn-Hilliard [33] or Navier-Stokes-Korteweg higher-order models [34].The IGA-FEM has been also successfully applied for solution of non-linear problems, such as wind turbine aerodynamics [35], incompressible hyper-elasticity [36], or blood flow simulations [37].This work's main motivation is to derive computational costs for a multi-frontal solver executed on different grids, either uniform or refined towards a point or edge singularity, with T-spline basis functions and analysis-suitable T-spline basis functions.Based on these estimations, we propose ordering algorithms resulting in a permutation of the sparse matrix such that the computational cost of multi-frontal solver algorithm factorization is linear for grids refined towards point or edge singularities when using analysis-suitable T-splines.
Input: Mesh, p (order of analysis suitable T-splines) Level = 1; Nodes = nodes from the top level of the mesh; Tree = build tree (Level, Nodes); Output: ordering = post-order traversal of Tree; build tree (Level, Nodes) { Node = create node (Nodes); Level = Level + 1; Node -> left = build tree (Level, first p + 1 nodes from Level neighboring to Nodes); Node -> right = build tree (Level, last p + 1 nodes from Level neighboring to Nodes); return Node; } The algorithm constructs a binary tree with nodes assigned to elements on levels of the refinement.The root is the largest element, and the leaves are the smallest elements.It orders the elements and T-splines assigned to them by post-order traversal of the tree.For grids refined towards a point, it is equivalent to browsing the mesh layer by layer; see Figure 5.
Computational grids used in CAD/CAE systems can be partitioned into sub-grids, either uniform (without local refinements) or refined towards a point or edge singularity, or they are a composition of such refinements towards several point or edge singularities.To estimate the computational costs for all kinds of grids, we take three representative grids: the point singularity, the edge singularity, and the uniform grid.Most other grids in two dimensions are constructed by combining these three representative grids, and so are the computational costs.Our ordering algorithm can be easily combined with the domain decomposition approach of complex grids.
Node -> right = build tree (Level, last p + 1 nodes from Level neighboring to Nodes); return Node; } The algorithm constructs a binary tree with nodes assigned to elements on levels of the refinement.The root is the largest element, and the leaves are the smallest elements.It orders the elements and T-splines assigned to them by post-order traversal of the tree.For grids refined towards a point, it is equivalent to browsing the mesh layer by layer; see Figure 5.
Computational grids used in CAD/CAE systems can be partitioned into sub-grids, either uniform (without local refinements) or refined towards a point or edge singularity, or they are a composition of such refinements towards several point or edge singularities.To estimate the computational costs for all kinds of grids, we take three representative grids: the point singularity, the edge singularity, and the uniform grid.Most other grids in two dimensions are constructed by combining these three representative grids, and so are the computational costs.Our ordering algorithm can be easily combined with the domain decomposition approach of complex grids.

Uniform Mesh
In this section, we estimate the computational cost of the multi-frontal solver on uniform grids.In this case, the B-splines are equivalent to T-splines and to the analysis-suitable T-splines.All these basis functions have identical supports on uniform grids.We assume that we have a uniform mesh that consists of 2 2s of blocks, each block being of size (p+ 1)(p+ 1).The number of mesh elements N is equal to (p+ 1)(p+ 1)2 2s and so is the number of T-splines.Some examples of meshes with p = 2 and 4 are presented in Figure 6.We assume even polynomial orders for the simplicity of the derivations, but our results remain valid also for odd cases.In the first step (I = 0), the multi-frontal algorithm will eliminate T-splines connected with each block's middle element.It will generate a frontal matrix with rows and columns related to the T-splines from this particular block.Only one row representing the T-splines from the center of the block can be eliminated here, since it is the support span inside the block.We will talk about the T-splines to be eliminated for simplicity of presentation, the ones having supports inside the processed blocks, and not about the rows and matrices.
these basis functions have identical supports on uniform grids.We assume that we have a uniform mesh that consists of 2 2s of blocks, each block being of size (p + 1)(p + 1).The number of mesh elements N is equal to (p + 1)(p + 1)2 2s and so is the number of T-splines.Some examples of meshes with p = 2 and 4 are presented in Figure 6.We assume even polynomial orders for the simplicity of the derivations, but our results remain valid also for odd cases.In the first step (I = 0), the multi-frontal algorithm will eliminate T-splines connected with each block's middle element.It will generate a frontal matrix with rows and columns related to the T-splines from this particular block.Only one row representing the T-splines from the center of the block can be eliminated here, since it is the support span inside the block.We will talk about the T-splines to be eliminated for simplicity of presentation, the ones having supports inside the processed blocks, and not about the rows and matrices.The algorithm will then join four neighboring blocks (33 for p = 2) into one bigger block and eliminate these T-splines that can be eliminated (T-splines located on the common edge, with support spans inside the two merged blocks).The algorithm will follow this pattern in the following steps: join four neighboring blocks into bigger blocks and eliminate T-splines that can be eliminated.Figures 7  and 8 present, step by step (iteration step i equal to 0, 1, and 2), the elements for which the T-splines can be eliminated (green color), for which the T-splines were eliminated in previous steps (white color), and the rest of the elements (blue color) for the mesh with p = 2 and s = 3. Figures 9 and 10 present, step by step (iteration step i equal to 0, 1, and 2 and 3), the elements for which the T-splines can be eliminated (green color), the elements for which the T-splines were eliminated in previous steps (white color), and the rest of elements (blue color).The pictures concern the mesh with p = 4 and the mesh with s = 3.For simplicity, only one block in each step of the algorithm is presented.The number of elements that can be eliminated in one block in each step i, for i> = 1, is proportional to the strip's thickness multiplied by the length of the strip (elements denoted by green color).The thickness of the strip is constant for each p and is equal to p.The thickness of the strip is equal to 2 for p = 2 and 4 for p = 4 (see Figures 7-10), and it will be 6 for p = 6.The length of the strip depends on the value of p and on the iteration step i (the bigger value of i, the longer the strip).The algorithm will then join four neighboring blocks (3•3 for p = 2) into one bigger block and eliminate these T-splines that can be eliminated (T-splines located on the common edge, with support spans inside the two merged blocks).The algorithm will follow this pattern in the following steps: join four neighboring blocks into bigger blocks and eliminate T-splines that can be eliminated.Figures 7  and 8 present, step by step (iteration step i equal to 0, 1, and 2), the elements for which the T-splines can be eliminated (green color), for which the T-splines were eliminated in previous steps (white color), and the rest of the elements (blue color) for the mesh with p = 2 and s = 3. Figures 9 and 10 present, step by step (iteration step i equal to 0, 1, and 2 and 3), the elements for which the T-splines can be eliminated (green color), the elements for which the T-splines were eliminated in previous steps (white color), and the rest of elements (blue color).The pictures concern the mesh with p = 4 and the mesh with s = 3.For simplicity, only one block in each step of the algorithm is presented.The number of elements that can be eliminated in one block in each step i, for i ≥ 1, is proportional to the strip's thickness multiplied by the length of the strip (elements denoted by green color).The thickness of the strip is constant for each p and is equal to p.The thickness of the strip is equal to 2 for p = 2 and 4 for p = 4 (see Figures 7-10), and it will be 6 for p = 6.The length of the strip depends on the value of p and on the iteration step i (the bigger value of i, the longer the strip).The solution of this recursive equation is the following: As such, the number of elements which can be eliminated in i-th step for a given p is proportional to p2 i p = 2 i p 2 (thickness multiplied by length of the strip) (elements denoted by green color).In one block, the number of all elements which were not eliminated before is proportional to the length of the strip multiplied by the thickness of the strip 2 i p 2 (elements denoted by green and blue).Summing up, for one block, for a given p and given iteration step I, we obtain a matrix of size bxb and can eliminate a rows, where a = p2 i p = 2 i p 2 and b = 2 i p 2 .So, the cost of elimination of a rows from the matrix of size bxb can be calculated as ab 2 = 2 i p 2 2 2i p 4 = 2 3i p 6 .We perform elimination in the following steps for i = 0, 1, …, s − 1.The cost of performing elimination in step i = 0 is 2 2s p 4 .The computational cost of elimination in all following steps is given by the formula: Summing up, the computational cost of performing eliminations step by step, merging blocks according to our algorithm, is equal to where N is the number of elements and N = O(2 2s p 2 ).Let L(i,p) denote the length of the strip for a given p in step i.The length of the strip can be given by following recursive equation:

Refined Meshes with T-Spline Basis Functions
The solution of this recursive equation is the following: As such, the number of elements which can be eliminated in i-th step for a given p is proportional to p2 i p = 2 i p 2 (thickness multiplied by length of the strip) (elements denoted by green color).In one block, the number of all elements which were not eliminated before is proportional to the length of the strip multiplied by the thickness of the strip 2 i p 2 (elements denoted by green and blue).Summing up, for one block, for a given p and given iteration step I, we obtain a matrix of size bxb and can eliminate a rows, where a = p2 i p = 2 i p 2 and b = 2 i p 2 .So, the cost of elimination of a rows from the matrix of size bxb can be calculated as ab 2 = 2 i p 2 2 2i p 4 = 2 3i p 6 .We perform elimination in the following steps for i = 0, 1, . . ., s − 1.The cost of performing elimination in step i = 0 is 2 2s p 4 .The computational cost of elimination in all following steps is given by the formula: Summing up, the computational cost of performing eliminations step by step, merging blocks according to our algorithm, is equal to O(N 1.5 p 3 ) (4 where N is the number of elements and N = O(2 2s p 2 ).

Refined Meshes with T-Spline Basis Functions
This section presents the estimations of the computational costs of adaptive grids with T-spline basis functions.The grids considered here do not have analysis-suitable T-junction extensions.

Mesh with Point Singularity
The mesh refined toward a point singularity contains layers of elements surrounding the singularity.This kind of mesh is common, e.g., in the case of point sources on the right-hand side.
In this case, each T-spline of order p "crosses" p corresponding edges.In Figure 11, we select one T-spline related to one finite element, for the case of p = 2.We color this element in blue.We also color in green all the elements with T-splines that overlap (have non-zero support intersections) with the T-spline corresponding to the blue element.The number of green elements represents the size of the frontal matrix.To eliminate the blue element (the row related to the T-spline over the blue element), we need to subtract this row from all the rows representing the T-splines denoted by green color.In Figure 12, we plot the analogous situation, this time for p = 4. Due to the fact that the support of elements now spans into four elements in both direction, we have to denote all the elements by green color.

Mesh with Point Singularity
The mesh refined toward a point singularity contains layers of elements surrounding the singularity.This kind of mesh is common, e.g., in the case of point sources on the right-hand side.
In this case, each T-spline of order p "crosses" p corresponding edges.In Figure 11, we select one T-spline related to one finite element, for the case of p = 2.We color this element in blue.We also color in green all the elements with T-splines that overlap (have non-zero support intersections) with the T-spline corresponding to the blue element.The number of green elements represents the size of the frontal matrix.To eliminate the blue element (the row related to the T-spline over the blue element), we need to subtract this row from all the rows representing the T-splines denoted by green color.In Figure 12, we plot the analogous situation, this time for p = 4. Due to the fact that the support of elements now spans into four elements in both direction, we have to denote all the elements by green color.The elimination algorithm is the following.In the first step (i = 0), the algorithm will eliminate the small corner element.In the second step (I = 1), the algorithm will eliminate 3 elements which are neighboring the corner element.In the third step (i = 2), the algorithm will eliminate 3 elements neighboring the 3 elements eliminated in the second step.The algorithm will follow the presented scheme step by step, eliminating, in each step, 3 elements neighboring the three elements eliminated in the previous step.So, for p = 2, in each step i = 0, 1, …, we will eliminate a rows from the matrix bxb, where a = 1 for i = 0 and a = 3 for i ≥ 1.The value of b is equal to 6 + (r − i), where r denotes the number of refinement levels (for the mesh from Figure 12, r equals 4).To sum up, for i ≥ 1, we have a = 3 and b = 6 + (r − i), so the cost of elimination of 3 elements in step i is ab 2 = 3(6 + r − i) 2 .We have r layers.The number of all elements in the mesh N is equal to 3*r + 1.The computational cost of elimination of elements of i-th layer is 3 (6 + r − i) 2 .The computational cost of the whole elimination process is: because N = 3r + 1.For one element, for a given p ≥ 4, we obtain a dense matrix of size bxb, where b is equal to N and can eliminate 1 row.So the computational cost of elimination of one T-spline for p ≥ 4 is O(N 2 ).Summing up, the computational cost of elimination of N elements for p ≥ 4 is The elimination algorithm is the following.In the first step (i = 0), the algorithm will eliminate the small corner element.In the second step (I = 1), the algorithm will eliminate 3 elements which are neighboring the corner element.In the third step (i = 2), the algorithm will eliminate 3 elements neighboring the 3 elements eliminated in the second step.The algorithm will follow the presented scheme step by step, eliminating, in each step, 3 elements neighboring the three elements eliminated in the previous step.So, for p = 2, in each step i = 0, 1, . . ., we will eliminate a rows from the matrix bxb, where a = 1 for i = 0 and a = 3 for i ≥ 1.The value of b is equal to 6 + (r − i), where r denotes the number of refinement levels (for the mesh from Figure 12, r equals 4).To sum up, for i ≥ 1, we have a = 3 and b = 6 + (r − i), so the cost of elimination of 3 elements in step i is ab 2 = 3(6 + r − i) 2 .We have r layers.The number of all elements in the mesh N is equal to 3*r + 1.The computational cost of elimination of elements of i-th layer is 3 (6 + r − i) 2 .The computational cost of the whole elimination process is: because N = 3r + 1.For one element, for a given p ≥ 4, we obtain a dense matrix of size bxb, where b is equal to N and can eliminate 1 row.So the computational cost of elimination of one T-spline for p ≥ 4 is O(N 2 ).Summing up, the computational cost of elimination of N elements for p ≥ 4 is Symmetry 2020, 12, 2070 13 of 26

Mesh with Edge Singularity
Now, we focus on meshes refined towards the edge.This kind of mesh is common in simulations involving, e.g., boundary layers with strong local gradients in one direction.Figure 13 presents an exemplary two-dimensional mesh with edge singularity.For this kind of mesh, each T-spline of order p "crosses" p/2 + 1 corresponding edges.Figure 13 presents the supports of corner T-splines for the cases of p = 2 and p = 4.

Mesh with Edge Singularity
Now, we focus on meshes refined towards the edge.This kind of mesh is common in simulations involving, e.g., boundary layers with strong local gradients in one direction.Figure 13 presents an exemplary two-dimensional mesh with edge singularity.For this kind of mesh, each T-spline of order p "crosses" p/2 + 1 corresponding edges.Figure 13 presents the supports of corner T-splines for the cases of p = 2 and p = 4.In order to eliminate one element, we have to consider p − 1 neighboring elements in each "direction".We introduce the refinement layers as presented in the picture.Notice that the first two layers have elements of identical sizes.For simplicity, we will omit the first layer in our consideration.Each T-spline of the even order can be associated with a single element, where it has its maximum value.When we build the linear equations system resulting from finite element method discretization, in the rows and columns, we have particular T-splines span over the mesh.The entries in the matrix correspond to integrals with multiplications of the T-spline pairs (and their derivatives, depending on the problem we solve).One from a row and one from a column.Thus, non-zero entries denote overlapping supports of pairs of T-splines.In the following figures, for each T-spline whose element is denoted by blue color, we color in green all the T-spline elements whose supports overlap with the blue T-spline (see Figure 14, central and right panel and Figure 15).In this way, we denote the sparsity pattern of the matrix's rows associated with the blue T-splines.
We have r + 1 layers (see Figure 14, where layers are denoted by different colors).The layer number 0 consists of 2 r elements.The i-th layer (i ≥ 1) consists of 2 r−i+1 elements.The number of all elements in the mesh N (equal to the number of T-splines) is equal to:  In order to eliminate one element, we have to consider p − 1 neighboring elements in each "direction".We introduce the refinement layers as presented in the picture.Notice that the first two layers have elements of identical sizes.For simplicity, we will omit the first layer in our consideration.Each T-spline of the even order can be associated with a single element, where it has its maximum value.When we build the linear equations system resulting from finite element method discretization, in the rows and columns, we have particular T-splines span over the mesh.The entries in the matrix correspond to integrals with multiplications of the T-spline pairs (and their derivatives, depending on the problem we solve).One from a row and one from a column.Thus, non-zero entries denote overlapping supports of pairs of T-splines.In the following figures, for each T-spline whose element is denoted by blue color, we color in green all the T-spline elements whose supports overlap with the blue T-spline (see Figure 14, central and right panel and Figure 15).In this way, we denote the sparsity pattern of the matrix's rows associated with the blue T-splines.
We have r + 1 layers (see Figure 14, where layers are denoted by different colors).The layer number 0 consists of 2 r elements.The i-th layer (i ≥ 1) consists of 2 r−i+1 elements.The number of all elements in the mesh N (equal to the number of T-splines) is equal to: Symmetry 2020, 12, x FOR PEER REVIEW 14 of 27

Mesh with Edge Singularity
Now, we focus on meshes refined towards the edge.This kind of mesh is common in simulations involving, e.g., boundary layers with strong local gradients in one direction.Figure 13 presents an exemplary two-dimensional mesh with edge singularity.For this kind of mesh, each T-spline of order p "crosses" p/2 + 1 corresponding edges.Figure 13 presents the supports of corner T-splines for the cases of p = 2 and p = 4.In order to eliminate one element, we have to consider p − 1 neighboring elements in each "direction".We introduce the refinement layers as presented in the picture.Notice that the first two layers have elements of identical sizes.For simplicity, we will omit the first layer in our consideration.Each T-spline of the even order can be associated with a single element, where it has its maximum value.When we build the linear equations system resulting from finite element method discretization, in the rows and columns, we have particular T-splines span over the mesh.The entries in the matrix correspond to integrals with multiplications of the T-spline pairs (and their derivatives, depending on the problem we solve).One from a row and one from a column.Thus, non-zero entries denote overlapping supports of pairs of T-splines.In the following figures, for each T-spline whose element is denoted by blue color, we color in green all the T-spline elements whose supports overlap with the blue T-spline (see Figure 14, central and right panel and Figure 15).In this way, we denote the sparsity pattern of the matrix's rows associated with the blue T-splines.
We have r + 1 layers (see Figure 14, where layers are denoted by different colors).The layer number 0 consists of 2 r elements.The i-th layer (i ≥ 1) consists of 2 r−i+1 elements.The number of all elements in the mesh N (equal to the number of T-splines) is equal to:   In the first step (i = 1), the algorithm will consider the number of blocks depending on the number of refinement levels in the mesh, denoted by r, and the order of approximation, denoted by p.The size of each block depends only on the order of approximation p.Then, there are the following steps for constructing blocks and eliminating selected T-splines from the blocks.We follow the layers of the mesh.For simplicity, the layer called i = 0 will be omitted.In the first step, we construct blocks related to the first layer of elements, where we group elements in the following way.We take p elements from the first layer and we denote them by blue color.Later, we search for the largest elements whose T-spline overlaps with the blue elements from the first layer.We construct a rectangular block spanning up to the largest blocks overlapping with our blue elements.In this way, we have all the "columns" of T-splines related to the non-zero entries in the rows on the blue elements in the block.For p = 2 and for i = 1, we will construct a block consisting of six elements from the first layer and three elements from the second layer.For p = 4 and for i = 1, we will construct a block consisting of 20 elements from the first layer, ten elements from the second layer, and five elements from the third layer.We noticed that the supports of the p central T-splines from two adjacent blocks do not overlap.In the next step, we construct a new, bigger rectangular block consisting of two adjacent blocks from the previous step and the next layer's corresponding elements.For p = 2 and for i = 2, we will construct a block consisting of 12 elements from the first layer, six elements from the second layer, and three elements from the third layer.For p = 4 and for i = 2, we will construct a block consisting of 40 elements from the first layer, 20 elements from the second layer, ten elements from the third layer, and five elements from the fourth layer (see Figure 16).The size of the block in the i-th step (the number of elements in the block from the first layer of elements) is equal to n1 = (p + 1) 2 (p/2) 2 (i−1) = (p + 1)2 (p/2+i−1) (where (p + 1) denotes the number of elements on the largest layer in the block, (p + 1)2 (p/2) is the number of elements in the first layer in the block from the first iteration i = 1, and (p + 1)2 (p/2) 2 (i−1) is the number of elements in the first layer in the block from iteration i).In the first step (i = 1), the algorithm will consider the number of blocks depending on the number of refinement levels in the mesh, denoted by r, and the order of approximation, denoted by p.The size of each block depends only on the order of approximation p.Then, there are the following steps for constructing blocks and eliminating selected T-splines from the blocks.We follow the layers of the mesh.For simplicity, the layer called i = 0 will be omitted.In the first step, we construct blocks related to the first layer of elements, where we group elements in the following way.We take p elements from the first layer and we denote them by blue color.Later, we search for the largest elements whose T-spline overlaps with the blue elements from the first layer.We construct a rectangular block spanning up to the largest blocks overlapping with our blue elements.In this way, we have all the "columns" of T-splines related to the non-zero entries in the rows on the blue elements in the block.For p = 2 and for i = 1, we will construct a block consisting of six elements from the first layer and three elements from the second layer.For p = 4 and for i = 1, we will construct a block consisting of 20 elements from the first layer, ten elements from the second layer, and five elements from the third layer.We noticed that the supports of the p central T-splines from two adjacent blocks do not overlap.In the next step, we construct a new, bigger rectangular block consisting of two adjacent blocks from the previous step and the next layer's corresponding elements.For p = 2 and for i = 2, we will construct a block consisting of 12 elements from the first layer, six elements from the second layer, and three elements from the third layer.For p = 4 and for i = 2, we will construct a block consisting of 40 elements from the first layer, 20 elements from the second layer, ten elements from the third layer, and five elements from the fourth layer (see Figure 16).The size of the block in the i-th step (the number of elements in the block from the first layer of elements) is equal to n1 = (p + 1) 2 (p/2) 2 (i−1) = (p + 1)2 (p/2+i−1) (where (p + 1) denotes the number of elements on the largest layer in the block, (p + 1)2 (p/2) is the number of elements in the first layer in the block from the first iteration i = 1, and (p + 1)2 (p/2) 2 (i−1) is the number of elements in the first layer in the block from iteration i).In the first step (i = 1), the algorithm will consider the number of blocks depending on the number of refinement levels in the mesh, denoted by r, and the order of approximation, denoted by p.The size of each block depends only on the order of approximation p.Then, there are the following steps for constructing blocks and eliminating selected T-splines from the blocks.We follow the layers of the mesh.For simplicity, the layer called i = 0 will be omitted.In the first step, we construct blocks related to the first layer of elements, where we group elements in the following way.We take p elements from the first layer and we denote them by blue color.Later, we search for the largest elements whose T-spline overlaps with the blue elements from the first layer.We construct a rectangular block spanning up to the largest blocks overlapping with our blue elements.In this way, we have all the "columns" of T-splines related to the non-zero entries in the rows on the blue elements in the block.For p = 2 and for i = 1, we will construct a block consisting of six elements from the first layer and three elements from the second layer.For p = 4 and for i = 1, we will construct a block consisting of 20 elements from the first layer, ten elements from the second layer, and five elements from the third layer.We noticed that the supports of the p central T-splines from two adjacent blocks do not overlap.In the next step, we construct a new, bigger rectangular block consisting of two adjacent blocks from the previous step and the next layer's corresponding elements.For p = 2 and for i = 2, we will construct a block consisting of 12 elements from the first layer, six elements from the second layer, and three elements from the third layer.For p = 4 and for i = 2, we will construct a block consisting of 40 elements from the first layer, 20 elements from the second layer, ten elements from the third layer, and five elements from the fourth layer (see Figure 16).The size of the block in the i-th step (the number of elements in the block from the first layer of elements) is equal to n1 = (p + 1) 2 (p/2) 2 (i−1) = (p + 1)2 (p/2+i−1) (where (p + 1) denotes the number of elements on the largest layer in the block, (p + 1)2 (p/2) is the number of elements in the first layer in the block from the first iteration i = 1, and (p + 1)2 (p/2) 2 (i−1) is the number of elements in the first layer in the block from iteration i).The number of elements in the first layer (i = 1) is equal to 2 r .The number of blocks constructed in the i-th step is equal to: The computational cost of elimination of T-splines from one block in the i-th step is equal to ab 2 , where a and b are defined as: (this is because in one considered layer, we have 2 (p/2) elements to eliminate, and in the other layers, we have a total of (2 (p/2) p(i-1)) elements to eliminate) (since we have p + 1 elements in one layer, and two times more in the previous layers, so the total number is (p + 1)(1 + 2 + 2 2 + … + 2 p/2 )).Thus, we have: The computational cost of the whole elimination process is: since the term Σi=1,...,r 2 −i i 3 converges to the constant.

Refined Meshes with Analysis Suitable T-Splines
This section presents the estimations of the computational costs of adaptive grids with analysis-suitable T-spline basis functions.The grids considered here have analysis-suitable T-junction extensions.The number of elements in the first layer (i = 1) is equal to 2 r .The number of blocks constructed in the i-th step is equal to: The computational cost of elimination of T-splines from one block in the i-th step is equal to ab 2 , where a and b are defined as: (this is because in one considered layer, we have 2 (p/2) elements to eliminate, and in the other layers, we have a total of (2 (p/2) p(i-1)) elements to eliminate) (since we have p + 1 elements in one layer, and two times more in the previous layers, so the total number is (p + 1)(1 + 2 + 2 2 + . . .+ 2 p/2 )).Thus, we have: The computational cost of the whole elimination process is: since the term Σ i=1, . . .,r 2 −i i 3 converges to the constant.

Refined Meshes with Analysis Suitable T-Splines
This section presents the estimations of the computational costs of adaptive grids with analysis-suitable T-spline basis functions.The grids considered here have analysis-suitable T-junction extensions.We color this element in blue.We also color in green all the elements with T-splines that overlap (have non-zero support intersections) with the T-spline corresponding to the blue element.The number of green elements represents the size of the frontal matrix.In Figures 22-25, we plot the analogous situation, this time for p = 4. Let us start with the mesh from Figure 17 and p = 2.In the first step (i = 0), the algorithm will eliminate the corner element.In the second step (i = 1), the algorithm will eliminate three elements which are neighboring the corner element.In the third step (i = 2), the algorithm will eliminate five elements neighboring the three elements eliminated in the second step.In the fourth step (i = 3), the algorithm will eliminate five elements neighboring the five elements eliminated in the third step.

Mesh with Point Singularity
Symmetry 2020, 12, x FOR PEER REVIEW 17 of 27

Mesh with Point Singularity
Figure 17 presents an exemplary 2D mesh with T-junctions with a point singularity for p = 2 and p = 4.In Figures 18-21, we select one T-spline related to one finite element, for the case of p = 2.We color this element in blue.We also color in green all the elements with T-splines that overlap (have non-zero support intersections) with the T-spline corresponding to the blue element.The number of green elements represents the size of the frontal matrix.In Figures 22-25, we plot the analogous situation, this time for p = 4. Let us start with the mesh from Figure 17 and p = 2.In the first step (i = 0), the algorithm will eliminate the corner element.In the second step (i = 1), the algorithm will eliminate three elements which are neighboring the corner element.In the third step (i = 2), the algorithm will eliminate five elements neighboring the three elements eliminated in the second step.In the fourth step (i = 3), the algorithm will eliminate five elements neighboring the five elements eliminated in the third step.

Mesh with Point Singularity
Figure 17 presents an exemplary 2D mesh with T-junctions with a point singularity for p = 2 and p = 4.In Figures 18-21, we select one T-spline related to one finite element, for the case of p = 2.We color this element in blue.We also color in green all the elements with T-splines that overlap (have non-zero support intersections) with the T-spline corresponding to the blue element.The number of green elements represents the size of the frontal matrix.In Figures 22-25, we plot the analogous situation, this time for p = 4. Let us start with the mesh from Figure 17 and p = 2.In the first step (i = 0), the algorithm will eliminate the corner element.In the second step (i = 1), the algorithm will eliminate three elements which are neighboring the corner element.In the third step (i = 2), the algorithm will eliminate five elements neighboring the three elements eliminated in the second step.In the fourth step (i = 3), the algorithm will eliminate five elements neighboring the five elements eliminated in the third step.

Mesh with Point Singularity
Figure 17 presents an exemplary 2D mesh with T-junctions with a point singularity for p = 2 and p = 4.In Figures 18-21, we select one T-spline related to one finite element, for the case of p = 2.We color this element in blue.We also color in green all the elements with T-splines that overlap (have non-zero support intersections) with the T-spline corresponding to the blue element.The number of green elements represents the size of the frontal matrix.In Figures 22-25, we plot the analogous situation, this time for p = 4. Let us start with the mesh from Figure 17 and p = 2.In the first step (i = 0), the algorithm will eliminate the corner element.In the second step (i = 1), the algorithm will eliminate three elements which are neighboring the corner element.In the third step (i = 2), the algorithm will eliminate five elements neighboring the three elements eliminated in the second step.In the fourth step (i = 3), the algorithm will eliminate five elements neighboring the five elements eliminated in the third step.The algorithm will follow the presented scheme step by step, eliminating, in each step, five elements neighboring the five elements eliminated in the previous step.In a similar way, for p = 4, the algorithm will eliminate one element for step i = 0, three elements for step i = 1, five elements for step i = 2, and seven elements in each step i ≥ 3.
To sum up, we will eliminate T-splines, layer by layer, starting from the corner element (left top).The cost of elimination of each element can be estimated by checking the number of "layers" of other elements which have assigned T-splines who overlap elements from the current layer, by determining the support of eliminated elements (p/2) and the number of elements in such a "layer" (p + 3).So, for one layer, we have: The computational cost of elimination of T-splines for one layer is: We have r layers (the refinement level is r).The number of all elements in the mesh N is equal to: The computational cost of elimination of elements of one layer is O(p 5 ).The computational cost of the whole elimination process is:

Mesh with Edge Singularity
Figure 26 presents an exemplary 2D mesh with T-junctions with an edge singularity for p = 2 and p = 4.There are the following steps for constructing blocks and eliminating selected T-splines from the blocks.For simplicity, the first p/2 + 1 layers will be omitted.In the first step, we construct blocks related to the first layer of elements, where we group elements in the following way.The algorithm will follow the presented scheme step by step, eliminating, in each step, five elements neighboring the five elements eliminated in the previous step.In a similar way, for p = 4, the algorithm will eliminate one element for step i = 0, three elements for step i = 1, five elements for step i = 2, and seven elements in each step i ≥ 3.
To sum up, we will eliminate T-splines, layer by layer, starting from the corner element (left top).The cost of elimination of each element can be estimated by checking the number of "layers" of other elements which have assigned T-splines who overlap elements from the current layer, by determining the support of eliminated elements (p/2) and the number of elements in such a "layer" (p + 3).So, for one layer, we have: The computational cost of elimination of T-splines for one layer is: We have r layers (the refinement level is r).The number of all elements in the mesh N is equal to: The computational cost of elimination of elements of one layer is O(p 5 ).The computational cost of the whole elimination process is:

Mesh with Edge Singularity
Figure 26 presents an exemplary 2D mesh with T-junctions with an edge singularity for p = 2 and p = 4.There are the following steps for constructing blocks and eliminating selected T-splines from the blocks.For simplicity, the first p/2 + 1 layers will be omitted.In the first step, we construct blocks related to the first layer of elements, where we group elements in the following way.
We take p elements from the first layer and we denote them by green color.Later, we search for the largest T-spline element overlaps with the green elements from the first layer.We construct a rectangular block spanning up to the largest blocks overlapping with our green elements.In this way, we have all the "columns" of T-splines related to the non-zero entries in the rows of the blue elements in the block.For p = 2 and for i = 1, we will construct a block consisting of six elements from the first layer and three elements from the second layer.For p = 4 and for i = 1, we will construct a block consisting of 20 elements from the first layer, ten elements from the second layer, and five elements from the third layer (see the right panel in Figure 26).We take p elements from the first layer and we denote them by green color.Later, we search for the largest T-spline element overlaps with the green elements from the first layer.We construct a rectangular block spanning up to the largest blocks overlapping with our green elements.In this way, we have all the "columns" of T-splines related to the non-zero entries in the rows of the blue elements in the block.For p = 2 and for i = 1, we will construct a block consisting of six elements from the first layer and three elements from the second layer.For p = 4 and for i = 1, we will construct a block consisting of 20 elements from the first layer, ten elements from the second layer, and five elements from the third layer (see the right panel in Figure 26).
We noticed that the supports of the p central T-splines from two adjacent blocks do not overlap.In the next step, we construct a new, bigger rectangular block consisting of two adjacent blocks from the previous step and the next layer's corresponding elements.For p = 2 and for i = 2, we will construct a block composed of 12 elements from the first layer, six elements from the second layer, and three elements from the third layer.For p = 4 and for i = 2, we will construct a block consisting of 40 elements from the first layer, 20 elements from the second layer, ten elements from the third layer, and five elements from the fourth layer (see the right panel in Figure 27).The size of the block in the i-th step (the number of elements in the first column of elements in the block) is equal to (p + 1): The number of all elements in the first layer (i = 1) is equal to 2 r .As such, the number of blocks constructed in the i-th step is equal to: We noticed that the supports of the p central T-splines from two adjacent blocks do not overlap.In the next step, we construct a new, bigger rectangular block consisting of two adjacent blocks from the previous step and the next layer's corresponding elements.For p = 2 and for i = 2, we will construct a block composed of 12 elements from the first layer, six elements from the second layer, and three elements from the third layer.For p = 4 and for i = 2, we will construct a block consisting of 40 elements from the first layer, 20 elements from the second layer, ten elements from the third layer, and five elements from the fourth layer (see the right panel in Figure 27).The size of the block in the i-th step (the number of elements in the first column of elements in the block) is equal to (p + 1): The number of all elements in the first layer (i = 1) is equal to 2 r .As such, the number of blocks constructed in the i-th step is equal to: We take p elements from the first layer and we denote them by green color.Later, we search for the largest T-spline element overlaps with the green elements from the first layer.We construct a rectangular block spanning up to the largest blocks overlapping with our green elements.In this way, we have all the "columns" of T-splines related to the non-zero entries in the rows of the blue elements in the block.For p = 2 and for i = 1, we will construct a block consisting of six elements from the first layer and three elements from the second layer.For p = 4 and for i = 1, we will construct a block consisting of 20 elements from the first layer, ten elements from the second layer, and five elements from the third layer (see the right panel in Figure 26).
We noticed that the supports of the p central T-splines from two adjacent blocks do not overlap.In the next step, we construct a new, bigger rectangular block consisting of two adjacent blocks from the previous step and the next layer's corresponding elements.For p = 2 and for i = 2, we will construct a block composed of 12 elements from the first layer, six elements from the second layer, and three elements from the third layer.For p = 4 and for i = 2, we will construct a block consisting of 40 elements from the first layer, 20 elements from the second layer, ten elements from the third layer, and five elements from the fourth layer (see the right panel in Figure 27).The size of the block in the i-th step (the number of elements in the first column of elements in the block) is equal to (p + 1): 2 (p/2) 2 (i−1) = (p + 1)2 2 (p/2+i−1) (17 The number of all elements in the first layer (i = 1) is equal to 2 r .As such, the number of blocks constructed in the i-th step is equal to:  The number of all elements in one layer (i-th layer) is equal to 2 r−i .Thus, the number of all elements in the considered layers is equal to: N = Σ i=1, . . .,r−p/2+1 (2 r−i ) = O(2 r ) ( The computational cost of elimination of T-splines from one block in the i-th step is equal to ab 2 , where a and b are defined as: a = 2 p/2 + 2 (p/2) p(i-1) = O(2 p/2 pi) (this is because in the actual layer, we have 2 (p/2) elements, and in the other layers, we have 2 (p/2) p(i-1) elements) b = (p + 1)(1 + 2 + . . .+ p/2) + 2 p/2 + 2 (p/2) 2p(i-1) = O(2 (p/2) pi) (since we have p + 1 elements in one layer, and this number multiplies by 2 in each following layer, we have (p + 1)(1 + 2 + 2 2 . . . 2 p/2 )).Thus.

Numerical Results
In this section we present the numerical results, obtained in Octave, confirming the estimated computational costs for T-splines and analysis-suitable T-splines with our ordering algorithms.They are summarized in Figures 27-32.For T-splines, the matrices are dense and cost of factorization is O(N 3 ).For analysis-suitable T-splines, for grids refined towards a point or an edge, our ordering delivers linear computational cost of factorization.For analysis-suitable T-splines, for grids refined towards an edge, our recursive ordering algorithm delivers a linear computational cost, which is up to 50-times faster than alternative Octave orderings.

Conclusions
In this paper, we estimated, for the first time, the computational costs of the multi-frontal solver when using T-spline and analysis-suitable T-splines on uniform and adaptive grids.We have shown that analysis-suitable T-splines have a linear computational cost on adaptive grids.We also proposed an ordering algorithm permuting the rows and columns of the resulting sparse matrix such that the computational cost of factorization is linear for grids refined towards point and edge singularities when using analysis-suitable T-splines.Our ordering algorithm constructs a binary tree, with nodes assigned to elements on levels of the refinement.The root is the largest element, and the leaves are the smallest elements.It orders the elements and T-splines assigned to them by post-order traversal of the tree.
In this paper, we compared computational costs of isogeometric analysis with T-splines and analysis-suitable T-splines.We focused on model two-dimensional symmetry-preserving grids, including a uniform grid, a grid refined to a point, and a grid refined to an edge.They are building blocks for adaptive two-dimensional grids.We conclude that the cost for a uniform grid is O(N 1.5 p 3 ) for both T-splines and analysis-suitable T-splines, as well as for B-splines (since supports of T-splines are equivalent to the supports of B-splines on uniform grids).The number of unknowns on the uniform grid is defined as N = O(2 2s p 2 ), where 2 2s is the number of elements and p denotes the T-splines (or B-splines) order.The cost for grids refined towards a point is O(N 3 ) for T-splines and O(Np 4 ) for analysis-suitable T-splines.For T-splines, the number of unknowns is N = 3r + 1 = O(r), where r is the number of layers in the mesh.For analysis-suitable T-splines, the number of unknowns is N = O(rp).Finally, the costs for grids refined towards an edge is O(N2 p p 2 ) for both T-splines and

Figure 1 .Figure 1 .
Figure 1.C 0 B-spline basis functions defined by tensor product of knot vectors [0 0 0 1 2 2 2] x [0 0 0 1 2 2 2]We enumerate the B-spline basis functions from 1 to 25, since we have five one-dimensional B-splines in each direction, and the tensor product results in 5*5=25 basis functions.We construct a matrix (see Figure2) where rows and columns correspond to the B-splines, and the non-zero values

Figure 3 .Figure 3 .
Figure 3. Left panel: The separators found by the nested dissections algorithm, using the peripheral node 23 as the starting node.Right panel: The elimination tree.

Figure 4 .
Figure 4. (Left panel) T-splines defined over a grid refined towards a point.(Right panel) Analysis-suitable T-splines defined over a grid refined towards a point.

Figure 4 .
Figure 4. (Left panel) T-splines defined over a grid refined towards a point.(Right panel) Analysis-suitable T-splines defined over a grid refined towards a point.

Figure 5 .
Figure 5. Ordering of elements = T-splines = rows of matrices for grids refined towards a point and an edge.

Figure 5 .
Figure 5. Ordering of elements = T-splines = rows of matrices for grids refined towards a point and an edge.

Figure 7 .
Figure 7. Left panel: Iteration step i = 0, p = 2, and s = 3.The elements for which the T-splines can be eliminated are denoted by green color.Right panel: Iteration step i = 1, p = 2, and s = 3.The elements for which the T-splines can be eliminated are denoted by light green color, elements for which the T-splines were eliminated in previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 8 .
Figure 8. Iteration step i = 2, p = 2, and s = 3.The elements for which the T-splines can be eliminated are denoted by dark green color, elements for which the T-splines were eliminated in the previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 7 . 27 Figure 7 .
Figure 7. (Left panel): Iteration step i = 0, p = 2, and s = 3.The elements for which the T-splines can be eliminated are denoted by green color.(Right panel): Iteration step i = 1, p = 2, and s = 3.The elements for which the T-splines can be eliminated are denoted by light green color, elements for which the T-splines were eliminated in previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 8 .
Figure 8. Iteration step i = 2, p = 2, and s = 3.The elements for which the T-splines can be eliminated are denoted by dark green color, elements for which the T-splines were eliminated in the previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 8 . 27 Figure 7 .
Figure 8. Iteration step i = 2, p = 2, and s = 3.The elements for which the T-splines can be eliminated are denoted by dark green color, elements for which the T-splines were eliminated in the previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 8 .
Figure 8. Iteration step i = 2, p = 2, and s = 3.The elements for which the T-splines can be eliminated are denoted by dark green color, elements for which the T-splines were eliminated in the previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 9 .
Figure 9. Iteration step i = 0 and I = 1, p = 4.The elements for which the T-splines can be eliminated are denoted by green color, elements for which the T-splines were eliminated in previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 9 .
Figure 9. Iteration step i = 0 and I = 1, p = 4.The elements for which the T-splines can be eliminated are denoted by green color, elements for which the T-splines were eliminated in previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 10 .
Figure 10.Iteration step i = 2 and I = 3, p = 4.The elements for which the T-splines can be eliminated are denoted by green color, elements for which the T-splines were eliminated in previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 10 .
Figure 10.Iteration step i = 2 and I = 3, p = 4.The elements for which the T-splines can be eliminated are denoted by green color, elements for which the T-splines were eliminated in previous steps are denoted by white color, and the rest of elements are denoted by blue color.

Figure 11 .
Figure 11.We select a single element and its T-spline and we denote it by blue color.Next, we browse all equal-size and larger elements, and we mark in green these elements whose T-splines have non-zero intersections with the T-spline span over the blue element.We consider quadratic T-splines here.

Figure 11 .
Figure 11.We select a single element and its T-spline and we denote it by blue color.Next, we browse all equal-size and larger elements, and we mark in green these elements whose T-splines have non-zero intersections with the T-spline span over the blue element.We consider quadratic T-splines here.

Figure 12 .
Figure 12.The elements with T-splines having non-zero intersection with the T-spline corresponding to the element, denoted by blue color, are denoted by green color, for p = 4.

Figure 12 .
Figure 12.The elements with T-splines having non-zero intersection with the T-spline corresponding to the element, denoted by blue color, are denoted by green color, for p = 4.

Figure 13 .
Figure 13.Two-dimensional mesh with edge singularity.Supports of T-splines of order p = 2 (left panel) and p = 4 (right panel) on meshes with edge singularity.

Figure 14 .
Figure 14.Left panel: Layers of mesh refined toward the edge.Central panel and right panel: The sparsity pattern of the rows of the matrix associated with the blue T-splines from the first and the second considered column, for p = 2.

Figure 13 .
Figure 13.Two-dimensional mesh with edge singularity.Supports of T-splines of order p = 2 (left panel) and p = 4 (right panel) on meshes with edge singularity.

Figure 13 .
Figure 13.Two-dimensional mesh with edge singularity.Supports of T-splines of order p = 2 (left panel) and p = 4 (right panel) on meshes with edge singularity.

Figure 14 .
Figure 14.Left panel: Layers of mesh refined toward the edge.Central panel and right panel: The sparsity pattern of the rows of the matrix associated with the blue T-splines from the first and the second considered column, for p = 2.

Figure 14 .
Figure 14.(Left panel): Layers of mesh refined toward the edge.Central panel and (right panel): The sparsity pattern of the rows of the matrix associated with the blue T-splines from the first and the second considered column, for p = 2.

Figure 15 .
Figure 15.The support of the T-spline related with an element from the first (left panel), second (central panel), or third layer (right panel), denoted by blue color, overlapping with supports of other T-splines, denoted by green color, for p = 2, 4, and 6, respectively.

Figure 15 .
Figure 15.The support of the T-spline related with an element from the first (left panel), second (central panel), or third layer (right panel), denoted by blue color, overlapping with supports of other T-splines, denoted by green color, for p = 2, 4, and 6, respectively.

Symmetry 2020 , 27 Figure 15 .
Figure 15.The support of the T-spline related with an element from the first (left panel), second (central panel), or third layer (right panel), denoted by blue color, overlapping with supports of other T-splines, denoted by green color, for p = 2, 4, and 6, respectively.

Figure 16 .
Figure 16.Eliminated blocks for p = 2 and 4 for the first, second, and third iteration steps.

Figure 16 .
Figure 16.Eliminated blocks for p = 2 and 4 for the first, second, and third iteration steps.

Figure 17
Figure 17 presents an exemplary 2D mesh with T-junctions with a point singularity for p = 2 and p = 4.In Figures 18-21, we select one T-spline related to one finite element, for the case of p = 2.We color this element in blue.We also color in green all the elements with T-splines that overlap (have non-zero support intersections) with the T-spline corresponding to the blue element.The number of green elements represents the size of the frontal matrix.In Figures22-25, we plot the analogous situation, this time for p = 4. Let us start with the mesh from Figure17and p = 2.In the first step (i = 0), the algorithm will eliminate the corner element.In the second step (i = 1), the algorithm will eliminate three elements which are neighboring the corner element.In the third step (i = 2), the algorithm will eliminate five elements neighboring the three elements eliminated in the second step.In the fourth step (i = 3), the algorithm will eliminate five elements neighboring the five elements eliminated in the third step.

Figure 17 .
Figure 17.The exemplary 2D meshes with a point singularity and T-junction extensions for p = 2 and p = 4 (left and right, respectively).

Figure 17 .
Figure 17.The exemplary 2D meshes with a point singularity and T-junction extensions for p = 2 and p = 4 (left and right, respectively).

Figure 17 .
Figure 17.The exemplary 2D meshes with a point singularity and T-junction extensions for p = 2 and p = 4 (left and right, respectively).

Figure 17 .
Figure 17.The exemplary 2D meshes with a point singularity and T-junction extensions for p = 2 and p = 4 (left and right, respectively).

2 18 ) 27 Figure 27 .
Figure 27.Exemplary blocks for p = 2 and p = 4 for the first iteration step (upper panel).Exemplary blocks for p = 2 and the first, second, and third iteration steps (bottom panel).The number of all elements in one layer (i-th layer) is equal to 2 r−i .Thus, the number of all elements in the considered layers is equal to: N = Σi=1,...,r−p/2+1(2 r−i ) = O(2 r )(19)

Figure 27 .
Figure 27.Exemplary blocks for p = 2 and p = 4 for the first iteration step (upper panel).Exemplary blocks for p = 2 and the first, second, and third iteration steps (bottom panel).

Symmetry 2020 , 27 Figure 28 .
Figure 28.Execution times and sparsity of a matrix for a mesh with T-splines, p = 2, refined towards a point.

Figure 28 .
Figure 28.Execution times and sparsity of a matrix for a mesh with T-splines, p = 2, refined towards a point.

Figure 28 .
Figure 28.Execution times and sparsity of a matrix for a mesh with T-splines, p = 2, refined towards a point.

Figure 29 .
Figure 29.Execution times and sparsity of a matrix for a mesh with T-splines, p = 4, refined towards a point.

Figure 29 .
Figure 29.Execution times and sparsity of a matrix for a mesh with T-splines, p = 4, refined towards a point.

Figure 28 .
Figure 28.Execution times and sparsity of a matrix for a mesh with T-splines, p = 2, refined towards a point.

Figure 29 .
Figure 29.Execution times and sparsity of a matrix for a mesh with T-splines, p = 4, refined towards a point.

Figure 30 . 27 Figure 30 .
Figure 30.Execution times and sparsity of a matrix for a mesh with analysis-suitable T-splines, p = 2, refined towards a point, permuted with our ordering.

Figure 31 .
Figure 31.Execution times and sparsity of a matrix for a mesh with analysis-suitable T-splines, p = 4, refined towards a point, permuted with our ordering.

Figure 31 .
Figure 31.Execution times and sparsity of a matrix for a mesh with analysis-suitable T-splines, p = 4, refined towards a point, permuted with our ordering.

Figure 31 .
Figure 31.Execution times and sparsity of a matrix for a mesh with analysis-suitable T-splines, p = 4, refined towards a point, permuted with our ordering.

Figure 32 .
Figure 32.Comparison of Approximate Minimum Degree (AMD) and Cuthill-McKee orderings, available in Octave, and our ordering for a matrix for a mesh with analysis-suitable T-splines refined towards an edge, p = 2, permuted with our ordering.

Figure 32 .
Figure 32.Comparison of Approximate Minimum Degree (AMD) and Cuthill-McKee orderings, available in Octave, and our ordering for a matrix for a mesh with analysis-suitable T-splines refined towards an edge, p = 2, permuted with our ordering.