Nelly Villamizar
email: n.y.villamizar@swansea.ac.uk

Angelos Mantzaflaris
email: angelos.mantzaflaris@inria.fr

Bert Jüttler
email: bert.juettler@jku.at

Completeness characterization of Type-I box splines

Introduction

Box splines are locally supported piecewise polynomial functions defined on uniform grids. They were first introduced by de Boor and DeVore in [START_REF] De Boor | Approximation by smooth multivariate splines[END_REF] and are considered a generalization of the univariate B-spline functions to the multivariate setting. From a geometric point of view, box splines can be seen as density functions of the shadows of higher dimensional boxes and half-boxes [START_REF] Prautzsch | Box splines, Handbook of computer aided geometric design[END_REF]. Box splines can also be introduced as a special case of the so-called simplex splines [START_REF] Dahmen | Recent progress in multivariate splines, Approximation theory[END_REF]. They possess a number of useful properties that make them well-suited for applications. For instance, it has been shown that box splines have small support (a few cells of the underlying grid), they are non-negative, form a partition of unity, and are refinable i.e., the box spline spaces on refined grids are nested [START_REF] Böhm | Subdividing multivariate splines[END_REF][START_REF] De Boor | Box splines revisited: Convergence and acceleration methods for the subdivision and the cascade algorithms[END_REF].

Box splines can be defined from an arbitrary set of directional vectors in R n , but of particular interest for Geometric Design are box splines surfaces which are defined on uniform triangulations of the plane. In this article, we focus on type-I box splines. They are splines defined on three directional meshes that are commonly known as type-I triangulations of R 2 .

From the rich literature on box splines, we mention a few monographs and survey articles which include [START_REF] Chui | Multivariate splines[END_REF][START_REF] Daehlen | Box splines and applications[END_REF][START_REF] De Boor | Box splines[END_REF][START_REF] De Boor | Box splines revisited: Convergence and acceleration methods for the subdivision and the cascade algorithms[END_REF][START_REF] Prautzsch | Box splines, Handbook of computer aided geometric design[END_REF], and a few representative publications on three specific topics. Firstly, a substantial number of results on the approximation power of box splines is described in the literature, e.g. [START_REF] Lyche | Quasi-interpolation projectors for box splines[END_REF][START_REF] Ron | The approximation order of box spline spaces[END_REF]. Secondly, several publications discuss techniques for the efficient manipulation of box spline bases. A general stable evaluation algorithm is devised in [START_REF] Kobbelt | Stable evaluation of box-splines[END_REF]. In [START_REF] Kim | Fast and stable evaluation of box-splines via the BB-form[END_REF] the problem of efficient evaluation of box splines is addressed by making use of the local Bernstein representation of basis functions on each triangle. Also, numerical integration schemes, which are important for applications, based on quasi-interpolation have been considered in [START_REF] Dagnino | Numerical integration based on trivariate C 2 quartic spline quasi-interpolants[END_REF][START_REF] Lamberti | Numerical integration based on bivariate quadratic spline quasiinterpolants on bounded domains[END_REF]. Recent applications of box splines include surface fitting [START_REF] Kang | Hierarchical B-splines on regular triangular partitions[END_REF], and solving linear elasticity problems in isogeometric analysis [START_REF] Giannelli | An immersed-isogeometric model: application to linear elasticity and implementation with THBox-splines[END_REF]. In other areas of mathematics, the theory of box splines has been proved useful to compute the volume of polytopes, and to deal with the integrate of continuous functions over polytopes [START_REF] Xu | Multivariate splines and polytopes[END_REF].

In this article, we are interested in the linear spaces of spline functions generated by the the translates of any fixed type-I box spline. These spline spaces share good approximation properties. For example, the set of translates of any type-I box spline form a partition of unity on R 2 , they are globally, and also locally, linearly independent. These properties of type-I box splines were studied by Dahmen and Micchelli in [START_REF]On the solution of certain systems of partial difference equations and linear dependence of translates of box splines[END_REF] and Jia in [START_REF] Jia | Linear independence of translates of a box spline[END_REF]. In particular, they investigated the linear independence of translates of a box spline in [START_REF]Local linear independence of the translates of a box spline[END_REF] and [START_REF]On the local linear independence of translates of a box spline[END_REF].

The linear independence property implies that the set of translates of any type-I box spline constitutes a basis for the spline space they span. In general this property is not satisfied for box spline functions associated to other uniform partitions, and that makes type-I box splines particularly relevant for applications. For instance, the set of translates of box splines defined from a set of four directional vectors, the so-called type-II box splines, are linearly dependent. For a concise treatment of type-II box splines, further references, and alternative proofs of box spline properties, see [START_REF] Lai | Spline functions on triangulations[END_REF]Chapter 12] and [START_REF] Chui | Multivariate splines[END_REF]Chapter 2].

A second interesting feature of type-I box splines is that, although for any fixed type-I box spline of degree d and order of smoothness r, the box spline translates form a basis, in general these translates do not generate all possible piecewise polynomial functions of degree d and global smoothness r over the three-directional mesh. More precisely, the box spline translates span a proper subspace of the space of C r -continuous spline functions S r d (G) of degree at most d on a three-directional mesh G (cf. Figure 1), for any d > 1. If the domain is taken as an infinite grid G, or as an infinite collection of triangles in G, then both spaces are infinite dimensional. However, if the domain is restricted to a finite collection of triangles Ω, as it is the usual setting in practice, then their finite dimension differ. Explicit dimension formulas in terms of the combinatorics of the domain Ω are well known for both spaces, dim S r d (Ω) can be computed using homological methods [START_REF] Mourrain | Homological techniques for the analysis of the dimensions of triangular splines spaces[END_REF], or Bernstein-Bézier methods as in [START_REF] Chui | Multivariate splines[END_REF]Chapter 2].

In this work, we provide a characterization for the space spanned by box splines translates based on supersmooth conditions across the edges of the underlying partition. We prove that for any fixed degree d and order of global smoothness r, the space of splines satisfying these extra local smoothness conditions is precisely the space spanned by the translates of the corresponding type-I box spline. From this we deduce that the type-I box spline spaces are complete with respect to the local polynomial space induced by the box spline translates. The proof of this result uses the Fourier transform of box splines, as well as the algebraic properties of the Bernstein-Bézier representations of type-I box splines. We generalise the classic definition of spline space in Definition 1 to link the spline polynomial pieces to specific polynomial subspaces of R[x, y], and present the main result in the paper in Section 4.

