Isogeometric Analysis and Shape Optimization of Holed Structures via the
Patch Removing Technique

In this study, a patch removing based Isogeometric analysis (PR-IGA) method is proposed to conduct the holed structural analysis with only one parametric domain, in which there are also no trimmed elements. The theoretical foundation of this novel patch removing approach is that any holed structure can be obtained by removing sub-patches (i.e., the holes) from an intact base patch. Since the parametric domains of these patches are all meshed by rectangular grids, the elements in the resulted holed structural parametric domain could all be untrimmed rectangles under certain mapping conditions. To achieve the special condition, a systematic technique consisting of T-spline local refinement and control points substitution/adjustment is provided. Due to the intactness of parametric elements, the analysis procedure of holed structures based on the proposed PRIGA is quite simplified and efficient compared to traditional multi-patch and trimming schemes. Moreover, after the deduction of analytical sensitivities related to structural mass and mechanical responses, the PR-IGA is directly employed in the holed structural shape optimization to successfully eliminate the need for model transformation during modeling, analysis and optimization processes. Numerical examples involving analysis and shape optimization of complex holed structures are presented to demonstrate the effectiveness of the proposed method.


Introduction
The Isogeometric analysis (IGA) approach originally proposed by Hughes et al. [1] is dedicated to integrate computer aided design (CAD) and finite element analysis (FEA) by invoking the isoparametric concept that the structural geometry and physical field share the same NURBS basis functions. This framework allows to model the structural geometry exactly at the coarsest level, to proceed with mesh refinement without interaction with the CAD system and to achieve high order smoothness in virtue of NURBS basis functions. Furthermore, the structural shape optimization based on this framework could be implemented straightforwardly since the coordinates and weights of NURBS control points are naturally the design variables [2][3][4][5][6].
IGA has attracted rapidly growing research interests for the above advantages [7][8][9][10]. However, due to the tensor-product form of NURBS basis functions, the complex holed structural analysis should be treated specially by the multi-patch IGA or by the trimming technique based IGA. With the former technique, the holed structure needs to be divided into NURBS patches artificially, and the corresponding modeling process is quite tedious and time consuming [11,12]; with the latter approach, there are inevitably many irregular trimmed elements in the holed structural parametric domain, of which the stiffness matrixes are calculated in a relatively complicated process [13,14].
Based on this observation, the present work is dedicated to performing the complex holed structural analysis by using a single patch with no trimmed elements. As shown in Figs. 1(a)-1(c), the complex holed structure Ω can be regarded as an intact base patch Ω base opened with a series of holes Ω hole_i (i = 1, 2, 3, 4). Since the base patch and hole-shaped sub-patches can all be represented by tensor-product NURBS surfaces, the holed structural parametric domain remove with only untrimmed elements could be obtained by removing regular parametric domain of hole-shaped sub-patches from the one of base patch, as plotted in Fig. 1(f). Figs. 1(d)-1(e) illustrate the main characteristics of IGA with multi-patch scheme and trimming technique, respectively. The former needs to partition the holed structure into many patches artificially; the latter could employ a single patch to analysis the holed structure, but with numerous irregular trimmed parametric elements.
Compared with the above two traditional approaches, the proposed patch removing based Isogeometric analysis (PR-IGA) technique can not only treat the problem of complex holed structure with only one patch, but also eliminate the trimmed elements in the parameter domain. With these characteristics, PR-IGA is very competitive for following reasons: The needlessness to deal with the problems caused by trimmed elements and irregular boundaries. Trimming technique is a powerful tool to describe the complex holed structure in CAD systems, but not an efficient way to realize the structural analysis by combining with IGA. The irregular trimmed elements not only bring several challenges to the stiffness calculation [13], but also make it difficult to impose essential and traction boundary conditions. The convenience to integrate with multi-patch IGA, and the resulted integration method could solve arbitrarily complicated engineering problems more flexibly. The engineering structure would sometimes be so complex that couldn't be modeled by one trimmed patch in CAD system or even be composed of different kinds of materials, so it is essential to use the multi-patch scheme to solve these problems. But the trimming technique can't work well with multiple patches for its irregular trimmed boundaries. The high efficiency to implement the shape optimization especially for complex holed structures. The control points of holed structural boundaries, which are ideal design variables for shape optimization, could be selected and adjusted directly in PR-IGA.
This paper is organized as follows. In Section 2, pertinent literature on NURBS, T-splines and the framework of isogeometric analysis is briefly reviewed. In Section 3, the PR-IGA strategies are explained, and the corresponding shape optimization schemes are discussed. In Section 4, several typical examples including numerical analysis and shape optimization are presented to verify the effectiveness and robustness of the proposed method. Finally, conclusions are drawn in Section 5.

