A High-Accuracy Single Patch Representation of Multi-Patch Geometries with
Applications to Isogeometric Analysis

This paper presents a novel approximating method to construct highprecision single-patch representation of B-spline surface from a multi-patch representation for isogeometric applications. In isogeometric analysis, multi-patch structure is not easy to achieve high continuity between neighboring patches which will reduce the advantage of isogeometric analysis in a sense. The proposed method can achieve high continuity at surface stitching region with low geometric error, and this technique exploits constructing the approximate surface with several control points are from original surfaces, which guarantees the local feature of the surface can be well-preserved with high precision. With the proposed approximating method, isogeometric analysis results using the new single-patch can be obtained efficiently compared with the original multi-patch structure. Several examples are presented to illustrate the effectiveness, accuracy and efficiency of the proposed method.


Introduction
Digital representation of object in real world is one of the basic problems in the field of Computer-Aided Design (CAD) and Computer Graphics (CG). Usually, we obtain the point clouds of real object through 3D laser scanning technology, then reconstruct the object by mesh, spline surface or implicit surface. As a standard representation in CAD industry, Non-Uniform Rational B-Splines (NURBS) or B-spline models are suitable to represent geometries because local and intuitive modification is allowed on these models. For objects with complex geometry, we usually use multi-patch representation to obtain approximate results with small errors. However, for the geometry with multi-patches, some geometric constraints should be added to achieve high-order continuity representation. Moreover, in the framework of isogeometric analysis [1] with multi-patches, more computation and operations will be involved to achieve the numerical simulation result over computational domain.
Interpolation and approximation of curves and surfaces through an array of points have been researched very early [2][3][4][5][6][7]. The algorithm of merging several B-spline curves to one B-spline curve has been proposed by Tai et al. [8]. They optimized the position of control points to merge two B-spline curves as precisely as possible, and adjusted the knots to avoid superfluous knots. Harish et al. [9] proposed an algorithm of merging two or more B-spline surfaces which is suitable for virtual reality applications. It's fast and robust, and even more robust than the common NURBS-based modeling software tool, Rhino. More B-spline curves and surfaces merging algorithms are referred in [10][11][12][13] and the references therein. The algorithm proposed in this paper focus on two or more B-spline surfaces, which combines least square method to increase the continuity at the stitching boundary and it is simpler than those existing methods. For surface approximation problem, least square method using B-splines is a simple approximation method, but it is sensitive to noisy data which may cause large errors. In order to reduce local error, Hoschek et al. [14] proposed an algorithm to construct a surface to fit an arbitrary point cloud and free boundary curves. It is found that in some cases they will not get the ideal fitting surface and might even fail, such as when the local density of sampling data points is high, or use high order polynomial to approximate, or the number of sampled data points is less than the number of control points of the B-spline. Weiss et al. [15] and Yang et al. [16] proposed weighted least square method according to importance of the sampling data points. They gave different weight factor to each different error term, and updated the weights by multiplying the resulted error on the weights in each iteration, thus the approximation is improved gradually. The selection and correction of the parameter during approximation is researched in [17][18][19]. Zhang et al. [20] optimized knot position to improve the quality of fitting surface, and their method can preserve geometric features. But the method takes scattered data points with parameterization as input. Shi et al. [21] proposed a scheme to construct G 1 smooth biquintic b-spline surfaces over arbitrary topology with interior single knots. This algorithm can deal with complex geometry for multi-patch representation. In the same year, a local scheme to construct convergent G 1 smooth bicubic b-spline surfaces is presented by Shi et al. [22]. The difference of this method with previous one is that the b-spline surfaces are not G 1 continuity between adjacent bicubic b-spline surfaces, but convergent to G 1 smoothness as number of control points increase.
After isogeometric analysis is introduced, many algorithms using multi-patch b-spline surfaces to represent geometry are proposed for isogeometric analysis [23][24][25][26][27][28][29][30][31][32][33]. Most of these methods consider the analysis-suitable multi-patch structure without high continuity between adjacent patches. Xu et al. [31] used global optimization to construct quadrilateral decomposition of the computational domain, and use local optimization to improve the quality of parameterization which result a G 1 continuous multi-patch parameterization. Kapl et al. [32] optimized a quadratic optimization problem with linear constraints to generate an analysis-suitable multi-patch parameterization.
In this paper, we propose a general approach to construct a high-accuracy single patch representation for an uncomplicated multi-patch spline structure. Information on the original multi-patch structure is partially kept in the final single-patch representation in order to keep as many geometric features as possible. In the proposed method, Newton iteration is applied to find the optimal parameters, and the final single-patch representation is obtained by using weighted least-square method with constraints. The results show that our method can improve the merging accuracy and the interior continuity of the geometry.
The rest of the paper is organized as follows. In Section 2, we briefly review the definition of B-spline surfaces. In Section 3, the procedure of our algorithm is described in details. Some experimental results and comparisons are presented to show the effectiveness of the proposed method in Section 4. We conclude this paper and describe the future work in Section 5.

