Bernard Mourrain

Keywords: spline constructions, multipatch domains, gluing data, extraordinary vertices, isogeometric analysis, approximation theory

In this thesis we investigate novel spline constructions over unstructured meshes to be applied for modeling purposes, approximation problems and in the numerical resolution of partial differential equations.

Being able to describe accurately a complex shape is not an easy task in geometric modeling, and it becomes even harder if we need the result to be suitable for numerical simulations. This is the challenge which motivated the topic of this work: explore new spline constructions that have potential to reproduce faithfully complicated geometries and, at the same time, are suitable to run isogeometric analysis experiments.

First we present the derivation of a globally G 1 smooth family of surfaces, defined by smoothing masks, approximating the well-known Catmull-Clark subdivision surface scheme. The resulting surface is a collection of Bézier patches, which are biquintic and join with G 1 smoothness around extraordinary vertices and bicubic elsewhere. Each Bézier point is computed using a locally defined mask which ensures, by means of quadratic gluing data, G 1 regularity around extraordinary vertices of the corresponding patches.

We continue with the description of a set of basis functions generating the space of biquintic G 1 spline over a quadrangular mesh. The basis is represented in terms of biquintic Bézier polynomials on each quadrilateral face. Starting from the equation defining the G 1 relations between two patches, achieved by quadratic gluing data functions, we perform an extraction procedure in order to obtain the values of the control point defining the different elements.

The latter basis functions, due to their definition, turn out not to be analysis-suitable. Driven by this, we investigate a new construction of G 1 spline basis functions with bidegree 5 suitable for isogeometric analysis simulations. The construction is made considering knot vectors composed of knots of multiplicity 5 and imposing C 1 regularity at the inner part of the resulting patches. Similarly to the earlier case, the coefficients defining the different functions are obtained by an extraction technique plus knot insertion.

The previous constructions are then used to solve two practical problems: the conversion of CAD models into smooth spline objects and the resolution of the shallow-water equation. The conversion is performed by fitting a point cloud representing a discretized model of the original CAD geometry, while the shallow-water equation is solved in the case of shallow lakes whose shape is faithfully approximated by a planar quad mesh.

Lastly, we present the definition of three cubic C 2 quasi-interpolation operators over arbitrary triangulations. The quasi-interpolants are locally generated by simplex spline basis functions defined on each triangle of the triangulation, which is subdivided according to the Wang-Shi split. The coefficients defining the operators in the simplex basis are computed easily by solving an Hermite problem, whose differential data is either given in input or reconstructed by using local cubic polynomials attached to the different features of the triangulation.

Résumé

Dans cette thèse, nous étudions de nouvelles constructions de splines sur des maillages non structurés à appliquer à des fins de modélisation, à des problèmes d'approximation et à la résolution numérique d'équations aux dérivées partielles.

Pouvoir décrire avec précision une forme complexe n'est pas une tâche facile en modélisation géométrique, et cela devient encore plus difficile si nous avons besoin que le résultat soit adapté aux simulations numériques. C'est le défi qui a motivé le sujet de ce travail : explorer de nouvelles constructions de splines qui ont le potentiel de reproduire fidèlement des géométries compliquées et qui, en même temps, sont adaptées à l'exécution d'expériences d'analyse isogéométrique.

Nous présentons tout d'abord la dérivation d'une famille de surfaces globalement G 1 lisses, définies par des masques de lissage, approchant le schéma bien connu de la surface de subdivision de Catmull-Clark. La surface résultante est une collection de points de Bézier, qui sont biquintiques G 1 autour de sommets extraordinaires et bicubiques ailleurs. Chaque point de Bézier est calculé à l'aide d'un masque défini localement qui assure, au moyen de données de collage quadratique, la régularité G 1 autour des sommets extraordinaires des patchs correspondants.

Nous poursuivons avec la description d'un ensemble de bases générant l'espace des splines G 1 biquintiques sur une maille quadrangulaire. La base est représentée en termes de polynômes de Bézier biquintiques sur chaque face du quadrilatère. sur chaque face du quadrilatère. A partir de l'équation définissant les relations G 1 entre deux patchs, obtenue par des fonctions de données de collage quadratique, nous effectuons une procédure d'extraction afin d'obtenir les valeurs du point de contrôle définissant les différentes bases.

Ces dernières bases, en raison de leur définition, s'avèrent ne pas être adaptées à l'analyse. C'est pourquoi nous étudions une nouvelle construction de fonctions de base spline G 1 avec un bidegree 5 convenant aux simulations d'analyse isogéométrique. La construction est réalisée en considérant des vecteurs de noeuds composés de noeuds de multiplicité 5 et en imposant une régularité C 1 dans la partie interne des patchs résultants. Comme pour la construction des bases précédentes, les coefficients définissant les différentes fonctions sont obtenus par une technique d'extraction et d'insertion de noeuds.

Les constructions précédentes sont ensuite utilisées pour résoudre deux problèmes pratiques: la conversion de modèles CAD en objets splines lisses et la résolution de l'équation des eaux peu profondes. La conversion est obtenue en ajustant un nuage de points représentant un modèle discrétisé de la géométrie CAD originale, tandis que l'équation des eaux peu profondes est résolue dans le cas de lacs peu profonds dont la forme est fidèlement approximée par un maillage quadratique planaire.

Enfin, nous présentons la définition de trois opérateurs quasi-interpolants cubiques C 2 sur des triangulations arbitraires. Les quasi-interpolants sont générés localement par des fonctions de base spline simplex définies sur chaque triangle de la triangulation, qui est subdivisée selon la division de Wang-Shi. Les coefficients définissant les quasi-interpolants dans la base simplex sont calculés facilement en résolvant un problème d'Hermite, dont les données différentielles sont soit données en entrée, soit reconstruites en utilisant des polynômes cubiques locaux attachés aux différentes caractéristiques de la triangulation.

Mots-clés: constructions de splines, domaines multipatch, données de collage, vertices extraordinaires, analyse isogéométrique, théorie de l'approximation

Introduction

Describing accurately a complex shape in terms of its main geometric features is a major challenge in geometric modeling. The latter aims at providing compact and efficient models for approximating or representing precisely geometric objects, favoring the exploitation of such models in simulation and optimization processes.

In this thesis, we investigate the construction of new smooth spline spaces over meshes and their applications in geometric modeling, approximation theory and IsoGeometric Analysis (IGA)-based numerical simulations.

Spline representations make a step towards this objective, by describing smooth shapes in terms of control points, which express the features of such parametrized surfaces. Thanks to their simple definition and malleability, splines are applied in a wide range of different topics; one of their main usage can be found in Computer-aided Design and Manufacturing (CAD-CAM) where, typically, surfaces are trimmed and glued together in patchworks to establish the entirety of a geometric object's shape. As a consequence, shape descriptions emerge that might lack precision, showing leaking patches, stretched components, and other deviations that don't align with the desired geometry of the objects.

An alternative approach explored to tackle this challenge is the utilization of subdivision surfaces. Beginning with a coarse mesh that governs the surface's geometric characteristics, a subdivision process is sequentially employed to achieve progressively finer meshes, ultimately converging toward a limit surface. Subdivision creates smooths objects but, unfortunately, they correspond to an infinite set of spline rings nearby an Extraordinary Vertex (EV) of the initial mesh. Moreover, the surface exhibits unwanted oscillations around EVs.

As another option to model smooth spline objects over meshes of any topology we have multipatch constructions. In this approach, the sought geometry is defined as a collection of single spline functions joining each other with a prescribed smoothness. The tangent plane continuity, or G 1 continuity, extensively investigated in the past few decades, is a proper example of commonly utilized smoothness in multipatch geometries.

Other types of regularity, e.g. standard C r , are more feasible when dealing with triangulations where the continuity can be simply achieved by using macro-elements split. This approach is common in the construction of quasi-interpolation operators or B-spline-like bases.

The investigation of novel smooth spline constructions over unstructured meshes is a revived topic due to its importance for isogeometric analysis discretizations. In fact, high quality results can be achieved by using G 1 basis functions when running numerical simulations to solve partial differential equations, also when dealing with challenging high-order equations. employee of Citroën, developed a numerically stable algorithm for evaluating the above curves.

In [START_REF] Curry | On Pólya frequency functions IV: The fundamental spline functions and their limits[END_REF] splines were described as linear combination of "fundamental spline functions"; these functions, subsequently renamed B-splines, were deeply investigated, for example, in [START_REF] Marsden | An identity for spline functions with applications to variationdiminishing spline approximation[END_REF][START_REF] De Boor | On calculating with B-splines[END_REF]. During these years there were already some results for tensor-product splines, but their most innovative application at that time was in the context of subdivision surfaces early introduced by Catmull and Clark in [START_REF] Catmull | Recursively generated B-spline surfaces on arbitrary topological meshes[END_REF], giving rise to a new fruitful line of research in CAD.

Also, further developments on triangular domains made their entrance: among the many works concerning the development on this topic we cite [START_REF] Frederickson | Triangular spline interpolation. Generalized triangular splines[END_REF][START_REF] Powell | Piecewise Quadratic Approximations on Triangles[END_REF][START_REF] Farin | Subsplines über Dreiecken[END_REF]. The age from 1980 to late 1990s had as main character multivariate spline constructions and their applications. A novel construction dated in those years were box-splines, about which we mention, for instance, [START_REF] Boehm | Calculating with box splines[END_REF][START_REF] Cohen | Discrete box splines and refinement algorithms[END_REF] and the classic monographs [START_REF] De Boor | Box Splines[END_REF]; moreover, in [START_REF] Loop | Smooth Subdivision Surfaces Based on Triangles[END_REF] Loop presented an innovative subdivision scheme whose limit surface is a box-spline.

The so-called geometric continuity began to be utilized in computer graphics, thanks to new works of the time such as [START_REF] Farin | Visually C 2 cubic splines[END_REF], as well as innovative spline-like constructions of which we quote [START_REF] Pottmann | The geometry of Tchebycheffian splines[END_REF][START_REF] Reif | TURBS-Topologically Unrestricted Rational B-Splines[END_REF].

Until then, geometric modeling and numerical simulations had little in common. In fact, engineers were implementing finite element methods for solving partial differential equations approximating the domain by polygonal meshes. Before 2005 there were already few works proposing numerical treatments of PDEs based on spline constructions, but the seminal paper proposing to unify the two can be identified with the work by Hughes, Cottrell and Bazilevs [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF], in which, for the first time, it was rigorously described how to use CAGD tools to efficiently solve differentials problems. Undoubdably this is, to date, one of the most prolific research area in numerical analysis: isogeometric analysis was born.

Context of the thesis

Dealing with smooth objects having complex geometry is not an easy task. One of the mostly used tools to perform geometric modeling are tensor-product B-splines. These functions enjoy good properties, such as local support, nonnegativity, partition of unity; moreover, they form a basis for the space of piecewise polynomials, on a certain rectangular region of the plane, with fixed maximal degree d. Similar constructions presenting analogous properties to the tensor-product case can be achieved on triangular partitions of the plane by using the so-called triangular Bernstein polynomials (see, e.g., [START_REF] Farin | Triangular Bernstein-Bézier patches[END_REF][START_REF] Lai | Spline Functions on Triangulations[END_REF]). Thanks to the interesting characteristics of these functions, they are widely involved in geometric modeling: in fact, in order to reproduce an object in terms of spline surfaces it is only necessary to trace out a piecewise linear approximation of the sought geometry, defining its control net. Nowadays, the literature presents high number of extensions of B-spline constructions in order to satisfy the most varied applications: among all the works we recall Hierarchical B-splines [Kra97; BC13], Generalized B-spline [KS99; MPS11], T-splines [Sed03; TMH21], Tchebycheffian B-spline [Maz04; RMS23], Truncated Hierarchical B-splines (THB) [START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF][START_REF] Evans | Hierarchical B-spline complexes of discrete differential forms[END_REF] and Locally refined B-splines (LR) [START_REF] Dokken | Polynomial splines over locally refined box-partitions[END_REF][START_REF] Bressan | Some properties of LR-splines[END_REF]. Many studies have also been focusing on analysing spline spaces, providing dimension formulas (or upper bounds for the latter) and basis computations. A recurring problem in applications is whether using constructions with high degree, which provide more freedom feasible to obtain smoother results but slower to evaluate numerically, or lower degree splines which are more computationally friendly but allow less flexibility. A possible trade off is the use of spline objects with inner respect to C 1 constructions in both visual quality and performances in numerical simulations. As suggested by the name, this continuity ensures that neighbouring patches share the same tangent plane when crossing their common edge. This result is achieved by using C 1 functions called gluing functions, which are defined along the edges of the considered geometry and solely depend on the valences of the two vertices identifying the edge; we have examples of different gluing data functions in [BMV17; BMX20; BH14].

Representing objects of complex shapes is common in CAD industry. Unfortunately, CAD representations present imperfections in the geometry which are unpleasant: in fact, such objects might present trimmed faces as well as gaps between neighbouring patches, making them non suitable for numerical experiments. In real applications, having proper representations of smooth objects which are also eligible for simulations is a must nowadays: that's a reason why one of the most prolific research topic in multipatch geometries is about isogeometric analysis suitable constructions.