B-splines and NURBS
B-spline basis functions can be defined by the Cox-de Boor iteration [15] and N i;p n ð Þ ¼ n À n i n iþp À n i N i;pÀ1 n ð Þ þ n iþpþ1 À n n iþpþ1 À n iþ1 N iþ1;pÀ1 n ð Þ for p ! 1 (2) in which ξ i is the ith knot, knot vector Ξ = {ξ 1 , ξ 2 , …, ξ n+p+1 } is a set of non-decreasing real numbers, p is the polynomial order and n is the number of B-spline basis functions. By using B-spline basis functions, a piecewise-polynomial B-spline curve is defined as where P i is the coordinate vector of the ith control point, and the control polygon is formed by {P i }. It should be noted that a standard knot vector Ξ in the CAD literature is an open knot vector, of which the first and last knot values appear p + 1 times. The B-spline curve with the standard knot vector can interpolate the control points at both ends.
The NURBS basis functions can be obtained by introducing weights and denominators into the B-spline basis functions where ω i is the weight of the ith control point. The NURBS curve that can exactly represent all ubiquitous shapes such as circles and ellipses is defined as The two dimensional NURBS basis functions can be constructed as R p;q i;j n; g ð Þ ¼ N i;p n ð ÞM j;q g ð Þx i;j P n i¼1 P m j¼1 N^i ;p n ð ÞM^j ;q g ð Þx^i ;ĵ (6) where N i,p (ξ) and M j,q (η) are univariate B-spline basis functions defined respectively on knot vectors Ξ = {ξ 1 , ξ 2 ,…, ξ n+p+1 } and η = {η 1 , η 2 ,…, η m+q+1 }, p and q are the polynomial orders, {ω i,j } are the weights corresponds to bidirectional control mesh {P i,j }. Similar to the NURBS curve representation (5), NURBS surface could be formulated as The NURBS surface is a tensor-product surface according to its definition (7), so that the introducing of a single control point will lead to the propagation of new control points along the entire row or/and column. To solve this problem and make the local refinement more convenient, some new spline forms such as T-splines [16][17][18] and THB-splines [19,20] are developed.

T-splines
T-splines introduce T-junctions into the NURBS framework and hence have the capacity of efficient local refinement [17]. Due to the existence of T-junctions, the T-spline blending functions, analogous to NURBS basis functions, cannot be directly obtained based on Ξ and η and should be calculated by using the T-mesh [21,22].
It can be seen from Fig. 2 that T-meshes are different with the NURBS parametric space because some inner intersection points are T-junctions. In order to infer the local knot vectors of blending functions in T-meshes, the anchors located at the centers of blending function supports need to be defined. Fig. 2(a) shows the cubic T-spline anchors P c1 and P c2, and Fig. 2(b) plots the quadratic T-spline anchors P q1 and P q2 . The local knot vectors of the above four anchors are given in Tab. 1 and marked by red dotted lines in Fig. 2. Since anchors, control points and blending functions of are one-to-one correspondence, the unidimensional blending functions N t i n ð Þ and M t i g ð Þ could be calculated according to the local knot vectors Ξ t i and η t i of the ith anchor. The two dimensional ith blending function is then given by where nt is the total number of the anchors in T-mesh, and ω i t is the weight of the ith control point. Supposing the Cartesian coordinate vector of ith control point is {P i t }, the T-spline surface expression can be finally written as

