Michelangelo Marsala
email: michelangelo.marsala@inria.fr

Angelos Mantzaflaris
email: angelos.mantzaflaris@inria.fr

Bernard Mourrain
email: bernard.mourrain@inria.fr

G 1 -smooth Biquintic Approximation of Catmull-Clark Subdivision Surfaces

Keywords: Bézier patches, gluing data, subdivision surface, extraordinary vertices, geometric continuity

Figure 1: Isophotes on approximate bicubic (ACC 3 ) patches and on our biquintic G 1 construction around a valence 6 vertex.

Introduction

Describing accurately a complex shape in terms of its main geometric features is a major challenge in geometric modeling. It 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.

B-Spline representations make a step toward this objective, by describing smooth shapes in terms of control points, which express the parametric features of such parametrized surfaces. However, for the representation of complex shapes, for instance those appearing in Computer-aided Design and Manufacturing (CAD-CAM), these parametric surfaces are usually trimmed and assembled in patchworks to define the complete geometric object. This results in shape descriptions, which may not be exact, that is, they exhibit leaking patches, stretched elements, etc., which do not adhere to the intended geometry of the objects.

Another approach that has been investigated to address this challenge involves subdivision surfaces. Starting from a coarse mesh, which controls the geometric features of the surface, a subdivision process is applied iteratively to obtain finer and finer meshes, which converge to the limit surface. The Catmull-Clark scheme [START_REF] Catmull | Recursively generated B-spline surfaces on arbitrary topological meshes[END_REF] is a well-known subdivision scheme, which operates on quadrilateral meshes. The limit surface is smooth, but it corresponds to an infinite number of bicubic B-Spline surfaces nearby an Extraordinary Vertex (EV) of the initial mesh. Moreover, the surface exhibits unwanted oscillations around EVs.

In [START_REF] Cashman | NURBS-compatible subdivision surfaces[END_REF] the author makes use of non-uniform B-spline knot insertion rules and spectral analysis to obtain subdivision schemes with bounded curvature around EVs. In [START_REF] Loop | Approximating Catmull-Clark subdivision surfaces with bicubic patches[END_REF], an approximation of Catmull-Clark subdivision surfaces by bicubic B-spline patches was proposed. It can be used to render efficiently such surfaces, but it also provides a compact geometric model description, composed of a collection of Bézier patches that can be directly used in CAD-CAM. The control points of the Bézier patches are obtained from the initial coarse mesh using explicit masks, which depend only on the valence of the vertices. Unfortunately, the bicubic Approximate Catmull-Clark surface (ACC 3 ) is not G 1 -smooth nearby an EV: the tangent planes of patches at a point of a common edge through an EV are distinct (see Fig. 1). For a mathematical definition of geometric continuity or G 1 -smoothness via connecting diffeomorphisms, we refer to [START_REF] Hahn | Geometric continuous patch complexes[END_REF].

The construction of G 1 -smooth B-spline surfaces has been extensively investigated over the last decades. To cite a few works, in [START_REF] Fan | On Smooth Bicubic Surfaces from Quad Meshes[END_REF] and [START_REF] Peters | On the complexity of smooth spline surfaces from quad meshes[END_REF], G 1 bicubic macro-patch surfaces for quad-meshes are studied. The minimal number of splits of these macro-patches is analyzed for linear gluing data between two adjacent patches. In [START_REF] Hahmann | Bicubic G1 Interpolation of Irregular Quad Meshes Using a 4-Split[END_REF], bicubic 4-split macro patches and linear gluing data are studied. In [START_REF] Blidia | G1-smooth splines on quad meshes with 4-split macro-patch elements[END_REF], the space G 1smooth 4-split splines on a quad-mesh with quadratic gluing data is analyzed, providing dimension formulae and basis constructions.

The article [START_REF] Peters | Biquartic C1 spline Surfaces over Irregular Meshes[END_REF] proposed a construction of a general G 1 surfaces with biquartic B-spline defined on the mid-point refined mesh. The control points are computed as explicit combinations of the initial vertices and of face points deduced from blending ratios of the face corners. Differently from our approach, this construction does not ensure the interpolation of the Catmull-Clark limit surface on the mesh's vertices.In [START_REF] Peters | Patching Catmull-Clark meshes[END_REF], a construction of bicubic B-spline with 6 internal knots along each direction is proposed to replace the neighborhood of an EV on a Catmull-Clark subdivision surface by a G 1 -smooth B-spline surface. The solution is given implicitly by describing and solving the linear constraints. Moreover, this construction requires more than one initial applications of Catmull-Clark subdivision on the input mesh; this requirement is not needed in our construction.

In [START_REF] Shi | A practical construction of G1 smooth biquintic B-spline surfaces over arbitrary topology[END_REF], arbitrary G 1 -smooth 4-split biquintic B-spline surfaces with a single interior knot on quadrilateral meshes are constructed by solving the so-called tangent and twist derivatives conditions at the vertices as well as some linear relations on the interior control points along an edge. The control points are computed as the solution of a linear system, in contrast to the present work, where a local and explicit computation for the control points is provided by smoothing masks. Indeed, instead of solving the linear system that expresses the regularity of the functions, we work directly with explicit masks for the Bézier control points, therefore removing the need of a linear solver.In [START_REF] Bonneau | Flexible G1 interpolation of quad meshes[END_REF], biquartic 4-split macro-patches are constructed using quadratic gluing data. In these works, the control points are explicitly described in terms of some free control points and some shape parameters; moreover, the G 1 conditions are obtained by solving the global linear system obtained imposing directly the G 1 constraints between adjacent patches. The remaining degrees of freedom are fixed by imposing the minimization of the thin-plate energy. Differently from that work, we provide a complete and explicit construction for the surface without solving any global system. In [START_REF] Karčiauskas | Improved shape for refinable surfaces with singularly parameterized irregularities[END_REF] a guided subdivision approach is developed to construct a bi-3, bi-4 or bi-5 C 1 Bézier surface. This approach, in which a face split is carried out, is intended to obtain a final surface with regular isophotes across extraordinary vertices; to do that, free parameters are introduced and fixed using minimization of a functional to obtain a better result, whereas we provide explicit mask formulas for the construction of the patches.

The research on smooth splines over unstructured meshes is a revived topic due to its importance for isogeometric analysis discretizations [START_REF] Hughes | Smooth multi-patch discretizations in isogeometric analysis[END_REF]. In [START_REF] Bercovier | Smooth Bezier Surfaces over Arbitrary Quadrilateral Meshes[END_REF], the space on G 1 -smooth Bézier patches of bidegree (n, n) is studied including dimensional analysis and basis constructions. The gluing data are built from bilinear parametrization of quadrilateral faces of the planar mesh. The idea of gluing data inferred from an initial C 0 multipatch B-spline model has also been explored for multipatch B-spline models in [START_REF] Collin | Analysis-suitable G1 multi-patch parametrizations for C1 isogeometric spaces[END_REF][START_REF] Kapl | An isogeometric C1 subspace on unstructured multi-patch planar domains[END_REF][START_REF] Kapl | Construction of analysis-suitable G1 planar multi-patch parameterizations[END_REF], and also in relation to locally refined splines [START_REF] Bracco | Isogeometric analysis with C1 hierarchical functions on planar two-patch geometries[END_REF][START_REF] Buchegger | Adaptively refined multi-patch b-splines with enhanced smoothness[END_REF]. In [START_REF] Scott | Isogeometric boundary element analysis using unstructured T-splines[END_REF] a construction of biquartic G 1 basis functions over T-meshes is presented. Te latter is devoted on basis computation for IGA simulation, rather than surface reconstruction. The basis functions are obtained solving a constrained optimization problem involving 20N + 1 constraints, with N the valence of the extraordinary vertex, which can be computationally expensive.Recent approaches also include the D-patch method [START_REF] Toshniwal | Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations[END_REF] or the approximate C 1 construction in [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF]. For the case of subdivision surfaces, analysis methodologies have been explored eg. in [START_REF] Barendrecht | Splines for engineers: with selected applications in numerical methods and computer graphics[END_REF][START_REF] Barendrecht | Efficient quadrature rules for subdivision surfaces in isogeometric analysis[END_REF][START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF][START_REF] Riffnaller-Schiefer | Isogeometric shell analysis with nurbs compatible subdivision surfaces[END_REF][START_REF] Zhang | Subdivision surfaces with isogeometric analysis adapted refinement weights[END_REF], as well as related issues, notably efficient quadrature [START_REF] Barendrecht | Efficient quadrature rules for subdivision surfaces in isogeometric analysis[END_REF][START_REF] Juettler | On numerical integration in isogeometric subdivision methods for PDEs on surfaces[END_REF].

