Tài liệu Assembly_matthias

Assembly strategies in isogeometric analysis Matthias M¨ller o March 16, 2015 Abstract These short notes shed some light on the practical assembly of bilinear and linear forms in Galerkin isogeometric finite element methods making use of B-spline basis functions. In particular, we present a quadrature-based algorithm for the efficient assembly of the stiffness matrix and the load vector of the Poisson problem in an arbitrary two-dimensional domain. The basic concepts are derived in a most general way and illustrated by a practical example which is repeatedly revisited. Version 2 1 Problem description Consider the two-dimensional Poisson problem −∆u = f in Ω ⊂ R2 , (1a) u=0 on Γ = ∂Ω, (1b) with a Lipschitz continuous domain Ω and assuming homogeneous Dirichlet boundary conditions to be prescribed along the entire boundary. Non-homogeneous Dirichlet boundary conditions can be easily realized via a so-called Dirichlet lift, whereas natural Neumann or Robin boundary conditions are implemented into the boundary integral terms that arise from performing integration by parts in the bilinear form (see Sections 3.3-3.4 in [1]). Here, purely homogeneous boundary conditions are assumed. The variational formulation reads: find the solution in the trial space u ∈ V := {w ∈ H 1 (Ω), w = 0 on Γ} (2a) a(w, u) = L(w) ∀w ∈ V, (2b) such that whereby the bilinear and linear forms are defined as follows: a(w, u) := ( w) u dx, (2c) Ω L(w) := wf dx. (2d) Ω Here, we adopt the notation w u = w · u for the dot-product. Let Ξ = {ξ1 , ξ2 , . . . , ξn+p+1 } and H = {η1 , η2 , . . . , ηm+q+1 } denote two open knot vectors (i.e., multiplicity of first and last knot is p + 1 and q + 1, respectively), where p and q stand for the degree of the B-spline basis functions. Example: As an example, that is used throughout these notes, consider the case p = q = 2 and let Ξ = {0, 0, 0, 0.5, 1, 1, 1}, (3a) H = {0, 0, 0, 1, 1, 1}. (3b) The resulting parametric domain Ω0 = (0, 1)2 is depicted in Figure 1. In order to represent the physical 1 η4 = η5 = η6 = 1 η1 = η2 = η3 = 0 ξ1 = ξ2 = ξ3 = 0 ξ4 = 0.5 ξ5 = ξ6 = ξ7 = 1 Figure 1: Parametric domain Ω0 with Ξ = {0, 0, 0, 0.5, 1, 1, 1} and H = {0, 0, 0, 1, 1, 1}. domain Ω from Figure 2 let the geometric mapping be defined as n m Ni,p (ξ)Mj,q (η)Bij . S(ξ) = (4) i=1 j=1 Here, ξ = (ξ, η) is the pair of parameter values in the parameter space Ω0 , and Ni,p and Mj,q stand for the univariate (=one-dimensional) B-spline basis functions of degree p and q, respectively. Moreover, Bij ∈ Rd are the control points in the d-dimensional physical space. For the two-dimensional (d = 2) y x physical domain depicted in Figure 2 the values of the control points Bij = (Bij , Bij ) are given in Table 1 (cf. Table 2.1 from [1]). It should be noted that using three-dimensional (d = 3) control points y x z Bij = (Bij , Bij , Bij ) in the bivariate geometric mapping (4) would create a two-dimensional manifold in R3 . This concept is commonly used to solve PDE-problems on surfaces (e.g., flow on the earth). (3,5) x = (x, y) (3,1.5) (-2,0) (0,0) Figure 2: Physical domain Ω: Biquadratic B-spline surface with Ξ = {0, 0, 0, 0.5, 1, 1, 1} and H = {0, 0, 0, 1, 1, 1} and control points defined in Table 1. 2 Transformation to parametric domain In order to solve problem (2b) numerically, the integrals needs to be transformed from the physical domain Ω into parametric domain Ω0 first. For this purpose, we will use the well-known integration rule w(S(ξ)) |det DS(ξ)| dξ, w(x) dx = Ω (5a) Ω0 where DS(ξ) = ∂S1 ∂ξ ∂S2 ∂ξ 2 ∂S1 ∂η ∂S2 ∂η (5b) i 1 1 1 2 2 2 3 3 3 4 4 4 j 1 2 3 1 2 3 1 2 3 1 2 3 Bij (0,0) (-1,0) (-2,0) (0,1) (-1,2) (-2,2) (1,1.5) (1,4) (1,5) (3,1.5) (3,4) (3,5) Table 1: Example: control points for the biquadratic B-spline surface depicted in Figure 2. is the Jacobian matrix of the geometric mapping (4). The individual entries are given by n m ∂S1 x = Ni,p (ξ)Mj,q (η)Bij , ∂ξ i=1 j=1 n m ∂S1 x = Ni,p (ξ)Mj,q (η)Bij , ∂η i=1 j=1 n (6b) m ∂S2 y Ni,p (ξ)Mj,q (η)Bij , = ∂ξ i=1 j=1 n (6a) (6c) m ∂S2 y = Ni,p (ξ)Mj,q (η)Bij , ∂η i=1 j=1 (6d) whereby the derivatives of the B-spline basis functions equal p p Ni,p−1 (ξ) + Ni+1,p−1 (ξ), ξi+p − ξi ξi+p+1 − ξi+1 q q Mj,q (η) = Mj,q−1 (η) + Mj+1,q−1 (η). ηj+q − ηj ηj+q+1 − ηj+1 Ni,p (ξ) = (7a) (7b) Hence, the entries of the Jacobian matrix of the geometric mapping read n m n m p ∂S1 p x x = Ni,p−1 (ξ)Mj,q (η)Bij + Ni+1,p−1 (ξ)Mj,q (η)Bij , ∂ξ ξi+p − ξi ξi+p+1 − ξi+1 i=1 j=1 i=1 j=1 n m n m ∂S1 q q x x = Ni,p (ξ)Mj,q−1 (η)Bij + Ni,p (ξ)Mj+1,q−1 (η)Bij , ∂η ηj+q − ηj ηj+q+1 − ηj+1 i=1 j=1 i=1 j=1 n m n m n (8b) m p p ∂S2 y x = Ni,p−1 (ξ)Mj,q (η)Bij + Ni+1,p−1 (ξ)Mj,q (η)Bij , ∂ξ ξi+p − ξi ξi+p+1 − ξi+1 i=1 j=1 i=1 j=1 n (8a) (8c) m ∂S2 q q y x = Ni,p (ξ)Mj,q−1 (η)Bij + Ni,p (ξ)Mj+1,q−1 (η)Bij . ∂η ηj+q − ηj ηj+q+1 − ηj+1 i=1 j=1 i=1 j=1 (8d) Moreover, we make use of the chain rule of differentiation. That is x u(x) = DS(ξ)− 3 ξ u(ξ). (9) Putting it altogether, the counterpart of the weak formulation (2b) in parametric space reads: find u ∈ V such that a(w, u) = L(w) for all w ∈ V with redefined (bi-)linear forms DS(ξ)− a(w, u) := DS(ξ)− w(ξ) u(ξ) |det DS(ξ)| dξ (10a) Ω0 w(ξ) |det DS(ξ)| DS(ξ)−1 DS(ξ)− = u(ξ) dξ. (10b) Ω0 G(ξ) w(ξ)f (S(ξ)) |det DS(ξ)| dξ. L(w) := (10c) Ω0 The 2 × 2 matrix function G(ξ) which involves the geometry map and its derivatives is termed the geometric factor [2] and it is reported to be symmetric uniformly positive definite on Ω0 [3]. An equivalent formulation based on the co-factor matrix Jc = DS(ξ)−1 det DS(ξ) reads [3] G(ξ) = − Jc S(ξ)Jc G11 = G21 det DS(ξ) G12 G22 (11) Remark: A similar transformation occurs in the parametric variant of classical finite elements when integrals over an individual physical element Tk (triangle/quadrilateral in 2D, tetrahedra/hexahedra in ˆ 3D) are replaced by integrals over the reference cell denoted by T . However, there are a few differences between the mapping in parametric finite elements and the geometry mapping in isogeometric analysis: 1. The mapping S(ξ) is defined globally, i.e. S : Ω → Ω0 rather than element-by-element as in finite ˆ elements, i.e. Fk : T → Tk , where Tk stands for an individual element of the triangulation. 2. The mapping S(ξ) is non-linear rather than the affine transformation Fk (ξ) = Aξ + b = x which is commonly used in classical standard finite elements. Remark: In so-called iso-parametric finite elements the element-by-element reference mapping Fk : ˆ T → Tk is defined in terms of the same finite element basis functions {ϕj }n that are used to represent j=1 the solution. Consider the quadratic P2 finite element depicted in Figure 3 with the degrees of freedom located at the vertices, edge midpoints and in the center. By defining the reference mapping as 7 x= ϕj (ξ)xj ˆ (12) j=1 ˆ letting ϕj denote the seven quadratic basis functions on the reference cell T , the so-defined iso-parametric ˆ finite element approach is able to represent quadratic boundaries exactly. y ˆ y T Fk 1 ˆ T 1 2 0 1 2 1 x ˆ 0 x Figure 3: Iso-parametric quadratic finite elements with quadratic reference mapping. 3 Quadrature-based assembly In what follows, we explain a simple idea for assembling the bilinear form (10b) efficiently by resorting to numerical cubature. This approach can be directly applied to the linear forms (10c). The main idea is to decompose the integral term over the entire parameter domain Ω0 into rectangular atomic contributions and apply standard cubature rules such as Gauss formulae on each rectangle afterwards. 4 Definition The index space is an equally-spaced representation of the tensor-product of the two knot vectors Ξ and H. It is spanned by the area [1, n + p + 1] in the i-direction and [1, m + q + 1] in the j-direction and is thus given by the rectangle [1, n + p + 1] × [1, m + q + 1]. The support of univariate B-spline basis functions is as follows: supp Ni,p (ξ) = [ξi , ξi+p+1 ), i = 1, 2, . . . , n, (13a) supp Mj,q (η) = [ηj , ηj+q+1 ), j = 1, 2, . . . , m. (13b) Thus, the support of the bivariate tensor-product B-spline functions reads supp ϕA (ξ) = supp Ni,p (ξ)Mj,q (η) = [ξi , ξi+p+1 ) × [ηj , ηj+q+1 ), A = 1, 2, . . . , n · m, (14) where A denotes a flat numbering of the degrees of freedom. Let us adopt row-wise/horizontal numbering. That is, given the one-dimensional data n, m, i, and j, the flat number of the (i, j)-th basis function is A = n(j − 1) + i. (15) On the other hand, the row and column indices can be recovered from the flat number A as follows: jA = A/n , iA = A − n(jA − 1). (16) It is therefore not necessary to store any mappings between the two-index numbering (i, j) and the flat numbering A since all required information can be computed efficiently on-the-fly. Example: Let p = q = 2 and iA = jA = 1 and iB = 3, jB = 2. The individual supports of the biquadratic B-spline basis functions ϕA=1 (ξ, η) = N1,2 (ξ)M1,2 (η) and ϕB=7 (ξ, η) = N3,2 (ξ)M2,2 (η) together with their virtual positions (’•’) in the index space are indicated in Figure 4. The green rectangle marks their common support in the index space. It should be noted that the common support can still be empty in the physical and parametric space if, e.g., ξ3 = ξ4 and/or η2 = η3 = η4 . η6 η5 η4 η3 η2 η1 ξ1 ξ2 ξ3 ξ4 ξ5 ξ6 ξ7 Figure 4: Index space with marked supports of ϕ1 (blue), ϕ7 (red) and ϕ1 ϕ7 (green). To develop an efficient strategy for assembling the stiffness matrix by numerical cubature consider the set of univariate B-spline functions depicted in Figure 5. The only non-zero basis functions on the interval [ξ2 , ξ3 ) of the index space are N1,2 , N2,2 , and N3,2 . This leads to Algorithm 1 for evaluating the values of all basis functions which are non-zero in a particular non-empty subinterval [ξk , ξk+1 ). 5 N1,2 ξ1 ξ2 N2,2 ξ3 N3,2 ξ4 ξ5 ξ6 ξ7 Figure 5: Quadratic B-spline fucntions with p = 2 and local support supp Ni,2 = [ξi , ξi+3 ). Algorithm 1 Evaluation of all B-splines of degree p that are non-zero in a given set of points. Let α = (α1 , α2 , . . . , αr ) be an array of points in [ξk , ξk+1 ). c c Initialize Bk,0 = 1 and Bj,0 = 0, j = k for c = 1, 2, . . . , r. for j = 1, . . . , p do for i = k − j, . . . , k do for c = 1, . . . , r do αc − ξi c ξi+j+1 − αc c c Bi,j = B B + ξi+j − ξi i,j−1 ξi+j+1 − ξi+1 i+1,j−1 end for end for end for c c return Bk−p,p , . . . , Bk,p for c = 1, 2, . . . , r. Let us consider the bilinear form (10b) for individual pairs (A, B) of basis functions in lieu of u and w and restrict it to a single cell [ξk , ξk+1 ] × [ηl , ηl+1 ] in the index space. That is ξk+1 ηl+1 ξk ξk+1 ηl ηl+1 S(ξ; k, l, p, q) := ϕA (ξ)G(ξ) ϕB (ξ) dξ (17a) ϕA,ξ G11 ϕB,ξ + ϕA,ξ G12 ϕB,η + ϕA,η G21 ϕB,ξ + ϕA,η G22 ϕB,η dξ (17b) = ξk ξk+1 ηl ηl+1 ξk ηl NiA (ξ)MjA (η) G11 (ξ, η)NiB (ξ)MjB (η) + (17c) G12 (ξ, η)NiB (ξ)MjB (η) + = (17d) 21 NiA (ξ)MjA (η) G (ξ, η)NiB (ξ)MjB (η) + 22 G (ξ, η)NiB (ξ)MjB (η) (17e) dξdη, (17f) where subscript ξ and η in (17b) denotes the first derivative with respect to the given variable. This term is obviously non-zero for all pairs (A, B) of basis functions which satisfy k − p ≤ iA , iB ≤ k and l − q ≤ jA , jB ≤ l. (18) Example: Let p = q = 2 and consider the indicated cell [ξk , ξk+1 ] × [ηl , ηl+1 ] in the index space. All possible locations for the basis functions ϕA and ϕB are marked by ’•’ in Figure 6. This allows us to device Algorithm 1 the efficient assembly of the global stiffness matrix corresponding to the discretised bilinear form (10b) by resorting to the local integral contribution (17c)–(17f). 6 Algorithm 2 Bilinear form assembly. Initialise the matrix S = {sAB } ≡ 0 for k = p + 1, . . . , n do for l = q + 1, . . . , m do 1) Choose one-dimensional Gauss cubature rules of order p and q and map the cubature points αc , c = 1, 2, . . . , r and βd , d = 1, 2, . . . , s both from the reference interval [−1, 1] to the parametric intervals [ξk , ξk+1 ] and [ηl , ηl+1 ], respectively, using linear mappings: ξk (1 − α) + ξk+1 (1 + α) 2 ηl (1 − β) + ηl+1 (1 + β) ˆ β= 2 α= ˆ 2) Compute the univariate basis functions Ni,p and Ni,p−1 in all cubature points α using the ˆ Algorithm 1 and store them locally as Ni,p,c and Ni,p−1,c , respectively. Apply the same ˆ procedure for Mi,q and Mi,q−1 using the cubature points β and store them as Mi,q,d and Mi,q−1,d . 3) Evaluate the geometric factor G(αc , βd ) for the tensor-product of the cubature points, i.e. [α1 , . . . , αr ] × [β1 , . . . , βs ] and store them locally as Gcd . 4) Loop over all pairs of basis functions (A, B), assemble the approximate integral (17c)–(17f) from precomputed data and store the result in the global matrix. Note that: – ωc and ωd denote the weights of the cubature rules of order p and q, respectively – the first derivatives of the basis functions are expressed by relations (7a)–(7b) – subscripts p and q in the basis functions are dropped to improve readability for iA = k − p, . . . , k do for jA = l − q, . . . , l do for iB = k − p, . . . , k do for jB = l − q, . . . , l do r s ωc ωd [NiA ,c MjA ,d (G11 NiB ,c MjB ,d + G12 NiB ,c MjB ,d )+ cd cd S= c=1 d=1 NiA ,c MjA ,d (G21 NiB ,c MjB ,d + G22 NiB ,c MjB ,d )] cd cd A = n(jA − 1) + iA B = n(jB − 1) + iB SAB := SAB + S end for end for end for end for end for end for 7 k, l + 1 k − 2, l k − 1, l k, l k − 2, l − 1 k−, l − 1 k + 1, l k, l − 1 k − 2, l − 2 k − 1, l − 2 k, l − 2 Figure 6: Possible positions of basis functions that are non-zero on the marked cell [ξk , ξk+1 ] × [ηl , ηl+1 ] for p = q = 2. 4 Quadrature-free assembly In [2, 3], a quadrature-free strategy for assembling the bilinear form (10b) efficient is proposed. In essence, the geometry factor (11) is projected onto a space that can be represented exactly by B-splines. That is, ΠG(ξ) = gk ϕk (ξ). (19) k This approximation is inserted into the bilinear form (10b). As a result, the stiffness matrix can be assembled efficiently using trivariate basis functions which are precomputed and tabulated. References [1] J. Austin Cottrell, Thomas J.R. Hughes, and Yuri Bazilevs. Isogeometric Analysis: Toward Integration of CAD and FEA. John Wiley & Sons, Ltd., 2009. [2] A. Mantzaflaris and B. J¨ttler. Exploring matrix generation strategies in isogeometric analysis. u Technical Report 1, G+S, April 2013. [3] A. Mantzaflaris and B. J¨ttler. Integration by interpolation and look-up for Galerkin-based isogeou metric analysis. Technical Report 14, G+S, June 2014. 8