Isogeometric Analysis (IGA)
The main idea of IGA is to adopt the same basis functions (i.e., shape functions) to express the geometry in CAD and approximate the physical fields in finite element analysis. In this work, the shape functions of IGA are composed of NURBS basis functions defined by Eq. (6) and T-spline blending functions defined by Eq. (8). The IGA equilibrium equation can be expressed as in which K and F are respectively the global stiffness matrix and the global load vector, B is the strain matrix, D is elasticity matrix, R is the shape function matrix and d denotes the unknown displacement vector. According to the core concept of isoparametric element, K can be calculated on the parametric region spanned by knot vectors Ξ and η as where the Jacobi matrix J is defined as 3 Patch Removing Based Isogeometric Analysis (PR-IGA) 3.1 Basic Principle As shown in Fig. 3, a holed structure Ω can be obtained by removing the sub-patch Ω sub from the base patch Ω base , and both patches could be represented by NURBS surfaces with intact parametric regionsΩ sub andΩ base as follows where R S k;l n; g ð Þand R B i;j n; g ð Þare basis functions, Q k,l and P i,j denote control points, superscripts/subscripts S and B represent the sub-patch and the base patch, respectively.
The precondition of the equality reflected in Fig. 3 is Based on the definitions of S S and S B in Eq. (15), their basis functions should satisfy and control points related toΩ sub need also to be the same, as i;j n; g ð Þ 6 ¼ 0 and n; g ð Þ 2Ω sub (18) where both the differences i-k and j-l are constant values that depend on the location ofΩ sub inΩ base .

Implementation Procedure
Step 1: Construct the base patch and sub-patch In order to build NURBS patches for Ω sub and Ω base , both the outer and inner boundaries of the holed structure Ω are divided into four segments, as shown in Fig. 4(a). The 4 split points (see the black dots in Fig. 4(a)) for each boundary should be selected appropriately, so that the base patch and sub-patch expressed by Eq. (15) could be constructed reasonably with uniform grids. For determining the sub-patch parametric domainΩ sub , 4 end-to-end grid lines (see the red thick lines in Fig. 4(a)) are selected nearby the hole boundary. And the sub-patch knot vectors Ξ sub and η sub could be obtained according to the knots values ξ s1 , ξ t1 , η s2 and η t2 related to the above 4 grid lines, as plotted in Fig. 4(b). Since Ξ sub and η sub should be open knot vectors, their first and last knots need to be repeated 3 times for the construction of quadratic NURBS patch, as Ξ sub ¼ n s1 ; n s1 ; n s1 ; n s1þ1 ; . . . ; n t1À1 ; n t1 ; n t1 ; n t1 f g ; η sub ¼ g s2 ; g s2 ; g s2 ; g s2þ1 ; . . . ; g t2À1 ; g t2 ; g t2 ; g t2 È É : ξ s1 , ξ t1 , η s2 , η t2 in knot vectors Ξ base and η base of the parametric domainΩ base should have at least 2 multiplicities to satisfy Eq. (17). And the T-spline local refinement technique illustrated in Section 2.2 is adopted here to realize the knots repetition to avoid the propagation of new control points along the entire row or column.
Step 2: Substitute sub-patch control points for the corresponding ones of base patch According to Eq. (18), the control points of base patch should be partially replaced with sub-patch ones. For the sake of carrying out the replacement process more intuitively, the concept of anchors [21] is introduced in this work, which are defined in the parametric domain and share a one-to-one mapping relationship with the control points, as drawn in Fig. 5.
The sub-patch anchors are all marked by black dots in Fig. 5(a), and in Fig. 5(b) only the anchors belonging to the domain [ξ s1 , ξ t1 ] × [η s2 , η t2 ] are marked. Since the number and distribution of these two groups of marked anchors are the same, the corresponding control points substitution could be   implemented directly. Moreover, considering that ξ s1 , ξ t1 , η s2 , η t2 are repeated twice in base patch parametric domain, the NURBS model of the holed structure could be obtained by deleting the control points related to the anchors marked by red dots in Fig. 5(b) after the control points substitution process, because the values of basis functions related to these anchors are all zero at the hole boundary.
Step 3: Adjust the control points of the holed structure Fig. 6 illustrates the construction process of the holed structural control mesh, and in Fig. 6(c) we can see that the control mesh is seriously distorted near the hole boundary after the previous substitution step.
In order to solve the mesh distortion problem, a classical mesh deformation technology proposed in [23] is adopted here to adjust the control points. To be specific, the displacements of the hole boundary control points are firstly predefined according to their variation caused by the substitution operation. Next, a virtual load vector is supposed to produce the given displacements, and the initial material elastic modulus is defined as E old . Then, the equivalent strain of each element could be obtained through the isogeometric analysis, and for the ith element its elastic modulus is revised to (20) in which " e i denotes the ith element strain and " e allow is the upper bound of the equivalent strain. Finally, a new isogeometric analysis is performed based on the updated element elastic moduli, and the resulted displacements are directly taken as the coordinate adjusting quantity of the control points. Fig. 6(d) shows the final control mesh of the holed structure adjusted by this technique.