Furthermore, we apply the completeness characterization of box splines into the construction of hierarchical spline spaces based on local refinements of a type-I triangulation. Hierarchical splines constitute a well-established approach to adaptive refinement in geometric modeling [START_REF] Forsey | Hierarchical B-spline refinement[END_REF] and numerical analysis [START_REF] Mustahsan | Finite element methods with hierarchical WEB-splines[END_REF][START_REF] Schillinger | An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces[END_REF][START_REF] Vuong | A hierarchical approach to adaptive local refinement in isogeometric analysis[END_REF]. Hierarchical tensor-product spline spaces were introduced by Kraft in [START_REF] Kraft | Adaptive und linear unabhängige Multilevel B-Splines und ihre Anwendungen[END_REF] using a selection mechanism for B-splines. The method has been refined leading to spline basis with better approximation properties, such as the partition of unity property, strong stability and full approximation power [17,[START_REF]Strongly stable bases for adaptively refined multilevel spline spaces[END_REF][START_REF] Speleers | Effortless quasi-interpolation in hierarchical spaces[END_REF][START_REF] Zeng | Dimensions of biquadratic and bicubic spline spaces over hierarchical T-meshes[END_REF]. It has also been adapted to Powell-Sabin splines [START_REF] Speleers | Quasi-hierarchical Powell-Sabin Bsplines[END_REF], Zwart-Powell elements and B-spline-type basis functions for cubic splines on regular grids [START_REF] Zore | Adaptively refined multilevel spline spaces from generating systems[END_REF]. In an earlier article, we constructed a hierarchical basis for quartic C 2 -continuous box splines [START_REF] Villamizar | Characterization of bivariate hierarchical quartic box splines on a three-directional grid[END_REF]. Quartic hierarchical box splines spaces have also been studied and used for surface fitting applications by Kang, Chen and Deng in [START_REF] Kang | Hierarchical B-splines on regular triangular partitions[END_REF] and [START_REF]Hierarchical box splines[END_REF]. Truncated hierarchical type-I box splines were considered in [START_REF] Kanduč | Adaptive isogeometric analysis with hierarchical box splines[END_REF] and [START_REF] Giannelli | An immersed-isogeometric model: application to linear elasticity and implementation with THBox-splines[END_REF] in connection to isogeometric analysis applications. Other subdivision schemes has been explored in [START_REF] Gérot | Elementary factorisation of box spline subdivision[END_REF]. A C 1continuous scheme based on cubic half-box splines was presented in [START_REF] Barendrecht | A bivariate C 1 subdivision scheme based on cubic half-box splines[END_REF].

The results we present in this article generalize our previous work [START_REF] Villamizar | Characterization of bivariate hierarchical quartic box splines on a three-directional grid[END_REF] on quartic box splines. Our results apply to type-I box splines of any polynomial degree with no restriction on the symmetry of their support.

The remainder of this paper is organized as follows. In Section 2 we introduce the relevant notation for type-I triangulations, spline functions and the directional derivatives. Section 3 concerns the definition and properties of type-I box spline spaces. We define the space of translates and recall existing results on local and global smoothness of these functions. In Section 4 we prove Lemma 2 which is the main result in the paper, and corresponds to the edge-contact characterization for type-I box splines. In Section 5 we construct the hierarchical type-I meshes and the corresponding hierarchical box spline spaces. This construction follows the approach presented in [START_REF] Mokriš | On the completeness of hierarchical tensorproduct B-splines[END_REF] and [START_REF] Villamizar | Characterization of bivariate hierarchical quartic box splines on a three-directional grid[END_REF]. We conclude the paper with some final remarks in Section 6.

Preliminaries

Throughout this article we assume that G is the uniform type-I triangulation of the real plane R 2 , see Figure 1. This triangulation is obtained by drawing in the north-east diagonals in the bi-infinity grid with grid lines at the integers. This triangulation of the plane is associated to three directional vectors, namely e 1 = (1, 0), e 2 = (0, 1) and e 3 = (1, 1), and therefore is also called a three directional mesh. Each line of G is parallel to one of these vectors and go through the points of the integer grid Z 2 .

The collection T of triangles in G are considered as open sets in, and we denote by E the set of all edges in T , and V = Z 2 the set of vertices. The set of edges is the disjoint union E = E 1 E 2 E 3 , where E i is the set of edges that are parallel to the vector e i . The edges in E i are called edges of type i. The combinatorial closure of a triangle ∈ T , denoted by ˆ , is the Fig. 1 Uniform type-I triangulation (or three directional mesh) of R 2 associated to the directional vectors e1 = (0, 1), e2 = (0, 1) and e3 = (1, 1). We denote this grid as G.

set consisting of the vertices and edges of , and itself. Analogously, ε of an edge ε ∈ E is the set consisting of the edge itself and its two vertices.

A multicell domain M is the triangulation in R 2 induced by a finite set of triangles { 1 , . . . , m } ⊂ T , i.e.

M = m i=1 ˆ i .
This means, that for every triangle i ∈ M , all the vertices and edges of i are considered as elements of M . The subspace of R 2 defined by the (topological) closure i of the triangles i defining the multicell domain M will be denoted by M * , namely

(1) M * = ∈M .
Given a multicell domain M , the diamond of an edge ε is defined as the union of all (at most two) triangles of M which have ε as an edge, that is

3(ε) = ∈M, ε∈ ˆ ˆ .
Similarly, the diamond of a vertex ν is defined by

3(ν) = ∈M, ν∈ ˆ ˆ ,
which is the union of the (at most six) triangles in M such that ν is a vertex of ˆ .

Notice that the diamond (•), of an edge or a vertex, depends on the multicell domain M . From the context, it will be clear the particular domain we are considering in each case.

We denote by R[x, y] the space of bivariate polynomials over the real numbers, and for d 0, P d ⊆ R[x, y] is the set of all bivariate polynomials in x and y of total degree d. In our presentation, the polynomial pieces that define the splines are taken from a finite vector subspace V of R[x, y]. This subspace V is not necessarily the same as P d for any polynomial degree d, it may be a proper linear subspace. In this setting, we define the space of continuous splines P(M, V) on a multicell domain M as follows.

Definition 1 Given a multicell domain M , and a vector subspace V ⊆ R[x, y], we define P(M, V) as the set of piecewise polynomials functions on M i.e.,

P(M, V) = f ∈ C 0 (M * ) : f | ∈ V| for each triangle ∈ M ,
where M * is as defined in Equation (1), and f | denotes the restriction of the function f to the triangle , and V| is the restriction to of the polynomials in V (seen as functions on R 2 ).

In particular, when V = P d , the space P(M, V) coincides with the usual space of C 0 -continuous splines (or piecewise polynomial functions) on M of degree at most d.

For any index, also called regularity vector, s = (s 1 , s 2 , s 3 ) ∈ Z 3 0 we consider the mixed directional derivative operator

D s : R[x, y] → R[x, y] p → (∇ • e 1 ) s1 (∇ • e 2 ) s2 (∇ • e 3 ) s3 (p) .
For a given multicell domain M and a polynomial vector space V, we extend the operator D s to elements f ∈ P(M, V) by applying D s to the restrictions f | , namely

(D s f )| = D s (f | ) , for each triangle ∈ M.
Definition 2 For a given index set I ⊂ Z 3 0 , and a vector space of functions V ⊆ R[x, y], we define the space of functions D I (M, V) on a multicell domain M by

D I (M, V) = f ∈ P(M, V) : D s f ∈ C 0 (M * ) for all s ∈ I ,
where M * is as defined in Equation (1).

Remark 1 Using the notation in Definition 2 we have

D ∅ (M, V) = P(M, V) .
If S r d (M ) denotes the space of globally C r -continuous spline functions on M of degree at most d, then S r d (M ) can be written as D I (M, P d ), where I = s ∈ Z 3 0 : s 1 + s 2 + s 3 r . In this paper, for a given multicell domain M , we shall consider piecewise polynomial functions, or splines, on M with a specific order of smoothness associated to each of the three directions associated to the grid G. Namely, for a regularity vector d = (d 1 , d 2 , d 3 ) ∈ Z 3 0 we define the index sets:

I d 1 = s ∈ Z 3 0 : s 2 + s 3 d 1 , I d 2 = s ∈ Z 3 0 : s 1 + s 3 d 2 , I d 3 = s ∈ Z 3 0 : s 1 + s 2 d 3 ,
and consider the spline space S d (M, V) defined as follows.

Definition 3 For a multicell domain M of the three directional grid G, a vector space V ⊆ R[x, y], and a vector d ∈ Z 3 , the spline space with edge smoothness d on M denoted S d (M, V) is defined as the set of piecewise polynomial functions on M such that the derivatives of order s ∈ D i are continuous across the edges of type i for i = 1, 2, 3. More precisely,

S d (M, V) = f ∈ P(M, V) : f | 3(ε) ∈ C di (3(ε) * ) for every ε ∈ E i ∩ M and i ∈ {1, 2, 3} , where 3(ε) * = ε∈ ˆ , ∈M ˆ , see Equation (1) above.
Later in this paper (see Definition 6 below), we shall introduce a spline space but with smoothness conditions at the vertices of the domain M , the notation in Definition 3 will be particularly convenient for that purpose.

In the following example we illustrate Definition 3 for a specific multicell domain in the grid G and a regularity vector d ∈ Z 3 0 . Example 1 Let M be the multicell domain in Figure 2. It is composed by four triangles denoted 1 , . . . , 4 , we take V = P 2 (the polynomials in R[x, y] of degree at most 2), and the regularity vector d = (0, 1, 0). Then

I d 1 = {(i, 0, 0) : i ∈ Z 0 }, I d 2 = {(0, j, 0), (1, j, 0), (0, j, 1) : j ∈ Z 0 } ,
and

I d 3 = {(0, 0, k) : k ∈ Z 0 } . If we define f ∈ P(M, P 2 ) by f | ˆ i = f i , where (2) f 1 = x 2 ; f 2 = (x -y + 1) 2 ; f 3 = (y -2x)(y -2); f 4 = 2(x -1)(y -x) + 2x -y 2 .
Then f is an element in S d (M, P 2 ). In fact, if we put

g 1,2 = (f 1 , f 2 ) = f | (ε1) , since f 1 -f 2 = (y -1)(2x -y + 1) then D s g 1,2 is a continuous function on (ε 1 ) for s ∈ I d 1 .
Similarly, if we put g i,i+1 = (f i , f i+1 ) it is easy to check that D s g 2,3 and D t g 3,4 are continuous functions for every s ∈ I d 2 and t ∈ I d 3 on (ε 2 ) and (ε 3 ), respectively.

In an analogous way as we defined a spline space associated to smoothness along the edges (Definition 3) of a multicell domain, we will introduce a space of splines with additional smoothness at the vertices of a given multivariate domain M . We prepare this definition by listing the possible vertex-vertex contact configurations ∩ = {ν} between any pair of triangles , ∈ T . First we need the following definitions.

Definition 4 Two triangles

and in the grid G, are said to be edgeconnected if there is a collection of triangles 0 , 1 , . . . , m ∈ T such that = 0 , = m and i-1 ∩ i ∈ E for every i = 1, . . . m. Such a collection of triangles 0 , 1 , . . . , m is called an edge-connected chain between and .

(0, 0)

(1, 1) (2, 2) (1, 2) ε 1 ε 2 ε 3 △ 1 △ 2 △ 3 △ 4 Fig. 2 Multicell domain M = ∪ 4 i=1 ˆ i , with edges ε i = i ∩ i+1 for i = 1, 2, 3.
Definition 5 If , ∈ T are triangles such that ∩ = ∅, we define the smoothness type ST( , ) ⊆ {1, 2, 3} as the set of edge-types that are in the shortest edge-connected chain in G between and . If = we define ST( , ) = ∅.

For any given pair of triangles ,

∈ T with a non-empty intersection, we can identify them with a pair of triangles from A to F in Figure 3, and their smoothness type ST( , ) becomes one of the subsets listed in the table on the left of Figure 3.

ST( ,

) A B C D E F A ∅ {1} {1,2} {1,2,3} {2,3} {3} B {1} ∅ {2} {2,3} {1,2,3} {1,3} C {1,2} {2} ∅ {3} {1,3} {1,2,3} D {1,2,3} {2,3} {3} ∅ {1} {1,2} E {2,3} {1,2,3} {1,3} {1} ∅ {2} F {3} {1,3} {1,2,3} {1,2} {2} ∅ △ A △ B △ C △ D △ E △ F 1 1 2 2 3 3

Fig. 3 Smoothness types of a pair of triangles ,

such that ∩ = ∅. We can identify , with two triangles in the picture on the right. The type ST( , ) is the corresponding index set shown in the table on the left side, which is constructed according to the shortest edge-connected chain between them (see Definition 4).

We now use Definitions 4 and 5 to introduce the space of strongly regular splines associated to a multicell domain M in the three directional grid G. Definition 6 For a multicell domain M ⊂ T , a vector space V ⊆ R[x, y], and a regularity vector d = (d 1 , d 2 , d 3 ), a spline f ∈ P(M, V) is said to be strongly regular if the for any pair of triangles , ∈ M such that ∩ = ∅, for every triple of indexes associated to the smoothness type ST( , ) the derivatives of f are C 0 -smooth. (Notice that the edge-connected chain are composed by triangles in the grid G and are not necessarily in M .) The set of all strongly regular splines on M will be denoted as Ŝd (M, V), it is given by

Ŝd (M, V) = f ∈ P(M, V) : f | U ∈ D I (U, V) for , ∈ M, ∩ = ∅ , U = ∪ , and I = i∈ST( , ) I d i .
The set Ŝd (M, V) is the linear space of splines with edge and vertex smoothness d on the multicell domain M .

Example 2 Let M be the multicell domain in Figure 2, V = P 2 and d = (0, 1, 0) as in Example 1. It is easy to check that the piecewise function f defined in Equation ( 2) is in Ŝd (M, V). For instance, if we take the triangles 1 and 4 , the smoothness type ST( 1 , 4 ) = {1, 2, 3}. Then I = 3 i=1 I d i = {(0, 0, 0)}, and in fact f 1 (1, 1) = f 4 (1, 1). Similarly, taking the triangles 2 and 3 , we get ST( 2 , 4 ) = {2, 3} and I =

3 i=2 I d i = {(0, 0, k) : k = 0, 1}.
The polynomials f 2 and f 4 and also their derivatives ∂(f i )/∂(x -y), for i = 2 and 3, have the same value at [START_REF] Barendrecht | A bivariate C 1 subdivision scheme based on cubic half-box splines[END_REF][START_REF] Barendrecht | A bivariate C 1 subdivision scheme based on cubic half-box splines[END_REF].

In contrast, if g is the function on M defined by g| i = g i with g 1 = 0, g 2 = y -1, g 3 = x 2 -2x + y and g 4 = x 2 -y, then g is also in S d (M, V), but it is not in Ŝd (M, V). In fact, for the triangles 2 and 4 , the derivatives

∂g 2 /∂(x -y) = -1 and ∂g 4 /∂(x -y) = 2x + 1. Then g| U ∈ D I (U, V) for U = ¯ 2 ∪ ¯ 4 .
3 . By Example 2 we also know that the set of strongly regular splines Ŝd (M, V) may be properly contained in the spline space S d (M, V).

In the following we give a sufficient condition for the equality S d (M, V) = Ŝd (M, V). For this we introduce the concept of over concave vertices and kissing triangles. Definition 7 For a multicell domain M in the three-directional grid G we say that a vertex ν ∈ V on the boundary of M is over-concave if star (ν) \ M consists of a single triangle, where star (ν) = ∈T, ν∈ ˆ ˆ .

Since we work exclusively in the three directional grid G then star (ν) consists of 6 triangles, 6 edges and v itself, see Figure 4.

Definition 8 Two triangles

and in M such that ∩ = {v} is a vertex v ∈ V will be called kissing triangles.

In Figure 3, for instance, the triangles {A, E}, {A, D} {A, C} are kissing triangles.

• • •

ν M Fig. [START_REF] Chui | Multivariate splines[END_REF] The vertex ν on the boundary of M is an over-concave vertex as in Definition 7.

Proposition 1 If M be a multicell domain in the three directional grid G, such that it does not have kissing triangles nor over-concave boundary vertices, then S d (M, V) = Ŝd (M, V) i.e., all splines with edge smoothness d on M are strongly regular.

Proof The statements follows directly from Definitions 3, 6 and 7.

Box splines on type-I triangulations

In this section we define box splines on the uniform type-I triangulation G defined in Section 2, Figure 1. This triangulation has vertices at all lattice points (i, j) ∈ Z 2 .

Definition 9

If β is a real-valued function on R 2 , we denote by supp (β) the support β, and it is defined as the set of points x ∈ R 2 such that β(x) = 0.

Recall from Definition 7, that the star of a vertex ν ∈ V , denoted star (ν), is defined by star

(ν) = ∈T, ν∈ ˆ ˆ .
It is the multicell domain composed by all the triangles ∈ T which have ν as one of their vertices, together with the edges and vertices of these triangles.

Definition 10

If n = (n 1 , n 2 , n 3 ) ∈ Z 3 is a triple of integers such that n i 1, the type-I box spline B n associated to n is defined recursively by

B n (x) = 1 0 B n-ei (x -te i )dt ,
for x ∈ R 2 and i ∈ {1, 2, 3} such that n-e i 1 = (1, 1, 1); the function B 1 is the classical Courant hat function with support on the star of the vertex (1, 1) given in Figure 5. More precisely, B 1 is the piecewise linear function on R 2 satisfying B 1 (1, 1) = 1 and B 1 (i, j) = 0 for (i, j) = (0, 0), (0, 1), (1, 2), (2, 1) and (2, 2).

(0, 0)

(1, 1) The coordinates n i of n denote the number of convolutions of B 1 along the directions e i . The support of the box spline B n is the zonotope in R 2 formed by the Minkowski sum of the direction vectors e i taken n i times, for i = 1, 2, 3, respectively (see Figure 5 for an example). Moreover, for every n, the box spline B n is strictly positive for all x in the interior of its support, and zero otherwise [START_REF] Lai | Spline functions on triangulations[END_REF]Theorem 12.2].