B-Spline Surfaces
A tensor product B-spline surface is defined as follows: where d u ; d v represent degrees in u direction and v direction respectively, n u ; n v are numbers of basis functions in each direction, N d u i ; N d v j ði ¼ 0; 1; . . . n u À 1; j ¼ 0; 1; . . . n v À 1Þ are B-spline basis functions determined by knot vectors U ¼ ½u 0 ; u 1 ; . . . ; u n u þd u and V ¼ ½v 0 ; v 1 ; . . . ; v n v þd v ; the control points of the B-spline surface is denoted as p ij ði ¼ 0; 1; . . . ; n u À 1; j ¼ 0; 1; . . . ; n v À 1Þ. And the i-th B-spline basis function N i;k ðuÞ of degree k is defined as B-spline surface in Eq. (1) can be rewritten in the following form: 3 High-Accuracy Conversion from Multi-Patch Structure to Single Patch Representation In the proposed method, the input geometry is a multi-patch structure, boundaries of geometry can be accurately preserved or preserved with high precision in the final single patch representation. The goal of the proposed algorithm is to construct one B-spline surface to approximate the original surface model as precisely as possible. Given the boundary condition of the final single patch, number of control points of the final surface is determined. We will represent the boundary with at least C 1 continuity, and optimize the interior control points with some conditions to preserve the features of surface.

Objective Function
In order to construct the single patch representation, knot vectors of two directions have to be set which determines the number of control points, and boundary curves are chosen according to the original multipatch representation, which preserves the boundaries as much as possible. A B-spline surface is determined by knot vectors, degrees and control points. In this paper, the control points will be optimized by weighted least-square method, hence a number of scattered points which is more than the number of unknown control points are sampled from the original surfaces. The objective function related to control points in the optimization framework is where q t ; t ¼ 1; Á Á Á ; m are sampling points and w t ; t ¼ 1; Á Á Á ; m are weights. ðu t ; v t Þ; t ¼ 1; Á Á Á ; m are parameters corresponding to q t ; and selection of initial value of ðu t ; v t Þ will be described in the next subsection.
. . . ; d nÀ1 Þ T and Q ¼ ðq 1 ; q 2 ; . . . ; q m Þ T are the vectors of control point and the sampled data points. This objective function is actually the sum of distance from each sampled point to the approximate surface. By optimizing the objective function in Eq. (3), optimal control points of the final B-spline surface will be obtained.
The objective function is quadratic, so the optimization can be simplified to a system of linear equations. By computing the derivatives with respect to control points, we get the following linear equations: Then it can be simplified as follows by using matrix representation Since the initial parameters ðu t ; v t Þ are not the best choices, we will modify them in each iteration until the target precision and quality requirements are satisfied.

Parameter Modification
Since parametric domain of spline surface is a rectangle, we firstly divide the parametric domain into n Â n rectangular elements uniformly, which correspond to n Â n quadrilaterals in the spline surface.
The nearest corner point s i of quadrilaterals to the sample point q t is chosen as initial point of Newton iteration. Vector q t s i ! can be expressed as where ðu i ; v i Þ are parameters of s i on surface Sðu; vÞ;ñ is the normal vector of Sðu; vÞ at ðu i ; v i Þ, and d is the distance from q t to s i . In order to find s i on Sðu; vÞ which is closest to q t and the corresponding parameters ðu i ; v i Þ, we can get ðDu; DvÞfrom Eq. (5), and s i will be substituted by Sðu iþ1 ; v iþ1 Þ in the next iteration. The iteration is terminated when Du < e and Dv < e where e is a tolerance.
2. Convert the surface Sðu; vÞ into a quad mesh with quadrilaterals. 3. Find the quadrilateral that closest to q t and set a corner point as initial point s i . 4. Update the parameters of s i which is closest to q t using Eq. (6).

