A new method for T-spline parameterization of complex 2D geometries

Engineering with Computers manuscript No. (will be inserted by the editor) A new method for T-spline parameterization of complex 2D geometries M. Brovka · J.I. L´ opez · J.M. Escobar · J.M. Casc´ on · R. Montenegro the date of receipt and acceptance should be inserted later Abstract We present a new strategy, based on the idea of the meccano method and a novel T-mesh optimization procedure, to construct a T-spline parameterization of 2D geometries for the application of isogeometric analysis. The proposed method only demands a boundary representation of the geometry as input data. The algorithm obtains, as a result, high quality parametric transformation between 2D objects and the parametric domain, the unit square. First, we define a parametric mapping between the input boundary of the object and the boundary of the parametric domain. Then, we build a T-mesh adapted to the geometric singularities of the domain in order to preserve the features of the object boundary with a desired tolerance. The key of the method lies in defining an isomorphic transformation between the parametric and physical T-mesh finding the optimal position of the interior nodes by applying a new T-mesh untangling and smoothing procedure. Bivariate T-spline representation is calculated by imposing the interpolation conditions on points sited both on the interior and on the boundary of the geometry. The efficiency of the proposed technique is shown in several examples. Also we present some results of the application of isogeometric analysis in a geometry parameterized with this technique. Keywords T-spline parameterization ·Simultaneous T-mesh untangling and smoothing ·Meccano method · Isogeometric analysis ·Quadtree M. Brovka ·J.I. L´ opez ·J.M. Escobar ·R. Montenegro University Institute for Intelligent Systems and Numerical Applications in Engineering (SIANI), University of Las Palmas de Gran Canaria, Spain E-mail: {mbrovka, jilopez,jmescobar,rmontenegro}@siani.es J.M. Casc´ on Department of Economics and History of Economics, Faculty of Sciences, University of Salamanca, Spain. E-mail: [email protected] 1 Introduction An open problem in the context of isogeometric analysis [1, 2, 7] is how to obtain a spline parameterization of a complex computational domain from the CAD description of its boundary. Parameterization is suitable for analysis if it does not have self-intersection, i.e., the Jacobian is strictly positive at any point of the computational domain. Moreover, in order to expect a high accuracy in numerical results it is necessary to obtain a good quality parameterization. Orthogonality and uniformity of isoparametric curves are desirable for the tensor-product structured parameterization. This task is not trivial and can be very time-consuming. For application of IGA it is essential to have an efficient method to construct T-spline parameterization. In the present work we investigate this problem for planar geometries. There are only a few works addressing this problem. In [32], the parameterization is found by solving a constraint optimization problem for the control points of a planar Bspline surface. Constrains are defined by imposing injectivity sufficient conditions in terms of control points, and the optimization consists in the minimization of some energy functions in order to reach a good orthogonality and uniformity of the parametric mapping. Another similar technique was proposed by these authors in [33]. They use a harmonic mapping obtained by solving an optimization problem for the control points. Additional term is added to the objective function in order to improve the quality where needed. The use of harmonic mapping is a common characteristic of several works dealing with 2D and 3D parameterization methods [20, 21, 22]. In this paper, we propose a different approach where the parameterization is accomplished by transforming isomorphically a T-mesh from the parametric domain to the physical one. The construction of this transformation is mainly based on a simultaneous T-mesh untangling and smoothing procedure. As far as we know, the only case of performing mesh untangling and smoothing procedure for a T-mesh, in order to construct T-spline representation of 3D domains, was described in [31, 34]. They remove tangled elements by maximizing the worst Jacobian. Smoothing is performed by moving each node towards the mass center of its neighboring elements. In our previous works [12, 14] we constructed the physical T-mesh of the solid using a volumetric parameterization obtained by deforming a tetrahedral mesh of the solid. In general, this approach does not provide an optimal Tmesh quality in the sense of its uniformity and orthogonality. Now, we propose a different approach where the optimization is applied directly to the T-mesh. We use an interpolation scheme to fit a T-spline object to the data instead of an 2 Fig. 1 General scheme of the T-spline parameterization method. approximation, as other authors do [32, 33]. This performs a more accurate adaptation of the T-spline to the input data. Our technique is simple and easy to implement. Satisfactory results are obtained with low computational effort for a variety of complex geometries. This paper is organized as follows. In next section we describe the main steps of the proposed algorithm. In section 3 we explain the process of boundary parameterization and the construction of T-mesh adapted to the singularities of the object boundary. Section 4 describes the simultaneous T-mesh untangling and smoothing procedure that leads to the construction of a high quality T-mesh of the object. The modeling of the geometry by means of bivariate T-splines is developed in section 5. A quality improvement strategy, based on adaptive refinement, is described in the section 6. In section 7 we illustrate the results of the application of the algorithm to some 2D domains and solve the Poisson equation in one of them using isogeometric analysis with Tsplines. Finally, in section 8 we present the conclusions and set out some challenges. 2 General scheme of the method In this section we summarize the proposed method in order to facilitate its understanding. Some ideas of the method are taken from our previous works on mesh untangling and smoothing and the meccano method [4, 5, 9, 11, 23, 24], but they have been adapted to the requirements of the present work. The algorithm includes the following stages: 1. Boundary parameterization and construction of an adapted T-mesh: A bijective correspondence between the input boundary of the object and the boundary of the parametric domain is defined. Then, an adapted T-mesh is generated by refining the initial mesh in order to approximate the geometry with a prescribed tolerance. During this process, the boundary nodes of the parametric domain are mapped to the boundary of the object. 2. T-Mesh optimization: We relocate the inner nodes of the T-mesh by applying a simultaneous mesh untangling and smoothing procedure. A previous relocation of the inner nodes is accomplished in order to facilitate this task. 3. Construction of a T-spline representation of the geometry: The T-spline parameterization is obtained by imposing interpolation conditions. As interpolation points, we take the vertices of the physical T-mesh obtained after the optimization process and other additional points. 4. Adaptive refinement to improve the mesh quality: If the quality of the mesh is not satisfactory, we apply an adaptive refinement in order to increase the degree of freedom in the areas with high distortion. Then, we return to step 2 and repeat the process until reaching a good T-spline parameterization. The T-mesh used in this work has a balanced quadtree structure [26]: all refinements are performed by dividing a cell into 4 equal cells and a procedure of balancing is applied in order to allow just one hanging node per edge. The input boundary representation is given by a polygonal, however it could be given by B-spline curves. A scheme of the algorithm is shown in Fig. 1. In the next sections we describe the details of each stage. 3 Boundary parameterization and construction of an adapted T-mesh In order to define a parametric mapping between the input boundary polygonal of the object and the boundary of the parametric domain, the unit square, we have to select four points of the polygonal that will correspond to the four corners of the square. These points divide the input polygonal into four parts that are mapped via chord-length parameterization into its corresponding edge of the square. The selection of the corners must be appropriate. Namely, the corners should be situated in convex areas of the input boundary polygonal, where the inner angle formed by the polygonal in these points is less than 180◦. Obviously, the optimal value of this angle is 90◦. Next, we construct an adapted T-mesh that approximates the input boundary with a pre-defined tolerance ε. To do that, an approximation error is calculated for each boundary cell and the cell is refined if this error is greater than 3 (a) (b) Fig. 2 (a) Illustration of the refinement criterion used for the construction of an adapted T-mesh in order to approximate the object boundary with a desired tolerance. (b) Example of the boundary approximation corresponding to the upper edge of the unit square. Red line represents the input boundary. ε. A cell refinement produces a new boundary point that is projected over the input boundary polygonal, obtaining a more accurate approximation of the geometry, as illustrated in Fig. 2(b). As approximation error, we evaluate the area of the triangles formed by the cell edge in the physical space and each vertex of the input polygonal approximated by this edge. A cell is refined if there is at least one triangle whose area exceeds ε. The refinement criterion is illustrated in Fig. 2(a). Another approximation error criterion can be used, such as the Hausdorff distance between the input boundary polygonal and the cell edge. Optionally, it is possible to start the boundary refinement and projection procedure with some initial refinement of the parametric domain in order to guarantee a sufficient number of free nodes in the interior of the geometry. This may facilitate the untangling procedure in complex geometries. As result of this stage, the position of the boundary nodes in the physical domain are known and the position of the inner nodes will be defined by means of the T-mesh optimization procedure developed in next section. Fig. 3(a) shows an example of the adapted parametric T-mesh constructed in this stage and, Fig. 3(b), illustrates the resulting tangled T-mesh of the physical domain that will be optimized. 4 T-mesh optimization The key of the proposed method lies in the optimization procedure that allows to obtain a high quality physical T-mesh, that is used to construct the T-spline representation of the object. 4.1 Previous relocation It is preferable to perform a previous relocation of the inner nodes in order to reduce the computational effort during the optimization process. In the present work, we have used for this purpose Coons patch [6, 15] to define a surface that interpolates given boundary curves. Let the given boundary curves be called x(ξ,0),x(ξ,1),x(0,η),x(1,η). The Coons function is defined as x(ξ,η) =(1−ξ)x(0,η)+ξx(1,η) +(1−η)x(ξ,0)+ηx(ξ,1) −1−ξ ξ x(0,0)x(0,1) x(1,0)x(1,1)1−η η This previous relocation procedure facilitates the untangling process, but in general does not obtain a satisfactory mesh quality and can produce self-intersections, as shown in Fig. 3(c). Therefore, it is essential to apply an efficient optimization algorithm. 4.2 Objective function The mesh optimization process is carried out by iterative relocation of each inner node of the mesh in such a way that 4 (a) (b) (c) (d) Fig. 3 Stages of T-mesh construction for the Spot test model geometry. (a) Parametric T-mesh adapted to the boundary of the geometry; (b) tangled physical mesh after boundary mapping (the color bold lines represent the correspondence between parametric and physical boundaries); (c) previous relocation by using Coons patch; (d) optimized physical T-mesh. the new position of the node improves the quality of the local submesh corresponding to this node. A local submesh is the set of all the elements connected with the movable or free node. The local objective function for a free node is based on algebraic shape quality metrics proposed by Knupp in [18, 19] for triangular and quadrilateral elements. Shape quality metric for a given triangle is defined in terms of the Jacobian matrix of the affine mapping from ideal triangle to the given one. This shape quality metric represents the deviation of the physical triangle from the ideal one. It attains its maximum value, 1, if the triangle is similar to the ideal one, and it equals 0 if the triangle is degenerated. The distortion metric of an element is defined as the inverse of its quality metric. In order to asses the quality of the local submesh for a given free node of a T-mesh, we have to decompose each neighboring cell into triangles and asses the quality of each triangle. For a T-mesh, this decomposition depends on the type of the free node. There are two types of free node: a regular node and a hanging node. The optimal position of each free node is determined by minimizing a local objective function. We define the objective function as a sum of shape distortion metrics of the triangles of the local submesh. For each triangle of the physical mesh, the corresponding triangle of the parametric mesh is used as its ideal element. Therefore, each cell of the physical mesh tends to have the same shape as its counterpart cell of the parametric mesh. Thus, repeating this procedure for all the inner nodes of the mesh, we achieve the physical mesh of the object be as similar to the parametric one as possible. A regular node is surrounded by four cells with equal or different sizes. In order to perform the mesh improvement, the local submesh is decomposed in twelve triangles, three triangles per cell whose qualities depends on the position 5 (a) (b) (c) (d) Fig. 4 Triangular decomposition of the free node local submesh. (a) Regular node case, where each cell is decomposed in three triangles; (b) hanging node case, where five triangles are formed in the cell where the node generates a T-juntion; (c) barriers and feasible region induced by the 12 triangles in the objective function for a regular node; (d) barriers and feasible region induced by the 11 triangles in the objective function for a hanging node. of the free node. Fig. 4(a) illustrates this decomposition and Fig. 4(c) shows the feasible region of the objective function. In a hanging node case, the free node is surrounded by three cells and the local submesh is decomposed in eleven triangles. The cell in which the node forms a T-junction is decomposed in five triangles whose qualities depend on the position of the free node. Each one of the other two cells is decomposed in three triangles, as was described in the case of a regular node. Fig. 4(b) shows the decomposition of a local submesh for hanging node case and the feasible region of the objective function. Note that, for the ideal case shown in Fig. 4(d), the feasible region induced by these eleven triangles is the same as the one obtained after a refinement of the T-junction cell, see Fig. 4(c). In order to define the objective function we introduce the following concepts. Let Tbe a triangle whose vertices are given by xk= (xk,yk)T∈R2,k=0,1,2 and TRbe the reference triangle with vertices u0= (0,0)T,u1= (1,0)Tand u2= (0,1)T. If we choose x0as the translation vector, the affine map that takes TRto Tis x=Au+x0, where Ais the Jacobian matrix of the affine map referenced to node x0, and expressed as A= (x1−x0,x2−x0). Let consider that TIis our ideal or target triangle whose vertices are v0,v1and v2. If we take v0= (0,0)T, the linear map that takes TRto TIis v=Wu, where W= (v1,v2)is its Jacobian matrix. As the parametric and real meshes are topologically identical, each triangle in the physical space has its counterpart in the parametric space. Affine map that takes TIto Tis given by x=AW−1v+x0, and its Jacobian matrix is S=AW−1. Note that this weighted matrix Sdepends on the node chosen as reference, so this node must be the same for Tand TI. Quality metrics of the triangle Tcab be defined in terms of the matrix S. For example, the mean ratio,q=2σ kSk2, is an easily computable algebraic quality metric of T, where σ=det(S)and kSkis the Frobenius norm of S. The maximum value of qis the unity, and it is reached when A=µRW , where µis a scalar and Ris a rotation matrix. In other words, qis maximum if and only if Tand TIare similar. Besides, any flat triangle has quality measure zero. We can derive an optimization function from this quality metric. Thus, let x= (x,y)Tbe the position of the free node, and let Smbe the weighted Jacobian matrix of the m-th triangle connected to this free node. We define the objective function of x, associated to an m-th triangle as ηm=kSmk2 2σm (1) The local objective function used for mesh quality improvement is defined by means of the inverse of mean ratio quality metric of each triangle of the local submesh. The function to be minimized is given by K(x) = M ∑ m=1 kSmk2 2σm (2) where Mis the number of triangles in the local submesh and Smis the Jacobian matrix associated to the affine mapping from the ideal triangle to the physical one. Objective function defined by Eq. 2 is appropriate to improve the quality of a valid mesh, but it does not work properly when there are inverted elements [16, 17]. In previous 6 works [9, 13] we have used a modified objective function K∗, where the untangling and smoothing are carried out in the same stage. This modified objective function K∗does not have singularities, it works as the original function K for the valid elements and tends to untangle the inverted and degenerated ones. This objective function is defined as K∗(x) = M ∑ m=1 kSmk2 2h(σm)(3) where h(σ) = 1 2(σ+√σ2+4δ2). Objective function K∗is smooth in R2, so the unconstrained optimization problem can be easily solved with any standard method (see for example [8]). In case of a conformal local submesh, the result obtained by minimizing the objective function K∗is, when possible, an orthogonal submesh, as shown in Fig. 5(a). However, not satisfactory result is obtained for a non-conformal submesh. In this case, two special situations can appear: a regular node surrounded by cells of different scales and a hanging node. In these situations, a variation in the position of the free node does not affect in the same way to the quality of the triangles of the local submesh. The objective function tends to form triangles as similar as possible to the reference ones, but the influence of the smaller cells are greater than the bigger ones. For example, in Fig. 5(b) it can be seen how the free node is moved toward the small cell and, therefore, the resulting mesh is not orthogonal. This problem is solved by a modification of the objective function K∗, namely multiplying the terms of the objective function by appropriate weights. 4.3 Weighted objective function The terms of the objective function K∗can be grouped according to the belonging to each cell of the local submesh. Each group is multiplied by an appropriated weight in order to avoid the problems mentioned in section 4.2. For a regular node, the weighted objective function is K∗ τ(x) =τ1 3 ∑ m=1 kSmk2 2h(σm)+τ2 6 ∑ m=4 kSmk2 2h(σm)+ +τ3 9 ∑ m=7 kSmk2 2h(σm)+τ4 12 ∑ m=10 kSmk2 2h(σm) (4) where each summation is the group associated to each cell and τiis the applied weight. This weight is equal to the scale factor of the cell in the parametric space. We assume that the smallest cells in the local submesh have scale factor τ=1 and the other cells can have scale factor τ=2 or τ=4, as illustrated in Fig. 6(a). Fig. 5(c) shows the resulting orthogonal mesh when these weights are applied. (a) (b) (c) Fig. 5 Resulting meshes after optimization with different objective functions. (a) Orthogonal mesh using K∗; (b) not satisfactory result for a non-conformal local submesh using K∗; (c) orthogonal mesh using weighted objective function K∗ τ. (a) (b) Fig. 6 Weights applied to the objective functions K∗ τ. (a) Regular node where τ1=1, τ2=τ4=2 and τ3=4; (b) hanging node where τ1= τ2=1 and τ3=8 5. A hanging node is surrounded by three cells as it was mentioned above. In this case, the weighted objective functions is K∗ τ(x) = τ1 3 ∑ m=1 kSmk2 2h(σm)+τ2 6 ∑ m=4 kSmk2 2h(σm)+τ3 11 ∑ m=7 kSmk2 2h(σm) (5) Hanging node is a more particular case because its local submesh is decomposed in different types of triangles. To guarantee the orthogonality in the local submesh after optimization, we have determined that the weights are τ3=8 5for the cell where the node forms a T-junction and τ1=τ2=1 for the other two cells, as shown in Fig. 6(b). The election of these weights is justified in appendix. 7 Fig. 7 One-to-one global mapping between parametric domain ˆ Ωand physical domain Ω. 5 Construction of a T-spline representation of the geometry In this section we describe the construction of T-spline representation for 2D geometries. We have to obtain a global one-to-one parametric transformation that maps the parametric domain into the physical domain S:ˆ Ω= [0,1]2→Ω, see Fig. 7. We start with a brief summary of the main concepts about B-splines and T-splines. 5.1 B-spline basis functions A set of B-spline basis functions Ni,p(i=1,2,...,n)of degree p, inferred from a knot vector Ξ={ξ1,ξ2,...,ξn+p+1} is defined by the Cox-de Boor recursion formula: Ni,0(ξ) = 1 if ξi≤ξ<ξi+1, 0 othewise. Ni,p(ξ) = ξ−ξi ξi+p−ξi Ni,p−1(ξ) + ξi+p+1−ξ ξi+p+1−ξi+1 Ni+1,p−1(ξ). A non-decreasing sequence Ξis called open knot vector if the first and the last knot are repeated p+1 times. At each knot of multiplicity kthe basis functions are Cp−k.In this paper we use B-spline functions of degree p=3. Bivariate B-splines are defined as a tensor product of two univariate B-spline functions: Ni,p(ξ) = Ni1,p(ξ)Ni2,p(η), where ξ= (ξ,η)and multi-index i= (i1,i2)∈I. The multiindex set is defined by I={1,2,...,n1}×{1,2,...,n2}. B-spline surface is defined as lineal combination of bivariate B-spline functions S(ξ) = ∑ i∈I PiNi,p(ξ), where Pi∈R3are the control points. In our case, Pi∈R2, so S(ξ)describes a planar domain. For more details about B-splines and NURBS see [25]. 5.2 T-splines In order to define T-spline basis functions of degree 3 over a given T-mesh in 2D, a local knot vector for both parametric directions should be assigned to each basis function Rα: Ξα=ξα1,ξα2,ξα3,ξα4,ξα5,Hα=ηα1,ηα2,ηα3,ηα4,ηα5. These knot vectors are inferred by traversing T-mesh edges. This basis function is associated to the central knot (ξα3,ηα3) that is called anchor. As we are using open knot vector structure, there are some blending functions that have the same anchor. The T-spline blending functions that we use in this work are rational B-spline functions defined as Rα(ξ) = Nα(ξ) ∑ β∈A Nβ(ξ)(6) being Nα(ξ) = N1 α(ξ)N2 α(η)the bivariate B-spline function defined over its local knot vectors Ξα={Ξα,Hα}, and Ais the index set of the basis spanned by T-mesh. A detailed report about T-splines and their relationship with isogeometric analysis can be found in [2]. 5.3 Interpolation We build bivariate T-spline surface representation of our physical domain as lineal combination of T-spline blending functions S(ξ) = ∑ α∈A PαRα(ξ)(7) where Pα∈R2is the control point corresponding to the α-th blending function. Control points Pαare found by imposing interpolation conditions. Assuming that the set of blending functions are linearly independent, we need as many interpolation points as blending functions. As interpolation points, first of all, we use the anchors of the blending functions. Each anchor coincides with a Tmesh vertex, ξv α, and its position in the physical space, xv α, was determined by the mesh optimization process. According to the strategy adopted in the section 5.2, that leads to an open knot vector structure along the boundary of the mesh, the total number of blending functions is greater than the number of vertices of the T-mesh. Therefore, we have to take additional interpolation points. These additional points are associated to functions whose knot vectors Ξαor Hαcontains exactly three repeated knots. We assign to each of these functions an interpolation point approximately sited where the function attains its maximum. The positions of these interpolation points in parametric and physical spaces are shown in Fig. 8(a). For each additional interpolation point ξα, its position in the physical space, xα, should be defined. 8 (a) (b) Fig. 8 Interpolation points in parametric and physical domains. (a) Parametric domain, ξvvertices of the mesh, ξeadditional interpolation points situated on the edges of the mesh, ξcadditional points situated in the four corner cells; (b) interpolation points in physical domain. The interpolation points ξe α, situated at the mid point of an edge of the parametric mesh, are mapped to the mid point of the same edge of the physical T-mesh. The positions in the physical space of the four interpolation points ξc αsituated in the center of the corner cells of the parametric mesh are determined by optimizing the local submesh formed after a fictitious quadtree partition of these cells. Finally we solve the linear system of equations xβ=Sξβ=∑ α∈A PαRαξβ,∀ξβ,β∈A(8) where ξβare interpolation points in parametric space and xβ are their images in the physical space. 6 Adaptive refinement to improve the mesh quality 6.1 Mean ratio Jacobian Our objective is to get high-quality geometry parameterization suitable for isogeometric analysis. The parametric Tspline mapping of Eq. 7 is suitable for analysis if it has positive Jacobian in all the domain. High distortion of the geometry can produce a large variation of the Jacobian that can lead to a poor accuracy in the numerical results. Therefore, a good uniformity and orthogonality of the isoparametric curves are desired for the parametric mapping S. A high quality of the optimized T-mesh is a necessary, but not sufficient, condition for a high quality of the T-spline mapping. It can happen that the Jacobian of the spline parameterization takes negative values even if all the cells of the T-mesh are valid. In order to assess the quality of the constructed parametric transformation we analyze the mean ratio Jacobian, given by Jr(ξ) = 2 det(J) kJk2,(9) where Jis the Jacobian matrix of the mapping Sat the point ξ= (ξ,η)and kJkis its Frobenius norm. The value of the mean ratio Jacobian at any point P0of the parametric domain is a shape quality metric for the infinitesimal triangle formed by two isoparametric curves of the physical domain passing through the point P0 0=S(P0), as illustrated in Fig. 9(a). In contrast to the scaled Jacobian, that represents a quality of the mapping Sin the sense of the orthogonality of its isoparametric curves, the mean ratio Jacobian represents both: a quality of the mapping in the sense of the orthogonality and uniformity of its isoparametric curves. Fig. 9(b) shows the comparison between the scaled Jacobian and the mean ratio Jacobian. Scaled Jacobian attains its maximum value 1 at the given point if the mapping conserves orthogonality of the isoparametric curves. Mean ratio Jacobian is equal 1 at the point P0if the mapping conserves orthogonality and produces the same length distortion in both parametric directions, i.e., the mapping is conformal at this point. It is easy to see that ∀ξ: 0 ≤ |Jr(ξ)| ≤ |Js(ξ)| ≤ 1, where Js=det(J) kSξkkSηkis scaled Jacobian. 6.2 Adaptive refinement Parameterization of complex geometries entails a severe distortion that can lead to appearance of low quality cells, even cells with negative Jacobian. This can be explained by the lack of degrees of freedom provided by the inner nodes. In order to improve the mesh quality in this case, we propose an adaptive strategy that refines all the cells with low quality. A similar idea was implemented for tetrahedral meshes in [10]. We proceed as follows. For each cell of the mesh, the mean ratio Jacobian is calculated at Gauss quadrature points. We use 16 =4×4 quadrature points per cell. A cell ˆ Ωeis marked to refine if, at least, one of its quadrature points has 9 (a) (b) Fig. 9 Mean ratio Jacobian. (a) A quality metric of the parametric mapping Sat any point P0in terms of the mean ratio of the triangle P0 0P0 1P0 2. (b) Comparison between the mean ratio Jacobian and the scaled Jacobian (a) (b) (c) (d) (e) (f) (g) Fig. 10 Adaptive refinement strategy to improve the parametric transformation quality in Gran Canaria island domain. (a) T-spline representation of the domain; (b) initial physical T-mesh; (c) final physical T-mesh; (d) initial T-spline parametrization with negative Jacobian; (e) resulting T-spline parameterization with no negative Jacobian after applying adaptive refinement; (f) mean ratio Jacobian of the initial parametrization; (g) mean ratio Jacobian of the final parametrization. mean ratio Jacobian less than a certain threshold δ. The refined T-mesh is optimized again and the process is repeated until a satisfactory mesh quality is obtained. Fig. 10 illustrates the efficiency of the proposed strategy. Additional refinements were applied to Gran Canaria Island domain with δ=0.2.The initial mesh with 3439 cells produces a Tspline parametric mapping with low quality in some areas and negative Jacobian in the North East part of the island. After adaptive refinement we have a mesh with 3577 cells and positive Jacobian in all the domain. Moreover, the minimum value of mean ratio Jacobian at the quadrature points is 0.21. 16 (a) (b) Fig. 20 Modification of the objective function in order to obtain, when possible, the orthogonality of the physical local submesh with hanging node. (a) Physical mesh and the optimal position of free node for the unweighted objective function. (b) The optimal position (x0,y0)of free node for the weighted objective function. Hanging node case In order to determine the weights of the objective function for the hanging node case, we pose the following problem. For the given physical local submesh, as shown in Fig. 20(a), we have to find the appropriate weights τ1and τ2 such that the optimal position of free node, according to the weighted objective function Kτ, will produce an orthogonal local mesh, as shown in Fig. 20(b). To do that, we evaluate the derivatives of the objective function Kτand enforce this derivatives to be equal to zero at (x0,y0). Objective function Kτis defined as Kτ(x,y) = τ1 3 ∑ i=1 ηi(x,y)+τ1 6 ∑ i=4 ηi(x,y)+τ2 11 ∑ i=7 ηi(x,y). The derivatives of the function η(S(x,y)) = kSk2 2det(S)are ∂xη=η(S)h∂xS,Si kSk2−∂xdet(S) 2det(S), ∂yη=η(S)h∂yS,Si kSk2−∂ydet(S) 2det(S), where the inner product h·,·iis defined as hA,Bi=Tr(AT,B). The derivative at (x0,y0)with respect to yis equal to zero due to the symmetry of the problem, and the derivative with respect to xis given by ∂xη(x0,y0) = −(a2−b2)(5τ2−8τ1) 2a2b. Thus, the point (x0,y0)will be the minimum, independently of values aand b, if τ2=8 5τ1. The weights we have used in this work are τ1=1 and τ2=8 5. It should be pointed out that the case of hanging node is a bit specific and its treatment is not completely analogous to the treatment of the regular node case. References 1. Bazilevs Y, Calo VM, Cottrell JA, Evans JA, Hughes TJR, Lipton S, Scott MA, Sederberg TW (2008) Isogeometric analysis: Toward unification of computer aided design and finite element analysis. In: Trends in Engineering Computational Technology, Saxe-Coburg Publications, Stirling, pp 1–16 2. Bazilevs Y, Calo VM, Cottrell JA, Evans JA, Hughes TJR, Lipton S, Scott MA, Sederberg TW (2010) Isogeometric analysis using T-splines. Comput Meth Appl Mech Eng 199:229–263 3. Bornemann PB, Cirak F (2013) A subdivision-based implementation of the hierarchical b-spline finite element method. Comput Meth Appl Mech Eng 253:584– 598 4. Casc´ on JM, Montenegro R, Escobar JM, Rodr´ ıguez E, Montero G (2007) A new meccano technique for adaptive 3-D triangulation. In: Proc. of the 16th International Meshing Roundtable, Springer, Berlin, pp 103–120 5. Casc´ on JM, Montenegro R, Escobar JM, Rodr´ ıguez E, Montero G (2009) The meccano method for automatic tetrahedral mesh generation of complex genuszero solids. In: Proc. of the 18th International Meshing Roundtable, Springer, Berlin, pp 463–480 6. Coons SA (1964) Surfaces for computer aided design. Springfield 7. Cottrell JA, Hughes TJR, Bazilevs Y (2009) Isogeometric Analysis: Toward Integration of CAD and FEA. John Wiley & Sons, Chichester 8. Dennis J, Schnabel R (1983) Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Classics in Applied Mathematics, Society for Industrial and Applied Mathematics, Englewood Cliffs, New Jersey, URL http://books.google.es/books?id= RtxcWd0eBD0C 9. Escobar JM, Rodr´ ıguez E, Montenegro R, Montero G, Gonz´ alez-Yuste JM (2003) Simultaneous untangling and smoothing of tetrahedral meshes. Comput Meth Appl Mech Eng 192:2775–2787 10. Escobar JM, Montenegro R, Montero G, Rodr´ ıguez E, Gonz´ alez-Yuste JM (2005) Smoothing and local refinement techniques for improving tetrahedral mesh quality. Computers & Structures 83:2423–2430 11. Escobar JM, Rodr´ ıguez E, Montenegro R, Montero G, Gonz´ alez-Yuste JM (2010) SUS Code: Simultaneous mesh untangling and smoothing code. http://www.dca.iusiani.ulpgc.es/proyecto20122014/html/Software.html 12. Escobar JM, Casc´ on JM, Rodr´ ıguez E, Montenegro R (2011) A new approach to solid modeling with trivariate T-splines based on mesh optimization. Comput Meth Appl Mech Eng 200:3210–3222 13. Escobar JM, Montenegro R, Rodr´ ıguez E, Montero G (2011) Simultaneous aligning and smoothing of surface triangulations. Engineering with Computers 27:17–29 14. Escobar JM, Montenegro R, Rodr´ ıguez E, Casc´ on JM (2012) The meccano method for isogeometric 17 solid modeling and applications. Engineering with Computers pp 1–13, published online, DOI 10.1007/ s00366-012-0300-z 15. Farin G, Hansford D (1999) Discrete Coons patches. Comput Aid Geom Design 16:691–700 16. Freitag LA, Knupp PM (2002) Tetrahedral mesh improvement via optimization of the element condition number. Int J Num Meth Eng 53:1377–1391 17. Freitag LA, Plassmann P (2000) Local optimizationbased simplicial mesh untangling and improvement. Int J Num Meth Eng 49:109–125 18. Knupp PM (2001) Algebraic mesh quality metrics. SIAM J Sci Comput 23:193–218 19. Knupp PM (2003) A method for hexahedral mesh shape optimization. Int J Num Meth Eng 58(2):319–332 20. Li B, Li X, Wang K (2010) Generalized polycube trivariate splines. In: SMI 2010, International Conference of Shape Modeling and Applications, pp 261–265 21. Li X, Guo X, Wang H, He Y, Gu X, Qin H (2007) Harmonic volumetric mapping for solid modeling applications. In: Proc. of ACM Solid and Physical Modeling Symposium, Association for Computing Machinery, Inc., pp 109–120 22. Martin T, Cohen E, Kirby R (2009) Volumetric parameterization and trivariate b-spline fitting using harmonic functions. Comput Aid Geom Design 26:648–664 23. Montenegro R, Casc´ on JM, Escobar JM, Rodr´ ıguez E, Montero G (2009) An automatic strategy for adaptive tetrahedral mesh generation. Appl Num Math 59:2203– 2217 24. Montenegro R, Casc´ on JM, Rodr´ ıguez E, Escobar JM, Montero G (2010) The meccano method for automatic three-dimensional triangulation and volume parametrization of complex solids. In: Developments and Applications in Engineering Computational Technology, Saxe-Coburg Publications, Stirling, pp 19–48 25. Piegl L, Tiller W (1997) The NURBS book. Springer, New York 26. Samet H (2006) Foundations of Multidimensional and Metric Data Structures. Morgan Kaufmann Publishers, Burlington, Massachusetts 27. Schillinger D, Deb´ e L, Scott M, Evans JA, Borden MJ, Rank E, Hughes TJR (2012) An isogeometric designthrough-analysis methodology based on adaptive hierarchical refinement of nurbs, immersed boundary methods, and T-spline cad surfaces. Comput Meth Appl Mech Eng 28. Scott MA, Li X, Sederberg TW, Hughes TJR (2012) Local refinement of analysis-suitable T-splines. Comput Meth Appl Mech Eng 213-216:206–222 29. Verf¨ urth R (1996) A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Technique. Wiley-Teubner, Chichester 30. Vuong AV, Giannelli C, Juttler B, Simeon B (2011) A hierarchical approach to adaptive local refinement in isogeometric analysis. Comput Meth Appl Mech Eng 200:3554–3567 31. Wang W, Zhang Y, Liu L, Hughes TJR (2013) Trivariate solid T-spline construction from boundary triangulations with arbitrary genus topology. Computer-Aided Design 45:351–360 32. Xu G, Mourrain B, Duvigneau R, Galligo A (2011) Parametrization of computational domain in isogeometric analysis: Methods and comparison. Comput Meth Appl Mech Eng 200:2021–2031 33. Xu G, Mourrain B, Duvigneau R, Galligo A (2011) Variational harmonic method for parameterization of computational domain in 2D isogeometric analysis. In: 12th International Conference on Computer-Aided Design and Computer Graphics, IEEE, Jinan, pp 223–228 34. Zhang Y, Wang W, Hughes TJR (2012) Solid T-spline construction from boundary representations for genuszero geometry. Comput Meth Appl Mech Eng 249252:185–197