(1, 0) (2, 1) (2, 2) (0, 1) (1, 2) (1, 1) (1, 1)
From the general theory of type-I box splines, it follows that the box spline B n is independent of the order in which of the vectors e i appear in the recursive construction of B n in Definition 10. This result follows immediately form the formula for the Fourier transform of a type-I box spline [START_REF] Lai | Spline functions on triangulations[END_REF]Theorem 12.6].

For any n, the box spline B n is a piecewise polynomial function on three directional triangulation G, and each polynomial is of total degree n = |n|-2,

where |n| = n 1 + n 2 + n 3 . It is also well-known that B n is a C r -continuous function on R 2 , where r = min i=1,2,3 |n| -n i -2 [28, Theorem 2.4].
There is a rich literature on type-I box splines, a detailed construction and the proof of structural and smoothness properties can be found for instance in [START_REF] Lai | Spline functions on triangulations[END_REF]Chapter 12] or [START_REF] De Boor | Box splines[END_REF]. In particular, it is known that each convolution along a direction e i increases the continuity of the box spline with respect to differentiation in that direction by one [START_REF] Lai | Spline functions on triangulations[END_REF]Theorem 12.3]. Thus, following the notation introduced in Section 2, d = (n

2 +n 3 -2, n 1 +n 3 -2, n 1 +n 2 -2) then (3) B n ∈ S d (G, P n ) .
We denote by B n (G) the set of integer translates of the box spline B n , which is defined as the set

(4) B n (G) = B n ( • -v) : v ∈ Z 2 .
The set B n (G) is also called the set of shifted box splines associated to the direction vector n.

Remark 3 Notice that translates in B n (G) have distinct support, it is the zonotope which is the support of B n shifted by v. In fact, the set B n (G) is (globally) linearly independent [START_REF] Lai | Spline functions on triangulations[END_REF]Theorem 12.19] i.e., if

v∈Z 2 a v B n (x -v) = 0, for all x ∈ R 2
then a v = 0 for all v ∈ Z 2 . Furthermore, it has been shown that the translates in B n (G) are also locally linearly independent i.e., if A is an open set, then the shifted box splines

B n (• -v) : supp B n (• -v) ∩ A = ∅
are linearly independent [START_REF]Line average algorithm: A method for the computer generation of smooth surfaces[END_REF][START_REF]Local linear independence of the translates of a box spline[END_REF]. Here, supp (β) denotes the support of the function β (Definition 9).

We now introduce the definition of the set of active box splines on a given multicell domain M in the type-I triangulation G.

Definition 11 If M ⊆ G is a multicell domain, we define Supp n (M ) = v ∈ Z 2 : supp B n (• -v) ∩ M * = ∅ ,
and the set of active box splines B n on M by

B n (M ) = B n (• -v) M : v ∈ Supp n (M ) ,
where M * is the closure of M in R 2 as defined in Equation (1).

In particular, if the multicell domain M = ˆ , for a triangle ∈ T , then B n ( ˆ ) is the set of translates B n ( • -v) ∈ B n (G) whose support contains . The number of elements in B n ( ˆ ) is given by ( 5)

φ(n) = n 1 n 2 + n 1 n 3 + n 2 n 3 ,
where n = (n 1 , n 2 , n 3 ).

Definition 12 If B n (• -v)
is the translate of the type-I box spline B n by v ∈ Z 2 , then we take (1, 1) -v as the point of reference of supp B 2 (• -v) in the lattice. For a triangle ∈ T , we define the 1-ring neighbourhood of as the set of reference lattice points (1, 1)

-v such that sup B n (• -v) ∩ = ∅.
For instance, if 2 = (2, 2, 2) then the elements in B 2 ( ˆ ) are the translates B 2 (•-v) associated to the φ(n) = 12 lattice points in a 1-ring neighbourhood of in the grid G, which are shown in Figure 6. Notice that to any type-I box spline B n and a triangle if ∈ T we can associate a linear space of polynomials. Namely, extending by linearity, we can take V n | as the space generated by the restriction of the active box splines B n to the triangle i.e.,

V n | = span B n ( ˆ ) .

Fig. 6 The 12 lattice points correspond to the 1-ring of a triangle ∈ T associated to the box spline B2 = B (2,2,2) . Each of these points is the reference point (1, 1) -v for the translate

B2(• -v) such that ⊆ supp (B2(• -v)).
By Equations ( 6) and ( 3), we see that

V n | is a linear subspace of P n ⊆ R[x, y], for any ∈ T and n = n 1 + n 2 + n 3 -2.
We now prove that for any n ∈ Z 3 , the subspace V n | ⊆ P n is independent of the choice of the triangle ∈ T . The proof of this result is a generalization of [START_REF] Villamizar | Characterization of bivariate hierarchical quartic box splines on a three-directional grid[END_REF]Proposition 27], where we consider the case 2 = (2, 2, 2) and the box spline B 2 .

Proposition 2 Let n = (n 1 , n 2 , n 3 ) ∈ Z 3 1 and f ∈ P n , with n = n 1 + n 2 + n 3 -2. Then, V n | = V n |

for any pair of triangles and in T .

Proof Let and be two triangles in the three-directional grid G. Denote by G 0 the triangulation of R 2 obtained by the lines parallel to the vectors e 1 = (a, 0), e 2 = (0, a) and e 1 = (a, a), for a fixed number a ∈ Z. Notice that for any a ∈ Z, we can see the grid G as a refinement of a grid G 0 . Denote by a the triangle in G 0 with vertices at (0, 0), (0, a) and (a, a). In particular, let us take ∈ Z + and a = 2 in such a way that the translate ˜ = a -v of a by a vector v ∈ Z 2 contains both and . Then, B n (1/2 •) is the correspondent box spline associated to n in the grid G 0 . By the refinement equation for box splines [28, Theorem 12.9], there exists a finite sequence

{c v } v∈Z 2 such that (7) B n (1/2 •) = v∈Z 2 c v B n (• -v) .
Let us denote by V 0 n | ˜ the span of the box spline translates B n (1/2 • -v) restricted to ˜ in G 0 . By the symmetry of the box splines supports, the number φ(n) of active translates on a triangle is independent of the grid and of the given triangle in the grid. Furthermore, these translates are linearly independent (see Remark 3). Thus, [START_REF]Line average algorithm: A method for the computer generation of smooth surfaces[END_REF] dim

V 0 n | ˜ = dim V n | = dim V n | .
Taking the restrictions to and , Equation ( 7) implies

V 0 n | ˜ | ⊆ V n | and V 0 n | ˜ | ⊆ V n | . Since V 0 n | ˜ is a polynomial subspace of P n , we have V 0 n | ˜ = (V 0 n | ˜ )| = (V 0 n | ˜ )| . In particular, both (V 0 n | ˜ )| and (V 0 n | ˜ )| have the same di- mension as V 0 n | ˜ .
The statement follows by applying Equation ( 8).

Proposition 2 implies that for a given n ∈ Z 3 1 the restrictions of the translates of the box spline B n to a triangle ∈ T define a polynomial space which is independent of . From now on, we will denote such polynomial space as V n .

Remark 4 For any

∈ T an n ∈ Z 3 1 , Remark 3 implies that B n ( ˆ ) ⊆ P n is a linearly independent set. Here P n is as before, the set of bivariate polynomials of degree at most n = n 1 + n 2 + n 3 -2, and n = (n 1 , n 2 , n 3 ). The number of elements in B n ( ˆ ) = φ(n) (Equation ( 5)) is the dimension of the space of polynomials dim V n associated to the type-I box B n . Thus, dim(

V n ) = n 1 n 2 + n 1 n 3 + n 2 n 3 n1+n2+n3 2
= dim P n . For n i 1, equality holds only for n = 1 = (1, 1, 1), which corresponds to polynomial space V 1 associated to the Courant hat function B 1 in Definition 10. For any other n ∈ Z 3 and the corresponding box spline B n , the polynomial space V n is a proper subspace of P n .

Characterization of box spline spaces

Let ∈ T , n ∈ Z 3 1 , and consider V n = span B n ( ˆ ). As observed in Remark 4, the set B n ( ˆ ) of the box spline translates with support on is linearly independent, and hence the restriction of a polynomial f ∈ V n to has a unique representation [START_REF]On the local linear independence of translates of a box spline[END_REF] f

| (x) = β∈Bn( ˆ ) λ β (f | )β(x), x ∈ ¯ , for coefficients λ β (f | ) ∈ R. Lemma 1 Let n = (n 1 , n 2 , n 3 ) ∈ Z 3 1 , f and f in V n , and M = ˆ ∪ ˆ , such that , ∈ T and ¯ ∩ ¯ ∈ E i for some i ∈ {1, 2, 3}. If f | , f | ∈ S d (M, V n ), then there exist polynomials g, h ∈ V n such that (10) f | , f | = g| , f | + h| , 0| , with λ β g| = λ β f | for every β ∈ B n ( ˆ ) ∩ B n ( ˆ ). Proof Define g| = β∈B( ˆ )∩B( ˆ ) λ β (f | )β| + β∈B( ˆ )\B( ˆ ) 0 • β| .
Then g| ∈ V n | , and extending by linearity we can see g| as the restriction to of a polynomial g in V n ⊆ R[x, y]. By construction, the pair (g| , f | ) satisfies the required condition, and taking h = f -g we obtain Equation [START_REF]On the solution of certain systems of partial difference equations and linear dependence of translates of box splines[END_REF].

Lemma 2 Let n = (n 1 , n 2 , n 3 ) ∈ Z 3 1 , f and f in V n , and , ∈ T , such that ¯ ∩ ¯ ∈ E i for some i ∈ {1, 2, 3}.
Then, for M = ˆ ∪ ˆ the following two statements are equivalent:

(i) The pair (f | , f | ) ∈ S d (M, V n ), with d = (n 2 +n 3 -2, n 1 +n 3 -2, n 1 + n 2 -2). (ii) For every β ∈ B n ( ˆ ) ∩ B n ( ˆ ), it holds λ β (f | ) = λ β (f | ). Proof (ii) ⇒ (i) Let (f, f ) ∈ V n be a pair of polynomials satisfying (ii). Then (11) (f | , f | ) = β∈Bn(M ) λ β (f | )β| ∪ + (h| , 0| ) + (0| , h | ), where h| = β∈Bn( ˆ )\Bn( ˆ ) λ β (f | )β| , and h | = β∈Bn( ˆ )\Bn( ˆ ) λ β (f | )β| .
Since β ∈ S d (M, V n ) for every β ∈ B n (M ), then the first term in Equation ( 11) is in S d (M, V n ). Now, if β ∈ B( ˆ ) \ B( ˆ ) then (β| , 0| ) = β| M , and similarly for β ∈ B( ˆ ) \ B( ˆ ). Hence, the last two terms in Equation [START_REF] De Boor | Approximation by smooth multivariate splines[END_REF] are also in S d (M, V n ), and (i) follows.