Contributions. We present a new construction of G 1 surfaces from quad-meshes; for our development we will assume that the input mesh features isolated EVs. The constructed surface, made of biquintic Bézier patches, converges quadratically to the Catmull-Clark subdivision surface of the quad-mesh, interpolating it at extraordinary vertices and reproducing it with C 2 (resp. C 1 ) regularity around regular vertices (resp. vertices connected to an EV). By analyzing the conditions of G 1 -regularity between two adjacent Bézier patches with the use of quadratic gluing functions, we present new explicit masks to compute the control points of the Bézier patches of the G 1 surface. We exploit the symmetry of the masks around (extraordinary) vertices to explicitly solve the linear systems obtained in the various steps of our construction, addressing the degeneracies of subsystems for even valence and introducing suitable degrees of freedom in the computation of the masks. Contrary to the classical subdivision approach, the method generates few number of control points (16 or 36 for each patch), which allows us to reconstruct directly the final surface; for this reason, the entire construction is computationally efficient. Moreover, a quantitative analysis of the G 1 surfaces corroborate the geometric quality of the construction. Outline. After introducing some definitions and notation in Section 2, we recall the work of Loop and Schaefer [START_REF] Loop | Approximating Catmull-Clark subdivision surfaces with bicubic patches[END_REF] in Section 3. By applying twice the degree elevation algorithm to Loop and Schaefer's bicubic surface and making use of gluing data functions we obtain sufficient degrees of freedom to impose G 1 condition introduced around extraordinary patches (Section 4). In Section 5, we analyze the linear system induced by the equations defining the G 1 regularity along edges and provide closed formulas for its inverse that yield masks associated to the control points of the Bézier representation. The system is treated by higher-order Taylor expansion around EVs and explicit solutions are presented in all cases. In Section 6, we analyze the quality of the families of solutions obtained by the above constructions. We present isophote and curvature analysis, as well as numerical results for four variants of our scheme. Furthermore, we compare with the CC limit surface and show that our construction converges quadratically to that limit. Section 7 concludes the present work.

Preliminaries

Before going into the construction of the masks of the G 1 surface, we need to introduce some definitions that we will use throughout the paper. First, we introduce the definition of Bézier patch: Definition 1 (Bézier Patch). A Bézier (quad) Patch is a surface defined as

Q(u, v) = n,m i,j=0 b i,j B i,j n,m (u, v) , (u, v) ∈ [0, 1] 2 ,
where b i,j ∈ R 3 are the Bézier control points and

B i,j n,m (u, v) = B i n (u)B j m (v) are the bivariate Bernstein polynomials of bidegree (n, m), that is, B k d (u) = d k u k (1 -u) d-k . Definition 2 (Mesh).
A mesh M is a collection of vertices, edges and faces that defines the shape of a polyhedral object. We call Bézier face point, a point which corresponds to the interior of a mesh face; a point which corresponds to the interior of an edge is a Bézier edge point and a Bézier vertex point is a point which is shared with more than one edge.

Let p be a vertex of a mesh: the valence N of the vertex is defined as the number of edges (or faces) to which it belongs: if N = 4 we will refer to p as a regular vertex, if N = 4 as an extraordinary vertex (EV). A patch presenting an EV among its vertices is called extraordinary patch. Moreover, we call first neighborhood of p the set of points directly connected to p by an edge and second neighborhood of p the set of points belonging to an extraordinary patch which are opposite to the EV (see Fig. 2).

Given an input mesh M, we will use a smoothing scheme to construct a smooth surface:

Definition 3 (Smoothing Scheme). Given a mesh M we define a smoothing scheme 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. The smoothing matrix contains in its columns the different masks M defining the scheme: we have vertex point masks, edge point masks and face point masks. A mask is a set of coefficients (weights) returning the rules with whom the control points are computed from the mesh vertices: the weights of a smoothing (or Bézier) mask must sum to 1 and they can be, possibly, negative.

Usually, the size of the smoothing matrix S depends solely on the valence of the vertex to whom it is attached and to its first and second neighborhood; 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. Here we will work with quad-meshes and we will provide masks which will return control points for each face of the mesh, therefore defining a collection of Bézier patches. Moreover, 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.

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

The bicubic approximate Catmull-Clark surface 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 (see. Fig. 1); moreover, the vertices generated by ACC 3 interpolate the Catmull-Clark limit surface. To recover regularity around these EVs without affecting the adjacent regular vertices, i.e. without compromising the interpolatory property and C 1 smoothness of the surrounding regular vertices, cubic or even quartic patches do not suffice. Consequently, we will work with biquintic patches, defined by 36 Bézier points each. As a starting point for our construction we use the degree elevated masks ACC 5 , which are obtained after applying twice the degree elevation algorithm 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

(a) N N +5 2 N +5 1 N +5 2 N +5 (b) N N +5 2 N +5 1 2(N +5) 1 N +5 1 N +5 1 2(N +5) (c) 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) (d)
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 the M i,j are the ACC 3 bicubic masks (see Fig. 2-(b-d)). Note that the degree elevation is applied here directly on the considered masks to obtain the new weights.

The masks for Bézier points can be written as vectors; 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 , (1)
where the first coordinate is the weight for the EV, followed by the weights for the vertices of first neighborhood and then the second neighborhood (cf. Fig. 2-(a)). We will use this vector notation for masks throughout the paper.

G 1 Constraints

The approach we use to recover regularity in the neighborhood of EVs is to compute new masks of a scheme returning Bézier control points (of a biquintic patch) s.t. the induced surface is G 1 -continuous. Definition 4 (Multipatch surface). A surface defined on a mesh M is a collection of patches f = (f σ ) σ∈M where σ is a face of the mesh and f σ is called the restriction of f on the face σ.

A surface is said to be G 1 if for every point on it there exists a unique tangent plane; using multipatch surfaces it is easy to understand that, without any constraints, this kind of regularity is not automatically ensured along a common edge shared by two adjacent patches. Therefore, to introduce the G 1 -regularity between two adjacent patches we use gluing data functions: From [START_REF] Hahn | Geometric continuous patch complexes[END_REF], a differentiable function f = (f σ ) σ∈M on a mesh M has to verify the following conditions:

   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 ) , u 1 ∈ [0, 1] , (2)
where f 0 , f 1 are the functions defined on faces σ 0 , σ 1 ∈ M sharing an edge e (see Fig. 3) and [a N,N , b N,N ] are the gluing data along the edge e. Since our construction is developed in the parametric domain, we take e = [0, 1]. The functions satisfying eq. ( 2) for all shared edges belong to the linear space of G 1 B-spline functions on M. In this work, we use quadratic symmetric gluing data [START_REF] Bonneau | Flexible G1 interpolation of quad meshes[END_REF][START_REF] Hahmann | Bicubic G1 Interpolation of Irregular Quad Meshes Using a 4-Split[END_REF][START_REF] Hahn | Geometric continuous patch complexes[END_REF] defined as

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 . (3)
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. Translating eq. ( 2) in terms of Bézier patches and using (3), 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) . (4)
If we call M (k) i,j the smoothing mask returning the control point b (k) i,j belonging to the patch k, solving eq. ( 4) 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 , we obtain after some reductions 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 , (5) 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 , (6) 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 , (7) 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 , (8)
M

(1)

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

(1)

5,1 + M (0) 1,5 = 2M (1) 5,0 , (10) 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 , (11)
with ā0 = 2 -a 0 . Masks verifying eq. ( 5)-( 11) 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 eq. ( 11): if we use quadratic gluing data, the degree of the left-hand side of ( 2) is ≤ 5 while the right-hand side has degree 4 + 2 = 6. Equation ( 11) reflects the constrain 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 eq. ( 9) and ( 10) are always satisfied by ACC 5 masks: 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 eq. ( 10). 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

Computing the Bézier masks

