Construction of polynomial spline spaces over T-meshes for its application in isogeometric analysis

UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA DOCTORAL DISSERTATION Construction of polynomial spline spaces over T-meshes for its application in Isogeometric Analysis Marina Brovka Las Palmas de Gran Canaria •June 2016 Discretization and Applications Division G,,, DE LAS P,l,W.S DE G"" CANAR" Instituto Universitario de Sistemas Inteligentes y Aplicaciones Numéricas en Ingeniería DAVID GREINER SÁNCHEZ, PROFESOR TITULAR DE UNIVERSIDAD Y SECRETARIO DEL INSTITUTO UNIVERSITARIO DE SISTEMAS INTELIGENTES Y APLICACIONES NUMÉRICAS EN INGENIERÍA (SIANI) DE LA UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA, CERTIFICA Que el Consejo de Doctores del Instituto Universitario de Sistemas Inteligentes y Aplicaciones Numéricas en Ingeniería (SIANI) en su sesión de fecha 27 de mayo de 2016 tomó el acuerdo de dar el consentimiento para su tramitación a la tesis doctoral titulada "Construction of polynomial spline spaces over T-meshes for its application in Isogeometric Analysis. Construcción de espacios spline polinómicos sobre T-meshes para su aplicación en Análisis Isogeométrico." presentada por la doctoranda Dña. Marina Brovka y dirigida por Dr. D. José María Escobar Sánchez y Dr. D. Rafael Montenegro Armas, a la vista de la idoneidad y calidad de su contenido, interés y relevancia del tema a nivel internacional. Para que así conste, y a los efectos oportunos se expide el correspondiente certificado a treinta y uno de mayo de dos mil dieciséis. Construction of polynomial spline spaces over T-meshes for its application in Isogeometric Analysis Programa de doctorado: Sistemas Inteligentes y Aplicaciones Num´ericas en Ingenier´ıa Instituto Universitario SIANI Autor: Marina Brovka Director: Jos´e Mar´ıa Escobar S´anchez Director: Rafael Alejandro Montenegro Armas Las Palmas de Gran Canaria, June 2016 To My Parents Acknowledgements I would like to express my deepest gratitude to my supervisors Prof. Rafael Montenegro Armas and Prof. Jos´e Mar´ıa Escobar S´anchez for their guidance and patience. My sincere thanks for giving me the opportunity to work in their research group, for their trust, continuous support and friendly attitude that contributed a lot to my motivation. I am very thankful to Prof. Jos´e Manuel Casc´on Barbero for his help, collaboration and numerous advices and suggestions. I am very happy and grateful for having my co-worker Jos´e Iv´an L´opez Gonz´alez, for our fruitful collaboration and for being a team, which made possible our research. Of cause, I would like to express my gratitude and appreciation for all the help, support and nice working atmosphere to all members of our research group: Dr. Albert Oliver Serra, Guillermo Socorro, Dr. Eduardo Rodr´ıguez Barrera, Jabel Ram´ırez Naranjo and Prof. Gustavo Montero Garc´ıa. I would like also to acknowledge The University of Las Palmas de Gran Canaria for their financial support (“Programa de FPI propio de la ULPGC”) and all the staff of the university institute SIANI for giving opportunities and facilitating the incorporation to scientific research. CHAPTER 1 Introduction 1.1 Isogeometric Analysis concept Computer-aided design (CAD) is the process of creating a design using computer software. After CAD has emerged (1960-1970) all engineering drawings, done manually using pencil and paper, was replaced by files describing the physical model. That entailed a huge advance in all engineering industry. Nowadays CAD is a widely used technology with numerous applications in automotive, shipbuilding, aerospace industries, industrial and architectural design, computer animation for special effects in cinema. The most commonly used mathematical tool for representing curves and surfaces in CAD are B-splines and NURBS (Non-Uniform Rational B-splines) [1, 2]. They offer a convenient free-form surface modelling and can exactly represent all conic sections. Efficient algorithms to handle spline object are available. Besides, B-splines possess some additional useful properties as higher smoothness and convex hull property. On the other hand, another important computational technique that has been under intensive development since its emerging (1950-1960) is Finite Element Method (FEM). This numerical method for solving partial differential equations has received a big amount of research effort and became a standard tool in engineering industry. The basic components of the method are the variational formulation of the physical problem and the discrete space used to approximate its solution. The approximation space is defined by its basis functions (Lagrange or Hermite interpolating polynomials). Each basis function is defined locally on its element. To perform the analysis, the computational domain of the problem should be decomposed in non-overlapping simple elements (triangles, quadrilaterals, tetrahedra, hexahedra). This step is referred to as mesh generation process. Taking into account that nowadays the geometry comes from CAD, there is a necessity to generate an analysis-suitable mesh of the object from its CAD representation. It is not a trivial task due to increasing complexity of engineering design. Finite Element Analysis (FEA) and CAD technologies have evolved separately and they use different geometrical tool. That impedes the communication between them 1 1 Introduction (a) (b) Figure 1.1: (a) Head NURBS model, (b) Rhinocero NURBS model. Images taken from http: // www. 3drender. com/ and make the process of data transferring very time-consuming. Creation of the analysis-suitable model and mesh generation process consumes the major part of all analysis time. Tight communication between CAD model and FEA geometry is also crucial to make possible design optimization. In addition, finite element mesh does not reproduce the exact model designed by CAD. For some problems geometry imprecision can lead to important error in numerical solution. To approximate a curved boundary with a desired precision an adaptive mesh refinement is required. This is possible only with an automatic interaction between exact CAD geometry and FEM mesh. With all these inconveniences it became evident the necessity to unify the design and the analysis in a unique engineering process. CAD has to provide directly a geometric model suitable for analysis, or it is necessary to have an automatic method to obtain analysis suitable model from its CAD representation. For that, the analysis process should be changed and adapted for its application on CAD geometries. The thought was proposed in 2005 by Tom Hughes and coworkers in [3]. The concept received the name Isogeometric Analysis (IGA). The idea is to use for the analysis the same basis functions (originally NURBS) that have been used for construction of the CAD model. Isogeometric Analysis can be seen as a generalization of FEM that uses basis functions with higher regularity. Since its introduction IGA attracted a lot of attention in research community and it has been an object of numerous investigation works in the last years. At present it is considered a promising tool for bridging the gap between CAD and FEA industries. Besides, it can offer some additional beneficial properties and possibilities compared with classical Finite Element Method. Here are some of them: •Working with the exact geometry leads to superior accuracy in the problems sensitive to geometry approximation (contact problems, boundary layer prob2 1 Introduction lems in aerodynamic and hydrodynamic). Even if the exact geometry is not achieved, smooth boundary approximation is necessary for some problems that require curvature continuity of the boundary. •Better accuracy of numerical solution due to higher regularity of the basis functions [4]. •Higher degree of continuity also allows the use of IGA for solving partial differential equations of order greater than two without the necessity to change variational formulation. For example, Cahn-Hilliard phase-field equation [5, 6, 7, 8]. •Possibility of collocation methods using strong formulation [9, 10, 11, 12]. However, to accomplish the ambitious idea of Isogeometric Analysis and make it work in practice, some issues should be resolved first. This question is revised in the next section. 1.2 State of the art IGA is relatively new method that needs a theoretical framework to be developed and some open problems to be solved. Clearly, a lot of question should be studied: error estimation theory for IGA [13, 14, 15], efficient quadrature rules [16, 17], efficiency of the direct and iterative solvers for the systems arising in the method [18, 19], efficient implementation of the structures and procedures arising in the method, imposing boundary conditions, the influence of parameterization quality on accuracy of the solution [20]. However, the most urgent and important problems of IGA are the local refinement problem and the parameterization of computational domain. 1.2.1 Local refinement Originally Isogeometric Analysis concept was proposed and tested using NURBS basis functions inasmuch as NURBS was the most common tool used by CAD software at that moment. However, it suffers from an important drawback: their tensor product structure does not allow for local refinement. NURBS surface is defined by a set of control points which lie, topologically, in a rectangular grid. A new control point insertion induces the insertion of an entire row of control points extended through the entire domain. That leads to a large number of superfluous control points. The representation of local features is inefficient and requires the use of several NURBS patches joined together. Frequently this join is discontinuous, i.e. it has gaps and overlaps that make the model not suitable for analysis. T-splines were proposed by Sederberg et al. [21] as an alternative to NURBS that permits local refinement. T-splines can be seen as a generalization 3 1 Introduction (a) “Head modeled (a) as a NURBS with 4712 control points and (b) as a T-spline with 1109 control points. The red NURBS control points are superfluous” (b) “NURBS head model, converted to a T-spline” Figure 1.2: Images taken from the original work of Sederberg et al. [22]: “T-spline simplification and local refinement”. of NURBS. T-splines allows the insertion of a partial row of control points that terminates in T-junction, which makes them locally refinable. T-spline plug-in are currently available for two NURBS-based CAD software Maya and Rhino, see http://www.tsplines.com/products/tsplines-for-rhino.html. T-splines offer a flexible tool for creating one-patch watertight surfaces with local detailed features and lower number of control points compared with NURBS, see Fig. 1.2. However, it appears that initially introduced T-splines are not suitable for their use in Isogeometric Analysis since they lack some properties essential for analysis and for proper convergence behaviour: linear independence, polynomial reproduction property and nesting behaviour of approximation spaces (this question will be illustrated in more details in Chapter 2. Besides, rational T-splines complicate and increase computational cost of the calculus of the derivatives. Therefore, it is still an open issue in the context of Isogeometric Analysis to find an alternative to NURBS that overcomes the local refinement problem and can be used in analysis. This issue has been the object of numerous research works in recent years. Analysis-suitable T-splines, proposed by Scott et al. in [23], are a class of T-splines defined over T-meshes that meet certain topological restrictions formulated in terms of T-junction extensions. Blending functions defined over an extended analysis-suitable T-mesh are linearly independent [24] and possess polynomial reproduction property. The refinement algorithm allows to accomplish highly localized refinements and construct nested T-spline spaces, but it presents an elevated implementation complexity and, as far as we know, the generalization of the 4 1 Introduction strategy to 3D cases is still an open question. Another well known approach for local enrichment of the approximation space is the hierarchical refinement, originally introduced by Forsey and Bartels in [25] and later developed in [26]. Recently, hierarchical refinement technique in the context of Isogeometric Analysis was described in [27, 28, 29]. This approach is based on a simple and natural idea to construct multilevel spaces by replacing coarse level functions with finer basis functions. Starting from an initial uniform mesh, hierarchical refinement scheme leads to sequential construction of nested spline spaces with linearly independent basis functions. Relatively simple implementation and straightforward generalization to 3D make it an attractive option for local refinement. However, a shortcoming of this strategy is the impossibility to define a spline space over a given arbitrary T-mesh, as well as the presence of redundant basis functions and excessive support overlapping. An interesting theoretical approach to the latter problem was given in [30]. The truncation technique is applied to redefine the function supports and reduce their overlapping. However, that elevates considerably the computational cost of the strategy. Other strategies for performing local refinement of spline spaces are C1continuous PHT-splines [31] and Locally Refined B-splines (LR-splines) [32]. It is worth mentioning another group of options for defining multivariate spline functions that do not have tensor product structure. Namely, bivariate simplex splines over triangulation [33]; simplex splines over Delaunay configurations, as a natural generalization of univariate B-splines, was proposed in [34, 35, 36]; box splines [37]; quadratic C1−continuous splines over Powell-Sabin triangle split [38]. Recent examples of isogeometric analysis using splines over triangulations can be found in [39, 40, 41]. Splines over triangulation could offer a better flexibility and adaptivity to irregular domains compared to more rigid tensor product splines. However, it is not a trivial task to define globally smooth spline spaces over triangulations. Currently available options presents a high complexity, and it seems that tensor product splines are still the most popular and standard option for CAD due to their simplicity. 1.2.2 Domain parameterization CAD models usually provide only the boundary surface of a solid. This spline representation of the boundary can be used directly in isogeometric shell analysis or isogeometric boundary integral method. But the application of Isogeometric Analysis requires a full volumetric representation of the geometry. So, another open problem of Isogeometric Analysis is how to generate a trivariate spline representation of a solid starting from the CAD description of its boundary. As it is pointed by Cottrell et al. in [42], “the most significant challenge facing Isogeometric Analysis is developing three-dimensional spline parameterizations from surfaces”. Parameterization is suitable for analysis if it does not have self-intersection, i.e. it is invertible. Moreover, in order to expect a high accuracy in numerical results it 5 1 Introduction is necessary to obtain a good quality volume parameterization. Orthogonality and uniformity of isoparametric curves are desirable for the tensor-product structured parameterization. It is not trivial task to obtain a good quality smooth global parametric mapping for complex domains and it can be very time-consuming. For the application of IGA it is essential to have an effective method to construct a spline parameterization. Here are some attempts to tackle the problem. In [43] the parameterization is found by solving a constraint optimization problem for a planar B-spline surface. Constraints are defined by imposing injectivity sufficient conditions in terms of the control points, and the optimization consists in minimizing some energy function in order to reach a good orthogonality and uniformity of the parametric mapping. The idea was extended for 3D in [44]. Another similar technique was proposed by these authors in [45, 46]. 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. For example, Li et al. [47] construct a harmonic volumetric mapping through a meshless procedure by using a boundary method. The algorithm can be applied to any genus data, but it is complex and requires placing some source and collocation points on an offset surface. Optimal results of source positions are unknown, and in practice they are chosen in a trial-and-error manner or with the help of human experience. Martin et al. [48] present a methodology based on discrete harmonic functions to parameterize a solid, where the input data are surface triangulation and a tetrahedral mesh of the solid. They solve several Laplace equations, first on the surface to establish surface parameterization and then on the complete 3D domain using FEM. The two obtained orthogonal harmonic functions are used to construct a structured hexahedral mesh of the solid, which is smoothed to remove irregularities. The user has to make an initial choice of two critical points to establish the surface parameterization and to fix a seed for generating the skeleton. Zhang et al. proposed in [49] a procedure to construct T-spline representation of a genus-zero solid from its boundary triangulation. First, a parametric mapping between the triangulation and the boundary of the parametric domain, a unit cube, is established using Floater parameterization method. Then an octree subdivision is carried out for the cube until the error between the T-mesh and the input triangle mesh is less than a threshold. During this process, the boundary nodes are mapped to the input triangle mesh and, then, the interior nodes are relocated via T-mesh optimization, which is maximizing the worst Jacobian of the trilinear map for each element. In [50] the method was extended for arbitrary genus topology. A smooth harmonic scalar field over the surface triangulation is computed, and saddle points are extracted to determine the topology of the object and to construct a polycube with the same topology that serves as the parametric domain for the trivariate T-spline. 6 1 Introduction 1.3 Goal and outline of the thesis This thesis is the result of our group research work on Isogeometric Analysis. As was mentioned before, it is relatively new method with many questions to explore. A lot of new interesting things about IGA have been studied and discovered since we have embarked on this journey. Some of these things will be reflected implicitly in the present work. But, mainly, our attention was focused on two open problems of IGA: construction of spline spaces with nice properties suitable for analysis and the problem of volumetric parameterization of computational domain from its boundary representation. Some preliminary results about parameterization method developed in our research are briefly included in this work, however this problem will be studied in depth and presented in another doctoral dissertation of the group. This thesis addresses the issue of construction of appropriate spline spaces over T-meshes for its application in Isogeometric Analysis. The main contribution of the dissertation is a new strategy for defining spline blending functions that span spaces with nice properties. The technique is designed for hierarchical T-meshes (multilevel meshes) with a quadand octree subdivision scheme. This type of meshes can be efficiently implemented with tree data structures [51], which are frequently used in engineering. Due to the elevated complexity of all current strategies, the main goal we pursue here is the simplicity and low computational cost, both in 2D and 3D. For that, we have to assume a restriction on the T-mesh. Namely, the T-mesh should fulfil the requirement of being a strongly balanced mesh. Assuming this reasonable and frequently used restriction over the T-mesh, we can define easily cubic spline functions that span spaces with desirable properties: linear independence, C2-continuous, polynomial reproduction property, nested spaces and a straightforward implementation. The manuscript is organized as follows. •Chapter 2 introduces some basic concepts about B-splines and NURBS and their use for curve and surface modelling, followed by the description of Isogeometric Analysis concept and its comparison with classical Finite Element Method. •In Chapter 3 the idea of T-splines is revised and some of their shortcomings are discussed and illustrated. Then, we discuss briefly the available solutions for the problem of local refinement and analyse their advantages and shortcomings that has motivated our research. •Chapter 4. The core of the dissertation is presented. We describe and illustrate the procedure for the construction of cubic spline space over a given T-mesh. The key of our technique is some simple rules used for inferring local knot vectors to define tensor product spline blending functions that span 7 1 Introduction spaces with desirable properties. The algorithms for an efficient implementation of the strategy for 2D and 3D meshes are given. •Properties and some characteristics of the constructed spaces are discussed in Chapter 5. Proof of some properties for 2D case is provided. •Chapter 6 presents a brief description of another result of our research: a method for spline parameterization of complex 2D geometries, which will be necessary for the computational examples. •Chapter 7. The approximation properties of the proposed spline spaces are tested in geometric modelling and Isogeometric Analysis. We apply IGA for different types of problems that involve adaptive refinement. Convergence behaviour is analysed. Some computational examples include a comparison with Finite Element Method. •Chapter 8 exposes the summary, concluding remarks and future research directions. 8 CHAPTER 2 Preliminaries In this chapter we introduce some basic concepts about spline theory relevant to this thesis, namely, B-splines and NURBS and their use for curve and surface modelling and Isogeometric Analysis. For more comprehensive review we refer the reader to [1, 52, 2]. Then, the Isogeometric Analysis concept is explained in details. 2.1 Spline theory. Basic concepts 2.1.1 Spline interpolation Data fitting is a fundamental concept for geometric design. A common problem in curve modelling is data interpolation: given data points piand their corresponding parameter values ξi, find a parametric curve that passes through the points pi. A well known method to solve this problem is polynomial interpolation, that is, find a n-order polynomial Pnthat satisfy interpolation conditions Pn(ξi) = pi, i = 0, . . . , n. However, a single segment polynomial interpolation have some shortcomings: a large number of data points requires a high degree polynomial which is inefficient and can present oscillations, known as Runge’s phenomenon, and numerical unstability (see [2] for an example of ill-conditioned Lagrange polynomial interpolation). A solution to this problem is piecewise polynomial interpolation. The idea is to construct a polynomial segments Ci(ξ) for each parameter interval [ξi−1, ξi] so that the resulting piecewise curve C(ξ) passes through the data points piand the polynomial segments join with a certain level of continuity, i.e., C(k) i(ξi) = C(k) i+1(ξi). This piecewise polynomial curve C(ξ) is called spline, see Fig. 2.1. The term spline comes from the name of flexible rulers that were used for technical drawings and were free to bent to pass through the specified points. Spline curve can be seen as a mathematical generalization of the physical spline. The set of all piecewise polynomial curves of degree nand some level of 9 2 Preliminaries A new knot ˆ ξ∈[ξk, ξk+1) is inserted, creating a new knot vector b Ξ = {ξ1, ξ2, . . . , ξk,ˆ ξ, ξk+1, . . . , ξn+p+1}. It is possible then to construct the same curve C(ξ) using the new basis {b Bi,p}n+1 i=1 defined on b Ξ: C(ξ) = n ∑ i=1 PiBi,p(ξ) = n+1 ∑ i=1 b Pib Bi,p(ξ).(2.2) The new control points b Piare calculated as b Pi=αiPi+ (1 −αi)Pi−1, where αi=       1,if i≤k−p ˆ ξ−ξi ξi+p−ξi ,if k−p+ 1 ≤i≤k 0,if i≥k+ 1 (2.3) Note that only some of the control points are changed with respect to the original control points Pi.When a new knot ˆ ξis inserted into a knot vector Ξ, the initial basis functions {Bi,p}n i=1 can be expressed in terms of the new basis functions {b Bi,p}n+1 i=1 defined on b Ξ.That is, a knot insertion produces nested spline spaces. Each function Bi,p is expressed as linear combination of the refined functions. Inserting these expressions into (2.2) leads to the formula (2.3) for the new control points. 2.1.3.1.2 Degree elevation. p-refinement When elevating the polynomial order of the basis functions, the inter-element regularity should be preserved. So the multiplicity of the original knots should be increased by one before degree elevation. No new knots are inserted. In this way, degree elevation also produces nested spaces and the original curve C(ξ) can be preserved. 16 2 Preliminaries 2.1.4 B-spline surfaces Bivariate B-splines are defined as a tensor product of univariate B-spline functions. Figure 2.8: Bivariate B-spline function. Let {Bi,p(ξ)}n1 i=1 and {Bj,q(η)}n2 j=1 be two sets of B-spline basis functions of degree pand q, defined over knot vectors Ξ = {ξ1, ξ2, . . . , ξn1+p+1}and H= {η1, η2, . . . , ηn2+q+1}, respectively, which form a Cartesian grid in parametric domain. Then, bivariate basis functions are defined as follows B(i,j),(p,q)(ξ, η) = Bi,p(ξ)Bj,q(η),(i, j)∈I, where the multi-index set Iis defined by I={1,2, . . . , n1} × {1,2, . . . , n2}. See an example of bivariate B-spline function in Fig. 2.8 and an example of a set of B-spline functions defined over a parametric grid in Fig. 2.9. (a) (b) Figure 2.9: Bivariate cubic B-spline functions defined over a parametric grid formed by the knot vectors Ξ = H={0,0,0,0,1/4,1/2,3/4,1,1,1,1}. Bivariate B-spline functions possess the same properties as the univariate Bsplines: 17 2 Preliminaries Figure 2.10: B-spline surface and its control net. •Local support:B(i,j),(p,q)(ξ, η) = 0 for ξoutside the rectangle [ξi, ξi+p+1)× [ηj, ηj+q+1). •Non-negativity:B(i,j),(p,q)(ξ, η)≥0. •Partition of unity:∑(i,j)∈IB(i,j),(p,q)(ξ, η) = 1,(ξ, η)∈[ξ1, ξn1+p+1]×[η1, ηn2+q+1). •Linear independence:∑(i,j)∈Ici,jB(i,j),(p,q)≡0⇒ci,j = 0,(i, j)∈I. •Except for the case p= 0 or q= 0,B(i,j),(p,q)attains exactly one maximum. •Each function B(i,j),(p,q)is a piecewise bivariate polynomial. It is (p−m) times differentiable in the ξ-direction, where mis the multiplicity of the ξ knot. Analogously for η-direction. A B-spline surface is defined as a linear combination of bivariate B-spline functions S(ξ, η) = ∑ (i,j)∈I Pi,j B(i,j),(p,q)(ξ, η), where the control points Pi∈R3form a control net. See an example of B-spline surface in Fig. 2.10. Analogously to B-spline curve, the B-spline surface does not interpolate the control points, except for the four corner control points. If triangulated, the control net forms a piecewise planar approximation to the surface. Also, B-spline surface conserves some properties of the spline curve: convex hull property, affine invariance and local control property, but no variation diminishing property. Remark 1. B-spline solids are defined in analogous way, as a linear combination of trivariate B-spline basis functions and inherit the properties of its univariate and bivariate counterpart. 18 2 Preliminaries 0 1 4 1 2 3 4 1 1 Figure 2.11: Example of NURBS functions derived from B-splines using the weights w1= 2, w2= 1, w3= 3, w4= 4, w5= 1, w6= 2 and w7= 1. 2.1.5 Non-Uniform Rational B-splines (NURBS) Although polynomial splines offer many advantages for the free-form modelling, there exist some important shapes which cannot be represented exactly with Bspline functions, e.g., circles, ellipses, cylinders, cones, spheres. This inconvenience can be solved using rational functions. A set of Non-Uniform Rational B-splines functions {Ri,p}for a given knot vector Ξ = {ξ1, ξ2, . . . , ξn+p+1}is defined from a set of B-spline basis functions {Bi,p}as follows Ri,p(ξ) = wiBi,p(ξ) ∑n j=1 wjBj,p(ξ), where {wi}is a set of npositive weights, see Fig. 2.11. Note that if wj= 1 for all j, then Ri,p =Bi,p. So, polynomial B-spline functions are special case of the NURBS functions {Ri,p}. In a similar fashion to B-spline curves, NURBS curve is defined as linear combination of the rational basis functions C(ξ) = ∑ i∈I PiRi,p(ξ), where Pi∈R2are the control points. A NURBS surface is defined as a linear combination of bivariate NURBS functions S(ξ, η) = ∑ (i,j)∈I Pi,j R(i,j),(p,q)(ξ, η), where R(i,j),(p,q)(ξ, η) are bivariate rational basis functions defined as R(i,j),(p,q)(ξ, η) = wi,jB(i,j),(p,q)(ξ, η) ∑(i,j)∈Iwi,j B(i,j),(p,q)(ξ, η). 19 2 Preliminaries parametric space index space Figure 2.12: Index space. Shaded zone corresponds to non-zero knot spans. Choosing an appropriate knot vector and values of the weights {wi}for NURBS enables to represent exactly conic sections. NURBS and NURBS curves (surfaces) inherit the properties of B-spline functions and B-spline curves(surfaces), respectively. For more details about NURBS we refer the reader to [1]. 2.1.6 Index space and parametric space Since some knots of the parametric space can have the same values, for some purpose it is useful to use the index space. In contrast to parametric space, where multiple knots are plotted as the same point, the index space is a uniform grid, where each knot is depicted as the unique knot, taking into account only the index of the knot. In the index space all knot intervals have the same length, while in parametric space some intervals can have zero length. Figure 2.12 illustrates an example of parametric grid and its corresponding index grid. Knot vector Ξ contains an interior multiple knot ξ6=ξ7= 0.5 and repeated boundary knots. Knot vector Hhas repeated knots only on the boundary. 20 2 Preliminaries 2.2 Isogeometric Analysis In this section we revise the idea of Isogeometric Analysis. Assuming that the reader is familiar with the Finite Element Method, we focus our attention on the similarities and the main differences between IGA and Finite Element Analysis. The main idea of Isogeometric Analysis is to use for the solution space the same basis functions that are used to model the geometry. This notion is called isoparametric concept, and it is common for IGA and FEA. However, “The fundamental difference between this new concept of Isogeometric Analysis and the old concept of isoparametric finite element analysis is that, in classical FEA, the basis chosen to approximate the unknown solution fields is then used to approximate known geometry. Isogeometric Analysis turns this idea around and selects a basis capable of exactly representing the known geometry and uses it as a basis for the fields we wish to approximate. In a sense, we are reversing the isoparametric arrow such that it points from the geometry toward the solution space, rather that vice versa.” (Cottrell et al. in [42].) 2.2.1 Variational formulation and Galerkin’s projection method Theoretical framework for both, Finite Element Analysis and Isogeometric Analysis, is based on two main steps: variational (weak) formulation of the problem and Galerkin’s projection method to approximate the solution of the weak problem. Next, we revise the method using the Poisson model problem. Let us consider a strong form of the boundary value problem −∆u=fin Ω, u=gon ΓD, ∇u·n=hon ΓN, (2.4) where Ω is a Lipschitz domain with boundary ∂Ω=Γ=ΓD∪ΓN, and the functions f: Ω →R, g : ΓD→R, h : ΓN→Rare given. Let Vbe a suitable Hilbert space V=H1(Ω) := W1 2(Ω),where the Sobolev space Wk 2(Ω) = {v∈L2(Ω) : ∂αv∈L2(Ω),∀ |α| ≤ k}. Then, the space V= H1(Ω) is equipped with the norm ∥v∥V= ∑ |α|≤1∫Ω (∂αv)2dΩ  1/2 and seminorm |v|V= ∑ |α|=1 ∫Ω (∂αv)2dΩ  1/2 . 21 2 Preliminaries The test function space V0(Ω) is defined as V0(Ω) = {v∈H1(Ω) : v|ΓD= 0}, and the solution space Vg(Ω) is the set of functions from Vwhich fulfil the essential boundary conditions, i.e., Vg(Ω) = {v∈H1(Ω) : v|ΓD=g}. If a function usatisfies the equation (2.4), it also satisfies ∫Ω−∆u v dΩ = ∫Ω f v dΩ,∀v∈V0(Ω). Applying the Green’s formula for integration by parts leads to ∫Ω∇u·∇vdΩ = ∫Ω f v dΩ + ∫ΓN h v dΓ + ∫ΓD ∂u ∂nvdΓ,∀v∈V0(Ω), where the last term vanishes. Then, variational formulation for the problem (2.4) consists in finding u∈ Vg(Ω) such that a(u, v) = F(v),∀v∈V0(Ω),(2.5) where a(·,·) is bilinear symmetric form a(u, v) = ∫Ω∇u·∇vdΩ and F(·) is a linear functional F(v) = ∫Ω f v dΩ + ∫ΓN h v dΓ. Under certain assumptions the variational problem (2.5) has a unique solution. Namely, Definition 1. A bilinear form a(·,·) on a normed vector space (V, ∥·∥) is coercive (elliptic), if there exist a constant α > 0 such that |a(u, u)| ≥ α∥u∥2,∀u∈V. Definition 2. A bilinear form a(·,·) is called bounded (continuous), if there exist a constant Msuch that |a(u, v)| ≤ M∥u∥∥v∥,∀u, v ∈V. Then, the existence and uniqueness of a solution for the problem (2.5) is guaranteed by the Lax-Milgram theorem. Theorem 1 (Lax-Milgram).Let Vbe a Hilbert space, a(·,·) : V×V→Ris a bounded, coercive bilinear form, and F(·) : V→Ris a continuous linear functional. Then, there exists a unique u∈Vsuch that a(u, v) = F(v),∀v∈V. 22 2 Preliminaries Next, we apply Galerkin method. To approximate the solution of the weak problem (2.5) we construct a finite-dimensional subspaces V0,h and Vg,h defined over the discretized domain Ωh.Let Vh= span{ϕ1, ϕ2, . . . , ϕN}be a finite-dimensional space spanned by basis functions defined over the discretization Ωh.Then, V0,h ⊂Vh is a subspace of functions that vanish on the boundary and Vgh,h is the subspace of functions that are equal to ghon the boundary, where ghis an interpolant of g. Then, finite-dimensional problem consists in finding uh∈Vgh,h such that a(uh, vh) = F(vh),∀vh∈V0,h(Ω).(2.6) Under the conditions of Theorem 1, there exists a unique uhthat solves the discrete problem (2.6). We search uh∈Vgh,h as uh= n ∑ i=1 uiϕi+gh,(2.7) where {ϕi}n i=1 ⊂V0,h.Substituting (2.7) into (2.6), Galerkin method gives rise to a linear algebraic system n ∑ i=1 uia(ϕi, ϕj) + a(gh, ϕj) = F(ϕj)∀j= 1, ..., n. That is, Ku =f,where K= [Kj,i]n j,i=1 is the stiffness matrix, being Kj,i =a(ϕi, ϕj); f= (F(ϕj)−a(gh, ϕj))n j=1 is the load vector, and u= (u1, u2, . . . , un)Tare the unknown coefficients called degrees of freedom. Under the assumptions of Theorem 1, the C´ea’s Lemma provides a bound of the error in norm ∥·∥Vbetween the exact solution and the approximate solution obtained by Galerkin projection: ∥u−uh∥V≤M αinf vh∈Vh∥u−vh∥V. The term infvh∈Vh∥u−vh∥Vmeasures the capacity of the finite dimensional space Vhto approximate the function u, and it is typically called an interpolation error. This error, in its turn, converges in terms of the mesh-size parameter hwith a certain order s, which depends on the smoothness of uand the degree of basis functions of the space Vh: inf vh∈Vh∥u−vh∥V≤C(u)hs. Thus, the Lemma asserts that the error of the finite element solution is of the same 23 2 Preliminaries (a) FEM basis. Lagrange functions. (b) IGA basis. B-spline functions. Figure 2.13: Comparison of basis functions for FEM and IGA. order that the interpolation error: ∥u−uh∥V≤M αC(u)hs. This bound on the error is called a priori error bound and it shows that finite element solutions uhconverges to uin H1-norm as h→0. 2.2.2 Basis functions The main difference between Finite Element Analysis and IGA lies in the type of basis functions used for discrete approximation space. In FEM the basis functions are defined locally on each element, which is key concept of the method. The discrete solution space Vhis spanned by the union of all basis functions defined on each element of the mesh, see Fig. 2.13(a). In IGA, however, the set of basis functions (B-splines, NURBS) is inferred from a global knot vector, which represents a discretization of the domain, Fig. 2.13(b). The knots partition the domain into the knot spans which is equivalent to the concept of element in FEA in the sense that it represents a discretization unit, where the functions are C∞-continuous. However, it is important to emphasize that a function support is not contained 24 2 Preliminaries Figure 2.14: Local parametric mapping for each element in Finite Element Method. just in one knot span, but it covers a several knot spans. The most important difference that provide spline basis functions, with respect to the functions used in classical Finite Element Method, is a higher continuity. In FEA any order basis functions are only C0across the element boundary. Meanwhile, pth-degree spline basis functions are Cp−1-continuous across the element boundary. Figure 2.13 illustrates this difference. For example, cubic B-splines possess C2-continuity (no repeated interior knots). It permits to construct numerical solution that is globally C2-continuous which is, in general, not possible with classical Lagrange type basis functions used in FEM. 2.2.3 Domain parameterization As a consequence of the global definition of the basis function in IGA, another difference, with respect to FEM, is the paremeterization of computational domain. In Finite Element Method each element of the physical domain has its own parametric mapping that transforms the reference element into the physical one. This local parameterization depends only on the nodes of the physical element and it is constructed with the basis functions defined on the reference element, see an illustration for 2D case in Fig. 2.14. For each physical element Ωka parametric mapping Sk:b Ω0→Ωkis defined as a linear combination Sk(ξ, η) = ∑ j Pk jˆ ϕj(ξ, η), where {ˆ ϕj}is the basis on the unique reference element b Ω0and Pk j∈R2.Then, the basis function ϕk ion the physical element Ωkis defined as ϕk i=ˆ ϕi◦S−1 k, where ˆ ϕiis the corresponding function on the reference element b Ω0. In Isogeometric Analysis, in contrast, the parameterization of computational domain is defined globally. A unique global parametric mapping S:b Ω→Ω 25 3 Local refinement problem this property. This can result in a poor approximation ability and convergence behaviour. Let us see an example illustrating this fact. We are going to use Tsplines to interpolate over a square domain Ω = [0,1]2a smooth function given by u(ξ, η) = exp −20((ξ−0.5)2+ (η−0.5)2),see Fig. 3.4(a). We define a set of rational bicubic T-spline functions {Rα}α∈ATfor a given mesh Tand interpolate the given function. The spline approximation is built as a linear combination of the blending functions uh(ξ) = ∑ α∈AT cαRα(ξ). The control values cαare found by imposing the interpolation conditions u(ξβ) = ∑ α∈AT cαRα(ξβ),∀ξβ, β ∈AT, As interpolation points we use the anchors of the functions, i.e., the vertices of the mesh. Adaptive refinement is performed according to the indicator based on the exact L2interpolation error: η(Ωe)2=∥u−uh∥2 L2(Ωe)=∫Ωe (u−uh)2dΩ. A cell Ωeis marked to be refined if η(Ωe)> γ maxi{η(Ωi)}, being γ∈[0,1]. Some steps of the adaptive refinement are shown in Fig. 3.4(b), and Fig. 3.4(c) illustrates the error convergence in L2-norm. As can be seen, no monotone convergence is obtained, since the error presents a significant jump for some mesh configurations. T-spline space corresponding to some meshes have poor approximation ability. To clarify the issue, the following study was carried out. Rational functions {Rα}are obtained from a set of polynomial spline blending functions {Bα}α∈ATby normalization. For each of these meshes we can check whether its blending functions {Bα}α∈ATspan a complete polynomial space over Ω. For that we interpolate some polynomial functions of degree 3 and check if the interpolation gains to represent exactly the polynomial (in the sense that the exact interpolation error ∥u−uh∥L2is zero). It was observed that some error increases on the Fig. 3.4(c) correspond to the meshes Twith “incomplete” space: P3(Ω) ⊂ ST= span{Bα}α∈AT.These meshes and their corresponding error are plotted in black colour in Fig. 3.4(b) and (c), respectively. The meshes that do span a complete polynomial space over Ω are marked in blue; and a proper convergence behaviour takes place for the sequence of these meshes. In Fig. 3.4(c) we also depict the error convergence for the same interpolation performed with polynomial spline functions proposed in this thesis, which, as can be seen, leads to monotone convergence. As another illustration of a poor approximation ability of T-splines, let us see an example of 2D parameterization. Suppose we are given a T-mesh T, a set of interpolation points {ξα}α∈AT∈R2in parametric space and their images in the 32 3 Local refinement problem (a) The function to approximate. (b) A sequence of T-meshes during the adaptive refinement. Polynomial Tmeshes (with polynomial reproduction property) are represented in blue. æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à ì ì ì ì ì 100 200 300 400 500 600 10-3 10-4 10-5 number of degrees of freedom L2error à our strategy æ T-splines (c) Error convergence, blue points corresponds to the meshes marked in blue. Figure 3.4: Example of using rational T-splines for interpolation. 33 3 Local refinement problem physical space {xα}α∈AT∈R2.We build a T-spline representation of the geometry as lineal combination of T-spline blending functions S(ξ) = ∑ α∈AT PαRα(ξ),(3.1) where Pα∈R2is the control point corresponding to the α-th blending function. The control points Pαare found by imposing the interpolation conditions xβ=S(ξβ)=∑ α∈AT PαRα(ξβ),∀ξβ, β ∈AT(3.2) The resulting T-spline parameterization is shown in Fig. 3.5(b). On the other hand, for a comparison we perform the interpolation of the same data points using polynomial spline blending functions, which is the main result of the thesis. Our blending functions {Nα}α∈ATdefined over the parametric T-mesh span a complete polynomial space over T, i.e., P3(Ω) ⊂span{Nα}α∈AT.The resulting spline parameterization is shown in Fig. 3.5(c). It can be appreciated that our polynomial parameterization leads to a more uniform parametric mapping, meanwhile rational T-spline parameterization presents some oscillation, see Fig. 3.5(d). Remark 4. Also, it worth mentioning that recently the issue of using rational spline functions, was discussed in [55] (It is Time to Drop the “R” from NURBS). The authors discuss the elevated complexity and computational cost of NURBS and that using of rational forms for parameterization also entails some disadvantages, regarding the quality. The use of NURBS in recent decades was motivated by their capacity of representing exactly conic sections. However, the authors suggest that there is no need for the rational form, since the conics can be approximated with sufficient accuracy in a robust and efficient manner with polynomial B-splines. 34 3 Local refinement problem (a) Parametric T-mesh. (b) T-spline parameterization. (c) Polynomial spline parameterization. (d) Some oscillations of rational T-spline parametrization. Figure 3.5: Example of using rational T-splines for domain parameterization. 35 3 Local refinement problem 3.4 Analysis-suitable T-splines T-splines was an important invent that overcame the drawback of tensor product NURBS spaces and provided a flexible tool for geometric modelling. However, for their use in analysis, T-splines lack some essential properties. As a solution to this problem, analysis-suitable T-splines were proposed in [23], which is a class of T-splines defined over T-meshes with certain topological restrictions formulated in terms of T-junction extensions. T-junction extensions include face and edge extensions, which are line segments originating at a T-junction, see Fig. 3.6. For cubic T-splines a face extension is created by marching from the T-junction, in the direction of a missing edge, until two perpendicular edges are intersected. An edge extension is then formed by marching in the opposite direction until intersect one edge. A T-mesh is called analysis-suitable if no T-junction extension intersects perpendicularly another T-junction extension. (a) (b) (c) (d) Figure 3.6: Analysis-suitable T-spline. (a) Initial, not analysis-suitable T-mesh. (b) Possible extensions to obtain analysis-suitable T-mesh. (c) T-junction extensions of the initial mesh, edge extension in red, face extension in green. (d) T-junction extensions of the final mesh do not intersect. T-spline functions defined over an analysis-suitable T-mesh are linearly independent [24] and possess polynomial reproduction property. The refinement algorithm allows to accomplish highly localized refinements and construct nested polynomial spline spaces with optimal approximation properties. However, to guarantee all these properties, it is necessary to work with T-meshes that fulfil the 36 3 Local refinement problem required topological restrictions. A given T-mesh should be extended by adding new edges until obtaining an analysis-suitable T-mesh, which, in general, is not trivial task that involves a global search of optimal T-mesh extension. Besides, the solution is not unique. The strategy presents an elevated implementation complexity and, as far as we know, the extension of the strategy to 3D cases is still an open question. 3.5 Hierarchical refinement scheme Another well known approach for local refinement is the hierarchical refinement scheme. This approach is based on a simple and natural idea to construct multilevel spaces by replacing coarse level basis functions by the finer functions. Starting from an initial uniform mesh, hierarchical refinement scheme leads to sequential construction of nested spline spaces with linearly independent basis functions. This approach is based on the refinability property of B-spline functions. Let Ξ0be the initial uniform knot vector. We can perform kuniform subdivisions to obtain a k-level knot vector Ξk. Then, a basis function Bk i,p defined on Ξkcan be expressed as a linear combination of p+ 2 basis functions defined on Ξk+1, i.e., Bk i,p = p+1 ∑ j=0 cjBk+1 j,p ,where cj= 2−p(p+1 j). The p+2 basis functions Bk+1 j,p of level k+1 are called children of the k-level function Bk i,p. For example, Fig. 3.7 shows a cubic B-spline function Bk 0,3defined on the knot vector Ξk={1,2,3,4,5}. An uniform subdivision of Ξkleads to the knot vector Ξk+1 ={1,1.5,2,2.5,3,3.5,4,4.5,5}and 5 children functions Bk+1 j,3(j= 0, . . . , 4) defined on it. Then, the function Bk 0,3is expressed as linear combination of its children Bk 0,3=1 8Bk+1 0,3+1 2Bk+1 1,3+3 4Bk+1 2,3+1 2Bk+1 3,3+1 8Bk+1 4,3. Figure 3.8 illustrates the idea of hierarchical refinement scheme. Figure 3.8(a) shows the initial (level 0) B-spline space, where the interval [4,9] is to be refined. The basis functions (blue discontinuous line), whose support is contained in the refined interval [4,9], are replaced by the finer basis functions of level 1 shown in Fig. 3.8(b). Next, the interval [5,8] is refined and all basis functions of level 1, whose support is contained in the refined area, are replaced by level 2 functions shown in Fig. 3.8(c). The final spline space, shown in Fig. 3.8(d), is composed by the basis functions from three figures represented in black continuous line. Note that all the functions are a scaled and translated version of each other. The idea of hierarchical refinement scheme for 2D is illustrated in Fig. 3.9. Relatively simple implementation and straightforward generalization to 3D make it an attractive option for local refinement. However, there are some drawbacks. Hierarchical refinement scheme permits to enrich a given spline space by replacing some functions by new functions of finer level, but it is not clear how to 37 3 Local refinement problem 1 2 3 4 5 123451 2 3 4 5 Figure 3.7: Refinability of B-spline functions. Rightmost figure shows the function Bk 0,3 and its children Bk+1 j,3scaled by their corresponding coefficients. define a spline space for a given T-mesh if no previous mesh (or space) is given. In other words, unlike T-splines, this strategy does not define a spline space from a given arbitrary mesh, but perform a sequential enrichment of spline spaces by adding new functions and eliminating the old ones to avoid linear dependency. In addition, hierarchical refinement scheme imposes a certain obligatory requirements on the refined area that extend the refined zone more than necessary. For example, for 2D mesh, in the worst case, refinement of only one cell leads to the refinement of totally 25 cells and introduction of 49 new basis functions. In 3D case the same situation leads to 125 new cells. And finally, another shortcoming of the strategy is an excessive support overlapping, which can affect the sparsity and condition number of the stiffness matrix. To solve this issue, a truncation technique was proposed in [30]. Truncated hierarchical splines are derived from the hierarchical basis by redefining the functions, namely, eliminating the contributions of the functions of finer levels to coarser level ones, and thus improving their locality. However, this complicates the strategy. 3.6 Motivation of our work We need to construct spline spaces with nice properties for their use in analysis and design. The main goal we pursue here is to find a simple practical solution to the problem of local refinement with straightforward implementation, both in 2D and 3D. Inspired in T-splines and hierarchical refinement scheme, the aim of our work was to elaborate another possible strategy to define spline spaces that combine only nice features of each strategy. The rest of the paper exposes the developed technique. 38 3 Local refinement problem (a) (b) 1 2 3 4 5 6 7 8 9 10 11 12 (c) 1 2 3 4 5 6 7 8 9 10 11 12 (d) Final space composed by the functions of three levels plotted in black. Figure 3.8: Hierarchical refinement scheme in 1D. 39 3 Local refinement problem level 0 level 1 level 0 level 1 Figure 3.9: Hierarchical refinement scheme in 2D. 40 CHAPTER 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes In this chapter we expose the main result of the thesis, that is, a new strategy to define tensor product spline functions over T-meshes. Due to the elevated complexity of all current strategies, the main goal we pursue here is the simplicity and low computational cost, both in 2D and 3D. For this, we have to assume a restriction on the T-mesh. Namely, the T-mesh should fulfil the requirement of being a0-balanced (strongly balanced) quadtree/octree mesh. Assuming this reasonable restriction over the T-mesh, we can define easily cubic spline functions that span spaces with desirable properties: linear independence, C2-continuous, cubic polynomial reproduction property, nestedness of spanned spaces and a straightforward implementation. The key of the strategy lies in some simple rules used for inferring local knot vectors for each blending function. 4.1 Main steps of the strategy. The strategy we propose here has some similarity with T-splines inasmuch as we define the blending functions from local knot vectors that are inferred by traversing T-mesh edges. Some additional rules and requirements are imposed for the local knot vectors in order to obtain spline spaces with desired properties. These additional rules were elaborated by studying and analysing different situations, when T-spline functions do not span a complete polynomial space. At the same time some features of the resulting rules and requirements are inspired by the hierarchical refinement scheme. The strategy is motivated by the idea of preserving nesting behaviour of the spaces under the mesh refinement. We have to prevent the situation, when a blending function cannot be reproduced with the basis of the new spline space after the refinement. This would guarantee the nestedness of the spline spaces. For example, in the initial mesh of the Fig. 4.1(a), we define a T-spline basis function Nαassociated to the vertex α. Then, we perform some cell refinements as shown in Fig. 4.1(b). It is easy to check that the new 41 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes Figure 4.6: Support notation for 2D case. modification for a bivariate function. Leftmost function support shown in Fig. 4.7(a) does not meet Condition 1. For the knot vector Ξ we have ∆ξ 3> ∆ξ 4, so the knot interval ∆ξ 4should be modified. Let us denote h= max(∆ξ 2, ∆ξ 3) = max(∆η 2, ∆η 3). Note that both maxima coincide due to the quadtree structure and the fact that the T-junctions are skipped. Then, to modify ∆ξ 4, we just double its size by redefining the fifth knot ξ5as ξ∗ 5←ξ3+ 2h. For the local knot vector H we have ∆η 2> ∆η 3, so the knot intervals ∆η 3and ∆η 4should be modified. Again, we just double their size. Namely, the knots η4and η5are redefined as η∗ 4←η3+h, η∗ 5←η3+ 2h. Note that basically we skip over some knots of the mesh in order to double the size of some knot intervals. Extension rule 2. Once Condition 1 is satisfied, in order to fulfil Condition 2, we check whether the vertices of the function support are situated over the mesh edges. If not, we modify the knot vectors by skipping over a knot for both parametric directions and placing this vertex over the mesh edges. Note that, again it implies to double the size of some knot intervals. In this case only exterior knot intervals are extended. An example of a function support violating Condition 2 is illustrated in Fig. 4.7(b). The corner vertex V5,5= (ξ5, η5) of this support is not situated over a mesh edge, so the intervals ∆ξ 4and ∆η 4are to be extended. For this the fifth knots for both parametric directions are redefined as ξ∗ 5←ξ3+ 3h, η∗ 5←η3+ 3h, and thus, the new vertex V5,5is placed over the mesh edges. The checking of Condition 2 and the extension is performed independently for each of the four quadrants of the function support. Note that for our 0-balanced quadtree we should make this checking only for some functions. For example, without loss of generality, the support vertex V5,5= (ξ5, η5) must be checked only if ∆ξ 3=∆ξ 4= ∆η 3=∆η 4. Remark 5. Note that, when a function support violates Condition 2, the extension of this support, in order to place its vertex over the mesh skeleton, can be made 48 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes (a) Condition 1 is not satisfied because ∆ξ 4<∆ξ 3and ∆η 3<∆η 2. Knot intervals to be changed are marked in blue. (b) Condition 2 is not satisfied because V5,5/∈skt(T). Figure 4.7: Extension rules. (a) An example of support modification by Extension rule 1. (b) An example of support modification by Extension rule 2. Figure 4.8: Extension rule 2. Other possibilities to place the corner vertex over the mesh skeleton. by modifying only one of the knot vectors instead of the both, see Fig. 4.8. This option leads to the loss of symmetry for the spline space. That is, a symmetric T-mesh would have a non symmetric function supports. For simplicity we extend the knot vectors in both parametric directions and conserve the symmetry. 49 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes Algorithm 1: Extension rule 1. Input: A knot vector Ξ = (ξ1, ξ2, ξ3, ξ4, ξ5), ξi∈[0,1]. 1Function Modify1(Ξ) 2Ξ∗←Ξ 3h= max(∆2, ∆3) 4if ∆2< ∆3and ξ2>0then 5ξ∗ 2←ξ3−h 6ξ∗ 1←ξ∗ 2 7if ξ∗ 1>0then ξ∗ 1←ξ3−2h 8 9if ∆1<∆2and ξ1>0then 10 ξ∗ 1←ξ3−2h 11 if ∆2>∆3and ξ4<1then 12 ξ∗ 4←ξ3+h 13 ξ∗ 5←ξ∗ 4 14 if ξ∗ 5<1then ξ∗ 5←ξ3+ 2h 15 16 if ∆4<∆3and ξ5<1then 17 ξ∗ 5←ξ3+ 2h 18 return Ξ∗ Output: A corrected knot vector Ξ∗that satisfies Condition 1. Algorithm 2: Extension rule 2 in 2D. Input: A 0-balanced mesh Tand a pair of local knot vectors S={Ξ,H}. 1Function Modify2(T, S) 2S∗←S 3h= max(∆ξ 2, ∆ξ 3) 4for n∈ {1,5}do 5for m∈ {1,5}do 6if (ξn, ηm)/∈skt(T)then 7ξ∗ n←ξ3+ 3hsgn(ξn−ξ3) 8η∗ m←η3+ 3hsgn(ηm−η3) 9return S∗ Output: A modified support S∗={Ξ∗,H∗}that satisfies Condition 2. The extension of any other function support is completely analogous to these two examples. In all possible cases, the extension of a function support implies to change some knot intervals by doubling its size. Detailed algorithms for Extension rule 1 and 2, used to modify a bivariate 50 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes function support according to Conditions 1and 2, are given in Algorithms 1 and 2. Figure 4.9 shows some examples of support extension procedure. Functions that satisfy both conditions and do not need be modified are given in Fig. 4.9(a). Examples of support extension according to Condition 1 are shown in Fig. 4.9(b), (c), (d) and (e). And Fig. 4.9(f), (g) and (h) illustrate support extension according to Condition 2 or both. Note that an exception for Condition 1 is a knot vector that contains a knot interval of length 0 due to the open knot vector structure along the boundary. In this case, a knot vector should fulfil the inequality (4.1) not taking into account the knot intervals of length 0. That is, Extension rule 1 is applied only to non-zero knot intervals, see Fig. 4.9(c), (d) and (e). Remark 6. It is important to highlight that application of the extension rules always place the redefined knots over the mesh edges, i.e., the extension rules just skip over some knots of the mesh, but do not invent new knots that are not induced by the T-mesh. Remark 7. It is important to underline that inferring of each function support does not depend on the rest of the functions, so the process can be parallelized. Remark 8. As was said at the beginning of the chapter, it is important take into account that the rules and algorithms we formulate here are designed for uniform quadtree subdivisions, where all cells of the same level have the same size. However, if for some reason, it is necessary to deal with non-uniform subdivision, the process of inferring of functions supports for such mesh can be translated to another, isomorphically equivalent, mesh (let us call it reference mesh), where all subdivisions are strictly uniform. 51 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes (a) Initial supports that satisfy both conditions and should not be modified. (b) Condition 1 is not satisfied because ∆ξ 1<∆ξ 2, so ξ∗ 1←ξ3−2h. (c) Condition 1 is not satisfied because ∆ξ 2<∆ξ 3, so ξ∗ 1←ξ3−2hand ξ∗ 2←ξ3−h. (d) Condition 1 is not satisfied because ∆η 1<∆η 2, so η∗ 1←η3−2h. (e) Condition 1 is not satisfied because ∆ξ 2<∆ξ 3and ∆η 2<∆η 3, so ξ∗ 2←ξ3−h, η∗ 2←η3−hand η∗ 1←η3−2h. (f) Condition 2 is not satisfied because V5,5/∈skt(T), so ξ∗ 5←ξ3+ 3hand η∗ 5←η3+ 3h. (g) Condition 2 is not satisfied because V1,1and V1,5/∈skt(T), so ξ∗ 1←ξ3−3h, η∗ 1←η3−3hand η∗ 5←η3+ 3h. (h) Conditions 1 and 2 are not satisfied because ∆ξ 2>∆ξ 3, ∆η 2>∆η 3 and V1,1/∈skt(T), so ξ∗ 4←ξ3+h, ξ∗ 5←ξ3+ 2h,η∗ 4←η3+h,η∗ 5←η3+ 2h, ξ∗ 1←ξ3−3h,η∗ 1←η3−3h. Figure 4.9: Examples of function support modification with Extension rule 1 and 2. Initial support is marked in blue and extended support in red. 52 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes Figure 4.10: Support notation for 3D case. 4.5.2 Support extension for 3D meshes In this section we give a description and illustration of the proposed strategy for defining trivariate spline functions over 0-balanced octree T-meshes. The skeleton skt(T) of a three-dimensional mesh Tis the union of all faces of the mesh. For a trivariate function let us denote the vertices of its support as Vn,m,k = (ξn, ηm, ζk), where n, m, k ∈ {1,5}. And the edge formed by two vertices Vn,m,k and Vp,q,r is denoted by E(n,m,k),(p,q,r). Then, the frame frm(supp Nα) of a trivariate function support is the union of its twelve edges. Figure 4.10 illustrates the introduced notation for the support of a trivariate blending function. The formulation of Condition 1 for the local knot vectors of a trivariate function is the same as for the 2D case. Condition 2 adapted to 3D meshes is stated as follows: Edges of the cuboidal function support should be situated over the mesh faces. The implementation of the strategy for 3D is similar to the 2D case. To satisfy Condition 1, Extension rule 1 is applied to each of the three local knot vectors of a function, analogously to the 2D case, using Algorithm 1. Extension rule 2. In order to fulfil Condition 2, we check whether the edges of a function support are situated over the mesh faces. If not, the two knot vectors perpendicular to this edge should be modified by skipping over a knot for both parametric directions and placing this edge over the mesh faces. Analogously to 2D, we double the size of some exterior knot intervals. The checking is performed independently for each of the eight quadrants of the function support and, in each quadrant, three edges should be checked. Figure 4.11 illustrates the support extension procedure for the quadrant of the vertex V5,1,1. Due to the octree structure only two cases can take place: (i) the quadrant contains one edge that is not situated over the mesh faces or (ii) the quadrant contains three edges and a vertex that are not situated over the mesh faces. Let see each case. (i) If a quadrant contains one edge that does not fulfil Condition 2, then two 53 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes (a) Only one edge (in blue) violating Condition 2. Node V5,1,1is sited in the center of the face of size 2h. The support is extended in two directions. (b) Three edges (in blue) violating Condition 2. Node V5,1,1is sited in the center of the cell of size 2h. The support is extended in three directions. Figure 4.11: Extension rule 2 for support modification of a trivariate function. knot vectors perpendicular to this edge are modified, see Fig. 4.11(a). For the function support shown in Fig. 4.11(a) left, the edge E(5,1,1),(5,1,5) is not situated over the mesh faces. Therefore, two knot vectors Ξ and H, perpendicular to this edge, are modified in order to place the edge over the mesh faces, namely, the knots ξ5and η1are redefined as ξ∗ 5←ξ3+ 3hand η∗ 1←η3−3h, where h= max(∆ξ 2, ∆ξ 3) = max(∆η 2, ∆η 3) = max(∆ζ 2, ∆ζ 3). (ii) If a quadrant contains three edges that are not situated over the mesh faces, then the three knot vectors are modified by skipping over a knot for each of the three parametric directions, see Fig. 4.11(b). Vertex V5,1,1and the three edges connected to it are not situated over the mesh skeleton, so all the three knot vectors are extended to place the three edges over the mesh faces: ξ∗ 5←ξ3+ 3h, η∗ 1←η3−3h, ζ∗ 1←ζ3−3h. 54 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes Algorithm 3: Extension rule 2 in 3D. Input: A 0-balanced T-mesh Tand three local knot vectors S={Ξ,H,Z}. 1Function Modify2 (T, S) 2S∗←S 3h= max(∆ξ 2, ∆ξ 3) 4mov(i) := i+ 4 sgn(3 −i) 5for n∈ {1,5}do 6for m∈ {1,5}do 7for k∈ {1,5}do 8if E(n,m,k),(mov(n),m,k)/∈skt(T)then 9η∗ m←η3+ 3hsgn(ηm−η3) 10 ζ∗ k←ζ3+ 3hsgn(ζk−ζ3) 11 if E(n,m,k),(n,mov(m),k)/∈skt(T)then 12 ξ∗ n←ξ3+ 3hsgn(ξn−ξ3) 13 ζ∗ k←ζ3+ 3hsgn(ζk−ζ3) 14 if E(n,m,k),(n,m,mov(k)) /∈skt(T)then 15 ξ∗ n←ξ3+ 3hsgn(ξn−ξ3) 16 η∗ m←η3+ 3hsgn(ηm−η3) 17 return S∗ Output: A modified support S∗={Ξ∗,H∗,Z∗}that satisfies Condition 2. Algorithm 3 explains the Extension rule 2 used to modify a trivariate function support according to Condition 2. 4.6 Classification of the functions supports (2D case) It is worth noting that the proposed rules and the 0-balanced restriction for the mesh lead to a very few number of all possible function supports that can be defined over a T-mesh. Let us consider any knot vector Ξ = (ξ1, ξ2, ξ3, ξ4, ξ5) and its knot intervals ∆ξ i=ξi+1 −ξi,i= 1, ..., 4.The set ∆ξ= (∆ξ 1, ∆ξ 2, ∆ξ 3, ∆ξ 4) determines a knot vector Ξ. Due to the Extension rule 1 and the balanced mesh condition, we have the following possible configurations for the set ∆ξ: (h, h, h, h), (h, h, h, 2h), (2h, h, h, h) and (2h, h, h, 2h), see Fig. 4.12(a). All possible bivariate function supports are obtained by combining all configurations for Ξ and H, i.e., combining two sets ∆ξ and ∆η. Note that, since we ignore the hanging nodes while inferring two knot vectors, the value of his the same for ∆ξand ∆η, i.e., h=∆ξ 2=∆ξ 3=∆η 2=∆η 3, see Fig. 4.12(b), and the value of hdetermines the level of a function. We say that 55 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes (a) (b) ... ... ... ... m m+1 m-1 m-2 (c) Figure 4.12: (a) All possible configurations for a knot vector (without repeated knots). (b) For any function support h=∆ξ 2=∆ξ 3=∆η 2=∆η 3.(c) A function support of level m, which is determined by the cell marked in red. a blending function is of level mif the largest cell contacting its anchor is of level m, see Fig. 4.12(b)(c). Thus, the functions we define according our strategy have a hierarchy. We can classify all function supports that can take place as follows. First, we can distinguish different function supports according to the underlying local mesh that leads to this support. Basically, it depends on the position of the anchor of the function, see Fig. 4.13. For scenario (a) the anchor of the m-level function coincides with the center of a (m-1)-level cell of the quadtree; scenario (b) the rest of the functions, i.e., the anchor of the m-level function does not correspond to any center of a (m-1)-level cell of the quadtree. Second, taking into account all the considerations about values of ∆ξand ∆η, we can classify all function supports in nine types that are illustrated in Fig. 4.13. First six types correspond to scenario (a) and the other three types correspond to scenario (b). Note that the supports in Fig. 4.13 are illustrated over the simplest (minimal) possible mesh configuration that gives rise to this support. However, it is important to notice that the same function support can take place over another (more refined) underlying mesh configuration, see for example Fig. 4.12(b)(c). 56 4 Strategy for construction of polynomial spline spaces over hierarchical T-meshes (a) (b) (c) Type 1 (d) Type 2 (e) Type 3 (f) Type 4 (g) Type 5 (h) Type 6 (i) Type 7 (j) Type 8 (k) Type 9 Figure 4.13: Two possible scenarios for a function support according to the underlying local mesh. (a) The anchor of the m-level function coincides with the center of a (m1)-level cell of the quadtree. (b) The rest of the functions. (c)-(k) Function support classification according to the underlying local mesh and configurations for ∆ξand ∆η. Remark 9. For convenience, we are going to put a name to the spline blending functions defined over T-mesh according to the proposed strategy. From now on, we are going to refer to the functions as Extended polynomial splines (EPsplines). 57 5 Properties of the constructed spline spaces 5.3 Linear independence For now we cannot provide a rigorous proof of this property. It is currently under study. However, for a better understanding of the strategy, some clues should be mentioned. For now, our strategy is designed for T-meshes without repeated interior knots. To prove the linear independence of our blending functions it is essential to take into account that each regular node of the T-mesh has only one function assigned to it. Some reasoning from [54] can be used to proof the linear independence of such functions. In particular, the authors show that T-splines without multiple interior knots can not be linearly dependent. Note that the example of linearly dependent T-splines given in [53] is for a T-mesh with multiple interior knots. The main reason of the linear independence of our spline functions lies in the fact that all the functions have different anchors (the third knot of its local knot vectors), i.e., there exist no pair of functions that have the same anchor. Thus, the proof of linear independence of our functions boils down to prove that a set of spline functions can be linearly dependent only if there are some functions with the same anchor. This hypothesis seems reasonable taking into account the knot insertion procedure, see Section 2.1.2.1 and Fig. 2.5. According to this formula, a spline function B0is split in two functions B0=αB1+βB2under a new knot insertion, where one of the functions has the same anchor as the original function B0. So, for the linearly dependent set {B0, B1, B2}there are two functions with the same anchor. It should be mentioned that the linear independence of our spline functions was verified in numerous computational experiments due to the non-singularity of the interpolation and stiffness matrices. 5.4 Non-negative weighted partition of unity Partition of unity for spline functions provides some useful properties as affine invariance and convex hull property for spline curves and surfaces, which are very important for geometric design. Taking into account the nesting behaviour of our spline spaces, the constant function 1 belongs to our spline space, i.e., 1 = ∑ α∈AT wαNα. Hence we can define weighted functions N′ α=wαNα.This new basis {N′ α}α∈AT forms a partition of unity: 1 = ∑ α∈AT N′ α. 64 5 Properties of the constructed spline spaces (a) (b) (c) Figure 5.5: Verification of the non-negativity for the coefficients of linear combination for 0-level function N0; (a) uniform mesh T0and its uniform support N0. (b)(c) some of the 216 possible mesh refinements of the support N0and the corresponding functions of T1necessary to represent N0. However, there are some issues that should be discussed. (i) Firstly, it is important that the weights are non-negative. In this way the basis {N′ α}α∈ATconserves the property of non-negativity. To verify this fact we reason as follows. Suppose that we start with the uniform mesh T0of level 0 whose functions form a strictly positive partition of unity: 1 = ∑ α∈AT0 N0 α,see Fig. 5.5(a). Then, let us consider a refined mesh T1of level 1. We can show that any function N0of the mesh T0is a linear combination of the functions from T1with non-negative coefficients, i.e., N0=∑ i∈AT1 ciN1 i,where ci≥0. For that we have to analyse all possible mesh configuration within the support of the function N0and verify for each of them that all coefficients of the linear combination to recover N0are non-negative. The support of the function N0 comprises 16 cells of level 0 that can be refined or not on the mesh T1,giving rise to different linear combinations to recuperate the function N0.Figure 5.5(b) and (c) illustrates some examples of such configurations and their functions necessary to recuperate the function N0. In total, there are 216 different mesh configurations to study. To do this, we perform a computational experiment. We generate all possible 216 mesh configurations and verify for each mesh T1that the function N0 is a linear combination of the functions of T1with non-negative coefficients. For that we interpolate the function N0using the approximation space of the refined mesh T1. This interpolation gains to reproduce exactly the function N0and all coefficients are non-negative. Having done this verification, we deduce that the 65 5 Properties of the constructed spline spaces (a) (b) Figure 5.6: Verification of the non-negativity for the coefficients of linear combination for a function N1of level 1. functions of the mesh T1form a non-negative partition of unity since 1 = ∑ α∈AT0 N0 α=∑ α∈AT0  ∑ i∈AT1 cα iN1 i =∑ i∈AT1 wiN1 i, wi≥0. Next, for a T-mesh T2of level 2 we repeat the reasoning and show that any function N1of the mesh T1is a linear combination of the functions of T2with non-negative coefficients, see Fig. 5.6. Indeed, due to 0-balanced condition, a function support N1is replaced by the new functions of level 2, Fig. 5.6(b), only if the support of N1is composed of the 16 cells of level 1, i.e., it is a scaled version of the support N0.So the entire procedure, described before for N0, can be repeated to show that the linear combination to represent the function N1includes only non-negative coefficients. Thus, assuming the non-negative partition for the mesh T1, we can deduce the non-negative partition for the mesh T2: 1 = ∑ i∈AT1 wiN1 i=∑ i∈AT1  wi∑ j∈AT2 ci jN2 j =∑ j∈AT2 w′ jN2 j, w′ j≥0. Repeating the reasoning we can show the non-negative partition of unity for a mesh of any level k. (ii) Second issue to discuss. There are mesh configurations for which some weights wican be zero. For example, Fig. 5.7(a) shows an initial mesh T0whose basis forms a strictly positive partition of unity, since it is uniform bivariate Bsplines, i.e., 1 = ∑ α∈AT0 Nα.After some refinements around the anchor of the function N0, shown in Fig. 5.7(b), the new space corresponding to T1consists of the same basis function of the mesh T0(due to Extension rule 1) plus 5 new functions of the finer level. These 5 new functions are not necessary for a partition of unity. Hence, they should be annulled in order to construct a strictly positive partition of unity. However, if the four cells that share the anchor of N0are refined, then the 66 5 Properties of the constructed spline spaces (a) (b) (c) Figure 5.7: (a) A mesh with partition of unity. (b) Refined mesh with some redundant functions. (c) Refined mesh without redundant functions and thus with strictly positive partition of unity. function N0is replaced by the finer ones, see Fig. 5.7(c). In this case, there are no redundant functions and the basis of the mesh T2form a strictly positive partition of unity. Also, a T-mesh can have some redundant functions due to Extension rule 2. Figure 5.8 shows such examples. The functions of the mesh of Fig. 5.8(a) form a strictly positive partition of unity. Then, it can be checked that the refined mesh of Fig. 5.8(b) has some new functions that are redundant for partition of unity. The anchors of active functions are marked by black circles, i.e., only these functions are necessary to form a partition of unity. (a) (b) (c) Figure 5.8: (a) Initial mesh with strictly positive partition of unity. (b) Refined mesh with redundant functions due to Extension rule 2. Active functions are marked with black circles, redundant functions with blue circles. (c) Mesh without redundant functions. Remark 10. Similar problem also happens in hierarchical refinement scheme. Therefore, an additional requirement is imposed on the refinement pattern, namely the refined zone should at least occupy the area of a function support of the current level in order to produce a replacement of this support by the finer ones. In this way, none of the new functions will be redundant for a partition of unity. 67 5 Properties of the constructed spline spaces It should be mentioned that in our computational examples presented in this manuscript we do not apply any additional treatment for our spline spaces, so the basis used for the analysis does not form a partition of unity. 5.5 Locality of the functions and possible excessive support overlapping It should be mentioned that the proposed strategy can lead, in some cases, to increased overlapping of the function supports of different refinement levels, which can affect the conditioning and sparsity of the stiffness matrix. Excessive accumulation of the functions can take place in some problems with a very sharp singularity, when the area marked to refine by the error indicator at each iteration is too small. According to our strategy, a cell refinement always adds at least one new blending function of the finer level, while some functions of the coarse level can stay unchanged after a refinement. Note that a coarse level function support is reduced only if the four cells that contact its anchor are refined, see Fig. 5.7(c). Thus, if the adaptive refinement produces a sequence of small nested refined areas, some functions of the coarse levels stay active (not null) in the refined zone, see Fig. 5.9(a). To avoid a possible accumulation of the functions of different levels, the following strategy can be adopted: if a function support intersects the cell marked to refine, then the coarsest cells sharing the anchor of the function are also refined. Application of this approach is illustrated in Fig. 5.9. In this example, let us suppose that the problem has a singularity at the center of the domain, and at each iteration, we refine only the four cells adjacent to the singularity point. Then at each iteration there are some functions whose support stays unchanged. However, if at each iteration we perform additional refinements, functions of the coarse level are replaced by the finer ones. Application of this additional refinement, where needed, can reduce an excessive function overlapping and improve the stiffness matrix conditioning and sparsity. However, due to the unnecessary extension of the refined zone, the optimal rate of convergence can be lost. In practice, we did not observe a significant advantage of this approach, and taking into account its computational cost, we do not apply this strategy in the computational examples presented in this work. Moreover, it was observed that for a large variety of problems, the excessive support overlapping is avoided naturally after various refinement steps due to the behaviour of the error indicator. Only for one computational example with a very strong singularity (section 7.2.8) the situation described in this paragraph takes place and additional refinements were applied. On the other hand, another possible solution for this problem is a more accurate and selective definition of the functions in order to obtain spaces with a better function supports locality. We believe that the strategy we propose in this work can be improved by including a little more sophisticated rules for inferring local 68 5 Properties of the constructed spline spaces (a) Three refinement steps around a singular point (black square) that lead to an accumulation of the function supports. Note that the red support is not changed during refinements. (b) The same three refinement steps applying additional refinements. Coarse level support is reduced. Figure 5.9: Additional refinements in order to avoid excessive overlapping of the function supports. knot vectors to avoid an unnecessary support extensions in some cases. Remark 11. It should be mentioned that the similar problem happens for hierarchical refinement scheme. To avoid this problem, analogously to our case, some additional refinements can be imposed, and these refinements are more extensive than in our case. As another solution to this issue, a truncation technique was proposed in [30]. Truncated hierarchical splines are derived from the hierarchical basis by redefining the functions, namely, eliminating the contributions of functions of finer levels to coarser level ones, and thus improving their locality. Truncated basis has smaller supports, non-negativity, partition of unity and local linear independence. Besides, the construction of spline spaces using truncated hierarchical B-splines is strongly stable versus weak stability of classical hierarchical refinement scheme, see [57]. Also, it is worth mentioning that a possible supports accumulation in a cell is inherent to the tensor product structure of the basis functions. Note that support accumulation also can take place with T-splines. Namely, it happens for the mesh refinements shown in Figure 5.9(a). However, for T-splines the support overlapping is considerably minor and it happens only for this particular situation. 69 Appendix 5.A Nesting behaviour of the spaces In this section we give a proof of the nestedness of the spline spaces constructed by means of our technique (2D case). 5.A.1 Step 1 of the proof First, we perform some reasoning in order to show that the proof can be reduced to a study of a few number of some basis cases. Lemma 1. For any two 0-balanced nested meshes T1⊂T2there exists a sequence of meshes {Te i}n i=0 such that (i) T1=Te 0⊂Te 1⊂Te 2⊂ ··· ⊂ Te n−1⊂Te n=T2. (ii) any mesh Te iis obtained from the mesh Te i−1by refining only one cell. (iii) all the meshes {Te i}n i=0 are 0-balanced. Proof. With abuse of notation, for a given T-mesh T, we denote its corresponding quadtree also as T. The term tree node is used exclusively in the context of the tree structures terminology, see section 4.2. Let us consider a 0-balanced quadtree mesh T1. And let T2be a 0-balanced refinement of T1. We denote k1= lev(T1) and k2= lev(T2), then k1⩽k2. The set of all internal nodes of the quadtree T2that do not belong to T1is denoted by N={Ni}n i=1. Elements of the set Ncan be sorted by increasing order of their levels: N:= {N1,N2, . . . , Nn}, where k1⩽lev(Ni)⩽lev(Ni+1)⩽k2for any i= 1, ..., n −1. Then, we start to refine the mesh T1, adding to the corresponding quadtree one by one the elements of the sorted set N. The addition of a new internal tree node corresponds to the refinement of one cell of the previous mesh. Note that internal node is a tree node that has children. So, the addition of the internal node Niof level mmeans that a m-level leaf node (cell) turns into the internal node after being subdivided. Thus, we construct a sequence of refinements that leads to the final mesh: T1=Te 0⊂Te 1⊂Te 2⊂ ··· ⊂ Te n−1⊂Te n=T2,see Fig. 5.A.1. 71 5 Properties of the constructed spline spaces Figure 5.A.1: Elemental refinements. A possible sequence of 0-balanced meshes {Te i}, where Te iis obtained from Te i−1by refining only one cell. Now, we have to assure that each of these meshes is 0-balanced. Indeed, let us suppose that the mesh Te iis the first one of the sequence that is not 0-balanced, i.e., the mesh Te i−1is 0-balanced and Te iis not. If Te iis obtained from Te i−1by refining a cell of level m, then, in order to obtain the final 0-balanced mesh T2, at least one neighbouring cell of level (m-1) must be refined. By construction (the nodes are sorted in increasing order), this does not occur, i.e., no cell of level (m-1) is refined later. Therefore, the rest of the meshes of the sequence, including T2, are also 0-unbalanced. This conclusion contradicts the assumption that the mesh T2 is 0-balanced. Definition 3. Refinements that fulfil the conditions of the Lemma 1 are called elemental refinement, i.e., a pair of 0-balanced meshes Te 1⊂Te 2is called elemental refinement if Te 2is obtained from Te 1by refining only one cell. Lemma 2. If any elemental refinement produces nested spaces: Te 1⊂Te 2⇒ STe 1⊂STe 2,then any two 0-balanced nested meshes produce nested spaces: T1⊂T2⇒ST1⊂ST2. Proof. For the sequence of elemental refinements constructed in Lemma 1 we have: ST1=STe 0⊂STe 1⊂STe 2⊂. . . ⊂STe n=ST2,and therefore ST1⊂ST2. Now, taking into account Lemma 2, our goal is to prove that any elemental refinement produces nested spaces. To proof this fact, it is crucial to take into account several considerations about our strategy. As was discussed in section 4.6, there are only few types of possible function supports. Each function is determined by its knot intervals ∆ξ= (∆ξ 1, ∆ξ 2, ∆ξ 3, ∆ξ 4) and ∆η= (∆η 1, ∆η 2, ∆η 3, ∆η 4).Due to the extension rules and the balanced mesh condition, we have the following possible 72 5 Properties of the constructed spline spaces (a) One-ring neighbourhood of the cell τ. (b) (m-1)-level function of the ring turns into m-level function. (c) m-level function of the ring turns into another m-level function. (d) m-level function, with the anchor out of the ring, is not affected. (e) (m-1)-level function, with the anchor out of the ring, is not affected. Figure 5.A.2: An example of elemental refinement Te 1⊂Te 2and the one-ring neighbourhood of the cell τ. (a) The contour of the one-ring neighbourhood and the anchors of its functions are marked in blue. Anchors, whose functions are changed after the cell refinement, are marked in black in the rightmost figure. (b) and (c) The functions of the one-ring neighbourhood that are changed after the refinement. (d) and (c) Functions whose anchors are out of the ring always stay unchanged. configurations for the set ∆: (h, h, h, h), (h, h, h, 2h), (2h, h, h, h) and (2h, h, h, 2h). It is important to highlight that the same function support can take place for different underlying meshes, since we have to skip some knots of the mesh in order to fulfil Condition 1 and Condition 2, see Section 4.4. Consequently, refinement of a cell does not affect some functions whose support is determined by the cells of a coarser level. Having said that, refinement of a cell affects only some of the function supports in the neighbourhood the cell. Namely, the anchors of the affected function supports belong to the one-ring neighbourhood of the cell. We are going to call one-ring neighbourhood of the cell τthe set of all cells that have contact with this 73 5 Properties of the constructed spline spaces scenario 2. The knot vector Ξ, associated to the anchor ξ3becomes larger after the refinement due to Extension rule 2. Namely, the vector Ξ = (ξ1, ξ2, ξ3, ξ4, ξ5) = (ξ1, ξ3−h, ξ3, ξ3+h, ξ3+2h), with knot intervals ∆= (∆1, ∆2, ∆3, ∆4) = (x, h, h, h), turns into the vector Ξ1= (ξ1, ξ2, ξ3, ξ4,ˆ ξ),where ˆ ξ=ξ3+3h, and its knot intervals becomes ∆1= (x, h, h, 2h),see Fig. 5.A.6(b). Then, Ξ2= (ξ2, ξ3, ξ4, ξ5,ˆ ξ) and the original function is split into two new functions B[Ξ] = B[Ξ1]−1 4B[Ξ2].(5.7) Analogously, if the knot vector is extended in the negative direction, i.e., Ξ = (ξ3−2h, ξ3−h, ξ3, ξ3+h, ξ5), Ξ1= (ξ3−3h, ξ3−2h, ξ3−h, ξ3, ξ3+h) and Ξ2= (ξ3−3h, ξ3−h, ξ3, ξ3+h, ξ3+ 2h), then B[Ξ] = −1 4B[Ξ1] + B[Ξ2].(5.8) And, if the knot vector is extended in both directions, i.e., Ξ = (ξ3−2h, ξ3−h, ξ3, ξ3+h, ξ3+ 2h), Ξ1= (ξ3−3h, ξ3−2h, ξ3−h, ξ3, ξ3+h), Ξ2= (ξ3−3h, ξ3−h, ξ3, ξ3+h, ξ3+ 3h) and Ξ3= (ξ3−h, ξ3, ξ3+h, ξ3+ 2h, ξ3+ 3h),then B[Ξ] = −1 4B[Ξ1] + B[Ξ2]−1 4B[Ξ3].(5.9) 80 5 Properties of the constructed spline spaces 5.A.2.1.2 Basic cases study Proposition 3. For a simple elemental refinement b Te 1⊂b Te 2any m-level function of Sb Te 1can be reproduced with functions from Sb Te 2. Proof. To study thoroughly all possible cases, m-level function supports can be divided in three types according to the position of its anchor with respect to the center of the cell τ, as shown in Fig. 5.A.7. Type 1 corresponds to the anchors aligned with an edge of the cell τ. Type 2 anchors are aligned with a diagonal of the cell and type 3 with a symmetry axis of the cell. For each of the three types we are going to study all possible situations when the function N0,associated to a certain anchor, undergoes some changes after refining the cell τ, and show that it can be recovered by the new space Sb Te 2. Figure 5.A.7: m-level functions from the one-ring neighbourhood of the τare classified in three types depending on the position of the anchor with respect to the cell τ. We adopt the following notation: ◦The function of the mesh b Te 1to study is denoted by N0. ◦The blending functions from the refined mesh b Te 2are denoted by b Ni. ◦The functions that belong to the new space Sb Te 2are denoted by ¯ Ni. That is, the function ¯ Nican be a blending function of b Te 2or a linear combination of those. Type 1. For this type of function only one of its local knot vectors can be changed after the refinement of the cell τ, and only the scenario 1 can take place. Without loss of generality, suppose that the initial knot vector of the function N0is H= (η1, η2, η3, η4, η5),and after the cell refinement the new function b N0associated to the same anchor, has its knot vector H2= (ˆη, η2, η3, η4, η5),see Fig. 5.A.8. Then, this situation can be seen as an insertion of the new knot ˆηbetween the first and the second knot of the vector H= (η1, η2, η3, η4, η5).Then, in virtue of (5.3) and (5.6) N0=B[Ξ]B[H] = 1 4B[Ξ]B[H1] + B[Ξ]B[H2] = 1 4b N1+b N0, 81 5 Properties of the constructed spline spaces (a) Modification of the function after the refinement. (b) N0=1 4b N1+b N0. Figure 5.A.8: Type 1 function N0is always split into two blending functions from the new space. where H1= (η1,ˆη, η2, η3, η4) and H2= (ˆη, η2, η3, η4, η5).It is easy to see that the functions b N1and b N0are new blending functions from Sb Te 2, see Fig. 5.A.8. Hence N0∈Sb Te 2. Note that, according to the formula (5.2), the value of the knot η5is irrelevant for the splitting of the function N0. It could be η5=η4+ 2hor η5=η4+hor η5=η4, see Fig. 5.A.9(a-b). Since the splitting is performed for ηdirection, the values of ξ1and ξ2are also irrelevant. Another possible case, besides the considered here for the Ξ knot vector, is ξ1=ξ2=ξ3,see Fig. 5.A.9(c). Combination of all possible values of η5and ξ1gives en total 6 = 3 ×2 possible situations for type 1 function N0. For all these cases the linear combination for the splitting of the function N0is the same and the functions b N1,b N0are blending functions of the mesh b Te 2.From now on, to denote which knot values are irrelevant for the subject of our study, we are going to mark them with red font in the figures. We are going to mention the other possible knot values, but we are not going to study mesh configurations corresponding to these cases, since they are equivalent to the studied case. Type 2. The only possible situation when this type of function can undergo some changes, after the refinement of the cell τ, is the following. For the mesh b Te 1the support of the function N0was extended due to Extension rule 2, because its corner 82 5 Properties of the constructed spline spaces (a) η5=η4+h. (b) η5=η4. (c) ξ1=ξ2=ξ3. Figure 5.A.9: Type 1 function N0. Other possible values for the irrelevant knots of the Fig. 5.A.8(b). was not situated over the mesh edges. However, after the refinement of the cell τ, this does not happen any more for the mesh b Te 2. Consequently the support of the function N0is reduced, because no extension is needed, see for example Fig. 5.A.10(a). Two cases are possible: (1) both knot vectors are changed or (2) only one of the knot vectors of N0is changed. For both cases only a reduction of the support (scenario 1) can take place. Let us consider each of two possible cases. Case 1. Both knot vectors of the function N0are modified. Suppose that for the function N0=B[Ξ]B[H] its knot vector Ξ = (ξ1, ξ2, ξ3, ξ4, ξ5) becomes (ξ1, ξ2, ξ3, ξ4,ˆ ξ) and H= (η1, η2, η3, η4, η5) becomes (ˆη, η2, η3, η4, η5) after the refinement of τ. That means that a new knot ˆ ξis inserted between the fourth and the fifth knot of Ξ and a new knot ˆηis inserted between the first and the second knot of H, see Fig. 83 5 Properties of the constructed spline spaces (a) Modification of the function after the refinement. (b) N0=b N0+1 4b N1+1 16 b N3+1 4b N2. (c) N0=b N0+1 4b N1+1 4b N2+1 16 b N3+1 16 b N4+1 16 b N5. Figure 5.A.10: Type 2 function N0, case 1: both knot vectors are modified. 5.A.10(a). Then, in virtue of (5.5), (5.6) and (5.4), the function N0is split into four functions B[Ξ,H] = (B[Ξ1] + 1 4B[Ξ2])(1 4B[H1] + B[H2]) = =1 4B[Ξ1]B[H1] + B[Ξ1]B[H2] + 1 16B[Ξ2]B[H1] + 1 4B[Ξ2]B[H2] = =1 4¯ N2+b N0+1 16 b N3+1 4¯ N1. (5.10) where Ξ1= (ξ1, ξ2, ξ3, ξ4,ˆ ξ), Ξ2= (ξ2, ξ3, ξ4,ˆ ξ, ξ5),H1= (η1,ˆη, η2, η3, η4) and H2= 84 5 Properties of the constructed spline spaces (ˆη, η2, η3, η4, η5).Now, to verify that N0∈Sb Te 2, we have to check that each of the four functions of (5.10) belongs to the new space Sb Te 2. The functions b N0and b N3are always blending functions of b Te 2for any underlying mesh configuration. And the functions ¯ N1and ¯ N2can be directly blending functions of b Te 2or a linear combination of those, depending on underlying mesh configuration. Namely, two different situations are possible. •For the mesh shown in Fig. 5.A.10(b), the functions ¯ N1and ¯ N2are directly blending functions of b Te 2.So, we denote them as b N1=¯ N1and b N2=¯ N2. Thus, from (5.10) we have N0=b N0+1 4b N1+1 16 b N3+1 4b N2. In this case the values of the knots η5and ξ1are irrelevant. Equivalent cases take place for other mesh configurations where η5=η4+hor η5=η4,and ξ1=ξ2−hor ξ1=ξ2.The linear combination for N0would be the same. •For the mesh shown in Fig. 5.A.10(c), the functions ¯ N1and ¯ N2are linear combinations of the blending functions of b Te 2, namely ¯ N1=1 4b N4+b N1and ¯ N2=1 4b N5+b N2.Therefore, from (5.10) follows N0=b N0+1 4b N1+1 4b N2+1 16 b N3+1 16 b N4+1 16 b N5, where b N1=B[Ξ2]B[(ˆη, η2, η3, η4,(η4+η5)/2)], b N4=B[Ξ2]B[(η2, η3, η4,(η4+η5)/2, η5)], b N2=B[((ξ1+ξ2)/2, ξ2, ξ3, ξ4,ˆ ξ)]B[H1], b N5=B[(ξ1,(ξ1+ξ2)/2, ξ2, ξ3, ξ4)]B[H1]. Case 2. Only one knot vector is changed. Suppose that for the function N0=B[Ξ]B[H] its knot vector Ξ = (ξ1, ξ2, ξ3, ξ4, ξ5) becomes (ξ1, ξ2, ξ3, ξ4,ˆ ξ) after the refinement of τ. That means that a new knot ˆ ξis inserted between the fourth and the fifth knot of Ξ, see Fig. 5.A.11(a). Then, in virtue of (5.3) and (5.5), the function N0is split into the two functions N0=B[Ξ]B[H] = B[Ξ1]B[H] + 1 4B[Ξ2]B[H] = b N0+1 4¯ N1, where Ξ1= (ξ1, ξ2, ξ3, ξ4,ˆ ξ) and Ξ2= (ξ2, ξ3, ξ4,ˆ ξ, ξ5).The function b N0is a blending function of b Te 2, and the other function ¯ N1can be expressed as linear combination of blending functions of b Te 2.Two situations are possible. 85 5 Properties of the constructed spline spaces (a) Modification of the function after the refinement. (b) N0=b N0+1 16 b N2+1 4b N1. (c) N0=b N0+1 16 b N2+1 4b N1+1 16 b N3. Figure 5.A.11: Type 2 function N0, case 2: only one knot vector is changed. •For the mesh shown in Fig. 5.A.11(b), the function ¯ N1is a linear combination of two blending functions of b Te 2, i.e., ¯ N1=1 4b N2+b N1, and, therefore, we have N0=b N0+1 4(1 4b N2+b N1)=b N0+1 16 b N2+1 4b N1, where b N0=B[Ξ1]B[H], b N2=B[Ξ2]B[η1,(η1+η2)/2, η2, η3, η4] and 86 5 Properties of the constructed spline spaces b N1=B[Ξ2]B[(η1+η2)/2, η2, η3, η4, η5],see Fig. 5.A.11(b). In this case the value of the knot η5is irrelevant. Other possible values are η5=η4+hor η5=η4. •For the mesh configuration shown in Fig. 5.A.11(c), the function ¯ N1is a linear combination of three blending functions of b Te 2; and the linear combination for the function N0is as follows N0=b N0+1 16 b N2+1 4b N1+1 16 b N3, where b N2=B[Ξ2]B[η1,(η1+η2)/2, η2, η3, η4], b N1=B[Ξ2]B[(η1+η2)/2, η2, η3, η4,(η4+η5)/2], b N3=B[Ξ2]B[η3, η4,(η4+η5)/2, η5],see Fig. 5.A.11(c). Type 3. For this type only one of the knot vectors of N0can be modified after the refinement of the cell τ. However, several different situations can arise: scenario 1 (the support is reduced after the refinement) and scenario 2 (the support is extended). Let us consider each of the two possible scenarios. Scenario 1. One of the knot vectors of the function N0is reduced after the refinement. Suppose that for the function N0=B[Ξ]B[H] its knot vector H= (η1, η2, η3, η4, η5) becomes (ˆη, η2, η3, η4, η5) after the refinement of τ, see Fig. 5.A.12(a). That means that a new knot ˆηis inserted between the first and the second knot of H, see Fig. 5.A.12(a). Then, in virtue of (5.3) and (5.6), the function N0is split into two functions N0=B[Ξ]B[H] = 1 4B[Ξ]B[H1] + B[Ξ]B[H2] = 1 4¯ N1+b N0, where H1= (η1,ˆη, η2, η3, η4) and H2= (ˆη, η2, η3, η4, η5). The function b N0is a blending function of b Te 2for any mesh configuration. And, analogously to anterior cases, the function ¯ N1can be a blending function of b Te 2or can be expressed as linear combination of those. Three different situations are possible. •For the mesh shown in Fig. 5.A.12(b) the function ¯ N1is a blending function of b Te 2.So, b N1=¯ N1, and we have N0=b N0+1 4b N1. For this case the values of the knots ξ1, ξ5and η5are irrelevant. Other possible values are η5=η4+hor η5=η4,and ξ1=ξ2−hor ξ1=ξ2,and ξ5=ξ4+hor ξ5=ξ4. 87 5 Properties of the constructed spline spaces (a) Modification of the function after the refinement. (b) N0=b N0+1 4b N1. (c) N0=b N0+1 4b N1+1 16 b N2, (d) N0=b N0+1 4b N1+1 16 b N2+1 16 b N3, Figure 5.A.12: Type 3 function N0, scenario 1. The knot vector His reduced after the refinement. 88 5 Properties of the constructed spline spaces •For the mesh shown in Fig. 5.A.12(c), the function ¯ N1is a linear combination of two blending functions of b Te 2:¯ N1=1 4b N2+b N1.Then N0=b N0+1 4b N1+1 16 b N2, where b N1=B[((ξ1+ξ2)/2, ξ2, ξ3, ξ4, ξ5)]B[H1], b N2=B[(ξ1,(ξ1+ξ2)/2, ξ2, ξ3, ξ4)]B[H1]. •For the mesh shown in Fig. 5.A.12(d), the function ¯ N1is a linear combination of three blending functions of b Te 2, namely ¯ N1=1 4b N2+b N1+1 4b N3,then N0=b N0+1 4b N1+1 16 b N2+1 16 b N3, where b N1=B[((ξ1+ξ2)/2, ξ2, ξ3, ξ4, ξ5)]B[H1], b N2=B[(ξ1,(ξ1+ξ2)/2, ξ2, ξ3, ξ4)]B[H1] b N3=B[(ξ2, ξ3, ξ4,(ξ4+ξ5)/2, ξ5)]B[H1]. Scenario 2. One of knot vectors of the function N0is extended after the refinement. Suppose that for the function N0=B[Ξ]B[H] its knot vector Ξ = (ξ1, ξ2, ξ3, ξ4, ξ5) becomes (ˆ ξ, ξ2, ξ3, ξ4, ξ5) after the refinement of τ, see Fig. 5.A.13(a). Then, the function N0can be expressed as linear combination of two functions N0=B[Ξ]B[H] = 1 4B[Ξ1]B[H] + B[Ξ2]B[H] = −1 4¯ N1+b N0, where Ξ1= (ˆ ξ, ξ1, ξ2, ξ3, ξ4) and Ξ2= (ˆ ξ, ξ2, ξ3, ξ4, ξ5). The function b N0is a blending functions of b Te 2for any mesh configuration. And, analogously to anterior cases, the function ¯ N1can be a blending function of b Te 2or can be expressed as linear combination of those. Two different situations are possible. •For the mesh shown in Fig. 5.A.13(b), the function ¯ N1is a blending function of b Te 2.So, we have N0=b N0−1 4b N1. For this case the values of the knots ξ5and η5are irrelevant. Several equivalent cases take place for other mesh configurations where η5=η4+hor η5=η4,and ξ5=ξ4+hor ξ5=ξ4. 89 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary parameterization of the patches with cube faces partition of the input triangulation in 6 patches (a) Surface paremeterization of the solid. adapted tetrahedral mesh (b) Adapted tetrahedral mesh of the cube (Kossaczk´y refinement). (c) Optimization of tetrahedral physical mesh. Figure 6.2: Main steps of the Meccano method. 96 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary Figure 6.3: T-spline representation of the solid obtained with Meccano method. 6.2 Parameterization method for 2D geometries Method for spline parameterization of 2D geometries is based on the idea of the Meccano method and a novel T-mesh optimization procedure. 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 spline representation is calculated by imposing the interpolation conditions on points sited both on the interior and on the boundary of the geometry. Next, we revise briefly the proposed method. A full version can be found in [63]. 6.2.1 General scheme of the method 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 97 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary 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 mesh untangling and smoothing procedure. 3. Construction of a T-spline representation of the geometry: The T-spline mapping is obtained by imposing interpolation conditions. As interpolation points, we take the vertices of the physical T-mesh obtained after the optimization process. 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 spline parameterization. 6.2.2 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. 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 ϵ. 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 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. Figure 6.1(a) shows an example of the adapted parametric T-mesh constructed in this stage; the tangled mesh obtained after projecting the boundary of parametric domain to the input boundary polygonal is shown in Fig. 6.1(b); and Fig. 6.1(c) illustrates the resulting optimized T-mesh. 6.2.3 T-mesh optimization The key of the proposed parameterization method lies in the optimization procedure that allows to obtain a high quality physical T-mesh, which is used to construct the T-spline representation of the object. 98 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary (a) (b) (c) Figure 6.1: Stages of T-mesh construction. (a) Parametric T-mesh adapted to the boundary of the geometry; (b) tangled physical mesh after projecting the boundary of parametric domain to the input boundary polygonal; (c) optimized physical T-mesh. The colors represent the correspondence between the parametric and physical boundaries. The mesh optimization process is carried out by iterative relocation of each inner node of the mesh in such a way that 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 node which is called free node. The local objective function for a free node is based on algebraic shape quality metrics proposed by Knupp in [64, 65] for quadrilateral elements. Shape quality metric for a triangle is defined in terms of the Jacobian matrix of the affine mapping from the 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 degenerate. The distortion metric of an element is defined as the inverse of its quality metric. Shape quality metric for a quadrilateral element is based on the set of shape quality metric of its simplicial elements that compose the given quadrilateral. Therefore, in order to asses the quality of the cells, we have to decompose each 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. A regular node is surrounded by four cells with equal or different sizes. Local submesh is decomposed in twelve triangles, three triangles per cell, whose qualities depends on the position of the free node. 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. 6.2 illustrates the decomposition of a local submesh for each case and the feasible region of the objective function. Feasible region is the subset of R2, where the free node can be placed, being the local submesh valid. The optimal position of each free node is determined by minimizing a local 99 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary (a) (b) (c) (d) Figure 6.2: Triangular decomposition of the 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 (red lines) and feasible region (light blue) 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. 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, during optimization process, 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 as similar to the parametric one as possible. Let Tbe a triangle whose vertices are given by xk= (xk, yk)T∈R2, k = 0,1,2. Consider that TIis our ideal or target triangle whose vertices are v0,v1and v2. And let the matrix Sbe the Jacobian matrix of the affine map that takes TIto T. Quality metrics of the triangle Tcan be defined in terms of the matrix S. For example, the mean ratio,q=2σ ∥S∥2, is an easily computable algebraic quality metric of T, where σ= det (S) and ∥S∥is the Frobenius norm of S. The maximum value of qis the unity, and it is reached when S=µR, where µis a scalar and Ris a rotation matrix. In other words, q= 1 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. Let x= (x, y)Tbe the position of the free node, and let Smbe the Jacobian matrix of the m-th triangle connected to this free node. We define the distortion measure for m-th triangle as ηm=∥Sm∥2 2σm . Then, the local objective function, used for mesh quality improvement, is defined 100 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary by means of distortion measures of each triangle of the local submesh K(x) = M ∑ m=1 ηm= M ∑ m=1 ∥Sm∥2 2σm , where Mis the number of triangles in the local submesh. The objective function Kbecomes discontinuous when the area of any triangle goes to zero. Due to these singularities, the function Kimproves the quality of valid elements but it does not work properly, when the mesh is tangled (σ≤0). In [61] we proposed a modification of Kby replacing σby the positive and increasing function h(σ) = 1 2(σ+√σ2+ 4δ2). Then, the modified distortion becomes η∗ m=∥Sm∥2 2h(σm). This modification eliminates the barriers associated with their singularities and the new objective function K∗=∑M m=1 η∗ mbecomes smooth in R2, so the unconstrained optimization problem can be easily solved with any standard method. In the feasible region the modified objective function K∗approximates the original function Kas δ→0, and then the minimum of the original and modified objective functions are nearly identical when δis small. When this region does not exist, the minimum of the modified objective function is located in such a way that it tends to untangle the local mesh. Thus, the modified objective function allows the simultaneous untangling and smoothing of the mesh. 6.2.4 Construction of a spline representation of the geometry We have to obtain a global one-to-one parametric transformation that maps the parametric domain into the physical domain S:b Ω = [0,1]2→Ω. Spline representation of physical domain is build as lineal combination of spline blending functions defined over the adapted parametric T-mesh T S(ξ) = ∑ α∈AT PαNα(ξ) where Pα∈R2is the control point corresponding to the α-th blending function. Control points Pαare found by imposing interpolation conditions. As interpolation points we use the anchors of the blending functions. Each anchor ξα coincides with a T-mesh vertex and its position in the physical space xαwas determined by the mesh optimization process, see Fig. 6.3. Note that for our EP-splines only regular vertices of the mesh are used for interpolation. For T-spline functions also hanging nodes (T-junctions) form a set of interpolation points. 101 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary Figure 6.3: Spot test model with 844 cells, 1174 control points. Parametric domain and spline representation of the physical domain. Finally we solve the linear system of equations xβ=S(ξβ)=∑ α∈AT PαNα(ξβ),∀ξβ, β ∈AT where ξβare interpolation points in parametric space and xβare their images in the physical space. Resulting spline parameterization of Spot test model is shown in Fig. 6.3. 6.2.5 Quality assessment and its improvement 6.2.5.1 Mean ratio Jacobian Our objective is a high-quality geometry parameterization suitable for isogeometric analysis. High distortion of the geometry can produce a large variation of the Jacobian that can lead to a poor accuracy in the numerical results. A good uniformity and orthogonality of the isoparametric curves are desired for the parametric mapping S. In order to assess the quality of the constructed parametric transformation we analyse the mean ratio Jacobian given by Jr(ξ) = 2 det (J) ∥J∥2, where Jis the Jacobian matrix of the mapping Sat the point ξ= (ξ, η) and ∥J∥ is 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 P′ 0=S(P0), as illustrated in Fig. 6.4. 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 102 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary Figure 6.4: Mean ratio Jacobian. A quality metric of the parametric mapping Sat any point P0in terms of the mean ratio of the triangle P′ 0P′ 1P′ 2. orthogonality and uniformity of its 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. Figure 6.5 shows the resulting spline representation of the Spot test model and the colormap of the mean ratio Jacobian. (a) parametric domain. (b) physical domain. Figure 6.5: Spot test model. Colormap of mean ratio Jacobian. 103 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary (a) (b) self-intersections. (c) no self-intersections. (d) negative Jacobian. (e) positive Jacobian. Figure 6.6: Adaptive refinement strategy to improve the parameterization quality in Gran Canaria Island domain. (a) Spline representation of the domain; (b) initial spline parametrization with negative Jacobian; (c) spline parameterization with no negative Jacobian after applying adaptive refinement; (d) mean ratio Jacobian of the initial parametrization; (e) mean ratio Jacobian of the final parametrization. 6.2.5.2 Adaptive refinement Parameterization of complex geometries entails a severe distortion that can lead to 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 have to increase the number of degrees of freedom in the area. We propose an adaptive strategy that refines the cells with low quality. 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 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. Figure 6.6 illustrates the proposed strategy. Additional refinements were applied to Gran Canaria Island domain with δ= 0.2.The initial mesh with 3439 cells produces a spline parametric mapping with low quality in some areas and selfintersections in the North East part of the island. After adaptive refinement we obtain 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. Figure 6.7 illustrates spline representation and the colormap of the mean ratio Jacobian of the Gran Canaria Island geometry. 104 6 Parameterization method for complex 2D and 3D geometries from representation of its boundary (a) Parametric domain. (b) Physical domain. Figure 6.7: Gran Canaria Island geometry with 3577 cells and 6054 control points. Colormap of the mean ratio Jacobian. 6.3 Parameterization method for 3D geometries To extend the parameterization method to 3D it is necessary to have an optimization technique for 3D hexahedral T-mesh, which is not a trivial task. This issue is studied in another doctoral dissertation of the group. For now we have implemented an optimization procedure for the regular conformal hexahedral meshes with no T-junctions. Using this procedure we can obtain volumetric parameterization for slightly distorted solids. For example, Fig. 6.1 shows a solid designed using Rhinoceros CAD software. The surface of the model is formed by six NURBS surfaces that can be exported and used as input data for construction of volumetric parameterization of the solid by means of our method. This slightly deformed shape can have one-patch volume parameterization using a unit cube as parametric domain. A uniform 8 ×8×8 hexahedral mesh is sufficient to represent well all the details of the model, since the input boundary parameterization of each face is defined over 4 ×4 parametric grid. The inner nodes of the physical mesh is relocated by applying our optimization procedure for hexahedral meshes, see [66]. Analogously to 2D case, the optimization is based on mean ratio quality metric for a tetrahedron given by q=3σ2/3 ∥S∥2, where σ= det (S) and ∥S∥is the Frobenius norm of S, being Sthe Jacobian matrix of the map from the target element to the physical one. Trivariate spline mapping Sis obtained by imposing interpolation conditions, where as interpolation points we take the vertices of the mesh. Figure 6.2 illustrates 105 7 Testing of the strategy. Computational examples (a) Input surface triangulation, 8650 nodes. (b) Spline representation of the surface. Figure 7.2: Igea’s face. Comparison of the input triangulation and spline approximation of the surface. 7.2.2 Adaptive refinement for interpolation problem In this example we interpolate the function with strong singularity u(r) = r1 2,(7.3) defined in the square domain [0,1]2, being r=√(x−0.5)2+ (y−0.5)2. Spline approximation of the function (7.3) is built as a linear combination of our blending functions uh(ξ) = ∑ α∈AT cαNα(ξ). And the control values cαare found by imposing the interpolation conditions u(ξβ) = ∑ α∈AT cαNα(ξβ),∀ξβ, β ∈AT, where the interpolation points ξβare the anchors of the functions, i.e., regular vertices of the mesh. A standard choice for an interpolation with B-splines is 112 7 Testing of the strategy. Computational examples (a) (b) (c) (d) æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à à à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò 101 102 103 104 105 100 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 10-9 10-10 numberof degrees of freedom error ò Uniform, H1 ì Uniform, L2 à Adaptive, H1 æ Adaptive, L2 -2 -1 (e) ∥uh−u∥L2(Ω) and |uh−u|H1(Ω). Figure 7.3: Adaptive refinement for interpolation with EP-splines. (a) Function to interpolate. (b) Initial mesh. (c) and (d) some steps of the adaptive refinement. (f) Error convergence of the adaptive refinement in L2-norm and H1-seminorm. Greville abscissae, which will be discussed in more details in the section 7.2.5. Generalization of this concept for T-mesh leads to so called Greville collocation points, which for our basis functions coincides with the anchors of the functions, except for the some functions with repeated boundary knots. 113 7 Testing of the strategy. Computational examples æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à à à à à à à ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò 101 102 103 104 105 100 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 10-9 10-10 numberof degrees of freedom error -2 -1 ò Lagrange, H1 ì Lagrange, L2 à Splines, H1 æ Splines, L2 (a) ∥uh−u∥L2(Ω) and |uh−u|H1(Ω). (b) Final mesh. Figure 7.4: Comparison with the interpolation using cubic Lagrange basis functions. (a) Error convergence. (b) Final triangular mesh using Kossaczk´y refinement algorithm. Adaptive refinement is performed according to the indicator based on the exact L2interpolation error: η(Ωe)2=∥u−uh∥2 L2(Ωe)=∫Ωe (u−uh)2dΩ. A cell Ωeis marked to be refined if η(Ωe)> γ maxi{η(Ωi)}, being γ∈[0,1]. Some steps of the adaptive refinement for the function interpolation (7.3) are shown in Fig. 7.3(b)-(d). Figure 7.3(e) illustrates the convergence in L2norm and H1-seminorm for the uniform and adaptive refinement. Optimal rate of convergence is observed for L2error and suboptimal rate for H1error, which can be considered normal, since the error indicator, used for refinement, is based on the exact L2error. Figure 7.4(a) shows the error convergence for the interpolation performed using cubic Lagrange basis functions defined over triangulation, where the adaptive refinement is performed using Kossaczk´y refinement algorithm. The resulting convergence behaviour is similar to one performed with spline functions. Figure 7.4(b) shows the final triangular mesh of the adaptive refinement with Lagrange functions. 114 7 Testing of the strategy. Computational examples 7.2.3 Poisson problem in a square domain In this subsection we present an example of the resolution of a Poisson problem in a square domain Ω = [0,1]2 {−∆u=fin Ω, u=gon ∂Ω.(7.4) In this case the parametric domain coincides with the physical one, and we use the parametric mapping Swhich is the identity. The problem (7.4) is set up in such a way that its exact solution is a function given by (a) (b) (c) (d) Figure 7.5: Isogeometric analysis with EP-splines for the Poisson problem (7.4). (a) Numerical solution corresponding to the final refinement. (b) Initial mesh. (c), (d) Several steps of the adaptive refinement. 115 7 Testing of the strategy. Computational examples æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à 101 102 103 104 105 10-1 10-2 10-3 10-4 10-5 10-6 10-7 numberof degrees of freedom L2error à uniform æ adaptive -2 (a) L2error convergence, ∥uh−u∥L2(Ω). æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à 101 102 103 104 105 100 10-1 10-2 10-3 10-4 10-5 numberof degrees of freedom H1error à uniform æ adaptive -3/2 (b) H1error convergence, |uh−u|H1(Ω). æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à 101 102 103 104 105 30 80 325 1300 4000 numberof degrees of freedom conditionnumber à uniform æ adaptive (c) Condition number of the stiffness matrix. æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à 101 102 103 104 105 0.002 0.01 0.04 0.1 0.4 1 numberof degrees of freedom density à uniform æ adaptive (d) Density of the stiffness matrix. Figure 7.6: Error convergence, condition number and density of the stiffness matrix for adaptive IGA for the Poisson problem (7.4). u(x, y) = exp (−7√(x−0.5)2+ (y−0.5)2)+ + exp (−7√(x−0.25)2+ (y−0.25)2)+ + exp (−7√(x−0.75)2+ (y−0.75)2). Some steps of adaptive refinement for the problem (7.4) and the numerical solution corresponding to the final refinement iteration are shown in Fig. 7.5. As expected, the error estimator has refined near the three singularity points. The convergence behaviour of the adaptive refinement in L2-norm and H1-seminorm are shown in Fig. 7.6(a) and (b). Suboptimal rates of convergence are obtained for the uniform refinement due to the presence of singularities, and the adaptive refinement gains the optimal rates both for L2and H1errors. The evolution of 116 7 Testing of the strategy. Computational examples the density of the stiffness matrix Aand its condition number are represented in Fig. 7.6(c) and (d), respectively. The density of the matrix Ais the fraction of non-zero elements and the condition number is defined as κ(A) =∥A∥2∥A−1∥2. It can be observed that the density and the condition number during the adaptive refinement evolve similarly to the uniform refinement. 7.2.4 Poisson problem in a complex domain In the next example, we present the result of solving a Poisson problem in a complex domain Ω, which is a puzzle piece, see Fig. 7.7. Parameterization of the computational domain is performed using the algorithm described in Section 6.2. This technique, based on a T-mesh optimization procedure, allows to obtain a good quality parametric mapping suitable for application of IGA. The mean ratio Jacobian Jr(ξ) = 2 det(J) ∥J∥2is used to evaluate the quality of the parameterization in the sense of its orthogonality and uniformity. Figure 7.7 shows the resulting parameterization and the colormap of its mean ratio Jacobian. Exact solution for the Poisson problem with Dirichlet boundary condition is a function with steep wave front given by u(r) = arctan(α(r−r0)), where r=√(x−xc)2+ (y−yc)2. Parameter αdetermines the steepness of the wave front and r0is its location. In this example α= 200 and r0= 0.6.The center of the wave front (xc, yc) = (0,0) is situated outside of the computational domain, so the function is smooth in Ω. The numerical solution of the problem and the mesh corresponding to the final refinement iteration are shown in Fig. 7.8(a)-(d). As expected, the error estimator has marked for refinement the zone of the wave front. The evolution of the exact error in L2-norm and H1-seminorm are shown in Fig. 7.8(e). Comparison with the uniform refinement is given in Fig. 7.8(f). Optimal rates of convergence are obtained for both adaptive and uniform refinement, due to the smoothness of the exact solution. 117 7 Testing of the strategy. Computational examples (a) Parametric domain. (b) Physical domain. (c) Colormap of the mean ratio Jacobian. (d) Colormap of the mean ratio Jacobian. Figure 7.7: Puzzle piece domain. Parameterization of the computational domain for the Poisson problem and its quality (mean ratio Jacobian). 118 7 Testing of the strategy. Computational examples (a) Final refinement in the parametric domain. (b) Final refinement in the physical domain. (c) Numerical solution in the parametric domain. (d) Numerical solution in the physical domain. æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì 104 3104 102 101 100 10-1 10-2 10-3 10-4 10-5 10-6 numberof degrees of freedom error -2 -3/2 ì H1error à L2error æ error estimator (e) Error convergence for the adaptive refinement. æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à ì ì ì ò ò ò 3103 104 3104 101 100 10-1 10-2 10-3 10-4 10-5 10-6 10-7 numberof degrees of freedom error ò uniform, H1 ì uniform, L2 à adaptive, H1 æ adaptive, L2 -3/2 -2 (f) Error convergence for the uniform and adaptive refinements. Figure 7.8: Results of the IGA adaptive refinement for the Poisson problem in the puzzle piece domain. 119 7 Testing of the strategy. Computational examples 7.2.5 Collocation method for Poisson problem in a complex domain Another advantage of the higher regularity of spline basis functions is the possibility to apply collocation method based on the strong formulation of the problem. The major benefit of collocation approach is the reduced number of evaluation points compared with the Galerkin method. Besides, the bandwidth and matrix sparsity are significantly reduced compared with a Galerkin method, which improves the performance of iterative solvers. Recently IGA collocation method was investigated in several works. A theoretical background for 1D case and some computational examples for 1D, 2D and 3D problems using NURBS-based collocation method was given in [9]. The efficiency of the method in comparison with IGA and standard FEM Galerkin approach was analysed in [10]. The authors report that IGA collocation outperforms in efficiency both Isogeometric Galerkin and Finite Element Methods. Adaptive isogeometric collocation using hierarchical refinement scheme was shown in [10]. Application of isogeometric collocation using analysis-suitable T-splines for secondand fourth-order problems was presented in [11]. Here, we solve the same Poisson problem as in the previous section. So, we have the strong form of the Poisson problem with Dirichlet boundary condition {−∆u=fin Ω, u=gon ∂Ω,(7.5) where Ω is the puzzle piece domain from the previous section. Numerical solution uhis searched as linear combination of the basis functions uh=∑ i∈I ciNi. Collocation method consists in enforce that uhsatisfies the differential equation and the boundary condition in a certain set of points {xi}that are called collocation points. That is, {−∆uh(xi) = f(xi)∀xi∈Ω, uh(xi) = g(xi)∀xi∈∂Ω.(7.6) Approximation space should be at least C2at collocation points, which is the case for our cubic splines. Analogously to Galerkin method, we split the uhin two parts: uh=∑ i∈I ciNi=∑ i∈I0 ciNi+∑ i∈I\I0 ciNi=∑ i∈I0 ciNi+gh. Coefficients {ci}i∈I\I0are found by interpolating Dirichlet boundary condition, i.e. 120 7 Testing of the strategy. Computational examples from the system ∑ i∈I\I0 ciNi(xj) = g(xj)∀xj∈∂Ω. And the coefficients {ci}i∈I0are found from the linear system −∑ i∈I0 ci∆Ni(xj)−∆gh(xj) = f(xj)∀xj∈Ω, j = 1, . . . , m. That is, Kc =f, where Kj,i =−∆Ni(xj), j = 1, . . . , m and i= 1, . . . , n = #I0. fj=f(xj) + ∆gh(xj), j = 1, . . . , m c={c1, c2, . . . , cn}are the unknown coefficients. The system is determined if the number of collocation points coincides with the number of basis functions. Otherwise, if m > n, the system is overdetermined and can be solved in a least-squares sense. 7.2.5.1 Greville abscissae as collocation points For a stability and good behaviour of the method is crucial to use a suitable set of collocation points. A standard option, working with spline, is Greville abscissae, which are the averages of the knots. Suppose we are given a global knot vector {ξ1, ξ2, . . . , ξn+p+1}that generates nB-spline functions, then its Greville abscissae are defined as ¯ ξi=ξi+1 +ξi+2 +···+ξi+p p, i = 1, . . . , n. These points correspond approximately to the points, where the B-spline functions attain their maximum and coincide exactly with them for an uniform knot vector, except for the functions with repeated knots on the boundary. Other possible options for collocation points are quadrature points, the maxima of spline basis functions, and the Demko abscissae. Some comparison between Greville abscissae and Demko abscissae can be found in [9], which shows that Greville abscissae may not lead to a proper convergence only in very unusual cases and, in general, it is the simplest and reliable option. Here, we choose Greville abscissae as collocation points. As we deal with T-mesh and we do not have global knot vector, we obtain a set of Greville points in parametric space as the knot average of each local knot vector, that is, for the local knot vectors Ξi={ξ1, ξ2, ξ3, ξ4, ξ5}and Hi= {η1, η2, η3, η4, η5}of the function b Niits Greville point is defined as ¯ ξi= (¯ ξi,¯ηi), where ¯ ξi=ξ2+ξ3+ξ4 3and ¯ηi=η2+η3+η4 3. 121 7 Testing of the strategy. Computational examples æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à ì ì ì ì ì ò ò ò ò ò 102 103 104 101 100 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 numberof degrees of freedom error -3/2 -2 ò uniform, H1 ì uniform, L2 à adaptive, H1 æ adaptive, L2 (a) Adaptive and uniform refinement, IGA. æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à ì ì ì ì ì ì ì ì ì ì ì ì ò ò ò ò ò ò ò ò ò ò ò ò 102 103 104 105 101 100 10-1 10-2 10-3 10-4 10-5 10-6 10-7 numberof degrees of freedom error ò FEM, H1 ì FEM, L2 à IGA, H1 æ IGA, L2 -3/2 -2 (b) Adaptive refinement, IGA and FEM. æ æ æ æ æ à à à à à ì ì ì ì ì ì ì ò ò ò ò ò ò ò 102 103 104 105 101 100 10-1 10-2 10-3 10-4 10-5 10-6 numberof degrees of freedom error ò FEM, H1 ì FEM, L2 à IGA, H1 æ IGA, L2 -2 -3/2 (c) Uniform refinement, IGA and FEM. Figure 7.14: Singularly perturbed elliptic equation (7.7). Error convergence and its comparison with FEM. 128 7 Testing of the strategy. Computational examples 7.2.7 Eigenvalue problem in a square domain Eigenvalue problem arises in different context of physics and engineering, e.g. structural vibration analysis, acoustic and electromagnetic wave phenomena, optics. In this section we test the behaviour of Isogeometric Analysis for eigenvalue problem in comparison with Finite Element Method. Considerably superior accuracy of IGA, comparing with Finite Element Analysis, for spectrum analysis was reported in [4]. We consider an eigenvalue problem in the domain Ω = [0,1]2that consist in finding λ∈Rand u= 0 such that      −∆u=λu in Ω, u= 0 on ∂Ω, ∥u∥= 1. (7.9) In this simple case the analytical solution of this problem is known. It has a countable sequence of eigenpairs given by (λn,m, un,m) = (π2(n2+m2),2 sin (πnx) sin (πmy)), n, m = 1,2, . . . . The variational formulation of the problem (7.9) consists in finding λ∈Rand u∈V0such that {a(u, v) = λ(u, v)∀v∈V0, ∥u∥= 1.(7.10) where a(u, v) = ∫Ω∇u·∇vdΩ,(u, v) = ∫Ω u v dΩ. For the finite dimensional space V0,T ={Ni}n i=1 the solution is searched as uh= ∑n i=1ciNi. The discretized eigenvalue problem for (7.10) consists in determine λh∈Rand uh∈V0,T such that {a(uh, Nj) = λh(uh, Nj)∀Nj∈V0,T , ∥u∥= 1.(7.11) Problem (7.11) leads to a generalized algebraic eigenvalue problem Kc =λhMc,(7.12) being Kj,i =a(Ni, Nj) the stiffness matrix, Mj,i = (Ni, Nj) the mass matrix, c={c1, c2, . . . , cn}is a eigenvector that defines a numerical eigenfunction uh, i.e. 129 7 Testing of the strategy. Computational examples 0 1000 2000 3000 4000 0 50000 100000 150000 n Λn numerical eigenvalues FEM numerical eigenvalues IGA exact eigenvalues (a) Sorted sequences of the exact eigenvalues λnand numerical eigenvalues λn,h obtained with IGA and FEA. 0 1000 2000 3000 4000 0.0 0.5 1.0 1.5 n relative error: HΛn,h-ΛnLΛn FEM IGA (b) Relative error of the numerical eigenvalues: (λn,h−λn)/λnfor IGA and FEM. Figure 7.15: Comparison of the numerical spectrum obtained via IGA and FEM with the exact spectrum of continuous problem. uh=∑n i=1ciNi.Matrices Kand Mare symmetric positive defined. The problem (7.12) have a finite sequence of real eigenvalues 0< λ1,h ≤λ2,h ≤λ3,h ≤ ···λn,h, which correspond to eigenfunctions u1,h, u2,h, u3,h,··· , un,h. Here we perform two numerical experiments. The purpose of the first one is to compare the numerical spectrum obtained with IGA and FEA. We solve the eigenvalue problem (7.9) using Isogeometric Analysis over a uniform 64 ×64 mesh, which corresponds to 4225 basis functions, and using Finite Element Method with 3rd order Lagrange functions over a uniform triangular mesh with a similar number of degrees of freedom (4513). We compare the sequence of numerical eigenvalues λn,h, sorted in increasing order, with the exact values λn. The corresponding sequences of eigenvalues are plotted in Fig. 7.15(a). Isogeometric result includes 4225 eigenvalues and FEM sequences contains 4513 values. Abscissa of each point represents the index of the eigenvalue. Figure 7.15(b) shows relative error for the same sorted numerical spectrum: (λn,h−λn)/λn. As can be seen, the FEM eigenvalues have considerably superior error and the difference becomes bigger for the highest eigenvalues. In the second experiment we perform an adaptive refinement in order to approximate well the 11-th eigenfunction u11. The index 11 corresponds to a simple 130 7 Testing of the strategy. Computational examples (a) Numerical solution for the eigenfunction u11 =u3,3= 2 sin 3πx sin 3πy. (b) (c) (d) (e) (f) (g) Figure 7.16: Eigenvalue problem. Adaptive refinement to approximate the eigenfunction u11 and its eigenvalue λ11. (a) Numerical solution for the final refinement step. (b)-(d) some steps of adaptive refinement with IGA. (e)-(g) some steps of adaptive refinement with FEA. (no repeated) eigenvalue with indexes n= 3, m = 3: (λ11, u11) = (λ3,3, u3,3) = (18π2,2 sin 3πx sin 3πy),see Fig. 7.16(a). 131 7 Testing of the strategy. Computational examples æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à ì ì ì ì ì ì ì ì ì ò ò ò ò ò ò ò ò ò 102 103 104 100 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 numberof degrees of freedom error ò FEM, H1 ì FEM, L2 à IGA, H1 æ IGA, L2 -2 -3/2 (a) Error convergence for the eigenfunction: ∥u11,h −u11∥L2and |u11,h −u11|H1. æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à 102 103 104 101 100 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 10-9 numberof degrees of freedom eigenvalue error: Λh-Λ -3 à FEM æ IGA (b) Error convergence for the eigenvalue: (λ11,h −λ11). Figure 7.17: Eigenvalue problem. Adaptive refinement to approximate the eigenpair (λ11, u11). Convergence behaviour. Residual-type error estimator for k-th eigenpair is given by η(Ωe)2=h2∥(∆uk,h +λk,huk,h)∥2 L2(Ωe)=h2∫Ωe (∆uk,h +λk,huk,h)2dΩ. Some steps of the adaptive refinement for Isogeometric Analysis and Finite Element Analysis are shown in Fig. 7.16. As can be observed, adaptive refinements lead to almost uniform refinement, since the exact solution is a smooth function. Error convergence for the eigenfunction u11 and its eigenvalue are shown in Fig. 7.17. It can be seen that both methods offer the optimal rates of convergence. However, IGA obtains a considerably better accuracy of the numerical solution for the eigenpair (λ11, u11). 7.2.8 Kellogg function This problem was designed by Kellogg while studying Poisson problems with intersecting interfaces. We consider elliptic boundary value problem with discontinuous diffusion coefficient in the domain Ω = [−1,1]2 −div(A(x)∇u) = fin Ω,(7.13) where the coefficient Ais piecewise constant: A=a1in the first and third quadrant, and A=a2in the second and fourth quadrant. For f≡0 the solution is given in polar coordinates by u(r, θ) = raµ(θ), 132 7 Testing of the strategy. Computational examples Figure 7.18: Kellogg function for a= 0.25. where µ(θ) =          cos ((π/2−c)a) cos ((θ−π/2 + b)a) if 0 ≤θ≤π/2, cos (ab) cos ((θ−π+c)a) if π/2≤θ≤π, cos (ac) cos ((θ−π−b)a) if π≤θ≤3π/2, cos ((π/2−b)a) cos ((θ−3π/2−c)a) if 3π/2≤θ≤2π, and the numbers a,band csatisfy the relations                    R:= a1/a2=−tan ((π/2−c)a) cot (ab), 1/R := −tan (ab) cot (ac), R=−tan (ac) cot (π/2−b)a, 0< a < 2, max{0, πa −π}<2ab < min{πa, π}, max{0, π −πa}<−2ac < min{π, 2π−πa}. (7.14) The exact solution has a discontinuous derivative along the axes y= 0, x= 0 and infinite derivative at the origin. It is known that the solution u∈H1+a. Here, we choose a= 0.25. Then, from (7.14), the following parameters for the solution of the problem (7.13) can be obtained: a1≈25.27414236908818, a2= 1, b=π/4 and c≈ −5.49778714378214.We have taken these data from [70]. The function ufor a= 0.25 is shown in Fig. 7.18. Due to the discontinuous diffusion coefficient, the error estimator for this problem includes the term for the gradient jump across the edges of the cell. Namely, the estimator is given by η(Ωe)2=h2∫Ωe (f+ div(A(x)∇uh))2dΩ + h∫∂Ωe[A(x)∂uh ∂n ]2 dΓ. This test problem is a challenging test for two reasons. 133 7 Testing of the strategy. Computational examples æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à à ì ì ì ì ì ò ò ò ò ò 102 103 104 100 10-1 10-2 10-3 10-4 numberof degrees of freedom error ò uniform, H1 ì uniform, L2 à adaptive, H1 æ adaptive, L2 -1 -1/2 (a) Uniform and adaptive refinement, C2IGA. æ æ æ æ æ à à à à à ì ì ì ì ì ò ò ò ò ò 102 103 104 101 100 10-1 10-2 10-3 10-4 10-5 numberof degrees of freedom error ò FEM, H1 ì FEM, L2 à IGA, H1 æ IGA, L2 -0.3 -0.1 (b) Uniform refinement, FEM and C2IGA. æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à à ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò 102 103 104 100 10-1 10-2 10-3 10-4 10-5 numberof degrees of freedom error ò FEM, H1 ì FEM, L2 à IGA, H1 æ IGA, L2 -1/2 -1 (c) Adaptive refinement, FEM and C2IGA. Figure 7.19: Kellogg function. Poor error convergence for the adaptive refinement using C2-continuous EP-spline spaces and its comparison with FEM solution. First, the problem has discontinuous diffusion coefficient, the solution has discontinuous derivative along the interfaces and infinite derivative at (0,0). The C2-continuous space seems to be not appropriate for approximating the solution of this problem. Very poor convergence order was obtained for the adaptive refinement using standard C2-continuous Isogeometric Analysis, see Fig. 7.19. Uniform refinement with C2spline spaces gives a similar to FEM error behaviour, Fig. 7.19(b); and for the adaptive refinement with C2spline spaces the convergence rates are considerably inferior to FEM’s rates, Fig. 7.19(c). None of the methods gained optimal rates with adaptive refinement. To tackle this problem we repeat the knots in the interior of the domain to define spline spaces that are C0continuous along the axes. The resulting error convergence for C0IGA is given in Fig. 7.20. For the uniform refinement with C0spline spaces the order of convergence coincides with the order of the C2case, and the accuracy is improved 134 7 Testing of the strategy. Computational examples slightly with respect to the C2case. However, for the adaptive refinement using C0-continuous spline space the resulting convergence is improved considerably, see Fig. 7.20(c). As can be seen, the convergence rates and accuracy are superior to Finite Element results. Suboptimal convergence order for H1error was obtained for the FEM adaptive refinement. Note that a stagnation of L2error, both for FEM and IGA method, takes place. The final meshes of adaptive refinement for IGA and FEM are shown in Fig. 7.20(d) and (e). Remark 15. Using C0-continuous spline spaces was beneficial for this problem. For that it was necessary to use multiple interior knots along the axes. In this particular case it was easy to create a repeated knots structure along the entire knot line of the domain. This implies to change knot intervals of some functions (repeating the knots) and add some additional functions to obtain a complete polynomial space. However, in general case, it is not so straightforward to define functions with reduced regularity along a certain interval (or point) of the domain, which can be necessary for some problems. For that we need to extend the strategy for T-meshes with repeated interior knots. Second, due to the strong singularity at the origin, during adaptive refinement the estimator marks to refine at each iteration a very reduced zone around the singular point (0,0). That leads to excessive accumulation of function support in this zone, which affects the condition number of the matrix and its sparsity. To avoid this phenomenon, we can apply the strategy explained in Section 5.5. The idea is to impose some additional refinements in order to extend the refined zone and not allow an excessive support accumulation during adaptive refinement. The results of the adaptive refinement using this approach are shown in Fig. 7.21. As can be seen, the application of this strategy allows to “delay” the stagnation of the error convergence and arrive to a higher accuracy of numerical solution. Final adapted mesh, which is more refined compared with the standard adaptive process, is shown in Fig. 7.21(b). It is worth mentioning that IGA refinement process for this problem leads to a very high mesh discretization level. Namely, the depth of the final quadtree is 42. This depth corresponds to the cell size of the order 10−13, which is close to the working precision. This fact, probably, can explain the stagnation of the errors. 135 7 Testing of the strategy. Computational examples æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à ì ì ì ì ì ò ò ò ò ò 102 103 104 100 10-1 10-2 10-3 10-4 10-5 10-6 10-7 numberof degrees of freedom error ò uniform, H1 ì uniform, L2 à adaptive, H1 æ adaptive, L2 -3/2 (a) Uniform and adaptive refinement, C0IGA. æ æ æ æ æ à à à à à ì ì ì ì ì ò ò ò ò ò 102 103 104 100 10-1 10-2 10-3 10-4 numberof degrees of freedom error ò FEM, H1 ì FEM, L2 à IGA, H1 æ IGA, L2 -0.1 -0.3 (b) Uniform refinement, C0IGA and FEM. æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò 102 103 104 100 10-1 10-2 10-3 10-4 10-5 10-6 numberof degrees of freedom error ò FEM, H1 ì FEM, L2 à IGA, H1 æ IGA, L2 -2 -3/2 (c) Adaptive refinement, C0IGA and FEM. (d) Final mesh, IGA. (e) Final mesh, FEM. Figure 7.20: Kellogg function. Isogeometric Analysis using C0-continuous EP-spline spaces and its comparison with FEM solution. 136 7 Testing of the strategy. Computational examples æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò 102 103 104 100 10-1 10-2 10-3 10-4 10-5 10-6 numberof degrees of freedom error ò FEM, H1 ì FEM, L2 à IGA, H1 æ IGA, L2 -3/2 -2 (a) Error convergence, C0IGA. (b) Final mesh of the adaptive process. Figure 7.21: Kellogg function. Isogeometric Analysis using C0-continuous EP-spline spaces. Adaptive refinement imposing additional refinements to avoid excessive support accumulation. 137