(i) ⇒ (ii) By Lemma 1 it is enough to consider (f | , 0| ) ∈ S d (M, V n ).
We need to show that λ β (f | ) = 0 for every β ∈ B( ˆ ) ∩ B( ˆ ). To simplify, by an abuse of notation we will denote Supp n ( ˆ ) ∩ Supp n ( ˆ ) simply by Supp n (ε 3 ).

Let ¯ ∩ ¯ = ε i ∈ E i ,

and denote by b i the barycentric coordinate relative to which vanishes at the edge ε

i . Let n = n 1 + n 2 + n 3 . Then (f | , 0| ) ∈ S d (M, V n ) if and only if the polynomial b n-ni-1 i divides f , see [2, Lemma 2.2].
Notice that the latter condition is satisfied if and only if the Bernstein-Bézier coefficients c jk of f | are zero for every 0 j n -n 1 -2.

Without loss of generality we can assume that i = 3 and is the triangle with vertices at (0, 0), (1, 0), and (1, 1). Thus, b 3 = x-y. We will prove that if

f ∈ V n and b n-n3-1 3 divides f then λ β (f | ) = 0 for every β ∈ B( ˆ )∩B( ˆ ).
We proceed by induction on n, with n i 1.

The induction base is n = 3, with n = (1, 1, 1), n 3 = 1 and n -n 3 -1 = 1. We have

B n ( ˆ ∪ ˆ ) = {β j } 3 j=0 ,
with β 0 = B n , and β i = B n (• + e i ) for i = 1, 2, 3. Thus, the only translates with support on are β 0 , β 2 and β 3 . Since Let us assume that the result is true for every m n for some n 3. Let n = (n 1 , n 2 , n 3 ) with n + 1 = n 1 + n 2 + n 3 . We consider three cases, in the first two cases only one of the indices n i is 2 and in the last case at least two of the indexes are 2.

β 0 = b 1 , β 2 = b 3 , and β 3 = b 2 , then f | = λ β0 (f | )b 1 | + λ β2 (f | )b 3 | + λ β3 (f | )b 2 | .
Case 1. Suppose n 1 2, and n 2 = n 3 = 1.

We have (x-y) n1 |f and D e1 f ≡ 0 mod (x-y) n1-1 . Notice that Equation ( 9) may be rewritten as ( 12)

f | = v∈Z 2 λ v (f )B n (• -v) ,
where λ v = 0 for every v / ∈ Supp n ( ), and

λ v (f ) = λ β (f | ) for β ∈ B( ˆ ) and β = B n (• -v).
Moreover, we know that for any linear combination of box splines [START_REF] De Boor | Box splines revisited: Convergence and acceleration methods for the subdivision and the cascade algorithms[END_REF] for a v ∈ R for every v ∈ Z 2 , and m = (n 1 -1, n 2 , n 3 ), see [START_REF] Lai | Spline functions on triangulations[END_REF]Lemma 12.3]. Thus,

D ei v∈Z 2 a v B n (• -v) = v∈Z 2 a v B m (• -v) -B m (• -v -e 1 ) = v∈Z 2 (a v+ei -a v )B m (• -v + e i ),
D e1 v∈Z 2 λ v (f )B n (• -v) = v∈Z 2 λ v+e1 (f ) -λ v (f ) B m (• -v + e 1 ) ≡ 0 mod (x -y) n1-1 ,
where m = (n 1 -1, 1, 1). By induction hypothesis λ v+e1 (f ) = λ v (f ) for every v ∈ Supp m (ε 3 ). Notice that, since n = (n 1 , 1, 1) then the elements v ∈ Supp v ( ˆ ) are either a multiple te 1 of e 1 or of the form e 2 + te 1 , for t ∈ Z. Thus, we have

f | = λ e1 (f ) t∈Z B n (• -te 1 )| + λ e2 (f ) t∈Z B n (• -e 2 + te 1 )| , for constants λ ei (f ) ∈ R, for i = 1, 2. Moreover, B n (• -te 1 )| ε1 = 0 for every t.
In particular, they are zero at the vertex (0, 0). But f (0, 0) = 0, and

1 = v∈Z 2 B n (• -v) = t∈Z B n (• -te 1 ) + B n (• -e 2 + te 1 ). Then f | (0,0) = λ e2 (f ) t∈Z B n (•-e 2 +te 1 )| (0,0) = λ e2 (f ) = 0.
We obtain λ e1 (f ) = 0 by considering the restriction at the vertex (1, 1).

Case 2. Suppose n 3 2 and n 1 = n 2 = 1.

By hypothesis (x -y)|f , and so D e3 f ≡ 0 mod (x -y). Following the same argument as above, λ v+e3 (f ) = λ v (f ) for every v ∈ Supp n (ε 3 ). Thus, there is a constant λ e3 which is equal to λ v (f ) for every v ∈ Supp n (ε 3 ), and

f | ε3 = λ e3 (f ) v∈Z 2 B n (• -v) ε3 .
Since f | ε3 = 0, and

v∈Z 2 B n (• -v) ε3 = 1, it follows λ e3 (f ) = 0.
Case 3. At least two of the indices n i are 2, say n 1 , n 2 2. By hypothesis (x -y) n1+n2-1 |f , and so D ei f ≡ 0 mod (x -y) n1+n2-2 for i = 1, 2. Similarly as before, rewrite f | as in Equation ( 12), and consider Equation ( 13) with i = 1. Thus,

D e1 v∈Z 2 λ v (f )B n (• -v) = v∈Z 2 λ v+e1 (f ) -λ v (f ) B m (• -v + e 1 )
≡ 0 mod (x -y) n1+n2-2 .

By induction hypothesis λ

v+e1 (f ) = λ v (f ) for every v ∈ Supp m (ε 3 ). Similarly, by considering D e2 f we get λ v+e2 (f ) = λ v (f ) for every v ∈ Supp m (ε 3 ), where m = (n 1 , n 2 -1, n 3 ). This implies that the coefficients λ v (f ) are equal to a constant λ, for every v ∈ Supp n (ε 3 ). Therefore, f | ε3 = λ v∈Z 2 B n (• -v)| ε3 . But f | ε3 = 0, and v∈Z 2 B n (• - v) = 1, in particular when taking the restriction to ε 3 . Hence, 0 = λ = λ β (f ), for every β ∈ B( ˆ ) ∩ B( ˆ ).
An analogous proof applies when the edge of intersection between ¯ and ¯ is parallel to one of the other two vectors e 1 or e 2 .

The edge-contact property plays a fundamental role for the construction of hierarchical spline spaces, we present that construction for any type-I box spline in Section 5. The following examples illustrate that importance of the properties of type-I box splines used to prove Lemma 2, we show that not every spline space possesses those properties.

Example 3 Let M = ˆ ∪ ˆ be the multicell domain in Figure 7. 7 The triangles and in G define the multicell domain M in Example 3, the triangles share the edge ε ∈ E2 that is parallel to the directional vector e2.

(-1, 0) (0, 0) (1, 1) (0, 1) ε △ △ ′ Fig.
Take the set B = {f 1 , . . . , f 4 } of spline functions f i on M defined as

f 1 = y 2 + 3x 2 + 4x, y 2 + x , f 2 = 4y 2 , 4y 2 , f 3 = x 2 , 0 , and f 4 = 0, x 2 .
Each pair of polynomials f i is defined on and , respectively. Then, B is contained in the space of splines S 0 (M, Q), where 0 = (0, 0, 0), and Q ⊆ P 2 is the linear subspace of polynomial spanned by x, x 2 , and y 2 .

We have

Q| = B| = y 2 + 3x 2 + 4x, 4y 2 , x 2
, and

Q| = B| = y 2 + x, 4y 2 , x 2 .
Then, g = 5 • (y 2 + 3x 2 + 4x) + 4y 2 , y 2 + x + 2 • (4y 2 ) ∈ S 0 (M, Q), but g is not an element in span (B). 3

The following example illustrates the importance of the spline space S d (M, V n ) (Definition 3) in the proof of Lemma 2.

Example 4 Let M = ˆ ∪ ˆ be the multicell domain defined from the triangles and such that ∩ = ε ∈ E 3 in Figure 8. Take n = (2, 1, 1), and consider the linear space V n = span B n ( ˆ ). The support and Bernstein coefficients of the box spline B (2,1,1) are displayed in Figure 8. The generators of V (2,1,1) can easily be described by restricting the translates

B (2,1,1) (• -v) ∈ B n (M ) to
and using the Bernstein coefficients of B (2,1,1) . Namely, V (2,1,1) is generated by

f 1 = 2(x -y)y + y 2 ; f 2 = 2(1 -x)y + y 2 ; f 3 = (x -y) 2 ; f 4 = (1 -x) 2 + 4(1 -x)(x -y) + (x -y) 2 + 2(1 -x)y + 2(x -y)y ; and f 5 = (1 -x) 2 .
In particular,

f = 4(x -y)y + 4(1 -x)(x -y) + (x -y) 2 = f 1 -f 2 + f 4 -f 5 ∈ V (2,1,1) . Let us notice that g = 0| , f | ∈ S 0 (M, V n ), but g / ∈ span B n (M ). In fact, for β(•) = B n (•) we have β| = f 1 , and λ β (f | ) = 1, but λ β (0| ) = 0.
On the other hand, ε ∈ E 3 i.e., ε is an edge parallel to the directional vector e 3 , and g is not a C 1 -continuous spline on M . Thus, g / ∈ S d (M, V n ) for d = (0, 1, 1).

3 In view of Lemma 2 we arrive to the following completeness result.

(0, 0) (1, 1) ε △

Corollary 1 (Completeness of type-I box splines) If M ⊆ G is an admissible domain, and n

= (n 1 , n 2 , n 3 ) ∈ Z 3 1 , then the generating set B n (M ) is complete for S d(n) (M, V n ), where d(n) = (n 2 + n 3 -2, n 1 + n 3 -2, n 1 + n 2 -2).

Hierarchical type-I box splines

In this section, to any given box spline B n we associate a space of special splines defined on a hierarchical grid. We construct a hierarchical basis for such space and prove that this basis is complete under certain assumptions on the domain hierarchy.

For an integer N 0, we recursively define a sequence of three-directional grids G 1 , G 2 , . . . , G N as follows. We take G 1 = G as the three-directional grid with vertices Z 2 introduced in Section 2, Figure 1. The hierarchical grids G for = 2, . . . , N are defined recursively, in such a way that G +1 is obtained from G by one global, uniform dyadic refinement step. More precisely, the grid G +1 = 1 2 G is the triangulation of R 2 with vertices at points 1/2 (k, k), obtained by drawing in the lines x = k/2 , y = k/2 , and x -y = k/2 , for all k ∈ Z. Thus, to construct G +1 , every triangle in the grid G is split into four smaller ones, as illustrated in Figure 9. The index will be called the level of the grid, and the number N specifies the number thereof. Each grid G is a uniform three-directional grid, similarly as for the grid G = G 1 before, we consider the triangles in this grid as open sets in R 2 .

G 1 G 2 G 3
Fig. 9 Three levels of three-directional hierarchical grids.

Let Ω be a domain of R 2 whose boundary ∂Ω is the union of edges from the grid G N . We define a hierarchical multicell domain H associated to the domain Ω as follows.

A nested sequence of subdomains of Ω is defined as a collection of domains

M such that ∅ = M 0 ⊆ M 1 ⊆ • • • ⊆ M N = Ω,
where M = ∈M , and M ⊆ G is a multicell domain in the grid G , for for each = 1, . . . , N . Thus, for each level , the boundary ∂M is a union of edges of the grid G . The difference between two successive subdomains, denoted D , is defined as the closure