In this section, we solve the system (5)- [START_REF] Catmull | Recursively generated B-spline surfaces on arbitrary topological meshes[END_REF] 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 [START_REF] Barendrecht | Splines for engineers: with selected applications in numerical methods and computer graphics[END_REF]. For Bézier points of order one we will solve eq. ( 5); since the latter involves Bézier points on different patches all around the EV, eq. ( 5) translates into a circulant system to be solved. Eq. ( 6) 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 eq. ( 5), 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, eq. ( 7)-( 8) 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 CC subdivision step to the mesh. 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 5.4.

Bézier masks of order one

We compute first the masks M 1,0 = M (1) 1,0 and M (0) 0,1 using eq. ( 5) as follows. If the valence of the EV is N = 3, eq. ( 5) 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 (5); 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 30) and [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF].

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) (12)
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 eq. ( 5) implies the following system to be solved:

(-a 0 I + C + C N -1 )M 1,0 = (2 -a 0 )M 0,0 , (13)
where I is the (2N + 1) × (2N + 1) identity matrix; Fig. 4-(a) shows the position of the control points involved in this construction.

Proposition 1. Let A = -a 0 I + C + C N -1 be as in eq. [START_REF] Collin | Analysis-suitable G1 multi-patch parametrizations for C1 isogeometric spaces[END_REF]. 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 . 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 2 -a 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) (14)
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 (13) can be therefore solved; the following theorem ensures its solvability.

Theorem 1. The system (13) admits a four-dimensional affine space of solutions.