Holed Structural Shape Optimization Based on PR-IGA
A typical shape optimization problem is to minimize the structural mass m with a given constraint on the maximum stress σ max , and its mathematical statement is written as where ρ is the material density, σ lim is the allowable von-Mises stress, δ is a small positive number, and the 2nd inequation is used to ensure det J ≥ δ to avoid the appearance of distorted elements during the whole optimization process.
The design variables are the control point coordinates and weight factors of the design model (optimization model). The sensitivity of the objective function m to an arbitrary design variable α S could be expressed as The von-Mises stress can be calculated as in whichr ¼ r x ; r y ; r x À r y ; ffiffi ffi 6 p s xy À Á T , so the sensitivity of von-Mises stress is [24] @r von @a s ¼ 1 2r vonr T @r @a s (24) where could be obtained by the derivation of r h ¼ DBd ¼ r x ; r y ; s xy Â Ã T @r h @a s ¼ @D @a s Bd þ D @B @a s d þ DB @d @a s (25) in which D just depends on the material property, so that the first item is zero. Besides, if the optimization problem considered in this work is design-independent with @F=@a s ¼ 0, the following sensitivity calculation formula of d could be derived from Eq. (10) [25] @d Based on the above analysis, Eq. (25) could be simplified as @r=@a s where A is the solution of the accompany equation KA ¼ DB ð Þ T . According to Eq. (13), the sensitivity of global stiffness matrix @K=@a s corresponds to where @B=@a s can be obtained by solving the sensitivity of @R=@x @R=@y ð Þ T as @ @a s @R=@x @R=@y ¼ @ @a s @x=@n @y=@n @x=@g @y=@g ! À1 @R=@n @R=@g ! ¼ @J À1 @a s @R=@n @R=@g (29) in which @J À1 @a s ¼ ÀJ À1 @J =@a s ð ÞJ À1 . It can be seen from Eqs. (22), (28) and (29) that the sensitivity analysis problem is reduced to calculating @J =@a s and @ det J ð Þ=@a s , which can be further converted to the solving of @x=@a s and @y=@a s according to Eq. (14). Based on definition equations of NURBS and Tspline (Eqs. (6)-(9)), the sensitivities computation related to point coordinates (x, y) T in physical domain are finally transformed into the solution of @P k =@a s and @x k =@a s , in which P k and ω k are respectively the coordinate vector and weight factor of the kth control point in the analysis model.
In this work, only the control point coordinates of the design model are used as design variables, so that the velocity field of the kth control point just contains @P k =@a s . The velocity fields of all the control points can be classified into two categories: the ones related to boundary control points and the ones related to interior control points. The former could be calculated based on the grid refinement scheme [26], and the latter could be obtained according to the grid updating strategy [27,28].
In order to clearly illustrate the shape optimization procedure based on the proposed PR-IGA, the detailed flowchart is provided in Fig. 7.

Numerical Examples
In this section, three typical examples for mechanical analysis and shape optimization are investigated to demonstrate the effectiveness of the proposed method.

The Holed Structure Analysis
The holed structure with geometrical sizes and boundary conditions shown in Fig. 8 has been studied in "multiple patches" Section of the IGA monograph [11]. F n = 5 N/mm denotes the uniformly distributed load applied on the two upper semicircular boundaries. Young's modulus and Poisson's ratio are 71.018 GPa and 0.3, respectively.
In order to facilitate the base patch modeling, the exterior shape of this holed structure is revised by introducing a new patch bounded by the red closed curve as shown in Fig. 9, which should be removed later as the sub-patches. After being filled by this new patch, the base patch in this example turns into a simple square (100 mm × 100 mm in dimension). Fig. 10(a) shows the initial control mesh of the square base patch with only 9 control points. Process ① is mainly to construct sub-patches (including the new patch used to fill the structure in Fig. 9) according to the refined control mesh of the base patch. Process ② includes the substitution and adjustment operations of control points that are explained in Steps 2 and 3 in Section 3.2. The resulted control mesh of the holed structure is shown in Fig. 10(c).  Figure 9: Illustration diagram of the filling operation Fig. 11(a) reflects the unique advantage of the present PR-IGA method that there is no trimmed elements in the parametric domain. The physical mesh in Fig. 11(b) is generated exactly along the holed structural boundaries, so that the proposed method could keep all the merits of the IGA method.
Since the parametric domain of this holed structure is meshed by rectangular grids, the analysis process of PR-IGA is quite simple and efficient compared to the multi-patch and trimming schemes. Fig. 12(a) plots the displacement contours obtained by PR-IGA with 1428 control points, and the analysis result is highly consistent with it of ANSYS software in which the finest mesh with 5442 nodes is adopted.