D = M \ M -1 .
The associated refined domain of level D ⊆ G is defined as D = T (D ), where T (•) is the triangulation operator which restrict the grid G to a given subset of the plane R 2 . More precisely,

T (Q) = ˆ ∈ G : ⊂ Q .
The hierarchical multicell domain H associated to Ω is then the collection of triangles form all levels of the refinement area

H = N =1 D .
Using this notation, the domain Ω can be written as the union Ω = ∈H .

Definition 14 Let H be a three-directional hierarchical multicell domain associated to a domain Ω ⊆ R 2 , and let B n be a box spline, for some triple n = (n 1 , n 2 , n 3 ) ∈ Z 3 1 . We define P(H, V n ) as the set of piecewise polynomial functions on H associated to B n i.e.,

P(H, V n ) = f ∈ C 0 (Ω) : f | ∈ V n | for each triangle ∈ H . If d = (n 2 + n 3 -2, n 1 + n 3 -2, n 1 + n 2 -2)
, the hierarchical box spline space with edge smoothness d on H is defined as the set

(14) S d (H, V n ) = f ∈ P(H, V n ) : f | 3(ε) ∈ C di (3(ε) * ) for every ε ∈ E i ∩ H and i ∈ {1, 2, 3} .
For d as above, we define the linear space of hierarchical box splines with edge and vertex smoothness d on the hierarchical multicell domain H as the set

(15) Ŝd (H, V n ) = f ∈ P(H, V n ) : f | U ∈ D I (U, V n ) for , ∈ M ,
for some , ∩ = ∅ , U = ∪ , and I = i∈ST( , ) I d i .

In Equation ( 14), the diamond of and edge ε ∈ H is taken over the multicell domain M such that ε is an edge of the a triangle ∈ M . Namely,

3(ε) = ∈M ,ε∈ ˆ ˆ .
Let us fix n = (n 1 , n 2 , n 3 ) ∈ Z 3 1 and the corresponding box spline B n . For each level = 1, . . . , N we denote by B n the set of translates of B n with respect to the grid G . In particular, we have B 1 n = B n (G) as defined in Equation ( 4), and recursively

B +1 n = β(2• ) : β ∈ B n .
In an analogous way as we introduce B n (M ) in Definition 11 for a multicell domain M in the three-directional grid G, we now define the set B n (M ) of active box splines on the multicell domain M in G as follows,

B n (M ) = β| M : β ∈ B n , and supp β ∩ M = ∅ .
A set of linearly independent box splines on a hierarchical multicell domain H can be constructed by a selection procedure analogous to that proposed by Kraft in [START_REF] Kraft | Adaptive und linear unabhängige Multilevel B-Splines und ihre Anwendungen[END_REF] in the context of tensor-product B-splines.

For all levels , we select box splines translates, and define the sets K as follows

K = β ∈ B n (M ) : supp β ∩ M -1 = ∅ .
The collection of these box splines translates in the levels = 1, . . . , N forms a hierarchical box splines basis given by ( 16)

K = N =1 K .
The linear independence of the functions in K is implied by the local linear independence of the box splines at each level, see [START_REF] Kraft | Adaptive und linear unabhängige Multilevel B-Splines und ihre Anwendungen[END_REF].

Then the question of completeness of hierarchical type-I box spline spaces can be stated as follows. For a given hierarchical multicell domain H, does the basis K in Equation ( 16) span the hierarchical box spline space S d (H, V n ) defined in Equation ( 14)? Does K span Ŝd (H, V n ) defined in Equation ( 15)?

In the following theorem we address the first question, and provide a sufficient condition for the completeness of the hierarchical spline basis.

Theorem 1 Let H be a three-directional hierarchical multicell domain, and let B n be a box spline, for some triple n = (n 1 , n 2 , n 3 ) ∈ Z 3 1 . The basis K in Equation ( 16) spans the hierarchical box spline space S d (H, V n ) if each multicell domain M of H is admissible (Definition 13) with respect to the grid level .

Proof The proof follows standard arguments already presented in [START_REF] Giannelli | Bases and dimensions of bivariate hierarchical tensorproduct splines[END_REF][START_REF] Mokriš | On the completeness of hierarchical tensorproduct B-splines[END_REF] for the case of hierarchical tensor B-spline bases. Let Ω ⊂ R 2 be the domain associated to H. We prove by induction on the levels that every spline function s

∈ S d (H, V n ) admits a representation (17) s = (h 1 + • • • + h N )| Ω ,
where h ∈ span B n (M ), and

(18) h | M = s| M -(h 1 + . . . + h -1 )| M .
for = 1, . . . , N . For any given level , all functions h k | M of lower levels k < are contained in S d (M , V n ). This follows from the relation span

B k n (M k )| M ⊆ S d (M , V n ).
It follows that s| M ∈ S d (M , V n ). Consequently, the right-hand side of Equation ( 18) is contained in S d (M , V n ). Since the multicell domain M is admissible, we conclude that h ∈ span B n (M ) according to Corollary 1. In particular, choosing = N in Equation ( 18) we get Equation (17). Moreover, the construction of the functions h ensures that

h | M -1 = 0| M -1 .
Since the box splines possess the property of local linear independence we can conclude that h ∈ span K . This completes the proof.

Theorem 1 is a consequence of Corollary 1, and therefore of Lemma 2. This result generalizes the completeness property of the space of translates of the quartic box spline B 2 proved in [39, Theorem 26].

Concluding remarks

Remark 5 The contact characterization property proved in Lemma 2 implies the completeness of the space spanned by the translates of type-I box splines B n on a multicell domain M , with respect to the spline space S d(n) (M, V n ), where d(n) = (n 2 +n 3 -2, n 1 +n 3 -2, n 1 +n 2 -2), and n = (n 1 , n 2 , n 3 ). This result holds whenever M is an admissible domain i.e., whenever M does not have over-concave vertices nor kissing triangles, and supp (β) is connected for all translates β of B n which have support on M . This is a sufficient condition, but it will also be interesting to prove necessary conditions to achieve this completeness property of the type-I box spline space. In this direction, a complete characterization of the vertex-vertex contact plays a crucial role. This calls for exploring the algebraic formulation of super-smoothness at vertices in order to proving an analogous to Lemma 2 for vertex-vertex contact of type-I box splines.

Remark 6

We partially undertake the vertex-vertex contact question raised in the previous remark, for low degrees. In particular, we have verified algorithmically the vertex-vertex contact (that is, Lemma 2 for the case ¯ ∩ ¯ ∈ V and the space Ŝd (M, V n )) for all type-I box splines up to n = (4, 4, 4). The approach is based on the fact that both statements (i) and (ii) of Lemma 2 yield linear relations on the coefficients λ β (f | ) and λ β (f | ). We can express each of the two statements as matrices A (i) and A (ii) . The latter matrix has the form A (ii) = I p 0 q -P I p 0 q , [START_REF] Giannelli | An immersed-isogeometric model: application to linear elasticity and implementation with THBox-splines[END_REF] where P is a permutation matrix, p = #B n ( ˆ ) ∩ B n ( ˆ ), and q = φ(n) -p.

The matrix1 A (i) can be computed as a product A (i) = C v L M . The factor C v describes the Bézier continuity conditions at the common vertex v, corresponding to the smoothness type ST( , ) (Definition 5) and to the regularity vector d = (n 2 + n 3 -2, n 1 + n 3 -2, n 1 + n 2 -2). The right factor L M is the Bernstein-Bézier representation of the translates of B n , considered independently on each of the two triangles of M :

L M = L 0 0 L (20)
or, equivalently, a change of basis matrix of the space C -1 (M * ).

Showing that Lemma 2 holds in this setting (and for a fixed n) is done by showing that A (i) and A (ii) are equivalent matrices. Indeed, in all our computations we obtain A (ii) as the reduced row echelon form of A (i) for type-I box splines up to total degree 12. This verifies that the matrices are equivalent and so the vertex-vertex contact lemma holds in these cases.

Remark 7

In the discussion of completeness of hierarchical type-I box splines spaces there are two main differences to the original approach, which was formulated for tensor-product splines.

First, the translates of a box spline do not span the whole space of bivariate polynomials of a given total degree. For this reason, this special polynomial subspace had to be identified. In some sense this situation generalizes the tensor-product case, where the B-splines span a polynomial space of a given (coordinate-wise) bi-degree, instead of the the space of bivariate polynomials of a given total degree.

Second, the constraints on the domains are entirely different, due to the differences in the characterization of contacts between polynomial pieces. For bivariate tensor-product splines, both edge-edge and vertex-vertex contacts could be characterized easily by the equality of spline coefficients. In the case of type-I box splines, we proved the characterization solely for edge-edge contacts. Consequently, the completeness of hierarchical splines requires more severe restrictions to the hierarchical grid.

We have alleviated these extra restrictions by proving algorithmically the vertex-vertex contact for small polynomial degrees (total degree up to 12). By proving a characterization of the vertex-vertex contact of box splines for arbitrary degree, as we indicated in Remark 5, we could relax these restrictions on the hierarchical grids and construct more general hierarchical type-I spline space. We leave this as a future research direction.

Remark 2 From

 2 Definitions 3 and 6, it is clear that both Ŝd (M, V) and S d (M, V) are contained in the space of splines that are globally C r -continuous on M , where r = min{d 1 , d 2 , d 3 }. Moreover, they are both contained in D I (M, V) for I = I d 1 ∩ I d 2 ∩ I d 3 (see Definition 2), and Ŝd (M, V) ⊆ S d (M, V)

Fig. 5

 5 Fig. 5 Support of the Courant hat function B1, it corresponds to star (ν) where ν is the vertex (1, 1) (left), support of the box splines B (2,1,1) (center) and B (2,2,1) (right).

Since the polynomials b 1

 1 , b 2 , b 3 are linearly independent, and by hypothesis b 1 |f , it follows λ βi (f | ) = 0 for i = 2, 3. This proves the statement for the case n = 3.

Fig. 8

 8 Fig. 8 Bernstein coefficients of the box spline 2B (2,1,1) .

Observe that ker A (i) = Ŝd (M, Vn).


Separations Engineering: Advances in Adsorption

By Gennady Gor and Benoit Coasne

This themed issue of Current Opinion in Chemical Engineering is dedicated to recent advances in adsorptionespecially in the context of fields relevant to chemical engineering. On this occasion, a series of 10 outstanding contributions, written by experts from both the academic and industrial worlds, have been gathered around key questions connected to adsorption phenomena in membrane and separation sciences, electrochemistry and energy conversion/storage, poromechanics and civil engineering, etc. All these contributions are short and focused mini-reviews -typically, about 2500 words -on the latest experimental and/or theoretical developments in emerging adsorption topics. About half of these mini-reviews deal with current advances in adsorption in chemical engineering fields: Adsorption in Membrane Science, Coupling adsorption/mechanics; Adsorption and coadsorption processes at the industrial scale; Adsorption effects in electrochemistry; Adsorptiondriven heat pumps/chillers. The other half of these mini-reviews focus on advanced techniques to probe adsorption: Spectroscopic tools for adsorption; Imaging techniques to probe adsorption; Probing adsorption using Nuclear Magnetic Resonance; Kinetic approaches for molecular modeling of adsorption; Free energy techniques for simulation of adsorption. Each of these contributions focuses on some of the most important recently published work and emphasizes the significance of these advancements towards the overall progress of the field. The authors were also asked to focus on a few key advances based on the recent literature (within the past 2 to 5 years) rather than survey an entire field and to provide their critical assessment of these developments and future prospects. An additional key feature of the Current Opinion journals is to provide an annotation for each of the most interesting articles within the bibliography of the reviews.