Error Computing
The procedure of single patch construction from the input multi-patch structure will stop when the merging error is small enough. The merging error is described by the distance between the sampling points q t and corresponding points on the approximating surface with respect to parameters obtained using Newton iteration. We define the maximal distance E max as follows, When E max < e, the single-patch representation is obtained.

Parameter Combination
A quadratic B-spline surface is determined by knot vector and control points. If there is a common boundary between two B-spline surfaces in one direction, suppose the V direction, then two surfaces can be merged exactly where knot vector of U direction is merged as follows, suppose U 1 ¼ ½0; 0; 0; u 1 ; u 2 ; . . . ; u n ; 1; 1; 1; From this procedure, the continuity along the common boundary is only C 0 on the merged surface. In particular, in order to increase the continuity of the final surface, we eliminate the duplicated knots In our algorithm, we will use this operation when we construct knot vectors for the final approximate surface.

Algorithm Framework
The flowchart of the algorithm is summarized as follows.
1. Obtain a set of sampling data points from the input multi-patch structure. 2. Determine boundary control points and interior control points which will be preserved in one-patch representation. 3. Parameterize the sampling data points using the method described in Section 3.2. 4. Compute the control points according to Eq. (4).

5.
Compute the maximal distance between the sampling data points and the approximating surface. Stop if the error is less than the given threshold; otherwise, return to Step 3.
In order to describe the method in details, we present the pseudo-code for the proposed framework.

Experimental Results
In this section, we will present some experimental results on single-patch approximation and two dimensional IGA applications to illustrate the effectiveness of the proposed method. Quadratic spline surfaces are used in every experiment. If there are more than two patches, they can be merged two by two with the proposed method. Therefore, only two surfaces being merged are illustrated. the maximal distance between the approximate surface and the data points avedist the average distance between the approximate surface and the data points pre_maxdist the maximal distance computed with previous approximate surface pre_avedist the average distance computed with previous approximate surface Fig. 2 shows a single-patch approximation using the proposed method, in which control points close to non-convex part of the input geometry boundary are preserved in the final single-patch representation. Better results are achieved compared with the method only using the geometric boundary information.

Single-Patch Approximation
In order to approximate the surface with given boundary information, we sampled 30 × 30 points on each surface. Fig. 2g shows the approximation results by weighted least-square method without any other constraints, in which some local self-intersections appear in the non-convex part. In our proposed method, some control points are preserved around the non-convex region, then we can obtain a much better result without self-intersection as shown in Fig. 2h, in which red control points are fixed points. Fig. 3 presents the color maps of distance error corresponding to the surface determined by control points shown in Figs. 2e and 2f after 20 iterations.
Tab. 1 lists the maximal and average distance between sampling points and the reconstructed single patch surface. For each iteration, we compare the merging results of approach with only boundary control points and the proposed method with boundary information and partial interior control points as constraints. The experimental results show that with few fixed interior control points, the merging error is reduced significantly. The next experimental example is shown in Fig. 4. Two surfaces that can be merged exactly, while the continuity on the common boundary is C 0 . If we use knot vector as shown in Eq. (8) to increase the continuity, the boundary control points around the common boundary will be approximated, not exactly preserved. The final iso-parametric curves of the merged surface are shown in Fig. 4h.
We also choose 30 × 30 points on each surface for weighted least-square method in this example. We can find that the nose on the face is a sharp feature, hence if we only preserve the boundary control points, there will be large merging error around this sharp region. In the proposed method, we fix some interior control points around the sharp region, and the merging error can be reduced immediately. Fig. 5 shows the color maps of the corresponding merging error after 20 iterations. The new proposed method improved the merging accuracy significantly and the continuity of merged surface along the common boundary is C 1 -continuity.
Tab. 2 lists the merging error for the 3D merging example in Fig. 5, and we can see that the maximal merging error by our method is reduced more than 80% than the maximal merging error of fitting surface by only fixing boundary curves. Furthermore, the average merging error is reduced more than 70%.
In the third example, we present a computational domain in which two patches can not be merged exactly because one boundary curve meets two boundary curves of another patch along the common boundary as shown in Fig. 6c. For this example, we use the knot vectors of the left patch as knot vectors of final merged surface. In order to use more information of the right patch, we insert some knots along the direction of common edge to increase the number of variables.
Tab. 3 shows that with the new proposed method, the merging error can be improved. Moreover, from the iso-parametric curves shown in Figs. 6g and 6h, we can see that the iso-parametric curves are more uniform by our merging method compared with the method only using boundary control points.