From the seminal paper [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] dated back in 2005 and the subsequent monograph [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF], isogeometric analysis has gained great interest. IGA is an efficient Galerkin-based method to numerically solve partial differential equations whose main idea lays on the isogeometric paradigm: it consists in reproducing in an exact way the domain of interest for our simulation via spline functions and use the same basis utilized in the representation of the domain to solve the target partial differential equation. Many papers provide challenging improvements in efficiency when compared to standard Finite Elements approaches, where the solution is obtained by using a set of local polynomial bases defined on a discretization of the selected domain. Particularly interesting are the constructions of G 1 -smooth B-spline basis functions to run numerical simulations on challenging multipatch geometries, that is the case when dealing with real-life applications. It is impossible to briefly give the proper credit to all the existing works on the topic; among all we refer to [JQ14; Wan18; ZQ19] for analysis-suitable G 1 bases on triangular meshes and to [CST16; KST18; KTC22; TT23] for basis functions defined over quadrilateral meshes. We also mention [START_REF] Hughes | Smooth multi-patch discretizations in Isogeometric Analysis[END_REF] for a more general description of multipatch discretizations and [START_REF] Verhelst | A comparison of smooth basis constructions for isogeometric analysis[END_REF] for a comparison of smooth basis constructions. Adaptive constructions on multipatch domains are also available, e.g. [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: Error estimator and convergence[END_REF][START_REF] Bracco | Isogeometric analysis with C 1 hierarchical functions on planar two-patch geometries[END_REF][START_REF] Bracco | Adaptive isogeometric methods with C 1 (truncated) hierarchical splines on planar multi-patch domains[END_REF][START_REF] Bracco | Adaptive isogeometric phase-field modeling of the Cahn-Hilliard equation: Suitably graded hierarchical refinement and coarsening on multi-patch geometries[END_REF], where the basis functions are locally refined around particular features of the considered geometry. Moreover, regarding G 1 spline basis functions on meshes with mixed topology we mention [START_REF] Mourrain | Dimension and bases for geometrically continuous splines on surfaces of arbitrary topology[END_REF][START_REF] Grošelj | C 1 -smooth isogeometric spline functions of general degree over planar mixed meshes: The case of two quadratic mesh elements[END_REF]. Lately, isogeometric analysis has been also used to solve partial differential equations using the so-called Boundary Element Method (BEM), where the search for a solution of a given equation in a certain domain is translated into some integral concerning only the boundary of the domain of interest. This approach gives rise to the IGA-BEM method, which can be more efficient than the standard approach in certain situations; examples are the works [Hel17; Kos18; Fal18] and references therein.

Motivated by the previous aspects, the objective of this thesis is to present novel smooth spline constructions on unstructured meshes that are both easy to construct and manipulate, and that produces high quality objects approximating complex shapes. We reach this goal by using multipatch spline geometries, whose control points are obtained from simple closed formulas ensuring smooth connections between patches. Being also interested in numerical simulations, we focus our interest on constructions providing spline objects that are also suitable for IGAbased numerical resolution of partial differential equations, i.e. possessing good approximation properties. In particular, this has been achieved investigating the construction of smooth basis functions that generate the sought spline space.

Overview

This thesis is organized in the following chapters:

• Chapter 1 recalls basics concepts and results which are necessary for the description of the novel constructions of this work; it also introduces the notation we will use throughout the thesis. In detail, the notions of B-spline curve and surface are introduced as well as a multivariate generalization given by simplex spline formulation. It also contains an overview on subdivision surfaces and on the numerical treatments of PDEs via isogeometric analysis approach.

• Chapter 2 is dedicated to the explicit construction of a G 1 Bézier surface obtained by means of smoothing masks on a quadrilateral mesh. After presenting the constraints ensuring the sought regularity, an explicit case-by-case solving strategy is deeply investigates. The chapter ends with several numerical experiments that demonstrate the quality of the novel proposal.

• Chapter 3 is devoted to the description and study of a set of basis functions generating the space of G 1 splines over a quad mesh. Following the topology of the given mesh, the different functions are obtained by solving simple linear systems derived from the equations defining the geometric continuity across two adjacent patches. Various numerical experiments confirm their power in point cloud fitting problems.

• Chapter 4 introduces an analysis-suitable extension of the basis functions developed in Chapter 3. This new formulation allows a space refinement through knot insertion which makes the construction suitable for the numerical solution of partial differential equations, as stated by the presented results.

• Chapter 5 gives a practical application of the constructions introduced in Chapters 2 to 4: it concerns a complete pipeline to convert CAD models to G 1 spline objects. The translation is achieved by fitting a point cloud, extracted from the input CAD model, with the G 1 basis functions we propose; this conversion also makes the resulting geometry suitable for numerical simulations. Several computational tests confirm the above.

• Chapter 6 presents another use of the tools developed in Chapters 2 to 4 regarding the free vibrations of a shallow lake. The equation strongly depends on the geometry of the lake's bathymetry, which we accurately reproduce by fitting a point cloud underlying its shape thanks to our G 1 basis functions. The cases of two real lakes are used to show the efficiency of this approach.

Chapter 1

Preliminaries

This chapter is dedicated to the basic concepts in the fields of approximation theory, Computer Aided Geometric Design and numerical simulation of partial differential equations. We start introducing B-spline basis functions, as well as B-spline curves, and then move on to the definition of bivariate tensor-product B-spline. Fundamental geometric algorithms to deal with B-spline object will be also provided, e.g. the knot insertion and the degree elevation algorithm. Later on, making use of the idea of geometric continuity, spline constructions on unstructured domains like quadrilateral mesh or triangulation are presented. Subdivision schemes, a powerful tool to reconstruct smooth surfaces starting from a coarse mesh, are also introduced together with an overview of classical schemes in the literature. Finally, in the context of numerical simulation of partial differential equation, the concepts of isogeometric analysis will be introduced.

B-spline functions

A basis spline function, better known as B-spline, is a piecewise polynomial function of a certain degree whose polynomial pieces join in order to achieve the desired global regularity.

There are several ways to define a B-spline, e.g. truncated power functions or blossoming; here we use the Cox-de Boor recurrence formula, which also gives an efficient algorithm for their implementation. We refer the reader to [Boo78; PT95; PBP02] for further details and proofs of the properties presented in this subsection.

Definitions and properties

Let t = {t 0 , . . . , t m } be a sequence of m+1 non-decreasing numbers, i.e. t i ≤ t i+1 , i = 0, . . . , m, called knot vector and whose elements are referred as knots. The i-th B-spline of degree d N i,d is defined as

N i,d (t) = t -t i t i+d -t i N i,d-1 (t) + t i+d+1 -t t i+d+1 -ti + 1 N i+1,d-1 (t), (1.1)
where

N i,0 (t) = 1 if t i ≤ t ≤ t i+1 , 0 otherwise.
The quotients in (1.1) may happen to be of the form 0/0, which we define to be 0. For brevity we often write N i,d instead of N i,d (t). Figures 1.1 and 1.2 present a set of quadratic and cubic B-spline basis obtained from an uniform and non-uniform knot vector, respectively. B-spline functions enjoy interesting properties:

• Local support: N i,d (t) = 0 if t / ∈ [t i , t i+d+1 ),
• Nonnegativity: N i,d (t) ≥ 0 ∀ i, d, t,

• Partition of unity: for an arbitrary knot span, i.e. the interval [t j , t j+1 ), we have

j ℓ=j-d N ℓ,d (t) = 1 ∀ t ∈ [t j , t j+1 ),
• In a given knot span [t j , t j+1 ) we have at most d + 1 non-zero B-spline basis functions, namely N j-d,d , . . . , N j,d ,

• Except for the case d = 0, N i,d (t) reaches exactly one maximum value,

• Derivatives formula: the k-th derivatives of a basis function is given by

N (k) i,d (t) = d t i+d -t i N (k-1) i,d-1 (t) - d t i+d+1 -t i+1 N (k-1) i+1,d-1 (t), (1.2) • At a knot N i,d (t) ∈ C d-µ
, where µ is the multiplicity of the knot.

B-spline curves

Let {P i } n i=0 be a set of n + 1 points in R 2 or R 3 called control points and {N i,d } the d-th degree B-spline functions defined on the knot vector t = {a, . . . , a d+1 , t d+1 , . . . , t m-d-1 , b, . . . , b d+1 } Figure 1.2. Cubic B-splines defined on the nonuniform knot vector t = {0, 0, 0, 0, 1, 2, 2, 3, 3, 3, 4, 5, 5, 5, 5}.

consisting of m + 1 knots. Unless stated otherwise, we assume that a = 0 and b = 1. A B-spline curve of degree d is defined by

C(t) = n i=0 P i N i,d (t), 0 ≤ t ≤ 1,
where the number of control points n + 1 is given by the relation n = md -1. The polygon formed by the {P i } is called the control polygon. Figure 1.3 shows an example of B-spline curves. The properties of B-spline curves, which will be now listed, follow directly from the properties of the basis functions N i,d presented in Section 1.1.1:

• In the case n = d and t = {0, . . . , 0, 1, . . . , 1}, the curve C(t) is a Bézier curve. This means that it can be rewritten as

C(t) = n i=0 P i B i d (t), with B i d (t) = d i t i (1 -t) d-i (1.3)
the Bernstein polynomials of degree d,

• Convex hull : the curve is contained in the convex hull of its control polygon,

• Endpoints interpolation: it results C(0) = P 0 , C(1) = P n ,

• Affine invariance: the affine transformation of a B-spline curve is obtained just transforming its control points,

• Local modification: moving a control point P i only affects C(t) in the interval [t i , t i+d+1 ),

• The control polygon provides a piecewise linear approximation of the curve,

• Variation diminishing: the curve crosses any plane at most as many times as it intersects its control polygon,

• Derivatives formula: using (1.2), the k-th derivatives of a B-spline curve can be computed by the formula

C (k) (t) = n-k i=0 P (k) i N i,d-k , where P (k) i =    P i k=0, d-k+1 t i+d+1 -t i+k P (k-1) i+1 -P (k-1) i k>0,
• C(t) ∈ C d-µ at a knot of multiplicity µ.

Tensor-product B-spline surfaces

The simplest bivariate B-spline construction is given by the tensor-product extension. It is obtained by taking a net of control points {P i,j } n,r i=0,j=0 , defining the control net, and the product of the univariate B-spline functions {N i,d } n i=0 , of degree d over the knot vector t and degree q over the knot vector s, i.e. {N j,q } r j=0 . More precisely, the B-spline surface S(t, s) of bidegree (d, q) is given by

S(t, s) = n i=0 r j=0 P i,j N i,d (t)N j,q (s), 0 ≤ t, s, ≤ 1,
with t and s defined similarly as in Section 1.1.2.

Figure 1.4 illustrates an example of bivariate tensor-product B-spline basis functions from which the B-spline surface in Figure 1.5 is obtained. It benefits from the same properties of B-spline curves, here shown for completeness: • if n = d, r = q, t = s = {0, . . . , 0, 1, . . . , 1}, then S(t, s) is a Bézier surface,

• Convex hull : if (t, s) ∈ [t i , t i+1 ) × [s j , s j+1
), then S(t, s) belongs to the convex hull of the points P i,j , id ≤ i ≤ i and jq ≤ j ≤ j,

• Corner points interpolation: it results S(0, 0) = P 0,0 , S(1, 0) = P n,0 , S(0, 1) = P 0,r , S(1, 1) = P n,r ,

• Affine invariance: an affine transformation is applied to the B-spline surface by applying it to the control net,

• Local modification: if the control point P i,j is moved, S(t, s) only changes in the rectangle

[t i , t i+d+1 ) × [s j , s j+q+1 ),
• The control net forms a piecewise planar approximation of the surface,

• Derivatives formula: the general formula returning the (k, l)-th derivative of a B-spline surface is given by

∂ k+l ∂ k t∂ l s S(t, s) = n-k i=0 r-l j=0 P (k,l) i,j N i,d-k (t)N j,q-l (s), with P (k,l) i,j = (q-l+1) P (k,l-1) i,j+1 -P (k,l-1) i,j s j+q+1 -s j+l ,
• S(t, s) is dµ (respectively qµ) times differentiable in the direction t (s) at a knot t (s) of multiplicity µ.

Smooth joints between B-spline patches

Let us consider here two adjacent B-spline surfaces of bidegree (d, q) We say that the patches S L (u, v) and S R (u, v) join C r continuous if the points

S L (u, v) = n,m i,j=0 P L i,j N i,d (u)N j,q (v) and S R (u, v) = n,m i,j=0 P R i,j N i,d (u)N j,q (v).
P L n-r,l , . . . , P L n,l = P R l,0 , . . . , P R l,r (1.4)
constitute the control polygon of some Bézier curve of degree r ∀ l = 0, . . . , m. In particular, C 0 continuity just requires that P L m,l = P R 0,l for all l = 0, . . . , m. Tensor-product schemes are very easy to use in presence of a regular rectangular net. However, in real life we have to deal with general quadrangular shapes. In this setting imposing higher smoothness conditions between adjacent patches e.g. C 1 smoothness is not straightforward; in fact, the C 1 relations between control points across an edge form a linear system of 3(d + 1) equations whose rank depends on the geometry of the quadrilaterals (see [START_REF] Bercovier | Smooth Bézier Surfaces over Unstructured Quadrilateral Meshes[END_REF][START_REF] Kapl | Isogeometric analysis with geometrically continuous functions on two-patch geometries[END_REF] for more details). A simpler type of continuity which can be used in these situations is the so called tangent plane continuity or G 1 continuity. In words, we say that two patches join G 1 smoothly if they can be reparametrized so that their first derivatives are identical along a common boundary curve; in Section 1.3.2 we will properly define this geometric continuity in the context of splines over topological quad meshes.

NURBS

Non Uniform Rational B-spline, or NURBS, curves and surfaces, are an extension of the concepts in Section 1.1.2 and Section 1.1.3 which allow more flexible. Their more interesting feature is the ability to exactly reproduce section of conics, which was not possible with B-spline objects.

Formally, a NURBS curve of degree d is defined by

C(t) = n i=0 w i P i N i,d (t) n i=0 w i N i,d (t) , 0 ≤ t ≤ 1,
where the basis functions {N i,d } n i=0 are defined over the knot vector t an in Section 1.1.2 and {w i } n i=0 , w i ∈ R, are the weights associated to the control points {P i } n i=0 . The higher the value of the weight, the closer the curve will be to the corresponding control point. In case of negative weights the curve results moved away from its associated control point.

Analogously, we define a NURBS surface of bidegree (d, q) on the knot vectors t and s as Section 1.1.2 by

S(t, s) = n i=0 r j=0 w i,j P i,j N i,d (t)N j,q (s) n i=0 r j=0 w i,j N i,d (t)N j,q (s) , 0 ≤ t, s ≤ 1.
Both NURBS curves and surfaces enjoy similar properties and formulas to the B-spline case; see, e.g., [START_REF] Piegl | The NURBS Book[END_REF] for a more in-depth look at the topic.

Fundamental algorithms

When dealing with B-splines it may be convenient, in some computation, to reformulate their algebraic expressions but without affecting the actual shapes. Two tools providing this kind of transformation are given by the knot insertion algorithm and degree elevation algorithm.

Knot insertion

Let C(t) = n i=0 P i N i,d (t) be a B-spline curve defined on the knot vector t = {t 0 , . . . , t m }. Suppose we want to insert a knot t ∈ [t k , t k+1 ) in the vector t defining the new knot vector t = { t0 = t 0 , . . . , tk = t k , tk+1 = t, tk+2 = t k+1 , . . . , tm+1 = t m }; if S t and St identify the vector spaces of B-spline curves on t and t respectively, clearly results that S t ⊂ St and dim (St) = dim (S t ) + 1. Hence C(t) has a representation on t of the form

C(t) = n+1 i=0 Q i Ni,d (t),
where Ni,d are the B-spline basis functions on t. The new control points Q i can be calculated either by direct computation, i.e. solving the linear system

n i=0 P i N i,d (t) = n+1 i=0 Q i Ni,d (t)
in the unknowns Q i , or using geometric identities given by the properties in Section 1.1.2. Both ways lead to the solution

Q i = α i P i + (1 -α i )P i-1 , where        1 i ≤ k -d, t-t i t i+d -t i k -d + 1 ≤ i ≤ k, 0 i ≥ k + 1.
When dealing with curves with no inner knots, the repeated knot insertion algorithm coincides with the de Casteljau algorithm for Bézier curves. In case of surfaces, the same algorithm can be applied separately first to the control points along the direction t and then to control points along the direction s, or vice versa.

Degree elevation

The goal of this algorithm is to identify a B-spline curve of degree d as an element of a curve space of higher degree d + 1. As for the knot insertion algorithm, this procedure do not affects the geometry of the curve. Now let C(t) = n i=0 P i N i,d (t) be a B-spline curve of degree d on the knot vector t. Since C(t) is a piecewise polynomial curve, it is possible to elevate its degree to d + 1 just by elevating the degree of its basis functions. Thus, there must exist control points Q i and a knot vector t such that

n i=0 P i N i,d (t) = n i=0 Q i N i,d+1 (t).
If the vector t has the form t = {0, . . . , 0

d+1 , t 1 , . . . , t 1 µ 1 , . . . , t m , . . . , t m µm , 1, . . . , 1 d+1 }, with µ 1 , . . . , µ m
the multiplicity of inner knots, in order to maintain the same regularity as the initial curve at the knots it must result that n = n+m+1 and t = {0, . . . , 0

d+2 , t 1 , . . . , t 1 µ 1 +1 , . . . , t m , . . . , t m µm+1 , 1, . . . , 1 d+2 }.
Again, the new control points can be computed by using the properties of B-spline curves in Section 1.1.2 thus obtaining

Q 1 = (1 -α i )P i + α i P i-1 , with α i = i d + 1 , i = 0, . . . , d + 1.
To degree elevate a B-spline surface it is sufficient to apply the algorithm separately to the control points of the net aligned to the two directions t and s.

Beyond tensor-product

The tensor-product structure is easy to deal with, but in real applications domains are not always rectangular. In this section we will introduce spline constructions over more general domains, that is triangulations and topological quad meshes. Moreover, we introduce the concept of geometric continuity, which is a widely used kind of regularity in multipatch approaches.

Splines on triangulations

We start this subsection introducing a local reference system, namely the barycentric coordinates system, that is much more useful than the usual Cartesian coordinates when dealing with triangles. Later, they will be used to define the so-called Triangular Bernstein polynomials and Bézier surfaces on a triangle. The notions presented in this subsection rely on [BS16; LS07].

Barycentric coordinates

Suppose △ is a nondegenerate triangle in R 2 with vertices p 1 , p 2 , p 3 . It will be useful to write △ = ⟨p 1 , p 2 , p 3 ⟩. Every point x ∈ R 2 has a unique representation in terms of its barycentric coordinates (β 1 , β 2 , β 3 ) with respect to △ such that If the point x lies inside the triangle △, then its barycentric coordinates are all positive; otherwise, some of them will result to be negative. Let △ 1 be the subtriangle △ 1 = ⟨x, p 2 , p 3 ⟩ and similarly △ 2 = ⟨x, p 3 , p 1 ⟩, △ 3 = ⟨x, p 1 , p 2 ⟩ (cf. Figure 1.6); it results that

x = β 1 p 1 + β 2 p 2 + β 3 p 3 , with β 1 + β 2 + β 3 = 1. p 1 p 2 p 3 x 3 1 2 (a) + + + + + - + -- -+ - --+ -+ + + - + p 1 p 2 p 3 (b)
β i = |△ i | |△| , i = 1, 2, 3,
where |△| indicates the signed area of the triangle △. Moreover, if x 1 and x 2 are two points in R 2 and △ their reference triangle, the barycentric directional coordinates (δ 1 , δ 2 , δ 3 ) of the vector x 2x 1 are defined as the difference of their corresponding barycentric coordinates.

Triangular Bernstein polynomials

Given a triangle △, we define the triangular Bernstein polynomials of degree d as

B i,j,k d (x) = d! i!j!k! β i 1 β j 2 β k 3 , ∀ i + j + k = d,
where (β 1 , β 2 , β 3 ) are the barycentric coordinates of the point x w.r.t. △. Figure 1.7 presents an example of these polynomials in the case d = 3. It is important to note that along each edge of a triangle, up to reordering of the indices, the triangular Bernstein polynomial coincides with the univariate Bernstein polynomial (1.3). They also possess properties which are similar to the univariate case, see Section 1.1.1. These properties are:

• Nonnegativity:

B i,j,k d (x) ≥ 0, x ∈ △,
• Partition of unity:

i+j+k=d B i,j,k d (x) = 1, ∀ x ∈ R 2 ,
• Recurrence relation: with (β 1 , β 2 , β 3 ) the barycentric coordinates of the point x w.r.t △ and B i,j,k d (x) = 0 in case of negative indices i, j, k,

B i,j,k d (x) = β 1 B i-1,j,k d-1 (x) + β 2 B i,j-1,k d-1 (x) + β 3 B i,j,k-1 d-1 (x),
• Degree elevation:

B i,j,k d (x) = i + 1 d + 1 B i-1,j,k d-1 (x) + j + 1 d + 1 B i,j-1,k d-1 (x) + k + 1 d + 1 B i,j,k-1 d-1 (x), (1.5)
• Directional derivatives: if u is a vector in R 2 with barycentric directional derivatives (δ 1 , δ 2 , δ 3 ) w.r.t. △, it holds that

D u B i,j,k d (x) = d δ 1 B i-1,j,k d-1 (x) + δ 2 B i,j-1,k d-1 (x) + δ 3 B i,j,k-1 d-1 (x) .
• Unimodal behavior : for i ≥ 1, j ≥ 1, k ≥ 1 it results that D u B i,j,k d (x) = 0 for any direction u and x ∈ △ if and only if

x = ip 1 + jp 2 + kp 3 d =: ξ i,j,k .
The points ξ i,j,k are called Greville points or domain points. The domain points are uniformly distributed over the triangle, and they can be triangulated by linking each pair of domain points

ξ i 1 ,j 1 ,k 1 and ξ i 2 ,j 2 ,k 2 such that |i 1 -i 2 | + |j 1 -j 2 | + |k 1 -k 2 | = 2. (1.6) 1.2.1.3 Bézier surfaces Let P i,j,k ∈ R 3 , i + j + k = d. A Bézier surface is defined as S(x) = i+j+k=d P i,j,k B i,j,k d (x),
where P i,j,k are its control points. Connecting these points as in (1.6) we obtain the control net of S(x). Figure 1.8 shows an example of cubic Bézier surface. They enjoy analogous properties to those of the tensor-product case: • Convex hull : S(x) is contained in the convex hull of its control points P i,j,k ,

• Corner points interpolation: S(1, 0, 0) = P 3,0,0 , S(0, 1, 0) = P 0,3,0 , S(0, 0, 1) = P 0,0,3 ,

• Directional derivatives formula: let u ∈ R 2 be a vector. The m-th directional derivative is given by

D (m) u S(x) = d! (d -m)! i+j+k=d-m ∆ m (P i,j,k ) B i,j,k d-m (x), m ≤ d,
where ∆ m (P i,j,k ) := P

[d-m]
i,j,k are the quantities obtained applying m steps of the triangular de Casteljau algorithm to the control points of S(x), i.e.

P [l-1] i,j,k = β 1 P [l] i+1,j,k + β 2 P [l] i,j+2,k + β 3 P [l]
i,j,k+1 , and (β 1 , β 2 , β 3 ) the barycentric coordinates of x.

• Degree elevation: by using (1.5), we obtain

i+j+k=d P i,j,k B i,j,k d (x) = i+j+k=d+1 P i,j,k B i,j,k d+1 (x),
where

P i,j,k = i d + 1 P i-1,j,k + j d + 1 P i,j-1,k + k d + 1 P i,j,k-1 and P i,j,k = 0 if i < 0 or j < 0 or k < 0.

A brief summary of simplex splines

Simplex splines represent an extension to the multivariate case of classical B-spline functions.

In fact, they benefit from the same malleability and similar properties to the univariate case. This subsection is based on [START_REF] Lyche | Construction of C 2 Cubic Splines on Arbitrary Triangulations[END_REF], while we point to, e.g. [START_REF] Prautzsch | Bézier and B-Spline Techniques[END_REF], for the proofs of the properties listed here.

To define formally a m-variate simplex spline of a given degree d, let Ξ = {ξ 1 , . . . , ξ n+1 } be a sequence of possibly repeated points (called knots) in R m , with n = d + m. We assume that the convex hull of Ξ is nondegenerate, i.e. vol m (⟨Ξ⟩) > 0. Consider now any simplex in R n , σ = { ξ1 , . . . , ξn+1 } with vol n (σ) > 0, satisfying the following projection property (cf. Figure 1.9):

π : R n -→ R m , π( ξi ) = ξ i , i = 1, . . . , n + 1.
Hence, the (normalized) simplex spline M Ξ can be defined as They benefit from properties analogous to the univariate case, such as

M Ξ : R m -→ R, M Ξ (x) = vol n-m (σ ∩ π -1 (x)) vol n (σ) .
• Local support: M Ξ has support ⟨Ξ⟩, • Nonnegativity: M Ξ (x) ≥ 0, ∀ x ∈ ⟨Ξ⟩,
• Normalization: M Ξ has unit integral,

• Knot dependence: M Ξ only depends on Ξ; in particular, the ordering of the knots is not affecting the simplex spline. It is also independent of the choice of the simplex σ,

• Derivatives formula: if u ∈ R m s.t. i a i ξ i = u, i a i = 0, then D u M Ξ = (d + m) d+m+1 i=1 a i M Ξ\ξ i , • Recurrence formula: for any x ∈ R m s.t. i b i ξ i = x, i b i = 1, it results M Ξ (x) = d + m d d+m+1 i=1 b i M Ξ\ξ i (x),

Splines on meshes

as well as the knot insertion algorithm (similarly to Section 1.1.5): for any y ∈ R m s.t.

i c i ξ i = y, i c i = 1, we have M Ξ = d+m+1 i=1 c i M Ξ∪y\ξ i .
If m = 1, then M Ξ is the normalized univariate B-spline of degree d with knot vector Ξ. If m = 2, i.e. the bivariate case, we define knot lines the lines in the complete graph of Ξ, which provide a polygonal partition of ⟨Ξ⟩. Furthermore, in this setting the simplex spline M Ξ is a polynomial of degree d = #Ξ -3 in each region of this partition, and across a knot line M Ξ ∈ C d+1-µ , with µ the number of knots that lie on the considered knot line (counted with multiplicity).

Splines on meshes

This section is devoted to present the construction of smooth spline objects over quad and triangular meshes: in particular, we present the constraints defining C 1 continuity across a triangular planar mesh and the more general G 1 smoothness for splines over quadrilateral meshes.

Splines on triangular planar meshes

Similarly to the case of rectangular B-spline patches, we present here the C r continuity conditions between adjacent triangular patches. Let △ L = ⟨p 1 , p 2 , p 3 ⟩ and △ R = ⟨p 4 , p 3 , p 2 ⟩ be two adjacent triangles on which the polynomials

S L (x) = i+j+k=d P L i,j,k B i,j,k d (x), S R (x) = i+j+k=d P R i,j,k B i,j,k d (x)
are defined. Then S L (x) and S R (x) join C r -smooth if and only if

P R m,j,k = i+ j+ k=m P L i,k+ j,j+ k B i, j, k m (p 4 ), j + k = d -m, m = 0, . . . , r.
(1.7)

The relation (1.7) is the triangular variant of (1.4) in the tensor-product setting. C 0 continuity only requires that P L i,j,0 = P R i,j,0 , for all i + j = d, while, in addiction, C 1 smoothness implies that the pair of triangles ⟨P L i,j,1 , P L i+1,j,0 , P L i,j+1,0 ⟩, ⟨P R i,j,1 , P R i+1,j,0 , P R i,j+1,0 ⟩, for each i + j = d -1, are coplanar. Even when dealing with triangulations imposing (global) C r continuity between adjacent triangles is not straightforward and resorting to geometric continuity can be a valid alternative (see [START_REF] Prautzsch | Bézier and B-Spline Techniques[END_REF]).

We conclude the subsection introducing the definition of spline space over a triangulation. Let T be a triangulation of a polygonal domain Ω ∈ R 2 , i.e. a partition of Ω consisting of non-overlapping triangles. The spline space of degree d on T with C r continuity is defined as

S r d (T ) := s ∈ C r (Ω) : s △ i ∈ P d , △ i ∈ T ,
where we identify with P d the space of polynomials in two variables and total degree d.

G 1 -smooth splines over topological quad meshes

Here we introduce the concept of splines over topological quad meshes, with a particular focus on geometrically smooth spline spaces. In order to define these objects, the concepts of quadrilateral mesh, together with its features, as well as transition map and gluing data functions will be introduced. We start introducing the concept of topological surface M (or mesh) as a collection of faces and edges together with their adjacency relations. An edge can be glued with at most one other edge and it cannot be glued with itself. Here we are interested in meshes which contain only quadrilateral faces and we shall refer to such a structure as quadrilateral mesh. An essential notion is the valence of a vertex, denoted with N , which is defined as the number of faces containing the vertex. An interior vertex is said to be regular (RV in short) if its valence is equal to 4, and it is called irregular (or extraordinary -EV) otherwise.

Let σ 0 , σ 1 be two faces of the mesh M sharing the common edge e; we can associate to each face a coordinate system, namely (u 0 , v 0 ) and (u 1 , v 1 ) on σ 0 and σ 1 , respectively (see Figure 1.10). Now, we can therefore define the (first order) transition map ϕ σ 0 ,σ 1 from the face σ 1 to σ 0 as the function

ϕ σ 0 ,σ 1 : R 2 -→ R 2 , (u 1 , v 1 ) -→ (u 0 , v 0 ) = v 1 b N,N ′ (u 1 ) + O(v 2 1 ) u 1 + v 1 a N,N ′ (u 1 ) + O(v 2 1 ) , (1.8)
which characterizes the change of coordinates from σ 1 to σ 0 . Moreover, the functions in (1.8)

a : e -→ R, b : e -→ R,
are C 1 functions called gluing data and they only depend on the two valences N and N ′ of the vertices identifying the edge e. We have all the ingredient to introduce the definition of G r continuity across an edge shared by two patches. Consider a multipatch function f , i.e. a collections of functions f = (f σ ) σ∈M defined over the faces of a mesh M; it is said to be G r -smooth if for any two functions

f σ 0 , f σ 1 J r e (f σ 1 )(u 1 , v 1 ) = J r e (f σ 0 • ϕ σ 0 ,σ 1 )(u 1 , v 1 ), (1.9)
where J r e indicates the jet (or Taylor expansion) of order r along the common edge e.

In particular, for the case r = 1 we are mostly interested in, previous (1.9) leads to the following two relations for each u 1 ∈ [0, 1]:

   f 1 (u 1 , 0) = f 0 (0, u 1 ), ∂f 1 ∂v 1 (u 1 , 0) = b N,N ′ (u 1 ) ∂f 0 ∂u 0 (0, u 1 ) + a N,N ′ (u 1 ) ∂f 0 ∂v 0 (0, u 1 ), (1.10)
where

f 0 = f | σ 0 , f 1 = f | σ 1
are the restrictions of f on the faces σ 0 , σ 1 . Geometrically, (1.10) means that the functions f 0 and f 1 share the same tangent plane along their interface, from where the name tangent plane continuity.

Finally, the spline space of bidegree (d, d) with regularity G 1 on a topological quad mesh M is defined as

S 1 d (M) := (f σ ) σ∈M : f i = f σ i ∈ P d , f i , f j join G 1 across e = σ i ∩ σ j ∀ i, j .
We point to [BMV17; MVV16; Man23] for a deeper insight on geometrically smooth spline spaces.

Subdivision surfaces

In this subsection we shortly recall the very basics of Subdivision surfaces: the definition of mask and subdivision matrix will be introduced as well as the concept of limit surface. Also, some classical schemes will be recalled and analyzed in order to have a general idea on how they behave. For a complete and accurate description of subdivision schemes, see, for example, [PR08; AS10], on which this paragraph is founded.

Definitions and properties

Let us consider a mesh M. The idea behind subdivision is to obtain smooth objects via iterative refinement of an input mesh; different refinement strategies lead to different limit objects which differ from both smoothness and shape. Each subdivision scheme is defined by its masks (or stencils): in general, masks relate to the vertices, edges and faces of the mesh. A mask is determined by a set of normalized coefficient (called weights), responsible for the definition of the recursive rule, which is applied to the considered vertex of the mesh and its first and second neighborhood. We call first neighborhood of a vertex v the set of points directly connected to v by an edge and, when dealing with quadrilateral meshes, second neighborhood of v the set of points belonging to the patches sharing v and opposite to it. Figure 1.11 presents an example of subdivision masks with underlined colors for the two neighborhoods.

In order to apply these subdivision operators (i.e. the masks) we make use of the so-called subdivision matrix, which is a matrix whose columns contains the subdivision masks and its size depends on the mesh size. If p = p (0) is a initial vector of points and S = S (0) is its related subdivision matrix, one step of subdivision is obtained by p (1) = p (0) S (0) , which can be generalized to an arbitrary subdivision step as p (k+1) = p (k) S (0) or p (k+1) = p (0) S (k) , with S (k) the subdivision matrix concerning the k-th subdivision step. The repeated subdivision converges to a surface called limit surface, whose regularity properties can be deduced carrying out a spectral analysis on the subdivision matrix defining the scheme.

An overview on classical subdivision surfaces

The aim of this subsection is to give a visual idea of how different schemes act on meshes, generating different limit surfaces. We will consider four classical subdivision schemes for surfaces, namely the Catmull-Clark [CC78], Doo-Sabin [DS78], Loop [START_REF] Loop | Smooth Subdivision Surfaces Based on Triangles[END_REF] and Butterfly [START_REF] Dyn | A Butterfly Subdivision Scheme for Surface Interpolation with Tension Control[END_REF] schemes.

Catmull-Clark scheme

The Catmull-Clark scheme acts on quad meshes. It is an approximating scheme, i.e. the vertices of the initial mesh are not maintained in the subdivision process, and its limit surface is devised as a generalization of bicubic B-spline surfaces with C 2 regularity everywhere except at EVs, where it is C 1 -smooth. Figure 1.12 shows the behaviour of this scheme on a quad mesh. Moreover, a direct evaluation of the limit surface is possible by using the formulas provided in [START_REF] Stam | Exact Evaluation of Catmull-Clark Subdivision Surfaces at Arbitrary Parameter Values[END_REF].

Figure 1.12. Three iterations of the Catmull-Clark scheme starting from the input mesh on the left.

Doo-Sabin scheme

Similarly to the Catmull-Clark scheme, Doo-Sabin rules are defined for quad meshes. It is an approximating scheme and it is based on a generalization of biquadratic uniform B-splines producing limit surfaces which are C 2 regular everywhere except at EVs, where we have C 1 smoothness; it is interesting to note that this scheme produces N -gons around EVs of valence N (see Figure 1.13). The limit surface evaluation proposed in [START_REF] Stam | Exact Evaluation of Catmull-Clark Subdivision Surfaces at Arbitrary Parameter Values[END_REF] is applicable also for this scheme.

Loop scheme

The Loop method for subdivision surfaces is an approximating scheme developed for triangular meshes. It is based on a quartic box spline and the limit surfaces it creates are C 2 continuous in regular regions while around EVs the smoothness drops to C 1 (EVs in triangular meshes are vertices with valence N ̸ = 6); Figure 1.14 presents an example of these surfaces. An explicit evaluation of the Loop limit surfaces can be achieved from [START_REF] Stam | Evaluation of Loop Subdivision Surfaces[END_REF].

The Butterfly scheme

An example of interpolatory scheme is the Butterfly scheme: in fact, as can be noticed in Figure 1.15, the vertices of the initial triangular mesh are also vertices of the resulting subdivided surfaces. The limit surface it generates, as established in [START_REF] Zorin | A method for analysis of C 1 -continuity of subdivision surfaces[END_REF] several years after the publication of the method itself, is C 1 -smooth everywhere but around vertices with valences N = 3 and N ≥ 8 where it is only continuous.

Isogeometric analysis

This last subsection is devoted to the introduction of isogeometric analysis. IsoGeometric Analysis (IGA) is a highly efficient technique for solving PDEs numerically. Its basic idea, presented in [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF], which unifies the Finite Elements Method (FEM in short) approach and Computer Aided Geometric Design, is to use the same basis functions for both reproducing exactly the computational domain, and for the numerical approximation of the PDE; more precisely, it is involved in the computation of the weak solution of a PDE. We start recalling some basics on PDEs, like their strong and weak formulation, and then continue introducing the so-called isogeometric paradigm. The content of this subsection has been extracted from [CHB09; BS16; Qua16], to which we also refer for the proofs of the results we recall.

PDEs: strong and weak formulation

Partial differential equations are differential equations containing the derivatives (spatial and/or temporal) of an unknown function. In particular, denote by u the unknown function in the d + 1 independent variables x = (x 1 , . . . , x d ) and t, a generic PDE can be denoted as

F t,
-∆u = f in Ω, (1.12)
where Ω ⊂ R 2 is a bounded domain and we denote by

∆u = d i=1 ∂ 2 u ∂x 2 i
the Laplace operator or laplacian. In case f = 0, (1.12) is know as Laplace equation. In order to obtain a unique solution for (1.12), we need to add some extra conditions regarding the behaviour of the solution on the boundary of the domain ∂Ω; the two standard boundary conditions are assigning the value of the solution on ∂Ω (Dirichlet condition) or the value of the normal derivative of the solution along ∂Ω (Neumann condition). These choices are used to define the general problem, stated as follow: find u such that We derive now the so-called weak form of the Poisson's equation with homogeneous Dirichlet conditions, that is h D = 0; the derivation of the weak formulation for other PDEs with eventually different boundary conditions is analogous. To do this, let v be a test function, i.e. a sufficiently smooth function which vanishes on ∂Ω = Γ D ; later on we will properly define this function. Multiplying both members of (1.12) by v and integrating over Ω we get

       -∆u = f in Ω, u = h D on Γ D , ∂u ∂ n = h N on Γ N , where ∂Ω = Γ D ∪ Γ N ,
- Ω ∆uv dΩ = Ω f v dΩ.
(1.13) By using the vectorial identity div (ψ∇A) = ψdiv (∇A) + ∇A • ∇ψ, with A a vector field and ψ a scalar function, and the divergence or Gauss theorem

Ω div(A) dΩ = ∂Ω A • n dγ,
where dγ represents the infinitesimal boundary element of ∂Ω, (1.13) becomes

Ω ∇u • ∇v dΩ = Ω f v dΩ,
where we used the fact that v is null on the boundary of the domain. Therefore, we arrive at the following weak formulation of the Poisson's problem:

find u ∈ H 1 0 (Ω) : Ω ∇u • ∇v dΩ = Ω f v dΩ ∀v ∈ H 1 0 (Ω), (1.14)
where, in order to have feasible operations, we take f ∈ L 2 (Ω) and we set

H 1 (Ω) := φ ∈ L 2 (Ω) : ∂φ ∂x i ∈ L 2 (Ω), i = 1, 2 , H 1 0 (Ω) := φ ∈ H 1 (Ω) : φ ∂Ω = 0 .
In general, we can define the Hilbert space H k (Ω) as

H k (Ω) := φ ∈ L 2 (Ω) : D α φ ∈ L 2 (Ω) ∀α, |α| ≤ k
with associated seminorms and norms

|φ| H k (Ω) = |α|=k Ω (D α φ) 2 dΩ, ∥φ∥ H k (Ω) = k m=0 |φ| 2 H k (Ω) ,
where α is a multi-index and the derivatives must be understood as distributional derivatives (cf. [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]). It is interesting to notice that the problem (1.14) is equivalent to the following variational problem:

find u ∈ H 1 0 (Ω) : J(u) = inf v∈H 1 0 (Ω) 1 2 Ω |∇v| 2 dΩ - Ω f v dΩ.
To conclude, we introduce now a more compact formulation of the weak form (1.14). Let V be an Hilbert space (V = H 1 0 (Ω) in our case). We can define the bilinear form

a : V × V -→ R, a(u, v) = Ω ∇u • ∇v dΩ
and the linear functional

F : V -→ R, F (v) = Ω f v dΩ such that the problem (1.14) becomes find u ∈ V : a(u, v) = F (v) ∀v ∈ V. (1.15)
We need the followings: a bilinear form a :

V × V -→ R is said to be • continuous if ∃ M > 0 such that |a(u, v)| ≤ M ∥u∥ V ∥v∥ V ∀u, v ∈ V ; • coercive if ∃ α > 0 such that a(v, v) ≥ α∥v∥ 2 V ∀v ∈ V.
The existence and uniqueness of a solution for the problem (1.15) is ensured by the following [BS08, Theorem 2.7.7] Theorem 1.5.1 (Lax-Milgram). Let V be an Hilbert space, a(•, •) : V × V -→ R a continuous and coercive bilinear form, F (•) : V -→ R a continuous linear functional. Then exist a unique

u ∈ V such that a(u, v) = F (v) ∀v ∈ V.
When dealing with nonhomogeneous Dirichlet conditions, i.e. h D ̸ = 0, the solution can be decomposed as

u = u 0 + u D ,
where u 0 ∈ H 1 0 (Ω) is obtained from (1.14), and u D ∈ H 1 D (Ω), with

H 1 D (Ω) := φ ∈ H 1 (Ω) : φ ∂Ω = h D ,
is given.

Strong and weak formulation of some classical PDEs

The aim of this paragraph is to introduce some classical PDEs, together with their weak formulations, that we are going to utilize later on in this work; ambient spaces, domains of definitions and boundary conditions will be properly defined when necessary.

• Biharmonic equation

∆ 2 u = f weak -----⇀ ↽ ----- strong Ω ∆u • ∆v dΩ = Ω f v dΩ, (1.16)
where ∆ 2 represents the biharmonic operator defined as

∆ 2 u = ∆(∆u) = d i,j=1 ∂ 4 u ∂x 2 i ∂x 2 j .
• Linear time-independent shallow-water equation

h∆u + ∇h • ∇u + κ 2 u = 0 weak -----⇀ ↽ ----- strong Ω h∇u • ∇v dΩ + Ω v∇h • ∇u dΩ + κ 2 Ω uv dΩ = 0,
(1.17) with h the bathymetry function and κ 2 = σ 2 /g ∈ R, where g is the gravitational acceleration and σ the frequency of the oscillations.

• Heat equation

∂u ∂t = c 2 ∆u weak -----⇀ ↽ ----- strong c 2 Ω ∇u • ∇v dΩ = Ω ∂u ∂t v dΩ, (1.18)
with c ∈ R the conductivity coefficient.

Chapter 1 Preliminaries

Galerkin's method and isogeometric spaces

Galerkin's approach for IGA is based on the decomposition of the variational problem into smaller ones defined onto subspaces of the initial space. Formally, let V be an Hilbert space and V h , h ∈ R, be a one-parameter family of spaces such that

V h ⊂ V, dim(V h ) < ∞ ∀h > 0.
The Galerkin's problem (or approximate problem) for (1.15) can be stated as

find u h ∈ V h : a(u h , v h ) = F (v h ) ∀v h ∈ V h . (1.19) If {ϕ i } i∈I , I = {1, . . . , dim(V h )}, is a basis of V h ,
a(u h , ϕ i ) = F (ϕ i ), i ∈ I. (1.20) Therefore, since u h ∈ V h , u h (x) = j∈I u j ϕ j (x),
with u j , j ∈ I unknown coefficients, (1.20) becomes j∈I u j a(ϕ j , ϕ i ) = F (ϕ i ), i ∈ I.

(1.21)

Denoting with K the stiffness matrix with entries

k i,j = a(ϕ j , ϕ i ),
f the vector of components f i = F (ϕ i ) and with u the vector of the unknown coefficients u j , (1.21) is equivalent to the linear system

Ku = f .
For completeness, we also introduce the mass matrix M , whose elements are given by

m i,j = m(ϕ j , ϕ i ), with m(u, v) = Ω uv dΩ;
this matrix may appear in the discretization of some variational formulation, e.g. (1.17). The existence and uniqueness of a solution of the Galerkin's problem is again ensured by the Lax-Milgram theorem (1.5.1), since V h is an Hilbert space (being a closed subspace of V ) and the form a(•, •) and functional F (•) are the same as in (1.15). Now, let assume that our domain Ω is given through a spline transformation of the parametric domain Ω = [0, 1] 2 . This means that there exist a map G, called geometry map, that is

G : Ω -→ R m , m = 2, 3, G(ξ) = i,j∈J c i,j ψ i,j (ξ),
where ξ = (ξ 1 , ξ 2 ), ψ i,j : Ω -→ R are spline basis functions of a certain degree d spanning the space S d and c i,j ∈ R 2 are their control points on a suitable indices set J = {1, . . . , dim(S d )},

ξ 2 ξ 1 [0, 1] 2 ξ 2 ξ 1 [0, 1] 2 G (1) G (2)
Ω (1) Ω (2)

Figure 1.16. A multipatch domain Ω consisting of two patches Ω (1) , Ω (2) with their associated geometry mappings G (1) , G (2) .

such that Ω = G( Ω). The founding idea behind isogeometric analysis is that the basis used to exactly describe the geometry will also serve as basis for the computation of the numerical solution. This concept brings us to the definition of the isogeometric function space

V h = v h ∈ L 2 (Ω) : v h ∈ S d • G -1
on which we solve our approximate problem. For an exhaustive description on mathematics of IGA we refer the reader to [START_REF] Beirão Da Veiga | Mathematical analysis of variational isogeometric methods[END_REF].

PDEs on manifolds

It is common in applications to deal with physical domains which are not planar. In this case, we need to reformulate the PDE in terms of operators on manifolds such that the Laplace-Beltrami operator. Manifolds are usually described as multipatch domains, i.e. as the union of single patch Ω (k) , each of them defined through a different geometry map G (k) , with k ∈ I Ω in some index set. Figure 1.16 shows this idea when dealing with a shape consisting of 2 patches.

Let us consider the patch Ω (k) of the 2-manifold domain Ω ⊂ R 3 . The Jacobian J (k) of the map G (k) , defined as

J (k) : Ω -→ R 3×2 , ξ -→ J (k) (ξ), J (k) (ξ) = ∂G (k) ∂ξ 1 (ξ) ∂G (k) ∂ξ 2 (ξ) ,
is used to define the first fundamental form of the mapping G (k)

G (k) : Ω -→ R 2×2 , ξ -→ G (k) (ξ), G(ξ) = J (k) (ξ) T J (k) (ξ),
Chapter 1 Preliminaries together with its determinant g (k) , which is

g (k) : Ω -→ R, ξ -→ g (k) (ξ), g (k) (ξ) = det G (k) (ξ) . Let φ (k) ∈ C 2 (Ω (k) ) and φ (k) (ξ) = φ (k) • G (k) (ξ)
. Finally, we have all the tools to define the gradient operator on the patch Ω (k) of the manifold Ω,

∇ Ω (k) φ (k) (x) = J (k) (ξ) G (k) -1 (ξ) ∇ φ (k) (ξ) • G (k) -1 (ξ), x ∈ Ω (k) , (1.22)
with ∇ representing the gradient operator in the parametric space, and the Laplace-Beltrami operator

∆ Ω (k) φ (k) (x) = 1 g (k) (ξ) ∇ • g (k) (ξ) G (k) -1 (ξ) ∇ φ (k) (ξ) • G (k) -1 (ξ), x ∈ Ω (k) .
(1.23) Hence, we can reformulate both strong and weak formulations of the PDEs (1.14) (1.16), (1.17) and (1.18) just replacing the gradient and laplacian operators with, respectively, the gradient operator on the manifold (1.22) and the Laplace-Beltrami operator (1.23). For further details on the description of PDEs on manifolds we refer to [START_REF] Dedè | Isogeometric Analysis for second order Partial Differential Equations on surfaces[END_REF][START_REF] Farahat | Isogeometric analysis with C 1 -smooth functions over multi-patch surfaces[END_REF].

Chapter 2

A G 1 approximation of Catmull-Clark surfaces Catmull-Clark subdivision scheme is undoubtedly the most famous and used scheme to obtain smooth surfaces from quad meshes. However, the limit surfaces that it generates may not be smooth enough for applications, e.g. in computer-aided design. The purpose of this chapter is to present a new, smooth and completely explicit construction of a globally G 1 family of Bézier surfaces, defined by smoothing masks, approximating the Catmull-Clark limit surface.

Starting from the C 0 scheme of Loop and Schaefer [START_REF] Loop | Approximating Catmull-Clark subdivision surfaces with bicubic patches[END_REF], we impose G 1 regularity around EVs by means of quadratic gluing data depending solely on their valences; this procedure leads to a system of equations of which we explore the space of solutions in order to achieve explicit formulas for the masks defining the patches' control points. Finally, we come up with four different schemes. In order to understand which scheme provides the most regular surface, we conduct isophotes and curvature analysis on an extensive benchmark, both visually and numerically. This chapter is based on [START_REF] Marsala | G 1 -Smooth biquintic approximation of Catmull-Clark subdivision surfaces[END_REF].

2.1 Bicubic Approximate Catmull-Clark (ACC 3 ) and degree elevated (ACC 5 )

The bicubic approximate Catmull-Clark scheme of [START_REF] Loop | Approximating Catmull-Clark subdivision surfaces with bicubic patches[END_REF] is a smoothing scheme reconstructing a bicubic surface starting from a quad-mesh which approximates the limit surface of the classical Catmull-Clark subdivision scheme; the application of this scheme to a mesh returns, on each face, 16 control points defining a bicubic Bézier patch. The surface obtained from ACC 3 is C 2 everywhere except in a neighborhood of the EVs, in which it is only C 0 ; moreover, the vertices generated by ACC 3 interpolate the Catmull-Clark limit surface. Usually, similarly to subdivision schemes, a smoothing scheme is defined through its smoothing matrix S = (M

(k) i,j ) as an operator such that b = p S, b (k) i,j = ℓ p ℓ M (k) i,j ℓ
, where p = (p 1 , . . . , p m ) ∈ R 3×m is the set of mesh vertices and b = (b (k) ), where b (k) is the matrix containing the control points of the k-th patch (see Figure 2.2 for their labeling). The smoothing matrix contains in its columns the different masks M defining the scheme and in this construction the weights of a smoothing (or Bézier) mask can be, possibly, negative. We will refer to smoothing masks using the notation M=[ω • , ω, ω ′ ] T , where ω • is the weight for the vertex point and ω = {ω i } N i=1 and ω ′ = {ω ′ i } N i=1 are the weights for its first and second neighborhoods, respectively. Lastly, in this construction we will consider meshes with isolated EVs, i.e. surrounded by regular vertices; this assumption is not restrictive since any linked EVs can be decoupled with one subdivision of the mesh. To recover regularity around these EVs without affecting the interpolatory property and C 1 (or C 2 ) smoothness of the surrounding regular vertices, linear gluing data and cubic or even quartic patches do not suffice. Indeed the G 1 relations for Bézier patches of bidegree (3, 3) or (4, 4) propagate up to control points neighboring the regular vertices, affecting then their properties. To avoid these effects, we will work with gluing data functions of degree two and biquintic patches, defined by 36 Bézier control points each. As a starting point of our construction we use the degree elevated masks ACC 5 , which are obtained after applying twice the degree elevation algorithm (Section 1.1.5.2) to ACC 3 masks; we will use underlined symbols to denote ACC 5 masks i.e. M = [ω • , ω, ω ′ ] T . For example, the ACC 5 masks M 2,2 and M 4,0 returning the Bézier points b 2,2 and b 4,0 , respectively, are given by

N N +5 2 N +5 1 N +5 2 N +5 (a) N N +5 2 N +5 1 2(N +5) 1 N +5 1 N +5 1 2(N +5) (b) N N +5 4 N (N +5) 1 N (N +5) 4 N (N +5) 4 N (N +5) 1 N (N +5) 4 N (N +5) 1 N (N +5) 4 N (N +5) 1 N (N +5) (c)
M 2,2 = 1 100 M 0,0 + 6 M 1,0 + 3 M 2,0 + 6 M 0,1 + 36 M 1,1 + 18 M 2,1 + 3 M 0,2 + 18 M 1,2 + 9 M 2,2 , M 4,0 = 1 5 3 M 2,0 + 2 M 3,0 ,
where M i,j are the ACC 3 bicubic masks (see Figure 2.1). Note that the degree elevation is directly applied on the considered masks to obtain the new weights. The masks for Bézier points can be written in vectorial form; for instance, the corner point mask M 0,0 is given as follows:

M 0,0 ∈ R 2N +1 , with M 0,0 = 1 N (N + 5) [N 2 , 4, . . . , 4, 1, . . . , 1] T = 1 N (N + 5) [N 2 , 4, 1] T ,
(2.1) where the first coordinate is the weight for the EV, followed by the weights for the vertices of the first neighborhood and then the second neighborhood; we will adopt this notation throughout the thesis.

G 1 constraints on Bézier patches

As we said in the previous section, our aim is to construct new smoothing masks defining Bézier patches that join G 1 at extraordinary regions. This can be obtained by the description in Section 1.3.2 as follow: let [a N,N ′ , b N,N ′ ] be quadratic symmetric gluing data functions [BH14; HBC08; Hah89]

a N,N ′ (u) = a 0 B 0 2 (u) + a 2 B 2 2 (u), where a 0 = 2 cos 2π N , a 2 = 2 cos 2π N ′ , b N,N ′ (u) = -1, (2.2)
with B i d the univariate Bernstein polynomials in (1.3) defined along the edge e shared by two adjacent patches. Since our construction is developed in the parametric domain, we take e = [0, 1]; moreover, in the cases we study we will assume that N ′ = 4 so that a 2 = 0, that is an EV is surrounded by regular vertices. Now, if we specialize (1.10) defining the G 1 constraints between two adjacent patches with f 0 and f 1 being two biquintic Bézier surfaces and using (2.2), we obtain

             5 i=0 b (1) i,0 B i 5 (u) = 5 i=0 b (0) 0,i B i 5 (u), 5 i=0 b (1) i,1 -b (1) i,0 + b (0) 1,i -b (0) 0,i B i 5 (u) = a 0 B 0 2 (u) 4 i=0 b (0) 0,i+1 -b (0) 0,i B i 4 (u) . (2.3) Calling with M (k)
i,j the smoothing mask returning the control point b

(k) i,j belonging to the patch k, solving (2.3) with respect to M (1) 0,1 , M (1) 1,1 , M (1) 2,1 , M (1) 3,1 , M (1) 4,1 , M (1) 5,0 , M (1)
5,1 , after some algebraic reductions we obtain the following relations to be satisfied:

M (1) 0,1 + M (0) 1,0 = ā0 M (1) 0,0 + a 0 M (1) 1,0 ,
(2.4)

5(M

(1)

1,1 + M (0) 1,1 ) = a 0 M (1) 0,0 + 5ā 0 M (1) 1,0 + 4a 0 M (1) 2,0 , (2.5
)

10(M (1) 2,1 + M (0) 1,2 ) = -a 0 M (1) 0,0 + 5a 0 M (1) 1,0 + 10ā 0 M (1) 2,0 + 6a 0 M (1) 3,0 , (2.6
)

10(M (1) 3,1 + M (0) 1,3 ) = a 0 M (1) 0,0 -5a 0 M (1) 1,0 + 10a 0 M (1) 2,0 + 10ā 0 M (1) 3,0 + 4a 0 M (1) 4,0 , (2.7)
M

(1)

4,1 + M (0) 1,4 = 2M (1) 4,0 , (2.8)
M

(1)

5,1 + M (0)
1,5 = 2M

(1) 5,0 , (2.9)

10(M (1) 3,0 -M (1) 2,0 ) = M (1) 0,0 -5M (1) 1,0 + 5M (1) 4,0 -M (1) 5,0 , (2.10)
with ā0 = 2a 0 . Masks verifying (2.4) to (2.10) generate a family of G 1 smoothing schemes. From now on the dependency of the mask from the patch will be specified only when it is different from patch 1. Particularly interesting is (2.10): if we use quadratic gluing data, the degree of the left-hand side of (1.10) is ≤ 5 while the right-hand side has degree 4 + 2 = 6. (2.10) reflects the constraint that the left-hand side has degree ≤ 5: in fact it is saying that even if we are working with biquintic patches, along the edge we have a quartic curve. Moreover, we notice that (2.8) and (2.9) are always satisfied by ACC 5 masks [LS08, Figure 3] : since we are in the vicinity of a regular vertex, we know that the ACC 5 surface is C 2 there; it means that the first derivatives along opposite directions across the regular vertex are the same, i.e. M (0)

1,5 -M (0) 0,5 = -(M (1) 5,1 -M (1) 5,0 ) as M (0) 0,5 = M (1)
5,0 , which is in fact (2.9). Then, since the limit surface is also at least C 1 across regular vertices (or C 2 in some cases), we have, additionally, equal second order derivatives along opposite directions i.e. M (0)

1,4 -M (0) 0,4 = -(M (1) 4,1 -M (1) 4,0 ) that is exactly (2.8) since M (0) 0,4 = M (1) 4,0 .

Explicit Bézier masks derivation

In this section we are going to solve the system (2.4) to (2.10) to obtain masks for the Bézier patches ensuring G 1 regularity of the resulting surface. The process is done incrementally with respect to the order of the Bézier mask where the order refers to the sum of the subscripts of the Bézier mask. In particular, for Bézier masks of order zero, i.e. vertex points, we take M 0,0 = M 0,0 as in (2.1) in order to not alterate the interpolating property we have at the vertices. For Bézier points of order one we will solve (2.4); since the latter involves Bézier points on different patches all around the EV, (2.4) translates into a circulant system to be solved. (2.5) presents constraints regarding masks of second order and it can be solved in two different ways; if we want to solve it with respect to M 1,1 , since it is related to points around the EV in different patches similarly to (2.4), it translates into a circulant system. On the other hand, if we want to solve it with respect to the mask M 2,0 it can be directly obtained from ACC 5 masks. Then, (2.6) and (2.7) regarding Bézier points of order three and four, respectively, can be solved by direct computation in either a symmetric or asymmetric approach. Moreover, the masks M 4,0 ,

M 5,0 , M 4,1 , M 5,1 , M 2,2 , M 3,2 , M 4,2 , M 5,2 , M 2,3 , M 3,3 , M 4,3 , M 5,3 as well as M 0,4 , M 1,4 , M 2,4 , M 3,4 , M 4,4 , M 5,4 , M 0,5 , M 1,5 , M 2,5 , M 3,5 , M 4,5 , M 5,5 are equal to those of ACC 5 .
Throughout the construction we assume that the first neighborhood of an EV is composed solely by regular vertices; this assumption is not restrictive, since it is satisfied by applying at most one step of any 4-split algorithm. Moreover, biquintic Bézier points will be computed only in the faces adjacent to EVs, otherwise we will use the standard bicubic patch of ACC 3 . First we focus on inner EVs and postpone the discussion of boundaries until Section 2.3.4.

Bézier masks of order one

We compute first the masks M 1,0 = M

1,0 and M (0) 0,1 using (2.4) as follows. If the valence of the EV is N = 3, (2.4) is automatically satisfied by ACC 5 masks. Then, we do not need to determine new masks for this relation. Otherwise, for valences N ≥ 5 we have to compute new masks verifying (2.4); given the mask for the vertex point M 0,0 = M 0,0 = [α • , α, α ′ ] T , we shall compute the mask M 1,0 = [β • , β, β ′ ] T . We observe that the masks M 0,1 = M

(1) 0,1 and M (0) 1,0 can be obtained by permuting the weights of M 1,0 . In particular, let C = diag(1, C, C) be a (2N + 1) × (2N + 1) block diagonal matrix where

C =        0 1 0 • • • 0 1 . . . 1 1 0 • • • 0        = Circ(0, 1, 0, . . . , 0) (2.11) is a N × N circulant matrix; it follows that M (0) 1,0 = C M 1,0 , M 0,1 = C -1 M 1,0 , where C -1 = diag(1, C -1 , C -1 ) and C -1 = C N -1
. Therefore (2.4) translates into the following system to be solved:

(-a 0 I + C + C N -1 )M 1,0 = (2 -a 0 )M 0,0 , (2.12)
where I is the (2N + 1) × (2N + 1) identity matrix; Figure 2.3-(a) shows the position of the control points involved in this construction. Our investigation towards explicit solutions leads to the following Proposition 2.3.1. Let A = -a 0 I + C + C N -1 be as in (2.12). The matrix A is symmetric and singular.

Proof. From the particular shape of A, its symmetry is straightforward since

C T = C -1 = C N -1
. Now, the matrix A can be decomposed as

A = -a 0 I + C + C N -1 = C N -1 (C 2 -a 0 C + I) = C N -1 (C -z 1 I)(C -z N -1 I) ,
where z 1 , z N -1 are the roots of the polynomial x 2a 0 x + 1, with a 0 = 2 cos (θ) , θ = 2π/N . If i = √ -1 is the imaginary unit, they are obtained as

z 1 = e iθ , z N -1 = e -iθ ,
i.e. they are the N -th roots of unity z k = e iθk corresponding to the angle θ for k = 1 and

k = N -1. Now, det(C -z i I) = (1 -z i ) det C -z i I 2 , i = 1, N -1,
where I is the N × N identity matrix; since C N -I = 0 we have

det C -λ I = (-1) N (λ N -1) (2.13) from which we deduce that det C -z i I = (-z i ) N + (-1) N +1 = 0, i = 1, N -1. This shows that A = C N -1 (C -z 1 I)(C -z N -1 I) is not invertible.
The system (2.12) can be therefore solved; the following theorem ensures its solvability.

Theorem 2.3.2. The system (2.12) admits a four-dimensional affine space of solutions.

Proof. Let b = (2a 0 )M 0,0 be the right-hand side of (2.12); we want to show that b ∈ Im(A) or, equivalently, W b = 0, where Im(W T ) = Leftker(A) = Ker(A) since A is symmetric.

Recalling the identity

0 = C N -I = N -1 j=0 (C -z j I),
with z k = e iθk the N -th roots of unity, it is easy to prove that we can take

W T = N -2 j=0 j̸ =1 (C -z j I) = diag(0, W T , W T ),
in which

W T = N -2 j=0 j̸ =1 ( C -z j I) (2.14)
and C as in (2.11): in fact

W A = (AW T ) T =     C N -1 (C -z 1 I)(C -z N -1 I) N -2 j=0 j̸ =1 (C -z j I)     T =   C N -1 N -1 j=0 (C -z j I)   T = 0
since the matrices (Cz j I) commute. Collecting the constants from the product we have

W b = b T W T T =   2 -a 0 N (N + 5) diag(N 2 , 4 I, I)(1 -z 0 , 1 -z 0 , 1 -z 0 ) N -2 j=2 (C -z j I)   T = 0
being z 0 = 1; due to the particular block structure of A, we deduce that corank(A) = 4. Thus the solutions of (2.12) form an affine space of dimension 4.

From the previous Theorem 2.3.2, in order to compute a solution of (2.12) we need to add 4 extra constraints to the coefficient matrix; then, since we arrive to an overdetermined system, the solution we get has to be meant in the least squares sense. To do that, let k 1 and k 2 be the two (column) vectors of one basis of W T : with this choice it results that

Im(W T ) = Span{k 1 , k 2 , k 3 , k 4 },
where

k 1 = (0, k T 1 , 0) T , k 2 = (0, k T 2 , 0) T , k 3 = (0, 0, k T 1 ) T , k 4 = (0, 0, k T 2 ) T .
Here we decided to impose the 4 extra constraints such that the new masks results to be as close as possible to M 1,0 ; this can be achieved projecting their difference onto Im(W T ); this projection, in principle, can be executed arbitrary, i.e. choosing any direction. Let M * be a particular solution (2.12); then, the general solution of (2.12) is given by

M 1,0 = M * + 4 i=1 µ i k i , µ i ∈ R. (2.15) Since the rows of A sum up to 2 -a 0 , i.e. A (1, 1, 1) T = (2 -a 0 ) (1, 1, 1) T , it is immediate
to notice that we can take M * = M 0,0 . Here we choose to perform an orthogonal projection

(M 1,0 -M 1,0 ) ⊥ Im(W T ), that is ⟨M 1,0 -M 1,0 | k i ⟩ = 0, i = 1, 2, 3, 4.
This leads to the solution obtained by solving

Kµ = b, (2.16)
where K = diag( K, K) and

K = ⟨ k i | k j ⟩ i,j=1,2
is the 2 × 2 orthogonal projection matrix and

b = k 1 k 2 k 3 k 4 T M 1,0 .
We also observe that if we solve (2.15) for N = 3, we obtain again the ACC 5 mask M 1,0 .

Example 2.3.3. In the case of an EV with valence N = 6, using the construction in Theorem 2.3.2 we obtain

W T = Span{k 1 , k 2 , k 3 , k 4 },
where

k T 1 = (1, 0, -1, -1, 0, 1), k T 2 = (-1, -1, 0, 1, 1, 0).
Hence, from (2.15) we compute the coefficients for the orthogonal projection which are

µ 1 = 3 110 , µ 2 = - 3 110 , µ 3 = 3 220 , µ 4 = 0,
and substituting them in (2.16) we obtain the weights for the mask M 1,0

M 1,0 = 1 660 [360,

Bézier masks of second order

Relation (2.5), i.e. the equation for the G 1 constraints among Bézier points concerning second order derivatives, can be solved both with respect to M 1,1 or M 2,0 . One approach, that is solving (2.5) w.r.t. the mask M 1,1 , will lead to a circulant linear system which presents different characteristics depending on the odd or even valence N of the EV to which it is attached. On the other hand, if we solve (2.5) w.r.t. M 2,0 , we end up with a direct formula for the solution.

Let us investigate the two cases.

2.3.2.1

Deriving M 1,1 by assigning M 2,0 : circulant system approach

We start by defining the mask M 2,0 , which is done differently for odd and even valence.

Similarly to (2.12), since M (0)

1,1 = C M 1,1 , (2.5) translates into the following linear system:

(I + C)M 1,1 = d , where d = 1 5 a 0 M 0,0 -5(a 0 -2)M 1,0 + 4a 0 M 2,0 .
(2.17)

The Bézier points involved in this relation are highlighted in Figure 2.3-(b). Let B = I + C be the matrix in (2.17). To solve the system we have to distinguish the cases N odd and N even.

Case N odd. In this case, from (2.10), replacing the masks we have we get to

M 2,0 = [δ • , δ, δ ′ ] T = 1 10 -M 0,0 + 5M 1,0 + 10M 3,0 -5M 4,0 + M 5,0 . (2.18)
We have the following:

Proposition 2.3.4. The matrix B = I + C is invertible with B -1 = diag 1 2 , B -1 , B -1 , B = I + C and B -1 = 1 2 N j=1 (-1) j-1 C j-1 . (2.19) Proof. It follows from   N j=1 (-1) j-1 C j-1   I + C = N j=1 (-1) j-1 C j-1 + N j=1 (-1) j-1 C j = N j=1 (-1) j-1 C j-1 + N +1 k=2 (-1) k-2 C k-1 = C 0 + C N = 2 I.
By using (2.19) are now able to compute explicitly M 1,1 .

Corollary 2.3.5. The weights for the mask

M 1,1 = [γ • , γ, γ ′ ] T in the odd case are given in terms of the masks M 0,0 = [α • , α, α ′ ] T , M 1,0 = [β • , β, β ′ ] T and M 2,0 = [δ • , δ, δ ′ ]
T by the following relations:

γ • = 1 5 (5 -2a 0 )α • + 2 5 a 0 δ • , γ i = 1 10 a 0 N j=1 (-1) j-1 α i+j-1 - 1 2 (a 0 -2) N j=1 (-1) j-1 β i+j-1 + 2 5 a 0 N j=1 (-1) j-1 δ i+j-1 , γ ′ i = 1 10 a 0 N j=1 (-1) j-1 α ′ i+j-1 - 1 2 (a 0 -2) N j=1 (-1) j-1 β ′ i+j-1 + 2 5 a 0 N j=1 (-1) j-1 δ ′ i+j-1 ,
with i = 1, . . . , N and where the indices i + j -1 are understood mod N when i + j -1 > N .

Proof. Since from (2.17) it results that M 1,1 = B -1 d, the proof is obtained by plugging in (2.19).

Case N even. If we are in presence of an EV with even valence, the matrix B is not invertible and results corank(B) = 2, due to corank( B) = 1. In fact, substituting λ = -1 in (2.13) it results that det( B) = 0 and B 1,1 = 1, where B 1,1 denotes the minor of B obtained removing the first row and column and since its cofactor is an upper triangular matrix with only ones on the diagonal. This implies corank(B) = 2. Hence, we must add two extra constraints in order to have a solution for the system; the obtained solution, also in this case, has to be meant in the least squares sense. First of all we must check if system (2.17) admits a solution.

Denoting Y = (w 1 w 2 ),
where

w 1 = (0, w T , 0) T , w 2 = (0, 0, w T ) T , w = (1, -1, 1, . . . , -1) T , it follows that Im(Y T ) = Span{w 1 , w 2 } = Leftker(B);
in fact, it easy to notice that Y T B = 0. Then, the system (2.17) admits a solution if and only if d ∈ Im(B) or, equivalently,

Y T d = 0. (2.20) Relation (2.20) is known as Vertex Enclosure Relation (V.E.R.). Let now denote M even 2,0 = [δ even • , δ even , δ ′ even ] T = 1 10 -M 0,0 + 5M 1,0 + 10M 3,0 -5M 4,0 + M 5,0 .
Checking property (2.20) for the masks M 0,0 , M 1,0 and M even 2,0 in (2.17), we notice that

Y T M 0,0 = 0, Y T M 1,0 = 0 but Y T M even 2,0 ̸ = 0.
Thus, in order to have a solution for system (2.17) we need to compute a new mask M 2,0 = [δ • , δ, δ ′ ] T ensuring (2.20). To do this we notice that

Im(B) = Leftker(Im(Y ));
in fact, it follows that BY = 0 due to the particular structure of B. Moreover, as for M 1,0 , we want that the mask M 2,0 results close to M even 2,0 . Therefore we compute the new mask M 2,0 by projecting its difference with M even 2,0 onto the left kernel B and forcing the V.E.R to be verified; since we decided to operate an orthogonal projection, the sought mask M 2,0 is obtained imposing ⟨M 2,0 -M even 2,0 | r i ⟩ = 0, i = 1, . . . , 2N + 1, where r i are the rows of B i.e. the vectors of one basis of Leftker(Im(Y )). Using the bilinearity of the scalar product we can rewrite the projection as

⟨M 2,0 | r i ⟩ = ⟨M even 2,0 | r i ⟩, which is in fact B M 2,0 = B M even 2,0 .
Finally, forcing the V.E.R. to be verified by the desired mask w T 1 M 2,0 = 0 and w T 2 M 2,0 = 0, we are able to condense the previous relations in the following linear system

B M 2,0 = f , where B =   B w T 1 w T 2   and f =   B M even 2,0 0 0   .
(2.21)

Reordering the rows of B we are able to achieve a diagonal block structure for the matrix in (2.21) with blocks (2, B, B) which brings to the following

Proposition 2.3.6. Let B = B w T . The pseudo-inverse matrix of B is a block matrix B+ = (B 1 |B 2 ) where B 1 = 1 2N N j=1 (-1) j-1 (N -2j + 1) C j-1 and B 2 = 1 N w, (2.22) with B 1 ∈ R N ×N and B 2 ∈ R N ×1 .
Proof. By definition of pseudo-invers matrix, B+ = ( BT B) -1 BT ; therefore it suffices to show that

( BT B) (B 1 |B 2 ) = BT .
Due to the block structure of B we have

BT B = B T B + w w T = 2 I + C + C -1 + w w T .
Noticing that

w w T = N j=1 (-1) j-1 j C j-1 , it results w w T B 1 = N j=1 N -2j + 1 = 0 since C w = C -1 w = -w. This implies that 2 I + C + C -1 w = 0. (2.23) Let H = N j=1 (-1) j-1 j C j-1 .
From (2.22) we deduce that

( BT B)B 1 = 2 I + C + C -1 + w w T B 1 = 2 I + C + C -1 1 2 w w T - 1 N H + 1 2N w w T = - 1 N 2H + CH + C -1 H ,
(2.24) and using the notation in (2.11), we have

H = Circ(1, -2, 3, . . . , -N ) , CH = Circ(-N, 1, -2, . . . , N -1) , C -1 H = Circ(-2, 3, -4, . . . , 1) , from which 2H + CH + C -1 H = -N Circ(1, 0, 0, . . . , 1) = -N B T .
(2.25) Substituting (2.25) in (2.24) we arrive at

( BT B)B 1 = B T .
In a similar way, using the identity (2.23),

( BT B)B 2 = 2 I + C + C -1 + w w T 1 N w = 1 N w w T w = w since w T w = N .
From Proposition 2.3.6, we are now able to compute explicitly the new mask M 2,0 ensuring the V.E.R as follow:

Corollary 2.3.7. The mask M 2,0 = [δ • , δ, δ ′ ] T ensuring (2.20) is obtained from M even 2,0 = [δ even • , δ even , δ
′ even ] T as follows:

δ • = δ even • , δ i = 1 2N N j=1 (-1) j-1 (N -2j + 1) δ even i+j-1 + δ even i+j , δ ′ i = 1 2N N j=1 (-1) j-1 (N -2j + 1) δ ′even i+j-1 + δ ′even i+j ,
with i = 1, . . . , N and where the indices i + j -1 and i + j have to be meant respectively mod N when i + j -1 > N and i + j > N .

Proof. Since M 2,0 = B + f the proof is obtained by the explicit formulas (2.21) and (2.22), using the block structure of B.

Hence, using mask M 2,0 we are sure to have a solution for (2.17). Now, since corank(B) = 2, to compute M 1,1 we need to add two extra constraints: once more we want to obtain a mask M 1,1 the closest possible to M 1,1 and it is ensure projecting their difference onto Im(Y T ) = Leftker(B). Like in the previous cases, we operate an orthogonal projection i.e.

⟨M 1,1 -M 1,1 | w i ⟩ = 0, i = 1, 2, which translates into ⟨M 1,1 | w i ⟩ = ⟨M 1,1 | w i ⟩.
Joining these two additional relations to system (2.17) and denoting

d =   d M 1,1 | w 1 M 1,1 | w 2   , (2.26)
the mask for M 1,1 satisfies the relation

B M 1,1 = d.
Similarly to Corollary 2.3.7, we have:

Corollary 2.3.8. The mask M 1,1 = [γ • , γ, γ ′ ] T in the even case is given in terms of M 0,0 = [α • , α, α ′ ] T , M 1,0 = [β • , β, β ′ ] T , M 1,1 = [γ • , γ, γ ′ ] T and M 2,0 = [δ • , δ, δ ′ ] T
by the following relations:

γ • = 1 5 (5 -2a 0 )α • + 2 5 a 0 δ • γ i = 1 10N a 0 N j=1 (-1) j-1 (N -2j + 1)α i+j-1 - 1 2N (a 0 -2) N j=1 (-1) j-1 (N -2j + 1)β i+j-1 + 2 5N a 0 N j=1 (-1) j-1 (N -2j + 1)δ i+j-1 + 1 N w i γ | w , γ ′ i = 1 10N a 0 N j=1 (-1) j-1 (N -2j + 1)α ′ i+j-1 - 1 2N (a 0 -2) N j=1 (-1) j-1 (N -2j + 1)β ′ i+j-1 + 2 5N a 0 N j=1 (-1) j-1 (N -2j + 1)δ ′ i+j-1 + 1 N w i γ ′ | w ,
with i = 1, . . . , N and where the indices i + j -1 have to be meant mod N when i + j -1 > N .

Proof. Because M 1,1 = B + d, the proof is obtained using (2.22), (2.26) and the block shape of B.

Example 2.3.9. For an EV with valence N = 6, and using the formulas in Corollary 2.3.7, the mask M 2,0 satisfying (2.20) is given by We notice that negative weights may occur.

M 2,0 = 1 1320 [592,

2.3.2.2

Deriving M 2,0 by assigning M 1,1 : direct approach based on ACC 5

If we assume the cross derivatives to be fixed, substituting the known masks in (2.5), we directly obtain

M 2,0 = 1 4a 0 -a 0 M 0,0 + 5(a 0 -2)M 1,0 + 5(M 1,1 + M (0) 1,1 ) , (2.28)
that is well defined since a 0 ̸ = 0 for N ̸ = 4, which is the case for the EVs we are analyzing. In Figure 2.3-(c) is figured the location of the treated control points.

Third and fourth order Bézier masks

First of all, depending on the previous solving strategies, from (2.10) we define the new mask

M 3,0 = 1 10 M 0,0 -5M 1,0 + 10M 2,0 + 5M 4,0 -M 5,0 .
When we modify the control point related to the mask M 3,0 the regularity relation for C 2 smoothness in Bézier surfaces across a vertex does not hold anymore. The resulting surface will be C 1 smooth around regular vertices linked to EVs instead of the initial C 2 regularity. From (2.6), it follows that

M 2,1 + M (0) 1,2 = r, where r = 1 10 -a 0 M 0,0 + 5a 0 M 1,0 -10(a 0 -2)M 2,0 + 6a 0 M 3,0 , (2.29)
where Figure 2.3-(d) presents the location of the involved Bézier points. By using the reparametrization M

(0) 1,2 -M (0) 1,2 = θ M 2,1 -M 2,
1 we obtain the one parameter family of solutions of (2.29)

     M 2,1 = 1 θ + 1 r + θM 2,1 -M (0) 1,2 , M (0) 1,2 = r θ θ + 1 - 1 θ + 1 θM 2,1 -M (0) 1,2 , (2.30)
where the parameter θ ∈ R \ {-1} has to be fixed. If we take θ = 1, (2.30) will return a symmetric solution; indeed, this choice reflects the fact that we are looking for a solution by projecting orthogonally the ACC 5 masks (M 2,1 , M

1,2 ) onto the G 1 solution space. For all the other values of θ \ {-1, 1} we will obtain a non symmetric solution for (2.30); in particular, for θ = 0, we obtain the asymmetric solution by projecting along the direction M (0) 1,2 . Moreover, with θ = ±∞ we mean the projection along the direction of M 2,1 ; as expected from (2.30), we don't have any solution for the value θ = -1 because it translates into a projection along a parallel direction to the G 1 solution space. For the fourth order Bézier masks, from (2.7) we have

M 3,1 +M (0) 1,3 = r, where r = 1 10 a 0 M 0,0 -5a 0 M 1,0 + 10a 0 M 2,0 -10(a 0 -2)M 3,0 + 4a 0 M 4,0 ,
(2.31) which can be solved as above, replacing r in (2.30) with r in (2.31); the positioning of the considered control points is shown in Figure 2.3-(e).

Treatment of boundaries

Meshes are not always closed and, therefore, they can possess boundaries in which EVs can eventually lay. Bézier control points returning G 1 smoothness around boundary regions can be obtained by applying (2.4) to (2.10) starting from the boundary masks presented in [START_REF] Loop | Approximating Catmull-Clark subdivision surfaces with bicubic patches[END_REF], after degree elevation. Since these masks ensure the interpolation of boundary edges, we are not going to modify them; we will recover the G 1 regularity by imposing G 1 relations only in the inner regions and using the ACC 5 boundary masks as initial and final value (see Figure 2.4). Let κ be the numbers of patches sharing the boundary vertex: for the Bézier masks of order one, imposing (2.4) (κ -1) times across the inner extraordinary edges we obtain, for κ ≥ 4, the following linear system:

A M 1,0 = h, where h =          (2 -a 0 )F (1) 0,0 -F (0) 1,0 (2 -a 0 )F (2) 0,0 . . . (2 -a 0 )F (κ-2) 0,0 (2 -a 0 )F (κ-1) 0,0 -F (κ-1) 0,1          , (2.32) M 1,0 = M (1) 1,0 , . . . , M (κ-1) 1,0 , F (1,...,κ-1) 0,0 , F (0)
1,0 and F , and A = Tridiag(-a 0 , 1, 1) ∈ R (κ-1)×(κ-1) , that is a tridiagonal matrix with entries -a 0 along the main diagonal and entries 1 on the first upper and lower diagonal. For κ = 3 it results that A = -a 0 I + C and the vector h has entries (2a 0 )F

(1)

0,0 -F (0) 1,0 and (2 -a 0 )F (κ-1) 0,0 -F (κ-1) 0,1
. System (2.32) can be solved with a similar strategy like in Section 2.3.1. Regarding the second order masks, even in this case, from (2.10) and as in (2.18), we have

M (i) 2,0 = 1 10 -F (i) 0,0 + 5M (i) 1,0 + 10M (i) 3,0 -5M (i) 4,0 + M (i) 5,0 , i = 1, . . . , κ -1.
Now, for the cross derivatives, imposing (2.5) (κ -1) times we get the system:

     1 1 0 . . . 0 1 1 . . . 0 . . . 0 1 1      M 1,1 = g, (2.33) where g = 1 5 a 0 F (i) 0,0 + 5(2 -a 0 )M (i) 1,0 + 4a 0 M (i) 2,0 , i = 1, . . . , κ -1, and M 1,1 = M (0) 1,1 , . . . , M (κ-1) 1,1
. Let R ∈ R (κ-1)×κ be the rectangular matrix in (2.33); since corank(R) = 1 we have a degree of freedom for the solution. If we want the masks to be close to ACC 5 , we can fix this free degree, for example, projecting orthogonally the solution of (2.33) onto the ACC 5 space; this translates into the solution

M 1,1 = M 1,1 + R T ρ,
where ρ is the vector containing the coefficients for the orthogonal projection given by

R R T ρ = -R M 1,1 + g.
The third and fourth order Bézier points are obtained in the same way as the case for an inner EV using (2.29) and (2.31) along the inner extraordinary edges.

Analysis of the solutions and numerical results

From the constructions presented in Section 2.3 it turns out that we can solve system (2.4)-(2.10) in multiple ways, obtaining a family of G 1 schemes. We are now going to analyze this family and identify the procedure returning the most regular surface. First, we will describe the results obtained using a set of test meshes [Pet] (cf. Figure 2.5) which contains shapes with EVs of different valence and particular geometric structure; from the results obtained by this procedure we will find out the element of the family returning the smoothest surface. Then, we will check the strength of the previous-found best strategy by applying it to the more complex meshes in Figure 2.12 and showing numerical results for the G 1 continuity.

Comparing the different schemes

Let us call CS (circulant system) the solving strategy in Section 2.3.2.1, NCS (no circulant system) the strategy in 2.3.2.2 and S (symmetric) and AS (asymmetric) the solutions for the higher derivatives obtained in 2.3.3 respectively with θ = 1 and θ = 0; then, we can define four global solving strategies as: NCS-S, NCS-AS, CS-AS and CS-S. To identify the smoothest surface we carried out an isophote and curvature analysis together with numerical evaluation of the G 2 errors on the surfaces obtained by each strategy. In Figures 2.6 and 2.7 the results for two meshes with EV of valence N = 5 and N = 8 are shown. In strategy NCS-S we notice that for both valences, isophotes are deviated away from the EV creating an unwanted white spot around it, as also underlined by the curvatures. For the asymmetric strategies NCS-AS and CS-AS we easily realize that the asymmetry of the Bézier points translates into a rotation of the patches around the EVs, which implies irregularities in the isophotes and unpleasant curvatures. Instead, in strategy CS-S, we see that both isophotes and curvatures related to the two valences are smooth and regular in a neighborhood of the EVs. To evaluate numerically the results we use the quantities

n := E ∥ n L -n R ∥ 2 2 , c G := max e∈E max u∈e c L G -c R G , c M := max e∈E max u∈e c L M -c R M
, where E := e∈E e with E the set of all edges e of the mesh M and n

L , c L G , c L M , n R , c R G , c R M
are respectively the unit normal vector, the Gauss curvature and the mean curvature of, respectively, the left and right patch sharing the common edge e. The G 1 quality is evaluated by n , which is a measure of the distance of our surface from the set of purely C 1 parametric surfaces. Indeed, this measure is zero on any C 1 interface, while it expresses the discrepancy or "jump" of the (non-unit) normal vector when there is a discontinuity across an interface. Moreover, the quantities c G and c M return, respectively, the maximum difference between the Gauss curvatures and mean curvatures on the patch interfaces: Table 3.25 shows the values for these quantities of each solving strategy carried out on the test meshes [Pet] in Figure 2.5. At this point, using both the numerical and visual approach, we can state that in any case the asymmetric approach AS is never returning very regular surfaces; on the other hand we see that in most cases the best approach is given by the strategy CS-S. From the numerical results presented in Table 2.1, we notice that for meshes with very quick and sudden change of slope of patches around EVs, strategy NCS-S returns slightly better surfaces in terms of curvature compared to CS-S. 1.2e-00 38.8e-00 0.3e-00 15.8e-00 6.0e-00 1.0e-00 47.5e-00 84.4e-00 0.6e-00

(a) el3 (b) hyp3 (c) linel3 (d) stairs3 (e) el5 (f ) house5 (g) beam6 (h) linel6 (i) saddle6 (j) smbeam6 (k) monk7 (l) el8 (m) bas8 (n) kpa8 (o) mitsu9 (p) bas10 (q) updwn10 (r) cube3456
Table 2.1. G 1 smoothness n , Gauss curvature difference cG and mean curvature difference cM obtained from strategies NCS-S, NCS-AS, CS-AS, CS-S for every test mesh in Figure 2.5.

Complex meshes

Using the best strategy CS-S obtained in Section 2.4.1, we developed a G 1 continuity analysis on the complex meshes in Figures 2.8 to 2.12; these meshes present a very complicate geometry together with a very high number of faces, vertices and, obviously, extraordinary vertices of several valence. In Table 2.2 we present a detailed description of each mesh: we show the number of EVs of any valence, the total number of EVs, the number of vertices and faces of the mesh and the numerical evaluation of the G 1 regularity. From the results in Table 3.25 and Table 2.2, we can affirm that this construction returns really smooth surfaces, with a G 1 error in the range of 10 -11 -10 -15 , i.e. numerically identical to zero. .12. Other surfaces obtained by complex meshes using the solving strategy CS-S.

Comparison with ACC 3 surface and Catmull-Clark limit surface

In order to investigate the accuracy of our G 1 approximation of the Catmull-Clark limit surface we analyse on the distance between the two surfaces; the two quantities we take into account are the geometry error and the normal error. For each point laying on the Catmull-Clark limit surface, to compute the geometry error we look for their closest points on the G 1 surface: formally, if P CC is a point on the Catmull-Clark surface, the corresponding geometry error geom err is obtaining by minimizing the distance dist(P CC , G 1 ) = min

P ∈G 1 ∥P CC -P ∥ 2 ,
where we define with G 1 to be our CS-S surface. To evaluate the normal error, we first compute the orthogonal projection P ⊥ of each point P CC onto G 1 . Then the normal error norm err is obtained evaluating the 2-norm of the corresponding unit normals attached to them

norm err = ∥n(P CC ) -n(P ⊥ )∥ 2 .
The same study is conducted for the ACC 3 surface in order to provide a comparison with our improved construction and the obtained result are shown in Figure 2.13. We also estimate the convergence rates of geom err and norm err for both the surfaces ACC 3 and G 1 . Starting from different Catmull-Clark subdivision levels of the same input mesh, we construct on each of them the two surfaces and we evaluate the maximum of the two errors over a dense sampling of the surface.

In accordance with the quadratic convergence of the Catmull-Clark subdivision scheme, the geometric error geom err of our G 1 construction converges quadratically, while norm err converges with a rate of 0.7.

Refinement step

Error We also observed a convergence rate of 1.5 for the normal vector obtained as cross product between the coordinate derivatives, which is in general a non-unit vector. On the other side, in the ACC 3 construction we observe that geom err converges with a rate of 1.5 while norm err presents an oscillatory non-convergent behavior. Figure 2.14 summarizes these results. Moreover, from the numerical results we also noticed that this construction returns unit normals at the vertices which are close to the real normal of the Catmull-Clark surface. Let T u , T v = C T u be the masks returning the tangents of the limit Catmull-Clark surface with weights described in [START_REF] Loop | Approximating Catmull-Clark subdivision surfaces with bicubic patches[END_REF]pag. 3,6] and Chapter 3

ACC 3 geom G 1 ACC geom ACC 3 norm G 1 ACC norm Slope 1.5 Slope 2 Slope 0.72
∂ u = M 1,0 -M 0,0 , ∂ v = M 0,1 -M 0,0 = C ∂

Geometrically smooth functions for point cloud fitting

This chapter is devoted to the construction of a G 1 set of basis function generating the space of geometrically smooth splines defined over a quad mesh with the aim of applying them in point cloud fitting problem. Point clouds arise in many different fields, from earth surveying to medical imaging, and being able to have a smooth representation of it is a crucial task. Here, we investigate this topic by using G 1 Bézier functions defined over a quadrilateral mesh.

After having recalled the equations generating a G 1 smooth spline space we will go further in the construction of a set of basis function spanning the targeted space; moreover, we also provide a combinatorial formula for the space dimension depending only on the main features of the given mesh. The numerical experiments we provide to show the quality of the resulting fitting are mainly obtained by using parametrized point cloud, i.e. points already equipped with a parametrization. Moreover, we present a full pipeline to deal with general point clouds in order to compute a quad mesh and a parametrization. This chapter is based on [START_REF] Marsala | G 1 spline functions for point cloud fitting[END_REF].

3.1 G 1 spline space on a mesh M Here we will recall the equation defining a G 1 joint between two adjacent Bézier patches nearby EVs on a quad mesh M: similarly to relations (2.4) to (2.10), if b (k) i,j are the control points laying on the k-th patch (see Figure 2.2 for their labelling), the geometric continuity constraints across an edge linked to an EV of valence N are: b

(1) 0,1 + b (0) 1,0 = ā0 b (1) 0,0 + a 0 b (1) 1,0 , (3.1) 5(b (1) 1,1 + b (0) 1,1 ) = a 0 b (1) 0,0 + 5ā 0 b (1) 1,0 + 4a 0 b (1) 2,0 , (3.2) 10(b (1) 2,1 + b (0) 1,2 ) = -a 0 b (1) 0,0 + 5a 0 b (1) 1,0 + 10ā 0 b (1) 2,0 + 6a 0 b (1) 3,0 , (3.3) 10(b (1) 3,1 + b (0) 1,3 ) = a 0 b (1) 0,0 -5a 0 b (1) 1,0 + 10a 0 b (1) 2,0 + 10ā 0 b (1) 3,0 + 4a 0 b (1) 4,0 , (3.4) b (1) 4,1 + b (0) 1,4 = 2b (1) 4,0 , (3.5) b (1) 5,1 + b (0) 1,5 = 2b (1) 5,0 , (3.6) 10(b (1) 3,0 -b (1) 2,0 ) = b (1) 0,0 -5b (1) 1,0 + 5b (1) 4,0 -b (1) 5,0 , (3.7)
with ā0 = 2a 0 . The gluing data function we use is the same as in (2.2). Starting from the above equations (3.1) to (3.7), we shall present an explicit construction of a set of basis functions generating the G 1 spline space that will follow the structure of the input mesh M. For this purpose, it is interesting to count its main features: we will refer to n V as the number of vertices of the mesh, which can be subdivided in n IEV inner EVs, n BEV boundary EVs and n RV regular vertices. The total number of edges is denoted by n E and is composed by the number of extraordinary and inner regular edges, n EE and n IRE respectively, that is, edges attached to extraordinary and regular vertices, as well as the number of boundary edges n BE . With n F we refer to the number of faces and lastly n C is the number of corners of the mesh.

Basis extraction

In this section we present an explicit construction of a set of basis functions generating the G 1 space over a quad mesh M, which we refer to as B. Following the topology of the input mesh M, we can distinguish the set of functions attached to the vertices B V (i.e. spanned by functions whose support lies on all the patch(es) sharing the vertex), set of basis functions attached to the edges B E (i.e. spanned by functions whose support lies on the patch(es) sharing the edge) and set of basis functions attached to the faces B F (i.e. spanned by functions whose support lies only on the interior of a single patch). This decomposition will be exploited in the proofs of Section 3.3. Thus, we arrive at a set of basis functions which can be decomposed as

B = n V i=1 B V i ∪ n E i=1 B E i ∪ n F i=1 B F i . (3.8)
Also in this construction we assume that all the EVs in M are isolated, that is, their one-ring neighborhood contains only RVs. Having EVs attached to regular vertices only, which translates into (3.5) and (3.6), implies that vertex, edge and face functions have value and derivative equal to zero on the boundary of their support. This is because the G 1 constraints, defined making use of the gluing functions (2.2), do not affect the final six points across the edge (i.e. the points responsible of the value and derivative of the basis functions at the endpoint of the edge) and they can therefore be set to zero. Thus, basis functions with local support on a collection of faces of M do not influence functions in the region surrounding its support. By combining (3.1) to (3.7) circularly around all the edges attached to an EV we can reformulate the G 1 constraints using a block staircase matrix, which is useful to better understand the basis extraction and their analysis which follows. Let b i,j = b (k) i,j , k = 1, . . . , N , be the vector containing all the points b (k) i,j attached to the neighborhood of the EV we are considering, u = (1, . . . , 1

N ) T , C ∈ R N ×N and the circulant matrix C introduced in (2.11) C =        0 1 0 • • • 0 1 . . . 1 1 0 • • • 0        = Circ(0, 1, 0, . . . , 0).
Note that b 4,0 = b 5,0 = 0, as a consequence of the locality of the G 1 constraints around the EV. The full system of G 1 relations around an EV can be written as

      -ā 0 u C 1 -a 0 u -5ā 0 I -4a 0 I 5C 2 -u 5I -10I 10I a 0 u -5a 0 I -10ā 0 I -6a 0 I 10I 10I -a 0 u 5a 0 I -10a 0 I 10ā 0 I 10I 10I                      b 0,0 b 1,0 b 2,0 b 1,1 b 3,0 b 2,1 b 1,2 b 3,1 b 1,3                = 0, (3.9)
where a 0 = 2 cos (2π/N ) is the first Bernstein coefficient of the gluing data function (2.2), ā0 = 2a 0 , I ∈ R N ×N is the identity matrix,

C 1 = -a 0 I + C + C N -1 , C 2 = I + C, (3.10) and b 1,2 = C b 1,2 , b 1,3 = C b 1,3 , with C in (2.11
). The strategy we exploit to obtain the Bézier coefficients defining the basis functions is the following: starting from (3.1) to (3.7) or eq. (3.9) we impose, one by one, the value one to each "free" coefficient involved in the G 1 constraints. Then, with this initial value we solve the G 1 relations (3.1) to (3.7) (or (3.9)), while we gradually set the value of any unconstrained coefficients that we encounter (i.e. coefficients which lead to an overdetermined equation) to zero.

The set B V of vertex basis functions

Here belong basis functions connected to inner and boundary vertices (extraordinary or regular) and corner vertices. We only provide the explicit construction for functions attached to inner extraordinary and inner regular vertices, since the construction for the remaining cases is analogous.

Construction of basis functions associated to an inner EV

Given an extraordinary vertex of valence N , we associate to it a basis function responsible of the value (of elements in the G 1 spline space) at the vertex and basis functions related to first and cross derivatives. We deduce during the construction that these basis functions are N + 3 in total.

Basis function related to the vertex value. To extract the basis function associated to the value at the vertex, we start solving the system (3.1) to (3.7) with initial value b 0,0 = 1 and continuing the construction fixing zero values for all the control points which are not constrained by any relation we will encounter during the construction. With this assumption, we can rewrite (3.1) in the form

C 1 b 1,0 = (2 -a 0 )u, (3.11)
with C 1 ∈ R N ×N defined in (3.10). The solution of (3.11) returns the control points involved in the definition of the first derivative of the basis with unit value at the vertex. From Section 2.3.1, in particular from Proposition 2.3.1 and Theorem 2.3.2, we can deduce that the matrix

C 1 = -a 0 I + C + C N -1
is singular and corank(C 1 ) = 2; for this reason, in order to obtain a unique solution we need to insert two extra constraints to the system (3.11). Let

Ker(C 1 ) = K = Span{k 1 , k 2 }
be the kernel of the matrix C 1 generated by the two vectors k 1 and k 2 ; we can choose our solution to be orthogonal to the space K, i.e.

⟨b 1,0 |k 1 ⟩ = ⟨b 1,0 |k 2 ⟩ = 0.
To compute the kernel K we can use formula (2.14). This procedure allows us to achieve a unique solution for this set of control points. Going further in the resolution of the system, using the solution of (3.1) and the circulant matrix C, and taking into account the constraint along the edge in (3.7) which becomes

b 2,0 = 1 2 b 1,0 - 1 10 u, (3.12)
we can reorder (3.2) as

C 2 b 1,1 = - 1 5 a 0 u + 5(2 + a 0 )b 1,0 , (3.13)
with C 2 as in (3.10). Following the analysis in previous Section 2.3.2.1, for odd values of the valence N the matrix C 2 is invertible, while for even occurrences we have corank(C 2 ) = 1.

To obtain a unique solution in the singular case we need to fix an extra constraint which we choose to be the orthogonality of the expected solution b 1,1 to Ker(C 2 ). Regarding the control points for the higher derivatives b 2,1 and b 3,1 , from (3.3) and (3.4) and using again (3.12) we come up with the relations

b 2,0 = 1 2 b 1,0 ,
are given by the equations

C 2 b 1,1 = 1 5 (2 + a 0 )b 1,0 , b 2,1 + b 1,2 = b 1,0 , b 3,1 + b 1,3 = 0,
which can be solved in the same way as (3.13) and (3.14). being well defined since a 0 ̸ = 0 ⇐⇒ N ̸ = 4, that is the case of an EV we are in fact investigating.

Regarding higher order derivatives, we have the points b

(k) 3,0 , b (k+1) 3,0 , b (k-1) 1,2 , b (k) 2,1 , b (k) 1,2 , b (k+1) 2,1 and b (k-1) 1,3 , b (k) 3,1 , b (k) 1,3 , b (k+1) 3,1
defined by the relations b

(k) 3,0 = b (k+1) 3,0 = 5 4a 0 , b (k-1) 1,2 + b (k) 2,1 = b (k) 2,1 + b (k+1) 1,2 = 1 2 5 a 0 -1 , b (k-1) 1,3 + b (k) 3,1 = b (k) 3,1 + b (k+1) 1,3 = 5 2a

Basis functions at an inner regular vertex

In presence of a regular vertex (RV), i.e. a vertex with valence N = 4, we expect to have by construction C 1 -smooth basis functions; this is in fact what (3.5) and (3.6) state. First we need to investigate how many RV basis function we have. To do that we first need to expand cyclically (3.5) and (3.6) to all the control points around the vertex; this procedure, using the notation in Figure 3.3-(a) leads to the following system:

            1 0 0 -2 0 0 1 0 0 0 1 0 0 -2 0 0 1 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 1 0 0 -2 0 0 1 0 0 0 1 0 0 -2 0 0 1 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1                           P Q R S T U V W X               = 0. (3.16)
Let D be the coefficient matrix in (3.16). The number of basis functions attached to the RV is given by corank(D) = 9rank(D) = 9 -5 = 4. Hence, a choice for the Bézier points returning linearly independent functions (by construction) verifying (3.16) is given by the coefficients in Figure 3.3-(b). In presence of RV connected to an EV it is also necessary to modify more control points for each extraordinary edge, according with (3.3) and (3.4).

Basis functions linked to extraordinary and regular boundary vertices and corners

The extraction of this type of functions is analogous to the constructions developed in the previous sections for the inner cases. Naming with κ the valence of an extraordinary boundary vertex, i.e. the number of patches attached to it, imitating the process in Section 3.2.1.1 we come up with κ + 3 basis functions, which is equivalent to N + 2 since κ = N -1. On the other hand, to obtain the functions connected to a regular boundary vertex and corner we need to copy the procedure shown in Section 3.2.1.2: in both cases, following the same strategy leading to (3.16), we come up again with four basis functions.

The set B E of edge basis functions

In this second set of functions, we have basis functions attached to inner and boundary edges, either extraordinary and regular ones. In the same way as the vertex functions we will present the explicit construction in the case of inner extraordinary and regular edge functions, whereas the construction for the remaining case is analogous.

Construction of basis functions connected to an extraordinary edge

These functions are obtained starting from (3.3),(3.4) and (3.7). Similarly to the construction in Section 3.2.1.1, to extract the Bézier coefficients for the functions, we need to set zero values at the free points appearing in the equations. In the construction of basis functions connected to extraordinary edges, the control points we need to nullify are all the points laying on the edge, i.

Basis functions belonging to an inner regular edge

These basis functions are obtained with a similar approach as in Section 3.2.1.2 for an inner regular vertex; are involved in this construction the six control points in the inner part of the edge, i.e. away from the influence of the vertices' equations. In order to determine the number of these functions, applying the C 1 constraints in (3.5) and (3.6) to the two layers of control points implicated in this analysis with the notation given by Figure 3.5-(a) we have:

1 -2 1 0 0 0 0 0 0 1 -2 1         P Q R S T Ū        = 0. (3.18) (a)
1 0 If D is the matrix in (3.18), the number of basis functions attached to an inner regular edge is given by corank( D) = 4 and a set of possible solutions verifying these constraints returning linearly independent functions is given by the configurations in Figure 3.5-(b).

1 2 0 0 0 0 0 1 2 0 1 0 0 1 0 1 2 0 0 0 0 0 1 2 0 1 (b)

Boundary edge basis functions

In presence of a boundary edge we have no smoothness constraints to impose. Hence, the basis functions in this case are the classical bivariate Bézier polynomials obtained assigning the unit value to the four control points involved in this setting, one at the time, to obtain the four basis functions we were looking for.

The set B F of face basis functions

To conclude the construction of our basis we miss the definition of functions belonging uniquely to a single patch. As for the case of boundary edge functions, here we have no regularity conditions to impose; thus the construction is the same as in Section. 3.2.2.3 returning four Bézier polynomial on each face.

Analysis of the basis and space dimension

The functions we built in Section 3.2 form a basis of the G 1 spline space over the mesh M.

Denoting by G 1 (M) such space, it results that G 1 (M) = ⟨B⟩. In this section we will provide a proof for this statement as well as a dimension formula for the underlined space. The first part of the proof, obtained by contradiction, focuses on the linear independence of the basis functions and it is carried out exploiting the property (3.8) and analysing their supports; the second half, instead, uses the equations defining the smoothness constraints of the basis functions in order to prove that they generate the desired spline space.

Theorem 3.3.1. The functions introduced in Section 3.2 form a basis B for the space G 1 (M) over a quad mesh M.

Proof. In order to prove that the set B is a basis of our space we need to prove that it is linear independent and the property of being a generating set. For the linear independence, we proceed by contradiction. Since we assume to have isolated EVs only, we have by construction that the supports of vertex basis functions attached to different vertices are disjoint; the same property holds for edge basis functions defined over different edges as well as face basis functions attached to different faces. Hence, we can restrict our proof on basis functions attached on elements sharing a common vertex. Therefore we consider an EV of valence N with its N extraordinary edges and N faces (the case with regular vertices and regular edges is analogous) and we assume that there exists a linear dependence between the basis functions supported on these elements, i.e. there exist non zero coefficients α, β, γ such that

N +3 i=1 α i B EV i + N j=1 2 i=1 β EE j i B EE j i + N j=1 4 i=1 γ F j i B F j i = 0. (3.19) From (3.19) it holds that supp   N +3 i=1 α i B EV i + N j=1 2 i=1 β EE j i B EE j i   = supp   - N j=1 4 i=1 γ F j i B i F j   ,
and in particular supp

  N +3 i=1 α i B EV i + N j=1 2 i=1 β EE j i B EE j i   ⊆ N +3 i=1 supp B EV i ∪   N i=1   2 j=1 supp B EE j i     supp   N +3 i=1 α i B EV i + N j=1 2 i=1 β EE j i B EE j i   ⊆   N i=1   4 j=1 supp B F j i     . It follows that supp   N +3 i=1 α i B EV i + N j=1 2 i=1 β EE j i B EE j i   ⊆   N +3 i=1 supp B EV i ∪   N i=1   2 j=1 supp B EE j i       ∩   N i=1   4 j=1 supp B F j i     = ∅.
(3.20)

Hence, (3.20) implies

N +3 i=1 α i B EV i + N j=1 2 i=1 β EE j i B EE j i = 0 = N j=1 4 i=1 γ F j i B F j i ,
that is γ

F j i = 0 ∀ i, j since B F j
i are linearly independent by construction. Now we have

N +3 i=1 α i B EV i = - N j=1 2 i=1 β EE j i B EE j i .
Let us fix two adjacent patches k and k -1 sharing the EV we are considering and let us focus on the control points b

(k) 0,0 , b (k) 1,0 , b (k) 0,1 , b (k) 1,1 , b (k-1) 1,0 , b (k-1)
1,1 . By definition, these control points are always equal to zero for all the basis functions B EE j i , while this is not the case for the B EV i ; hence, if we now look at the submatrices corresponding to the points we are considering for both B EV i and B EE j i , we deduce that in the first case the selected submatrix is invertible, while the second submatrix is the null matrix. Therefore it follows that ∀i, α i = 0 since the B EV i are linearly independent by construction and consequently, using the same argument, β

EE j i = 0 ∀i, j.
Regarding the generating property, let f ∈ G 1 (M). Since the basis functions B F j i attached to a face F j are C 1 smooth we can define the function

f ′ := f - n F j=1 4 i=1 c F j i B F j i , (3.21)
which is still G 1 and is such that the inner face coefficients vanish. The same procedure can be applied to the C 1 basis function attached to corners B

C j i , boundary edges B BE j i
and inner regular edges B IRE j i which can be used to define a new function starting from (3.21) as

f ′′ := f ′ - n C j=1 4 i=1 c C j i B C j i - n BE j=1 4 i=1 c BE j i B BE j i - n IRE j=1 4 i=1 c IRE j i B IRE j i ∈ G 1 ,
so that it has null coefficients in the above mentioned sectors.

Let s be the indices of the 9 control points in Figure 3.3 related to the basis functions B Hence

f ′′ [s] = 4 i=1 c RV j i B RV j i [s];
we proceed similarly with boundary EVs B

BEV j i
and boundary regular vertices B

BRV j i
to define respectively coefficients c

BEV j i , c
BRV j i so that the function

f ′′′ := f ′′ - n RV j=1 4 i=1 c RV j i B RV j i - n BRV j=1 4 i=1 c BRV j i B BRV j i - n BEV j=1 N BEV j +2 i=1 c BEV j i B BEV j i
is still a G 1 smooth function and it has zero coefficients in the previously treated regions. Now, the function f ′′′ has only non-zero coefficients around extraordinary vertices and their incident edges. We define t as the indices of the non-zero coefficients for an extraordinary vertex EV j (black points in Figure 3.1). Similarly to the case of regular vertices, from (3.9),

Āf ′′′ [t] = 0,
where Ā is the matrix in (3.9) containing the G 1 constraints. Then,

Ker( Ā) = Span B EV j i [t]
which leads to

f ′′′ [t] = 3N EV j +3 i=1 c EV j i B EV j i [t].
Finally, we have that

f ′′′′ := f ′′′ - n EV j=1 3N EV j +3 i=1 c EV j i B EV j i = 0,
which concludes the proof.

As consequence of the particular structure of the set B we have the following formula for the dimension of our spline space.

Corollary 3.3.2. The space G 1 (M) has dimension given by:

dim G 1 (M) = n V i=1 |B V i | + n E i=1 |B E i | + n F i=1 |B F i | = n EV i=1 N EV i + 3n IEV + 2(n BEV + n EE ) + 4(n IRE + n BE + n RV + n C + n F ).
Proof. The proof is obtained using the decomposition in (3.8) and summing up all the members of the basis functions for each feature of the mesh shown in Sections 3.2.1, 3.2.2 and 3.2.3.

Numerical experiments

We conclude this chapter presenting numerical experiments in which we assess the basis functions constructed in Section 3.2 in point cloud fitting problems. The setup of our investigation is the classical least squares fitting problem: given a point cloud P, i.e. a set of points P i ∈ R 3 , i = 1, . . . , n P , with associated parameters (u i , v i ) ∈ [0, 1] 2 on the patch ℓ i , we want to find the coefficients

c i ∈ R 3 of a G 1 surface G = k c k B k such that the quantity F = n P i=1 G (ℓ i ) (u i , v i ) -P i 2 2 + λE thin , λ ≥ 0, (3.22)
is minimal, where G (ℓ) , ℓ = 1, . . . , n F , is the geometry map

G (ℓ) : [0, 1] 2 -→ R 3 , (u (ℓ) i , v (ℓ) i ) -→ G (ℓ) (u (ℓ) i , v (ℓ) i ) = 36 i,j=1 b (ℓ) i,j B (ℓ) i,j (u (ℓ) i , v (ℓ)
(ℓ) = G (ℓ) ([0, 1] 2 ).
We also take into account in our minimization problem (3.22) an energy term given by the standard thin-plate energy

E thin = n F ℓ=1 [0,1] 2 G (ℓ) uu 2 2 + 2 G (ℓ) uv 2 2 + G (ℓ) vv 2 2 du dv, (3.24)
which is controlled by the parameter λ. The minimization of the functional in (3.24) is responsible for a regularisation effect on the final surface G (i.e. oscillations). Moreover, following the construction presented in [START_REF] Hoschek | Intrinsic parametrization for approximation[END_REF], we perform few iterations of parameter correction to further reduce the approximation error, if needed. Since our focus is on the quality of the G 1 basis for fitting, the parametrized data that we use in the following numerical experiments are obtained by evaluating certain input functions or surfaces to obtain parameters. In particular, starting from a planar mesh, in order to build point cloud data from a given analytic function f , we compute parameters (u i , v i ) by randomly sampling a certain number of points in the unit square [0, 1] 2 , for each patch. Therefore, the point cloud is obtained by the triple

(x i , y i , f (x i , y i )), with (x i , y i ) = G (ℓ i ) (u i , v i ).
A similar procedure is performed to obtain a point cloud from a quad mesh yielding an AC 3 surface [LS08]: considering the geometry map G ACC 3 , we obtain sample points

(x i , y i , z i ) = G (ℓ i ) ACC 3 (u i , v i )
randomly on the surface, with (u i , v i ) ∈ [0, 1] 2 . Furthermore, we present two examples of data fitting starting from point clouds without any given structure where both the quad mesh and a parametrization are computed just starting from the points data.

After having computed a least squares surface G from a point cloud P, let us define the array of errors err = {err i } nP i=1 whose entries are the quantities

err i = G (ℓ i ) (u i , v i ) -P i 2 ,
i.e. the Euclidean distance (ℓ 2 norm) from each point of the cloud and the corresponding value on the surface evaluated at its parameter. From this we define two quantities which will be used to asses the accuracy of the fitting, namely the maximum error and the root mean squared error (RMSE in short):

L ∞ = max i=1,...,n P err i , RMSE = 1 n P n P i=1 err 2 i .
(3.25)

Point cloud by analytic function evaluation

The experiments we propose here use point cloud data obtained by sampling an input function f (x, y) over a domain identified by a quad mesh M. In Example 3.4.1 we focus on the quality of the fitting showing how it improves when the number of basis functions increases, without the need to use extra smoothness constraints. On the other hand, Example 3.4.2 exhibit the power of our construction when thin plate energy is used in order to obtain an optimal result, but without increasing the dimension of the spline space. Finally, Example 3.4.3 provide an example of fitting from a point cloud derived by a trivariate function defined over the unit sphere.

Example 3.4.1. For this test, the point cloud is obtained evaluating the function

f T (x, y) = y 2 (cos(4(x 2 + y -1))) , (x, y) ∈ M T (3.26)
over the mesh M T , defining a triangle, formed by 3 patches which present an EV of valence N = 3; our sampling produced a point cloud formed by 1536000 points. This example uses no smoothing, which means we fix λ = 0. From the result in Figure 3.6-(a) can be noticed that constructing our spline space over a coarse mesh as M T leads to a fitting surface which is not approximating in a proper way our point set, due to its very oscillatory behavior. This issue can be solved by increasing the number of functions generating the spline space, i.e. increasing the number of patches defining the polygonal domain; in this example the mesh has been refined via Catmull-Clark subdivision. In Table 3.1 and Figure 3.7 are respectively listed and plotted the errors computed from (3.25) for 5 refinement levels. of the G 1 spline space. Moreover, we investigate the convergence of the fitting error and, as shown in Figure 3.7, we recovered the optimal convergence rate 6 using biquintic splines.

Example 3.4.2. Here, the cloud data is derived sampling the function

f E (x, y) = ν∈{-3,0,3} 2 3e √ (10x+ν) 2 +(10y+ν) 2 , (x, y) ∈ M E , (3.27)
where M E is a quad mesh tracing out an hexagon composed of 96 patches identifying an EV of valence N = 6 in its middle. The point cloud we obtain is formed of 153600 points. Here we show the power of the smoothing property of this construction: differently from Example 3.4.1, fixing the number of patches i.e. the number of basis function, which is equal to 1725 in our which presents much less irregularities, as can be noticed in Figure 3.8-(d). As expected, this procedure will produce at every iteration a smoother function than the previous, but on the other side this forced regularity constraint is reflected in an increase of the three errors in (3.25); this phenomenon is shown in Table 3.3 and graphically in Figure 3.9. λ 0 10 -3 10 -2 10 -1 L ∞ error 0.164e-00 0.165e-00 0.170e-00 0.196e-00 RMSE 0.308e-02 0.310e-02 0.312e-02 0.375e-02 Example 3.4.3. In this example, the point cloud is obtained by evaluating a trivariate function defined over the unit sphere S 2 . More precisely, the points are sampled from the function f S 2 (x, y, z) = max{0, sin (2πx) sin (2πy) sin (2πz)} + 1, (x, y, z) ∈ S 2 .

(3.28)

The data we get is composed of 540000 points, while the set of basis functions has been built over the quad mesh M S 2 approximating the unit sphere with 96 faces. We notice from the color plot in Figure 3.10-(d) that, understandably, the error is concentrated in the regions where the point cloud has peaks. This is due to the noticeable slopes created by the presence of the max function in (3.28). Table 3.4 shows the numerical results for the errors.

dim G 1 (M S 2 ) L ∞ error RMSE f S 2 1512 0.539e-01 0.107e-01

Point cloud from ACC 3 surfaces

We now present fitting examples obtained from big data sets. The point clouds utilized in this section are provided by randomly sampling the Approximate Catmull-Clark surfaces (ACC 3 ) obtained from the construction in [START_REF] Loop | Approximating Catmull-Clark subdivision surfaces with bicubic patches[END_REF]. Figures 3.11 to 3.16 show the data we use for our investigation obtained from the ACC 3 surfaces, starting from parameters sampling: the dimension of the clouds goes from a minimum of 549180 to a maximum of 968704 points. In all the experiments presented here we do not consider any smoothing parameter, i.e. λ = 0, but few iterations of parameter correction [START_REF] Hoschek | Intrinsic parametrization for approximation[END_REF] are performed to optimize the fitting error; the errors in Table 3.5 are obtained using the formulas in (3.25). Since ACC 3 surfaces are not G 1 at the EVs, it is not surprising to see that the maximal fitting error is located around the extraordinary vertices.

Quadrilateral mesh generation, parametrization and fitting

In the previous experiments, all the point clouds were already equipped with a parametrization. This is because the goal of this chapter is to show the quality of the fitting with the basis functions we propose, rather than improving upon parameter computation. Nevertheless, we now present few examples in which both the parametrization of an unorganised point cloud and the construction of a coarse quadrilateral mesh are derived. Starting from an input point cloud, we first reconstruct a triangular mesh by using the algorithm in [START_REF] Bernardini | The ballpivoting algorithm for surface reconstruction[END_REF]; then we apply the so-called 4-8 subdivision scheme [START_REF] Velho | 4-8 Subdivision[END_REF] to obtain a quadrangular representation of the previous triangular mesh. This procedure leads to a mesh, which is much more refined than we would need for our computation. Hence, we coarsen it using Rhino3D command QuadRemesh (see [START_REF] Mcneel | Rhinoceros 3D[END_REF]), which reduces the number of faces of the quad mesh without modifying its topology. This quad mesh M is the mesh on which we define the G 1 space and basis functions.

In order to associate parameters to a point P i = (x i , y i , z i ) of the point cloud, we compute its orthogonal projection

P ⊥ i = (x ⊥ i , y ⊥ i , z ⊥ i )
onto the G 1 spline surface G obtained from the construction in Chapter 2. Assume that the point P ⊥ i is laying on the ℓ-th patch of G. Let G (ℓ) be the geometry map associated with the patch ℓ defined in (3.23). We define the parameters (u i , v i ) associated with the point P i as

(u i , v i ) = G (ℓ) -1 (P ⊥ i );
This parameter computation is repeated for all the points of the point cloud. Figures 3.17 3.19, whose initial point cloud presents a very uneven distribution of points. Moreover, in all the numerical examples a regularization parameter λ = 10 -1 has been used in order to have a well-posed problem. The compression we obtain in these examples is smaller than the one in the experiments presented in Section 3.4.2; it could be increased doing more coarsening on the quad mesh but we didn't go further in this direction since the aim of the subsection is the general parametrization and mesh generation. Table 3

Comparison with C 0 fitting

Here we present a comparison between surfaces obtained by fitting a point cloud using our proposed G 1 basis functions and using the standard C 0 Bernstein basis on (the faces of) a quad mesh M, whose space is denoted by C 0 (M). In Table 3.7 we can notice that, despite the higher dimension of the continuous basis functions, the errors we get are comparable with the ones obtained from the G 1 construction using fewer basis functions; Chapter 4

dim C 0 (M) dim G 1 (M) L ∞ C 0 error L ∞ G 1 error RMSE C 0 RMSE G 1 f S 2
Analysis-suitable G 1 bases for isogeometric analysis

Complex geometries are commonplace when dealing with PDEs simulations arising from reallife applications; from this there is a need to have robust and analysis-suitable constructions to deal with such shapes. In this chapter we propose a new set of analysis-suitable G 1 spline basis functions for isogeometric analysis simulations on arbitrary domains Ω ⊂ R 3 . After recalling the generic constraints ensuring G 1 continuity between Bézier patches, we present the construction of vertex, edge and face functions and we provide a dimension formula for the spline space generated by such basis. Furthermore, we produce numerical experiments to ensure their suitability in IGA simulations on various PDEs. This chapter is based on [START_REF] Marsala | Analysis-suitable G 1 bases on quadrilateral meshes[END_REF].

Motivation

The basis functions introduced in Chapter 3 do not behave optimally when applied in the numerical resolution of PDEs on successive approximation of Catmull-Clark surfaces. This fact can be noticed from the plot in Figure 4.1, which reports the H 1 errors we obtain when solving the Poisson's equation (cf. (1.12) or following Section 4.4.2) over the domains obtained applying the G 1 construction introduced in Chapter 2 to the meshes in Figure 4.9. Indeed, except for the case of valence N = 4 that corresponds to the regular case, the other errors do not decrease, as expected in the optimal case, by a rate of 5. The reason behind this suboptimal result can be found in the refinement procedure we use during the simulations. By definition, the basis functions in Chapter 3 only allow the construction of a sequence of isogeometric spaces by iteratively splitting each quadrilateral mesh face regularly into four quadrilaterals (we opted for the Catmull-Clark split) and consequently reparametrizing via the G 1 framework in Chapter 2. It is easy to realize that the spaces obtained during this procedure are not nested. Moreover, this approach modifies the geometry, which changes at every refinement step. A rigorous explanation of why this suboptimal behaviour arises can be found in [START_REF] Takacs | Approximation properties over self-similar meshes of curved finite elements and applications to subdivision based isogeometric analysis[END_REF].

The goal of this chapter is to define a spline space and a basis over quad meshes allowing refinement techniques that don't change the geometry of the desired domain nor require a reparametrization at every step. This can be achieved, for example, and it is what we actually do, by equipping the construction with the possibility of using knot insertion. With this approach the resulting spline space will turn out to be nested and the domain's geometry is fixed.

General formulation of G 1 conditions

In order to describe our construction, we need to deal with G 1 constraints across Bézier patches defined over quad faces where, eventually, two linked EVs may occur. Actually, as will be explained in the next section, we will deal with EVs attached to valence 4 vertices whose edges are not all orthogonal to their adjacent ones; this vertex takes the name of non-crossing regular vertex (non-crossing RV in short) and, for our purposes, it will be treated like an EV. We assume again that the input quad mesh M presents isolated EVs. Similarly to the constructions in Chapter 2 and Chapter 3 we use here quadratic symmetric gluing data of the form

a N,N ′ (u) = a 0 B 0 2 (u) + a 1 B 1 2 (u) + a 2 B 2 2 (u), b N,N ′ (u) = -1,
where the B i d are the univariate Bernstein polynomials in (1.3) and the a i ∈ R, i = 1, 2, 3, are coefficients that we will define later on. As is well known, a Bézier patch is a particular case of tensor-product B-spline whose knot vectors are formed by zeros and ones only; for the biquintic case we are dealing with, the knot vector t B defining the grid t B × t B on which our patch lives is t B = {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1} = {0 6 , 1 6 }. i,j their control points defining the k-th patch labelled as in Figure 2.2, the system defining the G 1 constraints is defined by b

(1) 0,1 + b (0) 1,0 = (2 -a 0 )b (1) 0,0 + a 0 b (1) 1,0 , (4.2) 5(b (1) 1,1 + b (0)
1,1 ) = (a 0 -2a 1 )b

(1) 0,0 + (10 -5a 0 + 2a 1 )b

(1)

1,0 + 4a 0 b (1) 2,0 , (4.3) 10(b (1) 2,1 + b (0) 1,2 ) = (-a 0 + 2a 1 -a 2 )b (1) 0,0 + (5a 0 -10a 1 + a 2 )b (1)
1,0 + 2(-5a 0 + 4a 1 + 10)b

(1) 2,0

+ 3a 0 b (1)
3,0 , (4.4)

10(b (1) 3,1 + b (0) 1,3 ) = (a 0 -2a 1 + a 2 )b (1) 0,0 + 5(-a 0 + 2a 1 -a 2 )b (1) 1,0 + 2(5a 0 -10a 1 + 2a 2 )b (1) 2,0
+ 2(-5a 0 + 6a 1 + 10)b

(1) 3,0 + 4a 0 b

(1) 4,0 , (4.5)

5(b (1) 4,1 + b (0) 1,4 ) = (2a 1 -a 2 )b (1) 0,0 -(10a 1 -a 2 )b (1) 1,0 + 4(2a 1 -a 2 )b (1) 2,0 -2(10a 1 -3a 2 )b (1) 3,0 + 2(4a 1 + 5)b (1) 4,0 , (4.6) b (1) 5,1 + b (0) 1,5 = (2 + a 2 )b (1) 5,0 -a 2 b
(1) 4,0 , (4.7)

10(b (1) 3,0 -b (1) 2,0 ) = b (1) 0,0 -5b (1) 1,0 + 5b (1) 4,0 -b (1) 5,0 . (4.8)
These conditions, similarly to (3.9), can be reformulated in matrix form following the notation introduced in Section 3.2 as

M G 1 b = 0, (4.9)
where

M G 1 = (M 1 |M 2 ) is a block matrix with M 1 =          -(2 -a 0 )u C 1 -(a 0 -2a 1 )u (10 -5a 0 + 2a 1 )I -4a 0 I 5C 2 -u 5I 10I -10I (a 0 -2a 1 + a 2 )u -(5a 0 -10a 1 + a2)I 2(5a 0 -4a 1 -10)I -3a 0 I -(a 0 -2a 1 + a 2 )u 5(a 0 -2a 1 + a 2 )I -2(5a 0 -10a 1 + 2a 2 )I 2(5a 0 -6a 1 -10)I -(2a 1 -a 2 )u (10a 1 -a 2 )I -4(2a 1 -a 2 )I 2(10a 1 -3a 2 )I          M 2 =           10I 10I -4a 0 I 10I 10I -2(4a 1 + 5)I 5I 5I a 2 I -(2 + a 2 )I I I           and b the vector b = b 0,0 b 1,0 b 2,0 b 1,1 b 3,0 b 2,1 b 1,2 b 4,0 b 3,1 b 1,3 b 4,1 b 1,4 b 5,0 b 5,1 b 1,5 T , with b i,j = b (k) i,j , k = 1, . . . , N.
We notice that when a 1 = a 2 = 0, system (4.2) to (4.8) reduces to system (2.4) to (2.10). Given two biquintic patches defined over a generic tensor-product grid t × t, with t = {0 6 , t 1 ≤ t 2 ≤ . . . ≤ t s , 1 6 }, the equations defining the G 1 joins will depend also on the values of the knots (and their combinations), making them harder to be solved. In fact, calling with m the number of control points laying on the common edge shared by two patches, the G 1 equations on a general knot vector can be sketched as b

(1) 0,1 + b (0) 1,0 = (2 -a 0 )b (1) 0,0 + a 0 b (1) 1,0 , (4.10) b (1) 1,1 + b (0) 1,1 = L 1,1 (b (1) 0,0 , . . . , b (1)
m-1,0 ; t i , a 0 , a 1 , a 2 ), (4.11) b

(1)

2,1 + b (0) 1,2 = L 2,1 (b (1) 0,0 , . . . , b (1)
m-1,0 ; t i , a 0 , a 1 , a 2 ), (4.12) . . .

b

(1)

m-3,1 + b (0) 1,m-3 = L m-3,1 (b (1) 0,0 , . . . , b (1) m-1,0 ; t i , a 0 , a 1 , a 2 ), (4.13) b (1) m-2,1 + b (0) 1,m-2 = L m-2,1 (b (1) 0,0 , . . . , b (1)
m-1,0 ; t i , a 0 , a 1 , a 2 ), (4.14) b

(1)

m-1,1 + b (0) 1,m-1 = L m-1,1 (b (1) m-2,0 , b (1)
m-1,0 ; t i , a 2 ), (4.15) with t i ∈ t, i = 1, . . . , s, and where we indicate with L j,1 (• ; •) the linear combination of the variables present in the first argument by the coefficients present in the second argument concerning the j-th pair of control points. In addiction to (4.10) to (4.15) we also need to consider edge constraints, i.e. linear combinations of control points on the edge, which are as many as the knot spans present in t.

Let now assume that our knot vector is only formed by s uniform inner knots with multiplicity 5 (e.g. equal to the degree of the patch), that is

t = {0 6 , t 5 1 < t 5 2 < • • • < t 5 s , 1 6 }, (4.16)
where we followed the notation in (4.1). Under this hypothesis, the spline spaces associated to the knot vectors will be nested and, as consequence of the knot insertion algorithm, the patch will be split in several biquintic Bézier components which joins C 0 between them; moreover this split procedure gives rise to valence 4 points on the edge which can be seen as extraordinary vertices that are the so called non-crossing RV, in which the gluing data functions are nonzero. We will refer to them as virtual EVs. Finally, in order to have enough regularity in our space, we will impose extra C 1 continuity across the split lines of each patch. Figure 4.2 shows an example of patch with no split and with one split. Using (4.2) to (4.8) (e.g. system (4.9)) and the construction presented in Section 3.2, we provide in the next section the derivation of a set of basis function spanning the G 1 spline space which result to be analysis-suitable for IGA simulations.

Construction of the basis

We present here the construction of a set of basis functions generating the G 1 spline space over a quad mesh M associated to the knot vector t as in (4.16). We will refer to this set as B t , to stress its dependency on the knot vector. Similarly to what we did in Section 3.2, we will follow the topology of the input mesh M in order to separate the set of functions attached to (b): control points in a patch with a 4-split given by the knots t = {0 6 , 1/2 5 , 1 6 }. The red dot identifies the EV, the blue square the virtual EVs while the blue dots outline the C 0 split lines.

the vertices B t V , set of basis functions attached to the edges B t E and set of functions attached to the faces B t F . Hence, similarly to (3.8) we can decompose the set of basis functions as

B t = n V i=1 B t V i ∪ n E i=1 B t E i ∪ n F i=1 B t F i . (4.17)
Depending on the basis functions we want to compute, the strategy we exploit is different.

Regarding the vertex functions, their control points are obtained via knot insertion starting from the vertex basis functions defined on a Bézier knot vector, while the edge basis is computed concatenating around a virtual EV the relations (4.2) to (4.8), thanks to the particular shape of the patch we obtained under the assumption (4.16). The functions related to the faces are again the standard C 1 Bézier basis. Let investigate these constructions more in details.

Vertex basis functions: the set B t

V

To this set belong basis functions connected to inner and boundary EVs, regular vertices as well as corners. We only present the construction for functions supported around an EV since the other cases are analogous or already presented in Section 3.2.1.

Construction of basis functions corresponding to an inner EV

For an inner extraordinary vertex of valence N , also in this construction we have attached to it N + 3 different basis functions; this because the degrees of freedom of system (4.2) to (4.8) around an EV are the same as system (2.4) to (2.10). The computation of the basis functions attached to an inner EV defined on the knot vector t are obtained via knot insertion from the functions presented in Section 3.2.1.1. Let B EV k , k = 1, . . . , N + 3 be any of these basis functions (defined on t B ); if t (as in (4.16)) is our target knot vector, we can repeatedly apply the knot insertion algorithm to all the patches defining the basis function until we reach t. The output function B EV t k we get from this procedure presents a support which is bigger than desired; in fact, as we already noticed in the previous two sections, the G 1 constraints across an edge only involve the points along the edge and its above and below layers. In order to reduce the support of the latter function we apply truncation, i.e. we substract to the elevated basis function a linear combination of the C 1 functions living in the inner part of the patch. This procedure do not affect the G 1 regularity across edges. Figure 4.3 shows graphically the steps of this procedure on a patch. To conclude, we can summarize the entire procedure as follow: given an EV of valence N and any related vertex function B EV k , k = 1, . . . , N + 3 defined on t B , for suitable sets of indices I, I and I t , in order to derive the vertex function B EV t k defined on t we do:

B EV k = i,j∈I b i,j B i,j    knot insertion B EV t k = i,j∈ I b t i,j B t i,j    truncation B EV t k = i,j∈I t b t i,j B t i,j , where b t i,j = 0 if i, j ≥ 2, b t i,j = b t i,j
otherwise.

Figure 4.4 presents the set of vertex basis functions for an EV of valence N = 5 defined on the knot vector t = {0 6 , 1/2 5 , 1 6 }.

Bases linked to boundary EVs

In presence of a boundary EV of valence κ, with κ the number of patches sharing the vertex, similarly to the case of an inner EV we come up with κ + 3 basis functions. Since κ = N -1, it is equivalent to N + 2 and we can compute them applying the same procedure of repeated knot insertion as in Section 4.2.1.1.

Basis functions at regular vertices and corners

When dealing with a regular vertex, by construction we have attached to them standard C 1 basis functions; same fact for a corner. Both of them can be computed following the strategy in Section 3.2.1.2 and Section 3.2.1.3.

Edge basis functions: the set B t

E

We find in this set those functions whose support lays on the two patches shared by an edge, inner or boundary either extraordinary or regular. Here we present the explicit construction in the case of inner extraordinary since the remaining cases are analogous to Section 3.2.2.2 and Section 3.2.2.3.

Extraordinary edge basis functions

For a given knot vector t as in (4.16) and an extraordinary edge, we will investigate separately the construction of the edge functions which vanish along the edge from those that are nonzero along it. We assume at this step that our knot vector has only one multiple inner knot, i.e.

t = {0 6 , 1/2 5 , 1 6 }. (4.18)
The general situation will be treated later. To compute these edge basis functions we exploit the fact that knot vectors as (4.18) define an internal split of the patch in four Bézier subpatches. This split also creates across the edge what we call a virtual EV, i.e. a valence 4 point with no vanishing gluing data on it; hence we can circularly chain around this virtual EV, two patches by two, the relations (4.2) to (4.8), using each time the proper gluing data defined on the subedge identified by the split. In fact, along the extraordinary edge, at each split of the face corresponds a split of the gluing function defined on it: the coefficients identifying the new pair of subdivided (quadratic) gluing data can be computed by using the de Casteljau algorithm (cf. In addiction, since we want C 1 smoothness in the inner part of the face, when concatenating the equations between subpatches that do not cross the extraordinary edge we use (4.2) to (4.8) with null values for the gluing data coefficients that is in fact imposing C 1 continuity. The labelling for the coefficients of the new pair of gluing data function, as well as the subpatches, is made starting from the right outgoing edge of the virtual EV, in a counterclockwise ordering: we have a R 0 , a R 1 , a R 2 for the outgoing right edge and a L 0 , a L 1 , a L 2 for the outgoing left edge. It is straightforward that

a L 0 = -a R 0 , a R 0 ̸ = 0.
As we already explained, these coefficients can be computed applying the de Casteljau algoritmh to the initial gluing data function defined along the original edge. Figure 4.6-(c) illustrates the labeling of the new gluing data and subpatches orientation. Moreover, in each subpatch we order the control point following the Bézier structure in Figure 2.2 starting from the virtual EV which represent the point ḃ(k) 0,0 , k = 0, . . . , 3, where we introduced the over dot to differentiate the notation from the classical EV case. Figure 4.5 illustrates this labelling.

After building the system, the strategy we exploit to obtain the Bézier coefficients defining the basis functions is analogous to the strategy we used in Section 3.2.1: starting from (4.2) to (4.8) (or (4.9)) we impose, one at the time, the value one to each free coefficient appearing in the G 1 system. Consequently, we solve (4.2) to (4.8) with this initial values while gradually setting the value of any unconstrained coefficients that we encounter to zero.

Since we are dealing with edge functions, we can set to zero the values of the points neighbouring the vertices defining the edge e.g. the points ḃ(0)

4,0 , ḃ(0) 5,0 , ḃ(0) 4,1 , ḃ(0) 5,1 , ḃ(1) 0,4 , ḃ(1) 0,5 , ḃ(1) 1,4 , ḃ(1) 1,5 , ḃ(2) 4,1 , ḃ(2) 5,1 , ḃ(3) 1,4 , ḃ(3) 1,5
. The derivation of both types of basis functions is founded on two circulant systems, which come up after the concatenating procedure explained before, whose solutions return

0 1 2 3 ḃ(0) 0,0 ḃ(0) 1,0 ḃ(0) 2,0 ḃ(0) 3,0 ḃ(0) 3,1 ḃ(0) 2,1 ḃ(0) 1,1 ḃ(0) 0,1 ḃ(1) 1,1 ḃ(1) 1,2 ḃ(1) 1,3 ḃ(1) 0,1 ḃ(1) 0,2 ḃ(1) 0,3 ḃ(2) 1,1 ḃ(2) 2,1 ḃ(2) 3,1 ḃ(3) 1,0 ḃ(3) 1,1 ḃ(3) 1,2 ḃ(3) 1,3 Figure 4.5.
Labelling of the control points involved in the construction of the edge basis functions. The ordering is performed counterclockwise, starting from the patch zero, around the virtual EV (blue square). the points ḃ1,0 and ḃ1,1 (the bold refers to the vectorial notation introduced in Section 3.2 and (4.9)). These systems are:

(a) (b) 1 a L 2 a L 1 a L 0 2 a R 0 a R 1 a R 2 0 3 (c)
    -a R 0 1 0 1 1 0 1 0 0 1 a R 0 1 1 0 1 0     ḃ1,0 =     2 -a R 0 2 2 + a R 0 2     ḃ0,0 (4.19) and     1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 1     ḃ1,1 = d, (4.20)
where

d =          1 5 a R 0 -2a R 1 ḃ(0) 0,0 + 2 -a R 0 + 2 5 a R 1 ḃ(0) 1,0 + 4 5 a R 0 ḃ(0) 2,0 2 ḃ(1) 1,0 - 1 5 a R 0 + 2a L 1 ḃ(2) 0,0 + 2 + a R 0 + 2 5 a L 1 ḃ(2) 1,0 - 4 5 a R 0 ḃ(2) 2,0 2 ḃ(3) 1,0          . (4.21)
Let analyze their properties. If we call A ∈ R 4×4 the matrix in (4.19) it is trivial to notice that det(A) = 0 and corank(A) = 2.

In fact, it results that

Ker(A) = Span{ k1 , k2 }, where k1 = (-1, -a R 0 , 1, 0) T , k2 = (0, -1, 0, 1) T since A k1 = 0, A k2 = 0 and rank( k1 k2 ) = 2.
Next, if C is the matrix in (4.20), since we are in the case of even valence N = 4, from the analysis in Section 2.3.2.1 we know that

corank( C) = 1 and Ker( C) = Span{ w}, w = (1, -1, 1, -1) T .
From the previous investigation we can conclude that in order to obtain a solution, similarly to the constructions presented in Chapter 2 and Chapter 3, we need to add extra conditions to the systems in (4.19) and (4.20); again, these additional constraints will be the orthogonality of the desired solutions to the kernel of the corresponding matrix plus the verification of the Vertex Enclosure Relation (V.E.R.) (2.20). We underline that all the solutions we achieve must be meant in a least square sense, since we will deal with overdetermined systems. Following paragraphs show how to compute explicit solutions for problems (4.19) and (4.20).

Basis functions vanishing along the edge. To obtain basis functions which are zero along the edge we need to impose ḃ(0) i,0 = ḃ(2) i,0 = 0, i = 0, . . . , 5. (4.22) Specializing (4.9) for a virtual EV, that is valence N = 4 and properly subdivided gluing data, and substituting in it the assumption (4.22), we get that corank(M G 1 ) = 6, which means there are only six edge functions vanishing along the edge being G 1 across it and C 1 internally the faces for the knot vector t in (4.18). Four of them are the basis living in the inner part of the two subedges which derive from (4.4) to (4.5) and, alike to Section 3.2.2.1, we have

ḃ(0) 2,1 + ḃ(3) 1,2 = 0, ḃ(0) 3,1 + ḃ(3) 1,3 = 0, ḃ(1) 1,2 + ḃ(2) 2,1 = 0, ḃ(1) 1,3 + ḃ(2) 3,1 = 0,
that can be solved by setting, for example,

ḃ(0) 2,1 = ḃ(0) 3,1 = ḃ(1) 1,2 = ḃ(1) 1,3 = 1, ḃ(3) 1,2 = ḃ(3) 1,3 = ḃ(2) 2,1 = ḃ(2) 3,1 = -1.
Their graph is illustrated in Figure 3.4. Consequently the remaining two basis functions will be obtained by solving systems (4.19) and (4.20). If we substitute (4.22) in (4.19), it reduces to

1 1 1 1 ḃ(1) 1,0 ḃ(3) 1,0
= 0 from which we get the two (linearly independent) solutions

ḃ(1) 1,0 = ḃ(3) 1,0 = 0 (4.23) and ḃ(1) 1,0 = 1, ḃ(3) 1,0 = -1. (4.24)
To deduce the remaining control points we need to substitute, one at the time, previous solutions in (4.20). First, replacing (4.23) in (4.20) one has

C ḃ1,1 = 0 whence ḃ1,1 ∈ Ker( C) =⇒ ḃ1,1 = w,
which identifies the nonzero coefficients of the desired basis function. Going on, filling (4.24) in (4.20) we additionally have Nonzero edge basis functions. To extract this second set of functions we need to assume that there exists at least one nonvanishing control point along the edge. This assumption reflects an extra degree of freedom on (4.4) and (4.5), which can be solved symmetrically or not; in the light of Section 2.4.1 we opt for a symmetric solution. Under these hypothesis, from (4.9) adapted for a virtual EV, we obtain

      1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 1 1 -1 1 -1       ḃ1,1 = 2     0 -1 0 1     , (
corank(M G 1 ) = 2,
meaning that only two basis functions which don't vanish at the edge exist. The computation of the first basis function begins assuming that ḃ0,0 = 1 (and consequently ḃ(0) 3,0 = ḃ(2) 3,0 = 0); from (4.8) follows ḃ(0) 2,0 = -

1 10 + 1 2 ḃ(0) 1,0 , ḃ(2) 2,0 = - 1 10 + 1 2 ḃ(2) 1,0 . (4.26)
The analysis of (4.19) reveals that two extra conditions are required to ensure the uniqueness of its solution; the first additional constrain we impose is the feasibility of the points ḃ1,0 for the V.E.R. To do that, we derive the equation of the V.E.R from (4.21) forcing

⟨d| w⟩ = 0 which translates into 2 - 3 5 a R 0 + 2 5 a R 1 ḃ(0) 1,0 -2 ḃ(1) 1,0 + 2 + 3 5 a R 0 + 2 5 a L 1 ḃ(2) 1,0 -2 ḃ(3) 1,0 = 2 5 a R 1 + a L 1 , (4.27)
where we used (4.26). Adding (4.27) to (4.19) results

       -a R 0 1 0 1 1 0 1 0 0 1 a R 0 1 1 0 1 0 2 - 3 5 a R 0 + 2 5 a R 1 -2 2 + 3 5 a R 0 + 2 5 a L 1 -2        ḃ1,0 =        2 -a R 0 2 2 + a R 0 2 2 5 a R 1 + a L 1        (4.28)
and, if we name A V the matrix in (4.28), turns out that

corank(A V ) = 1, Ker(A V ) = Span{ kV } with kV = (0, -1, 0, 1) T .
Finally, the single solution is obtained after setting the extra condition

⟨ ḃ1,0 | kV ⟩ = 0 giving rise to          -a R 0 1 0 1 1 0 1 0 0 1 a R 0 1 1 0 1 0 2 - 3 5 a R 0 + 2 5 a R 1 -2 2 + 3 5 a R 0 + 2 5 a L 1 -2 0 -1 0 1          ḃ1,0 =          2 -a R 0 2 2 + a R 0 2 2 5 a R 1 + a L 1 0          (4.29) whose solution is ḃ1,0 = (1, 1, 1, 1) T . (4.30)
To obtain the following coefficients defining the nonzero basis function we substitute (4.30) in (4.26) and (4.21); then, adding the extra relation ⟨ ḃ1,1 | w⟩ = 0 to increase the rank of the matrix in (4.20) we get

      1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 1 1 -1 1 -1       ḃ1,1 =          2 - 12 25 a R 0 2 2 + 12 25 a R 0 2 0          obtaining the solution ḃ1,1 =          1 - 6 25 a R 0 1 + 6 25 a R 0 1 + 6 25 a R 0 1 - 6 25 a R 0         
.

The last coefficients concerning the higher order derivatives are extracted from (4.4) and (4.5), while after substituting in them the previous solutions we got they become

ḃ(0) 2,1 + ḃ(3) 1,2 = 4 5 1 - 3 5 a R 1 , ḃ(1) 1,2 + ḃ(2) 2,1 = 4 5 1 - 3 5 a L 1 and ḃ(0) 3,1 + ḃ(3) 1,3 = - 6 25 a R 2 , ḃ(1) 1,3 + ḃ(2) 3,1 = - 6 25 a L 2 .
To define the second nonzero basis functions we retrace the steps we made for the previous function by setting ḃ0,0 = ḃ(2) 3,0 = 0 but ḃ(0) 3,0 = 1; this choice updates the edge constraint (4.8) as ḃ(0)

2,0 = 1 2 ḃ(0) 1,0 + 1, ḃ(2) 2,0 = 1 2 ḃ(2) 1,0 . (4.31)
In this framework, the points ḃ1,0 needs to satisfy the V.E.R.

-

3 5 a R 0 + 2 5 a R 1 ḃ(0) 1,0 -2 ḃ(1) 1,0 + 2 + 3 5 a R 0 + 2 5 a L 1 ḃ(2) 1,0 -2 ḃ(3) 1,0 = - 4 5 a R 0 . (4.32)
Substituting the previous V.E.R. (4.32) in (4.29) and taking into account the initial values we fixed to extract this second basis function, we obtain

         -a R 0 1 0 1 1 0 1 0 0 1 a R 0 1 1 0 1 0 2 - 3 5 a R 0 + 2 5 a R 1 -2 2 + 3 5 a R 0 + 2 5 a L 1 -2 0 -1 0 1          ḃ1,0 =          0 0 0 0 - 4 5 a R 0 0          , from which ḃ1,0 = h     2 a R 0 -2 a R 0     , h = a R 0 8a R 0 -a R 1 + a L 1 , (4.33)
where h is well defined under the hypothesis of separated EV in the input mesh. Replacing the solution (4.33) and (4.31) in (4.21), and adding the additional constraint ⟨ ḃ1,1 | w⟩ = 0 as for the previous case, for the points ḃ1,1 we get from (4.20) the new system

      1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 1 1 -1 1 -1       ḃ1,1 = 2h           2 5 5 + 13 2 a R 0 + a L 1 a R 0 - 2 5 5 + 3 2 a R 0 + a L 1 a R 0 0           returning the vector ḃ1,1 = h 5     10 + 13a R 0 + 2a L 1 -(10 + 3a R 0 + 2a L 1 ) -(10 + 3a R 0 + 2a L 1 ) 10 + 13a R 0 + 2a L 1     .
Lastly, the values for the points linked to the higher derivatives can be deduced substituting (4.33), (4.31) and the initial assumption ḃ(0) 3,0 = 1 in (4.4) and (4.5) coming to Extraordinary edge functions for knot vectors with multiple inner knots. The previous paragraph shows the explicit construction of basis functions living across the two patches shared by an extraordinary edge when the knot vector present a single inner knot with multiplicity 5 as (4.18); the entire procedure led to reasonably simple closed formulas defining the control points of the different functions in Bézier form. In principle, the same pipeline can be applied for any knot vector with an arbitrary number of inner knots (4.16) by concatenating the G 1 constraints around every virtual EV appearing plus forcing C 1 smoothness in the inner part of the patch, but unfortunately the results are not obtained as easily as in the previous case. To explain why, let us consider a knot vector of type t = {0 6 , 1/3 5 , 2/3 5 , 1 6 }.

ḃ(0) 2,1 + ḃ(3) 1,2 = h a R 0 (90 -16a R 0 + 28a R 1 + a R 2 -2a L 1 )a R 0 + 4 a R 1 + 5 2 a L 1 -a R 1 , ḃ(1) 1,2 + ḃ(2) 2,1 = -h 10 -6a L 1 + a L 2 and ḃ(0) 3,1 + ḃ(3) 1,3 = h a R 0 -2(5 + 16a R 0 -2a R 1 + a R 2 + 2a L 1 )a R 1 + 2(5 + a R 2 )a L 1 + (80 + 13a R 2 )a R 0 , ḃ(1) 1,3 + ḃ(2) 3,1 = 3a L 2 h.
(4.34)

The case with s ≥ 3 inner knots follows directly. In presence of two inner knots our patch will be subdivided in 9 Bézier subpatches while the edge will be split in three subedges, giving rise to two virtual EVs; as consequence, also the gluing data defined over the original edge has to be redefined in each subedge via the de Casteljau algorithm leading to a triple of quadratic gluing data that is globally identified by 9 coefficients. That is why, being chained by these points, the two systems must be solved simultaneously and consequently each degree of freedom relative to a virtual vertex we need to fix to begin the computation will inevitability influences the points around the second virtual EV, and vice versa. Hence, the basis functions will have a support which lays all along the control points involved in the 8-times-chained system depending therefore on all the coefficients of the local gluing data functions defined along the edge. Despite all these dependencies (from both the elevated number of involved control points and gluing data coefficients), we are able to derive explicit solutions for the desired basis functions which, however, are not elegant and easy to deal with for implementation purposes. For this reason we opt for a numerical computation of these basis functions by numerically evaluating the kernel of the resulting chained system: this is made, as in the previous explicit construction, separating the computing of basis functions vanishing across the edge from by those that are nonzero on it. Additionally, as expected, computing the corank of the chained system we find that the zero functions along the edge are 10, while the nonzero are 4.

(a) (b) (c) (d)
The same chained dependency appears when dealing with knot vectors having an arbitrary number of inner knots: in this case, the pairing of the equations is made concatenating, twoby-two, the G 1 systems of two consecutive virtual EVs as well as the numerical computation. Moreover, if say our knot vector has s inner knots, we will end up with 4s + 2 zero edge basis functions and 2s nonzero edge functions along an extraordinary edge.

Face basis functions: the set B t

F

The only missing functions in this construction are those whose support only lay on a single patch. The concerned basis functions are C 1 and belong to this set also the ones obtained forcing the C 1 continuity across the split lines defined by the Bézier split. For a knot vector t with s inner knots, by using combinatorial formulas, we have (4s + 2) 2 basis functions in this set.

Analysis of the basis and space dimension

The functions we provided in previous Section 4.2 define a set of basis functions for the G 1 spline space over the mesh M and the knot vector t. We denote such space as G 1 t (M) = B t . The proofs of the following results demonstrating this statement, based on the decomposition in (4.17), are identical to those of Theorem 3.3.1 and Corollary 3.3.2 and for this they will be left out.

Theorem 4.3.1. The functions introduced in Section 4.2 form a basis B t for the space G 1 t (M) over a quad mesh M.

Corollary 4.3.2. The space G 1 t (M) has dimension given by:

dim G 1 t (M) = n V i=1 |B t V i | + n E i=1 |B t E i | + n F i=1 |B t F i | = n EV i=1 N EV i + 3n IEV + 2n BEV + 4(n RV + n C ) + 4(2s + 1)(n IRE + n BE ) + 2(3s + 1)n EE + (4s + 2) 2 n F ,
with s the number of inner knots in t, counted without multiplicity.

We observe that, as should be, for s = 0 we find again the dimension formula in Corollary 3.3.2.

Numerical experiments

Here we present some numerical experiment concerning the quality of the functions we propose when dealing with isogeometric analysis simulation. The spaces involved in the approximation process are a series of nested spaces defined over subsequent finer knot vector. Given an initial set of basis functions, the refined ones are obtained inserting 2 L -1, L ∈ N, uniform inner knots of multiplicity 5 in both parameter directions, where L represent the level of refinement.

Let us consider a mesh M composed of n F quad faces σ ℓ . The space of isogeometric functions we consider is given by

V := v ∈ L 2 (Ω) : v • G (ℓ) ∈ G 1 t (M) σ ℓ , ℓ = 1, . . . , n F , (4.35)
where G (ℓ) is the geometry map defining the G 1 ACC spline space introduced in Chapter 2 used to compute the ℓ-th patch of the physical domain

Ω = n F ℓ=1 Ω (ℓ) , Ω (ℓ) = G (ℓ) (σ ℓ ) . (4.36)
We parametrize each mesh's face such that σ ℓ = [0, 1] 2 ; moreover it results that V ⊂ H 2 (Ω).

The refined spaces will be denoted V h , where h is the size of the elements which will be properly specified in each example; generally we have h = O(2 -L ).

First we use the proposed functions for numerical experiments of L 2 projection for both planar and non planar domains and then continue in solving Poisson's equation and biharmonic equation.

L 2 projection

Let Ω be a multipatch domain as in (4.36), g : Ω -→ R a smooth function defined on it and let {ϕ i } i∈I , I = {1, . . . , dim( V h )} be a set of G 1 -smooth isogeometric functions forming a basis of the subspace V h ⊂ H 2 (Ω). We want to approximate the function g by the function

u h (x) = i∈I c i ϕ i (x), c i ∈ R, (4.37)
in the least square sense i.e. we compute the coefficients {c i } i∈I such that

∥ u h -g∥ 2 L 2 = Ω ( u h (x) -g(x)) 2 dx -→ min c i . (4.38)
Introducing the mass matrix M = (m i,j ) i,j∈I and the vector g = (g i ) i∈I whose entries are given by over three different domains obtained applying the G 1 ACC construction on the meshes in Figure 4.9 with valences N = 3, 4, 5. Table 4.1 shows the numerical errors for each step of subdivision while Figure 4.10 represents their plot, indicating that the convergence rate is optimal with respect to the L 2 norm, i.e. O(h 6 ). 4.12 present the numerical results we gained. We observe that also in this case we achieve the optimal convergence rate O(h 6 ) for the L 2 error.

m i,j = Ω ϕ i (x)ϕ j (x) dx and g i = Ω g(x)ϕ i (x) dx,
N = 3 N = 4 N = 5 L dim( V h ) L 2 error dim( V h ) L 2 error dim( V h ) L 2 error 0 72 1.586e-
Sphere Y -shape L dim( V h ) L 2 error dim( V h ) L 2 error 0 360 1.008e-

Poisson's equation

Let consider again a multipatch domain Ω. Here we test our basis functions on the following Poisson's problem

-∆u(x) = f (x), x ∈ Ω, u(x) = h D (x), x ∈ ∂Ω, (4.40)
with f, h D ∈ L 2 (Ω). If {ϕ} i∈I , I = {1, . . . , dim( V h,D )} is a set of basis functions for the subspace V h,D ⊂ H 1 (Ω), that is the subspace of functions in (4.35) verifying the Dirichlet boundary condition (4.40), using the Galerkin's approach in Section 1.5.2 we can translate the problem (4.40) in a system of linear equations

Kc = f
for the unknown coefficients c = (c i ) i∈I , where the entries of the stiffness matrix K = (k i,j ) i,j∈I and the vector f = (f i ) i∈I are defined, respectively, as

k i,j = Ω (∇ϕ i (x)) T ∇ϕ j (x) dx and f i = Ω f (x)ϕ i (x) dx,
such that our solution can be written as in (4.37). 4.3 together with their plot in Figure 4.13, indicating that we achieve optimal convergence rate of O(h 5 ) in H 1 norm.

N = 3 N = 4 N = 5 L dim( V h,D ) H 1 error dim( V h,D ) H 1 error dim( V h,D )

Biharmonic equation

We conclude this chapter presenting the numerical approach to a fourth-order problem, that is the biharmonic equation. Given a multipatch domain Ω, we can introduce the Cauchy problem for the biharmonic equation as

       ∆ 2 u(x) = f (x), x ∈ Ω, u(x) = h D (x), x ∈ ∂Ω, ∂u ∂ n = h N (x), x ∈ ∂Ω, (4.41)
where f, h D , h N are functions belonging to L 2 (Ω). As for the Poisson's problem, using the Galerkin's approach we can translate (4.41) in a linear system. To do that, let {ϕ i } i∈I be a set of basis functions for the subspace of functions in (4.35) verifying the Dirichlet-Neumann conditions (4.41) V h,D,N ⊂ H 2 (Ω). Introducing the matrix K = (k i,j ) i,j∈I and the vector f = (f i ) i∈I whose entries are the quantities

k i,j = Ω ∆ϕ i (x)∆ϕ j (x) dx and f i = Ω f (x)ϕ i (x) dx + ∂Ω h N (x)ϕ i (x) dγ,
we end up solving the linear system with unknowns c = (c i ) i∈I ,

Kc = f ,
allowing us to write the sought solution in the form of (4.37).

Example 4.4.4. We numerically solve the biharmonic equation (4.41) on the domains described by the meshes in Figure 4.9, defining the input functions f, h D and h N such that, similarly to the Poisson's equation, the exact solution is given by the function in (4.39). Table 4.4 resumes the H 2 errors while Figure 4.14 shows their convergence plot; also in this case we gain optimal convergence rate of O(h 4 ) in H 2 norm. As we already mentioned in Example 4.4.2, the proposed functions can be used to solve PDEs over nonplanar domains; in next Chapter 5 we will show some numerical experiments regarding this topic.

N = 3 N = 4 N = 5 L dim( V h,D,N ) H 2 error dim( V h,D,N ) H 2 error dim( V h,D,

Summary

This chapter has been devoted to the description of a set of analysis-suitable G 1 biquintic spline basis functions. Following the idea of Chapter 3, we proceeded with an extraction strategy starting from the equations describing the G 1 constraints between two adjacent spline patches defined on a knot vector whose knots presents multiplicity 5; this assumption allowed us to perform an easy and almost fully explicit construction for such functions. Several IGA simulations on classical PDEs have been provided to show the optimality of our novel proposal.

Chapter 5

A pipeline from CAD models to spline representation Shape modelling and analysis are crucial operations which directly impact engineering and industrial processes in many sectors of our society. The last several decades have witnessed the development of many powerful tools for computer-aided design (CAD), computer-aided engineering (CAE), and computer-aided manufacturing (CAM). These tools assist in handling the complex computations required to convert from the digital model of a shape to its actual production. Also boundary representation (B-rep) is a widely used tool when representing manufacturable geometry in Mechanical CAD (MCAD) processes. Such computations can include digital shape description, model reparation, meshing, numerical simulations, and optimisation. Currently, they require specific engineering efforts, are time consuming, and prone to errors and approximations: this explains why alternative approaches are preferred. Moreover, such model are not suitable for numerical simulations.

Using the constructions developed in Chapters 2 to 4 we present here a complete pipeline to convert CAD models into smooth G 1 spline representations, which are suitable for isogeometric analysis. Starting from a CAD boundary representation of a mechanical object, we perform an automatic control cage extraction by means of quadrangular faces, such that its limit Catmull-Clark subdivision surface approximates accurately the input model. Then we compute a basis of the G 1 spline space over the quad mesh in order to carry out least squares fitting over a point cloud, acquired by sampling the original CAD geometry. Finally, we use the basis functions to perform isogeometric analysis simulations of realistic PDEs on the reconstructed G 1 model. This chapter is based on [START_REF] Marsala | From CAD to Representations Suitable for Isogeometric Analysis: a Complete Pipeline[END_REF], and the work has been developed in collaboration with Sam Whyman and Mark Gammon from ITI.

Control cage generation from MCAD geometry

In this section, we detail the construction of the control cage from the boundary representation of a model. The foundation to our approach is to represent each MCAD edge with a cubic B-spline with a multiplicity-4 knot at each end. This allows the end curvature to be controlled by so-called slope control points which are not themselves part of the control cage topology (Figure 5.1). These can equally be thought of as tangent vectors stored at the ends of the spline. The limit surface is then defined to be the tensor product surface of the boundary B-splines. Since the slope control points affect only the first two knot spans, it follows that the first two layers of patches depart from the usual behaviour of the regular regions of a Catmull-Clark subdivision surface (i.e., away from the EVs). This modification constitutes the addition of Bézier edge conditions to the standard Catmull-Clark scheme [START_REF] Shen | Conversion of trimmed NURBS surfaces to Catmull-Clark subdivision surfaces[END_REF], and is illustrated by Figure 5.2. The 3D location of the limit surface at the boundary is influenced only by the control points on the boundary. Therefore, neighbouring control cages sharing the same boundary control points will maintain at least C 0 continuity between their limit surfaces, regardless of the positions of the interior control points, or the slope control points. The use of Bézier edge conditions imposes strict topological requirements upon the control cage. Namely, no EV may be placed on the boundary, nor within the first two layers of control cage faces, due to the tensor product nature of the limit surface within these regions. These topological constraints, to which engineering-grade subdivision surfaces must adhere, pose challenges when generating a suitable control cage. Each MCAD vertex must be treated as a corner, i.e., it must be associated with a 2-valent control cage vertex (such as in Figure 5.2). The appearance of any EV is restricted solely to the interior of the control cage. Meeting these requirements through traditional quad-dominant meshing techniques is challenging, and therefore we present a novel automatic approach referred to as SubD layering, which is tailored to satisfying the prescribed topological constraints. The process involves partitioning an MCAD face into regions of structured and unstructured mesh. The structured regions form a boundary layer from which EVs are fully excluded. These are formed by constructing 4-sided blocks around each MCAD vertex (referred to as corner blocks), and then connecting these to form 4-sided edge blocks associated with each MCAD edge. The remaining interior region constitutes a block with topology identical to the original MCAD face, and this is filled with unstructured mesh. This process is illustrated in Figure 5 the corner blocks are determined using the two-dimensional medial axis [START_REF] Ali | Block Topology Generation for Structured Multi-block Meshing with Hierarchical Geometry Handling[END_REF] as a guide to ensure that the blocks do not intersect. An example of this is given in Figure 5.4. Once the domain partitioning has been completed, we generate a coarse mesh of the layer faces. The boundary layer is meshed using a transfinite interpolation technique [START_REF] Gordon | Construction of curvilinear coordinate systems and application to mesh generation[END_REF] and is specified to be one element thick. The interior region is meshed using quad-dominant meshing technology [START_REF] Ali | Block Topology Generation for Structured Multi-block Meshing with Hierarchical Geometry Handling[END_REF]. The coarse mesh is then subdivided twice, and this ensures that there are no EVs within the first two layers of control cage faces from the boundary, thus satisfying the topological requirements for the Bézier edge conditions. This refinement of a coarse mesh is illustrated in Figure 5.5. The SubD layering process is applied to each MCAD face in the model in turn, such that the resulting control cage meets the topological requirements dictated by all MCAD edges, both external and internal (i.e., connectivity 1 and 2, respectively). By stipulating that the common vertices along the edges between contiguous faces are shared, a single limit surface is able to represent the entire model while remaining fully watertight (i.e. C 0 continuous). This guarantee is not maintained for MCAD B-rep geometry as a network of NURBS patches. An example of a watertight limit surface is given in Figure 5.6. An example of the control cage

Control cage adjustment

The SubD layering process outlined in Section 5.1 is primarily focused on achieving the correct topology for the control cage. The second step of the process is to adjust the control cage so that its limit surface coincides with the target MCAD geometry. To achieve this, we use the ability of our engineering-grade subdivision surfaces to accurately represent geometry; the behaviour of the edges is governed by the B-spline edge conditions, which can be made to respect the MCAD edges via a least squares fitting process. For the interior, we found that 1000 rounds of iterative control cage adjustment give a sufficiently good approximation for each of the four models. This is demonstrated in Figure 5.9, which shows the front section of the NASA CRM.

Point cloud sampling

The fitting procedure for our multipatch G 1 spline representation combines the generation of points on the geometry with a standard regression between the sample points and the parametric points on the subdivision surface. We require a set of samples points which lie exactly on the geometry, including the edges. Each point must be mapped to a corresponding patch of the subdivision surface, together with a local patch parameter coordinate. This is achieved by sampling each patch of the limit surface at a predetermined set of parameter values, using explicit evaluation [START_REF] Stam | Exact Evaluation of Catmull-Clark Subdivision Surfaces at Arbitrary Parameter Values[END_REF]. The corresponding 3D positions on the limit surface (which approximates the geometry) are then moved exactly onto the MCAD geometry. Given that the limit surface samples already lie very close to the target geometry, they may be projected easily onto the embedding NURBS surfaces of the MCAD B-rep faces. We also apply additional measures aimed at preventing creases and folds appearing in highly curved regions. The result is a uniform non-folding structured grid of sample points lying on the target geometry, for each patch of the limit surface.

From CAD to G 1

Here we report the numerical experiments concerning the different steps of the proposed construction. First we present the fitting procedure to reconstruct a G 1 surface from an MCAD point cloud using the technique introduced in Chapter 3, then we will use the previous result as the geometric domain over which to solve the heat equation using geometrically smooth basis functions in the IGA environment. The CAD models used for the numerical investigation are standard target examples for this type of problem. These are: Car model, Dream Chaser shuttle model, KCS hull model and the NASA CRM. All of them present special features and sharp edges which are a notable characteristic to be recovered in the fitted surface. The samplings of the original models are obtained following the procedure explained in Section 5.3. In order to obtain a precise result, the target point clouds contain large amounts of data. For the same reason, the control cage obtained from the MCAD (Section 5.1) presents a significant quantity of faces. The models are represented in their original scale, i.e., 1 unit = 1 m. All of the numerical experiments have been performed on two different machines: the first has been devoted to the control cage generation and the MCAD sampling, while the second ran the basis computation and consequent spline fitting and IGA simulation. Their specifications are: Windows 10, Intel(R) Xeon(R) CPU E3-1240 v6 @ 3.70GHz, 16.0 GB RAM, 4 cores, and Windows 10, Intel(R) Core(TM) CPU i7-9850 @ 2.60GHz, 16.0 GB RAM, 6 cores, respectively. Figures 5.12 to 5.15 show for each of the four models, in order, the quad mesh extracted from the input MCAD model by using the technique in Section 5.1, the point cloud acquired from the initial geometries following the idea in Section 5.3, the fitted surface in solid and multipatch coloring and a color plot representing the approximation error in Eucledian norm. Moreover, Table 5.1 summarises the dimensions of the starting point cloud, that is the number of points it contains indicated with n P , the number n F of faces forming the control cage and the number of basis functions n b generating the spline space together with the approximation errors evaluated with the formulas in (3.25). From the L ∞ errors represented in Table 5.1, as well as in the color map of the models' color plot, it can be noticed that the highest errors are located, as expected, around the EVs and near the sharp regions of the CAD models. This is because we are fitting sharp edges with high smoothness basis functions which cannot properly recreate the actual shape of the model in these regions (see Figure 5.11). In order to increase the quality of the fitting (i.e., decrease the error) and faithfully reproduce the characteristics of the input model, our construction allows us to identify the sharp edges of the CAD model which are to be preserved: this can be done by labelling as sharp the concerned edges and defining locally on those edges only C 0 basis functions. In this manner our output surface will manifest the desired sharp features. Figure 5.11 shows a detail of the KCS hull model where the top edge is first computed with G 1 smoothness, and then recomputed as a C 0 sharp feature. Regarding IGA simulations, various experiments with the heat equation have been performed: considering the following Cauchy problem

       ∂u ∂t (x, t) = c 2 ∆ Ω u(x, t) , (x, t) ∈ Ω × (0, T ] , u(x, 0) = u 0 (x) , x ∈ Ω , u(x, t) = u D (x, t) , (x, t) ∈ ∂Ω × (0, T ] ,
(5.1)

with Ω the physical domain we are looking at, ∆ Ω the Laplace-Beltrami operator defined in (1.23), c, T > 0 and u 0 (x), u D (x, t) given initial data, we run the simulations for each of the four CAD models, considering as the final time the instances T = 0.1, T = 0.2, T = 0.3 and T = 0.4 minutes. The Galerkin's discretization of (5.1) can be obtained in a similar way as (4.40). These time dependent integrations have been carried out using 20 time steps

Summary

In this chapter we provided a practical application of the constructions described in Chapters 2 to 4 to convert CAD models into G 1 objects. Starting from an input CAD model, using tools from the software CADfix have been possible to build a quad mesh such that its limit Catmull-Clark surface faithfully approximates the CAD model and also a point cloud discretization of the model itself. With this information, we were able to define the basis functions we previously developed on the given quad mesh and we used them to fit the cloud and obtained a smooth and analysis-suitable model on which we run IGA simulations.

Chapter 6

Free natural vibrations of a shallow lake

In this chapter we present another application of the novel constructions described in Chapters 2 to 4 concerning the IGA simulation of shallow-water equation for real case studies on lakes. Shallow-water equation is an important and popular tool to investigate the behaviour of water surfaces, i.e. their oscillations, in the shallow framework; in particular it is widely used to analyze the surface movement of rivers and lakes nearby their shores as well as coastal tsunami. As can be noticed from its formulation (1.17), the solution of the shallow-water equation strongly depends of the bathymetry function h which describes seabed's shape of the considered environment: for this reason, in order to have an accurate solution, having a precise representation of the bathymetry function is essential. The focus of our investigation is on the numerical simulation of the so-called free vibrations of a shallow lake, that are the oscillations of a lake's surface due only to the force of gravity. Since there are no analytic solutions for the general shallow-water equation, we first make sure the code is correct by analyzing a particular case where an explicit solution is available; then we move further to the general setting considering domains representing real lakes whose bathymetry function is given in terms of point cloud obtained from a GPS rastering. The steps to achieve a truthful solution are the following: starting from cloud data returning the bathymetry of a selected lake, we extract from it the curve delineating the shape of the lake. Then, using Rhino3D modeler we reconstruct a quad mesh (that discretize the lake's surface) on which we build our spline space using the bases in Chapter 3; these basis functions will be used both to reconstruct the bathymetry function by fitting the input point cloud and to solve the shallow-water equation. This is the reason why, in order to obtain a very precise solution, we are going to use fine quad meshes for both the fitting and the simulation.

Problem statement

Let Ω be our physical domain identifying the surface of a lake living on the plane x = (x, z) and h(x) the bathymetry function returning the bathymetry value y = -h(x). The Cauchy problem for the shallow-water equation can be formulated as

   h(x)∆u(x) + ∇h(x) • ∇u(x) + κ 2 u(x) = 0, x ∈ Ω, ∂u ∂ n (x) = 0, x ∈ ∂Ω, (6.1)
where κ 2 = σ 2 /g, g the gravitational acceleration and σ the frequency of the free oscillations.

The homogeneous Neumann condition is a natural choice for this kind of equations. We observe that (6.1) is actually describing an eigenvalue problem: introducing the differential operator D h , which depends on the function h, defined as

D h φ = h∆φ + ∇h • ∇φ, (6.2)
(6.1) can be reformulated as

D h u = -κ 2 u,
meaning that the solution we seek is an eigenfunction of the operator D h with eigenvalue -κ 2 . In general, there is no explicit solution for the problem (6.1) unless special cases such as very regular shape for the domain Ω and constant depth h; in order to test the correctness of our implementation we first focus on the particular case where the lake Ω is a square domain of length ℓ with depth h = constant.

The test case

Assuming constant bathymetry, (6.2) reduces to

D h φ = h∆φ, from which (6.1) becomes    ∆u(x) = -κ 2 u(x), x ∈ Ω, ∂u ∂ n (x) = 0, x ∈ ∂Ω, (6.3)
with κ2 = κ 2 /h, that is a well-known equation in acoustic applications named Helmholtz equation. In the additional case of a square lake, an exact solution of (6.3) can be simply obtained via separation of variables, leading to:

u(x, z) = Λ cos kπx ℓ cos jπz ℓ , (x, z) ∈ Ω, (6.4)
with Λ ∈ R and k, j, ℓ any integers such that κ2 = π 2 k 2 + j 2 ℓ 2 .

Figure 6.1 shows the exact solution (6.4) in presence of a unit square (i.e. ℓ = 1) for k = j = 1, k = 2, j = 3 and Λ = 1. Similarly to the problems in Section 4.4, by using the basis functions introduced in Chapter 4, the Galerkin's discretization of (6.3) brings to

Kc = -κ 2 M, (6.5)
where the stiffness K = (k i,j ) i,j∈I and mass M = (m i,j ) i,j∈I matrices have entries

k i,j = Ω (∇ϕ i (x)) T ∇ϕ j (x) dx, m i,j = Ω ϕ i (x)ϕ j (x) dx (6.6)
and c = (c i ) i∈I is the vector containing the unknown coefficients of the desired solution expressed as linear combination of the bases {ϕ i } i∈I as in (4.37). Hence, our numerical investigation is devoted to the computation of the eigenvalues of the generalized eigenvalue problem defined in (6.5) whose related eigenfunctions are the solution of (6.3).

The general case

Once we make sure that our implementation is correct, we can deal with the general problem (6.1) defined over a domain Ω with arbitrary shape and any bathymetry function h. As we said above, in this setting the problem (6.1) has no explicit solution. Its Galerkin's formulation, as for the test case, leads to the generalized eigenvalue problem Kc = -κ 2 M, with M as in (6.6) and stiffness matrix K = ( k i,j ) i,j∈I

k i,j = Ω h(x) (∇ϕ i (x)) T ∇ϕ j (x) dx + Ω ϕ j (x) (∇h(x)) T ∇ϕ i (x) dx
In the next section we present two numerical simulation for the shallow-water equations based on real data: starting from the 3D point cloud of a lake's seabed obtained from a GPS rastering we are able to delineate the shape of the lake (defining our domain Ω) as well as the h function via fitting the cloud using the basis functions in Chapter 3 defined over a quad meshing of Ω.

Simulations on real lake data

The data used in these experiments is taken from the dataset analyzed in [START_REF] Khazaei | GLOBathy, the global lakes bathymetry dataset[END_REF]. In particular we use data from the Rogagua Lake, in Bolivia and from the Orta lake, in Italy; the two lakes have been chosen based on their very different shape and bathymetry conformation.

In order to solve the shallow-water problem (6.1) using the technique introduced in Chapter 4 we need to define a mesh domain delineating the shape of the lake as well as a function returning the bathymetry of the lake. This goal can be achieved by following the next steps: first, starting from the point cloud acquired from a GPS rastering of the lake's seabed (represented as triangular mesh), using the Rhino3D modeler command ExtractMeshEdges we extract the boundary edges and we get a curve identifying the shape of the lake. Since the extrapolated curve is obtained from a set of points it is a piecewise linear curve; therefore, to achieve a smoother and more realistic result we apply few iterations of Laplacian smoothing [START_REF] Sorkine | Laplacian Mesh Processing[END_REF]. This curve will identify our physical domain Ω. A triangular meshing of the internal area can be made by the command PlanarMesh of Rhino3D, and its conversion to quads by using the QuadRemesh option. As accurately replicating the lake's shape is essential for a precise simulation, the resulting mesh will exhibit a high number of faces. Finally, on this mesh we build our spline bases which will be used to both reconstruct the bathymetry function h by fitting the input point cloud and solving the shallow-water equation.

Rogagua lake

The first experiment is performed on the Rogagua lake model from the repository analyzed in [START_REF] Khazaei | GLOBathy, the global lakes bathymetry dataset[END_REF]. Figure 6.2 illustrates the meshing and fitting steps properly explained in Section 6.2 while Table 6.2 summarizes the features of the input point cloud, i.e. its number of points n P , the number of quad faces n F of the extrapolated quad mesh and the number of bases n b defined on it as well as the errors (defined in (3.25)) while fitting the bathymetry function. It also contains the dimensions of the lake's model bounding box, expressed in meters. The error color plot in Figure 6.2-(e) points-out that the (very few) regions where the fitting error is localized corresponding to the boundary of the domain, and this is due to its slight modification induced by the Laplacian smoothing in order to get a proper outline. Finally, we computed the first 15 solutions for the shallow-water problem (6.1), shown in Figure 6.3, that correspond to the eigenfunctions with associated eigenvalues (and frequencies) listed in Table 6.1.

Eigenvalues

Frequencies 2.872e-08 2.802e-05 1.601e-04 5.307e-04 1.658e-02 3.962e-02 3.713e-04 5.037e-04 5.087e-04 6.034e-02 7.028e-02 7.063e-02 6.527e-04 8.495e-04 9.009e-04 8.001e-02 9.128e-02 9.399e-02 1.150e-03 1.294e-03 1.489e-03 1.060e-01 1.127e-01 1.208e-01 1.573e-03 1.838e-03 1.878e-03 1.242e-01 1.343e-01 1.357e-01

Table 6.1. First 15 eigenvalues, and corresponding frequencies, for the shallow-water simulation on Rogagua lake.

Orta lake

The second investigation concerns the Orta lake, whose model belongs to the same dataset as before. Its irregular shape and seabed is very different from the previous Rogagua lake, and together with the presence of an island make this case very interesting. Similarly to the previous experiment, Table 6.2 presents the point cloud features and the attributes of the quad mesh we achieved from it. Also in this case the fitting error is localized nearby the boundary of the domain. Table 6.3 lists the first 15 eigenvalues and frequencies identifying the 15 solutions of the shallow-water equation depicted in Figure 6.5. Chapter 7

Cubic C 2 spline quasi-interpolants on arbitrary triangulations

Quasi-interpolation is a general approach to construct efficient local approximants in a prescribed space to a given function f or just a given set of data. It is often preferred over interpolation because it is a local method and does not require the solution of possibly large linear systems. A quasi-interpolation operator is usually obtained as a linear combination of a suitable set of functions, which are positive, stable, and locally supported, or at least have a strong fast decay, in order to achieve local control.

In this chapter we present the construction of C 2 cubic spline quasi-interpolants on a given arbitrary triangulation where each triangle is refined according to the cubic Wang-Shi (WS) split. The main tool in our construction is the simplex spline basis for the local space of C 2 cubics on the cubic WS split presented in [START_REF] Lyche | Construction of C 2 Cubic Splines on Arbitrary Triangulations[END_REF]. In order to guarantee global C 2 smoothness between the local simplex spline representations, our construction starts from a Hermite interpolation problem that uniquely identifies elements of the cubic WS spline space on a prescribed triangulation. For a given function f , different C 2 cubic spline quasi-interpolants are then obtained by feeding different sets of Hermite data to this Hermite interpolation problem. Finally, various numerical examples illustrating the performance of the introduced quasi-interpolants are presented. This chapter is based on [START_REF] Marsala | Maximally smooth cubic spline quasiinterpolants on arbitrary triangulations[END_REF].

7.1 The spline space S 2 3 (△ WS 3 )

Let us consider a reference (macro-)triangle △ = ⟨p 1 , p 2 , p 3 ⟩ identified by three noncollinear points p 1 , p 2 , p 3 in R 2 . We divide each edge of △ into three equal segments, respectively, resulting into nine boundary points. Then, we refine △ into a number of subelements delineated by the complete graph connecting those boundary points, which consists of the three edges of the triangle and 18 interior lines. The resulting partition is called the WS 3 (or cubic Wang-Shi) split of △ as it is the third member of a family of splits originally proposed by Wang and Shi [START_REF] Wang | S µ µ+1 surface interpolations over triangulations[END_REF]; see Figure 7.1 for an illustration. We denote the obtained mesh structure by △ WS 3 : it has 58 vertices and 75 polygonal regions. We consider the space of C 2 cubic splines on the partition △ WS 3 , i.e., S 3]: sequences of knots for a set of simplex spline basis functions for S 2 3 (△WS 3 ). Each black disc shows the position of a knot, and the number inside indicates its multiplicity.

Hermite interpolation in S 2

3 (△ WS 3 )

An efficient way to identify the elements of the space S 2 3 (△ WS 3 ) is to solve a suitable Hermite interpolation problem. Taking into account the dimension of the space, a natural choice is to consider six Hermite data values associated with each vertex of △, three Hermite data values associated with each edge of △, and one function value associated with the triangle △ itself. Due to its importance later on, we explicitly state the following result which is a byproduct of [LMS22, Theorem 3]. Proposition 7.2.1. Let f j ∈ R, j = 1, . . . , 28, be given. There exists a unique spline s ∈ S 2 3 (△ WS 3 ) such that ρ j (s) = f j , j = 1, . . . , 28, (7.2)

where the operators ρ 1 , . . . , ρ 18 related to the vertices of △ are defined as

ρ 1 (φ) = φ(p 1 ), ρ 2 (φ) = φ(p 2 ), ρ 3 (φ) = φ(p 3 ), ρ 4 (φ) = D p 1 p 2 φ(p 1 ), ρ 5 (φ) = D p 1 p 3 φ(p 1 ), ρ 6 (φ) = D p 2 p 3 φ(p 2 ), ρ 7 (φ) = D p 2 p 1 φ(p 2 ), ρ 8 (φ) = D p 3 p 1 φ(p 3 ), ρ 9 (φ) = D p 3 p 2 φ(p 3 ), ρ 10 (φ) = D 2 p 1 p 2 φ(p 1 ), ρ 11 (φ) = D 2 p 1 p 3 φ(p 1 ), ρ 12 (φ) = D 2 p 2 p 3 φ(p 2 ), ρ 13 (φ) = D 2 p 2 p 1 φ(p 2 ), ρ 14 (φ) = D 2 p 3 p 1 φ(p 3 ), ρ 15 (φ) = D 2 p 3 p 2 φ(p 3 ), ρ 16 (φ) = D p 1 p 2 D p 1 p 3 φ(p 1 ), ρ 17 (φ) = D p 2 p 3 D p 2 p 1 φ(p 2 ), ρ 18 (φ) = D p 3 p 1 D p 3 p 2 φ(p 3 ), (7.
3) the operators ρ 19 , . . . , ρ 27 related to the edges of △ are defined as

ρ 19 (φ) = D q 3 p 3 φ(q 3 ), ρ 20 (φ) = D q 1 p 1 φ(q 1 ), ρ 21 (φ) = D q 2 p 2 φ(q 2 ), ρ 22 (φ) = D 2 p 3,1 p 3 φ(p 3,1 ), ρ 23 (φ) = D 2 p 2,1 p 2 φ(p 2,1 ), ρ 24 (φ) = D 2 p 1,2 p 1 φ(p 1,2 ), ρ 25 (φ) = D 2 p 3,2 p 3 φ(p 3,2 ), ρ 26 (φ) = D 2 p 2,3 p 2 φ(p 2,3 ), ρ 27 (φ) = D 2 p 1,3 p 1 φ(p 1,3 ), (7.4)
and finally the operator ρ 28 related to the triangle △ is defined as

ρ 28 (φ) = φ(q), (7.5)
with (see also Figure 7.3)

p 1,2 = 2 3 p 2 + 1 3 p 3 , p 1,3 = 1 3 p 2 + 2 3 p 3 , q 1 = 1 2 p 2 + 1 2 p 3 , p 2,1 = 2 3 p 1 + 1 3 p 3 , p 2,3 = 1 3 p 1 + 2 3 p 3 , q 2 = 1 2 p 1 + 1 2 p 3 , p 3,1 = 2 3 p 1 + 1 3 p 2 , p 3,2 = 1 3 p 1 + 2 3 p 2 , q 3 = 1 2 p 1 + 1 2 p 2 , (7.6
) For the sake of completeness, we report in Table 7.1 the collocation matrix of the operators ρ j applied to the basis functions B i , i, j = 1, . . . , 28; see also [START_REF] Lyche | Construction of C 2 Cubic Splines on Arbitrary Triangulations[END_REF]Section A.3]. From [LMS22, Corollary 1] we know the following related result.

and q = 1 3 p 1 + 1 3 p 2 + 1 3 p 3 . (7.7) p 1 p 3,1 p 3,2 p 2 p 1,2 p 1,3 p 3 p 2,3 p
Corollary 7.2.2. For given data f k,α,β , g k , g k,l , and h 0 , there exists a unique spline s ∈ S 2 3 (△ WS 3 ) such that

D α x D β y s(p k ) = f k,α,β , 0 ≤ α + β ≤ 2, k = 1, 2, 3, D n k s(q k ) = g k , D 2 n k s(p k,l ) = g k,l , k, l = 1, 2, 3, k ̸ = l, s(q) = h 0 , (7.8)
where n k is the outgoing normal direction of the edge opposite to the vertex p k , and the points p k,l , q k and q are defined in (7.6) and (7.7).

7.3

The spline space S 2 3 (T WS 3 )

Let T be a triangulation of a polygonal domain Ω ∈ R 2 and let T WS 3 denote its refinement obtained by taking the WS 3 split of each of its triangles. We consider the space of C 2 cubic splines on T WS 3 , i.e.,

S 2 3 (T WS 3 ) := s ∈ C 2 (Ω) : s τ ∈ P 3 , τ is polygon in T WS 3 .
The unisolvency of the Hermite interpolation problem stated in Corollary 7.2.2 implies that

dim(S 2 3 (T WS 3 )) = 6n V + 3n E + n T ,
where n V , n E , n T are the number of vertices, edges, and triangles of T , respectively; see [LMS22; Wan01; WS90]. Moreover, any spline function of S 2 3 (T WS 3 ) can be locally constructed on each (macro-)triangle △ i of T via the Hermite data (7.8), and the corresponding spline piece on △ i can be represented in the form

s △ i = 28 j=1 c i j B i j , c i j ∈ R, (7.9)
where {B i 1 , . . . , B i 28 } is the scaled simplex spline local basis (7.1) for the triangle △ i . Any function, which is represented locally in the form (7.9) on each Then, to obtain a consistent value for f 19 , we only need to compute the value D p 1 p 2 s(q 3 ). This can be achieved by observing that, since s ∈ S 2 3 (△ WS 3 ), its restriction to the edge p 1 p 2 is a univariate C 2 cubic spline with two inner knots at p 3,1 and p 3,2 , respectively. Thus, it is uniquely determined by its values and first and second derivatives at the two ends of the edge, i.e., by the values of ρ 1 (s), ρ 2 (s), ρ 4 (s), ρ 7 (s), ρ 10 (s), ρ 13 (s) or, in other words, by the values f 1 , f 2 , f 4 , f 7 , f 10 , f 13 previously computed. This allows us to compute D p 1 p 2 s(q 3 ) and so to get a consistent value for f 19 . Similarly, we can compute f 20 and f 21 . Now, let us consider f 22 . The direction p 3,1 p 3 involved in the definition of the operator ρ 22 (see Figure 7.5) can be written as p 3,1 p 3 = ηp 1 p 2 + ζn 3 , for some real values η, ζ. Thus, we have

△ i of T , is C 2 -smooth over △ i but not necessary C 2 -
ρ 22 (s) = D 2 p 3,1 p 3 s(p 3,1 ) = η 2 D 2 p 1 p 2 s(p 3,1 ) + 2ηζD p 1 p 2 D n 3 s(p 3,1 ) + ζ 2 D 2 n 3 s(p 3,1 ).
Hence, from (7.8) we set

f 22 = η 2 D 2 p 1 p 2 s(p 3,1 ) + 2ηζD p 1 p 2 D n 3 s(p 3,1 ) + ζ 2 g 3,1 .
The value D 2 p 1 p 2 s(p 3,1 ) can be computed from the analytical expression of s along the edge p 1 p 2 previously obtained. Finally, the restriction of D n 3 s along the edge p 1 p 2 is a univariate C 1 quadratic spline with two inner knots at p 3,1 and p 3,2 , respectively. Therefore, such univariate spline is uniquely determined by its values and first derivatives at the two ends of the edge, which can be computed from

f k,α,β , 1 ≤ α + β ≤ 2, k = 1, 2,
and by its value at the midpoint q 3 , i.e., by g 3 . Similarly, we can compute f j , j = 23, . . . , 27.

It is worth to notice that the computed values f j , j = 1, . . . , 28, are linear expressions in terms of the given Hermite data (7.8).

Quasi-interpolant from exact Hermite data

For this construction, we assume that for each triangle △ i of T , i = 1, . . . , n T , the Hermite data f i k,α,β , g i k , g i k,l , and h i 0 , as in (7.8), are sampled from a given (sufficiently smooth) function f . These data allow us to compute the values f i exa = (f i exa,1 , . . . , f i exa,28 ) T (7.10) according to the strategy outlined in Section 7.4.1. Then, to find the spline satisfying the problem stated in (7.2), we solve the linear system

B µ i exa (f ) = f i exa ,
where the matrix B ∈ R 28×28 is given by

B = (b l,k ) 28 l,k=1 = ρ l (B i k ) 28 l,k=1 . (7.11)
Note that, due to the definition of ρ l , l = 1, . . . 28, and B i k , k = 1, . . . , 28, the matrix B is independent of i. Actually, it is just the transpose of the matrix reported in Table 7.1. It can be directly checked that

∥B∥ ∞ = 366, ∥B -1 ∥ ∞ = 3523 729 ≤ 5. (7.12) Let µ i exa (f ) = (µ i exa,1 (f ), . . . , µ i exa,28 (f )) T , we define Q exa : C 2 (Ω) -→ S 2 3 (T WS 3 ) by Q exa f △ i = 28 j=1 µ i exa,j (f )B i j , i = 1, . . . , n T . (7.13)
By construction, the quasi-interpolant Q exa f enjoys the global smoothness, cubic polynomial reproduction, and locality properties stated at the beginning of Section 7.4. Moreover, we have the following result as a consequence from Corollary 7.2.2. Proposition 7.4.1. The quasi-interpolant in (7.13) is a projector onto S 2 3 (T WS 3 ), i.e.,

Q exa s = s, ∀s ∈ S 2 3 (T WS 3 ).
Let ∥f ∥ ∞,Ω be the standard L ∞ -norm of the function f on the domain Ω,

∥D k f ∥ ∞,Ω = max α+β=k ∥D α x D β y f ∥ ∞,Ω ,
and let h T be the length of the longest edge of the triangulation T with θ T its smallest angle. Then, from its construction described in Section 7.4.1, we know that f i exa in (7.10) can be bounded by

∥f i exa ∥ ∞ ≤ C ∥f ∥ ∞,Ω + h T ∥Df ∥ ∞,Ω + h 2 T ∥D 2 f ∥ ∞,Ω , (7.14)
for some constant C, which is independent of f and h T but depends on θ T .

be the known values of a given function f . As in Section 7.4.3, we build different polynomials related to the features of the considered triangulation, namely its vertices, edges, and faces. Let v k , k ∈ {1, . . . , n V }, be a vertex of T and denote by |v k | the number of triangles attached to v k . We define the weighted vertex polynomial V k ε f as

V k ε f = 0≤α+β≤3 c k α,β x α y β , (7.22)
where c k = (c k 0,0 , c k 0,1 . . . , c k 3,0 ) T is obtained by solving the overdetermined linear system

W k V k c k = W k z k
in least-squares sense. Here V k ∈ R 10|v k |×10 is the Vandermonde matrix in terms of the monomial basis of P 3 evaluated at all the modified Bézier domain points of the triangles sharing v k , and z k is the vector containing the values

z k = { z i ℓ , i : v k ∈ △ i , ℓ = 1, . . . , 10}.
Moreover, W k ∈ R 10|v k |×10|v k | is the diagonal matrix containing the weights corresponding to the points in (7.21). More precisely, for a given point b i ℓ that belongs to a triangle △ i attached to v k , its weight ω( b i ℓ ; v k ) is chosen as

ω( b i ℓ ; v k ) = 1 1 + 10 dist( b i ℓ ,v k ) h △ i 2 , (7.23)
where dist( b i ℓ , v k ) = ∥ b i ℓv k ∥ 2 is the standard Eucledian distance and h △ i is the length of the longest edge of △ i . A similar procedure can be followed to define the weighted edge polynomial E j ε f attached to the edge e j of T , j ∈ {1, . . . , n E }, E j ε f = 0≤α+β≤3 c j α,β x α y β , (7.24)

where c j = (c j 0,0 , c j 0,1 . . . , c j 3,0 ) T is determined by solving the overdetermined linear system W j V j c j = W j z j in least-squares sense. Here V j ∈ R 20×10 is the Vandermonde matrix in terms of the monomial basis of P 3 evaluated at all the modified Bézier domain points of the triangles sharing e j , and z j is the vector containing the values z j = { z i ℓ , i : e j ∈ △ i , ℓ = 1, . . . , 10}.

The weights in the diagonal matrix W j ∈ R 20×20 are calculated similar to (7.23). For a given point b i ℓ that belongs to a triangle △ i attached to e j , we choose the weight ω( b i ℓ ; e j ) as ω( b i ℓ ; e j ) = 1 1 + 10 dist( b i ℓ ,e j ) h △ i 2 .

Lastly, the triangle (or face) polynomial F i ε f attached to the triangle △ i , i ∈ {1, . . . , n T },

F i ε f = 0≤α+β≤3 c i α,β x α y β , (7.25)
is defined as the cubic polynomial interpolating at the points in (7.21), which can be built in a manner similar to (7.16). Figure 7.7-(b) illustrates the positioning of the above weighted polynomials defined in (7.22), (7.24), and (7.25). These polynomials can take the same role as the averaged polynomials introduced in Section 7.4.3 and can be used to recover the Hermite values involved in Corollary 7.2.2. The new quasi-interpolant Q wls f is built as follows. Consider the triangle △ i of T , i ∈ {1, . . . , n T }. With reference to (7.8), from the weighted vertex polynomials V k ε f in (7.22) we sample the values f i k,α,β ; from the weighted edge polynomials E j ε f in (7.24) we sample the values g i k and g i k,l ; and finally, the value h i 0 is sampled from the triangle polynomial F i ε f in (7.25). Given these Hermite data we follow the same construction as in Section 7.4.2. More precisely, applying the conversion strategy explained in Section 7.4.1, we find the values f i wls = (f i wls,1 , . . . , f i wls,28 ) T , from which we compute the coefficients µ i wls (f ) = (µ i wls,1 (f ), . . . , µ i wls,28 (f )) T by solving B µ i wls (f ) = f i wls , with B as in (7.11). The quasi-interpolant Q wls : C 0 (Ω) -→ S 2 3 (T WS 3 ) is then identified by

Q wls f △ i = 28 j=1
µ i wls,j (f )B i j , i = 1, . . . , n T . (7.26)

Once more, by construction, the quasi-interpolant Q ave f enjoys the global smoothness, cubic polynomial reproduction, and locality properties stated at the beginning of Section 7.4.

Error estimates

In this section we analyze the approximation order of the proposed quasi-interpolants to a sufficiently smooth function f . To this end, we recall a classical approximation result about polynomials on triangles; see, e.g., [LS07, Theorem 1.3]. Given a triangle △ with maximal edge length h △ , for every f ∈ C d+1 (△) there exists a polynomial p d ∈ P d such that

∥D k (f -p d )∥ ∞,△ ≤ K d h d+1-k △ ∥D d+1 f ∥ ∞,△ , (7.27)
for all 0 ≤ k ≤ d. The constant K d depends only on d. Note that a similar estimate in L q -norm can be shown for any function f belonging to the Sobolev space W d+1,q (△); see [START_REF] Lai | Spline Functions on Triangulations[END_REF].

Proposition 7.4.3. Given a function f ∈ C 4 (Ω), the quasi-interpolant in (7.13) satisfies ∥f -Q exa f ∥ ∞,Ω ≤ K exa h 4 T ∥D 4 f ∥ ∞,Ω , (7.28)

where the constant K exa is independent of f and h T .

Proof. Consider a triangle △ i of T . Then, using a triangle inequality,

∥f -Q exa f ∥ ∞,△ i ≤ ∥f -p 3 ∥ ∞,△ i + ∥p 3 -Q exa f ∥ ∞,△ i ,
where p 3 is the same polynomial as in (7.27) for d = 3. Moreover, the polynomial reproduction property of the quasi-interpolant implies ∥p 3 -Q exa f ∥ ∞,△ i = ∥Q exa (p 3f )∥ ∞,△ i and from Proposition 7.4.2 we deduce

∥p 3 -Q exa f ∥ ∞,△ i ≤ 5C ∥f -p 3 ∥ ∞,△ i + h T ∥D(f -p 3 )∥ ∞,△ i + h 2 T ∥D 2 (f -p 3 )∥ ∞,△ i .
Combining the above estimates with (7.27) results in (7.28), where K exa = K 3 (1 + 15C).

Similar approximation estimates can be derived for the quasi-interpolants Q ave f in (7.20) and Q wls f in (7.26) as well.

Numerical experiments

This section contains several numerical experiments illustrating the quality of the proposed quasi-interpolants Q exa f , Q ave f , and Q wls f , described in Section 7.4. We start by applying our constructions to a degree four polynomial and the well-known Franke's function, both defined on a square domain. We consider a sequence of nested tridirectional meshes to verify the approximation order of the quasi-interpolants (see Section 7.4.5). Next, we apply our constructions to a sigmoid function and a function with three peaks, again defined on a square domain, but now we use adaptively refined triangular meshes in order to capture more accurately the local features. Lastly, we present an example dealing with a non-square domain. In all the examples we take the value ε = 1/2 in (7.21) for the construction of Q wls f .

Being able to evaluate any spline s ∈ S 2 3 (△ WS 3 ) represented in the basis (7.1) is an important aspect in our constructions. Since the basis functions are (scaled) simplex splines, they could be evaluated through recurrence or simply by means of a lookup-table process; see [LMS22, Section 5.1] and in particular [LMS22, Tables 1 and2] for their explicit expressions. Any point p ∈ △ can be uniquely characterized and searched via a Boolean vector. Indeed, thanks to the cross-cut nature of the WS 3 split, each of the 75 polygonal regions in △ WS 3 is uniquely identified by the sign of the linear equations defining the 18 interior lines of the split; see [LMS22, Eq. ( 28)]. Hence, in order to detect in which polygonal region of the macro-triangle a given point p belongs to, we only need to evaluate all the 18 interior lines at p and to collect the resulting signs in a Boolean vector.

Polynomial of degree four

We start our numerical investigation considering as given function a polynomial of degree four. From Section 7.4.5 we expect that the L ∞ norm of error we commit when approximating a polynomial of degree four with our three QIs decreases as O(h 4 T ), with h T the maximum length of the macroelements in the triangulation. In order to demonstrate the optimality of our proposed approximations, we consider five nested tridirectional meshes in [0, 1] 2 with axis-aligned edge lengths from 2 0 to 2 -4 ; this experiment has been developed by using the bivariate quartic polynomial f ∈ P 4 f (x, y) = 1 4 x 4 + x 3 y. (7.29) In Figure 7.8 are presented, respectively, the graph of the polynomial f and its corresponding quasi-interpolant Q ave f , while Figure 7.10 contains the 5 nested meshes we used to test the optimal convergence rate shown in Figure 7.9. As we can notice from the convergence plot, for the polynomial case the errors obtained with the averaging and weighted least-squares are almost the same, while with the exact data the error is smaller.

Franke's function

As first non-polynomial example we consider the well-known Franke's function, which is a standard test in approximation methods. Franke's function presents two Gaussian peaks of different heights and a smaller dip; it is analytically defined as: e -((9x-4) 2 +(9y-7) 2 ) .

(7.30) Imitating Section 7.5.1, we test the error convergence of our proposed constructions over tridirectional meshes in [0, 1] 2 with axis-aligned edge lengths 2 -i , i = 0, . . . , 4 (see Figure 7.10). Figure 7.11 presents the graph of the Franke's function and its corresponding quasi-interpolant Q ave f . In this case, we notice from the error plot in Figure 7.12 that the quasi-interpolant Q wls f returns smaller errors than Q ave f , when the size of the mesh becomes small enough; the construction with exact data Q exa f returns the most accurate representation. The optimal convergence rate is once again achieved. We restrict its definition to the squared domain [-1, 1] 2 . It easy to notice that in the region neighboring the line of equation y = x the sigmoid presents a very steep jump between its maximum and minimum value. Hence, in order to construct precise quasi-interpolants preserving this feature, we will use in this example non-uniform triangulations which are finer, at every step, in the regions of the domain close to the slope. The graph of the sigmoid function and its Q exa f quasi-interpolant are contained in Figure 7.13, while Figure 7.14 shows these adaptative triangulations. Figure 7.15 presents the error plot for this investigation. Similarly to the previous example, the best approximation is achieved with the quasi-interpolant Q exa f , while when dealing with recreated Hermit data the Q wls f quasi-interpolant does a better job. 10 1.4 10 1.5 10 1.6 10 1.7 10 1.8 10 1.9 10 -3

10 -2 10 -1 √ dof L ∞ error Q exa f Q ave f Q wls f

Function with three peaks

In real applications it is common to deal with functions which presents some discontinuities: in these situations, most of the usual numerical tools have to be adapted in order to treat them. Here we show that our construction works perfectly also in this challenging setting without any ah-hoc tuning. Let f be the function defined as From standard computation it is easy to notice that the function f presents three cusps in correspondence of the pairs (-3, -3), (0, 0), (3, 3) for x and y, respectively. Thus, in these points its first (and second) derivative will result discontinuous; as consequence, in this experiment we won't be able to construct the quasi-interpolant Q exa f , while there are no problems for the other two approaches Q ave f and Q wls f . Furthermore, an interesting aspect of our construction is to regularize discontinuous functions with faithful C 2 approximations. Figure 7.16-(a) shows the plot of the function f defined in (7.32), while Figure 7.16-(b) the graph of its QI Q ave f . Due to the presence of the peaks we decided to use adaptative meshes which present an increasing number of triangles nearby the regions of the domain where the irregularities are; Figure 7.17 shows these locally refined triangulations while in Figure 7.18 is illustrated the behavior of the L ∞ error for Q ave f and Q wls f . Once again, we notice that the approximation achieved with Hermite data from local weighted least-square polynomials produces a smaller error.

f

Function on a pentagon domain

In the previous examples we set our approximation problem over triangulations of a squared domain; moreover, this is not a restrictive assumption. In fact, another key point of our construction is its generality over arbitrary triangular meshes. In the following experiment we consider the polygonal domain ⊂ [0, 1] 2 identified by the set of vertices {(0.2, 0), (0, 0.55), (0.5, 1), (1, 0.55), (0.8, 0)} (7.33) on which we define a function f that, for a matter of graphic convenience, we take such that it has zero value along the boundaries of . Let r i (x, y) = 0, i = 1, . . . , 5 be the implicit equations of the five lines identifying the pentagonal domain (7.33); so we have the polynomial q(x, y) = 5 i=1 r i (x, y), which is used to define the objective function f of our investigation f (x, y) = e -((x-0.5) 2 +(y-0.45) 2 ) q(x, y).

Conclusion and perspectives

In this thesis, we have investigated new smooth spline structures suitable for modelling shapes, as well for approximation problems and numerical simulations. We began in Chapter 2 presenting the construction of a G 1 family of surfaces, described by mean of smoothing masks, having the property of approximating the well-known Catmull-Clark subdivision surface. The coefficients of the above masks are explicitly defined through explicit formulas generated from the equations ensuring the geometric continuity between Bézier patches. The solving procedure led to multiple solutions which have been individually analyzed in order to select the one returning the smoothest result, as has been proven in the several proposed experiments.

Then we moved to the investigation of a set of basis function generating the G 1 spline space introduced in Chapter 2. Chapter 3 was devoted to the creation of such functions through a constructive procedure together with a deep analysis of their properties and a dimension formula for the space they generate which only depends on the combinatorial features of the mesh. We used these objects to solve point cloud fitting problems, whose high-quality result is provided by the various numerical experiment we proposed.

Chapter 4 was dedicated to the investigation of analysis-suitable basis functions for IGA simulations. This construction is essentially led by mimicking the building process made for the bases in Chapter 3, but generalizing the previous definitions for more general knot vectors. Also in this case a closed formula for the spline space dimension has been given. The suitability of the proposed bases in numerical IGA simulations of PDEs is shown in some classical tests.

Chapters 5 and 6 presented two interesting applications of the novel constructions in Chapters 2 to 4. The first one concerned the conversion, via point cloud fitting, of CAD models into smooth spline representations. We proposed a complete pipeline which describes the several steps of this translations combined with practical examples. Moreover, this approach turns CAD geometries (that are not appropriate for simulations) to analysis-suitable structures. Equally interesting, the second application focused on the solution of the shallow-water equation when dealing with lakes. Due to its strong dependence to the bathymetry conformation, a crucial point for a realistic study is to gain a precise representation of it: in our case, it is obtained via fitting a cloud describing a real lake seabed with our technique, to which follows the numerical resolution of the equation.

Finally, Chapter 7 was centered on the construction of cubic C 2 quasi-interpolant operators supported on a generic triangulation refined according to the cubic Wang-Shi split. The quasi-interpolation operators are computed by means of a local simplex spline basis which allows a simple derivation of coefficients for their representation in the simplex basis. Their approximation power and quality is confirmed by the examples we proposed which treat several challenging situations.

Several interesting perspectives can be considered. Concerning the construction of the G 1 ACC surface, it might be generalized to patches of arbitrary degree; in fact, the choice to use bidegree 5 Bézier patches, in order to have sufficient degrees of freedom, can be substituted with spline patches of any degree defined on an appropriate knot vector. Also interesting is the generalization to quad meshes presenting connected EVs. In this framework, the systems defining the geometric constraints around each vertex result to be interconnected with each other, making the solution of the equations more difficult. It would also be interesting to extend of the ACC scheme for the construction of G 2 surfaces; actually, starting from the G 1 scheme presented in this thesis, we believe that the G 2 smoothness can be reached imitating the solving strategy used for the tangent plane continuity, but applied to the proper equations governing the desired higher continuity. More challenging, but not impossible, is the extension of the G 1 ACC construction to the volumetric case. It is not a straightforward task, since some theoretical tools are not well defined yet, as for example the definition of gluing data function across faces of a volumetric mesh and the compatibility condition around a vertex.

A similar outlook may be applied to the construction of the basis functions. More precisely, the definition of G 1 functions can be generalized to patches of arbitrary degree and knot vector as well as to topologies with connected EVs. Also the G 2 case is interesting to be studied and the more demanding volumetric generalization.

Also of interest is to develop a pipeline converting CAD objects to IGA suitable models that takes into account the features of the input geometry, since it is a sensitive aspect in practical applications.

Lastly, regarding the quasi-interpolation operators, a similar approach may be used to define quartic C 3 quasi-interpolants by using the simplex spline basis defined in [START_REF] Lyche | A local simplex spline basis for C 3 quartic splines on arbitrary triangulations[END_REF]. It is also interesting to investigate their application to the isogeometric analysis environment, e.g. in the definition of ad-hoc quadrature rules for triangular domains.

Figure 1 . 1 .

 11 Figure 1.1. Quadratic B-splines defined on the uniform knot vector t = {0, 1, 2, 3, 4, 5}.

Figure 1 . 3 .

 13 Figure 1.3. Cubic B-spline curve on t = {0, 0, 0, 0, 1/3, 1/2, 1/2, 2/3, 1, 1, 1, 1} (a) and degree 7 Bézier curve (b) obtained from the same control polygon with coinciding endpoints (red).

Figure 1 . 4 .

 14 Figure 1.4. Single biquadratic B-spline basis function (a) and entire bases set (b) defined on t × s, with t = s = {0, 0, 0, 1/2, 1, 1, 1}.

Figure 1

 1 Figure 1.5. Biquadratic B-spline surface and its control net (red) obtained by using the basis functions in Figure 1.4.

Figure 1 . 6 .

 16 Figure 1.6. Geometric interpretation of barycentric coordinates (a) and subregions of R 2 defined by their signs (b).

Figure 1

 1 Figure 1.7. A cubic triangular Bernstein polynomial (a) belonging to the set of 10 cubic triangular Bernstein polynomials (b).

Figure 1 . 8 .

 18 Figure 1.8. Cubic Bézier surface and its control net (red).

Figure 1

 1 Figure 1.9. Univariate M Ξ as projection of a tetrahedron over a line.

Figure 1 . 10 .

 110 Figure 1.10. Reference systems on two adjacent patches.

Figure 1

 1 Figure 1.11. (a): first neighborhood (green) and second neighborhood (blue) of a vertex (red) with valence N = 5. Masks for a new face point (b), edge point (c) and vertex (d) defining the Catmull-Clark algorithm. α = 1 -7/4N, β = 3/2N, γ = 1/4N .

Figure 1 .

 1 Figure 1.13. Three iterations of the Doo-Sabin scheme starting from the input mesh on the left.

Figure 1 .

 1 Figure 1.14. Three iterations of the Loop scheme starting from the input mesh on the left.

Figure 1 .

 1 Figure 1.15. Three iterations of the Butterfly scheme starting from the input mesh on the left.

Figure 2

 2 Figure 2.1. From (a) to (c): ACC3 masks proposed in [LS08] for face ( b1,1), edge ( b1,0 and b0,1) and vertex ( b1,1) Bézier points, respectively. The coloring refers to Figure 1.11-(a).

Figure 2 . 2 .

 22 Figure 2.2. Control points labeling in a biquintic patch.

Figure 2 . 3 .

 23 Figure 2.3. Bézier points involved in G 1 constraints. (a): disposal of the control points involved in the left-hand side (circles) and right-hand side (squares) of (2.12) involved in the computation of the first derivatives. (b)-(c): disposal of Bézier points regarding the left-hand side and right-hand side of (2.17) and (2.28), respectively, involved in the computation of the second derivatives. (d)-(e): disposal of the control points utilized, respectively, in the left-hand side and right-hand side of (2.29) and (2.31) for the computation of the third and fourth derivatives.

Figure 2 . 4 .

 24 Figure 2.4. Example of a boundary EV with four incident faces. Red square: Bézier point for the border EV b (0,...,3) 1,0 . Red circles: initial (b (0) 1,0 ) and final (b (3) 0,1 ) Bézier points for the G 1 relations around a boundary EV.

  , the ACC 5 masks for the boundary (or frontier) control points b (1,...,κ-1

Figure 2 . 5 .

 25 Figure 2.5. Benchmark meshes in [Pet]: the numbers in the names refers to the EVs valences of the meshes.

Figure 2 . 6 .

 26 Figure 2.6. Mesh in Figure 2.5-(e). From left to right: solving strategies NCS-S, NCS-AS, CS-AS, CS-S. Upper row: isophotes. Middle row: Gauss curvature. Lower row: mean curvature.

Figure 2 . 7 .

 27 Figure 2.7. Mesh in Figure 2.5-(l). From left to right: solving strategies NCS-S, NCS-AS, CS-AS, CS-S. Upper row: isophotes. Middle row: Gauss curvature. Lower row: mean curvature.

Figure 2

 2 Figure 2.8. Alien. (a): multipatch color. (b): solid color. (c): isophotes. (d): zoom on an EV of valence 7.The surface was obtained using the solving strategy CS-S.

Figure 2

 2 Figure 2.9. Hand. (a): multipatch color. (b): solid color. (c): isophotes. (d): zoom on an EV of valence 3.The surface was obtained using the solving strategy CS-S.

Figure 2 Figure 2

 22 Figure 2.10. Rabbit. (a): multipatch color. (b): solid color. (c): isophotes. (d): zoom on an EV of valence 6.The surface was obtained using the solving strategy CS-S.

  (a) bird (b) dinosaur (c) disk (d) dolphin (e) gear (f ) hammer (g) head (h) venus

Figure 2

 2 Figure2.12. Other surfaces obtained by complex meshes using the solving strategy CS-S.

Figure 2

 2 Figure 2.13. Catmull-Clark limit surfaces with error colormap representing the geometry error (a)-(c) and normal error (b)-(d) between the ACC3 surface (left objects) and our G 1 surface (right objects) using the data in [LS08, p. 8].

Figure 2 .

 2 Figure 2.14. Experimental convergence rates of the geometry and normal errors for ACC3 and G 1 CS-S under successive Catmull-Clark subdivisions.

  u the masks returning the derivative directions of our surface obtained from (2.15). Denoting D the matrix with columns T u , T v , ∂ u , ∂ v , we observe that two trailing singular values σ i (D) of D are small, which means that the four vectors are almost coplanar i.e. the tangent planes generated by Span{T u , T v } and Span{∂ u , ∂ v } are approaching each other, also verified by the converge of the unit normal (cf. Figure2.14).

b 2, 1 Figure 3 . 1 .

 131 Figure 3.1. Coefficients for a basis function attached to the value at the vertex (a), attached to the first derivatives at the vertex (b)-(c) and attached to the cross derivatives at the vertex (d) for an EV of valence N = 3.

  Figure 3.1-(b) and (c) show the graph of the resulting functions for an EV with N = 3 while Figure 3.2-(b) and (c) presents their plot when in presence of an EV with valence N = 5.

Figure 3 . 2 .

 32 Figure 3.2. Basis function attached to the value of the vertex (a), value of first derivatives (b)-(c) and value of cross derivatives (d)-(e)-(f)-(g)-(h) for an EV of valence N = 5.

  0 , which are obtained making use of (3.3), (3.4) and (3.7). Previous equations can be solved as in (3.14). Repeating the same construction for all the patches in the ring we come up with the N basis functions attached to the cross derivatives. The result in the case N = 3 is shown in Figure 3.1-(d), while Figure 3.2-(d) to (h) presents the plot of the set of cross derivative basis functions for an EV of valence N = 5.

Figure 3

 3 Figure 3.3. (a): labelling of the control points around an inner regular vertex. (b): values for the Bézier points of the four basis functions.

  e. b 0,0 , b 1,0 , b 2,0 , b 3,0 , b 4,0 and b 5,0 . This assumption transforms (3.3) and (3.4) into b2,1 + b 1,2 = 0, b 3,1 + b 1,3 = 0, (3.17)which define the two basis functions living on an extraordinary edge. The simplest solution satisfying (3.17) is to takeb 2,1 = b 3,1 = 1, b 1,2 = b 1,3 = -1,or vice-versa. Figure3.4 presents their plot.

Figure 3 . 4 .

 34 Figure 3.4. The two edge basis functions defined on an extraordinary edge.

Figure 3 . 5 .

 35 Figure 3.5. (a): labelling of the control points across a regular inner edge. (b): values for the Bézier points of the four basis functions.

  RV j i attached to a regular vertex RV j and let denote by f [s] the corresponding control coefficients of the function f . Since from (3.16) Af [s] = 0, it follows that Ker(A) = Span B RV j i [s] .

  the G 1 (M) space, and b (ℓ) i,j the control points associated with the bivariate Bernstein polynomials B (ℓ) i,j defining the ℓ-th patch Ω

Figure 3

 3 Figure 3.6. (a): final point cloud obtained sampling (3.26) at the last refinement step. (b): resulting surface. (c): multipatch representation of the surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 2.

  Figure 3.6 shows the input point cloud (a) and the results at the last step of subdivision (b)-(c) together with its error color plot (d), while Table 3.2 and Figure 3.7 present the resulting errors, when the number of points increases along with the dimension dim G 1 (M T ) -00 0.122e-00 0.834e-01 0.178e-01 0.790e-02 RMSE 0.491e-01 0.149e-01 0.335e-02 0.353e-03 0.570e-04 Table 3.1. Maximal error L ∞ and RMSE for the surfaces in Example 3.4.1 obtained under 5 Catmull-Clark subdivision steps from a point cloud of 150528 elements. dim G 1 (M T ) 0.225e-00 0.844e-01 0.340e-01 0.385e-02 0.187e-03 0.351e-05 RMSE 0.349e-01 0.140e-01 0.338e-02 0.256e-03 0.723e-05 0.119e-06

Figure 3 . 7 .

 37 Figure 3.7. Experimental behaviors of errors for Example 3.4.1. (a): errors obtained when increasing the number of basis functions but keeping fix the size of the point cloud. (b): errors obtained when increasing both the number of basis functions but and the size of the point cloud at every refinement step.

Figure 3 . 9 .

 39 Figure 3.9. Experimental behavior obtained from the errors for Example 3.4.2 when increasing the smoothing parameter λ.

Figure 3

 3 Figure 3.10. (a): point cloud obtained sampling the trivariate function defined in (3.28) over the unit sphere S 2 . (b): quad mesh used to define the G 1 space. (c): final least squares surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 3.

Figure 3

 3 Figure 3.11. Bird. (a): point cloud obtained sampling an ACC3 surface. (b): quad mesh used to define the G 1 space. (c): final least squares surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 100.

Figure 3

 3 Figure 3.12. Dinosaur. (a): point cloud obtained sampling an ACC3 surface. (b): quad mesh used to define the G 1 space. (c): final least squares surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 100.

Figure 3 Figure 3

 33 Figure 3.13. Hammer. (a): point cloud obtained sampling an ACC3 surface. (b): quad mesh used to define the G 1 space. (c): final least squares surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 100.

Figure 3

 3 Figure 3.15. Rabbit. (a): point cloud obtained sampling an ACC3 surface. (b): quad mesh used to define the G 1 space. (c): final least squares surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 100.

Figure 3

 3 Figure 3.16. Venus. (a): point cloud obtained sampling an ACC3 surface. (b): quad mesh used to define the G 1 space. (c): final least squares surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 100.

  , 3.19 and 3.21 present three examples of point cloud fitting obtained using the pipeline presented in this section. Particular interest goes for the example in Figure

Figure 3 . 17 .

 317 Figure 3.17. Uniform toy. (a): input point cloud. (b): reconstructed quad mesh used to define the G 1 space. (c): final least squares surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 2.

Figure 3 . 18 .

 318 Figure 3.18. Reflection lines on the C 0 (a) and G 1 (c) surfaces obtained from the point cloud in Figure 3.10-(a) and relative zooms, (b) and (d), around an EV.

  Figures 3.18 and 3.20 present the reflection lines for the treated examples.

Figure 3 . 19 .

 319 Figure 3.19. Nonuniform toy. (a): input point cloud. (b): reconstructed quad mesh used to define the G 1 space. (c): final least squares surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 2.

  Comparison between C 0 and G 1 dimensions and errors for the examples in Figure3.10 and 3.17

  .

Figure 3 . 20 .

 320 Figure 3.20. Reflection lines on the C 0 (a) and G 1 (c) surfaces obtained from the point cloud in Fig. 3.17-(a) and relative zooms, (b) and (d), around an EV.

Figure 3

 3 Figure 3.21. Icosahedron. (a): input point cloud. (b): reconstructed quad mesh used to define the G 1 space. (c): final least squares surface. (d): error color plot of the ℓ2 distances. The cloud is contained in a bounding box whose longest length is 2.

Figure 4 . 1 .

 41 Figure 4.1. Convergence plot of the error for the Poisson's problem when using the basis functions introduced in Chapter 3. Except for the regular case N = 4, we obtain a suboptimal rate.

  , specializing (1.10) for two Bézier functions with b (k)

Figure 4

 4 Figure 4.2. (a): control points in a Bézier patch with knots t = {0 6 , 1 6 }. The red dot identifies the EV.(b): control points in a patch with a 4-split given by the knots t = {0 6 , 1/2 5 , 1 6 }. The red dot identifies the EV, the blue square the virtual EVs while the blue dots outline the C 0 split lines.

Figure 4 . 3 .

 43 Figure 4.3. Control points involved in the computation of the new vertex basis function. The point surrounded in red refers to the EV. (a): nonzero control points of a basis function defined on tB. (b): nonzero control points after degree elevation to reach the knot vector t = {0 6 , 1/2 5 , 1 6 }. The light green points refer to the C 1 functions laying on the face. (c): control points defining the new vertex function on t after truncation.

Figure 4 . 4 .

 44 Figure 4.4. Set of truncated vertex basis functions for an EV of valence N = 5 defined on the knot vector t = {0 6 , 1/2 5 , 1 6 }. Their non truncated version is shown in Figure 3.2.

  Figure 4.6-(a)-(b)).

Figure 4 . 6 .

 46 Figure 4.6. Extraordinary vertex (red) and regular vertex (green). (a): extraordinary edge where it is defined a single gluing data function (red line). (b): virtual EV (blue square) which splits in two the extraordinary edge, defining a pair of new gluing functions (orange and green lines). The dashed blue lines identify the regions where we impose C 1 continuity. (c): labelling of the subdivided gluing data and subpatches around a virtual EV.

  4.25) where we already imposed the extra condition ⟨ ḃ1,1 | w⟩ = 0 since rank( C) = 3. Finally, (4.25) returns the solution ḃ1,1 = (-1, -1, 1, 1) T that together with (4.24) defines the other sought function. In Figure 4.7-(a) and (b) are shown these two zero basis functions.

Figure 4

 4 Figure 4.7-(c) and (d) presents these two nonzero edge basis functions.

Figure 4 . 7 .

 47 Figure 4.7. Nonzero edge basis functions (a)-(b) and zero edge basis functions (c)-(d) for an extraordinary edge on the knot vector t = {0 6 , 1/2 5 , 1 6 }.

  Figure 4.8-(a) underlines with different colors the gluing data functions and virtual EVs living on the splitted edge. If now we chain the G 1 equations (4.2) to (4.8), first around a virtual EV and then around the next one, we realize that the two systems share the pairs of control points referring to the highest order derivatives, i.e. related to (4.4) and (4.5) (cf. Figure 4.8).

Figure 4

 4 Figure 4.8. (a): subdivided gluing data (colored lines) and virtual EVs (blue and violet squares) when considering the knot vector t = {0 6 , 1/3 5 , 2/3 5 , 1 6 }. (b): control points involved in the G 1 constraints around two virtual EVs. In orange are highlighted the common control points responsible of the concatenating of the two G 1 systems.

Figure 4 . 9 .

 49 Figure 4.9. Quad meshes presenting EVs of different valences used to define our physical domain Ω involved in the IGA simulations.

  the previous minimization problem (4.38) can be reformulated as M c = g, where c = (c i ) i∈I is the vector containing the unknown coefficients. Example 4.4.1. Here we investigate the L 2 projection problem (4.38) for the function g(x, y) = cos(2πx) sin(2πy) (4.39)

Figure 4 . 10 .

 410 Figure 4.10. Convergence plot of the error for the L 2 projection problem in Example 4.4.1

  Figure 4.11. Nonplanar domains utilized in the L 2 projection experiment

Figure 4 . 12 .

 412 Figure 4.12. Convergence plot of the error for the L 2 projection problem over nonplanar domains in Example 4.4.2

Figure 4 . 13 .

 413 Figure 4.13. Convergence plot of the error for the Poisson's problem in Example 4.4.3

Figure 4 . 14 .

 414 Figure 4.14. Convergence plot of the error for the biharmonic equation in Example 4.4.4

Figure 5 . 1 .

 51 Figure 5.1. Schematic of a cubic B-spline with knot vector {0, 0, 0, 0, 1, 2, 3, 4, . . . }. The slope control point is shown in red, and all other control points in blue. The location of the red point influences the shape of the spline only within the first two knot spans.

Figure 5

 5 Figure 5.2. (a): schematic showing the application of Bézier edge conditions to a subdivision surface. The red lines represent the rows of control points which approximate the MCAD edges as cubic B-splines. The blue dashed lines represent the rows of slope control points. Shaded patches are defined using the tensor product of the boundary B-splines, and the unshaded patches are the usual regular bicubic B-spline patches. (b): 3D control cage with vectors pointing to the implied positions of the slope control points. Each boundary control point stores one vector, whereas the corners store three.

Figure 5 . 3 .

 53 Figure 5.3. Schematic outlining the SubD layering process. (a): MCAD B-rep face with 5 vertices. (b): corner blocks are constructed around the MCAD vertices. (c): corners blocks are connected to form edge blocks. (d): Domain is partitioned into corner (red) and edge (green) blocks which can receive structured mesh, and one interior (grey) region which receives unstructured mesh.

Figure 5 . 4 .

 54 Figure 5.4. SubD layering for the wheel arch of the Car model, showing corner (red), edge (green) andunstructured regions (grey). The 2D medial axis (blue) determines the shape of the blocks so that they do not intersect.

Figure 5

 5 Figure 5.5. (a): initial coarse mesh for a 5-sided face. (b): one level of refinement. (c): two levels of refinement (with mesh smoothing applied). The EVs are highlighted in red. After two levels of refinement the location of the EVs meet the topological requirements for the Bézier edge conditions.

Figure 5

 5 Figure 5.6. A single subdivision surface representing the wing tip of the SNC Dream Chaser model, which is comprised of multiple MCAD B-rep faces.The control cage is shown in green, and the limit surface is shaded in orange. Bézier edge conditions (with slope control vectors shown in grey) are applied such that contiguous regions of the limit surface meet with the same C 1 discontinuity as the original B-rep. However, the limit surface is guaranteed to be exactly C 0 continuous along the join reflecting the MCAD edges (blue).

Figure 5 . 7 .

 57 Figure 5.7. Subdivision surface approximation of the KCS ship hull. The control cage is shown in wireframe, and the limit surface shaded.

Figure 5

 5 Figure 5.8. (a): shaded MCAD B-rep faces of the KCS ship hull model, with edges shown in black wireframe. (b): embedding NURBS patches for each face.

Figure 5 . 9 .

 59 Figure5.9. Heat map of the distance from the limit surface to the target MCAD geometry, for the forward section of the NASA CRM. A maximum fitting error of approximately 0.005% of the aircraft length is achieved after 1000 rounds of iterative control cage adjustment.

  Figure 5.10 shows an example of the typical distribution of the point cloud sampling. The procedures introduced here and in previous Section 5.1 have been carried out by using the software CADfix [CAD].

Figure 5 . 10 .

 510 Figure 5.10. An example of the distribution of points in the sampling of the MCAD geometry, for the tail of the NASA CRM. Each patch of the subdivision surface receives a uniform grid of samples, which are projected onto the embedding surfaces of the MCAD B-rep faces.

Figure 5 Figure 5 . 12 .

 5512 Figure 5.11. Zoom-in of the top edge of the KCS hull model. (a): smoothed edge with consequent oscillations due to the G 1 continuity. (b): sharp edge achieved by imposing C 0 regularity along the sharp edge.

Figure 5

 5 Figure 5.13. Dream Chaser shuttle model. (a): input MCAD model. (b): quad mesh extrapolated from the MCAD geometry. (c): point cloud sampling of the original MCAD model. (d): surface in solid color. (e): surface in multipatch color. (f): error color plot representing the ℓ2 distance between the point cloud and the resulting surface.

Figure 5

 5 Figure 5.14. KCS hull model. (a): quad mesh extrapolated from the MCAD geometry. (b): point cloud sampling of the original MCAD model. (c): surface in solid color. (d): surface in multipatch color. (e): error color plot representing the ℓ2 distance between the point cloud and the resulting surface.

Figure 5

 5 Figure 5.15. NASA CRM. (a): quad mesh extrapolated from the MCAD geometry. (b): point cloud sampling of the original MCAD model. (c): surface in solid color. (d): surface in multipatch color. (e): error color plot representing the ℓ2 distance between the point cloud and the resulting surface.

Figure 5 .

 5 Figure 5.16. Solution for the heat equation on the Car model at the instants from T = 0 to T = 0.4 min.

Figure 5 . 17 .

 517 Figure 5.17. Solution for the heat equation on the Dream Chaser model at the instants from T = 0 to T = 0.4 min.

Figure 5 .

 5 Figure 5.18. Solution for the heat equation on the KCS ship model at the instants from T = 0 to T = 0.4 min.

Figure 5 . 19 .

 519 Figure 5.19. Solution for the heat equation on the NASA CRM at the instants from T = 0 to T = 0.4 min.

Figure 6 . 1 .

 61 Figure 6.1. Exact solution (6.4) for the problem (6.3) when Ω is the unit square. (a): solution for Λ = k = j = 1. (b): solution for Λ = 1, k = 2, j = 3.

Figure 6 . 2 .

 62 Figure 6.2. Rogagua lake. (a)-(b): point cloud representing the lake's seabed. (c): curve defining the lake's shape. (d): quad mesh of the domain identified by the boundary curve (c). (e): ℓ2 fitting error color plot. (f): fitted surface determining the bathymetry function h.

Figure 6 . 3 .

 63 Figure 6.3. Solutions of the shallow-water equation for the Rogagua lake corresponding to the first 15 eigenvalues, ordered in increasing magnitude.

Figure 6 . 4 .

 64 Figure 6.4. Orta lake. (a)-(b): point cloud representing the lake's seabed. (c): curve defining the lake's shape. (d): quad mesh of the domain identified by the boundary curve (c). (e): ℓ2 fitting error color plot. (f): fitted surface determining the bathymetry function h.

Figure 6 . 5 .

 65 Figure 6.5. Solutions of the shallow-water equation for the Orta lake corresponding to the first 15 eigenvalues, ordered in increasing magnitude.

2, 1 Figure 7 . 3 .

 173 Figure 7.3. Labeling of the knots on the boundary of the triangle △.

7. 2 Figure 7

 27 Figure 7.4. From (a) to (g): the simplex spline basis function B1, B4, B10, B16, B19, B22, B28.

Figure 7 Figure 7 . 9 .

 779 Figure 7.8. (a): plot of the polynomial f (x, y) in (7.29). (b): plot of Qavef .

Figure 7

 7 Figure 7.10. From (a) to (e): triangulations of the domain [0, 1] 2 with axis-aligned edge lengths 2 -i , i = 0, . . . , 4 involved in the examples in Sections 7.5.1 and 7.5.2.

Figure 7 Figure 7

 77 Figure 7.11. (a): graph of the Franke's function in (7.30). (b): plot of Qavef .

Figure 7

 7 Figure 7.13. (a): plot of the sigmoid in (7.31). (b): plot of Qavef .

Figure 7

 7 Figure 7.14. From (a) to (e): adaptative triangulations of the domain [-1, 1] 2 utilized for the construction of Qexaf, Qavef and Q wls f in the example in Section 7.5.3.

Figure 7 .

 7 Figure 7.15. Error plot in logarithmic scale of the L ∞ error for the quasi-interpolants of the sigmoid function (7.31) over the adaptative meshes in Figure 7.14.

Figure 7

 7 Figure 7.16. (a): graph of the function in (7.32). (b): plot of Qavef .

Figure 7 Figure 7

 77 Figure 7.17. From (a) to (f): adaptative triangulations of the domain [-1, 1] 2 used in the construction of Qavef and Q wls f in the example in Section 7.5.4.

  Figure 7.19 shows the plot of the function f and its quasi-interpolant Q ave f , respectively, while in Figure 7.20 are reported the consequent triangulations of the domain which have been used to investigate the behavior of the L ∞ error presented in Figure 7.21.

Figure 7

 7 Figure 7.20. From (a) to (d): successive triangulations of the pentagonal domain ⊂ [0, 1] 2 used in the example in Section 7.5.5 for the assembly of Qexaf, Qavef and Q wls f .

Figure 7 Figure 7

 77 Figure 7.19. (a): plot of the function in (7.34). (b): plot of Qavef .

Table 3 .

 3 2. Maximal error L ∞ and RMSE for the surfaces in Example 3.4.1 obtained under 6 Catmull-Clark subdivision steps with increasing number of elements in the point cloud to estimate the rate of convergence. The point clouds have respectively 1500, 6000, 24000, 96000, 384000 and 1536000 points.

Table 3

 3

.3. Maximal L ∞ error and RMSE for the surfaces in Example 3.4.1 for the surfaces in Example 3.4.2 computed making use of 1725 basis functions and progressively bigger smoothing parameter λ. The maximum length of the box containing the model in Figure 3.8 is equal to two.

Table 3 .

 3 4. Spline space dimension, maximal error L ∞ and RMSE for the fitting presented in Example 3.4.3.

The point cloud in Figure

3

.10 is surrounded by a box with longest size 3.

  .6 summarizes the key points of these experiments. Table 3.6. Points in the input cloud, quad mesh and corresponding spline space details, maximal error L ∞ and RMSE for the fitting examples in Figures 3.17, 3.19 and 3.21. All the models are contained in a bounding box whose longest length is 2.

		n P	n F coarse quad mesh dim G 1 (M)	L ∞ error RMSE
	uniform toy	113846	2104	33280	0.190e-01 0.248e-02
	nonuniform toy 30908	1860	29391	0.250e-01 0.277e-02
	icosahedron	208708	7892	125205	0.189e-01 0.142e-02

  Table 4.1. Spline space dimension and L 2 error for each level of subdivision in Example 4.4.1Example 4.4.2. The construction of the basis functions in Section 4.2 holds without further adjustments also for nonplanar domains. Here we present again an example of L 2 projection performed on the two multipatch domains in Figure4.11 containing various EVs of different valences: the objective function we consider in this case is

			01	100	3.021e-02	118	5.315e-02
	1	234	5.518e-03	324	4.355e-04	388	2.188e-03
	2	846	1.705e-04	1156	1.359e-05	1408	1.268e-04
	3	3222	4.055e-06	4356	2.892e-07	5368	3.379e-06
	4	12582	7.010e-08	16900	5.035e-09	20968	5.300e-09
	5	49734	1.195e-09	66564	8.157e-11	82888	8.675e-11
			g(x, y, z) = cos	πx 10	sin	πx 10	.
	Table 4.2 and Figure

  Table 4.2. Spline space dimension and L 2 error for each subdivision level of the two domains in Example 4.4.2

			06	846	5.347e-06
	1	1464	4.769e-08	3222	2.232e-07
	2	5976	1.856e-09	12582	8.449e-09
	3	24264	3.474e-11	49758	9.808e-11

  Table 4.3. Spline space dimension and H 1 error for each level of subdivision in Example 4.4.3

							H 1 error
	0	72	5.757e-01	100	1.283e-01	118	2.035e-01
	1	234	3.234e-02	324	3.485e-03	388	1.091e-02
	2	846	1.772e-03	1156	1.282e-04	1408	9.592e-04
	3	3222	5.826e-05	4356	4.330e-06	5368	4.660e-05
	4	12582	1.830e-06	16900	1.370e-07	20968	2.063e-07
	5	49734	5.806e-08	66564	4.265e-09	82888	5.257e-09

  Table 4.4. Spline space dimension and H 2 error for each level of subdivision in Example 4.

							N ) H 2 error
	0	72	6.184e 01	100	1.428e 01	118	2.081e 01
	1	234	8.675e-00	324	1.091e-00	388	2.219e-00
	2	846	1.011e-00	1156	6.770e-02	1408	3.216e-01
	3	3222	5.705e-02	4356	4.221e-03	5368	3.283e-02
	4	12582	3.244e-03	16900	2.637e-04	20968	4.675e-04
	5	49734	1.911e-04	66564	1.747e-05	82888	2.374e-05

  t step under the Crank-Nicolson method. Figures 5.16 to 5.19 present the results of the IGA simulations. These are obtained by setting, as initial conditions, an heat source at the likely location of the engines within the various vehicles represented. This demonstrates a realistic analysis of the thermal behaviour of such models. Table5.2 presents the running times required to compute each fragment of the pipeline.

		Car	Dream Chaser	KCS	NASA CRM
	n b	166198	168187	166881	116607
	n F	10432	10576	10456	7336
	n P	1262272	1279696	1265179	887656
	L ∞ error	7.965e-03	4.610e-02	4.644e-01	2.946e-03
	RMSE	1.292e-04	1.660e-03	5.890e-03	8.264e-05
	Maximal length	3	9	230	1.70

Table 5 .

 5 1. Fitting errors, spline space and CAD model features for the experiments in Figures 5.12 to 5.15.

  Table 5.2. Detailed elapsed time for each step of the procedure.

		Car	Dream Chaser	KCS	NASA CRM
	Cage generation	1 min 45 sec	1 min 4 sec	55.25 sec	57.88 sec
	MCAD sampling	22 min 28 sec	26 min 24 sec	22 min 23 sec	22 min 53 sec
	Basis computation	11.39 sec	10.98 sec	10.04 sec	6.03 sec
	Fitting	2 min 43 sec	3 min 11 sec	2 min 51 sec	1 min 42 sec
	IGA simulation	15.65 sec/t step	9.39 sec/t step	7.39 sec/t step 1 min 17 sec/t step

Table 6 .

 6 2. Spline space, mesh, point cloud and model features concerning the reconstruction of the bathymetry functions in Figures 6.2 and 6.4.

Table 6

 6

.3. First 15 eigenvalues, and corresponding frequencies, for the shallow-water simulation on Orta lake.

  smooth across the edges of T . Conditions on the local spline coefficients c i j to ensure global C 2 smoothness are discussed in [LMS22, Section 4]. The construction of a stable global basis with local support for S 2 3 (T WS 3 ) is also presented in [LMS22, Section 4].

	7.3 The spline space S 2 3 (T WS 3 )	137
	19
	15 ρ 16 ρ 17 ρ 18 ρ
	10 ρ 11 ρ 12 ρ 13 ρ 14 ρ
	8 ρ 9 ρ
	7 ρ
	6 ρ
	3 ρ 4 ρ 5 ρ
	ρ 1 ρ 2 ρ

Table 7 .

 7 Figure 7.5. Example of directions attached to the edge p1p2 defining the operators ρ19, ρ22, ρ25 and its outgoing normal.

	7.3), (7.4) and (7.5).

1. Values of ρj(Bi) for i, j = 1, . . . , 28, where ρj is defined in (

Acknowledgements

This PhD thesis has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 860843.

Publications

The contributions of the thesis are based on the following works: • [START_REF] Marsala | Analysis-suitable G 1 bases on quadrilateral meshes[END_REF] Analysis-suitable G 1 bases on quadrilateral meshes, with Angelos Mantzaflaris and Bernard Mourrain, in preparation.

• [START_REF] Marsala | From CAD to Representations Suitable for Isogeometric Analysis: a Complete Pipeline[END_REF] From CAD to representations suitable for isogeometric analysis: a complete pipeline, with Angelos Mantzaflaris, Bernard Mourrain, Sam Whyman and Mark Gammon, preprint. https://hal.science/hal-04185850.

• [START_REF] Marsala | Maximally smooth cubic spline quasiinterpolants on arbitrary triangulations[END_REF] Maximally smooth cubic spline quasi-interpolants on arbitrary triangulations, with Carla Manni and Hendrik Speleers, submitted for publication.

Implementations

The results present in this thesis have been obtained by using the following languages/libraries:

• The construction of the G 1 surfaces, as well as the basis functions, were obtained using Julia language [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF]. The developed codes are available in the repository https://gitlab.inria.fr/AlgebraicGeometricModeling/G1Splines.jl.

• The least squares fitting and the IGA simulations have been carried out using the G+Smo library [Man20; Lan15; Jüt14]. The developed codes are available at https://github.com/gismo/gismo.

• The quasi-interpolation operators have been implemented in MATLAB language [START_REF]MATLAB. version 9.9.0[END_REF].

Summary

In this chapter we presented a novel construction of G 1 surfaces defined through masks over quad meshes which approximate the well-known Catmull-Clark limit surface obtained from the same initial geometry. Starting from the C 0 scheme in [START_REF] Loop | Approximating Catmull-Clark subdivision surfaces with bicubic patches[END_REF], by degree elevating the extraordinary patches and using of quadratic gluing data functions we achieved explicit solutions for the weights identifying the different masks, after having adequately fixed the degrees of freedom that occurred during the construction; these masks would return the control points forming the net of a biquintic Bézier patch. Several numerical tests on a benchmark of challenging meshes ensured the quality of the latter construction.

case, we compute the fitting surface increasing the smoothing parameter λ from 10 -3 to 10 -1 with a step of 10 -1 . From Figure 3.8-(b) we notice that the fitting surface presents several wrinkles around the middle peak; by increasing the smoothing factor λ we recover regularity in the output function

Summary

In the present chapter we described an easy procedure to build biquintic G 1 smooth Bézier basis functions over quadrangular meshes to be utilized in point cloud fitting problems. Starting from the equations expressing the geometric continuity constraints obtained making use of quadratic gluing data, we performed an extraction approach, based on the different degrees of freedom present in the system, leading to linear systems returning the control point defining the sought functions. Moreover, we provided an analysis of the latter basis functions as well as a dimension formula for the spline space they generate. We concluded providing numerical experiments establishing their quality in point data fitting.

Summary

This chapter has presented another practical application of the tools introduced in Chapters 2 to 4 to solve the shallow-water equation when dealing with shallow lakes. Given a point cloud reproducing the seabed of a certain lake, we provided a pipeline to obtain from it a quad mesh reproducing the lake's shape; having defined on the latter mesh the G 1 basis functions, we were able, using the same basis set, to both fit the point cloud obtaining a smooth description of the bathymetry function and solve the shallow-water equation using the IGA formulation. Two examples on real lake data have been provided showing the power of the proposed approach.

Chapter 7 Cubic C 2 spline quasi-interpolants on arbitrary triangulations see, e.g., [LMS22, Theorem 2]. The partition △ WS 3 looks (and actually is) quite complicated, not only due to the large number of regions but also because of their shape as it is a mixture of triangles, quadrilaterals, and pentagons. Therefore, dealing in practice with the elements of the space S 2 3 (△ WS 3 ) is prohibitive if we just rely on a direct approach by considering a polynomial representation on each region of the split and by taking into account the interelement smoothness. However, the split has a particularly beautiful feature: it is a cross-cut partition of △, formed by the complete graph of the nine points selected on the boundary of △. Hence, it is natural to consider bivariate cubic simplex splines (see Section 1.2.2) to bypass the complex geometric structure of the split and to obtain an elegant representation of S 2 3 (△ WS 3 ). Following [LMS22, Theorem 3], a basis of S 2 3 (△ WS 3 ) consisting of (scaled) simplex splines, {B 1 , . . . , B 28 }, (7.1) can be constructed. More precisely, each B i is obtained by scaling a cubic simplex spline defined by six (counting multiplicity) knots taken among the boundary points of the WS 3 split; the 28 knot sequences are reported in Figure 7.2. There are seven different types of simplex in the set (7.1): Figure 7.4 shows a representative for each type. The basis functions in (7.1) enjoy several interesting properties:

• the scaling factors ensure that they form a nonnegative partition of unity;

• they inherit recurrence relations and differentiation formulas from the simplex spline structure;

• they admit simple conditions for C 2 joins to neighboring triangles in a given triangulation T ;

• cubic polynomials can be represented through a Marsden-like identity;

7.2 Hermite interpolation in S 2 3 (△ WS 3 ) 133

• they lead to well-conditioned collocation matrices for Lagrange and Hermite interpolation using certain sites;

• a control net can be formed that mimics the shape of the spline function represented in this basis.

Construction of quasi-interpolants

In this section we present the construction of three C 2 cubic quasi-interpolants in S 2 3 (T WS 3 ) for an arbitrary triangulation T . For a given (sufficiently smooth) function f , the common idea of the three approaches is to build a quasi-interpolant Qf by using the local simplex spline basis (7.1) on each macro-triangle △ i of the triangulation refined according to the WS 3 split. More formally,

where {µ i j (f )} n T ,28 i=1,j=1 are linear functionals to be constructed so that the following properties hold:

• global smoothness: Qf is globally C 2 ;

• polynomial reproduction: Qp = p, ∀p ∈ P 3 ;

• locality: µ i j (f ) depends on values (and/or derivatives) of f only on triangles

In the spirit of finite (macro-)elements, an efficient way to obtain linear functionals µ i j enjoying the above properties is to construct Qf | △ i as the solution of the Hermite problem in Corollary 7.2.2 where the Hermite data are directly taken (see Section 7.4.2) or accurately estimated (see Sections 7.4.3 and 7.4.4) from f .

A direct solution of such a Hermite problem, however, would require a triangle-dependent collocation matrix and thus would lead to the computation of n T different collocation matrices in T . To circumvent this issue, we consider the Hermite problem stated in (7.2), whose collocation matrix is the same for all triangles in T . The values f j at the right-hand side have to be properly deduced from those in (7.8) in order to get global C 2 smoothness. This is briefly described in the next section.

Consistent local Hermite data

Suppose the Hermite data (7.8) are given for the triangle △. In this section we outline how to compute the values f j , j = 1, . . . , 28 so that the solution s of the problem stated in Proposition 7.2.1 solves the Hermite interpolation problem stated in Corollary 7.2.2 as well.

Of course, f 28 = h 0 . Moreover, the values f j , j = 1, . . . , 18, can be immediately obtained by standard calculus from the data

Let us now focus on f 19 . The direction q 3 p 3 involved in the definition of the operator ρ 19 (see Figure 7.5) can be written as

for some real values γ, δ. Therefore, from (7.8) we deduce ρ 19 (s) = D q 3 p 3 s(q 3 ) = γD p 1 p 2 s(q 3 ) + δD n 3 s(q 3 ) = γD p 1 p 2 s(q 3 ) + δg 3 . Proposition 7.4.2. The quasi-interpolant in (7.13) satisfies

where the constant C is defined in (7.14).

Proof. Consider a triangle △ i of T . Since the basis functions in (7.1) form a nonnegative partition of unity on △ i , we have

The bounds in (7.12) and (7.14) conclude the proof.

Construction with Hermite data from averaged polynomials

In practical problems it is often unrealistic to assume the availability of data describing a bivariate function up to second order derivatives. Hence, in order to solve the problem stated in Corollary 7.2.2, we need to be able to generate reasonable values up to second derivatives for our spline approximation. We describe a first strategy in this section, making use of local cubic polynomial interpolants defined on each triangle of T . Let us start by defining the Bézier domain points b i ℓ , ℓ = 1, . . . , 10, related to the triangle △ i of T , i ∈ {1, . . . , n T }. For simplicity, we define them in barycentric coordinates:

, b i 10 = (0, 0, 1) . (7.15) See Figure 7.6-(a) for an illustration of these points. Suppose we know the values of a function f at them, i.e., z i ℓ = f (b i ℓ ), ℓ = 1, . . . , 10, and set z i = (z i 1 , . . . , z i 10 ) T . We can then build the local cubic polynomial interpolant I i f on △ i , defined by

where c i = (c i 0,0 , c i 0,1 . . . , c i 3,0 ) T is obtained by solving the linear system

and V i ∈ R 10×10 is the Vandermonde matrix in terms of the monomial basis of P 3 evaluated at the Bézier domain points of △ i . The idea behind this construction is to use these local polynomials to produce reasonable derivative values to be used in our approximation process. Since the Hermite data in (7.3) to (7.5) refer to values attached to the vertices, edges, and triangles of T , we can orient our investigation towards the same features. Let v k be a vertex of T , k ∈ {1, . . . , n V }. Making use of the local interpolants defined over the triangles sharing the vertex v k , we define the averaged vertex polynomial V k f as

where |△| stands for the area of the triangle △. The choice of the weights ω l in the above averaging favors the polynomials defined over triangles with small area. Similarly to the vertex case, for each edge e j of T , j ∈ {1, . . . , n E }, we introduce the averaged edge polynomial E j f as

Furthermore, to define the triangle (or face) polynomial F i f attached to the triangle △ i , i ∈ {1, . . . , n T }, we simply take the local cubic interpolant I i f , i.e.,

The positioning of these vertex, edge, and triangle polynomials is visualized in Figure 7.7-(a).

We have now all the ingredients at hand to build the quasi-interpolant Q ave f . Consider the triangle △ i of T , i ∈ {1, . . . , n T }. With reference to (7.8), from the vertex polynomials V k f in (7.17) we sample the values f i k,α,β ; from the edge polynomials E j f in (7.18) we sample the values g i k and g i k,l ; and finally, the value h i 0 is sampled from the triangle polynomial F i f in (7.19). Given these Hermite data we follow the same construction as in Section 7.4.2.

More precisely, according to the strategy outlined in Section 7.4.1, we first compute the values f i ave = (f i ave,1 , . . . , f i ave,28 ) T , and solve the linear system B µ i ave (f ) = f i ave , with B defined in (7.11). Let we then define the quasi-interpolant Q ave : C 0 (Ω) -→ S 2 3 (T WS 3 ) by

From its construction it follows that the quasi-interpolant Q ave f enjoys the global smoothness, cubic polynomial reproduction, and locality properties stated at the beginning of Section 7.4.

Construction with Hermite data from local weighted least-squares polynomials

In this section we propose a different method to reconstruct the Hermite data required in Corollary 7.2.2 by building local cubic polynomials through weighted least squares. Let us assume, as in Section 7.4.3, that the function values of f are known in certain points, which we select to be a modification of the standard Bézier domain points (7.15). Fix ε > 0. For a given triangle △ i , i ∈ {1, . . . , n T }, we define the modified Bézier domain points b i ℓ , ℓ = 1, . . . , 10, in barycentric coordinates as:

.21) Note that the standard Bézier domain points are a special case with ε = 1. The smaller the value of the parameter ε, the closer the points will be to the vertices of the triangle. Figure 7.6-(b) shows the positioning of these points for ε = 1/2. The idea here is to build a single cubic polynomial defined over the neighborhood of each vertex or edge of the triangulation obtained through weighted least squares, where the weights are defined by measuring the distance of each modified Bézier domain point to the target vertex or edge.

More precisely, let z i ℓ = f ( b i ℓ ), ℓ = 1, . . . , 10, i = 1, . . . , n T , Chapter 7 Cubic C 2 spline quasi-interpolants on arbitrary triangulations

The quasi-interpolant Q wls f , among the two defined from reconstructed Hermite data, shows a slightly smaller approximation error. We reach again the expected rate of convergence.

Summary

In this chapter we provided the construction of three cubic C 2 quasi-interpolants over arbitrary triangulations. The quasi-interpolants are locally spanned by a simplex spline basis defined on cubic Wang-Shi refinement of the polygonal domain, making the whole construction very elegant. Thanks to their locality, we are able to compute the coefficients defining each quasiinterpolant just solving a small linear system whose right-hand side contains proper Hermite data ensuring the global C 2 continuity. If not given as input, we presented two different techniques to reconstruct consistent Hermite informations from local cubic polynomials defined thanks to a set of known values of a given function. Several numerical experiments are provided to state the actual approximation quality of the proposed quasi-interpolants.