Adsorption is a process of adhesion of atoms, molecules or ions in either gas or liquid phase, at the surface of solids or liquids. Most of the papers in this issue focus on physisorption, i.e. the process in which the adsorbed species do not form chemical bonds with the surface but only bind through physical forces (dispersion, electrostatic, etc.). Also the papers focus on adsorption on solid surfaces, and do not discuss adsorption on liquid phases. This however does not mean that chemisorption or adsorption on liquid surface is of lesser interest. The research of various aspects of adsorption is very broad, and this issue could have covered many more topics. We decided to limit them to those which have not been covered in recent review papers. There is a number of fields relevant to chemical engineering in which adsorption has recently emerged as a key ingredient. While in some applications adsorption effects have been recognized for a long time as important, there are over domains where they have been overlooked or only considered in a very effective fashion. This includes applications belonging to membrane science, poromechanics i.e. the field of mechanics of porous media which accounts for the coupling with the thermodynamics of confined fluids, gas or liquid adsorption in soft matter physics, phase transformation such as for heat pumps and chillers, etc. In fact, developments in these fields have often been done without specifically accounting for adsorption effects until recently. The present themed issue on Advances in Adsorption from Current Opinion in Chemical Engineering intends to highlight these recent developments. This is a timely issue as these efforts have been quite impressive in the last 10-20 years with significant breakthroughs reported -with or without direct technological transfer already made.

As indicated above, the mini-reviews in this themed issue subdivides into (1) five contributions focusing on advanced techniques to probe adsorption and (2) five contributions dealing with current advances in adsorption relevant to chemical engineering fields. In the two following paragraphs, we list and comment each of these minireviews.

Advanced techniques to probe adsorption • The paper by Ustinov presents an overview of the recent advances in development of kinetic Monte Carlo (kMC) method for application to modeling equilibrium phenomena [1]. As it is obvious from the name, the kMC method has been used for a number of years for modeling kinetics. However, in the recent years, it was shown to be a viable alternative to conventional (Metropolis) Monte Carlo to predict thermodynamic equilibrium, in particular adsorption equilibrium. This review outlines the main features of the method, and shows examples of application of this method for gas adsorption on solid surfaces. Importantly, the kMC method provides a straightforward way to calculate the Gibbs free energy, while performing simulations in the canonical ensemble • Another review focusing on free energy calculations -in the specific context of the flexibility of compliant MOF materials -is the work by Tanaka and Miyahara [2]. Such advanced techniques, which are being more and more frequently used in modern research, allow going much beyond standard Monte Carlo and molecular dynamics simulations. This review shows how, under given thermodynamic conditions, the stability and the molecular mechanism responsible for MOF flexibility can be assessed using different thermodynamic integration schemes (using either the adsorption isotherm or the hydrostatic pressure as data to be integrated). Beyond the fundamental insights gained into the different flexibility modes of this class of porous materials, such advanced molecular simulation techniques are receiving increased attention to screen and predict in silico different materials for gas storage and separation applications.

• The review by Pini and Joss [3] discusses how imaging techniques applied to adsorption in microporous materials provide a mean to "see the unseen". Recent progress made and on-going developments in this field aim at designing in situ techniques that probe transient and localized phenomena involved in porous adsorbents and catalysts. The review shows that, owing to the range of techniques available, involving microscopy, spectroscopy and tomography-based apparatuses, the rich and complex phenomena spanning several time and length scales can be monitored to gain insights into the adsorption but also the kinetics at play. With this respect, a strength of these methods, which can be applied to almost any system under static or kinetic conditions, is that they allow closing the gap by probing adsorption systems from the nanometer to the macroscopic scales. • The review by the Allgeier's group to some extent continues the possibilities to "see the unseen", discussing the application of nuclear-magnetic resonance (NMR) relaxometry for probing the porous media [4]. Starting from the physics of NMR, the paper describes how the measured relaxation times of nuclear spins can be related to the surface chemistry and porosity. This relation makes NMR relaxometry a promising tool for characterization of range of various geological porous media (such as shales), as well as materials used in catalysis and gas separations. • Tsyganenko et al discuss the process of controlled and selective desorption of molecules using resonance infrared laser [5]. The idea is based on direct excitation of vibrational states of adsorbed species. In this process only those adsorbed molecules whose vibrational frequency equals to the exciting radiation will be activated and desorbed from the surface. The review discusses the basic studies of these processes, highlighting potential applications for isotope separations, which have not been succeeded yet.

Current advances in adsorption relevant to chemical engineering fields • Two reviews discuss the phenomena related to the coupling between adsorption and mechanics, approaching that from totally different perspectives, and discussing literature in different fields. The paper of Vandamme introduces the phenomenon of adsorptioninduced deformation -appearance of stresses and strains due to adsorption of fluids in porous materials [6]. It proceeds with discussion of various theoretical approaches to predict those effects. Then the reciprocal effects are discussed: when materials are compliant, how this compliancy affects their adsorption properties. The review concludes with the discussion of specific challenges in theories of mechanics of adsorption: effects of pore morphology, pore-size distribution, anisotropy and non-linear responses. • While also discussing the mechanicsadsorption coupling, the review by Weissmüller [7] focuses mainly on electrosorption phenomenon. In this case instead of a gas or a liquid, described by its chemical potential, the adsorbing species is an electrolyte, described by the electric potential. The paper starts from the fundamental thermodynamic relations of the surface mechanics of electrodes, and proceeds to the overview of recent experimental measurements of electro-chemomechanical coupling, relating them to the calculations based on the density functional theory. The review also demonstrates the similarity between the mechanical effects of electrosorption and adsorption. • The review by Julbe and coworkers focuses on the role of adsorption phenomena and surface reactions in the efficiency and control of membrane-based processes [8]. Such adsorption effects in this context are discussed by considering both gas-phase and liquid-phase separations for membranes ranging from inorganic to composite (i.e. organic/inorganic) membranes. In addition to reviewing the importance of adsorption for the design of separation processes using membranes, their contribution discusses other adsorption-related issues to be further explored: fouling phenomena, novel surface treatments and the possibility to use adsorption as a tool to characterize the microstructure of membranes and its evolution when used in real process conditions. • The review by Pullumbi and coworkers deals with the important topic of gas separation using adsorption-based processes [9]. While these technologies are already in use for many processes such as purification, separation, etc., the authors have identified opportunities to further improve their efficiency and expand the scope of their applications. In particular, the role of numerical simulation aspects in boosting technological transfer towards the industrial world is highlighted. This includes research in the fields of multiscale modeling and simulation tools, process design and complex optimization, synthesis scale up/processing together with materials characterization but also knowledge discovery in already existing databases • The review by Hastürk et al. focuses on the design of optimal adsorbents for adsorption heat transformation applications -either for cooling or heating processes. This important field of research keeps receiving very broad attention as it is a very promising environmentally friendly application [10]. In this contribution, after presenting the fundamentals of the thermodynamics of adsorption heat transformation, the authors adopt a material scientist point of view by considering different types of adsorbents. This includes samples as different as activated carbons, silica gels, zeolites and metal organic frameworks. Promising materials such as composite materials made of zeolites or silica gel with salts like LiCl or CaCl2 are also envisaged owing to their large uptakes and transformation enthalpies.


Mingyang Chen

Benoit Coasne
email: benoit.coasne@univ-grenoble-alpes

Robert Guyer

Dominique Derome

Jan Carmeliet

Benoit Coasne

Molecular Simulation of Sorption-Induced Deformation in Atomistic Nanoporous Materials

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1.

Introduction

When using superscirpt referernces, there is no space betweeb the reference and the word before Sorption-induced deformation is a common phenomenon in various types of porous materials [START_REF] Gor | Adsorption-Induced Deformation of Nanoporous Materials-A Review[END_REF] . This phenomenon arises at the molecular level because of the strong intermolecular forces at the sorption surface induced by the adsorbent-adsorbate interactions [START_REF] Suzuki | Grand Canonical Monte Carlo Simulation for Nitrogen Adsorption in Graphitic Slit Micropores: Effect of Interlayer Distance[END_REF] . Specifically, for nanoporous materials, experimental and theoretical studies find sorption-induced deformation generally associated with first order phase transitions which correspond to capillary condensation and evaporation [START_REF] Neimark | Capillary Condensation in MMS and Pore Structure Characterization[END_REF] . As a result, the capillary hysteresis usually observed in sorption isotherms is accompanied by hysteresis in strain isotherms [START_REF] Shao | Different Measures for the Capillarity-Driven Deformation of a Nanoporous Metal[END_REF] . Moreover, it is found that the deformation has distinct characteristics before and after phase transition [START_REF] Prass | Capillarity-Driven Deformation of Ordered Nanoporous Silica[END_REF][START_REF] Balzer | Relationship between Pore Structure and Sorption-Induced Deformation in Hierarchical Silica-Based Monoliths[END_REF] , which indicates that the sorption-induced deformation is governed by different mechanisms at different adsorption stages.

When you have to rferences; they should be separated by commas. For instance 9,10 not 910. From a mechanical perspective, sorption-induced deformation of nanoporous materials is generally attributed to either disjoining pressure [START_REF] Gor | Adsorption of N-Pentane on Mesoporous Silica and Adsorbent Deformation[END_REF] , surface stress 910 or Laplace pressure [START_REF] Prass | Capillarity-Driven Deformation of Ordered Nanoporous Silica[END_REF] effects. In the formulat it should be W(e) not W alone. Also, w should be w(e) etc. The disjoining pressure is defined as , where W(e) is the surface potential that stabilizes the adsorbed film of a thickness e. Note that this definition indicates that disjoining pressure acts in the direction perpendicular to the surface. However, it must be recalled that the surface potential term scales with the film surface area (which pertains to the direction parallel to the pore surface). The concept of disjoining pressure is often used to explain sorption-induced deformation. Gor and Neimark [START_REF] Gor | Adsorption-Induced Deformation of Mesoporous Solids: Macroscopic Approach and Density Functional Theory[END_REF] use the disjoining pressure concept as formalized in the Derjaguin-Broekhoff-de Boer (DBdB) theory to explain the deformation of cylindrical pores during sorption before the occurrence of capillary condensation. A similar approach is adopted by Kolesnikov et al [START_REF] Kolesnikov | Effects of Enhanced Flexibility and Pore Size Distribution on Adsorption-Induced Deformation of Mesoporous Materials[END_REF] when considering sorption in a pore network. The surface free energy of the liquid-vapor interface is written as dependent on the liquid film thickness through the use of the disjoining pressure term.

On similar systems, the concept of disjoining pressure is also used in computational modeling work using methods such as Grand Canonical Monte Carlo (GCMC) and density =d / d We  functional theory DFT. For example, conducting GCMC simulations, Cornette et al [START_REF] Cornette | Binary Gas Mixture Adsorption-Induced Deformation of Microporous Carbons by Monte Carlo Simulation[END_REF] find that adsorption-induced deformation in a slit pore is driven by the solvation pressure which is equivalent to the disjoining pressure. Similar observation is found in Grégoire et al. who used DFT calculations to describe the mechanics and physics of nanoconfined fluids [START_REF] Grégoire | Estimation of Adsorption-Induced Pore Pressure and Confinement in a Nanoscopic Slit Pore by a Density Functional Theory[END_REF] .