Proof. Let b = (2 -a 0 )M 0,0 be the right-hand side of (13); we want to show that b ∈ Im(A), i.e. 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) and C as in eq. ( 12): 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 (C -z 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
since z 0 = 1; due to the particular block structure of A, we deduce that corank(A) = 4. Thus the solutions of ( 13) form an affine space of dim 4.

In order to compute a solution of eq. ( 13), we add 4 extra constraints to the coefficient matrix; then, since we arrive to an overdetermined system, the solution can be obtained applying the Least Square (LS) algorithm. 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 .
To obtain a new mask the closest possible to M 1,0 we project 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 of (13); then, the general solution of ( 13) is given by

M 1,0 = M * + 4 i=1 µ i k i , µ i ∈ R . ( 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 ) i.e. M 1,0 -M 1,0 | k i = 0, i = 1, 2, 3, 4, where •|• represents the usual Euclidean scalar product. This leads to the solution

µ = K -1 b , (16)
where

K -1 = diag( K -1 , K -1 ), 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 [START_REF] Hahmann | Bicubic G1 Interpolation of Irregular Quad Meshes Using a 4-Split[END_REF] for N = 3, we obtain again the ACC 5 mask M 1,0 .

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

W T = Span{k 1 , k 2 , k 3 , k 4 }, k T 1 = (1, 0, -1, -1, 0, 1), k T 2 = (-1, -1, 0, 1, 1, 0)
. Hence, from eq. ( 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 eq. ( 16) we obtain the weights for the mask M

Bézier masks of second order

Equation ( 6) can be solved both with respect to M 1,1 or M 2,0 . Let us investigate the two cases.

5.2.1.

Deriving M 1,1 by assigning M 2,0 : circulant system approach In this section we use a circulant system approach to obtain M 1,1 . We start by defining the mask M 2,0 , which is done differently for odd and even valence. Since M (0) 1,1 = C M 1,1 , eq. ( 6) 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 . (17)
The involved Bézier points are highlighted in Fig. 4-(b). Let B = I + C be the matrix in [START_REF] Hughes | Smooth multi-patch discretizations in isogeometric analysis[END_REF]. To solve the system we have to distinguish the cases N odd and N even.

5.2.1.1. Case N odd. In this case, from eq. ( 11), substituting known masks we get

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

Proposition 2. The matrix B = I + C is invertible with B -1 = diag 1 2 , B -1 , B -1 where B = I + C and B -1 = 1 2 N j=1 (-1) j-1 C j-1 . (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 ,
that we have relation [START_REF] Kapl | Construction of analysis-suitable G1 planar multi-patch parameterizations[END_REF].

We are now able to compute explicitly M 1,1 .

Corollary 1. 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, . . . , N , (20)
where the indices i + j -1 are understood mod N when i + j -1 > N ; a similar relation for γ is obtained by substituting the values α , β , δ in (20) .

Proof. The proof is obtained by substituting M 1,1 = B -1 d in (19).

5.2.1.2.

Case N even. If we have an even valence, the matrix B is not invertible and in fact corank(B) = 2, due to corank( B) = 1 since from eq. ( 14) with λ = -1 we have det( B) = 0 and B 1,1 = 1 where B 1,1 is the minor of B obtained removing the 1-st row and column, since its cofactor is an upper triangular matrix with only ones on the diagonal. This implies corank(B) = 2. It means that we must add two extra constraints in order to have a solution for the system; the obtained solution has to be meant in the LS sense. First of all we must check if system (17) admits a solution: denoting w 1 = (0, w T , 0) T , w 2 = (0, 0, w T ) T , where w = (1, -1, 1, . . . , -1) T , and Y = (w 1 w 2 ), it follows that Im(Y T ) = Span{w 1 , w 2 } = Leftker(B); in fact, one can see that Y T B = 0. Then, the system (17) admits a solution if and only if d ∈ Im(B) or, equivalently,

Y T d = 0 . (21)
Relation ( 21) is known as Vertex Enclosure Relation (V.E.R.).

In this section 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 [START_REF] Karčiauskas | Improved shape for refinable surfaces with singularly parameterized irregularities[END_REF] for the masks M 0,0 , M 1,0 and M even 2,0 , we notice that [START_REF] Karčiauskas | Improved shape for refinable surfaces with singularly parameterized irregularities[END_REF] in order to have a solution for system [START_REF] Hughes | Smooth multi-patch discretizations in isogeometric analysis[END_REF]. 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, we want that the mask M 2,0 is 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

Y T M 0,0 = 0, Y T M 1,0 = 0 but Y T M even 2,0 = 0. Thus, we need to compute a new mask M 2,0 = [δ • , δ, δ ] T ensuring
B = B w T 1 w T 2 and f = B M even 2,0 0 0 . ( 22
)
Reordering the rows of B we are able to achieve a diagonal block structure with blocks (2, B, B) which brings to the following

Proposition 3. 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 , (23)
with

B 1 ∈ R N ×N and B 2 ∈ R N ×1 . Proof. By definition 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 . We notice that w w T = N j=1 (-1) j-1 j C j-1 . Thus, we have 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. Let H = N j=1 (-1) j-1 j C j-1 . Using (23) 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 . ( 24
)
Using the notation in [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF], 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 .
Substituting this in [START_REF] Peters | SurfLab Shape Gallery Page[END_REF] we arrive at ( BT B)B 1 = B T . Similarly, using the same identity as above,

( 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 3, we have the following:

Corollary 2. The mask M 2,0 = [δ • , δ, δ ] T ensuring (21) 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, . . . , N , (25)
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 ; a similar relation for δ is obtained substituting in [START_REF] Peters | Biquartic C1 spline Surfaces over Irregular Meshes[END_REF] the values δ even .

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

Hence, using mask M 2,0 we are sure to have a solution for [START_REF] Hughes | Smooth multi-patch discretizations in isogeometric analysis[END_REF]. To compute M 1,1 we add two extra constraints because corank(B) = 2: since we want to obtain a mask M 1,1 the closest possible to M 1,1 we project 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 [START_REF] Hughes | Smooth multi-patch discretizations in isogeometric analysis[END_REF] and denoting

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

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

Corollary 3. The the mask M 1,1 = [γ • , γ, γ ] T in the even case is given in terms of 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 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, . . . , N , (27)
where the indices i + j -1 have to be meant mod N when i + j -1 > N ; a similar relation for γ is obtained substituting in [START_REF] Peters | On the complexity of smooth spline surfaces from quad meshes[END_REF] the values α , β , δ .

Proof. Because M 1,1 = B + d and the block shape of B, using eq. ( 23) and ( 26) we deduce relation [START_REF] Peters | On the complexity of smooth spline surfaces from quad meshes[END_REF].

Example 2. For an EV with valence N = 6, and using the formulas in Corollary 2, the mask M 2,0 satisfying ( 21) is given by

M 2,0 = 1 1320
[592, 294, 130, -8, 18, -8, 130, 79, 8, -1, -1, 8, 79] T .

Thus, using eq. ( 28) in Corollary 3, the cross derivative mask is M In this case, substituting the known masks in (6), 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 ) , (29)
that is well defined since a 0 = 0 for N = 4, which is the case for the EVs we are analyzing. In Fig. 4-(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 eq. ( 11) 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 eq. ( 7), 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 , (30)
where Fig. 4-(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 eq. ( 30)

     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 , (31)
where θ ∈ R \ {-1} has to be fixed. If we take θ = 1, eq. ( 31) will return a symmetric solution; indeed, this choice reflect the fact that we are obtaining the solution by orthogonal projection of 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 [START_REF] Toshniwal | Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations[END_REF]; 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 eq. ( 31), we don't have any solution for the value θ = -1 because it translates into a projection along a direction parallel to the G 1 solution space. For the fourth order Bézier masks, from eq. ( 8) 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 , (32)
which can be solved as above, replacing r in [START_REF] Toshniwal | Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations[END_REF] with r in [START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF]; the positioning of the considered control points is shown in Fig. 4-(e).

0,1 ) Bézier points for the G 1 relations around a boundary EV.

Treatment of boundaries

Meshes are not always closed and, therefore, they can presents boundaries in which EVs can eventually lay. Bézier control points returning G 1 smoothness around border situations can be obtained by applying equations ( 5)-( 11) 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 Fig. 5). Let κ be the numbers of patches sharing the boundary vertex: for the Bézier masks of order one, imposing eq. ( 5) (κ -1) times across the inner extraordinary edges we obtain the following linear system:

(-a 0 I + C + C N -1 ) M 1,0 = h , where h =         (2 -a 0 )F (0) 0,0 -F (0) 1,0 (2 -a 0 )F (1) 0,0 . . . (2 -a 0 )F (κ-2) 0,0 (2 -a 0 )F (κ-1) 0,0 -F (κ-1) 0,1         , (33)
M 1,0 = M (1)
1,0 , . . . , M (κ-1) 1,0 , F (0,...,κ-1) 0,0 , F , and -a 0 I + C + C N -1 = A ∈ R (κ-1)×(κ-1) as in eq. [START_REF] Collin | Analysis-suitable G1 multi-patch parametrizations for C1 isogeometric spaces[END_REF]. System (33) can be solved with a similar strategy like in Section 5.1. Regarding the second order masks, even in this case, from eq. ( 11) and as in eq. ( 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 eq. ( 6) (κ -1) times we get the system:

    1 1 0 . . . 0 1 1 . . . 0 . . . 0 1 1     M 1,1 = g , (34)
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 [START_REF] Zhang | Subdivision surfaces with isogeometric analysis adapted refinement weights[END_REF]; 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 (34) 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 eq. ( 30) and eq. ( 32) along the inner extraordinary edges.

Analysis of the solutions and numerical results

From the constructions presented in Section 5 it turns out that we can solve system ( 5)- [START_REF] Catmull | Recursively generated B-spline surfaces on arbitrary topological meshes[END_REF] in more than one way, 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 [START_REF] Peters | SurfLab Shape Gallery Page[END_REF] (cf. Fig 6) which contains shapes with EVs of different valence and particular geometric structure; from the results obtained by this procedure we will find out the solving strategy, i.e. 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 Fig. 12 and showing numerical results for the G 1 continuity. Figure 6: Test meshes in [START_REF] Peters | SurfLab Shape Gallery Page[END_REF]: the numbers in the names refers to the EVs valences of the meshes.

Comparing the different schemes

Let us call CS (circulant system) the solving strategy in Section 5.2.1, NCS (no circulant system) the strategy in 5.2.2 and S (symmetric) and AS (asymmetric) the solutions for the higher derivatives obtained in 5.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 Fig. 7-8 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 . 1 shows the values for these quantities of each solving strategies carried out on the test meshes [START_REF] Peters | SurfLab Shape Gallery Page[END_REF] in Fig. 6. 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 . 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.

Complex meshes

Using the best strategy CS-S obtained in Section 6.1, we developed a G 1 continuity analysis on the complex meshes in Fig. 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 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 6 and Table 12, 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. 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 correspondent 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 Fig. 9. 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. Interestingly, geom err converges quadratically while norm err converges with a rate of 0.7 (we also observed a convergence rate of 1.5 for the non-normalized normal 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. Fig. 10 summarizes these results.

Conclusion

In this work we presented an explicit construction of a G 1 smoothing scheme providing smooth globally biquintic surfaces. Starting from the work of Loop and Schaefer [START_REF] Loop | Approximating Catmull-Clark subdivision surfaces with bicubic patches[END_REF] and making use of Bézier degree elevation algorithm we impose G 1 constraints in order to obtain a set of equations defining a family of surfaces with the desired smoothness; the solutions of these equations are analyzed case by case with respect to the order of Bézier masks and presented in such a way as to have explicit smoothing masks for any valence of the EVs. Moreover, we exhibit a numerical analysis of the new surfaces obtained with this construction to be able to classify their regularity: all the results presented throughout the paper have been obtained by implementing codes using Julia language for the surface reconstruction and G+Smo library [START_REF] Mantzaflaris | An Overview of Geometry Plus Simulation Modules[END_REF] for obtaining the numerical results. It has also been noticed that the whole construction is computationally efficient, even in presence of complex geometries. Two interesting future research directions are envisaged. First, we would like to advance our constructions by providing explicit basis functions and use them in the isogeometric setting, eg. for analysis of shell structures, cf. [START_REF] Riffnaller-Schiefer | Isogeometric shell analysis with nurbs compatible subdivision surfaces[END_REF]. Second, in the spirit of G 1 -smooth volumetric splines (cf. [START_REF] Birner | Bases and dimensions of C1-smooth isogeometric splines on volumetric two-patch domains[END_REF][START_REF] Birner | The space of C1-smooth isogeometric spline functions on trilinearly parameterized volumetric two-patch domains[END_REF][START_REF] Xie | Interpolatory catmull-clark volumetric subdivision over unstructured hexahedral meshes for modeling and simulation applications[END_REF]), we would like to extend the construction to subdivision solids.

Figure 2 :

 2 Figure 2: (a): first neighborhood (green) and second neighborhood (blue) of a vertex (red) with valence N = 5. (b-d): ACC 3 masks proposed in [22] for face ( b 1,1 ), edge ( b 1,0 and b 0,1 ) and vertex ( b 1,1 ) Bézier points, respectively.

Figure 3 :

 3 Figure 3: Left: local coordinate systems between two patches. Right: control point labeling in a biquintic patch.

  ) that is exactly eq. (9) since M

Figure 4 :

 4 Figure 4: Bézier points involved in G 1 constraints. (a) first derivatives: disposal of the control points involved in the left-hand side (circles) and right-hand side (squares) of (13). (b)-(c) second derivatives: disposal of Bézier points regarding the left-hand side and right-hand side of (17) and (29), respectively. (d)-(e) third and fourth derivatives: disposal of the control points utilized, respectively, in the left-hand side and right-hand side of (30) and[START_REF] Weinmüller | Construction of approximate C1 bases for isogeometric analysis on two-patch domains[END_REF].

Figure 5 :

 5 Figure 5: Example of a boundary EV with four incident faces. Red square: Bézier point for the border EV b (0,...,3) 1,0. Red

  , the ACC 5 masks for the boundary (or frontier) control points b

Figure 7 :

 7 Figure 7: Mesh in Fig. 6-(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 8 :

 8 Figure 8: Mesh in Fig. 6-(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 9 :

 9 Figure 9: Catmull-Clark limit surfaces with error colormap representing the geometry error (a)-(c) and normal error (b)-(d) between the ACC 3 surface (left meshes) and our G 1 surface (right meshes) using the data in [22, p. 8].

Figure 10 :

 10 Figure 10: Experimental convergence rates of the geometry and normal errora for ACC 3 and G 1 CS-S under successive Catmull-Clark subdivisions.

Figure 11 :

 11 Figure 11: From left to right: multipatch color, solid color, isophotes and zoom on EVs of valences, respectively, 7, 3, 6 and 8, of some surfaces obtained by complex meshes using the solving strategy CS-S.

Figure 12 :

 12 Figure 12: Surfaces obtained by complex meshes using the solving strategy CS-S.

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

	1,1 =	1 9900	[5016, 1512, 1512, 258, 144, 144, 258, 604, 185, 43, -4, 43, 185] T .
	5.2.2.

Table 1 :

 1 G 1 smoothness n , Gauss curvature difference c G and mean curvature difference c

			el3	hyp3	linel3	stairs3	el5	house5	beam6	linel6	saddle6
		n	4.4e-15	3.4e-15	8.5e-15	3.3e-15	9.9e-15	7.0e-16	2.2e-15	1.5e-14	3.9e-15
	NCS-S	c G	1.9e-00	3.1e-00	3.4e-00	2.7e-00	1.2e-00	11.3e-00 84.7e-00	8.1e-00	1.3e-00
		c M	0.9e-00	2.5e-00	0.3e-00	3.5e-00	0.7e-00	1.4e-00	2.4e-00	2.8e-00	0.5e-00
		n	4.4e-15	3.4e-15	8.9e-15	3.3e-15	9.9e-15	7.5e-16	2.2e-15	1.5e-14	3.9e-15
	NCS-AS	c G	13.6e-00	3.4e-00	132.8e-00 11.0e-00	7.9e-00	16.5e-00 92.1e-00	75.2e-00	1.6e-00
		c M	4.2e-00	2.7e-00	23.7e-00	7.0e-00	4.1e-00	4.9e-00	7.3e-00	23.7e-00	1.4e-00
		n	4.4e-15	3.4e-15	7.0e-15	3.3e-15	9.9e-15	1.0e-15	2.2e-15	1.5e-14	4.0e-15
	CS-AS	c G	5.0e-00	1.3e-00	33.0e-00	4.0e-00	1.3e-00	3.7e-00 43.3e-00	16.0e-00	4.1e-00
		c M	2.0e-00	1.4e-00	8.2e-00	4.2e-00	0.9e-00	1.5e-00 16.7e-00	5.7e-00	3.5e-00
		n	4.4e-15	3.4e-15	6.9e-15	3.3e-15	9.9e-15	9.5e-16	2.6e-15	1.5e-14	4.4e-15
	CS-S	c G	0.8e-00	1.3e-00	6.8e-00	3.9e-00	0.1e-00	3.4e-00 17.4e-00	0.9e-00	1.5e-00
		c M	0.4e-00	1.4e-00	0.6e-00	4.0e-00	0.1e-00	0.8e-00	0.9e-00	0.3e-00	0.6e-00
			smbeam6	monk7	el8	bas8	kpa8	mitsu9	bas10	updwn10 cube3456
		n	2.9e-15	6.2e-15	3.0e-14	5.6e-15	1.0e-14	1.9e-14	1.8e-14	1.0e-14	3.6e-14
	NCS-S	c G	45.6e-00	8.3e-00	7.5e-00	7.7e-00 56.6e-00	1.6e-00 50.0e-00	70.4e-00	0.7e-00
		c M	4.3e-00	6.4e-00	2.9e-00	5.5e-00	6.5e-00	1.4e-00 11.8e-00	19.1e-00	0.5e-00
		n	2.9e-15	6.2e-15	3.0e-14	5.6e-15	1.0e-14	1.9e-14	1.8e-14	1.0e-14	3.2e-14
	NCS-AS	c G	49.5e-00	38.3e-00 107.5e-00 12.1e-00 63.0e-00 23.4e-00 91.8e-00	106.2e-00	3.3e-00
		c M	5.1e-00	10.6e-00	40.8e-00 13.6e-00 17.6e-00 12.8e-00 27.6e-00	29.7e-00	1.9e-00
		n	2.9e-15	6.1e-15	2.6e-14	5.7e-15	1.1e-14	1.9e-14	1.8e-14	1.0e-14	4.0e-14
	CS-AS	c G	35.1e-00	92.0e-00	43.6e-00 16.1e-00 27.6e-00 14.6e-00 87.9e-00	25.3e-00	1.3e-00
		c M	7.8e-00	39.0e-00	16.6e-00 16.8e-00 10.7e-00	8.3e-00 49.5e-00	87.4e-00	1.9e-00
		n	2.9e-15	6.1e-15	2.6e-14	5.7e-15	1.1e-14	1.9e-14	1.8e-14	1.1e-14	4.5e-14
	CS-S	c G	6.0e-00	91.4e-00	0.9e-00 15.1e-00 27.0e-00	1.0e-00 83.3e-00	24.5e-00	0.7e-00
		c M	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

M obtained from strategies NCS-S, NCS-AS, CS-AS, CS-S for every test mesh in Fig.

6

.

Acknowledgments

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


Ludovic Leleu

Thomas Martzel

Arona Fall

Morgane Sanselme

Vincent Levacher

Sylvain Oudeyer
email: sylvain.oudeyer@univ-rouen.fr

Jean-François Brière
email: jean-francois.briere@insa-rouen.fr

Diastereoselective addition of redox active esters to azomethine imines by electrosynthesis

Thanks to metal-and catalyst-free electrochemical conditions in an undivided cell, a series of readily available redox-active N-(acyloxy)phthalimide esters led to an efficient and highly stereoselective addition (85:15 to 95:5 dr) of putative radical species to chiral (racemic and enantioenriched) C5-substituted azomethine imines to provide an array of 31 polyaminated hydrazine derivatives as a single diastereoisomer.

Radical-based redox processes rapidly expanding in modern organic chemistry, thanks to the power of the renewal technology-driven synthesis by means of visible-light photoredox strategies and electrosynthesis. 1,2 Beside ionicbased reagents approaches, the radical-promoted transformations have markedly enlarged the toolbox of organic chemists, by affording complementary selectivities and an exquisite functional group tolerance, although the transformation of important building blocks remain to be addressed. 3 Azomethine imines (AI) have a versatile reactivity profile both for cycloadditions and nucleophilic additions of organometallic anionic-reagents (Scheme 1A), which opened routes to the construction of biologically important hydrazine derived architectures. 4 Very recently, the radical-based transformation of these dipoles was also developed upon visible-light photoredox catalysis (Scheme 1B). 5,6,[START_REF]For analogous radical-based transformation of N-Aminopyridinium Ylides[END_REF] Xiang, Chen and Yang and colleagues demonstrated the single-electron reduction of Dorn-Otto azomethine imines (AI) by an iridiumderived photocatalyst, followed by a radical-radical crosscoupling reaction of transient difluoroalkyl radicals originated from BrCF2C(O)R partners. 5 On the other hand, the group of Paixão developed radical α-amino alkylation and carbamoylation reactions of similar azomethine imines AI, triggered by the photocatalytic oxidative SET of either dimethylamine or dihydropyridine (DHP) precursors. 6 Then, the obtained N-centered radical, or isomer derived thereof, undergoes a reductive SET by the reduced photocatalyst securing a redox-neutral process. In these pioneering and elegant achievements, a single example of a chiral C5substituted azomethine imine was described but, unfortunately, this led to low or not mentioned diastereoisomeric ratios. Thus, the complete potential of azomethine imines in radical and stereoselective transformations is still in its infancy. In that context, we devised an alternative plan from readily available chiral C5-susbtituted azomethine imines 2 4 based on an unprecedented electrochemical approach (Scheme 1C). The discovery in the 80's of N-(acyloxy)phthalimides (NHPs) 1 as potent redox-active esters (RAEs) has elicited a large array of developments in photoredox chemistry, capitalizing on the wide availability of acid starting materials. 8,[START_REF] Murarka | For recent selected reviews[END_REF] More recently, the alkyl NHP ester derivatives 1 emerged as versatile radical precursors upon a decarboxylative cathodic reductive SET event. 9d,10,11 We reasoned that a sequential cathodic electrolysis, 1,2c advantageously occurring at the same electrode would be suited to our strategy. 12 Accordingly, the first addition of the radical species (stemming from the reduction of RAE 1) to azomethine imine 2 would occur at the cathode-solution interface facilitating thereby the subsequent reductive SET of obtained N-centered radical followed by the final protonation. Based on this Electrochemical-Chemical-Electrochemical-Chemical (ECEC) mechanism, we also assumed that the C5chiral pyrazolidinone moiety could allow the stericdiscrimination of the two diastereotopic faces of 2 during the approach of the radical species. 3.0 F.mol -1 (3h13) at 5 mA 70 14

No current [b] 0

Reaction conditions: carried out with 2a (0.2 mmol), 1a (2 equiv) in DMF (0.1 M) at RT in an undivided cell. Yield of the major diastereoisomer determined on the crude mixture by 1 H NMR with Bn2O as internal standard.

[a] High electrical resistance of the solution was observed without supporting electrolyte leading to a much longer reaction time to deliver 2.6 F.mol -1 .

[b] 20% of 3aa and 65% of 2a estimated by 1 H NMR after 30 minutes of electrolysis and a stirring maintain during 2h17 afterwards.

To our delight, and based on working hypotheses, a mixture of N-Boc glycine-derived NHP 1a (2 equiv) and azomethine imine 2a, upon constant current electrolysis with Glassy Carbon (GC) electrodes, led rapidly (less than 3 hours at room temperature) to the corresponding product 3aa in 80% NMR yield for the major diastereoisomer (determined with an internal standard, Table 1, entry 1). 9b Thanks to 1 H NMR and HPLC/MS analyses, a diastereoisomeric ratio (dr) of 90:10 have been estimated on the crude reaction mixture. Interestingly, among the screening of electrodes (entries 2-4, see SI for the complete investigation), the cheap Stainless Steel (SS) electrodes could be used as well with a slight decreased in yield (72%, entry 4). The reaction performed at best (1) in polar aprotic solvent such as DMF (entries 5-6), (2) in the presence of both 0.5 equivalent of n-Bu4NBF4 as supporting electrolyte (entries 7-8), and (3) with 2.5 equivalents of Hunïg base as a sacrificial oxidable species (entries 9-10). An excess of NHP ester 1a (2 equiv, entry 1) was required in order to improve the yield (1.5 equiv of 1a gave 47% yield, entries 1 versus 11), likely due to a competitive dimerization of the transient α-amino radical species (see ESI). Eventually, the delivery of 2.6 F.mol -1 proved to be the optimal charge at 5 mA to get the maximum yield of 80% (Entry 1 versus entries 12-14).

Scheme 2. Scope and limitation. Reaction conditions: carried out with 2 (0.2 mmol), RAEs 1 (2 equiv) in DMF (0.1 M) at RT in an undivided cell. Isolated yields (%) of the major diastereoisomer after purification by column chromatography; diastereoisomeric ratio determined by 1 H NMR, and HPLC/MS if required, on the crude product. [a] Yield determined by 1 H NMR on the crude mixture with Bn2O as internal standard. [b] Isolated yields of the two major diastereoisomers. [c] Structures 3fa (CCDC2154728), 2g-(Z) (CCDC2154729) were ascertained by X-ray analyses. [d] Starting from 2a-(S) (0.2 mmol, 99% ee) giving 3fa in 57% yield and >98% ee. [e] Product 3ia unstable on silica gel and was purified as a stable derivative after N-Boc formation (3ia', 27% over two steps). [f] The reaction was performed in diluted conditions (0.05 M) for solubility issue.

With these efficacious electrosynthesis conditions in hand, we investigated the scope of this reaction with regard to various NHP esters 1a-i (Scheme 2a). Beside the model glycine derived REA 1a, providing product 3aa in 69% isolated yield of the major diatereoisomer and >90:10 dr, the homologated β-alanine derivative 1b also led to the corresponding addition product 3ba in a 84:16 dr and a modest yield of 18% due to purification issue (55% yield estimated by NMR). The addition of the secondary radical intermediate derived from phenyl alanine derivative 1c occurred smoothly to the dipole 2a in a likely similar diastereoinduction, but the relative configuration of the second stereocenter α to R 2 was not controlled affording 4 diastereoisomers (3ca, 40:14:38:8 dr) in 31% (3ca') and 29% (3ca'') isolated yields respectively of the two major stereoisomers. Pleasingly, disubstituted acyclic (R 1 = R 2 = Me, 1d) and cyclic (cyclopropyl to cyclohexyl, 1e-h) N-Boc glycine derivatives, via a putative tertiary radical species, gave rise to the formation of the original products 3da-3ha in 45-64% yields and an excellent stereoinduction (> 95:5 dr, minor stereoisomer hardly seen by 1 H NMR or HPLC/MS). The reaction with the cyclobutyl precursor was uneventfully performed to 1.2 mmol scale with a slight improvement of the isolated yield (from 64% of 72%) or the faraday efficiency (from 50% of 55%, see SI for details) for 3fa. Most importantly, the reaction starting from readily available enantioenriched azomethine imine 2a (99% ee) 4,13 led to the corresponding addition product 3fa in similar enantiomeric excess (>98% ee). Thereby, no significant racemization event takes place during this radical process. Beside amino acid derivatives, this electrosynthesis reaction could also be applied to various alkyl residues. For instance, isopropyl NHP ester 1i furnished 3ia in 57% NMR yield (> 90:10 dr). However, the obtained product turned out to be unstable during purification on silica gel, but was eventually isolated after a N-Boc protection in moderate 27% yield over two-steps. Although further investigation would be required to explain this phenomenon, deshydrogenated products were identified by MS-analyses of the crude mixture, which could stem from the formation of an azomethine imine after loosing of hydrogen. Furthermore, a crystal of the diastereoisomer 3fa suitable to Xray diffraction analysis revealed an intramolecular hydrogen bonding between the NH of the pyrrolidinone and the N-Boc group which might account for the higher stability of addition products derived from aminoacids. Subsequently, making use of the cyclobutyl-glycine precursor 1f (Scheme 2b), the addition reaction was performed on various azomethine imines either unsubstituted derivatives (3fd, 57%) or flanked on the pyrazolidinone part with a methyl (3fb, 63%) and a bulkier isopropyl moiety (3fc, 60%) with high dr in both cases. A limitation was reached with the more sterically hindered dimethyl azomethine imine affording product 3fe in a low estimated yield (<20%). As long as the iminium moiety of the azomethine imine is concerned (Table 1c), variously ortho, meta and parasubstituted aryl (3ff-3fn, 49-81%), a naphthyl (3fo, 88%) and heterocycles such as thiophene and indole derivatives (3fp-3fq, 68-71%) underwent highly stereoselective C-C bond formation (> 95:5 dr). Pyridine-derived precursor 2r turned out to be poorly soluble in the reaction media, and even in diluted conditions led to 3fr in a modest yield of 29%. Eventually, the azomethine imines derived from linear (2s) or branched aliphatic aldehydes (2t-2u) proved to be compatible with these electrochemical conditions leading to the corresponding products 3fs-3fu in high stereoselectivities (> 95:5 dr) but slightly lower yields (46-49%). Being interested in extrapolating these electrochemical conditions to other substrates, without further optimization, it was also demonstrated that the NHP ester 1f could be added to isoquinoline derived C,N-cyclic azomethine imine 4 to form 5 in 50% yield (Scheme 2d). Eventually, the hydrazone 6b, having an ester pendant, was also transformed into the products 7ib and 7fb in 47% and 67% yields respectively as stable products. Most notably, the analogous hydrazone 6a derived from a benzaldehyde was resistant to these conditions which highlights the complementary reactivity with azomethine imines 2. Mindful of taking advantage of the obtained and original chiral products such as 3fa (Scheme 3), we were pleased to observe that an uneventfully N-N bond reductive cleavage took place in the presence of Raney Nickel to afford the polyaminated derivative 8 in 61% yield (as a single diastereoisomer). Next, owing to the significant of hydrazine derived fused-bicyclic heterocycles as bioactive compounds in the literature, 14 we undertook a cyclization strategy of 3fa. The acetal derived from benzaldehyde in the presence of Bi(OTf)3, furnished the original spirocyclic heterocycle 9 in 51%, due to partial conversion of 57%, with more than 95:5 dr along with an in-situ N-Boc deprotection. We carried out preliminary model reactions to get insight into this electrochemical process. The TEMPO (3 equivalents), used as a radical scavenger, largely inhibited the formation of product 3fa in the optimized conditions. Next, N-Boc glycinederived NHP 1a without azomethine imine yielded the (BocNHCH2)2 dimer likely stemming from the recombination of the formed radical species. Accordingly, based on both cyclic voltammetry (CV) experiments (see ESI for further details) and a background literature on the use of NHP ester as RAE in electrosynthesis, [START_REF] Li | For metal-based cathodic reduction of NHP esters[END_REF][START_REF] Liu | For metal-and catalyst-free cathodic reduction of NHP esters[END_REF] we propose the possible mechanism depicted in Scheme 4. The transformation would start by the SET from the cathode to the NHP ester 1, leading to the radical species 10, thus triggering the subsequent fragmentation/decarboxylation events. The thus obtained amino radical would add diastereoselectively to the less shielded face of the azomethine imine 2, at the opposite to substituent at C5 (as depicted in Scheme 4). For this induction model, the initial Z-configuration of 2g was proven from X-ray analyses. Actually, CV experiments clearly showed that NHP esters 1 (1a, E = -1.25 V; 1f, E = -1.27 V vs SCE in DMF) is more easily reduced than the azomethine imine partner such as 2a (Ered = -1.66 V), ruling out a priori a radical-radical couplingmechanism similar to the one reported by Yang and colleague. 5 Then, the obtained radical 11, or an isomer derived thereof, 6a would undergo a second reductive SET, followed by a protonation, along with a sequential cathodic electrolysis. This cathodic sequence was balanced at the counter electrode by the two electrons oxidation of the Hunïg base (E = +0.88 V), as a sacrificial reaction, while furnishing the proton atoms required to obtain 3. Furthermore, it is believed that i-Pr2EtN protects the obtained product 3 or azomethine imine 2 (2a and 3fa, Eox = +1.56 V) to anodic oxidative events, albeit some decomposition at the end of the electrolysis cannot be ruled out, especially when more than 2.6 F.mol -1 is applied (see Table 1, entries 9-10 and 13).

In conclusion, we demonstrated that a likely sequential cathodic electrolysis process leads to an efficient way to add alkyl acid-derived redox active N-(acyloxy)phthalimide esters to chiral C5-substituted azomethine imines with excellent diastereoselectivities, to provide versatile and original hydrazine heterocycles. This strategy opens new opportunities for radical-based transformation of these dipoles upon electrochemical metal-and catalyst-free conditions, while expanding the scope of the growing cathodic-based transformations. 2c This work has been partially supported by University of Rouen Normandy, INSA Rouen Normandy, the Centre National de la Recherche Scientifique (CNRS), European Regional Development Fund (ERDF), Labex SynOrg (ANR-11-LABX-0029), Carnot Institute I2C, the graduate school for research XL-Chem (ANR-18-EURE-0020 XL CHEM), and by Region Normandy through the "ElectROrg" project. This work is part of the EFLUX project funded by the European Union through the European Regional Development Fund (ERDF) and was also supported by the French National Research Agency (ANR) as part of the ANR-16-CE07-0011-01 project OMaChem. We are grateful to Prof. Hassan Oulyadi for further NMR analyses.

Scheme 1 .

 1 Scheme 1. Reactivity of azomethine imine derivatives and radical species.

Scheme 3 .

 3 Scheme 3. Chemical transformations. [a] Relative configuration proposed from NOE experiment.

Scheme 4 .

 4 Scheme 4. Mechanistic proposal.

Table 1 .

 1 Proof of principle and optimization. a

	Entry	Deviation from the standard conditions	Yield [%] 3aa
	1	None	80
	2	Cgraph(+)/Cgraph(-) as electrodes	27
	3	Mg(+)/GC(-) as electrodes	60
	4	SS(+)/SS(+) as electrodes	72
	5	DMA or MeCN or THF or NMP as solvents	61-68
	6	i-PrOH as solvent	0
	7	No n-Bu4NBF4 [a]	58
	8	n-Bu4NBr instead of n-Bu4NBF4	67
	9	No i-Pr2EtN	16
	10	Use of 1.5 equiv of i-Pr2NEt	62
	11	1.

5 equiv of NHP ester 1a 47 12 2.2 F.mol -1 (2h21) at 5 mA 69 13

Conflicts of interest

There are no conflicts to declare.

Notes and references


Cl) is currently the only in situ cosmogenic nuclide applicable in carbonates, Caand K-rich feldspars and aphyric silicate rocks. Because the production reactions of 36 Cl are more numerous and complex than those of other cosmogenic nuclides (e.g. 10 Be, 3 He), comprehensive and user-friendly calculators are essential for routine application of 36 Cl to Earth surface research questions. However, the existing, commonly used calculators are affected by several drawbacks, such as oversimplification of the muon production model, or the absence of choice regarding production rate parameters. Here, we present the novel 36 Cl exposure age calculation tool in development, "CREp 36 Cl dev", which we have implemented as a new functionality in the existing cosmic-ray-exposure program CREp for 10 Be and 3 He exposure ages. Taking advantage of the most recently developed scaling, muon and lowenergy neutron models, "CREp 36 Cl dev" allows calculating surface exposure ages, providing several choices between scaling and muon models and production parameters.

Introduction

The so far existing CREp online calculator for in situ cosmogenic 10 Be and 3 He exposure age calculations [START_REF] Martin | The CREp program and the ICE-D production rate calibration database: a fully parameterizable and updated online tool to compute cosmic-ray exposure ages[END_REF]; https://crep.otelo.univ-lorraine.fr/) is a versatile tool with advantageous characteristics. It has a user-friendly web interface combined with a high computation speed and a permanent link to the ICE-D production rate database, permitting to account for the most recent production rate calibrations (http://calibration.ice-d.org/). It also allows the choice between global, regional or local production rates, several scaling schemes, atmospheric pressure models, and Virtual Dipole Moment database. Chlorine-36 ( 36 Cl) is currently the only in situ cosmogenic nuclide applicable in carbonate environments, and it is now quasi-routinely analyzed in Ca-and K-rich silicate rocks and minerals. However, given the complexity of the 36 Cl production systematics, computation of 36 Cl exposure ages consists in intricate and dense procedures. Unlike most of the other nuclides, multiple production pathways and target elements have to be considered, with their relative contributions depending on the chemical composition of the sample: spallation of the target elements K, Ca, Ti and Fe; muon capture by Ca and K; capture of low-energy neutrons (derived both from cosmic ray particles and from nucleogenic, i.e. non-cosmogenic reactions in the rock) by target element 35 Cl. The flux of these low-energy neutrons within the rock sensitively depends on the concentration of elements with certain neutron-absorption and moderation characteristics, and the nucleogenic 36 Cl production pathway is in addition controlled by the concentrations of U and Th in the rock. Also, the fraction of stopped negative muons that are captured by Ca and K to produce 36 Cl is dependent on the bulk composition (Fabryka [START_REF] Fabryka-Martin | Production of Radionuclides in the Earth and Their Hydrogeologic Significance, with Emphasis on Chlorine-36 and Iodine-129[END_REF]. Therefore, accurate 36 Cl age calculations require input not only of the 36 Cl concentrations deduced from accelerator mass spectrometry (AMS) measurements and site-and sample-specific data collected during sampling, but also of detailed chemical compositions of the bulk rock and of the analyzed material. Two comprehensive and user-friendly 36 Cl calculators have been published and routinely used over the past years, the excel® spreadsheet for 36 Cl exposure ages and profiles by [START_REF] Schimmelpfennig | Sources of in-situ 36 Cl in basaltic rocks. Implications for calibration of production rates[END_REF] and the CRONUS-Earth online calculator [START_REF] Marrero | Cosmogenic nuclide systematics and the CRONUScalc program[END_REF]. They account for all the complex reaction mechanisms, mainly based on the calculations presented in [START_REF] Phillips | An improved approach to calculating low-energy cosmicray neutron fluxes near the land/atmosphere interface[END_REF] and [START_REF] Gosse | Terrestrial in situ cosmogenic nuclides: theory and application[END_REF]. These calculators have considerably facilitated the application of 36 Cl, leading to a growing number of studies in the fields of geomorphology, paleoseismology and paleoclimatology that are based on the use of this cosmogenic nuclide. However, these existing calculators are affected by several drawbacks, such as the need to calculate scaling factors apart and oversimplification of the muon production model, potentially leading to inaccurate calculations at depth [START_REF] Schimmelpfennig | Sources of in-situ 36 Cl in basaltic rocks. Implications for calibration of production rates[END_REF], and the absence of choice regarding production rate parameters by imposition of a global default production rate [START_REF] Marrero | Cosmogenic nuclide systematics and the CRONUScalc program[END_REF]. To overcome these weaknesses, we are developing the new CREp 36 Cl exposure age calculator, based on the Excel® spreadsheet of [START_REF] Schimmelpfennig | Sources of in-situ 36 Cl in basaltic rocks. Implications for calibration of production rates[END_REF], (see Appendix A and B in [START_REF] Schimmelpfennig | Sources of in-situ 36 Cl in basaltic rocks. Implications for calibration of production rates[END_REF] for equations adapted from [START_REF] Gosse | Terrestrial in situ cosmogenic nuclides: theory and application[END_REF]Phillips, 2001, and[START_REF] Phillips | An improved approach to calculating low-energy cosmicray neutron fluxes near the land/atmosphere interface[END_REF], and we included the following improvements.

-The site-specific, time-dependent scaling is automatically calculated. Users have the same choices as for 10 Be and 3 He calculations, i.e. they can select between the two scaling methods Lm (Lal, 1990;[START_REF] Stone | Isotope production[END_REF][START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10 Be and 26 Al measurements[END_REF] and LSD [START_REF] Lifton | Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes[END_REF] combined with the choices between two atmospheric pressure models (US standard atmosphere; NOAA, 1976;and ERA-40, Uppala et al., 2005) and several Virtual Dipole Moment databases. The scaling method LSD is nuclide-specific andin the case of 36 Cl -target element-specific, and was implemented following the model integrated in the CRONUS-Earth online calculator [START_REF] Marrero | Cosmogenic nuclide systematics and the CRONUScalc program[END_REF]. Note that the original study that first conceptualized the LSD model [START_REF] Lifton | Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes[END_REF] did not account for the scaling of 36 Cl. -Muogenic 36 Cl production is not anymore approximated by an exponential function [START_REF] Gosse | Terrestrial in situ cosmogenic nuclides: theory and application[END_REF], but computed after integration of the muon flux and stopping rates as a function of depth (Heisinger et al., 2002a, b), following the implementation of [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10 Be and 26 Al measurements[END_REF] for the Lm scaling model and that of [START_REF] Lifton | Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes[END_REF] for the LSD scaling model.

In the current development version of the calculator, sea level/high latitude (SLHL) production rates for spallation of Ca, K, Ti and Fe and the parameter for cosmogenic 36 Cl production from low-energy neutrons Pf(0) can either be selected from published values in dropdown menus, or users can manually type in their own values. Note that some of the proposed published production rates were not calibrated with neither of the two scaling methods available in CREp and are therefore not always coherent with the user's choice of scaling. This incoherence is planned to be tackled in a future version of CREp. Note that the systematic errors arising from this currently employed simplification are generally within the uncertainties of the calibrated spallation production rates.

Instructions for use and necessary inputs

The use of "CREp 36 Cl dev" follows a very similar principle to the current CREp online calculator for 10 Be and 3 He exposure age calculations. Users need to go through four tabs to choose scaling-related parameters, load their data and choose productions rates, before obtaining the exposure age results.

Tab 1 Parameters (Fig. 1): After ticking the nuclide of interest, users can choose between:

• two scaling models: (1) the empirical Lal-Stone time-dependent model [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10 Be and 26 Al measurements[END_REF][START_REF] Lal | Cosmic ray labeling of erosion surfaces: in situ nuclide production rates and erosion models[END_REF][START_REF] Stone | Isotope production[END_REF] with the muon implementation of [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10 Be and 26 Al measurements[END_REF], and (2) the Lifton-This technical note is a non-peer reviewed preprint submitted to EarthArXiv Sato-Dunai (LSD) theoretical model [START_REF] Lifton | Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes[END_REF] and the associated muon implementation.

• two atmosphere models: (1) the one based on the ERA-40 data-base [START_REF] Uppala | The ERA-40 re-analysis[END_REF], and

(2) the standard atmosphere (N.O.A.A, 1976). • three geomagnetic databases according to (1) [START_REF] Muscheler | Geomagnetic field intensity during the last 60,000 years based on 10 Be and 36 Cl from the Summit ice cores and 14[END_REF], (2) Lifton et al. (2014), (3) Lifton ( 2016), or import their own geomagnetic database for paleomagnetic corrections. Tab 2 Samples (Fig. 2): By clicking on "Download template", users download an excel file that has to be filled with all sample-specific information before loading it on the CREp web page. All sample input is then displayed for verification. Note that there are 51 columns for data input per sample (see Appendix), which can be verified on the display by scrolling to the right. This technical note is a non-peer reviewed preprint submitted to EarthArXiv Tab 3 Production Rates (Fig. 3): In this "CREp 36 Cl dev" calculator, the production rates for spallation of Ca, K, Fe and Ti and the production parameter for the low-energy neutron capture pathway (Pf(0)) can either be selected from dropdown menus with published values or users may type manually their own values. Note that the proposed published production rates are sometimes not scaled with the same parameters as the ones chosen in "1 Parameters". Tab 4 Results (Fig. 4): The calculated exposure ages are given with full uncertainties, i.e. including the production rate errors as well as analytical errors, and with their analytical errors only (without PR error). The scaling factors for each production reaction is displayed (sp: spallation; th: thermal neutron capture; epi: epithermal neutron capture; mu: muon capture), and a summary of the associated chosen parameters is shown. Ticking "PDF (graphic)" produces a graph with gaussian distributions of all ages.

Clicking on "PDF (tab) download" provides an excel file with these results including the probability values of all ages, and the references of the chosen parameters.

Work in progress

The existing 36 Cl production rate calibration data is currently being implemented in the ICE-D calibration data base with the objective to connect it to CREp 36 Cl, similar to the current functioning of 10 Be and 3 He calculations. This will allow for 1) regularly updates on global elemental production rates, whenever new calibration date will be added, and 2) potentially for the choice of local production rates. However, it must be noted that 36 Cl production rate calibrations are still much scarcer than those of the other nuclides, in particular with regard to mono-elemental production rate material. Further work in progress concerns the implementation of an updated muon production models that are relevant for 10 Be and 3 He calculations at depth (<20 cm) and the addition of a new Virtual Dipole Moment database, which will be applicable to all nuclides in CREp.

Fig. 1 :

 1 Fig. 1: Screen shoot of the Parameter tab of "CREp 36 Cl dev".

Fig. 2 :

 2 Fig. 2: Screen shoot of the Samples tab of "CREp 36 Cl dev".

Fig. 3 :

 3 Fig. 3: Screen shot of the Production Rate tab of "CREp 36 Cl dev".

Fig. 4 :

 4 Fig. 4: Screen shot of the Results Rate tab of "CREp 36 Cl dev".

APPENDIX

Instructions for 36 Cl data input in the sample file:

Fill numbers in all cells. If you do not know a sample-specific value, either use a standard value or zero if applicable.

All yellow columns refer to the composition of the "target" material, i.e. the sample material that was totally dissolved for 36 Cl extraction. Usually this is either an etched mineral separate or rigorously etched whole rock grains.

All green columns refer to the major element concentrations in the "bulk" material, i.e. the chemically untreated and homogenized bulk sample, the composition of which can be significantly different from the target material. If it can be assumed that target and bulk material of the same sample have similar compositions (e.g. most limestones or aphyric volcanics that have only been leached with HNO3), the same sample concentrations can probably be used for both the target and bulk materials.

All blue columns refer to the trace element concentrations in the "bulk" material, i.e. the chemically untreated and homogenized bulk sample, the composition of which can be significantly different from the target material. If it can be assumed that target and bulk material of the same sample have similar compositions (e.g. most limestones or aphyric volcanics that have only been leached with HNO3), the same sample concentrations can probably be used for both the target and bulk materials.

"Ztop" refers to the depth of upper side of the sample "Zbottom" refers to the depth of the lower side of the sample (the sample thickness results from the "Ztop" and "Zbottom" and does not have to be provided separately) "formation age" is the age when the rock was formed (in years), e.g. eruption age for volcanic rocks. This is necessary for the estimation of the nucleogenic 36 Cl contribution.