The Multi-Material Structural Analysis
The proposed PR-IGA could also be applied to the analysis problem of multi-material structure [29], which is here considered to be made up by the holed structure and sub-patches (i.e., the holes) with different materials. Since no trimmed parametric elements are adopted in PR-IGA, there is only one material in each element of the resulted multi-material structure. Therefore, the corresponding solution procedure is quite simple and efficient.  According to the implementation procedure in Section 3.2, we can easily obtain the parametric mesh and control mesh of the multi-material structure, as shown in Fig. 14. And the only difference between the modeling processes of the holed structure and multi-material structure is whether the inner control points of sub-patches should be deleted. As the material of the sub-patch (i.e., the deep-colored circular region) differs from it of the holed structure (i.e., the light-colored region), we marked the sub-patch meshes in red color in Fig. 14. Fig. 15 plots the von-Mises stress contours obtained by the present method and the ANSYS software, and it could be observed that the maximum values (326.126 MPa and 327.413 MPa) of the von-Mises stress in these two contour plots are very close to each other.

The Holed Structural Optimization
The holed structure shown in Fig. 16 is a bracket for the aero-space application, and its shape optimization has already been considered as a standard benchmark in the structural optimization community [30][31][32][33]. Young's modulus and Poisson's ratio are 207.4 GPa and 0.3, respectively. The optimization problem expressed by Eq. (21) is investigated here with the upper bound of the von-Mises stress σ lim = 800 MPa. The method of moving asymptotes (MMA) [34] is chosen as the optimization algorithm, and stopping criterion is that the change rate of objective function (structural mass m) is less than 10 −6 .
Although the multi-patch based IGA has been successfully applied in the bracket shape optimization [35], there are still some shortcomings in their work, such as the difficulties of modeling the irregular patches and choosing the variable control points, the loss of smoothness and high-order continuity between adjacent patches, as shown in Fig. 17(a). It is worth mentioning that these shortcomings could be completely overcome by the following PR-IGA based shape optimization method. Fig. 18 draws the initial control mesh of the base patch, in which the control points are classified into 3 categories: variable control points, conjunct control points and linked control points. Only the coordinates of variable control points could be adopted as design variables. The movement rules related to latter two  categories of control points are decided by variable control points and structural geometric features (arcs and straight lines) that should be preserved. Fig. 19 shows the refined meshes of the bracket analysis model. It can be seen from Fig. 19(a) that there are no trimmed elements in the parametric domain. The physical mesh and the von-Mises stress contours of the optimized bracket are plotted in Fig. 20, and the maximum von-Mises stress of the design result is 799.3205 MPa that is less than σ lim . Also, it can be seen that the optimized geometry in Fig. 20 is in great agreement with the result in Fig. 17(b). Fig. 21 draws the convergence curve of the bracket mass, which is reduced from 0.341 kg to 0.113 kg (about a drop of 67%).

Conclusions
In this work, the PR-IGA method is proposed to realize the holed structural analysis and shape optimization by using a single patch with no trimmed elements. The basic principle, implementation procedure as well as application manners in shape optimization related to the proposed novel method are provided in this paper in detail. Numerical examples show that it is quite simplified and efficient to carry out mechanical analysis and shape optimization of complex holed structures by using the PR-IGA method.
It is worth noting that the PR-IGA method could also be applied to solving the multi-material structural analysis problem in the similar way to analyzing the holed structure. Besides, due to the adoption of highorder NURBS basis functions, the present IGA-based method have many merits in handling the beam and shell problem. Whereas, by using only the proposed PR-IGA, it is inefficient or even impossible to analyze the engineering structures with quite complex external profiles and topologies. And for this type of problem, it is necessary to combine the PR-IGA and multi-patch IGA to achieve the complementary advantages.