A second explanation of sorption-induced deformation stems from the work of Bangham [START_REF] Bangham | The Expansion of Charcoal Accompanying Sorption of Gases and Vapours[END_REF] which was extended by Shuttleworth [START_REF] Shuttleworth | The Surface Tension of Solids[END_REF] . In this approach, the surface stress can be expressed according to Shuttleworth equation as , where is the surface energy and is the surface strain. Surface stress arises from the spreading effect of free energy of the solid-liquid interface and the dependence on strain. Unlike disjoining pressure, surface stress acts within the interface and drives the deformation along the inplane direction of the surface. Zhang [START_REF] Zhang | Mechanics of Adsorption-deformation Coupling in Porous Media[END_REF] showed that the swelling of coal subjected to pure CO2 and mixtures of CO2/CH4 injections can be modeled assuming that all deformation is due solely to surface stress. Furthermore, as indicated by Gor et al. [START_REF] Gor | Revisiting Bangham's Law of Adsorption-Induced Deformation: Changes of Surface Energy and Surface Stress[END_REF] , the surface stress can be either positive or negative so that both swelling and shrinkage of the solid can be observed. A third explanation of sorption-induced deformation relies on the Laplace pressure, which characterizes the pressure difference across a curved boundary between a vapor and a liquid. The Laplace pressure is generally invoked to explain the deformation occurring after capillary condensation and in the presence of a liquid meniscus at the open end of the pore.

The three different mechanisms mentioned above can have combined effects on the deformation. As a result, a comprehensive description of adsorption-induced deformation should address all these mechanisms. A reasonable and straightforward approach to conduct such analysis is to use a simple pore geometry such as the slit pore model. The major advantage of this pore model over other geometries is that one can easily separate contributions from the directions perpendicular and parallel to the adsorption surface (which involve disjoining pressure and surface stress, respectively). Studies of sorption behavior using the slit pore model date back to the 1950s when Steele and Halsey 18

A S  d / d A A A S S S A      A S  A 
modified the theory of the interaction of gas molecules with a structureless plane to treat capillary space. In the same spirit, Stoeckli [START_REF] Stoeckli | The Gas-Solid Interface Calculations of Adsorption Potentials in Slot-Like Pores of Molecular Dimensions[END_REF] calculated the adsorption potential for a slit pore. Derjaguin and Churaev 20 addressed capillary condensation in slit pores with the concept of disjoining pressure. However, disjoining pressure is here used solely to model phase transition, while the related mechanical effect is not considered. In addition to such analytical studies, sorption behavior is also extensively studied using numerical simulations. These studies [START_REF] Coasne | Multiscale Adsorption and Transport in Hierarchical Porous Materials[END_REF][START_REF] Balbuena | Classification of Adsorption Behavior: Simple Fluids in Pores of Slit-Shaped Geometry[END_REF][START_REF] Suzuki | Grand Canonical Monte Carlo Simulation for Nitrogen Adsorption in Graphitic Slit Micropores: Effect of Interlayer Distance[END_REF][START_REF] Jagiello | 2D-NLDFT Adsorption Models for Carbon Slit-Shaped Pores with Surface Energetical Heterogeneity and Geometrical Corrugation[END_REF][START_REF] Striolo | Water Adsorption in Carbon-Slit Nanopores[END_REF] show that the hysteresis observed in the sorption isotherm strongly depends on the pore width and surface energy heterogeneity. Few studies focus on the sorption-induced deformation in slit pores by looking at the in-pore pressure. It was shown by Long et al [START_REF] Long | Pressure Enhancement in Carbon Nanopores: A Major Confinement Effect[END_REF][START_REF] Long | Under Pressure: Quasi-High Pressure Effects in Nanopores[END_REF][START_REF] Long | On the Molecular Origin of High-Pressure Effects in Nanoconfinement: The Role of Surface Chemistry and Roughness[END_REF] with GCMC method that the in-pore pressure is strongly enhanced by confinement of nanopores. As the result, the in-pore pressure, in both normal and tangential directions, can be several orders of magnitude larger than the bulk phase pressure in equilibrium. This high in-pore pressure indicates that deformation induced by sorption in nanoporous material can be significant, and the high sensitivity of the in-pore pressure over the bulk pressure shows that the variation of deformation can be large during the sorption process. Moreover, the results by Long et al also show that the pressure tensor differs in normal direction and longitudinal direction, indicating an anisotropic deformation pattern. Similar high in-pore pressure and related deformation are also given by other works employing either GCMC method 2930 or density functional theory 14 31 .