Application in Isogeometric Analysis
In this section, we apply the proposed merging method in isogeometric analysis. The input multi-patch structure is merged into a single B-spline patch and the specified PDE problem is solved on the single patch. Compared with the IGA method on the multi-patch structure, higher accuracy can be achieved with the proposed single patch representation. Consider the following Poisson equation the weak form of Eq. (8) is given as wðxÞ 2 H 1 ðÞ; w @ j ¼ g; such that for all v 2 H 1 0 ðÞ, aðw; vÞ ¼ f ðvÞ; where aðw; vÞ ¼ R rwrvd; f ðvÞ ¼    where P i denotes the control points and N i ðs; tÞ are B-spline basis functions. The finite dimensional function space V h can then be defined as where R i ðx; yÞ ¼ N i G À1 ðx; yÞ and n is number of basis functions which are zero on the boundary of domain .
The weak form of the isogeometric approximation becomes: find u h 2 V h ; such that for all v h 2 V h ; Due to finite dimension of V h , the approximate solution can be written as u h ¼ P n i¼1 u i R i ðx; yÞ. And Eq. (11) becomes a linear equation system, If the computational domain consists of a set of connected sub-domains 1 ; 2 ; Á Á Á ; s and each sub-domain is represented as    where R ij ðx; yÞ ¼ N ij G À1 i ðx; yÞ;m i is number of basis in i . On the common boundary C ii 0 of two connected sub-domain i and i 0 , u h i and u h i 0 are set to be identical to ensure that solution u h is continuous over . And C ii 0 has the same representation for both sub-domains i and i 0 .
From the above analysis, numerical solution space V h of Poisson equation defined on can be chosen as V h ¼ fR ij ; j 2 I i ; i ¼ 1; . . . ; s; R ij þ R i 0 j ; j 2 K ii 0 ; for each pair of connected sub-domains i and i 0 g Here, I i is the index set such that R ij ðx; yÞ is a basis function which is non-zero in the interior of i for each j 2 I i . K ii 0 is the index set such that for each j 2 K ii 0 , R ij ðx; yÞ and R i 0 j ðx; yÞ are basis functions of interior boundary C ii 0 in i and i 0 respectively. We implement the above framework on a two-dimensional surface shown in Fig. 2, where u has an exact solution x 2 þ y 2 on the given domain. Fig. 9 shows the numerical error, the L 2 error is 3.09214 × 10 −4 with DOF (Degree of Freedom) is 180 for the original multi-patch representation while the L 2 error is 2.51835 × 10 −4 with DOF is 168 for the new single-patch representation.  In summary, by the proposed patch-merging method, we can get more accurate numerical results with the new single-patch compared with the original multi-patch structure for incomplicated domains.

Conclusion
In this paper, we propose a new surface-merging method to construct a single patch representation from a multi-patch structure for isogeometric analysis. Besides the sampled points of geometry and the given boundary information, part of interior B-spline information is also preserved to construct the final surface. Compared with the method only using sampled points of geometry and boundary information, our method can achieve high-accuracy merging results. The resulting single patch is suitable for isogeometric analysis on the input 2D computational domain. On the other hand, the merged representation is C 1 continuous in the interior part, which improved the merging result directly. But if the multi-patch structure is very complex, such as many features need to be preserved, our algorithm may not produce a one-patch representation as good as we expect. This is a major limitation that we need to improve.  In the future, we will investigate on automatic selection of the fixed control points. The generalization to trivariate volume-merging problem is also a part of future work [26,28,29].