However, a great limitation of works stated above is that the pore systems built are not deformable, thus one can only deduce the deformation by assuming a mechanical constitutive relationship (Hooke's law for instance). A deformable slit pore model was considered by Schoen et al [START_REF] Günther | Sorption Strain as a Packing Phenomenon[END_REF][START_REF] Günther | Sorption Strains and Their Consequences for Capillary Condensation in Nanoconfinement[END_REF][START_REF] Schoen | Phase Transitions in Nanoconfined Fluids: Synergistic Coupling between Soft and Hard Matter[END_REF][START_REF] Schoen | Capillary Condensation in Deformable Mesopores: Wetting versus Nanomechanics[END_REF] where the porosity is located between two Einstein solid substrates (i.e. each atom is attached to its original position through a harmonic potential).

While giving good results for the normal strain, Einstein solids cannot model the deformation along the longitudinal direction properly. Indeed, as each solid atom is connected to a fixed point instead of neighboring atoms, the longitudinal stress depends on the absolute displacement instead of the relative displacement between neighboring atoms so non-physical strain fields are obtained.

While many experimental and theoretical studies have been conducted on the sorptioninduced deformation of nanoporous materials, most research has been limited to a single deformation mechanism out of the three at play in real materials. In this paper, we aim at characterizing the contributions of these different mechanisms in sorption-induced deformation with the help of a realistic slit pore model described at the molecular scale.

The hysteresis observed in sorption and deformation at the single pore level is also discussed.

Models and Methods

If for langmuir (or equivalent) you need to have subsections numbering as well

Preparation of the solid model

In this paper, we build a slit pore model composed of two atomic solid substrates (see Figure 1). Each substrate, which is made up of two layers of solid atoms, is first generated in a simple cubic crystal with a lattice parameter a = 0.3822 nm. To avoid crystallization of the confined liquid that would be induced by the use of a crystalline solid surface, amorphization of the subtrates is achieved by applying a random spatial displacement in the x-y plane to each solid atom. Note that the solid surface remains plane as no extra roughness is introduced. Then, a harmonic bond is added between each pair of nearest neighboring solid atoms to mimic a covalent chemical bond as in real solids. To further reduce the problem complexity, we classify all the bond lengths modified into 20 length bins and replace each with the average length of its corresponding bin. This allows describing the system by specifying just 20 different harmonic bonds instead of specifying parameters for every bond length. However, the replacement introduces pre-stresses into the substrates which need to be relaxed before sorption simulation. Specifically, the solid model is relaxed at constant temperature 87 K using a Molecular Dynamics (MD) run in the NVT ensemble. Using the strategy above, one can produce any desired stiffness of the substrate by tuning the bond strength. In this work, the stiffness of the solid substrates is set to 10 GPa.

As for the longitudinal direction, the substrates are continuous solid materials with a stiffness which solely depends on the bond strengths. I don't undertstand the previous sample: why do you say that it s continuous? If it s made of atoms, it s not really continuuous. As for the normal direction, the stiffness arises from the surrounding nanoporous void which provides a certain stiffness during deformation. As the two substrates of the slit-pore model are separated by the pore space, they can move freely in the normal direction so that we add a linear spring of a force constant k connecting the two solid substrates to account for the stiffness in the normal direction. Specifically, the spring connects the centers of mass (COM) of the two outer layers of solid atoms and only accounts for the relative displacement in the normal direction. One advantage in doing this is that one can reproduce any target modulus in the normal direction by tuning the spring constant k. In fact, k can be related to Young's modulus in the normal direction Ezz:

(0)
where Lx, Ly, Lz are the sizes of the solid model in the x, y, z directions, respectively.

Simulation details use subsection numbering

As both sorption and deformation arise at the molecular scale, we use hybrid molecular simulations including both Grand Canonical Monte Carlo (GCMC) and Molecular Dynamics (MD). The GCMC method, which is often used to determine sorption isotherms in various materials, fails to model the sorption-induced deformation directly as it relies on the Grand Canonical ensemble which imposes a constant volume (and therefore does not allow deformation). On the other hand, MD allows the deformation of the material but cannot model the sorption process. Thus, we employ a hybrid GCMC/MD method, which runs iterations of GCMC and MD simulations until an equilibrium state is reached. Specifically, each block of GCMC and MD includes 1000 trial insertion/deletion attempts of the fluid molecule and 1000 steps of Verlet integration of MD. In this work, 10 5 blocks are carried out to reach equilibrium.

Argon is taken as a simple prototypical adsorbate model and is described using the 12-6 Lennard-Jones (LJ) potential with = 0.3405 nm and = 119.8 K (kB is the Boltzmann constant). The solid-fluid interaction is also modeled using the LJ potential with = 0.45 nm and . The mechanical properties of the solid are solely controlled by the harmonic bonds and the spring described above. The sorption process is carried out at a constant temperature T = 87 K. The simulation box boundaries are made periodic in all the three directions of space. In the y and z directions, some void is left between the substrates and the box boundary to allow the deformation. However, this additional space is not accessible for the argon atoms so that it cannot be considered that the stress imposed in this direction corresponds to the argon gas pressure imposed through the use of constant chemical potential simulations. I modified the previous sentence, please check that it is ok.

To prevent any tangential drift of the upper and lower solid plates during the simulation, the center-of-mass position of each solid plate is constrained along the tangential direction (y and z direction) by adjusting the coordinates of the atoms every timestep while maintaining relative coordinates of any pair of atoms in the solid plates.

Results and discussion

Sorption and strain isotherms

The sorption isotherms are shown in Figure 2 Sorption results in deformations in the normal and longitudinal directions as shown in Figure 2 (a) (center and bottom). In the normal direction, no significant deformation is observed in the film region during adsorption until condensation occurs where an abrupt shrinkage occurs along the phase transition. In the filled region, as the pressure keeps increasing, the shrinkage reverts back to zero strain as the menisci flatten. Upon desorption, shrinkage increases in a reversible fashion as the relative pressure decreases until evaporation occurs. Afterwards, as the pressure keeps decreasing, the strain goes back to zero in the film region. In the longitudinal direction, shrinkage is observed even at very low relative pressures in the film region and reaches a maximum at a point around 0.15 p/p0. Then, shrinkage decreases as the relative pressure keeps increasing in the film region (the origin of the maximum in the longitudinal shrinkage will be discussed below). A small yet abrupt shrinkage is observed at the condensation point. In the filled region, the longitudinal strain increases back to zero in a fashion similar to that observed for the normal strain. Upon desorption, the longitudinal strain decreases beyond the condensation point and jumps at the evaporation point back to the strain observed upon adsorption in the film region.

Altogether, the results above show that the system deforms quite differently in the normal and longitudinal directions in the film region while both normal and longitudinal strains follow the same trend in the filled region. A pronounced hysteresis is observed in both the normal strain and longitudinal strain isotherms. Obviously, this hysteresis is concomittant with the hysteresis observed in the sorption isotherm. However, no hysteresis is found when plotting the strain as a function of the adsorbed amount in Figure 2(c), therefore

showing that the strain can be unambiguously defined for a given fluid configuration.

Deformation mechanisms

For the normal direction, the deformation is driven by the internal normal pressure of the fluid, pzz. In the film region, the two films are separated by a vapor phase which is in contact with the vapor reservoir located outside the porosity. Considering the force balance in the z-direction, we get pzz ~ p where p is the vapor pressure in the reservoir. Thus, the maximum vapor pressure that can be achieved in the film region is around 6000 Pa, which corresponds to the condensation pressure. Given the modulus Ezz = 5 GPa used to maintain the two pore surfaces at a distance corresponding to the nominal pore size (see the definition of the spring used in the section method), the expected strain in the film region please use scientific notation for Langmuir not engineering notation:

1E-6 should be 10 -6 is negligible in agreement with the values obtained in our molecular simulations. However, in the filled region, the pore space is filled with a closely packed liquid phase so that two mechanical effects arise: (1) the meniscus formed at the open ends imposes a Laplace pressure on the liquid and ( 2) the packing behavior of the atoms in the liquid phase can introduce a disjoining pressure term. The pressure due to the meniscus is given by the Laplace pressure:

(0

)
where Rg is the ideal gas constant, VL is the molar volume of the condensed phase, T is the temperature, pV and pL are the pressures of the vapor and liquid phases, respectively. Under the condition of linear elasticity and small deformation, the strain in the normal direction caused by the Laplace pressure can be calculated as:

(0)
As shown in Figure 2 (a), this estimated is in reasonable agreement with the simulated normal strain but with a constant offset. This indicates that the variation of strain with the relative pressure in the filled region is mainly driven by the variation of the Laplace pressure. The constant offset in Figure 2 (a), which corresponds to an additional term arising from the liquid pressure, is thought to be due to the disjoining pressure. Generally the disjoining pressure is strongly pore size dependent and decays while oscillating upon increasing the pore size. It indicates that in our simulation the disjoining pressure is negative so that further shrinkage with respect to the Laplace pressure effect applies. Note / 1.2E 6

zz pE     0 ln g V VL L RT p pp Vp  ln g L V L V zz zz L zz SAT RT p p p E V E p    L zz
 that this disjoining pressure only persists in the filled region as the surfaces see each other through the liquid phase. In contrast, in the film region, the two films are separated by the vapor phase so that they do not interact with each other and the disjoining pressure becomes negligible (we note that the pore size is larger than the interaction cutoff used in our simulations so that no disjoining effect is expected in the gas phase).

For the longitudinal direction, significant deformation is observed in the film region. This is generally explained by the surface stress. However, a detailed analysis of the strain isotherm at very low relative pressure indicates that shrinkage starts at the very beginning of the adsorption processi.e. even before a single-layer film is formed. As shown in Figure 3a, the strain decreases (increasing shrinkage) with the adsorbed amount with a steep initial slope, therefore indicating that the effect of scattered fluid atoms on the deformation of the substrates can be very significant. To further illustrate the longitudinal shrinkage observed in the film region, we consider the simple following model where adsorbed atoms interact with solid substrate as shown in Figure 4(a). We use in this model the same solid-liquid Lennard-Jones interaction potential used in the molecular simulations. The initial distance between nearest neighbor solid atoms is taken to be 0.38 nm. Because of the additivity of the energy for liquid-solid pair interaction potentials, the energy of the system can be expressed as the sum of solid-solid, solid-liquid and liquid-liquid energy terms:

Here the solid-solid, solid-liquid and liquid-liquid energy are USS, USLand ULL. The solid is assumed to behave in a linear elastic fashion so that USS is expressed as a quadratic strain term, , where C is the elastic constant. Both USL and ULL are expressed in terms of a LJ potential. The distance between the i-th solid atom and the liquid atom is Ri while the distance between the neighboring solid atoms is a. For the sake of simplicity, the distance between the neighboring gas atoms is written a where  is the distance ratio between gas atom distance and solid atom distance. By minimizing U, one can get at equilibrium.

Let us now consider the very simple situation where only one liquid atom is adsorbed at the solid surface. As shown in Figure 4 (b), the solid-liquid interaction energy (black curve) increases with the strain, which indicates that the system can reach lower energy level upon shrinking. Now if we take into account the contribution of the elastic energy (blue curve), we find an energy minimum at a negative longitudinal strain, which indicates that the solid undergoes shrinkage because of the solid-fluid interaction. Another important factor is the 2D internal film pressure which arises from the fluid-fluid interaction energy Uff. At very low relative pressure, the argon atoms are far from each other which results in small attractive inter-molecular forces and therefore a relatively small negative pressure contribution. Considering the force balance in the longitudinal direction, the negative fluid pressure will cause additional shrinkage of the solid substrate as indicated by the solid red However, if the film density keeps increasing, the attractive inter-molecular force decreases and eventually becomes repulsive, resulting into a swelling strain of the solid substrate (= 1, dotted red line in Figure 4(b)).

yy i ii U U U U C R R a a                   2
In the filled region, the longitudinal strain is found to follow a similar trend with that of the normal direction: abrupt shrinkage at condensation and evaporation points. The magnitude of shrinkage decreases as the relative pressure increases, which indicates that the Laplace pressure also takes effect in the longitudinal direction (Eq. ( 5) also applies by substituting

Ezz with Eyy).

Isosteric enthalpy of adsorption Again use subsection numbering

The sorption-induced deformation of the slit pore can also be investigated by checking the isosteric enthalpy of adsorption of the system along the adsorption/strain isotherm:

(0

)
where E is the potential energy of the system and N is the adsorbed amount. Because of the additivity of the energy, the isosteric heat of adsorption can be split into three contributions:

The The deformation of the system can be checked by looking at since it characterizes the change of deformation energy with adsorbed amount. According to Figure 5, is relatively small compared to the other two contributions, which indicates that the effect of deformation on the sorption process is not significant. On the other hand, we find a larger

() 1 () SL LL SS st st st st q E N E N q q q              SL st q LL st q SS st q SL st q LL st q SS st q SS st q contribution of
in the filled region as shown in the inset of Figure 5. This can be attributed to the larger dependence of the strain on the adsorbed amount (larger /N) .

Considering the linear elasticity of the solid substrates, should satisfy (0)

As shown in Figure 2 (c), the slope is steeper in the filled region which results in a larger .

SS SS ii st ii y z U q V E NN          SS st q st q SL st q LL st q SS st q 2
respectively. The inset shows a zoom-in the data for . The open and closed symbols represent the adsorption and desorption data, respectively.

Effect of temperature and solid/fluid interaction strength

We further discuss the influence of two factors on the sorption-induced deformation: the solid-fluid interaction strength and the temperature. The solid-fluid interaction strength is tuned by changing the solid-fluid bonding energy in the LJ potential. As shown in Figure 6(a), the solid-fluid interaction does not strongly affect the normal strain. This can be attributed to the fact the normal strain is mainly governed by the Laplace pressure, which is independent of solid-liquid interaction as a first order approximation. On the contrary, the solid-liquid interaction significantly influences the longitudinal strain in the film region.

More in detail, a stronger solid-liquid strength leads to larger shrinkage magnitude. As the shrinkage is caused by the attractive force between the solid and liquid atoms, a stronger solid-liquid strength induces larger attractive forces, resulting in larger shrinkages in the film region. On the other hand, as the longitudinal strain is governed by Laplace pressure in the filled region, the influence of the solid-fluid interaction becomes negligible in the bulk phase region. The temperature plays an important role on the normal strain by influencing the Laplace pressure (Figure 6b). We find that the simulated strain is consistent with the values predicted from Eq. (3). For the longitudinal strain, the temperature does not have a significant influence in the film region. Although the change of temperature influences the liquid pressure, the strain is mainly governed by the solid-fluid interaction and the virial part of the liquid pressure. On the other hand, the influence of the temperature on the longitudinal strain in the bulk phase region is similar to that on the normal strain as they are both governed by the Laplace pressure. predicted by Eq.( 5) for T = 100 K and 87 K, respectively.

Conclusion

In this work, the sorption-induced deformation in an atomistic slit pore is simulated. We find that the sorption-induced deformation is governed by different mechanisms depending on the adsorption region considered. In the liquid-filled region, the Laplace pressure is dominating both the longitudinal and normal strains. In contrast, in the film region, the normal strain can be neglected while the longitudinal strain is driven by both the solidliquid attractive forces and the 2D pressure in the liquid film.

Figure 1 .

 1 Figure 1. Deformable atomistic slit pore model used in the present work. The slit pore model (right) is assumed to be taken from a linear elastic porous medium containing slit pores (left). The mechanical effect of surrounding material is simply described using a spring connecting the two substrates. The longitudinal and normal directions are defined as the y and z directions, respectively.

  (a,top). Starting from an empty pore, upon increasing the pressure, adsorption starts with a film forming at the solid wall at low relative pressure (Figure2b(I)) -this domain is referred to as "film region" in what follows. Capillary condensation (Figure 2 b(II)) occurs at a certain chemical potential at which the vapor phase between the liquid films is replaced by a liquid phase forming menisci -this domain is referred to as "filled region". Afterwards, the adsorption amount increases gradually as the menisci at the two open ends flatten. On the other hand, desorption starts with the retreating of menisci (Figure 2 b(III)) until the fluid separates at a certain relative pressure and only adsorbed films are left on the solid walls. This process is referred to as ff evaporation. The difference in the chemical potential at which capillary condensation and evaporation occur results in a hysteresis loop in the sorption isotherm. The simulated sorption behavior is consistent with experiments and theories conducted by other researchers 3637 .

Figure 2 .

 2 Figure 2. (a) Sorption isotherms (top), normal strain (center) and longitudinal strain (bottom) as observed from molecular simulation. The adsorbed amount is given in pore density n which is defined as n = N/V, where N and V are the number of adsorbed molecules and pore volume, respectively. The open and closed symbols represent the adsorption and desorption branches, respectively. The dashed lines are the strain isotherms in the filled region predicted using Eq. (3). (b) Fluid density map corresponding to points I, II, III in (a); (c) Strains zz and yy as a function of adsorbed amount n. The black and red symbols correspond to the normal strain zz and longitudinal strain yy, respectively. The open and closed symbols represent the adsorption and desorption branches, respectively.

  Figure 3b illustrates the film formation process of the first layer, showing the densification of the adsorbed atoms at the solid surface. More specifically, the adsorbed fluid atoms make the solid shrink to maximize the density of fluid-solid interactions.

Figure 3 .

 3 Figure 3. (a) Longitudinal strain as a function of pore fluid density before the completion of first adsorbed layer. (b) Typical molecular configuration showing the adsorbed film at points (I) (II) (III) in (a).

  5) in Figure4(b). As the relative pressure increases, the adsorbed argon film becomes denser ( decreases) and the amplitude of the attractive inter-molecular forces becomes larger; as a result, the liquid pressure become more negative which leads to a larger shrinkage strain in the solid substrate (dashed red line, = 3 in Figure4 (b)).

Figure 4 .

 4 Figure 4. (a) Schematic view of a system consisting of a solid substrate and adsorbed gas atoms. G and S represent gas and solid atoms respectively. The solid-solid, solid-gas, gasgas distances are a, Ri, and a. (b) Relationship between the system energies and the tangential strain at different film density. The tangential strain at equilibrium corresponds to the minimal energy at each density. As the film density increases, the strain can turn

  solid-fluid , fluid-fluid and solid-solid contributions. There s something really wrong here. You can simplify by qst = -<E> which is obviously not correct. Is this what you did? You need to use the following formula: According to Fig. 5 decreases and increases with increasing the adsorbed amount until a plateau is reached. These results are consistent with the data reported by other researchers. This can be attributed to the fact that at low adsorbed amount the enthalpy change associated with adsorption mainly arises from the contribution of solid-fluid interaction as the fluid molecules are adsorbed at the solid substrate. On the other hand, at larger adsorbed amounts, the enthalpy change mainly arises from the fluid-fluid interaction.

Figure 5 .

 5 Figure 5. Different contributions to the isosteric enthalpy of adsorption as a function of the adsorbed amount n. The black, purple, blue, red symbols represent , , ,

Figure 6 .

 6 Figure 6. Effect of different factors on the strain isotherm. (a) Effect of solid-fluid interaction strength: the red, blue and black symbols represent different interaction strengths: εsf = εff, 1.5 εsf and 2.0 εsf, respectively. The open and closed symbols represent the adsorption and desorption data, respectively. The black solid line is the strain isotherm as predicted by Eq. (5). (b) Effect of temperature: the red, blue and black symbols represent data for T = 100 K 92.5 K and 87 K. The open and closed symbols represent the adsorption and desorption data, respectively. The red and black solid lines are the strain isotherms as

Acknowledgements

The authors acknowledge the support of the Swiss National Science Foundation (SNSF) (No.143601).